title: Comparing logistic regression models

text: I have a 2(truthfulness) x 2 (immediate test) study where I am looking to compare 2 logistic regression models. Both are classifying the same outcome variable: truthfulness (truthful versus deceptive -- an experimental manipulation where they were told to tell the truth or lie). The independent variable in both is a score on a delayed test. The first model uses scores for people who did the immediate test, the other model uses scores for people who did not do the immediate test. I want to know whether truthfulness can be classified **better** using scores on the delayed test when people have previously done the immediate test. However, I am confused as to how to compare logistic regression models. Is it possible to do that even though models are estimated from different samples? Is it possible to interpret differences between classification accuracy tables as meaningful? And is it possible to statistically test the differences between models?


title: Confusion related to precision and recall

text: I was reading this paper. However, I have this confusion about the results reported by them. They have given the results in terms of precision and recall. However, their results show that the precision and recall value is greater than 1. I am a bit confused why it is so? Any insights. ![enter image description here](http://i.stack.imgur.com/Nh3GZ.png)


title: Beyond the basics: intermediate medical statistics textbooks suggestions

text: I am a soon-to-be physician. During my studies I have taken a class in biostatistics. I own Martin Bland's "An introduction to medical statistics", which was a required textbooks at the time, and Harvey Motulsky's "Intuitive biostatistics", which I've purchased on my own initiative based on positive reviews on Amazon. They served me well, but I'm interested in clinical research and I feel that I need to stop scratching the surface and delve deeper to understand statistics. I'm looking for one or more textbooks I can use to self study. Ideally the textbook(s) should go beyond the basics, and explain modern statical techniques used in biomedical research (e.g reading these forums I've seen a lot of mentions of bayesian statistics, which seems to be the Next Big Thing). If possible, it should feature examples, and not assume a particular statistical package (unless it is R, that is the software I use). I should probably mentions that my math skills are... rusty. So feel free to list the kind of math skills required to understand the book(s) you are suggesting.


title: Limiting the range of possible values for distribution fitting with GAMLSS in R

text: I'm trying to fit ST5 distribution for reaction times with GAMLSS. Sometimes mu parameter of the distribution is negative, which is implausible. Is it possible to limit the range of parameter values somehow?


title: Metric to compute structural similarity between two directed graphs

text: I'm working on a small project in which I try to compare directed a-cyclic graphs. Say I have (directed) three graphs: 1)                                   X                        /   \                      /     \     START - X - X - X - X - END              \                      \                    X        2)                        X - X - X - X            /              \           /                \     START -  X - X - X - X - END           \                /            \              /              X - X - X - X       3)                                   X - X           X - X - X                        /               /          \                      /               /            \     START - X - X - X - X - X - X - X - X - X - X - END      Note, that the labels (X) are of no importance; I'm only interested in the structural properties and similarities between the graphs. Do you have any pointers to literature or any clever ideas how to compute the structural distance or similarity between such graphs? Many thanks!


title: Detecting changes in distribution of multiple variables

text: I am a bit new to this field. So I needed help in finding out which topic should I focus on for achieving this. Suppose I have N dependent random variables. I have n samples of each of these random variables. Now I want to check if the kth sample follows the pattern that is suggested by the previous (k-1) samples. For example if there are 2 variables - X1 and X2. For the first 10 samples lets say X1 > X2. But if for the next sample X1 < X2, then this sample is anomalous. I need to detect such an anomalous sample. Can you tell me what statistical machine learning concepts should I study in order to understand which techniques are applicable in such a scenario? From my limited knowledge this appears to be a classification problem that must be attacked with unsupervised technique.


title: Cost functions for weighting sensitivity and specificity in binary classification problem

text: I'm searching for a combination of `sensitivity` and `specificity` cost function because i want have more weight for sensitivity ( sensitivity is more impotent for me rather than specificity). After searching i found this :               Final_Cost = ( (Cb/Cg)/( 1+(Cb/Cg) )*Bg + ( 1/( 1+(Cb/Cg) ) )*Gb      `Cb` is misclassification cost of positive and `Cg` is misclassification cost of negative. `Bg` is number of false positive detected and `Gg`is number of false negative detected. We should specify `Cb/Cg`. Is this a good function for calculating cost? Is there any other better functions?


title: Bayesian estimation of $N$ of a binomial distribution

text: This question is a technical follow-up of this question. I have trouble understanding and replicating the model presented in Raftery (1988): _Inference for the binomial $N$ parameter: a hierarchical Bayes approach_ in WinBUGS/OpenBUGS/JAGS. It is not _only_ about code though so it should be on-topic here. ## Background Let $x=(x_{1},\ldots,x_{n})$ be a set of success counts from a binomial distribution with unknown $N$ and $\theta$. Further, I assume that $N$ follows a Poisson distribution with parameter $\mu$ (as discussed in the paper). Then, each $x_{i}$ has a Poisson distribution with mean $\lambda = \mu \theta$. I want to specify the priors in terms of $\lambda$ and $\theta$. Assuming that I don't have any good prior knowledge about $N$ or $\theta$, I want to assign non-informative priors to both $\lambda$ and $\theta$. Say, my priors are $\lambda\sim \mathrm{Gamma}(0.001, 0.001)$ and $\theta\sim \mathrm{Uniform}(0, 1)$. The author uses an improper prior of $p(N,\theta)\propto N^{-1}$ but WinBUGS does not accept improper priors. ## Example In the paper (page 226), the following success counts of observed waterbucks are provided: $53, 57, 66, 67, 72$. I want to estimate $N$, the size of the population. Here is how I tried to work out the example in WinBUGS ( **updated** after @Stéphane Laurent's comment):               model {          # Likelihood       for (i in 1:N) {         x[i] ~ dbin(theta, n)       }          # Priors          n ~ dpois(mu)     lambda ~ dgamma(0.001, 0.001)     theta ~ dunif(0, 1)     mu <- lambda/theta          }          # Data          list(x = c(53, 57, 66, 67, 72), N = 5)          # Initial values          list(n = 100, lambda = 100, theta  = 0.5)     list(n = 1000, lambda = 1000, theta  = 0.8)     list(n = 5000, lambda = 10, theta  = 0.2)      The model does _sill_ not converge nicely after 500'000 samples with 20'000 burn-in samples. Here is the output of a JAGS run:               Inference for Bugs model at "jags_model_binomial.txt", fit using jags,      5 chains, each with 5e+05 iterations (first 20000 discarded), n.thin = 5      n.sims = 480000 iterations saved              mu.vect  sd.vect   2.5%     25%     50%     75%    97.5%  Rhat  n.eff     lambda    63.081    5.222 53.135  59.609  62.938  66.385   73.856 1.001 480000     mu       542.917 1040.975 91.322 147.231 231.805 462.539 3484.324 1.018    300     n        542.906 1040.762 95.000 147.000 231.000 462.000 3484.000 1.018    300     theta      0.292    0.185  0.018   0.136   0.272   0.428    0.668 1.018    300     deviance  34.907    1.554 33.633  33.859  34.354  35.376   39.213 1.001  43000      ## Questions Clearly, I am missing something, but I can't see what exactly. I think my formulation of the model is wrong somewhere. So my questions are:   * Why does my model and its implementation not work?   * How could the model given by Raftery (1988) be formulated and implemented correctly? Thanks for your help.


title: Compressed Population Complexity in Minimum Description Length (MDL)

text: I am studying the MDL and found it is sum of model complexity and compressed population complexity. To my understanding, model complexity refers to number of bits to encode the model, which can be computed by number of bits per sample (i.e. log N) multiply by number of permutations. On the other hand, the compressed population complexity is the number of data its can compressed. However, it is usually estimated by the average number of bits for encoding the data multiply by the number of population (i.e. N), which can be computed by a entropy function. However, can anyone help me to understand how can the compressed population complexity can be estimated by the average no. of bits for data encoding. thanks.


title: Markov Switching Regimes

text: I have one pending question about finding the transition probabilities in the Markov Switching method. I know that Eviews 8 has this function, but since the latest version I have is Eviews 7, I'm trying to do this method in excel (my God i hope this is possible). My question is, how does one or on what basis does one compute or get the values for mean for regime 1 and mean for regime 2, their standard deviations, and p11 as well as p22? A lot of the excel files that I've been able to download online (for reference purposes) show these parameters as inputs but do not specify how they get them. I'm attaching as well the paper that I'm currently referring to just in case you're curious: "Analysis of US Business Cycle with a Vector Markov Switching Model" by Zenon G. Kontolemis. enter link description here


title: Sample size calculation for paired ordinal data

text: What formula should be used for sample size calculations for paired ordinal data? Imagine we have an intervention and we want to measure the effect based on before and after measurements of an ordered outcome (say a Likert scale). I know how to calculate for a t-test for paired continuous data and how to calculate for a Mann-Whitney test for ordinal data (Whitehead 1993), but what about paired ordinal? Either a nonparametric solution or a parametric solution through a GLM for ordinal response would work I think.


title: Multivariate gaussian log-likelihood

text: Let's say that we have a 2x5 matrix A where rows correspond to observations and columns correspond to variables from a p-variate Gaussian distribution, and we want to learn the inverse covariance matrix structure. Since n is smaller than p, we cannot invert the sample covariance matrix, and let's say that we have a 5x5 matrix "theta" as an estimate of the inverse covariance matrix of the data A. Then, the log-likelihood of this data under inverse covariance estimation "theta", using the log-likelihood formulation for multivariate gaussian is (I am using MATLAB notation):                   -n/2*p*log(2*pi) + n/2*log(det(theta))-n/2*trace(cov(A)*theta)      Then, let's say, we have one more 2x5 matrix B, and we are using the same estimate "theta" of the inverse covariance matrix for this data. Then, I want to calculate the log-likelihood of all 4 observations under estimation "theta", and I was expecting that the two ways loglik(A)+loglik(B) and loglik([A;B]) would give us the total likelihood result and would be exactly equal to each other. Log-likelihood of a dataset is always calculated by multiplying likelihoods for each of the observations (and in the log scale this corresponds to additions), right? Then, log-likelihood of the augmented matrix with all observations should be equal to sum of the log-likelihoods of independent matrices. But they are not, when I try on MATLAB. Do you have an idea about why? Below is the MATLAB code I am running:                   function loglik = exactMVGloglikelihood(X, theta)             n = size(X, 1);             p = size(X, 2);             loglik = -n/2*p*log(2*pi) + n/2*log(det(theta))-n/2*trace(cov(X)*theta);         end              A = randn(2,5);         B = randn(2,5);         theta = randn(5,5);         loglik1 = exactMVGloglikelihood(A, theta) + exactMVGloglikelihood(B, theta);         loglik2 = exactMVGloglikelihood([A;B], theta);      I expect loglik1 and loglik2 to be exactly equal, but they are not. Any help is welcome. Thanks.


title: Reference for distribution of sample coefficient correlation

text: Some one please share to me the paper to understand how to derive the $E(r)$ and $Var(r)$ for sample coefficient correlation (R) thanks..


title: adf.test in R, define hypothesis?

text: How do I define a hypothesis using adf.test (tseries) in R? When I use adf.test, I would like to add a intercept or a trend in the hypothesis, since ADF has three different hypotheses. I would like to test the following hypotheses: H0: Random walk H0: Random walk with drift H0: Random walk around a trend with drift Is it possible?


title: Modifying decision function in LibSVM

text: The decision function for C-support vector classification is $sign(wT\phi(x)+b)=sgn(\sum_{i=1}^{l}y_i\alpha_iK(x_i,x)+b)$ How can I modify it to $sgn(wT\phi(x)+b)=sgn(\sum_{i=1}^{l}y_i\alpha_iK(x_i,x)+b+aconstant)$ that aconstant is a ARBITRARY number


title: Calculating the probability of an inequality with two random variables

text: I am analyzing a timing circuit I designed, and I need to calculate the probability of a certain event (bit error). For example, I have derived this equation: $(1 + x) / d < 1 / M$, where x is a random variable with a normal distribution, d is a random variable with uniform distribution, and M is a constant. I would like to be able to calculate the probability that this inequality is true, for a given distribution for x, a given distribution for d, and a given M. For example, let's say x has a normal distribution with mean = 0 and SD = 0.05, d has a uniform distribution between 0 and 1, and M = 3. How would I calculate the probability that the inequality is true? If there was a single random variable, I could do it. But I am lost and don't know how to solve it. I tried searching and reading a statistics textbook, but I don't know what to look for.


title: What is the reason the $\log$ transformation is used with right-skewed distributions?

text: I once heard that > log transformation is the most popular one for right-skewed distributions in > linear regression or quantile regression I would like to know is there any reason underlying this fact? Why log transformation is fitful to the right-skewed distribution? How about left- skewed distribution?


title: Marginalization in gaussian

text: If $\theta|\mu \sim N(\mu,\sigma_o^2)$ and $\mu \sim N(0, \sigma_1^2)$ what is the marginalized $P(\theta)$. $\theta$ and $\mu$ both are nx1 vectors $P(\theta) = \int P(\theta|\mu)P(\mu)d\mu$ Is it $N(0,\sigma_o^2I+\sigma_1^2A)$? where A is nxn matrix with all ones?


title: How many samples are needed for a paired trial with multiple tries?

text: I'm having problem understanding how big sample size is needed for a situation where there are pairs in the before/after groups but these have all multiple data points. ## Problem Pose a situation where you have X dart throwers. Each dart thrower have thrown $n_i$ darts for some average $m_i$ points with a "bad" dart. Now an expected better dart is introduced, so the dart throwers all start throwing with the new dart. How many throws are needed in order to be able to say that the new darts are better at say 5% significance with 80% power? (Note that this will be a function of the actual difference in average points.) ## Thoughts on solution I am aware of the `power.t.test`, which can be used for to compare all the points before and after. However, as the dart throwing expertise could vary substantially amongst the throwers, it should probably be a paired test. But, when reading up on the `power.t.test(..., type = "paired")`, it tells you how many pairs (throwers) are needed (and presumably assumes only one data point with the bad and good dart). So, I'm considering to do a `power.t.test` for each of the authors but how would I assemble the results into one? Is there a paired power test which could answer how big sample is needed in the situation where you have pairs but also multiple data points for each pair? Does the original distribution matter?


title: How is cross validation different from data snooping?

text: I just finished "An Introduction to Statistical Learning". I wondered whether using cross-validation to find the best tuning parameters for various machine learning techniques is different from data snooping? We are repeatedly checking which value of the tuning parameter results in a best predictive result in the test set. What if the tuning parameter we arrive at just happens to fit this particular test set by chance, and won't perform well on some future test set? Please excuse my novice understanding of machine learning, and I'm eager to be educated. EDIT: Please see @AdamO answer on the definition of "data snooping". I used the term very inaccurately in my question.


title: What is a good way of testing for a relationship between two count variables?

text: I have counts of occurrences of two types of words (A and B) in several texts. What I would like to test is whether the frequencies of occurrence of both types of words across texts is 'correlated'. However, using Pearson's correlation is probably not correct, because my data is not continuous, and in addition the counts are often quite low (sometimes zero). What is a good way to test my hypothesis?


title: What is the estimation techniques used in lm() in R?

text: I'm wondering what is the estimation techniques used in lm(). If it's OLS, how could we perform a log likelihood test by logLik()? What's the difference between lm() and ols(), mle() and other functions?


title: Latent Class Models

text: What is difference(Or relation) between "Latent Class Analysis"and Structural Equation models" and "Latent Growth Curve Models"?


title: WinBUGS multinomial all parameters are random variables

text: everyone here. I want to use multinomial distribution, but p[] and N are all random. model is correct and data loaded , model compilied, but it can't still run , it says:value of binomial f[1] must be an integer. Is there anyone here know about it? thank you very much . this is my coda:               model {        for(s in 1:S) {        for(i in 1:E) {              # stochastic part of the likelihood              Sample[s,i] ~ dpois(lambda[s,i]);              # link function              lambda[s,i] <- (w[s]*l[i])*inprod(A[i,1:T],f[1:T]);              }         }       for(m in 1:T) {        delta[m] ~ dgamma(alpha[m],1);         p[m] <- delta[m] / sum(delta[]);        }       #  p[1:T] ~ ddirch(alpha[1:T])            N ~ dpois(lambda1);       q[1]<- p[1];       for(k in 2:T){            q[k] <- p[k]/(1-sum(p[1:(k-1)]));        }      f[1] ~ dbin(q[1],N);      for(j in 2:(T-1)) {           tmp[j] <- N-sum(f[1:(j-1)])           f[j] ~ dbin(q[j],tmp[j])      }      f[T]<-N-sum(f[1:(T-1)])          }


title: WinBUGS multinomial with unknown parameter and size

text: I want to use multinomial distribution, but i don't know the `p[]` and `N`. So `p[]` and `N` are all random variables. This is my winBUGS code:               model {       for(s in 1:S) {        for(i in 1:E) {              # stochastic part of the likelihood              Sample[s,i] ~ dpois(lambda[s,i]);              # link function              lambda[s,i] <- (w[s]*l[i])*inprod(A[i,1:T],f[1:T]);              }          }              for(m in 1:T) {                 delta[m] ~ dgamma(alpha[m],1);                 p[m] <- delta[m] / sum(delta[]);         }        #  p[1:T] ~ ddirch(alpha[1:T])             N ~ dpois(lambda1);        q[1]<- p[1];        for(k in 2:T){             q[k] <- p[k]/(1-sum(p[1:(k-1)]));          }         f[1] ~ dbin(q[1],N);         for(j in 2:(T-1)) {                tmp[j] <- N-sum(f[1:(j-1)])                f[j] ~ dbin(q[j],tmp[j])         }        f[T]<-N-sum(f[1:(T-1)])           }      The model is correct, the data loads, and model compiles, but it won't run. It says: > value of binomial f[1] must be an integer. Does anyone here know about this? Thank you very much.


title: Conditional on binomial question

text: Here goes the question from textbook: A small restaurant seats 20 diners, and is full every night. The chef knows from previous experience that 40% of the diners order steaks, so she always has 12 in her fridge at the beginning of the evening. If a customer orders steak, what is the probability that he/she will receive one? Apparently, the solution in here https://nz.answers.yahoo.com/question/index?qid=20120626203618AAsWm24 talks about correction factor. What is the correction factor? and where may I apply it in this context? Thanks in advance for any suggestion related to the theory of this question.


title: Regression: Comparing Two Slope Coefficients

text: Let's say we have two different variables (stress and attitude). I am analyzing how the influence confidence levels. When I determined the regression equation, the slope coefficient of stress on confidence is -.15. Then I did the regression equation for attitude, and it had the slope coefficient of .7. Based on these two slope coefficients, can I determine if one of these two variables is more useful in predicting confidence, or do I need to do other steps, like calculating the correlation coefficient? Thanks!


title: How to deal with unequal sample sizes while fully embrace a dataset?

text: **_Imagine the situation:_** Mythical Seafolk use holes in the seabed as their burrows. Each hole has two parameters - diameter and depth. **Majority of holes are unoccupied** due to their surplus (n = 235). Occupied holes (n = 15) are (generally) expected to be much deeper and with larger diameter than random.               # generate data     set.seed(1234)     x <- runif(250, min=0, max=10)     y <- runif(250, min=0, max=10)     rbPal   <- colorRampPalette(c("deepskyblue","darkblue"))     my.data <- data.frame(x_coor = x, y_coor = y,                           diameter = c(abs(rnorm(15)+2.5), abs(rnorm(235))),                           depth = c(abs(rnorm(15)+2.5), abs(rnorm(235))),                           usage = rep(c(1,0), times = c(15, 235)))     my.data$col <- rbPal(10)[as.numeric(cut(my.data$depth,breaks = 10))]          # look at the situation     # occupied holes are marked by red circles     plot(my.data$x, my.data$y, cex = my.data$diameter, col = my.data$col, pch = 20)     grid(5, 5, lwd = 0.75, lty = 2, col = "grey")      points(my.data$x[1:25], my.data$y[1:25], pch = 1, cex = 2,            col = "red", lwd = 1.75)      ![enter image description here](http://i.stack.imgur.com/8O1Xg.png) Although this analysis seem to be rather straightforward, **I have some doubts about comparing unequal samples (15 vs 235)**. Such problem do not occur in other _Seafolk burrow selection studies_ because **mapping of seabed is very expensive and time demanding** and when researchers find 10 used holes, they continue to map the seabed only till they get find additonal 10 nonused holes (which are most probably located in close vicinity of used holes). Due to this data collection approach the sample sizes of compared groups are usually almost equal. **However, we were able to find ALL holes in seabed - which can be now considered as a "handicap".** **_We do not want to do random sampling of 15 holes from 235 to obtain equal sizes (15 vs 15). Why to throw away such large amount of data!? Is there any solution/approach to this situation which is able to fully embrace such unique dataset?_** _I have in mind two possible ways (please comment or add others):_ **1.** Take the sample of unoccupied holes (n = 235) and randomly resample them to (for example 50) smaller groups (all with n = 15). Subsequently, each of these reduced groups will be compared with sample of occupied holes by binomial logistic regression. By this I will obtain 50 logistic curves which will be "averaged" into just one curve - final model. **2.** Take the sample of unoccupied holes (n = 235) and randomly resample them to (for example 50) smaller groups (all with n = 15). Subsequently, from each of these reduced groups (and from sample of used holes also) we gain (mean, sd, sample size). Thus, we will have six numbers for each comparison and we calculate effect size - standardized mean difference (SMD). Effect sizes will be inserted into meta-analytic model and tested if there is a heterogeneity present among effect sizes ( _to see if the "randomnes" of creating unused holes sample has the significant impact on test result_ ).


title: Determining the confidence interval for a non-normal distribution

text: I have a set of data I ran on a simulation using R with a population size (N) of 1000, sample mean of 64.93, and a standard deviation of 27.61. The distribution is positively skewed and a Wilks-Shapiro test done in JMP shows that the data is not normally distributed. I need the confidence intervals of this data to be able to test some experimental data I have. I have learned that using the `boot.ci()` function in R can give me confidence intervals using bootstrapping. As shown in http://www.statmethods.net/advstats/bootstrapping.html, the boot function requires: `bootobject <- boot(data= , statistic= , R=)`, where data is the data of interest, statistic is a function that produces the statistic to be bootstrapped, and R is the number of times to do it. My data is simply a vector of 1000 samples. I am not sure what statistic I should do. All the examples I have read online require two data sets to make a correlation statistic or a linear regression statistic. The data can be downloaded from: http://ktdllc.com/data/simp2.csv Any help on this issue would be greatly appreciated.


title: Building a flexible predictive model to put together fashionable outfits

text: I have a rather open-ended question that I'm hoping people can help me brainstorm for. Imagine you wanted to design an application to help you put together outfits from your closet.   1. To simplify the problem, imagine each outfit needs two components, a top and a bottom.   2. Each top and bottom has a set of features: i.e., color, fabric thickness, etc.    3. There are certain guiding principles that predict whether a given top will go with a given bottom.  For example, colors on opposite sides of the color wheel go well together; tops made out of thick fabric go with bottoms made out of thick fabric. One way is to implement a formula that predicts how fashionable the outfit will be given the guiding principles described above. Now, what if you want to allow your "formula" to adapt based on feedback? How would you go about implementing that? The simpler the model, the better. I've heard something about using hierarchical Bayesian networks (which frankly, I don't understand well) but was wondering what the good folks at Cross Validated would suggest. * * * By the way, I'm not actually going to build the application described above. I'm more interested in it as a thought experiment and a way to learn more about predictive models.


title: If I compare conditions within two preexisting groups is that a quasi-experimental design?

text: I have two groups of participants, one group of psychosis patients and one group of controls. All of the participants complete a computer task that has three different conditions. So, I want to make both within- and between-groups comparisons. Since the differentiating factor between the groups is illness condition (something that can't be manipulated), I figured I run a mixed design ANOVA (split-plot ANOVA). Thus, labeling psychosis and healthy as two independent variables. Is this analysis right? If it is right, what do I call this type of design? Would it be appropriate to call it "Quasi-experimental independent groups repeated measures design"?


title: Different Parameter Estimates of SAS and Minitab for Unreplicated Factorial Designs

text: I am using famous example from Myers and Montgomery (1995), which is $2^{4}$ factorial design to compare some methods for unreplicated factorial designs concerning active effects and mean squared errors of the predicted values. However, in order to run the analysis, I need to get parameter estimates for full model, which gives zero degrees of freedom for error in ANOVA. ![enter image description here](http://i.stack.imgur.com/1Jdgz.gif) I should use SAS and Minitab. You can see the SAS code I use below                   data rate;               do a=-1,1; do b=-1,1; do c=-1,1; do d=-1,1;                  input y @@;                  ab=a*b; ac=a*c; ad=a*d; bc=b*c; bd=b*d; cd=c*d;                  abc=a*b*c; abd=a*b*d; acd=a*c*d; bcd=b*c*d;                  abcd=a*b*c*d;                  output;               end; end; end; end;               datalines;               45 71 48 65 68 60 80 65 43 100 45 104 75 86 70 96               ;            run;         ----------------------------         proc reg data=rate outest=regout;               model y=a b c d ab ac ad bc bd cd abc abd acd bcd abcd;            run;      In Minitab, I am using "Stat > DOE > Factorial > Analyze Factorial Design" section Here are the outputs for SAS and Minitab, respectively. ![SAS output](http://i.stack.imgur.com/ahn8Y.png) ![Minitab Output 1](http://i.stack.imgur.com/mI2YE.png) ![Minitab Output 2](http://i.stack.imgur.com/lkEUs.png) As you can see, both programs give different parameter estimates. I am pretty sure that SAS uses Ordinary Least Squares (OLS) for regression. For minitab, I am not certain about it uses OLS or Maximum Likelihood in order to fit a general linear model for factorial designs. But there's no official document I could find. My first question is: Do you have any official document for Minitab regarding the usage of OLS for parameter estimations of a general linear model? My second question is: How can I deal with that kind of problem? For example, I would like to find active effects using Lenth's method. For that, I need to get effects (which is two times of parameter estimates), as I can see these two programs give different contrasts.


title: how to shape the dataframe for mlogit in R

text: I have a question concerning `mlogit` function. In my experimental design I collect 80 trials for each participant . The subjects are asked to provide an answer by choosing among 4 different alternatives (not ordered). Now, I would like to model the choices of my participants as a function of nine continuous predictors, however I can’t understand how the data frame must be re-shaped in order to be used with `mlogit`. Does anybody know what would be the proper way to arrange the dataset? Here is an extract of the dataset:               # Subject Trial Choice  X1      X2      X3      X4      X5      X6      X7      X8      X9     Subj1       1   2       0.00    0.62    0.04    -0.84   0.12    -0.06   0.58    -0.05   0.04     Subj1       2   2       0.00    0.62    0.04    -0.84   0.12    -0.06   0.58    -0.05   0.04     Subj1       3   4       0.02    -0.01   0.02    -0.23   0.01    0.30    0.12    -0.01   0.01     Subj1       4   3       0.01    -0.02   0.03    -0.24   0.03    0.37    0.12    -0.02   0.03     Subj1       5   3       0.01    -0.02   0.03    -0.24   0.03    0.37    0.12    -0.02   0.03     …           …   …        …        …       …       …      …       …       …        …       …     Subj1       80  3       0.01    -0.02   0.03    -0.24   0.03    0.37    0.12    -0.02   0.03     Subj2       1   3       0.01    -0.04   0.03    -0.22   0.03    0.33    0.10    -0.02   0.03     Subj2       2   2       0.02    0.01    0.02    -0.12   0.01    0.26    0.07    -0.01   0.01     Subj2       3   1       0.03    -0.01   0.02    -0.12   0.02    0.23    0.07    -0.02   0.02      Subj2       4   2       -0.02   0.00    0.00    -0.19   -0.01   0.35    0.06    -0.04   0.01      Subj2       5   2       0.01    -0.01   -0.01   -0.10   0.00    0.27    0.05    -0.01   0.00     …           …   …        …        …       …       …      …       …       …        …       …     Subj2       80  2       0.01    -0.01   -0.01   -0.10   0.00    0.27    0.05    -0.01   0.00     Subj3       1   4       0.02    0.31    0.01    -0.20   0.01    0.04    0.12    -0.01   0.01     …           …   …        …        …       …       …      …       …       …        …       …


title: Expert forecasting software evaluation

text: I have question on evaluating forecasting software with expert systems: Are there any objective assessment (not from the manufacturer) on expert forecasting software in the literature? I found two references but both are old:   1. Principles of Forecasting (2001) Armstrong   2. Forecasting methods and applications (1998) Makridakis et al. SCA and AUTOBOX seem to be two expert systems in the market, are there any comparisons made on these two software? Are there any other software that have expert systems? Any guidance in this regard would be greatly appreciated.


title: Selection of regressors

text: I have a question about the package **leaps** which I am using for model selection. I would like to compare 4 different selection methods: forward, backward, stepwise and best subset. I used the code below:               library(leaps)     forward <- regsubsets(Response ~.,data = mydata, method = "forward", nbest=1)       backward <- regsubsets(Response ~.,data = mydata, method = "backward", nbest=1)     stepwise <- regsubsets(Response ~., data = mydata, method = "seqrep", nbest=1)     best subset <- regsubsets(Response ~.,data = mydata, method = "exhaustive", nbest=1)     # adjusted R2     opt = par (mfrow =c(2,2))     plot(forward, scale = "adjr2", main = "Forward Selection")     plot(backward, scale = "adjr2", main = "Backward Selection")     plot(stepwise, scale = "adjr2", main = "Stepwise selection")     plot(best subset, scale = "adjr2", main = "Best subset selection")      Using these commands I obtained figures below: ![enter image description here](http://i.stack.imgur.com/KlDWC.png) I am wondering why figure A and D are similar to each other (and also figure B and C). I would expect different algorithms to select models in a different way. For instance models selected with forward selection method should be chosen based on the significance level/ AIC value. On the other hand models selected with best subset selection method should be chosen based on the sample statistics. I am also wondering why forward selection does not choose one variable at the time adding it to the existing model? Also Fig B shows that backward selection starts with eight variables in the model. Why it does not start with all the variables and excludes one at the time? I would be very grateful for these clarifications.


title: Implementing Neural Network for time series

text: I am currently working on neural networks for time series forecasting, my doubt is do we need to account for issues like trend,non stationarity and seasonality while using neural networks as opposed to the box jenkins methodology. If yes, then how do we account for it . For instance do we need to change the cost function etc.I am currently using nnet and neuralnet package in R. But I don't find mention of these issues in their documentation. If no , then does it imply that the model learns about seasonality from data itself, because for instance their is a time series with an upward trend, along with some variations, and if we fit a neural network on such series then it might learn just an interpolation in that case


title: Bias items and probability estimates in LibSVM

text: I have two questions in using LIBSVM   1. The decision function for C-support vector classification is  $$\text{sgn}\left(w^T\phi(x)+b\right)=\text{sgn}\left(\sum_{i=1}^ly_i\alpha_iK(x_i,x)+b\right)$$ How can I understand the physical meaning of $b$? How can I enable LIBSVM to output these two types of information, $w$ and $b$?   2. LIBSVM is able to output probability estimates. When reading the articles of SVM, I also read similar term like confidence? I am not sure whether they are talking the same thing?


title: How to calculate number of participants required to compare mean scores on a questionnaire between two groups?

text: I have to create a questionnaire that compares the sportiness of Facebook users vs. non-Facebook users. I need to calculate the number of participants that are required so that the result is statistically relevant. Unfortunately I have no idea how I can calculate that. Could you help me?


title: F Test: $F_{4, 81}(5\%) \approx F_{4,60}(5\%) - (81 - 60)(F_{4, 60}(5\%) - F_{4,120}(5\%))/60$. Where does this come from?

text: I had to compute the F-test in my ANOVA question and use interpolation if necassary. The first time, I had to work out the values of $F_{4,90}(5\%)$ which I said was approximatley $F_{4, 60}(5 \%)$ and then did my test and got the right answer. In the second bit, in my 3-way ANOVA table, when I'm checking to see if two of the factors interact, I need to get the values for $F_{4, 81}(5 \%)$ and so I thought I'd do the same thing and so its roughly $F_{4, 60}(5 \%)$ and then do it like this. However it turns out you get something like: $$ F_{4, 81}(5\%) \approx F_{4,60}(5\%) - (81 - 60)(F_{4, 60}(5\%) - F_{4,120}(5\%))/60$$ Where does this come from and why is it like this?


title: Calculating factor scores in SPSS

text: I have done a factor analysis on different attributes of different cigarette brands (e.g. how expensive consumers find it, etc.). Three factors were calculated out of the many attributes after reduction. Now I would like to do a cluster analysis (euclidean distance and Ward method) on factor scores calculated for each person. THe question is how can I calculate factor scores per person in SPSS or Excel?


title: Do we need to apply the same transformation of predictors on a test dataset?

text: This maybe a stupid question. But I hope someone could help me. I first divide the dataset into training (75%) and test (25%). Then fit a logistic regression model on training data set. When fitting the model, I did some modification on independent variable, such as normalization, transformation stuff. Now I have the model and ready to test it using test data set. But do I need to do the same transformation to the independent variable in the test data set?


title: Determine significant transitions in data values

text: I don't have much experience with statistics and suspect this is a common problem with a common solution. Normally I would just eye-ball the data and call it a day, but at some point you ought to learn the right way to do things. I have a dataset of about 17,000 areas with a cumulative distribution:               gid          st_area2d           cume_dist      9810     0.0000280288363734472   0.0000568052715291979      10804    0.0000387003730889654   0.0001136105430583960      3087     0.0000492940355178313   0.0001704158145875940      7037     0.0000504638710612963   0.0002272210861167920     ...      When I make a scatter plot of the area and percentile there are two clear transitions- one where the percentiles were very small and another where they become very large (relatively): ![CSUM Chart](http://i.stack.imgur.com/ZDOBW.png) How can I best determine exactly which percentile marks those transitions? So I want a way to determine that somewhere around the 95th percentile my areas get much larger very quickly and somewhere are 98 it levels off again. My guess is comparing the slope between points; at both transitions the slope will become significantly less steep. I suspect there is standard method, but I don't know the terminology or have the background to defend one method over another.


title: Writing out the log likelihood for a mixed effects model with multiple variance components

text: For a given model: $y_{ik} = \mu_k + \beta X + u_k X + b_i + \epsilon_{ik}$ where $u \sim N(0, \sigma_{u}^{2} I)$, $b \sim N(0, \sigma_{b}^2 I)$ and $\epsilon \sim N(0,\sigma_{e}^2 I)$. Assuming the $\epsilon_i$ are normally distributed, I am trying to write out the log likelihood for the above equation. However, given that there are multiple variance components in the equation, how do I go about doing this? Here's what I have so far. $l(Y; X, \beta, \sigma_{u}, \sigma_{b}, \sigma_{e}) = - \frac{n}{2} log(2\pi + V^2) - \frac{1}{2V^2} \sum_{i=1}^{n} (y_{ik} - \mu_k - \beta X - u_k X - b_i)^2$ where -- $V = \sigma_u^2 + \sigma_b^2 + \sigma_e^2$


title: How to testing the probability of two factors being found together?

text: I know how many times a particular factors has been found: `values<-cbind(f=c("S1","S10","S3","S4","S5","S6","S7","S8","S9"),count=c(9,8,4,24,20,4,8,21,5))` I also have information on the number of times factors are found alone, in pairs, triplets ...               together<-cbind(count=c(25,11,11,1,1,0,2,0),fnb=c(1,2,3,4,5,6,7,8))      S4 and S5 are found together in 13 cases which seems very high compared to any of the other associations. I would like to work out whether these two factors are more frequently found together than you would expect if the factors were randomly associated. How do you formulate this into a statistical test? Any help is very much appreciated.


title: how to assess multilevel model assumptions using residual plots

text: I am not clear as to how to assess if a multilevel model fit using `lmer` satisfies the assumptions of Normality and Homoscedasticity? I have used the following `r` code and i find that histogram of residuals and qq plot satisfy the assumption that residuals are normally distributed. But as multilevel models have residual both at individual level as well as at the group level, should the residuals at both the levels be normally distributed? Any help is greatly appreciated.               residuals<-resid(results)     summary(residuals)     hist(residuals)     qqnorm(residuals)     qqline(residuals)     xyplot(resid(results) ~ fitted(results))


title: Fisher's exact test in 3x2 contingency table

text: I have two groups of patients (A and B) with a congenital malformation which might present itself in 3 forms (a or b or c). Sample sizes are small as you can see, so I think the best test to check whether there's a statistically significant difference between the 3 forms in the 2 groups is a Fisher's exact test in a following 3x2 contingency table:               group A: 2 a, 12 b, 1 c     group B: 5 a,  3 b, 1 c      My question is: how should I interpret the p value? I don't understand what is that referred to. Having the p value, how can I say that one of the three forms is statistically significantly more represented than the others (if true)?


title: Issue of multiple scaling method

text: I plan to use multiple scaling in R and start with a toy example. There are two matrices. The first contains 10 observations and the second one contains two replicated rows. It appears that the repeated rows are projected to different positions. Intuitively I expect the replicates are projected to (approximately) the same position in lower dimensional plot but seems not. Below is the code:                aa <- matrix(sample(1:10,1000,replace=T),nrow=10)      bb <- rbind(aa,aa[c(3,5),])      dist1 <- dist(aa)      dist2 <- dist(bb)      fit1 <- cmdscale(dist1,eig=TRUE, k=2)      fit2 <- cmdscale(dist2,eig=TRUE, k=2)      par(mfrow=c(2,1))      x1 <- fit1$points[,1]          y1 <- fit1$points[,2]      plot(x1,y1)           x2 <- fit2$points[,1]          y2 <- fit2$points[,2]      plot(x2,y2)      So I am not sure if anything wrong or that is the nature of this method. In this case, what would be the explanation or solution (project to a higher dimension, choose other methods)?


title: Characterization of normal distribution

text: How can one prove that the random sample is normally distributed if the sample mean and sample variance are independent?


title: Does rpart use multivariate splits by default?

text: I know that R's `rpart` function keeps the data it would need to implement multivariate split, but I don't know if it's actually performing multivariate splits. I've tried researching it online looking at the `rpart` docs, but I don't see any information that it can do it or is doing it. Anyone know for sure?


title: How to conduct community division of a social network with R?

text: I am trying to use R to conduct community division within my weighted network (based from an association matrix). I tried with igraph but I encountered some problems. I usually use the program Socprog (Whitehead 2009) for my analysis, but as I would like to conduct a community division with Newman modularity (2006) on 1000 bootstraps of my data in order to calculate the comembership matrice resulting from all the community division of the bootstraps. So I am looking for a script that would conduct a community division (Newman 2006) from an association matrice and that would give a list of individuals with their cluster. Does anyone have a R script that works for this analysis? Any advice to implement this analysis would be much apreciated. Best regards


title: Fitting GAMM model in R

text: when fitting GAMM with R, I would like to know why when the smooth function is linear, the confidence interval is zero around the middle (the dotted lines are crossing each other at the middle)


title: Statistically correct to compare two Spearman's $\rho$-Vaules?

text: I have different variables (A, B, C, D) which are all ordinal scaled. I have another variable (Z) which is ordinal scaled, too. My sample size is approx. 1.500. Is it in the statistical sense correct, to calculate the Spearman's $\rho$ for each relationship - (1) A & Z, (2) B & Z, (3) C & Z and so on - and compare these values subsequently? Maybe using the Fisher Z-Test? Possible result: None of the relationships differs significantly from another one but (1) is slightly stronger than (2) but not as strong as (3)?


title: DBSCAN & OPTICS

text: I'm looking for DBSCAN and OPTICS matlab codes for clustering 'm' pattern in 3D space: the imput is a matrix where the 'm' rows are vectors a the 'n'=3 columns are X,Y,Z axes. I would also like the result was plotted in 3D graph. Someone can help me ? Thanks a lot ! Thank, useful! But I don't know how to do plot of results ...


title: Smoothing spline with dependent coordinates in R

text: I have a series of patients in whom I have measured the value of a certain blood quantity at several time points. However, the time points vary a lot and the number of measurements range between 2 and 5. I would like to fit a smoothing spline to a graph were I have each patient's blood values on the Y-axis and time on the X-axis. I would like to have a graph like the one shown below ![enter image description here](http://i.stack.imgur.com/ZVLsf.png) I am using R, and with the `smooth.spline` function I can have a smoothing spline with all data-points used, but I think I should consider each patient's data points as dependent and hence I cannot put them just to a two vector. How could I produce a smoothing as shown above?


title: Repeated multiple regression for LASSO significance testing

text: I'm currently working on a statistical modelling problem in biology. We have cellular measurements of proteins in every cell in a tissue, and I'm using regression analysis to see if a given protein is affected by the protein content in nearest neighbours. Overall, I'd like a general linear model of the form $Y = \beta X + B$ where $Y$ is the matrix of measurements of multiple proteins across cells, and $X$ is likewise but the average over nearest neighbours of a given cell. However, in each cell we believe there is a network that internally affects protein levels, so I've used LASSO (we believe the underlying matrix is sparse, and have many proteins so many noise variables) within each cell to regress a given protein off all the others, leaving a matrix of residuals $\tilde{Y}$ (which naturally has a covariance of near 0), to give a linear model $\tilde{Y} = \beta X + B$. As we believe the interactions in this model to also be sparse, I've been using LASSO regression with the significance testing found in               Bühlmann, Peter, Markus Kalisch, and Lukas Meier. "High-Dimensional Statistics with a View Toward Applications in Biology." Annual Review of Statistics and Its Application 1  (2014): 255-278.      However, as $\Sigma = 0$ for $\tilde{Y}$ and we want to select different predictor variables for each response, I've been repeating multiple regression for each response variable, rather than performing general linear regression on all the response variables at once. The question I'd like to ask is, is this a valid approach? Are there any obvious flaws with the method? Any other comments generally? Thanks in advance.


title: Long-run cointegration and I(2) variables

text: I'd like to run a test for long-run cointegration of two time series: CO2 emissions/capita and energy use/capita in China. My data includes 23 observations for years 1990-2012. Now as far as I know I have to have _I(1)_ series to test for their cointegration. Nevertheless, unit root tests point out that these two variables might not be a result of _I(1)_ processes, but rather of _I(2)_. Have a look at the test statistics of ADF: ![enter image description here](http://i.stack.imgur.com/KiHbL.png) And ERS test: ![enter image description here](http://i.stack.imgur.com/X7aYZ.png) The autocorrelation functions of first differences possess some autoregressivity: ![enter image description here](http://i.stack.imgur.com/CkMRg.png) ![enter image description here](http://i.stack.imgur.com/jFYQE.png) Now the stats of ERS test for first differences are a bit ambigous, but it may be due to misspecification of the model. It seems that these time series are indeed _I(2)_. About long-run cointegration testing; I know that for instance ARDL test allows to have mixed I(0) and I(1) series included, but its statistics are not valid for I(2) variables. Same for Johansen test, its statistics are not valied for I(2). Am I right to this point? If I take first differences of the series and test for their cointegration, then the results are not valid regarding long-run cointegration. And finally, how do tests work in cases including structural breaks? For instance have a look at the plots of both original series: ![enter image description here](http://i.stack.imgur.com/gGtbQ.png) ![enter image description here](http://i.stack.imgur.com/Gnr8s.png) There are obviously breaks in the series. How do they affect the outcome?


title: Given time series data, how to model the frequency of someone changes his job?

text: I am given a time series data vector (ordered by months and years),which contains only `0`s and `1`s. `1` s represent a person changes his job at a particular a month. **Questions:** What model can i use to determine model how frequently this person change his job ? In addition, this model should be able to predict the probability of this person changing his in the next 6 months. A poisson process ? (I have studied poisson process before however I have no idea when and how to apply it). Any assumptions that data need to meet before applying the poisson process ? Would love to gather more information on how to model something like this. Thanks


title: Linear and semi-log regression model

text: Is this equation: $$\log{(y)} = a + bx$$ semi-log or log-linear mode (or it is the same thing)? I have two models: linear (1) and semi-log (2). The values of $R^{2}$, adjusted $R^{2}$, and Standard Error are:   * Linear: $R^{2} = 0.6780,~\mathrm{adj.}~R^{2} = 0.6513,~~\mathrm{SE}=94.101$    * Semi-log: $R^{2} = 0.5803,~\mathrm{adj.}~R^{2} = 0.5455,~~\mathrm{SE}=0.5493$ How to interpret this values especially from the second model?


title: Gibbs Sampling Inserting Some Known Predictors

text: Imagine you would like to use a simple Gibbs sampling to resample from a joint probability distribution which is difficult to model (but you know all the conditionals $Pr\left(X_i|X_1,...,X_{i-1},X_{i+1},...,X_N\right)$, for $i=1,...,N$, hence the Gibbs sampling). At the end, you would have a set of many drawn vectors of variables $\left(x_1^{(k)},x_2^{(k)},...,x_N^{(k)}\right)$, for $k=1,...,K$. Imagine now you would like to assign a vector of variables $\left(y_1^{(k)},y_2^{(k)},...,y_M^{(k)}\right)$ to each of the already drawn vectors. You have some conditional probabilities like $Pr\left(Y_j|Y_1,...,Y_{j-1},Y_{j+1},...,Y_M,X_1,...,X_N\right)$. Again, modelling joint distributions is difficult. Do you think it can work a Gibbs sampling with these conditionals by inserting sequentially the previously drawn $x_i^{(k)}$ and, after discarding some steps, collecting the sampled $y_j^{(k)}$? Otherwise, do you have a strategy to apply in these cases? Thanks!


title: What are the limitations of ANOVA vis a vis multiple regrssion analysis"

text: It seems that multiple regression is a better way to do analysis. For example,the ANOVA simply gives t-value (or F-value for two or more variables) and does not indicate magnitude of effect. Contrary-wise, the regression analysis indicates both the t-value(s) as well as magnitude of regression coefficients. Does that mean we should discard ANOVA completely.


title: Designing a probability function with a desired kurtosis and skewness

text: Is it possible to design a probability function P(x) such that the distribution will have a specified kurtosis, K, and skewness, S? I am not otherwise interested in any other properties of this distribution. P(x) = f(K,S) ...f is some function of kurtosis, skewness.


title: How to calculate type II error between glm fits in R?

text: I'm using a GLM with logistic link function to try to predict Y (0 or 1) as a function of a ton of predictor variables (A, B, C, etc.). Some of the predictor variables (A*, B*, C*, etc.) have been shown in other studies to be significant predictors. I want to show essentially that Y is unrelated to all of the other predictor variables, and I thought the simplest way to do this would be to run the full model (Y ~ .) and the null model (Y ~ A* + B* + C* + ...), and then use anova() to compare the two and show that they aren't different (i.e. have equal predictive power). However, anova() only outputs p-values (type I error), but I need a type II error rate here (since I want to show that the models are the same, I need a false negative rate for that). Any ideas on how to approach this?


title: How to interpret this binomial logistic regression output in R

text: I have run my logistic regression model to find out whether the gender of a test name is a predictor of the gender of the word assigned to the test name.               q2.glm=glm(item.gender ~ test.gender, data=q2data, family="binomial")     summary(q2.glm)      R has returned the following coefficients:               Coefficients:                      Estimate Std. Error z value Pr(>|z|)           (Intercept)  0.05365    0.12384   0.433  0.66484            test.gender  0.48794    0.17836   2.736  0.00623 **      I understand the P value for test.gender means I can reject the null hypothesis that the gender of the test name has no influence on the gender of the word assigned to it, however I don't understand how to interpret the other values. Please can someone explain what the other values mean? What does the intercept category mean? many thanks - I have a big project due soon and am struggling to understand it!


title: Basic question regarding n of 1 study. Four weeks intervention

text: I initially will apologise for any lack of clarity with regards to this question. I'm posing this on behalf of my sibling who is in the biomedical sciences. They have done an intervention on a patient. They are measuring the effect of the intervention by looking at test scores (6 questions), and repeating this test over four weeks. The data looks like this Week 1: 4/6 Week 2: 5/6 Week 3: 5/6 Week 4: 6/6 My sibling wants to see if the increase in the test scores is significant or not. I didn't really know how to go about this. Any help would be greatly appreciated


title: Meaning of Chi-squared test output

text: I have an interaction table with plants species (column) and the number of time each plant has been found to live in each environment (row). I'd like to know if there is a differential "preference" for plants species to live in different environment. I guess I should run a chi squared test. I'd need some help to make sense of the output. The greatest is the residual, the biggest is the "preference" of one plant species to its environment. Is it correct? The p.value gives the proportion of times that we would find such extreme pattern if all species would have no "preference" of environment. Is it correct? What if I'd like to have such p.value for each plant species. Should I run a chi-squared test for each plant species independently? Or is there some kind of post-hoc test of the chi-squared test?


title: Cronbach's Alpha with missing data

text: In `R`, there are several packages that calculate Cronbach's alpha. If the data contain no missing values, all packages converge to the same value. With missing data, listwise deletion is a possible way to go (the only option in SPSS or packages `MBESS` and `psy` btw). However, listwise deletion might lead to dropping a lot of data and therefore something like pairwise deletion might seem more appealing in some situations (let's say data are MCAR).   Some of these packages have something like `na.rm=T`, but then the results differ, and they differ more the more missings a data set contains. The reason for this is the calculation of the variance-covariance matrix $S$. For example, the package `coefficientalpha` employs a robust $S$ matrix naturally leading to different results. However, other packages use _pairwise_ $S$ matrices and still the results differ.   Let's make a specific example. The formula for alpha is $$\alpha = \frac{p}{p-1}(1-\frac{\sum_{i=1}^p \sigma_{y_i}^2}{\sigma_x^2})$$ with the number of items $p$, the variance of the $i$th item $\sigma_{y_i}^2$ and the variance of the total score $\sigma_x^2$. For the total variance in the denominator, the package `ltm` uses something like `var(rowSums(dat, na.rm=T))`, whereas the package `psych` uses something like `sum(cov(dat, use = "pairwise"))`. These two estimates are the same for complete data, but differ with missing data. (Package `epicalc` seems to use a third way of calculating the variance(s), but I have not yet mastered their code.) **What is the reason for different calculations of $\sigma_x^2$ and when should I use which one?**               > require("MASS", "psych", "ltm")     > n <- 200     > it <- 10     > V <- matrix(.4, ncol=it, nrow=it)     > diag(V) <- 1     > dat <- mvrnorm(n, rep(0,it), V)   # complete data     > m1 <- matrix(rbinom(n*it,1,.1), nrow=n, ncol=it)   # 10% missings     > dat[m1==1] <- NA   # missing data     > alpha(as.data.frame(dat), na.rm=T)$total[[1]]   # psych package       [1] 0.8595489     > cronbach.alpha(dat, na.rm=T)$alpha   # ltm package     [1] 0.8105867


title: What is a high dimensional multivariate data set?

text: I'm new to data analysis and data mining. Often in the papers I'm reading, they use the term "high dimensional multivariate data set." Currently, my task is to detect an outlier and visualize the same from a large complex data set. But how does one find out whether I have a multivariate high dimensional data set or not?


title: How to value comparative Likert results

text: ![enter image description here](http://i.stack.imgur.com/DXzxU.png) I have absolutely no knowledge of statistics. My question is how to make sense of some evaluation results I received. The attached graph represents my results of a Likert question compared to my department, division, and college. My question is, how should I interpret these results? While I have more students who "strongly agree" on this question, I have fewer who merely "agree." What interpretations of my performance are authorized here? Am I performing better? Similarly? Is there some empirically valid way of making assertions about this data? There are twelve other such graphs, and I have to write a summary of the findings. But it is hard for me to know what kind of claims I can make based on the data.


title: Help me calculate how many people will come to my wedding! Can I attribute a percentage to each person and add them?

text: I am planning my wedding. I wish to estimate how many people will come to my wedding. I have created a list of people and the chance that they will attend in percentage. For example `Dad 100% Mom 100% Bob 50% Marc 10% Jacob 25% Joseph 30%` I have a list of about 230 people with percentages. How can I estimate how many people will attend my wedding? Can I simply add up the percentages and divide it by 100? For example, if I invite 10 people with each a 10% chance of coming, I can expect 1 person? If I invite 20 people with a 50% chance of coming, can I expect 10 people?


title: Constant forecasts in SPSS

text: I have weekly data for the last four years. I am using SPSS to do forecasting. I am getting a constant value in the forecast period. What could be the reason behind it? Is it due to defining weekly dates? I am not able to figure it out. I have defined the week dates using following code: `DATE WEEK 1 52` ![enter image description here](http://i.stack.imgur.com/2Wqhy.jpg) Thanks in advance.


title: What does it take to learn a Model?

text: In Machine Learning, there are various models. I tried to learn few probabilistic models of Machine Learning. I read the theory, worked on problems and tested my results and could analyze my data according to the need like classification, etc. I am a Computational Linguist and I am trying to post many questions on the models(Naive Bayes and HMM) I learnt and going to learn (CRF,MEM). I was thinking if there are any set of questions that I should be able to answer after learning a model, upon answering whom I may say I know a model? Is there any set of questionairre available in the web or in any text book( preferably with the answers)? If any one of the esteemed members of the group may kindly let me know? Regards, Subhabrata Banerjee.


title: $E(\overline{x})$ vs $E(x)$?

text: If I have a problem down to $$k(E(x^2) - E(x)^2) = E(\overline{x}^2) - E(\overline{x})^2$$ ...can I say $k = 1$? $X$ has an unbiased mean and a variance. Each $x$ in $X$ is independent. Or is it more complex than that? If it's not clear, $\overline{x}$ signifies the average (or mean or whatever, I always thought they were the same until this homework assignment...)


title: Does A independent of B, and B correlated with C imply that C is independent of A?

text: Assume 3 random variables, $A, B, C$. If $A \perp B$, but $Cov(B,C) \neq 0$, can we say anything about $Cov(A,C)$? I think it can either be 0 or not 0, but it seems like there should be more general results than that.


title: "Zero-inflated continuous covariates", Can they cause problems in logistic regression?

text: I pose a very similar question to this, although I felt the advice given does not apply to my particular situation; I am using logistic regression models for an animal habitat occupancy study, and all the predictor variables I am interested in contain >50% zeros (although they have a decent range of values in the higher percentiles). Can this cause bias or influence how I should interpret the estimated coefficients? A 2-stage analysis, as suggested in the linked question, doesn't seem to make sense because all the predictors share this zero-inflated distribution. Thanks for any insights **EDIT Clarifications suggested by Peter Flom;** Sample size ~ 500 (300 "0"s, 200 "1"s) There are 5 IV's; a typical five-number summary looks like this; min= 0.000 lower= 0.000 median=0.000 upper= 0.289 max= 16.887 Also, Mean= 0.468, SD= 1.467 correlations between the 5 IV's all absolute r < 0.3 The IV's are hectares of specific habitat types. Every sample has >0 hectare(s) for at least 1 of the IV's. An example run of the model in R;                   Call:         glm(formula = use ~ x.1 + x.2 + x.3 + x.4, family = binomial,              data = mydata)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -2.2338  -0.9312  -0.8679   1.3231   1.6432            Coefficients:             Estimate Std. Error z value Pr(>|z|)         (Intercept) -0.78412    0.12814  -6.119 9.41e-10 ***     x.1          0.19866    0.06366   3.121  0.00181 **      x.2          0.06956    0.02618   2.657  0.00788 **      x.3          0.05238    0.02265   2.313  0.02074 *       x.4         -0.09995    0.13814  -0.724  0.46935         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for binomial family taken to be 1)              Null deviance: 634.10  on 473  degrees of freedom     Residual deviance: 611.18  on 469  degrees of freedom     AIC: 621.18          Number of Fisher Scoring iterations: 4


title: How to visualize the effect of a regression parameter

text: I am arguing that I can control error vs. coverage by modifying a certain parameter. After running an experiment with leave-many-out validation I have a set of errors along with the parameter value for each error. I want to visualize how I can trade coverage for accuracy by placing restrictions on my parameter. If I didn't have a parameter to modify I could use a regression error characteristic curve. These curves plot the error tolerance on the $x$-axis vs. the fraction of test points within that tolerance on the $y$-axis. Drawing several REC curves using different values of my parameter leads to a plot with many overcrossing lines which is hard to interpret. I've also found some work on "ROC surfaces" but I don't know how widespread it is. Right now I'm using a scatter of mean error vs. test points retained with points colored by the parameter value: ![enter image description here](http://i.stack.imgur.com/OUQCp.png) It's clear to me what is going on in my figure (e.g. I can set the parameter to 10^2 and get 80% coverage at mean error 200). However, I still feel like I've run into a very common problem that has standardized solutions. Is there is another visualization technique I should be looking into?


title: Combining Multiple Principal Components?

text: Background: I am an evolutionary biologist interested in how morphological variation within a population of individuals contributes to variation in their competitive ability. To get at this in the tadpoles that I work on, I want to create a 'morphological index' that incorporates four different traits that are thought to be involved in competitive ability. Given that two of the traits are continuous traits and two are ordinal, I used a PCA using a correlation matrix to create this index. I then want to see how this morphological index influences various measures of competitive ability using multiple regression. Problem: PC 1 varies in how much variation it explains depending on what population I look at. It ranges from 34% to 64% and it probably averages about 45%, which is lower than I would like. At the least, it seems that PC 2 is explaining quite a bit of the variation. I am relatively new to statistics, so I am wondering if there is a simple way to combine PC 1 and PC 2 into a single morphological index? I'll apologize right now for my ignorance of statistics that likely shines through in this post. I'm happy to provide any important details that I likely omitted. Thanks for your help!


title: Methodology for finding source of MANCOVA interaction

text: I am trying to find out the best method for determining the source of a significant MANCOVA interaction. I have 17 traits measured for four locations. Is there an alternate to univariate ANCOVAs investigating each trait independently? These are body measurements which aren't independent of one another, so it would be good to investigate all traits at once.


title: Multivariate Bayesian formula

text: I got there example graphs bishop's PRML (8.2.1)   1\. **a <\- c -> b**   $$ p(a,b,c) = p(a|c)p(b|c)p(c) --(1)\\\ p(a,b) = \sum_c p(a|c)p(b|c)p(c) --(2) $$ **Q1** : Can I use a new graph to represent the p(a,b) in formula (2) ? 2. **a->c->b** $$ p(a,b,c) = p(a)p(c|a)p(b|c)--(3)\\\ p(a,b) = p(a)\sum_c{p(c|a)p(b|c)} = p(a)p(b|a) --(4) $$   **Q2** : $$\sum_c{p(c|a)p(b|c)} = \sum_c{p(b,c|a)} ? why$$


title: Area Under ROC Curve for Multiple Classes

text: I am working with a highly class-skewed three class classification problem. The class percentages are A = 1.8%, B = 17.5% and C = 80.7%. According to this paper, the following definition of multi-class AUC is insensitive to class distributions and therefore, I am using it:![enter image description here](http://i.stack.imgur.com/oqdEN.png) Now I am testing 10 different classfiers on the dataset and only one of them has $AUC_{total} > 0.5$. Am I right in assuming that $AUC_{total} = 0.5$ for a random classifier? If yes, are 9 of the classifiers worse than a random classifier (on this data) just because they have a lower $AUC_{total}$? If yes again, can you suggest possible reasons for such poor performance and what can be done. P.S.: Even the best performing classifier has $AUC_{total}$ of about 0.6.


title: Mathematics behind classification and regression trees

text: Can anyone help explain some of the mathematics behind classification in CART? I'm looking to understand how two main stages happen. For instance I trained a CART classifier on a dataset and used a testing dataset to mark its predictive performance but:   1. How is the initial root of the tree chosen?   2. Why and how is each branch formed? My dataset being 400 thousand records with 15 columns and 23 classes achieves a 100% accuracy from a confusion matrix, I use 10-fold crossvalidation on the dataset. I would be really greatful if anyone could help explain the stages of CART classification?


title: Prediction interval with measurement error

text: My data is $(Y,X)_i,i=1...N$ and i want to fit a regression model to relate $Y$ and $X$. It is known that $X$ and $Y$ are erroneous measurement of true variable $(U,V)$, but the standard errors $(\sigma_X,\sigma_Y)_i, i=1...N$ are known to me. I also know that the standard error is correlated with the true value (and thus the observed value as well). i.e. $cov(\sigma_X,U)\neq 0$ Now, if I want to do inference on the model parameters, I know I need to use measurement error models (Deming, Ripley's MLFR etc.) However, if I only want point prediction, I can just use OLS naively. (Fuller (1987), p.74) My question is, what if I want to get a prediction interval? I shouldn't use the naive OLS prediction interval because there is measurement error. So what should i do? If answer is not available I will appreciate pointer to reference too. Thanks in advance.


title: How to report permutation ANOVA results in a scientific publication?

text: I have been using `lmp` function from lmPerm package to calculate differences between treatment levels (2-ANOVA). Now, I am wondering how to report `lmp` results in a scientific publication? I am aiming for a biological or geochemical journal, but have been using APA style in reporting statistics previously. I am not sure whether this is the desired way either, but at least their rules seem well defined.


title: Difference between theoretical model and fit

text: I have two sets of experimental data for exponential decays of isotopes and I calculated the ratio of the two isotopes as a function of time. The theoretical model and the fit of the experimental data are respectively : $1.04\text{e}^{-0.000856008t}$ $1.03154\text{e}^{-0.000457549t}$ The two curves are near of each other (I'm only interested in comparing the [0-110] x-range) but is there a way to calculate an error between the two curves (like a standard deviation)? I need to determine if the data are normally distributed so I thought getting the standard deviation from the theory/fit would be a correct way... Any ideas?


title: Use of nlme for mixed models of repeated measures with two classifying factors and a covariate

text: Suppose we have a mixed model for repeated measures with response being $y$ and a covariate being $x$(say, age). There are two classification factors $a$ and $b$. Say, 5 psychological patients($a$) in 8 conditions($b$) are measured. I want to know how I get the estimates of intercepts and slopes (and their standard errors) of the random effects model when:   1. $a$ fixed, $b$ random    2. $b$ fixed, $a$ random    3. Both random How should the appropriate code look like? Besides, can you tell me something about the meaning of the following error message?                 nlminb problem, convergence error code = 1       message = iteration limit reached without convergence (10)


title: Which search range for determining SVM optimal parameters?

text: I am using SVM for classification and I am trying to determine the optimal parameters for linear and RBF kernels. For the linear kernel I use cross- validated parameter selection to determine C and for the RBF kernel I use grid search to determine C and gamma. I have 20 (numeric) features and 70 training examples that should be classified into 7 classes. My question is which search range I should use for determining the optimal values for the C and gamma parameters.


title: Collaborative filtering through matrix factorization with logistic loss function

text: Consider collaborative filtering problem. We have matrix $M$ of size #users * #items. $M_{i,j} = 1$ if user i likes item j, $M_{i,j} = 0$ if user i dislikes item j and $M_{i,j}=?$ if there is no data about (i,j) pair. We want to predict $M_{i,j}$ for future user, item pairs. Standard collaborative filtering approach is to represent M as product of 2 matrices $U \times V$ such that $||M - U \times V||_2$ is minimal (e.g. minimizing mean square error for known elements of $M$). To me logistic loss function seems more suitable, why are all algorithms using MSE?


title: Likelihood-ratio test in the big picture

text: I am trying to understand the big picture of statistical tests. That is how to choose the correct statistical test under different situations. I found this table where choosing a statistical test is based on the types of the dependent and independent variables. But I didn't find the likelihood-ratio test in this table, so I wonder which category in this table the likelihood-ratio test belongs to or where is the place of the likelihood-ratio test in the big picture of all statistical tests.


title: Statistical significance in canonical correlation analysis

text: I do canonical correlation analysis between two multivariate datasets $X$ and $Y$. For each pair of canonical variates (x-y pair) I get the canonical correlation coefficient. How can I test its statistical significance?


title: Confusion related to intractability in topic models

text: ![enter image description here](http://i.stack.imgur.com/tb2va.png) I was reading this paper related to topic models. I am a bit confused why the marginal likelihood is not tractable and how converting the graphical model into the new one actually helps. First I don't understand how the coupling between $\beta$ and $\theta$ result in making it intractable. It points to a reference. But I cannot access the article (Dickey 1983). So I cannot understand why it is intractable. can anyone please provide me some information?


title: Beyond the normal distribution: what if a particular distribution can't be assumed

text: In a random sample of 150 community college students, the mean number of hours spent studying per week is 11.7 hours and the standard deviation is 4 hours. Without assuming anything about the distribution of the number of hours community college students study per week, at least what percentage (approximately) of the students study between 5.3 and 18.1 hours per week?


title: Issue in graph construction

text: I have a symbolic representation of time series obtained from SAX toolbox. I was wondering if it is possible to construct a graph where each node represents a unique symbol and the edges represent the transition providing that there is no transition to itself. For ex, let the time series T of n=20 data points be represented as `T=[1 1 2 1 2 1 3 1 1 1 2 2 3 3 3 1 1 2 3 3 1]'` where number of alphabets used to symbolize = 3 and they are (1,2,3). I have combined co-occuring symbols together so, the compressed time series becomes               T' = [1 2 1 2 1 3 1 2 3 1 2 3 1]'      In order to construct a graph (Esp. fuzzy cognitive map ) with fuzzy membership values from T' where the nodes will be (1,2,3), there will be an edge from $Node_i$ to $Node_j$ and it will have a weight $W_{ji}$. How do I find the weights ? I do not know which theory to search for this kind of problem and so if there are any ideas as to what can form the weights. Thank you


title: Error implemented a specific factorial design in Minitab (default generators)

text: I'm intending to implement the following factorial design I wish to obtain this alias structure - with a 2^(5-2) factorial design.               I+ACE+BDE+ABCD     A+CE+ABDE+BCD     B+ABCE+DE+ACD     C+AE+BCDE+ABD     D+ACDE+BE+ABC     E+AC+BD+ABCDE     AB+BCE+ADE+CD     AD+CDE+ABE+BC      I've obtained this (the alias structure) by hand, using the design generators I=ACE and I=BDE. I've controlled my alias structure and confirmed it is the correct one. However, using Stat->DOE->Create factorial design-> 2 level design (specify generators) and using 5 factors And subsequently using Default generators: ACE, BDE Gives the error "Block are aliased with the mean". Why doesn't this work?


Suppose a business is operating at the three-sigma quality level. If projects have an average improvement rate of 50% annually, how many years will it take to achieve six-sigma quality?


title: Interpretation of the final cluster centers (cluster analysis)

text: I have a question concerning the interpretation of the final cluster centers. I performed a cluster analysis based on a pca (the variables are based on a five point Likert-scale). I got the following result for one factor:                       cluster_1        cluster_2       cluster_3      cluster_4       cluster_5               0,31            0,39               -0,82          0,63            0,35      Is this factor also interesting for the description of cluster 1, 2 and 5? Or should I only mention its influence for the clusters 3 and 4? Thanks a lot!


title: Basic questions about discrete time survival analysis

text: I am attempting to carry out a discrete time survival analysis using a logistic regression model, and I'm not sure I completely understand the process. I would greatly appreciate assistance with a few basic questions. Here is the set up: I'm looking at membership in a group within a five year time window. Each member has a monthly record of membership for each month that member is in the group. I'm considering all the members whose membership began during the five year window (to avoid "left censorship" issues with members who joined earlier). Each record will be indexed by time, with time one being the month the member joined. So, a member who stays for two and a half years will have thirty monthly records, numbered from one up to thirty. Each record will also be given a binary variable, which will have a value of one for the last month of membership, and zero otherwise; a value of one for the binary variable marks the event that the member has left the group. For each member whose membership continues beyond the five year analysis window, all the binary variable values will be zero (these are the right-censored individuals in the survival analysis). So, the logistic regression model is built to predict the values of the binary event variable. So far, so good. One of the typical ways to evaluate a binary predictive model is to measure the lift on a holdout sample. For the logistic regression model I have built to predict the membership ending event, I have computed the lift on a holdout data set with a five to one ratio of non-events to events. I ranked the predicted values into deciles. The decile with the highest predicted values contains seventy percent ones, a lift of more than four. The first two deciles combined contain sixty-five percent of all the ones in the holdout. In certain contexts this would be considered a fairly decent predictive model, but I wonder whether it's good enough to carry out a survival analysis. Let $h[j,k]$ be the hazard function for individual $j$ in month $k$, and let $S[j,k]$ be the probability that individual $j$ survives through month $k$. Here are my fundamental questions:   1. Is the discrete hazard function, $h[j,k]$, the conditional probability of non-survival (leaving the group) in each month?   2. Are the predicted values from the logistic regression model estimates of the hazard function? (i.e., is $h[j,k]$ equal to the model predicted value for individual $j$ in month $k$, or does something more need to be done to obtain hazard function estimates?)   3. Is the probability of survival up to month q for individual $j$ equal to the product of one minus the hazard function from month one up to $q$, that is, does $S[j,q] = (1 - h[j,1]) \cdot (1 - h[j,2]) \cdot \ldots \cdot (1 - h[j,q])$?   4. Is the mean value of $S[j,k]$ over all individuals $j$ for each time $k$ a reasonable estimate of the overall population mean survival probability?   5. Should a plot of the overall population mean survival probability by month resemble the monthly Kaplan-Meier graph? If the answer to any of these questions is no, then I have a serious misunderstanding, and could really use some assistance / explanation. Also, is there any rule of thumb for how good the binary predictive model needs to be in order to produce an accurate survival profile? Thanks! \-- TMK --


title: Model Matrix and GBM

text: I am using the GBM package to create a prediction model. Can GBM handle categorical predictors, or do I need to create a model matrix before I input ?


title: Higher r-squared value on test data than training data?

text: I am trying to create a linear regression model. I split my data into training and testing data, and built a model. The R-squared value on the training data is 0.840. Then I ran the model on the test data. When I calculate R-squared, I get 0.982:               y.predicted <- predict(lm1, newdata=test)     y.actual <- donation_test$yval     errors <- (y.actual - y.predicted)     1 - sum(errors^2)/sum(y.actual^2)     [1] 0.9823576        What I am doing wrong? It seems very unlikely that my model fits my test data better than my training data.


title: Logistic regression and odds ratio

text: If the odds ratio is greater than one with an insignificant p value for a variable in logistic regression should the variable be kept in the model? Can I select the variable with odds close to 1? Question 2: In the validation and testing dataset do not have same mean how to handle it?


title: How to compare means of repeated measures between groups in SPSS

text: I'm having trouble regarding which test to run on my data in SPSS. I have two groups (condition1 and condition2) which is also my independent variable (categorical). Both groups were asked to rate their brand attitude (dependent variable 1) and purchase intention (dependent variable 2) at two points in time (before and after an intervention). So I have repeated measures at only two points in time. Hypothesis: the increase in brand attitude and purchase intention is greater for people in group1 (condition1) than in group2 (condition2). I already did a paired t test for each group to see the change in brand attitude and purchase intention, but now I want to compare the groups and see if the difference in means is significant. Can you tell me which test to run?


title: ARIMA model selection in SAS

text: I want to fit an ARIMA model to my data, however I don't know how to select the best ARIMA model using the `proc ARIMA` command. First of all I know I want to transform my data with the cubic root transformation.               data robberies_cubic_root;     set robberies;     robberiescubicroot=robberies**(1/3);     run;      Now I also want to take the first difference of the transformed data, thus the `var=robberiescubicroot(1)`. At this point I just want to fit ARIMA models where p=0,...,4, and q=0,...,4. (i.e, I want to fit ARIMA(0,1,0),...,(0,1,4), (1,1,0),...,(1,1,4),(2,1,0),...,(2,1,4)....all the way to (4,1,4).) And then choose the best ARIMA model. However when I use `estimate p=(0,1,2,3,4) q=(0,1,2,3,4)`, it only seems to fit a ARIMA(4,1,4) model for me.               proc arima data=robberies_cubic_root;     identify var=robberiescubicroot(1);     estimate p=(0,1,2,3,4) q=(0,1,2,3,4) plot;     forecast lead=5 id=month out=forecast;     run;


title: Canonical reference for Naive Bayes classifier

text: I'm looking for a reference for the Naive bayse classifier to put in my work. Not sure what I'm missing but a scholar search didn't yield any meaningful results. any idea ? **update** : I meant for a reference to put in the bibliography...


title: intuitive understanding of max-min in optimization

text: I have problem in understanding max-min in optimization. I know what does it mean to maximize obejctive function or minimize it, I don't really understand what does it mean to have objective function like max min at the same time. I put example below: ![enter image description here](http://i.stack.imgur.com/kFL9C.png) Can anybody explain it to me ?


title: Logistic regression with mutually exclusive predictor variables

text: I am working on a logistic regression approach to predict the clinical status of patients (No disease vs Disease). I already have quite strong evidence indicating that the number of genes hit by a certain kind of variation is significant in predicting disease. Among the genes that I am considering, one of them (gene A) apparently plays a central role: I would like to investigate if the effect caused by the number of genes involved is driven by "gene A + others" combinations or "any multiple genes combination". I developed a model with the 3 independent dummy variables _geneA.alone_ (T if a sample has variation only in gene A), _geneA.plus.others_ (T if a sample has variation in gene A and other gene(s)), and _only.others_ (T if a sample has variation in one or more genes, but not in gene A).   I performed two nested model comparisons using the _anova()_ function in R:               anova(glm(Status~geneA.alone,...),glm(Status~geneA.alone+geneA.plus.others,...),      test="Chisq")          anova(glm(Status~geneA.alone+geneA.plus.others,...),        glm(Status~geneA.alone+geneA.plus.others+only.others,...),test="Chisq")      In both cases I observe an improvement of the model fit, with a significant p-value (in the order of e-5 / e-6). So, my first conclusion was that, when multiple genes are involved, it is not only the combinations "geneA+other(s)" that are important in determining disease, but also those not involving gene A. However, I am really doubting that this is the right way to investigate it:   * Does it make sense to have three mutually exclusive variables as predictors? Can they be considered as independent variables? My feeling is that the significant p-values are a natural consequence of the fact that, by using mutually exclusive predictors, I deliberately subtract part of the information from a variable, and shift it to another...hance the significance.   * Would it be more appropriate to use only the variables _number.of.genes_ (numeric) and _geneA.involved_ (T/F), and compare the model `Status~n.of.genes+geneA.involved` with `Status~n.of.genes*geneA.involved`? My thought, in this case, is that if the second model results in a significantly better fit, it may be a demonstration that the _n.of.genes_ is important only if geneA is involved...does it make sense? Any help would be greatly appreciated! Thank you!


title: Weakly or strongly identified parameters?

text: I was reading this paper. Authors talk about "weakly identified parameters" or "strongly identified parameters" continuously. Rigorously, what does it mean when one say a parameter is strongly or weakly identified? Thanks in advance!


title: Inferring an unmeasured value

text: Suppose we have a table of finishing times for an number of 100 metre races where each competition has a different mix of entrants. We suspect that some races are slower than others due to, say, a headwind, but this wasn't measured. How do I find the underlying time for each runner and the correction for each race? I have devised my own crude technique to solve this problem. I assume               Tobs(e,c) = Tund(e) + Tcor(c) + E(e,c)          where           Tobs(e,c) = observed time of entrant e in competition c     Tund(e) = underlying time of entrant e     Tcor(c) = time correction for competition c     E(e,c) = error matrix           sum of all Tcor = 0     sum of all rows in E(e,c) = 0     sum of all columns in E(e,c) = 0           I then iteratively choose values for Tcor until the above relations are satisfied. Is there a better way? Using R? Please be gentle. As you might have guessed by now, I don't much much about stats or maths. Thanks!


title: Comparing smoke and hormones level in two groups of people. Which test?

text: I'm totally not a pro in statistics so i hope the question is not too simple! I have two groups of patients. One group are smokers, the second are non- smokers. I want to compare the level of some hormones in the two groups to see if they are influenced by smoking. On many articles I found researchers using Mann- Whitney test to see if the means are significantly different among the two groups. But reading around t-test i'm not sure anymore about which test should I use... unpaired t-test? unpaired t-test with Welch correlation? Mann- Whitney? one-tailed, two-tailed? Can someone explain me which test to use is this situation and why? (in a way comprehensible by normal people?) Many Thanks!!!! **Some more info:** The two groups are N = 102 (smokers) and 194 (non- smokers), and i red that for N > 30 unpaired t-test and mann-whitney should be similar but applying three tests i got three different p values:   * Unpaired t-test: p = 0.0245   * Unpaired t-test with Welch correction: p = 0.0318   * Mann-Whithney: p = 0.0620 As you can see p values varies a lot and with Mann-Whitney significance is even lost (>0.05)! **Update 2:** Also correlation between number of daily cigarettes and hormone level gives strange results: Spearman test gives r = 0.1786 and p = 0.0725, while Pearson gives r = 0.2472 and p = 0.0123. So again with a parametric test the result is significative while with a non-parametric is not!!! I suppose that in this case non-parametric is more correct since we cannot be sure to have the right number of cigarettes, especially since this is a retrospective study.


title: Response reviewer on a repeated measures logistic regression

text: I got a paper rejected at a conference with the following comment: > The used method for the analysis is correct, but I find it strange that the > authors evaluate the results as a 4-condition experiment instead of the 2x2 > experiment they set up. This makes their analysis of the interaction effect > a bit clumsy, and they end up testing all possible combinations of > conditions (which means they should have done a Bonferroni correction). A > factorial analysis would have been much more straightforward. As I do not have the possibility to respond, can anybody explain what the reviewer meant with: > the authors evaluate the results as a 4-condition experiment instead of the > 2x2 experiment they set up > > A factorial analysis would have been much more straightforward. Can anybody explain what a factorial analysis is, and how a 4-condition experiment differ from a 2x2 experiment. My experiment consist of 2x2 within-subjects design. For example A (yes/no) x B (yes/no), where AxB (no/no) is the control condition. I chose to perform a repeated measures logistic regression analyses as I have a:   * Within-subjects design   * A four level categorical IV (yes/yes, yes/no, no/yes, no/no)   * A two level categorical DV (yes/no) To test significance I contrasted the categories against a reference. I varied the reference to asses differences between categories. I don't see the "clumsiness" the reviewer implies, as this is the way I've learned it using a regression.


title: Given two biased samplers, create a less biased sampler

text: Suppose you have two biased sampling algorithms, and you don't know the bias of the algorithms. Can you combine the algorithms somehow and reduce the bias? Partly my motivation is discretizing stochastic differential equations. Say you sample using Euler's method and the Milstein method. Both are biased, but can you combine the algorithms somehow are get a less biased sample? Another example: suppose you have two MCMC algorithms that are having a hard time converging to the stationary distribution, can you combine the MCMC algos and produce a better algo? Further motivation: In machine learning, one almost always gets a further improvement by averaging predictions with (ideally) different biases. Now, this is because in ML often our target is a expected value ($E[Y |X]$ or $P(Y|X) = E[ 1_{Y} | X]$), and averages are a great approximation to expected values. Obviously averaging samples does not work, but it there an analogous approach?


title: Wild cluster bootstrap seems really simple. Too simple. Am I missing someting?

text: I've been dealing with the problem of how to construct confidence intervals on penalized spline estimators in the presence of cluster-wise auto-correlation and heteroskedasticity. My previous thread on this is here: Penalized spline confidence intervals based on cluster-sandwich VCV One of the commenters suggested using a wild bootstrap. This is a bit tricky in the presence of clustered data, but this paper by Cameron et al. at UC Davis proposes what seems to be a relatively simple approach that accommodates the clustered context: http://ideas.repec.org/a/tpr/restat/v90y2008i3p414-427.html Basically, the residual vector for each cluster has an equal chance of being left alone, or being multiplied by -1. Thats it. No resampling. Over many iterations, you get many different combinations of cluster-wise residual vectors, but your combinations are restricted to the set of residual vectors that you already have, potentially multiplied by -1. (I'm aware that there are alternative things that you might multiply your residual vector by, summarized on wikipedia: http://en.wikipedia.org/wiki/Bootstrapping_%28statistics%29#Wild_bootstrap) Basically this all seems too simple. Am I missing some complicating wrinkle? Given that I'm programming my own estimators in R, I need to program my own bootstraps as well. I want to make sure I get the procedure right before I go and code. Thanks in advance.


title: How to interpret factor scores saved using the regression method in SPSS?

text: I have 30 attribute variables with response on a 1 to 5 scale (e.g., extremely important, very important, important, not that important, etc.). I ran factor analysis in SPSS and extracted 6 factors and saved the values using the regression approach. Now I want to understand them. I can make sense of the 6 factors themselves, but not the individual scores.   * What is the intuition behind the factor saved scores for each respondent?   * What is the best way to display them visually?


title: Lasso for GEE model

text: 1. Can a LASSO be applied for predictor selection in a logistic GEE (generalized estimating equations) model for longitudinal data?   2. Is there an implementation of LASSO for a logistic GEE model for longitudinal data in either R or MATLAB? Thanks!


title: Panel data with two time factors

text: I've got some data that come from a economic experiment. I'm trying to fit a model using plm() function (plm package). I'm following a very good tutorial (Panel data econometrics in R: The plm package). The only problem is the function allow the following data structure:               individuals / time / x1...xp / y       and I've got the following structure:               individuals / time 1 (session) / time 2 (round) / ....      I was reading the function gls() (nlme package), so may be I can explore the autocorrelation structure, add it to the model and treat session and rounds as dummy variables. But what about individuals? They are the same within each session, but change between sessions. So individuals is a random factor, isn't it?. Any advice will be appreciated


title: Raw Prices vs. Daily Return vs. Price Ratio - ADF unit-root test

text: As the title suggests, I'd like to ask what the difference between using Raw Prices vs. Daily Return vs. Price Ratio in the Augmented Dickey-Fuller unit- root test. The context is that I am trying to find cointegrated pairs of stocks from a given time series in R. Take a look at these three different resources, all using different metrics (Raw, Daily and Ratio) for the ADF test:   1\. Raw Prices   2\. Daily Returns   3\. Price Ratio (A UBS Research monograph, sorry couldn't find link but here is an excerpt): > Adjusted Dickey-Fuller (ADF) test – a test of stationarity of the price > ratio with null hypothesis of a unit root. Also worth noting that #3 doesn't perform linear regression (OLS), performing the ADF test on the price ratio itself, while #1 and #2 do it on the spread (residuals) of OLS. Any reason for this? So if anyone could shed some light on the issue, that would be fantastic. Maybe it's just a matter of preference, but hopefully there are reasons for choosing them, with solid evidence backing it up to boot!


title: How to interpret GARCH volatility forecast?

text: I am trying to model the volatility of gold using GARCH models and here is the forecast that I obtained. ![Gold Volatility](http://i.stack.imgur.com/ImrB8.jpg) The plot seems to indicate that the volatility decreases as time passes but I was expecting to see something like an increased volatility over time. Can anyone explain why this is the case? Please don't pay attention to the x-axis don't know how to change that, it's like that by default, the real date is 2009-2013.


title: significant difference in preference

text: I have the following data                    X black brown     1 women    25    20     2   men    41    44      How can I see if there is a significant difference in preference between men and women using R