title: Deep (C)RBM - How to generate data?

text: I want to create a CRBM that learns sequences of data. Since the data can be quite complex, I want to use several layers, but I am struggling how to use several layers. I found this article about Deep RBM. It is about how to train several layers. The problem I have with this paper is in equation 19: It says that h1 depends on v and h2, but when I just start giving data into the RBM, there are no values for h2, so how do I generate h1 then? Also, how does this algorithm is used with the connections used for history values like there are in CRBMs? If nobody knows how to answer the question relating to the linked article, I would be happy if someone could show me how deep CRBMs are being trained and used, my one-layer CRBM works fine.


Can a moving average value result be higher than the highest constituent value on which the moving average is based? I realize that somewhere here there maybe an answer, but I have been unable to find one specific to this question. I do not believe that a "moving average value result" based on a set of values over the moving time period specified (such as the previous 4 weeks sales) can be higher than the individual highest constituent sales total that is one of the values considered in the moving average calculation. Am I right here? Thank you.


title: Appropriate statistical measure for rank correlation

text: In the course of my work, something of a problem has arisen. It can be summarized like this: there's a large pool of job applicants, let's say N of them. What needs to be done is that these applicants be ranked in order of desirability to the company: 1, 2, 3,... N. Until now, this has been done rather laboriously by discussions (read: arguments) between multiple people responsible for the decision. The idea is NOT to eliminate said discussions, but to reduce the amount of work necessary. So, the new idea being batted about is that one could introduce a sort of "score" for each applicant... For the sake of argument, let's say that each applicant was rated from a low of 1 to a high of 5 for each of several different dimensions (e.g., raw intelligence, sociability, etc.), and the ratings were simply added up to generate a final score. Then the applicants could be ordered by score. Of course, there will be many many ties with this particular system, but it would reduce the problem of ranking each applicant by creating "bins" that could then be sorted through individually, i.e., independent of those outside a given bin. But how to convince people that the "score" had any correlation with rankings as done by the old, laborious method? It seems to me that one of the common nonparametric rank correlation coefficients -- Spearman rho, Kendall tau -- should be suited to this purpose. (In particular, Kendall tau-b, as it is designed to handle ties.) But I am not a statistician, and I wonder if something more appropriate exists. After all, using a nonparametric statistical test ignores the fact that the one of the variables being compared is simply the list (1, 2, 3,..., N)! Or perhaps my understanding of what "nonparametric" means in this context is flawed. Any insight into what methods I can use to attack this problem would be appreciated.


title: Comparing two groups with many zeros

text: I am comparing the difference in time-activity-budgets of two populations of seabirds, those in the presence of ship disturbance and those not in the presence of ship disturbance. Focal animals were chosen haphazardly and observed for up to 10 minutes. Not all animals were observed for 10 minutes but the purpose of this analysis we have chosen to include all observations that exceed 3 mins. For this reason I have chosen to model proportion of time spent rather than actual time. I do have time (secs) spent in each category and total observation time if that makes the analysis easier to accomplish. This would need to be weighted somehow though. Three activities (flight, loafing, foraging) were defined and proportions sum to 1. I am interested in 1. determining if the time spent in flight between the two groups differs and 2. if there is a difference, do both the other activities decrease/increase or only one. Flight is not a common activity to observe in this species and therefor there are a lot of zeros in the data. Here is the distribution of proportion time spent flying for the two populations ![histogram of proportion time spent flying per observation](http://i.stack.imgur.com/G01jR.jpg) I think that I need to use generalized linear modeling with a zero inflated Poisson distribution but I am unclear on how to make that work correctly. Also it appears that the data is over-dispersed. I am working in R. Any help or advice for papers to read is appreciated. Thanks. **UPDATE: 2/6/14** I ran a zero-inflated negative binomial model and ship presence and time of day came out as significant predictors. Since my data was in proportional format I needed to transform the data to integers without loosing the weighting that proportions offered my. I multiplied all the proportions by 600 (max observation time) allowing each observation to have equal weight when using seconds flying. The mean % time flying was best predicted with presence of ship (0,1) and time of day (morning, mid-day, evening). The predicted times matched closely with the observed times. Interpretation: These birds appear to spend about 3 times more time flying when ships are present then absent for all three time periods. I then wanted to figure out if they were which activity(s) they were doing less of when ships were present (I measured three mutually exclusive events: flying, foraging, or loafing). Foraging time was also heavily dominated by zeros so I modeled it in the same way. This time when comparing the best fit model with predictors (just presence of ship this time) it did not improve on the fit of the null model (p = 0.10). Could this be that the true distribution is not zero-inflated negative binomial or can I interpret this as there is not a significant difference in the mean time spent foraging when ships are present and absent? This is how the raw data appears to me. The Mann-Whitney U test also comes up non-significant for foraging (although it is suggestive of a difference, p = 0.069). Assuming that there is no true difference in foraging and a marked difference in flying am I safe to assume that birds in this specific region (I am aware of scope of inference) compensate for extra flying time by loafing less rather than foraging less when ships are present (raw data supports this claim). Thanks for all the input so far.


title: Calculate probability distribution for n of participants in common to two parties given their size

text: I have a set of parties $V=\left\{ {V_i, V_j,..., V_n}\right\}$. Each party has a corresponding attribute $S=\left\{ {S_i, S_j,..., S_n}\right\}$ indicating the number of people participating to the event. For each pair of events $V_i$ and $V_j$ I can infer the number of participants in common as $E_{ij}$, with $E_{ij}=0$ meaning no one participate to both parties. Now I want to use my data to estimate, given $S_i$ and $S_j$, the probability distribution for $E_{ij}$ so to highlight which pair of $V_i$ and $V_j$ have a higher/lower $E_{ij}$ than expected given the data. In other words, I would like to see for each $E_{ij}$ a probability distribution like this ![enter image description here](http://i.stack.imgur.com/C8mxs.jpg) which (I guess) should tell me that given $S_i$ and $S_j$, $E_{ij}$ as probability .2 of being 3, .23 of being 4, and so on. Then If the _real_ value of $E_{ij}$ is, let's say, 10 I am in the position to tell that $E_{ij}$ has an unexpected value given the rest of the data. Does it make sense doing it? How should I calculate the probability distribution for $E_{ij}$? Should I use a linear regression with $S_i$ and $S_j$ as independent variables? NOTE: The figure was generated in R with the following code:               x = rbinom(100000, 20, 0.2)     plot(prop.table(table(x)),type="l")


title: Alternative to Shannon's entropy when probabilty equal to zero

text: I have a serie of objects for whom I know the probability of belonging to 10 classes. This probability can be null (see example below with 4 classes).                  A    B    C    D     1  0.4  0.0  0.2  0.4     2  0.1  0.3  0.4  0.2     3  0.0  0.0  0.0  1.0      In order to get for each object an information about the quality of the classification, I wanted to calculate the Shannon's entropy but it does not work when one of the classes has a probability equal to zero (log(0)=-Inf). My question: Is there a parameter similar with Shannon's entropy (or an adaptation) which can manage probability equal to zero ?


title: Non-significant p-values for factor levels with only 0s in negative binomial glm using glm.nb() in R

text: I am trying to fit a negative binomial GLM to fish catch data with month of the year (factor) as my explanatory variable. I have selected the month with the greatest number of catches as my reference level and fitted the model using the `glm.nb()` function from the MASS package.               model1<-glm.nb(catch~month+offset(log(Duration_hours*NumberHooksDeployed)),      catchdata)      Using `anova(model1)` shows that the factor `month` is significant and my pseudo $R^2$ is about 55%. The problem I have is that months with no catches at all have p-values greater than 0.05 (close to 1), even though I would expect them to be significantly lower than the reference level. Months with only very few catches are significant on the other hand. The output from the `predict()` function makes sense, so I think the model should be fine. Is there a problem with the likelihood test R uses, or am I fitting the wrong model to my data? From other posts I figured that I might have to use a likelihood ratio rather than Wald's test; any suggestions on how to do that in R? Any help would be greatly appreciated. See model output below, the months with 0 catches are May and June:               Call:     glm.nb(formula = NumCaught ~ factor(month) + offset(log(Effort_mins *      Nhooks)), data = f[f$Area == "PW", ], init.theta = 1.181267049,      link = log)          Deviance Residuals:      Min       1Q   Median       3Q      Max       -1.8445  -0.7766  -0.2148   0.0000   1.5386            Coefficients:                          Estimate Std. Error z value Pr(>|z|)         (Intercept)            -2.647e+00  3.192e-01  -8.292  < 2e-16 ***     factor(month)May       -2.754e+01  1.024e+05   0.000 0.999785         factor(month)June      -2.790e+01  3.122e+05   0.000 0.999929         factor(month)September -2.516e+00  6.516e-01  -3.861 0.000113 ***     factor(month)December  -1.112e-01  4.403e-01  -0.253 0.800543         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           (Dispersion parameter for Negative Binomial(1.1813) family taken to be 1)          Null deviance: 80.420  on 45  degrees of freedom     Residual deviance: 36.178  on 41  degrees of freedom     AIC: 140.81          Number of Fisher Scoring iterations: 1                         Theta:  1.181            Std. Err.:  0.525           2 x log-likelihood:  -128.815


title: Training a Logistic Regression Model

text: I created a quick fun Excel Spreadsheet tonight to try and predict which video games I'll enjoy if I buy them. I'm wondering if this quick example makes sense from a Logistic Regression perspective and if I am computing all of the values correctly. Unfortunately, if I did everything correctly I doubt I have much to look forward to on my XBOX or PS3 ;) I laid out a few categories and weighted them like so (Real spreadsheet lists twice as many or so): ` 4                   4             3         1 Visually Stunning   Exhilirating  Artistic  Sporty ` Then I went through some games I have and rated them in each category (ratings of 0-4). I then set a separate cell to be the value of Beta_0 and tuned that until the resulting percentages all looked about right. Next I entered in my expected ratings for the new games I **was** looking forward to and got percentages for those. Example: Beta_0 := -35 ` 4                   4             3         1 Visually Stunning   Exhilirating  Artistic  Sporty 4                  4             0         1 ` Would be calculated as P = 1 / [1 + e^(-35 + (4*4 + 4*4 + 3*0 + 1*1)] P = 88.1% If I were to automate the regression am I correct in thinking I'd be tuning Beta_0 to make it so the positive training examples come out high and the negative training examples come out low? I'm completely new to this (just started today thanks to this site actually!) so please have no concern about bruising my ego, I'm eager to learn more. Thanks!


title: Assessing reliability of a questionnaire:  dimensionality, problematic items, and whether to use alpha, lambda6 or some other index?

text: I am analyzing scores given by participants attending an experiment. I want to estimate the reliability of my questionnaire which is composed of 6 items aimed at estimating the attitude of the participants towards a product. I computed Cronbach's alpha treating all items as a single scale (alpha was about 0.6) and deleting one item at a time (max alpha was about 0.72). I know that alpha can be underestimated and overestimated depending on the number of items and the dimensionality of the underlying construct. So I also performed a PCA. This analysis revealed that there were three principal components explaining about 80% of the variance. So, my questions are all about how can I proceed now?   * Do I need to perform alpha computation on each of these dimension?    * Do I have remove the items affecting reliability? Further, searching on the web I found there is another measure of reliability: the lambda6 of guttman.   * What are the main differences between this measure and alpha?    * What is a good value of lambda?


title: About goodness-of-fit p-value and Pearson's $\chi^2$ statistic

text: In a simulation, I generate data from distribution A, and then I hope to test if each of some other distributions B,C,D... can fit the data equally well. Obviously there is an initial fit of the model under distribution A, so I can get the Pearson's chi^2 statistic as a measure of goodness-of-fit. My question: is it possible that under distributions B,C,D... I may get a even smaller chi^2 statistic than the one obtained from the true model (i.e. distn A)? The reason for this question is that I am considering whether a 1 or 2-sided test should be done, i.e. H0: model A fits the data; Ha: not H0. Intuitively (maybe I am wrong), a more general model (I mean, more general than A, such as A=Poisson, B=negative binomial) might provide an equally well or better fit for the data generated from A, so the Pearson's statistic is expected to be smaller if I use distn B to model data simulated from distn A. In this case, I should get a large p-value indicating no lack-of-fit, right? Any suggestions on choosing alternative statistics for evaluating GOF? Thank you very much!


title: Independence of a single sample?

text: Is it possible to speak of the independence of a single sample (of multiple data points) relative to the underlying distribution? Or is independence only relevant _between_ samples?


title: When should I use lasso vs ridge?

text: Say I want to estimate a large number of parameters, and I want to penalize some of them because I believe they should have little effect compared to the others. How do I decide what penalization scheme to use? When is ridge regression more appropriate? When should I use lasso?


title: Bayesian Updating Without Conjugate Prior

text: I'm trying to create some kind of iterative Bayesian algorithm, which continuously updates as more data is gathered. However, the distribution of my data is such that there does not exist a conjugate prior, so I'm thinking of using something like Gibbs sampling to generate my posterior distribution. However, if I want to do another iteration of this Bayesian algorithm with new data, is it possible without a conjugate prior? Because I only have the samples of the posterior distribution, not the actual distribution. Thanks!


title: Help on Statistical analysis of Profits

text: We have been working on a project where we have 2 or more trading agents on a market stock and for the outcome we have the profits that each trader generated after running 1000 times or more each market session. We need to statistically compare the profits of this traders to know which one is better, we have the idea of using confidence intervals on the mean, and we also where thinking about using the U-Test but we are not sure if its suitable to use it for such a large samples and comparing more than 2 traders. Any help will be appreciated. Thanks


title: How do I prepare data which has a trend for use in a Copula model?

text: I want to use a set of daily water quality data including 3 parameters in a Copula model. Somebody told me these data do not have a condition of a random variable to use in copula, and I should do some work on it (my guess is that it involves detrending the data). Can anybody help me in finding methods to prepare the data to use in a copula method? I would appreciate if you can refer me to good journal papers or books.


title: Likert scale and yes/no calculation of one score

text: My survey is based on observations and interviews. The observation questions are yes/no type questions and interview questions are 4 points. > $1=$ never, $2=$ sometimes, $3=$ usually, and $4=$ always.   >  Yes $=1$, No $=0$. If I add up all these with highest possible score and divide it by the total score received, I can then retrieve a total score for the survey. For example, let's say I have 5 responses for the Likert and 5 for the yes/no. $5\times4$ $=20$ (with 4 being the highest possible score you can get on the Likert scale) and $5\times1$ $=5$ $(1$ being the highest possible score for yes/no question). I then add up both to find the highest possible score $=25$. If I get a combined score of $20$, do I divide it by the possible score, $20\div25=80$%? Can I build my questionnaire and calculations like this? My main objective is to get one percentage / score for the entire questionnaire. I thought of accomplishing this in 2 ways:   1. Get the sum of earned score from both scales divided by total possible scores for both scales.    2. Get a percentage by using the above calculation for individual scales, which will give me two sets of percentages, and then take an average for those percentages.  I am not sure which way I should go with, or if this even makes sense doing it this way.


title: Effect sizes from LMM

text: My question is, how to get effect sizes for a linear mixed model? I am using the following model in SPSS:               MIXED       transfer  BY distance training rotation sequence  WITH pretest       /CRITERIA = CIN(95) MXITER(100) MXSTEP(5) SCORING(1) SINGULAR(0.000000000001)           HCONVERGE(0, ABSOLUTE) LCONVERGE(0, ABSOLUTE) PCONVERGE(0.000001, ABSOLUTE)       /FIXED = pretest rotation sequence distance training distance*training  | SSTYPE(3)       /METHOD = REML       /REPEATED = distance | SUBJECT(ResponseID) COVTYPE(CS) .      I've also done it in R, where it looks like this:               library( nlme )     options(contrasts=c("contr.sum","contr.poly"))     lm1 <- lme(transfer ~ training * distance + rotation + sequence + pretest,                 random=~1|ResponseID, method="REML", data=wide.data )     fit <- anova.lme(lm1,type='marginal')     print( fit )      The dv, transfer, measures pretest to posttest improvement for a given level of distance. Distance is a categorical within-subjects variable with either 2 or 3 levels depending on the dataset. Training is a categorical between- subjects variable with either 2 or 4 levels depending on the dataset. Rotation and sequence are binary categorical between-subjects variables. Pretest is a covariate with different values for each level of distance. I've already submitted the results of the analysis for publication, and got a reviewer comment asking for odds ratios. Ideally I'd like to get odds ratios for each effect in the model. If I can't get odds ratios, some other measure of effect sizes might be OK. Most important, I need to get it for the training * distance interaction and the covariate (pretest) because these are the only significant effects. Ideally I'd get it straight from SPSS but as far as I can see, SPSS can't do it. I couldn't figure out how to do it in R either. Second best would be to put the output of one of these into some other software that could calculate it for me. I found some free software that could compute eta-squared for regular mixed ANOVA but not for linear mixed model. Other questions about effect sizes seem to relate to other models, not LMMs. I did see another similar question about LMMs which is currently unanswered. Any ideas?


title: Difference between p values of interaction

text: I was wondering which p-value to use when i am testing the significance of an interaction between two continuous variables. Here is my model with the interaction: Coefficients:                                 Estimate Std. Error t value Pr(>|t|)         (Intercept)       110.2584    23.9969   4.595 0.000417 ***     lg_hag              36.0716    18.2778   1.974 0.068510 .           raceblack         -12.9419     1.8198  -7.112 5.24e-06 ***     racemexican        -9.0254     2.3207  -3.889 0.001636 **      racemulti/other    -9.8012     2.0640  -4.749 0.000311 ***     sexfemale           3.9612     1.3405   2.955 0.010442 *       age_yr             -1.1009     0.1691  -6.511 1.38e-05 ***     as.factor(educa)2   6.7660     2.2988   2.943 0.010686 *       as.factor(educa)3  14.1302     2.2394   6.310 1.92e-05 ***     as.factor(educa)4  15.4776     1.8657   8.296 8.95e-07 ***     as.factor(educa)5  19.8650     2.0524   9.679 1.40e-07 ***     lg_hag:age_yr       -0.5236     0.2555  -2.049 0.07862 .       So in this model, the p value for the interaction is 0.07862, limit to the significance However, when I run the model without interaction, and then i do a anova (with a Chi square test as the models are GLM) between the two models (with and without interaction), i got a p-value of 0.04 (which is supposed to be the real p value of inteaction, in my view)               Working (Rao-Scott) LRT for lg_hag:age_yr     Working 2logLR =  4.199194 p= 0.041893      Please if someone could give me an advice on what is the difference between the two p value, and which one to use to attest for the significance of the interaction. Thank you in advance


title: Using the normal equations to calculate coefficients in multiple linear regression

text: I am trying to understand how to get the coefficient of a multiple linear regression. The formula is: $b = (X'X)^{-1}(X')Y$ I try to calculate $b$ without package and with the `lm` package inside R. Doing so, I got different results. I want to know why. Did I made a mistake? Or does the `lm` package calculate differently because of the intercept?               > y <-  c(1,2,3,4,5)     > x1 <- c(1,2,3,4,5)     > x2 <- c(1,4,5,7,9)     > Y <- as.matrix(y)     > X <- as.matrix(cbind(x1,x2))     > beta = solve(t(X) %*% X) %*% (t(X) %*% Y) ; beta                 [,1]     x1  1.000000e+00     x2 -1.421085e-14     > model <- lm(y~x1+x2) ; model$coefficients      (Intercept)           x1           x2      1.191616e-15 1.000000e+00 1.192934e-15       ![3d](http://i.stack.imgur.com/KHD2q.png) # Update As Alex and the other told me, it was a question of roundoff error. Therefore, I decided to take another data from the book "Essential Statistics for business and economics" by Anderson and all. In this case, the coefficients are the same in both `lm` function and in my own matrix.               > y <- c(9.3, 4.8, 8.9, 6.5, 4.2, 6.2, 7.4, 6, 7.6, 6.1)     > x0 <- c(1,1,1,1,1,1,1,1,1,1)      > x1 <-  c(100,50,100,100,50,80,75,65,90,90)     > x2 <- c(4,3,4,2,2,2,3,4,3,2)     > Y <- as.matrix(y)     > X <- as.matrix(cbind(x0,x1,x2))          > beta = solve(t(X) %*% X) %*% (t(X) %*% Y);beta              [,1]     x0 -0.8687015     x1  0.0611346     x2  0.9234254     > model <- lm(y~+x1+x2) ; model$coefficients     (Intercept)          x1          x2       -0.8687015   0.0611346   0.9234254


title: Beta Distribution with Random Shape Parameters (Binomial or Poisson Binomial)

text: I've run into an interesting model into some of my research and I was hoping somebody with more mileage can point me to any resources or other applications. I'm familiar with the Beta-Binomial, but now I'm trying to do some analysis where the shape parameters associated with the Beta prior are RVs distributed as Binomial or Poisson-Binomial. Has anybody seen this or studied it before?


title: Mantel test vs. Pearson's correlation coefficient

text: Why is it correct not to measure the correlation between two distance matrices by calculating the Pearson's correlation coefficient between all distance pairs but to use the Mantel test instead? Edit: Wikipedia says that because "distances are not independent of each other – since changing the "position" of one object would change n-1 of these distances (the distance from that object to each of the others) – we can not assess the relationship between the two matrices by simply evaluating the correlation coefficient between the two sets of distances and testing its statistical significance. The Mantel test deals with this problem." This argument is not clear to me. Why is dependence between the distances a problem?


title: Generalization of multivariate normal distribution and classification

text: I am interested in a family of multivariate distributions that can be seen as a generalization of the multivariate normal distribution, insofar as they are defined by an expectation value $\vec \mu$ and a covariance matrix $\Sigma$, plus a monotonously decreasing function $g(d)$ such that the density is $$ p(\vec x) \propto g \left ( \Delta(\vec x, \vec \mu) \right ) $$ where $$ \Delta(\vec a, \vec b) = \sqrt { (\vec a - \vec b)^T \Sigma^{-1} (\vec a - \vec b) } $$ is the Mahalanobis distance. The multivariate normal is of course recovered by $g(d) = \exp (- \frac12 d^2 )$. My **first question** is: What is the name of this family of distributions? It is simple to show that for classification of a given data point to one of two or more classes, each of which is described by such a density with different $\mu$ but identical $\Sigma$ and $g(d)$, optimal classification boundaries are piecewise linear (hyperplanar). My **second question** is: Is this a standard result, and if yes, what is the standard literature (textbook) reference for it?


title: Farlie-Gumbel-Morgenstern Bivariate Gamma Distirbution

text: Given the variables $X$ and $Y$, which are correlated, $X\ge0$, $Y\ge0$ and each follow a gamma distribution with different shape parameters, i.e.,$X\sim Gamma(a_1,\alpha)$ and $Y\sim Gamma(a_2,\alpha)$. I understood that the Joint PDF $f_{X,Y}(x,y)$ can be obtained by making use of the Farlie Gumbel Morgenstern (FGM) Bivariate Gamma Copula which is applicable for the case of different shape parameters. The joint PDF is given by $$h_{X,Y}(x,y) = f_X(x)f_Y(y) \left[1+\lambda\left(2F_X(x)-1\right)\left(2F_Y(y)-1\right)\right],\quad |\lambda|\le 1 $$ where $f(.)$ denotes the marginal PDF and $F(.)$ denotes the marginal CDF. My major questions are:   * How to determine the value of $\lambda$?   * Is it true that the FGM copula is used only for weak correlation?    * Is there another method of calculating Pearson's correlation coefficient rather than the equation given by $$ \rho = \lambda K(a_1)K(a_2)$$ where $K(a)=\displaystyle\frac{1}{2^{2a-1} \beta(a,a) \sqrt{a}}$ and $\beta(.,.)$ is the beta function.


title: Coefficients of bootstrapping in logistic regression

text: I have seen several articles and CrossValidated questions on bootstrapping ( this, this or this for example); there are a lot of theoretical and statistical explanations, however since they are so theory based, I am afraid I might be understanding the use wrongly. Hence my questions: 1) When I make a non-parametric bootstrapping (changing the sample for every run) with logistic regression on my data, I basically will end up with several different coefficients for each predictor for each run. Eventually I'll have the confidence interval for each predictor as well. I understand until that point. My question is; assuming that the distribution is normal, when I want to come up with a final model on practice, can I just take the mean of the confidence intervals for each predictor and consider this as my final model coefficient? 2) If the answer to question #1 is yes, is this the only way of choosing coefficients while bootstrapping? If not, what else? I encountered in a few more articles a method called "bagging". This seems to be my main purpose. 3) This one is more of a curiosity question: Can above methodology be applied to the categorical predictors when they are assigned with Weight Of Evidences? I know we can split the categorical predictors into dummy variables; but how would I treat each coefficient if I want to use WOE methodology?


title: Neural networks creates negative output

text: I am using a simple feedforward neural network in MATLAB to predict the output for inputs in the range [1e-5, 0.3]. (These are the activations of another network.) I am using a sigmoid function for the hidden layer, and a linear function for the output layer. Input units are 6, hidden units are 4, and the output unit consists of one neuron. The outputs range between [58, 1696]. I normalized the outputs too and turned off the mapminmax function in order to avoid over-normalizing! Weirdly, it creates negative outputs. Could it be because of the input range? I would appreciate if anyone could tell me what is happening here. Any thoughts? First I need to make an update : I changed the number of hidden neurons to 10. Now one time it gives me negative outputs the other time positive. I cannot find the answer except relating it to the random initialization of the weights. Here is the excerpt of my code :               clear all;     load features;load labeldata;           net=feedforwardnet(10);       IPF={'fixunknowns','remconstantrows'} ; OPF={'remconstantrows'};     net.trainParam.lr=0.01;     net.trainParam.max_fail=10;     NN = train(net,features,y); %y is normalized value of labeldata so it is in range[0,1]     wb = getwb(NN);     net=NN;          %% TEST THE NET (NOT WITH NEW TEST DATA BUT WITH THE TRAINING DATA, SO WE EXPECT GOOD RESULTS)     [pred_learnedFeatures]=net(features);     % scaled_out=pred_learnedFeatures*(max(labeldata)-min(labeldata))+min(labeldata);     fprintf('MSE w/o scaling %f%\n',sum(( pred_learnedFeatures- y).^2)/size(labeldata,2));     fprintf('\n MSE w scaling %f%\n',sum(( scaled_out-labeldata).^2)/size(labeldata,2));


title: How do you calculate simple statistics for left censored data in R?

text: I have left censored data where the distribution is known (it's near enough lognormal, at least in theory). I'd like to calculate some simple summary stats: geometric mean and standard deviation in this case. I've previously used R's `NADA` package for this but it is no longer on CRAN. Is there an alternative available? EDIT: I contacted Lopaka Lee, the package maintainer, and he says that he's > working on updating the package for the upcoming release of R so hopefully the package's absence is only temporary.


title: Discerning between two different linear regression models in one sample

text: Suppose I observe a sample $(y_i,x_i)$, $i=1,...,n$. Suppose that I know the following: $y_i=\alpha_0+\alpha_1x_i+\varepsilon_i$, $i \in J\subset\{1,...,n\}$ $y_i=\beta_0+\beta_1x_i+\varepsilon_i$, $i \in J^c$ where $\varepsilon_i$ are i. i. d. and $J$ is not known in advance. Is it possible to estimate $\alpha_0,\alpha_1,\beta_0,\beta_1$? Or at least test the hypothesis that $J=\varnothing$? If $J$ is known the problem is very easy to solve. Going through all the subsets is not feasible, since we have $2^n$ possible combinations. If we assume $J=\{1,...,k\}$ with unknown $k=1,...,n$, it is the classical change- point problem, for which many tests are available. I suspect that this maybe ill-posed problem, so I wanted to check before trying to solve it. Here is a simple illustration of the problem:               N <- 200     s1 <- sample(1:N,N %/% 2)     s2 <- (1:N)[!(1:N) %in% s1]          x <- rnorm(N)     eps <- rnorm(N)          ind <- 1:N          y <- rep(NA,N*T)     y[ind %in% s1] <- 2+0.5*x[ind %in% s1]+eps[ind %in% s1]/5     y[ind %in% s2] <- 1+1*x[ind %in% s2]+eps[ind %in% s2]/5     y          sal1 <- ind %in% s1          plot(x, y)     points(x[sal1], y[sal1], col=2)     abline(2, 0.5, col=2)     abline(1, 1)      Graphically it is more or less obvious that we have two different models. Maybe it is possible to use some classification or data-mining techniques for solving this problem?


title: Basic Question: how does feed-forward neural network solve regression?

text: This is a fairly basic question but I can't seem to find an answer on the net (perhaps I'm searching the wrong things). Regression is trying to predict continuous outputs. Since a neural network uses a clamping function (typically giving a value between 0 and 1) before the output.. if the output can only be between 0 and 1, how can a neural network learn a regression function? Thanks in advance


title: R, arima() with include.mean=TRUE, still has no mean reported

text: I have a regression with ARMA errors, which I am fitting with `arima()`. I know that the ARMA model is being fit on these residuals from the regression. My problem is that when I use `include.mean=TRUE` the output does not return any estimate of the mean (of the residuals, which yes, I know is always zero). The documentation says that ARMA model will include a mean if this option is set to `TRUE`. I do see an intercept reported, which is the intercept from the regression of y on x, not the intercept from the ARMA model on the residuals. How do I get the `include.mean` argument to work? Is there any other way to get the intercept from the ARMA model on residuals? Is R just censoring this because residuals are always mean zero (if intercept is in regression)? MWE is below, thanks in advance. `set.seed(123) y = as.xts(ts(rnorm(20)+3, start=c(1980,1), freq=12)) x = as.xts(ts(rnorm(20)+5, start=c(1980,1), freq=12))` `fit1 = arima(y, xreg=as.data.frame(x), order=c(1,0,0), include.mean=TRUE) fit1 # intercept shown is from the regression, 3-(-.0935)*(5), not for the residual ARMA model`


title: Textbook deriving Metropolis-Hastings and Gibbs Sampling

text: I have fairly good practical experience with Metropolis-Hastings and Gibbs sampling, but I want to get a better mathematical understanding of these algorithms. What are some good textbooks or articles that prove the correctness of these samplers (more algorithms would also be great)?


title: Why is the true (test) error rate of any classifier 50%?

text: In section 7.10.2 of Elements of Statistical Learning, it says that the true (test) error rate of any classifier is 50%. I'm having trouble understanding the intuition behind this. If you have a binary class (1 or 0) and your classifier is a die where if you roll a 1-5, the classification is 1 and if you roll a 6, the classification is 0. Then suppose that the true value of your binary class is 1. Then I would think that the error rate would be converge to 1/6 over time. The excerpt from the text is below. > Consider a classification problem with a large number of predictors, as may > arise, for example, in genomic or proteomic applications. A typical strategy > for analysis might be as follows: 1. Screen the predictors: find a subset of > “good” predictors that show fairly strong (univariate) correlation with the > class labels 2. Using just this subset of predictors, build a multivariate > classifier. 3. Use cross-validation to estimate the unknown tuning > parameters and to estimate the prediction error of the final model. Is this > a correct application of cross-validation? Consider a scenario with N = 50 > samples in two equal-sized classes, and p = 5000 quantitative predictors > (standard Gaussian) that are independent of the class labels. The true > (test) error rate of any classifier is 50%. We carried out the above recipe, > choosing in step (1) the 100 predictors having highest correlation with the > class labels, and then using a 1-nearest neighbor classifier, based on just > these 100 predictors, in step (2). Over 50 simulations from this setting, > the average CV error rate was 3%. This is far lower than the true error rate > of 50%.


title: Question regarding parameters and variable selection in Mahout algorithm for logistic regression

text: Below is the list of parameters in Mahout logistic regression.   1. What does "passes" mean? In detail please  > --passes passes the number of times to pass over the input data   2. What does "features" mean? Is it the number of variables I am taking as predictor? But I see people entering more or less at times for the number of variables.   3. Is there any way of selecting variables with good predictive power in machine learning something like ginie coeffecient?    4. How do we deal with multicollinearity in machine learning? Is it by regularisation? In statistics we remove the variable using variance inflation factor. Is it computed in machine learning?              Command: mahout trainlogistic     Flag       Any argument to flag     Description     --help                              print this list     --quiet                             be extra quiet     --input       input                 where to get training data     --output      output                where to get training data     --target      target                the name of the target variable     --categories  number                the number of target categories to be considered     --predictors p [p ...]              a list of predictor variables     --types      t [t ...]       a list of predictor variables types (numeric, word or text)     --passes     passes                 the number of times to pass over the input data     --lambda     lambda                 the amount of coeffiecient decay to use     --rate       learningRate           the learning rate     --noBias                            do not include a bias term     --features   numFeatures            the number of internal hashed features to use


title: Friedman chi-squared = 0. What does it mean?

text: My data was paired and do not fit normality. So I used a Friedman test.However, the result turn out to be Friedman chi-squared = 0, df = 1, p-value = 1. How do this happen? My data was in this way:                   I    II     1   2    1     2   4    1     3   5    0     4   0    2     5  14    3     6  11    3     7   7    6     8   8    9     9  12   19     10 15   13     11 19   32     12 22   74     13 15    1     14  2    1     15  2    4     16  1    6     17  5    9     18  1    8     19  6    6     20  8   11     21 12   17     22 26   22     23 38   29     24 43   58     25 71   20


title: I've been trying to wrap my head around the use of eigenvalues in cluster analysis. What does it tell me about my clustering behavior?

text: In a typical hierarchical cluster output from using SAS, the first table given lists all of the eigenvalues. From what I understand, eigenvalues are derived from covariance between the variables. What I fail to understand is how this would assist me in understanding the underlying clustering behavior. Anyone know a simple answer of why the eigenvalues are important to know for cluster analysis? Here's the output: ![enter image description here](http://i.stack.imgur.com/uIIhe.png)


title: Transforming confusion matrix to distance matrix

text: I have an annotation task which I need to analyze for accuracy, category fine tuning and reliability , there are many annotators assigning multiple set valued categories to items. And I have come up with confusion matrices for each pair of annotators like below:               +-----------------------------------------+     |  |5|  .   1   .   .   .   .   2   .   . |     |   7  |.|  .   .   .   .   .   5   .   1 |     |  34   .  |.|  .   .   .   1   4   .   1 |     |   7   1   .  |.|  .   .   1  46   .   . |     |   .   .   .   .  |.|  1   .   .   .   . |     |   5   .   .   .   .  |.|  .   1   .   . |     |   5   .   .   .   .   .  |.|  .   .   . |     |   5   .   .   .   .   .   1  |2|  .   5 |     | 198   .   1   .   .   2   1  11  |.|  1 |     |   1   .   .   .   .   .   2   .   .  |.||     ------------+-----------------------------+       To transform it into an undirected distance matrix , I transposed the above and averaged it to get the following symmetric matrix which I need to represent as distance matrix:               |   |5|   3   17    3    .    2    2    3   99    0 |     |    3   |.|   .    0    .    .    .    2    .    0 |     |   17    .   |.|   .    .    .    0    2    0    0 |     |    3    0    .   |.|   .    .    0   23    .    . |     |    .    .    .    .   |.|   0    .    .    .    . |     |    2    .    .    .    0   |.|   .    0    1    . |     |    2    .    0    0    .    .   |.|   0    0    1 |     |    3    2    2   23    .    0    0   |2|   5    2 |     |   99    .    0    .    .    1    0    5   |.|   0 |     |    0    0    0    .    .    .    1    2    0   |.||      I have about 190 such matrices from multiple annotator pairs and there is nothing like reference and test categories as the annotator category assignment needs to be analyzed for fine tuning, merging the categories. So is the above approach correct to start with and how to transform multiple such matrices (dimensions)into a single visualization allowing to analyze the categories and inter annotator reliability apart from PCA and the usual kappa and related statistics? Also would such a distance matrix give a correct indication of the real components in PCA for cluster analysis to base further decisions on?


title: Inferring regression coefficient of two variables through a third

text: I have 3 time series X, Y, Z If the regression coefficient (slope) of X and Y is **a** , Y and Z is **b** , then what can I say about the regression coefficient of X and Z?


title: Log prices vs log returns

text: I am doing a basic analysis for my research, regarding the relationship between the global oil & gas price and the Australian aggregate and sectoral equity markets. My supervisor has me using EVIEWS 8 courtesy of its availability at university. My data comprises of price indexes for World O&G price, Aus financial market, Aus O&G market, etc. My supervisor instructed me to use returns in order to prevent serial correlation in error term, but I'm not entirely sure I have it correct. Could someone clarify please? I don't really understand the difference between taking ln(Price) and ln(First Difference), and whether ln(First difference) can be taken as the returns of that particular market.


title: Getting spline coefficients in R

text: I'm fitting a natural basis spline on a data set of the form:               splineModel=lm(dist~bs(speed, df=3), data=cars)      using `bs` function of `spline` package. I would like to know whether is possible to print the actual spline equations fitted (the tracing) within each knot. I would like to have the coefficients of all piecewise fitted polynomials.


title: (very basic) One-sample test for binary data

text: I've repeatedly measured a continuous variable and each measure has been assigned a populational percentile range it falls into (percentile ranges were estimated for general population in another study). The exemplary barplot depicts the share of each group in my sample during different time points. I'd like to test, whether the percentage, that specific range constitutes in certain time point (eg. share of 90-100th percentile in the 6th time point) differs between my group and general population. The simplest way would be to perform a one-sample t-test (eg. against the mean of 0.1 in case of 90-100th percentile, as 10% of values will fall above 90th percentile in the population). But is there any alternative if my sample gets small? It would be ideal to perform Fisher's exact test, but I have no reference group - I can only assume, that in general population 1/10 samples will fall into that percentile range. ![Imgur](http://i.imgur.com/Kvc9R48.png)


title: How to get the difference of two numbers (standard deviation)?

text: I have many pairs like:               A  B     10 20     15.3 19.5     23 13     45 32      etc etc, I create those pairs calculating the standard deviation of many vectors (R vectors) and then group them by 2(A - B). Now, I have to filter those pairs by their percentage difference. I set a percentage like `30%` then I filter all the pairs to get those who have a difference **BELOW** or **EQUAL** 30% My question is, is this formula correct to get the difference (in % of two numbers)               ((A / B) -1) * 100      My doubt is that I need their difference to understand if i can take them or not but if I do: ((B / A) -1) *100 obviously I get a different result, example:               A = 10     B = 20     A-B % = -50 %  ((10 / 20) -1) *100     B-A % = 100 %  ((20 / 10) -1) *100      So, How can I better study the real difference of two numbers(standard deviation difference)? **EDIT (@Henry)**               > a <- log(sd(rnorm(250)))     > b <- log(sd(rnorm(250)))     > a/b     [1] 0.6049963      **0.60%** of change ? Thank you!


title: What happens if the full model has the smallest AIC?

text: Suppose I am modelling a binary outcome $y$ with $n$ covariates. If this model has the lowest AIC among all the other ones...would this be the one to use? Also in dealing with categorical variables, should you delete all levels or none? For example suppose a categorical variable $z$ has the following levels: $z_1, z_2$ and $z_3$. If $z_1$ is the reference level then I include $z_2$ and $z_3$ in the model. If I don't want to include the variable $z$ in the model....would I delete both $z_2$ and $z_3$? Or could I include one of them?


title: Histogram/distribution fitting for this dataset with unequal and open-ended intervals?

text: I have this income distribution data for various groups: https://docs.google.com/spreadsheet/ccc?key=0Akwg3n_e05cCdEdtT0VZYU5keW5DVkNoNmpBWmdzeUE As you can see, I have intervals/bins with varying widths. I also have an open-ended interval and the mean income in that open interval varies a lot between groups. **For pedagogical purposes, I'd like to be able to calculate how many people and/or percent is in a given interval, given an equal bin width.** Example:   On an x-axis from 0-2000 with a bin width of 10, I'd like to be able to say how many people is in the 590-600 group. 1) Is this doable?   2) Do I need to fit a distribution to a histogram or how do I do it?   2b) What would you base this histogram / distribution on (mean, median etc.)? (As mentioned, the mean in the open-ended interval varies a lot between groups, something I'd like to take into consideration = show truthfully graphically.)   3) Can you please explain this in simple mathematical operations based on the data given?


title: How to compare two power law distributions?

text: Disclaimer: I am a biologist, so please don't hesitate to correct me if I am making unwise assumptions. I am looking at the distribution of a unique 6mer in a series of 2100mer (I have 415 independent 2100mer). More precisely, I am looking into how many times that 6mer is present in each 2100mer. As I deal with DNA sequence (A, T, G or C), each 6mer is expected to be found in random DNA every 4096 bp or 0.51 times in each 2100mer. Clearly a non-normal distribution (qqplot). Thus, I believe I am looking at a power law distribution or what should be a power law. Now, I want to compare that observed distribution to a simulated one and see if they are the same. Here are the raw counts:               RawObserved = c(87, 112, 111, 60, 20, 18, 7)       RawExpected = c(44.6, 28.7, 19.0, 7.7, 2.0, 1.5, 0.5)        Here is my data rendered as proportion:               NormObserved = c(0.21, 0.27, 0.27, 0.14, 0.05, 0.04, 0.02)      NormExpected = c(0.11, 0.07, 0.05, 0.02, 0.00, 0.00, 0.00)        I clearly cannot use a t-test but I was surprised to see that a chi-square test returned a non-significant result. To me those distributions look very different, qualitatively and I would like to statistically test that. I hope my explanation is sufficient and makes sense. Thanks!


title: What distribution to use for this QQ plot?

text: I have a dataset and I made the QQ-plot against the $N(0,1)$ distribution. The plot is included below. ![](http://i.stack.imgur.com/MlXsY.png) My statistics is rusty to say the least (meaning what little knowledge I had is now rusted away!) Clearly the normal distribution is not a good fit -- the tails in my data are heavier than those in a normal distribution. But the data appears not to be skewed relative to the normal distribution. So, which distribution might be a better fit? Is it something a transformation like Box-Cox could cure to make it fit better? Edit: My data isn't strictly positive so Box-Cox is out... But there may be another transformation that works. Edit 2: I have increasingly large datasets and I need to see what, if anything, the datasets converge to. This is exploratory analysis, I would rather find a distribution that fits the data rather than transforming the data to a distribution. All of this is done with SciPy which reports the asymptotic value of the biased kurtosis as -1 and skewness as 0. But I don't know how to use that information to determine which distribution this might be, aside from checking all of them to see which ones give a better $r^2$. Edit 3: Based on a comment from gung, I checked it against a uniform distribution: ![](http://i.stack.imgur.com/vNETx.png) Sure enough, that's considerably better although still shows a difference in the tails.


title: Bootstrapping vs. permutation tests

text: There are several popular resampling techniques, which are often used in practice, such as bootstrapping, permutation test, jackknife, etc. There are numerous articles & books discuss these techniques, for example Philip I Good (2010) Permutation, Parametric, and Bootstrap Tests of Hypotheses My question is which resampling technique has gained the more popularity and easier to implement? Bootstrapping or permutation tests? Many Thanks


title: Is it true that $F_U(U(\bar\omega)) = F_X(z)$ implies $z=X(\bar\omega)$?

text: For a statistics class, I have to prove a result which leads me to the following question. If I can show that it is true, my proof is done. So here is the question. Suppose that $U: \Omega \rightarrow [0,1]$ is a random variable uniformly distributed on $[0, 1]$ and $X: \Omega \rightarrow \mathbb{R}$ is continuously distributed as some distribution $F_X(·)$ (where $\Omega$ is the sample space). Is it true that > $F_U(U(\bar{\omega})) = F_X(z_{\bar{\omega}})$ for all $\bar{\omega} \in > \Omega \implies z_{\bar{\omega}}=X(\bar{\omega})$ or equivalently, > $P \{\omega \in \Omega~|~ U(\omega) \leq U(\bar{\omega})\} = P \{\omega \in > \Omega~|~ X(\omega) \leq z_{\bar{\omega}}\}$ for all $\bar{\omega} \in > \Omega \implies z_{\bar{\omega}}=X(\bar{\omega})$ ? I suspect that it is true, but I haven't been able to prove or to disprove it.


title: Neural network  model to predict treatment outcome

text: I want to predict the outcome of a particular treatment (remitted or not) using demographic, plasma biomarker, genetic, and clinical data. IS a neural network model the best way of doing this? What advantages does this have over traditional logistic regression model building? How limited am I with only 120 cases and up to 40 covariates, depending on collinearity? How do I pare these down? I would normally tend towards factor analysis but will a neural net combine collinear variables sensibly? Any ideas on combining multimodal data like that would be helpful, or a starting point for reading - already have Ripley's MASS.


title: Relative importance weight with cforest

text: I am new in using RF. I want to use it to compute the relative importance of the features. I found the weight is very small ("party" package, cforest). Is there anyway to get these weights in a range of 0-1? The total weight would sum up to 1? For example, if $x_1$, $x_2$ and $x_3$ are the features, the relative importance for these features are something like, $r_1: 0.5, r_2:0.4$ and $r_3:0.1$. Thanks in advance.


title: How to shape a distribution?

text: I have a data set of teams. Each team has a sum of points. I want to divide the teams in ranges of points (0 - 999; 1000 - 1999; etc) such as that the top range (with the most points) is smaller (has less teams) than the next and so on. Suppose the number of ranges is 7 and the angle is X degrees like in this figure: ![a](http://i1198.photobucket.com/albums/aa451/KakaOriginal/album/distribuicao.png) How do I find the point ranges' limits to form that figure? EDIT: The bars have the same base size regardless of the respective range size.


title: Calculating forecast error with time series cross-validation

text: I have a forecasting model for a time series and I want to calculate its out- of-sample prediction error. At the moment the strategy I'm following is the one suggested on Rob Hyndman's blog (near the bottom of the page) which goes like this (assuming a time series $y_1,\dots,y_n$ and a training set of size $k$)   1. Fit the model to the data $y_t,\dots,y_{t+k-1}$ and let $\hat{y}_{t+k}$ be the forecast for the next observation.   2. Compute the forecast error as $e_{t} = \hat{y}_{t+k} - y_{t+k}$.   3. Repeat for $t=1,\dots,n-k$   4. Compute the mean square error as $\textrm{MSE}=\frac{1}{n-k}\sum_{t=1}^{n-k} e_t^2$ My question is how much I have to worry about correlations because of my overlapping training sets. In particular, say I want to forecast not only the next value, but the next $m$ values, so that I have predictions $\hat{y}_{t+k},\dots,\hat{y}_{t+k+m-1}$ and errors $e_{t,1},\dots,e_{t,m}$, and I want to construct a term-structure of prediction errors. Can I still roll the window of the training set forward by 1 each time, or should I roll it forward by $m$? How do the answers to these questions change if there is significant autocorrelation in the series that I'm predicting (conceivably it is a long-memory process, i.e. the autocorrelation function decays as a power law rather than exponentially.) I'd appreciate either an explanation here, or links to somewhere where I can find theoretical results about the confidence intervals around the MSE (or other error measures).


title: Maximum Likelihood Estimation question - minimum log likelihood

text: I know the formula for the likelihood of some parameters given the data. The result has to be maximised and I can avoid multiplication using the log. How can I make this a minimisation problem (i.e. rewrite the given formula). I presume it is not as easy as taking the negative log and than adding things - or is it? Thanks. Christian


title: Is there a way to measure the information content (or difference) in two matrices?

text: **Prologue** I have a matrix $N$ with dimensions $I\times J$, $I >> J$, and $n_{ij}$ is assumed to follow a Poisson distribution with mean $\lambda_{ij}$. Using a method (say _NtoM_ ), I reduce $N$ to a matrix $M$ with dimensions $I \times K$ (rows correspond to that of $N$, however columns are not representing the same quantities as in $N$), $I >> K$, and the entries of $M$ are positive real numbers. _Edit_ (thanks to @whuber's comment) This is the relation between $N$ and $M$: For any $i \in \{1,\ldots, I\}$ and a particular $j$, $\lambda_{ij}$ will be on of the $m_{ik}$, $k \in \{1,\ldots, K\}$. However, the key is that for a particular $i$, $\lambda_{ij_1}$ and $\lambda_{ij_2}$ need not be the same if $j_1\neq j_2$. $K$ can be assumed to be less than $J$. **Question** Is there a way to measure the difference of information content (gain or loss) in reducing matrix $N$ to $M$ using the method _NtoM_?


title: What is a minimum sample size for a paired t-test and what is a non-parametric equivalent if data is non-normal?

text: I have 4 sample pairs, let's say (x1, y1),...,(x4, y4).   * What is the minimum sample size for a paired t-test?   * Which assumption should I check for paired t-test?    * If my data is non-normal, what is an alternative non-parametric test?


title: Highly irregular time series

text: I have data for the population of a number of different fish, sampled over a period of about 5 years, but in a very irregular pattern. Sometimes there are months between samples, sometimes there are several samples in one month. There are also many 0 counts How to deal with such data? I can graph it easily enough in R, but the graphs are not particularly illuminating, because they are very bumpy. In terms of modeling - with species modeled as a function of various things - maybe a mixed model (aka multilevel model). Any references or ideas welcome Some details in response to comments There are about 15 species. I am trying to both get an idea of any trends or seasonality in each fish, and look at how the species are related to each other (my client originally wanted a simple table of correlations) The goal is descriptive and analytic, not predictive Further edits: I did find this paper by K. Rehfield et al., which suggests using Gaussian kernels to estimate the ACF for highly irregular time series http://www.nonlin-processes-geophys.net/18/389/2011/npg-18-389-2011.pdf


title: Sample size via confidence intervals

text: I'm looking for a R package to compute sample size for the differences of two proportions via confidence inetrval (CI). Note that there are several methods to compute CI for the difference in proportions, and for some of these methods one can compute sample size with pencil and paper, but I was wondering if there is package invoving several of these methods. Thanks


title: Implementation for Co-Clustering

text: I am looking for existing implementations for co-clustering (aka biclustering). I came up with `biclust` function available in MATLAB, but still I am wondering if there are more implementations for this subject, either in MATLAB or R.


title: Finding the equation for a data set

text: I asked this question on StackOverflow and someone pointed for me to ask this here. Lets assume that I have 3 variables: x, y, z which are integers. I have 5000 values of each and I want to predict the 5001th (or other) value given this data set. The thing is, the values oscillate a lot so I was thinking that I could possibly find the equation of this data set computationally, and I could approximate the 5001th value using this equation. What is the preferred way for me to do this?


title: Dropping term for correlation between random effects in lme and interpretting summary output

text: I want to fit a model without a correlation term between the random effects with `lme`. In `lmer` this is fairly straighforward....               # lmer without correlation term     m1 <- lmer(distance ~ (1|Subject) + age + (0+age|Subject) + Sex, data = Orthodont)     VarCorr(m1)     # Groups    Name        Std.Dev.     # Subject   (Intercept) 1.474105     # Subject.1 age         0.099979     # Residual              1.402591      With `lme` I think I can drop the correlation term using the following specification...               # lme without correlation term?     m2 <- lme(distance ~ age + Sex, data = Orthodont, random = list(~ 1 | Subject, ~-1+ age | Subject))     VarCorr(m2)     #             Variance            StdDev         # Subject =   pdLogChol(1)                       # (Intercept) 2.172946296         1.47409169     # Subject =   pdLogChol(-1 + age)                # age         0.009996006         0.09998003     # Residual    1.967260819         1.40259075      I am not entirely convinced that these are the same models though, partly because I can not find any resources that details how to specify this specific form and partly because the output from `print` is a little mystifying to me...               m2      # Linear mixed-effects model fit by REML     #   Data: Orthodont      #   Log-restricted-likelihood: -218.3227     #   Fixed: distance ~ age + Sex      # (Intercept)         age   SexFemale      #  17.5806928   0.6601852  -2.0117005      #      # Random effects:     #  Formula: ~1 | Subject     #         (Intercept)     # StdDev:    1.474092     #      #  Formula: ~-1 + age | Subject %in% Subject     #                age Residual     # StdDev: 0.09998003 1.402591     #      # Number of Observations: 108     # Number of Groups:      #                Subject Subject.1 %in% Subject      #                     27                     27       In particular, what is `Subject %in% Subject` referring to? and why are the residual considered as part of the second random effect term?


title: fitting LMEMs for repeated measures with no correlation between intercept and slope

text: For a simulation study, I contrast the power of different LMEMs for repeated measures. To get p-values, I use likelihood ratio tests where I compare a model including a fixed treatment effect with one that has the same random effects structure but without having the fixed treatment effect. I want to specify a model in which random intercept and slope are allowed to correlate and one in which it is not allowed. But, when I extract the p-values of both models, it appears that they are exactly the same. Why? I am not an expert for mixed models but that seems weird. Here's my code using lmer syntax:               # correlation allowed       ml3 <- lmer(rt ~ treatment + (1+treatment|subject),data=df)       ml0 <- lmer(rt ~ 1 + (1+treatment|subject),data=df)       lrt <- anova(ml3,ml0)       pVal <- lrt$"Pr(>Chisq)"[2]          # no correlation       ml2 <- lmer(rt ~ treatment + (1|subject) + (0+treatment|subject) ,data=df)       ml0 <- lmer(rt ~ 1 + (1|subject) + (0 + treatment|subject),data=df)       lrt <- anova(ml2,ml0)       pVal <- lrt$"Pr(>Chisq)"[2]      Data:               dput(DF)     structure(list(subject = structure(c(1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,      3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,      3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,      4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,      4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,      5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,      6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,      6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L,      7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L,      7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L,      8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L,      8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L,      9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L,      9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 10L, 10L,      10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L,      10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 11L, 11L,      11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L,      11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L,      11L, 11L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,      12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,      12L, 12L, 12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L, 13L, 13L, 13L,      13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L,      13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 14L, 14L, 14L,      14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L,      14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L,      14L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L,      15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L,      15L, 15L, 15L, 15L, 15L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L,      16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L,      16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 17L, 17L, 17L, 17L,      17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L,      17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L,      18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L,      18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L,      18L, 18L, 18L, 18L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L,      19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L,      19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 20L, 20L, 20L, 20L, 20L,      20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L,      20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 21L,      21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L,      21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L,      21L, 21L, 21L, 22L, 22L, 22L, 22L, 22L, 22L, 22L, 22L, 22L, 22L,      22L, 22L, 22L, 22L, 22L, 22L, 22L, 22L, 22L, 22L, 22L, 22L, 22L,      22L, 22L, 22L, 22L, 22L, 22L, 22L, 23L, 23L, 23L, 23L, 23L, 23L,      23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L,      23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 24L, 24L,      24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L,      24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L,      24L, 24L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L,      25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L,      25L, 25L, 25L, 25L, 25L, 25L, 26L, 26L, 26L, 26L, 26L, 26L, 26L,      26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L,      26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 27L, 27L, 27L,      27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L,      27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L,      27L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L,      28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L,      28L, 28L, 28L, 28L, 28L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L,      29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L,      29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 30L, 30L, 30L, 30L,      30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L,      30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L     ), .Label = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10",      "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21",      "22", "23", "24", "25", "26", "27", "28", "29", "30"), class = "factor"),          treatment = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L), .Label = c("1", "2", "3"), class = "factor"),          rt = c(551.798792586772, 693.014255128461, 715.599061613616,          670.119831344829, 777.748610260388, 736.018489208224, 636.791011800404,          864.593711496912, 604.529352905588, 596.673178487122, 860.858066937491,          717.975814131377, 531.672833100059, 571.150454430927, 644.315598150879,          601.914697283216, 583.92746647402, 702.714068138085, 660.346853172676,          541.292786332608, 608.233103066463, 740.593415976325, 686.059921551164,          706.522723402261, 567.648255604935, 596.111352599386, 625.779084220279,          752.776987343973, 922.314285125596, 720.736074757203, 768.585671134519,          539.657760625667, 431.193030969184, 739.341430343149, 581.505474510558,          485.905153431116, 524.085545405872, 876.566370460358, 631.259754679943,          587.887351105621, 624.365050240473, 642.528438460209, 440.661792577731,          517.142782023978, 705.840003729944, 557.122142924839, 645.711579229236,          477.292943229673, 578.522058679457, 623.879658296107, 480.855063147831,          622.295733392922, 611.623490658329, 594.974733982977, 546.239019853272,          551.638287622872, 567.791819285002, 539.239628365136, 541.328446070423,          609.931976806498, 549.492601324081, 479.862984098331, 592.411150981731,          466.224011597179, 489.388878789762, 565.187127159354, 806.196199699478,          565.001071713299, 449.529036961143, 446.824243314547, 357.993777663337,          370.054045045062, 546.443479822161, 473.894296409884, 335.821704077378,          370.498649523398, 486.052525038318, 436.53292033153, 359.637460079864,          333.146018287273, 597.894114487158, 551.993792800734, 518.563432886515,          513.629383189428, 572.062676720248, 366.611317255576, 353.934207291842,          492.273303938824, 414.632984933654, 456.987565377718, 585.524844348671,          453.577328112778, 665.072248688078, 459.204631254183, 452.028605442515,          422.731506299078, 522.84363619892, 771.586286956136, 478.422495080758,          530.925498291748, 457.030618882822, 446.313635696342, 537.708665959068,          815.929601138346, 460.420404065423, 603.027278932425, 538.526470664698,          571.491835856551, 567.234631106499, 450.878624452358, 650.340961680322,          598.829083718722, 620.85411026516, 573.498196791879, 519.953442801483,          1143.14393274202, 505.028670926264, 685.875665196364, 605.316954852204,          645.269429082978, 678.192056130499, 605.671978651269, 650.564580984954,          641.331733928499, 687.164180542278, 613.873849203194, 789.829709495785,          560.793473918547, 707.169378961089, 680.753196215641, 659.262985906231,          600.712008959484, 662.275074291484, 661.346206480403, 568.31000899618,          661.439508442242, 761.227769640367, 699.901658463283, 631.422448673388,          734.257735977184, 585.776345181453, 714.587957176744, 893.931699334816,          586.343993838929, 664.205207596859, 828.003782888565, 906.448165648461,          584.196768113385, 747.575564348236, 687.698668648395, 398.221092516595,          490.332613905338, 502.683377386602, 451.168200674477, 620.606534311108,          369.820458042713, 429.483129392912, 628.153937257066, 476.31856841443,          608.016880363378, 402.588700424079, 460.302916138576, 341.209753425326,          600.165531243784, 454.003777748405, 589.089266888531, 504.033320854066,          399.871492203846, 421.426563579218, 423.375093277487, 587.614013919312,          689.18637317583, 652.10069672704, 553.995320740249, 570.86533170596,          636.399559100471, 801.517490092303, 653.425223465164, 684.914139340214,          639.418654954543, 555.718100869331, 677.062768399616, 579.433200999322,          561.757369869387, 672.316124102021, 701.108131071079, 588.129426947175,          438.090900053591, 562.520435558598, 610.177372103278, 564.672192806652,          552.305226838045, 586.912866128373, 872.521433158083, 654.971253063189,          575.068762782096, 784.098108527601, 648.265348029739, 590.541840435637,          552.569131260877, 554.840084955354, 582.798864712891, 573.196470707737,          512.123960183202, 579.838037093289, 710.216216611067, 779.786949219207,          615.995650564573, 549.096392807351, 600.781394864656, 415.016144351654,          765.924387691343, 401.541419432177, 436.050367769487, 536.508634405116,          445.112952169149, 478.003493101049, 509.819044087032, 490.265270275681,          594.667876389766, 781.844411855516, 827.832351086729, 379.401116898897,          469.280230588986, 397.115839743604, 874.524377877882, 612.130504039819,          802.270319490186, 651.842161968928, 581.489774054855, 686.457677143518,          570.172663082147, 566.996565453736, 577.947675356248, 494.016772046721,          546.065861910691, 506.178677541412, 527.653822550596, 470.043764013502,          595.080116592997, 464.590366280242, 684.362069491853, 534.310814471562,          545.7046301149, 452.141529834992, 619.652055160652, 568.61376011316,          576.847350527713, 514.248803061826, 585.909312171032, 687.892034205561,          907.133281713537, 549.603068658537, 617.860688444804, 423.424246676122,          524.28348263976, 593.203848577403, 431.733188413523, 476.284708033659,          588.88583865225, 437.275988819986, 733.45270912684, 592.366412341047,          606.958434204909, 721.61902078205, 604.596941234802, 534.65440311647,          526.29928462654, 655.076689084035, 560.740728878595, 591.083376633783,          533.331301213643, 750.32841350028, 547.366173885661, 602.313382446308,          787.158938523176, 534.80106549099, 454.37886245909, 599.535859565986,          607.126697674517, 668.173760533712, 589.060272311024, 590.188448587092,          711.910534337354, 528.634489779135, 600.468858484032, 580.666817624455,          659.907090614686, 596.395917159692, 994.163737779338, 662.059444540888,          637.256716085147, 714.353436812361, 587.212427691626, 676.527668439672,          609.004414569998, 667.364031145608, 788.145350832559, 725.891539439069,          561.397498270981, 430.027498616446, 500.437956195847, 463.801763917305,          458.358780003907, 383.304386810731, 598.957692241571, 409.89510543858,          390.650415086637, 552.072907469115, 388.554580583084, 671.244783776164,          433.841093781351, 423.187562794827, 502.566122911232, 469.869008810394,          547.610270179268, 501.091740213331, 417.336826115574, 500.284514580019,          460.835882303962, 650.071068396249, 523.313503950421, 861.366681123829,          879.241985731583, 673.655630620254, 448.199583227711, 578.587129494665,          654.742597624172, 623.62363768736, 665.526175470944, 942.738238156293,          1006.15443845549, 667.679153732234, 711.686114156855, 642.929069350516,          685.862290196822, 1062.63097632175, 758.162511396556, 827.547233897549,          668.688764986398, 791.497544557741, 838.143090686178, 681.935257212825,          758.732997665222, 661.724656922782, 793.560801116029, 896.416624866383,          642.617709357462, 633.832129070135, 751.515360586321, 616.151652306802,          684.496510560379, 655.310039878885, 710.298048482024, 606.373767619465,          754.268924528687, 822.582103710613, 820.556840434073, 733.785411148237,          584.830824784288, 588.316573524589, 572.95505735157, 559.402915982595,          640.891735376065, 482.407652048448, 569.682285396545, 517.277707765673,          698.102946480301, 651.001615070688, 665.691471843539, 511.440973330271,          504.930464361447, 613.891964397534, 454.7073692139, 513.19138352863,          422.708112768038, 347.049510934991, 523.980248957572, 480.301125161823,          633.307276361827, 799.987010744151, 533.354042715484, 410.150445477125,          809.568249688128, 531.41976915349, 792.355614308461, 747.208014043674,          607.571115317023, 546.485408007754, 633.55875460818, 767.73368427773,          676.492693414302, 751.649529779836, 984.189814104173, 684.929427919003,          615.787024482925, 567.942282464503, 571.041675151281, 614.028930539252,          779.839834490734, 630.179209124113, 651.603035032816, 788.591687415382,          799.38918489533, 701.842888543902, 693.932887722425, 624.800556024233,          659.981040765196, 542.243217216484, 721.703181143723, 607.818766172507,          586.813797432694, 610.206652108693, 837.694469363876, 763.535995041537,          758.89830766469, 616.838182390385, 578.107924042397, 628.314074464124,          676.917384461922, 635.824489980127, 495.143374853889, 815.582155744321,          534.740299502999, 568.853739307473, 771.28519095763, 673.064072347686,          713.558399193608, 599.316767121742, 689.796480377465, 673.175516507747,          518.390229271871, 784.139459988779, 536.808895738866, 591.342581848355,          563.762291009613, 679.413099342014, 490.762928348403, 575.612328735691,          504.631070884374, 689.919729220693, 545.809277581445, 641.095314483731,          618.332256267043, 641.913937397485, 953.129874375365, 646.628853366556,          631.881991258933, 671.647395089865, 503.290799393907, 506.506064370266,          529.718096437596, 484.255438291713, 861.643688089666, 625.018895601203,          768.279996867868, 596.708155034627, 671.714642028838, 844.072568247028,          578.694918479722, 422.588349061727, 594.493346157147, 520.812331184257,          741.876339265066, 555.516494731537, 579.386393427601, 622.316950052304,          523.488853303438, 595.901305255518, 553.512680895547, 557.643582245011,          624.788623102115, 479.715363427417, 574.354660431126, 524.472350214463,          660.841590121958, 608.621321258764, 631.743182107793, 711.470012104646,          617.432792370567, 573.363544694191, 617.449976333406, 563.127159530709,          583.25391667852, 678.936105477067, 572.153554376884, 609.829503412847,          606.069768210344, 693.198276061625, 647.952198803514, 508.091779167254,          654.226813385831, 530.92016927824, 504.963210966908, 512.488303835862,          763.325818301401, 664.628862733417, 477.385861339593, 566.148674353306,          578.70068655976, 606.967024346421, 697.460752784057, 662.304772796768,          520.905460930742, 629.14344808993, 815.023764792173, 680.359748369552,          628.317980877129, 679.505810999772, 512.999611466799, 656.728389486035,          548.409794219861, 619.925877003775, 581.949057396067, 663.545400676099,          666.518874913722, 691.701483159255, 616.896649470106, 595.504424960074,          574.172251324537, 552.787259430621, 678.018976276998, 607.329759814185,          611.581207725935, 690.981992177989, 564.832150097104, 540.336710300887,          533.631062681699, 553.612294126468, 479.592612257575, 805.559491265258,          528.304765655223, 496.66528049325, 503.82305630743, 696.245302300331,          566.070769246181, 659.391688013324, 697.528902380277, 524.903347139913,          599.821891499886, 605.43297053286, 663.035359384042, 714.444647793395,          578.371129029246, 615.320894052349, 586.420779403222, 611.255799029828,          583.666658817928, 585.113768358993, 590.122958856932, 629.219469590256,          538.928053998428, 686.894125956954, 600.89266348967, 577.953187882193,          554.413905790583, 579.409330807908, 745.040084235899, 891.087132273406,          552.988614856682, 820.955690300634, 634.469321378978, 673.975047013567,          756.507601731563, 639.35059215201, 722.557986588015, 586.447409643988,          656.432427481585, 682.783381677787, 644.716232116734, 572.306442663379,          517.147784935371, 455.957276558869, 581.819706567048, 484.576715810217,          481.058650198264, 769.846887958231, 614.717393882487, 512.55467514312,          562.454770697369, 470.842224095898, 473.821865893767, 525.60888531351,          504.615217687803, 633.463711580414, 697.478798243637, 542.901502870182,          554.078075963646, 552.734146037028, 505.049122827182, 486.831379133217,          784.789844765716, 625.603196289942, 670.995369035953, 685.579926259636,          540.482850768361, 439.117039069522, 502.605387735171, 683.149979103402,          516.322278257158, 642.3671240847, 552.631029671279, 647.736853458454,          560.328629000192, 730.565596415312, 618.355157539931, 663.62525661894,          762.952777957374, 628.608584740535, 624.968635247218, 661.70575556195,          574.702252033339, 612.175432188694, 517.854558597715, 675.314287039473,          542.173699486536, 627.693153783529, 692.014742774091, 739.358457751625,          594.351746882543, 487.81864701434, 584.028504991851, 772.039898984639,          467.800944704621, 699.587386648698, 711.892383008835, 615.312970618784,          675.467367812567, 504.86811313757, 549.182867476271, 399.814725143066,          558.16544073586, 422.858340371991, 636.197179849367, 452.640724786824,          620.185492727861, 463.138266913543, 513.777642749675, 562.005709606924,          536.148107655772, 609.263894686695, 569.500530985324, 431.43349765191,          584.398797922899, 499.315449753743, 619.406287942484, 530.179301913325,          581.286972325074, 609.256211854971, 603.54936265609, 631.687639186526,          657.738657905408, 662.587176694764, 566.971357910094, 858.130855520899,          636.509916564228, 961.258701590037, 654.287721552112, 722.951283166332,          632.132002720104, 487.318833408862, 568.846556179602, 618.040023485574,          496.276276900436, 575.584711170303, 581.733646148308, 429.189732200854,          542.010860963286, 494.264804962282, 422.488093162063, 458.026827419797,          453.276318995818, 469.307056127784, 569.3897375524, 579.968164508765,          499.796867133562, 434.308570876294, 715.325075398682, 716.351002032214,          628.027210297141, 500.606473063414, 481.830926575354, 613.318935569666,          976.743117685661, 644.275338785824, 478.107026795071, 557.000656140104,          705.452526296914, 485.589962586432, 499.786618070234, 558.712600821937,          590.922333630969, 542.410481740128, 489.828649613243, 588.761989730902,          546.978243344109, 606.196264556647, 564.236942812287, 536.823507783435,          579.521762239388, 515.499050687035, 630.900495035976, 656.332012380511,          587.856990348358, 522.305185738772, 739.045222536055, 730.745657623434,          924.844562132056, 784.35778551794, 602.551997131679, 647.756594982111,          698.734409940058, 836.546703053691, 698.539340777214, 1026.87197547421,          629.904387631378, 1256.38776880042, 550.201102905894, 602.768212743634,          721.187968335008, 978.957227830475, 720.414204972345, 674.653879707098,          788.5848329187, 719.806764303146, 717.541641935441, 692.033638361742,          697.745018761048, 666.284578038868, 772.479906582772, 659.065480010219,          820.683359002167, 676.970146466469, 874.368371442289, 727.177048657942,          698.30875424695, 779.956903736863, 845.624757018358, 713.683982932567,          726.314529163246, 619.106987740244, 602.275574812157, 1012.31009024481,          900.357542354599, 644.181000629658, 696.237940779481, 723.968447714559,          728.433707295704, 738.277840656323, 624.899240530917, 602.62262026982,          741.637155120584, 771.780154658976, 655.858038042129, 746.047425940232,          620.210165076071, 727.335465168732, 575.696484591354, 664.112900876799,          635.568414156208, 813.389693290661, 782.803606862186, 611.54036820155,          626.182943686362, 629.796912199246, 850.589128120044, 751.317837909887,          698.370084520991, 684.731100139964, 675.610584876072, 746.655183232173,          700.916031879367, 672.999065959173, 754.165715908428, 1000.39209579409,          711.494061908177, 641.566680754116, 804.073919378063, 684.705472972499,          703.074833865838, 461.242335766073, 427.275150289403, 448.904300292751,          504.467553858542, 468.445048327652, 492.18239230431, 441.918536235364,          516.654435536123, 552.509951287832, 409.431285285276, 435.346043325067,          547.743500224329, 531.812309405229, 364.578903523756, 508.390034089605,          379.665384727272, 509.53781675453, 357.585528744983, 441.32777083335,          443.431829371385, 606.59281583625, 686.319960511486, 625.7347791126,          738.463790041674, 563.274182025531, 464.717827380926, 507.784383381725,          588.175980498022, 413.597103546393, 718.241261054521, 748.95432032953,          753.3423119569, 624.488936342165, 685.583685719862, 556.272453690569,          553.85940984926, 547.648759204925, 600.572261526898, 798.618024234413,          726.941645252511, 571.72160960877, 550.034219181198, 627.264280360843,          678.989621824996, 643.317665763982, 562.10080966811, 612.803383559254,          681.139896626253, 602.318525361381, 725.533189662524, 624.989203686985,          757.711095633453, 672.222798550503, 814.596048431427, 661.778896020412,          598.861275115565, 720.798422223033, 607.767472082705, 786.918370564599,          609.319159908796, 654.545951651791, 590.381966203786, 580.331158503126,          515.973182460188, 652.661059652561, 545.988014818838, 695.948607980764,          541.582595061958, 589.514921567389, 543.237248780138, 629.781187086802,          495.114814971547, 713.705169086826, 604.33689023605, 515.485770936165,          519.257552917307, 699.375113218082, 538.736700962025, 475.293688282428,          575.221421714611, 504.071893399439, 480.754098260713, 529.183372881182,          600.773688635732, 528.753221287108, 623.467259878089, 524.506017554373,          578.950778202312, 558.484848311201, 548.077100964434, 628.137735553388,          536.855310919075, 671.288977759216, 656.468309447081, 554.138278777839,          652.139385673766, 576.644938824018, 623.229206449598, 690.260009908557,          908.122945156817, 512.745373098672, 512.90318083329, 732.257651656802,          542.85407478119, 497.700590782599, 523.882208542902, 448.976579619777,          563.466660067041, 505.398939062326, 668.814284148356, 512.653186957994,          503.040433891059, 433.619712591384, 464.197386550985, 408.88198638402,          681.920233753602, 504.514813136438, 570.133904166935, 491.416987899975,          440.029552147731)), .Names = c("subject", "treatment", "rt"     ), class = "data.frame", row.names = c(NA, -900L))      Thanks in advance


title: The Bayesian approach to computing estimator bias and variance

text: From what I understand, jackknife and bootstrapping are frequentist methods for computing statistics (bias, variance, etc.) of an estimator. Given a sample of my data and an estimator, and assuming little about the generative process, how can I compute the same statistics of my estimator using a Bayesian approach?


title: What is the "equivalent" of normal distribution in an interval?

text: What is the most "natural" family of distributions on an interval [0,1] indexed by their mean $\mu$ and standard deviation $\sigma$? I am looking for something occurring in "nature", like normal distributions, except that it has to take values from 0 to 1.


title: Interpreting Cox regression

text: My regression results are as follows.               Variable    B          exp(B)     age        0.004575    1.004586     capital    0.2250      1.2524     year      -0.1026      0.9024      Can someone please help me to interpret these? (eg: one year increase in age would increase/decrease the hazard by Xx%)


title: Correlation of Fixed Effects in lmer

text: I'm a psycholinguistics student with few knowledge in statistics and I have some doubts about a Correlation of Fixed Effects in lmer function (lme4 package). So, if my question is stupid… hummm… I'm sorry! My response variable is RT (Reaction Time in a self-paced reading experiment) and my independent variables are Ant (PP, NP) and Verbo (SG, PL). I have modeled the data with intercepts for Sujeitos (the people who are doing the task) and Item (the sentences i've used), asking for principal effects and interactions between the variables. Here is the model `lmer(RT~Ant*Verbo+(1|Sujeitos)+(1|Item))` and that's the coefficients for fixed effects: ![Table Fixed Effects](http://i.stack.imgur.com/MM27e.jpg) So, I have made a table of the coefficients for the interactions. My problem is: that -0.705 and -0.716 correlation effects are a problem for my interaction terms? I'm saying this because the coefficients for the condition NP:SG came from the coefficients of SG only: ![NP:SG](http://i.stack.imgur.com/vsYZo.jpg) and the coefficients for PP:PL came from the PP only: ![PP:PL](http://i.stack.imgur.com/t0f8G.jpg) So, to me, there is no problem here, because I'm not contrasting PP:SG x PP (correlation = -0.705) and I'm not contrasting PP:SG x SG ((correlation = -0.716)). But I do that when getting the coefficients for PP:SG: ￼![PP:SG](http://i.stack.imgur.com/xBlx7.jpg) In the last case, there is a contrast between PP:SG x PP and SG. So, the correlation could be a problem. Is this correct? And, if so, how can I deal with this? I've read some thinks in Jaeger's blog and in this book: `Howell, 2010. Statistical Methods for Psychology`, but it doesn't help much. Thank you.


title: Using Real-Time Forecast Errors to Improve Forecast

text: I'm trying to forecast customer orders, and want to incorporate real-time forecast errors into my forecast. Say it's Monday, and I forecast that the customer will order 100 units on Wednesday. To make things as easy as possible, suppose that he'll just make one order this week. There are a lot of cases I want to cover:   * The order occurs early: If I see an order for 120 units on Tuesday, then I want to adjust my forecast to be 0 units on Wednesday.    * Some of the order occurs early: If I see an order for 60 units on Tuesday, then I want to adjust my forecast to be 40 units on Wednesday.   * The order occurs late: If Wednesday passes, and there are no orders, then I want to adjust my forecast to be 100 units on Thursday.   * etc. I've never seen a problem like this, where the future forecast should somehow depend on past forecast errors. An additional problem is that there's a lot of noise in the underlying process, so even if I could get the timing 100% correct, I'd see a lot of forecast errors. Also, there could be multiple orders per week. I am, however, 100% sure that the order will happen. Has anyone seen a way to deal with problems like this?


title: Detrending a time series regression model

text: In my text book is says: > Regress each of $y_t, x_{t1}$, and $x_{t2}$ on a constant and the time trend > $t$ and save the residuals... What do they mean regress each on a constant? What constant specifically? This is coming from a section in the textbook (introductory econometrics to be exact) dealing with detrending a time series.


title: Kolmogorov-Smirnov test in likelihood function

text: I want to test how well my data fits a uniform distribution and use this as one factor in a likelihood function I am constructing. Unfortunately, I have no solid basis in statistics. So far, I figured out that a good criterion to use is the Kolmogorov-Smirnov test. I can do the test, but I do not know what the resulting number I get means. What I want to have is the value of the PDF describing whether or not my sample is coming from the given distribution, so that I can plug this into my likelihood function. How do I get it? Thank you!


title: When dealing with time series data must the time periods be equal intervals or annual?

text: So like the title states, do the time periods in time series data have to be of equal intervals? For example, if we wanted to have data on the level of education recorded by year, would we need observations for 1996, 1997, 1998,...,2002 or could we have data for 1996, 2000, and 2004? It seems to me that it would need to be of equal spacing intervals, but I just want to make sure that a time interval of say 1997, 1998, 2000, 2003 wouldn't be acceptable.


title: Cointegration difficulties using Stata

text: I am working on a time series analysis with 52 quarterly data concerning a variety of possible determinants of CO$_2$ emissions by transport (CO$_2$ taxes, GDP, load factor, transport volume). This means I will have at least 3 independent variables. I am using Stata and I want to see if there is a long and short term relationship between them so I have started a cointegration process.   1. I have conducted an ADF test for unit roots test and 3 out of 5 variables are I(1) but the other 2 are I(0). I have read here in another answer that in a multivariate model it is possible to continue with the cointegration even if you have I(0) and I(1). Is that true and how does this change the interpretation of the results?   2. Also I have a major problem to separate VAR from VECM. What is the difference between them?   3. Do I have to do first a VAR and afterwards a VECM or just one of them? Which is more accurate?


title: How to interpret basic three variable mediation model?

text: I have run a multiple regression to determine which factors predict the willingness to pay for members of another group. The multiple regression determined that symbolic threat (fear of a cultural and traditional attack on group values by others) was the only significant predictor of how much one would pay for members of another groups. I did however, in my research find that attitudes towards the outgroup in question entirely mediates this effect. The final results are as shown: (I can't post a picture because my reputation isn't 10 yet so I'll write it as best I can, all ---- lines should be solid)                      Attitude toward outgroup              ㅅ                 \              /                   \       .633***/                     \.482***            /                       \           /                         v       Symbolic ------------------> Amount of Tax Willing        Threat   .495*** (.190 ns)    To pay on outgroup.      I am however having difficulty in interpreting what this mediation means.   * Does this mediation mean that without any attitudes towards the outgroup, symbolic threat would no longer predict the amount one would pay for the other group?    * Does this mean that existing attitudes combined with symbolic threat faithfully predict the amount one will pay for the other group, while separately no prediction can be made? Any ideas or useful references would be really helpful thank you.


title: Purposeful selection and confounding

text: I conducted purposeful selection as outlined in Jewell's _Statistics for Epidemiology_. The log likelihood tests showed covariates, which I considered to be confounding though not significant in the modeling, as not significant when excluded. In the end, however, goodness of fit test was significant. When re-introduced the "confounding" covariates, goodness of fit was not significant. Unfortunately, I am unable to find much research in this subject as it is more so for the purpose of policy and program practice so a reference for which variables to hold is lacking. Should I still keep covariates believed to be confounders regardless of the log likelihood test results?


title: Balanced incomplete block design

text: In a BIBD, let $a$ denote the number of treatments; $b$ denote the number of blocks; $k$ denote the number of treatments exactly in each block; $\lambda$ is the number of times each pair of treatments appears in the same block. Each treatment occurs $r$ times. How are the relations $r=\binom{a-1}{k-1}$ and $\lambda=\binom{a-2}{k-2}$ obtained?