title: Difference between null distribution and sampling distribution

text: I am asking this questions for clarifications on the terms between "null distribution" and "sampling distribution". I find that when some one says **null distribution** , actually they mean the same thing as when others say **sampling distribution**. In this hypothesis testing article [1], you can see the following description of an example _We are considering a random variable Y which is normally distributed. (This is one of the model assumptions.)_ _Our null hypothesis is: The population mean µ of the random variable Y is a certain value µ0. For simplicity, we will discuss a one-sided alternative hypothesis: The population mean µ of the random variable Y is greater than µ0. (i.e., µ > µ0)_ _Another model assumption says that samples are simple random samples. We have data in the form of a simple random sample of size n._ _To understand the idea behind the hypothesis test, we need to put our sample of data on hold for a while and consider all possible simple random samples of the same size n from the random variable Y._ _For any such sample, we could calculate its sample mean ȳ and its sample standard deviation s. We could then use ȳ and s to calculate the t-statistic t = (ȳ - µ0)/(s/√n)_ _Doing this for all possible simple random samples of size n from Y gives us a new random variable, Tn. Its distribution is called a **sampling distribution**._ _The mathematical theorem associated with this inference procedure (one-sided t-test for population mean) tells us that if the null hypothesis is true, then the sampling distribution has what is called the t-distribution with n degrees of freedom._ I have no problem with understanding the content, my main concern is about the term "sampling distribution". Here the so-called sampling distribution refers to the test statistic distribution, given the null hypothesis is true. It is a theoretical distribution. According to wikipedia [2], null hypothesis seems to mean the same thing. I have read many lectures notes on statistics, I find both terms coexist. But if you search for sampling distribution, you get many more results. Could someone clarify my doubts-- **Do null distribution and sampling distribution mean the same thing?** Reference: [1] http://www.ma.utexas.edu/users/mks/statmistakes/hyptest.html   [2] http://en.wikipedia.org/wiki/Null_distribution


title: Usage of LDA with more than two classes

text: I'm reading about the Linear Discriminant Analysis by Fisher and I have a couple of questions about its usage.   1. If you have k>2 classes in a two-dimensional space you find k−1 vectors that you need to use to project the sample data. Is it possible that one sample is closer to different means along different vectors?   2. Suppose you have, for example, 3 classes in a 3-dimensional space. You get two vectors, but now you need planes to be able to project your samples. Do you create those by combining the vectors obtained by LDA?


title: What's the best way to transform this vector to be normal-like?

text: Assume I have an outcome like this:               x<-c(rep(2.83,3),          rep(3.17,4),          rep(3.83,4),          rep(4.17,5),          rep(4.83,8),          rep(5.5,3),          rep(7.17,5),          rep(8.17,7),          rep(8.83,12),          rep(9.5, 12),          rep(9.83,17),          rep(10.17,50))      hist(x)     hist(log(x))      I would use it as a continuous dependent variable in a regression analysis but have no idea how to make it normal-like. I have tried few algorithms like log, quadratic, or boxcox, but it seems neither way is making better. Does anybody here have some suggestions?


title: What do I gain if I consider the outcome as ordinal instead of categorical?

text: There are different methods for prediction of ordinal and categorical variables. What I do not understand, is how this distinction matters. Is there a simple example which can make clear what goes wrong if I drop the order? Under what circumstances does it not matter? For instance, if the independent variables are all categorical/ordinal, too, would there be a difference? This related question focuses on the type of the independent variables. Here I am asking about outcome variables. **Edit:** I see the point that using the order structure reduces the number of model parameters, but I am still not really convinced. Here is an example (taken from an introduction to ordered logistic regression where as far as I can see ordinal logistic regression does not perform better than multinomial logistic regression:               library(nnet)     library(MASS)     gradapply <- read.csv(url("http://www.ats.ucla.edu/stat/r/dae/ologit.csv"), colClasses=c("factor", "factor", "factor", "numeric"))          ordered_result <- function() {       train_rows <- sample(nrow(gradapply), round(nrow(gradapply)*0.9))       train_data <- gradapply[train_rows,]       test_data <- gradapply[setdiff(1:nrow(gradapply), train_rows),]       m <- polr(apply~pared+gpa, data=train_data)       pred <- predict(m, test_data)       return(sum(pred==test_data$apply))     }          multinomial_result <- function() {       train_rows <- sample(nrow(gradapply), round(nrow(gradapply)*0.9))       train_data <- gradapply[train_rows,]       test_data <- gradapply[setdiff(1:nrow(gradapply), train_rows),]       m <- multinom(apply~pared+gpa, data=train_data)       pred <- predict(m, test_data)       return(sum(pred==test_data$apply))     }          n <- 100          polr_res <- replicate(n, ordered_result())     multinom_res <- replicate(n, multinomial_result())     boxplot(data.frame(polr=polr_res, multinom=multinom_res))      which shows the distribution of the number of right guesses (out of 40) of both algorithms. ![polr_vs_multinom](http://i.stack.imgur.com/fUKhS.png) **Edit2:** When I use as scoring method the following               return(sum(abs(as.numeric(pred)-as.numeric(test_data$apply)))      and penalize "very wrong" predictions, polr still looks bad, i.e. the plot above does not change very much.


title: Including a variable which is identical for each cross section member

text: I am doing a static panel analysis: four independent variables, seven cross sections and $T = 13$ years. I want to include a variable (exchange rate) which varies over time but is identical (has the same value) for each of the cross section members. ### Question: **Is it possible to use a variable in the regression which has an identical value sequence for all of the cross section members?**


title: Which measure of model fit to report when performing likelihood based regression: AIC, BIC, Pseudo R-square?

text: I'd like to hear your opinions on the following:   * What parameters would you report when estimating different likelihood based regression? AIC, BIC, Pseudo $R^2$?    * What is the standard to report?  It should be a parameter which answers the question of how good the specified model is.


title: How to visualize a GraphML multitree?

text: I have a multitree that represents the lineages of all the fishgroups in a breeding program. It's stored as a double adjacency table with `fish_id`, `sire_id` and `dam_id`. This means a particular fishgroup is only directly aware of its parents (if any) and knows nothing about descendents. I can conveniently output this data as a GraphML doc, but I'm having difficulty finding a viewer that can display a multitree given simple input. InfoVis can, but requires a verbose and redundant description of the structure. Cytoscape Web can take GraphML objects, but can't automatically layout a multitree (it does it's best, but I'd like a hierarchical layout). Does anyone know of another option? My GraphML looks like this:               <?xml version="1.0" encoding="UTF-8"?>     <graphml xmlns="http://graphml.graphdrawing.org/xmlns"         xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"         xsi:schemaLocation="http://graphml.graphdrawing.org/xmlns              http://graphml.graphdrawing.org/xmlns/1.0/graphml.xsd">         <key id="currentfish" for="node" attr.name="cfish" attr.type="string">             <default>#ffff00</default>         </key>         <key id="parcount" for="edge" attr.name="parcount" attr.type="int"></key>         <key id="parenttype" for="edge" attr.name="parenttype" attr.type="string"></key>         <key id="fish" for="node" attr.name="fish" attr.type='string'></key>         <graph id="G" edgedefault="directed">             <node id='f22'>                 <data key='fishname'>ZS6181</data>             </node>             <edge id='es_23' source='f57' target='f23'>                 <data key='parcount'>0</data>                 <data key='parenttype'>sire</data>             </edge>             <node id='f23'>                 <data key='fishname'>90716</data>             </node>             <edge id='ed_23' source='f42' target='f23'>                 <data key='parcount'>0</data>                 <data key='parenttype'>dam</data>             </edge>             <node id='f24'>                 <data key='fishname'>ZS6377</data>             </node>             <node id='f25'>                 <data key='fishname'>ZS6375</data>             </node>             <edge id='ed_26' source='f25' target='f26'>                 <data key='parcount'>10</data>                 <data key='parenttype'>dam</data>             </edge>             <node id='f26'>                 <data key='fishname'>ZS6375 F1</data>             </node>             <edge id='es_26' source='f25' target='f26'>                 <data key='parcount'>10</data>                 <data key='parenttype'>sire</data>             </edge>             <node id='f27'>                 <data key='fishname'>ZS56181 F1</data>             </node>             <edge id='es_27' source='f43' target='f27'>                 <data key='parcount'>9</data>                 <data key='parenttype'>sire</data>             </edge>             <edge id='ed_27' source='f43' target='f27'>                 <data key='parcount'>9</data>                 <data key='parenttype'>dam</data>             </edge>             <node id='f28'>                 <data key='fishname'>ZS6377 F1</data>             </node>             <edge id='es_28' source='f44' target='f28'>                 <data key='parcount'>7</data>                 <data key='parenttype'>sire</data>             </edge>             <edge id='ed_28' source='f44' target='f28'>                 <data key='parcount'>7</data>                 <data key='parenttype'>dam</data>             </edge>             <node id='f29'>                 <data key='fishname'>100128 AB</data>             </node>             <edge id='es_30' source='f25' target='f30'>                 <data key='parcount'>10</data>                 <data key='parenttype'>sire</data>             </edge>             <node id='f30'>                 <data key='fishname'>ZS6375 F1A</data>             </node>             <edge id='ed_30' source='f25' target='f30'>                 <data key='parcount'>10</data>                 <data key='parenttype'>dam</data>             </edge>             <edge id='ed_31' source='f45' target='f31'>                 <data key='parcount'>0</data>                 <data key='parenttype'>dam</data>             </edge>             <node id='f31'>                 <data key='fishname'>AB 100223</data>             </node>             <edge id='es_31' source='f25' target='f31'>                 <data key='parcount'>0</data>                 <data key='parenttype'>sire</data>             </edge>             <edge id='es_32' source='f45' target='f32'>                 <data key='parcount'>0</data>                 <data key='parenttype'>sire</data>             </edge>             <edge id='ed_32' source='f45' target='f32'>                 <data key='parcount'>0</data>                 <data key='parenttype'>dam</data>             </edge>             <node id='f32'>                 <data key='fishname'>AB 100319</data>             </node>             <edge id='es_33' source='f32' target='f33'>                 <data key='parcount'>8</data>                 <data key='parenttype'>sire</data>             </edge>             <edge id='ed_33' source='f31' target='f33'>                 <data key='parcount'>8</data>                 <data key='parenttype'>dam</data>             </edge>             <node id='f33'>                 <data key='fishname'>AB 100714</data>             </node>             <edge id='ed_34' source='f28' target='f34'>                 <data key='parcount'>10</data>                 <data key='parenttype'>dam</data>             </edge>             <node id='f34'>                 <data key='fishname'>AB 100715</data>             </node>             <edge id='es_34' source='f59' target='f34'>                 <data key='parcount'>10</data>                 <data key='parenttype'>sire</data>             </edge>             <edge id='ed_35' source='f27' target='f35'>                 <data key='parcount'>4</data>                 <data key='parenttype'>dam</data>             </edge>             <node id='f35'>                 <data key='fishname'>AB 100722</data>             </node>             <edge id='es_35' source='f28' target='f35'>                 <data key='parcount'>5</data>                 <data key='parenttype'>sire</data>             </edge>             <node id='f36'>                 <data key='fishname'>ZS6377 F2</data>             </node>             <edge id='es_36' source='f28' target='f36'>                 <data key='parcount'>8</data>                 <data key='parenttype'>sire</data>             </edge>             <edge id='ed_36' source='f28' target='f36'>                 <data key='parcount'>6</data>                 <data key='parenttype'>dam</data>             </edge>             <node id='f37'>                 <data key='fishname'>ZS6375 F2</data>             </node>             <node id='f38'>                 <data key='fishname'>100730 AB</data>             </node>             <edge id='es_39' source='f34' target='f39'>                 <data key='parcount'>6</data>                 <data key='parenttype'>sire</data>             </edge>             <edge id='ed_39' source='f38' target='f39'>                 <data key='parcount'>6</data>                 <data key='parenttype'>dam</data>             </edge>             <node id='f39'>                 <data key='fishname'>110208</data>             </node>             <node id='f40'>                 <data key='fishname'>110412</data>             </node>             <edge id='ed_41' source='f27' target='f41'>                 <data key='parcount'>2</data>                 <data key='parenttype'>dam</data>             </edge>             <node id='f41'>                 <data key='fishname'>110413</data>             </node>             <edge id='es_41' source='f24' target='f41'>                 <data key='parcount'>5</data>                 <data key='parenttype'>sire</data>             </edge>             <node id='f42'>                 <data key='fishname'>90318</data>             </node>             <node id='f43'>                 <data key='fishname'>ZS56181-90705</data>             </node>             <node id='f44'>                 <data key='fishname'>ZS6377-90913</data>             </node>             <node id='f45'>                 <data key='fishname'>ZS56181</data>             </node>             <node id='f57'>                 <data key='fishname'>90317</data>             </node>             <node id='f59'>                 <data key='fishname'>ZS63775 F1</data>             </node>         </graph>     </graphml>


title: Interpretation of weights in non-linear least squares regression

text: I am conducting a non-linear least squares regression fit using the python scipy.optimize.curve_fit function, and am trying to better understand the weights that go into this method. I have a distribution of raw data points that I wish to fit to a Gaussian cumulative distribution function. I created a function for this that takes three parameters: an average, a standard deviation, and a scale factor for if my distribution doesn't quite approach 1. My confidence in each raw data point is based on a separate instrumental count that doesn't necessarily have anything to do with the value of the data point, so I'm trying to include this in my fit using weights. Specifically, small and large values of x have less certainty, so I want them to matter less in the regression. When I conduct a fit by passing these raw counts as weights, the fit is not particularly good, whereas if I pass them as (1/counts) the fit improves. I have plotted the raw data, fits, and normalized weights for these two options. What I am trying to understand is how to interpret the weights. I would have thought that higher values for the weights imply more importance in the regression. However, it seems that it is actually correct to use the weights where the "bad" raw data points have a higher weight. Why is this and how should I be interpreting the weights? Also, is there a better resource for understanding weights in non-linear regression? Most of the sources I have found have not explained weighting in a way I can comprehend. Edit: I added a second plot that now shows actual (non-normalized) counts, along with the corrected fit (weighted according to counts) according the proper fitting technique shown below. ![enter image description here](http://i.stack.imgur.com/k7irq.png) ![enter image description here](http://i.stack.imgur.com/bfgbo.png)


title: SPSS: "Weight cases" by multiple factors, by taking the product of weights?

text: My understanding is that a sample can be adjusted based on a known population using the "Weight Cases" function in SPSS. What is the proper procedure for weighting by multiple factors? Let's say I have WEIGHT_SEX and WEIGHT_ETHNICITY. Can I compute a WEIGHT_COMPOSITE by taking their product? Are there any downsides to doing this that I should be aware of?


title: Why can weather prediction be so correct?

text: I've done a test using ARMA model on some financial series. It turns out the prediction rate is really very bad – close to half time correct and half time wrong… I am new to ARMA model so what I tried is very simple, following the textbook of deciding (p,q) for ARMA model by ACF and PACF first, and then use half data to do regression and half data to do test. I often heard that weather forecasting is used time series model. And in my daily experience, I feel it is pretty good forecasting. I wonder why the weather forecast could be done so great? Also, if I want to push my toy ARMA model to next level, what is the direction that I should put into effort ? Following up:   There are good answers below, so I am summarize it a little up :   1. From a systematic point of view, market as a system is evolving, while weather is more stable from one year to another (truth doesn't change).    2. Market is more complicated in some sense, while weather is ruled by physics laws (more understood structure).    3. Weather is more periodic, which adds to the predictability.


title: How to classify data having sub-instance features?

text: I am trying to use machine learning on some peculiar (at least for me) data. Usually, when I do machine learning I am use to have the data in this format:                   Feat1 Feat2 Feat3 Feat4  Class     I1 | .0  | .0  | .2  | .3  | TRUE  |     I2 | .5  | .0  | .1  | .3  | TRUE  |     I3 | .0  | .0  | .1  | .3  | FALSE |      In this case, my data is in this format:                       F1   F2   F3   Class      I1 a | .0 | .2 | .3 |         b | .0 | .2 | .3 |         c | .0 | .2 | .3 | TRUE  |           I2 a | .9 | .1 | .1 |         b | .0 | .0 | .0 |         c | .0 | .0 | .0 | TRUE  |           I3 a | .0 | .0 | .1 |         b | .0 | .2 | .3 |         c | .0 | .2 | .1 | FALSE |      So I have multiple rows concurring to the final classification of an instance, and I can compute the features only at this sub-level. One way to aggregate the features in a single row would be to, for example, using things such as 'mean', 'min/max', etc. But there is no clear choice for these aggregations in our dataset. The number of rows is not known a-priori, so I cannot create a single row aggregating the multiple rows to model the instance. For example, I cannot have:                      F1-a  F2-a  F3-a  F1-b  F2-b  F3-b  F1-c  F2-c  F3-c  Class      I1 a | .0  | .2  | .3  | .0  | .2  | .3  | .0  | .2  | .3 | TRUE  |      Moreover, this does not seem to be a multiple-instance learning problem, because I am not trying to classify the single row, but the instances (which are made of features from multiple rows). Thank you!


title: How do I get a a p-value for the output of an lme model with lme4?

text: My problem can be summarized very simply: I'm using a linear mixed-effects model and I am trying to get p-values using pvals.fnc(). The problem is that this function seems to have some trouble estimating p-values directly from the t-values associated with model coefficients (Baayen et al., 2008), and I don't know what is going wrong with the way I do it (i.e. according to what I have read, it should work). So, I'm explaining my model below, and if you can point out what I am doing wrong and suggest changes I would really appreciate it! **DESCRIPTION** : I have a 2 by 2 within subjects design, fully crossing two _categorical_ factors, "Gram" and "Number", each with two levels. This is the command I used to run the model:               >m <- lmer(RT ~ Gram*Number + (1|Subject) + (0+Gram+Number|Subject) + (1|Item),data= data)      If I understand this code, I am getting coefficients for the two fixed effects (Gram and Number) and their interaction, and I am fitting a model that has by- subject intercepts and slopes for the two fixed effects, and a by-item intercept for them. Following Barr et al. (2013), I thought that this code gets rid of the correlation parameters. I don't want estimate the correlations because I want to get the p-values using pvals.fnc (), and I read that this function doesn't work if there are correlations in the model. The command seems to work:               >m     Linear mixed model fit by REML      Formula: RT ~ Gram * Number + (1 | Subject) + (0 + Gram + Number | Subject) + (1 |Item)         Data: mverb[mverb$Region == "06v1", ]         AIC   BIC logLik deviance REMLdev      20134 20204 -10054    20138   20108     Random effects:      Groups      Name        Variance  Std.Dev. Corr                Item       (Intercept)   273.508  16.5381                     Subject     Gramgram        0.000   0.0000                                 Gramungram   3717.213  60.9689    NaN                          Number1        59.361   7.7046    NaN -1.000       Subject     (Intercept) 14075.240 118.6391                     Residual                35758.311 189.0987                    Number of obs: 1502, groups: Item, 48; Subject, 32          Fixed effects:                  Estimate Std. Error  t value     (Intercept)    402.520     22.321  18.033     Gram1          -57.788     14.545  -3.973     Number1         -4.191      9.858  -0.425     Gram1:Number1   15.693     19.527   0.804          Correlation of Fixed Effects:                 (Intr) Gram1  Numbr1     Gram1       -0.181                   Number1     -0.034  0.104            Gram1:Nmbr1  0.000 -0.002 -0.011      However, when I try to calculate the p-values I still get an error message:               >pvals.fnc(m, withMCMC=T)$fixed     Error in pvals.fnc(m, withMCMC = T) :      MCMC sampling is not implemented in recent versions of lme4       for models with random correlation parameters      Am I making a mistake when I specify my model? Shouldn't pvals.fnc() work if I removed the correlations? Thanks for your help!


title: Is bootstrapping appropriate for estimating a multivariate normal covariance matrix using a small sample size?

text: Let $p$ be a matrix, where each row represents an observation of a 4-variate normally distributed random variable $\mathcal{N}_4(\mu,\Sigma)$.   * Is there any Bootstrap methode to get a good estimation for Σ?   * If not, is the number of the following sample enough to bootstrap the distribution of any statistic T which operates on the population where Y comes from?  ### Here's what I tried so far:               Y<-data.frame(response=c(10,19,27,28,9,13,25,29,4,10,20,18,5,6,12,17),                    treatment=factor(rep(1:4,4)),                    subject=factor(rep(1:4,each=4))                    )          p<-matrix(Y$response,4,4,byrow=T)     B<-1000     sampleB<-sample(1:4,4*B,replace=T)     fit<-lm(p[sampleB,]~1)     cov(residuals(fit))      I also tried               require(nlme)     require(mgcv)     nSubj <- 20     sampleB<-sample(1:4,nSubj,replace=T)     y<-data.frame(response=c(t(p[sampleB,])),                treatment=factor(rep(1:4,nSubj)),                subject=factor(rep(1:nSubj,each=4))                )          fit <- lme(response~-1+treatment,y,random=~1|subject,correlation=corSymm())     extract.lme.cov(fit,y)[1:4,1:4]      but I get the error code:               Error in lme.formula(response ~ -1 + treatment, y, random = ~1 | subject,  :      nlminb problem, convergence error code = 1     message = iteration limit reached without convergence (10)


title: Gamma distribution from pessimistic/most-likely/optimistic values

text: Is it possible to obtain a Gamma distribution from knowing the pessimistic, most-likely and optimistic value of a random variable? The values are known with a 90% confidence, i.e. nine times out of ten, the actual value of the random variable will fall within the pessimistic-optimistic bound while it will be skewed to the most-likely value. **Question:** How can I compute the shape and scale parameters for Gamma distribution from the lower and upper limits and the mode?


title: Logistic regression with categorical predictor

text: I'm facing a problem with a binomial `glmer` model. I want to find if differences in flower presence in pine trees is due to procedence of the tree. My model is as follows: `FlorMas ~ Proc + (1|Blq)`. Proc is a factor with nine levels, one of it (`TAMR`) presents no flower at all (variable value for all `TAMR` trees is 0). This model gives me this output:               Generalized linear mixed model fit by the Laplace approximation      Formula: FlorMas ~ Proc + (1 | Blq)         Data: flower.data       AIC   BIC logLik deviance      593 647.7 -285.5      571     Random effects:      Groups Name        Variance Std.Dev.      Blq    (Intercept) 0.18476  0.42983      Number of obs: 1067, groups: Blq, 8          Fixed effects:                  Estimate Std. Error z value Pr(>|z|)         (Intercept)   -1.7668     0.2958  -5.974 2.32e-09 ***     ProcTAMR     -16.8758  1080.5608  -0.016  0.98754         ProcARMY      -0.3543     0.3910  -0.906  0.36490         ProcASPE      -1.4891     0.5260  -2.831  0.00464 **      ProcCOCA      -2.4947     0.7619  -3.274  0.00106 **      ProcMIMI      -1.2040     0.4930  -2.442  0.01459 *       ProcORIA      -1.5360     0.5739  -2.676  0.00744 **      ProcPLEU      -1.9437     1.0538  -1.845  0.06511 .       ProcPTOV       0.1693     0.3508   0.483  0.62945         ProcSCRI       0.5060     0.3346   1.512  0.13050              Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      I don't understand that values for `TAMR` procedence, as if it has all zero values it should be different from the others. Any help will be appreciated.


title: Comparing cases based on multiple binary variables

text: A researcher wants to compare 50 websites. There is 80 binary variables for each website that belong to 4 main measures. Actually There is 50*80 matrix that each row represent a website and each column is a binary variable. He asked me a statistical method for comparing (and ranking) the websites based on these data. Maybe we can use the main 4 measures (as sum of related variables) instead of 80 variables. Thanks for your suggestions...


title: Sum of exponential random variables follows Gamma, confused by the parameters

text: I've learned sum of exponential random variables follows Gamma distribution. But everywhere I read the parametrization is different. For instance, Wiki describes the relationship, but don't say what their parameters actually mean? Shape, scale, rate, 1/rate? Expoenetial distribution: $x$~$exp(\lambda)$ $$f(x|\lambda )=\lambda {{e}^{-\lambda x}}$$ $$E[x]=1/ \lambda$$ $$var(x)=1/\lambda$$ Gamma distribution: $\Gamma(\text{shape}=\alpha, \text{scale}=\beta)$ $$f(x|\alpha ,\beta )=\frac{1}{{{\beta }^{\alpha }}}\frac{1}{\Gamma (\alpha )}{{x}^{\alpha -1}}{{e}^{-\frac{x}{\beta }}}$$ $$E[x]=\alpha\beta$$ $$var[x]=\alpha{\beta}^{2}$$ In this setting, what is $\sum\limits_{i=1}^{n}{{{x}_{i}}}$? What would the correct parametrization be? How about extending this to chi-square? Thanks!


title: Weighting a control group

text: In reviewing some market research, I came across a study which decided to weight the results/outcomes of the control group so the distribution of participants in the control sub-groups (0,1,2,3,4,5 - a measure of their pre- experiment activity of a certain activity) would equate with that of the treatment group. The experiment was designed to introduce an incentive to encourage purchase activity. Both treatment and control groups were supposed to be selected at random (N is in the 1000s for each group) - what would be the proper way(s) to address differences in the pre-distributions across groups when analyzing the results/outcomes of the test?


title: Using total score from multi-scale instrument in structural equation modeling

text: I wonder if you can help me to get the right answer to a question about structural equation modeling. Imagine someone trying to validate a questionnaire using component factor analysis that highlighted 5 subscales. Is it right to use the total score of the questionnaire in a model created by structural equation modeling? That is, would it be correct to consider that the questionnaire is 'validated' and use its total score? Would you think that this model is right or valid or should the subscales be considered as observed variables and total score as latent variable in the model?


title: Interpretation of mixed model with repeated measures

text: I have performed the following experiment: Two genotype groups of two different ages were monitored in the presence of three drugs one after the other. So, each subject after a control period received drug A, then drug B and then drug C. I want to know what is the effect of every drug compared to control and between each of them and also check for drug*genotype*age interactions. As far as I understand, I have a mixed model with one within- subjects factor (drug) and two between-subjects factors (genotype and age). I'm using SPSS and I'm running repeated measures test (responses follow normal distribution and Levene's test of equality of error variances is not significant). **My question is:** How many levels should be included in the within-subjects factor? Should I run one test with 4 within-subjects levels or run a separate test for each comparison between drugs with 2 levels in each test? I think that the correct procedure would be to run one test with 4 levels, but, I am concerned about the question that such a comparison asks. The reason is, that none of the drugs are expected to work in the same way. Some are expected to show reduced response and others increased response. So, if I run such a comparison with all 4 levels I feel I'm asking if _all_ of the drugs affect the two genotypes and ages in a similar way, which is different to my original question of finding differences between _each_ of the drugs. Many thanks


title: How to fit the model for crossed and nested design using lme function in R?

text: Suppose A is a fixed factor and B is a random factor. A and B are crossed. Another random factor C is nested in B. Taking into account the AB and AC interactions, how should I fit the model using lme function in the nlme package? Thanks.


title: Dummy variable that is a sub-category of another dummy variable

text: I'm wondering if there are any statistical issues with using a dummy variable that is a subcategory of another dummy variable. For example: Let's say I am trying to predict a binary outcome of whether a child goes to college. If I have one dummy variable "NYC" for "lives in New York City" and another "Bronx" for "lives in the Bronx" (which is part of NYC) then "NYC" will always be 1 when "Bronx" is 1 but not vice versa. As I understand it this shouldn't create a "dummy variable trap" but is this a problem in some other way?


title: Avoid rater tendency in confidence questionaires

text: For the development of a new chemical process, a handful experts should be asked a few questions (4 to 6), what values they expect for several experimental parameters, i.e. temperature ranges and how confident they are regarding their guesses on a scale from 1 to 10. I would like to weight their responses corresponding to the confidence. Some people tend to be very confident with average ratings of 8, other tend to be less confident with ratings of 3 to 4. To avoid such things as the Dunning-Kruger effect, I would like to get a rating of them, what confidence level means what and normalize their responses. I thought of giving them a questionnaire with questions, they know for pretty sure, and others, where they most likely have to guess and take the confidence level their as well. If people give high confidence scores for questions, they could not possibly know, I can assume, that an 8 is more something other would rate with a 5. I now wonder, if such a questionnaire to get the rater’s tendency in confidence scores already exists, possibly already widely evaluated? I am not sure, if Cross Validated is the right corner of Stack Exchange to post this question, if not I am very sorry; feel free to delete this question then.


title: Time Series: Seasonality and trend

text: I am interested in financial time series and I have a small question regarding the use of the forecast package. The time series I am interested in is a monthly one and present clear evidences of seasonality. Thus there is a need to use a lag 12 differentiation. After this step, I apply a full collection of KPSS and ADF tests on the resulting time series and the conclusion is that it is stationnary up to a linear trend. However if I perform those tests before applying the lag 12 diff operator, the conclusion is that the time series is not stationnary and needs to be differentiated. Here are my questions: 1) Is it better checking the stationnarity through the tests before or after differentiating at lag 12 ? 2) If it is better (as I do for the moment) to differentiate at lag 12 priorly, what should I do if the outcomes of the tests are that the resulting time series is stationnary up to a linear trend ? I guess I cannot simply use auto.arima with D=1 and xreg=1:N since the linear regression would then be done before differentiating at lag 12. Thanks.


title: conditional mean of a geometric RV

text: Say, there are three nodes: $S$, $R$, $D$. $S$ transmits to $R$, $R$ stores the packets, and later transmits to $D$. At any time, either $S$ or $R$ is selected to transmit according to some random selection rule. The probability that $S$ is selected to transmit is $P_S$, and that of R is $P_R$. For e.g., if $S$ is selected repeatedly, $R$ has to store the packets in its queue which it transmits later when it's turn arrives. We assume that the queue size is $L$. Let us model a geometric RV $X$ which refers to the number of time $R$ selected before $S$ is selected. In that sense, $P_S$ is the probability of success and $P_R = 1-P_S$ is the probability of failure, according to conventional terminology. The average value of $X$ is given by $P_R/P_S$. However, if the $R$ is selected to transmit repeatedly, it is possible that the queue at $R$ becomes empty. So until $S$ is selected again, there is no transmission. The vice-versa is also possible, but I'm not interested in it. I'm trying to figure out the average value of the number of time $R$ is selected when the queue is empty. I have an expression of probability of empty queue which is a function of $P_R$ and $L$. So basically I have to find a conditional mean of $X$. Any help or direction will be highly appreciated.


title: Variance of predictions from lme4 models at existing levels?

text: I know that this has been widely discussed, but even after having read many discussions on this matter, I'm still not sure if I'm understanding properly. So I have a dataset of studies reporting prevalence in several settings, here is an exemple:               data<-data.frame(id_study=c("UK1","UK2","UK3","FRA1","FRA2","BEL1","GER1","GER2","GER3"),                  country=c("UK","UK","UK","FRANCE","FRANCE","BELGIUM","GERMANY","GERMANY","GERMANY"),                  n_events=c(130,49,18,10,50,6,34,20,51),                  n_total=c(3041,580,252,480,887,256,400,206,300),                  study_median_age=c(25,50,58,30,42,26,27,28,36),                  pop_type=c(2,0,1,0,0,0,1,0,2))     data$prevalence<-data$n_events/data$n_total      Since the prevalence is linked to the sample age and type, the objective is to estimate the predicted prevalence for each study for a sample with standardized median age (30) and type ("A"). Then these predictions will be pooled by country using a Dersimonian & Laird method. This is the model I chose:               mod<-glmer(cbind(n_events, n_total-n_events) ~ study_median_age + pop_type +          (1|country) + (1|id_study), data=data, family="binomial"(link=logit))      I tried to adapt the method from the r-sig-mixed-models FAQ in order to get the variance of the estimate, which will be used in the weighting of the Dersimonian & Laird meta-analysis:               newdat=data.frame(id_study=data$id_study,                       country=data$country,                       study_median_age=30,                       pop_type="A",                       n_events=0,                       n_total=data$n_total)          pred<-predict(mod,newdat[,1:4],type="link",REform=~(1|country) + (1|id_study))     mm <- model.matrix(terms(mod),newdat)     newdat$distance <- mm %*% fixef(mod)          pvar1 <- diag(mm %*% tcrossprod(vcov(mod),mm))     tvar1 <- pvar1 + VarCorr(mod)$country[1] + VarCorr(mod)$id_study[1]        If I understood correctly, tvar1 is an estimation of the variance of the prediction, taking into account the uncertainty on the two random intercepts. I can consider tvar1 as the variance of the prediction for the further meta- analysis (with inverse variance weighting), and also for calculation of the confidence interval of the prediction.               data$prediction<-mod@resp$family$linkinv(pred)     data$predictionSE2<-tvar1     data$predictionCI<-paste0("[",round(mod@resp$family$linkinv(pred-1.96*tvar1),4),";",                               round(mod@resp$family$linkinv(pred+1.96*tvar1),4),"]")     data      So my questions are:   * Did I properly adapt the method for the calculation of the variance to this model with 2 random intercepts?   * Can I use this variance in the meta-analysis? Thank you!


title: Predict Coefficients with glmnet()

text: ridge.fit1 <- glmnet(x,y, alpha = 0)            ridge.fitt <- coef(glmnet(x, y, alpha = 0, lambda = 3))     ridge.pred <- predict(ridge.fit1, s = 3, type = 'coefficients')      What's the difference between ridge.fitt and ridge.pred?


title: Mean of covariance matrices

text: I'm trying to generalise a formula that takes the mean of some variances to it work with vectors. I'm not sure it makes sense to take the variance between a bunch of vectors, rather it is more suited to generate a covariance matrix. So, is there a way to take the mean of a bunch of covariance matrices?


title: Reading only two out of three columns with read.csv

text: I have an ascii dataset which consists of three columns, but only the last two are actual data. Now I want to create a dotchart of the data by using `read.csv(file = "result1", sep= " ")`. R reads all three columns. How do I avoid this?


title: How are the standard errors computed for the fitted values from a logistic regression with parameter 'response' in R?

text: I am looking for the equation for the se.fit values when using logistic regression in R. I have seen this answer - How are the standard errors computed for the fitted values from a logistic regression? but in my case, I'm calling "predict" function with 'type' parameter set to be "response". In this case, the equation given in the linked I attached doesn't hold. Here is an example of the predict function I'm calling:               predicted.resutls <- predict(glm.model, train.data, type = "response", se.fit=TRUE)


title: How can I perform joint Schur decomposition or joint spectral decomposition?

text: I have a collection of matrices $A_1,\ldots,A_n$, and I want to do joint Schur decomposition on them (they all have the same unitary matrices in their decomposition). I couldn't find any implementation of this (I was looking for Python, but I would be happy to find anything at the moment). The closest I found was an implementation that works for two matrices ($n=2$) for Python. Anyone is familiar with how to do that or has a Python/other implementation? Also, is there more stable algorithms to do joint decomposition to identify the eigenvalues of $n$ matrices simultaneously? EDIT: My matrices are diagonalizable, so a joint spectral decomposition into eigenvectors and eigenvalues (when the eigenvectors are all the same) would be fine as well.


What is a Gaussian Discriminant Analysis (GDA)? What materials should one read to understand how a GDA works and where it comes from? Try to explain this for someone at a high-school level.


title: How to compare two groups on a set of dichotomous variables?

text: I have two groups (G1, n=10; G2, n = 10) each representing a separate condition. Participants in each group answered 20 questions and each question is a dichotomous variable coded 0 and 1 (VDD). I want to compare the group 1 with group 2. I am having some trouble understanding if I have it right, for every participants of both group, to mean their answer (since the variable is dichotomous). It would give me a probability to get an answer more than the other one I guess, but I don't know if I have the right to do that.


title: Can you recommend an online survey platform for 5k+ participants?

text: I am planning to perform an online behavioral survey across a nationwide sample, and I expect several thousand responses. I expect to have not that many questions (perhaps 3 pages, 8 qs each), require skip logic, and would like it to redirect back to a page of my choice. I also need it to be hosted, as local webhosting is not an appealing option for me. I do not need reports or fancy user interface. Has anyone been involved in a survey of this scale, and do you have a platform you recommend?   * I am aware of REDCap but my institution does not have access to it.   * SurveyMonkey's Gold plan appears to meet my needs, but I would like to review alternatives   * AFAIK Google Docs Forms doesn't have skip logic   * LimeSurvey was recommended in a related answer Thanks for any help.


title: Density estimation with a truncated distribution?

text: I have some data which is clearly truncated on the left. I wish to fit it with a density estimation that will handle it in some way instead of trying to smooth it down. What known methods (as usual, in R) can address this? Sample code:               set.seed(1341)     x <- c(runif(30, 0, 0.01), rnorm(100,3))     hist(x, br = 10, freq = F)     lines(density(x), col = 3, lwd = 3)      ![enter image description here](http://i.stack.imgur.com/MDPOL.png) Thanks :)


title: How do I explain that software implemented model selection procedures should not be used unsupervised?

text: I know that people generally say that procedures which select a model based on information criterion lead to inconsistent model selections. I read a paper by Leeb and Potscher (2005), MODEL SELECTION AND INFERENCE: FACTS AND FICTION which supports this theory. Regardless of the kind of information criterion we use (a consistent one like BIC or a conservative criterion like AIC) we will sometimes get seriously biased variables. As there is no way to actually check if our model is biased in this way, since the examples provided by Leeb and Potscher depend on their sample design (even the design is a very general one), we could be very dull and use the model selection procedure nevertheless. Another more common issue is that an omitted variable biases the coefficients as well, if a regressor is correlated with the error term. So I would say that one should always check the changes in coefficients estimated if a regressor gets kicked out. Even if the relationship between an omitted regressor and the dependent variable is very low (which is the reason why it got kicked out) it would still be a viable source of bias. So one must not only check the usual significance criterion, but also check the changes in coefficients, to minimize the risk of biased coefficients?


title: Calculating demand size and intervals

text: could you explain how to determine a the mean and standard deviation of the following: demand size, demand interval and demand per period, with respect to minimum, 25%ile, median, 75%ile and maximum from the data below:   Weeks: Orders: Inventory               0 :               : 112224         04/08/2009: 2000    :110224         18/08/2009: 2000    :108224         01/09/2009: 8000    :100224         15/09/2009: 2040:   98184         06/10/2009; 20160;  78024         13/10/2009; 10080 ; 67944         20/10/2009; 15960;  51984         27/10/2009; 10080;  41904         03/11/2009; 18000   ;136128         10/11/2009; 13440;  122688         17/11/2009; 7200;   115488         24/11/2009; 2040;   113448         08/12/2009; 20880;  92568         05/01/2010; 25000;  67568         12/01/2010; 16000;  51568         19/01/2010; 20600;  30968         26/01/2010; 16120;  14848         02/02/2010; 20000;  -5152         09/02/2010; 20000;  -25152         16/02/2010; 20000;  -45152         02/03/2010; 18000;  49072         09/03/2010; 18000;  31072         16/03/2010; 8040;   23032         23/03/2010; 10080;  12952         30/03/2010; 12000;  952         06/04/2010; 15960;  -15008         13/04/2010; 10080;  -25088         20/04/2010; 12000;  -37088         27/04/2010; 15960;  59176         04/05/2010; 12960;  46216         11/05/2010; 4400;   41816         18/05/2010; 14000;  27816         25/05/2010; 11200;  16616         01/06/2010; 13240;  3376         08/06/2010; 3340;   36         15/06/2010; 8140;   -8104         22/06/2010; 8260;   95860         29/06/2010; 12000;  83860         06/07/2010; 4060;   79800         13/07/2010; 9000;   70800         20/07/2010; 14000;  56800         27/07/2010; 6840;   49960         03/08/2010; 8040;   41920         10/08/2010; 15240;  26680         17/08/2010; 15240;  123664         24/08/2010; 8040;   115624         31/08/2010; 16040;  99584         07/09/2010; 5080;   94504         14/09/2010; 8320;   86184         21/09/2010; 9880;   76304         28/09/2010; 8040;   68264         05/10/2010; 9280;   58984         12/10/2010; 3880;   55104         19/10/2010; 8040;   47064         26/10/2010; 16520;  142768         02/11/2010; 15000;  127768         09/11/2010; 960;    126808         16/11/2010; 17040;  109768         23/11/2010; 17040;  92728         30/11/2010; 17040;  75688         07/12/2010; 17040;  58648         14/12/2010; 17040;  41608         21/12/2010; 16520;  25088         28/12/2010; 14760;  10328         04/01/2011; 10080;    248         11/01/2011; 11160;  101312         18/01/2011; 5760;   95552         25/01/2011; 16520;  79032         01/02/2011; 16480;  62552         08/02/2011; 16520;  46032         15/02/2011; 1080;   44952         01/03/2011; 16520;  28432         08/03/2011; 17400;  11032         15/03/2011; 17400;  105856         22/03/2011; 3240;   102616         29/03/2011; 17520;  85096         05/04/2011; 14400;  70696         12/04/2011; 14400;  56296         19/04/2011; 13560;  42736         26/04/2011; 16080;  26656         03/05/2011; 13680;  12976         10/05/2011; 17280;  107920      Demand sizes is the sizes of demand when demand occurs, Demand per period is all periods includes demand occurring and nondemand occurring. The time bucket is 1 month, the history length is 22 months, and the lead time is 8 weeks


title: How to model continuous predictor variable, $x\ge0$, showing continuity "jump" at $x=0$?

text: I'm modeling hospital inpatient stays with a logistic regression model, and have encountered a number of predictor variables measuring # of outpatient care episodes which have an inflated number of observations at $x=0$ (sometimes over 90%). In addition, there's a discontinuous "jump" in observed y values comparing $x=0$ with $x>0$. Plotting a decile chart with $x=0$ to $10$ and looking at the proportion of people with inpatient stays within each decile reveals the following pattern as a typical example: ![enter image description here](http://i.stack.imgur.com/LI1Pc.jpg) One solution I'm tinkering with is to describe x with two component variables,   * a binary variable $x_1$ where $x_1=0$ if $x=0$, $x_1=1$ if $x>0$, and    * a continuous variable $x_2=x$.  Then describe the model equation as: $\text{logit}(y)=\beta_0 + (\beta_1 + \beta_2x_2)x_1 = \beta_0 + \beta_1x_1 + \beta_2x_1x_2$ `If x1=0 then logit(y) = beta0. Else if x1=1 then logit(y) = (beta0 + beta1) + beta2*x2`. However I'm unsure if this is a standard procedure for modeling predictor variables that exhibit this behavior?


title: Evaluation on Binary Outcome data

text: Let us say that we have an experiment where people are shown images of tyres from car/bike and have to tell what it really belongs to. The objective of the experiment is to validate the hypothesis that people from age group 25-35 are much better at recognizing that the tyre images are from car or bikes. So the final data that is obtained is as follows Sample data for the same would be : Bike = 0 Car = 1               -------------------------------------     AGE    Actual_Image  Predicted_Image      -------------------------------------         15         0              1     27         1              1     34         0              0     55         1              0      What is the correct way to see if the hypothesis is correct or not ? Since I am a beginner in this field some links would be appreciated. I have looked at `chi-square test` for the same but am unable to understand how to apply it in this problem. As a secondary question I would also like to know if `d-prime` value can help in this ? For example if I calculate d prime value for each group and say that the group having the highest d-prime value is the best at recognizing the images ? P.S : Not any homework/test. Looking to learn more about validating hypothesis/ significance value etc.


title: How important is domain knowledge in our profession?

text: _**or:_** _Does selecting a domain when entering a job narrows your future options for domains and hence jobs ?_ To make this question as broadly applicable as possible ...   * profession refers to all kinds of data analysts, from statisticians over machine learner programmers to data miners.    * imagine you were asked to give an advice to an audience containing both students and professionals of different age classes **Maybe a starting point:** The competitions on Kaggle have shown that outsiders can outperform the models created by company employees (see e.g. here). On the other hand, my (limited) work experience has led me to the conclusion, that understanding how and where the data has been generated is absolutely mandatory to create an abstract environment where something like a Kaggle competition can happen. Furthermore, without domain knowledge, I find it hard to report the results to other layers / departments. Some associate the last skill as key to the "new" profession "Data Science" (see e.g. here or here).


title: Analysis of forest inventory data - non-random samples

text: I apologise, this isn't a single question. It is more of a general problem on which I am working and am seeking guidance for how to proceed. I have been provided with an inventory dataset of plot measurements (tree height, diameters etc) within a non-random selection of n forest stands (from a total N in the complete forest estate) that are dispersed spatially over a large area. I have point estimates for the forecasted yield per hectare for each of these stands and I can calculate the observed yield per hectare based on the inventory data. The 0.02 ha plots were located at random within each forest stand aiming at a sampling intensity of approximately 1 plot per 6 hectares. As such I believe that these plots represent a simple random sample and can be analysed accordingly to give me **stand** specific information such as mean yield per hectare. I believe that I could do this for each forest stand without too many problems. I am trying to make sure I understand the limitations of the data with a view to being able to make valid statements at the stand level, the aggregated stands and the complete forest estate. I would like to be able to statements such as:   1. The mean yield for stand n(i) is 200m3/ha with an upper and lower confidence interval bound (a,b).   2. The predicted yield for stand n(i) overestimates/equals/underestimates (lies outside the confidence interval or some more appropriate criteria) of the observed yield by x%.   3. On average the predicted yield across **all** stands overestimates/equals/underestimates   the observed average yield by x%.   4. Region appears to be an important factor in whether the forecasted yield overestimates/equals/underestimates the observed value.   5. Stand age appears to be an important factor in whether the forecasted yield overestimates/equals/underestimates the observed value. And so on... While the n forest stands were a non-random selection, I believe that I could combine the forest stand estimates that I calcuated and then assess (or at least produce) the distribution of (say) the mean yield per hectare and the variation across this n stands or the distribution of the mean difference between the observed yield and the forecasted yield. This would give me some insight into the variability across the selection of stands and perhaps give me some basis on which to suggest whether the forecasted yields are appropriate. However, due to the non-random stand selection procedure, any conclusions that I make are not applicable to the estate level. Is this a correct assumption or am I already making too many assumptions about my ability to combine the estimates I calculated at the stand level. One side effect with calculating the stand level averages is that the resulting number of data points that I have available reduces markedly. This means that comparing factors that might influence differences between the stands becomes difficult. I am thinking that it may be useful to take the entire set of plots (i.e. across all the stands) and treat them as a simple random sample but I do not understand what the implications of this approach are. I am erring towards re-sampling without replacement from all the plots available to me because I think that would at least put me on firm ground for making inference about the set of stands that have been included in the inventory. Appreciate your time and comments.


title: Naive Bayes Classifier - measure accuracy after training

text: I have built a prototype Naïve Bayes Classifier in an Excel spreadsheet. My data is a transaction (an order) with 13 parameters. This translates directly to a feature vector of (feature_1, feature_2, …, feature_13). Each parameter of an order, if it is incorrectly set, can be the cause of an error in the order. I am trying to classify each order according to the cause of the error in the order. So I’m modelling 14 classes: each of the parameters, plus an ‘other’ class for some unknown error that doesn’t fit into any of the other classes. My question is: how do I measure the increase in accuracy of the system after it is trained with each new order? And on a related note, should I expect the average increase in accuracy (after each new order) to drop over time? I.e., the system will get less and less value out of each new feature vector that is fed into it?


title: Convergence in Distribution\CLT

text: Given that $N = n$, the conditional distr. of $Y$ is $\chi ^2(2n)$. $N$ has marginal distr. of Poisson($\theta$), $\theta$ is a positive constant. Show that, as $\theta \rightarrow \infty$, $\space \space (Y - E(Y))/ \sqrt{\operatorname{Var}(Y)} \rightarrow N(0,1)$ in distribution. Could anyone suggest strategies to solve this. It seems like we need to use CLT (Central Limit Theorem) but it looks tough to get any information on $Y$ on it's own. Is there a rv that can be introduced to take a sample of, to generate $Y$? This is homework so _hints_ appreciated.


title: Implementing a Logit Model With Multiple Predictors

text: Let's say I have the following equation:               Won = B0 + B1*(Bid)      Once I know B0 and B1, I can generate the probability curve and find the probability of "won" for each "bid." However, let's say I add a variable zipcode, and add it to the logit model. In R, I save it as a factor and end up with:               Won = B0 + B1*(Bid) + B2*(zip5001-6000) + B3*(zip6001-7000) + B4*(zip7001-8000)      Once I find all the coefficients, I want to find the probability curve for each zipcode, such that I can say, "when the bid is x and the zipcode is y, the probability of winning is f." What I assumed is that for B2, I would take:               Won = B0 + B1*(Bid) + B2*(zip5001-6000) + 1*(zip6001-7000) + 1*(zip7001-8000)      However, this doesn't give me the probability of "won" when bid is x and B2 is 5001-6000. I'm doing this by adding the coefficients to the general equation for the probability curve. How would I be able to generate a curve for the desired probabilities and how should the equation look? The probability equation I am using is = $$\text{Prob} = \frac{1}{1 + e^{-z}}$$ where $$z = B_0 + B_1X_1 + \dots + B_nX_n.$$


title: Evaluation of Autocorr/Part Autocorr values

text: I am practicing MA and AR modelling by using autocorrelation and partial autocorrelation values. My data is in the image below; I can see that only at lag 12 there is a value that might be considered but I have a dilemma: what model is this if any?   1. To be a MA of higher order the autocorr should be 0 After lag 12, but not before.   2. To be an AR I should see value oscillating and decaying up to lag 12. Am I correct in my thinking or am I wrong? I performed the Box-Jenkins test (`lbqtest`) with Matlab and only at lag 12 I have the rejection of the null hypothesis ( h =1 ) for an Alpha of 0.10. ![enter image description here](http://i.stack.imgur.com/b9wAG.png)


title: Sampling with replacement in R randomForest

text: The randomForest implementation does not allow sampling beyond the number of observations, even when sampling with replacement. Why is this? Works fine:               rf <- randomForest(Species ~ ., iris, sampsize=c(1, 1, 1), replace=TRUE)     rf <- randomForest(Species ~ ., iris, sampsize=3, replace=TRUE)      What I want to do:               rf <- randomForest(Species ~ ., iris, sampsize=c(51, 1, 1), replace=TRUE)     Error in randomForest.default(m, y, ...) :        sampsize can not be larger than class frequency      Similar error without stratified sample:               rf <- randomForest(Species ~ ., iris, sampsize=151, replace=TRUE)     Error in randomForest.default(m, y, ...) : sampsize too large      Since I was expecting the method to take bootstrap samples when given replace=TRUE in both cases, I was not expecting this limit. My objective is to use this with the stratified sampling option, in order to draw a sufficiently large sample from a relatively rare class.


title: Covariance of $\mathbf{x}$ and $\mathbf{x}^\prime A \mathbf{x}$ when $\mathbf{x} \sim N_{m}(\mu, \Sigma)$

text: This is part of problem 5.23 in A First Course in Linear Model Theory, Dey and Ravishanker. It was on a previous midterm and I didn't know how to do it, but now I am studying for the final and would like to figure it out. Suppose $\mathbf{x} \sim N_{m}(\mu, \Sigma)$ and $A$ is a symmetric idempotent matrix with rank $p<m$. Edit: Show that $Cov(\mathbf{x}, \mathbf{x}^\prime A \mathbf{x})=2\Sigma A \mu$ I can find $E(\mathbf{x}^\prime A \mathbf{x})$ and $Var(\mathbf{x}^\prime A \mathbf{x})$ using the formula for the cumulants of $\mathbf{x}^\prime A \mathbf{x}$, $$ \kappa (\mathbf{x}^\prime A \mathbf{x})=2^{r-1}(r-1)![tr(A\Sigma)^r+r\mu^\prime A(\Sigma A)^{r-1}\mu] $$ but I don't see how I can get the covariance above using the cumulants without more information. Somehow I need to calculate $E(\mathbf{x} \mathbf{x}^\prime A \mathbf{x})$ but I can't figure out how to find this quantity. Hints at how to get started would be great.


title: Evaluating methods for news recommendation from an online experiment

text: I have data from an experiment that evaluated 3 methods for recommendation of news articles based on their content. The experiment showed an article that the user read and 6 recommendations presented in random order below it (screenshot: http://cl.ly/image/2R10073Z0i39). The task of the participants was to choose one recommendation that they found the most interesting. The 6 recommendations were generated by 3 methods, each creating 2 recommendations. I evaluated precision, recall and the F1 measure. But I received feedback that since the task was to choose one recommendation from 6, the users didn't have a chance to respond to every recommendation and so precision can't be applied. I also calculated the R-score (Expected Utility) of the methods (the problem here is that the user saw only 2 recommendations from each method at once). My question is, what would be the best way to evaluate the 3 methods based on the data from this experiment? I have data about user clicks, read articles and the recommendations that the users saw (actually I have top 10 recommendations by each method at each point - including the ones the users didn't see). Thank you.


title: Measures of autocorrelation in categorical values of a Markov Chain?

text: **Direct Question:** Are there any measures of auto-correlation for a sequence of observations of an (unordered) categorical variable? **Background:** I'm using MCMC to sample from a categorical variable and I'd like a measure of how well the sampling method I've developed is mixing across the posterior distribution. I'm familiar with acf plots and auto-correlation for continuous variables, but I've been stuck looking at the transition probability matrix for this categorical variable... Any thoughts?


title: R Prefmod for paired comparisons

text: I'm trying to use the prefmod package to analyzed paired comparisons, but unfortunately the documentation is a bit poor and the author is no longer with us. I have a survey with 32 questions, split into two blocks (16 questions per respondent). Here's a sample of the first 10 observations:                   V24 V27 V29 V31 V33 V35 V37 V39 V41 V43 V45 V47 V49 V51 V53 V55 V57 V59 V61 V63 V65 V67 V69 V71 V73 V75 V77 V79 V81 V83 V85 V87 V89 block     1    2   2   2   2   2   2   2   2   2   1   2   1   2   1   2   1   1  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA     1     2    2  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA   2   1   1   1   1   2   2   2  NA   2   1   1   2   2   1   1     2     3    2   2   2   1   2   2   2   2   1   1   2   1   1   1   2   1   1  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA     1     4    3   2   2   2   2   2   2   2   2   1   2   1   2   1   2   1   1  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA     1     5    3  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA   2   1   1   1   1   2   2   2   2   1  NA   2   1   2   1   1     2     6    2  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA   2   2   1   1  NA   2   2   2   1   2   1   1   2   2   1   1     2     7    1   2   2   2   2   2   2   2   2   1   2   1   2   2   2   1   1  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA     1     8    3   2   2   2   1   2   1   2   2   1   2   1   2   1   2   1   1  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA     1     9    2   2   2   2   2   2   2   2   2   1   2   1   2   1   2   1   1  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA     1     10   2   2   2   2   2   2   2   2   2   1   2   1   2   1   2   1   1  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA     1      V24 is an individual-specific variable, block is what block they were randomly assigned to (1 or 2), but the rest are their choices over the pairs. How do I get this to run in prefmod? I have attempted to do use the package's llbt.design command, but the "number of items" doesn't make sense because unlike the examples, there aren't a set number of objects here - the choices have 3 attributes with 4 levels each, for a total of 64 possible scenarios (48 appear amongst all the questions.


title: How to interpret an inverse covariance or precision matrix?

text: I was wondering whether anyone could point me to some references that discuss the interpretation of the elements of the inverse covariance matrix, also known as the concentration matrix or the precision matrix. I have access to Cox and Wermuth's _Multivariate Dependencies_ , but what I'm looking for is an interpretation of each element in the inverse matrix. Wikipedia states: "The elements of the precision matrix have an interpretation in terms of partial correlations and partial variances," which leads me to this page. Is there an interpretation without using linear regression? IE, in terms of covariances or geometry?


title: Test for accuracy of prediction of financial time series

text: I would like to test a proprietary method of financial time series prediction that says, in effect, that on a certain day/week in the future, to an accuracy of +/- 2 to 3 days/weeks, a financial time series will make a relative high or low on this predicted date, and furthermore, within a defined period in the future there will only be a certain number of such relative highs and lows. What statistical tests exist that can ascertain whether such a predictive methodology is non-random and actually has predictive ability in as much as it accurately calls the highs and lows? For the purpose of such test(s) the years 2007 to 2011 inclusive can be considered unseen, out of sample data. One idea I have is to plot the data and count how far the actual highs and lows are from their predicted dates of occurrence and devise an error metric and then randomly permuting the data to obtain new time series with random highs and lows, applying the same error metric and then comparing the real error metric with the distribution of the random error metrics.


title: R - Clustering and Multidimensional Scaling

text: I am trying to get to grips with Clustering and Visualisation. I have a set of data (a matrix) that I want to cluster (using R) and then visualise (HTML5 Canvas). So, I can use MDS to get the coordinates of a matrix, for example:               cells <- c(1, 1, 2, 1, 4, 3, 5, 4)     rnames <- c("A", "B", "C", "D")     cnames <- c("X", "Y")     x <- matrix(cells, nrow=4, ncol=2, byrow=TRUE, dimnames=list(rnames, cnames))     d <-dist(x)     m <- cmdscale(d,eig=TRUE, k=3)     m <- cmdscale(d,eig=TRUE, k=2)     print(m)      But, should I first cluster this data (using something like k-means) or should I cluster the output of MDS? Or can I just cluster (and get the coordinates using some other method) and ignore MDS? What is the relationship between MDS and K-means, if any? I am unsure what is the best way to approach? Any suggestions? Thanks, s


title: Probabilities and Linear Regression

text: I am having some trouble understanding the relationship between the weight provided from a multiple linear regression and calculating probabilities. The weights given create an equation whose values are not part of a probability distribution. Rather, they predict the y value. How can this be transformed into a probability distribution? If I have my weights, and all of the attributes which are weighted occur, then I should have a probability of 1 of the event happening (likewise a p of 0 if all the weighted attributes do not occur). I am more concerned with calculating those probabilities when only some of the attributes occur, given that the more heavily weighted attributes should be given more, well, weight.


title: General linear hypothesis testing in repeated-measures multivariate linear model

text: In a multivariate linear model $Y = XB + U$, where within-subject (or repeated-measures) factors are coded as simultaneous response variables in $Y$, a general linear hypothesis can be formulated as $H_0: L B M = 0$, where $L$ specifies the contrasts among the levels of a between-subjects variable while $M$ codes the contrasts among the levels of a within-subject factor, $B$ and $U$ ares the regression coefficient and residual matrices respectively. The traditional univariate testing (e.g., ANOVA) for an effect associated with a within-subject factor can be computed through the sum of squares and cross products (SSP) for the above hypothesis ($H$) and for the errors ($E$), assuming sphericity. From the references I have access to (e.g., SAS documentation), it seems that the computations are done through an orthogonalization process of $M$, and then each $F$-statistic is calculated as $trace(H)/trace(E)$. Unfortunately all the references I've read so far only provide some basic descriptions without any formulas. Thus I have two questions: 1) Is there a formula in which the $F$-statistic of the trace ratio can be expressed directly in terms of the $M$ matrix without the orthogonalization? I vaguely remember seeing somewhere a formula like $F = trace[H (M' M)^{-1}] / trace[E (M' M)^{-1}]$, but I can't recall the literature. Could someone provide a reference for the formula? 2) If there is no such explicit formula, how is the orthogonalization exactly done in SAS/SPSS?


title: Oversampling in logistic regression

text: I was trying to find out whether an oversampling can really make a model better. On this blog page, it says that it can improve a decision tree, but it shouldn't improve a logistic regression. Quotation below: > _Standard statistical techniques are insensitive to the original density of > the data. So, a logistic regression run on oversampled data should produce > essentially the same model as on the original data._ **Question:**   Is it really true that using an oversampling in a logistic regression will not improve results?


title: Python + R: How much R do I need to know to use random R libraries in python Using Rpy2

text: I'm just getting into data science from a qualitative background, and I want to be able to do the most stuff with least amount of overhead. (mostly machine learning and agent based modeling). I'm currently learning python, and would hate to switch languages in the middle of my coursera class mitx 6.00x. Also, I like that python is better for other stuff like scripting, data munging, NLP etc. However, I am eventually going to need to use some R packages (probably) it seems. Can I get away with learning only a minimum of R (how much and what exactly), if I learn pyhton + packages and use Rpy2 to access random R libraries, or would i still have to learn both languages fully? Thanks!


title: Matthews correlation coefficient - how much is a good classifier

text: I'm comparing performance of different classifiers on the data sets derived from financial markets, getting different accuracy and precision measures but Matthews correlation coefficient and Kappa statistic seldom exceeds 0.2. My data sets are quite big like 20000x170 so it should not be a problem with not enough of data. So my question is: what value of MCC and Kappa can be considered as a 'good classifier' on such data.


title: Determine individual consistency from mixed model?

text: To investigate the difference between three tasks I tested 30 individuals on each of the three tasks and ran the linear model as follows:               model <- lme(behaviour ~ task, random = ~ 1|subject, data=dat, method="ML")      With contrasts using `glht` to investigate if task 1 was different from task 2 and task 3. To furthermore investigate if individuals are consistent in their behaviour across the three different tasks I now simply test this by running three correlation tests such as `cor.test(behaviour.task1, behaviour.task2)`. However, I was wondering if it would not possible to get this information directly from the model, i.e. by looking at the intercept. I have found a couple papers that have done something similar like this by reporting if the intercept is significant, however they don't explain how they performed their analyses nor can I see how showing that the intercept is significantly different from 0 in a mixed model means that individuals were consistent. **Update based on Wolfgangs answer:**   First, to calculate the overall ICC of the model I ran               as.numeric(VarCorr(model)[1,1]) / (as.numeric(VarCorr(model)[1,1]) + model$sigma^2)      This shows the ICC = 0.336. To determine if the consistency of individuals is significant, I calculated if the ICC is different from 0 as follows:               model0 <- gls(behaviour ~ task, data=dat, method="ML")     anova(model0, model)      The output shows:               Model df      AIC      BIC    logLik   Test  L.Ratio p-value     model0     1  4 1211.739 1221.648 -601.8695                             model      2  5 1205.301 1217.688 -597.6505 1 vs 2 8.438081  0.0037      Thus the ICC is significantly larger than 0 ( _P_ < 0.005). Now, as I expect the tasks will differently affect the consistency of individuals, I run a new model as follows:               model2 <- gls(behaviour ~ task, correlation = corSymm(form = ~ 1 | subject), data=dat, method="ML")     summary(model2)      And get the following output:               Generalized least squares fit by maximum likelihood     Model: t.notcovered ~ session      Data: cot.sum      Subset: maleSub         AIC      BIC   logLik     1203.922 1221.263 -594.961          Correlation Structure: General      Formula: ~1 | rat       Parameter estimate(s):      Correlation:        1     2         2 0.233           3 0.453 0.584      thus showing that the correlation between task 1 and 2 is 0.233, between task 1 and 3 is 0.453 and between task 2 and 3 is 0.584. Now to see if this model provides a better fit than the model that assumes a single common correlation I ran the following:               anova(model, model2)            Model df      AIC      BIC    logLik   Test  L.Ratio p-value     model      1  5 1205.301 1217.688 -597.6505                             model2     2  7 1203.922 1221.264 -594.9610 1 vs 2 5.378838  0.0679      showing that the more complex model does not yield a significantly better fit, suggesting that the consistency does not differ significantly between the tasks, but as we know from the simpler model, individuals were overall consistent in their behaviour. _I added the above information with the idea that other people with a similar question may learn from this as well._


title: Mean of residuals in quantile regression are significantly differ from 0

text: Is it necessary to have mean of residuals which is equal to 0 in Quantile regression?


title: fitting a function

text: I have a data consists of three independent variables and dependent variable, I want to fit a model to these data, p=f(x1,x2,x3), how can I find the best fit for these data? it might be linear, quadratic,... could you please guide me to find the optimum way to find the model.


title: Stationarity - assumptions and examination

text: I am examining rodent captures on six permanent rodent trapping grids measuring 150 x 150 meters and consisting of 121 trap stations evenly spaced 15 meters apart. There are six such trapping grids on the study site that is < 1000 hectares in size. I would like to interpolate the capture data to create a Kriged surface of rodent activity. An assumption of interpolation is that the data are stationary. As Fortin & Dale (2005) state > stationarity is required for making inferences from a model that > characterizes the process of the spatial structure of data at locations that > are not sampled. From what I understand, a process can be described as stationary when its statistical properties (mean and variance) do not vary across space. But isn't variation across space why we conduct spatial analysis in the first place? Stationarity is very often introduced in the spatial/geostatistical analysis literature but, I have yet to find solid direction and information about   1. what scale, or for which types of studies, it is reasonable to assume your data are stationary,   2. how to examine and verify data are stationary, and lastly,   3. once quantified in some way just how much difference from one area to the next area qualifies your data as non-stationary? Thus far after reviewing the literature the concept and the examination of stationarity seems highly subjective, arbitrary and/or obfuscated. If anyone can provide some practical advice with this problem I'd greatly appreciate it!


title: Concordance test and corresponding plot

text: I have a list of genes identified using algorithm_A, then using algorithm_B and finally using algorithm_C; then a second list of genes identified using algorithm_A, then using algorithm_B and then algorithm_C, and so on for a total of 100 gene lists. To compare the overlap of the identified gene list using the three methods, I used the Cohen test for each gene lists. I would like to know first if it is the best choice and second how can I graphically plot the Cohen test results for the 100 gene lists.


title: How to extract/compute leverage and Cook's distances for linear mixed effects models

text: Does anyone know how to compute (or extract) leverage and Cook's distances for a `mer` class object (obtained through `lme4` package)? I'd like to plot these for a residuals analysis.


title: Gaussian Kernel function vs normal distribution function

text: I read from this link and thought that kernel density functions are used for solving the unrealistic normal distributions or specification errors. But when I read the description of kernel density in MATLAB, the estimate is based on a normal kernel function. Then what is the point? I mean it still assumes normality and the problem is unsolved. Why use kernel density instead of using normal distribution when it is still based on assumption?


title: A simple distribution problem

text: Suppose that $x \sim U(0, 255)$ and and $y \sim U(0, 255)$, and that $x$ and $y$ are independent. What is the distribution of $x - y$? If I were to draw 1,000 samples from this distribution, how many times should I expect to see each possible value? NB. These are all _discrete_ random variables, in case that wasn't clear... * * * **Update:** I've written a computer program to work out the expected count for each difference value. However, when I ask the computer to sum them all, it doesn't add up to 1,000. Clearly I've done _something_ wrong somewhere, but I'm not sure precisely what. I'm probably just being dumb; if somebody could walk me through this slowly, it would be helpful. The reason I'm doing this is that I want to do a chi-squared test to see if the data I'm getting actually does follow the expected distribution. Clearly this doesn't work if you calculate your expectations wrong. (!) In addition, I think I may have over-simplified the program description. What my program _actually_ does is obtain $x_1 \ldots x_{1000}$, and then compute $d_n = x_{n+1} - x_n$. Would that mean that the $d$ values are _not_ statistically independent? If so, presumably a chi-squared test won't work at all. (I suppose I could throw away all the even-numbered $d_n$ so I'm left with non-overlapping pairs...)


title: Does the average increase of a "re-roll" increase with the size of the dice?

text: I'm having trouble modeling the following scenario: > Roll two dice and take the best result. What is the proper way to model the average result from this? I understand intuitively that the average should go up, but I'm having trouble understanding exactly what plays into this increase. Furthermore, is the increase greater or simply proportional in dice of larger size (more sides)? For comparison's sake, consider a normal 6-sided die, and a 12-sided die.


title: Confusion related to decorrelation of PCA

text: I have this confusion related to PCA. PCA assumes that the variables given have corresponding latent variables which are uncorrelated. What PCA does it project the data into the eigen space which preserves maximum variance and decorrelates the dimensions. Why does PCA makes such assumption that the underlying latent variables are uncorrelated and all that PCA does is get those latent variables. Insights?


title: How to intuitively understand formula for estimate of pooled variance when testing  differences between group  means?

text: Suppose I want to compare the difference between means of samples selected from two populations (the treatment and control). Assume both groups have normally distributed observations. Then $$Z = \frac{(\bar{X}_{t}- \bar{X}_{c})-(\mu_{t}-\mu_{c})}{\sqrt{\left(\frac{\sigma^{2}_{t}}{n_t}+ \frac{\sigma^{2}_{c}}{n_c} \right)}}$$ Suppose that $\sigma_{t}^{2}$ and $\sigma_{c}^{2}$ are unknown but can be assumed equal to $\sigma^2$. Why is the pooled estimate $S_{p}^{2}$ for $\sigma^2$ equal to $$S_{p}^{2} = \frac{S_{t}^{2}(n_{t}-1)+ S_{c}^{2}(n_{c}-1)}{[n_t+n_c-2]}$$ where $S_{t}^2$ and $S_{c}^2$ are the sample estimates of the treatment and control groups. I know this has something to do with degrees of freedom. But I never could really "grok" its definition. In short, how do we get the pooled estimate and what are degrees of freedom intuitively?


title: Should I use Fisher's exact test for this?

text: I'm relatively new to statistics and I'm a little confused about which test I should be using for this problem. Let’s say I am collecting data about responses to a book under two different conditions (C1 and C2). In each condition, I am looking at several attributes about the book and asking for ratings using a 5 point Likert scale. The example tables below show the responses for the “enjoyability” attribute of the book, out of sample size of, say, 22.               C1     Rating:        1   2   3   4   5     Frequency:     2   1   2   9   8          C2     Rating:        1   2   3   4   5     Frequency:     10  6   3   3   0      I would like to compare the “enjoyability” ratings under each condition and test for significance (is the book significantly more enjoyable under C1?). Is Fisher’s exact test what I should be using? Is there another test which might be more suitable?


title: Using mle2() for age-period-cohort models

text: I've followed guidelines for comparing models in Chapter 6 of Bolker's _Ecological Models and Data in R_ , applying code used in this section to cancer count data. The models include parameters for age group (a) and birth cohort (c), and are as follows:               # Value of lambda is the same for all age groups and birth cohorts     poisfit.0 <- mle2(counts ~ dpois(lambda = a*c*pyrs), start=list(a=3,c=0.7),data=x)        Warning message: In dpois(x, lambda, log) : NaNs produced Coefficients: Estimate Std. Error z value Pr(z)   a 3.1932757 0.0022079 1446.28 < 2.2e-16   c 0.8429316 0.0083643 100.78 < 2.2e-16 -2 log L: 17550.32                # Use coefficients from simplest model to build more complex models       start.ab <- coef(poisfit.0)          # Model with varying parameter values for age groups and birth cohorts     poisfit.ab = mle2(counts ~ dpois(lambda = a*c*pyrs), start = start.ab, data = x, parameters=list(a~x\$ageGroup, c~x\$bCohort))          # Same as above, but without intercept     poisfit.ab1 = mle2(counts ~ dpois(lambda = a*c*pyrs), start = start.ab, data = x, parameters=list(a~x\$ageGroup-1, c~x$bCohort-1))      The data look like this, and go on for 110 rows (14 age groups, 16 birth cohorts):   ageGroup bCohort pyrs counts   1 13 1.4006 0   1 14 32.6925 1   1 15 49.7632 1   1 16 49.7059 0   2 12 1.2528 0   2 13 30.3879 3   2 14 47.7709 1   2 15 52.3699 0   2 16 55.7669 2   3 11 1.0257 0 A sample of coefficients for the model lacking intercept:               Estimate    Std.    Error   z   value   Pr(z)        a.x\$ageGroup1  0.101115    0.072032    1.4037  0.1603938        a.x\$ageGroup2  0.189811    0.078609    2.4146  0.0157517            a.x\$ageGroup3  0.335937    0.086470    3.8850  0.0001023        a.x\$ageGroup4  0.760921    0.117493    6.4763  9.400e-11        a.x\$ageGroup5  1.282271    0.141785    9.0438  <   2.2e-16     a.x\$ageGroup6  2.221877    0.180402    12.3162 <   2.2e-16     a.x\$ageGroup7  3.694786    0.227830    16.2173 <   2.2e-16     a.x\$ageGroup8  6.629915    0.313676    21.1362 <   2.2e-16     a.x\$ageGroup9  10.103253   0.400364    25.2351 <   2.2e-16     a.x\$ageGroup10 16.198128   0.537537    30.1340 <   2.2e-16     a.x\$ageGroup11 24.141987   0.686824    35.1502 <   2.2e-16     a.x\$ageGroup12 32.166317   0.811856    39.6207 <   2.2e-16     a.x\$ageGroup13 39.877367   0.941045    42.3756 <   2.2e-16     a.x\$ageGroup14 36.982237   1.169403    31.6249 <   2.2e-16     c.x\$bCohort1   0.058183    0.026073    2.2315  0.0256476        c.x\$bCohort2   0.135923    0.023366    5.8171  5.988e-09         There are two questions I have about this model:   1. The age group parameters blow out of proportion starting with age group 8. Why is this happening? Where did I go wrong?    2. The starting parameter values in the simple model poisfit.0 are both positive. They must both have the same sign in order for mle2() to accept them. If one is positive and the other negative, mle2() issues this error: Error in optim(par = c(-0.2, 3.4), fn = function (p) : initial value in 'vmmin' is not finite I don't know enough about what goes on under the hood to understand why this happens. I'd appreciate insight into these questions. Thank you, Susana


title: Friedman's Anova on qualitative data in SPSS

text: I have some qualitative data and having trouble preparing the data before I perform the analysis in SPSS. My understanding of Friedman's anova is that the data must be ranked from lowest to highest. My experimental design was repeated measures using a single variable (interaction technique) with three levels (a,b,c for example). The experiment involved participants ranking there agreement to 13 questions (from the IBM Usability Satisfaction Questionnaire). Effectil my SPSS spreadsheet has 39 columns (3 levels * 13 questions). Here participants had to rank there agreement on a scale between 1 - 5 with 1 being Strongly Agree and 5 being Strongly Disagree. Hence questions ranked with a value of 1 are more positive. Considering this I am confused how I enter the data into SPSS. If 1 is better should I invert this value and represent 1 as a ranking of 5 when entering into SPSS, and invert 5 -> 1 etc etc? I have tried both approaches and get different results. Which is correct? Any help would be greatly appreciated as I'm rather confused.


title: Power analysis with R

text: I have two questions regarding the R script. If I want to know the size of the sampling when I do a t-test with R, I can use this R script :               power.t.test(delta=1,sd=1.5,sig.level=0.05,power=0.8)      If I want to know the size of the sampling when I do the test of Pearson correlation coefficient, what is the R script? And for a chi-square test what is the R script?


title: Comparing sets of parameters extracted from different exponential fits to each other for significance

text: I take two simultaneous measurements from two different regions of a biological system, where both signals can be modelled by simple exponential rise and decay: $$ I(t)=(1-e^{-t/\tau_{rise}})(e^{-t/\tau_{decay}}) $$ First, by kernel matching to get initial coefficients, and later by nonlinear regression (`nlinfit` in MATLAB), I am able to get the fits for the two signals.The system is relatively noisy. Also, since every acquisition is different and each should be used independently for analysis, I cannot have means of repetitions. So, strictly, I only have 50 time points for each curve. My aim is to be able to compare these two fits. In short, I would like to be able to say whether the two are statistically different. One would perhaps, right this moment, suggest ANOVA or alike (I'm not the statistics expert at all!) to tell whether these two populations are different from each other. I request more: From the fits, I am able to extract their peak value, and when this peak occurs. I would like to have a measure on how reliable these extracted values are, considering the noise. For a noisier acquisition, my extracted values should have a lower confidence. My **first question** is: What is the best method to assess the confidence to my fits, and to the values that I can extract from these fits? I believe one cannot simply look at $R^2$ for nonlinear regression. **Second** and the **real question** is related to the fit comparison, or rather peak value comparison: So, I have the peak information. What matters most for me is if the peak of _fit A_ is greater and appears earlier than of _fit B_. In cases when I have strong signal and low noise, even if the peak values are not so different, by intuition, I think I should be able to compare the peaks for difference. In cases when noise is high (for any of the two signals or both) and signal is low, again by intuition, the values need to be separated more to be able to have sound a statement, otherwise I shouldn't be able to have a confident remark. How do I do this? I have a hunch that this is related to goodness of fit. Over a quick dinner conversation in a meeting, someone hinted me towards bootstrapping. However, I'm not sure whether I was concise enough to deliver my needs, nor have any idea on how to apply bootstrapping. Your help is greatly appreciated. I hope I was clear.


title: Modeling membership function given some survey data or empirical distribution

text: For example, I have a set of numbers (say 0 to 10) that are presented to 100 subjects. Each subject is asked whether the number is a small or a large number. The results are that 100 people think zero is a small number, 70 people think one is a small number, etc. Now I use a certain distribution, say, exponential to estimate the parameter of the distribution based on these data. I assume people give subjective opinions and that most of them think 0 is a small number so it may be assumed that the degree of truth of zero being a small number is 1. Finally, by multiplying the exponential probability density function (PDF) by a constant, I can scale the (PDF) such that the peak of it hits one and making it a membership function. Does this sound right?


title: $\sigma$-algebra intersection of infinite subsets

text: I am working out a book on Lebesgue measure by Bartle, and would like to see the steps that go into the construction of a proof for the following: > Show that any $\sigma$-algebra of subsets of $\mathbb{R}$ which contains all > open intervals also contains all closed intervals. That is, $[a,b]=\cap_{n=1}^{\infty}(a-\frac{1}{n},b+\frac{1}{n})$


title: How to decide if to do dimensionality reduction before clustering?

text: Is there any agreement on when to reduce data dimension before clustering in order to avoid curse of dimensionality? My intuition is that if I have say 1000 points and data dimension is 10 then it is OK to cluster. But if dimension is 50 then it is not OK because data points become sparse and hard to cluster (as a result expect to obtain "too much" clusters).


title: Dynamic factor analysis vs factor analysis on differences

text: I'm trying to wrap my head around dynamic factor analysis. So far, my understanding is that DFA is just factor analysis plus a time series model on the scores (the loadings remain fixed). However, in the cases that I've seen, the model on the scores is just a random walk with a diagonal correlation matrix. This seems identical to normal factor analysis applied to the differences. What am I missing? If you know of any good references to get me started, I'd appreciate them. I'd actually like to find something that allows the loadings to be slowly-varying; my context for thinking about that is West&Harrison-style DLMs, which hasn't got me far.


title: ANOVA mixed design in R shows incorrect Df

text: I'm running mixed design ANOVA using R. Somehow, I got incorrect values for the Df. There are two factors: `time` (within subjects), and `task` (between subjects). `Time` has four levels and `task` has three levels. I expected to find Df 3 for `time` and Df 2 for `task`. However, the ANOVA result I got showed 1 for the task Df under "Error:subject". Why does this happen? The dataset looks like this:               > data     subject  response  time  task       1        1.0076499 time1 task1       2        0.9939471 time1 task1       3        0.9981924 time1 task1     4        0.9960073 time1 task1     5        1.0064204 time1 task1     6        0.9990992 time1 task1     7        0.9959636 time1 task1     8        0.9922857 time1 task2     9        1.0007059 time1 task2     10       0.9952433 time1 task2     11       0.9976965 time1 task2     12       0.9957487 time1 task2     13       1.0013376 time1 task2     14       0.9894985 time1 task2     15       1.0049207 time1 task3     16       0.9904670 time1 task3     17       1.0015606 time1 task3     18       1.0086701 time1 task3     19       0.9970305 time1 task3     20       0.9847713 time1 task3     21       0.9990247 time1 task3     1        1.0047615 time2 task1     2        0.9930625 time2 task1     3        1.0449758 time2 task1     4        1.0332383 time2 task1     5        1.0117085 time2 task1     6        0.9849849 time2 task1     7        0.9486615 time2 task1     8        0.9795459 time2 task2     9        0.9803583 time2 task2      The script I wrote is here:               data = read.csv("anova_data.csv", header=T)     anova_mixed = aov(response~(task*time)+Error(subject/time)+task, data)     summary(anova_mixed)      Then, the result:               Error: subject          Df   Sum Sq  Mean Sq     task  1 0.001163 0.001163          Error: subject:time          Df   Sum Sq   Mean Sq     time  3 0.002001 0.0006671          Error: Within               Df   Sum Sq   Mean Sq F value  Pr(>F)                  task       2 0.005760 0.0028798   6.628 0.00234 **     time       3 0.000273 0.0000909   0.209 0.88976        task:time  6 0.001679 0.0002798   0.644 0.69468        Residuals 68 0.029544 0.0004345                        ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


title: Intra class correlation and experimental design

text: I am trying to understand how to apply the intraclass correlation coefficient. The situation we have is the following: we have $k$ judges who can assign numerical ratings to targets, and we want to test the reliability and consistency of these judges. One of the judges is considered to be an expert. So we pick $n$ targets, and each of the $k$ judges judges all $n$ targets twice, so we have $2nk$ measurements per experiment. Each judgement is somewhat expensive and takes a long time, so we do not want to make a judge do any more work than strictly necessary. The situation is as in case 3 as explained in section _Models for Reliability Studies_ in Shrout and Fleiss. Then we compare each judge to herself when rating the same target, so that's $n$ intra-ratings per judge, and we compare each judge to the expert judge, which is $2n$ inter-ratings per judge. There is no inter-rating for the expert. Then we compute an inter and an intra ICC using the following numpy code:               from numpy import isnan, nan, logical_not, logical_or          def exclude_nan(x,y):         """         Exclude NaN values if either entry in a pair of vectors has NaN         """         idx = logical_not(logical_or(isnan(x), isnan(y)))         x = x[idx]         y = y[idx]         n = len(x)         return [x, y, n]          def compute_icc(x, y):         """         This function computes the inter-class correlation (ICC) of the         two classes represented by the x and y numpy vectors.         """              if all(x == y):             return 1              [x, y, n] = exclude_nan(x,y)              ## Need at least 3 data points to compute this         if n < 3:             return nan              Sx = sum(x); Sy = sum(y);         Sxx = sum(x*x); Sxy = sum( (x+y)**2 )/2; Syy = sum(y*y)              fact = ((Sx + Sy)**2)/(n*2)         SS_tot = Sxx + Syy - fact         SS_among = Sxy - fact         SS_error = SS_tot - SS_among              MS_error = SS_error/n         MS_among = SS_among/(n-1)              ICC = (MS_among - MS_error) / (MS_among + MS_error)              return ICC      This corresponds to $\operatorname{ICC}(3,1)$ in Shrout & Fleiss linked above. The question now is, what $n$ should we pick? Making the same assumptions as in Case 3 of Shrout & Fleiss, what are the statistical properties of the ICC estimator?


title: What is the correct procedure to choose the lag when performing Johansen cointegration test?

text: When preforming Johansen Cointegration test for 2 time series (the simple case) you need to decide the lag you want to use. Doing the test for different lags return different results: for some lag levels the null hypothesis can be rejected but for others it can't. My question is what is the right method based on the input data to decide what lag I need to use when preforming the Johansen Test? p.s. I submitted this question to quant.stackexchange but some suggested it is better fit to this group.


title: Question about aggregated data and the effect on specific contributors

text: I have a set of data about the incidence of certain types of infections that occur in hospitals. The analysis has looked at the incidence trend for all aggregated results in a time series. This uses a GAM fitted spline to identify the incidence trend for aggregated data. Over the past year the incidence trend has increased sharply for the aggregated data set. The corresponding incidence trends for each participating hospital do not look much different. Is there a way to test a hypothesis that a specific individual hospital or hospitals drive the increase in the aggregated incidence?