title: Advice on how to analyze these data

text: In my experiment, I have 312 separate samples, each of which is made up of several different food items. I multiply each food by a specific multiplier to convert it to energy, and then sum the energy in each sample to get a total energy value for the sample. There are several different multipliers possible to select from for each food, resulting in different energy values. I want to test the energy output per sample using five different methods of selecting multipliers (i.e. max value, min value, etc.). To me, it seemed like multiple paired t-tests made sense in this situation, because the energy values are all dependent on the original food mass in each sample. But I am worried about family-wise error when performing several tests. Is multiple t-tests the best solution? Would repeated-measures ANOVA work better? Thank you!


title: Calculating the BIC for multidimensional, clustered data

text: I have some microarray data (~15 samples) which I've clustered via pam, with a range of cluster sizes and I want to find out the optimal k with BIC. I basically want to re-implement the BIC score from the x-means paper and this stat.stackexchange post answered some basic questions. But it seems that their definition of sigma is for the unidimensional case. How would I calculate the Covariance matrix for my multidimensional dataset to plug into the multivariate Gaussian log-likelihood function? I could be missing something obvious, but I can't seem to find a reference to explain the multivariate case for cluster models. I can add a reproducible example if needed. update: Here's the formula for variance: $$ \sigma^2 = \frac{1}{R-K}\sum_{i}(x_i - \mu_{(i)})^2 $$ Here, $x_i$ is the sample point and $\mu_{(i)}$ is the cluster center for the cluster which the sample belongs to. In the multivariate case, a point is defined by a vector of size $n$ (for example, a row $i$ in the data matrix) so the mean $\mu_{(i)}$ should also be 1 by n vector. How then do they get a single number for the variance? * * * 1 X-means: extending K-means with efficient estimation of the number of clusters, Pelleg & Moore


title: Testing data-generating hypotheses for reaction times

text: My dependent variable are reaction times (RTs) in a psychological experiment. The experiment features a control condition and an experimental condition and for each subject, I have several hundred realizations of each condition (the variance in this experiment is pretty large, both the intra- and inter- individual). I calculate differences by subtracting the RTs in the experimental condition from those in the control condition and in the mean over trials and subjects, I get a value significantly different from zero. Now my problem: From a theoretical perspective, I derive the hypothesis that the experimental condition differs only in some of the trials significantly from the control condition and that the differences in the means comes from this proportion, only. In other words, I suspect that the RT from some trials in the experimental condition is drawn from the same distribution as the control trials and only some are from a different distribution. I could possibly come up with a theoretical estimate of the size of this proportion but an approach that would allow to estimate it from the data would be even better. Is there a way to statistically test this hypothesis against the hypothesis that all experimental trials are from the experimental distribution? **[EDIT]** I guess, I can rephrase the problem as follows: Can I test whether the data from the experimental condition is from a bimodal distribution rather than from a unimodal one and, furthermore, is one of the modes identical to the mode of the control-distribution? **[/EDIT]** thanks a lot, matthias


title: d20: effect of "face-up" value

text: I was wondering whether you have a suggestion on how to program a specific R-program for (a lot of) d20 testing. I have the null hypothesis that the face-up value of any dice (I'm mostly interested in d20, but the mechanic is the same) does not have any effect, positive or negative bias, on the coming roll. The way I would like to make the experiment:   * I create a table with 10 small boxes on it with walls that all can be removed at the same time. I will make a mechanic tipping of the table to simulate a human throw.   * I will write down every result for each dice, input it a spreadsheet and later `scan()` it into R. Here is my problem now: I want R to tell me which face-up value gave the highest value afterwards. This means I need to program functions with twenty variables and then somehow, probably by using `lm()`, make it possible for me to analyze whether the p-value is lower/higher than 0.05. Of course, I could force the face-up value to be the same every time, but I would rather test all face-up values at the same time if possible.


title: Is it a problem that the sizes of randomized groups can differ?

text: My website has advertising for various companies - say 10,000 companies. Each company has it's own page but all the pages are currently similar except by which company is on the page. A success/response/dependent variable is a visitor clicking on the advertisement. Each company's page has a different success rate - whereby 50% of visitors click for some company pages, and 1% for others. I am testing a new page layout that will be for every company. We have 4 different potential new pages, for a total of 5 pages being tested including the control. A randomization algorithm will assign each visitor into one of the 5 buckets, however this does not guarantee that each bucket will get an equal representation of stores. Whereby, one bucket could get more visitors to stores with the 50% success rate, or vice versa. How big of an issue is this? A possible solution is to serially (not randomly) assign visitors into buckets based on the store. So for Store A, visitor 1 to that store gets in bucket 1, visitor 2 bucket 2, and all the way to 5, whereby visitor 6 goes to bucket 1 again. Is this necessary or even statistically valid or does it create bias? Thanks - i am happy to explain more or even to summarize this in a more universal format if that is preferred.


title: Intelligence Squared Scoring and Winner Determination

text: There is an NPR podcast called Intelligence Squared. Each episode is a broadcasting of a live debate on some contentious statement such as "The 2nd amendment is no longer relevant" or "Affirmative action on college campuses does more harm than good". Four representatives debate-- two for the motion and two against. To determine which side wins, the audience is polled both before and after the debate. The side that gained more in terms of absolute percentage is deemed the winner. For example:                         For    Against  Undecided      Before   18%      42%       40%      After    23%      49%       28%           Winner: Against team -- The motion is rejected.      Intuitively, I think this measure of success is biased and I am wondering how one would poll the audience to determine the winner in a fair way. Three issues I immediately see with the current method:   * At the extremes, if one side starts with 100% agreement, they can only tie or lose.   * If there are no undecided, then the side with less initial agreement can be viewed as having a larger sample size from which to draw.   * The undecided side is not likely to be truly undecided. If we assume that the two sides are equally polarized, it seems our prior belief about the undecided population should be $\text{Beta}(\text{# For}, \text{# Against})$ if each was forced to take a side. Given that we have to rely on audience polling, is there a more fair way to judge who wins?


title: Can I write estimate $\pm$ its standard error?

text: Suppose I have an MLE estimate $\hat{\theta}$ for a parameter $\theta$, and $\hat{\sigma}$ is the sqrt of the inverse of the negative of the Hessian of the log likelihood at $\hat{\theta}$. Can I write $\hat{\theta} \pm \hat{\sigma}$? In general, if I have an estimate $\hat{\theta}$ for a parameter $\theta$, and an estimate $\hat{\sigma}$ of the standard deviation of $\hat{\theta}$, can I write $\hat{\theta} \pm \hat{\sigma}$? Will others understand what I mean? Is it also what others write?


title: What Gaussian process covariance function captures an affine mean?

text: In the documentation for GPML, the author trained a GP with an affine mean function and isotropic squared exponential covariance function. Then there is an exercise to the reader: > Try training a GP where the affine part of the function is captured by the > covariance function instead of the mean function. That is, use a GP with no > explicit mean function, but further additive contributions to the > covariance. My question is: **How do I design the covariance function in order to capture affine mean?**


title: Why is Poisson regression different with glmer and gamlss?

text: I have a set of count data that seems to fit "Poisson" = not overdispersed, alpha = 0. The problem is, I get different results using `gamlss` vs `glmer`. Any help explaining the difference would be appreciated: `Treatment` is a factor with 6 levels. `Trial` is a random intercept factor with 3 levels.               Glmer.POI <- glmer( Y ~ Treatment + offset(log(Total)) + (1|Trial),      family=poisson))     Gamlss.PO <- gamlss(Y ~ Treatment + offset(log(Total)), random=~1|Trial, family=PO)          # Glmer result:  AIC = 284. deviance = 270     # Gamlss result: AIC = 388. deviance = 376      I get different coefficients and different p-values with the 2 models. When I look at the diagnostic plots, I see residuals that increase with fitted values in the gamlss model "heteroscedasticity", but not with the glmer. My intuition tells me this is related to how random effects are handled, but how? I think I had similar problems with another dataset using negative binomial. I pulled the question from this site to concentrate on this one: If gamlss and glmer behaviour with random effects explains both problems, Great! But I'd like to know how a biologist (non statistician) "poor soul" like me can understand how to use these tools and choose the right one! ;-)


title: Best practices for extrapolating data

text: I have a set of variables that parameterize a logistic equation bacterial growth model. The parameters change based on temperature (e.g., growth speeds up at higher temperatures) and so it is necessary to fit some secondary models in order to be able to estimate growth under a variety of temperatures. Only five temperatures (16, 20, 25, 30 and 35 C) have been measured under laboratory conditions, and so the secondary equations are fit from these five sets of parameters. I have noticed that the growth model is very sensitive to slight changes in the parameters - even the small difference between the observed and predicted parameters resulted in a large end difference in growth. So it is very desirable to include the real observed set of parameters in the model. I can do this using a higher order polynomial and I am comfortable with the values that are interpolated in between (because at the very least they are probably close, and the model then includes the exact values at temp=16, 20, etc). However, to make a long story long, once I extrapolate below the lowest observed values (with this higher order polynomial) things get crazy. I know extrapolation is always kind of dicey in terms of uncertainty, but I do need to have values for parameters down to 10 C. Are there some methods that tend to work better than others? Or what could I even search for? (This is a new subject to me and wikipedia's entry on extrapolation mostly suggests "don't do it".) If all else fails I can hold the parameters constant from 10 - 16 C (which is conservative) but it would be nice to use something better. Edits: In order to elaborate, the logistic equation I am using is: $y(t)= \frac{\alpha }{1+e^{\beta-rt}}$ $\alpha ,\beta$ and a $r$ are my parameters that change with temperature. Here are the parameters that were fit to growth curves measured in a laboratory:               Temperature  α        β         r     16          3.2366    3.4698    0.0297     20          3.5557    3.3770    0.1256     25          3.6941    3.9876    0.3300     30          3.9081    4.0331    0.5422     35        3.9868      3.8333    0.8405      And I would like to fit curves to these to know (for instance) what these three will be at 18 C. It is standard in the microbiological literature to use certain types of relationships but I find that they're incredibly inaccurate and I trust the laboratory observed values the most so I'm willing to use something ad hoc. Here is an image of what it looks like for $\beta$ using a polynomial ![enter image description here](http://i.stack.imgur.com/EU8sG.png) And, here is a truncated image from the paper where I got the original data just to give you an idea (Yang et al. 2009. Food microbiology 26: 606-614). ![enter image description here](http://i.stack.imgur.com/y6ejl.png)


title: How To Fix a Non-Representative Sample Using Ordinal Logistic Regression and Predict Appropriately?

text: After searching for this question, I did find this -- but it didn't seem to be asking the same question and I'd like to extrapolate on it if possible to get more into fixing it generally and the prediction side of things. If it's a duplicate, I'm happy to take it down if somebody can direct me elsewhere. In my dataset, I have four assessments I'm trying to predict: 1 [Good] to 4 [Bad]. In my sample, I have:   * 5% '1's   * 59% '2's   * 14% '3's   * 22% '4's My model seems to be working using the `polr` function -- though I'm not convinced of its accuracy (besides the fact that when I redo it, I should probably remove `sex`). Here's the output:               > y <- as.factor(assessment)     > x <- cbind(race, partisanship, sex, age)     >      > olr <- polr(y ~ x, mydata)     > summary(olr)          Re-fitting to get Hessian          Call:     polr(formula = y ~ x, data = mydata)          Coefficients:                      Value Std. Error t value     xrace          0.49485   0.214426  2.3078     xpartisanship -0.00990   0.002942 -3.3654     xsex          -0.21304   0.299763 -0.7107     xage           0.01486   0.006812  2.1819          Intercepts:         Value   Std. Error t value     1|2 -1.4763  0.8253    -1.7887     2|3  1.8049  0.8237     2.1913     3|4  2.4739  0.8290     2.9842          Residual Deviance: 667.1306      AIC: 681.1306      (1401 observations deleted due to missingness)     >      > (ct <- coef(summary(olr)))          Re-fitting to get Hessian                               Value  Std. Error    t value     xrace          0.494853619 0.214426447  2.3078012     xpartisanship -0.009899843 0.002941635 -3.3654220     xsex          -0.213044564 0.299762649 -0.7107108     xage           0.014863269 0.006812136  2.1818808     1|2           -1.476276983 0.825333331 -1.7887039     2|3            1.804900438 0.823656309  2.1913272     3|4            2.473925644 0.829013589  2.9841798     > pval <- pnorm(abs(ct[, "t value"]), lower.tail = FALSE) * 2     > (ct <- cbind(ct, `p value` = pval))                          Value  Std. Error    t value      p value     xrace          0.494853619 0.214426447  2.3078012 0.0210101944     xpartisanship -0.009899843 0.002941635 -3.3654220 0.0007642666     xsex          -0.213044564 0.299762649 -0.7107108 0.4772634401     xage           0.014863269 0.006812136  2.1818808 0.0291183306     1|2           -1.476276983 0.825333331 -1.7887039 0.0736625078     2|3            1.804900438 0.823656309  2.1913272 0.0284281263     3|4            2.473925644 0.829013589  2.9841798 0.0028433959      When I look at the prediction output for this dataset using the model, it gives me probabilities very similar to the percentages I mentioned above with slight variation:               > head(predict(olr, mydata, type = "p"), 10)             1         2         3         4     1  0.07743232 0.6132735 0.1227231 0.1865711     2  0.08238534 0.6225283 0.1185267 0.1765596     3  0.05479841 0.5518965 0.1440269 0.2492782     4  0.05490456 0.5522788 0.1439214 0.2488952     5  0.05753405 0.5614067 0.1413146 0.2397447     6  0.03093548 0.4283424 0.1645377 0.3761845     7          NA        NA        NA        NA     8  0.04295889 0.5013185 0.1555760 0.3001466     9  0.01934162 0.3248235 0.1618940 0.4939409     10         NA        NA        NA        NA      **Question 1:** The main problem with this is that the data I imported is not fit to any representative sample; it's just all the data I have. My initial thought process was that it didn't have to be representative of what I was trying to predict, because I was just trying to give the model more data to accurately predict '1's,'2's, '3's, and '4's. But I was apparently wrong in thinking this, since its predictions are mostly all '2's, because 59% of what I loaded in were '2's. Is there a better way to do this? **Question 1.5:** Is this happening because of the non-representative sample, or is it more because I'm trying to predict the same dataset that I used to make the model with? Thanks for the help!


title: Comparing Percentages Between 2 Groups When Each Group's Total Percentage Exceeds 100

text: I am comparing the rate of specific Response Patterns A through L among students in Group 1 vs. Group 2 (rates for the first five patterns are exemplified in the table). These are response patterns that a research assistant and I coded from student interviews. I also created a code for Multiple Patterns, as I wanted to also compare the rate at which each group of students demonstrates more than one pattern. The problem is that I did not want to make this category mutually exclusive of the others; hence, the percentages across Patterns A through L for each group have different sums, both exceeding 100% (117.23 for Group 1 and 121.94 for Group 2). Each percentage in the table was calculated from the number of times that the pattern showed up in the group divided by the total number of interview items answered by the group as a whole. I calculated it this way because although my interview consisted of up to 11 items per student, some students were administered less (e.g., 8 or 9 items) due to interruptions. Each student's response to each item was coded as demonstrating one or more of the 12 response patterns. **My question is: Would it be valid for me to make statements like "Group 2 exhibited Response Pattern E more frequently than Group 1" given these complications? If not, what is another way to make comparisons between the groups for each pattern?** To obtain the percent for each cell and have each group add up to 100, the proportion could be taken out of the total number of codes assigned to each group instead of total number of items administered across each group, but I'm not sure if it makes sense to do this. I know that I cannot make statements about whether such differences are significant using my current set-up, but I am okay with that for now because this is a more qualitative element of my dissertation research, which already contains other quantitative elements.                          % of Group's Responses Showing Pattern                       Group 1     Group 2                      (n = 22)    (n = 21)      Pattern:   A _______________15.79______15.82   B ________________6.70_______6.63   C _______________16.75______20.92   D _______________12.44_______5.31   E _______________6.22________23.98   F   G   H   I   J   K   L   Multiple Patterns ___16.75_______20.92


title: Equivalent mixed models yielding different results in SAS

text:
Below I show two equivalent ways to write a mixed model in R and SAS. The two R models as well as the two SAS models yield the same estimates of the random and fixed effects and the same standard errors of the estimates of the fixed effects. But the two SAS models do not give the same confidence intervals of the fixed effects. My question is: _which are the correct confidence intervals ?_ Here are simulated data:               library(mvtnorm)     I <- 3      J <- 6      K <- 5      n <- I*J*K     set.seed(666)     tube <- rep(1:J, each=I)     position <- rep(LETTERS[1:I], times=J)     Mu_i <- 3*(1:I)     Mu_ij <- c(t(rmvnorm(J, mean=Mu_i, sigma=diag(I)+2)) )       tube <- rep(tube, each=K)     position <- rep(position, each=K)     Mu_ij <- rep(Mu_ij, each=K)     dat <- data.frame(tube, position)     sigmaw <- 2      dat$y <- rnorm(n, Mu_ij, sigmaw)     dat$tube <- factor(dat$tube)     > str(dat)     'data.frame':   90 obs. of  3 variables:      $ tube    : Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
          $ position: Factor w/ 3 levels "A","B","C": 1 1 1 1 1 2 2 2 2 2 ...      $ y       : num  2.76 5.5 2.54 1.56 6.46 ...     > head(dat)       tube position        y     1    1        A 2.759443     2    1        A 5.496689     3    1        A 2.540150     4    1        A 1.558261     5    1        A 6.462050     6    1        B 4.239749      The corresponding nlme model is the following:               > # firstly set position C as the "intercept" for concordance with SAS     > dat$position <- relevel(dat$position, "C")     > # load nlme     > library(nlme)     > # fit the model     > ( fit1 <- lme(y ~ position, data=dat, random= list(tube = pdCompSymm(~ 0+position ))) )     Linear mixed-effects model fit by REML       Data: dat        Log-restricted-likelihood: -199.0196       Fixed: y ~ position      (Intercept)   positionA   positionB         8.526544   -4.800535   -3.322507           Random effects:      Formula: ~0 + position | tube      Structure: Compound Symmetry               StdDev   Corr            positionC 1.892433                 positionA 1.892433 0.082           positionB 1.892433 0.082 0.082     Residual  1.932750                      Number of Observations: 90     Number of Groups: 6       This model is equivalent to a 2-way ANOVA with mixed effets (in the sense that the marginal models are the same), which is more easy to fit with lme4:               > library(lme4)     > lmer(y ~ position + (1|tube) + (1|tube:position), data=dat)     Linear mixed model fit by REML      Formula: y ~ position + (1 | tube) + (1 | tube:position)         Data: dat       AIC BIC logLik deviance REMLdev      410 425   -199    402.5     398     Random effects:      Groups        Name        Variance Std.Dev.      tube:position (Intercept) 3.28587  1.81270       tube          (Intercept) 0.29543  0.54354       Residual                  3.73552  1.93275      Number of obs: 90, groups: tube:position, 18; tube, 6          Fixed effects:                 Estimate Std. Error t value     (Intercept)   8.5265     0.8493  10.039     positionA    -4.8005     1.1595  -4.140     positionB    -3.3225     1.1595  -2.866          Correlation of Fixed Effects:               (Intr) postnA     positionA -0.683            positionB -0.683  0.500      Below I check that the results are indeed the same:               > # same standard erros     > sqrt(diag(fit1$varFix))     (Intercept)   positionA   positionB        0.8493533   1.1594505   1.1594505      > # the total variance in the second model is given in the first model:     > sqrt(1.81270^2+ 0.54354^2)     [1] 1.892437      Well. Now here are the two equivalent SAS models:               /* First model */     PROC MIXED DATA=dat ;     CLASS POSITION TUBE ;     MODEL y = POSITION ;     RANDOM POSITION / subject=TUBE type=CS ;     RUN; QUIT;          /* Second model */     PROC MIXED DATA=dat ;     CLASS POSITION TUBE ;     MODEL y = POSITION ;     RANDOM TUBE TUBE*POSITION ;     RUN; QUIT;      Results are identical to the R results. But SAS assigns different degrees of freedom for the fixed effects and consequently gives different confidence intervals, as shown below. The first model gives degrees of freedom 5, 10, 10:               Effect  position    Estimate    StandardError   DF  Alpha   Lower   Upper     Intercept       8.5265  0.8494  5   0.05    6.3432  10.7099     position    A   -4.8005 1.1595  10  0.05    -7.384  -2.2171     position    B   -3.3225 1.1595  10  0.05    -5.9059 -0.7391     position    C   0   .   .   .   .   .      whereas the second model gives degrees of freedom 15, 15, 15:               Effect  position    Estimate    StandardError   DF  Alpha   Lower   Upper     Intercept       8.5265  0.8494  15  0.05    6.7162  10.3369     position    A   -4.8005 1.1595  15  0.05    -7.2718 -2.3292     position    B   -3.3225 1.1595  15  0.05    -5.7938 -0.8512     position    C   0   .   .   .   .   .


title: Consequences of identical distribution of random variables

text: Consider two real-valued random variables $Y_1: \Omega \rightarrow \mathbb{R}$ and $Y_2:\Omega \rightarrow \mathbb{R}$. Even if they have the same domain and codomain, $Y_1 $ and $Y_2$ generate different $\sigma$-algebras, $\sigma(Y_1)$ and $\sigma(Y_2)$, and are defined in different probability spaces, $(\Omega, \mathcal{F}_{Y_1}, \mathbb{P}_{Y_1})$ and $(\Omega, \mathcal{F}_{Y_2}, \mathbb{P}_{Y_2})$. Suppose that $Y_1$ and $Y_2$ are identically distributed: $\mathbb{P}_{Y_1}(Y_1 \in E)=\mathbb{P}_{Y_2}(Y_2 \in E)$ $\forall E \in \mathcal{B}(\mathbb{R})$, i.e. $\mathbb{P}_{Y_1}(\omega \in \Omega | Y_1(\omega) \in E)=\mathbb{P}_{Y_2}(\omega \in \Omega | Y_2(\omega) \in E)$ $\forall E \in \mathcal{B}(\mathbb{R})$. Could you help me in saying whether the following statements are true? (1) Identical distribution does not imply that $\mathbb{P}_{Y_1}$ and $\mathbb{P}_{Y_2}$ are the same functions because the sets $\{\omega \in \Omega | Y_1(\omega) \in E \}$ and $\{\omega \in \Omega | Y_2(\omega) \in E \}$ could be different. (2) Identical distribution does not imply that $Y_1$ and $Y_2$ are the same function, i.e. $Y_1(\omega)=Y_2(\omega)$ $\forall \omega \in \Omega$. Therefore, even if $Y_1$ and $Y_2$ are identically distributed, $\sigma(Y_1)$ and $\sigma(Y_2)$ are different. (3) Identical distribution implies that the image of $Y_1$ and $Y_2$ is the same.


title: Clustering algorithm for my situation?

text: Here is my situation. I have a corpus of over 500,000 news. Now I need to cluster the news based on closeness in time and cosine similarity, using vector-space model and TF-IDF weights. I want to cluster the news that report the same event. Clustering performance is beyond time efficiency. I think about using a slide window, only clustering news "week per unit"(making an assumption that one event is reported within a week). Any ideas for this problem? Looking forward for the a suggestions:) Thanks in advance.


title: Test for stationarity in an unbalanced panel

text: I have an unbalanced panel model and I need to check it for stationarity. So, I need to perform a Unit Root test (I think I will use a Fisher Type Test?). But I am a bit confused whether (1) I need to perform a panel unit root test separately for each variable I have in my model; or (2) I need to perform a kind of a combined panel unit root test for all variables in the model (maybe including them into a group or something like that...). (3) How I determine the lagged differences in the command for Fisher type test in Stata? (4) And then if I have a unit root, how can I make the model stationary?


title: Regression to the mean - indication and correction

text: I analyze the correlation between a change score in a variable A (A_t2-A_t1) that has been measured twice (A_t1; A_t2) and another variable B that has been measured once (at t1). 1.) Is it OK to assess whether there might be regression towards the mean by correlating the change score (A_t2-A_t1) with the baseline value (A_t1)? (and proceed with 2 if there is a negative correlation) 2.) To rule out that the significant correlation between B and change score (A_t2-A_t1) is only due to regression to the mean, which of the following tests is correct a) partial correlation between B and change score (A_t2-A_t1), correcting for A_t1 b) partial correlation between B and A_t2, correcting for A_t1 The results seem to be identical, but is one more correct than the other? Could you maybe provide a citation of a paper?


title: one-way MANOVA effect size interpretation

text: I'm having problem understanding effect sizes in one-way MANOVA model. In my case, I'm having a 13 variables which I used for clustering and then the produced `clusterID` variable with 5 levels serves as IV together with other four DVs(count measures) in one-way MANOVA model. I'm using R so my model is implemented something like               manova.fit <- manova(cbind(V1, V2, V3, V4) ~ ClusterID, data=mydata)     summary(manova.fit, test="Pillai")        MANOVA          Df      Pillai    approx F      num Df      den Df      Pr(>F)        cluster           4       0.761        4.29          16         292           0      Residuals          73                                                                   Then, I'm calculating effect size as               etasq <- etasq(manova.fit, anova=F, partial=T)                          eta^2     ClusterID        0.1986109      Now, I have several questions:   1. **First of all, are partial $\eta^2$ values the same as $\eta^2$ in the case of oneway MANOVA?** In the code above, it didn't matter whether I've put `partial=T` or `partial=F`, it always produced the same result. I know it is true for oneway ANOVA, but here they say that _"Kondo-Brown and Fukuda (2008) correctly chose to use partial eta2 because their design was a MANOVA, which by definition involves non-independent or repeated measures"_.   2. Also, I'm not sure how to interpret it. In [1], pg 619, there is explained that _"when there is one degree of freedom in the numerator, the square root of $\eta^2$ equals $r$..."_ and in addition that _"When there is more than one degree of freedom in the numerator, $\eta^2$ equals $R^2$"_. In my case R code returns that I have num DF = 16, so the question is **can I interpret it as $R^2$?**   3. If I can interpret it as $R^2$, **would it be considered as small, moderate or large effect size, according to [2]?** There is a formula on page 284 $$ f = \sqrt{\frac{\eta^2}{(1-\eta^2)}} $$ and if I plug in the numbers I get $f=0.497$ which is in [2] considered a large effect size which is _"not typical for social science research"_ which scares me and makes me question my data analysis :D Thanks, Vita [1] Levine, T. R. and Hullett, C. R. (2002), Eta Squared, Partial Eta Squared, and Misreporting of Effect Size in Communication Research. Human Communication Research, 28: 612–625. pdf [2] Cohen, J. (1988), Statistical power analysis for the behavioral sciences.


title: Estimating the non- parametric conditional probability

text: I have a set of observation from two parameters, let say $x$ and $y$ and then I want to make the conditional probability of $x$ for the given $y$, $p(x|y)$. So first I use `scipy.stats.gaussian_kde` to obtain the joint probability distribution of $p(x,y)$ and use **normalized histogram** and **one dimensional interpolation** to obtain $p(y)$ and then in the grid points of $x$ and $y$ plane, I compute $\frac{p(x,y)}{p(y)}$ in order to get $p(x|y)$. I also want to measure cumulative distribution of the conditional probability in the grids. Is it the correct approach to estimate the **_conditional probability_** or there are other efficient and faster approaches?


title: kNN randomization test in R?

text: I am trying to determine significant differences between groups of data using a k-nearest neighbor randomization test in R. This test basically looks for similarity amongst homogenous groups and separates them out using random clustering. In the literature, this test is called a "K-nearest neighbor (kNN) randomization test," however, I'm not certain if it called by other names elsewhere. For my specific data, I have isotopic ratios given for various prey items. I have already grouped together ecologically similar prey item types, and now want to see if those groups differ from one another in their isotopic signature. This is where the kNN test would come in. Thanks for all of your answers - I'm new to this site, so I'll address your inquiries as applicable in the comments section.


title: Binomial logistic regression in SPSS using survey weights

text: I am running a logistic regression in SPSS with a sample that uses survey weights. The sample size is 1000 and the weights are along the lines of .86 or 1.23 depending on the case. I am using the weights option in the `Data` command to weight my sample before I run my logistic regression but then I get the unweighted frequencies in the `Categorical Variables Coding` part of the output. At the beginning my output gives me a `Case Processing Summary`, which says `Unweighted cases` with a superscript "a" which then says "If weight is in effect see classification table for the total number of cases." My classification tables all are properly saying 1000 cases which is right. I have 7 predictor variables (all categorical) and one outcome variable (dichotomous). My output seems to make sense, but I am troubled that the frequencies are all the unweighted frequencies. Are my results still legitimate? I am afraid that I may not understand the weighting quite properly.


title: Standards on presenting graphs and charts in vector format?

text: In his up-voted answer to _Software needed to scrape data from graph_ , Alexey writes (with respect to the need for data scraping off of raster graphs) "But nowadays the good practice is to publish graphs in vector form." To rip off my comment to Alexey from that page: **Is there a good reference for best practices around which vector format(s) to present graphs of data in?** For example, ought I use an eps encapsulation of an svg file in my LaTeX manuscripts, or am I supposed to output graph in LaTeX directly? Or something else?


title: Differences between Multi-layer NN, Hopfield, Helmholtz and Boltzmann machines

text: There are several resources describing these types of networks, but none of them explicitly address the differences between them. What are the differences between theses models?


title: Is the sign test appropriate to test whether values in one distribution are generally lower than the other?

text: Generally, the sign test is used to test the hypothesis that the difference median is zero between two continuous distributions (Sign test). I am trying to understand whether it can be adapted to prove that the values from one distribution are generally lower than the other. Or are there **better tests** which can be designed for this purpose? Data available: Suppose that I have the average GMAT scores from 50 universities. There are 2 averages per university, the first one being the average of those candidates who have studied mathematics at a graduate level while the second one is for those who have not. A simple plot shows that students who have studied math score higher. However, I would like to prove this fact statistically using a test (and a p-value?). Conversely, I may also like to prove in another case that there is no difference. ![Plot of average scores by universities separated by Math NoMath](http://i.stack.imgur.com/xCFKR.png)


title: REML use in GLMM

text: I would like to make an enquiry regarding the use of ML or REML for GLMM. At the moment, I'm focusing on doing model selection to evaluate the effect of four fixed effects and one random effect on a response variable (continuous data) in a study of behavioural ecology. My global model is close to a Gaussian distribution after log-transformation, but does not completely satisfy the assumption of normality of residuals. Therefore, I'm using GLMM (family = Gaussian) rather than LMM. My question is: should I use ML or REML as a method of model fitting? I'm quite new to R and model selection, but from what I have read in several books (e.g., Generalized, Linear and Mixed Models, by McCulloch and Searle), I understand that REML should only be used to select random effects for a constant set of fixed effects, so I'm inclined to use ML for model estimation. I'm not 100% sure if that's the right thing to do. The default in LMM is set to REML = TRUE, but I'm not sure what is the best method of estimation for GLMM.


title: Is there a "generalized least norm" equivalent to generalized least squares?

text: In a standard regression problem \begin{equation} \mathbf{y} = \mathbf{X} \beta + \mathbf{e} \ , \end{equation} the solution to $\beta$ when the system is overdetermined is $\hat{\beta}= \left(\mathbf{X}\Sigma^{-1}\mathbf{X}\right)^{-1}\mathbf{X}^T\Sigma^{-1}\mathbf{y}$, where $\Sigma = \sigma^{-2} \operatorname{Var}(\mathbf{e})$. The solution to an underdetermined system when $\Sigma\equiv \mathbf{I}$ is $\hat{\beta} = \mathbf{X}^T\left(\mathbf{X}\mathbf{X}^T\right)^{-1}\mathbf{y}$. Can a similar formula to the least squares case be written when $\Sigma \ne \mathbf{I}$ such that $\hat{\beta} = \mathbf{X}^T\Sigma^{-1}\left(\mathbf{X}\mathbf{X}^T\Sigma^{-1}\right)^{-1} \mathbf{y}$? I have not been able to find any good references on this system and would appreciate if you could suggest any also.


title: Estimate the second moment of a latent variable using a conditionally unbiased proxy

text: **The Setup:** Let $X_t$ denote an unobservable stochastic sequence where the first two _unconditional_ moments are finite constants; ie $\mathbb{E} X_t = \mu < \infty$ and $\mathbb{E} X_t^2 = \gamma < \infty$. Assume that $X_t$ obeys some minimum set of regularity conditions such that a weak law of large numbers applies to its first two sample moments; ie conditions that limit the dependence in the sequence and bound the higher order moments. Let $Y_t$ denote a conditionally unbiased proxy for $X_t$. That is, $\mathbb{E} [Y_t | X_t] = X_t$, and hence $\mathbb{E} Y_t = \mathbb{E} X_t = \mu$. Assume that $Y_t$ also obeys a minimum set of regularity conditions such that a weak law of large numbers applies to its first two sample moments. Let $S_T = T^{-1} \sum_{t=1}^T Y_t$. Given the assumptions on $X_t$ and $Y_t$, it follows that $\mathbb{E} S_T = \mu$ and $S_T \overset{\mathbb{P}}{\rightarrow} \mu$, hence $S_T$ is a consistent, unbiased estimator for $\mu$ **My Question:** Given the setup, how can I estimate $\mathbb{E} X_t^2 = \gamma$? **Things I Know:** From Slutsky's theorem, $(S_T)^2 \overset{\mathbb{P}}{\rightarrow} (\mathbb{E} X_t)^2 \neq \mathbb{E} X_t^2 = \gamma$, via Jensen's inequality, so in general $(S_T)^2$ will be a poor estimator for $\gamma$. Also, $T^{-1} \sum_{t=1}^T Y_t^2$ may not be a good estimator since the assumption set does not guarantee that $\mathbb{E} X_t^2 = \mathbb{E} Y_t^2$. **Final Note:** I would be happy with any answer that improves upon the properties of $(S_T)^2$. So an approximation type estimator obtained via some sort of series expansion would be great. That's what I've mainly been looking into thus far, but with little luck. Even suggestions regarding expansions I might look into would be most welcome. Thanks in advance to all responders.


title: Difference between -c and -wi paramters of LIBSVM

text: I have an unbalanced data and I want to use LIBSVM from MATLAB. LIBSVM has two different parameters for cost: -c and -wi. What is the difference between them?


title: Testing statistically significant difference between samples

text: I have 1 variable measured for 8 samples. so, for e.g. height of 8 men. I want to say something about: if height of men is significantly different amongst them. How can I do this? Without any reference theoretical mean, I can not do a simple t-test.


title: Definition of simultaneous equation modelling

text: I have been asked to comment on a statistical analysis proposal that suggests using simultaneous equation modelling as a method. I have a basic general knowledge of statistics (Masters level) but very little experience with econometrics-specific modelling. Does this technique go by any other name? If it is solely used in econometrics, I'm looking for a somewhat non-theoretical applied definition. I'm also looking for an indication of which software it can be implemented in out of a choice between Stata/SAS/SPSS. I have looked through various online sources and am either not convinced of their accuracy or they present the information too theoretically.


title: LibSVM cost weights for unbalanced data doesn't work

text: I have a dataset where the number of negative labeled values is 163 times the number of positive labeled values. That is, I have an unbalanced data set. I tried:               model = svmtrain(trainLabels, trainFeatures, '-h 0 -b 1 -s 0 -c 10 -w1 163 -w-1 1');     [predicted_label, accuracy, prob_estimates] = svmpredict(testLabels, testFeatures,                                                               model, '-b 1');      Accuracy was nearly 99%. I searched and found a post (#7) [dead link?] where it says: > have you tried weighting on a smaller scale (ie: <1) I changed my model to:               model = svmtrain(trainLabels, trainFeatures, '-h 0 -b 1 -s 0 -c 10 -w1 0.5 -w-1 0.003');     [predicted_label, accuracy, prob_estimates] = svmpredict(testLabels, testFeatures,                                                               model, '-b 1');      I have still high accuracy every time because of unbalanced data. Any ideas? **PS:** I am trying to implement the first challenge of KDD Cup 2008 - Breast Cancer. **I want to rank the candidates by decreasing order.**


title: Variable Selection for a binary classification problem

text: I have a binary response variable (0,1) with 300 (approximately) normally distributed variables as predictors all on (0,1). The training set has only 250 observations and the test set has approx. 20k observations. The goal is to make a model to best classify the test set observations as either 0 or 1. (it's an approximately balanced training set, with around 125 each 0 and 1). With 300 variables and only 250 observations, I'm attempting to do a little bit of variable reduction. First I simply checked the correlation between the predictor and each of the response variable, to see if there were any variables that were not terribly correlated with the response (relative to the others): ![enter image description here](http://i.stack.imgur.com/xgubK.png) there doesn't seem to be any real break where I could reasonably say these variables are much more correlated with the response and i'll toss the rest. So I tried a PCA analysis:               train_pca <- prcomp(xVar)     plot(sort(train_pca$sdev), main = "PCA$sdev graph")     sort(apply(abs(train_pca$rotation[,1:25]), 1, sum))         plot(sort(apply(abs(train_pca$rotation[,1:25]), MARGIN = 1, sum)), ylim = c(0.7,2.1), main = "PCA Rotation Graph")      I graph the rotation of each of the vectors to see again if there are any that are much greater that explain the variation than the others. ![enter image description here](http://i.stack.imgur.com/MA2yE.png) Again not any big break where I can decide that certain variables are better than the others. Lastly I built a random forest with the importance control on:               randy_1 <- randomForest(y = yVar, x = train[,5:304], importance = TRUE, ntree = 1000)      and made an importance graph to look at the variables: ![enter image description here](http://i.stack.imgur.com/6EJBk.png) There appear to a few variables that appear better than others, but if you look at the mean gini reduction for all the variables, there really is just a uniform descent. So basically at this point I'm more or less looking to either (1) use all the variables which isn't right (I don't think), b/c the question I've been posed suggests that variable reduction was necessary, or (2) pick a point that seems sort of arbitrary to cut off which variables to use and which to not consider (either the top half of the PCA or top half of random forest importance variables for example). Anyone have any other ideas for variable selection that I haven't considered but should be? Thanks.


title: Why is square root taken for sample count "N" in standard deviation formula?

text: I'm trying to understand a very basic concept of standard deviation. From the formula $\sigma= \sqrt{ \dfrac{ \sum\limits_{i=1}^n (x_i-\mu)^2} N } $ I can't understand why should we halve the population "N" i.e why do we want to take $\sqrt{N}$ when we didnt do ${N^2}$? Doesn't that skew the population that we are considering? Shouldn't be the formula be $\sigma= \dfrac{ \sqrt{ \sum\limits_{i=1}^n (x_i-\mu)^2} } {N} $ Appreciate your help.


title: Type of data missingness in R

text: I am working on a short panel of 3 periods with a few hundred subject, and as the question and as the question suggests, I have some blanks. I know that there are not due to attrition since it is rating data of microfinance institutions and ratings have always been issued (each rating takes 3 years of accounting information, so some info related to past years is not reported). How can I use R to know whether data is MAR, MCAR or MNAR? I need this info to check whether I can apply GEE to the data. Thanks.


title: Capturing changing time trends

text: Suppose that you have ten years of monthly sales data and are interested in forecasting future sales. Consider the simple model $$\begin{equation*} y_t = \alpha_m + \beta a_t + y_{t-1} + \epsilon_t, \end{equation*}$$ where $y_t$ is log monthly sales, $\alpha_m$ is a monthly fixed effect and $a_t$ is the year at time $t$ (2011, say). The time trend $\beta$ is not constant, however; it is getting smaller over time and heading toward 0. **How can I model this?** I tried adding time-squared and log-time terms, but these send forecasts crashing out-of-sample (that is, the growth rates become negative, rather than limiting to 0). I could create discrete blocks of time and have time period-year interactions, but   * This seems like an ad hoc approach and   * There aren't obvious breaks in the time trend, but rather a more-or-less linear decline that then slows and asymptotes at 0. I could model the year-on-year growth rates instead of log sales. But the year-on-year difference is strongly correlated with a 12 period lag of itself (because log sales is not a unit root process). I would have to throw away two years of data to estimate this model with the necessary lags, which I am disinclined to do. Are there any other suggestions that you could offer? Thanks!


title: Is it necessary for a distance measure used in clustering to correspond to some valid vector space?

text: I have defined an distance measure based on some properties of points. But I'm not even sure that it corresponds to a valid distance in some vector space. Is this a necessary condition for clustering ? If yes, how do I check that it's a valid distance in some vector space. Just like Mercer's theorem can verify that the kernel is an valid cross product in some valid vector space, are there any such tests for distance ?


title: Derivation of IV estimator?

text: Assume, the model we are trying to estimate is: $Y=\beta_ß+\beta_1X+U$ where x and u are correlated: $Cov(X,U)\neq 0$ Then OLS is inconsistent: $\beta_1^{OLS}=\frac{Cov(Y,X)}{Var(X)}=\frac{Cov(\beta_0+\beta_1X+U,X)}{Var(X)}=\beta_1\frac{Cov(X,X)}{Var(X)}+\frac{Cov(X,U)}{Var(X)}=\beta_1+\frac{Cov(X,U)}{VaR(X)}\neq 0$ If we introduce an instrumental variable, for which it holds: $Cov(Z,X)\neq 0$ $Cov(Z,U) = 0$ One can show, that OLS is consitent, so that $\beta_1^{IV}=\beta_1+\frac{Cov(Z,U)}{Cov(X,Z)}=\beta_1$ I wanted to proof this. Therefore I tried the following: $\beta_1^{OLS}=\frac{Cov(Y,Z)}{Var(Z)}=\frac{Cov(\beta_0+\beta_1X+U,Z)}{Var(Z)}=\beta_1\frac{Cov(Z,X)}{Var(Z)}+\frac{Cov(Z,U)}{Var(Z)}=?$ This does not lead to a nice solution (first of all $\frac{Cov(Z,X)}{Var(Z)}\neq 1$ and I have $Var(Z)$ in the denominator and not $Cov(X,Z)$)? Why, where is my mistake and how can I proof this?


title: Algorithm for price optimization

text: I'm trying to figure out a way for calculating price optimization in a commerce environment. In other words, I'm trying to analyze how a company can increase revenue and profitability by analyzing costs, customer buying behaviors, competitive activity, demand signals and historic data. How can I come up with a strong algorithm to do this? Performance is not an issue since I'm looking to do analysis in a big data environment.


title: How to show independence and interpret the result

text: > If $X$ is symmetrically distributed about zero , then show that $U=|X|$ and > $$V = \left\{ \begin{array}{} +1, & X \geq0 \\ -1, & X<0 \end{array} \right > .$$ are independently distributed and interpret the result. Here distribution is not given only it is known that it symmetrically distributed . I think this result holds for continuous and discrete distribution . How can I show independence and interpret the result . Please help **Edit:** Here independence means statistical independence. > **Defn1:** Two random variables $X$ and $Y$ are independent if and only if > $$ F_{XY}(x,y)=F_X(x).F_Y(y) \forall x,y \in R$$ where $ F_{XY}(x,y)$ is > their joint distribution function and $F_X(x)$ and $F_Y(y)$ are their > marginal distribution functions. I also know the independence with PDF/PMF. > **Defn2:** Two random variables $X$ and $Y$ , forming an absolutely > continuous random vector,are independent if and only if > $$f_{XY}(x,y)=f_X(x).f_Y(y) \forall x,y \in R$$ where $ f_{XY}(x,y)$ is > their joint probability density function and $f_X(x)$ and $f_Y(y)$ are their > marginal probability density functions. For discrete case $f$ is replaced by $p$ .


title: R plot time indexed densities

text: This is pretty close from what I would like to do. ![enter image description here](http://i.stack.imgur.com/vr9ET.jpg) More precisely, I would like to represent graphically a sequence of inferences on time-indexed RVs $X_{t_i}$ by plotting at each time $t_i$ the density for $X_{t_i}$ (the horizontal axes representing the values taken by the $(X_{t_i})_i$ and all variables being homogeneous) at the vertical level $t_i$ (the vertical axis representing the time). I note that the inference relies on MCMC methods and thus my data are set of samples. I spent many times to find a solution on the web but did not manage to find it. If somebody knows how to do that and could help me, I will really enjoy it.


title: Delta method in a multinomial logistic regression

text: I am estimating a multinomial logistic model with R (package "mlogit"). I use the estimated coefficients to get the estimated odds and then I apply the Delta Method (package "deltamethod") to get robust standard errors. So far, I am following the instructions here. Then, I divide the estimated coefficients by the new standard errors to get new t-values and apply `2*(1-pnorm(abs(tvalues), mean = 0, sd = 1))` to get p-values. The problem is that after applying the Delta Method, some variables that have been highly significant before the transformation become insignificant. What is the cause? Did I make a mistake somewhere? Thank you very much for your help!


title: Average of root mean square error

text: Is taking the average of different rmse valid? for example average rmse = (rmse1+rmse2+rmse3)/3 Thank you for your help!


title: Kullback-Leibler divergence Pareto Distribution

text: What is the Kullback-Leibler divergence for a Pareto Distribution? Given $p(x)$ = $ \alpha$ $\frac{x^{\alpha}_{min,1}}{x^{a+1}}$.


title: Modelling risk aversion

text: How would you model risk aversion in an automated planner? This section in Planning Algorithms, briefly talks about this, but doesn't really describe a solution. However, example 4 is an interesting scenario:               Suppose that you have performed the task and are about to win the prize.     Just to add to the drama, the host offers you a gambling opportunity. You     can select action u1 and receive $10,000, or be a gambler by selecting u2      and have probability 1/2 of winning $25,000 by the tossing of a fair coin.     In terms of the expected reward, the clear choice is u2 . However, you just     completed the unpleasant task and expect to earn money. The risk of losing     it all may be intolerable. Different people will have different preferences     in this situation.      So in this case, the agent has choice u1=10000*1.0=10000, or u2=25000*0.5=12500 Now supposed the agent accumulated $500 in expenses in order to arrive at this scenario, and needs to at least break even. How would you model two different agents, one that's extremely risk-averse and so would always choose u1, since it gaurantees it'll cover its costs, and another agent that's more reckless and so would always choose u2, even though there's a greater chance of it getting nothing? This is some very crude Python to model out the classical utility-based approach, and one that I think answers my question:               #!/usr/bin/env python          expenses = 500     u1 = (10000, 1.0)     u2 = (25000, 0.5)     actions = [u1, u2]          def decide_utility(actions):         choices = [(r*p, i) for i,(r,p) in enumerate(actions)]         best = max(choices)         return best[1]+1          def decide_risk(actions, expenses, risk_factor):         choices = [((p if r>=expenses else 0.0)*risk_factor, i) for i,(r,p) in enumerate(actions)]         best = max(choices)         return best[1]+1          print decide_utility(actions) # chooses u2     print decide_risk(actions, expenses, risk_factor=1.0) # risk-averse, chooses u1     print decide_risk(actions, expenses, risk_factor=0.0) # risk-prone, chooses u2      Does this make sense? Is there a better way to do this?


title: Problems estimating covariance matrices with small $n$, smaller $p$?

text: It is well known that estimating large covariance matrices from small samples is problematic. For instance, the $p \times p$ sample covariance matrix $\Sigma_n$, estimated from $n$ samples, is not necessarily a consistent estimator of the true covariance matrix $\Sigma$ when $p > n$. In contrast, I am interested in the case where $p < n$, but with the caveat that $n$ is still relatively small; for instance, $n = 30$ and $p = 3$. Is there a reference that details estimation errors that may arise using the sample covariance matrix in this case? In particular, if the data are heavy- tailed? Are there pathological cases? Google brings up plenty of results for the $p > n$ or $p \approx n$ cases, but not the one I'm interested in.


title: Bayesian statistics for the working scientist

text: What is a good resource that explains Bayesian statistics for an experimental scientists who already use Frequentist statistics? What are the Bayesian analogues for all the tools working scientists use, like ANOVA and linear regression (and how do the Bayesian analogues compare), and what software can be used to perform the Bayesian methods? I've heard people arguing that we should publish Bayesian "likelihood ratios" instead of doing significance testing. How would one do this?


title: If a sequence of distributions converges to a degenerate, does that imply the variance strictly decreases?

text: If $F_i = G(F_{i-1}), F_0 = x$ is a sequence of distributions that converges to a degenerate distribution as $i \to \infty$, does that imply that the variance of $F_i$ decreases with $i$? Specifically, I am interested in the inverse Kumaraswamy distribution: $G(x) = (1-(1-x)^a)^b$ when $a=b$... Note that $G$ produces the cdf of the Max of b draws from the distribution of the Min of a draws from its argument (which is going to be a cdf). So if it's not true for the general case, but you have insights on how/why it would be true in this case, I would appreciate it. I can intuit (and graph for specific values of a) that the mean and variance decrease in this case -- taking the Min of $a$ draws from $F$ decreases the mean, then taking the max of $a=b$ draws from the new distribution increases the mean, but not by as much as the first operation lowered it, so the final mean is less than that of the original. But I am shaky on how to prove it (and my other questions today have basically been trying to get at this point). PS -- sorry to overload with questions today! As I mentioned, they've all been in the goal of this question here.


title: How to estimate measurement error on  spatially autocorrelated gridded data when aggregating grid cells

text: I would like to estimate the measurement error when aggregating (via arithmetic mean) gridded spatial data. The goal is to come up with the mean elevation (or some other spatially continuous variable) +/- some measure of uncertainty for each aggregated region. However, spatial data is more often than not spatially autocorrelated which may influence the measurement error. This leads to the following questions: 1/ How can I assess the magnitude of the influence of spatial autocorrelation on the standard deviation of the aggregated data? 2/ Are there other statistics that can or should be used to estimate measurement error? I offer the following example as a starting point of discussion. For simplicity I’m presenting a one dimensional dataset.               x <-runif(100,0,100)     xa <- x     # Simulate 1-D spatial autocorrelation     for (i in 4:100) {       x[i] <-  mean(xa[(i-4):i])     }      The black text in the figure is the arithmetic mean and the green text is the standard deviation. ![enter image description here](http://i.stack.imgur.com/e7Dlx.jpg)


title: How do I adjust standard errors in a research study in which the control group is constructed via matching with replacement?

text: I have a treatment sample of 200 firms. I'm using propensity score matching to pair each treatment observation with one control (sampling with replacement, in order to minimize bias associated with poor matches). This results in a sample of only 102 unique control firms (many of them are matched to multiple treatment observations). My final sample consists of 400 observations, approximately 100 of which are duplicates. My question is, what's the proper way to adjust my standard errors in a regression framework to account for the (complete) lack of independence in these duplicate observations?


title: Asymptotics of the survival function for Anderson Darling distribution?

text: I am using the `ADinf` procedure of Marsaglia & Marsaglia to compute the CDF of the Anderson Darling statistic. I am interested in the survival function, 1 minus the CDF, for large values of the test statistic. The implementation given by Marsaglia & Marsaglia does well up until roundoff issues creep in around $z=30$: ![survival function versus A2 stat](http://i.stack.imgur.com/7vpYD.png) This is somewhat to be expected. I was wondering, however, if the asymptotics of this function are known. If so, I could just approximate the tail for $z \ge 30$. Alternatively, if there is an implementation of the CDF that gives the upper tail, I would be happy with that.


title: local ordinal regressions

text: I have found no tool to estimate an "oprobit"-like ordinal model using the local regression (local likelihood) framework. It's not built into R's locfit. What is a good way to proceed? Example: I have a predictor $x\in [-1,1]$, and outcome $y\in\{1,2,3,4,5\}$. I would like to non parametrically predict (trace out, smooth) $Pr(y=n \mid x)$ for all $n \in \{1,2,3,4,5\}$ and x from $(x,y)$ observations. A latent index (with latent cutoffs) could drive the result, as there is a monotonic ranking in levels of y in $1<2<3<4<5$.


title: Regression model selection when there are more variables than cases

text: I have a database with 200+ variables and less then 50 cases. I need to choose an optimal model that predicts one dependent variable. Are stepwise/lasso regressions still appropriate methods to build such model?


title: Is there a simple equivalence test version of the Kolmogorov–Smirnov test?

text: Has two one-sided tests for equivalence (TOST) been framed for the Kolmogorov–Smirnov test to test the negativist null hypothesis that two distributions _differ_ by at least some researcher-specified level? If not TOST, then some other form of equivalence test? Nick Stauner wisely points out that (I should already know ;) that there are other nonparametric TOST equivalence tests for null hypotheses for stochastic equivalence, and, with more restrictive assumptions, for median equivalence.


title: microRNA analysis statistical methodology review

text: I am doing a review of statistical methodology used for microRNA data obtained from Affymetrix platform. I have data from 178 patients and their prognosis information as well as recurrence of disease and survival outcome (dead,alive). I learnt that I can treat recurrence rate to separate patients into two groups (recurrence 1 case and recurrence 0 control). Can you please give me a review of statistical methodologies used for such analysis? The first aim is to find a prognostic model ie correlations between microRNA and clinical outcomes and then validate it on other datasets. Thanks A. Ulgen


title: Conditional inference trees vs traditional decision trees

text: Can anyone explain the primary differences between conditional inference trees (`ctree` from `party` package in R) compared to the more traditional decision tree algorithms (such as `rpart` in R)?   1. What makes CI trees different?   2. Strengths and weaknesses? Update: I have looked at the paper by Horthorn et al that Chi refers to in the comments. I was not able to follow it completely - can anyone explain how variables are selected using permutations (e.g. what is an influence function)? Thanks!


title: PC Transformation for multinormal

text: Homework problem asks to compute the principal components of a p-variate multinormal vector $X$ ~ $N(0,\Sigma)$ where $\Sigma$ = $\delta^{2}I_{p} + \sigma ^{2} I_{p}I_{p}^{T}$. I use the PC transformation from our multivariate text where the rows of Y are the PCs and $\Gamma ^{T}$ is matrix of eigenvectors for the decomposition of $\Sigma$: Y = $\Gamma ^{T} (X - \mu)$. I'm stuck in trying to compute the spectral decomposition for $\Sigma$. I know that to compute eigenvalues, it's necessary to take the determinant of $\Sigma$ - $\lambda I_{p}$, where $\Sigma$ is expanded from above. How do I do this?


title: Looking for a name for a "mean-influencing" statistic

text: I have been using a statistic to measure "influence" or some analogous concept in the following way: First calculate the mean for a given sample, then calculate the mean for the given sample _excluding_ each observation in turn. Subtract the latter from the former. The result is a vector measuring the amount by which the sample mean is increased by the inclusion of each observation. Something like this (more or less): $$ \frac{\sum_{i = 1}^{n}}{n} - \frac{\sum_{i \neq 1}^{n}}{n - 1} $$ To me, this recalls the concept of "leverage" in regression, but I am wondering whether   1. This specific measure has been used elsewhere -- especially in the statistical literature.   2. This specific measure has a name Thanks in advance for your time.


title: Log odds as prior

text: I have this problem: In betting situations one is often interested in odds, referring in the thumbtack tossing $\theta / (1 - \theta)$. Alternatively one may consider the log-odds: $\lambda=log(\theta/(1-\theta))$ Show that a uniform distribution for $\lambda$ implies the following distribution for $\theta$: $p(\theta)=\theta^{-1}(1-\theta)^{-1}$ What problems are associated with this distribution if it is used as a prior in the thumbtack tossing problem? Now I have the thumbtack thingy, I know what it will happen if we use that p as prior, but the uniform for $\lambda$ implies the p($\theta$) I don't get it, I have tried to get to the p($\theta$) from both sides but I always get to $log(\theta)-log(1-\theta)$ and can't get out, any help??


title: ordinal regression and calcuating effect size when predictor is categorial

text: I've conducted an ordinal regression with two categorical variable (attitude toward fraud or probability of being caught=punishment). attitude range is 1:no wrong to 4:very wrong. probability of being caught range is 1:unlikely to 4: very likely. the research question is whether attitude toward fraud or probability of being caught is a stronger predictor of occurrence of fraud? The results are as below: How can I obtain an effect size or something like beta in regression when IV's are categorical? do you suggest another test? what if I wanted enter the age and gender as control? Model Fitting Information   Model -2 Log Likelihood Chi-Square df Sig. Intercept Only 1558,485   Final 359,231 1199,254 6 ,000 Link function: Logit. Parameter Estimates   Estimate Std. Error Wald df Sig. 95% Confidence Interval Lower Bound Upper Bound Threshold               [occurence = ,00]   5,364   ,119    2046,148    1   ,000    5,132   5,596     [occurence = 1,00]  6,704   ,129    2713,403    1   ,000    6,452   6,957     [occurence = 2,00]  7,942   ,157    2570,944    1   ,000    7,635   8,249     [occurence = 3,00]  9,175   ,227    1627,175    1   ,000    8,729   9,621      Location               [attitude=1,00] 2,563   ,118    474,966 1   ,000    2,333   2,794     [attitude=2,00] 2,155   ,093    537,802 1   ,000    1,973   2,337     [attitude=3,00] ,794    ,094    71,764  1   ,000    ,610    ,977     [attitude=4,00] 0a  .   .   0   .   .   .     [punish=1,00]   ,444    ,137    10,430  1   ,001    ,174    ,713     [punish=2,00]   1,182   ,118    99,799  1   ,000    ,950    1,414     [punish=3,00]   ,514    ,121    18,118  1   ,000    ,277    ,751     [punish=4,00]   0a  .   .   0   .   .   .      Link function: Logit.   a This parameter is set to zero because it is redundant.


title: Multiple comparisons for correlation matrix?

text: Does the computation of an N x N correlation matrix for N unrelated variables require multiple comparisons correction for all the computed pairwise correlations (assuming each computed correlation is a 'comparison' in the sense of being a 'statistical test')? If so, would FWER be a more appropriate measure than Bonferroni, or do any assumptions about the variables have to be made before using FWER? For instance, does it make a difference in what correction should be used whether the variables are a-priori likely to be correlated (e.g. test scores from the same group of subjects) vs uncorrelated (e.g. each variable represents the test score of a different sample)?


title: Is there an R package for nonlinear mixed effects model with spatial autocorrelation?

text: I want to fit a nonlinear mixed effects model for repeated measures data. The subjects on which the measures were repeated are spatially structured. Is there a package in R for this? I also have access to Stata, SAS, Mplus, and Mathematica. Alternatively, is it acceptable to just include random effects for membership in a neighborhood of radius R, and maybe row and column effects (the spatial structure is 2D). Thanks.


title: Logistic model, what is more important: Anova Chi-sq test or significance of coefficients

text: I have a logistic model with 8 variables. I ran a chi square test in R (anova(glm.model,test='Chisq')) and 2 of the variables turn out to be predictive when ordered at the top of the test and not so much when ordered at the bottom. the summary(glm.model) suggests that their coefficients are insignificant (high p-value). In this case it seems that the variables are not significant. I wanted to ask which is a better test of variables significance - the coefficient significance in the model summary or the chi square test from anova(). Also - when is either one better over the other? I guess it's a broad question but any pointers on what to consider will be appreciate. Thanks


title: Maximum Likelihood Estimation

text: Can someone please explain Maximum Likelihood Estimation in very lay-man terms? I understand that the Likelihood function is obtained by multiplying all the probabilities in the distribution. Or if I am wrong, please correct.


title: Analysing Several Existing Data Sets

text: I have several data sets that gathered similar data on a similar population group in one setting. Can I simply combine the data on common variables or do I need to engage in some pooled analysis?


title: Daylight saving time in time series modelling (e.g. load data)

text: I would like to model intraday electricity load data. I have quarter-hourly data for each day of the year. In many countries daylight saving time is practive. This has the consequence that once a year a day is 1 hour longer (it has 25 hours) and once a year a day is shorter (only 23 hours). I don't think that changing the time (e.g. to UTC) is a solution here because people go to work at 8:00 a.m. in their "local" time not in UTC. What I have found in the literature is e.g. here MODELING AND FORECASTING SHORT-TERMELECTRICITY LOAD USING REGRESSION ANALYSIS. There the extra hour is discarded and the missing hours is filled with average values. What is the most used procedure here? How do you cope with this probem?


title: Testing two trendlines for statistical significance

text: **Background** I am comparing data on how ATPase activity from flight and leg muscle from locusts reacts to different concentrations of calcium. I plotted a scatter graph of each data set on the same graph (ATPase activity on the Y axis, calcium concentration on the X axis) and have a trend line for each. The trend line for leg muscle has a steeper gradient, suggesting it is more responsive to changes in calcium. I carried out 9 repeats for each muscle type at each concentration of calcium. **Question** How do I test if there is a statistically significant difference between the trends of each data set? I am a novice user of R and excel, though any general pointers would be appreciated.


title: A model with binomial data?

text: I have a set of data to analyze and I would welcome some ideas. I have a series of indepedent observations when a call was produced by individuals. Each observation includes the ID of the caller, and information about other individuals in the vicinity (is there an ally? is there a friend?). There are about 400 of such observations, for about 25 individuals (that is, each of them contributed a certain amount of observations, and not the same amount). ID Called? Dominants? Friends? x yes yes no x no no yes y yes yes yes and so on... I am trying to relate call production to the presence of either dominant or ally presence (or both if they interact), controlling for the ID of the caller. Given that the dependent variable (call) and the two fixed factors (dominant, friends) have a binomial structure, what could I possibly run to show that one factor in particular predict the presence of a call? Thank you by advance


title: Regarding redundant training data in building SVM-based classifier

text: To build a SVM-based classifier, I have a training data set consisting of N data points. Some of them are redundant. For instance, there have 50 data points which are exactly the same, and there have other 100 data points which are exactly the same. I have two choices, remove the redundant ones and construct the reduced data set; keep the original data set. Will the resulting classifier be different after applying these two different choices?


title: Confidence interval for difference of means in regression

text: Suppose I have a quadratic regression model $$ Y = \beta_0 + \beta_1 X + \beta_2 X^2 + \epsilon $$ with the errors $\epsilon$ satisfying the usual assumptions (independent, normal, independent of the $X$ values). Let $b_0, b_1, b_2$ be the least squares estimates. I have two new $X$ values $x_1$ and $x_2$, and I'm interested in getting a confidence interval for $v = E(Y|X = x_2) - E(Y|X=x_1) = \beta_1 (x_2 - x_1) + \beta_2 (x_2^2 - x_1^2)$. The point estimate is $\hat{v} = b_1 (x_2 - x_1) + b_2 (x_2^2 - x_1^2)$, and (correct me if I'm wrong) I can estimate the variance by $$\hat{s}^2 = (x_2 - x_1)^2 \text{Var}(b_1) + (x_2^2 - x_1^2)^2 \text{Var}(b_2) + 2 (x_2 - x_1)(x^2 - x_1^2)\text{Cov}(b_1, b_2)$$ using the variance and covariance estimates of the coefficients provided by the software. I could use a normal approximation and take $\hat{v} \pm 1.96 \hat{s}$ as a 95% confidence interval for $v$, or I could use a bootstrap confidence interval, but is there a way to work out the exact distribution and use that? Thank you.


title: Creating auto-correlated random values in R

text: We are trying to create auto-correlated random values which will be used as timeseries. We have no existing data we refer to and just want to create the vector from scratch. On the one hand we need of course a random process with distribution and its SD. On the other hand the autocorrelation influencing the random process has to be described. The values of the vector are autocorrelated with decreasing strength over several timelags. e.g. lag1 has 0.5, lag2 0.3, lag1 0.1 etc. So in the end the vector should look something that: 2, 4, 7, 11, 10 , 8 , 5, 4, 2, -1, 2, 5, 9, 12, 13, 10, 8, 4, 3, 1, -2, -5 and so on.


title: Bandwidth selection for smooth reliability diagram

text: Following up on "How to evaluate quality of probability estimator for Bernoulli experiments?", I want to visualize the quality of an estimator for probability forecasting using a Reliability Diagram. I do not however want to put the observations into "bins", which seems to me an unprecise estimation, but rather estimate the reliability using a smoothing process similar to Kernel Smoothing. So, given $n$ observations $o_i \in \{0,1\}$ and $n$ predictions $\hat{p}_i \in [0,1]$, I want to estimate $f(x)=\frac{\sum_{i}o_i \cdot 1_{\{x\}}(\hat{p}_{i}(x))}{\sum_{i} 1_{\{x\}}(\hat{p}_{i}(x))}, x \in [0,1]$, that is the number of successes divided by the number of predictions for a probability value $x$. For a perfect predictor, $f(x)=x$. What I am doing right now is trying to estimate $f$ by calculating $\hat{f}(x)=\frac{\sum_{i} K(\frac{x-x_i}{h}) \cdot o_i}{\sum_{i} K(\frac{x-x_i}{h})}$, where for $K$ I currently use the Epanechnikov Kernel. The remaining question is: how do I properly set the bandwidth parameter $h$? I am aware that there are methods to select an optimal bandwidth for Kernel Density Estimation, but since this is not KDE, I am not sure on how to proceed here.


title: Modelling the effect of a 2 by 4 mixed design on a three-level nominal dependent variable

text: A colleague just asked me this question: ### Context: A psychological study had   * 2 `group`s of participants (between subjects)   * 4 `context`s (within subjects))   * each participant provided a `response` in each of the four contexts   * there were three categorically distinct response options (lets call them $A$, $B$, and $C$) ### Question   * What is an appropriate statistical model for analysing the effect of `group` and `context` on `response`?   * What would you tell someone who analysed this data as a 2 by 4 by 3 log-linear model? (e.g., problems with assuming independence of observations)


title: The boundary proprerty of conjugate exponential family

text: Consider a conjugate prior in the k-parameter canonical exponential family: $$p(w|n_0,y_0)=c(n_0,y_0) \exp(n_0(y_0 w'-b(w)))$$ where $n_0>0, y_0 \in Y$ is the pseudo observation vector, $Y$ is an open subset of $\mathbb{R}^k$. I have to prove that the normalization constant $c(n_0,y_0)$ goes to zero when $y_0$ tends to the boundary of $Y$. Any suggestion? I would like to find a proof that works both when $Y$ is bounded and unbounded (e.g. $Y=\mathbb{R}^k$).


title: Calculating mean age from grouped census data

text: I hope you can help me with a question regarding calculating mean age from grouped census data. If the age categories used were [0–4], [5–9], [10–14], and [15–19] years, how would you calculate the midpoints? I initially assumed the midpoints would be 2, 7, and so on. However, I read in a worked example that the midpoint should be 2.5 when age range is 0 to 4. I am assuming this has something to do with the babies not actually being zero years, but I am not exactly sure why the midpoint would be 2.5. Can anyone assist? Many thanks


title: Categorization/Segmentation techniques

text: First, let me say that I am a bit out of my depth here, so if this question needs to be re-phrased or closed as a duplicate, please let me know. It may simply be that I don't have the proper vocabulary to express my question. I am working on an image processing task in which I identify features in an image, and then classify them based on their properties, including shape, size, darkness, etc. I'm quite experienced with the image processing portion of this, but think I could improve the methods I use for classification of the features. Right now, I set thresholds for each of the parameters measured, and then classify features according to some simple logic based on which thresholds the feature has crossed. For example (the actual properties and groupings are more complex, but I'm trying to simplify irrelevant portions of my project for this question), lets say I'm grouping features into the groups "Big and Dark," "Big and Light" and "Small". Then a feature $A$ will be in "Big and Dark" iff (size($A$)>sizeThreshold) & (darkness($A$)>darknessThreshold). The goal is for the classification to agree with the classification done by an expert-level human, so I can set the thresholds to produce the best match between the groupings made by human and computer on some test set, and then hope that the classification works well with new data. This is already working pretty well, but I see one particular failure mode which I think may be fixable. Let's say feature $A$ is known to belong to "Big and Dark." The human classified it in this way because, while is was just barely big enough, it was very very dark, which made up somewhat for the lack of "bigness." My algorithm would fail to classify this feature properly, because classification is based on rigid binary logic, and requires all thresholds to be crossed. I would like to improve this failure by making my algorithm better mimic the human guided process, in which a deficiency in one parameter can be compensated by an abundance of another. To do this, I would like to take each of the base properties of my features, and convert them into some sort of score which would be a predictor of the group to which the feature belongs. I have thought of many ways of doing this, but they are mostly ad hoc ideas, based on my background in vector calculus and physics. For example, I've considered treating each feature as a vector in the N-D space of feature properties, and calculating the projection of each feature along certain vectors, each of which would measure the degree to which a feature belongs in the group. I am sure there is a more rigorous and better established technique for doing this sort of thing, but my background is relatively weak in statistical analysis, so I'm looking for a shove in the right direction. Even the name of a technique, or a link to a textbook would be helpful. **TL;DR:** What techniques are useful in classifying objects based on a large number of descriptive parameters?


title: What is the best way to measure goodness of fit between data and functions?

text: I am not that familiar with statistics and was wondering about the best way to quantify the difference between a set of data and a function. My scenario is that I have a set with maximum heart rate data from around 60 people, with their age. I'm trying to see which of all "mathematical models" out there that describe the relation between maximum heart rate and age, is the least wrong. I'm much better at understanding the biology behind this, not the mathematics. Thank you for any assistance and pointers!


title: Interpretation of an equation

text: I'm trying to describe the equation in the attached screenshot in my own words. $\lvert(\overline\chi_1-\overline\chi_2)\rvert > s_1(t_{1 @ 99 \%}) + s_2(t_{2 @ 75 \%})$ The equation tests whether 99% of population 1 lies outside of 75% of population 2. The equation is taken from Donegan, T., 2012. Geographical variation in Immaculate Antbird Myrmeciza immaculata, with a new subspecies from the Central Andes of Colombia. Bull. Brit. Orn. Club, 132, pp.1–40. My current interpretation is as follows: The difference between the means of the population 1 and 2 must be greater than 99% of a standard deviation in population 1 plus 75% of a standard deviation in population 2. I realise my interpretation may not be quite right. Can anybody give a better interpretation?


title: Measuring what's 'lost' in PCA dimensionality reduction?

text: I'm working with multiple time series signals ${\{X_i\}}$ and one method to remove suspected noise is with PCA. But unlike most methods which remove the components with least variance, we remove much more substantial components (often the second or third). This is often done because empirically it's been shown to improve results in certain analyses (like simple pairwise correlation) but has been detrimental to other analysis like clustering. So I was curious of general metrics to compare before and after dimensionality reduction. Would looking at things like entropy be useful?


title: Interpretation of interaction and main effect in a mixed model when one is significant and the other is not

text: I have repeated (x4) measurements on 90 subjects. The outcome is zero inflated. The model output gives these estimates:               Coefficients:                              Estimate Std. Error z value   Pr(>|z|)       (Intercept)               -2.2151     1.6998   -1.30    0.193      time                       0.5283     0.2167    2.44    0.015     X                          0.0791     0.3267    0.24    0.811     time:X                    -0.1525     0.0611   -2.49    0.013     Y                          1.0348     0.4488    2.31    0.021     time:Y                    -0.0037     0.0583   -0.06    0.950          Number of observations: total=342, Subject=90     Random effect variance(s):     Group=Subject                 Variance StdDev     (Intercept)    6.302   2.51     Negative binomial dispersion parameter: 3.9265 (std. err.: 0.69723)     Zero-inflation: 0.13169  (std. err.:  0.027756 )      `X` and `Y` both vary at the subject-level only, not at the time-level. As you see, here the main effect of `X` is not significant, while the `time:X` interaction is significant. On the other hand, the main effect of `Y` is significant, while the `time:Y` interaction is not. I believe that in the case of Y, this means that subjects with higher values of `Y` have higher trajectories of the outcome. Is this correct ? However, how is the `time:X` interaction to be interpreted ?


title: Interpretation of McLeod-Li test for ARCH effect

text: I want to know if I need to use ARCH/GARCH model for my time series. How can I use `McLeod.Li.test` in R to do that and how should I interpret the result? the description in R help was not clear for me because in the example of the R help they used difference of the log of the time series and I do not know why? Also, how should I interpret the resulted plot? R help example is like this:               data(CREF)     r.cref=diff(log(CREF))*100     McLeod.Li.test(y=r.cref)      Thanks


title: Sample weighting

text: **How do I work out an adequate/representative sample in the following scenario?** My intention is to use correlation and regression models to test the relationship between income level and money spent on a particular product by store type. **Scenario** There are 200 stores, with a total of 10,000 registered customers. The stores can be grouped in terms of their floor space as large (100 stores), medium (70 stores) and small (30 stores). The large stores account for 60% of the registered customers, the medium stores account for 30% of the customers and the small stores account for 10% of the customers. **Sampling Method** Using the cluster sampling technique, I have randomly chosen 2 large, 2 medium and 2 small stores. (This is on the assumption that the stores are similiar as per their group size.) Then, using stratified random sampling, I chose 50 customers from each of the six stores such that I have 100 large store customers, 100 medium store customers and 100 small store customers. (This was to ensure gender and age balance.) My final sample is then 300 customers. **Is this a representative sample (given the scenario)?** Or, do I have to use some sort of weighting to ensure the final sample reflects:   * the store size i.e. 100 are large stores, 70 are medium stores and 30 are small stores.    * customer distribution i.e. 60% of the customers are from large stores, 30% are from medium stores and 10% are from small stores) Here a method I propose to use. It is **Probability Proportionate to Size (PPS) Sampling Method.** Please see my comment below.


title: How can I perform weibull analysis on monthly recorded data of wind speeds?

text: I've read several articles about how to perform Weibull distribution but they all did it with wind speed data in a time series manner (example: recorded data every 10mins and then averaged to every hour). My prob is that I've got only mean monthly wind speed for a year (from January to December). So how should I deal with this situation to perform Weibull distribution on these data? I have data for 5 years namely 2006, 2007, 2008, 2009, 2010 and like i said, its monthly wind speed from January to December. I'd appreciate if you could help out. Thank you in advacnce...


title: Canonical Correlation Definitions

text: Suppose that we have two sets of variables $\textbf{X}$ and $\textbf{Y}$. Let $$ \textbf{X} = \begin{pmatrix} X_1 \\ X_2\\ \vdots \\X_p \end{pmatrix}$$ and $$ \textbf{Y} = \begin{pmatrix} Y_1 \\ Y_2\\ \vdots \\Y_q \end{pmatrix}$$ where $p \leq q$. Define the following linear combinations: $$U_i = a_{i1}X_1 + a_{i2}X_2 + \cdots + a_{ip}X_p$$ $$V_i = b_{i1}Y_1 + b_{i2}Y_2 + \cdots + b_{iq}Y_q$$ Does it matter what we call the first canonical variate pair? For example is the first canonical variate pair $(U_1, V_1)$? Or would we define the first canonical variate pair as $(U_1, V_3)$?


title: Correlation between growths variables

text: I'm performing correlations between growth variables (e.g. debt growth and net income growth) and the correlation is rather strong. Nevertheless, I was wondering why when I perform the same correlation for debt growth and net income (e.g. net income 2008), so not comparing growth and growth I don't get any correlation. Finally I also tried correlation between debt and net income for my sample and didn't get any correlation neither. Could somebody explain me why there are correlations between growths variables and not between the variables themselves? does it make sense to perform correlations between growth variables? Thank you very much for your help!


title: Reference book for ecological studies

text: I'm analyzing a data base from an ecological study and I'd like to read a good reference book about this kind of studies as there are some limitations concerning their analysis; for example, they are susceptible to ecological fallacy.


title: Bayesian estimation of tridiagonal covariance

text: I want to estimate covariance of a multivariate normal distribution from data using a Bayesian method. I want to force the result to be tridiagonal. I am looking for an appropriate prior or method. Is there any prior or method for this?


title: Summing up ordered items to a sum-score: implications?

text: It is common practice in psychology to add together ordered items to sum- scores, and then use these sum-scores for subsequent statistics. This is usually done without explicitly mentioning measurement level, and I guess people simply assume items to be interval or ratio. I wonder about this practice, and the implications. An example are screening instruments for mental disorders. Items are mostly ordered (and not interval or ratio), and a sum-score is used as threshold to determine whether someone has a disease or not. Isn't it correct that for a sum to be meaningful, equidistant points between each of the scale elements would be needed, allowing for an interpretation of differences in the _degree_ of order? Adding up the numbers 2, 2, and 3 to a sum-score only makes sense if the first 2 is equal to the second 2, and if the 3 is 1.5 times larger than each of the 2s (which is highly problematic with ordered items, obviously). Is there (critical) statistical literature regarding that topic? What would be your opinion? I have looked for a while now but probably wouldn't even know the technical terms to look properly.


title: Event Prediction through Machine Learning

text: I have a large data set consisting of ca. 40 categorical data items and a few interval data items (real numbers, less than 5 such items). Most categories should have a lot of values that repeat themselves over and over and very few that appear very rarely. Some categories are also overcategories of others (like country and city). The outcome of each data is either 1 if the event occured or 0 if it did not occur. The idea is to calibrate a machine learning model or a staistical model that can predict for every given data row the probability that the outcome is 1. The data set I will use will have at least 1 million rows. What machine learning approaches and statistical models will perform well on such a task? My initial thoughts are logarithmic regression and support vector machines (with extensions like random forest) How do I deal with the interval data items? The easiest approach is obiously to convert different ranges into categories, which I think will not be a problem. What libraries and tools can I use, when my data set has a size of 10 GB? I am interested in tools/libraries that include machine learning algorithms but also statistical tools to help me find attributes with significant influence on the outcome. I can code in Java and C++ at the moment. I looked into Root, a data analysis tool from CERN for lare data sets and its machine learning module TMVA but it can only handle real numbers and integers as input as far as I know.


title: Simple linear regression: confidence interval?

text: Given data $(x_{1}, y_{2}),..., (x_{n}, y_{n})$, consider the simple linear regression model $E(Y | X=x) =3.2 + B_{1}x, Var(Y|X=x)=\sigma^{2}$. Construct a 95% confidence interval for $B_{1}$. I know that we first have to estimate the $\sigma$, but I was told that the degrees of the confidence interval is n-1, which I don't get because doesn't the simple linear regression have 2 degrees of freedom since there are 2 parameters? Also, I calculated the variance of $\hat{B_{1}}$ to be $\frac{\sigma^{2}}{\Sigma {x_{i}}^{2}}$ but I don't know why in the confidence interval, the $se(\hat{B_{1}}) = \frac{\hat{\sigma}}{\Sigma{x_{i}}^{2}}$. I would really appreciate it if someone can explain my two issues to me. Thanks in advance.


title: Generate synthetic data from a graph (e.g., bipartite)

text: Does anyone know of literature describing methods for generating synthetic data from a graphical structure that has been learned from observations of real data? The graph, depending on the type of data or the questions being asked of it, could be bipartite (as in topic modeling), or directed (as in a link graph-type application). I have spent quite some time researching this, and the closest I've been able to find on this is research on generating synthetic data from statistical distributions that have been learned from training data.


title: Random Forest and cluster-level bootstrapping

text: I'm working with cluster-correlated data (individuals nested into households). I would like to modify the bootstrap sampling process to make a household- level bootstrap instead of the subject-level bootstrap. My aim is to ensure that the entire household is included in the sample used to grow each tree of the forest. I'm working with R package 'randomForest'.


title: How can I improve a simple regression here?

text: I'm trying to solve the following issue : Let's say that I have three normal random variables a,b,c, non correlated. Let's say that I only have two observations of these, M and N where :               M = a + b     N = a + c      And another variable I'm trying to explain :               R = a + b + c + err      I can do a simple regression on this data, `lm(R ~ M + N + 0)`, I get the same loading for both `M` and `N`, equal to `0.66`. My final model is `1.33a + 0.66b + 0.66c`. I know that this model is the best linear model here, but what could I use to improve it. Are there methods to properly get `a, b, c` from `M`and `N` and do the regression on them ? I've tried PCA but I only get two vectors, not three. Basically, I'd need a method that would project `M` and `N` into the three interesting components here. **EDIT** I realise I haven't been clear here. $M$, $N$, $R$ are three vectors, with a high number of observations. If they had the relations above with three (unknown) factors $a$, $b$, $c$, then a simple regression would give me : $\hat{R} = 0.66 M + 0.66 N = 1.33 a + 0.66 b + 0.66 c$. I am wondering if there is an alternative technique for modelling $R$, that would give a closer result to the actual $R$. Something that would somehow reconstruct the common term $a$ from $M$ and $N$ automatically, and take that into account to finally get : $\hat{R} = 1 a + 1 * (M - a) + 1 * (N - a)$ I hope I'm a bit more clear in my explanations ...


title: How to run a sensitivity analysis with dependent variables?

text: I am trying to run a sensitivity analysis on an efficiency indicator: the ratio between nitrogen outputs (milk, meat, crops) and nitrogen inputs (fertilizers, cattle feed, symbiotic fixation…). I have 4 output variables and 12 input variables. However, they are not independent. For instance, fertilizer input is correlated with crop output. It seems that the methods I read about (ANOVA, Morris, Sobol) are based on the hypothesis that all variables are independent. This allows testing the effect of the variation of each variable while keeping all other variables constant. Do you have any suggestion of how I could take these correlations into account in my sensitivity analysis ?


title: How can I determine weibull parameters from data?

text: I have a histogram of wind speed data which is often represented using a weibull distribution. I would like to calculate the weibull shape and scale factors which give the best fit to the histogram. I need a numerical solution (as opposed to graphic solutions) because the goal is to determine the weibull form programmatically. **Edit:** Samples are collected every 10 minutes, the wind speed is averaged over the 10 minutes. Samples also include the maximum and minimum wind speed recorded during each interval which are ignored at present but I would like to incorporate later. Bin width is 0.5 m/s ![Histogram for 1 month of data](http://i.stack.imgur.com/C0WCk.png)


title: How to analyse reading data with two word types and three types of dependent variables?

text: In my study, I have one group encoding 2 different types of words (within subjects) and 3 different dependent measures namely recall, recognition, rating. I think that I can conduct repeated measure analysis. But I am confused about what should I write for factor name. I think I have one factor with two levels because all subjects encoded 2 different types of words in the same list. Is that right? And I defined three measure names to define DVs.Is that OK?


title: How to simulate head to head competition based on winning percentages?

text: If I have two players that are about to engage in a head to head competition, how do I simulate a winner based only on win percentage? One idea I had was to create a "weighted winning pool" where I add the two percentages and then randomly pick a number in that pool. If it's in the percentage range of a team then that team wins. e.g. 100% and 100% =200. If the random number is less than 100 then team A wins, if not then team b wins   e.g. #2 100% and 50% = 150. If the random number is less than 100 then team A wins if not team B wins. Would this work or is there a more accepted way of doing this? Thanks in advance.


title: Datasets for clustering algorithms

text: I am asked to give a lecture on clustering algorithms for an audience that is not very technical. With that in mind, I wanted to do a simple exercise where I will ask the audience to identify groups from a dataset. However, I cannot find good datasets that could be usable for this purpose. Is there a dataset of customers and some products they have bought that I can use for this purpose? Or any other dataset that might look suitable!


title: Calculating the Variance using Delta Method

text: I'm trying to find the variance of $L$, $Var(L)$, using the delta method (I want to find a closed form). $L$ is defined as: $$L = \frac{A}{B} + \frac{C}{D}$$ All $A$, $B$, $C$, and $D$ are dependent. But I'm not sure how to proceed.


title: Questions about seqtreedisplay and GraphViz

text: After I installed Graphviz, I used the following command to plot a tree resulting from a tree-structured analysis of my sequences: `> seqtreedisplay(mytree, filename=mytree, seqdata=mysequences, type="I", border=NA)` However, I had this error message: > 'dot' is not recognized as an internal or external command, operable program > or batch file. Error in disstreedisplayInternal(tree = tree, filename = > filename, tmpdisstree = tmpdisstree, : You should install GraphViz to use > this function: see http://www.graphviz.org In addition: Warning message: > running command 'C:\Windows\system32\cmd.exe /c dot -Tpng > -otmpseqtree1b6cc27794e.png tmpseqtree1b6cc27794e.dot' had status 1 I would appreciate your suggestions!


title: Are there known statistical tests to determine the parameters of a machine learning algorithm?

text: What i am trying to do is to take a machine learning algorithm that is already trained.Run the algorithm several times and collect the output. Conduct a statistical analysis on the output and try to create an inference regarding the parameters of the algorithm Is this possible ? Are there known statistical tests which does this or which does something similar to this ? Adding from the comment I was not sure whether to add statistical bias or not.I ll try to explain my thought flow when i added that tag. If in fact we can learn something about the algorithm (it could be the parameters used or something else) from the output, then it should be(i am not sure) because the output is biased in some way. My aim to find out if there is some way to learn something about the algorithm used, using the properties of this output.Have you heard of such a thing ?