title: Clustering using correlation of residuals

text: I have worked out a naive approach for clustering different variables, and I wonder what better approaches are there. I have data with many more dimensions (variables) than observations (samples), and I wish to cluster the variables by their correlation. However, the observations correspond to different groups, for example treatment and control. If I just calculate the correlation coefficient, all variables that go up in one group will have a positive correlation, and that is not what I intend; I want to cluster them by within- group effects and disregard between-group effects. To give you a specific example: variables are abundances of certain substances in the serum of patients. Serum samples are the observations, and each serum sample is associated with a number of abundances (for example, abundances of amino acids, cytokines etc.). There are, say a hundred observations (serum samples), each associated with the same 1000 variables (abundances). I want to see which abundances are correlated: for example, which cytokines are correlated with which amino acids. The approach that I am taking is to fit a simple linear model which explains the effect of various groups and correlate the residuals from the model. This heuristic works in the sense that it gives "reasonable" (in terms of biology) grouping of the variables, for example it clusters genes with similar functions. Consider this example in R.               groups <- factor( rep( c( "a", "b" ), each= 50 ) )     set.seed( 123 )     x1 <- rnorm( 100 )     x2 <- rnorm( 100 ) # not correlated with x1 yet     x3 <- x1 + rnorm( 100 ) # positively correlated with x1     # simulate group effect     x1 <- x1 + rep( c( 0, 2 ), each= 50 )     x2 <- x2 - rep( c( 0, 2 ), each= 50 )     x3 <- x3 - rep( c( 0, 2 ), each= 50 )      Now, x1 and x2 appear negatively correlated (although within each group there is no correlation), and x1 and x3 appear not correlated at all (although within each group there is a positive correlation): ![enter image description here](http://i.stack.imgur.com/q1lkX.png) Then, I calculate the residuals:               x1r <- lm( x1 ~ groups )$residuals         x2r <- lm( x2 ~ groups )$residuals     x3r <- lm( x3 ~ groups )$residuals      The situation now is very much different and much more what I intend: ![enter image description here](http://i.stack.imgur.com/Ig7FJ.png) My question is this: what would be a statistically more appropriate way of achieving the same? Note that neither of x1, x2 or x3 can be treated as independent variables or response variables -- they are merely correlated, without a given as to which one is determines which other. My ultimate goal is to create clusters of variables that are correlated after removing the effects of the groups.


title: How to calculate ${\rm Pr}(F_A|F_B)$ or ${\rm Pr}(F_A| F_B, F_C, F_D)$ in a four node network?

text: I have just started learning these things, and I would really appreciate any help on this. There are 4 nodes in a network (say $A,B,C,D$), and each node has certain failure probabilities, namely ${\rm Pr}(F_A)$, ${\rm Pr}(F_B)$, ${\rm Pr}(F_C)$, ${\rm Pr}(F_D)$. Then, I have a $4 \times 4$ correlation Matrix which contains the correlation coefficients for every pair of nodes. **How can I calculate ${\rm Pr}(F_A|F_B)$ or ${\rm Pr}(F_A| F_B, F_C, F_D)$?**


title: How do I add this P-value to normal probability plot in Minitab?

text: ![Probability plot of standardized residuals.](http://i.stack.imgur.com/gVwaj.png) I go to Stat from main menu then Regression -> Regression... Although I have selected almost everything I cannot plot the extra info (p-value, mean, N, etc ) above in circle. I use Minitab 16.2.0.


title: Datasets for exploring Cognitive Diagnosis Modelling

text: I recently became interested in _Cognitive Diagnosis Modelling_ , which classifies respondents to skill profiles based on their responses and a hypothesized relation between skills and items (Q- or Transfer-Matrix). I had a look at at https://pslcdatashop.web.cmu.edu/ but the datasets there are quite big. Does anyone have a suggestion, where I can find _records of testing data with reasonable size, where a Q-Matrix is already specified_?


title: Multilevel models including random slopes: how to calculate variance

text: In a linear mixed model, you take the covariance between data into account by adding a random intercept per cluster. For example, you measure the effect of a drug campaign over time on students, and add a random intercept per school and a random intercept per student within a school: $Y_{ijk}=\beta_{0}+\beta_{1}*Time+b_k+b_{jk}+\epsilon_{ijk}$ With Y the outcome variable "how much does he knows with drugs. We assume that the random intercepts are normally distributed with mean and variance: $b_k \sim N(0, \sigma_3^2)$, $b_{jk} \sim N(0, \sigma_2^2)$ and $\epsilon_{ijk} \sim N(0, \sigma_1^2)$ With these variances of random intercepts we can calculate the variance of the treatment effect on the different levels, for example on the school level. But what if I want to model that the campaign has not the same average increase over time in every school? I add a random slope per school: $Y_{ijk}=\beta_{0}+\beta_{1}*Time+b_{0k}+b_{1k}*Time+b_{jk}+\epsilon_{ijk}$ Is in this model the variation of $b_{0k}$ still meaningful? Is my intuition correct that the variation between different schools increases over time, as the random slopes drive the schools further apart? And as a second question, without random slope we also know that the covariance between students within the same school is given by: $Cov_{school}=\frac{sigma_3^2}{sigma_1^2+sigma_2^2+sigma_3^2}$ What is this covariance if a random intercept is included? Thank you in advance! UPDATE: included a plot of the data I am working with, level 1 are the 8 repeated measurements at each measurement occasion, level 2 are the difference measurement occasions, level 3 are the patients. ![enter image description here](http://i.stack.imgur.com/MKt8e.jpg)


title: Why is posterior density proportional to prior density times likelihood function?

text: According to Bayes' theorem, $P(y|\theta)P(\theta) = P(\theta|y)P(y)$. But according to my econometric text, it says that $P(\theta|y) \propto P(y|\theta)P(\theta)$. Why is it like this? I don't get why $P(y)$ is ignored.


title: Is it always necessary to correct the sample correlations for a meta-analysis?

text: Some authors eg Hunter et al (1982) propose that we must correct sample correlations for the measurement error.It seems measurement error need not be corrected for when we take a recourse to the meta-analysis of sample- correlations In a meta-analytic perspective, measurement error(here it is construct bias) has a little effect on the sampling error variance.


title: What does it mean to scale random variables?

text: We just started learning asymptotic theory, and to prove the lindberg-levy central limit theorem, weak law of large numbers etc, we 'scale and standardize' the RVs so it ends up having a standard normal distribution when sample size T tends to $\infty$. The math behind it is straightforward but can someone give a quick intuition or an example of what it means to 'scale' a RV, by $T^{1/2}$ for example?


title: Best possible regression models

text: I am looking for different types of regression models for such data in which the dependent variable is always decreasing with the independent variables. The two's, I know and have used so for are               Y = C + C_a A + C_b B      and               Y = C + C_a A + C_b B + C_a2 A^2 + C_b2 B^2 + C_ab A B      The second one give me the best fit so far. I want to try some other ones too but don't know. Please suggest some. Thanks


title: How to condition a point process to the actual sampling data?

text: **Figure A** demonstrates a point process (object=rectangle) with marks as the length and width of the rectangle. **Figure B** shows a realization of the point process regarding the information extracted from the input data.   The **question** is:   How to condition (using pre-conditioning or post-conditioning) the realization to the existing data points with marks? ![enter image description here](http://i.stack.imgur.com/8koYz.jpg) **More Explanation:**   Figure **E** demonstrates point process which has location, and width and height marks added.   Figure **A** is our point process. Figure **B** added marks as width and height variables. Figure **C** is our realization of the point process by means of information extracted from Figure **A**. Figure **D** is final marked point process, a realization. We seek a method to condition this realization to Figure **B**. ![enter image description here](http://i.stack.imgur.com/TEghN.jpg)


title: Linear regression with weekly spending data

text: I want to make a linear regression where the dependent variable is the amount spent on some good and the independent variables are some characteristics of the consumer (age, gender, wage, etc). Specifically, I want to know how much this characteristics influence the annual spending on this good. But the way that the survey was made is asking how much was spent on the good in the week of reference. So a lot of zero consumption was recorded (more than 90% of the records are zero). So, every regression that I make to explain the spending as a function of the consumer characteristics have no significance. What is expected if the good has some sporadic way of being consumed. Is there any work around for this problem?


title: How to do Simple Confirmatory Factory Analysis/SEM in R?

text: I have a data set consisting of 166 observations on 24 variables (`a1`, `a2`, `a3`, `a4`, `b1`, ..., `b4`, ..., `f4`). Responses where made on a 6-point scale from 1 to 6:               d <- read.table("http://pastebin.com/raw.php?i=m1ZJuKLH")     str(d)     ## 'data.frame':   166 obs. of  24 variables:     ##  $ a1: int  7 7 7 1 1 7 7 7 7 1 ...     ##  $ a2: int  7 7 4 7 5 5 1 7 7 4 ...     ##  $ a3: int  7 7 5 1 1 7 7 3 7 1 ...     ##  $ a4: int  7 7 6 7 1 7 1 7 7 1 ...     ##  $ b1: int  1 2 5 7 1 4 1 7 4 2 ...     ## [...]     ##  $ f4: int  6 6 4 1 7 7 7 7 7 1 ...      Variables with a common letter (i.e., all `a`s, all `b`s, ..., all `f`s) are of the same type. I would like to test whether the latent factors behind the variable types are correlated with each other (or whether you can set the correlation to 0). An additional complication is that the latent variables `c` to `f` share a further latent variable (termed `c.to.f`) and I only want to test the correlation between `a`, `b` and `c.to.f`. I would like to run a confirmatory factor analysis (which essentially is a structural equation model) in R testing this. There are at least two mature packages of doing so `sem` and `openMX`.   **I am interested in opinions/code on which package would be the best or perhaps easiest to specify such a model.** * * * Edit: I would like to accept an answer which includes code samples.


title: What are the steps for Simulated Annealing (SA)?

text: **Simulated Annealing** is a method of simulation in which each entity in the system is perturbed in some characteristics e.g., location due to increasing and decreasing the heat (the energy) of the entire system.   My question is about the application of simulated annealing approach to conditioning simulations. For instance what are the steps to condition a simulation (for example see conditioning point process) to the existing sampling data?


title: Number of observations used for ARIMA modeling

text: Our professor keeps on writing in his slides that when you test different ARIMA models on your time series, one always has to keep T fixed. I assume he is talking about fitting your model using the same number of observations. If you try a model with a higher lag, you loose one observation to fit the data on. Is this correct? Is this why he writes restrict your sample? We have 214 observations, should I for example to be on the safe side only use 205 observations? And how does Eviews know to always use the same number of obersations? Ref screenshot: ![enter image description here](http://i.stack.imgur.com/93Zjo.png) Thanks in advance.


title: Theoretical error bounds of classification and regression trees

text: So, some algorithms were motivated by theoretical work, such as in the case of boosting. Adaboost was introduced as an algorithm for solving the hypothesis boosting problem. The bounds on the training error and generalization error were formulated along many other extensions and explanations about how and why the algorithm works. My questions: 1.) Regarding classification and regression tress (CART), would it be correct to say that learning a decision tree was motivated by data and not by theory? 2.) Regarding classification and regression tress, are there theoretical bounds on the training and testing errors? If so a reference would be great.


title: How to add two random variable's pdf?

text: X,Y are independent random variables. X's pdf = f(x) Y's pdf = g(x) if Z= X+Y what is the Z's pdf? Can it be calculated?


title: Inner correlation of occurrences (burstiness?) in R

text: I want to measure the inner correlation of the occurrences of events i.e. I want to distinguish between the two (drawn) and say "in the second sample events occur more conglomerate compared to the first": ![enter image description here](http://i.stack.imgur.com/NxRBP.png) Isn't this different from burstiness (i.e. easier to compute as it doesn't invole traffic-values)? I have drawn a possible representation of the expected result (for S2) in the image below the samples (just as it is in my mind now), but don't hesitate to suggest different proposals. I looked out and tried different build-in functions and googled/ searched CrossValidated for R and burstiness or "inner correlation", but I didn't make progress. Maybe I am looking for the wrong search terms. If I wanted to really measure burstiness, it seemed to me (by reading papers about burstiness) that there would be no common known measure (but many different one). I we would pick one it would be cool to justify why choosing this specific one. This is how the time series is represented in R, currently: (written by dput in file)               c(3.256861, 3.377142, 3.941173, 4.304236, 4.485358, 4.606512,      4.707296, 5.473004, 5.714746, 5.815394, 5.835405, 5.936067, 5.957008,      6.964611, 7.045158, 7.065171, 7.165824, 7.669618, 8.17324, 8.273692,     9.503988, 9.604991, 9.624853, 9.725522, 10.237766, 10.954529,      11.378399, 12.687714, 13.291919, 13.41258, 13.67527, 14.380529,      14.743638, 15.247138, 15.851832, 15.952875, 15.972497, 16.456259,      16.476052, 17.201506, 17.463708, 18.068535, 18.309645, 18.390292,      18.410299, 18.430323, 18.531736, 18.652921, 18.793662, 19.297076,      19.639692, 19.760698, 20.768096, 20.868441, 20.990499, 21.494412,      21.856368, 22.199341, 22.219143, 22.440472, 22.481118, 23.327013,      23.447678, 23.811188, 23.843, 24.113302)


title: How to calculate a ROC curve? What kind of inputs do you need for that?

text: I have a microarray data which I ran through a continuous covariate (say "X"). I did this using 4 different methods. For the results obtained from each of the 4 methods, I have the following from each of the 4:   1. p-value,    2. FDR significant value,    3. "Bonferonni significance" (saying "true" or "false")   4. "Holm significance" (saying "true" or "false") I wish to obtain an ROC curve showing lines for sensitivity and specificity from each of the 4 methods. I kind of know what is an ROC curve, but even after reading through a few links, I dont feel clear of the concept of what kind of input it needs to create it. So I apologise if this seems to be a dummy question. So my questions are:   1. Can a ROC curve be created from the 4 values (p value, FDR, Bonferoni & Holm) which I have mentioned above? Or is it that I have to calculate the FP and FN first?   2. Can you suggest me an easy to use package which could do this for me, by giving me a ROC curve with lines of 4 methods using the input above? Your help appreciated. Thanking you in advance.


title: How to create artificial data with one binary response variable?

text: I want to check various classification model like random forest, tree, knn,etc. I used some bench marking data set but now I need to simulate my own data set with a binary response variable.


title: ANOVA post hoc with concerned about the sphericity

text: According to this post the best way to do post hoc analysis in ANOVA in R is t-test + ptukey. Because lme + glht "essentially assumes the sphericity assumption holds". How are t-test + ptukey "concerned about the sphericity"? What are the disadvantages of this method?


title: How to compute the PDF of a sum of a discrete and a continuous random variable?

text: I have a problem with this exercise in probability and statistics: Calculate the probability density function (PDF) of $$Z=X+Y$$ where $Y$ is discrete random variable which is equal to $-1$ and $1$ with equal probability; $X$ is a standard gaussian random variable independent of $Y$. I know that the PDF of the sum of two continuous independent variables is given by the convolution of the marginal PDFs $$f_z(z)=f_x*f_y,$$ but if one of the two variables is discrete what should I do?


title: Understanding the Polya urn model

text: I am having trouble understanding the plot below (taken from Edwin Chen's blog).   * What is the **x-axis** supposed to represent? Shouldn't color be a categorical variable? Does the **x-axis** have to be on the real line for the Polya urn model?   * Are we supposed to see significant changes across runs for the same alpha? ![enter image description here](http://i.stack.imgur.com/7G5hD.png)


title: Factorial experiment with non-random sample – is it possible?

text: Factorial experiments have several advantages when one wants to examine the effect of several independent variables (factors) on an outcome. Perhaps most importantly, they often allow additional factors without major increases in sample size. However, sampling is often the most difficult thing to begin with, and in some cases, even a 2x2 factorial design can require a sample that is practically impossible to obtain. It is not necessarily the sample size that is the problem, but the feasibility of obtaining a _random_ sample. I'm considering a case where, for practical and ethical reasons, I have no options but to compare a preselected group of people to another group which I can select randomly from a larger population. The purpose is to examine whether a certain treatment has an effect on a behaviour which everyone in the population (including the preselected group) will carry out continuously in collaboration with each other. Everyone in the population has gone through an integration process in which they have learned, one way or the other, the behaviour. The question is whether the treatment is effective for quickly learning the desired behaviour in the preselected group – is the behaviour learned in a shorter time and with more accuracy than what previous newcomers have been able to learn without the treatment. There will also be other factors in the design except the treatment. In an ideal case, I could randomly select two groups that are part of the population, and give the treatment to one of them. However, it is impossible in this case, because I cannot give the treatment to anyone but the preselected group. In addition, due to ethical reasons, I cannot randomly split the preselected group in half and give the treatment only to one half: everyone in the group must get the treatment. Thus I am forced to give the treatment to everyone in the preselected group and use a random sample of the existing population as a proxy for the non-treatment group. Are there any valid ways of conducting a factorial experiment under these conditions where I have to deal with a non-random sample? If not, are there any other ways to assess the effect of the treatment under these conditions?


title: Is there a reason to score on subsets of your data?

text: If you're training a model using cross validation on one set of data, and then you're scoring that model on a separate set of data, is it advised to score, say, 1000 times on subsets of your data? I'm thinking taking bootstrap samples or something, and then scoring the model, then another boot strap, etc. It seems like this way you'd get an idea of the stability of model performance, but I'm not sure if that's legit or not...


title: How to find the joint density of 3 normal variables?

text: Suppose $U, W, V, S$ are four independent normal random variables with mean $0$ and variance $1$. Let $X=W+U$, $Y=2W+S$, $Z=3W+V$. What is $f(X, Y, Z)$? Thanks!


title: Need help understanding response from Metis

text: I was wondering if any of you could help me understand the response I got from this clustering algorithm (Metis). As you probably can see, I'm trying to cluster IP addresses based on common records accessed. The graph file contains one node per IP, and weighted edges, where the weight is equal to the number of common accesses. A couple of specific questions:   * This algorithm requires a pre-defined number of clusters. How should I, based on this response, decide on the optimal number? (The algorithm is really fast, so repeated runs is possible).   * How should I go about evaluating this clustering? I don't have a training/test set as of now. Should I create one? How would I do that with a binary access log? If you have any other tips regarding this algorithm, they are more than welcome! Thanks for taking a look! Here comes the response:               ******************************************************************************     METIS 5.0 Copyright 1998-13, Regents of the University of Minnesota     (HEAD: , Built on: Jul  1 2014, 14:23:44)     size of idx_t: 32bits, real_t: 32bits, idx_t *: 64bits          Graph Information -----------------------------------------------------------      Name: gt2rec_gt2ip_b_rec_sorted_nobots_pdf_full.graph, #Vertices: 103924, #Edges: 34678560, #Parts: 10          Options ---------------------------------------------------------------------     ptype=kway, objtype=cut, ctype=shem, rtype=greedy, iptype=metisrb     dbglvl=0, ufactor=1.030, no2hop=NO, minconn=NO, contig=NO, nooutput=NO     seed=-1, niter=10, ncuts=1          Direct k-way Partitioning ---------------------------------------------------      - Edgecut: 16083558, communication volume: 662037.           - Balance: constraint #0:  1.030 out of 0.000           - Most overweight partition:      pid: 4, actual: 10704, desired: 10392, ratio: 1.03.           - Subdomain connectivity: max: 9, min: 9, avg: 9.00           - The original graph had 1243 connected components and the resulting         partitioning after removing the cut edges has 2362 components.     Timing Information ----------------------------------------------------------      I/O:                  7.162 sec      Partitioning:        10.511 sec   (METIS time)      Reporting:            2.118 sec          Memory Information ----------------------------------------------------------       Max memory used:       1021.427 MB     ******************************************************************************


title: When to transform predictors?

text: I am analyzing data from an experiment in which treatment levels increase quadratically, e.g. the treatment levels are $0, 1, 4, 9$. When analyzing the response using regression, would it make sense to use the square root of the treatment level as a predictor? If so, how would this affect interpretation?


title: What is the best site to hire contract or permanent stats folks

text: I really enjoy this resource -- its been super helpful to learn and read along as I explore some concepts for a new venture. But now I need real help -- in the short term contract, and in the long term statistical pros. What are the best sites were folks look for side projects and gigs? What is the best advice you have for someone hiring stats folks?


title: For classification with Random Forests in R, how should one adjust for imbalanced class sizes?

text: I am exploring different classification methods for a project I am working on, and am interested in trying Random Forests. I am trying to educate myself as I go along, and would appreciate any help provided by the CV community. I have split my data into training/test sets. From experimentation with random forests in R (using the randomForest package), I have been having trouble with a high misclassification rate for my smaller class. I have read this paper concerning the performance of random forests on imbalanced data, and the authors presented two methods with dealing with class imbalance when using random forests. **1. Weighted Random Forests** **2. Balanced Random Forests** The R package does not allow weighting of the classes (from the R help forums, I have read the classwt parameter is not performing properly and is scheduled as a future bug fix), so I am left with option 2. I am able to specify the number of objects sampled from each class for each iteration of the random forest. I feel uneasy about setting equal sample sizes for random forests, as I feel like I would be losing too much information about the larger class leading to poor performance with future data. The misclassification rates when downsampling the larger class has shown to improve, but I was wondering if there were other ways to deal with imbalanced class sizes in random forests?


title: Standard deviation of means over two dimensions

text: (Apologies if the notations are "unusual", I'm not sure what the correct notations should be. I'm putting an example at the end of the question.) Let's assume there was an initial dataset of an n by m matrix $M=(x_{ij})$ with $1<=i<=n$ and $1<=j<=m$, from which the following two vectors have been calculated:   * the vector of the column-wise means, that is, a vector of length $n$ where each element of index $j$ is the mean of the $m$ elements in column $j$ of the matrix: $(\bar{x}_{*j})$ with $1<=i<=n$, where $\bar{x}_{*j} = \frac{1}{m}\sum_{i=1}^m{x_{ij}}$.   * the corresponding vector of the columns-wise standard deviations ($s_{*j}$). I'd like to get the mean and standard deviation of the row-wise vector of means, that is, the column vector made of the mean of each row: $(\bar{x}_{i*})$ with $1<=j<=m$, where $\bar{x}_{i*} = \frac{1}{n}\sum_{j=1}^n{x_{ij}}$, assuming that the initial matrix has been lost. The mean of the elements in the row-wise vector is also the mean of the elements in the column-wise vector (which is also the mean of all the elements in the matrix, with equal weight): $\frac{1}{m}\sum_{i=1}^m{\bar{x}_{i*}} = \frac{1}{m}\sum_{i=1}^m({\frac{1}{n}\sum_{j=1}^n{x_{ij}}}) = \frac{1}{n\times m}\sum_{i=1}^m\sum_{j=1}^n{x_{ij}}$ Is there any way to get the standard deviation within the row-wise vector of means without having the original matrix? For example: $$ M = \left(\begin{matrix} 1 & 2 & 3\\ 7 & 5 & 4\\ 8 & 2 & 3\\ 5 & 2 & 4 \end{matrix}\right) $$ $$ (\bar{x}_{*j}) = \left(\begin{matrix}5.25 & 2.75 & 3.5\end{matrix}\right) $$ $$ (s_{*j}) = \left(\begin{matrix} 2.68 & 1.29 & 0.5\end{matrix}\right) $$ $$ (\bar{x}_{i*}) = \left(\begin{matrix} 2\\ 5.33 \\ 4.33 \\ 3.66\end{matrix}\right) $$ The mean of $\{2, 5.33, 4.33, 3.66\}$ is $3.83$ (also the mean of ${5.25, 2.75, 3.5}$). The standard deviation of $\{2, 5.33, 4.33, 3.66\}$ is $1.21$. Would there be any way of calculating this standard deviation, knowing $(\bar{x}_{*j})$, $(s_{*j})$, but without knowing $M$ nor $(\bar{x}_{i*})$? Is there anything else that could have been calculated "column-wise" (that is, independently for each column, thus excluding the "row-wise" means themselves) that would help find out the (preferably sample) standard deviation of the "row-wise" means? Thank you.


title: Why are random variables defined as functions?

text: I'm having problems understanding the concept of a random variable as a function. I understand the mechanics (I think) but I do not understand the motivation... Say $(\Omega, B, P) $ is a probability triple, where $\Omega = [0,1]$, $B$ is the borel-$\sigma$-algebra on that interval and $P$ is the regular lebesgue measure. Let $X$ be a random variable from $B$ to $\{1,2,3,4,5,6\}$ such that $X([0,1/6)) = 1$, $X([1/6,2/6)) = 2$, ..., $X([5/6,1]) = 6$, so $X$ has a discrete uniform distribution on the values 1 through 6. That's all good, but I do not understand the necessity of the original probability triple... we could have directly constructed something equivalent as $(\{1,2,3,4,5,6\}, S, P_x)$ where $S$ is all the appropriate $\sigma$-algebra of the space, and $P_x$ is a measure that assigns to each subset the measure (# of elements)/6. Also, the choice of $\Omega=[0,1]$ was arbitrary-- it could've been $[0,2]$, or any other set. So my question is, why bother constructing an arbitrary $\Omega$ with a $\sigma$-algebra and a measure, and define a random variable as a map from the $\sigma$-algebra to the real line?


title: Methods to best test lead/lag relationships

text: I was wondering if you can share your experiences on what you feel is the best method to test lead / lag relationships between I(1) time series variables (i.e stock prices) and advantages and disadvantages of your proposed method(s). Also if you have links to academic papers that further describe these methods I would greatly appreciate them. I have read several papers that speak of VECMs, IRSUR, simple OLS, threshold regression etc.. but I'm not sure which to use for my study. I am trying to establish lead/lag relationships in intraday stock data returns (using 1 minute price time series) between stock prices of ~50 companies in one particular sector. Appreciate the help.


title: Chow Structural Test in R

text: I am sitting on a pile of data concerning wages at a local company and other information, such as the gender, whether the person in question belongs to a minority group etc. What I would like to investigate is whether an additional year of education gives the same relative increase in wages for lower and higher education. For this purpose, I have divided the original data into two subcategories; one group with at most 12 years of education and the other with 13 or more years of education. What I would like to perform is a Chow test of structural change of the form \begin{equation} F=\frac{\frac{(RSS-(RSS_{\text{lower}}+RSS_{\text{higher}}))}{2}}{\frac{RSS_{\text{lower}}+RSS_{\text{higher}}}{n-2k}}=\frac{(n-2k)((RSS-(RSS_{\text{lower}}+RSS_{\text{higher}}))}{2(RSS_{\text{lower}}+RSS_{\text{higher}})} \sim F(k,n-2k) \end{equation} where $n$ is the total number of observations and $k$ is the number of explanatory variables. Clearly, I could simply calculate the $RSS:$s directly and then construct $F$ explicitly.                total<-lm(SALARY~EDUC+GENDER+MINORITY)         lower<-lm(SALARY1~EDUC1+GENDER1+MINORITY1)         higher<-lm(SALARY2~EDUC2+GENDER2+MINORITY2)         RSStot<-sum(residuals(total)^2)         RSSlow<-sum(residuals(lower)^2)         RSShig<-sum(residuals(higher)^2)         ((576-6)((RSStot-(RSSlow+RSShig)))/(2(RSSlow+RSShig))      > Nevertheless, I am sure there must be a way to do this directly in R. How > exactly would that be? I would truly appreciate any enlightenment from any > kindhearted spirit. Cheers to all!


title: Statistically back-calculating: Markov Chain?

text: I would like to calculate the value of bacteria on 4 surfaces $i=\{1..4\}$. A person touches some of those 4 surfaces at random and a count is made on their finger after each surface contact ($x_i$). Someone lost the bacteria count ($x_i$) after each surface but I do know the total count (X) on a persons's finger after they've touched a number of surfaces. I also know which ones and in which order. What I know:   1. Final bacteria count on a person's finger: $X$   2. Transfer efficiency from surface to finger: $PT_i=\displaystyle \frac{\text{Finger contact area}}{\text{Area of surface}_i}\frac{1}{\gamma_i}$ where $\gamma$ is a surface dependent constant.   3. The number of times the person touched a particular surface: $h_i$. If I had surface counts $C_i$, the summation of bacteria is linear: ie $\begin{eqnarray} h_1C_1PT_1&=&x_1\\ h_2C_2PT_2&=&x_2\\ \vdots\quad &=& \quad \vdots\\ h_iC_iPT_i&=&x_i \end{eqnarray}$ such that summing over all surfaces $i$ the total count x is: $\displaystyle \sum_i h_iC_iPT_i=\sum_i x_i=X$. Can I back calculate $C_i$, without $x_i$ even statistically or probabilistically?   Best regards.


title: How to show that answers to a set of 16 questions are different based on gender.

text: I have set of 16 questions that have been answered on a numerical scale, from 1 to 7. They are answered by 200 couples, and I want to show there is a difference in the answers on gender. On a casual glance at the data, I see there is little difference between the genders for 6 of the questions, but is for the other. If I run a series of t-tests for each, and show significant differences for 10 of them; I should be able to conclude that on the set as a whole, there are differences. Is there a better way to show that for this series, answers differ based on gender? Mu


title: Power calculations/sample size for biomarker study

text: We have a potential biomarker for predicting whether a patient has cancer or not. The biomarker test result is binary being positive or negative. We want to get some sense of the amount of patients that need to be tested to determine whether this biomarker is a good predictor or not. From reading on the internet it seems that the way to go is to look at the sensitivity (for the number of cases) and specificity (for the number of controls). It is suggested that you should treat this situation as a one- sample proportion test, but it remains unclear how you should go about estimating what the sensitivity is and the range you are prepared to except. If say I consider any biomarker with a sensitivity of greater than 0.8 to be "good", how would you set the two variables up? I would like my null hypothesis to be the biomarker is no better than a random assignment i.e. a sensitivity of 0.5. Could anyone give an example of the best way to do this (especially if it is in R).


title: Checking assumptions LMM: residual plot with diamond shape

text: I am running a linear mixed model and want to check its assumptions. The model I run is comparing if males and females behave differently over time (timeclass=1,2,3,4):                x <- lme(response ~ timeclass*sex, random = ~ 1|subject,  method="ML", data=dat)      The code to create the plots:               plot(x)     plot(x, response ~ fitted(.) | sex, abline = c(0,1))     qqnorm(x,~resid(.)|timeclass)     qqnorm(x,~resid(.)|sex)     hist((resid(x) - mean(resid(x), na.rm=T)) / sd(resid(x), na.rm=T), freq=F); curve(dnorm, add = TRUE)      This is the output I get: ![residual plot](http://i.stack.imgur.com/4gYnC.png) ![linear relationship males and females](http://i.stack.imgur.com/Pykhd.png) ![qqplot for sex](http://i.stack.imgur.com/nm6aJ.png) ![qqplot for timeclass](http://i.stack.imgur.com/IyXNj.png) ![histogram](http://i.stack.imgur.com/66naf.png) The plots except the residual plot seem to show that there seems to be a linear relationship between the fitted values per sex and the response and that the errors are relatively normally distributed in both sex and across the four time classes. However, the first plot has a very clear diamond shape, almost so clear that it seems values are cut off in the four triangular angles. Can someone help me understand what is going on and what this means? Is the model I run valid based on the plots for assumptions?


title: Can one obtain a p-value for all levels of factor when doing a GLM?

text: Please help! I have recently been criticized for using pairwise comparisons to explain all three levels of a factor within a negative binomial GLM rather than all levels at once. I was told that it is "long-winded" and "uneccessary". I was under the impression that in GLMs one cannot bulk all levels of a factor together to obtain a test statistic and corresponding p-value. Obviously if a factor is "insignificant" at any level then carrying out a post-hoc analysis is pointless. My levels all have there own p-values therefore I discussed these values from the below global model. I was told to do an ANOVA instead which I don't believe is suitable for overdispersed, zero- inflated data. p-value for all levels of a factor anyone? (Below, lower field layer 0, upper field layer 1 and change1 is in intercept)                   Deviance Residuals:          Min       1Q   Median       3Q      Max          -2.4284  -0.7956  -0.3862   0.4045   2.4233                Coefficients:                                 Estimate Std. Error z value Pr(>|z|)            (Intercept)                    4.3410884  1.8219786   2.383  0.01719 *          Height                         0.0373584  0.0119929   3.115  0.00184 **         Width                         -0.0007891  0.0008246  -0.957  0.33859            MeanMin                       -0.1731877  0.1404434  -1.233  0.21752            as.factor(Site_Treat)2        -0.4080256  0.2480438  -1.645  0.09998 .          as.factor(Change)2            -0.4940398  0.1755487  -2.814  0.00489 **         as.factor(Change)3            -0.1613766  0.1763677  -0.915  0.36019            as.factor(Lower_Field_Layer)1  0.4873488  0.2931585   1.662  0.09643 .          as.factor(Lower_Field_Layer)2 -0.3292409  0.3717863  -0.886  0.37585            as.factor(Upper_Field_Layer)2 -0.0081040  0.3257734  -0.025  0.98015            ---         Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1              (Dispersion parameter for Negative Binomial(4.7795) family taken to be 1)              Null deviance: 96.392  on 46  degrees of freedom         Residual deviance: 47.968  on 37  degrees of freedom         AIC: 403.94      Best wishes, Platypezid


title: Non-overlapping aggregation and autocorrelation

text: I have a time series that follow an MA(1) process , I do non-overlapping aggregation with aggregation level $m=2,3,\dots$, and then I calculate autocorrelation of basic and aggregated series, for negative value of moving average parameter, aggregated autocorrelation becomes close to zero, means that the aggregated series want to be like white noise, I was wondering if anybody can explain the reasons behind this behavior? ![enter image description here](http://i.stack.imgur.com/D0Xun.png)


title: Bounding mutual information given bounds on pointwise mutual information

text: Suppose I have two sets $X$ and $Y$ and a joint probability distribution over these sets $p(x,y)$. Let $p(x)$ and $p(y)$ denote the marginal distributions over $X$ and $Y$ respectively. The mutual information between $X$ and $Y$ is defined to be: $$I(X; Y) = \sum_{x,y}p(x,y)\cdot\log\left(\frac{p(x,y)}{p(x)p(y)}\right)$$ i.e. it is the average value of the pointwise mutual information pmi$(x,y) \equiv \log\left(\frac{p(x,y)}{p(x)p(y)}\right)$. Suppose I know upper and lower bounds on pmi$(x,y)$: i.e. I know that for all $x,y$ the following holds: $$-k \leq \log\left(\frac{p(x,y)}{p(x)p(y)}\right) \leq k$$ What upper bound does this imply on $I(X; Y)$. Of course it implies $I(X; Y) \leq k$, but I would like a tighter bound if possible. This seems plausible to me because p defines a probability distribution, and pmi$(x,y)$ cannot take its maximum value (or even be non-negative) for every value of $x$ and $y$.


title: Interpreting confidence intervals

text: Suppose I have a $95 \%$ confidence interval for $\mu$ with population variance $\sigma^{2}$ known. It will be of the form $$[\bar{X}-1.96 \sigma/\sqrt{n}, \bar{X}+1.96 \sigma/\sqrt{n}]$$ The statistical interpretation of this is that $95\%$ of the confidence intervals will contain the true mean? I know individually, there is a $50\%$ chance for a confidence interval to contain the mean. But collectively, my above sentence is valid?


title: Is the mean squared error used to assess relative superiority of one estimator over another?

text: Suppose we have two estimators $\alpha_1$ and $\alpha_2$ for some parameter $x$. To determine which estimator is "better" do we look at the MSE (mean squared error)? In other words we look at $$MSE = \beta^2+ \sigma^2$$ where $\beta$ is the bias of the estimator and $\sigma^2$ is the variance of the estimator? Whichever has a greater MSE is a worse estimator?


title: Limited Information Maximum Likelihood (LIML) estimation in R?

text: Curious whether anyone knows a package, or has written an implementation themselves, for conducting instrumental variables regressions using LIML in R. All of the R packages I have seen for IV regressions seem to resort to 2SLS/GMM as opposed to LIML, which may have more desirable finite sample qualities compared to 2SLS (Hayashi 2000) For additional context, stata's ivregress command includes options to use LIML estimation, and hoping someone has already implemented something similar in R so I don't have to write it myself.


title: Applying Recency In Logistic Regression

text: Are there any statistical concepts or theories on how to effectively measure recency to where recent events are given more weight than older ones. I'm creating a logistic regression model and would like to apply an adjustment to various factors based on the events recency. ...or, is it purely up to me to come up with an arbitrary formula? For example: One of my projects is predicting via logistic regression the performance of pro golfers in upcoming tournaments. Their recent form (how they played last week) is generally more important than how they played six months ago. Are there any specific techniques/approaches that tap into this concept?


title: Do we need a Bonferonni correction when running numerous univariate regression analyses?

text: For my data, I need to run numerous linear regression analyses, separately for each of my two groups of subjects as well as different conditions. Together, I would have to run approximately 40 regression analyses across the 3 conditions for each group. I'm aware that we should correct for multiple comparisons for the regression coefficients within a model. But in my case, should I also do a similar correction for the p-value associated with the regression model itself (ie: the F test)? For example, if in one condition, I have to run 25 regression analyses for each group, would it make sense to divide alpha=.05 by 25? I've been going through different websites and textbooks, but can't seem to find an answer. Any suggestions/help would be appreciated!


title: Alternatives to glmnet for feature selection on data with lots of NAs

text: I have a surgical database in which there are approximately 100,000 observations and 200 features. Each observation corresponds to a separate patient while the features correspond to either preoperative, perioperative, or postoperative variables (e.g. preoperative labs, length of operation, and days until death). As such some of the features are factors, while others are continuous data types. The database is not "clean" in that many of the features contain NAs because that data was not collected. Each patient has some feature with NAs in it, so when I try to select only observations without NAs that results in 0 observations. That being all said my end goal is to see whether any of these preoperative variables might be predictive of mortality. Initially I wanted to use L1 regularization (lasso) via glmnet to perform feature selection, however because of the NAs I can't run it. Are there any alternatives or techniques to bypass this problem? I assume that it's more common than not to have a database with NAs, so I wanted to see if you guys could fill me in. Thanks!


title: What could cause different conversion rates for the same content in separate split tests?

text: We have an events based tracking system for our website, with split testing built-in and we are using ABBA for the calculations. The problem comes up when we are doing consecutive split tests. For example, we test X and Y, Y wins, so the next test is Y and Z. We are seeing, most of the time, that the conversion rate of Y differs in the two tests... Sometimes by as much as 50%. We have a tool to track conversion rate over time, and I've included two graphs of conversion rate between two different tests of new and baseline. You can see that between the two tests, the previous new and the subsequent baseline differ by a significant amount even though they should theoretically be the same. Both experiments have about 200,000 runs for new and baseline - 400k total - 800k over two experiments, and the P values on both are .0001. I'm just unable to explain why the conversion rate of "new" in test1 was 1.28% but the conversion rate of "baseline" in test2 was .85% - when they are exactly the same content. I can certainly think of some things that would contribute to something like this: only catching part of a user cycle/time of day variations/seasonal variations/holidays, but my gut tells me that it shouldn't really make a 50% difference. We have, however, eliminated the possibility of Googlebot and other bots affecting the results, as we don't even record events for bots. Unfortunately, we don't have anything we can compare with. We do use Google analytics, but we only record our key metrics there, so there's no apples – apples. The experiments were run immediately after each other, and I'm not aware of any special marketing campaigns (We don't do very much marketing, and no CPC). We are seeing some pretty good changes in traffic from certain countries, and I can definitely look into that to see if it could contribute... But, again, 50%? Test 1: "New" wins ![Test 1 conversion](http://i.stack.imgur.com/1WSDF.png) Test2: The old "New" is now "Baseline"...note the mismatch in conversion ![Test 2 conversion](http://i.stack.imgur.com/l1LFv.png)


title: How important it is to have an intuition of why various statistical tests work?

text: I sometimes have difficulty understanding the intuition behind some tests. Is it important to know why some test works in a particular situation, or is it enough to learn things at the level "OK, now I can check homoscedasticity" or "I have no idea why we divide by 12 in Mann–Whitney U-test, but trust me, it works"? Or are there books where the authors have explained the ideas behind every test and provided some reasons to "prove" that the tests work?


title: How to find principal components without matrix algebra?

text: Principal component analysis _can use_ matrix decomposition, but that is just a tool to get there. How would you find the principal components without the use of matrix algebra? What is your objective function (goal), and what are your constraints? Thanks!


title: Correlation between two time series

text: What is the easiest way / method to compute the correlation between two time series that are exactly the same size? I thought of multiplying $(x[t]-\mu_x)$ and $(y[t] - \mu_y)$, and adding up the multiplication. So if this single number was positive, can we say these two series are correlated? I can think of some examples however where a linearly another exponentially growing time series would have no relation to eachother, but the computation above would report they were correlated. Any thoughts?


title: Goodness-of-fit for linear, mixed effects models

text: Can someone tell me why this is wrong: I want to estimate the goodness-of-fit for a linear mixed effects model. I regress the fitted values against the raw values and calculate an $R^2$. I suppose an accepted alternative would be to run the model with just the random effect and no fixed effects and compare AIC scores, but I'm specifically interested in the scenario above.


title: Uncertainties from likelihood surface

text: I have just been reading a paper which uses the following method to calculate parameter uncertainties from a likelihood surface of two parameters. Step 1) Calculate the likelihood for finding a given pair of parameters given the data at each point on the grid. Step 2) Normalise the surface to 1. The Likelihood value now give the probability to the parameter given the as L(a|x) is proportional to p(x|a). I'm not totally convinced by this - not least as the extent of the grid is arbitrary. I'm probably going to ticck the 'forward it on to a statician' box as it's a review and I'm not familiar with this method but I was curious to whether this works under any circumstances. I always assumed Wilkes' theorem was about as well as you could do this.


title: Two (or multiple) t-test replaced with a single pairwise.t.test()?

text: This is an artificial example. I want to compare the Income of the _same_ group of people (i.e. "pairing") while they were working in US, Spain, and Nigeria. In other words, I use the "US" income as a reference and simply conduct 2 Two-sample Paired T-tests:               t.test(dat[Country=="US",]$Income, dat[Country=="Spain",]$Income, paired=TRUE)     t.test(dat[Country=="US",]$Income, dat[Country=="Nigeria",]$Income, paired=TRUE)      But how is this different from the following?               pairwise.t.test(dat$Income, dat$Country, paired=TRUE)      The help for pairwise.t.test() states, "...with corrections for multiple testing". Is that the only difference?


title: Simulating a Gaussian process with an exponentially decaying covariance function

text: I'm trying to generate many draws (i.e., realizations) of a Gaussian process $e_i(t)$, $1\leq t \leq T$ with mean 0 and covariance function $\gamma(s,t)=\exp(-|t-s|)$. Is there an efficient way to do this that wouldn't involve computing the square root of a $T \times T$ covariance matrix? Alternatively can anyone recommend an `R` package to do this?


title: Regression through the origin

text: We have the following points: $$ (0,0)(1,51.8)(1.9,101.3)(2.8,148.4)(3.7,201.5)(4.7,251.1)(5.6,302.3)(6.6,350.9)(7.5,397.1)(8.5,452.5)(9.3,496.3)$$ How can we find the best fitting line $y=ax$ through the points? My calculator has the option to find the best fitting line $y=ax+b$ through these points, which is: $$y = 53.28x + 0.37$$ How can I find the best fitting $y=ax$? It seems to me we can't just remove the $0.37$ without compensating in the $a$?


title: Likelihood ratio tests on linear mixed effect models

text: I am currently running some analyses on a linguistic data set with a mixed effect model. The problem is, I think that one random factor should be excluded while my colleague thinks it should be included. The two options are:               lmer(intdiff ~ stress * vowel_group + (1|speaker) + (1|word), data)     lmer(intdiff ~ stress * vowel_group + (1|speaker), data)      How do we check which model best fits our data set? It was suggested that I use a likelihood ratio test, but as far as I can tell, there isn't a function in R that can be applied to 2 linear mixed effects models. Is there another way to tell which model is more predictive? Thanks


title: How to find the MLE for two Dirichlet distributions with same mean

text: I have samples from two groups. Samples from each group follow a Dirichlet distribution. The two Dirichlet distributions have the same mean but different precision. How to find the MLE of mean and two precisions?


title: Trade-off between using explicit inputs and using fewer inputs for neural networks

text: Suppose you have 2 distinct features x1 and x2 with discrete values between 1 and 10. When building a neural network, I have 2 options:   1. Create a new feature x3 = 10*x1 + x2 and use this as the only input.   2. Use both inputs x1 and x2 for the neural network. How can I decide which is the best option? Here is a practical example inspired from poker: You have 7 cards. The first 2 cards are from player 1. The next 2 cards are from player 2 and the last 3 cards are the flop. You want to determine who has the strongest hand. The two options are:   1. Treat each card as a number from 1 to 52.    2. Make two inputs: One for the rank (2,3,4,5,6...J,Q,K,A) and one for the suit (club,diamond,heart,spade). More precisely, I want to understand how choosing each option will affect how much training data is needed, how many neurons in the hidden layer I should use and how much time will be needed for training.


title: Covariance matrix for Gaussian Process and Wishart distribution

text: I'm reading through this paper on Generalised Wishart Processes (GWP). The paper calculates the covariances between different random variables (following Gaussian Process) using squared exponential covariance function, i.e., $K(x,x') = \exp\left(-\frac{|(x-x')|^2}{2l^2}\right)$. It then says that this covariance matrix follows GWP. I used to think that a covariance matrix computed from **linear covariance function ($K(x,x') = x^Tx'$)** , follows Wishart Distribution with appropriate parameters. My question is, how can we still assume the covariance to follow a Wishart distribution with squared exponential covariance function? Also, in general, what is the necessary condition for a covariance function to produce a Wishart distributed covariance matrix?


title: Medical Insurance Fraud Detection: Text analysis

text: I'm trying to analyse a dataset to detect fraudulent insurance claims. Unfortunately, other than basic demographics the rest of the claim is a free format OCR scanned text file made from documents submitted by the patient. These documents include lab tests, hospital bills etc., patient summaries. Right now, manual coding is one option but expensive. Is there any chance that applying textual analysis / text-mining might yield some predictive strength? It need not be perfect but even if it can reliably identify high fraud potential claims, that's a plus. Any ideas how I could model this?


title: Fitting Mancova model on factor scores

text: I have two matrices with multiple variables. The two matrices come from same sample.I would call them matrix A and Matrix B.Each of them have many correlated variables . I would want to evaluate the relationship between matrix A and matrix B taking matrix B as explanatory and matrix A as dependent . When i looked at the structure of Matrix A .I was able to describe it into 4 factors and matrix B into 3 factors. I fitted a mancova model using factor scores of matrix A and factor scores of matrix B. Is this a valid way to do the analysis.


title: Learning curve shows decreasing accuracy

text: I'm working on a random forest classifier with 10-folds CV to aestimate the hyperparameter 'mtry' (chosen by maximizing AUROC). I decided to pre-split the training set in 8 samples equals in size (about 100 observations each). I iterated the train process 8 times. At each iteration I added one sample to my actual training set. So, the first train model was fitted using a training set of 100 observations, the second had 200...the last one has the entire original training set (about 800 observations). At each iteration I measured the missclassification error (FP + FN) for both training and test set in %. I finally plot this learning curve. On the x axis I have the training set size, on the y axis the missclassification error in % (black for the training set, red for the test set). The gap between train and test performance suggests me I have high variance in the data (overfitting) but I really cannot understand...why does the training error start so high, then suddenly drop, then start to rise again as training set size increases? If I have a overfitting problem, how can I address it (close the gap) if I cannot obtain more training examples? ![enter image description here](http://i.stack.imgur.com/FKsnK.png)


title: Strange behavior of location-tests for different datasets

text: I have two datasets, $D1$ and $D2$, each containing many samples $S_i$ of 100 numbers, each sample from a different (unknown) population. For example: $S_1=\{0.125, 0.122, 0.053, \dots\} \\ S_2=\{0.001, 0.013, 0.034, \dots\} \\ \dots $ The are both quite large. $D1$ contains almost 6M samples, and $D2$ contains almost 500K. I randomly split each sample $S_i$ in two disjoint subsamples of size 50 each, let them be $S_i'$ and $S_i''$. For each subsample I run a t-test for the null hypothesis of mean equal zero, so I have two p-values: $p_i'$ and $p_i''$. My question of interest is: how often do I observe significance in $S_i''$ if I observed significance in $S_i'$. So, for different significance levels, I plot this proportion: $\displaystyle\frac {\#\{p_i' \leq \alpha~and~p_i''\leq\alpha\}} {\#\{p_i' \leq \alpha\}}$ With dataset $D1$ there is nothing rare: the lower the significance level the smaller the proportion because it is harder to observe significance. But when I plot it with dataset $D2$, I have this weird behavior in the plot: it kind of goes up with alpha, but then at some point it goes down with large alphas. I can't figure out why this is happening. The same happens with other tests such as Wilcoxon and Sign. Any ideas? ![enter image description here](http://i.stack.imgur.com/KPUOV.png)


title: How to interpret the coefficients from a beta regression?

text: I have some data that is bounded between 0 and 1. I have used the `betareg` package in R to fit a regression model with the bounded data as the dependent variable. My question is: how do I interpret the coefficients from the regression?


title: Centralization measures in weighted graphs

text: I'm using the `igraph` package in `R` to analyze network data. I'm currently trying to calculate some centrality measures for the vertices of my graph and the corresponding centralization scores. My network is both directed and weighted.               require(igraph)     set.seed(12152)          m <- expand.grid(from = 1:4, to = 1:4)     m <- m[m$from != m$to, ]     m$weight <- sample(1:7, 12, replace = T)          g <- graph.data.frame(m)      I have no trouble using the `closeness` function to obtain the closeness centrality for each vertex:               closeness(g, mode = "in")     closeness(g, mode = "out")     closeness(g, mode = "total")      However, it appears that the `centralization.closeness` function from `igraph` does not work for directed graphs. `igraph` does include a way to calculate a custom centralization from the individual centrality scores in a graph (the `centralize.scores` function), but that function requires the user to specify the theoretical maximum of the centrality measure, and it's not obvious to me what that would be in this weighted example (I believe the built-in `centralization.closeness.tmax` function in `igraph` assumes an unweighted graph). Does anyone know how to calculate a centralization score in a weighted graph? Is there a good way to accomplish this in R with `igraph` or some other package?


title: What does a significant intercept mean in ANOVA?

text: Consider a model with a continuous response variable and a categorical explanatory variable. I appreciate that in R, a summary.lm output of an anova on this data gives you rows that represent the mean value of each factor level. The significance stars represent the significance of the difference between the mean of each level and the "intercept", which represents the mean of the first level of the factor. What I am wondering is what do significance stars on this intercept term represent? Simply that the mean of this particular factor level is different from zero?


title: When combining p values, why not just averaging?

text: I recently learned about Fisher's method to combine p-values. This is based on the fact that p-value under the null follows a uniform distribution, and that $$-2\sum_{i=1}^n{\log X_i} \sim \chi^2(2n), \text{ given } X \sim \text{Unif}(0,1)$$ which I think is genius. But my question is why going this convoluted way? and why not (what is wrong with) just using mean of p-values and use central limit theorem? or median? I am trying to understand the genius of RA Fisher behind this grand scheme.


title: Minimum training sample size required for a classifier

text: What is the best method to determine the minimum number of training samples required for a classifier? I am only comparing one classifier (four class problem), discriminant function analysis (DFA) with different training sample sizes from different data sets. However, I could be considered to be comparing more than one classifier as I am also comparing DFA using the within matrix and separate matrix methods. The data is accelerometer data used to classify animal behaviours. The different data sets are all from the same original data but have had treatments applied, eg smoothed, outlier removed. I have run 100 interations of randomly selected increasing subsample sizes (10 - 2550) from each data set through each classifier. Averaging the classification accuracy or error rate for each subsample size separately for each data set gives me a learning curve for the classifier on each of the four classes. Using the area at which the learning curve levels out, I would like to suggest a possible minimum training sample size for each class required for a future study. Rather than visually looking at the learning curve and making a recommendation, is there an appropriate statistical test. Am I talking about data-driven feature selection? I dont need to extrapolate as it looks like I have enough training samples. I dont think ROC curves are applicable as I only have the classification rates from each run of each sub sample sizes rather than predicted membership probabilities which is what I think i need for ROC curves. I have looked over several posts here which have been useful in advancing my thinking. This post (How large a training set is needed?) relates to my requirement but I dont see it suggesting an appropriate test. This paper (http://jmlr.org/papers/volume7/demsar06a/demsar06a.pdf) suggests the Friedman test with the corresponding post-hoc tests for comparison of more classifiers over multiple data sets. However, it doesn't appear to take into account the increasing sub sample size from the same training data set, rather it uses completely independent data sets.


title: Distribution of deviances from logistic regression

text: I would like to know what is the theoretical distribution of the vector of (signed) deviances $d_i$ from a logistic regression. Failing that, is there a way to to transform them unto approximately normally distributed r.v.? My model for the data generating process is the simplest one, e.g. $$y_i\sim\mathcal{B}(p_i)$$ where $$p_i=\pmb x_i'\pmb\theta$$ and $\mathcal{B}$ is the Bernoulli distribution. so in R code:               n<-1000     p<-5     x<-matrix(rnorm(n*(p-1)),nc=p-1);     y<-cbind(1,x)%*%rep(1,p)     y0<-1/(exp(-y)+1);     y0<-rbinom(n,1,y0)          fit<-glm(y0~x,family="binomial")     y<-residuals(fit,type="deviance")     qqnorm(y)     qqline(y)


title: Adjusting regression for correlated errors-in-variables

text: I have a set of data points $\{x_i\}$. These data points are grouped so that (say) $i\in\{1,2,3\}$ is group $A$, $i\in\{4,5,6,7\}$ is group $B$, etc. I would like to test the null hypothesis of no linear relationship between group means and group variances. A naive approach would be to calculate, for each group, the sample mean $\overline{x}_g$ and sample variance $S^2_g$, then regress $$ \overline{x}_g = \beta S^2_g + \varepsilon_g $$ This doesn't work, since in small samples, $\hat{\beta}$ will be a biased estimate of $\beta$. For example, if $\beta=0$ but the $\{x_i\}$ are very positively skewed, a spuriously high $\hat{\beta}$ will be estimated. This is because   * Groups that do not contain positive outliers will have low $S^2_g$ and low $\overline{x}_g$   * Groups that do contain these outliers will have high $S^2_g$ and high $\overline{x}_g$.  Is there a general technique for solving this problem? (Cross-posted from Math Stack Exchange.)


title: Financial instruments ranking based on price movement (outliers, consistency)

text: I have a lot of financial instruments and I need to rank them based on price movement "stability", i.e. instruments with with big outliers should have the lowest ranks and instruments with stable and consistent moves should have the highest ranks. To solve this issue I use the ratio of the `Arithmetic Mean` to the `Standard Deviation` for every instrument and then sort them by this ratio. But there is a problem: this method works only roughly and some "bad" instruments have higher ranks than "good" instruments. Are there are different approaches for ranking? For data points I use $\ln(\text{High} / \text{Low})$ and period is 750 data points for every instrument.


title: Survey questions with many different response formats: Is it a problem?

text: A co-worker has asked me to help setup a spreadsheet that is used in the analysis of customer care emails. In particular, we are looking at three categories - the assignment of meta-data, content and communication skills. A typical example question may be: "How well is the email tailored towards the customers questions?" "How well did the employee answer the customers questions?" I'm not a statistician but have conducted a few surveys in the past. After seeing the initial survey, which was basically an ad-hoc collection of questions combined with varying scaling mechanisms for each question. I figured the first step would be to iterate the goals and work back from there; otherwise I'd be basically devising my own methods from scratch, which feels a little flakey. All this led me to re-design the survey using the Likert scale due to it's simplicity in implementation and analysis. However, without a deep statistical background I'm finding it hard to validate why it is better to use a proven methodology than a collection of different scales and measurements aggregated into a single number. How can I convince them of the benefits if using a 'methodology'/scale such as Likert rather than various choice of measurements?


title: When to use weight of evidence in a predictive model?

text: Suppose I am building a binary logistic regression model with around 400 variables then how weight of evidence(WOE) can be used in here? Also, I have known that WOE is helpful in dimension reduction in a model then why not go for factor analysis rather going for WOE? What should be the size of the number of variables in order to use WOE? Is there any other use of WOE in developing a model? What is the most appropriate scenario where WOE would be best suited to use?


title: Deriving posterior pdf in classical linear normal regression model under noninformative prior

text: Question: Assume the following classical linear normal regression model: \begin{gather*} y_{i} = \beta_1 x_{1i} + \beta_2 x_{2i} + \cdots + \beta_K x_{Ki} + e_i \\ \underbrace{\boldsymbol{y}}_{n \times 1} = \underbrace{\boldsymbol{X}}_{n \times K}\underbrace{\boldsymbol{\beta}}_{K \times 1} + \underbrace{\boldsymbol{e}}_{n \times 1} \end{gather*} where $e_i \sim iid \ N(0, \sigma^2)$ for $i = 1, 2, \cdots, n$ and $x_{1i} = 1$ for $i = 1, 2, \cdots, n$ Define $\boldsymbol{\beta} = (\beta_1, \beta_2, \cdots, \beta_K)'$ and assume a noninformative prior of the form $p(\boldsymbol{\beta}, \sigma) \propto \frac{1}{\sigma}$, derive the conditional posterior pdfs for $\boldsymbol{\beta}$ and $\sigma$, that is, $p(\boldsymbol{\beta}|\sigma, \mathbf{y})$ and $p(\sigma|\boldsymbol{\beta}, \mathbf{y})$. * * * I can derive the conditional posterior pdf for $p(\boldsymbol{\beta}|\sigma, \mathbf{y})$ but I am not sure how to derive it for $p(\sigma|\boldsymbol{\beta}, \mathbf{y})$. The answer to the latter is: $p(\sigma|\boldsymbol{\beta}, \mathbf{y})=\frac{2}{\Gamma\left(\frac{n}{2}\right)}\left(\frac{na}{2}\right)^{\frac{n}{2}} \frac{1}{\sigma^{n+1}}\exp\left[-\frac{na}{2\sigma^2}\right] $ where $a = \frac{1}{n}\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)'\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)$ * * * The following is my working to derive the posterior pdf for $p(\boldsymbol{\beta}|\sigma, \mathbf{y})$ We have: \begin{align*} p(\boldsymbol{e}|\sigma) & = \prod_{i=1}^n p(e_i) \\ & = \prod_{i=1}^n (2\pi \sigma^2)^{-\frac{1}{2}}\exp\left[-\frac{1}{2\sigma^2}e_i^2\right] \\ & = (2\pi \sigma^2)^{-\frac{n}{2}}\exp\left[-\frac{1}{2\sigma^2}\sum_{i=1}^n e_i^2\right] \\ & = (2\pi \sigma^2)^{-\frac{n}{2}} \exp\left[-\frac{1}{2\sigma^2} \boldsymbol{e'}\boldsymbol{e}\right] \end{align*} Since $p(\boldsymbol{y}) = p(\boldsymbol{e})\left|\frac{\partial \boldsymbol{e'}}{\partial \boldsymbol{y}} \right|$ and $\boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{e} \Longleftrightarrow \boldsymbol{e} = \boldsymbol{y} - \boldsymbol{X}\boldsymbol{\beta}$, then $\frac{\partial \boldsymbol{e'}}{\partial \boldsymbol{y}} = \frac{\partial \boldsymbol{y'}}{\partial \boldsymbol{y}} = \boldsymbol{I}_{n}$. Hence, $\left|\frac{\partial \boldsymbol{e'}}{\partial \boldsymbol{y}} \right| = \left| \boldsymbol{I}_{n}\right| = 1$. So: \begin{align*} p(\boldsymbol{y}|\boldsymbol{\beta}, \sigma) & = p(\boldsymbol{e}) \\ & = (2\pi \sigma^2)^{-\frac{n}{2}} \exp\left[-\frac{1}{2\sigma^2}\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)'\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)\right] \ \ \text{since} \ \ \boldsymbol{e} = \boldsymbol{y} - \boldsymbol{X}\boldsymbol{\beta} \end{align*} Thus we have: $$L(\boldsymbol{\beta}, \sigma|\boldsymbol{y}) \propto (\sigma^2)^{-\frac{n}{2}} \exp\left[-\frac{1}{2\sigma^2}\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)'\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)\right]$$ Given $p(\boldsymbol{\beta}, \sigma) \propto \frac{1}{\sigma}$, we have: \begin{align*} p(\boldsymbol{\beta}, \sigma|\boldsymbol{y}) & \propto L(\boldsymbol{\beta}, \sigma|\boldsymbol{y}) \times p(\boldsymbol{\beta}, \sigma) \\ & \propto \frac{1}{\sigma^n}\exp\left[-\frac{1}{2\sigma^2}\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)'\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)\right] \times \frac{1}{\sigma} \\ & = \frac{1}{\sigma^{n+1}} \exp\left\{-\frac{1}{2\sigma^2}\left[(n-K)\widehat{\sigma}^2 + (\boldsymbol{\beta} - \boldsymbol{b})'\boldsymbol{X}'\boldsymbol{X}(\boldsymbol{\beta}-\boldsymbol{b})\right] \right\} \\ & \text{where} \ \ \widehat{\sigma}^2 = \frac{1}{n-K} (\boldsymbol{y} - \boldsymbol{X}\boldsymbol{b})'(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{b}) \ \ \text{and} \ \ \boldsymbol{b} = (\boldsymbol{X}'\boldsymbol{X})^{-1}\boldsymbol{X}'\boldsymbol{y} \\ & = \underbrace{\frac{1}{\sigma^K} \exp\left[-\frac{1}{2\sigma^2}(\boldsymbol{\beta} - \boldsymbol{b})'\boldsymbol{X}'\boldsymbol{X}(\boldsymbol{\beta}-\boldsymbol{b})\right]}_{f^n(\boldsymbol{\beta}|\sigma, \boldsymbol{y})} \underbrace{\frac{1}{\sigma^{n-K+1}} \exp\left[-\frac{(n-K)\widehat{\sigma}^2}{2\sigma^2}\right]}_{f^n(\sigma|\mathbf{y})} \end{align*} Note that $f^n(\boldsymbol{\beta}|\sigma, \boldsymbol{y})$ is a kernel for $\boldsymbol{\beta}$, where $\boldsymbol{\beta} \sim MVN(\boldsymbol{b}, \sigma^2(\boldsymbol{X}'\boldsymbol{X})^{-1})$ and $f^n(\sigma|\mathbf{y})$ is a kernel for $\sigma$ where $\sigma \sim IG(v, \widehat{\sigma})$ and $v = n-K$, $\widehat{\sigma}^2 = \frac{1}{v} \left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)'\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)$. Thus we have: $$p(\boldsymbol{\beta}, \sigma|\mathbf{y}) = p(\boldsymbol{\beta}|\sigma, \boldsymbol{y})p(\sigma|\mathbf{y})$$ where: \begin{gather*} p(\boldsymbol{\beta}|\sigma, \boldsymbol{y}) = (2\pi)^{-\frac{K}{2}}|\sigma^2(\boldsymbol{X}'\boldsymbol{X})^{-1}|^{-\frac{1}{2}}\exp\left[-\frac{1}{2\sigma^2}(\boldsymbol{\beta} - \boldsymbol{b})'\boldsymbol{X}'\boldsymbol{X}(\boldsymbol{\beta}-\boldsymbol{b})\right] \\ p(\sigma|\mathbf{y}) = \frac{2}{\Gamma\left(\frac{v}{2}\right)} \left(\frac{v\widehat{\sigma}^2}{2}\right)^{\frac{v}{2}}\frac{1}{\sigma^{v+1}}\exp\left[-\frac{v\widehat{\sigma}^2}{2\sigma^2}\right] \end{gather*}


title: Should I include interactions in ridge regression?

text: I'm running the a ridge regression broadly analogous to the one at http://www.mathworks.com.au/help/stats/ridge.html. However, I have a lot more predictors (12) than the example (which only has 3). Thus when I include the interactions the number of parameters blows out. What should I do? I don't have any particular reason for thinking there would/would not be interactions.


title: Finding a sufficient statistic

text: Consider an i.i.d. sample $(X_{1},\ldots, X_{n})$ where the $X_{i}$ have density $f(x) = k \cdot \exp(−(x − θ)^4)$ with $x$ and $\theta$ real, obtain the sufficient statistic and its dimension. What does this example show? Is the statistic that you have found also complete?


title: is there any toy similarity matrix to test my random walk code?

text: I am looking for a toy similarity matrix of any kind of entity class to test my random walk code. Is there any set like this or how can I create my own set? I need to interpret the results of the random walk that I run to find the dominant instances of the data.


title: AR(1) process with heteroscedastic measurement errors

text: **1. The problem** I have some measurements of a variable $y_t$, where $t=1,2,..,n$, for which I have a distribution $f_{y_t}(y_t)$ obtained via MCMC, which for simplicity I'll assume is a gaussian of mean $\mu_t$ and variance $\sigma_t^2$. I have a physical model for those observations, say $g(t)$, but the residuals $r_t = \mu_t-g(t)$ appear to be correlated; in particular, I have physical reasons to think that an $AR(1)$ process will suffice to take into account the correlation, and I plan on obtaining the coefficients of the fit via MCMC, for which I need **the likelihood**. I think the solution is rather simple, but I'm not pretty sure (it seems so simple, that I think I'm missing something). **2. Deriving the likelihood** A zero-mean $AR(1)$ process can be written as: $$X_t = \phi X_{t-1}+\varepsilon_t,\ \ \ (1)$$ where I'll assume $\varepsilon_t\sim N(0,\sigma_w^2)$. The parameters to be estimated are, therefore, $\theta = \{\phi,\sigma_w^2\}$ (in my case, I also have to add the parameters of the model $g(t)$, but that's not the problem). What I observe, however, is the variable $$R_t = X_t+\eta_t,\ \ \ (2)$$ where I'm assuming $\eta_t\sim N(0,\sigma_t^2)$, and the $\sigma_t^2$ are known (the measurement errors). Because $X_t$ is a gaussian process, $R_t$ is also. In particular, I know that $$X_1 \sim N(0,\sigma_w^2/[1-\phi^2]),$$ therefore, $$R_1 \sim N(0,\sigma_w^2/[1-\phi^2]+\sigma_t^2).$$ The next challenge is to obtain $R_t|R_{t-1}$ for $t\neq 1$. To derive the distribution of this random variable, note that, using eq. $(2)$ I can write $$X_{t-1} = R_{t-1}-\eta_{t-1}.\ \ \ (3)$$ Using eq. $(2)$, and using the definition of eq. $(1)$, I can write, $$R_{t} = X_t+\eta_t = \phi X_{t-1}+\varepsilon_{t}+\eta_t.$$ Using eq. $(3)$ in this last expression, then, I obtain, $$R_{t} = \phi (R_{t-1}-\eta_{t-1})+\varepsilon_{t}+\eta_t,$$ thus, $$R_t|R_{t-1} = \phi (r_{t-1}-\eta_{t-1})+\varepsilon_{t}+\eta_t,$$ and, therefore, $$R_t|R_{t-1} \sim N(\phi r_{t-1},\sigma_w^2+\sigma_t^2-\phi^2\sigma^2_{t-1}).$$ Finally, I can write the likelihood function as $$L(\theta) = f_{R_1}(R_1=r_1) \prod_{t=2}^{n} f_{R_{t}|R_{t-1}}(R_t=r_t|R_{t-1}=r_{t-1}),$$ where the $f(\cdot)$ are the distributions of the variables that I just defined, .i.e., defining $\sigma'^2 = \sigma_w^2/[1-\phi^2]+\sigma_t^2,$ $$f_{R_1}(R_1=r_1) = \frac{1}{\sqrt{2\pi \sigma'^2}}\text{exp}\left(-\frac{r_1^2}{2\sigma'^2}\right),$$ and defining $\sigma^2(t) = \sigma_w^2+\sigma_t^2-\phi^2\sigma^2_{t-1}$, $$f_{R_{t}|R_{t-1}}(R_t=r_t|R_{t-1}=r_{t-1})=\frac{1}{\sqrt{2\pi \sigma^2(t)}}\text{exp}\left(-\frac{(r_t-\phi r_{t-1})^2}{2\sigma^2(t)}\right)$$ **3. Questions**   1. **Is my derivation ok?** I don't have any resources to compare other than simulations (which seem to agree), and I'm not a statistician!   2. **Are there any derivation of this kind of things in the literature for $MA(1)$ proccesses or $ARMA(1,1)$ proccesses?** A study for $ARMA(p,q)$ proccesses in general that could be particularized to this case would be nice.


title: Is the logistic regression model integrable into the potential outcome framework?

text: (1) I preprocessed data with a matching method (CEM) (2) On this data (binary dependent variable, explanatory variables) , I ran a logit model and estimate predicted probabilities or odds ratios. I would like to integrate my estimates of the logit model into the potential outcome framework (I would like get an average treatment effect for the treated (ATT) of my explanatory variable) I know it's possible for OLS (I found it in Gelman & Hill 2007, Chapter 9), but I did not find anything on logit.


title: How to calculate a correlation between variables which are proportions of observations varying in sample size

text: Is it possible to calculate the correlation between two proportions, while accounting for the individual sample size that is used to calculate the proportion per individual? For example: How do I calculate the correlation between the proportion of heads from one coin, and the proportion of heads from another coin, when tossed by the same individual. Not every individual tosses every coin as often (some individuals even toss one of the two coins once). Is this even possible?


title: In SPSS two factor ANOVA, how to compare control groups across other factors?

text: I have two factors: **Brainregion (4 levels)** and **Group (5 levels)**. Glutamate is my dependent variable. I am running General Linear Models in SPSS (two factor ANOVA). Can I get the output in such a way that: I can compare means (post-hoc test) of control with other levels such as alcohol, EW etc. (Group) across CC, CB etc. (Brainregion)?


title: Research using basic OLS methods

text: I am looking for papers, preferably in economics, where simple OLS is used. The reason is that to convey the workings of linear regression, I want to have real data from real papers. Most research uses simultaneous models or other bells and whistles. Also few researchers release their data. I could generate my own data, but as I said, I want to have the whole story of the data and research behind it. Wooldridge (Introductory Econometrics) does a lot of simple analyses from past papers, but as far as I know, these are just simplified version of the real models.


title: Bayesian inference on a sum of iid real-valued random variables

text: Let $X_1$, $X_2$, ..., $X_n$ be iid RV's with range $[0,1]$ but unknown distribution. (I'm OK with assuming that the distribution is continuous, etc., if necessary.) Define $S_n = X_1 + \cdots + X_n$. I am given $S_k$, and ask: What can I infer, in a Bayesian manner, about $S_n$? That is, I am given the sum of a sample of size $k$ of the RV's, and I would like to know what I can infer about the distribution of the sum of all the RV's, using a Bayesian approach (and assuming reasonable priors about the distribution). If the support were $\{0,1\}$ instead of $[0,1]$, then this problem is well- studied, and (with uniform priors) you get beta-binomial compound distributions for the inferred distribution on $S_n$. But I'm not sure how to approach it with $[0,1]$ as the range... **Full disclosure** : I already posted this on MathOverflow, but was told it would be better posted here, so this is a re-post.


title: To what degree is FWE correction required for multiple permutation based t-tests

text: I am performing multiple t-tests to compare conditions of data. I have had to resort to permutation based t-tests (based on difference of means) due to high variance between samples. My question is if I am using permutation based t-tests is FWE correction still required? My understanding is that because is it non-parametric the probability of a false negative is not contingent on the number of tests applied to the same set of data. So FWE correction should not be required? * * * I will expand on my specific issue. I am comparing scan-path data for subjects looking at two different types of faces (inverted and upright). I need to compare a number of measures for different face areas (eyes, nose , mouth, ears and 'other') between the areas. For this I have set up a t-test for each area, for each measure, and then correct for multiple comparisons. For example, I am testing the number of fixations, and fixation duration for each face area. I am constrained these tests in order to be able to compare to previous work. I have not been able to rely on any omnibus test thus far because of the non normality of my data greatly reduced the power of any omnibus test I have tried.


title: Hierarchical modeling when data not missing at random

text: Gelman & Hill (2006) explain how a multilevel regression model can compensate for data sparsity by pooling parameter estimates. But what happens if the data are not sparse at random? For example, suppose we are evaluating a new blood pressure medication. On every visit we measure a patient's blood pressure, administer the treatment, and measure blood pressure again. We may have fewer opportunities to observe the effect of a treatment in some patients than in others, so we may want to use pooling. But the number of observations per patient could be related to the fact that the treatment is ineffective for some patients (so they do not want to come in), or causes side effects two hours later (so they do not want to come in), or the medication puts patients in a coma and they drop out of the study. Is pooling no longer appropriate in such a study?


title: How do you calculate the standard deviation on a multiplicative scale for a distribution that has been transformed logarithmically?

text: I know the value for the 16% quartile, so I know the additive deviation for the given distribution. How do I find the deviation of the log of the given distribution on a multiplicative scale?


title: Correlated errors between variables [between / within subjects, within variables, everything]!

text: I am trying to figure out good model to analyze social dynamics using a model. I have commenting data from a website. Such as 500 articles, 10000 users, and their 75000 comments. Every user does not comment on every article. The opinions posed by each user is linked to the opinion of the previous users, as well as the opinion of the article itself. The opinions of previous users are also linked to the previous users and the article. I want to estimate the impact of opinion of previous user, average user opinion, article opinion and variation of opinions on the current opinion. Can you suggest a methodology for such a problem? Everything is correlated. I also would like to ask, what software and code would you use to model a problem like this? Codes are very very appreciated.


title: Comparing logistic regression models of winning percentage

text: I previously asked this question: link I've generated the two logistic regression models described in the link using all available data (8 seasons for all 30 teams), and now I want to decide which is a better predictor of win percentage in the coming season, winning percentage in the previous season or Pythagorean expectation from the previous season. I know that if you want to evaluate the predictive abilities of a model, it's often good to evaluate on data different from that which generated the model. I wonder how important that is here, given that the models are of the same type, created from the same data. It seems as though we could say that the model which was a better fit with respect to the data which created both models is better in this case. So, the question is: 1) do I need to see how the models compare on a different data set, and 2) what's the standard measure for model comparison in a situation like that I've described. Basically, how do I know which is better?


title: Have I computed these likelihood ratios correctly?

text: I'm the author of the ez package for R, and I'm working on an update to include automatic computation of likelihood ratios (LRs) in the output of ANOVAs. The idea is to provide a LR for each effect that is analogous to the test of that effect that the ANOVA achieves. For example, the LR for main effect represents the comparison of a null model to a model that includes the main effect, the LR for an interaction represents the comparison of a model that includes both component main effects versus a model that includes both main effects _and_ their interaction, etc. Now, my understanding of LR computation comes from Glover & Dixon (PDF), which covers basic computations as well as corrections for complexity, and the appendix to Bortolussi & Dixon (appendix PDF), which covers computations involving repeated-measures variables. To test my understanding, I developed this spreadsheet, which takes the dfs & SSs from an example ANOVA (generated from a 2*2*3*4 design using fake data) and steps through the computation of the LR for each effect. I would really appreciate it if someone with a little more confidence with such computation could take a look and make sure I did everything correctly. For those that prefer abstract code, here is the R code implementing the update to ezANOVA() (see esp. lines 15-95).


title: OLS is BLUE. But what if I don't care about unbiasedness and linearity?

text: The Gauss-Markov theorem tells us that the OLS estimator is the best linear unbiased estimator for the linear regression model. But suppose I don't care about linearity and unbiasedness. Then is there some other (possible nonlinear/biased) estimator for the linear regression model which is the most efficient under the Gauss-Markov assumptions or some other general set of assumptions? There is of course one standard result: OLS itself is the best unbiased estimator if in addition to the Gauss-Markov assumptions we also assume that the errors are normally distributed. For some other particular distribution of errors I could compute the corresponding maximum-likelihood estimator. But I was wondering if there is some estimator which is better-than-OLS in some relatively general set of circumstances?


title: Is there a test in R whose null hypothesis is that there is a difference between the means of two populations?

text: I want to conduct a statistical test that will allow me to say that two populations are similar rather than different. As statisticians know, I can't do this with tests like the T test because the null hypothesis of the test is that the populations' means are equal. Hence my question.


title: Calculate boundary for MAE given RMSE

text: E.g. from the Netflix prize I know that the best RMSE = 0.8563 where the test dataset has a size of n=1,408,789. Can I calculate a boundary for the MAE. If not, why can't I calculate a boundary? I started like that: RMSE = $\sqrt{\frac{\sum_{i=1}^{n}(y_{i}-\hat y_{i})^2}{n}}=0.8563$ Now I want to calculate the MAE or a boundary for the MSE: MAE = $\frac{1}{n}\sum_{i=1}^{n}{|y_{i}-\hat y_{i}|}$ I know that: $\sum_{i=1}^{n}(y_{i}-\hat y_{i})^2=RMSE^2*n$ from the first formula. Is it possible to determine a boundary for the MAE e.g. MAE < 1?


title: How to construct a reasonable prior and likelihood for Bayes modelling?

text: To apply Bayes inference for data analysis or machine learning, we have to construct _prior_ and _likelihood_ , right? But if I fail to come up with a reasonable _prior_ and _likelihood_ , then the Bayes model will not be meaningful, right? I wonder is there any technique could be utilized to construct reasonable Bayes models? Let me make it concrete by an example, given a dataest with features $X$ and target $y$, where $y=\{0,1\}$ and $X$ is composed of $p$ variables. This is an ordinary binary classification problem, to do Bayes inference, what and how should I specify for the prior and likelihood? One more question, if we could not specify a reasonable prior, then what's the point of using Bayes modelling?


title: Statistical significance and winner of multivalent test

text: This should be a fairly straight forward, but I can't seem to figure out the answer. I'm doing some A/B(C/D/E...) testing on a website and measuring impressions and clicks. What method should I be using to determine the statistically significant winner? This is sample data that I would have. There could be any number of tests, all being displayed in relatively-equal distribution and measuring clicks as a successful result. Such as:               Test #  |  Impressions  | Clicks       1     |      50       |   12       2     |      55       |   15       3     |      53       |   30       4     |      49       |   22      What algorithm should I be using to determine the winner in a statistically significant way?


title: How to compare two matrices?

text: I am working on Markov transition matrices. I would like to find a statistical test to compare them. The first matrix is considered the population transition matrix and the second one is obtained by estimation from data (counting the transitions). My first idea was to use a Chi square test on each line of the two matrices. But I have an issue: some of the transitions might be null. When this happens the Chi square test does not work any more. If someone has an idea, it would be great. In addition, if you have some academic paper on the subject I would like to have the references.


title: Using lmer for repeated-measures linear mixed-effect model

text: EDIT 2: I originally thought I needed to run a two-factor ANOVA with repeated measures on one factor, but I now think a linear mixed-effect model will work better for my data. I think I nearly know what needs to happen, but am still confused by few points. The experiments I need to analyze look like this:   * Subjects were assigned to one of several treatment groups   * Measurements of each subject were taken on multiple days   * So:      * Subject is nested within treatment     * Treatment is crossed with day (each subject is assigned to only one treatment, and measurements are taken on each subject on each day) My dataset contains the following information:   * Subject = blocking factor (random factor)   * Day = within subject or repeated measures factor (fixed factor)   * Treatment = between subject factor (fixed factor)   * Obs = measured (dependent) variable **UPDATE** OK, so I went and talked to a statistician, but he's an SAS user. He thinks that the model should be: **Treatment + Day + Subject(Treatment) + Day*Subject(Treatment)** Obviously his notation is different from the R syntax, but this model is supposed to account for:   * Treatment (fixed)   * Day (fixed)   * the Treatment*Day interaction   * Subject nested within Treatment (random)   * Day crossed with "Subject within Treatment" (random) So, is this the correct syntax to use?               m4 <- lmer(Obs~Treatment*Day + (1+Treatment/Subject) + (1+Day*Treatment/Subject), mydata)      I'm particularly concerned about whether the Day crossed with "Subject within Treatment" part is right. Is anyone familiar with SAS, or confident that they understand what's going on in his model, able to comment on whether my sad attempt at R syntax matches? Here are my previous attempts at building a model and writing syntax (discussed in answers & comments):               m1 <- lmer(Obs ~ Treatment * Day + (1 | Subject), mydata)      How do I deal with the fact that subject is nested within treatment? How does `m1` differ from:               m2 <- lmer(Obs ~ Treatment * Day + (Treatment|Subject), mydata)     m3 <- lmer(Obs ~ Treatment * Day + (Treatment:Subject), mydata)      and are `m2` and `m3` equivalent (and if not, why)? Also, do I need to be using nlme instead of lme4 if I want to specify the correlation structure (like `correlation = corAR1`)? According to Repeated Measures, for a repeated-measures analysis with repeated measures on one factor, the covariance structure (the nature of the correlations between measurements of the same subject) is important. When I was trying to do a repeated-measures ANOVA, I'd decided to use a Type II SS; is this still relevant, and if so, how do I go about specifying that? Here's an example of what the data look like:               mydata <- data.frame(       Subject  = c(13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 29, 30, 31, 32, 33,                     34, 35, 36, 37, 38, 39, 40, 62, 63, 64, 65, 13, 14, 15, 16, 17, 18,                     19, 20, 21, 22, 23, 24, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39,                     40, 62, 63, 64, 65, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24,                     29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 62, 63, 64, 65),        Day       = c(rep(c("Day1", "Day3", "Day6"), each=28)),        Treatment = c(rep(c("B", "A", "C", "B", "C", "A", "A", "B", "A", "C", "B", "C",                            "A", "A", "B", "A", "C", "B", "C", "A", "A"), each = 4)),        Obs       = c(6.472687, 7.017110, 6.200715, 6.613928, 6.829968, 7.387583, 7.367293,                      8.018853, 7.527408, 6.746739, 7.296910, 6.983360, 6.816621, 6.571689,                      5.911261, 6.954988, 7.624122, 7.669865, 7.676225, 7.263593, 7.704737,                      7.328716, 7.295610, 5.964180, 6.880814, 6.926342, 6.926342, 7.562293,                      6.677607, 7.023526, 6.441864, 7.020875, 7.478931, 7.495336, 7.427709,                      7.633020, 7.382091, 7.359731, 7.285889, 7.496863, 6.632403, 6.171196,                      6.306012, 7.253833, 7.594852, 6.915225, 7.220147, 7.298227, 7.573612,                      7.366550, 7.560513, 7.289078, 7.287802, 7.155336, 7.394452, 7.465383,                      6.976048, 7.222966, 6.584153, 7.013223, 7.569905, 7.459185, 7.504068,                      7.801867, 7.598728, 7.475841, 7.511873, 7.518384, 6.618589, 5.854754,                      6.125749, 6.962720, 7.540600, 7.379861, 7.344189, 7.362815, 7.805802,                      7.764172, 7.789844, 7.616437, NA, NA, NA, NA))


title: How can I generate uniformly distributed points on a circle?

text: I am looking to generate 450 data points in R. There are three distinct sets 150 of each distributed in a circular band with different radii (at 1, 2.8 and 5). In particular, I'm looking to reproduce the first graph on p546 of The Elements of Statistical Learning. ![enter image description here](http://i.stack.imgur.com/rMILV.png) I would be very grateful for some assistance in the R code that could generate this data. Thanks!


title: How to build a relevant text classifier?

text: I would like to build a message classification system which classifies a given message into either of 2 class - Relevant/Not. I don't have any labelled dataset. I only have certain keywords which defines relevance. I can include more keywords (synonyms, nearly same words, etc). For example: If I want to define relevance for Automobile the relevant keywords I'll be having are - Cars, bikes, automobiles, gas, tires, engine, oil, ... Now, the classifier should check a new message (email/SMS) and classify them whether the message belongs to this category or not.               Message1: Always use abc car oil.     Message2: My dish had too much of oil.      In this case, _Message1_ is relevant to automobile domain rather _Message2_. Please, advice me how to build a training dataset and how do I generalize this model for any domain?


title: Testing time homogeneous Markov chains

text: I am working with different transition diagrams and want to calculate the likelihood ratio statistic for testing time-homogeneous. I saw that there are already some comparable questions, but I still don't understand the method. For example, if I have the following transition diagram (probabilities) for 1 year:               From state 1 to state 1,2,3: 0,3 0,4 0,3   <- (Based on 80 companies)     From state 2 to state 1,2,3: 0,2 0,7 0,1   <-(Based on 500 companies)     From state 3 to state 1,2,3: 0,5 0,2 0,3   <- (Based on 300 companies).      And the following diagram for 2 years:               From state 1 to state 1,2,3: 0,2 0,5 0,3   <-(Based on 100 companies)     From state 2 to state 1,2,3: 0,1 0,2 0,7   <-(Based on 1000 companies)     From state 3 to state 1,2,3: 0,6 0,2 0,2   <-(Based on 500 companies)      How can I calculate the likelihood ratio statistic with the same method as done in the article Testing simple Markov structures for Credit Rating transitions (http://www.occ.gov/publications/publications-by-type/occ-working- papers/2008-2000/wp2004-3.pdf). I hope that someone can help me with this question based on this example!