title: How do I determine which functional form is correct?

text: I have an assignment that asks me to explore different series of data with respect to the effects that an incinerator built in 1981 in Massachusetts, has on the price of houses in that area. I have around 30 different sets of data with regards to the houses ranging from 'number of bathrooms', 'distance from the incinerator', 'square footage of house' to things like 'price before incinerator was built' and the 'age of the house'. I've started off with which variables I believe explain the changes (or no changes) in the house prices in the area and have around 6 explanatory variables with 'price' as the dependent. I know the Ramsey RESET test can be used to confirm functional form is correct; however, I'm not sure on when I should log a variable, square a variable or just keep it a linear functional form? Should I adopt the attitude of trying out different functional forms on different variables until I come to a point where my regression passes the RESET test, White test for heteroscedasticity and the Normality tests, or should I try and think about what the data looks like before I start taking the log of variables randomly. Thanks for any responses in advance, sorry if this question seems strikingly simple, I'm new to all this.


title: The operation of chance in a deterministic world

text: In Steven Pinker's book Better Angels of Our Nature, he notes that > Probability is a matter of perspective. Viewed at sufficiently close range, > individual events have determinate causes. Even a coin flip can be predicted > from the starting conditions and the laws of physics, and a skilled magician > can exploit those laws to throw heads every time. Yet when we zoom out to > take a wide-angle view of a large number of these events, we are seeing the > sum of a vast number of causes that sometimes cancel each other out and > sometimes align in the same direction. **The physicist and philosopher Henri > Poincare explained that we see the operation of chance in a deterministic > world either when a large number of puny causes add up to a formidable > effect, or when a small cause that escapes our notice determines a large > effect that we cannot miss.** In the case of organized violence, someone may > want to start a war; he waits for the opportune moment, which may or may not > come; his enemy decides to engage or retreat; bullets fly; bombs burst; > people die. Every event may be determined by the laws of neuroscience and > physics and physiology. But in the aggregate, the many causes that go into > this matrix can sometimes be shuffled into extreme combinations. (p. 209) I am particularly interested in the bolded sentence, but I give the rest for context. My question: are there statistical ways of describing the two processes that Poincare described? Here are my guesses: 1) "A large number of puny causes add up to a formidable effect." The "large number of causes" and "add up" sound to me like the central limit theorem. But in (the classical definition of) the CLT, the causes need to be random variables, not deterministic effects. Is the standard method here to approximate these deterministic effects as some sort of random variable? 2) "A small cause that escapes our notice determines a large effect that we cannot miss." It seems to me like you could think of this as some sort of hidden Markov model. But the (unobservable) state transition probabilities in an HMM are just that, probabilities, which is by definition once again not deterministic.


title: Item analysis for a likert-type questionnaire - item discrimination, point-biserial, factor analysis, cronbach's alpha, other diagnostics?

text: I'm working on a questionnaire that is intended to measure two constructs. There are presently 26 items in the questionnaire -- eleven for construct 1 and 15 for construct 2 -- and each item has five available responses: "strongly disagree" "disagree" "neutral" "agree" "strongly agree". In our first iteration of the questionnaire, we used ten items for each construct, and had about 60 respondents. After the results came in, we made some judgments based on item discrimination (dichotomized), point-biserial correlation, cronbach's alpha after item deletion, and factor analysis, and dropped several of these questions and replaced them with others. We now have 68 responses to the present questionnaire (the one with 26 items). I'd like to make sure I cover my bases as we do the reliability and validity analysis, and am specifically wondering if we should carry out additional diagnoses or approach the ones we've already calculated in different ways. Our approach to identifying potentially problematic items is to look for:   1. item discriminations below .3 on the dichotomized responses - where agree/strongly is coded as 1 and the others are coded as zero, and we take the difference between the top 3'rd of respondents and the bottom 3'rd of respondents.   2. low point-biserial (spearman's) correlations between items and the total score for the associated construct (total score = mean response for that construct). We don't have a threshold defining "low" ... and is spearman's or pearson's the appropriate measure?   3. items whose inclusion causes cronbach's alpha to decrease -- identified by running tests with "cronbach's alpha after item is deleted." Would G6 or Omega be more appropriate here?   4. Factor analysis -- we've been looking for items with low loadings in either dimension. Here, we're wondering about two thngs: a) whether to "swap" factors from their intended construct to the one which they appear to load with more strongly. Is that a defensible approach, even if the measure doesn't have face validity with respect to its intended construct? b) if an item loads strongly in both dimensions (where it appears that there's little separation) is that a problem? I didn't think it was but I thought it might be worthwhile to add. If anyone has any thoughts on these questions, or if an important procedure comes to mind that we haven't consider, I'd love to hear it. Thanks in advance, and I apologize if this post is not properly directed!


title: Sign of the unnormalized log likelihood in Ising model

text: Here is a section of _Machine Learning: a Probabilistic Perspective_ by Kevin Patrick Murphy ![enter image description here](http://i.stack.imgur.com/jatxV.png) I don't understand in (19.18) why there is a negative sign. For me, $\log \tilde{p}(\mathbf{y})=\sum_{s\sim t}\log\psi_{st}(y_s,y_t)$ holds. When $y_s$ and $y_t$ agree, $\log\psi_{st}(y_s,y_t)=w_{st}$, otherwise $\log\psi_{st}(y_s,y_t)=-w_{st}$. So shouldn't it be $\log \tilde{p}(\mathbf{y})=\sum_{s\sim t} y_sw_{st}y_t$? Also, when all entries of $\mathbf{y}$ agree, $\sum_{s\sim t} y_sw_{st}y_t$ is maximized because all summands are positive. So where is the problem? Do I miss something? Thank you.


title: Are $\mathbb{F}_2$-linear combinations of random variables in an i.i.d. Bernoulli process again an i.i.d. Bernoulli process?

text: I'm having trouble understanding how certain combinations of random variables can correlate. The problem is as follows: I have a binary $m \times n$ matrix $A$ of full rank (over the finite field $\mathbb{F}_2$ with two elements), where each row has exactly $w$ $1$'s (i.e., their weights are all $w$). I randomly pick an $n$-dimensional binary vector $\boldsymbol{v} = (v_0,\dots,v_{n-1}) \in \mathbb{F}_2^n$, where each $v_i$ is $1$ with probability $p$ and $0$ with probability $1-p$. The probability $p$ is assumed to be strictly smaller than $\frac{1}{2}$. Here, $v_i$ are all independent, so that picking $\boldsymbol{v}$ is seen as an i.i.d. Bernoulli process of probability $p$. I am going to hash this vector $\boldsymbol{v}$ into an $m$-dimensional vector $\boldsymbol{v}' = (v'_0,\dots,v'_{m-1})$ by taking $\boldsymbol{v}' = A\boldsymbol{v}^T$ over $\mathbb{F}_2$. So, if I write the $i$th row of $A$ as $\boldsymbol{r}_i$, each $v'_i$ is the product $v'_i = \boldsymbol{r}_i\cdot\boldsymbol{v}^T$ between the $i$th row of $A$ and $\boldsymbol{v}$. In other words, I'm just taking the mod $2$ sum of some bits of $\boldsymbol{v}$. Because the weight of $\boldsymbol{r}_i$ is assumed to be $w$, the probability $P(v'_i = 1)$ that the $i$th bit $v_i'$ of the hashed vector $\boldsymbol{v}'$ is $1$ is $$P(v'_i = 1) = \sum_{x: \text{odd}}^{w}\binom{w}{x}p^x(1-p)^{w-x},$$ which is a constant if we fix $p$ and $w$, and is also uniform across all $v_i'$. Now, if $A$ had linearly dependent rows, the hashed bits $v_i'$ are clearly not independent. My question is: _**What if all rows of $A$ are linearly independent? Is $\boldsymbol{v}'$ an i.i.d. Bernoulli process with probability $P(v'_i = 1)$? If not, can I approximate it as one if $\boldsymbol{r}_i$ and $\boldsymbol{r}_j$ have only few $1$'s at the same columns for all $i \not=j$?_** I remember I read in some research paper in electrical engineering that $\boldsymbol{v}'$ is i.i.d. Bernoulli if rows of $A$ are linearly independent, ~~although I can't seem to remember where it was.~~ (I found one recent paper that says this is the case: V. Toto-Zarasoa, A. Roumy, C. Guillemot, _Maximum likelihood BSC parameter estimation for the Slepian-Wolf problem,_ IEEE Comm. Lett. **15** (2011) 232–234. It's Lemma 1 on the second page.) But now I think about it, this is counterintuitive (to me) because if two rows of $A$ have $1$ at the same column, that means that I took at least one same $v_i$ when hashing $\boldsymbol{v}$ for those two rows, so the corresponding $v_i'$ look correlated. If my reasoning and Alecos Papadopoulos's answer are correct, the last part of the above question in boldface becomes essential. For instance, are there standard methods for evaluating how close or similar a given set of random variables like $\boldsymbol{v}'$ is to i.i.d. Bernoulli?


title: Regression with some associated ordinal dependent variables

text: I have some associated categorical dependent variables that are ordinal in nature (with 4 or 5 categories). If I want to see the effect of a set of independent variables (which can be both continuous and categorical) what is the right way to perform the regression? Should that be a multivariate generalized linear model? Actually if I performed separate regression for each DV, I would use ordinal logistic regression. But as the variables are associated, I am afraid, that will be misleading to use in my case! What should be the right procedure? Can SPSS handle this? A reference that I can follow is cordially expected.


title: Correlation coefficient between two arrays of 2D points?

text: I have two arrays of 2D points and I need to estimate their correlation. What formula should I use? Thank you. Example of arrays: X: [[1,5],[2,5],[1,7],[4,1]] Y: [[3,4],[1,6],[4,6],[4,3]]


title: eta-squared in GLM

text: If I am right, you can calculate eta-squared in an ANOVA by: SS (covariate of interest) divided by SS (total) I was wondering if the calculation and interpretation is the same in a General Linear Model (GLM).


title: Samples realizing conditioning on intersections of events

text: I have a data set that is essentially of the form of a sequence of pairs: $$\left(\cap T_i, P'(W|\cap T_i)\right)$$ The sequence might have duplicate $\cap T_i$'s with different $P'(W|\cap T_i)$'s (where $P'$ denotes an empirical probability). The events $T_i$ are not independent. I would like estimate the probability $$P(W|T_i)$$ for each of the $i$'s. Is this feasible? What "kind" of technique do I need to use?


title: machine learning applications in number theory

text: > is there any research into or applications of machine learning in number > theory? am also looking for (leading examples of) statistical/empirical analysis of number theory questions. also wondering if genetic algorithms in particular have ever been used in these areas.   * roughly related question on other site: Why can machine learning not recognize prime numbers?   * an area in number theory that seems to have had some statistical analysis, the Collatz conjecture.   * possibly somewhat related, automated theorem proving.


title: Can I use Wilcoxon Signed-Rank Test to test pairs of data in a time series?

text: I have a dataset of measurements (non-normally distributed). Each line of this dataset is composed by two measurements of the same variable in a specific time (i.e. XYY [time, meas1, meas2]). This dataset contains 24 lines. Is it correct to use Wilcoxon Signed-Rank Test to compare these two measurements? Is it correct to use this test to compare the mean value of these measurements between different experiments? Thanks


title: SVMlight support vectors inside the margin

text: I am trying to understand a linear model obtained with SVMlight. The model file contains a number of support vectors and their corresponding $\alpha_i y_i$. Since all $y_i \in \lbrace -1,+1 \rbrace$, I know that $\alpha_i = |\alpha_i y_i|$. From the constraints in the dual problem, it is given that $- \leq \alpha_i \leq C$ for all $i$. From the literature, I understand the following interpretation of $\alpha_i$ values:   * $\alpha_i = 0$ : a training instance has sufficient distance to the hyperplane (i.e., it lies beyond the margin); therefore these instances vanish from the model.   * $\alpha_i < C$ : a training instance lies on the margin, i.e., $y[w \cdot x - b] = 1$ for that instance.   * $\alpha_i = C$ : a training instance lies in the margin, i.e., $y[w \cdot x - b] = 1 - \xi_i$ with positive $\xi_i$ Furthermore, the function $f(x) = w \cdot x - b = \sum_i y_i \alpha_i (x_i \cdot x) - b$ gives a value that is proportional to the distance of $x$ to the margin. In my SVMlight model, I obtain values that seem inconsistent to me. If I calculate the $f(x)$ for all support vectors in the model, I obtain values in between $-1$ and $+1$. So all support vectors where $|f(x_i)| < 1$ have to lie inside the margin. Since this is only the case if $\alpha_i = C$, this implies that all of their $\alpha_i$ should be the same: $$ \forall i : |f(x_i)| < 1 \Rightarrow \alpha_i = C $$ However, this is not the case. In fact, there are many support vectors with $|f(x_i)| < 1$ and $\alpha_i < C$. Is there some severe misunderstanding in my interpretation of the model? Why are there support vectors in the margin with $\alpha_i < C$?


title: How to estimate true concentration level of antibodies according to optical density in ELISA?

text: In ELISA Is it necessary to obtain true concentration level of antibodies for diagnosis? Is there a model allowing us to do so? I know there is a 4-parameters logistic model, but will people conducting tests use it? How about other test, like NAA test?


title: Feature Importance in each fold and repeat after repeated cross validation in caret

text: this is my first post on Cross Validated so I apologize in advance if I'm not yet familiar with any conventions regarding forum posts. Currently, I'm working on a feature selection task using elastic net in caret and I would like to visualize the feature importance for each model trained in each cross validation step. However I don't seem to find a way to access the coefficients, other than those of the final model. Here is more or less a minimal example of what I'm doing.               library(caret)     library(doMC)          registerDoMC(16) # register 16 cores      Define tuning grid of elastic net parameters.               grid <- expand.grid(.lambda = seq(0, 1, length=20),                         .alpha = seq(0, 1, length = 11))      Provide control object defining repeated cross validation resulting in 5-fold cross validation each repeated 5 times.               ctrl <- trainControl(method = "repeatedcv",  # cross-validation method                          number = 5,             # number of folds                          repeats = 5,            # number of complete sets of folds                          allowParallel = TRUE)   # utilize parallelization      Train models with defined cross validation scheme and parameter grid Furthermore the featuers will be centered and scaled.               model <- train(x = iris[,-5],                    y = iris$Species,                    method = "glmnet",                    type.gaussian = "naive",                    tuneGrid = grid,                    trControl = ctrl,                    preProc = c("center", "scale"))      Alright, now I can get some information about the test performance after repeated cross validation.               model$resample[with(model$resample, order(Resample)), ]          Accuracy Kappa   Resample     12 1.0000000  1.00 Fold1.Rep1     19 1.0000000  1.00 Fold1.Rep2     25 1.0000000  1.00 Fold1.Rep3     2  0.9666667  0.95 Fold1.Rep4     9  0.9333333  0.90 Fold1.Rep5     1  0.9000000  0.85 Fold2.Rep1     8  0.9333333  0.90 Fold2.Rep2     15 0.9666667  0.95 Fold2.Rep3     22 0.9333333  0.90 Fold2.Rep4     16 0.9333333  0.90 Fold2.Rep5     18 1.0000000  1.00 Fold3.Rep1     11 0.9666667  0.95 Fold3.Rep2     5  1.0000000  1.00 Fold3.Rep3     3  1.0000000  1.00 Fold3.Rep4     6  1.0000000  1.00 Fold3.Rep5     23 1.0000000  1.00 Fold4.Rep1     7  1.0000000  1.00 Fold4.Rep2     10 1.0000000  1.00 Fold4.Rep3     17 0.9333333  0.90 Fold4.Rep4     24 1.0000000  1.00 Fold4.Rep5     13 0.9333333  0.90 Fold5.Rep1     20 0.9666667  0.95 Fold5.Rep2     14 0.9333333  0.90 Fold5.Rep3     21 1.0000000  1.00 Fold5.Rep4     4  1.0000000  1.00 Fold5.Rep5      However, I don't see how to access the coefficients for the models generating the respective cv accuracies, to visualize the variable importance in the same way it is possible for the final model.               plot(varImp(model))      ![Variable Importance for final model](http://i.stack.imgur.com/EsFUk.png) I would very much appreciate your help.


title: Chosing the reference group in a regression with a non-ordinal predictor

text: When running a model with a categorical predictor, the output contains n-1 p.values where n is the number of levels in the factor. One level is taken to be the reference group. For example: I observed plants in different environments and want to know if the substrate1, substrate2 and substrate3 (categorical non-ordinal predictors) influence the height, number of leaves and weight (response variables) of the plants. My question is: Which group should be the reference? Note: None of the group might be consider as a control group or as a normal group while the others are differ from what is normal. For example I will run the three models. Heigth ~ substrate1 + ... weigth ~ substrate1 + ... Number_of_leaves ~ substrate1 + ... `substrate1` is made of 5 groups: 'a','b','c','d','e'. It non-ordinal! Which group of `substrate1` would you suggest me to choose as the reference group for each model? - One of the two extreme group (even if it is not the same group in each model) - The group which is the closest to the overall mean (even if it is not the same group in each model) - The group which contain the greatest number of observations - Any group but it is better to always chose the same group Here are some graphs to make more sense of my question. Should I always choose group 'c' or should I chose 'c' for the first model, 'd' for the second model and 'e' for the last model. Or other... ![enter image description here](http://i.stack.imgur.com/MFddZ.png) ![enter image description here](http://i.stack.imgur.com/dZRRM.png) ![enter image description here](http://i.stack.imgur.com/GUy6g.png)


title: Calculate the likelihood of two consecutive events

text: I have a sequence of events for a number of users. Now I want to see if there is any correlation or a likelihood score between two consecutive events. For example, I have 10 types of events (A1,A2,A3,A4,A5,A6,A7,A8,A9,A10) and the time between two consecutive events. Now I want to find out how likely it is (or what is the correlation) that two events come after each over time to predict whether the two events are related to each other or not.


title: Should I compute both Pearson's and Spearman's coefficients of correlation?

text: I have been reading over some statistics notes, and certain things are introduced without motivation (note: I am a mathematician and not a statistician, so my stats knowledge is woefully inadequate, so I know little about the motivation). One of these oddly-introduced things is Spearman's correlation coefficient. So, I have been doing some reading, and I now understand that this coefficient measures how monotonic your data looks, while Pearson's measures how linear the data looks. And so on. I also found this old question from here, and one of the answers points out that if Spearman's coefficient has absolute value greater than the absolute value of Pearson's coefficient then you might be better looking for something monotonic and not just linear. However, Spearman's coefficient was introduced in these rather odd notes in a section on linear regression. I was wondering if it has any business being there? In the sense that, if your data is linear then it is monotonic. So if your data has Pearson's coefficient has absolute value basically equal to 1 is there any point in working out Spearman's coefficient? Even if Spearman's coefficient had absolute value 1 on the nose then...so what?


title: Analyzing differences in preference for proportional data

text: We have conducted an experiment using a two-alternative forced-choice task where listeners must indicate whether they prefer the sound quality of an auditory stimulus in interval A or interval B. These preference judgments are made for different types of signal processing, which are randomized. We have three signal processing conditions for which we would like to compare preferences (A, B, C). Our data end up being the proportion of the time that an individual preferred one type of processing over the other, so we have the proportions A~B, B~C, and A~C. We had considered doing an ANOVA with processing as a within-subject factor, but this does not seem appropriate given the dependence of preferences that occurs because of our forced-choice approach. Any suggestions are appreciated.


title: Repeating questions based on giving correct answers

text: I am making a web-app that is supposed to help people memorize Japanese kana symbols using flashcard. How it works at the moment is that on each new flashcard, user is presented with one symbol which gets chosen randomly from a pool of X symbols. User than has to recognize that symbol and if he is successful, correct count and total number of tries for that symbol both go up by one. If the user failed to recognize the symbol only the total number of tries goes up by one. I would like to make the symbols that user had lower success recognizing appear more frequently to make the learning more efficient. Is there a distribution, or a statistical method (using some kind of weight factors or something like that) that is used to do something like this which I could use? I don't even have a clue what to google for.. How do I properly model weight factors based on success rate (percentage)?


title: Grid search on k-fold cross validation

text: I've a dataset of 120 samples in a 10-fold cross validation setting. Currently, I pick the training data of the first holdout and do a 5-fold cross-validation on it to pick the values of gamma and C by grid search. I'm using SVM with RBF kernel. Since I'm doing a ten 10 cross-validation to report precision,recall, do I perform this grid search in the training data of each holdout (there are 10 holdouts, each having 10% test and 90% training data) ? Wouldn't that be too time consuming ? If I use the gamma and C of the first holdout and use it for the rest of 9 holdouts of the k-fold cross-validation, is that a violation because I would have used the train data to get gamma and C and again use the portion of train data as test in second holdout ?


title: Searching for a model like ZIP, but to model three data-generating processes, not just two

text: What type of test would fit my data in this situation? I'm measuring a count variable. Like in a zero-inflated poisson, I have numerous observations that are zeros because of a separate process. But I also have observations that have extremely high outcomes in the count variable because of a distinct process. Can I model a "zero-inflated & high-outcome inflated" poisson or something like it?


title: Calculating Bias from Cox PH model

text: I am using COX PH model to fit a lifetime data set and estimating parameters. By the by I am also simulating data and trying to find out estimates of parameters. Now I want to calculate bias and standard deviation to determine the size of the bias How can I do it? Can anybody give me suggestions and some articles name related to this problem?


title: Intuition behind auto/cross-correlation

text: I am really having trouble understanding the intuition behind autocorrelation. I mean why are we even calculating a correlation of some series with itself? Can anybody explain it in lay terms with an everyday example from life? Also, I see that while calculating the ACF, we have to take lags. But how does one know in advance how many lags to take so as to make them correlate better?


title: How do I set up a single continuous independent variable for a repeated measures linear model?

text: I'm doing an education study and I'm trying to see the effect of website usage on quiz scores in a class of college students. There were 30 students in the class and they took 10 weekly quizzes over the course of the semester. They also had a website where they could log in and use a study tool. The website kept track of the time they were online. So my dependent variable is quiz score (continuous) and my predictor variable is website usage over the previous seven days (also continuous). My first instinct was to make a 300-line database, matching up each quiz score with the amount of website usage in the preceding week, and then run a regression. But I can't really do that, can I? Instead of 300 independent measurements, I really have 10 repeated measurements on 30 subjects. So I tried using "General Linear Model > Repeated Measures" in SPSS, but I can't figure out how to tell the program that all those columns for website usage are a single, continuous predictor variable. Any guidance? Am I on the right track here? Or should I be using a different analysis altogether?


title: Is there a specific name for this formula?

text: I was wondering if there was a specific name for this formula and if it is involved in an established statistical process (like how chi-squared tests are for comparing observation classes or t tests are for distributions without a known population stddev), or if this is an original development from the referenced article author. I wish to understand its usage a bit more is why. How it is referenced in the article:   > A further test, based on the work of Ellegard and Thomson, was devised to > determine "Distinctive Numbers" ... The "distinctive numbers" which result > are therefore given the symbol D. This test is in fact finding where > x-epsilon is the expected occurrence, x-0 the observed occurrence and x the > average occurrence in the grand total. The effect of this is to increase the > denominator when the Stowe manuscript uses a word either more often or less > often than the average, giving a lower D value and weighting the result so > that increased notice is taken of occurrences when abnormalities in the > texts coincide. The formula in question: $$ D=\sum\frac{[(x_\epsilon-\bar{x}) - (x_0-\bar{x})]^2}{|x_\epsilon-\bar{x}|} $$ The article where I found it: Harvey, P. (1970). Stylistic Analysis


title: Given count data with many zero observations, what is a reasonable amount of zero observations in the data?

text: I have sales data which records at what time (by second) and how many were sold. Therefore, the data is count data with around 90-150 sales during a 3-day period. If I agggregated it to the 10 minute interval, I will have around 75% zero observations because for example 100/(24*3*6) = 0.2314815. But I think this is still leads to too high a zeros/total ratio even when applying the data to zero-inflated Poisson or similar types of models. What is a reasonable amount of zero observations contained in the dataset to apply zero- inflated Poisson like models? According to the literature or examples elsewhere, I think 40% zero observations is acceptable. But I don't know if 50% okay.


title: Can a deep belief network (stacked RBMS) be used solely as a dataset generator?

text: I have a large dataset (tens of thousands of predictors) on which I would like to perform feature reduction with the intent of better model-building for prediction. Deep Belief Networks seem to address the problem in such a way that they produce an output that bears some resemblance to that of PCA, i.e. a series of 'virtual' dimensions that can be used to explain variation in the data. Unlike PCA however, it seems that every example of DBNs follow the original example of Hinton, being used to generate a series of starting weights for a feed-forward neural network. Is there a reason that nobody appears to use stacked RBMs as, say, input for linear regression or k-nearest neighbor?


title: Inconsistency between 'read.csv' and 'write.csv' in R

text: The R function `read.csv` works as the following as stated in the manual: "If there is a header and the first row contains one fewer field than the number of columns, the first column in the input is used for the row names." That's good. However, when it comes to the function `write.csv`, I cannot find a way to write the csv file in a similar way. So, if I have a `file.txt` as below:               Column_1,Column_2     Row_1,2,3     Row_2,4,5      Then when I read it using `a = read.csv('file.txt')`, the row and column names are `Row_x` and `Column_x` as expected. However, when I write the matrix `a` to a csv file again, then what I get as a result from `write.csv('file2.txt', quote=F)` is as below:               ,Column_1,Column_2     Row_1,2,3     Row_2,4,5      So, there is a comma in the beginning of this file. And if I would read this file again using `a2 = read.csv('file2.txt')`, then resulting `a2` will not be the same as the previous matrix `a`. The row names of the matrix `a2` will not be `Row_x`. That's, I do not want a comma in the beginning of the file. How can I get rid of this comma while using `write.csv`?


title: When using linear regression analysis to get the fitted values of an outcome, why do the more extreme values tend to be predicted closer to the mean?

text: I am working on a project in which I am using several independent variables to "predict" the values of an outcome using linear regression. In R this is done quite simply as               model  <- lm(outcome ~ predictor1 + predictor2 + predictor3)     fitted <- model$fitted.values      I am interested in the _difference_ between the predicted values and the actual values - i.e. how accurate the predictors are.               residuals <- model$residuals      My question relates to the relationship between `residuals` and `outcome`. Samples with lower values of `outcome` tend to have negative values for `residuals`, and vice versa for samples with high `outcome` values. Plotting the values against one another is the simplest way to see this: ![Various plots](http://i.stack.imgur.com/U5dPa.png) The $R^2$ for the original LM (outcome ~ predictors) is 0.42, the $R^2$ between `residuals` and `outcome` is 0.58, and the $R^2$ between `fitted` and `outcome` is 0.39. What could explain the phenomenon? Why would samples with high `outcome` tend to be predicted lower than they actually are, and vice versa for lower values of `outcome`? Or indeed, am I missing something conceptually here? Many thanks for your input * * * Edited (13.08.20) to include an updated plots and terminology (now use "residuals" rather than "difference") - but in essence the questions remains the same. Thanks all for the input so far.


title: Representing repeated measurements within sample plots in linear mixed effects model in R

text: I have collected data on gas fluxes from plots of soil subjected to 5 different treatments ("D2", "K2", "M", "N", and "O2"), which also possessed variable clay contents. The experiment was laid out in a randomized complete block design, with 4 replications. Within each plot, two separate measurements of flux were performed. The resulting data.frame resembles the following:               block   treatment   subsample       flux        clay     1           D2          1           112068.6003 14.8     1           D2          2           129223.1641 14.8     1           K2          1           256712.4712 15.5     1           K2          2           113343.9756 15.5     1           M2          1           85794.47834 16.4     1           M2          2           -33620.6990 16.4     1           N           1           70283.98133 18.2     1           N           2           49569.84621 18.2     1           O2          1           100553.1116 13.4     1           O2          2           38885.99674 13.4     2           D2          1           96968.58451 15.8      I want to build a linear mixed effect model that takes account of this subsampling, and have come up with:               flux.lme <- lme(flux ~ treatment + block, random = ~1|subsample)      Which produces an ANOVA table:               > anova(flux.lme)                    numDF denDF   F-value p-value     (Intercept)        1    26 158.15781  <.0001     treatment          4    26   8.88691  0.0001     clay               1    26   1.72640  0.2003     block              3    26   1.59188  0.2153     treatment:clay     4    26   1.73011  0.1736      This output seems a bit strange to me, as the denominator degrees of freedom should not be 26, which is taking each repeated sampling as independent experimental unit. It should instead be based on the number of “main plots”, which is 20. In my case the denominator d.f. should be 20-1-3-1-4-4=7. Is is possible to instruct the lme() function to use this value?


title: Markov chains that do not contain all the states in the model

text: I am trying to understand Mixture Markov Models in order to cluster a set of sequences that do not necessarily all have the same states occurring in them. If I have several sequences that I am trying to fit a Mixture Markov Model to, how do I go about the fact that some of the chains contain a part of all possible states while the others might contain other states? For example: sequence 1: 1 2 4 2 4 1 sequence 2: 3 2 3 3 3 3 sequence 3: 1 2 3 1 1 2 How would I define a transition probability matrix for each of the sequences, given that the set of possible states is {1,2,3,4}? In the context of Mixture Markov Modeling, shouldn't a transition probability matrix be 4x4 and contain transfer probabilities for each state from the collection -> each state of the collection? Take for instance sequence1: My understanding is that the TPM of sequence1 should be:                  1   2   3   4          1  0.5 0.5 0   0          2  1   0   0   0          3  0   0   0   0          4  0.5 0.5 0   0      This seems to be incorrect, since TPM is a stochastic matrix whose rows sum up to 1, and the third row is just all 0s. Should I completely remove state 3 if I am building a transition probability matrix of a specific sequence? The reason I need to have a Markov chain for each sequence is that I am trying to calculate log-likelihood distance between each two sequences in order to build a distance matrix and get medoids for initializing the Mixture Markov Model. The formula for distance that I am trying to use is: $D(seq_i,seq_j)=1/2* log(p(seq_i|v_j)+p(seq_j|v_i))$, where $v_i,v_j$ are the Markov Chains representing the sequences $seq_i$ and $seq_j$ correspondingly. How would I define sequences 1,2,3 in terms of their transition probability matrices and single ML-estimated Markov chains representing each sequence?


title: Is logistic regression still valid where there are repeated measurements and generalized estimating equations is the real deal?

text: I have a model with a binary dependent variable (DV) and 5 independent variables, all of which are matched (each person, twice). I think since these matched INDEPENDENT variables can be considered "repeated- measures", _generalized estimating equations_ (GEE) is the best approach here. However, I have used binary logistic regression. Do you think this analysis is valid? Or is it just of lower power or less elegant than GEE? I can handle a lower power or less elegance, but not invalid. By "invalid" I mean where some assumptions are not met and the test is incorrect. My guess is that since binary logit discards the correlations between the repeated measures, it is a special case of repeated measures with zero correlation. So it might still be valid, but less useful than GEE. Am I right? Besides, I doubt if matching the independent variables is considered repeated measures. I am confused. Please guide me. Many thanks


title: How to compare two related ordinal variables?

text: I would like to compare two measurements of a variable (anxiety) on the same subjects at different times. I was planning on using the paired samples t-test in SPSS, but I believe the variables are ordinal, since the questions were scaled 0,1,2,3. How can I compare them?


title: behavioral data on contrast sensitivity

text: I have done behavioral experiments,and have no clue even how to start making statistics on the data. Here is the description of the experiment:   Mice are subjected to an experimental task where they have to distinguish stripe patterns of various contrasts and width. It is known that at a certain spatial frequency (stripe width) the contrast sensitivity is best, and with wider or narrower stripes the sensitivity becomes worse (see e.g. here) So: For every mouse (n=7) we tested at a **given stripe width** and determined the contrast at which it was just able to see the stripes. That is,we get 7 data points per width. We did this for 3 widths: the "optimal" width, and a wider and a narrower width. As result, we get some V-shaped curve (connecting the means for each stripe width).   Next, we tested at a **given contrast** (relatively high contrast, i.e. stripes can be seen well) and checked how narrow the stripes can be made. Again, we get 1 data point for each mouse. We did this for 2 different contrasts. Taken together, for each mouse we have 5 data points in the form {stripe_width,contrast}, where 3 times the width was fixed, and contrast may vary from mouse to mouse, and 2 times the contrast was fixed, and the width may vary from mouse to mouse. Now comes the tricky part: These mice are known to become worse as they become older. So we repeated the same test, with the same mice, 3 times within a couple of weeks. In the plot below, black is the original measurement, red after 1 week, green after 2 weeks. ![Mathematica graphics](http://i.stack.imgur.com/Ujrst.png) I am now looking for a statistical test to show that they really see worse with age (they do - I want to quantify it). Some even become blind, so that I do not have numerical values for {stripe_width,contrast}, but one could imagine (maybe) that one can set the needed contrast to Infinity (or some other high value compared to the real data) for the experimental tests where one keeps the widths fixed. The blind ones are not plotted in the graph above, and are not included in calculating the mean. I hope someone out there can help me. I very much appreciate it.


title: Spherical Gaussian Sigma dimension

text: I think I am confused with this thing. If we have a 3 dimension Gaussian then the MLE estimate for $\mu$ is a vector with 3 element $$\mu(1)' = \frac{1}{n}\sum_{j = 1} ^ n x_j\text{ and so on }\mu(2)\text{ for }y\text{ and }\mu(3)\text{ for }z$$ also $\sigma$ $$ \sigma(1)' = \sqrt{\frac{1}{n}\sum_{j = 1} ^ n(x_j - \mu')}\text{ and so on }\sigma(2)\text{ for }y\text{ and }\sigma(3)\text{ for }z$$ When we know that the $x,y,z$ are not independent, in which case the covarience matrix looks like this: $\sigma\cdot I$ where $\sigma$ is a number and $I$ an $3\times 3$ identity matrix? How can we compute $\sigma$? Thanks for the answers an


title: What is the best test for validating trends in time-series data?

text: I have annual temperature data from a variety of weather stations in the Caribbean and I want to be able to show statistically that the trends for each station are significant, either positive or negative or no trend. I can obviosly plot the data and fit a line, but I want to make a stronger argument than just eyeballing it. I was thinking about a runs's test or a Mann-Kendall test. Most everything I've read talks about ARIMA models, concernes about autocorrelation, etc. but I feel like those are for econometric predictions. I'm not trying to predict anything. I just want to be able show whether or not the increase or decrease in temperature is due to random chance, or not. (my n=49) I may have larger n values for other stations, but for right now, working with 1 station, I have 49 observations. It would be nice to be give some idea of how far above or below random chance it is. I don't know this helps, but I actually have monthly data, but to avoid periodicity issues I've so far been playing around with July values to represent N. Hemisphere summers. In terms of hypothesis testing, if the Caribbean is not warming, then none of the weather stations would show increasing trends in temperature.


title: How is the intercept calculated in a generalized linear model and why is it different from a linear model?

text: I have a set of count data that I have fitted both a linear model to and a Poisson generalized linear model. The mean of the raw data is 233.375 and the standard deviation 279.983. I have been surprised when I fit a Poisson glm in R that the intercept = 5.53865. A negative bionomial glm in R gives an intercept of 5.4526. Is this because the data are sparse and/or over-dispersed? The data and R code are:               data <-c(2, 25,1121,361,251,123,123,81,25,215,4,196,0,353,968,336,179,229,92,204,35,299,8,371)     lm(data ~ 1)     glm(data ~ 1, family = poisson)     glm.nb(data ~1, link = log)


title: What is the correct sample size when calculating CI for a subset of results?

text: Suppose I have a study of 3031 people that obtains responses to various questions (95% CL). One of the questions (Q1a) gets a yes response from 616 of the people. Of those 616, only 34 have a QualityB (the remaining 582 don't) that is dependent on a yes answer to (Q1a). When determining the margin of error for QualityB (Confidence Interval) using an online calculator, do I use 34 or 616 for the sample size? I think I should use 616 (and 3031 population, 5% percentage), but I am not sure. Result +/- 1.54 Using 34 for sample size result is +/- 7.29 http://www.surveysystem.com/sscalc.htm


title: Omitted variable bias in linear regression

text: I have a philosophical question regarding omitted variable bias. We have the typical regression model (population model) $$ Y= \beta_0 + \beta_1X_1 + ... + \beta_nX_n + \upsilon, $$ where the samples are from $(Y,X_1,...,X_n)$, and then a bunch of conditions by which the OLS estimates behave quite well. Then we know that, if we omit one of the main variables, $X_k$, this might bias the estimates of $\beta_0, \beta_1, ..., \beta_{k-1}, \beta_{k+1}, ..., \beta_n$. This would affect, at least, the estimated effect of the rest of the variables on $Y$, and also the hypothesis tests about $beta_1,beta_2, ...$, as the predicted values are not reliable. The thing is, we don't know which variables are in the true population model. Instead, we have a bunch of candidates from which we should analyze and find out the most appropriate subset. This process of variable selection uses OLS estimates and hypothesis tests again. Based on that, we reject or include different variables. But since each candidate model is omitting relevant variables (you will never be able to find the true model), wouldn't these decisions be based on biased results? Why then, should we trust them? (I'm thinking of forward stepwise method, for instance, where you pick 1 variable then add the rest. You compare the models doing inference, and I'm thinking that omitted variables may be disturbing everything.) I was never too worried about this topic until I began thinking of it, and I'm sure I'm wrong somewhere.


title: Understanding the whiskers of a boxplot

text: I have a question regarding the interpretation of the whiskers of a boxplot. I have read the following: "On the top and bottom of the rectangle, the “whiskers” show the range of 1.5 times the distance between the 0.25- and 0.75-quantiles"¨, but do not entirely understand what is meant with "distance". It can't be that the probability mass is meant, since between the 0.25 and 0.75 quantile we obviously have always the same percentage of data. What then is exactly the idea?


title: Which "mean" to use and when?

text: So we have arithmetic mean (AM), geometric mean (GM) and harmonic mean (HM). Their mathematical formulation is also well known along with their associated stereotypical examples (e.g., Harmonic mean and it's application to 'speed' related problems). However, a question that has always intrigued me is "how do I decide which mean is the most appropriate to use in a given context?" There must be _at least some rule of thumb_ to help understand the applicability and yet the most common answer I've come across is: _"It depends"_ (but on what?). This may seem to be a rather trivial question but even high-school texts failed to explain this -- they only provide mathematical definitions! I prefer an English explanation over a mathematical one -- simple test would be "would your mom/child understand it?"


title: R knn variable selection

text: I have a data set that's 200k rows X 50 columns. I'm trying to use a `knn` model on it but there is huge variance in performance depending on which variables are used (i.e., `rsqd` ranges from .01 (using all variables) to .98 (using only 5 variables)). This kind of compounds my problem as now I need to determine `k` _and_ which variables to use. Is there a package in R that helps with selecting variables for a `knn` model, while tuning `k`? I've looked at `rfe()` in `caret` but it seems to only be built for linear regression, `randomforest`, naive bayes, etc but no `knn`. As an aside, I've tried manually building a loop to use the caret train function like this:               for(i in 2:50){     knnFit <- train(x[,i],y,...) ## trains model using single variable     }      My problem is that `knnFit$results` prints all of the results and `knnFit$bestTune` only prints the final parameter of `k`.               > data1 <- data.frame(col1=runif(20), col2=runif(20), col3=runif(20), col4=runif(20), col5=runif(20))     > bootControl <- trainControl(number = 1)     > knnGrid <- expand.grid(.k=c(2:5))     > set.seed(2)     > knnFit1 <- train(data1[,-c(1)], data1[,1]     + , method = "knn", trControl = bootControl, verbose = FALSE,     + tuneGrid = knnGrid )     > knnFit1      20 samples      4 predictors          No pre-processing     Resampling: Bootstrap (1 reps)           Summary of sample sizes: 20           Resampling results across tuning parameters:            k  RMSE   Rsquared       2  0.485  0.124          3  0.54   0.369          4  0.52   0.241          5  0.528  0.232             RMSE was used to select the optimal model using  the smallest value.     The final value used for the model was k = 2.           > knnFit1$results           k      RMSE  Rsquared RMSESD RsquaredSD         1 2 0.4845428 0.1241031     NA         NA         2 3 0.5401009 0.3690569     NA         NA         3 4 0.5197262 0.2410814     NA         NA         4 5 0.5277939 0.2317607     NA         NA          > knnFit1$bestTune           .k         1  2      I need some way to print the RMSE/rsqd/other metric for the best single performing model (i.e., just "R-Squared: .91"). Any suggestions?


title: Error propagation with linear regression

text: I'm trying to obtain an estimation of the uncertainty related to an analytical method: my function is just a linear regression $f: y=ax+b+ \epsilon$ with $y_i=\frac{R_i}{C}$, both $R$ and $C$ are affected by a measured uncertainty, $s_R$ and $s_C$ respectively. I want to calculate the total uncertainty with the generic formula for error propagation $$s_f^2=\sum_{i=1}^{n}\Big[{\Big(\frac{\partial f}{\partial x_i}\Big)s_{x_i}}\Big]^2$$ My $x_i$ are $R_i$ and $C$, so $$s^2_f=\sum_{i=1}^{n}\Big[{\Big(\frac{\partial f}{\partial R_i}\Big)s_{R_i}}\Big]^2+\Big[\Big(\frac{\partial f}{\partial C}\Big)s_{C}\Big]^2$$ with $$\Big(\frac{\partial f}{\partial R_i}\Big)^2= \Big(\frac{\partial a}{\partial R_i}\Big)^2+\Big(\frac{\partial b}{\partial R_i}\Big)^2+\Big(\frac{\partial \epsilon}{\partial R_i}\Big)^2$$ and a similar expression for $\Big(\frac{\partial f}{\partial C}\Big)^2$. The problem for which I need your help is the calculation of $\Big(\frac{\partial \epsilon}{\partial R_i}\Big)^2$. Trying to solve the problem without taking into account the error term, lead me to the solution $$s^2_f=\frac{1}{n}\sqrt{s^2_R+n^2\bar y^2 s_C^2}$$ so, an hyperbolic branch that, increasing $n$, approaches asymptotically to $0$. I think it has sense but it doesn't take into account the probabilistic part. What do you think about it? Thank you.


title: When building a regression model using separate modeling/validation sets, is it appropriate to "recirculate" the validation data?

text: Suppose I've got an 80/20 split between modeling/validation observations. I've fit a model to the modeling data set, and I'm comfortable with the error that I'm seeing on the validation data set. Before I roll out my model for scoring future observations, is it appropriate to combine the validation back with the modeling data to get updated parameter estimates on the 100% data? I've heard two perspectives on this:   1. The validation that I've performed was a validation of the model structure, ie the set of predictors and transformations that I've applied. It doesn't make sense to leave 20% of my data on the table when I could use that to update my estimates.   2. The validation that I've performed was partly a validation of the parameter estimates that I calculated on my modeling data set. If I update the model fit, I've changed the estimates and I no longer have an objective way to test the performance of the updated model. I've always followed argument #1, but more recently I've heard several people argue #2. I wanted to see what others think about this. Have you seen any good discussions in the literature or elsewhere on this topic?


title: Fit Linear Model Using Generalized Least Squares with Unknown Covariance Matrix

text: I want to fit a linear model of $y$ on $x$ with respect to a grouping factor and with respect to a covariance structure for the dependent variable within groups. I have a dataset similar to the following:               y1 <- rnorm(25)     y2 <- y1 + rnorm(25)     y3 <- y2 + rnorm(25)     y4 <- y3 + rnorm(25)     y <- c(rbind(y1, y2, y3, y4)) # this is vector of the form: y1[1], y2[1], .., y4[1], y1[2], .., y4[2], ..., y1[25], .., y4[25]     x <- runif(100)     # Dataset:     df <- data.frame(y = y, x = x, group = factor(rep(1:25, rep(4,25))))       I need to estimate the coefficients for the $x$-variable in each _group_ (these effects are nested in _group_ and they are not random) and to do statistical inference for these coefficients. Of course, variables in each _group_ are dependent:               cor(cbind(y1, y2, y3, y4))      so I should use a method that considers the latter (unknown) covariance structure of the dependent variable. I found the function _gls_ in the _nlme_ library which is supposed to do that:               library(nlme)     # covariance structure:     cs <- corSymm(form = ~ 1 | group) # or maybe something different than group?     cs <- Initialize(cs, data = df)     # variance structure:     df$group2 <- factor(rep(1:4, 25)) # We make a new grouping factor for our data. Is that righr?     vs <- varIdent(form = ~ 1 | group2 ) # or maybe something different instead of factor(rep(1:4, 25))?     vs <- Initialize(vs, data = df)      But I am not really sure if the grouping factors for _cs_ and _vs_ is right.               # Model:     gls(y ~ x %in% group, correlation = cs, weights = vs, data = df)       The help for _gls_ is very poor and I don't have the book from Pinheiro, Bates. So I am not really sure if this is doing what I want. Also I don't really understand why the estimated covariance structure by _gls_ is so different from _cor(cbind(y1, y2, y3, y4))_.


title: Weibull MLE: what is the method/algorithm used to perform the optimization?

text: Regarding fitting a Weibull two parameters probability density function (pdf), with the `fitdistr` from MASS package, in R: 1 - What is the iteration algorithm/method used for scale and shape Maximum Likelihood Estimation (MLE)? I used to think it was Newton-Raphson, but after looking here and here I got confused. There are Nelder–Mead, quasi-Newton and conjugate-gradient algorithms, but I am not sure which one is used for the iteration process on Weibull 2 parameters pdf MLE. 2 - Based on the answer to question 1, what would be a good reference to study the method?


title: Generalization of Chernoff bound to Bernoulli schemes

text: I have this (I think) simple problem, which I guess is just a generalization of the Chernoff bound to Bernoulli schemes. Assume I have a bowl with N balls with R different colors. The balls are not necessarily evenly distributed. How many balls do I need to sample (with replacement) in order to decide with confidence % C which color occurs most often in the bowl? Well, that depends on the distribution of course. Let us assume that $0 < \text{prop(most_often)} \leq 1$ denotes the proportion of the most often occuring color. Now, for $R = 2$ (and thus $\text{prop(most_often)} \geq \frac{1}{2}$) we can plug these values into the Chernoff bound. But can that be generalized to $R \geq 3$?


title: Adstock v/s Polynomial Distributed lag

text: I am running a Market Mix model. There is a general practice in the industry to create an ad stock variable for media variables (like TV). I came across polynomial distributed lag model. I was thinking to used polynomial distributed lag instead of creating Adstock variables. Can anyone tell me whether should I move ahead with polynomial distributed lag or to continue with adstock variables?


title: Help with understanding inferential concept

text: I am reading an old article on Information Value Theory and to properly understand the material, I need to get my head around the following inferential concept where **_u_** , **_v_** are random variables and **_S_** is the state of information. > A most important inferential concept is expansion; it allows us to encode > our knowledge in a problem in the most convenient form. The concept of > expansion permits us to introduce a new consideration into the problem. > Suppose, for example, that we must assign a probability distribution to a > random variable _u_. We may find it easier to assign a probability > distribution to _u_ if we had previously specified the value of another > random variable _v_ and we may also find it easy to assign a probability > distribution directly to _v_. In this case the expansion equation, > > { _u|S_ } = ∫ _v_ { _u|vS_ } { _v|S_ } > > shows that all we have to do to find { _u|S_ } is multiply { _u|vS_ } by { > _v|S_ } and sum over all possible values of the random variable _v_. I will be grateful if someone could explain what { _u|vS_ } means. I am no good statistician and will appreciate an intuitive explanation on a concrete example.


title: Discriminative learning and generative learning

text: I'm wondering why in Generative learning algorithm, they try to maximize the probability $\prod_{i=1}^np(x^{(i)}, y^{(i)})$ While in Discriminative learning algorithm, it is $\prod_{i=1}^np(y^{(i)} | x^{(i)})$   1. What is the reason of maximizing these terms?   2. Why are these optimization objectives different? What is the meaning of this difference?


title: Discrete scenarios from an ARMA forecast

text: Given an ARMA model and a historical time series, I'm trying to create a set of $n$ forecast scenarios, where each scenario $s$ is a potential future hourly time series $x_s[t]$, with a given probability, $p_s$. I'm familiar with how to generate an ARMA point forecast for a given horizon - how do you extend this theory to forecast a series?


title: Training models on data that may be incorrectly classified?

text: I am working with some data that has been classified by domain experts. However, the classification they use is not 100% accurate. How can I deal with data that may not be correctly classified? Are there any models which are more resilient to data that are not classified with 100% accuracy? Does such a classification bias entail limitations on the performance of a trained model?


title: How to organize the data for this repeated-measure within-subject setup?

text: This question is a follow up to the question I asked yesterday. Here is the current situation again. I have data for a within subject measurement, where each variable is measured repeatedly. The repetitions are used to cancel out random fluctuation of the measured variable. Each subject was exposed to each of the conditions, so I have multiple datapoints per condition and subject. Currently the data is organized like this:               Subject; Cond1_rep1; Cond1_rep2; ... Cond1_rep20; Cond1_mean; Cond2_rep1; Cond2_rep2; ... Cond2_rep20; Cond2_mean ... Cond8_rep1; Cond8_rep2 ... Con8_rep20; Cond8_mean     1      ; ....     2      ; ....      So I have eight Conditions, and around 20 repetitions per condition per subject. Currently there are about 20 subjects to be included in the analysis (no removal of outliers and other data cleanup is done yet). Now for analysis the first attempt would be to do a multivariate ANOVA on all the means for all different conditions to see if there is any difference in the distributions. However if I use means like this I will use the fact, that each value probably has a very low error and the likelihood of a false result is actualy lower than estimated by the ANOVA. So I should somehow include the repetitions in the analysis. The suggestion I was given in the previous question was to use a generalized linear model with repetitions and include all the repetitions in the analysis. However in this case I would need an additional variable identifying the conditions to be used as an independent variable in the generalized linear model. So I would have to somehow add this in my data. What I would try would be first to split up the eight conditions. The eight conditions come from variations of three independent factors, so I can easily split them up. Then I would separate my data into multiple rows, so I can introduce a condition variable. This way I would organize like this:               Subject; Factor1; Factor2; Factor3; Rep1; Rep2; ... Rep20     1      ; 0      ; 0      ; 0      ; ...     1      ; 1      ; 0      ; 0      ; ...     1      ; 0      ; 1      ; 0      ; ...     1      ; 1      ; 1      ; 0      ; ...     ...     1      ; 1      ; 1      ; 1      ; ...     2      ; 0      ; 0      ; 0      ; ...     2      ; 1      ; 0      ; 0      ; ...      This way I would get eight rows per subject with a different combination of factors. In this case I can do a repeated measures GLM with the factors as independent variables and the 20 repeated measurements as repeated dependent variables. Would this setup be correct or are there any drawbacks? If there are, is there any better way of analysing this data without loosing the fact that each measurement was repeated multiple time (as when you take the mean first). **EDIT** : I am not so much interested in a way to restructure the data like this, but my main question if there could be a problem when I do a GLM like this, without taking into account that some data points where gathered from the same person. The final analysis in SPSS should not be a problem.


title: Interpretation of GEE coefficients

text: Suppose blood pressure is a continuous outcome variable and you run a linear GEE with following predictors: age (years), weight (lbs), and smoking (yes/no). How would you interpret the coefficients for these predictors? Would it be the same as in "regular" linear regression?


title: Optimization of an objective function with noise

text: I want to estimate the parameters that minimize a multivariate function with noise based on realized observations. As a simple example, say I have observed the univariate sample below: ![enter image description here](http://i.stack.imgur.com/drzJH.png) The samples were generated from a quadratic function plus noise, but in general the function will be unknown and I don't want to make any assumptions about it other than it is "smooth". In this case, the function at $x \approx 5.8$ happens to be the minimum, but I want to make a statistical argument that the function is minimized at $x=5$. How should I approach this? Should I fit some kind of nonparametric model to the sample data, and then optimize the fitted, smoothed values? Or are there algorithms that tackle the optimization of noisy function directly?


title: Should I randomize survey questions if a respondent will fill in identically structured forms multiple times?

text: I'm having an argument with a colleague (neither of us have a strong statistical / research methods background) on whether we should randomize questionnaire layout. The survey in question is a feedback forms on teacher performance which are distributed for high school pupils. So one pupil gets up to 20 surveys (on each teacher he has contact with) which have identical, multiple-choice questions. Most of the questions also have identical scale. The pupil fills the surveys in online and does not necessarily need to do everything at once (but he definitely can). We do not have control how / when they choose to do it. My colleague argues that by randomizing the order of questions we reduce the possibility that pupils will fill in the questionnaires without thinking too much - e.g. if one fills 5 forms at once, the evaluations of teachers will be very similar. In my opinion, though, we would actually gain from making sure that the survey questions are grouped together based on the area of question (e.g. teacher communication skills, usage of IT, etc) - this would ensure that pupils provide consistent answers. I guess there's a point in both approaches - but maybe there are some best practices how to deal with that? Oh, and if somebody could dig up some references, that would be ultra-cool. My colleague likes facts :)


title: Two years of data describing occurence of violence- testing association with number of patients on ward

text: I have two years of data which looks basically like this: Date **___ __ _**Violence Y/N? __ Number of patients 1/1/2008 **___ ___ _**0 **___ ___ __ ___ ____** 11 2/1/2008 **___ __ ___** 0 **___ __ ___ __ ___ __** 11 3/1/2008 **___ ____** _1 **___ ___ __ ___ ____** 12 4/1/2008 **___ ____** _0 **___ ___ __ ___ ____** 12 ... 31/12/2009_ ** ___ __** 0_ ** ___ __ ___ __ ___ __** 14 i.e. two years of observations, one per day, of a psychiatric ward, which indicate whether there was a violence incident on that day (1 is yes, 0 no) as well as the number of patients on the ward. The hypothesis that we wish to test is that more patients on the ward is associated with an increased probability of violence on the ward. We realise, of course, that we will have to adjust for the fact that when there are more patients on the ward, violence is more likely because there are just more of them- we are interested in whether each individual’s probability of violence goes up when there are more patients on the ward. I've seen several papers which just use logistic regression, but I think that is wrong because there is an autoregressive structure (although, looking at the autocorrelation function, it doesn’t get above .1 at any lag, although this is above the “significant” blue dashed line that R draws for me). Just to make things more complicated, I can if I wish to break down the results into individual patients, so the data would look just as it does above, except I would have the data for each patient, 1/1/2008, 2/1/2008 etc. and an ID code going down the side so the data would show the whole history of incidents for each patient separately (although not all patients are present for all days, not sure whether that matters). I would like to use lme4 in R to model the autoregressive structure within each patient, but some Googling comes up with the quotation “lme4 is not set up to deal with autoregressive structures”. Even if it were, I’m not sure I grasp how to write the code anyway. Just in case anyone notices, I asked a question like this a while ago, they are different datasets with different problems, although actually solving this problem will help with that one (someone suggested I use mixed methods previously, but this autoregression thing has made me unsure how to do this). So I’m a bit stuck and lost to be honest. Any help gratefully received!


title: Automatic test measuring dissimilarity between two time series

text: I have two time series, series1 and series2. My aim is to find **how much Series2 is different from Series1, automatically/quantitatively**.   Image can be seen in its original size by clicking here. ![alt text](http://i.stack.imgur.com/0N5fN.gif) Series1 is the expected result. Series2 is the test/incoming series. I am providing a histogram plot, where Series2 is represented in dark brown colour. You can also note in the x-axis between 221 and 353 there is a significant variation. ie Series2 is less than Series1. I am coding using C++. I think, crosscorrelation will help, but produces a value based on similarity rather than dissimilarity. I see people talk about Kolmogorov-Smirnov Test. Is this the test which i should be performing? UPDATE 1: I am trying to perform a template matching. I have divided my template image in to 8x8 blocks as well as my incoming test image. I am trying to compare one block in template image with the same block(based on the spatial pixel positions) in the test image. I calculate the intensity sum within each block. I obtain series1 for the Template image and have Series2 for the test image.


title: Confusion related to conditional Gaussian distribution

text: I have a certain confusion. I refer to this paper. Let's say I have $p$ variables $x_1, x_2, \dots, x_p$ which follow a multivariate Gaussian distribution. Now suppose I have $N$ examples or samples from this distribution. However, some of the values in each example are missing. Suppose I want to impute the values. Let's say using only the available examples I calculate the $\mu$ and $\Sigma$ vector which are the mean and covariances (initial estimate). Let $K = \Sigma^{-1}$. Now I found this formula $E[x_{ij}|x_{obs,i}, \mu,K]$ = $x_{ij}$ if $x_{ij}$ is observed $=c_j$ if $x_{ij}$ is not observed where $c = \mu-(K_{mis,mis})^{-1}*K_{mis,obs}(x_{obs,i}-\mu_{obs})$ $E[x_{ij}x_{ij'}|x_{obs,i},\mu,K]$ = $(K_{mis,mis})^{-1}_{jj'} + c_jc_{j'}$ if both $x_{ij}$ and $x_{ij'}$ are unobserved. However, $E[x_{ij}x_{ij'}|x_{obs,i},\mu,K]$ = $x_{ij}x_{ij'}$ if both $x_{ij}$ and $x_{ij'}$ are observed. I didn't get how the term $(K_{mis,mis})^{-1}_{jj'}$ disappears if either $x_{ij}$ or $x_{ij'}$ is observed.


title: Forecasting unemployment rate with plm

text: Please see question for the background. Following the advice of @kwak and @Andy W, I have decided to use the package **plm** in `R` to fit my model. Here an excerpt of the data `df` (the numbers are made up, not real data!):                   reg year          ur         mur cl     1    1 2001 0.000698717 0.012483361  1     2    2 2001 0.008283762 0.011899896  1     3    3 2001 0.001863817 0.012393738  1     4    4 2001 0.005344206 0.012126016  1     5    5 2001 0.007475083 0.011962102  1     6    6 2001 0.002785111 0.012322869  1      where **reg** is the region indicator, **year** is the year of measurement, **ur** is the unemployment rate, **mur** is the average unemployment rate in the neighboring regions (given by **cl** ) excluding the current region (see @kwak suggestion). Below is the code that I used in `R`               fm <- pgmm(log(ur)~lag(log(ur),1)+log(mur)|lag(log(ur),2:40),data=df)      I have several question regarding this model:   * I guess I should choose `effect="individual"` to avoid having time dependent intercepts (?). But doing so the code is crashing (Error in names(coefficients) <- c(namesX, namest))!   * Which `model` should I choose, `onestep` or `twosteps`?   * How do I decide on the numbers of instruments (40 is just a guess)? Assuming that I have fitted this model successfully, how do I simulate from it? Regards


title: Calculating standard deviation associated with percentage change in left ventricular mass

text: I have two averages +/- strandard deviations of left ventricular mass before and after of physical activity. Before physical activity, left ventricular mass is 150 +/- 20 grams After physical activity, left ventricular mass is 120 +/- 10 grams The left ventricular mass percentage variation, after physical activity is: (120-150)/150*100= - 20% I want to find out the standard deviation of a percentage change of left ventricular mass after physical activity compared to left ventricular mass before physical activity. Now my question is, how to take into account the standard deviations for this percentage value? You can write the formula with the numbers in my example?


title: How can I format and export PCA output in Stata?

text: Are there any tools that can format and write to external file output of principal component analysis in Stata? I'm thinking about something that will work in similar vein to [excellent] family of estout commands.


title: Why squaring $R$ gives explained variance?

text: This may be a basic question, but I was wondering why an $R$ value in a regression model can simply be squared to give a figure of explained variance? I understand that $R$ coefficient can give the strength of a relationship, but I don't understand how simply squaring this value gives a measure of explained variance. Any easy explanation of this? Thanks very much for helping with this!


title: Comparing two sets of pixels to determine whether they belong to the same object

text: I have two sets of data, and I want to know if the second set is sufficiently different from the first to be considered different. More specifically, I have a sample set A from a number of pixels in the vicinity of point X in an image, and another sample B of pixels in the vicinity of point Y in the same image. I want to know if point X and point Y could be part of the same object in the image (based on color values). For example, if A and B are both mostly blue then the answer is yes, but if A is red and pink and B is blue then the answer is no. My only idea so far (based on vague memories of a statistics class I took years ago) is to calculate the standard deviation in A and B, and from that calculate a threshold that gives you 95% certainty, and see if they are within a certain distance of each other. Is this correct? Otherwise, what's the best way to do this?


title: important features in RandomForest - Sklearn

text: 1) How to find important features in RandomForest classifier (in sklearn) with high statistical significant? 2) The input data I have is unbalanced which I simply repeat data to compensate that. When I do shuffling, some of these repeated data go to training and some goes to test which definitely increase the prediction accuracy. I know in reality it is not correct but what about if I only want to find important feature? 3) What is the meaning of a feature value (for example, 0.05) in randomforest? in Feature importance it is said that this value means 5% of data is classified correctly by this feature? To find the most important features can I sum up their values till become 0.9 and then I say 90% of data is classified correctly by this features?


title: Using a survival model to predict two different types of "death"

text: I am using a survival model to predict membership counts. Other than predicting signups that are being added, we are trying to use a coxph model to predict accounts that will get canceled. In our problem, there are two types of "deaths" that can happen. The company can terminate the account because the customer didn't pay, or the customer can cancel his account on his own. Right now, we are building two separate models to predict the two types of death - it is important for us to predict how many of each type there will be. People who had their account terminated are right censored in the account canceled model, and people who canceled their account are right censored from the account terminated model. To do the prediction, we want to use the hazard function to model the probability a person will cancel - or terminate - on day $t$, given that they have stayed with our service for $n$ days. Since these two models are not independent - if we predict your account will be terminated, then you cannot be canceled - how can I predict the number of terminated accounts and number of canceled accounts? Edit regarding confusion between censoring and "death" In this problem, a customer may have forgotten to pay and his account gets terminated. In the account terminated model, we would count this as a death, but in the account canceled model we right censored this customer. This may not have been the right thing to do, but it is easy to change the analysis so please suggest improvements if you have any. I believed that he should be right censored because maybe the customer did intend to cancel his account on his own accord later down the line, but we wouldn't know because his account simply expired and disappeared from our system. Therefore, if a customer did cancel on his own accord, we counted that as a death in the account cancel model, if his account expired we right censored the account, and if the account has not been canceled at all and is still with us as of today, we right censored it as well


title: Non negative lasso implementation in R

text: I am looking for some open source or an existing library I can use. As far as I tell the glmnet package isn't very easily extensible to cover the non negative case. I may be wrong, Any one with any ideas much appreciated. By non-negative I mean that all coefficients are constrained to be positive (> 0).


title: Random forest like techniques (bagging, random feature subset) for SGD methods

text: Are there any well-known results/tools/literature on using bagging and random feature subset selection for regression or SGD-based methods?


title: Linear mixed effects vs. Repeated Measures

text: I have 46 replicate trees I sampled in a 2 year study. My objective is: To estimate the proportion of variation among-trees w/in year, and proportion variation due to among year and tree interaction. Response: 12, 15, 32, 42, 13, 26, 10, 15 Year: A, A, A, A, B, B, B, B, Tree: 1, 2, 3, 4, 1, 2, 3, 4 The above is an example of my dataset I planned on running a repeated measures ANOVA (because I sampled the same tree in 2 years) but because circularity was violated, I decided to go with a mixed effects model. With trees as random and year as fixed. I'm then running a variance components analysis. Q1. Is the Mixed Effect Model a good choose to answer my initial questions? My model is this (in R code): lme <-lme(response~year,random=~1|tree/year,df) Also, I am not getting p-values for the interaction term after running Linear Mixed Effects Model (nor do I when I run a normal aov (2 way ANOVA) model using tree*year as a factor.) Wh Q2. A variance comp. analysis says that most of the variation (55%) is due to year and tree interaction. So why am I not getting the p-values for the interaction term? Thanks!!


title: Computing and comparing semantic networks generated by grouped individuals

text: I have 2 groups, each containing 20 individuals. For each individual I have a dissimilarity matrix of a number of judged items. I'd like to evaluate whether there is an overall difference between the groups in the semantic networks underlying their judgements. As a first step, I computed 2D MDS spaces for each group separately using the `smacofIndDiff()` function from the `smacof` package in R, which I understand permits one to characterize the overall space for a group of individuals who may vary slightly from one to another. This yields two 2D spaces (really, two 2xK matricies, where K is the number of items evaluated for dissimilarity), but I'm not sure how to proceed with regards to evaluating the degree to which the resulting spaces differ and whether this difference might reasonably be expected by chance. Any suggestions? Note that I'm not wed to the 2D MDS idea, so feel free to suggest alternative ways to represent the judgement data as a semantic space.


title: High-dimensional Regression Datasets

text: Am looking for pointers to publicly(online) available high-dimensional regression datasets for evaluating my research work. By high-dimensional, am looking for regression datsets with the number of predictor variables (features) ranging anywhere from 50 to 2000. Most of the regression datasets on statLib(CMU) and UCI-ML repository have datasets with the number of predictor variables $<$ 50.


title: Compare the difference of two probabilities or a ratio of probabilities?

text: In an experiment with continuous-value random variable, e.g. measuring length of cucumbers. I would like to compare the probability of getting a particular length range in two different conditions. e.g. $P1=P(0.1<L<0.2, fertilizer A)$, $P2=P(0.1<L<0.2, fertilizer B)$. So I create 2 histograms for the lengths corresponding to each of the 2 fertilizers, $H_A$ and $H_B$. From this histograms, probabilities over ranges of L are calculated from bin counts, etc. My aim is to show that only for some ranges (or even 1 range) of L the probabilities (or counts of the above histogram) differ while in most other length ranges the probabilities do not differ. 1) what is a good way to express this difference? The usual ratio $\frac{P1}{P2}$ suffers when $P1$ or $P2$ is zero. a good one is $P1-P2$ but then I would like this to be normalised over P so as to have something like a percentage change for comparison, I thought something along: $\frac{P1-P2}{P1+P2}$ 2) I would like to plot these probability differences for various L and if possible to use a statistical test which tests the significance of probability differences for just one range and not overall (what a t-test does). I guess the significance of difference depends also on the counts for that particular range of lengths.


title: Formula for marginal probability in CRF++

text: On the website for CRF++ http://crfpp.sourceforge.net/ they mention that marginal probabilities can be output for each possible label. My question is, in CRF theory, what's the formula for this marginal probability? I think it's: $$p_{_\theta}(Y_{_j} = y | X ) = \frac{\alpha_{_j}(y|X)\beta_{_j}(y|X)}{Z_{_\theta}(X)}$$ where Yj is a sebsequence of Y, which is the sequence of labels, conditioned on the observed variable X, and Z(X) is a normalization function of X, and α and β are the forward and backward probabilities used in forward-backward algorithms. Is this correct?


title: Evaluation and Testsets for NNMF

text: I am trying to evaluate my recommender system which uses Non-negative Matrix Factorization. Some things that I evaluate are   * How does the size of the feature matrix affect the recommendations   * How does the amount of "rated movies" affect the recommendation   * ... Should I evaluate these things with a usual trainingset-testset separation? Like 60% Trainingset and 40% Testset? What is about the fact that every factorization of NNMF leads a bit different results, because of non-convexity and random initialization? I would use a random 60/40 separation, but I thought I´d better ask before I do it wrong. And if I would decide to make several runs of NNMF in order to select the best results, then I would need to do cross validation right? Or can I use the same training/testset for every run?


title: Reciprocal roots and eigenvalues relationship in time series

text: I came across a result in a time series textbook the other day and have not been able to understand why it is true (the authors don't give a proof but just state it as true). I want to show that the eigenvalues of the matrix $\mathbf{G}$ given by $$G= \begin{pmatrix} \phi_1&\phi_2 &\phi_3 &...&\phi_{p-1} & \phi_p\\ 1 & 0 &0 &...& 0 &0\\ 0 & 1 & 0 &... &0 &0\\ \vdots & & & \ddots&0&0\\ 0 & 0 &...&...&1 &0 \end{pmatrix} $$ correspond to the reciprocal roots of the $AR(p)$ characteristic polynomial $$\Phi(u)=1-\phi_1u-\phi_2u^2-...-\phi_pu^p$$ The one thing i was able to deduce is that the eigenvalues of $\mathbf{G}$ must satisfy $$\lambda^p-\phi_1\lambda^{p-1}-\phi_2\lambda^{p-2}-...-\phi_{p-1}-\phi_p=0$$


title: How to calculate relative error?

text: Let's say I have a sensor that measures pressure in a range from 10 - 60 mmHg. This sensor has an error of +/- 0.8 mmHg. Is there a way to quantify how large this error is with respect to the range of values the sensor can sense? I want to figure out what the relative size of this error is. Thanks!


title: What type of regression to use with negative values?

text: If possible, please explain things like I'm 5. I know very little about this subject, but would like to learn more. I have a data frame (in R) containing `player_id`, `points`, `away`, `opponent_fact_1`, `opponent_fact_2`. `points` can be negative. `away` lets us know if the game was at home or away (0 or 1). `opponent_fact_1` and `opponent_fact_2` gives us a stat about the opponent. `opponent_fact_1` is on a scale of about 0.0-5.5. `opponent_fact_2` is on a scale of about 70.0-95.0. The issue with the two facts is that there are fewer opponents as you reach the upper and lower bounds, so fewer data points exist at those levels. How can I determine how much of an influence `away`, `opponent_fact_1`, and `opponent_fact_2` has on a player's `points`? I asked someone online how to do this and he said use poisson regression, but didn't go into detail. Why would regression be helpful here? What is it? And I read that you can't use poisson regression with negative values? Also, how do I deal with the fewer data points around the upper and lower bounds? I'm using R, so any examples in R would be awesome. Explaining the output would be even better. I hope this isn't asking for too much. **EDIT: Added sample data**                 player_id opponent_team_id away  points opponent_fact_1 opponent_fact_2      1       695               22    0     0.0        2.888889           81.58      2       695               30    1     1.2        2.750000           81.58      3       695                4    1     3.0        3.714286           69.57      4       695               20    0    -3.0        3.000000           84.09      5       695               14    0     0.0        2.444444           72.97


title: Inconsistency in model diagnostics multinomial regression

text: I am using ordinal regression (parallel) on a total dataset of 15706 points, for 10 categories, and testing whether additional predictors significantly improve the model fit. I started with 2 predictors, obtained a deviance of 16530 and plotted expected vs observed counts: ![Counts - 2 predictors](http://i.stack.imgur.com/6R6Oh.png) Then I added a 3rd predictor (which I logically believe should help improve the fit), obtained a deviance of 15241, also compared this model with the first one (2 predictors) using the likelihood ratio test and found that the more complex model fits the data better. However, when I plotted the counts for this case, it looks worse than the fit using 2 predictors. ![Counts - 3 predictors](http://i.stack.imgur.com/Bc13v.png) What could be possible reasons for this? For information, I am using matlab and counting unique combinations of predictors for each response category, which means I will actually have more points in the case with 3 predictors (including more points with only one instance). I cannot see why this would be a problem but this might affect the results in a way I am not aware of.


title: Generating survival times for a piecewise constant hazard model with two change points

text: When there are two change points in a piecewise constant hazard model then the density function becomes some triangle exponential distribution. In this situation I can't generate the survival time from the CDF using probability integral transformation. Can any one help me to generate the survival time from this model?


title: Which test? Categorical data, not independent groups.

text: I would really appreciate it if people could give me advice on which statistical test that I should be using. I have a dataset which comprises of categorical data: presence or absence of different types of behaviour. Specifically, I am investigating, out of those who engage in a given behaviour, how many engage in different topographies of that behaviour. Participants have been assessed for multiple topographies and may engage in multiple topographies (therefore the groups aren't necessarily independent). For this latter reason, I am unsure whether a chi-squared test is appropriate? As an example: in total 619 individuals engaged in self-injury. Of these, 261 self-hit, 425 self-bit, 146 skin-picked and, finally, 174 self-scratched. All participants were assessed for all behaviours however we do not have data (and, unfortunately, no way to access data) for each individual- only overall numbers. Therefore, I cannot say for a given participant whether, for instance, they engaged in self-scratching picking but not biting or hitting. I am keen to see if there is a statistically significant difference between the prevalence of each of these topographies of behaviour. At present, I have coded the data in SPSS with topography in the first column (1=hitting, 2=biting etc) and presence/absence in second column (1= present, 2=absent). I did a chi-squared test preliminarily which came out as significant (X(3, N=2468)=317.678, p≤.001) however, as mentioned previously, I am unsure whether this is appropriate. Any advice would be greatly appreciated!


title: Expected value of minimum order statistic from a normal sample

text: **UPDATE Jan 10th 2014:** _the mistake was found - a math typo in one of the sources used. Preparing correction..._ **UPDATE Jan 25th 2014:** _the mistake is now corrected. Please ignore the calculated values of the Expected Value in the image uploaded - they are wrong- I don't delete the image because it has generated an answer to this question ._ The density of the minimum order statistic from a collection of $n$ i.i.d continuous random variables with cdf $F_X(x)$ and pdf $f_X(x)$ is $$f_{X_{(1)}}(x_{(1)}) = nf_X(x_{(1)})\left[1-F_X(x_{(1)})\right]^{n-1} \qquad [1]$$ If these random variables are standard normal, then $$f_{X_{(1)}}(x_{(1)}) = n\phi(x_{(1)})\left[1-\Phi(x_{(1)})\right]^{n-1} = n\phi(x_{(1)})\left[\Phi(-x_{(1)})\right]^{n-1}\qquad [2]$$ and so its expected value is $$E\left(X_{(1)}\right) = n\int_{-\infty}^{\infty}x_{(1)}\phi(x_{(1)})\left[\Phi(-x_{(1)})\right]^{n-1}dx_{(1)}\qquad [3]$$ where we have used the symmetric properties of the standard normal. In Owen 1980, p.402, eq.[ **n,011** ] we find that $$\int_{-\infty}^{\infty}z\phi(z)\left[\Phi(az)\right]^m dz= \frac {am}{(\sqrt {a^2+1})(\sqrt {2\pi})}\int_{-\infty}^{\infty}\phi(z)\left[\Phi\left(\frac {az}{\sqrt {a^2+1}}\right)\right]^{m-1}dz\qquad [4]$$ Matching parameters between eqs $[3]$ and $[4]$ ($a=-1$, $m=n-1$) we obtain $$E\left(X_{(1)}\right) = -\frac {n(n-1)}{2\sqrt{\pi}}\int_{-\infty}^{\infty}\phi(x_{(1)})\left[\Phi\left(\frac {-x_{(1)}}{\sqrt 2}\right)\right]^{n-2}dx_{(1)}\qquad [5]$$ Again in Owen 1980, p. 409, eq [ **n0,010.2** ] we find that $$\int_{-\infty}^{\infty}\left[\prod_{i=1}^{m}\Phi\left(\frac{h_i-d_iz}{\sqrt {1-d_i^2}}\right)\right] \phi(z)dz= \mathcal Z_m(h_1,...,h_m;\{\rho_{ij}\})\qquad [6]$$ where $\mathcal Z_m()$ is the standard multivariate normal, $\rho_{ij}=d_id_j,\; i\neq j$ is the pair-wise correlation coefficients and $-1\le d_i\le 1$. Matching $[5]$ and $[6]$ we have, $m=n-2$, $h_i=0,\;\forall i$, and $$\frac{d_i}{\sqrt {1-d_i^2}} = \frac 1{\sqrt 2} \Rightarrow d_i = \pm \frac 1{\sqrt 3} \forall i \Rightarrow \rho_{ij} = \rho = 1/3$$ Using these results, eq $[5]$ becomes $$E\left(X_{(1)}\right) = -\frac {n(n-1)}{2\sqrt{\pi}}\mathcal Z_{n-2}(0,...,0;\rho=1/3)\qquad [7]$$ **This multivaririate standard normal probability integral of equi-correlated variables, all evaluated at zero** , has seen enough investigation, and various ways to approximate and compute it have been derived. An extensive review (related to the computation of multivariate normal probability integrals in general) is Gupta (1963). Gupta provides explicit values for various correlation coefficients, and for up to 12 variables (so it covers a collection of 14 variables). The results are **(THE LAST COLUMN IS WRONG)** : ![enter image description here](http://i.stack.imgur.com/RCJKt.jpg) Now if we graph how the value of $\mathcal Z_{n-2}(0,...,0;\rho=1/3)$ changes with $n$, we will obtain ![enter image description here](http://i.stack.imgur.com/XDy7G.jpg) **So I arrive at my three questions/requests:**   1) Could somebody check analytically and/or verify by simulation that the results for the expected value are correct (i.e. check the validity of eq $[7]$)? 2) Assuming that the approach is correct, could somebody give the solution for normals with non-zero mean and non-unitary variance? With all the transformations I feel really dizzy. 3) The value of the probability integral seems to evolve smoothly. How about approximating it with some function of $n$?


title: How to prove that a t-distribution can be written as a ratio distribution?

text: If $X \sim N(0,1)$ and $Y \sim \chi^2(n),$ then it's "known" that $Z = X/\sqrt{Y/n}$ is $t$ distributed. Is there anywhere a proof for this? That in the end one can see the $t$ distribution?


title: Need normal distribution for kappa Fleiss and ICC (2, 1)?

text: I'm analysing a study of 4 raters and 34 subjects. I'm looking at the inter rater reliability. For the ordinal results, I did a weighted kappa of Fleiss and for the other continuous results I calculated an ICC (2,1). My question is if the data has to follow a normal distribution to be able to use the Kappa of Fleiss and the ICC (2, 1) correctly?


title: How to correctly apply the Nemenyi post-hoc test after the Friedman test

text: I'm comparing the performance of multiple algorithms on multiple data sets. Since those performance measurements are not guaranteed to be normally distributed, I chose the Friedman Test with the Nemenyi post-hoc test based on Demšar (2006). I then found another paper that, aside from suggesting other methods like the Quade test with subsequent Shaffer post-hoc test, they apply the Nemenyi test differently. # How do I apply the Nemenyi post-hoc test correctly? **1. Using the Studentized range statistic?** In Demšar's paper it says to reject the null hypothesis (no performance difference of two algorithms) if the average rank difference is greater than the critical distance CD with $$ CD = q_{\alpha}\sqrt{{k(k+1)}\over{6N}} $$ > "where critical values qα are based on the Studentized range statistic > divided by $\sqrt{2}.$" After some digging I've found that you those "critical values" can be looked up for certain alphas, for example in a table for $\alpha = 0.05$, for infinite degrees of freedom (at the bottom of each table). **2. or using the normal distribution?** Just when I thought I knew what to do, I found another paper that confused me again, because they were only using the normal distribution. Demšar is stating a similar thing at page 12: > The test statistics for comparing the i-th and j-th classiﬁer using these > methods is $$ z = {{(R_i − R_j)}\over{\sqrt{{k(k +1)}\over{6N}}}} $$ The z > value is used to ﬁnd the corresponding probability from the table of normal > distribution, which is then compared with an appropriate $\alpha$. The tests > differ in the way they adjust the value of $\alpha$ to compensate for > multiple comparisons. At this paragraph he was talking about comparing all algorithms to a control algorithm, but the remark "differ in the way they adjust ... to compensate for multiple comparisons" suggests that this should also hold for the Nemenyi test. So what seems logical to me is to calculate the _p_ -value based on the test statistic $z$, which is normally distributed, and correct that one by dividing through $k(k-1)/2$. However, that yields completely different rank differences at which to reject the null hypothesis. And now I'm stuck and don't know which method to apply. I'm strongly leaning towards the one using the _normal distribution_ , because it is simpler and more logical to me. I also don't need to look up values in tables and I'm not bound to certain significance values. Then again, I've never worked with the _studentized range statistic_ and I don't understand it.


title: Significance of independent variables in linear regression models

text: I am trying to make some sense out of the results of a linear regression model. I have a dependent variable X, and, say, 3 independent variables Y1 Y2 Y3. I set up 5 models : (m1) X ~ Y1 (m2) X ~ Y1 + Y2 (m3) X ~ Y1 + Y3 (m4) X ~ Y1 + Y2 + Y3 (m5) X ~ Y1 + Y2*Y3 In models 1 to 3, all independent variables are significant. In model 4, when I include both Y2 and Y3 in the model, they both turn un- significant. However, when I include an interaction between Y2 and Y3 (m5), all variables main effects and the interaction turn significant. I am wondering whether I am overfitting the model, or if there might be one logical way to interpret the changes of significances for Y2 and Y3 in these analyses (m4 vs m5). With regards [edit] Y2 is numeric, Y3 is factor (3 levels). The parameters' sign are the following : (m5) main Y2 : +, main Y3 : +, Y2*Y3[2] : -, Y2*Y3[3] : -,


title: How many stickers do I need to complete my FIFA Panini album?

text: I'm playing the FIFA Panini Online Sticker Album, which is an Internet adaption of the classic Panini albums that are usually published for the soccer world cup, European championship, and possibly other tournaments. The album has placeholders for 424 different stickers. The purpose of the game is to collect all 424. The stickers come in packs of 5, which can be obtained through codes found online (or, in case of the classic printed album, bought at your local newsstand). I make the following assumptions:   * All stickers are published in the same quantity.    * One pack of stickers does not contain duplicates.  How can I find out how many packs of stickers I need to acquire in order to be _reasonably sure_ (let's say 90%) that I have all 424 unique stickers?


title: How to find MLE when samples depend on the estimated parameter

text: Can you show me what I'm doing wrong here? This is the homework problem: Consider a random sample $Y_1, \ldots , Y_n$ from the pdf $f_Y(y;\theta) = 2y\theta^2$ where $0\le y \le \frac{1}{\theta}$. Find the maximum likelihood estimator for $\theta$. The likelihood function for this pdf is $\prod_{i=1}^{n} f_Y(y;\theta) = \prod_{i=1}^{n} 2y_{i}\theta^2 = 2^n \theta^{2n} \prod_{i=1}^{n} y_i$. So to me, it looks like I can maximize the likelihood function by choosing $y_{max}$ for the MLE. But then what do I do about the condition $0\le y \le \frac{1}{\theta}$? I tried looking up how to do this but the only example I can find is the one about the uniform distribution, but that one is much more straightforward because the condition on y is just $0 \le y \le \theta$ and not $\frac{1}{\theta}$. To maximize $\frac{1}{\theta}$ I would choose my MLE to be $y_{min}$ but that wouldn't maximize the likelihood function. Any help would be GREATLY appreciated.


title: Starting Steps for building logistic regression model?

text: I have a categorical response variable. It is binary and represents the win or loss of a deal. Some of the independent variables used to predict the response are also categorical (like Geo, Region, and others...). These categorical variables have more than 3 categories. The rest of the variable are counts (like #face2face activities, #of CXO/VP meet, Business development activities). Should I use logistic regression to predict the response variable? If yes, please specify the steps needed to come up with the best model for prediction. How should I check the quality of the model, so as to decide which one is the best?


title: Is it helpful to have monotonic features when using a random forest for classification?

text: I am training a random forest for binary classification. Here is a plot of one of my features, which is an integer giving the number of months since an event. The y-axis gives the proportion of cases ($y_i=1$) for each value of the feature. Error bars indicate the $95\%$ binomial confidence interval. ![enter image description here](http://i.stack.imgur.com/uJW4T.png) The value on the extreme left hand side occurs when there has been no such event, hence the number of months since the event is undefined. It has currently been encoded with the value of $-999997$ but this is entirely arbitrary. My question is whether it is advantageous when training the random forest to re-encode this so that it fits in with the (roughly) monotonic relationship between the probability that $y_i=1$ and the number of months since the event. I have a vague understanding of how the decision trees that make up the random forest are constructed, and intuitively it seems to me that leaving it as $-999997$ requires a more complex rule. How much of a problem is this?