title: Why the analysis of deviance of glm gives df=0 for one covariate?

text: I am applying a logistic regression on the effect of dose,age, PS, menopausal and pairID on the response variable. The data come from a case-control study where controls were matched/paired based on age, PS and menopausal. The pair ID is then recorded as PairID and is also a covariate in the model. The mode is:               mod2 <- glm(response~ dose+ Age + PS + menopausal + PairID, family = binomial(link = "logit"), data=dat)      The model output shows the estimated coefficient:               summary(mod2)          Deviance Residuals:       Min        1Q    Median        3Q       Max       -0.89847  -0.36858  -0.05885  -0.02615   3.00978            Coefficients: (2 not defined because of singularities)                      Estimate Std. Error z value Pr(>|z|)               (Intercept)     -32.73273   11.42767  -2.864 0.004179 **      MaxDose           0.07984    0.02287   3.491 0.000482 ***     Age               0.62496    0.27268   2.292 0.021910 *       PS1              -0.39054    0.89232  -0.438 0.661628         PS2              17.09581 1268.45566   0.013 0.989247         PS3              15.76341 1833.21516   0.009 0.993139         menopausalpost  -12.55850    5.28287  -2.377 0.017444 *       menopausalmale  -23.87464   10.51685  -2.270 0.023200 *       PairID2           5.14332    3.02191   1.702 0.088754 .       PairID3           3.60030    2.20299   1.634 0.102200         PairID4         -13.42706    6.47253  -2.074 0.038036 *       PairID5          -0.16041    1.50785  -0.106 0.915276         PairID6         -13.20065 1268.45572  -0.010 0.991697         PairID7         -19.62178 1833.21665  -0.011 0.991460         PairID8          -7.08125    3.55520  -1.992 0.046393 *       PairID9          -0.05779    1.84367  -0.031 0.974994         PairID10         -6.59761    3.44665  -1.914 0.055593 .       PairID11          3.13177    1.55304   2.017 0.043744 *       PairID12          9.04193    4.16930   2.169 0.030106 *       PairID13         15.00443    7.03528   2.133 0.032946 *       PairID14          1.00987    1.58307   0.638 0.523527         PairID15          5.45367    2.91475   1.871 0.061336 .       PairID16         -3.13399    2.18068  -1.437 0.150672         PairID17          6.85244    3.63025   1.888 0.059081 .       PairID18          5.60300    3.02755   1.851 0.064217 .       PairID19               NA         NA      NA       NA         PairID20        -11.53556    5.23541  -2.203 0.027569 *       PairID21        -11.09677    5.23554  -2.120 0.034048 *       PairID22          6.15099    3.16279   1.945 0.051799 .       PairID23               NA         NA      NA       NA         PairID24        -18.62866 1268.45749  -0.015 0.988283         PairID25          4.83989    1.81209   2.671 0.007565 **      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           (Dispersion parameter for binomial family taken to be 1)              Null deviance: 226.13  on 540  degrees of freedom     Residual deviance: 167.62  on 511  degrees of freedom     AIC: 227.62          Number of Fisher Scoring iterations: 17      From this output, seems that the covariate menopausal has an significant effect. However, If I apply analysis of deviance to test the significance of each covariate, covariate menopausal shows df=0. Does anyone know why?               drop1(mod2, test = "Chisq")          Single term deletions          Model:     collapsed ~ MaxDose + Age + PS + menopausal + PairID                Df Deviance    AIC    LRT  Pr(>Chi)         <none>          167.62 227.62                          MaxDose     1   213.08 271.08 45.462 1.557e-11 ***     Age         1   175.22 233.22  7.602  0.005829 **      PS          3   173.49 227.49  5.868  0.118227         menopausal  0   167.62 227.62  0.000                   PairID     22   186.31 202.31 18.684  0.664773


title: Kernel interpretation of Manhattan distance / Lp norms

text: Euclidean distance corresponds to the linear kernel similarity via: $$d_\text{Euclidean}(x,y)^2=\sum_i (x_i-y_i)^2=\sum_i x_i^2 + y_i^2 - 2 x_i y_i = \left<x,x\right>+\left<y,y\right>-2\left<x,y\right>$$ Is there any kernel-correspondence of Manhattan distance, and other $L_p$ norms (Minkowski metrics)? Or don't they yield a kernel function somehow, too? I'm trying to understand how these popular family of metrics relates to kernel similarity.


title: How to test whether a time series is significantly non-zero at any point in time?

text: I have some data that varies over time, and cannot be described by a simple model. In fact we don't care what the model is at all, we just want to know whether it is significantly non-zero at any point in time. Ideally, the result would also tell us at which time points the data is significantly different from 0. It is possible to do bootstrapping, but I am wondering if there is a more specific way to assess significance (a specific type of ANOVA, or t-test with corrections for multiple comparisons - but how to determine number of comparisons since the data is clearly correlated in time.). Thanks in advance - I'm sure this is a common problem, but after hours of reading I still can't find a clear answer. ![Example](http://i.stack.imgur.com/ncV1a.jpg)


title: Comparing very different models

text: I'm running an analysis and I've tried very different approaches, some frequentist and some Bayesian. Now I would like to compare the performance or fit of the models, and I was wondering how I can do it. If I only had Bayesian models, I could compare models by DIC, Bayes Factors, and/or posterior predictive checks. If I only had frequentist models, I could use AIC or BIC, for instance. Is leave one out cross validation suitable here?


title: What is the English name for a statistics term that I'm looking for?

text: I'm helping my father to translate a text from my own native Latvian language to English. Even though I'm a fair English speaker and have a bit of mathematics in my past, I never did have a firm grasp of statistics. :( The term I'm looking for has something to do with error calculation. A GPS device is measuring coordinates and then trying to estimate the error. In wikipedia I've come across two different likely terms - Mean squared error and Root mean square error. From what I can understand, they are not one and the same thing, but what each of them means is beyond me. A direct translation from Latvian language would be "average (mean?) squared error estimate". No idea what that means either, I'm afraid. Perhaps someone here has a better idea?


title: Using the log versus the level variable in Heckman second stage

text: Consider the Heckman selection model where my dependent variable is in log form for the second stage. I want to find the difference in computing derivative of unconditional expectation when using the log form versus the level form. \begin{equation} \frac{\partial E(\log y)}{\partial x}- \frac{\partial E(y)}{\partial x} \end{equation} I also want to know whether \begin{equation} \frac{\partial E(\log y)}{\partial x} \end{equation} is same as \begin{equation} \frac{\partial \log E(y)}{\partial x} \end{equation} A quick hint or reference to a paper would be appreciated.


title: Visualize movie/actor relationships

text: First off, I'm a programmer but my experience with true statistics ended at A-Level so I'm looking to all of you for help with a little side project I've been tinkering with. At home I use Plex Media Center to display all of my movies. I built an export tool for this to generate a HTML file containing information on your library so that others can view it online. After I made this tool I realised I now had access to a wealth of data about films and the actors in them. And this is where you guys (and gals) hopefully come in. I want to visualize the relationships between actors and movies somehow. Initially I just used a node graph library to map all actors who have been in more than one movie to all their movies and ended up with this: http://www.flickr.com/photos/dachande663/5574979625/ [section of a 5000x2500px image] The problem is, with anything more than 250 movies it just turns into a mess of spaghetti that's impossible to follow. I've looked into arc diagrams but think it would just be even more confusing. My question therefore is: how do I visualize this? Size isn't too much of an issue as I'd love to print this out on a large canvas and actually hang it up. Also, I'll eventually replace the text with images of the respective movies and actors. What I'm trying to avoid is having a million lines snaking everywhere. I've tried to find the most important movies and place them more centrally but at the moment that's more guess work than actual logic. Are there libraries that can do a better job of this, or even a better way of displaying the data (dropping actors as nodes and adding them as edge labels)? I'm currently using Dracula graph, which provides an okay-starting point but can change as needed. Any input will be much appreciated. Cheers.


title: Applying principal component analysis on variables of different dimensions?

text: Is it possible to perform principal component analysis on variables that have different number of rows of data? For instance, if there were just three variables - Var1, Var2, Var3 - where Var1 has 100 rows of data, Var2 has 50 rows, and Var3 has 20 rows. Can PCA be performed on such a data set? If PCA is possible, do I simply pad the empty ones with zero or is some sort of transformation needed? Had all variables been the same dimension (e.g. Var1, Var2, Var3 all has 100 rows of data), then there would be no issue. Unfortunately my data set is not so trivial. Please advise. Thanks!


title: Is the absolute value of distance covariance a metric?

text: I'm reasonably certain the absolute value of the distance covariance satisfies   1. d(x, y) ≥ 0 (non-negativity, or separation axiom)   2. d(x, y) = 0 if and only if x = y (identity of indiscernibles, or coincidence axiom)   3. d(x, y) = d(y, x) (symmetry) But I'm not sure about: 4) d(x, z) ≤ d(x, y) + d(y, z) (subadditivity / triangle inequality). I'm thinking in particular of a space of timeseries, but I don't think that matters too much.


title: Interpreting R's ur.df (Dickey-Fuller unit root test) results

text: I am running the following unit root test (Dickey-Fuller) on a time series using the `ur.df()` function in the `urca` package. The command is:               summary(ur.df(d.Aus, type = "drift", 6))      The output is:               ###############################################      # Augmented Dickey-Fuller Test Unit Root Test #      ###############################################           Test regression drift                Call:     lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)          Residuals:           Min        1Q    Median        3Q       Max      -0.266372 -0.036882 -0.002716  0.036644  0.230738           Coefficients:                  Estimate Std. Error t value Pr(>|t|)        (Intercept)  0.001114   0.003238   0.344  0.73089        z.lag.1     -0.010656   0.006080  -1.753  0.08031 .      z.diff.lag1  0.071471   0.044908   1.592  0.11214        z.diff.lag2  0.086806   0.044714   1.941  0.05279 .      z.diff.lag3  0.029537   0.044781   0.660  0.50983        z.diff.lag4  0.056348   0.044792   1.258  0.20899        z.diff.lag5  0.119487   0.044949   2.658  0.00811 **     z.diff.lag6 -0.082519   0.045237  -1.824  0.06874 .      ---     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1           Residual standard error: 0.06636 on 491 degrees of freedom     Multiple R-squared: 0.04211,    Adjusted R-squared: 0.02845      F-statistic: 3.083 on 7 and 491 DF,  p-value: 0.003445                Value of test-statistic is: -1.7525 1.6091           Critical values for test statistics:            1pct  5pct 10pct     tau2 -3.43 -2.86 -2.57     phi1  6.43  4.59  3.78        1. What do the significance codes (Signif. codes) mean? I noticed that some of them where written against: z.lag.1, z.diff.lag.2, z.diff.lag.3 (the "." significance code) and z.diff.lag.5 (the "**" significance code).   2. The output gives me two (2) values of test statistic: -1.7525 and 1.6091. I know that the ADF test statistic is the first one (i.e. -1.7525). What is the second one then?   3. Finally, in order to test the hypothesis for unit root at the 95% significance level, I need to compare my ADF test statistic (i.e. -1.7525) to a critical value, which I normally get from a table. The output here seems to give me the critical values through. However, the question is: which critical value between "tau2" and "phi1" should I use. Thank you for your response.


title: Non-negative coefficient restriction in stacked regression

text: I've been reading about stacked regression, as described, for example, here. It seems it's important that when you regress against the first predictors, you require that the regression coefficients are non-negative. How does one do that in R?


title: How does one cross-validate a model of binomially distributed data?

text: I'm ultimately looking to apply cross-validation to a much more complicated scenario, but for now consider the simple task of leave-one-out-cross- validating logistic regression. When you compare the left-out value against its respective model's predicted value, I presume you must convert the predicted value to the probability scale, but then do you also round the predicted value so that its either 0 or 1, or leave it unrounded? Any ideas on which is the better approach (and maybe why)? Example data and by-hand loocv of a logistic model:               myData = data.frame(      x = seq(-2,2,length.out = 10)     )     myData$y = rbinom(nrow(myData),1,plogis(myData$x))          myData$preds = NA     for(i in 1:nrow(myData)){         thisFit = glm(             data = myData[(1:nrow(myData)!=i),]             , formula = y~x             , family = binomial         )         myData$preds[i] = predict(thisFit,newdata=myData[(1:nrow(myData)==i),])     }          myData$predsProb = plogis(myData$preds)          #one approach:     myData$squaredError1 = (myData$y - myData$predsProb)^2     SSE1 = sum(myData$squaredError1)          #alternatively:     myData$squaredError2 = (myData$y - round(myData$predsProb))^2     SSE2 = sum(myData$squaredError2)          #which is a "better" estimate of future prediction error, SSE1 or SSE2? Why?


title: How to choose between ANOVA and ANCOVA in a designed experiment?

text: I am conducting an experiment which has the following:   * DV: Slice consumption (continuous or could be categorical)   * IV: Healthy message, unhealthy message, no message (control) (3 groups in which people are randomly assigned - categorical) This is a manipulated message about the healthiness of the slice. The following IV's could be considered as individual difference variables:   * Impulsivity (this could be categorical ie. high versus low or continuous and is measured by a scale)   * Sweet taste preference (this is also measured by a questionnaire which is 3 options to choose for each question)   * BMI - participants will be weighed measured accordingly (this could also be considered either categorical or continuous). As the groups will be randomly assigned to one of 3 groups I assume I am doing an ANOVA of some sort and would possibly use Factorial ANOVA as I am interested in which IV effects the DV the most but also the interactions between the IV's as research indicates that there are relationships between some combinations. But I am not completely sure of this due to needing to know whether it's best to have the IV's all categorical or continous or mixed. Or is ANCOVA a possibility or even regression but I am not sure about that given they are assigned to groups then categorised based on their answers to surveys. I hope this makes sense and look forward to hearing from someone about my query.


title: What's a good tool to create sankey diagrams?

text: How can I create sankey diagrams?


title: Most confusing statistical terms

text: We statisticians use many words in ways that are slightly different from the way everyone else uses them. This cause lots of problems when we teach or explain what we are doing. I'll start a list (and now I'll add some definitions, per comments):   * Power is the ability to correctly reject a false null hypothesis. Usually, this means correctly saying "something is happening".   * Bias - a statistic is biased if it is systematically different from the population parameter associated with it.   * Significance - results are statistically significant at some percent (often 5%) in the following situation: If the population which the sample comes from has a true effect of 0, a statistic at least as extreme as the one gotten from the sample would only occur 5% of the time.   * Interaction - Two independent variables interact if the relationship between the dependent variable and one independent variable is different at different levels of the other independent variable But there have to be many others!


title: Linear mixed model, negative information criteria values and Hessian matrix not positive definite

text: I am analyzing (in SPSS 19) the data from a field experiment using a Linear Mixed Model with Repeated Measures. When checking the model, I always obtain negative values in the information criteria table. I know the rule of thumb “smaller is better”, but how should these negative values be interpreted? As example, when I select a type of covariance “Unstructured” the Bayesian information criterion (BIC) reads $-642.868$ and with a “Compound symmetry” covariance it reads $-497.270$. So, which type of covariance, “Unstructured” or “Compound symmetry”, is better? The second question is that I receive the following warning message: "The Final Hessian matrix is not positive definite although all convergence criteria are satisfied. The MIXED procedure continues despite the warning. Validity of subsequent cannot be ascertained". May I ignore it or I have a real trouble? Is there some way to fix it?


title: What does accuracy mean in SPSS's automatic linear modeling output?

text: When running the Automatic Linear Modeling function of SPSS version 20, it will return an accuracy value. The higher the accuracy, the more predictive the model (I assume). However, what does this accuracy really indiciate? Is it the R2? Adjusted R2? Something else?


title: Multivariate regression for 100 binary output variables

text: I was asked to do some regression analysis on anonymised data, which I don't have yet. It will have a few categorical inputs and ~100 boolean outputs, something like this:               Category   Type         Size       Blood pressure     Y1   Y2    Y3  ....   Y99       -------------------------------------------------------------------------------     cat        persian      small      high               1    0     0   ....   1     cat        persian      big        low                0    1     1   ....   0     dog        wolfhound    big        normal             1    1     0   ....   0     ...        ...          ...        ...                      ......        duck       scoter       small      above normal       0    1     1   ....   1        (I just made up the meanings here). The number of samples in the real data set will be about a million, without missing values. The task is to predict the output, eg `P(Y86 = 1 | {cat, siamese, big, low})` and generally understand the relationship between inputs and outputs. My questions:   1. which approaches are worth trying, what are the pros/cons?   2. which `R` packages can be of help?


title: Performing "all-possible regressions" in R

text: I am trying to implement _all-possible regressions_ in order to select the best predictors of stock returns from an exhaustive list of potential economic/fundamental variables. My response variable _y_ (i.e. stock returns) is a panel of 3000 securities (cross-section), each having 384 observations (time-series). Would anyone please suggest the best way to handle this procedure in R, in the context of panel data? I came across the package `leaps`, but it addresses the case of _y_ as a response **vector** rather than a response **matrix**. Thank you very much,


title: Needle-in-a-haystack Regularized Regression

text: I'm in a setting where I am trying to model a continuous output variable given ~100 variables and ~100k datapoints. The signal-to-noise ratio is extremely low, and colinearity is very high. Among the variables are many "needle-in-a- haystack" binary-valued features. A "needle-in-a-haystack" binary-valued feature, $f$, is one where $Pr[f==1]$ is small (~0.01), but where it is important for our model to be unbiased when $f==1$. When I use OLS, the resultant model is properly unbiased when $f==1$. However, the model has undesirable characteristics stemming from noise and colinearity. When I try elastic-net regularization, the noise/colinearity problems go away. However, it appears that the act of regularizing causes the model to disregard bias for the needle-in-a-haystack features. Even when $f$ is selected by the model, the model generates unacceptably large residuals when $f==1$. I'm wondering how I can get the best of both worlds. I am currently training an elastic-net regularized model first, and then training a second OLS model to predict the residuals from the needle-in-haystack features. This seems to work decently, but I'm wondering if there is a more standard way.


title: Feature selection for "final" model when performing cross-validation in machine learning

text: I am getting a bit confused about feature selection and machine learning and I was wondering if you could help me out. I have a microarray dataset that is classified into two groups and has 1000s of features. My aim is to get a small number of genes (my features) (10-20) in a signature that I will in theory be able to apply to other datasets to optimally classify those samples. As I do not have that many samples (<100), I am not using a test and training set but using Leave-one-out cross-validation to help determine the robustness. I have read that one should perform feature selection for each split of the samples i.e.   1. Select one sample as the test set   2. On the remaining samples perform feature selection   3. Apply machine learning algorithm to remaining samples using the features selected   4. Test whether the test set is correctly classified   5. Go to 1. If you do this, you might get different genes each time, so how do you get your "final" optimal gene classifier? i.e. what is step 6. What I mean by optimal is the collection of genes that any further studies should use. For example, say I have a cancer/normal dataset and I want to find the top 10 genes that will classify the tumour type according to an SVM. I would like to know the set of genes plus SVM parameters that could be used in further experiments to see if it could be used as a diagnostic test.


title: Can you use the paired t-test in what seems like an unpaired t-test situation?

text: I came across a situation where a large retailer is testing a new retail program. It selected specific pilot retail stores to test this program. Next, it selected a control group of retail stores so that it would be as similar as possible as the pilot retail stores. It matched each pilot store to one control store based on numerous sociodemographic and business metrics. The selection criteria was to generate a P value from a student t test as high as possible, indicating a good match between two stores. This all made good sense. On the other hand, they told me they used the paired t test. This was a surprise. I thought this would be a classic unpaired t test situation.


title: What kind of test should I use for two different objects?

text: Say I'm trying to compare the speeds of 2 different cars in 5 different climates and roads. I've got the mean of each one and the standard deviation. Assuming all required assumptions are made, what test should I use to conclude if there's either enough or not enough evidence that 1 car is faster than the other?


title: How to deal with "not applicable" values in categorical variables

text: My situation:   * small sample size: 116   * binary outcome variable   * long list of explanatory variables: 44   * explanatory variables did not come from the top of my head; their choice was based on the literature.   * most cases in the sample and most variables have missing values. Approach to feature selection chosen: LASSO R's glmnet package won't let me run the glmnet routine, apparently due to the existence of missing values in my data set. It became clear from some of the answers to my previous question that I have do deal with more basic issues before considering imputation methods. I would like to add here new questions regarding that. On the the answer suggesting the coding as constant value and the creation of a new variable in order to deal with 'not applicable' values in continuous variables and the usage of group lasso:   * Would you say that if I use group LASSO, I would be able to use the approach suggested to continuous predictors also to categorical predictors? If so, I assume it would be equivalent to creating a new category - I am wary that this may introduce bias. If that is not advisable, what would be in the case of categorical variables?   * Does anyone know if R's glmnet package supports group LASSO? If not, would anyone suggest another one that does that in combination with logistic regression? Several options mentioning group LASSO can be found in CRAN repository, any suggestions of the most appropriate for my case? Maybe SGL? OBS: I am not a statistician.


title: Can you make population inferences from two different samples of the same population?

text: Say that 100 police-trained dogs were chosen and randomly divided into two groups. One group would be kept with trainers who don't expect the dogs to get too trained and the second group is kept with trainers that believe in the dogs and know they can get trained even further. After a month both groups were tested and the group of dogs with the trainers who believed in them were much better at doing what they were trained for. Could you make population inferences? Why or why not? I think that you could as you used the 100 and that's your population but I'm not sure. edit: The dogs were randomly chosen as well. Then the 100 dogs were randomly divided into groups of 50. So 2 groups in total.


title: How to do multivariate outlier detection in mixed data with category?

text: I have a data table where the entries are in the following format. The first column is `category`, which represent the product category. I have 5 such categories. Feature 1,2 and 3 are different features associated with each of these categories. These features may be either categorical or numerical. I am planning to apply multivariate outlier detection `mvoutlier` in R package to explore this data. I have following questions that I am trying to answer:   1. Should I consider each of these categories separately and perform the analyses independent of other categories?   2. How to handle the categorical values in this data? Thanks               category    f1  f2  f3     a1  1   33.4    333     a1  0   23  444     a1  0   30  333     a1  0   34  300     a2  1   56  600     a2  1   60  609     a2  1   64  630     a2  1   66  650     a3  0   99  900     a1  0   30  320     a3  0   99  1000     a1  0   30  340     a2  0   59  600


title: Linear combination of processes yielding higher order of integrarion

text: Is it possible to have an integrated random walk process by linearly combining finite number of random walks? Is it possible to have a random walk process by linearly combining I(0) processes? Thanks


title: Running simulated annealing optimization multiple times and averaging results?

text: I am using a simulated annealing algorithm to optimize a cost function. Given the fact that simulated annealing is a stochastic algorithm and will not give the same results each time you run it, I was wondering whether it is reasonable to run it multiple times and take an average of the results. If not, what would increase my chances of finding the global minimum? From what I have read adaptive simulated annealing is a good way. Thanks for your time.


title: Which test to choose when there is only one duplicate

text: I was wondering which test I can use to test if there is a significant difference in phenotype between cells given 6 different treatments against untreated cells. The problem is that I only have one counting from each treatment, but a large sample size for each treatment (n= 1000).


title: Is it valid to log-transform percentages?

text: My thesis is to examine the validity of the relationship between tourism represented by tourist exports (real growth) and economic growth using a dummy variable and force workers. I tested the normality of the variables and they are not normal, so I converted the variables using a logarithmic transformation in order to make them normal. My question is: **it is valid to log-transform a percentage variable?**


title: What are real strengths of association when using Pearson correlation?

text: I have located a table that breaks downs Pearson's correlation values into 3 categories. See below: ![enter image description here](http://i.stack.imgur.com/dJ3Ka.png) My question is this. One could have a large positive correlation of 1 as well as a perfect correlation. Is this table valid?


title: Feature selection: permutation test Vs deleting a variable

text: In feature selection for predictive models, it is usually applied a permutation test. In this test, all the values of one variable are randomly permuted and the prediction accuracy is extracted for each permutation. For example, if we have 3 variables (features), the ACC with all variables is 0.9. And, in the permutation test, we get ACC of 0.2, 0.1 and 0.9, respectively, when the first, second and third variables are permuted. Thus, the third variable can be dropped because it does not help for the prediction. However, if we delete the variable or change all the values to zero, instead of randomly permuting the values, do we get the same result in the permutation test?


title: Is it possible in R (or in general) to force regression coefficients to be a certain sign?

text: I'm working with some real world data and the regression models are yielding some counterintuitive results. Normally I trust the statistics but in reality some of these things can not be true. The main problem that I am seeing is that an increase in one variable is causing an increase in the response when, in fact in reality, they must be negatively correlated. Is there a way to force a specific sign for each of the regression coefficients? Any R code to do this would be appreciated as well. Thanks for any and all help!


title: pairing algorithm in R

text: I have two sets of elements _M_ and _N_ , and a scalar-valued distance/similarity function between one element from _M_ and one from _N_. The problem is to generate a set of pairs (one item from _M_ and one from _N_ ) which minimizes the sum of the distance function. Caveats   * Using an existing package in R would be ideal   * _M_ and/or _N_ may contain "extra" elements that don't have a partner; therefore _M_ and _N_ may be of unequal length   * We have a starting configuration which is reasonably close My first thought in addressing this is a monte carlo or genetic algorithm. So perhaps as a solution someone could show how to use an R genetic algorithm package for this problem. Also see this unanswered question.


title: Is there a test for omitted variable bias in OLS?

text: I am aware of the Ramsey Reset test which may detect nonlinear dependencies. However, if you just throw out one of the regression coefficients (merely linear dependencies), you may get a bias, depending on the correlations. This is obviously not detected by the Reset test. I did not find a test for this case, but this statement: "You cannot test for OVB except by including potential omitted variables". It is probably a reasonable statement, isn't it?


title: why feature scaling or weighting is important in surpervised learning?

text: I can understand feature scaling or weighting is important in unsurpervised learning case, because we want an good representation of "similarity". But why it is also important in surpervised learning case like in svm? Surpervised learning are supposed to learn those weights. Thanks!


title: Evaluation of a semi-supervised ranking model

text: I learned ranking model with graph-based semi-supervised approach, while labeled (just positive) and unlabeled (positive and negative) data is both used in training. With the model, all of the data could be ranked with a score after the training processing. Is there any way to evaluate the system just by the positive labeled data if I do not want to label the unlabeled data? Actually I would like to randomly split the positive labeled into two parts, and keep half as positive and half as unlabeled for the further evaluation. But I did not find a good metric.


title: Test of stationarity in Dickey and Fuller test

text: Why is that in DF (Dickey and Fuller) test, the test verifies that $Tst = 0$ or $Tst < 0$. Why not test $\Phi = 1$ or $\Phi < 1$? Need some guidance on this


title: What characteristics of the distribution of a test statistic can be inferred using a bootstrap?

text: **Setup:** Let $p_n(x) = \mathbb{P}(S_n \leq x)$ and let $p_n^*(x) = \mathbb{P}^*(S_n^* - S_n \leq x)$, where $S_n$ is some zero-mean test statistic that can validly be bootstrapped, eg a sample mean, and $S_n^*$ is a bootstrap re-sample of $S_n$. Note that $\mathbb{P}^*$ implies the probability is conditional on the underlying data. Assuming the conditions of some bootstrap theorem, eg stationary bootstrap, are fulfilled, we have the main bootstrap result: \begin{equation} \sup_{x \in \mathbb{R}} | p_n(x) - p_n^*(x)| \overset{\mathbb{P}}{\rightarrow} 0, \: \textrm{as} \: n \rightarrow \infty \end{equation} My question has two parts: **Question 1:** The main bootstrap result is usually used to assert that bootstrapped confidence intervals are asymptotically valid. Can we make this assertion because of Slutsky's theorem? If so, why isn't continuity an issue? Currently, my understanding is that the argument works as follows: 1) For some fixed point $x \in \mathbb{R}$, we use the above result to assert that $p_n^*(x) \overset{\mathbb{P}}{\rightarrow} p_n(x)$, as $n \rightarrow \infty$. 2)Use Slutsky's theorem to assert that $g(p_n^*(x)) \overset{\mathbb{P}}{\rightarrow} g(p_n(x))$, as $n \rightarrow \infty$, for some Borel function $g(a)$ that is continuous at $a$. 3) Choose $g(a)$ to be the quantile function, ie $g(F(x)) = \inf \\{ x \in \mathbb{R} : \lambda \leq F(x) \\}$ where $F(x)$ is a cdf and $\lambda \in (0, 1)$. The result then follows. The problem with this argument as it stands is that, for general cdf functions, the quantile function is only guaranteed to be left continuous. I'm assuming that we get around this problem by exploiting the fact that $p_n(x)$ and $p_n^*(x)$ are both converging to the Normal, which has a continuous cdf, but I've never seen this stated formally anywhere, so I am wondering if my reasoning is actually correct. **Question 2:** Having established that we can build asymptotically valid confidence intervals, I'm wondering what other characteristics of the distribution of $S_n$ we can consistently estimate? I know that most bootstrap theorems also include a result of the form $\mathbb{V}^*S_n^* \overset{\mathbb{P}}{\rightarrow} \mathbb{V} S_n$, as $n \rightarrow \infty$, so that takes care of the variance. What about other moments? For example, can we infer from the main bootstrap result that $\mathbb{E}^* S_n^{*4} \overset{\mathbb{P}}{\rightarrow} \mathbb{E} S_n^4$, as $n \rightarrow \infty$? As in the case of confidence intervals, I tried to prove this using Slutsky's theorem, but was not comfortable with the resulting expression: \begin{equation} \int x^4 dp_n^*(x) \overset{\mathbb{P}}{\rightarrow} \int x^4 dp_n(x), \: \mathrm{as} \: n \rightarrow \infty , \end{equation} since it is not clear to me that we are not accumulating lots of "small" errors in the integral. **EDIT:** I asked the question here since it is about the bootstrap. However, since it is also about probability theory, if users feel it would be more appropriate in stack exchange mathematics, please let me know.


title: Interpret interactions and logarithms in linear regression

text: I'm currently playing around with linear regression in R, and I've come up with a regression that fits data quite well. I'm just having some problems with interpreting the coeffiecents of my model. I know how to interpret log- log models in a simpler form, but when I have interactions I'm not quite sure how to interpret them. Here's my output from R:                   Call:         lm(formula = log(y) ~ log(x1) + x2 * log(x1) + x3 * log(x1) +          I(x3^2), data = Data)              Residuals:              Min       1Q   Median       3Q      Max          -0.56943 -0.12082  0.00012  0.11123  0.54579               Coefficients:                         Estimate   Std. Error t value          Pr(>|t|)             (Intercept) -2.393889950  0.545879641  -4.385 0.000025149470154 ***         log(x1)      0.497477722  0.056113496   8.866 0.000000000000009 ***         x2          -0.000264760  0.000055476  -4.773 0.000005220020368 ***         x3           0.041126987  0.017930934   2.294           0.02357 *           I(x3^2)     -0.000688879  0.000231778  -2.972           0.00358 **          log(x1):x2   0.000031580  0.000006691   4.720 0.000006494076511 ***         log(x1):x3   0.003145219  0.001277909   2.461           0.01528 *           ---         Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1               Residual standard error: 0.1932 on 119 degrees of freedom         Multiple R-squared: 0.9865,     Adjusted R-squared: 0.9859          F-statistic:  1454 on 6 and 119 DF,  p-value: < 2.2e-16       I've been googling for the past hour, but I can only find answers to some simpler models like the answer given here: Interpretation of log transformed predictor or http://www.ats.ucla.edu/stat/sas/faq/sas_interpret_log.htm I hope someone out there can help me with interpreting the interaction terms and the polynomial term in my model.


title: Multiple imputation for missing values

text: I would like to use imputation for replacing missing values in my data set. I have some constraints , for example I dont not want imputed variable `x1` be less than sum of my two other variables let say `x2 and x3`. I also want `x3` to be imputed by `0 or 14 or >= 14 and` and x2 to be imputed by `0 or 16 or >= 16`. I tried to define these constraints in spss for example for multiple imputation, but in spss I can only define maximum and minimum. Is there any way to define furthur constraint in spss or do you know any R package which let me define these constraints for imputation of missing value?                  x1 =c(21, 50, 31, 15, 36, 82, 14, 14, 19, 18, 16, 36, 583, NA,NA,NA, 50, 52, 26, 24)        x2 = c(0, NA, 18,0, 19, 0, NA, 0, 0, 0, 0, 0, 0,NA,NA, NA, 22, NA, 0, 0)        x3 = c(0, 0, 0, 0, 0, 54, 0 ,0, 0, 0, 0, 0, 0, NA, NA, NA, NA, 0, 0, 0)        dat=data.frame(x1=x1, x2=x2, x3=x3)        > dat            x1 x2 x3        1   21  0  0        2   50 NA  0        3   31 18  0        4   15  0  0        5   36 19  0        6   82  0 54        7   14 NA  0        8   14  0  0        9   19  0  0        10  18  0  0        11  16  0  0        12  36  0  0        13 583  0  0        14  NA NA NA        15  NA NA NA        16  NA NA NA        17  50 22 NA        18  52 NA  0        19  26  0  0        20  24  0  0


title: What to do with confounding variables?

text: I need to do an experiment. First let me describe present situation. The company that I work for is a cinema. It has a gaming section where people who are waiting for movies can pass time by playing games. People can pay only by using prepaid membership card. Unfortunately this gaming section is not generating enough sales. We are trying to find the cause(s). My hypothesis is if we accept cash as payment, sales will increase. My plan is to have experimental group and control group. The experimental group will accept cash payment, the control group doesn't. The sales of both groups are tallied before and after the experiment. The difficult thing about this is that I can't find a way to isolate the 'cash payment' factor from other factors:   * When the movie playing in the cinema is good, more people will come and sales will also increase   * Each cinema only has one gaming section, I can't split it into two sections (one accepts cash, the other doesn't)   * If several sites accept cash and several others don't, I don't think I can compare the results directly because the visitors are different, the number of gaming units are different I'm looking for suggestions to isolate this 'cash payment' variable, or maybe another approach altogether.


title: State space form of time varying AR(1)

text: I would like to implement the model proposed in Dynamic modeling of mean- reverting spreads (Kostas Triantafyllopoulos, Giovanni Montana). They propose to model a time serie Y_t with the following equations:               (1) Y_t = A_t + B_t * Y_(t-1) + e_t       (2) A_t = Phi1 * A_(t-1) + nu1_t       (3) B_t = Phi2 * B_(t-1) + nu2_t      That can be expressed in a state space form:               (1') Y_t = F_t * theta_t + e_t       (2') theta_t = Phi * theta_(t-1) + nu_t      with               F_t = (1, Y_t-1)       Phi = diag(Phi1, Phi2)      I would like to use R to perform a bayesian update of this model. I have studied the package 'dlm' but in the book Dynamic linear models with R, written by Giovanni Petris (author of the 'dlm' package), it is written (page 113) > "The matrix F_t of a DLM cannot depend on past values of the observations". However, it seams that it is the case in the model proposed above. Can someone understand this last sentence and eventually help me perform this implementation in R? Thank you Fred


title: What is a good use of the 'comment' function in R?

text: I just discovered the `comment` function in R. Example:               x <- matrix(1:12, 3,4)     comment(x) <- c("This is my very important data from experiment #0234",                     "Jun 5, 1998")     x     comment(x)      This is the first time I came by this function and was wondering what are common/useful uses of it. Since it is quite difficult to search "R comment" in google and find relevant results, I was hoping someone here might share with his experience.


title: Does efficiency imply unbiased and consistency?

text: If I can prove that for an estimator $\hat{k}( \theta)$ I can write: $$\frac{\partial l(X_1, \dots , X_n)}{\partial \theta} = a(n, \theta)(\hat{\theta} - \theta)$$ Am i sure that the estimator is unbiased? and consistent? NB:   * $l$: is the log likelihood   * $X_1$ is generated from a regular model   * $\hat{\theta}$ is the estimator for $\theta$   * $a(\cdot,\cdot)$ is a function of $n$ and $\theta$ (without any particular meaning i guess)


title: Uniform Distribution Test

text: I've got a data-set which I assume is uniformly distributed. Say I've got `N=20000` samples and a suspected `p=0.25`. This means that I would expect each option to show up roughly `5000` times. How do I calculate the following interval `[5000 - x, 5000 + x]` such that I can say with a certain confidence that the data-set is probably _NOT_ uniformly distributed since the number of times an option shows up falls outside of the interval? **EDIT** ABCDBCDADBCDA, BDCAADBCDADBA, ADCDBDACDBDAD, CDBDACDBDACDA, That's some sample data. A sample is one cookie string! Now I want for each position in that cookie string determine if a character there is too rare or too common at that position. So I count, for all samples, the number of A's on positon 0, the number of B's, C's and D's. Suppose I get a count of 5 A's on position 0 and I would expect a count of roughly 50 A's then the character A is too rare at position 0. That's what I want to do for each character position.


title: How to apply Shapiro test in R?

text: I'm pretty new to statistics and I need your help. I just installed the R software and I have no idea how to work with it. I have a small sample looking as follows:               Group A : 10, 12, 14, 19, 20, 23, 34, 41, 12, 13     Group B :  8, 12, 14, 15, 15, 16, 21, 36, 14, 19      I want to apply t-test but before that I would like to apply Shapiro test to know whether my sample comes from a population which has a normal distribution. I know there is a function shapiro.test() but how can I give my numbers as an input to this function? Can I simply enter shapiro.test(10,12,14,19,20,23,34,41,12,13, 8,12, 14,15,15,16,21,36,14,19)?


title: Error in estimating the size of a set?

text: Suppose we have a set A and a subset B. If we know |A|, then we can calculate |B| by finding the probability p that an element chosen uniformly at random from A belongs to B. Specifically |A|p=|B|. Suppose we generate n elements of A uniformly at random and use this data to estimate p (number of elements in B divided by n) and hence estimate |B|. > How reliable is this estimate? I.e. how can we compute the error? As a side question, is there a name for this technique? (it seems to be a mathematical version of the mark-and-recapture technique)


title: How to treat "Missing Value Analysis" test results (problems)

text: I have a problem with Missing Value Analysis. I am using SPSS version 20. I am trying to test whether missing values are at complete random. As I know in order to ensure missing values are completely random (MCAR: missing completely at random), Sig. value should be greater than 0.05. I have got two cases. First case, I have got Sig value less than 0.05 as illustrated by the picture ![EM \(Missing Value Analysis by SPSS\)](http://i.stack.imgur.com/lVIdV.png) Second case: I have got degree of freedom is zero and Sig is unspecified for other variable. (See, picture below) ![EM \(Missing Value Analysis by SPSS](http://i.stack.imgur.com/E2Bip.png) Please, what should I do? Thank you


title: Pearson Correlation for Clustered data

text: I have a sample that was drawn using cluster sampling. Let say I have variables x and y, and I would like to check the significance of correlation between those two variables in SPSS 20 taking into account cluster sampling. The problem is that SPSS "complex samples" dialog doesn't offer such an option. Would it be correct to define sampling plan (with "Complex Sampling" -> "Prepare for analysis") and use Complex Sample GLM to estimate simple model y = x + e with this sampling plan and then check if coefficient for x variable is significant?


title: Conservativeness of tests based on a discrete random variables

text: For discrete test statistics, the distribution of the corresponding $p$-value is discrete and stochastically larger than the uniform distribution. Hence the corresponding hypothesis test based on the p-value (reject if p-value is less than 0.05, for instance) is always conservative in the sense that the probability of making a Type I error will be smaller than 0.05. I know sometimes it is recommended to use the mid-pvalue. But I think there is no proof that using the mid-pvalue still controls the Type I error. Is there any other way to reduce the conservativeness? Can anyone who is familiar with this field give some hint or point out some existing literature on this?


title: When subject-based analysis is better than answer-based (and vice versa)?

text: In experimental psychology one of the generally accepted practices in data analysis is the aggregation on subject level. For example, several measurements of reaction time is gathered for each subject using different stimuli, then mean reaction time is computed for each subject. Those scores are further analyzed by ANOVA or other methods. My questions are:   1. When this procedure is more appropriate than a simple analysis of raw data (i.e. the analysis without any aggregation)? And when analysis of raw data is more appropriate?   2. Has anyone met the discussion of these question in literature? UPD: To give an example, suppose I'm interested in the effect of a drug on a speed of face detection. So I for each participant I have the data for drug and placebo condition (within subject design). I can analyze a) answers for each face for each subject for each condition, b) answers for each subject aggregated by condition, and, perhaps c) answers for each face aggregated by condition. When (a) is better than (b) and when it is otherwise?


title: Difference between different methods of outlier determination?

text: how to calculate Cooks distance. How is this method of determining outlier different from the quartile method.


title: Derive minimum positive rate of change for co2 data

text: I have CO2 (in parts per million) data of a closed room. The CO2 data is recorded along with timestamps. Typical difference between two samples is around five minutes. **My aim is to find occupancy of the room.** Theoretically, while the room is occupied the CO2 levels increase. And, when the room is unoccupied - CO2 levels decrease. So, I take 1st derivative of CO2 data with respect to time which gives me rate of change in CO2. At this stage, I judge the room to be occupied with positive rate of change (rate of change > 0) in CO2 and not occupied otherwise. I have some other CO2 data as well where I don't have validation head-counts. I use the same above approach to analyze this data and I believe there are many false positives in results (I have validation data of the room with actual Boolean head-counts) corresponding lower positive rate of change. The thought behind false positives being certainly there is because I know for sure the room was unoccupied for those particular duration. I want to derive a positive cutoff from training data (the one with validation set) so that I can apply the same approach on new data without any validation set being present. The room will be judged to be occupied when ( **rate of change** > **positive_cutoff** ). This can decrease the number of false positives in results. Have a look sample graphs: Smooth CO2 (Data with validation set) :![enter image description here](http://i.stack.imgur.com/fRy6c.jpg) Results after rate of change analysis: ![enter image description here](http://i.stack.imgur.com/hXeVS.jpg) CO2 data without validation set : ![enter image description here](http://i.stack.imgur.com/4P9CV.jpg) Results after rate of change analysis on above data : ![enter image description here](http://i.stack.imgur.com/etZLE.jpg) Now, I realize that for all actual occupancies, the rate of change is significantly higher than zero. Hence, the idea of deriving a positive threshold from training data (co2 values + actual headcounts) and use that threshold for predicting occupancy of room in test data. How should I approach the problem from here on?


title: Can the use of dummy variables reduce measurement error?

text: If the continuous variables are measured with error, can the use of dummy variables mitigate the problem? For instance, IQ measures intelligence with error. So will using a dummy of high, medium, low IQ mitigate the measurement error problem? Thanks!


title: Two independent non-random samples

text: Is there any statistical tool available to test the hypothesis regarding differences in means of two independent **non-randomly** selected samples. Is **Mann-Whitney-Wilcoxon test** suitable for such conditions?


title: ICA unreasonable result

text: I am using the following code for Independent Component Analysis of some EEG data:               def ICA(data):         import mdp              data0 = data.mean(0)         FastICA = mdp.nodes.FastICANode()         ics = FastICA(data - data0)         #W = FastICA.get_projmatrix()         A = FastICA.get_recmatrix()              return ics, A, data0      This function generally works and I am satisfied with the results, but for some cases, it produces a zero value for some components. Does this mean anything special? I am pretty sure that these data are similar to other normal ones (at least in shape) and there isn't too much noise or any extraordinary difference (at least in the first looking). I should mention hat i even do the whitening procedure on the data, but it doesnt help. Basically what would zero component mean? Any idea about the initial constraints that the data should have? Do I need to apply first any other kind of procedure on the data?


title: How to explain the value of y-axis in Partial Dependency plot

text: sample numbers = n   number of classes = 2   The function is $$F(x) = \sum_{i = 1}^nf(x,x_{ic}) $$ where $$f(x,x_{ic}) = \log p_1(x,x_{ic}) - \frac{(\log p_1( x,x_{ic} )+\log p_2( x,x_{ic}))}{2}$$ ![enter image description here](http://i.stack.imgur.com/ema8x.png) In the image, as x = 0.25, y is about 0.78.   How to I explain the value 0.78   Thanks in advance for your help.


title: expressing this probability distribution over different variables

text: I have a likelihood function as follows: $$ P(y|x,w, \phi) = \frac{\phi}{2\pi} \exp ^{-0.5 (y-t(x, w)'\phi (y-t(x,w)) } $$ Here $y$ and $x$ are two observed values. $\phi$ is also some given parameter and $t$ is a non-linear transformation parameterised by $w$. $w$ is a two dimensional position vector. $t$ is basically a spatial transformation of $x$. So as usual, the likelihood is the data term and is telling me the probability of observing $y$ for a given $x$, $w$ and $phi$. So the cost function becomes the sum of square difference between $y$ and $t(x, w)$ under the Gaussian noise model where the noise precision is given by $\phi$. So, I am trying to use Expectation Propagation (EP) to get a Gaussian approximation to the posterior on $w$. Let us call this distribution $q(w)$. EP approximates each of the likelihood terms with factors where each of the factors is a Gaussian. So, in my case it is (ignoring normalisation for the moment): $$ q(w) = \prod_i f_i(w_i) \approx \prod_i P(y|x,w, \phi) $$ So, the likelihood factorises as a product over the terms and each of these ith terms is approximated by the corresponding factor $f_i(w_i)$. Each of these factor terms is a low dimensional Gaussian on $w_i$. Now, EP is an iterative algorithm where at each iteration we need to do the following computation: $$ \prod_{i\neq j} f_i(w_i) P(y_i|x_i,w_i, \phi) $$ So, I need to do this multiplication where this factor $f_i$ directly depends on $w_i$ with the exact likelihood i.e. $P(y_i|x_i,w_i, \phi)$. The $f_i$ terms are 2-dimensional Gaussians and so, I need to transform this likelihood so that it is some 2-dimensional distribution over $w_i$ (I realise it will not be a Gaussian, of course). $t$ is a non-linear transformation that is parameterised using first order spline i.e. simple linear interpolation along the two dimensions. So, t(x, w) can be written as: $$ t(x, w) = \sum_{d} x_d \; B_{d}(w) $$ Note that $w$ is two dimensional and $B_d$ is the dth first degree B-spline coefficient. So, we are interpolating the observed $x$ at spatial locations $w_d$ using linear interpolation.


title: Find distribution of the data from Q-Q plot

text: With 10 thousands, Monte Carlo simulation, I have generated the Q-Q plot ![enter image description here](http://i.stack.imgur.com/1DPbO.png) Is it possible to infer what distribution follows my sample? I'm new to Q-Q plot.   As far as I understand, the sample is not normal distributed because the dot are not on the line.   But, then, which distribution behave like that? This is related to this question. **Edit**   To give an idea of the data, here is the empirical cumulative distribution function of the sample ![enter image description here](http://i.stack.imgur.com/wI4FD.png)


title: Is it statistically sound to use Lasso for variable selection even when $n\gg p$?

text: I have a classification task where $n\gg p$ (like 440000 vs. 23). I want to use Lasso (`glmnet` in R) to select the variables first, then use techniques like random forest or boosting for the actual classification. I read this from _Statistics for High-Dimensional Data: Methods, Theory and Applications_: ![enter image description here](http://i.stack.imgur.com/09CU7.png) I am not sure whether Lasso will still do variable selection when $n\gg p$. Any suggestions? Thank you.


title: Why do we assign fractional ranks in rank tests?

text: Both Mann-Whitney U test and Wilcoxon signed rank test follow this procedure:   1. sort the data   2. assign ranks, **Ties receive a rank equal to the average of the ranks they span**   3. compute the statistic using a **sum the ranks**   4. compare statistic against critical values and make verdict For me the 3rd step seems like an **unnecessary complication**. If we intend to sum the ranks anyway, there is no difference if we sum 2 and 3 or 2.5 and 2.5. Am I missing something?


title: How to determine good $\alpha$ for elastic net?

text: Using the paramater names from What is lambda in an elastic net model (penalized regression)? I wonder about best practices or suggestions for how to determine $\alpha$, i.e. the paramater determining the weighting of the Lasso and the Ridge. My guess is that one does a 1 dimensional grid search (is there another term for that btw?) starting with 0.5 and then going of towards 0 and 1. My feeling here is that we want to not go off with the same step size all the time. I mean if one of these values is much bigger than the other perhaps we want a value really close to 0 or 1 to compensate? But how close to 0 and 1 do we test? What is the standard approach here? Is there such a thing?


title: Question about Good Turing Discounting

text: I've got a question about Good Turing discounting. I understand the how and the why, but I'm having trouble wrapping my head around it. Say we're discounting the probability of each n-gram that occurs less than 5 times, so for every n-gram we use GT-discounting if it occurs less than 5 times, and the probability mass given to unseen events is `N_1 / N`. If the training set follows a Zipfian distribution and has a large number of n-grams occur once, `N_1/N` is going to be a fairly large number, and the discounted probability mass of the other n-grams plus the mass for the unseen events will be greater than one. Is this correct, or am I missing something? What is the probability of a single unseen event? It can't be the `N_1/N`, that's far too large of a probability for a single unseen event.


title: Confidence interval for predicted probabilities

text: I fitted the following multinomial regression:               library(car)     p1<-c(1,2,3,4,3,4,3,4,3,2,1,2,1,2,1,2,3,4,3,2,3,4,3,2,2,2,3,4,3,3,4,3,4)          d1<-c(1,2,3,4,3,4,3,4,3,2,1,2,1,2,1,2,3,4,3,2,3,4,3,2,1,2,3,4,3,2,2,2,1)          d1<-as.ordered(d1)     test<-multinom(p1~d1)     predi<-expand.grid(d1=c("1","2","3","4"))          pre<-predict(test,predi,type="probs")      The output is a table of the predicted probabilities for every coefficient. I can also order the results for the confidence interval of the coefficents with:               confint(test)      My question is: is it possible to get the results for the confidence interval for the predicted probabilities? It means for every amount in the "pre" output! PS: I found a similar question here in ["plotting confidence intervals"][1]Plotting confidence intervals for the predicted probabilities from a logistic regression The main answer is perfect for my question, but I do not know how to combine with multinomial regression. I hope you understand my bad english :) Thank you for your help


title: Weighted Least Squares for NMF on concatenated images of different sizes

text: I have a large set of images decomposed into 5 spatial frequency bands, with each band downsampled appropriately. So for each image I have 64x64, 32x32, 16x16, 8x8 and 4x4 for each spatial frequency. Now I want to do NMF on the combined image set, considering all spatial frequencies together. The first option would be to upsample all spatial frequencies to 64x64, but then I have ~20,000 features and I do not have a machine withe enough RAM to do that. Doing NMF on the 5456 downsampled vectors work, but obviously its biased very much towards to the high spatial frequency. I saw NMF uses alternating least squares algorithm, so I thought if I could weight the components of the different spatial frequencies appropriately, this could correct the imbalance in number of pixels and result in equal overall weight to each SF component. The key part of the ALS algorithm consists of the following lines:               h = max(0, w0\a);     w = max(0, a/h);      So I modified this to:               h = max(0, lscov(w0, a, weight));     w = max(0, a/h);      Based on my first thought about how to weight to normalise the number of elements I thought I would weight each SF (high to low) by [1 4 16 64 256] respectively (because each lower SF has 4x less elements). But this produced solutions weighted far too much to the low SF. Instead weights of [1 2 4 8 16] give solutions that look great - but obviously I have to come to this in a completely ad hoc way so I am a bit worried if it is justifiable. So the questions:   * Is my modification to the NMF ALS algorithm valid? I couldn't find a way to add weights to the second iteration for w - because it is along the other axis so needs weights along the number of different images direction.   * Why does 2x weights seem to be correct rather than 4x weights (is there a square somewhere that I am missing - I thought weights were applied like XWX so would not be squared)?


title: Which confidence intervals should be used?

text: I would like to plot 95% confidence intervals for my data. It seems there are multiple options for calculating them, e.g.: (examples shown using R code)               library(psych)          # Make some data     set.seed(1)     N <- 10     data <- data.frame(A.1=rnorm(N, 1, 1), A.2=rnorm(N,2,1), B.1=rnorm(N,2,1),     B.2=rnorm(N, 1, 1))     data          desc <- describe(data)     alpha <- 0.05               # Either:     # (1) using 1.96 (the quantile from the normal distribution)     ci <- desc$se * 1.96         # or ci <- desc$se * qnorm(1-alpha/2)          # (2) using the quantile from the t distribution     ci <- desc$se * qt(1-alpha/2, desc$n - 1)          # Note that 1 and 2 would be used like this:     # upper <- desc$mean + ci         # lower <- desc$mean - ci          # (3) using R's confint function:     ci <- confint(lm(values ~ 0 + ind, stack(data)))          # (4) using the trimmed mean:     library(WRS)     ci <- apply(data, 2, trimci)          # (5) various bootstrap methods, e.g.:     ci <- apply(data, 2, trimcibt)      I am a bit lost as to which method I should be using. I understand bootstrap methods should be more accurate because they do not assume the distribution is normal, but I rarely see them applied in practice. Does anyone have any guidelines or advice?


title: Exploratory factor analysis and latent variable measuring actual use

text: I have a technology acceptance model where I assume that 10 independent variables (each latent variable is composed of three items, based on a 1-5 likert scale) should influence the actual use (dependent variable) of a contactless payment system. Only the dependent variable is composed of two items, the first one measure the frequency on a 1-5 likert scale (1.never - 2. few times per year - 3. few times per month - 4. few times per week - 5. every day) while the second one measure the number of uses (so it is an open ended question, linked to the previous item). I need to conduct a SEM analysis so I should start with an exploratory factor analysis. My sample is composed of 350 cases and I've found that for the two items measuring actual use the distributions are strictly not normal. In particular, only 48% declared to use the system and usage ranged from 1 to 730 per user over a period of one year. The first item has a Skewness of 1,535 and a Kurtosis of 2,180 and the second item has a Skewness of 6,934 and a Kurtosis of 54,817. So I would like to ask if it is possible to conduct the analysis. Should I first standardize all the variables, because of the different scales of measure, and then transform both the two actual use variables to normalize the distributions? Thank you so much, it is for my thesis and I have big big troubles


title: Linear regression with fat-tailed errors

text: I'm testing a linear model that explains stock returns with some contemporaneous factors; the model is assumed to satisfy OLS assumptions except that the errors (i.e., unexplained stock returns) have fat tails: let's say they are described by an alpha stable or Student's t distribution. How should I estimate the model? I assume I should transform the data, or use a generalized linear model? I'm willing to be liberal with the distribution assumptions, as long as the calibration of the distribution and the estimation itself are not too complex.


title: Why would R return NA as a lm() coefficient?

text: I am fitting an lm() model to a data set that include indicators for the financial quarter (Q1, Q2, Q3, making Q4 a default). Using lm(Y~., data = data) I get a NA as the coefficient for Q3, and a warning that one variable was exclude because of singularities. Do I need to add a Q4 column?


title: How to analyze mood data over time?

text: I've collected data for individuals on a team about their daily mood. Each day, individuals rate their mood on a scale that is assigned the following values: 0 (bad), 5 (so-so), 10 (good). Using the data from the individuals, I can come up with an average mood for the team for that day. I have several months worth of data for around 10 teams, with team size varying from 3 to 12 people. I would like to be able to do some analysis with this data to answer questions like "Is this team in a better or worse mood than normal?" instead of just being able to report a numeric value. What sort of statistical models or calculations should I be looking into? Here is what I have implemented so far:   * Use the past N days of data to compute a "baseline" average and standard deviation. (N is ~30-90)   * Use the past M days of data to get capture the "current feelings". (M is ~3-5)   * Compute the z-score and use different "bands" to categorize (more than 10% worse than the mean => team is not doing well, more than 10% better => team is doing great, inbetween => team is normal) While doing my own research, I've run into things like "time-series analysis" and (seemingly) advanced techniques like ARIMA. Would these approaches be worthwhile to investigate further? Is my method valid and reasonable? I am looking for some guidance. My background is in software engineering; my last statistics class was in high school so if it is more complicated than standard deviation, please try to explain it like I'm still in the 11th grade :-)


title: What are some common tools, initial approaches to data in a prediction problem when facing too many predictors?

text: If one is given several hundred features (of both categorical and continuous type) what are some approaches to determining which features to keep or even drop? Data as such is difficult to visualize and other than PCA (which I think qualifies because it throws away information deemed not as important by variance) I'm not sure what other common data pruning/cleanup methods I can look into?


title: What is the difference between Visualization and BI Visualization?

text: I am currently doing visualization on my project and I need to know whether the visualization I have done is BI visualization or not. So is there any difference between Visualization and BI Visualization? If yes, then what are the differences?


title: Ordinal logistic regression to find support for inverse proportionality of a continuous predictor on an ordinal dependent variable

text: I want find statistical support that dependence is inversely proportional to power. To do so, I have   * ~260 cases,   * with four questions about dependence, and   * one question about power The questions about dependence are on a continuous scale, whereas the question about power only allows three ordered answers (I am powerful, equilibrium, the other one is powerful). To support the inverse proportionality of power and dependence (in this application), is the right way to do a ordinal logistic regression? I have plotted grouped error bars for the data (`error.bars.by()`) and they show the case I want to prove quite clearly; however I suppose I need the right figures on top of that as well. Thanks for any suggestions and advice. * * * I had a look at some threads on Cross validated (1,2,3) about similar questions, and as far as I understand there is not a clear answer to how the above described _mutual_ inverse correlation/proportionality could be tested. Correct me if I'm wrong. * * * Thanks for answers and comments so far. As far as I understand ordinal logistic regression helps me to find support for a relation of power with dependence, if I use dependence as predictor and power as dependent variable.


title: Validating logistic regression - formulas needed for simulation

text: This is related to a data classification problem having a Boolean output variable. Summary: Once I perform the ML task using Logistic regression, I get the required coefficients. I use the multivariate model $\theta_0 + \theta_1X_1 + .. = Y$. Having obtained the $\theta$'s, am looking for some straightforward formulas that I can use to get the following:   1. For the coefficients:      * Standard Error: I think it is $(X^TX)^{-1} * (X^TY)$. Please correct me if wrong     * Wald Stats: Is it square of $\text{Coefficient} / \text{StdError}$?     * 95 %CI: Is it $(\text{lower, upper}) = (\text{coefficient}_i - 1.96 SE_i,~ \text{coefficient}_i + 1.96 SE_i)$?     * p-value and significance   2. For the whole model      * Null model -2 Log likelihood     * Full model -2 Log Likelihood     * Chi-squared     * DF     * Significance level   3. ROC Curve      * Area under the curve (AUC) I am new to ML and trying out C# based simulation, so any vector-based formulas would be highly appreciated.


title: What is the minimum historical data/sample data required for a time series forecasting analysis?

text: Are there any statistical power analysis/sample size deteminations methods for time series data analysis/forecasting? For example if I have time series of 30 data points, how can I with _confidence_ use a particular statistical methods like exponential smoothing or arima for predict the future ? I have seen in some textbooks that have a cursary mention on historical data points required for ARIMA would be 50 or 60. But I have not encountered a formal approach on how much history is required for a a particular time series forecasting method. I did a thorough search on major time series textbooks and the internet, I'm unable to find any literature on this topic. Any guidance would be helpful.


title: Best way to simply store data for statistical analysis in R

text: I have been using text files to store my data for R without any problem for some time now. But for a recent project the files' sizes are getting too big for raw text files to handle. What's the best simple alternative?


title: How to determine a good simulation time?

text: So I'm doing network simulations, and I'd like to know how long does each simulation run needs to be. My network is quite simple: its composed of **multiple** M/M/1/H queues (Markovian processes + finite waiting queues):   * clients arrive at some node(s), are put into waiting queues   * each node process one client at a time   * a processed client is diriged to another node, or goes out of the system   * all random processes follow a Poisson distribution The finite waiting queues prevents it of being a Poisson process: the output of a node is no longer Poisson, so neither will be the inputs of the nodes after it. I thought I could try to plot the variance and mean of the waiting time (inside the queues) and see when it becomes "stable"? Would that be a good solution? How else could I do that? Many thanks,


title: Moderated regression using PROCESS (Hayes) with categorical predictors (dummy coding)

text: Has anyone tried to run a moderated regression using a dummy code with `PROCESS`? With moderated regression I thought you had to entered all the independent variables together, but process allows only one independent variable at a time. Thanks


title: Prediction using multinomial logit regression in R

text: I've got three variables, a factored (c) and two ordinal independent (a,b). Every variable has five categories (1,2,3,4,5). Thus, I fitted a multinomial logit regression (testus, see below) with the car package. Now I would like to predict the probabilities.               > a<-sample(5,100,TRUE)     > b<-sample(5,100,TRUE)     > c<-sample(5,100,TRUE)     > required(car)     > a<-as.ordered(a)     > b<-as.ordered(b)     > c<-as.factor(c)     > testus<-multinom(c~a+b)     > predictors <-     > expand.grid(b=c("1","2","3","4","5"),a=c("1","2","3","4","5"))     > p.fit <- predict(testus, predictors, type='probs')     > probabilities<-data.frame(predictors,p.fit)      Now I got the predicted probabilities for a under b and c.               `head(probabilities)       b a         X1         X2         X3         X4        X5     1 1 1 0.10609054 0.22599152 0.20107167 0.21953158 0.2473147     2 2 1 0.20886614 0.27207108 0.08613633 0.18276394 0.2501625     3 3 1 0.17041268 0.24995975 0.16234240 0.13111518 0.2861700     4 4 1 0.23704078 0.21179521 0.08493274 0.03135092 0.4348804     5 5 1 0.09494071 0.09659144 0.24162612 0.21812449 0.3487172     6 1 2 0.14059489 0.17793438 0.29272452 0.26104833 0.1276979`      The first two cols shows the categories of the independent variables a and b. the next five columns show the conditional probabilities (p.e. P(c=1|b==1&&a==1)=0,10609). But now I would like to know only the predicted probabilities for c under a or the predicted probabilities for c under b. Is this possible?


title: Numerical example to understand Expectation-Maximization

text: I am trying to get a good grasp on the EM algorithm, to be able to implement and use it. I spent a full day reading the theory and a paper where EM is used to track an aircraft using the position information coming from a radar. Honestly, I don't think I fully understand the underlying idea. Can someone point me to a numerical example showing a few iterations (3-4) of the EM for a simpler problem (like estimating the parameters of a Gaussian distribution or a sequence of a sinusoidal series or fitting a line). Even if someone can point me to a piece of code (with synthetic data), I can try to step through the code.


title: Plotting relationship between three variables in Matlab

text: I am new to Matlab and data visualization in general. I am plotting an A vs. B relation based on 100 values of (A,B) using `plot(A,B)` which gives me a nice 2-D line plot. Now I want to run 1,000 trials of this and see how the A vs. B relation varies across the trials. (A has a random component and B is dependent on A so the trials will yield different results.) What is the best way in matlab to visualize/analyze this? (I presume drawing 1,000 2-D line plots and looking at them isn't an option.)


title: Transforming distance matrix to univariate data

text: Is there a valid mathematical or statistical approach to converting site by site distance objects/matrices to a univariate site vector or column in a dataframe? I have a Euclidean distance matrix for distances between all sites, a cost distance matrix for distances between all sites, and an environmental distance for distances between all sites. I need to find some way of transforming each of these distance matrices to a single value per site, i.e. one euclidean distance value per site that encompasses the variation of all sites. I have used the first axis of variation in Principal Coordinates Analysis (PCoA) and I'm not sure if I like this approach. Several literature searches have not revealed any potential solutions, hence I am putting it to you. Any ideas would be immensely helpful.


title: Looking for estimates for my data using cumulative beta distribution

text: I'm trying to find the estimates for $alpha$ and $beta$ of the beta distribution (red curve) that will fit the black curve (empirical) better. ![trying to find parameters to make the red curve fit the black curve better](http://i.stack.imgur.com/DS11t.jpg) I have a function that finds "non-parametric" inverse cumulative function of my data (`boot.mean`). I would like to find a distribution that would fit. So far, I think cumulative beta distribution (alpha = beta = 2) could somehow be used. I would love to see a distribution that fits the bill better, though...               # raw data produced by a function (inverse cumulative distribution)     boot.mean <- c(37.021, 35.051, 29.091, 27.094, 22.058, 18.994, 16.944, 12.897, 7.903, 4.926, 3.939, 1.94, 1.94, 0.968)          #"fidge" (not sensu Climategate) boot.mean for comparison to qbeta     scaled <- 1 - (boot.mean/boot.mean[1])     scaled[1] <- 0.01 #dreader zero be gone     range(scaled)      [1] 0.0100000 0.9738527          # this is the theoretical curve     x <- seq(0, 1, length = 100)     y <- qbeta(x, shape1 = 2, shape2 = 2)          # all along the x axis     x.axis <- seq(from = 0, to = 1, length = length(scaled))          # plot empirical and the theoretical values     plot(x.axis, scaled, type = "l")     lines(y, x, col = "red")          # I'm just an x-con trying to fit a distribution to my data     (beta.fitted <- MASS::fitdistr(x = scaled, densfun = qbeta, start = list(shape1 = 2, shape2 = 2)))      Error in optim(x = c(0.01, 0.0532130412468599, 0.214202749790659, 0.268145106831258,  :         non-finite finite-difference value [2]      In addition: There were 50 or more warnings (use warnings() to see the first 50)


title: How to follow up a 3-way interaction with one continuous variable in F-test fixed effects in mixed model

text: In the Estimates of Fixed Effects SPSS-table, I don't understand why certain parameters are set to zero an just a dotted line appears (because they are redundant according to SPSS). Moreover, I don't get what the t-tests in this table mean. What is being tested here? Parameter Estimate Std. Error df t Sig. 95% Confidence Interval [treatment=,00] * [stimuluscode=-1,00] * QuestionnaireScore ,371900 ,133225 1462,000 2,792 ,005 ,110567 ,633233 [treatment=,00] * [stimuluscode=,00] * QuestionnaireScore ,379735 ,133225 1462,000 2,850 ,004 ,118402 ,641068 [treatment=,00] * [stimuluscode=1,00] * QuestionnaireScore 0b 0 . . . . . [treatment=1,00] * [stimuluscode=-1,00] * QuestionnaireScore 0b 0 . . . . . [treatment=1,00] * [stimuluscode=,00] * QuestionnaireScore 0b 0 . . . . . [treatment=1,00] * [stimuluscode=1,00] * QuestionnaireScore 0b 0 . . . . . a. Dependent Variable: inclusion. b. This parameter is set to zero because it is redundant. Any help or comments are very much appreciated! best wishes, Mariska


title: Create positive-definite 3x3 covariance matrix given specified correlation values

text: Let's say I have 3 variables, `A`, `B`, and `C` and I want to generate data where `A` and `B` are correlated at `r=x`, `A` and `C` are correlated at `r=y`, and `B` and `C` are correlated at `r=z`.   1. Is there any algorithm that will tell me, given specified values for `x`, `y`, and `z`, if any set of variances for `A`, `B` and `C` will yield a positive definite covariance matrix?    2. Given the existence of such an algorithm, is there another algorithm that, in cases where a positive definite covariance matrix _is_ possible, will give me a set of variances that achieve a positive definite matrix?


title: BY and WITH command in mixed model SPSS

text: I coded a variable as 0/1 (species 1 vs. species 2). I assumed that a dichotomous independent variable can be put in the SPSS mixed model as either a "factor" (two levels) ("by") or, if coded 0/1 also as a covariate (predictor) ("with") and that it would not make a difference. But I was wrong. The statistics do change when I use WITH vs. BY. What is the difference between the two and when should I use them? **between subject factor:** species subject: 0-1 **within subject factors:** same_different (stimulus characteristic, do subjects look at their own or the other species): 0-1 stimuluscondition: 0-1 luminance stimulus: continuous variable that I want to control for (it is not of interest) **dependent variable:** pupilsize of subject (continuous) The 'noint' keyword is used to suppress the intercept, because "with int" already defines an intercept.               mixed pupilsize **with** **int** luminance_of_stimulus lin quadr species_subject same_different stimuluscondition       / fixed int luminance_of_stimulus                 lin quadr                 species_subject  same_different  stimuluscondition                 species_subject*lin same_different*lin stimuluscondition*lin                 species_subject*quadr same_different*quadr stimuluscondition*quadr                 species_subject*same_different                 species_subject*stimuluscondition                 same_different*stimuluscondition                 species_subject *same_different*stimuluscondition                 species_subject*same_different*lin                 species_subject*stimuluscondition*lin                 same_different*stimuluscondition*lin                 species_subject*same_different*quadr                 species_subject*stimuluscondition*quadr                 same_different*stimuluscondition*quadr                 species_subject *same_different*stimuluscondition*lin                 species_subject *same_different*stimuluscondition*quadr             | noint       / random int lin | subject(id)  covtype(DIAG)       / random int lin | subject(id*condition)  covtype(DIAG)       / repeated = time | subject(id*condition)  covtype(AR1)       / print solution testcov r.               mixed pupilsize **by** species_subject **with** **int** luminance_of_stimulus lin quadr same_different stimuluscondition       / fixed int luminance_of_stimulus                 lin quadr                 species_subject  same_different  stimuluscondition                 species_subject*lin same_different*lin stimuluscondition*lin                 species_subject*quadr same_different*quadr stimuluscondition*quadr                 species_subject*same_different                 species_subject*stimuluscondition                 same_different*stimuluscondition                 species_subject *same_different*stimuluscondition                 species_subject*same_different*lin                 species_subject*stimuluscondition*lin                 same_different*stimuluscondition*lin                 species_subject*same_different*quadr                 species_subject*stimuluscondition*quadr                 same_different*stimuluscondition*quadr                 species_subject *same_different*stimuluscondition*lin                 species_subject *same_different*stimuluscondition*quadr             | noint       / random int lin | subject(id)  covtype(DIAG)       / random int lin | subject(id*condition)  covtype(DIAG)       / repeated = time | subject(id*condition)  covtype(AR1)       / print solution testcov r.


title: Least Squares Dummy Variable Regression dropping states/variables?

text: So I am running a Least Squares Dummy Variable Regression (LSDV1) involving data from 21 states observed 3 times (2007, 2008, 2009) and dropping one of the dummy value states. I have 8 independent variables that I am seeking to use (ENTER method). When I regress I find two strange things: 1) 7 of the 8 independent variables are excluded 2) 1 of my 20 states (21 - 1 dummy dropped) is excluded, too Huh UPDATE: Changed IV to independent variable               LOGS          Model Summary                    Model   R   R Square    Adjusted R Square   Std. Error of the Estimate     1   .883a   .780    .676    .07702               ANOVAa                           Model       Sum of Squares  df  Mean Square F   Sig.     1   Regression  .886    20  .044    7.466   .000b         Residual    .249    42  .006                 Total   1.135   62                    Coefficientsa                                    Model       Unstandardized Coefficients     Standardized Coefficients   t   Sig.    Collinearity Statistics              B   Std. Error  Beta            Tolerance   VIF     1   (Constant)  .854    .080        10.694  .000                 EducationandTrainingActivities  .024    .021    .202    1.126   .267    .162    6.180         StateIsAlaska   .024    .076    .039    .322    .749    .363    2.751         StateIsArizona  .109    .055    .173    1.971   .055    .675    1.481         StateIsArkansas -.333   .063    -.529   -5.302  .000    .525    1.905         StateIsCalifornia   -.108   .055    -.172   -1.956  .057    .675    1.481         StateIsDelaware -.307   .055    -.487   -5.532  .000    .675    1.481         StateIsFlorida  -.063   .055    -.101   -1.144  .259    .675    1.481         StateIsIdaho    .060    .063    .095    .950    .347    .525    1.905         StateIsIllinois .049    .055    .078    .890    .379    .675    1.481         StateIsIndiana  .061    .063    .097    .967    .339    .525    1.905         StateIsKansas   .045    .063    .071    .715    .479    .525    1.905         StateIsKentucky .019    .063    .030    .300    .765    .525    1.905         StateIsMichigan .080    .055    .127    1.440   .157    .675    1.481         StateIsMissouri -.031   .055    -.049   -.561   .578    .675    1.481         StateIsNewJersey    .115    .063    .182    1.823   .075    .525    1.905         StateIsNorthCarolina    .043    .055    .068    .772    .444    .675    1.481         StateIsUtah .047    .076    .075    .626    .535    .363    2.751         StateIsVermont  .107    .055    .170    1.932   .060    .675    1.481         StateIsVirginia .020    .076    .032    .266    .791    .363    2.751         StateIsWashington   .083    .055    .132    1.495   .142    .675    1.481               Excluded Variablesa                          Model       Beta In t   Sig.    Partial Correlation Collinearity Statistics                             Tolerance     1   ManagementActivities    .b  .   .   .   .000         PersonnelServices   .b  .   .   .   .000         PlanningandResearchActivities   .b  .   .   .   .000         InformationSystemActivities .b  .   .   .   .000         SupportServices .b  .   .   .   .000         FinanceandBudgetActivities  .b  .   .   .   .000         PublicInformationandLiaisonActivities   .b  .   .   .   .000         StateIsTexas    .b  .   .   .   .000


title: What are the chances of something *not* occurring?

text: Suppose I'm gambling using a strategy of doubling my bet whenever I lose to recoup my losses from previous bets. If my initial bet is 1/2048 of my capital, I can bet 10 times before I run out of money. Statistics says that the chances of this happening are ~1/718. As the law of large numbers states, though, this is in no way a guarantee that such an outcome will manifest in any certain number of iterations. Is it possible, therefore, to calculate the chances that this _doesn't_ happen? For example, is there an equation we can put x (the number of bets) into and determine the likelihood of this 1/718 chance _not_ occurring? My chances of losing within the first ten bets would be extremely small, and my chances after ten thousand would be pretty high, so is there any way to calculate the middle ground?


title: Calculated breeding values using markers using animal model in R

text: Animal model (frequently used in animal science and sometime in human or plants) is mixed model with: $y$ = $X$$b$ + $Z$$u$ + $e$ y is observed values for any quantitative variable, $Xb$ is fixed effects terms while $Zu$ is random effects, and $e$ is residual. Covariance between terms of $u$ can be defined as: $G$ = $A$$\sigma$$^2$$u$ . $A$ is additive relationship matrix whereas $\sigma$$^2$$u$ is additive genetic variance. $A$ can estimated from pedigree or from markers most commonly as $M$$M'$ where $M$ is marker scores (-1,0,1) x observation matrix. ![enter image description here](http://i.stack.imgur.com/ofgea.jpg) If there are no fixed effects in the model above model can be re-written as: $y$ = $Z$$u$ + $e$ As $A$ is estimated from large number of markers with multi-colinarity it is logical to shrink the $u$ toward 0 using ridge shrinkage. From the above information we can calculate **Best Linear Predictors (BLUPs).** ![enter image description here](http://i.stack.imgur.com/wWdnE.jpg) Here for individual i, the BLUP estimate is $yi$, which is summation of $ai$ is additive genetic effect and $ei$ is residual. Is there any R software package that can do this ? Examples appreciated.


title: N-sigma curves for a non-linear least square curve fit

text: I'm using `python`'s scipy.optimize.curve_fit routine (which uses a non-linear least squares) to fit an exponential function of the form:               f(x) = a * exp(b*x) + c      to a set of data. The result looks like this: ![enter image description here](http://i.stack.imgur.com/aTMpL.png) where the black triangles are the data set and the blue curve is the `f(x)` fitted by the routine. The output of the process includes the optimal value for each parameter `a, b, c` along with a covariance matrix from where I can obtain the standard deviation for each of them, ie I get: $a = p_a \pm e_a$ $b = p_b \pm e_b$ $c = p_c \pm e_c$ where $p_i$ are the optimal values and $e_i$ its standard deviations. What I want is to generate the N-$\sigma$ curves (where N is an integer) for the fit. As an example, here's an article where the authors show a similar exponential best fit (although a decaying one) along with its $1\sigma$ upper and lower curves (see caption): ![enter image description here](http://i.stack.imgur.com/1TJez.png) This site gives an example of how to do this but I'm not sure the procedure is entirely correct. How can I correctly obtain such curves with the data I have available?


title: Select most different observations

text: I have a polar chart representing eight dimensional observations as 'rings' around the center. Furthermore, I am limited to Excel for graphing purposes, which means I can only display 255 series at once. My goal is for illustrative purposes only - it will be printed, so details will get lost anyways. I would like to maximize the information content in the graph. So, given the 255 series limit, how should I select the observations to show? Is there some kind of measure to determine the amount of added information to the data, which is not affected by their order, to then drop those with the least info, until only 255 remain? I appreciate the help!


title: predict category by using K-NN algorithm having text features

text: I would like to predict the category of the provided data by using K-NN algorithm. Here is an example of the training data set               category sex      age     1        male     19-20-21     5        female   40-41-42     x        male     18        <- provided data      the features of the training data set could be numbers such as age but especially additionally text data. How could you suggest me to calculate distance since I have not data numbers but text data. By using hamming distance perhaps?


title: Large sample size for t-test

text: I am working on a physics experiment, and I ran 500 runs to calculate a parameter. I want to use the t-test to see if my mean from 500 runs matches the actual value.   * u0: mean is actual value    * u1: mean is not actual value When I tried with 100 data first, it returned a high p-value which did not reject the null hypothesis. So I was very happy with my results. However, when I tried with 500 data, it provided me low p-value (<0.05). So I have to reject my hypothesis. I find this very strange because my mean from 100 data and 500 data are very close. I believe the higher sample size changed that insignificant difference to significant difference. So what sample size is recommended? Any article that I can look for?


title: Estimating total number of people from an observed sample

text: The well known "German tank problem" shows how to answer the question: "If I have tanks which have an increasing serial number, and I see a sample of tanks and record their serial numbers, what is the likely total number of tanks". This question is analogous but is where there is no ordering to the observations, eg with people. Here's a hypothetical example (and the one I am most interested in). Suppose you go to a company website and they provide a number of CVs for staff of a particular job title (eg analyst or whatever). The question is, given this knowledge, how many staff are there likely to be with that title? To formalise this, let the number of observed people be $m$ and the total number of people be $N$. The question is then: What is $p(N|m)$. I appreciate that there may be company policies at work here, eg they may want to show all of the people on a particular level, or some representative sample. Clearly, $N \geq m$. Bayes' theorem gives that $p(N|m) = P(m|N) p(N)$. Let's ignore the prior $p(N)$ for now (or equivalent assume that it is flat), giving $p(N|m) \sim p(m|n)$. Using combinatorics, the number of ways you can get the observed $m$ people from a larger set is $N \choose m$. So immediately this implies that $p(m|N) \propto 1/ {N \choose m}$. There is a normalisation constant which is obtainable by requiring that $\sum_m^\infty C/{N \choose m} = 1$. The problem is that $1/ {N \choose m}$ is a very steeply declining function for even moderate values of $m$. For example, using $m=5$ then $p(5|5) \sim 0.8$, $p(6|5) \sim 0.13$, $p(7|5) \sim 0.038$, $p(8|5) \sim 0.014$ etc. My intuition is that if you observe five people you shouldn't conclude there is a 97% chance that there are between 5 and 7 people. What is going here? I suppose I could change my prior, but this conflicts with the fact that I want to be indifferent to the number of people (i.e. I shouldn't have to assume that 7 people are more likely than 5 people to get what I want). Help please? * * * It may not be possible to answer this question at all. From a maximum entropy perspective, all I know is that $N \geq m$. The maximum entropy distribution which describes this is just a uniform distribution on $[m, \infty)$. So perhaps I am hoping there is some trick here that will give me info that doesn't exist...


title: I have 2 problems with conditional probability, involving coins and dice

text: Hi I have these problems:               A coin is tossed three times. What is the probability that exactly two heads     occur, given that     (a) the ﬁrst outcome was a head?     (b) the ﬁrst outcome was a tail?     (c) the ﬁrst two outcomes were heads?     (d) the ﬁrst two outcomes were tails?     (e) the ﬁrst outcome was a head and the third outcome was a head?      Which I have done it like:               T1 T2 T3          H  H  H     H  H  T     H  T  H     H  T  T     T  H  H     T  H  T     T  T  H      T  T  T      I put all the outcomes, then by counting them by hand:               p(2H|T1=H)=2/4     p(2H|T1=T)=1/4     p(2H|T1=H,T2=H)=1/2     p(2H|T1=T,T2=T)=0     p(2H|T1=H,T3=H)=1/2      Is this correct? Then, is there another way of doing it? Counting by hand just feel wrong, there must be a formal way. Then:               A die is rolled twice. What is the probability that the sum of the faces is     greater than 7, given that     (a) the ﬁrst outcome was a 4?     (b) the ﬁrst outcome was greater than 3?     (c) the ﬁrst outcome was a 1?     (d) the ﬁrst outcome was less than 5?      Again counting, this really feels so wrong, I'm expecting some formula to get them.               p(sum>7|D1=4)=3/6     p(sum>7|D1>3)=12/18     p(sum>7|D1=1)=0/6     p(sum>7|D1<5)=6/24      Am I correct? Thanks.


title: pairwise variance ratio in Stata

text: Suppose, I want to compute the pairwise correlation in Stata, then I do the following:               sysuse auto     pwcorr      Now, I want to compute the variance for each variables (numeric) and then obtain variance ratios of all possible combinations. For example, I have :                  sum           Variable |       Obs        Mean    Std. Dev.       Min        Max     -------------+--------------------------------------------------------             make |         0            price |        74    6165.257    2949.496       3291      15906              mpg |        74     21.2973    5.785503         12         41            rep78 |        69    3.405797    .9899323          1          5         headroom |        74    2.993243    .8459948        1.5          5     -------------+--------------------------------------------------------            trunk |        74    13.75676    4.277404          5         23           weight |        74    3019.459    777.1936       1760       4840           length |        74    187.9324    22.26634        142        233             turn |        74    39.64865    4.399354         31         51     displacement |        74    197.2973    91.83722         79        425     -------------+--------------------------------------------------------       gear_ratio |        74    3.014865    .4562871       2.19       3.89          foreign |        74    .2972973    .4601885          0          1          .       This gives the standard deviation for all variables. Now if I want to compute the variance ratio of price and mpg, then this will be given by: var_ratio_pm=variance of price/variance of mpg and var_ratio_mp=variance of mpg/variance of price. In sum, I want a matrix of variance ratios. I wonder whether there is a command in Stata like `pwcorr` for that.


title: ICC - unclear comment from reviewer

text: I'm new here and not a statistician, so I apologize in advance and I sure hope someone can help! I'm working on a study where we have 10 locations $(5 \ intervention, 5 \ control)$. We have two measurements of activity for each participant $(n= 700)$ at each location. I've received the following comment from a reviewer: "Typically an intracluster correlation coefficient (ICC) is calculated to adjust for the correlation between participants at a single location" I know how to calculate ICC using SPSS, but because there are different numbers of participants at each location, I'm a little confused... can anyone help? Thanks so much in advance!