title: Stata "Regress" Command and Finding Average Difference

text: I want to determine the average difference between 1998 treatment schools (group1=1) and comparison schools (group1=0) in the following two characteristics:Proportion of female children(“female”),Average child year of birth (“yob”) I am new to Stata, so I don’t know what to type in the command window. What should I type? I've learned how to calculate the regression between an x and y variable where it's like               regress y x      but I don't know how to calculate regression between group 1 and the comparison group when compared to another characteristic. I need 2 regressions, so would set would I type?               regress group1=1 yob female     regress group1=0 yob female          regress group1=0 group2 yob     regress group1=1 group2 female      or something else? and is the average difference given to me or is that something I need to calculate on my own? All I know is that I have to use the `regress` command. Here is the data set:               schid   group1  group2  group3  part98  female  yob     115 1   0   0   .7439648    .3834144    1987.234     117 1   0   0   .7443243    .3705711    1986.428     142 1   0   0   .7687114    .3864889    1987.516     143 1   0   0   .699026 .3519214    1987.588     158 1   0   0   .5439946    .3648172    1987.836     165 1   0   0   .7515087    .3641604    1987.808     191 1   0   0   .6320185    .3460737    1987.057     206 1   0   0   .8539617    .3606257    1987.215     207 1   0   0   .8054298    .3759444    1986.204     210 1   0   0   .8510154    .3543401    1986.542     215 1   0   0   .6107619    .387303 1987.144     216 1   0   0   1   .4003763    1986.481     219 1   0   0   .8350134    .3801926    1986.603     220 1   0   0   .6586895    .3857815    1986.905     265 1   0   0   1   .3786209    1986.815     267 1   0   0   .9961755    .3797742    1986.82     268 1   0   0   .9998046    .3439624    1986.505     274 1   0   0   .9886107    .448456 1986.628     279 1   0   0   .9877396    .387717 1986.717     281 1   0   0   .994603 .4124718    1985.863     286 1   0   0   .9647818    .3563205    1986.472     288 1   0   0   .9494266    .3386317    1985.715     290 1   0   0   .9806792    .4015853    1986.362     296 1   0   0   .9931837    .4222261    1986.637     297 1   0   0   1   .3925168    1986.104     109 0   1   0   .7627826    .3316738    1986.306     116 0   1   0   .6412753    .3338825    1987.693     118 0   1   0   .5040858    .3230579    1986.826     119 0   1   0   .6373736    .3749958    1986.878     133 0   1   0   .5865115    .4493426    1988.716     134 0   1   0   .3969593    .4464827    1987.21     141 0   1   0   .619329 .3765736    1987.549     164 0   1   0   .7135595    .2955894    1987.438     166 0   1   0   .7807278    .3349105    1987.209     182 0   1   0   .6992276    .3615404    1987.088     183 0   1   0   .7401822    .3925673    1988.551     192 0   1   0   .5432901    .3736948    1987.006     203 0   1   0   .6527734    .3902833    1986.182     208 0   1   0   .6212902    .4109006    1986.253     221 0   1   0   .8271556    .9911665    1988.778     262 0   1   0   .9702736    .3943461    1987.15     263 0   1   0   .6658556    .413118 1987.252     264 0   1   0   .9434038    .3477578    1986.569     269 0   1   0   .9521447    .3844051    1986.523     273 0   1   0   .8159674    .4032201    1987.036     277 0   1   0   .74529  .4287328    1986.644     284 0   1   0   .943206 .334808 1987.396     289 0   1   0   .9922523    .3737274    1986.665     292 0   1   0   .992666 .421622 1986.812     293 0   1   0   .8902047    .4432617    1985.719     108 0   0   1   .6897349    .2809297    1986.995     140 0   0   1   .7615737    .3513303    1987.509     144 0   0   1   .585315 .3754444    1986.324     167 0   0   1   .5862138    .3253592    1987.472     168 0   0   1   .7034473    .3261905    1987.37     189 0   0   1   .5167804    .2934488    1986.522     190 0   0   1   .6437579    .2906108    1986.846     193 0   0   1   .8333147    .3409633    1986.912     201 0   0   1   .7961279    .4090333    1986.359     202 0   0   1   .8110762    .4032104    1986.049     205 0   0   1   .8138164    .4128731    1986.051     211 0   0   1   .9859384    .3808163    1985.383     212 0   0   1   .8145156    .4102123    1986.173     213 0   0   1   .9905288    .3850739    1985.755     218 0   0   1   .9398516    .3266357    1986.124


title: Predict 2 responses from two co-variates

text: I'm not quite sure how I should fit a model that has two responses. The data consists of target (x,y) co-ordinates and actual (x,y) co-ordinates. I would like to fit a model to predict a new set of points from the target, however having two variables is putting me off. There also appears to be no relationship between points so I'm not even sure what kind of modelling to do. Would a multivariate Gaussian be appropriate? edit: here is the link to the data: https://dl.dropbox.com/u/12945652/data.csv               sessionID       targetX         targetY         touchX          touchY          name    phone     253             0.034375        0.564606742     0.034375        0.584269663     Subject15       Iphone4     253             0.60625         0.901685393     0.625           0.907303371     Subject15       Iphone4     253             0.2875          0.764044944     0.253125        0.783707865     Subject15       Iphone4     253             0.25            0.814606742     0.228125        0.823033708     Subject15       Iphone4     253             0.059375        0.657303371     0.04375         0.679775281     Subject15       Iphone4     253             0.06875         0.91011236      0.05            0.963483146     Subject15       Iphone4     253             0.625           0.603932584     0.68125         0.643258427     Subject15       Iphone4     253             0.265625        0.530898876     0.26875         0.533707865     Subject15       Iphone4     253             0.596875        0.974719101     0.603125        0.963483146     Subject15       Iphone4     253             0.290625        0.533707865     0.25            0.530898876     Subject15       Iphone4     253             0.946875        0.575842697     0.98125         0.564606742     Subject15       Iphone4     253             0.825           0.54494382      0.909375        0.542134831     Subject15       Iphone4     253             0.4625          0.556179775     0.4875          0.570224719     Subject15       Iphone4     253             0.6875          0.814606742     0.70625         0.828651685     Subject15       Iphone4     253             0.734375        0.873595506     0.775           0.91011236      Subject15       Iphone4


title: Kernel logistic regression

text: I heard Kernel Logistic Regression is a classical combination of kernel methods and Logistic regression, but I cannot find any major reference (book, or paper) on this topic. Can you give me any suggestions? Thanks.


title: How to predict logistic model by accounting for the error distribution?

text: Suppose I have a logistic model such that:               set.seed(123)     df <- data.frame(y=rbinom(100,1,0.5),                      x1=rnorm(100,10,2),                      x2=rbinom(100,20,0.6))          fit <- glm(y~x1+x2, data=df, family="binomial")      I would predict $\hat{y}$ for m times based on the `fit` by accounting the error distribution $Error \sim N(0,\sigma)$, which means for every set of $\hat{y_{i}}, i=1\cdots m$: $\hat{y_{i}}=fixed(\hat{y_{i}})+error_{i}$. Would somebody here know how to realize that?


title: How to use a meta-analysis to assess the equivalence of two different procedures?

text: I'm making a comparison of two surgical procedures, and I want to know how I can analyze the data I have for a meta-analysis. The overall surgical procedure is the same, but some authors choose to modify the procedure a bit. So basically I have a large number of studies with the standard procedure and a smaller number with a modified procedure. The outcome is percent complications, e.g., 20% of the study population had mild complications and 10% had severe complications. What would be the best way to show if the modified procedure is as good as the standard procedure?


Are goodness of fit and lack of fit the same? If they are tests, is one's null the other's alternative?


title: What is EVIR's dgpd output?

text: I am trying to evaluate the output of the "dgpd" function in the "evir" package in R. I am able to pass a vector of exceedences as well as gpd parameter estimates, code below, but the vector of output has large numbers. What unit of measure is ssociated with the dgpd output?               exceedence.vector <- subset(negRets, negRets > 0.02)     exceedence.vector          densityOfExceedence <- dgpd(exceedence.vector,xi=tau,beta=sigma)     densityOfExceedence


title: Generating non-homogeneous spatial Gaussian data

text: I want to generate a spatial data following multivariate Gaussian distribution. However, I don't want it to be homogeneous, meaning I don't want the correlation/covariance to be homogeneous. I want it to be heterogeneous. Any suggestions on how to generate such data? Here is an example of what I mean by heterogeneous covariance matrix                  25.8621   25.1207   24.6305   24.3867        25.1207   25.3719   24.8768   24.6305        24.6305   24.8768   25.3719   25.1207        24.3867   24.6305   25.1207   25.8621      There are four variables. As you can see the marginal variance is different which means C(0) is dependent on location. So does C(h) for any h. This is one method of generating heterogeneous data. But I want some spatial data consisting of like land, lakes which have very less covariance between each other.


title: Sample size calculation for univariate logistic regression

text: How does one calculate the sample size needed for a study in which a cohort of subjects will have a single continuous variable measured at the time of a surgery and then two years later they will be classified as functional outcome or impaired outcome. We would like to see if that measurement could have predicted the bad outcome. At some point we may want to derive a cut point in the continuous variable above which we would try to intervene to diminish the probability of the impaired outcome. Any ideas? Any R implementation.


title: Gentler approach to Bayesian statistics

text: I recently started reading "Introduction to Bayesian Statistics" 2nd Edition by Bolstad. I've had an introductory stats class that covered mainly statistical tests and am almost through a class in regression analysis. What other books can I use to supplement my understanding this one? I've made it through the first 100-125 pages fine. Afterwards the book begins to talk hypothesis testing which is what I'm very excited to cover but there a couple of things throwing me:   * The use of probability density functions in calculations. In other words how to evaluate such equations.   * This whole sentence: "Suppose we use a beta(1,1) prior for pi. Then given y=8, the posterior density is beta(9,3). The posterior probability of the null hypothesis is..." I _believe_ beta(1,1) refers to a PDF where the mean is 1 and the stdev is 1? I don't get how it would change to a beta(9,3) as a posterior density function. I do get the concept of priors vs posteriors and understand how to apply them using a table manually. I get (I think!) that pi represents the supposed population proportion or probability. I don't get how to connect this together with data I would run into on a day to day basis and get results.


title: How to apply neural networks on multi-label classification problems?

text: **Description:** Let the problem domain be document classification where there exists a set of feature vectors, each belonging to 1 or more classes. For example, a document `doc_1` might belong to `Sports` and `English` categories. **Question:** Using neural network for classification, what would the label be for a feature vector? would it be a vector constituting of all the classes such that 0 value is given to non-relevant classes and 1 for relevant classes? So if the class labels' list is `[Sports, News, Action, English, Japanese]`, then for document `doc_1` the label would be `[1, 0, 0, 1, 0]` ?


title: How to take account for industry and year effects?

text: Controlling for variables like, firm size etc, is quite obvious. But I really have no idea on how to control for industry effects and year effects. I've read that you must insert dummy variables, with two digits SIC of GICS codes. But how does one do that in SPSS (or Stata). Can you please help me with this?


title: Soft-thresholding vs. Lasso penalization

text: I am trying to summarize what I understood so far in penalized multivariate analysis with high-dimensional data sets, and I still struggle through getting a proper definition of _soft-thresholding_ vs. _Lasso_ (or $L_1$) penalization. More precisely, I used sparse PLS regression to analyze 2-block data structure including genomic data (single nucleotide polymorphisms, where we consider the frequency of the minor allele in the range {0,1,2}, considered as a numerical variable) and continuous phenotypes (scores quantifying personality traits or cerebral asymmetry, also treated as continuous variables). The idea was to isolate the most influential predictors (here, the genetic variations on the DNA sequence) to explain inter-individual phenotypic variations. I initially used the mixOmics R package (formerly `integrOmics`) which features penalized PLS regression and regularized CCA. Looking at the R code, we found that the "sparsity" in the predictors is simply induced by selecting the top $k$ variables with highest loadings (in absolute value) on the $i$th component, $i=1,\dots, k$ (the algorithm is iterative and compute variables loadings on $k$ components, deflating the predictors block at each iteration, see Sparse PLS: Variable Selection when Integrating Omics data for an overview). On the contrary, the spls package co-authored by S. Keleş (see Sparse Partial Least Squares Regression for Simultaneous Dimension Reduction and Variable Selection, for a more formal description of the approach undertaken by these authors) implements $L_1$-penalization for variable penalization. It is not obvious to me whether there is a strict "bijection", so to say, between iterative feature selection based on soft-thresholding and $L_1$ regularization. So my question is: **Is there any mathematical connection between the two?** **References**   1. Chun, H. and Kele ̧s, S. (2010), Sparse partial least squares for simultaneous dimension reduction and variable selection. _Journal of the Royal Statistical Society: Series B_ , _72_ , 3–25.   2. Le Cao, K.-A., Rossouw, D., Robert-Granie, C., and Besse, P. (2008), A Sparse PLS for Variable Selection when Integrating Omics Data. _Statistical Applications in Genetics and Molecular Biology_ , _7_ , Article 35.


title: Interpreting output of igraph's fastgreedy.community clustering method

text: With the help of several people in this community I have been wetting my feet in clustering some social network data using igraph's implementation of modularity-based clustering. I am having some trouble interpreting the output of this routine and how to use it to generate lists of members of each community detected. This routine outputs a two-column matrix and a list of modularity values. From the docs: > merges: A matrix with two column, this represents a dendogram and contains > all the merges the algorithm performed. Each line is one merge and it is > given by the ids of the two communities merged. The community ids are > integer numbers starting from zero and the communities between zero and the > number of vertices (N) minus one belong to individual vertices. The first > line of the matrix gives the first merge, this merge creates community N, > the number of vertices, the second merge creates community N+1, etc. > > modularity: A numeric vector containing the modularity value of the > community structure after performing every merge. Working with this explanation and looking at the example at the bottom of the man page, I think the communities in this graph are               1st community: 0 1 2 3 4     2nd community: 10 11 12 13 14     3rd community: 5 6 7 8 9      Can someone who has used this method before confirm whether the right approach produces this result? I have basically (i) ignored the last two merges and (ii) gone over each row in the 'merges' matrix, combining each pair of vertices into a set while watching out for vertex values that are larger than the number of vertices (and therefore refer to another row in the 'merges' matrix). Many thanks! ~l


title: Sequential experiment: controlled design and metrics

text: This is a generic question. Let me put forth an example scenario. Say, I had 2 techniques for allocating my daily budget within 4 stocks. Upon allocation, I get the data about the stock's performance the next day. (Am lazy, and go to a movie after allocation). Problem 1:I have to evaluate or come up with an experimental design to test the two techniques in this sequential experiment over a given number of days. Problem 2 with caveat: If I split my budget into half and have "one" of the two techniques allocate it amongst the four stocks in the same quantities in two parallel experiments at the same time-I find the following issue on the second day. Though everything from the budget to the allocation was the same. Both the experiments gave different performances -the next day- due to the randomness in the system. Under this situation- where performances are different even under "one" technique- How would I evaluate or design an experiment for comparing two techniques over a given number of days? Problem 3: If instead of getting the performances the next day- I get minute by minute or hourly performances, and would like to evaluate the two technqiues-what would be your line of thought?


title: Does the log likelihood become unimodal when the sample size goes to infinity?

text: I know that, under the usual regularity conditions, the MLE converges to the true parameter values as the sample size gets large. And the scaled MLE tends to being normally distributed. However, in a number of real world cases, finding the global maximum is difficult due to the presence of local minima in the likelihood. Intuitively, I would expect that the number and depth of the local minima would decrease as the sample size increases, leading to a unimodal likelihood. But is this true? Has it been proven?


title: Calibrating a household survey to household-level and person-level control totals

text: Imagine a (reasonably large) household survey where all persons in every household have been questioned. For the purpose of microsimulation, this survey needs to be expanded to a full population. In a first step, weights are attached to each observation so that external control totals are obeyed (calibration). If we only have control totals that describe how many households of this-and- that type are in a zone, we can use IPF (also known as raking) which gives a maximum-likelihood estimate of the weights. Minimizing the relative entropy is equivalent to raking/IPF. **EDIT** : But what if we have control totals at person and household level? Like, telling us how many households of which type _and_ how many persons of which sex/age/education level/... there are. I was unable to find a "standard" approach here. Is raking/IPF the "correct" approach from a statistical point of view? Are there other options? What would be, from a statistical point of view, the most reasonable approach to calibrate the weights in the presence of control totals at household and person level? See the original question for more context. (It was probably too big, I'm splitting it into parts.)


title: Fit VAR model with unknown order in Matlab

text: I have a multivariate observed time series $Y_t$ and I want to find the best fitting VAR process for it. I have the econometric toolbox in Matlab and can use 'vgxvarx' if I pre- specify an order for the VAR process. However I want the order of the process to be estimated as well. Is there a function to do this in Matlab?


title: Girsanov Theorem

text: I am interested in understanding Kallenberg's Proof of Girsanov's theorem. The Theorem says: If $Q=Z_{t}P$ and $Q$,$P$ are probability measures and $Z$ is a.s. continuous, then for all local P-martingale $X - X':=X-Z^{-1}[X,Z]$ is a local $Q$-martingale. I do not get the parts which say that $-[X,Z]+[X',Z]=0$ and that $X'Z$ is a $Q$-martingale because of $X'Z-(X'Z)_{0}=X'*M+Z*m$, where $[.,.]$ is the quadratic covariation, and $X'*Z$ is a stochastic integral.


title: Significance of adjacency in correlation matrix with ordered variables

text: The following values are from the Discrete module in BayesTraits. Each value in the matrix is the difference in log-likelihood values between a 4-parameter model of independent evolution and an 8-parameter model of dependent evolution. I apologize if I used the terms "correlation" or "correlation matrix" incorrectly. The documentation is here . The data here is the result of an analysis with two inputs: A phylogeny, or evolutionary tree, and a matrix of presence/absence data (1/0) for 2 traits, for each species in the tree. The matrix looks like this: ![enter image description here](http://i.stack.imgur.com/2YYzW.png) Normally, I could use this to get p-values for correlations between 2 variables. The p-value for the correlation between 4 and 5 would be the p-value of 2 (85.83) in a $\chi^2$ table with four degrees of freedom. Also, I should mention that the order of the variables is very important. Each number (1-9) corresponds to a serially homologous trait on an animal, and the serial homologues appear in the same order in which they are presented here. 1 comes first, then 2, etc. So here's what I want to do: I want to get a p-value **_for the whole matrix_** that tells me whether adjacency significantly impacts correlation values. A p-value below 0.05 would indicate that the correlation values between a trait and the traits to either side of it (e.g., 4 compared to 3 and 4 compared to 5) are higher than would be expected by chance. I see this in many matrices but would like to confirm my suspicions. Is this doable in R?


title: Can someone give a simple guide of Dirichlet process clustering?

text: I devoted my last 3 days to understand Dirichlet process and Dirichlet process clustering but because of my lacking knowledge of stochastic process and measure theory (I did not take any course on those) I cannot get it what Dirichlet process exactly is. At least after 3 days I know what is Dirichlet distribution :) but cannot step forward. If some one give me a naive introduction for Dirichlet process and clustering I will be appreciated. Basically my intuition confused in;   * what is the difference between (Dirichlet) distribution and (Dirichlet) process   * How Poly's url or stick breaking emphasis Dirichlet process,   * Most importantly how Gibbs sampling based clustering of Dirichlet mixture model utilizing Dirichlet process. How you extend my capability on those axes, I am grateful as this much.


title: Programmatically check if two classes are linear separable

text: How to programmatically check if two classes are linear separable? (only yes/no answer) Maybe I can train linear SVM on data(use parameters to overfit) and then run test on same data and if accuracy will be 100% dataset is linear separable.


title: Sum of sample mean and sample variance sampling distribution

text: Let $X_1, X_2, \cdots, X_n$ be an identical and independently distributed sample from $N(\mu, \sigma^2)$, define: $$D = \frac{1}{t}\left[\overline{X} + \frac{1-\rho}{2} S^2\right]$$ where: $t$ and $p$ are fixed constants $\overline{X} = \frac{1}{n}\sum_{i=1}^n X_i$, i.e, the sample mean $S^2 = \frac{1}{n} \sum_{i=1}^n \left(X_i-\overline{X}\right)^2$, i.e, the biased sample variance Then what is the sampling distribution of $D$? * * * I know that $\overline{X} \sim N\left(\mu, \frac{\sigma^2}{n}\right)$ and $n\frac{S^2}{\sigma^2} \sim \chi^2(n-1)$, but how do I derive the sampling distribution of $D$?


title: Is this an application of Frisch-Waugh-Lovell theorem?

text: Firstly, the regression of Y on $$X_1, X_2, X_3, X_4$$ yields: $$Y=\beta_0+\beta_1X_1+\beta_2X_2+\beta_3X_3+\beta_4X_4+e (1)$$ Then, regress Y on $$X_2, X_3, X_4$$ $$Y=\beta'_0+\beta'_2X_2+\beta'_3X_3+\beta'_4X_4+e'(2) $$ Lastly, regress $$X_1$$ on $$X_2, X_3, X_4$$ $$X_1=\beta''_0+\beta''_2X_2+\beta''_3X_3+\beta''_4X_4+e'' (3)$$ Now I regress e' on e'' $$e'= \alpha_1+\alpha_2e''+\epsilon$$ Is there a simple method to prove that $$\alpha_2=\beta_2$$ Can I plug (3) into (1) (which yields equation (4) ) and then identity the residuals of (2) and (4) ? p/s: I'm sorry but can anyone tell me how to add comment? I try the button "add comment" but I still can't. I'm new here.


title: Determining the correct sample size

text: Which is the best formula to use when determining a the sample size in a customer and employee satisfaction surveys?


title: How to get the values of the probability density of a given data using R? The goal is to get the values of the function not the graph

text: The problem is the following: If I have data from 2002 to 2005 of operational risk then I want to get the value of density associated with this.               X # is the data from 2002 to 2005, X takes values [-200,200]      density(X) # the density of X       The question is how to get the values of this density? For example, if I take any value X0 in X. say 109, then how to get this value?


title: Significance: p-value or t-value?

text: I'm writing an economics thesis and I've found many econometric papers in which I found regression tables without the p-value. They only show the coefficient of the independent variable and its t-value, like in the example image. ![enter image description here](http://i.stack.imgur.com/Jgc8u.png) My econometrics knowledge is very elementary, but I thought that the most important thing to evaluate the significance of a variable was the p-value. Searching on the web, I've found that sometimes $t-value>2$ is considered as a good value for significance. Is that right? Can I say that when $t-value>2$ the coefficient is statistically significant? Edit: after Rob's suggestion, here's a link to the paper


title: Survival analysis cycle time dependent variable

text: I've a survival analysis study going on and I'm having trouble with dealing with time-dependent variables. To deal with time-dependent variables I have to create intervals when the variable changes. One of the variables is "season" wich has values of breeding, post-breeding or winter. Since all the observations start at a different time (different days and different years) it's not easy to o the math for the intervals because there are three changes per year. The final table would be huge, since I have 154 observations and some of them have 5 years of following... I'm using the R survival package. Can anyone help me to figure an easy way to deal with this?


title: Quadratic models with R. The use of poly(..) and I(..) functions (R-language)

text: What causes the different results below?               var1 = c(0.04875,0.13725,0.28350,0.50975,0.77425,0.94700,0.05325,0.14050,0.29725,0.51525,0.79000,0.95400,0.04625,0.15250,0.29000,0.53300,0.79825,0.95225,0.05025,0.14625,0.28800,0.52625,0.78200,0.95925,0.04700,0.14225,0.30325,0.53500,0.79325,0.95875,0.04775,0.13850,0.28675,0.54250,0.78300,0.95175,0.05150,0.12725,0.30175,0.54725,0.79475,0.96275,0.05375,0.14100,0.30050,0.53275,0.78100,0.96175,0.05450,0.15300,0.29650,0.52850,0.80100,0.95675,0.05425,0.13975,0.30875,0.56025,0.80575,0.96100,0.05100,0.15350,0.31175,0.53300,0.78900,0.96000,0.04650,0.13525,0.29600,0.53625,0.78475,0.96375,0.05375,0.13900,0.29600,0.53725,0.78700,0.95800,0.05075,0.14350,0.29225,0.54525,0.80275,0.95800,0.05050,0.13200,0.29850,0.52700,0.80525,0.96150,0.05150,0.14050,0.29450,0.54375,0.79450,0.96375,0.05375,0.13525,0.30475,0.55250,0.79425,0.96025,0.04950,0.14500,0.29425,0.52250,0.78475,0.95650,0.05225,0.14425,0.29225,0.53150,0.80425,0.95375)          var2 = c(1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6)          var3 = c(4,4,4,4,4,4,6,6,6,6,6,6,8,8,8,8,8,8,10,10,10,10,10,10,12,12,12,12,12,12,14,14,14,14,14,14,16,16,16,16,16,16,18,18,18,18,18,18,20,20,20,20,20,20,22,22,22,22,22,22,24,24,24,24,24,24,26,26,26,26,26,26,28,28,28,28,28,28,30,30,30,30,30,30,32,32,32,32,32,32,34,34,34,34,34,34,36,36,36,36,36,36,38,38,38,38,38,38,40,40,40,40,40,40)               summary(aov(var1~as.factor(var2)*var3))     summary(aov(var1~as.factor(poly(var2,1))*poly(var3,1)))      I might not understand well how poly(..) and I(..) work. Could you give me some hints about these two functions? (I aim to create models of different degree (quadratic, cubic, etc..) and compare their BIC (or AIC)).


title: Testing slopes in multivariate adaptive regression splines (MARS/earth)

text: I am using the `earth` package in `R` to estimate the number of breakpoints in a curve. There is only a single predictor. I was hoping there was a sensible way to compare a slope for some segment against 0. Normally I would use the estimate and standard error, but it's not obvious (to me) what to do for segments that are combinations of multiple slopes:                       estimate    SE     (intercept) 2.78        .015     h(age-4)    -.01        .0029     h(4-age)    -.35        .015     h(age-11)   -.076       .0055      So the slope after age 11 is -.086. How do I quantify my confidence that this slope really is negative?


title: Evaluating predictive models

text: **I am looking for ways of evaluating the performance/success of predictive (classification) models for economic purposes.** I know of:   * Direct accuracy percentage   * AUC score   * Net profit   * Rate of return   * Sharpe ratio **Could you inform me of other possible ways of evaluating the performance/success of such models?** Thank you. PS: In case it's relevant to your answer: there's 2 possible categories the model will classify data objects in.


title: References for Survival Analysis

text: I am looking for a good book/tutorial to learn the survival analysis. I am also interested in references on doing survival analysis in R


title: What test(s) for my data?

text: I'm trying to analyse data for my psych undergrad dissertation but struggling with which tests to use because my project is unlike anything I've had experience with previously (and also I don't have a mind for statistics). I have 200 cases of patients who've attended therapy, and the dates of every session they attended. I made extra variables for the number of days wait between each session but I want to know whether different diagnoses (e.g. depression, anxiety...) had an effect on how often patients came for therapy. The trouble is that the diagnostic groups have very different n in each (there are 5 main diagnostic groups, ranging from n = 109 to n = 3!). I'm not sure if there's anyway to do this? Any advice much appreciated. Thanks!


title: Typical Dimensions of a Doc Term Matrix

text: I'm working on a text classification task with about 4.5 million short documents (abt 100-150 words per doc). After performing the standard preprocessing steps (stemming, stripping punctuation and numbers, lower casing, etc), I created a doc term matrix with approx. 500,000 dimensions. When creating the doc term matrix, I ignore terms with term frequencies less than 5 and terms that appear in more than 50% of the documents. I have a feeling that the number of dimensions I have is way too big and I'm doing something wrong - can someone with experience in the area of text classification tell me if the number of dimensions I have is reasonable? (Note: in the next step of my analysis, I am using LSA to compress things down, but I want to amke sure each step is done correctly).


title: Estimating individual-level outcome probabilities given from group-level training data

text: Given multi-level data, I'm interested in estimating member-level outcomes when only group-level outcomes, along with both group- and member-level predictors, are known. I'm lacking in statistical terminology, so I'll explain with an example: We'd like to estimate the likelihood that a given high school student will be granted college admission. In this case, "groups" will be schools and "members" will be students. We have:   1. Aggregated, group-level data on the dependent variable (total number of accepted/rejected applications from each high school)    2. Aggregated, group-level independent variables for each school (number of students who have taken at least three advanced/honors courses, distribution of GPA, number of students in the top 10% of standardized test scores, etc.)   3. The same independent variables at the member/individual level (for each student, whether he/she has taken at least three advanced/honors courses, actual GPA, whether they scored in the top 10% on a standardized test, etc.) What we're lacking, however, is member-level outcome data on college admission (for each student, whether he/she was admitted to college). Intuitively, I get the sense that some estimation may be made for group members even in absence of this individual-level outcome data: we may see, for instance, that schools with high proportions of students taking at least three advanced/honors courses have correspondingly high acceptance rates. It makes sense that, knowing this, a student's likelihood of being granted college admission should be positively influenced by whether he/she has taken at least three advanced/honors courses (possibly along with other, group-level effects). Is it possible, then, to model group-level outcomes and apply those models to individual cases to estimate individual-level probabilities? If so, what terminology describes this type of model? Any pointers to literature on this topic would be appreciated as well.


title: Fitting series to a distribution

text: I am reading this note on Mean-CVAR optimization. The authors argue against the prevailing assumption that asset classes are normally distributed and propose using a truncated levy flight distribution to describe the behavior of asset classes. Now, here is my comprehension of fitting a historical time series to a new distribution, please kindly correct me if my thinking is wrong.   1. You determine what parameters are needed for your new distribution: In the case of the levy flight distribution, you require the tail index, a skewness parameter, a scale parameter and a location parameter   2. You estimate those parameters from your historical time series: In this case, the tail index is estimated as Kurtosis, the skewness parameter is estimated as historical skewness, the scale parameter as standard deviation and the location parameter as the historical average.   3. You simulate your new variables by applying your estimates to the inverse of the cumulative distribution function of the levy flight distribution. **Can someone explain how this is done?** More questions:   * Once you fit a new distribution, does that imply that your old estimates are no longer relevant?   * What are some good ways to test for the goodness of fit for a levy flight distribution?


title: Difference in chi square for factor analysis

text: A couple of questions relating to factor analysis I am conducting. Parallel analysis suggests 5 factors (or possibly 6) for my model and both make clinical sense. However, the seventh factor only has 2 items - it is interesting that they came out together as they are slightly more extreme items and it does make sense for them to be together. They cross-load in the 5 factor model. There is very little difference between the 5 and 6 factor models in fit statistics (TLI,CFI RMSEA). However the difference in chi-square is highly significant (90 on 30 degrees of freedom). I am inclined to stick with the 5 factor model, removing the two items (but mentioning them in the report) but I am unsure how to justify this due to the significant chi-square. While the chi-square test for the model is strongly dependent on sample size is the difference in chi-square? If not, how to justify ignoring it? Do you think parallel analysis to select the number of factors takes precedence over the chi-square difference? *sneaky edit


title: How can I measure which function is a better fit to a set of datapoints

text: The problem is from here. Namely, I have set of data points               data = [[90.00, 2.0], [97.40, 5.0], [104.8, 14.0], [112.2, 12.0], [119.6, 11.0], [127.0, 6.0], [134.4, 3.0], [141.8, 1.0], [149.2, 2.0], [156.6, 1.0]]      I have to fit a curve of the form $$\frac{m}{\sigma\sqrt{2\pi}}\cdot e^{-(x-\mu)^2/(2\sigma^2)}$$ to this data. In the link, Volker Brown got the values $$m = 405.75796954829985, \mu = 111.86913960269014, \sigma = 11.968861052746961$$ and as I tried my version of Sage to the same to that set of data, I got $$m = 405.7572223708457, \mu = 111.86914614035226, \sigma = 11.96881745593664$$ What kind of ways there are to measure which one is better fit to that set of data? I heard that sometimes one can use log-likelihoods to select the best model and sometimes Akaike's information criterion is suitable for that. But which one should I choose in this situation? Thusfar, I haven't studied the theory behind the AIC. Another thing that confuses me is Brown's comment "You need to give some hint for the initial values." Does that mean that there is no standard way to measure how good a particular curve fits to the data? Or is that just a bug in a previous version of Sage that has been fixed in the version 5.6, 21 January 2013?


Is the skew normal distribution log-concave?


title: Multinomial logistic regression?

text: How can I look for collinearity diagnostics, particularly condition indexes and proportions of variance in Multinomial logistic regression?


title: Test for IID sampling

text: How would you test or check that sampling is IID (Independent and Identically Distributed)? Note that I do not mean Gaussian and Identically Distributed, just IID. And idea that comes to my mind is to repeatedly split the sample in two sub- samples of equal size, perform the Kolmogorov-Smirnov test and check that the distribution of the p-values is uniform. Any comment on that approach, and any suggestion is welcome. **Clarification after starting bounty:** I am looking for a general test that can be applied to non time series data.


title: What are good references to jump from discrete to continuous time series models?

text: I've been studying discrete time series models for a while. I'm familiar with the usual time series models (ARMA, ARIMA, GARCH, etc.), but I want to study their continuous counterparts (CAR and CARMA processes seem very interesting) in order to implement some models to unevely spaced time series. Despite of the above, I've been trying to find some introductions to continuous time series theory and methods as good as Brockwell & Davis (1991) book, but I couldn't find any, just some reviews here and there (specially the ones by Brockwell). Do you have any suggestions for a good book on this topic? Thanks in advance! PS: Please take into consideration that my background doesn't have much measure theory, and in order to read the discrete Time Series literature I had to pretty much teach myself from lecture notes and reviews I found on the web (suggestions on this topic are also welcomed!).


title: On the corrections for multiple comparisons

text: I've used the Mann-Whitney U-test in order to check for possible differences between two relatively small groups (12 subject in each) on 15 different items. I didn't have any pre-defined assumption about which items would be different. However this simple Mann-Whitney U-test yielded a significant different (P<0.05) for **ALL** items tested. **My question** : if/what corrections for multiple comparisons should I use here? Does the fact that I have found a significant different in all items make any difference?


title: Proofs of the central limit theroem

text: I know there are different versions of the central limit theorem and consequently there are different proofs of it. The one I am most familiar with is in the context of a sequence of identically distributed random variables, and the proof is based on an integral transform (eg. characteristic function, moment generating function), followed by first order approximations to obtain a function to which the inverse transform can be applied. I am interested to know if there are any flaws in this approach - I have been told informally that it is not completely rigorous - but why ?


title: What is the distribution of the ratio of sums of squared normal random variables?

text: Given that $Y_1, \ldots, Y_6$ are independent normal random variables with mean $0$ and variance $\sigma^2$. Find the distribution of the following $ T=\displaystyle\frac{Y^2_{1}+Y^2_{2}+Y^2_{3}}{Y^2_{4}+Y^2_{5}+Y^2_{6}} $


title: Archetype path clustering

text: I have a sample of customers who took various paths to an end state. Between the start and end, they may have completed various interactions in between. Some of these interactions are dependent on previous actions; others, don't. Thus, I have a number of trinary categorical variables representing a finite number of interaction types (didn't do a thing, did a thing prior to another event, did a thing after another event); and I have time stamps associated with most of the interactions. Moreover, I have other demographic attributes that can represent confounding factors. I would like to reduce the sample into a smaller number of path archetypes. To me, this reads a little like each of path analysis, time-to-event analysis, latent class analysis, and multiple correspondence analysis. I have some experience in clustering, but not when there is interdependence or semi-strong correlation between variables. And certainly not when there is a time element involved. Any suggestions on how to approach? I figure it'll be an evolutionary process where I may start out with a less powerful technique, and perhaps select certain archetypes for more powerful deeper analysis. I would like to conduct, preferably, in R or SAS. Thanks for any advice.


title: Determing whether a moderation effect exists

text: I am currently trying to interpret the results of a regression analysis where I am looking for a moderation effect. While the interaction term is significant, the main effect of the potential moderator with the independent variable is not significant. Is a moderation effect only dependent on the interaction effect, or does the main effect of the potential moderator also have to be significant with the same direction?


title: Confidence interval for Bernoulli sampling

text: I have a random sample of Bernoulli random variables $X_1 ... X_N$, where $X_i$ are i.i.d. r.v. and $P(X_i = 1) = p$, and $p$ is an unknown parameter. Obviously, one can find an estimate for $p$: $\hat{p}:=(X_1+\dots+X_N)/N$. My question is how can I build a confidence interval for $p$?


title: Confidence Score for a Single Proportion

text: I was asked to calculate a confidence score for a single proportion representing a conversion rate. I saw a lot of documents that calculate confidence _intervals_ for proportions based on the binomial distribution; however I did not see any material for a single conversion (or proportion). If $X \sim Bin(n,p)$, with large sample size normal approx to binomial is valid, then $$ \hat{p} \leadsto \mathcal{N}\bigg( p, \sqrt{\frac{p(1-p)}{n}} \bigg) $$ Using this approximation, I saw statistical significance involving a second ratio (through $z$-scores) can be calculated, etc. However, can any error metric be calculated for a single proportion? Ideally this metric would be a number between 0% and 100%.


title: Relevance of opinion poll (500 people from 10.5 millions)

text: In local news [1] comes with information that 58% of people is for something. This information is based on option pool that ask 500 people from 10.5 million about that question. What is confidence interval for people that are for that something? To solve this I try following. The $n$ people from $10^7$ is for something. I randomly select 500 of them and I will try to find $n_1,n_2$ for which: $$ n_1: \frac{\binom{n_1}{290}\binom{10^7-n_1}{210}}{\binom{10^7}{500}}<0.05 $$ $$ n_2: \frac{\binom{n_2}{290}\binom{10^7-n_2}{210}}{\binom{10^7}{500}}>0.95 $$ But this is not easy to solve. So I am looking for some solvable algorithm. [1] link


title: Percent error for linear regression model

text: Suppose I fit a linear regression $y = \beta x + \rm error$. In this situation, $x > 0$, $\alpha > 0$, and therefore $y > 0$. Moreover, the $\rm error$ is normally distributed with mean $0$ and standard deviation $\sigma$. Now, if I calculate another error metric, percent error, which is $100*\rm predicted / actual$, then the error decreases with the value of $y$. That is, the model looks better with increasing values of $y$. This seem paradoxical. Any thoughts? An example is suppose we are predicting sales as a function of amount of discount. Sales increases with the amount of discount. The model, on percent error basis, appears to perform well for higher amount of discount (higher sales). Any insights are appreciated.


title: what's the pdf and covariance for this distribution?

text: I am stuck on a problem and wonder if anyone can give me some suggestions. $X_1, X_2, X_3$ all follow a $\text{Uniform}[0,1]$ distribution and are subject to the constraint $X_1+X_2+X_3\leq 1$. What's the joint distribution for $(X_1, X_2, X_3)$, that is, what's $p(X_1, X_2, X_3)$, and what's the variance-covariance matrix for it? I get the joint distribution by geometrical way, that the pdf should be $1/6$. However, I can't calculate the variance-covariance matrix for it. I wonder how to get it?


title: Bayesian inference and curve fitting

text: I'm a silly scientist, and I've got a process I'm interested in. It gives me lots of time-series data. I've got an explicit, analytical model that is firmly rooted in reality, and basically boils down to fitting exponentials. Unfortunately, the model is under-determined with the experiments I can do; two parameters often end up highly correlated. As ever, it's the parameters I really care about. At the moment, I'm using matlab's `lsqcurvefit` to fit my model to data. I estimate appropriate starting values from other means, and fit parts of my data to different models, then use those extracted values on other (later) data to get other parameters. It does a cracking job, and usually converges on sensible values pretty quickly. I can construct appropriate _a priori_ estimates for the likely distribution of all these parameters -- it's all just chemical kinetics, and I can go and measure some of these parameters in analogous environments _in vitro_. I'd really love to do something more intelligent than just using `lsqcurvefit` like a muppet. Am I right in thinking that there are Bayseian methods to tackle this nut? What should I do? Is this a solved problem in bayseian inference? If you could possibly point me to an appropriate introduction, I'd really love to read it. I'm a physicist working in a biological situation, if that helps you to think of any texts.


title: Interpretation of Kolmogorov-Smirnov Critical Value Generated Distributions

text: As a non-statistician, I need help in interpreting a customer specified two- part reliability requirement that I think involves KS. Requirement Part 1 R[4 years] must be greater than or equal to 0.95 and R[8 years] must be greater than or equal to 0.85 I have plotted the reliability (survival) function of a 2-parameter Weibull distribution that meets the above requirement in Plot A below. The shape parameter is 1.664 and the characteristic life is 23.844 for this distribution. ![Plot A\]!\[enter image description here](http://i.stack.imgur.com/bY5JN.png) Requirement Part 2 The confidence level shall be 90% when demonstrating the Part 1 requirement via product life testing. It’s the Part 2 that I’m a bit shaky on. On page 8-54 of MIL-HDBK-338B (http://www.sre.org/pubs/Mil-Hdbk-338B.pdf) there is a table showing KS critical “d” values as a function samples size, N and significance level, alpha (also note the plot on page 8-57). From this table I took a d value of 0.264 based on a signficance value of 0.10 and a sample size of 20. Plot B below shows my result. My interpretation of Plot B is that after running a life test on 20 samples that if the resulting reliability plot does not fall below the lower boundary shown in Plot B then we have met the requirements. ![Plot B\]!\[enter image description here](http://i.stack.imgur.com/ilZqf.png) I have two questions:   1. Did I translate the Part 2 requirement properly when I used an alpha of 0.10 to obtain the KS critical value of 0.264 ? In other words, does a 90% confidence equal a 0.10 significance within the KS context ? If not, can someone provide guidance ?   2. How would _you_ interpret Plot B ? Many thanks. * * * Response to owensmartin's answer Thank you for your reply. I'll take your "your reasoning is not incorrect" statement as a big vote of confidence. I just have a few items below that I'd appreciate comments from anyone on. a. You are right in that my reliability functions are also survival functions. I believe the Nominal curve is simply the compliment of the nominal CDF. b. Although I haven't computed the statistical power, I'm not that surprised when you say it is "very low at this sample size". However our customer is comfortable with the sample size of 20 and by extension the resulting low power. c. Our customer is also comfortable with using the referenced Military Handbook and the associated KS critical values shown therein even though they may be approximate. As time permits I'm hoping to learn how to compute these critical values "exactly" so as not to have to rely on the handbook values. d. To answer your question as to whether I really need the Weibull fit ? The short answer is no as there is nothing particularly "magical" about the Weibull distribution. The slightly longer answer is that I'm not sure how else to produce a nominal survival curve that I can then apply the d value of 0.264 to. The nominal curve shown in Plot B was constructed by solving the two simultaneous equations for the shape parameter and characteristic life knowing that R[4] =0.95 and R[8]=0.85 meet the requirements. This is shown in Plot A. I suspect that other commonly used distribution types may not be able to meet both of these requirements simultaneously. But because of the inherent flexibility of the Weibull, it _is_ able. e. As for being sure about the 90% confidence, I think as long as the fitted distribution resulting from the 20 sample life test does not drop below the lower boundary shown in Plot B would you agree that we can say we are 90% confident we meet the stated requirements? The only exception I can see to this is if when fitting the test data to a distribution that if the degree of fit is so poor as to raise question to its validity. But in solving engineering problems, this is always present as a concern that needs to be dealt with. f. Regarding failed vs didn't fail test results, we normally strive for having each of the 20 samples fail so that we can avoid dealing with the added uncertainty associated with censored or suspended data points. We achieve this via accelerated-life testing methods which essentially compresses time by increasing either the duty cycle and/or the stress value(s). Thanks for any further insight into this.


title: Observations with very high or very low residuals in regression

text: If a regression model is applied and there exist residuals that is very high or very low (meaning outliers compared to the others), is it good practice to get rid of those observations and then do the regression again, particularly if you have a very big sample of data?


title: My data lies on a linear plane

text: My purpose is to classify 3 classes from an EEG data. When I plotted my data on feature space in order to visualize, I found they lie on such a linear plane (please see my figures). Before plotting, I had applied several methods for feature extraction (ex: Common Spatial Pattern, ...) I will explain my procedures here briefly. I used Common Spatial Pattern. My data are EEG signals consisting of 2 channels. I apply pass-band filter to extract 8 more sub-bands for each original channel. The sample rate is 128Hz, and I chose Windows size is 1 second, so my data at 1 second can be denoted as $\\{E_n\\}_{i=1}^M \in R^{\rm{channel} \times 128}$ where $M$ is a number of trial. Each trial will be associated with a specific class. I'm trying to find $W$ such that $W^\intercal E = S$ and $S_n \in R^{d \times 128}$ ($d<\rm{channel}$). Each trial's feature is $(x_1 ... x_d)$ where $x_i = \log \bigl(\frac{var(s_i)}{\sum{s_i}}\bigr)$. I chose $d=2$ for figure 1 and $d=3$ for figure 2 and plotted $S$. My question is: what is the phenomenon? Is the problem in my data or feature extraction? ![I projected onto 2D plane](http://i.stack.imgur.com/Q1yFz.jpg) ![I projected onto 3D plane, each color correspond with one class](http://i.stack.imgur.com/3KMxJ.jpg)


title: problems using arima model

text: Forgive me for my lack of knowledge. but i am having troubles using arima model in e views. I don't know how to include independent(control) variable in my regression expression. th model i am using is ls Dy c beta ar(1) ar(2) ma(1) ma(2) Dx1 Dx2 Dx3 can anyone guide if this model is correct or not. i have searched many websites and papers but couldn't find how to regress with control variable. thank you in anticipation


title: Deriving the optimum value of a function

text: I have a function $f(t) = \sum_{i=1}^{N} |y_i-t|$. What will be the optimum value of t that will minimize it. How to derive it? Similarly what is the optimal value of t which minimizes $f(t) = \sum_{i=1}^{N} |y_i-t|^{\infty}$?


title: Should a Poisson regression be carried out with only 3 data points?

text: I'm trying to test the relationship between the number of adults counted and the percentage heather cover over 3 areas. The data looks like this:               Area      Cover    Count     Field1     53.7      216     Field2     19.5        2     Field3     39.8     6106      Given that count data is involved I was going to carry out a Poisson regression. I know a relationship is unlikely given this data but my question is, is 3 data points enough for a regression analysis? This is important as I also want to test the relationship between adult counts and nest number, adult counts and vegetation height etc. but all of which only have 3 data points as I was only studying 3 areas on a site. With this type of data, I don't know if any other test could be used given that it is count data. I'm still a beginner when it comes to stats so any help would be much appreciated. I'm using R to carry out all my analyses.


title: Who is the best writer among statisticians?

text: I'd like to learn their writing style. Could you recommend a good writer statistician?


title: How Does a Disparity in Number of Documents (Training Data Points) Affect Text Classification?

text: I have collected a fairly clean set of data (5,410 documents) to train a text classifier. I am now attempting to improve my classification success. (Note: When I trained/tested the classifier from in-sample data, I scored really well. Now that I'm using out-of-sample data, my success is down about 25%. Bad news, but not unexpected for a total ML newb.) Does having an equal (or more equal) distribution of training documents improve classification success? For instance, although I have only six categories -- Labels A,B,C,D,E, F -- the number of training docs I have for each label varies widely, as Label A = 328; Label B = 1,973; Label C = 148; Label D = 822; Label E = 40; and Label F = 1,242. Any thoughts?? (Note: I don't actually know if this is a representative distribution of what I'm likely to to see out of sample. Is there an ideal algorithm to help with this uncertainty?) Presently I am using Scikit's `svm.LinearSVC()` to classify. Thank you for any insight!


title: Beta Vector from Scaled PLS Model not Reproducing Predicted Responses

text: I am running R 3.0.2 and using the PLS package to build partial least squares models. The problem I am having is that when I apply feature-scaling to my design matrix the resulting coefficient vector does not reproduce the predicted response values on my training and cross-validation samples. Here is my code:               model.train <- cppls(y ~ ., data = data.train, ncomp = 20, scale = TRUE,                           validation = "none", weights = w)     yhat.train <- fitted(model.train)     B <- coef(model.train, intercept = TRUE, ncomp = 20)     yhat.cv <- predict(model.train, newdata = data.cv, ncomp = 20, type = "response")      If I then take the B vector and multiply by my design matrix, the resulting values do not equal the values generated by the `fitted()` and `predict()` functions; in other words, `Y != BX`. However if I run the exact same code with `scale = FALSE` in the `cppls()` function, `Y = BX`. My question is whether the `coef()` function is simply failing to re-scale the betas, or if I'm making some error in my routine. I would really like to be able to apply feature-scaling to my design matrix, but I need the correct beta vector as an output. Would greatly appreciate any advice.


title: Basic GARCH (1,1) question

text: Background to question: Hi, I was trying to fit a GARCH(1,1) model to the variance of log returns of a series, and ARMA(0,0) for the mean. I was using the `fGarch` package in R to do this. The aim of the modeling is to generate a predicted volatility number to feed into the Black-Scholes model to an generate option price and therefore option deltas. I plan to backtest the delta from GARCH volatility to hedge my option positions (as opposed to deltas derived from implied vol prices). Questions: A) I used the `predict` function in the package to generate a 'n-day ahead' volatility forecast. As I understand GARCH, these numbers are annualized standard deviation numbers. To hedge a 1 month option I want to forecast 30 day volatility. I can simply put 'n-days ahead = 30 to get the numbers, but **how do I combine those 30 numbers to get a annualized vol number?** B) Could anyone also please explain how to use the `nroll` argument in the package? Basically I want rolling GARCH estimates of volatility. For example, at day 10, I want to use the past 10 days of data to get a vol prediction for day 11, at day 50 I want to use 50 days of data for vol prediction of day 51 etc. Any help would be greatly appreciated.


title: ANOVA using R - unsure of if analysis is appropriate, and if variables need to be numeric

text: I'm doing an ANOVA to compare the results of two trials. Individuals had to walk towards a bucket blindfolded, with their deviation from the target recorded. I want to see if individuals get better with practice (are trial-2 deviations lower than trial-1) and is there a relationship between the trials regarding gender and handedness (left handed/right handed). **Is an ANOVA appropriate?** Also, I'm having trouble in R because I've put gender as M and F, and Handedness as L or R. Do I need to make these numerical or not?


title: Sales forecasting

text: I have weekly sales data for the previous year of short-term foods. I want to forecast sales for the next 30 days. I am able to forecast sales figure for a week. But the problem I am facing is how can I generate a daily forecast? Is there any method for doing this or I have to distribute a weekly forecast for all seven days?


title: Calculating primary groups

text: I am preparing to analyze some data and I have a question about how to go about it. It is a 39 question survey with 5-item Likert scale responses. I will be using the responses to separate the participants into the category that they score highest in for further analysis. The problem is that Category A corresponds to 14 questions, Category B has 11 questions, and Category C has 14 questions. What is the best way to go about determining which one the participant scores highest in, since the # of questions is not equal? Will it work to take their total possible points in each category and divide it by the number of points possible, then compare their percentages across the three categories to see which is the highest?


title: Confusion related to kruskal wallis test

text: I have some confusion related to Kruskal wallis test. I have an example lets say               X=[2 2 35 10 9 8 11 12];     Y=[1 1 1 2 2 2 2 2];      Y is the group variable Now when I ran the kruskalwallis test               p = kruskalwallis(X,Y,'off')      I got p values of around 0.4. I was assuming the Kruskal wallis test takes the median. So it should have been robust when I added an outlier with value 35 in the third position. Why isn't it robust to that. Is it because I have very few samples. Can anyone explain?


title: Augmented Dickey Fuller output conflicting in Stata

text: I am required to perform unit root testing on a given time series. The output obtained in Stata is somewhat confusing me. To the best of my knowledge I am obtaining two conflicting results, Stata indicating that the time series fits a unit root process while also seemingly saying that the coefficient is significantly different from zero, hence contradicting the idea that the time series indeed follows a unit root process. This is the output obtained in Stata:               . dfuller series3, noconstant regress lags(0)          Dickey-Fuller test for unit root                   Number of obs   =        48                                     ---------- Interpolated Dickey-Fuller ---------                   Test         1% Critical       5% Critical      10% Critical                 Statistic         Value              Value            Value          ------------------------------------------------------------------------------      Z(t)              8.943            -2.623            -1.950            -1.609     ------------------------------------------------------------------------------      D.series3 |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]     -------------+----------------------------------------------------------------      series3 |          L1. |   .1314633   .0147005     8.94   0.000     .1018897    .1610369     -----------------------------------------------------------------------------      Clearly the test statistic is indicating that the null hypothesis is accepted and hence the time series fits a unit root process. However, the p-value obtained for L1 is 0.000 indicating that it is significantly different from zero (0.1314633 cannot be treated as being 0). I read somewhere that this might be related to the Augmented Dickey Fuller test itself; however, I cannot seem to find this link again. I also ran the Augmented Dickey Fuller test using lags as well as the Philips-Perron test and obtained similar results. Would anyone have any idea why the ADF is giving such results and what might be the cause of this anomaly?


title: Support vector machine for text classification

text: I am currently having a data set, class 1 with about 8000 short text files and class 2 with about 3000 short text files. I applied LibSVM and tried a couple of parameter combinations in the cross-validation experiment. Generally the class 1 precision falls into the range of (85%, 90%); the class 2 precision falls into the range of (70% , 75%); the recall of both class 1 and class 2 fall into the range of (80% , 85%). For the text classification purposes, I built text feature space following the common approaches, tokening the document, filtering the stopwords and building the word vector using tf-idf or binary frequency, etc. I also tried n-gram model to build the feature space. But these approaches did not improve the performance a lot. I would like to know are there any other ways that may help tune the LibSVM to improve the performance. LibSVM provides grid search for parameter setting up, but it runs pretty slow.


title: How to fit a model for a censored binary outcomes?

text: Suppose I have a data frame such that:               set.seed(2014)     df<-data.frame(y=rbinom(100,1,0.3),futime=as.integer(rnorm(100,100,10)),          age=rnorm(100,50,5),sex=rbinom(100,1,0.5))     df[df>100]<-100       y futime   age sex     1 0    100 49.72   1     2 0    100 53.63   0     3 0     91 55.03   1     4 0    100 50.58   1     5 0     98 51.32   1     6 0    100 54.91   1      where `y` is the outcomes and _futime < 100_ are right censored of follow up time. Difference in survival analysis where we model the survival time given their censoring status, here I'm interested in whether a subject may have experienced an event or not given its censoring status or followup time. I am wondering which is the best model for this kind of situation. Or, is it enough to just use the `futime` as a regular predictor in a normal logistic model? I'd think a special situation should be accounted here. Would somebody have some suggestions?


title: Plotting and presenting longitudinal data, options?

text: **Hypothetical Scenario** : A few continuous variables, each measured repeatedly at, say 12 time points, each with say, 150 observations. There are small fluctuations from one time point to another (i.e. change is not that smooth), but the variables are showing overall time-related changes. **Question** : What are different options and their pros and cons in plotting and presenting change in such data? Would be excellent if example graphs/codes could be attached. This is probably a very vague and general question, and I understand that it very much depends on what one wants to highlight, the kind of analyses conducted, the number of time points, the number of observations, etc. But it would be very helpful to have an idea how people usually plot their longitudinal data and what (nice) options are available. I'm particularly interested in hearing the options to highlight time-related changes integrated with individual differences (e.g. variability at each time point and variability in growth curve trajectories) in these changes. For those interested, here is some randomly selected 100 observations from one segment of 12 time points to play with. Missings were coded as NA and data were structured wide.


title: Multiple linear regression

text: I am new to regression analysis. Therefore, I would be grateful if you could kindly help me to set a model for my dataset. Briefly, I have images of diseased eyes (10) and control eye images (25). I generated 167 statistical features of all these eyes to find the most discriminative features or the degrees to which I can discriminate all of the features. I also have the thickness layer where there is an assumption that illness occurs as thickness decreases. I figured out that I need to do a multiple regression analysis: using the thickness values of the eye as dependent variable and 167 features as independent variable. I scaled all the values between 0 to 1. I don't know how to start with in Matlab. I used the `regress` command but the answer was not quite right to me. Could you please advise me if the model is right and what to do next, how to present the results.


title: Using log-linear models for presence/absence data in wildlife

text: I'm working on a project wherein I compare the presence/absence of a number of bird and herptile species between wetlands that have received three different treatments. The populations were surveyed across two different years. So the response variable is a binary categorical variable, and the predictor variables (wetland treatment and year of sampling) are also categorical. The situation gets a bit more complicated in that different sampling schedules were used between the two years, resulting in different sample sizes between the two years. I'm looking at modeling animal presence/absence using the glm() function in R, but I'm not sure if there might be a more appropriate approach?


title: Significance test for large sample sizes

text: This is more of a theoretical question. Super large sample sizes will almost always show a significance when a $\chi^2$ test is done. Is there any other statistical test of significance (an alternative to $\chi^2$) which is good for testing independence when the sample size is very large? This is the context of my problem: I have 2 large datasets of phrases. Set1 corresponds to the Google n-grams set and set2 is a smaller set corresponding to the phrases found in one single website. Now consider a phrase: say ('Technology') found in Set2. I want to test and see if this phrase is specific to this website (it could be if it is a Technology website) or if it is a general phrase. So I am performing a $\chi^2$ test between the frequency of phrases in the two sets as follows:                                       Set1                        Set2     Not_Technology   (set1) 2,674,797,869,255    (set2) 46,168,477.00      Technology       (set1) 1710231              (set2) 1991      I understand that this might not be the best method to test whether a phrase is a general phrase or not, so if you have any suggestions or criticisms I am happy to listen to them.


title: Is multicollinearity implicit in categorical variables?

text: I noticed while tinkering with a multivariate regression model there was a small but noticeable multicollinearity effect, as measured by variance inflation factors, _within_ the categories of a categorical variable (after excluding the reference category, of course). For example, say we have a dataset with continuous variable y and one nominal categorical variable x which has k possible mutually exclusive values. We code those $k$ possible values as 0/1 dummy variables $x_1, x_2,\dots ,x_k$. Then we run a regression model $y = b_0 + b_1x_1 + b_2x_2 + \dots + b_{k-1}x_{k-1}$. The VIF scores for the $k-1$ dummy variables turn out to be non-zero. In fact, as the number of categories increases, the VIFs increase. Centering the dummy variables doesn't appear to change the VIFs. The intuitive explanation seems to be that the mutually exclusive condition of the categories within the categorical variable causes this slight multicollinearity. Is this a trivial finding or is it an issue to consider when building regression models with categorical variables?


title: Are integer results from random number generators unlikely?

text: If I generate a random float value between 0…1, say to 40 digits, or _n_ digits, **aren't the chances of getting a true zero (0) or a true one (1) incredibly small?** On the zero condition, every 0–9 digit has to be zero, and on the one condition, the first integer must be one and the rest should default to zero. Does that logic of how a number is represented factor in to how random numbers are generated, or otherwise influence how unlikely that result is? _I don't have great insight into the scope of how random numbers are generated. If you ask, "What is your definition of a random number?" assume a generated digit could exhibitstatistical randomness_. * * * Part 2 If I generated numbers within a range, say 5…8, aren't all integer-value results [5,6,7,8] very unlikely as well? _Note: This whole question has no intended application; I am just curious. This question was partly influenced by the mathematical elaborations ofViHart on Youtube, specifically Proof some infinities are bigger than other infinities, and extra-specifically Cantor's diagonal argument_. _P.S. If any moderators want to move this question elsewhere, that's okay with me!_


title: Count data correlation

text: My dataset is basically a time series - each case is a one day and several events are being observed for each case. I am not experienced in time series analysis but I still would like to get at least some information out of that data. As I understand, all of the variables are counts that express how many times at that particular day the event occurred. In addition to that, none of the variables is normally distributed, but rather Poisson distributed and the vast majority of the values is 0 so I cannot run a simple linear regression. I tried to use the Generalized Linear Model in SPSS, but unfortunately don't know hot to interpret the output :(. I would like to know at least whether there is some relationship between those two variables, so I ran Pearson correlation, is it ok, or have I violated some assumptions? If I would lag one of the variables by 1 day and then run the correlation, would the output mean that occurrence of the event1 (variable a) has an effect on occurrence of the event in variable 2 the day after?


title: SVM Training: Working Set Selection

text: This is related to Joachims's 1998 paper on training SVMs (link to paper). In 11.3, I understand how the term $V(\mathbf d)$ arises as a result of a first order approximation, and why it needs to be minimized. > ![reference to the paper](https://i.imgur.com/W6SPMQD.png) I do not understand how the method in 11.3.2 arises from the description in 11.3. > ![reference to the paper](https://i.imgur.com/WtnyiwT.png) How does the quantity $\omega_i$ represent $V(\mathbf d)$? And how does picking the top and bottom $q/2$ terms imply minimizing $V(\mathbf d)$?


title: About Identification in a 3 equation SEM

text: I got this example and I was wondering about a certain statement: $$ \begin{aligned} y_1 &= \alpha_{12}y_2 + \alpha_{13}y_3 + \beta_{11}z_1 + u_1 \\\ y_2 &= \alpha_{21}y_1 + \beta_{21}z_1 + \beta_{22}z_2 + \beta_{23}z_3 u_2 \\\ y_3 &= \alpha_{32}y_2 + \beta_{31}z_1 + \beta_{32}z_2 + \beta_{33}z_3 + \beta_{34}z_4 + u_3 \end{aligned} $$ It is written that in the first equation of the model we could use all the excluded exogenous variables i.e. $z_2, z_3, z_4$ as instruments for the two endogenous regressors $y_2,y_3$. But as I remembered I can't use $z_2$ nor $z_3$ for $y_2$ because these variables already appear in equation 2. In the same sense I cannot use any of those z's for $y_3$. In my understanding I could use $z_2, z_3, z_4$ for $y_1$ but not for $y_2, y_3$ right?


title: How to calculate accuracy of each feature

text: I read some paper (* about predicting user retention in StumbledUpon) and saw the authors provide a list of features with accuracy of each feature with the following explanation: > ![enter image description here](http://i.stack.imgur.com/PLwYO.png) As > decision trees use information gain to find the best feature, we instead use > decision trees with single features and report their accuracies on the > decision tree classification for retained and non-retained URLs. Table 2 > shows the ranking based on chi-square measure and the decision tree accuracy > on individual features. I am confused about how they can do it. I can just imagine using chi-squared to make a rank and to know the most important features (features evaluation), but in this case I don't know how to calculate the accuracy of each feature. Thank you.


title: Regression adjustment propensity score

text: I have read many articles, but they are hard to understand. Please explain regression adjustment for propensity scores in a mathematical way. My understanding of regression adjustment is to include propensity score as an independent variable into the model (where _Y_ is outcome variable): _Y_ = intercept + coefficent$_a\cdot$Propensity + coefficient$_b\cdot$treatment When we fit the model, if coefficent$_a$ (which is the coefficient of propensity score term in the equation) has a value greater than 0.05, then we can say treatment effect is not significant. Am I correct? However, some articles said treat effect = $Y_t-Y_c-B(X_t-X_c)$. What does that mean? I thought treatment effect is equal to $Y_t-Y_c$, what does the $B(X_t-X_c)$ do here? Does $X_t$ include all the covariances and treatment term, or only covariance terms?


title: Convert an interval variable to ordinal in SPSS

text: I am currently conducting a study on reverse mortgages for my master thesis, and I'd like to convert one interval variable into ordinal in order to run a ordinal probit regression analysis. My interval variable is: "A reverse mortgage is generally a good deal" --> 1-7 Likert scale (where 1 = totally disagree, 4 =neutral, 7 = totally agree). How can I how transform this variable?


title: Understand the reasons of using Kernel method in SVM

text: I understand that one can use kernel functions (i.e. radial kernel) to create non-linear decision boundary. However, there is something with my logic and I am sure there is something that I clearly misunderstood: I understand that Kernel functions operate in a high-dimensional, implicit feature space without ever computing the coordinates of the data in that space, but rather by simply computing the inner products between the images of all pairs of data in the feature space. This operation is often computationally cheaper than the explicit computation of the coordinates. Here is where my logic went wrong: So, I believe $K: \mathbb{R}^N \rightarrow \mathbb{R}$, i.e. it maps input from high-dimensional space to 1-dimensional space. (Not sure if this is correct. ) However, I watched Andrew Ng's lecture videos on SVM, he mentioned that Kernel can also take original data in $\mathbb{R}^1$ and maps it to very high- dimensional feature space $\mathbb{R}^N$. This becomes a contradiction and is very confusing. Please correct my misunderstanding. Thanks.


title: Prediction using Naive Bayes of klaR package fails

text: I am trying to replicate a example that I found in Tom Mitchell's book _Machine Learning_ (1997), using R. It is a example from chapter 6. There are 14 training examples (shown below) of the target concept PlayTennis, where each day is described by the attributes Outlook, Temperature, Humidity, and Windy. Training examples:               Outlook,Temperature,Humidity,Windy,Play     overcast,cool,normal,true,yes     overcast,hot,high,false,yes     overcast,hot,normal,false,yes     overcast,mild,high,true,yes     rainy,cool,normal,false,yes     rainy,mild,high,false,yes     rainy,mild,normal,false,yes     sunny,cool,normal,false,yes     sunny,mild,normal,true,yes     rainy,cool,normal,true,no     rainy,mild,high,true,no     sunny,hot,high,false,no     sunny,hot,high,true,no     sunny,mild,high,false,no      Here's my code:               library("klaR")     library("caret")          data = read.csv("example.csv")          x = data[,-5]     y = data$Play          model = train(x,y,'nb',trControl=trainControl(method='cv',number=10))          Outlook <- "sunny"     Temperature <- "cool"     Humidity <- "high"     Windy <- "true"          instance <- data.frame(Outlook,Temperature,Humidity,Windy)          predict(model$finalModel,instance)      The example tries to predict the outcome for               Outlook=sunny, Temperature=cool,Humidity=high and Wind=strong      The problem is that I am getting a different prediction from the one in the book. Here are the probabilities I've got from my code:               no          yes     0.001078835 0.9989212      Here are the book's probabilities:               no     yes     0.0206 0.0053      My code classifies the unseen data as Yes and the book's classifier classifies it as No. Shouldn't both give the same answer since we are using the same naive Bayes classifier? EDIT: I replicated the example using scikit-learn MultinomialNB classifier and I have got the following probabilities               no    yes     0.769  0.231      which are similar to the **normalized** probabilities of the book. Normalized probabilities of the book               no     yes     0.795  0.205


title: Probability that a continuous random variable assumes a fixed point

text: I'm in an introductory statistics class in which the probability density function for continuous random variables has been defined as $P\left\\{X\in B\right\\}=\int_B f\left(x\right)dx$. I understand that the integral of $\int\limits_a^af(x)dx=0$ but I can't rectify this with my intuition of a continuous random variable. Say X is the random variable equal to the number of minutes from time t that the train arrives. How do I calculate the probability that the train arrives exactly 5 minutes from now? How can this probability be zero? Is it not possible? What if the train _does_ arrive exactly 5 minutes from now, how could it occur if it had probability 0? Thanks.


title: Can the selling price be treated as a random variable?

text: I have some data from my company that I have been looking at. I have been comparing deals that we have "lost" to deals that we have "won". I have been comparing the product averages and standard deviations, and have found that the average and standard deviation (for some industries) are lower for the "lost" deals. We sell each product at set prices, but products can have special discounts when requested which causes "random prices." Some products will be sold higher or lower depending on the customer. When I look at the data, the prices are slightly negatively skewed but seem to be normally distributed. Does this mean that I can treat the price as random and apply the central limit theorem? Also, I will be suggesting more aggressive lower pricing if I can treat the range as normally distributed. Is this logical?


title: Convergence errors in parametric bootstraps (PBmodcomp) of lmer models

text: I am using PBmodcomp from the pbkrtest to perform a parametric bootstrap model comparison. However, for some of the comparisons a warning message stating that the models failed to converge appear. A example of the R script for my lemr models can be found here: http://stackoverflow.com/questions/25111939/error-runnning-parametric- bootstrap-pbmodcomp-on-lmer-objects               modelfit.04b[[1]] <- PBmodcomp(output.04b[[1]], output.04a[[1]])     Warning messages:     1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :       Model failed to converge: degenerate  Hessian with 1 negative eigenvalues     2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :       Model failed to converge: degenerate  Hessian with 1 negative eigenvalues          modelfit.04b[[1]]     Parametric bootstrap test; time: 611.59 sec; samples: 1000 extremes: 156;     large : y ~ age_c + gender_R2 + ibphdtdep + iyeareducc + apoegeno + age_c *          apoegeno + ABCA7_carrier + age_c * ABCA7_carrier + ABCA7_carrier *          apoegeno + age_c * ABCA7_carrier * apoegeno + (age_c | pathid)     small : y ~ age_c + gender_R2 + ibphdtdep + iyeareducc + apoegeno + age_c *          apoegeno + ABCA7_carrier + age_c * ABCA7_carrier + (age_c |          pathid)              stat df p.value     LRT    6.5422  4  0.1621     PBtest 6.5422     0.1568      My models have no issues with convergence outside of the bootstrap. So my question is does this affect the p.value I get from the PBmodcomp function? Is it still valid? Thanks


title: Are there problems with inference using linear regression on observational data with highly skewed distributions of predictor values?

text: I am using a linear regression model to perform inference on some observational data. The samples are from an observational study and highly skewed along some of the dummy variables in the regression. My understanding is this is okay as long as the the errors are uncorrelated and have unequal variance. The predictions for the poorly sampled regions of the independent variable may not be very accurate, but it shouldn't affect the validity of the inference. Is this correct?


title: Correlation of the result of an algorithm to the coefficient of variation of the input

text: Let $A$ and $B$ be two algorithms to solve the same problem and $f(.)$ be the evaluation criteria of the solutions. I use $\frac{f(A(I))}{f(B(I))}$ to find the relative performance of the two algorithms. Both $A$ and $B$ take a list of positive numbers as input. I want to show that algorithm $A$ works much better (compared to $B$) if the input $I$ has a higher coefficient of variation. Assume $X$ is the vector of coefficient of variation of the inputs and $Y$ is the relative performance of the algorithms as described above. Here is a plot of "coefficient of variation" vs "performance ratio" for 2000 randomly generated inputs: ![enter image description here](http://i.stack.imgur.com/Wod3W.png) Correlation of the "coefficient of variation" and "performance ratio" for the sample data is about $0.87$. Now, I have two questions:   1. Based on the above data, can I say the variables are strongly correlated?   2. What sort of other analysis/plot can I use to show the correlation if them clearly? Edit: I asked quite a similar question at Mathematica.SE and got an answer. I'll put the answer below.


title: Grid search and tolerance in libsvm

text: 1. I'm using libsvm and the 3-fold cross validation to select the best C and gamma, but I'm not sure for the range to use in the grid search. Is there any standard way to choose this range? I used: `log2c =(-10:3:4)` and `log2g=exp(-10:3:4)`   2. If I change the `tolerance -e`, the result change, can any one clarify tis point for me? thank you.


title: Statistical test for a value being significantly further from the population mean: is it a Z-test or a T-test?

text: How significant is a value compared to a list of values? In most cases statistical testing involves comparing a sample set to a population. In my case the sample is made by one value and we compare it to the population. I am a dilettante in statistical hypothesis testing confronted with perhaps the most basic problem. It is not just one test but hundreds of them. I have a parameter space, and must do a significance test for every point. Both value and background list (population) are generated for each parameter combination. Then I am ordering this by p-value and find interesting parameter combinations. In fact, the finding of parameter combinations where this p-val is high (nonsignificance) is also important. So let's take one single test: I have a computed value generated from a selected set and a background set of values computed by choosing a random training set. The computed value is 0.35 and the background set is (probably?) normally distributed with a mean of 0.25 and a very narrow std (e-7). I actually don't have knowledge on the distribution, because the samples are computed from something else, they are not random numbers samples from some distribution, so background is the correct word for it. The null hypothesis would be that "the mean of the sample test equals my computed value, of 0.35". When should I consider this to be a Z-test or a T-test? I want the value to be significantly higher than the population mean, therefore it is a single-tailed test. I am a bit confused as to what to consider as a sample: I either have a sample of one (the observation) and the background list as the population OR my sample is the background list and I am comparing that to the whole (unsampled) population which according to the null hypothesis should have the same mean. Once this is decided, the test goes to different directions I guess. If it is a T-test, how do I compute its p-value? I would like to compute it myself rather than using an R/Python/Excel function (I already know how to do that) therefore I must establish the correct formula first. To begin with, I suspect a T-test is a bit too general, since in my case the T-test would be linked to the sample size and would have the form: $$T=Z/s,$$ where $$Z=\frac{\bar{X}}{\frac{\sigma}{\sqrt{n}}}$$ and s is $$s=\hat{\sigma}/\sigma$$, the sample std versus the population std. So I have two cases: either my sample size is the size of the population, which I "guess" would mean I am dealing with a Z-test, or the population statistics (n and std) are unknown but the distribution can be in some way approximated and I am really dealing with a T-test. In any case my following questions are:   1. How do I compute a p-value? (i.e. not using an R/Python/Excel function or p-value table look-up but actually compute it based on a formula, because I want to know what I am doing)   2. How do I decide a significance threshold based on my sample size? (a formula would be nice)


title: Assumptions in general linear model and in multivariate linear regression model?

text: I am reading Wikipedia and some notes I found online, and still not very clear about their definitions. Both the general linear model and the multivariate linear regression model assume $$ Y = X\beta + e $$ where $Y$ is $n \times q$, $X$ is $n \times p$, $\beta$ is $p \times q$, and $e$ is $n \times q$. $n$ is sample size, $p$ is the number of input variables, and $q$ is the number of output variables. What does the general linear model further assume? What does the multivariate linear regression model further assume? For example, following are some assumptions that I can think of on the model, more specifically, on the error matrix $e$.   * the mean of every entry of $e$ is zero   * Any two entries of $e$ in the same column but different rows are uncorrelated   * Any two entries of $e$ in different columns and different rows are uncorrelated   * All entries of $e$ in the same column but different rows have the same variance. If all of above assumptions hold, can the estimation be done for each output variable separately? If not all of the above assumptions hold, is it good to estimate for each output variable separately? Thanks!


title: Model to use in assessing the responsiveness of a categorical variable to policy changes over time?

text: I am interested in determining if the decisionmaking of a particular government body was responsive to policy and statutory changes that occurred at known points in time. I can classify the decisions of this body as 1) rejected proposal; 2) accepted modified proposal; 3) accepted proposal (without change). I could classify the policy or statutory changes dichotomously -- "with deregulation" and "without deregulation," for example. The intervention model discussed in this previous question - "Quantifying effect of a categorical variable in time series analysis" \- seems highly relevant, but I think there are important differences between my dataset and the dataset described there:   1. In the linked question, the OP had a continuous response variable. In my case, I could take the categorical "decision" variable and convert it to a continuous "rate accepted" or "rate rejected" variable, BUT   2. The timing of the response variable may not be as easily coded. There is a lag between the time at which a decision was requested from the government body and the time at which the body made this decision, and this lag varies considerably from observation to observation. Further, decisions generally occur in bunches (a phenomenon I can't readily explain) -- that is, a scheme of intervals used to measure "rate accepted" would have considerable variation in N from interval to interval. In general, I'm interested in whether there's a modification to the intervention model described in the link above that I should investigate, or any other model that might be relevant? A colleague suggested ARIMAX, of which I know nothing. Thanks in advance for any help you may be able to provide.


title: Within group correlation of pearson's residuals

text: I have social network data in which an "ego" names a friend "alter". I am running a regression in R in which attributes of alter are predictors of outcomes for ego. So each observation is dyadic with variable measures for both ego and alter. There are multiple observations for each ego which are accounted for by using a gee model, clustering on ego. The problem is that i have been asked to also account for multiple observations of alter, or least to demonstrate that interdependence among the multiple alters is not impacting the final results. There are multiples of the same alter in the dataset as well as multiples of the same ego. The two options seem to be some kind of cross clustering and I am not sure if that is possible in R. Another option which was suggested was to run a within- group correlation of some sort on the pearson's residuals, with the groups being the alters for each observation. I had considered some sort of ICC but the number of times any individual alter shows up in the dataset ranges from 1-7. As far as I can tell, ICCs expect that the number of measures for each group in the dataset be the same. Does anyone know how to do a within group correlation which can handle groups within which there are differing numbers of measures? I have looked online and have not come across anything that seems to address this. Thanks in advance for any suggestions!


title: Conceptual explanation of regression coefficient?

text: I understand, mathematically, how to get the estimates of coefficients in ordinary least squares. What I am struggling with is coming up with a conceptual, geometric explanation for the correlation part of the simple linear regression coefficient. Say we have the following equation: $$ y = \hat\beta_1x + \hat\beta_0 $$ where $$ \hat\beta_1 = \frac{S_{xy}}{S^2_x} = r_{xy}\frac{S_y}{S_x} $$ and $$ \hat\beta_0 = \overline{y} - \hat\beta_1\overline{x} $$ If I understand correctly, we are essentially just getting estimates of spread ($\hat\beta_1$) and position ($\hat\beta_0$) shifts that can turn $x$ into $y$ as best as possible (least squares). So from what I understand, the ratio of standard deviations ($\frac{S_y}{S_x}$) makes it so $x$ has the same spread as $y$, and then the estimate of $\beta_0$ moves $x$ to the same position as $y$. What I am missing is what the correlation $r_{xy}$ is doing in this linear transformation. If we have shifted the spread and position of $x$ as above, what is $r_{xy}$ representing?


title: How to compute a binomial probability?

text: The percentage of female americans 25 years old and older who have completed 4 years of college or more is 26.1. In a random sample of 200 American women who are at least 25 what is the probability that at most 50 have completed 4 years of college or more?


title: What's in an adjective? Efron's "significant" --> "interesting"

text: This footnote by Efron appears in the chapter on false discovery rate control: > I am trying to avoid the term "significant" for the rejected cases as > dubious terminology even in single-case testing, and worse in the false > discovery rate context, preferring instead "interesting". \--Bradley Efron > "Large Scale Inference" 2010, p.47 What are the arguments pro and con adopting this change in nomenclature? Is it significant that such a prominent statistician is avoiding "significance"? Biomedicine and healthcare in particular are domains where quantitative evidence is largely based on association studies, whether RCTs or observational. "Links" between therapies and outcomes, or biomarkers and disease states, depend on hypothesis testing framework. Could adopting "interesting" on a large scale induce reconsideration of the strength of evidence associated with of such linkages, even if p-value or z-scores still quantify the results?


title: Bayesian Analysis in R/WinBUGS and SAS (Proc MCMC)

text: I have a Bayesian model here in R/WinBUGS. This data "pumps" has two columns, time and fail. Each observation is a pump and fail[1] tells how many failures there were in time[1]. I want to know how I can write this same model in SAS using Proc MCMC. Thank you!               pump <- function() {         for(i in 1:10){             fail[i] ~ dpois(lam[i])             theta[i] ~ dgamma(a,b)             lam[i] <- theta[i]*time[i]         }         a ~ dgamma(1.5,.25)          b ~ dgamma(1.5,.25)     }          filename <- file.path(tempdir(),'pump.bug')      write.model(pump,filename)     time <- pumps[,1]     fail <- pumps[,2]          data <- c('time','fail')     parameters <- c('fail','lam','theta','a','b')     pump.sim <- bugs(data, inits=NULL, parameters, model.file='pump.bug',         n.iter=20000, n.burnin=1000, n.chains=1, n.thin=1, debug=T)


title: MLE vs least squares in fitting probability distributions

text: The impression that I got, based on several papers, books and articles that I've read, is that the recommended way of fitting a probability distribution on a set of data is by using maximum likelihood estimation (MLE). However, as a physicist, a more intuitive way is to just fit the pdf of the model to the empirical pdf of the data using least squares. Why then is MLE better than least squares in fitting probability distributions? Could someone please point me to a scientific paper/book that answers this question? My hunch is because MLE does not assume a noise model and the "noise" in the empirical pdf is heteroscedastic and not normal.


title: Why is the expected value of the number of trials before the first success larger than the number of trials with a 50% probability of success?

text: For a Bernoulli distribution with parameter p, the number of trials with a 50% probability of at least one success is about (1/p) * ln(2). But the expected value of the corresponding geometric distribution is 1/p. Intuitively I would expect them to be the same (1/p). What am I missing?