title: Contribution of each covariate to a single prediction in a logistic regression model

text: Say, for example, that we have a logistic regression model which outputs the probability that a patient will develop a particular disease based on many covariates. We can get an idea of the magnitude and direction of the effect of each covariate in general by examining the coefficients of the model and considering the change in the odds-ratio. What if we want to know for a single patient what his or her biggest risk factors / the biggest factors in his or her favor. I'm particularly interested in those which the patient could actually do something about. What is the best way to do this? The way I'm currently considering is captured in the following R code (taken from this thread):               #Derived from Collett 'Modelling Binary Data' 2nd Edition p.98-99     #Need reproducible "random" numbers.     seed <- 67          num.students <- 1000     which.student <- 1          #Generate data frame with made-up data from students:     set.seed(seed) #reset seed     v1 <- rbinom(num.students,1,0.7)     v2 <- rnorm(length(v1),0.7,0.3)     v3 <- rpois(length(v1),1)          #Create df representing students     students <- data.frame(         intercept = rep(1,length(v1)),         outcome = v1,         score1 = v2,         score2 = v3      )      print(head(students))          predict.and.append <- function(input){         #Create a vanilla logistic model as a function of score1 and score2         data.model <- glm(outcome ~ score1 + score2, data=input, family=binomial)              #Calculate predictions and SE.fit with the R package's internal method         # These are in logits.         predictions <- as.data.frame(predict(data.model, se.fit=TRUE,      type='link'))              predictions$actual <- input$outcome         predictions$lower <- plogis(predictions$fit - 1.96 * predictions$se.fit)             predictions$prediction <- plogis(predictions$fit)             predictions$upper <- plogis(predictions$fit + 1.96 * predictions$se.fit)                   return (list(data.model, predictions))     }          output <- predict.and.append(students)          data.model <- output[[1]]          #summary(data.model)          #Export vcov matrix      model.vcov <- vcov(data.model)          # Now our goal is to reproduce 'predictions' and the se.fit manually using the      vcov matrix     this.student.predictors <- as.matrix(students[which.student,c(1,3,4)])          #Prediction:     this.student.prediction <- sum(this.student.predictors * coef(data.model))     square.student <- t(this.student.predictors) %*% this.student.predictors     se.student <- sqrt(sum(model.vcov * square.student))          manual.prediction <- data.frame(lower = plogis(this.student.prediction -    1.96*se.student),          prediction = plogis(this.student.prediction),          upper = plogis(this.student.prediction + 1.96*se.student))          print("Data preview:")     print(head(students))     print(paste("Point estimate of the outcome probability for student",     which.student,"(2.5%, point prediction, 97.5%) by Collett's procedure:"))     manual.prediction     print(paste("Point estimate of the outcome probability for student",     which.student,"(2.5%, point prediction, 97.5%) by R's predict.glm:"))         print(output[[2]][which.student,c('lower','prediction','upper')])      I'm considering looking additionally at               this.student.prediction.list <- this.student.predictors * coef(data.model)      and trying to get the information out of the individual addends of the sum which is the probability estimate, but I'm not sure how to do it. I could look at   * Which variables make the largest absolute contribution to the probability estimate and take those to be the largest risk factors.   * Which variables differ by the greatest amount from their mean proportion, i.e. see what proportion each variable contributes to the probability estimate on average and see which variables differ from this proportion by the greatest amount in this particular observation    * A combination thereof: weight the absolute difference between mean proportion and observed proportion by the mean proportion and take those variables with the greatest weighted values Which of these make the most sense? Would any of these approaches be a reasonable way to answer the question? Additionally, I'd like to know how I could get confidence intervals for the additive contributions of individual covariates to the probability estimate.


title: significant difference between coefficients

text: What is the correct way to test for significant differences in parameter estimates in the following case: I have a dependent variable (height for age) and a independent variable (household state pension benefits). These are observed in two years. To see if there is an endogeneity problem I want to include pension benefits in year two in the regression of year 1. The equation I want to estimate is thus: $$y_1 = b_0 + b_1 * x_1 + b_2 * x_2 + e$$ I then want to test if $b_1$ is significantly different from $b_2$ (in stead different from zero). I found this previous question, but I am not sure if this is exactly the same problem. Also, is there a standard way of implementing this in R? This is my first post, so sorry if I didn't post correctly, please feel free to edit.


title: How can I determine if values are 'significantly' similar / close?

text: I have an experimental value and I want to find out which of the theoretical values are 'significantly' close to the experimental value. Say for example my experimental value is:               3000      and my theoretical values are:               2500,3005,3300,2750,2900,3100,2995,3000,3010,3150,2850      How would I find out which of my theoretical values are 'significantly' close to my experimental values. I was assuming I could calculate this using standard deviation but I would have to calculate the mean, but my mean isn't necessarily going to be at my experimental value. What process/statistical test would be best appropriate to find the values significantly close to my experimental value?


title: Non-Normal residuals after Box-Cox

text: I applied Box-Cox on a regression model with heteroscedasticity. Although the transformed model looked better than the original, it still showed signs of heteroscedasticity. What else can I do after Box-Cox to almost eliminate heteroscedasticity? The original model showed a skewness (regarding normality) of 1.4 where the transformed showed a skewness in the residuals of 0.8. What other transformations can I use and/or should I look more closely to the original observations and take out observations that are outliers?


title: ARIMAX model or ARDL?

text: I would like to study the impact the advertising of a product on its sales (weekly data for 5 years). As the final aim is to forecast what would be the impact on sales of a change in the advertising presence in the media, I was considering including also other variables (such as competitors prices, macroeconomics variables,…) As I know that advertising has a long lasting effect and that, for example, at time t I would be only able to know prices at t-1 I need to include lags. Would ARIMAX model be appropriate? Shall I just try to use a trend+seasonality+lagged variables? Thank you.


title: Classification of Huge number of classes

text: I have a dataset of samples belonging to >100 classes. I want to classify and/or cluster these classes. I have the following questions: 1) Is one classifier efficient for such problem? or one classifier for each one/subset of classes? (From my point of view: the efficient solution is to discover the features discriminate each class from all others and solve the problem as 1-to-all classification problem. Any suggestion on that?) 2) About 60% of these classes have 1 or 2 samples at maximum!. How can I create new samples from these 1-sample classes. Do you think any of SMOTE (synthetic minority oversampling technique) techniques are workable in this case. Regards,


title: How to find the expected distance between two uniformly distributed points?

text: If I were to define the coordinates $(X_{1},Y_{1})$ and $(X_{2},Y_{2})$ where $$X_{1},X_{2} \sim \text{Unif}(0,30)\text{ and }Y_{1},Y_{2} \sim \text{Unif}(0,40).$$ How would I find the expected value of the distance between them?      I was thinking, since the distance is calculated by $\sqrt{(X_{1}-X_{2})^{2} + (Y_{1}-Y_{2})^{2}})$ would the expected value just be $(1/30 + 1/30)^2 + (1/40+1/40)^2$?


title: What does "principled" mean, as in "principled Bayesian analysis"?

text: I am generally curious what the term "principled" means. It was used in the title of an unpublished manuscript, "Combining Computer Models in a Principled Bayesian Analysis". In addition, Zhang 2004 says "There are several advantages to a principled Bayesian analysis" here: ![enter image description here](http://i.stack.imgur.com/QybsP.png) Otherwise, I find relatively few google hits for the phrase. What does the term "principled" mean in this context? Zhang, 2004 Causal Inference with Instrumental Variables. Ch. 8 in Gelman and Meng, Applied Bayesian Modeling and Causal Inference from Incomplete-Data Perspectives. Wiley.


title: p-values of Mann-Whitney U test identical for raw and log-transformed data

text: I am new to Stats and came across this problem while running my analyses on SPSS which I cannot explain. How is it that even after transforming my data by logging it, it still has the same p-value as the raw data set that was not transformed? Sorry for sounding really ignorant here! Thank you in advance for your help.


title: Transforming data: correlate regardless of distribution

text: Is there a way to correlate data regardless of distribution? I know the Choleksy transformation is used for normally distributed data, but is there a general method that applies to any case? To include more information:   * I want to be able to decide any distribution shape for any variable.   * I would like to correlate 2 or more variables.   * In the end, I want the variables scaled (this can easily be done afterwards) I now realized, after running some tests, that the answer probably lies in reordering the data (I do not care about maintaining "pairs", so I am _not_ talking about observed x and y _within_ the same person, so these values can be shuffled around). I found out that even for a small n=10, reshuffling 100.000 times has a correlation of 0.3 included with a small tolerance.               a <- runif(10)     b <- rnorm(10)          cors <- 1:100000     for (i in 1:100000) {         cors[i] <- cor(a, sample(b, length(b))     }     Output:     > sort(cors)[89981]     > 0.2999949      Would there be an easy way to get to this reordering quickly (just finding out that this correlation exists, takes 7 seconds to run)? Also, this might not be the best way to go for the 3 variable case.


title: Is there an R equivalent of SAS PROC FREQ?

text: Does anyone know of an R equivalent to SAS `PROC FREQ`? I am trying to generate summary descriptive statistics for multiple variables at once.


title: What is the preliminary scale in M-estimates for regression

text: I'm trying to get aquainted with robust regression methods and there's something about M-estimators that I don't understand. In "Robust statistics" (Maronna, Martin, Yohai) it is said that if both our x's and y's are random, the M-estimator of regression is given as a solution to $$ \sum \rho (\frac{r_i (\hat{\beta})}{\hat{\sigma}}) = min $$ where $\rho$ is a bounded $\rho$ function and $\hat{\sigma}$ is a "high breakdown preliminary scale". What is that "preliminary scale" supposed to be? How do I understand preliminary in that context?


title: Student t as mixture of gaussian

text: How to show the following: Using the student t-distribution with $k > 0$ degrees of freedom, location parameter $l$ and scale parameter $s$ having density $ \frac{\Gamma ((k+1)/2)}{\Gamma ((k/2) \sqrt{k \pi s^2}} \\{ 1 + k^{-1}( \frac{x-l}{s} )\\}^{-(k+1)/2}$ How to show that $t$-distribution can be written as a mixture of Gaussian distributions by by letting $X$~ $N(\mu,\sigma^2)$, $\tau = 1/\sigma^2$~$\Gamma(\alpha,\beta)$ and integrating the joint density $f(x,\tau|\mu)$ to get marginal density $f(x|\mu)$. What are the parameters of the resulting $t$-distribution, as functions of $\mu,\alpha,\beta$ I got lost in calculus. By integrating the joint conditional density with the Gamma distribution. Please help! Thanks a lot!


title: Software for working with exponential random graphs

text: I am looking for software packages for working with **exponential random graph models** (fitting/generating them and sampling from the graph distributions). I have only found two packages so far, both using R: ergm/statnet and RSiena, with the former being much more popular. Are there any other packages available? I'm particularly interested in methods that are both published/documented and have an available implementation to play with.


title: What book is recommendable to start learning statistics using R at the same time?

text: # Books to Learn Statistics using R ## What exactly is the book I'm looking for. What I am looking for is a book that teaches you statistics while using R to give you hands-on experience and thus end up helping you learn R together. I've seen on amazon many books that attempts to do that, but not with R. Examples are Minitab and SAS. ## Are the R Book and Statistical Computing an option? - _Still not answered_. The R Book and Statistical Computing: An Introduction to Data Analysis using S-Plus seems viable, but a reader opinion here would be helpful and welcome. ## How the book relate to statistics courses? To be even more precise on what I was looking for, consider these two courses learning outcomes on statistics from a math department at the university Im currently a student: Intermediate Statistics and Probability & Statistics, that is, I'm looking in a book a normal statistics course going to intermediate level but rather than just board and paper having you learning and using R instead. That also means I am looking for a book that assume I want to learn statistics from the beginning. ## This book is for researchers too. I am also a software engineer researcher, but I guess the current situation where you are found with mountains of data and want to learn statistics to go on writing code to automate that is pretty much applicable to many other fields. That means I'm am not interested on learning every single detail of every single property for every single curve, but am more concerned on making sense of data for my research domain, although I would not mind if the book wanted to go deep on that. As a final motivation, I find myself reading scientific papers in different sort of communities that claim results based on statistical inference while there is no readable proof if the statistics assumptions/constraints are being violated or not. A R book that is not much about statistics won't ensure I am not following up on this practice, which is also why I decided looking for a book that is akin to a statistics course using R rather than playing around with a overview book. ## Related questions in Cross Validated.   * What books provide an overview of computational statistics as it applies to computer science? \- Differs that the question looks for an overview while this is to learn statistics using R.   * Open source statistical books gives a list of open source (open books) available online.  ## Answers and feedback for this question. ### @Julie Suggested books were few I already come across but are an example that unfortunately doesn't suits me: Introductory Statistics with R, Using R for Introductory Statistics, Statistics: An Introduction using R are few of the books that I already looked on amazon but are about an statistics overview or make assumptions that requires previous statistics knowledge. The problem with overview books is mostly about not calling attention to the assumptions, constraints and provide enough explanation to result in make sense of the information. If you believe there is no book that could fit on this needing as well or think the R book or the Statistical Computing: An Introduction to Data Analysis using S-Plus would fit this, I would also appreciate this type of answer. ### @Christopher Aden Introduction to Probability and Statistics Using R seems to be the closest one but still broad general to what I was looking for. _What I was expecting for is a book such asDavid S. Moore, The Basics of Statistics because:_   * It covers all statistics subjects.    * It uses two tools, miniTab and other to give hands-on learning on the just explained method.   * It very much highlight assumptions and constraints. This is very important for a researcher who has not taken a in depth statistics course and want to use statistics. Hardly overview books will cover them, which is dangerous for researchers.      * You can see the book table of contents here. Notice how the focus is statistics and the tool usage is to improve understanding and get the student to know how to use tools to do the statistics after learning in an easier way. Its not about the tool, its about statistics!  _I want exactly the same thing, but using R._ ### @Gregory Demin It uses R as pedagogy examples, assumes you want to learn statistics and best of all, it is open source. Unfortunately, does not cover ANOVA nor ANCOVA, or more advanced subjects. ### @Peter Ellis Good suggestion for a textbook that covers what is wanted in this question. ## Books in the asker opinion that answer the question. @Peter Ellis and @Gregory Demin. ## Collection of R Books on Amazon Amazon discussion about R books for different students background may be found here. ## Video Lectures teaching Statistics using R Google Tech Talks from 2007 that also motivated this question and covers more about Data Mining rather than statistics but using R together here.


title: Hidden Markov model implementation

text: I've been looking for the internet explanation about the hidden Markov model (HMM) and its implementation but I think there's not quite good explanation. Does anyone know and guide me through this like tutorial and its implementation. There will be feedback and not do this as charity. So I hope there will be someone interested with this subject please contact me very further so we can talk soon. Thank you very much.


title: What location parameter is modelled by robust regression?

text: There is quite some number of ways how to robustly fit a linear regression model, e.g. using M-estimation based on Tukey's biweight loss or on Huber's loss, see e.g. Wikipedia. I got two questions about the two mentioned cases:   1. Is it correct that these are modelling the mean (more precisely _the center_ ) of the conditional response distribution as long as the error distribution is symmetric?   2. What location parameter are they modelling in general? Is it a specially weighted average of the response variable?


title: Abusing Linear Models under Multicollinearity: Simulation for 'realistic' movement of predictors

text: I have a reasonable understanding of why multicollinearity is a problem is regression models, along the lines of this excellent post. To summarise my understanding, for a regression model of $y = \alpha + \beta_1x + \beta_2z$ (where $x$ and $z$ are correlated), beta coefficient estimates (as well as being unstable) are difficult to interpret, as a situation where you might increase $z$ without increasing $x$ is unlikely to occur, and not supported by the data. I understand multicollinearity is less harmful to purely _predictive_ as opposed to explanatory or descriptive models. I'm interested in another interpretation: _If I decided to increase $z$, and let $x$ vary as it pleases in reaction, what would I see happen to $y$, accounting for the fact that $x$ is likely to move with $z$, and also have it's own effect?_ In other words, accepting the causal interpretation that $x$ and $z$ both cause $y$, and are themselves correlated to some extent (.7 say), how would all three variables move if $z$ is (linearly) increased by some amount? I've tried to model this sort of thing before, fitting $y = \alpha + \beta_1x + \beta_2z$ (model 1), and $x = \alpha + \beta_1z$ (model 2). Hypothetical increased $z$ values are produced, and resulting $x$ values are predicted with model 2. The hypothetical $x$ and $z$ values are used to predict $y$ using model 1. However this feels very unsatisfactory, complicated simulations are required to capture uncertainty (I used `sim` in `arm`). Additionally, my gut tells me that apart from being painfully inelegant, it's a bad idea for other reasons I can't put my finger on.   * Is such an 'observational'/conditional-when-I-feel-like-it interpretation possible?   * Does anyone know of a better method for this interpretation?   * Can anyone recommend a paper or `R` package along these lines?   * Is the above multi-model mess at-all valid? I'm aware that a model along the lines of $y = \alpha + \beta_1z$ would yield a similar answer to the two-stage mess above, but would lose information in $x$. I understand that these ideas are similar to structural equation modelling, but apart from having scant knowledge of SEM, I'm yet to find an `R` package which allows flexibly extending these models with different link functions for proportional odds models, etc.


title: How do MCMC methods allow the estimation of the posterior distribution in this example?

text: I am reading a book example (diagram from p10) in which a person scores 9/10 on which we assumed a uniform prior. The posterior distribution could be easily worked out analytically, but the book gives an example of estimating the posterior distribution by MCMC methods. I can see from the bottom part of the diagram that about about 95% of samples had a theta of between 0.59 and 0.98, where theta represents 'ability' in the sense of accuracy rate on these tests. However, I don't understand how a theta value is obtained for each sample in the chain, and what the relationship is between a sample in a chain and the proceeding sample. The book also does not explain why three chains are chosen or where the initial values for these chains come from. ![enter image description here](http://i.stack.imgur.com/qGmjR.png)


title: How to compare two distance matrices?

text: Suppose that I have two distance matrices for the same set of items. By a distance matrix I mean a square matrix whose (i,j)th entry holds the distance (in terms of cosine similarity) between ith and jth items. The ith and jth items are the same items in both matrices. Such a situation might happen when we gather information about a set of items from two different sources. What I want to do is to compare these two distance matrices. Whether they are similar or not with respect to the distance relations among items. One idea is to find the correlation between the elements of the two matrices (only the upper triangular elements since these are symmetric matrices). This can be done by Mantel's test. Another thing comes to mind is to build clusters out of these distance matrices and compare the resultant two clusterings. But does this give any additional information over Mantel's test? Can we use other measures to understand the similarity between these two distance matrices or the above two methods are enough? To make things concrete, the items are documents. One set consists of wikipedia documents written in English and the other set consists of the same documents written in another language (say German). Documents are encoded as tf*idf vectors and their similarity values are measured using cosine similarity. So one distance matrix hold the similarities of the English documents and the other one holds the similarities of the German documents Thanks Ahmet


title: 22 word frequency Variable has been selected for predicting whether tomorrow is up or down day (binary target), how to interpret

text: First of all, I appreciate lot of creator of this web and any attention on this question. I asked related question here: 52 variables after backward variable selection on logistic regression on 160 variable at beginning, whether it is illusion or good modeling . The Statistic Guru @PeterFlom pointed out my naive choice to dichotomize the positive value of volatility whose self-correlated characteristics also violated the assumption that residuals should be independent. I totally agree with his opinion and thank him again for the precious time and patience from him. I really learned lot. Now I changed my target and modelling way, little bit long description again below:   1. Target Variable: Binary Target of day T + 1: if the close price is higher than the open price of EURUSD, I mark it as 1, and if the close price is lower than the open price of EURUSD< I mark it as 0. Little assessment: this time I believe here is not violation of assumption about independent residual like the one we had in volatility. We know volatility has kind of self-correlated characteristics, But in my humble opinion, here is no similar characteristics for the up and down of the Price.The distribution of the target looks like this pic: ![enter image description here](http://i.stack.imgur.com/kbmp9.jpg)   1. Predictor Variable: again, frequency of 160 trading related key words from famous trading online forum. For example, the word “trader” in 2008-01-01 we found it occurred for 50 times and we found in 2008-01-01 we have totally 500 posts on that forum. Therefore, we kind of normalize the appearance of "trader" by dividing 50 of "trader" by 500 of "number of posts". We got 0.1 here for 2008-01-01. We have 160 words dictionary and the training data is from 2008-01-01 to 2011-12-31. Each row has 160 column with number indicating the frequency of specific words. That is all. Nothing variable of price, volatility here. Kind of pure sentiment data based on frequency counting. @gung I normalize the frequency by the reply number of day T, do you think I should still use log on the normalized frequency which is usually less than 1?    2. The goal and story line: whether we could use the combination of frequency of words could predict whether tomorrow Day T + 1 is up or down. The reason why I do not use linear regression to predict exact log return or percentage return is because the most words normalized frequency distributed in this way and not easy to transform them in to regular distribution shape, like normal distribution, that is why I give up on predicting the EXACT change of price, Please have a look on the distribution of "limit" and "Correction": ![enter image description here](http://i.stack.imgur.com/m8Gkv.png) ![enter image description here](http://i.stack.imgur.com/OAXgv.png)   3. Again: Software I use: SAS Enterprese Miner 3 in SAS 9.1.3. Here is definition and details of modeling: A. I select logistic regression B. No data partition because I only have 1000days as my training data so I select to use complete data which means no sampling process and only training data. C. I select the backward as we want to see what are those significant predictor out of 160, here is some shortage if using forward. D. There is no validation data set so I chose the Crossvalidation-error as my ultimate criteria. Significant level for variable to stay is 0.05. 5.Result: thankful for the sas easy usage, I got the following graph directly from the the result and modeling manager. This time I got 22 word frequency variables left as predictor, below is the result: ![enter image description here](http://i.stack.imgur.com/fZ8TQ.png) ![enter image description here](http://i.stack.imgur.com/vzb2a.png) ![enter image description here](http://i.stack.imgur.com/gxV0w.png)   1. my naive interpretation about the three picture above. The goodness of fit is looking good. We got almost 76% accuracy of predicting on whether tomorrow is up or not and much better than the 50% baseline as random trader. The cumulative response chart also looks smooth. And we got 22 variable and it looks less scaring but still kind of heavy degree of freedom?    2. @PeterFlom, should I use proc glmmix here as you recommended over the shortage of backward selection on variables?   3. I pick up the top 10% day with highest score(predicted probability) from model above and do LONG, and bottom 10% day with lowest score from model above and do SHORT. The percentage gain is 50% from these 200 trading days out of 1000 days as dataset. Is this a good validation for model?


title: Library for using plates in bayesian networks

text: In my Bayesian network there are plenty of repetitive variables leading to the use of plates(http://en.wikipedia.org/wiki/Plate_notation). I do not want the exponential space complexity in implementing the bayesian network. But the way traditional libraries(like dlib on c++) offers bayesian network, it is always exponential. Is it possible to actually use the plates and reduce the space complexity? If yes, please help with some libraries that does this(Preferably c++)


title: Is it true that decision trees can handle dependent features very well?

text: I have once heard that decision trees can handle dependent features very well, compared to e.g. Naive Bayes classifiers. However, I can't find any scientific source for that. So my question: is it true and if yes, is there any source?


title: Interactions in multiple regression

text: I would like to predict BMR or Basal Metabolic Rate (energy expended by the body while awake) using the predictors age, weight, height and gender. Please suggest the interactions that I should be considering here. I am using spss to do this regression. Shouldn't the significance of interactions between continuous variables and that with gender be tested?


title: Reparametrisation of a model when an interaction is significant to facilitate the interpretation

text: It is admitted that it is complex to interpret main effects when they are involved in an interaction. Lets take a regular linear model, with two categorical 2 level variables A and B who are interacting together. The model can be written:               lm(Response ~ A + B + A:B)      Lets call I the intercept, a2 the estimate associated to level 2 of variable A, b2 the estimate associated to level 2 of variable B, and c22 the interaction estimate. For me, there is 2 cases:   1. The interaction estimate is small compared to main effects estimates. In that case, I see no issue in interpreting main effects and then specify the interaction is significant but small and explain what it changes.    2. The interaction estimate is big compared to main effect, and may even change the direction of effects. In that case, it is very difficult to interpret estimates, and it is better to look the effect of B for each level of A, or the effect of A for each level of B, depending on our underlying understanding of the interaction between A and B.  We focus on case 2. Let's say B is a modifier of A effect, and I want to look what is the effect of A in each level of B. To look what the effect of A is for each level of B, one would sum the estimates:   * Effect of A for level 1 of B is a2    * Effect of A for level 2 of B is a2+c22 My issue is that when doing this, we do not know the confidence interval of the effect of A for level 2 of B (i.e. of a2+c22). Therefore, I wondered whether a good solution to interpret interactions in that case would be to reparametrise the model as:               lm(Response ~ B + A:B)      This returns two estimates for the interaction instead of one (c21=a2 and c22=a2+previous_c22), so that we directly have the estimated effect of A for each level of B, and moreover we have their confidence intervals. Is it a good solution? If it is, why don't people do that? Thanks!


title: Can someone explain (in simplest words possible) minimum statistics, ancillary,completeness and MUVE

text: So, I understand what sufficient statistics is but I am unable to extend those ideas to minimal statistics, completeness, ancillary and eventually connect it with UMVU estimators? Can someone explain it to me or point me the (easy to understand) material? Thanks


title: Proving conditional independence from mutual independence

text: Suppose I have a Bayesian inference problem where $E$ is the result of a chemistry experiment. The chemistry experiment consists of a controlled sequential digestion of one end of the protein using a specific enzyme that has known statistical properties. The sequence of the protein is unknown at the beginning of the experiment. The data collected is the number of reactions and the times in which they occur. The goal is to infer sequence information using the enzyme. $P\left(x_{1},\, x_{2},\,\ldots,\, x_{n}|E\right)$ = $\dfrac{P\left(x_{1},\, x_{2},\,\ldots,\, x_{n}\right)P\left(E|x_{1},\, x_{2},\,\ldots,\, x_{n}\right)}{P\left(E\right)}$ Suppose also that the discrete random variables $X_{1},\, X_{2},\,\ldots,\, X_{n}$ of the posterior are known to be mutually independent a priori. These random variables of the posterior give sequence information (the 20 amino acids) for the first n positions in the protein. This Bayesian inference calculation is performed for all amino acid sequences of length n. The marginal distributions are then taken to infer the unknown sequence. How do I prove the conditional independence of the posterior random variables given the evidence: $P\left(x_{1},\, x_{2},\,\ldots,\, x_{n}|E\right)$ = $P\left(x_{1}|E\right)P\left(x_{2}|E\right)\cdots P\left(x_{n}|E\right)$ .


title: Prove $\mathbb E(X|X>0)=\mathbb E(X)/\mathbb P(X>0)$

text: The following is a lemma I read $$\mathbb E(X|X>0)=\frac{\mathbb E(X)}{\mathbb P(X>0)}.$$ Could anyone tell me how to prove it, please? Thank you! **Update** Using the conditional expectation formula I got the following $$\mathbb E(X|X>0)=\sum_x x\mathbb P(X=x|X>0)=\sum_x x \frac{\mathbb P(X=x, X>0)}{\mathbb P(X>0)}.$$ What is $\mathbb P(X=x, X>0)$, please?


title: Normal distribution is better than Poisson for large sample set?

text: I'm at a beginner level, so please bear with me. This is a call center use case. For every week, certain number of calls are received. The average is about 20. This seemed like a Poisson distribution to me (rate per interval). I took about 100 data points. I tried to answer the question "what is the probability of receiving more than 5,10,20,30 etc. calls per week". I computed the corresponding Poisson results. I created another column and answered the same question using Normal distribution. I then generated the probability by directly querying from the database (ie # of weeks where there were more than x calls/total number of weeks). I found something interesting, the actual data followed Normal distribution and not Poisson distribution. Why is that ? I assumed that "number of calls" is discrete & so it should be Poisson. The question I'm getting at is   1. When is it appropriate to use normal distribution when you have discrete variables ?   2. Where is Poisson distribution appropriate ?    3. What is the distribution to use to model discrete variables ? I use Libre Office and SQL Server.


title: Collinearity in time dummies, fixed effect regression

text: ![enter image description here](http://i.stack.imgur.com/7xQQ3.png) I am running a fixed effect panel regression with 81 groups x 20 periods, so approx 1620 (unbalanced) observations. I use the following to create dummies:               *create timedummy       tabulate refper, generate(refdummy)        which to visual inspection looks fine. But when running the `xtreg`, it removes 6 time dummies due to collinearity. Can anyone explain what is happening here?


title: Need help selecting test for nominal IV with categorical DP

text: I am going to attempt to give you as much background as possible.   My **independent variable** is nominal and describes the number of days after which an insect was exposed to new conditions (0,1,3,5,7).   The **dependent variable** is the insect's condition after the experimental period (Healthy, Sick, Dead, Diapause).   From exploratory analysis, there are very clear connections between Days and Condition. When I graph the percentage of each outcome against days, there are very obvious correlations. I.e., Dead is a very rare outcome at Day=0 but become steadily more common as Day increases.   What I am missing at this moment, however, is the appropriate stat test to demonstrate this trend is statistically relevant. Any help would be greatly appreciated. I am using SPSS. 1) All insects are individuals and were recorded as such, so yes I can say #42 has a discrete outcome.   2) It's a 0,1,3,5,7 scale because those are the only days such manipulations were possible for a combination of logistical and physiological reasons.   3) No outcomes are reversible and all observations were made on the same day REGARDLESS of what day the environmental changes were initiated, 28 days after adult emergence. So it not that the 7 day data set was older than the 0 day data set, it's just that the environmental change happened at day 7 instead of day 0.


title: conditional probability, chemical reactor residence times

text: The following scenario is given in Dekking's "A Modern Introduction to Probability and Statistics" to illustrate conditional probability: > "Consider a continuously stirred reactor vessel where a chemical reaction > takes place. On one side fluid or gas flows in, and it flows out on the > other side. > > Residence time: How long particles stay inside the vessel. Rt={residence > time > t) = e^-t, assuming continuous stirring > > ![enter image description here](http://i.stack.imgur.com/6Qjk5.png) **Please explain how the author equates numerators in the second and third expressions, P(R4 ∩ R3) and P(R4)**.


title: Matrix Factorization Recommendation Systems with Only "Like Ratings"

text: I'm trying to build a recommendation system, but I only have data on what my user's have "liked" i.e. all non-missing data has the same numeric value. Is it possible for me to using matrix factorization methods without actually having "ratings?" (multiple numeric values for user ratings rather than just an indicator that the user has "liked" the item?) If so, how?


title: Tuning the inflation parameter in mcl

text: I read on the mcl documentation that the inflation parameter can be used to tune the granularity of the clusters. I am not very familiar with graph theory. What is the granularity of the clusters? Can I use this inflation parameter to choose to keep the clusters for which the nodes are very interconnected? (If, for instance, I have 4 nodes a,b,c and d that are all interconnected and a 5th node e connected to a only, then is it possible to play with the inflation parameter to get two clusters {a,b,c,d} and {e}?)


title: Radial distribution to cartesian coordinates

text: I'm trying to go from a radial distribution in polar coordinates to cartesian coordinates. I'm using Matlab and have tried pol2cart, but are having some trouble getting a matrix with probabilities. The distribution is a log-normal along the radius, homogenous along theta.What I'd like to have is a 2D matrix with the probabilities so I could plot a heat map or similar, getting something like this: ![enter image description here](http://i.stack.imgur.com/2CNkT.jpg) Thank you in advance, Cheers!


title: Variance-covariance matrix for ridge regression with stochastic $\lambda$

text: In ridge regression with design matrix $X$, outcomes $y$, fixed regularization parameter $\lambda$, and errors $\epsilon\sim\mathcal{N}(0, \sigma^2I)$, the computations for the ridge regression coefficients $\hat\beta$ (aka the solution to $\arg\min_b \big[(y-Xb)'(y-Xb) + \lambda b'b\big]$) and their variance-covariance matrix $var(\hat\beta)$ are: $$ \begin{align*} M &:= (X'X + \lambda I)^{-1}X' \\\ \hat\beta &= My \\\ var(\hat\beta) &= \sigma^2MM' \end{align*} $$ The computation on $var(\hat\beta)$ relies on $\lambda$ being treated as a constant value, which makes $M$ constant. Thus, we can apply the identity $var(My) = Mvar(y)M'$ and the fact that $var(y) = \sigma^2I$ to derive $var(\hat\beta)$. However when $\lambda$ is selected using $X$ and $y$ (for my application via cross-validation), $\lambda$ and $M$ become stochastic. How can $var(\hat\beta)$ be updated to deal with a stochastic $\lambda$ selected via cross-validation?


title: How much time is reasonable for training Restricted Boltzmann Machine?

text: I have implemented Binary RBM in Matlab. I am using 60000 images as an input to train RBM. It takes approximately 11.3 minutes. I used tic and toc functions to evaluate above mentioned elapsed time. Is that, looks reasonable time? or it is slow?


title: Cluster analysis of open-question responses

text: I am currently analysing a 1100-observation dataset of open-question responses. I was wondering whether I could really do a segmentation based on these responses. What I've done is so far is conducted a qualitative coding of each response, and then produced a binomial variable for each of the codes, like this:               Response                     Donuts     Cakes        I like donuts and cakes      1          1     I like cakes                 0          1     Hello world                  0          0      The number of binomial variables in my real dataset is 22, with some categories containing 9 cases all the way up to 322. Do you think it would be valid to conduct some sort of cluster analysis based on these binomial variables? What method should I better use (LCA, cluster, principal components - I haven't done much of those, so need for-dummies advice)? Additionally, my kinda philosophical concern was - are zeroes in this table really zeroes, or are they **missings**? It might be different from explicit tick-box questions, where the absence of a tick **really means** that a person's not having this attribute. Here, he might have expressed that he likes cakes, but **maybe he also likes donuts, just forgot to write about it** , i.e. he would have told us if we asked explicitly. I presume, in the latter case, there isn't much clustering to be done? Thanks a lot for your help guys - I am not a statistician, but this will be a big step forward to my dissertation.


title: Forecasting the impressions a specific ad campaign will recieve

text: I am currently working on forecasting the impressions an ad campaign will get given certain constraints specified by the ad provider (eg. only target men in California on websites with certain capabilities). I have forecast the total views all ad spaces will get using an ARIMA model in R. I now need to find some way to include the constraints. The data I have lists all of the properties of the views each ad space gets and I have processed this also. How can I determine the correlations between the different constraints using the data I have and somehow hit my ARIMA forecast with a time-dependent dampener to forecast the actual impressions a certain ad campaign will get given the constraints?


title: What is the difference between simple linear model and loess model?

text: Can somebody explain to me the difference between linear model and loess model in statistics? I need to explain this to non-math people.


title: Geometric mean of two numbers

text: I need to find the geometric mean of 0 and 600. Would it just be the square root of 600? I know that the formula is $$ \bar x_{\mathrm{GM}} = \left(x_1 x_2 \cdots x_n\right)^{1/n}, $$ but in this case one of my values is 0, so I'm not sure.


title: Confusion related to derivation of a formula in gaussian process

text: I was reading this paper related to online Gaussian processes. ![enter image description here](http://i.stack.imgur.com/3MWl9.png) ![enter image description here](http://i.stack.imgur.com/liJHA.png) I didn't get how equation set 2 was derived. Any suggestions?


title: Mutual information with a Dirac delta type pdf

text: What does the $MI(X,Y)$ convey about $Y$, when one of the probability distributions, $X$ is trivial and has all the probability concentrated at a single point?


title: unstandardized coefficients where mean-centered X variables have +/- values from 0

text: > **Possible Duplicate:**   >  How to interpret regression coefficients for a variable with takes positive > and negative values? I'm performing negative binomial regression and one of my centered-mean X variables (with 200 observations) ranges between -0.2 and +0.2 around the 0. Do regression coefficients only account for a positive increase from 0 on its effect on the dependent, or can/do/how they account for a negative shift from zero? How should I interpret my unstandardized B? Its value is negative. FYI. I'm testing to see if the increase in violence is influenced by a positive or negative change in vegetation cover (mean centered).


title: Alternative to Pearson's chi-square goodness of fit test, when expected counts < 5

text: In a part of my research, I am fitting probability distribution models for the count data(and binomial data) using Poisson, Binomial, Negative Binomial and Beta Binomial models. I have few data sets in which the observed frequency and hence the expected modeled frequencies are too small(less than 5).. I want to compare the goodness of fits using a proper method. I know, Pearson's Chi- Square test will not be applicable for my purposes as it has some limitations on expected counts. Here is a similar situation which I also deal with.. ![enter image description here](http://i.stack.imgur.com/zS6P9.png) I am looking for an alternative Goodness of fit test measure which can be applied to this kind of situations. Any suggestion is greatly appreciated.. Thank you.


title: R Error in manova(classification ~ ., data) : need multiple responses

text: I'm trying to get a MANOVA running to get Wilk's Lambda from the summary but I get a constant "need multiple responses" error. Looking at past answers I figured that one of the variables might be missing from the data frame but it doesn't seem like it.. here's head(data):               head(data)           ror penalty  last30 classification     1 4.97782 1.35074 4.83647           Good     2 4.65229 1.17034 4.99381           Good     3 4.91005 2.33616  4.8811           Good     4 4.85295 1.90597 4.82852           Good     5 4.79894 1.67042 4.81287           Good     6 4.50941 0.96402 4.88748           Good      head(classification)                head(classification)     [1] "Good" "Good" "Good" "Good" "Good" "Good"      code:               raw <- read.csv("~/Downloads/rp_r/draft_b.csv")     params <- cbind(raw$scaled, raw$penalty_scaled, raw$scale)         classification <- c(1:length(raw$final_rank))      cutoff = 20       for (i in 1:length(classification)) {         if (raw$final_rank[i] < cutoff) {           classification[i] <- "Bad"         } else {           classification[i] <- "Good"         }       }            data=data.frame(cbind(ror=params[,1], penalty=params[,2], last30=params[,3], classification))     m <- manova(classification ~., data)         summary(m,test="Wilks")      Thanks heaps!


title: What does it mean to have negative values for both lower and upper bounds on a 95% Confidence Interval on ANOVA?

text: My data shows 95% CI = [-1.07, -0.40], what does the negative sign indicate and how do I interpret the data? as I think it is significant for not capturing the value of 0.


title: Integrating an empirical CDF

text: I have an empirical distribution $G(x)$. I calculate it as follows                   x <- seq(0, 1000, 0.1)         g <- ecdf(var1)         G <- g(x)      I denote $h(x) = dG/dx$, i.e., $h$ is the pdf while $G$ is the cdf. I now want to solve an equation for the upper limit of integration (say, $a$), such that the expected value of $x$ is some $k$. That is, integrating from $0$ to $b$, I should have $\int xh(x)dx = k$. I want to solve for $b$. Integrating by parts, I can rewrite the equation as $bG(b) - \int_0^b G(x)dx = k$, where the integral is from $0$ to $b$ ------- (1) I think I can calculate the integral as follows                   intgrl <- function(b) {             z <- seq(0, b, 0.01)             G <- g(z)             return(mean(G))          }      But when I try to use this function with                   library(rootSolve)         root <- uniroot.All(fun, c(0, 1000))      where fun is eq (1), I get the following error                   Error in seq.default(0, b, by = 0.01) : 'to' must be of length 1        I think the issue is that my function `intgrl` is evaluated at a numeric value, while `uniroot.All` is passing the interval `c(0,1000)` How should I solve for $b$ in this situation in R?


title: How to predict given Poisson regression?

text: I run some Poisson regressions with the following results: (with number of associations an individual belongs to as the dependent variable)               Coefficients:                      Estimate Std. Error z value       (Intercept)       -0.92      0.11     -8.43        age(decades)       0.07      0.01      5.64             female            -0.10      0.04     -2.57           education(degree)  0.18      0.02     10.55         income(in 1000     0.09      0.01      8.38          east              -0.46      0.05     -9.27     TV                -0.08      0.02     -4.98      I'm trying to make a nice hook in my work however I have struggles to capture the following: I would like to say something about a 25 year west male who holds a particular degree(education=5) earns 1000 .- per month and watches 1 hour TV per day.


title: Linear model with categorical dependent variable

text: I need to identify relationships between multiple 'counting' independent variables (representing numbers of amino acids in a protein sequence) and one ordered categorical dependent variable (an experimental outcome) which assumes values between 1 and 7. The distribution of the dependent variable looks the following: ![Distribution of dependent variable.](http://i.stack.imgur.com/Zukwm.jpg) I know that for such cases linear regression models are not appropriate, but I'm under big time pressure and it would be a big plus if I can get something out of a linear model in a fast and dirty way. Hence my questions: **1) Do I have to expect biased coefficients even in the case I did not miss any relevant independent variable?** **2) Is the 'real' standard error larger, smaller or unrelated to the one given during the model estimation?** Below are the two diagnostic plots created in R. Any help would be greatly appreciated! ![QQ plot.](http://i.stack.imgur.com/cuiGc.jpg) ![Standardized residuals vs. fitted values.](http://i.stack.imgur.com/RblGQ.jpg)


title: Multilevel model where an interaction is a varying slope

text: I have a multilevel problem where I want to have a random intercept and a random slope. However the random slope is the interaction of two predictors. In this case, do I also have to allow random slopes for the individual predictors? For example, with lme4 is this OK? lmer(y ~ x + w + x:w + (1 + x:w|group), data=mydata) Or should i use lmer(y ~ x + w + x:w + (1 + x + w + x:w|group), data=mydata) even if I am not interested in either x or w to be random?


title: Benchmarking model in speech recognition with different language

text: My supervisor asked me to benchmark my method in classifying speech signal with other language. I am doing Malay language speech recognition. To benchmark my method/feature used, I need to test English speech. I am wondering, while doing testing (on the English data), should I use English data as training and testing using the same method/feature in Malay, or using Malay trained data to test the English data. I am really confused. Please help me.


title: p value for the test statistic in python using rtool

text: I have the following data,               These are the x and y values to find correlation      x=115778,171235,1     y_a=31920,49327,0     y_b=83858,121908,1     cor_a= 0.99947573036   // These are the correlation values using numpy (x,y_a)     cor_b= 0.999915675755   // (x,y_b)     test_statistic_a= 30.870014264  // test statistic value using formula     test_statistic_b= 76.9983294567     d_f=1  // degrees of freedom       How to find the p_value for the above data and test the hypothesis? using rpy2


title: Understanding ROC curve[edited]

text: I'm having trouble understanding ROC curve. Is there any advantage/improvement in area under the ROC curve if I build different models from each unique subset of train set and use it to produce probability? For example, if y has values of $\\{a, a, a, a, b, b, b, b\\}$, I build model $A$ by using $a$ from 1st/4th values of $y$ and 8th/9th values of $y$ and build model $B$ by using remained train data. Finally, generate probability. Any thoughts/comments will be much appreciated. Here is r code for better explanation for my question               Y= factor(0,0,0,0,1,1,1,1)     X=matirx(rnorm(16,8,2)     ind=c(1,4,8,9)     ind2=-ind          mod_A=rpart(Y[ind]~X[ind,])     mod_B=rpart(Y[-ind]~X[-ind,])     mod_full=rpart(Y~X)          pred=numeric(8)     pred_combine[ind]=predict(mod_A,type='prob')     pred_combine[-ind]=predict(mod_B,type='prob')     pred_full=predict(mod_full,type='prob')      so my question is, area under ROC curve of pred_combine VS pred_full


title: Analyzing neuralnet functions in R

text: The 'neuralnet' package in R allows us to use neural network algorithm with backpropagation. I want to use the function for prediction. I saw a tutorial on neuralnet in which predictions on the iris dataset were made. I am still confused as to how the parameters in neuralnet function affects the output, and how to make the prediction more accurate. I am using compute method for prediction, and I want to know what is the ratio of training and test data set; and is there any minimum of records that need to be added for good training of the neuralnet function?               neuralnet(formula, data, hidden = 1, threshold = 0.01,                       stepmax = 1e+05, rep = 1, startweights = NULL,                learningrate.limit = NULL,                learningrate.factor = list(minus = 0.5, plus = 1.2),                learningrate=NULL, lifesign = "none",                lifesign.step = 1000, algorithm = "rprop+",                err.fct = "sse", act.fct = "logistic",                linear.output = TRUE, exclude = NULL,                constant.weights = NULL, likelihood = FALSE)      How do I choose the value of hidden, threshold, and when to choose which activation function? I have read the paper on neuralnet and CRAN documentation too, but if you guys have some interesting read, please share, I am here to learn.


title: Does covariance stationarity lead to mean stationarity necessarily?

text: Traditionally a weak stationary process is also called covariance stationary, but those 3 properties are exposed: $$E[Xt] = μ , \forall t$$ $$var(Xt) = \sigma^2, \forall t$$ $$cov(Xt, Xt−j) = \gamma_j, \forall t$$ that are respectively mean, variance and covariance stationarity. We can merge 2 and 3 by extending to $j=0$ and $\gamma_0 = \sigma^2$, and let's call it covariance-only stationarity But seeing the name **covariance stationarity** and also trying to figure out examples where the variance is stationary and not the mean, or where covariances (j=0 too) are stationary but not the mean, I was wondering if maybe covariance-only stationarity could lead to mean stationarity.   1. Is there a proof for covariance-only stationarity → mean stationarity?   2. or are both covariance-only stationarity and mean stationarity required for having a covariance stationary process? note: if 2., in my opinion 'covariance stationary' is a very misleading term, weak stationary would be better to use


title: What is data augmented by the additive inverse?

text: I am reading Biclustering of expression data (Cheng and Church, 2000) The paper is about the Cheng and Church biclustering algorithm and its main metric, the mean squared residue (MSR). It is said that > A particular score that applies to expression data transformed by a > logarithm and augmented by the additive inverse is the mean squared residue. To me this means that MSR can identify scale biclusters when data is log transformed (because they become shift biclusters and have additive increments instead of scale biclusters that have multiplicative increments). But what does it mean that it applies to data augmented by the additive inverse? I was thinking that this score is minimized (good) also when elements of bicluster show a inverse match, like in the following image: ![Credits to www.kemaleren.com](http://i.stack.imgur.com/RTjvC.png) (This image is taken from the CC post on www.kemaleren.com) This can be interesting because, when applied to gene expression data, it clusters together genes that behave in the same (or in the opposite way) in some conditions.


title: Is the p value in a permutation test the same as the p value in the original test?

text: This is an extension from my last post. It seems to be too long to discuss there.   * When testing some null $H$ versus alternative $K$ by a test statistic $U(X)$, the p-value for $U$ on a sample $X$ can be defined as the infinimum of size $\alpha$ of the test rule, with varying critical value $c$ while still rejecting the null on $X$: $$p_U(X) := \sup_{F \in H} P(U(Y) \geq U(X) | Y \sim F), \quad i.e. \quad \alpha(U(X)),$$ where $\alpha(c)$ means the size of the test rule $(U, c)$.   * When the permutation test with the set $G$ of permutations can be applied to the above test (the null, $U$ and $X$), the p-value of the permutation test is said to be $$ p_P(X) := \frac{\\# \\{\pi \in G: U(\pi X) \geq U(X)\\}}{|G|} \\\ \stackrel{?}{=} P(U(Y) \geq U(X) | Y \sim F_0) $$ where $F_0 \in H$ is a common distribution of $X$ and $\pi X$ under null, $\pi X$ means applying the permutation $\pi$ to the sample $X = \\{x_1, ..., x_n\\}$, and the part after $\stackrel{?}{=} $ is what I thought of myself. Questions:   1. I was wondering if $p_U(X) = p_P(X)$? I.e., is the p-value of the permutation test equal to the p-value based on $U(X)$, for sample $X$?   2. If yes, how are they equal when the null is composite?      * $p_U(X)$ considers all the null distributions $F$'s one by one and then take $\sup_{F\in H}$,     * $p_P(X)$ seems to equal $P(U(Y) \geq U(X) | Y \sim F_0)$ for a single null distribution $F_0$ of $X$, and doesn't consider other null distributions of $X$. Thanks!


title: Progression between two exams under three conditions

text: I want to test if three types of learning mehods have different effects on the progression between two exams. The data was collected like this:   * A first exam was done with all participants. They had various level of knowledge of the subject, ranging from expert to novice.   * The participants were randomly separated into three groups which were given a lecture based on text, pictures or animations.   * They all passed a second exam (the same for everyone, but not the same as the first one) on the same subject. I would like to know which learning method was the best one in that case. I was thinking of calculating the "progression" like this: $$ \text{score on the second test} - \text{score on the first test} $$ And to run an ANOVA on this. However I realized that participants with very good previous knowledge had a low progression because they scored high on both tests while participant with very low previous knowledge had a much better progression even if the score at the second test was below average. I thought of removing the participants with high scores on the first test, but I am not sure it is a very good solution. I also thought of doing something with the ranking: average the ranking of the user under each condition before and after the lecture and see if one condition allowed a greater increase in average ranking. What do you think?


title: Gaussian process regression toy problem

text: I was trying to gain some intuition for Gaussian Process regression, so I made a simple 1D toy problem to try out. I took $x_i=\\{1,2,3\\}$ as the inputs, and $y_i=\\{1,4,9\\}$ as the responses. ('Inspired' from $y=x^2$) For the regression I used a standard squared exponential kernel function: $$k(x_p,x_q)=\sigma_f^2 \exp \left( - \frac{1}{2l^2} \left|x_p-x_q\right|^2 \right)$$ I assumed that there was noise with standard deviation $\sigma_n$, so that the covariance matrix became: $$K_{pq} = k(x_p,x_q) + \sigma_n^2 \delta_{pq}$$ The hyperparameters $(\sigma_n,l,\sigma_f)$ were estimated by maximizing the log likelihood of the data. To make a prediction at a point $x_\star$, I found the mean and variance respectively by the following $$\mu_{x_\star} = k_\star^T (\mathbf{K}+\sigma_n^2\mathbf{I})^{-1} y$$ $$\sigma_{x_\star}^2 = k(x_\star,x_\star)-k_\star^T(\mathbf{K}+\sigma_n^2\mathbf{I})^{-1} k_\star$$ where $k_\star$ is the vector of the covariance between $x_\star$ and inputs, and $y$ is a vector of the outputs. My results for $1<x<3$ are shown below. The blue line is the mean and red lines mark the standard deviation intervals. ![The results](http://i.stack.imgur.com/vNc1r.png) I'm not sure if this is right though; my inputs (marked by 'X's) do not lie on the blue line. Most examples I see have the mean intersecting the inputs. Is this a general feature to be expected?


title: Is R output reliable (specially IRT package ltm)

text: This may seem like a foolish question, but it is critical for me to have some expert advice on this. I am doing some research on students' learning levels in India. I am using R software (`ltm` package) to calculate student scores (ability estimates as per IRT model). Features of the data: No of students: 15,000 Test papers: 3 different forms No of questions in one form: 50 My concern is whether R outputs are reliable enough to publish the results in research community or the results from R software can be questioned, specifically with respect to `ltm` package. Also if R can handle such large data comfortably without compromising the quality. I would be very grateful if you could share your inputs and any documentation/research paper/reference you may have.


title: Approximating a Poisson distribution using a partially observed Gaussian

text: _Note that the problem specification has changed since the original posting. Thanks to whuber for helping me better specify the question: in an attempt to make the question general, I had left out several important constraints._ **Problem** : Let $y$ be distributed Poisson with rate parameter $\lambda$: $p(y | \lambda) = \frac{\lambda^y}{y!}e^{-\lambda}$ Let $z^*$ be distributed Gaussian with mean $\mu$ and variance $1$: $p(z^* | \mu) = \frac{1}{\sqrt{2 \pi}}e^{-\frac{(z^* - \mu)^2}{2}}$ Define $z$ as follows: $z = k$ if $\gamma_{k} \le z^* < \gamma_{k+1} $, where $\gamma_0 = -\infty$. Let $\mu = \log\lambda$. How do we choose a single sequence of $\gamma_k$ for $k = 1, 2, \ldots, \infty$, that does not depend on $\lambda$, in order to best approximate the distribution of $y$ with the distribution of $z$, across the full range of possible values of $\lambda \in (0,\infty)$. It is not clear to me exactly what the right criterion is here for operationalizing "best approximating" across the full range of $\lambda$. **Background** : The idea here is to do something similar to the Albert and Chib (1993) approach to Gibbs sampling binary and categorical models. Instead of having a single or small number of cutpoints, here we would have an infinite series of cutpoints. The question is how to choose the cutpoints to best approximate a Poisson distribution. As a practical matter, it is pretty easy to draw doubly-truncated normal random variates (Robert 2009). There is an existing approach to using data augmentation to fit a Poisson regression (Frühwirth-Schnatter and Wagner 2005), but it seems to involve rapidly increasing computation with increasing values of the count variable, which is a very unattractive feature for my application. I am working on a Gibbs sampler for a hierarchical model to be applied to a very large dataset of count data. I have no particular reason to believe that the distribution of counts are particularly close to Poisson, I am mostly interested in computational efficiency and not having a probability distribution for the counts that is very far from a standard distribution. **References:** James H. Albert and Siddhartha Chib, "Bayesian Analysis of Binary and Polychotomous Response Data" Sylvia Frühwirth-Schnatter and Helga Wagner, "Data Augmentation and Gibbs Sampling for Regression Models of Small Counts" Christian Robert, "Simulation of truncated normal variables"


title: Identifying parameters in BUGS linear regression

text: With a linear regression defined in BUGS, how should one implement model identification constraints such as having the mean of a group of parameters be zero, or having the group of parameters sum to zero? EDIT: To clarify, I'm thinking of a model with, say, varying slopes, e.g.,               for (i in 1:N) {       y[i] ~ dnorm(mu[i], tau)       mu[i] <- beta0 + beta1[group[i]]*x1[i]     }      What's the right way to specify that the mean of all the beta1[] is zero? I _don't_ want to say that each beta1[group[i]] is drawn from a distribution with zero mean; I do want to say that all the beta1[] are drawn from a distribution the mean of which is fixed at zero. In pseudocode that is probably not legal BUGS or even JAGS, something like               for (j in 1:nGroups) {       beta1[j] ~ dnorm(0, beta1.tau)     }     mean(beta1.mu) <- 0


title: Gaussian Process covariance matrix gets zero determinant

text: I have a Gaussian process regression implementation and developed some example data to test the capabilities of those methods. In the posterior calculation one gets the covariance matrix $K$. For some sample data this matrix has a 0 determinant and thus it is not invertible. Can someone see a problem in the covariance matrix composition that leads to such behaviour? My Covariance matrix looks like this: $$ \begin{pmatrix} K(X,X) & K(X_*,X) \\\ K(X, X_*) & K(X_*,X_*) \end{pmatrix} $$


title: Sampling from an unknown distribution

text: I am using the principle of maximum entropy to fit a density to a given set of samples. I want to generate new set of samples from the approximated density. Is there any way to draw samples from a density obtained by maximum entropy approach?


title: Predicting dependent variables where independent variables are samples with differing sample size

text: Is there a standard method of dealing with independent variables, where they are samples from a known distribution, and the sample sizes differ from sample to sample? I'll give an example of what I mean. Suppose your sample consists in $n$ individuals, who have been given 3 different types of test. Let $Z_{i,j}$ be number of successes for person $i$, on test $j$; $n_{i,j}$ is the number of times person $i$ takes test $j$ (the sample size). We further assume that $Z_{i,j} \sim B(n_{i,j}, p_{i_j})$, where $p_{i,j}$ represents the ability of person $i$, at test $j$. The different individuals have different abilities at the test. We want to use the sample to construct a model that let us use the results of test 1 and 2 for an individual to predict their result on test 3. We could fit a generalized linear model, for example, where $Z_{i,3}$ is the independent variable, and $\frac{Z_{i,1}}{n_{i,1}}$ and $\frac{Z_{i,2}}{n_{i,2}}$ are dependent variables, but we wouldn't take into account the differing sample sizes. Ideally, we want to construct a model that accounts for the relationship between an individual's ability on the different tests ($p_{i,1}, p_{i,2},p_{i,3}$), and use that to predict his ability on test 3, given his ability on tests 1 and 2, but we only know the results of the tests, the $Z_{i,j}$'s. ADDITION: Suppose we have a model of the relationship between each of the $p_{i,\cdot}$'s, for example $p_{i,3}=f(\beta_0+\beta_1g(p_{i,1})+\beta_2g(p_{i,2}))$, where $f$ and $g$ are, say, logit functions. How would I obtain those parameters? Is there a standard algorithm for this, or a function in R? Is maximum likelihood an appropriate method of fitting? Also, when fitting the model, what pitfalls are there? Are there metrics or tests for comparing different models, or determining goodness-of-fit? Is there a term for this form of modelling?


title: Data dependency detection software

text: I need to have **software for detecting dependencies** among data attributes. The simplest example: we have a database table representing results of some social quiz. Columns represent questions; rows, people who passed the quiz. We are searching for different interesting dependencies in the data. _Examples:_   * 60% of men are introverts;    * 90% of people in USA aged 7-15 like milk chocolate;    * passport number identifies people aged 16+;    * children document number identifies people aged <16; etc, etc, etc. There are _many types of possible dependencies_ and a lot of methods for finding them (actually I am interested in methods too, so good links and books are appreciated). Valuable information for me would also include statistics of the importance of different types of dependencies: how to determine which data dependency method is better if we have sets of dependencies found by those methods? We need to somehow mark every dependency (surely marks will be different for different subjects). Does anybody know of software (ideally in .NET and free) that can analyze data that way? I need real products for data dependency detection, not software for creating such products (so please, do not offer Matlab, C#, C++ etc. as answers).


title: Simple, robust and fast algorithm for stochastic gradient descent

text: I'm trying to solve an optimisation problem with stochastic gradient descent with the following properties:   * It has a very large (1,000,000+ element) parameter vector.    * Empirically, there seems to be a single maximum (though I can't prove this) so hill climbing is fine, however the problem is definitely not convex   * I can get gradient samples at any point, but the samples have quite a bit of noise and getting a large number of samples is expensive.   * It needs to be an online algorithm Currently, I'm using a simple online gradient descent with momentum. It works, but has two problems:   * It seems quite slow to converge   * It requires quite careful hand-tuning of the learning rate Is there a better algorithm that I could use for this situation?


title: Probabilistic graphical models textbook

text: Is Koller's "Probabilistic Graphical Models" suitable as a textbook? Or is there another book which is more recommendable as textbook for a master- course? Disclaimer: cross-posted from quora.com, where I got no answer.


title: Formulating a hierarchical model in educational measurement

text: I'm having some trouble positing a model in levels for the question below. **Question** : How does the percentage of students born in the United States at the school level moderate the effect of the use of applications or models while learning science (denoted as SCAPPLY) on PISA science scores (PISA - outcome variable)? **Level-1 Predictors:** ESCS = Index of economic, social and cultural status (or socioeconomic status (SES)); SCAPPLY; SCINVEST = An index measuring the degree to which students report that they conduct investigations when learning science; SCINTACT = Student Interaction. **Level-2 Predictor** : NATIVE = the percentage of students born in the US for each school in the sample In the model, we allow the effect of ESCS to vary at level-2; all other slopes at level-1 remained fixed at level-2. Also, we allow the school-level intercepts to vary and the variability in the ESCS slopes to vary with the school-level intercepts. So far I have the following. \----------------------------------------------------EDIT------------------------------------------------------- Level 1: $$PISA_{ij} = \beta_{0j} + \beta_{1j}(ESCS_{ij}) + \beta_{2j}(SCAPPLY_{ij}) + \beta_{3j}(SCINTACT_{ij}) + \beta_{4j}(SCINVEST_{ij}) + e_{ij}$$ Level 2: \begin{align} \beta_{0j} &= \gamma_{00} + \mu_{0j}\\\ \beta_{1j} &= \gamma_{10} + \mu_{1j}\\\ \beta_{2j} &= \gamma_{20} + \gamma_{21}(NATIVE_j)\\\ \beta_{3j} &= \gamma_{30}\\\ \beta_{4j} &= \gamma_{40} \end{align} Did I posit my models in levels correctly?


title: Difference-in-difference-in-differences estimator

text: Two towns, $X$ and $Y$. In each town:   * Pool cross-sections of male and female hourly wages, one from the year before a wage-discrimination policy took effect and one from the year after. Consider the following model:   * $wage_i=\beta_0+\beta_1after_i+\beta_2female_i+\beta_3X_i +\beta_4after_iX_i +\beta_5after_ifemale_i +\beta_6X_ifemale_i +\beta_7after_ifemale_iX_i+u_i$ where   * $after_i=1$ if date after gender-wage discrimination policy; $0$ if date before gender-wage discrimination policy. **The gender wage discrimination policy only applies to females in town $X$.**   * $female_i=1$ if female; $0$ if male.   * $X_i=1$ if town $X$; $0$ if town $Y$. The treatment group is females from town $X$. **Does this model measure the effect of a gender-wage discrimination policy on the average wage of women relative to men in town $X$ compared to the average wage of women relative to men in town $Y$?** In this model, the affect of the policy is captured by $\beta_7$:the 'difference-in-difference-in-differences' estimator. For $\hat{\beta_7}$ to be interpreted as the causal effect of the policy on wages in town $X$, $E(u_i|female_i, after_i, X_i)=0$ has to be assumed. **Does this mean that once gender, date of policy, and town have been controlled for, there exist no unobserved factors which _change over both town and time_ that affect wages?**


title: How do I test statistical significance between normally distributed data and non normally distributed data?

text: I have observed values (80 in total) in different measured circumstances (80 in total). I want to test whether my observed values are the result of chance or do measured values affect my observed values. My problem is that measured values are non-normally distributed and observed values are normally distributed according to Shapiro-Wilk test which I run on SPSS. Measured values' p-value was 0,001 in that test and observed values' p-value was 0,102 in that test. What test should I use that I would know that are my observed values statistically significant? All I need to test is that are my observed values result of chance or did measured values affect them. If I use t-test, are both measured and observed values required to be normally distributed? Thank you for reading, any help is appreciated


title: Where did the frequentist-Bayesian debate go?

text: The world of statistics was divided between frequentists and Bayesians. These days it seems everyone does a bit of both. How can this be? If the different approaches are suitable for different problems, why did the founding fathers of statistics did not see this? Alternatively, has the debate been won by Frequentists and the true subjective Bayesians moved to decision theory?


title: Can I use a confidence interval of poisson mean for its variance

text: In a Poisson distribution the mean equals the variance. I would like to find a confidence interval of the variance. Is my reasoning below correct?   Using the central limit theorem I construct a 95% confidence interval for the mean $\mu$   $L \leq \mu \leq U$   $\mu=\sigma^2$   Therefore   $L \leq \sigma^2 \leq U$   It would seem to me that the inequality should work like any other inequality in mathematics but statistics can throw a curve ball sometimes so I'm not certain about it. I can't find any papers discussing if this approach is valid.   Another good example of this is a confidence interval for the mean and median of a normal distribution. The mean confidence interval is smaller but the median confidence interval is more robust so either might be preferred as an estimate of the other.


title: How can I compare two classifications?

text: My job is to develop a classification that groups together diagnoses of similar cost and procedures of similar cost, by similar anatomy. For instance a hip replacement would not be in the same group as a heart transplant as they cost different amounts and are not clinically similar. We revise this huge classification yearly and I've been tasked to prove that changes we make, make it 'better'. I know :-( Numerically we used to look at reduction in variance, does this way of grouping explain more variance than this way of grouping? However, a lot of the clustering decisions are based on clinical reasoning, not numerical. Two procedures may cost the same but you may not put them in the same group as one requires anaesthetic and one doesn't! This obviously cannot be measured in the same way. I'm not wanting to assess the accuracy of the algorithm (ie has the patient dropped into the right group for what is wrong with them) but a method of comparing the way groups are defined. Thanks!! I couldn't add a comment to your question so I've added here: The classification is used to collect costs from hospitals. Each hospital calculates the average cost to them to deliver the service described by each group within the classification (there are around 1500 groups within the classification such as Heart transplants, disorders of spleen, procedures on lymphatic system etc..). These hospital level costs are collated to provide a national average for each group which in turn is used to calculate a tariff for the government to pay those hospitals for that service. The variance I described is the variance in costs reported - one hospital may say it costs them £2000 to do a hip replacment, another may say £8000. If we say there are 20 different ways to define a hip replacement (different kinds of surgery), 19 of the 20 may cost around 2K while one may cost £20k. It would be appropriate to move that 20K hip replacement into a different group with other 20K procedures. How could we show statistically that making this change has improved the quality of the costs within each group by reducing the cost variance in each group? Hope that makes sense!


title: "On/off problem" resources

text: I need to study in depth the on/off problem with different approaches: frequentist, bayesian - frequentist hybrid, profile likelihood. The on/off problem is a counting experiment where you observe $n_{on}$ events when the source is "on" in the time $T_{on}$ and $n_{off}$ events when the source is "off" in time $T_{off}$ and you want to make inferences about the source power (the simpler case power > 0 ?) Do you know some good books? Some references on the problem:   * Tucker 2008   * Cousins Linemann & Tucker 2008   * Cranmer 2011 (page 36)


title: How can the F distribution be used, other than for hypothesis testing and confidence interval estimation?

text: I am trying to fit informed prior distributions to data using MLE, and F occasionally provides a best fit (lowest AIC value). I am starting with only very basic knowledge of probability theory, so I am not even confident that it is appropriate to fit an F in these cases. One of the interesting aspects of probability theory is that distributions often represent specific data generating processes - and are therefore useful for specific data types. However, I can not find any such justification for the F distribution. For example, the NIST handbook states: > Since the F distribution is typically used to develop hypothesis tests and > confidence intervals and rarely for modeling applications, we omit any > discussion of parameter estimation. However, relationships between the F and Beta distributions suggest that F might be appropriate in general for ratios. Are there theoretical (not computational) reasons to use or not use the F distribution as an informed prior distribution, or is having a continuous parameter with values > 0 sufficient justification?


title: How to determine whether a classifier like adaboost is weak?

text: I run the cross-validation experiment for a given data set, and tried two different approaches: one is based on SVM, another is based on SVM plus Adaboost. But the confusion matrix for two experiments are exactly the same. I am confused on how to explain this kind of result. Adaboost is supposed to start with a weak classifier, but how to determine whether a classifier is weak?


title: What type of data analysis (ANOVA) should I run on this experiment?

text: I am not sure how I would analyse this data: The experiment is looking at the accuracy of emotion recognition in children with ADHD. The children will be asked to identify what emotion a face is portraying (out of 6 expressions; happiness, sadness, fear, anger, disgust and neutral). They will firstly be shown the photographs and asked to identify the emotion, then they will be instructed to focus on the eyes and asked to identify the emotions and then asked to focus on the mouth. (three conditions- free gaze, eye gaze and mouth gaze). The experiment involves 2 groups of participants, ADHD children and matched controls. The DV is accuracy of emotion recognition. How would I analyse this, would this be a repeated measures ANOVA to compare accuracy of emotion recognition by participant group (ADHD or controls) gaze condition (free, eye, mouth) and emotion type? I think this will determine whether there was a significant difference in the accuracy of identifying emotions in the difference gaze conditions and types of emotion portrayed. And if significant, follow up Bonferroni tests will be conducted to show which emotions were recognised the least accurately? t tests? planned comparisons? Im not too sure what design this is or what type of ANOVA should be conducted. Would I do seperate ANOVAs for each emotion? I am a bit confused!


title: ANOVA for binomially-distributed data

text: Is ANOVA valid for binomial data? I have 7 binomial groups and I want to test whether the distributions of all of the groups are equal or whether the distribution of the first group is significantly different than the distribution of groups 2-7. What tests would accomplish this? I was thinking of just combining groups 2-7 and then performing a chi-square test on group 1 vs groups 2-7. Is this valid? Thank you!


title: Significance of varibles after stepwise regression

text: I did stepwise regression with my multiple regression model and using AIC as a measure of fit with the `step` function in R. Afterwards some variables that the stepwise regression did not eliminate was not significant (> 0.05 p-value). Does this mean i have to take out those variables with large p-values or what is a normal procedure?


title: Formula symbols for mixed model using lme4

text: **Edit Note:** Since I posted this question, it has been suggested that I read some documents and I am still thinking on the subject. I have added some new understanding as marked between _*_ * asterisks: please do correct them if I am wrong and add remaining questions.* I often confuse four symbols ":", "|", "/" and "*" in the formula, particularly while doing mixed model. Can somebody explain clearly the differences between them? I have the following working example:               require(lme4)     mydata <- expand.grid(xvar1 =factor(1:10), xvar2 =factor(1:3),replication = factor(1:3))     mydata$yvar <- rnorm(nrow(mydata), 10, 5)          fm1 <- lmer(yvar ~ 1 + (1|xvar1) + (1|xvar2) + (1|xvar1:xvar2), mydata)     ***My understanding: Cross classification between xvar1, xvar2***            fm2 <- lmer(yvar ~ 1 + (1|xvar1) + (1|xvar2) + (1|xvar1/xvar2), mydata)     ***is not correct model, should remove 1|xvar1, 1|xvar2***      ***** has fixed intercept 1 , I do not know if that is correct technically**      **fm2 <- lmer(yvar ~ 1 + (1|xvar1/xvar2), mydata)*****       *is xvar1 is nested in xvar2 and is essentially same as fm1*          fm3 <- lmer(yvar ~ 1 + (1|xvar1) + (1|xvar2) + (xvar1|xvar2), mydata)     Warning message:     In mer_finalize(ans) : singular convergence (7)     ***I do not have idea on "|",***           fm4 <- lmer(yvar ~ 1 + (1|xvar1) + (1|xvar2) + xvar1:xvar2, mydata)     Error in mer_finalize(ans) : Downdated X'X is not positive definite, 31 ;     ***xvar1:xvar2 interaction term, technically correct !***          fm5 <- lmer(yvar ~ 1 + (1|xvar1) + (1|xvar2) + (1|xvar1*xvar2), mydata)     Error in validObject(.Object) :        invalid class "ngTMatrix" object: all row indices must be between 0 and nrow-1     In addition: Warning message:     In Ops.factor(xvar1, xvar2) : * not meaningful for factors         **again (1|xvar1) and (1|xvar2) are not necessary**     ***fm5 <- lmer(yvar ~  1 + (1|xvar1*xvar2), mydata) is equal to***      ***fm5 <- lmer(yvar ~ 1 + (1|xvar1) + (1|xvar2)+ (1|xvar1*xvar2)***               fm6 <- lmer(yvar ~ 1 + (1|xvar1) + (1|xvar2) + (xvar1/xvar2), mydata)      ***same as fm2***          fm7 <- lmer(yvar ~ 1 + (1|xvar1) + (1|xvar2) + (0 + xvar1|xvar2), mydata)        Warning message:        In mer_finalize(ans) : singular convergence (7)        ***I have no idea on xvar1|xvar2, 0 means no intercept I believe***       **Edit: Here is something I learned from R documentation on linear model formula**               (1) yvar ~ xvar1 + xvar2 + xvar1:xvar2  - cross over classification          is same as yvar ~ xvar1 * xvar2           (2) yvar ~ xvar1 / xvar2  - nested classification (         however means xvar1 + xvar2 + xvar1:xvar2)          is same as yvar ~ xvar1 %in% xvar2      I am not sure those applicable to lm model hold true for mixed model. I am still not on track on use of "|" which is I believe unique in mixed models.


title: Detecting outliers using 95% PI around a natural spline fit

text: Please see the picture below: ![enter image description here](http://i.stack.imgur.com/pjfW9.png) I wanted to mark the points that are not consistent with their adjacent points as outlier. What I did was to fit a natural spline fit to 1000 observations (the purple line in the picture is the fitted line). I then consider 95% prediction interval and 99% prediction interval around the fit. I consider any point outside of the 99% PI as outlier. If you look at the picture above, you would see that the red dots are the points that were selected as outliers. In particular look at the outliers after x = 150 (x axis). These are very consistent and should have not been selected as outlier. As you noticed the PI interval is very narrow and I think this is because of two things:   1. n = 1000 is large   2. data is highly correlated Question is how can I take this correlation into account and get a realistic 95% PI? Do I need a time series approach? What is your suggestion? Your help is greatly appreciated


title: distance measure of angles between two vectors, taking magnitude into account

text: Suppose I have two vectors, `v1` and `v2`, from which I can calculate the angle between these two vectors as a measure of their "distance", using the arccos function, say. For example:               v1 = c(100,200, 500,600)         v2 = c( 50, 30,  10,5)           v3 = c( 10,  7,  30,40)           # pairwise angles     acos( as.numeric((v1 %*% v2) / (norm_vec(v1) * norm_vec(v2))) ) * 180 / pi  # 66.8017     acos( as.numeric((v2 %*% v3) / (norm_vec(v2) * norm_vec(v3))) ) * 180 / pi  # 66.67337     acos( as.numeric((v1 %*% v3) / (norm_vec(v1) * norm_vec(v3))) ) * 180 / pi  # 8.061138      This kind of measure will give similar distances (angles) **regardless of the magnitude of the elements in the vectors**. For instance, the distance between `v1` and `v2` is 66.80, and the distance between `v2` and `v3` is similarly 66.67, but clearly the **magnitude** of `v2` and `v3` are more similar than `v1`, so I am thinking of a measure that will also take the "magnitude" into account when calculating the dissimilarity. In other words, dist(v1, v2) should be greater than dist(v2,v3), but this result is obtained by **still using the pairwise angle idea**. Thanks! **UPDATE** Thank you all for your replies! For the equivalence of Euclidean and angle distances, I use the same vectors as above to calculate the Euclidean distances:               # Euclidean distance     ed <- function(x1,x2)       sqrt(sum((x1 - x2) ^ 2))     ed(v1,v2)       # 790.9014     ed(v2,v3)       # 61.26989  --> imply v2, v3 are "closer"     ed(v1,v3)       # 761.4782  --> imply v1, v3 are "far apart"      As you can see, `v1, v2` and `v1, v3` have similar magnitude of distances in Euclidean case, but for arccosine, they are different (66.8 vs. 8.06). Is there anything particular that is revealed by the angle distance but not the Euclidean distance? I think the **orientation** information is emphasized in the angle distance.


title: What is the nonparametric variant of multivariate ANOVA with more than one one factor (and with interactions between factors)?

text: I analyze real data where assuptions of classical multiway (more than one factor assumed) MANOVA (more than one response variable) are not met (normality, equal variances). Is there any nonparametric variant? Note that Friedman test in only alternative to repeated measures ANOVA, not to more general multiway MANOVA with interactions.


title: Simple technique to identify number of clusters in dataset

text: I have a survey app (programmed using Ruby On Rails), and I am required to cluster the responses. I am using a Ruby library called AI4R and my code (in the event it is useful...) looks like the below (example code filched from AI4R)               # 5 Questions on a post training survey     questions = [   "The material covered was appropriate for someone with my level of knowledge of the subject.",                      "The material was presented in a clear and logical fashion",                      "There was sufficient time in the session to cover the material that was presented",                      "The instructor was respectful of students",                      "The instructor provided good examples"]          # Answers to each question go from 1 (bad) to 5 (excellent)     # The answers array has an element per survey complemented.      # Each survey completed is in turn an array with the answer of each question.     answers = [                  [ 1, 2, 3, 2, 2],   # Answers of person 1                 [ 5, 5, 3, 2, 2],   # Answers of person 2               ]          data_set = DataSet.new(:data_items => answers, :data_labels => questions)          # Let's group answers in 4 groups     clusterer = Diana.new.build(data_set, 4)      This in turn lets me create graphs like this (the survey has questions which are linked to themes/axes). ![enter image description here](http://i.stack.imgur.com/Vsus8.png) The problem is that right now you have to pick (read guess) the number of clusters to pass into AI4R. I saw on Wikipedia that there is a technique called the Elbow Method (illustrative picture from Wikipedia), ![enter image description here](http://i.stack.imgur.com/OfGnA.jpg) which compares the number of clusters with the variance that they explain. This technique would be perfect for my needs, but I don't know how to implement it in Ruby (or with pen and paper). What stats technique can I use to calculate the number of clusters vs the percentage variance that they explain?


title: Specifying mixed-linear model in lme() for pre-/post-treatment design

text: Dear crossvalidated community! I am analyzing a data of a study in which every participant's performance in a memory task was measured twice (once under placebo, once under drug treatment in random order). Moreover every subject was tested in a blood test measuring the level of stable blood marker. I would like to analyze the data with respect to drug-effects and blood marker effects using lme() from the nlme() package. 'treatment' codes placebo/drug 'session' codes the first or second testing in order to account for learning/re-test effects 'blood' codes for the blood marker level (continous variable), which should be entered as a covariate in the analysis 'id' codes for the subjects Would the following equation correct for the described data:               # some data     data <- structure(list(id = structure(c(1L, 12L, 14L, 15L, 16L, 17L,      18L, 19L, 20L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 13L,      1L, 12L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 2L, 3L, 4L, 5L, 6L,      7L, 8L, 9L, 10L, 11L, 13L), .Label = c("1", "10", "11", "12",      "13", "14", "15", "16", "17", "18", "19", "2", "20", "3", "4",      "5", "6", "7", "8", "9"), class = "factor"), sess = structure(c(1L,      2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L,      2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L,      2L, 2L, 1L, 1L, 2L, 2L, 1L), .Label = c("a", "b"), class = "factor"),          treat = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,          2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L         ), .Label = c("drug", "placebo"), class = "factor"), memory = structure(c(2L,          12L, 17L, 18L, 1L, 20L, 19L, 15L, 14L, 16L, 8L, 7L, 5L, 9L,          4L, 6L, 10L, 3L, 11L, 13L, 2L, 12L, 17L, 18L, 1L, 20L, 19L,          15L, 14L, 16L, 8L, 7L, 5L, 9L, 4L, 6L, 10L, 3L, 11L, 13L), .Label = c("1",          "10", "11", "12", "13", "14", "15", "16", "17", "18", "19",          "2", "20", "3", "4", "5", "6", "7", "8", "9"), class = "factor"),          blood = structure(c(3L, 11L, 17L, 18L, 20L, 19L, 9L, 10L,          6L, 16L, 5L, 2L, 12L, 15L, 4L, 13L, 14L, 8L, 7L, 1L, 3L,          11L, 17L, 18L, 20L, 19L, 9L, 10L, 6L, 16L, 5L, 2L, 12L, 15L,          4L, 13L, 14L, 8L, 7L, 1L), .Label = c("1", "100", "17", "25",          "38", "39", "40", "44", "55", "6", "62", "66", "70", "77",          "81", "83", "85", "86", "89", "94"), class = "factor")), .Names = c("id",      "sess", "treat", "memory", "blood"), row.names = c(NA, -40L), class = "data.frame")          # the model     m1 <- lme(memory~sess*treat*blood, random=~1|id/treat, data=data)     summary(m1)


title: Which high performance methods exist to determine the statistical significance of web analytics metrics?

text: In the field of web analytics, how can one calculate the statistical significance of a metric (for instance Repeat Visitor's). Imagine that we have millions of visitors daily and we are registering each visit, through cookies maybe. We store data per visit/visitor in a relational database and want to show a number for interesting metrics on demand in a report. Performance is an issue here, so we need fast queries and computations. Also we don't want to make any assumptions, about the data distribution (definitely not binomial for all metrics). What are the best relevant techniques to solve this problem?


title: Get frequency of a noisy sample

text: I'm writing you because i mostly need a general algorithm, or even a hint, which i will afterwards port on my system (an embedded Arduino microcontroller); I have a set of around 650 samples covering 5 seconds representing a periodic signal but with quite a lot of noise; I need the frequency of the main/dominant signal. Samples look similar to this https://www.dropbox.com/s/fw196r6yf1awhrh/untitled2.bmp Here you have a dump with around 5k samples https://www.dropbox.com/s/efwvyn5oec7ixgg/samples.txt In glad to get any hints on how can i achieve this. Thank you, Alex


title: Does "x per y" always refer to the arithmetic mean?

text: I recognize that there are many types of averages. That said, does "x per y" (e.g. visits per day) always refer specifically to the **arithmetic mean**? For some context, I saw this recently where someone claimed the cited figure was actually in reference to a median. While I think it might have been OK to use the median when declaring something as an average, it seems to me that it is clearly wrong to denote anything besides the arithmetic mean with "x per y." This is highlighted by the common usage of "x/y" to indicate "x per y." Am I correct to always assume usage of the arithmetic mean here? **EDIT** : Apparently my original post was not entirely clear. I am specifically referring to situations where y is unit of time, so clearly "x per y" is a frequency or a rate. I would always interpret it as the (arithmetic) mean frequency or mean rate. However, I saw somebody else interpret it as the median frequency (the median of the y discrete time intervals).


title: Using MAD as a way of defining a threshold for significance testing

text: If I have a set of terms each term having a particular frequency associated with it (the number of the times the term has appeared in fixed corpus of papers), then is the following method of significance testing valid?   1. calculate the median absolute deviation (MAD) of the GO term frequencies in the given corpus, for sample $S$ : ${\rm MAD}(S) = 1.4826 \times {\rm median}(|x_{i} - {\rm median}(S) |)$   2. get ${\rm thresh} = 2.7\times MAD(S) + {\rm median}(S)$   3. use ${\rm thresh}$ as a threshold above which the GO terms are deemed significantly associated with the given corpus and below which the GO terms are deemed non-siginificant.


title: Generating volcano scatterplot

text: This might look a very silly question to many of you but please answer. I am interested in generating a volcano plot for my dataset, which has four columns and all values are in log2. column 1 has name while column 2, 3, 4 has values in treated, untreated states and fold change.               DDR1    7.8519007358    8.0207450402    0.1688443044     RFC2    8.1756591822    7.7784602732    -0.397198909     HSPA6       7.8186806878    7.274639204 -0.5440414838     PAX8    8.0207450402    3.8519007358    -4.1688443044     CALR    7.7784602732    1.8519007358    -5.9265595374     MAPK1   4.8519007358    7.7784602732    2.9265595374     MAPK1   4.8819007358    7.2784602732    2.3965595374     KRAS    3.8519007358    7.8519007358    4      I am intrested in generating scatter plot from these values with marking outliers at 2 fold. I am bit confused and puzzled with the command in R as it is asking more than one parameter like p-value and things which I do not have.                 drawVolcanoPlot(M,p,m=1,p.cut=0.05,p.transform=log10,ylab=NULL,colramp=NULL,na.rm=TRUE,...)      Kindly help me in generating this scatterplot with my dataset. Solutions in python are also welcome. Thank you for your time and consideration.