title: Interpretation of Coefficient

text: In Poisson Regression, suppose we have the following regression equation: $\ln(E(Y_i|X_1)) = \beta_0+ \beta_{1}X_{1}$ where $Y$ is the number of heart attacks. Also suppose $X_1$ is a binary variable (e.g. smokers versus non- smokers). How would you interpret the coefficient $\beta_1$? Smokers, on average, have $\exp(\beta_1)$ more heart attacks than non-smokers?


title: Constructing a covariance matrix and finding first principal component

text: So far in the topic of principal component analysis, I seem to be jumping around in understanding these various topics, and have a few questions that I need to help solidify my overall understanding. What I know how to do: calculate eigenvalues from a covariance matrix. My first question: Given a 2×2 matrix, for example, let's say: $$\left[\begin{matrix}2&5\\\4&2\end{matrix}\right]$$ How can I compute the covariance matrix to move forward in finding the eigenvalues, etc.? And in the problem I am trying to solve, I am required to find the value of first principal component given a covariance matrix. It says the first principal component corresponds to the eigenvector that captures the maximum variability. Second question: I'm trying to understand the process behind this, so do I calculate the eigenvalues first, then use that to get the eigenvector? Then from there it says something about converting it to a unit vector in order to find the PC1 (first principal component). I don't want direct answers, but more of the idea and process that I am missing when trying to solve these. Thanks.


title: Normalized correlation with a constant vector

text: I am confused how to interpret the result of performing a normalized correlation with a constant vector. Since you have to divide by the standard deviation of both vectors (reference: http://en.wikipedia.org/wiki/Cross- correlation ), if one of them is constant (say a vector of all 5's, which has standard deviation of zero), then the correlation is infinity, but in fact the correlation should be zero right? This isn't just a corner case, in general if the standard deviation of one of the vectors is small, the correlation to _any_ other vector is very high, which obviously doesn't make sense. Can anyone explain my misinterpretation?


title: What is the Hotelling $T^2$ used for?

text: Hotelling's $T^2$ distribution arises in testing differences between means of different populations. But is it often used? Can it be implemented in a modeling procedure, let's say, logistic regression?


title: Chi-square of a sample with respect to a known population's marginal distribution?

text: This is a question to reflect some very basic understanding of the logic of the chi-square for a planned introductory essay. _[update] I've tried to improve the question with clearer examples and better focusing of the problems upon which I stumbled. It might still be a bit weak but I can't do better at the moment_ I'm used to the chi-square as a measure for the deviance from an expected discrete distribution, where the expected frequencies are computed via the marginal frequencies of a crosstab over the empirical data of two categorical items. Let's discuss the 2x2-table here which has one degree of freedom. Assume the marginal row frequencies being [53,47] and the marginal column frequencies being [40,60] giving the expected frequencies as $$ {\bf A}: \text{ expectation for cell-frequencies by equal ratios} \\\ \small \begin{array} {r|rr|r} &0&1&all\\\ \hline 0&21.2&18.8&40\\\ 1&31.8&28.2&60\\\ \hline all&53&47&100 \end{array} $$ First I began to think what it means, that the empirical contingency table is bounded by the minimal entry of the marginal frequencies: in this case the cell [0,0] can vary between _0_ and _40_ only, so we have at most 41 possible outcomes depending on the possible frequencies in that cell. If we assume a normal random process, which generates that frequencies, the (expected) frequency in this cell should be centered around the mean of them _20_ : $$ {\bf B}: \text{ means of possible ranges for frequencies in cells}\\\ \small \begin{array} {r|rr|r} &0&1&all\\\ \hline 0&20&20&40\\\ 1&33&27&60\\\ \hline all&53&47&100 \end{array} $$ This expectation in cell[0,0] is not equal to the expected frequency of 21.2 . First question: can the latter concept be brought/translated into the first (and standard) one? How could the difference (and relation or possible non- relation) between that two concepts be best explained? To answer this myself I went back one further step and asked, what does it mean at all to base the chi-square computation on the sample's marginal distribution if this distribution is itself subject of random... So I generated a "population" data set with _N=10000_ cases distributed exactly like our table **A** and took _1000_ random samples each with _n=100_. Over the _1000_ samples I've got different marginal frequencies and consequently, each sample has different parameters for its expected frequencies and so for its possible chi-square.   The table of the deviances of the empirical marginal frequencies from the population marginal frequencies was $$ {\bf C}: \text{ deviations of empirical marginal frequencies from population} \\\ \small \begin{array} {r|rrrr} \text{dev } & \text{dev from}& \text{dev from}& \text{dev from}& \text{dev from}\\\ \text{value} & 53 & 47 & 40 & 60 \\\ \hline -18&0&0&0&1\\\ -15&1&1&1&1\\\ -14&1&1&0&0\\\ -13&2&4&1&1\\\ -12&2&6&5&4\\\ -11&10&9&5&5\\\ -10&8&12&7&14\\\ -9&13&15&17&17\\\ -8&27&20&21&23\\\ -7&34&38&26&31\\\ -6&45&35&31&37\\\ -5&53&43&62&50\\\ -4&46&46&50&49\\\ -3&81&67&67&63\\\ -2&55&83&84&67\\\ -1&85&71&97&80\\\ 0&86&86&83&83\\\ 1&71&85&80&97\\\ 2&83&55&67&84\\\ 3&67&81&63&67\\\ 4&46&46&49&50\\\ 5&43&53&50&62\\\ 6&35&45&37&31\\\ 7&38&34&31&26\\\ 8&20&27&23&21\\\ 9&15&13&17&17\\\10&12&8&14&7\\\11&9&10&5&5\\\12&6&2&4&5\\\ 13&4&2&1&1\\\14&1&1&0&0\\\15&1&1&1&1\\\18&0&0&1&0 \end{array} $$ The table of occurrences of chi-square for those deviations of the empirical marginal frequencies from the population's is $$ {\bf D}: \text{ deviations of the empirical marginal frequencies} \\\ \text{ from that of the population }\\\ \text{(in terms of}\ \chi^2\ \text{values)} \\\ \begin{array} {r|rr} & & \text{backwards} \\\ \chi^2 \text{-value} & \text{freq} & \text{cum freq} \\\ \hline 0.0&12.8&100.0\\\ 0.5&22.0&87.2\\\ 1.0&12.7&65.2\\\ 1.5&11.2&52.5\\\ 2.0&10.5&41.3\\\ 2.5&7.8&30.8\\\ 3.0&3.3&23.0\\\ 3.5&4.8&19.7\\\ 4.0&3.2&14.9\\\ 4.5&2.5&11.7\\\ 5.0&2.3&9.2\\\ 5.5&1.0&6.9\\\ 6.0&1.7&5.9\\\ 6.5&1.2&4.2\\\ 7.0&0.9&3.0\\\ 7.5&0.6&2.1\\\ 8.0&0.5&1.5\\\ 9.0&0.1&1.0\\\ 9.5&0.5&0.9\\\ 10.0&0.1&0.4\\\ 10.5&0.2&0.3\\\ 17.0&0.1&0.1 \end{array} $$ I'm stumbling at the simple fact that we compute the chi-square as deviation from an expected frequency - but where the expected frequency is based on an _empirical_ marginal frequency which is _itself_ subject of a random process - and for instance has some confidence-interval when I infer from the sample onto the population.   So - in reverse - having an empirical marginal distribution in our single empirical sample the conclusion to the population's marginal distribution is a matter of confidence intervals. _(Is here also a maximum-likelihood aspect lurking around anywhere ?)_   This reminds me of the praxis, that we use the sample's variation as estimate of the population's variation-parameter, and do tests based on this assumption. Q: In the justification/formula for the chi-square-distribution as basis for the significance test - can we find some point, where the randomness of the marginal frequencies is reflected in the formulae?


title: Power for two sample t test

text: I am trying to understand power calculation for the case of the two independent sample t-test (not assuming equal variances so I used Satterthwaite). Here is a diagram that I found to help understand the process: ![enter image description here](http://i.stack.imgur.com/S6P3Y.jpg) So I assumed that given the following about the two populations and given the sample sizes:               mu1<-5     mu2<-6     sd1<-3     sd2<-2     n1<-20     n2<-20      I could compute the critical value under the null relating to having 0.05 upper tail probability:               df<-(((sd1^2/n1)+(sd2^2/n2)^2)^2) / ( ((sd1^2/n1)^2)/(n1-1) + ((sd2^2/n2)^2)/(n2-1)  )     CV<- qt(0.95,df) #equals 1.730018      and then calculate the alternative hypothesis (which for this case I learned is a "non central t distribution"). I calculated beta in the diagram above using the non central distribution and the critical value found above. Here is the full script in R:               #under alternative     mu1<-5     mu2<-6     sd1<-3     sd2<-2     n1<-20     n2<-20               #Under null     Sp<-sqrt(((n1-1)*sd1^2+(n2-1)*sd2^2)/(n1+n2-2))     df<-(((sd1^2/n1)+(sd2^2/n2)^2)^2) / ( ((sd1^2/n1)^2)/(n1-1) + ((sd2^2/n2)^2)/(n2-1)  )     CV<- qt(0.95,df)               #under alternative     diff<-mu1-mu2     t<-(diff)/sqrt((sd1^2/n1)+ (sd2^2/n2))     ncp<-(diff/sqrt((sd1^2/n1)+(sd2^2/n2)))               #power     1-pt(t, df, ncp)      This gives a power value of 0.4935132. Is this the correct approach? I find that if I use other power calculation software (like SAS, which I think I have set up equivalently to my problem below) I get another answer (from SAS it is 0.33). SAS CODE:               proc power;           twosamplemeans test=diff_satt              meandiff = 1              groupstddevs = 3 | 2              groupweights = (1 1)              ntotal = 40              power = .             sides=1;        run;      Ultimately, I would like to get an understanding that would allow me to look at simulations for more complicated procedures. EDIT: I found my error. should have been 1-pt(CV, df, ncp) NOT 1-pt(t, df, ncp)


title: Comparison negative binomial model and quasi-Poisson

text: I have run negative binomial and quasi-Poisson models based on an hypothesis testing approach. My final models using both methods have different covariates and interactions. It seems that there are no patterns when I plot my residuals in both cases. Thus, I was wondering which test I could use to see which model fits my data better as the quasi-Poisson does not have any likelihood or AIC… Also, I have a lot of overdispersion which makes me think that the negative binomial would be more appropriate, but I don't know if I can choose my model based on common sense…


title: Questions about the order statistics of uniform distributions

text: I refer to the Simes (1986) paper found here. In this setting, $P_{(1)}$ through $P_{(n)}$ are the order statistics of $n$ independent Uniform$[0,1]$ random variables and, for $0\le \alpha \le n$, $$A_n(\alpha) = \Pr\\{P_{(j)}\gt j\alpha/n; j = 1, 2, \ldots, n\\}.$$ On page 2, I'm unsure of how the integral in the proof is obtained. It asserts that because (a) the $P_{(j)}/P_{(n)}$ are the order statistics of $n-1$ independent uniform variables and (b) the distribution function of $P_{(n)}$ is $p^n$ ($0\le p \le 1$), then $$A_n(\alpha) = \int_\alpha^1 A_{n-1}\left(\frac{\alpha(n-1)}{p n}\right) n p^{n-1} dp.$$ Simes mentioned that he thanks the referee for a shorter version of the proof and I'm assuming that it cut off a corner that would lead to my better understanding. I suppose that my two specific questions are these… 1: Why consider the $P_{(1)}/P_{(n)} ..... P_{(n-1)}/P_{(n)}$ order statistics of $n-1$ independent uniform random variables on $(0,1)$ independent of $P_{(n)}$? 2: How is the $A_{n-1}$ {$\alpha$$(n-1)/np$} obtained in the integral? Many thanks if you are able to help me figure this out!


title: Market Basket Analysis: comparing rules between two models

text: Given two independent MBA models 1 and 2 (each model is a set of rules with calcualted support, confidence and lift metrics) that were generated on subsets of large population of transactions, how to effectively compare rules between models? In particular, how to detect the biggest movers (up and down) from one model to the other. So far, my approach has been as follows:   * set thresholds for support, confidence, and lift (e.g. given rule $ A: support(A) > 0.001, confidence(A) > 0.1, lift(A) > 10 $) and filter out rules below these thresholds from both models;   * for each rule $ A $ compare its log-lifts between models 1 and 2: $$ \mid log(lift_1(A)) - log(lift_2(A)) \mid $$ and select rules with largest difference The goal is to identify set of rules that change the most between two models. I am looking for both validation of given approach and better alternatives - thank you!


title: Linear regression forecast underestimation

text: I have the following multiple linear regression model:               Call:     lm(formula = Y ~ X1 + X2 + X2 + X3 + X4 + X5 + X6 + X7,          data = my.model, na.action = na.omit)          Residuals:         Min      1Q  Median      3Q     Max      -43.836  -1.507   0.010   1.485  46.231           Coefficients:                    Estimate Std. Error t value Pr(>|t|)         (Intercept)  -0.0244927  0.0245157  -0.999    0.318         X1           -0.3484619  0.0134383 -25.931   <2e-16 ***     X2            0.1195273  0.0106940  11.177   <2e-16 ***     X3            0.1224587  0.0108849  11.250   <2e-16 ***     X4           -0.0010173  0.0028247  -0.360    0.719         X5            0.5496942  0.0156319  35.165   <2e-16 ***     X6           -0.2287941  0.0145018 -15.777   <2e-16 ***     X7           -0.2315801  0.0146361 -15.823   <2e-16 ***     X8            0.0005465  0.0003595   1.520    0.128         ---     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1           Residual standard error: 2.936 on 35849 degrees of freedom       (12534 observations deleted due to missingness)     Multiple R-squared: 0.05968,    Adjusted R-squared: 0.05947      F-statistic: 284.4 on 8 and 35849 DF,  p-value: < 2.2e-16       The model is affected by multicollinearity but my question is about the forecast, so this shouldn't be an issue. I checked the absolute values of my model forecast and compared against the actual Y absolute values. The average of the absolute predicted values is significantly lower than the absolute observed values mean:               > lm1.predict = predict(lm1, mydata)     > mean(abs(lm1.predict))     [1] 0.3294776     > mean(abs(mydata$Y))     [1] 1.206954      Does this mean that the linear regression variables I am using tend to underestimate the outcomes? Can any other conclusion be derived from this simple comparison? **EDIT** Another way to look at this is to calculate the absolute difference between each observation and the relative outcome:               > mean(abs(mydata$Y - lm1.predict))     [1] 1.208378      These are the diagnostic from the regression: ![enter image description here](http://i.stack.imgur.com/Y7zmf.jpg)


title: Statistical testing of stochastic processes or multivariate testing?

text: I've two groups of data, each containing an equal number of tuples. Each tuple is filled from observations of a phenomenon at regular intervals and are drawn from a discrete set. For each group we are observing independent phenomena, even though in each group they have the same characteristics. I want to test whether these two groups are created from "different" categories of phenomena. To make it clearer, let say that group A and group B contains temperature observations over a specified period of time and coming from different places of a country, with locations from group A coming from the north, while those in group B coming from the south. I want to test whether group A and group B shows a different "behaviour". Is that possible? I tried to search the Internet and tried to figure out the correct set of terms to search for in order to get some hints, but to no avail. Thank you for any help.


title: Confidence intervals using standard deviation - mistake or misunderstanding?

text: Hedges et al. 1999, _Ecology_ **80** : 1150-1156 reintroduces the old concept of using the natural logarithm of response ratio for ecologists as preferred statistics over p-values in comparative experiments and meta-analyses. Response ratio, $R = \bar{X}_{E =treatment}/\bar{X}_{C=control}$; $L = ln(R) = ln(\bar{X}_{E}) - ln(\bar{X}_{C})$. They indicate that taking a logarithm of the response ratio helps to tone down the violations against assumptions of normal distribution and homoscedasticity, which are evident in most experimental data sets. Authors write: "If $\bar{X}_{E}$ and $\bar{X}_{C}$ are normally distributed and $\bar{X}_{C}$ is unlikely to be negative, then L is approximately normally distributed with mean approximately equal to the true log response ratio and variance, v, approximately equal to" $\frac{(SD_{E})^2}{n_{E}\bar{X}_{E}^2}$ + $\frac{(SD_{C})^2}{n_{C}\bar{X}_{C}^2}$ They continue: "An approximate 100(1-$\alpha$)% confidence interval for the individual log response ratio parameter $\lambda$ is given by" $L - z_{\alpha/2}\sqrt{v}\leq \lambda \leq L + z_{\alpha/2}\sqrt{v}$ "where $z_{\alpha/2}$ is the 100(1-$\alpha/2$)% point of the standard normal distribution, and the corresponding confidence interval for the (unlogged) response ratio $\rho$ is obtained by taking the antilogs of the confidence limits for the log response ratio." In here, as I understand it, the authors indicate that standard deviation instead of standard error should be used to calculate the confidence intervals. Is this my misunderstanding or a mistake in the article?


title: What type of statistics to use to show correlation between two discrete sets of data?

text: I have two variables, each with the same set of 5 unranked possible values (let's call them A/B/C/D/E), and a set of data such as the following (A/A, A/D, B/B, D/D, E/E, B/A etc), where the first letter is the value of the first variable and second letter is second variable. How would I show correlation between the statistics i.e. to show high probability of the two variables being the same letter? Sorry if this is not clear, and thanks in advance for any help! :)


title: Why does this criterion characterize the median of a continuous random variable?

text: If $X$ is a continuous random variable, then let $$\min_a{\mathbb{E}\:| X - a |} = \mathbb{E}\: | X - m |$$ Why is $m$ the median of $X$?


title: Testing classifier with binomial test when group sizes are unequal

text: I have data from 50 human subjects, who are divided into groups A and B (30 participants are in group A and 20 participants in group B). I also have a range of measurements from each subject. I have used a machine learning algorithm (SVM) to predict the subject groups from the features using leave- one-out cross-validation. I would like to have a p-value for my classifier. If I have understood correctly, I can treat each prediction as a Bernoulli trial, and test for significance using a binomial test. However, I'm a little confused about the details. The paper (pdf) I was following, in the section 'result significance', assumes that the subject groups are of equal size. That's not the case here. What do I need to do to make the test work?


title: Relationship between correlation of x vs. x+y and y vs. x+y

text: I am working with some data in which two component scores X and Y are added together to make a total score X + Y. Is there any known relationship between the correlation R for X vs. X + Y and the R for Y vs. X + Y? I've played around a little bit, but I am not seeing anything except in the case where X and Y themselves are uncorrelated. (In that case, the R2 values sum to 1.) ETA: With all due respect, this question is NOT a duplicate of the one already cited. I made abundantly clear that I am not talking about cases in which the variables have no correlation, which is a key feature of the other question. I am asking a much broader question.


title: Is there a reliable algorithm to identify a low-dimensional log-linear model that fit high-dimensional data (if one exists)?

text: I am working with data from a complex survey, with over 100 variables per observation. From this I will be selecting 30 to 40 variables. All variables are categorical, with number of levels ranging from 2 to 23, with a mean of about five. I would like to summarize my data with a log-linear model that highlights the important correlations between these variables without prematurely imposing a causal structure on the model. I am still doing coding and data cleanup, but I am pretty sure that I will not be able to implement the complete saturated model and then work my way down by testing whether I can knock off high-level interaction terms, as most of the texts seem to advise me to do. If I am understanding this correctly, the saturated model would contain about 900 quintillion coefficients. I'm using R, and I don't think all those coefficients are going to fit in my 4 gig of RAM. Also, in the smallest unit of geography that I am examining, there are only about 50,000 observations. 19 quadrillion coefficients per observation will leave me short a few degrees of freedom. So if this project is to be feasible, I would need to find a model that fits reasonably well without needing more than a 5th- or 6th-order interaction for any combination of variables. My hope was to build up some blocks of mutually independent variables that I could count on excluding from groups of variables with higher-order interactions. However, considerable cold water was thrown on this notion by the answer to my last question, that mutually independent variables can nonetheless be correlated conditional on additional variables not in the original set. See: Does independence imply conditional independence?. This implies that, for example, three mutually independent variables could nonetheless participate in obligatory 4- or 5-way interactions So my question is, does anybody know of good algorithms for finding reasonably well-fitting log-linear models for high-dimensional data that are not egregiously invalid and that do not require that you start with a fully saturated model and work down? Here is one that I was wondering about, but I am not too attached to it, and you should skip the rest of these ramblings if you are aware of an algorithm that is known to work. Is it reasonable to assume that models will not completely skip two or more levels? For instance, if five variables are not nested in any models that require 6th- or 7th-order interactions, is it reasonable to assume that, in real data not artificially designed to display odd patterns of correlation, that they will not be nested in any models that require interactions of 8th order or higher? Suppose, for example, I end up with 30 variables. Then I could examine all the 4000-odd 3-variable combinations and identify those that require a three-way interaction. For those that do not, I could test all the 4-way combinations that contain those three variables, and if there are none, I can truncate the model at the two- way level for those three variables, as a 5-way interaction would involve skipping two levels. Similarly, for the three-variable combinations that do require a three-way interaction, I can test all the five-variable combinations that contain them. If there are none, truncate. Whenever I find an n-way interaction that is required, I search the n+2 variable space, adding new highest-required-interaction-levels wherever I find them, and truncating wherever I do not. In the end, I have a single structure of hierarchical models where we always know for any interaction term not contained in a higher-level interaction that there are no interactions that are one or two levels higher, and so, by assumption, that there are no required higher-level interactions. Of course, I would have to do all this on just half my data, so I could use the other half to test whatever hypothetical structure I might find. All that would be rather clever if I actually had any reason to believe hierarchical log-linear models can not skip two steps. I suppose that if none of the triples that do not need a 3-way interaction term are embedded in 5-ples that require a 5-way interaction without there being one or more intervening required 4-way interaction terms, that would lend some support to the hypothesis. I could actually test that. I would be very grateful for whatever thoughts people could offer. Respectfully yours, andrewH


title: One group, four conditions and six tests. Testing the effect of altitude

text: I'm testing the effect of four different altitudes (0m, 1000m, 2000m and 3000m) on repeated sprint performance (6 sprints at each altitude). I have 8 test subjects who will be tested at each altitude. That means I'm getting 8x4x6=192 sets of data. What I want to investigate is the difference in performance at the 4 different altitudes and wether or not there actually IS a difference. How do I best represent this and what would be the relevant statistical tests to run? To me it is a bit confusing since I have multiple conditions, and I want to compare each condition with the other conditions (but I'm not sure about the best way to do this). Thanks in advance for your help. Let me know if I can clarify.


title: Question of Holt-Winters, parameter chosen

text: I am using Holt-Winters to do a time-series forecasting. The package chose gamma equal to 1 for me. I am wondering what that means. The prediction works pretty well overall. When will you use this method and how do you choose the parameters of Holt- Winter?


title: Detecting variants of unknown significance

text: I have two groups of people, a healthy population which has carriers of disease and non-carriers, and I have a sick population. I know the frequency of genetic variants in both groups. I want to know if I can say anything about the variants with any statistical certainty and how I would do that.               Variant      Frequency in Healthy Population    Frequency in Sick Population     X            30                                 25     Y            5                                  700     Z            600                                600      Each group (Healthy and Sick) has 1000 people in it. Can I say anything about each variant?


title: Best practices for dealing with shifting, inconsistent seasonality

text: This question is related to a previous post I've looked at (Calculation of seasonality indexes for complex seasonality), but deals with more granular data (daily instead of weekly), and transforming holiday seasons (instead of just holiday) creating a larger question about the reasoning behind dummy variables. As an introduction, I'm looking to forecast daily retail sales for a highly cyclical e-commerce company. Each year, buying patterns shift for several reasons:   1. Calendar shift (where, for example, 5/17/13 is a Friday, while 5/17/12 is a Thursday);    2. Shifting of holidays set by the day of the week (i.e. Labor Day, Thanksgiving, Mother's Day);   3. Shifts in the time and length of shopping seasons (i.e. in 2012, there were 32 calendar days, 21 workdays, and the season started on Nov. 28 in the holiday shopping season, while this year, there are 26 calendar days, 17 workdays, and the season starts on Nov. 29). These changes affect the overall shape and size of the seasonal curve, and not just the behavior at or near the holiday.  Among the retail forecasters I've encountered, there's a few options for dealing with forecasting shifting seasonality:   * One is to only use exact calendar equivalent years to parcel out seasonality -- for example, 2013 exactly matched 2002 in terms of the placing of holidays. Then, you apply a 365-day ARIMA on the data. But, in e-commerce (and overall) a lot has happened since 2002 that affects how seasonality looks.    * The second is exponential smoothing with dummy variables. But, as consumer behavior doesn't just change on the holiday (or the days approaching it), but relates to how far the holiday is in relation to other shifting dates (and the day of the week), you hit the problems of either creating a boatload of dummy variables for the distance from and to each holiday that can over-specify your model (and not make it easily generalizable to the year you want to forecast), or get a "general" seasonality that usually under-preforms. You could also do a fractional polynomial to approximate the curve by days until (or between) one holiday or another, but this creates the usual problems with fracpolys -- notably over-fitting and wildly strange results at certain points in and out of the sample.    * The third option I've seen is to adjust the data from previous years to match the size of each holiday season by, in the case when the holiday shopping season is a week shorter than previous years, taking out a week and then adding back that revenue across the remaining season, or doing the reverse when the season is a week longer. After that's done, you can do a seasonal decomposition to parcel out the seasonality. The issue with this is that it involves a ton of subjectivity, and assumes that shoppers will evenly spread the lost week, rather than "back load" their behavior as the holiday more quickly approaches. Fixing the "back load" would then create even more subjectivity of when shoppers will start responding to the proximity to the holiday (such as Christmas).  So, I'm wondering the experience of other forecasters (especially those with retail experience) with this issue of shifting seasonality by day. Looking forward to your thoughts. I can also add some dummy data if you want to see exactly what I'm talking about if you'd like.


title: ARIMA and external regressors in SAS and R

text: So I remember reading somewhere that when we have external regressors, `auto.arima` cannot make correct predictions for the order of difference for either seasonality or the main time series itself (correct me if I'm wrong!) Now, I'd like to know whether we'd need to difference the external regressors as well? Also, in the case of having external regressors (a few time series and a few dummies for seasonal patterns in those time series), can `auto.arima` even calculate the optimal MA and AR? Also, I have weekly seasonality as well as quarterly and yearly seasonality; since I can't specify that many seasonalities in auto.arima, I'm inputting a lot of dummy variables for quarters and months; will that yield mathematically correct results? Further, for those of you who have worked with SAS, when using the forecast procedure and estimating the input variables (the external regressors), does it automatically calculate the MA and AR for each external regressor?


title: What is delta in the maximization step in this EM algorithm?

text: ![EM Algorithm](http://i.stack.imgur.com/Oxvip.png) The algorithm is used to classify english vs non-english tweets from unlabeled data. Given n observed tweets (x1 ... xn) where each tweet xi is a collection of d words (xi1 ... xid). y is the class (english/non-english). We use EM to estimate the parameters (thetas). It is used in this paper: http://www2013.wwwconference.org/companion/p593.pdf


title: Importance sampling of finite path of stochastic difference equation

text: Before passing to question, let me briefly recap what's importance sampling of random variables is about. Suppose $\xi$ is a real-valued random variable with density $f$, and let $g:\Bbb R\to \Bbb R$ be some function. The task is to use Monte-Carlo to compute the integral $$ \mathsf E[g(\xi)] = \int_\Bbb R g(x)f(x)\mathrm dx \approx \frac1N\sum_{i=1}^N g(\xi^i) $$ where $\xi^i$ are iid random variable distributed as $\xi$. Since the approximation error may have a big variance, the idea of the importance sampling is to change a density of $\xi$ to some $\hat f$ and use that $$ \int_\Bbb R g(x)f(x)\mathrm dx = \int_\Bbb R g(x)w(x)\hat f(x)\mathrm dx $$ where the weighting function is given by the fraction $w = f/\hat f$. As a result, $$ \mathsf E[g(\xi)] = \mathsf E[g(\hat \xi)w(\hat \xi)] \approx \frac1N\sum_{i=1}^N g(\hat\xi^i)w(\hat\xi^i) $$ where $\hat\xi^i$ are iid with densities $\hat f$. Then one runs the optimization problem to find the best choice of $\hat f$, that is the one which minimizes the variance. * * * In my case I have a similar problem. Let us consider a discrete-time stochastic process $X$ with a state space $E$, given by the stochastic difference equation of the form $$ X_{k+1} = r(X_k,\eta_k), \quad X(0) = x\in E, \tag{1} $$ where $\eta_k$ is a sequence of iid real-valued random variables with some density $h$, and $r$ is a jointly measurable function. Let $\mathsf P$ be the induced probability measure on $E^{n+1}$ and let $A$ be measurable subset of $E^{n+1}$. I am interested in using importance sampling to evaluate $\mathsf P(A)$. The method of importance sampling described above easily extends to the case of random elements with a range $\Bbb R^m$, when their density is precisely known. In my case $E$ is a subset of $\Bbb R^m$, and I know the function $h$, but it would be almost impossible to get an expression of the density of $X = (X_0,X_1,\dots,X_n)$ since the function $r$ may have an extremely complicated shape. Due to this reason, I restrict myself just to change of $\mathsf P$ which is induced by the change of the distribution of $\eta$: $$ \begin{align} \mathsf P(A) &= \int_{E^{n+1}} 1_A(x_0,\dots,x_n)\mathsf P(\mathrm dx_0\times\dots\times\mathrm dx_n) \\\ & = \int_{\Bbb R^n}1_A(R(y_0,\dots,y_n))h(y_0)\dots h(y_{n-1})\mathrm dy_0\times\dots\times \mathrm dy_{n-1} \\\ & = \int_{\Bbb R^n}1_A(R(y_0,\dots,y_n))w(y_0,\dots,y_{n-1})\hat h(y_0)\dots \hat h(y_{n-1})\mathrm dy_0\times\dots\times \mathrm dy_{n-1} \\\ & = \mathsf E[1_A(\hat X_0,\dots,\hat X_n)w(\hat \eta_0,\dots,\hat \eta_{n-1})] \end{align} $$ where $R$ is a function which for every noise realization $(\eta_0,\dots,\eta_{n-1})$ gives the path of the process $(X_0,\dots,X_n)$ according to $(1)$, and further $\hat h$ is a new density, $\hat \eta_k$ are distributed according to $\hat h$, $$ (\hat X_0,\dots,\hat X_n) \sim R(\hat \eta_0,\dots,\hat\eta_{n-1}) $$ is distributed according to the new distribution of the noise. * * * My question is the following: in case I need to do importance sampling of $(1)$, I need to perform the change of measure from $\mathsf P$ to some $\hat{\mathsf P}$ in such a way, that I get an explicit shape of the Radon- Nikodym derivative $w = \mathrm d\mathsf P/\mathrm d\hat{\mathsf P}$. Is the way I described the only one which guarantees the knowledge of explicit shape of $w$? I am pretty sure, there shall be some literature on the importance sampling of stochastic difference equations that deals with such problems, but I didn't find anything particularly appropriate yet. * * * Was asked on MSE but didn't get much attention.


title: Multidimensional scaling

text: I am running into some problems performing a multidimensional scaling image. First of my dataset is quite large (330.000 fields: 33000 rows, 10 columns). The output image needs to contain 10 dots, 1 for each column. The fields in the dataset contain intensity values which first need to be converted to a distance matrix, however when doing so I run into a memory limit in R. I was hoping you could help solve this problem or perhaps point me to a program that can handle these kind of data sizes.


title: State space models for time series forecasting

text: I am new to time series forecasting and have been slowly working my way through the different approaches available. I've so far mainly been using ets and arima models available in the R forecast package. As I need to forecast at a daily level taking into consideration weekends and public holidays, I have also been exploring arimax models with xreg values. As far as I can tell, it does not seem possible to use xreg values with ets models. With my sample data ets also can't forecast the step function that occurs at weekends accurately - the smoothing means the only saturday is forecast correctly but sunday is also incorrect. I cam across a series of slides from Nate Derby - http://www.sas.com/offices/NA/canada/downloads/presentations/Victoria2008/Time.pdf The final slides elude to the fact that state space models might make it easier to model timeseries. How powerful and implementable are state space models? Is this something I need to seriously consider?


title: Chi Square Test for Independence in R and Python

text: Consider the following R code and output:               row1 = c(0,23,0,0)     row2 = c(0,1797,0,0)     data.table = rbind(row1, row2)     chisq.test(data.table)              Pearson's Chi-squared test          data:  data.table     X-squared = NaN, df = 3, p-value = NA      Now consider the same in Python:               import scipy.stats     scipy.stats.chi2_contingency([[0,23,0,0], [0,1797,0,0]])          Traceback (most recent call last):       File "<stdin>", line 1, in <module>       File "/usr/lib/python2.7/dist-packages/scipy/stats/contingency.py", line 236, in          chi2_contingency         "frequencies has a zero element at %s." % zeropos)     ValueError: The internally computed table of expected frequencies has a zero element at [0, 0, 0, 1, 1, 1].      Is this expected behaviour? Should I just trap for the error in Python. A search for the message "The internally computed table of expected frequencies has a zero element at" did not reveal anything useful.


title: Modelling dependency between random variables

text: Assume we have 2 random variables $(X,Y)$ with same distribution $G(.)$. Suppose we are interested in a joint probability of these r.v.'s without any information about dependence or correlation between these random variables. Would it help to assume $X$ as a linear (or non-linear) function of $Y$? And lastly, if we knew $corr(X,Y)$ and $F_X(.),F_Y(.)$ could we derive $F_{X,Y}(.)$? Thank you in advance...


title: Assessing normality of distribution

text: I'm having trouble understanding how to assess normality for discrete data. I understand the the K-S test is available to test the normality distribution for continuous data. However, what would be the equivalent for discrete data that have a Poisson or binomial shape? Would these be the chi-square and binomial test respectively? I'm not sure if these are testing significance rather than normality. Any help is greatly appreciated!


title: Appropriate clustering techniques for temporal data?

text: I have temporal data of activity frequencies. I want to identify clusters in the data that indicate distinct periods of time with similar activity levels. Ideally I want to identify the clusters _without_ specifying the number of clusters a priori. What are appropriate clustering techniques? If my question does not contain enough information to answer, what are the pieces of information that I need to supply to determine appropriate clustering techniques? Below is an illustration of the kind of data/clustering I am imagining: ![clustering across time](http://i.stack.imgur.com/6ipJx.png)


title: model comparison when alternatives are not all nested within one another

text: I am running a glmm with three fixed effects: opponent 1 size ("1")   opponent 2 size ("2")   opponent 1 size - opponent 2 size ("diff")   I am unable to run all three variables in the model at once because of the "diff" variable being correlated with the "1" and "2" variables. How, then, do I decide which is the best model, when all the different combinations I can test are not nested within one another? Or can I consider the first two variables to be nested within the third one, since they are both used to calculate "diff"? I have these combinations of variables and their corresponding AICs (variable(s) followed by AIC in parentheses): variable(s) (AIC):   diff (223)   1 (231)   2 (262)   1,2 (265)   1,diff (265)   2,diff (265)


title: Calculation of probability of selling a Poisson number of products

text: You sell 10 product on average every day. The model is based on a Poisson distribution. To be enable to cover the indirect cost, the manager needs to sell 8 products per day. Two questions: >   1. What is the probability to sell between 0 - 5 products in a day? > >   2. What is the probability to have at least 6 products in a single day? > >


title: MCMC sampling from a function with a conditional distribution as Dirac delta function

text: Suppose the following distribution $p(x,y)$, $p:[0,1]^2 \to [0,+\infty]$ and: $p(x,y) = \left\\{ \begin{array}{ll} 0 & \mbox{if } y>0 \\\ +\infty & \mbox{if } y = 0 \end{array} \right.$ i.e. we have: $$p(y|x) =\delta_0$$ and $$p(x)=1$$ If this is the distribution that I want to sample from using MCMC methods, what do I have to do when I get $+\infty$ as the value of $p(x,y)$? * * * I am not sure but I realised in such a case the sample from joint this distribution is gonna be obvious! $y=0$ in all sample pairs. And $x$ uniformly chosen from $[0,1]$. I would be happy if somebody could confirm that


title: Sample size analysis for a study

text: I am a clinical researcher and orthopedic surgeon and one of my PhD student is starting on the last part of her thesis. Her work is centered around the potential benefits of interprofessional teaching at an orthopedic ward. The study is a retrospective cohort study on register data for hip fracture patients. Our hypothesis is that an interprofessional ward (with medical students and nurse students guided by real doctors and nurses) does not put patient at higher risk for adverse events than a control ward. As the primary outcome variable we have a proxy variable for adverse events; readmission rate at 3 months. The number of patients treated at the interprofessional ward is about 1:4 compared to the control ward. My question: Assuming a 3-month readmission rate of about 20% (about normal for this patient group) in the control ward, how large a sample size do I need when I aim at detecting a 5% difference in readmission rates? This is for 80% power and $p \le 0.05$. When performing this calculation as a superiority calculation for proportions in for instance SamplePower 3, I wind up at approximately 500 and 2000 in the interprofessional and control ward, respectively. Is this calculation correct, or should I do a non-inferiority analysis for proportions instead where the minimum clinical difference I am interested in is 5%? I cant find any such non-inferiority calculators where the number of subjects in the two groups are of different numbers, as in this study.


title: What is the best way to analyze data if your experimental design changes while running the experiment?

text: I imagine this is a somewhat common situation in practice. I am thinking mainly in terms of preclinical drug trials. 1) During the course of the study a new technique is learned or some time/money is freed up to be able to use a new technique as a secondary outcome. 2) Some key piece of equipment breaks, the lab member with the skill to perform a technique leaves, or a planned technique is looking like it is yielding inconsistent results for unknown reasons so that outcome is dropped for financial reasons. 3) When time comes to analyze data it comes to the knowledge of the researcher that the planned methods are inadequate/inappropriate and they wish to use a different approach. 4) Experiment 1 (e.g. test drug A) does not look like it will pan out, so a financial decision is made to stop that experiment and instead test drug B. Alternatively some other part of the experiment may be modified such as cell culture or animal strain. 5) Study is stopped prematurely or additional funds are freed up to increase sample size. I can think of more, but you get the idea. Are p values of any use under these circumstances? Should corrections for multiple comparisons be made in the case of #4? At what point does deviation from the plan invalidate the hypothesis testing procedure? The best way to me would be to simply plot the data and describe it in the hopes it may be useful for someone. However this behaviour is discouraged in favor of performing significance/hypothesis tests and the data may not be published if these are not included. **EDIT:** I was thinking about it and the actual best thing to do is to take the data and come up with a model to explain it that makes a precise prediction that can then be tested. It seems like this should always be the best thing to do, though.


title: Predict Failure time/ Weibull analysis

text: Suppose I want to predict how many machine are going to fail in the next three years. The data collected are in days, so we want to predict no. of fails in next 1095 days. All the machines (100 of them) started from, say, 1975 and until today, only 13 have failed. Now, I want to predict how many are going to fail in the next three years. Some of the issues I encountered, in my opinion:   1. Ideally, we have the fixed time to run the tests on these machines, BUT in this case, we don't have the fixed time. I can assume today's date as the final date, that is, my testing time interval is [1975, Today's date]   2. So, I have to put all the other machines that didn't fail as Censored observation.    3. BUT when I do that, the data does not follow, I believe, Weibull distribution, a standard process to analyze failure times (Weibull analysis).    4. Do we have enough data to come up with an appropriate model? And if this is what we have to work with, how can we do that?   5. Can we still use Weibull analysis approach for this case, and estimate the scale and shape parameter and use it to find the expected no. of failures given by:                     N(t)=(1-e^(t/a)^b)n      where _a_ is the scale parameter and _b_ is the shape parameter; _n_ is the population of machines which began operating continuously at the same time t=0 * * * Here are some of the arbitrary failure times ( I cannot share the data due to confidentiality issues):               10162, 8300, 11110, 11520, 11520, 8460, 7320, 11424, 11112, 11321, 11584, 10436, 9560      These are given in days from, say Jan 1, 1975. These are the machines that are built to last, so it is not uncommon to see these kind of minimal failure rates.


title: Physical significance of heavy tail distributions

text: I have come across several usage of the term "heavy tail distributions", but unable to find some good resources in the Internet to answer few of my questions:   1. Compared to what the tail of the distribution is said to be "heavy"?   2. Given some data, and their CCDF plots (a log-log plot, say) how do we say its heavy? Simply by curve-fitting with a power-law CCDF?   3. What is the physical significance (or effect) of the tail being heavy? For example, say the inter-arrival times follow a heavy tail distribution -- what are its implications?


title: ARIMA modeling: can seasonal data be seasonally stationary?

text: I am new to ARIMA modeling and currently encountering a weird situation with time series of count data. The time plot shows clear seasonal patterns.ACF also hints on presence of seasonality. However, seasonal unit root test in R shows that series is seasonally stationary. If I include seasonal differencing (D=1), I cannot find a single model where residuals satisfy normality assumption, even if I perform log- or - square root- or Box-Cox transformations of original series. If I do not include seasonal differencing (keeping seasonal AR and MA parameters is in the model), I easily identify a model with great diagnostics (residuals are white noise and normally distributed). Having hard time solving the puzzle whether seasonal-looking data can be seasonally stationary. Will appreciate any suggestion.


title: graphical and statistical approach to evaluate the skewness and kurtosis of a data set

text: In practice, we may always be asked to check the skewness and kurtosis of a data set. I have two questions. Given a probability distribution, how can we determine/evaluate the skewness and kurtosis of this data set? If we have the data sample itself, which kind of statistics can help us evaluate its skewness and kurtosis?


title: Regressing latent variables on observed variable using `sem` package in R

text: I'm having difficulty with some syntax in R using the `sem` package. I'm trying to do the following:   1. Create four latent variables with three (attitudinal) indicators each   2. Regress these latent variables on an observed (behavioral) variable   3. Compare this to a model where I regress four latent variables with three indicators each on a higher-order latent variable Using the `sem` package in R, my code for doing steps 1-2 is:               mydata.cov<-cov(mydata)     model.mydata <- specify.model()      F1 ->  X1, lam1, NA     F1 ->  X2, lam2, NA      F1 ->  X3, lam3, NA     F2 ->  X4, lam4, NA      F2 ->  X5, lam5, NA     F2 ->  X6, lam6, NA     F3 ->  X7, lam7, NA     F3 ->  X8, lam8, NA     F3 ->  X9, lam9, NA      F4 ->  X10, lam10, NA      F4 ->  X11, lam11, NA     F4 ->  X12, lam12, NA     behav -> F1, lam13, NA     behav -> F2, lam14, NA     behav -> F3, lam15, NA     behav -> F4, lam16, NA     X1 <-> X1, e1,   NA      X2 <-> X2, e2,   NA      X3 <-> X3, e3,   NA      X4 <-> X4, e4,   NA      X5 <-> X5, e5,   NA      X6 <-> X6, e6,   NA      X7 <-> X7, e7,   NA      X8 <-> X8, e8,   NA      X9 <-> X9, e9,   NA      X10 <-> X10, e10,   NA      X11 <-> X11, e11,   NA      X12 <-> X12, e12,   NA      F1 <-> F1, e13, NA      F2 <-> F2, e14, NA      F3 <-> F3, e15, NA      F4 <-> F4, e16,  NA      behav <-> behav, NA, 1      Presumably, the second model would be similar code:               model2.mydata <- specify.model()      F1 ->  X1, lam1, NA     F1 ->  X2, lam2, NA      F1 ->  X3, lam3, NA     F2 ->  X4, lam4, NA      F2 ->  X5, lam5, NA     F2 ->  X6, lam6, NA     F3 ->  X7, lam7, NA     F3 ->  X8, lam8, NA     F3 ->  X9, lam9, NA      F4 ->  X10, lam10, NA      F4 ->  X11, lam11, NA     F4 ->  X12, lam12, NA     F5 -> F1, lam13, NA     F5 -> F2, lam14, NA     F5 -> F3, lam15, NA     F5 -> F4, lam16, NA     X1 <-> X1, e1,   NA      X2 <-> X2, e2,   NA      X3 <-> X3, e3,   NA      X4 <-> X4, e4,   NA      X5 <-> X5, e5,   NA      X6 <-> X6, e6,   NA      X7 <-> X7, e7,   NA      X8 <-> X8, e8,   NA      X9 <-> X9, e9,   NA      X10 <-> X10, e10,   NA      X11 <-> X11, e11,   NA      X12 <-> X12, e12,   NA      F1 <-> F1, e13, NA      F2 <-> F2, e14, NA      F3 <-> F3, e15, NA      F4 <-> F4, e16,  NA      F5 <-> F5, NA, 1          model.sem2 <- sem(model2.mydata, mydata.cov, nrow(labor_coef))     # print results (fit indices, parameters, hypothesis tests)      summary(model2.sem)     stdCoef(model2.sem)      Unfortunately, neither syntax appears to work. Both models should be identified, and R just gives me a generic error message when I run the summary command: Updated: Now, the model will run, but summary returns the following: `Error in summary.objectiveML(data.sem) : coefficient covariances cannot be computed`. Thoughts?


title: How to reduce number of items using factor analysis, internal consistency, and item response theory in conjunction?

text: I am in the process of empirically developing a questionnaire and I will be using arbitrary numbers in this example to illustrate. For context, I am developing a psychological questionnaire aimed at assessing thought patterns commonly identified in individuals who have anxiety disorders. An item could look like "I _need to check the oven repeatedly because I can't be sure its off_ ". I have 20 questions (5-point Likert) which may be comprised of one or two factors (note that in reality I have closer to 200 questions, comprised of 10 scales, and each scale may be comprised of two factors). I am willing to erase about half the items, leaving 10 questions on one of two factors. I am familiar with exploratory factor analysis (EFA), internal consistency (Cronbach's alpha), and item characteristic curves in item response theory (IRT). I can see how I would use any single of these methods to determine which items are the "worse" within any single scale. I appreciate that each method also answers different questions, although they may lead to similar results and I am not sure what "question" is most important. **Before we start, lets make sure I know what I am doing with each of these methods individually.**   * Using EFA, I would identify the number of factors, and remove the items that load the least (lets say <.30) on their respective factor or that cross-load substantially across factors.   * Using internal consistency, I would remove items that have the worse "alpha if item deleted". I could do so assuming one factor in my scale, or do it after an initial EFA to identify the number of factors and subsequently run my alpha for each factor.   * Using IRT, I would remove items that do not assess the factor of interest along their (5 Likert) response options. I would be eyeballing item characteristic curves. I would basically be looking for a line on a 45 degree angle going from option 1 on the Likert scale all the way up to 5 along the latent score. I could do so assuming one factor, or do it after an initial   EFA to identify number of factors, and subsequently run the curves for each factor. I am unsure which of these methods to use in order to best identify which items are the "worst". I use worst in a broad sense such that the item would be of detriment to the measure, either in terms of reliability or validity, both of which are equally important to me. Presumably I can use them in conjunction, but I am not sure how. If I was to go ahead with what I know now and give it my best shot I would do the following:   1. Do an EFA to identify number of factors. Also delete items with bad loadings on their respective factors, since I don't want items that load badly regardless of how they would do in other analyses.   2. Do IRT and remove bad items judged by that analysis as well, if any remain from the EFA.   3. Simply report Cronbach's Alpha and don't use that metric as a means to delete items. **Any general guidelines would be greatly appreciated!** Here is also a list of specific questions that you can perhaps answer:   1. What is the practical difference between removing items based on factor loadings and removing items based on Chronbach's alpha (assuming you use the same factor layout for both analyses)?    2. Which should I do first? Assuming I do EFA and IRT with one factor, and both identify different items that should be removed, which analysis should have priority?  I am not hard set on doing all of these analyses, although I will report Chronbach's alpha regardless. I feel like doing just IRT would leave something missing, and likewise for just EFA. Thank you!


title: Word association

text: I'm doing an experiment with moodcongruency (influence of mood on word retrieval). Prior to inducing a mood, the participants will memorize a list of words containing an equal amount of positive, negative or neural words. Before the actual experiment, I want to make sure that the words are indeed associated (by a sample) with a positive, negative or neutral mood. So I've summarized a list containing 37 words and asked 24 respondents to rate their association with either (1) negative, (2) neutral of (3) positive meaning (nominal). From these results, I'm looking to distil the most reliable words into a list for each association (pos, neg, neut). How would I perform this in SPSS?


title: Standard deviation in SERVQUAL

text: I'm doing a project on measuring customer satisfaction using SERVQUAL. How do you calculate the standard deviation of each statement? For example, if for one of the questions, there were: >   * 2 respondents as strongly disagree (SD) >   * 5 disagree (D) >   * 20 neutral (N) >   * 21 agree (A) >   * 2 strongly agree (SA) (sample being 50) > How would you calculate the standard deviation if the scale were `SD=1`, `D=2`, `N=3`, `A=4` and `SA=5`? Also, how do you calculate the standard deviation for one dimension on the whole: e.g., the standard deviation for one of the dimensions in SERVQUAL say, Tangibles which has 4 statements or questions, and the standard deviation for Reliability which has 5 questions, etc.?


title: Confusion about scatter matrices

text: I am learning about evaluating clustering outcome and am confused about the scatter matrices. Hoping to get some help here. The within-cluster scatter matrix $S_W$is defined as: $$ S_W=\sum _{ k=1 }^{ K }{ \sum _{ x\in { C }_{ k } }^{ }{ \left( x-{ \mu }_{ k } \right) { \left( x-{ \mu }_{ k } \right) }^{ T } } } $$ The between-cluster matrix $S_B$ is defined as: $$ S_B=\sum _{ k=1 }^{ K }{ N_k }{ \left( { \mu }_{ k }-{ \mu } \right) { \left( { \mu }_{ k }-{ \mu } \right) }^{ T } } $$ where $K$ is number of clusters, $x$ is a member in cluster $C_k$, $\mu_k$ is centroids of cluster $C_k$, $N_k$ is number of members in cluster $C_k$, $\mu$ is the mean of the whole dataset. My dataset has the form of $m$ by $d$, i.e. $m$ data points of $d$ dimensions (features). After clustering, each cluster $C_k$ has the form of $N_k$ by $d$. So, a point $x$ has the dimension of $1\times d$, likewise $\mu_k$ has the dimension of $1\times d$. And, $S_W$ and $S_B$ are scalar. Why are they not matrices? What exactly should I expect in the elements of the matrices? With that result, the use of $trace \left( S_W \right)$ and $trace \left( S_B \right)$ become irrelevant. I sure has not understood the subject correctly. I will appreciate any help here. **Update 1:** The scatter matrix for each cluster is given as: $$ S_k=\sum _{ x\in { C }_{ k } }^{ }{ \left( x-{ \mu }_{ k } \right) { \left( x-{ \mu }_{ k } \right) }^{ T } } $$ which gives me a scalar quantity for $x$ and $\mu_k$ of size $1 \times d$. What should $S_{ k(i,j) }$ be? Given a dataset of $n$ observations with $d$ variables/features (i.e. $n \times d$), what should be the size of $S_k$? To further clarify my problem, suppose one of the clusters has 2 members (rows) with 3 variables/features (columns): $$ C_k = \begin{bmatrix} 1 & 2 & 3 \\\ 4 & 5 & 6 \end{bmatrix} $$ Then $$ { x }_{ 1 }=\left[ \begin{matrix} 1 & 2 & 3 \end{matrix} \right] \\\ { x }_{ 2 }=\left[ \begin{matrix} 4 & 5 & 6 \end{matrix} \right] \\\ { \mu }_{ k }=mean\left( { x }_{ 1 },{ x }_{ 2 } \right) =\left[ \begin{matrix} 2.5 & 3.5 & 4.5 \end{matrix} \right] \\\ \left( { x }_{ 1 }-{ \mu }_{ k } \right) { \left( { x }_{ 1 }-{ \mu }_{ k } \right) }^{ T }=\left[ \begin{matrix} -1.5 & -1.5 & -1.5 \end{matrix} \right] \left[ \begin{matrix} -1.5 \\\ -1.5 \\\ -1.5 \end{matrix} \right] =6.75\\\ \left( { x }_{ 2 }-{ \mu }_{ k } \right) { \left( { x }_{ 2 }-{ \mu }_{ k } \right) }^{ T }=\left[ \begin{matrix} 1.5 & 1.5 & 1.5 \end{matrix} \right] \left[ \begin{matrix} 1.5 \\\ 1.5 \\\ 1.5 \end{matrix} \right] =6.75\\\ { S }_{ k }=\sum _{ x\in { C }_{ k } }^{ }{ \left( x-{ \mu }_{ k } \right) { \left( x-{ \mu }_{ k } \right) }^{ T } } =6.75+6.75=13.5 $$ So, $S_k$ is scalar and consequently $S_W$ will be scalar. Where did I go wrong in the above computation of scatter matrix for one cluster? **Update 2:** So, the data points should be column vectors. $$ { x }_{ 1 }=\left[ \begin{matrix} 1 \\\ 2 \\\ 3 \end{matrix} \right] \\\ { x }_{ 2 }=\left[ \begin{matrix} 4 \\\ 5 \\\ 6 \end{matrix} \right] \\\ { \mu }_{ k }=mean\left( { x }_{ 1 },{ x }_{ 2 } \right) =\left[ \begin{matrix} 2.5 \\\ 3.5 \\\ 4.5 \end{matrix} \right] \\\ \left( { x }_{ 1 }-{ \mu }_{ k } \right) { \left( { x }_{ 1 }-{ \mu }_{ k } \right) }^{ T }=\left[ \begin{matrix} -1.5 \\\ -1.5 \\\ -1.5 \end{matrix} \right] \left[ \begin{matrix} -1.5 & -1.5 & -1.5 \end{matrix} \right] =\left[ \begin{matrix} 2.25 & 2.25 & 2.25 \\\ 2.25 & 2.25 & 2.25 \\\ 2.25 & 2.25 & 2.25 \end{matrix} \right]\\\ \left( { x }_{ 2 }-{ \mu }_{ k } \right) { \left( { x }_{ 2 }-{ \mu }_{ k } \right) }^{ T }=\left[ \begin{matrix} 1.5 \\\ 1.5 \\\ 1.5 \end{matrix} \right] \left[ \begin{matrix} 1.5 & 1.5 & 1.5 \end{matrix} \right] =\left[ \begin{matrix} 2.25 & 2.25 & 2.25 \\\ 2.25 & 2.25 & 2.25 \\\ 2.25 & 2.25 & 2.25 \end{matrix} \right]\\\ { S }_{ k }=\sum _{ x\in { C }_{ k } }^{ }{ \left( x-{ \mu }_{ k } \right) { \left( x-{ \mu }_{ k } \right) }^{ T } } =\left[ \begin{matrix} 2.25 & 2.25 & 2.25 \\\ 2.25 & 2.25 & 2.25 \\\ 2.25 & 2.25 & 2.25 \end{matrix} \right]+\left[ \begin{matrix} 2.25 & 2.25 & 2.25 \\\ 2.25 & 2.25 & 2.25 \\\ 2.25 & 2.25 & 2.25 \end{matrix} \right]=\left[ \begin{matrix} 4.5 & 4.5 & 4.5 \\\ 4.5 & 4.5 & 4.5 \\\ 4.5 & 4.5 & 4.5 \end{matrix} \right] $$ And finally, $S_W=\sum _{ k=1 }^{ K } S_k$. Yes, this is matrix! Do I get it right this time? The size of $S_W$ is $d \times d$ ($d$ being dimension/no. of features of the data points). $trace(S_W)$ is then the sum-of-squared-error. **Update 3:** Using the approach given by @ttnphns, the data matrix for cluster k can be arranged in rows (while the equations above have data in columns): $$ X_k = \begin{bmatrix} x_1^T \\\ x_2^T \end{bmatrix} = \begin{bmatrix} 1 & 2 & 3 \\\ 4 & 5 & 6 \end{bmatrix} $$ Each columns has 2 elements, hence to center the matrix columns, the required center matrix is $$ C_2 = I(2) - \frac{1}{n}O(2) = \begin{bmatrix} 0.5 & -0.5 \\\ -0.5 & 0.5 \end{bmatrix} $$ where $I(2)$ is identify matrix of size 2, $O(2)$ is 2-by-2 matrix of all 1's. Center columns of matrix $X_k$, $$ X_k^c = C_2X_k = \begin{bmatrix} 0.5 & -0.5 \\\ -0.5 & 0.5 \end{bmatrix}\begin{bmatrix} 1 & 2 & 3 \\\ 4 & 5 & 6 \end{bmatrix}\\\=\begin{bmatrix} -1.5 & -1.5 & -1.5 \\\ 1.5 & 1.5 & 1.5 \end{bmatrix} $$ Do $X_k^{c~T}X_k^c$ gives $S_k$ $$ S_k=X_k^{c~T}X_k^c = \begin{bmatrix} -1.5 & 1.5 \\\ -1.5 & 1.5 \\\ -1.5 & 1.5 \end{bmatrix}\begin{bmatrix} -1.5 & -1.5 & -1.5 \\\ 1.5 & 1.5 & 1.5 \end{bmatrix}\\\=\begin{bmatrix} 4.5 & 4.5 & 4.5 \\\ 4.5 & 4.5 & 4.5 \\\ 4.5 & 4.5 & 4.5 \end{bmatrix} $$ This is the same as in Update 2.


title: How to interpret coefficients in a vector autoregressive model?

text: Can I interpret the coefficients in a VAR model in the same way as I do in a normal OLS regression?


title: Is it valid to apply a t-test on scores obtained by cross validation?

text: I have a 12% holdout cross validation and I've done a bunch of repetitions which gives me a distribution of scores for each method. I'd like to do some sort of hypothesis testing for comparing methods. Each score I have is unitary value between 0 and 1. A view at a histogram yields an approximately normal distribution. Could a t-test be valid in this scenario?


title: Metric warning using caret's rfe

text: I am using the `caret` package to do feature selection with `rfe` while training a `knn` classifier. I want to tune both the `k` parameter and the `size` of the variable subset by taking the best value of the metric `ROC` from the `twoClassSummary()` function. My code to accomplish that is:               caretFuncs$summary <- twoClassSummary     knn.rfeC <- rfeControl(functions = caretFuncs, method = "repeatedcv", number=10, repeats=5,                            verbose = TRUE, returnResamp = "final")     knn.grid <- expand.grid(.k=seq(1, 20, 2))     knn.trainC <- trainControl(method = "boot", number=25, verboseIter=TRUE, returnResamp="final",                            classProbs=T, summaryFunction = twoClassSummary)     sizes <- 2:(ncol(train.x)-1)     set.seed(96)     knnR <- rfe(train.x, train.y, sizes = sizes, rfeControl = knn.rfeC, method="knn",                tuneGrid=knn.grid, metric="ROC", trControl=knn.trainC)      However, when `rfe` finishes, I get the following warning message:               1: In train.default(x, y, ...) :     The metric "Accuracy" was not in the result set. Sens will be used instead.      which I gather it means that the `metric` argument of `rfe` is only working on the subset size selection, but not the `k` parameter optimization, even though the `twoClassSummary` function is passed to the `summaryFunction` argument of `trainControl`. Am I missing something here? I would appreciate any help you could provide. EDIT: The output for `sessionInfo()` is:               > sessionInfo()     R version 2.13.1 (2011-07-08)     Platform: i386-pc-mingw32/i386 (32-bit)          locale:     [1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United Kingdom.1252        [3] LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C                                [5] LC_TIME=English_United Kingdom.1252              attached base packages:     [1] stats     graphics  grDevices utils     datasets  methods   base               other attached packages:     [1] caret_4.92      cluster_1.14.0  reshape_0.8.4   plyr_1.5.2      lattice_0.19-30          loaded via a namespace (and not attached):     [1] grid_2.13.1  tools_2.13.1


title: Interpretation of the regression coefficient of a proportion II

text: Following question: Interpretation of the regression coefficient of a proportion type independent variable In my model I have a log dependet variable. As indenpendet variables I have one proportion $X_1$ that lies within [0.001, 0.30] and a second proportion $X_2$ that lies within [0.12, 0.99]. $E(\log(y))$ = $\alpha + \beta_1X_1$ + $\beta_2X_2$ How do I interpret the coefficients? Let´s say, the estimated $\beta_1$ is 1.5 , the estimated $\beta_1$ is 0.8. Or do I have to adjust the range before running the regression?


title: Correct use of AIC

text: It is well known that the AIC can be used to compare nested models. Additionally, I believe I am correct in saying that you can also use the AIC to compare non-nested models on the same dataset (please correct me if I am in fact wrong). However, it is not correct to use the AIC to compare between different data sets. In my scenario, I have 5 measurements on individuals over time plus an outcome variable. If I was to regress the outcome variable linearly on all the measurements over time, I would be able to obtain an AIC for this model which uses the entire dataset. Now I want to consider only 2 of the measurements of all individuals plus the outcome variable. Technically, I am now using a subset of the original dataset as I have lost information on the other three measurements. However, isn't this the same as fitting a nested model since I kept all individuals but lost 3 explanatory variables? So is it justified to compare the AIC I get from this model with that obtained from the full model?


title: Simulate forecast sample paths from tbats model

text: Using the excellent forecast package by Rob Hyndman, I came across the need to not only have prediction intervals, but to simulate a number of future paths, given past observations of a time series with complex seasonalities. There is something for less complex time series with one or two seasonalities only (simulate.ets() in the forecast package), but in my case, I would require the equivalent of simulate.ets() for the more complex tbats model. I assume that the data necessary for creating such paths is already present in the fit object, yet the possibility to create sample paths seems to be not directly accessible. Therefore, I have come up with a naive solution and would like to know, whether this approach is correct.               require(forecast)     fit = bats(test,use.parallel=T,use.damped.trend=F,use.trend=T,seasonal.periods=seasonal.periods)      Naively, I imagine that sample paths can be constructed by using the point forecast from               fit           > forecast(fit)              Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95     1960.016       24.48576 23.82518 25.14633 23.47550 25.49602     1960.032       24.79870 23.88004 25.71735 23.39374 26.20366     1960.048       25.31743 24.39878 26.23608 23.91247 26.72239     1960.065       25.69254 24.77389 26.61120 24.28759 27.09750      1960.081       26.06863 25.14998 26.98729 24.66367 27.47359     1960.097       26.43215 25.51350 27.35080 25.02719 27.83711     1960.113       26.77674 25.85809 27.69540 25.37179 28.18170      and simply adding randomly drawn values from the model fitting procedure.               > fit$errors     Time Series:     Start = c(1959, 2)      End = c(1960, 1)      Frequency = 365        [1]  0.140656913 -0.455335141 -0.558989185  1.697532911 -0.114406022  0.366182718 -0.377056927  0.396144296      Therfore, with               prediction = forecast(fit)     errors = fit$errors          path = prediction$mean + sample(errors, size = length(prediction$mean))     plot(ts(path))      one sample path can be constructed. http://i.stack.imgur.com/mO83W.png Is this a valid way of constructing sample paths? If not, what would be a correct way? Thanks a lot for any help!


title: Seasonality criterion?

text: I'm working on a model of saving account amounts. It shows a seasonality: On the 1st day of every month the amount always decreases. This is reasonable as people probably need draw money to pay rental or other everyday bills. However, before seasonality adjustment is applied, I'm wondering: Is there a criterion of seasonality? For example, is there some “ratio” or “index” that could suggest the existence or obvious character of seasonality?


title: Excel xy chart with unequal x values in series

text: I am a senior doing a science fair project on the efficiency of fuels. For each fuel I tested, I have a time series of temperature. I want to plot these into an excel xy chart, how do I plot them if the constants (time) are not the same for each set of data?


title: Combining probabilities of nuclear accidents

text: The recent events in Japan have made me think about the following. Nuclear Plants are usually designed to limit risk of serious accidents to a 'design basis probability' for example, say, 10E-6/year. This is the criteria for a single plant. However, when there is a population of hundreds of reactors, how do we combine the individual probabilities of a serious accident ? I know I could probably research this myself but having found this site I a sure there is someone that will be able to answer this question quite easily. Thanks


title: categorical feature ranking

text: I would like to rank categorical features by the order or importance in a classification/regression setting. **Input** There are two features, which are survey questions:   1. "how is your mood?": four possible answers are crappy, mediocre, good and great!   2. "are you happy with your life?": three possible answers are no, maybe and yes. **Output** I am trying to predict the likelihood of a person being satisfied in a relationship, which can either be a probability on [0,1] scale (regression), or a set of 5 categories: extremely dissatisfied, dissatisfied, don't care, happy, extremely happy. **Problem** I am trying to rank which question is more appropriate in predicting how satisfied the person is in a relationship (if at all), and score the relevancy of the question on a categorical scale: not useful, somewhat useful, sorta useful, useful, extremely useful. The user can answer 0, 1 or both questions (so, sometimes the user votes on their relationship happiness and I don't have any other information on the user). **Approaches Tried So Far** I tried the following approaches:   1. LASSO: turned the questions and answers into a binary design dummy matrix, and mapped cases with 0 questions answered to the intercept. So, in the case of no answered questions, the row of the design matrix looks like [1 0 0 0 0 0 0 0] (4 zeros for first question and 3 zeros for second question). I then try to infer the likelihood of a question being relevant using a scheme on the non-zero LASSO regression coefficients (which is a bit hacky) and then bin the likelihood into each of the relevancy categories.   2. Partykit::ctree. Fit a decision stump, see which question the conditional inference framework picks first, exclude the data for that question (drop feature), regress onto the next question. If at any stage the decision stump is just a node, then stop the algorithm. Relevancy is ranked by p-value.   3. Similar idea to 2, but using CHAID multiway split trees.   4. Forward/backward feature selection, rank question relevancy by the magnitude of the cross-validation error.  **Question** Is there an algorithm or another approach which is better-suited for what I'm trying to do?


title: Can you use $R^2$ and regression to estimate cause and effect?

text: Does $R^2 * \mbox{slope}$ provide an estimation of cause and effect? For example, if $R^2$ of miles driven and coffee consumed is $x\%$ and the slope of the coffee-and-miles dataset is $y$ (miles/coffee), can you say "$z$ of each cup of coffee drunk is related to the number of miles driven"?


title: Is the product of multivariate lognormal distributions is multivariate lognormal distribution

text: A. Mulitivariate normal distribution case with latent variable $X$ Mapping from the low-dimensional space $X$ in Q-dimensional space (Q=2) to the high-dimensional space of $Y$ in D-dimensional space (say D=10) can be considered as $$Y=WX+\mu+\epsilon$$ Distribution of $Y$ given $X$ taking isotropic noise $\epsilon$ approximately as $N(0,\sigma^2I)$ can be defined as $$P(Y|X)=(2\pi\sigma^2)^{-\frac{D}{2}} exp(-\frac{1}{2\sigma^2}||Y-WX-\mu||^2),$$ Were $W$ is taken as $D$ by $Q$ weight matrix. If I take the Gaussian Prior over the latent variable $X$ with zero mean and unit variance $$P(X)=(2\pi)^{-Q/2}exp(-\frac{1}{2}X^TX)$$ Then the marginal distribution over $Y$ can be define as $$P(Y)=\int P(Y|X)P(X)dX$$ $$=(2\pi)^{-\frac{D}{2}} |C|^{-\frac{1}{2}}exp(-\frac{1}{2}(Y-\mu)^TC^{-1}(Y-\mu)),$$ where $C=\sigma^2I+WW^T$ B. Mulitivariate normal distribution case with latent variable $X$ Now here comes my question if I define distribution of $P(Y|X)$ as log normal when taking isotropic noise $\epsilon$ approximately $LogN(0,\sigma^2I)$ $$P(Y|X)=(2\pi\sigma^2)^{-\frac{D}{2}}(\prod_{i=1}^{D}Y_{i}^{-1}) exp(-\frac{1}{2\sigma^2}||ln(Y)-WX-\mu||^2),$$ Where $W$ is taken as $D$ by $Q$ weight matrix. If I take the lognormal Prior over the latent variable $X$ with zero mean and unit variance $$P(X)=(2\pi)^{-Q/2}(\prod_{i=1}^{Q}X_{i}^{-1})exp(-\frac{1}{2}ln(X)^Tln(X))$$ Then the marginal distribution over $Y$ can be define as $$P(Y)=\int P(Y|X)P(X)dX$$ $$=(2\pi)^{-\frac{D}{2}} (\prod_{i=1}^{D}Y_{i}^{-1})|C|^{-\frac{1}{2}}exp(-\frac{1}{2}(ln(Y)-\mu)^TC^{-1}(ln(Y)-\mu)),$$ where $C=\sigma^2I+WW^T$ Questions: I have two parts of my question: The first "is the product of two multivariate lognormal distribution is also multivariate lognormal distribution (as seen in equation of $p(Y)$ in the lognormal case)" If it is how can I prove this property? and the second part of the question is "how can I perform marginalisation over $Y$ in case of lognormal case and will the marginalising out the latent variable $X$ will also lead to the lognormal distribution" and how can I prove this property Assuming the same properties to that of normal distribution I have re-written the lognormal distribution is the expression for the PDF of $Y$ written above is true.


title: Generalized method of moments versus standard least squares estimation

text: I was thinking that in a very standard case such as a simple linear model with iid errors and no endogeneity, I would get the same results using the a simple least square estimate (such as provided by `lm()` in R) and a simple GMM estimator with an identity matrix as the weighting matrix and the matrix of regressors as instruments. I've just run the following simulation in R to compare the results of the `lm()` function and of the `gmm()` function :               > library("gmm")     > set.seed(1234567)     > N  <- 1000     > dd  <- data.frame(id = 1:N)     > dd$u  <- rnorm(N)         > dd$x  <- 1 + rnorm(N)     > dd$y  <- 1 + dd$x + dd$u     > m1  <- lm(y ~ x, data = dd)     > m2  <- gmm(y ~ x, x = ~ x,  wmatrix = "ident", data = dd)      I've got the same coefficients but I don't have exactly the same standard errors :               > coefficients(m1)     (Intercept)           x        1.0273856   0.9690455      > coefficients(m2)     (Intercept)           x        1.0273856   0.9690455      > sqrt(diag(vcov(m1)))     (Intercept)           x        0.044285    0.031726      > sqrt(diag(vcov(m2)))     (Intercept)           x       0.04367438  0.03127408       I'm not sure to understand why the standard errors are different in this case. Is this due to statistical theory or to the `gmm()` function ? See my gist file if you want to run the example.


title: R - Fitting two shapes by count of close values

text: I have two sets of data:               > dput(a)     c(15984.24, 17359.7, 17341.47, 18461.18, 24924.5, 46039.49, 41595.55,      19403.54, 21827.9, 21043.39, 22136.64, 26229.21, 57258.43, 42960.91,      17596.96, 16145.35, 16555.11, 17413.94, 18926.74, 14773.04)     > dput(b)     c(1.16776468870609, 1.25440938833205, 1.37597370494326, 1.4006623420762,      1.92223911868668, 2.8959286303691, 2.89756460878735, 1.53090508809756,      1.6460907811141, 1.68787497981851, 1.8009924919931, 1.88472999563212,      2.7976160336472, 3.16461975637253, 2.17047399047576, 2.30305266750383,      1.86211958973053, 1.37116485177612, 1.43838336037203, 1.33379850762213     )      And I would like to fit a and b such that I have the highest count of close value (by a threshold of 2% for example). So doing it manually, the fit would be something like:               > plot(a)     > lines(b*12000)      ![a and b superimposed](http://i.stack.imgur.com/TpERW.png) I wrote a function to count the close points:               CountClosePoints = function(x,y) {              count = 0         for(i in(1:20))         {             if( (!is.na(x[i])) & (!is.na(y[i])) )             {                 error = abs( x[i] - y[i] )                 threshold = 2*x[i]/100 # 2%                 if (error < threshold) count = count + 1 # one more point close to the target shape             }         }              return (count)     }      Now I would like to find K to maximize: CountClosePoints(a,K*b) for example: CountClosePoints(a,13000*b) = 4 CountClosePoints(a,12500*b) = 5 etc I do not know how to find K.


title: Does this de-noising algorithm have a name?

text: I recently was asked to help characterize some surfaces (raised ink dots on paper) that were scanned with a profilometer, yielding height data across a grid of x & y coordinates. The critical question was to characterize the average dot height. The data were rather noisy as the paper was not perfectly flat, so I played around with various approaches to removing this noise and came upon an approach that appears to do very well (to the eye at least). See below for pseudocode and R code. The logic is that every row and column definitely has an area with no dots on it, so ideally the minimum for each row and column should be around zero. However, one obtains different surfaces depending on whether row minimums are shifted to zero or whether column minimums are shifted to zero. To solve this ambiguity, I simply do both and average the results. I then noticed that repeating this process many times appeared to converge on a very nicely de-noised final result. Is there a name for this approach or something similar?               #Pseudocode ("data" is a RxC matrix of height values)     repeat the following 1000 times:         create a copy of data, call it "copy1"         for every row of height data in copy1:             subtract the minimum of that row from each value in the row (i.e. shift the row's data so that the minimum is zero)         create a second copy of data, call it "copy2"         for every column of height data in copy2:             subtract the minimum of that column from each value in the column (i.e. shift the column's data so that the minimum is zero)         replace data with the cell-wise average of copy1 and copy2      _               #R-code (data is in "a")     for(i in 1:1e3){         c1 = a         for(i in 1:ncol(c1)){             c1[,i] = c1[,i]-min(c1[,i])         }         c2 = a         for(i in 1:nrow(c2)){             c2[i,] = c2[i,]-min(c2[i,])         }         a = (c1+c2)/2     }


title: How can I determine the probability that a sample has a certain average given details about the total population?

text: Imagine that I have goal scoring data for a hockey team. They average 2.43 goals scored per game with a standard deviation of 1.63 for this data set. They have played 82 games. I have two subsets of this data set; each subset represents data for the team when a specific goalie is on the ice (Goalie A and Goalie B). Each goalie has played 41 games in no particular order. When Goalie A has been on the ice, the team averages 2.92 goals per game. When Goalie B has been on the ice, the team averages 2.00 goals per game. Assume that the goalies have no influence on their own team's goal scoring rates. What I am trying to understand is this: is there a way to compute the probability that a team could have data subsets that differ so greatly from one another? How likely is it that the team scores almost a full goal more with one goalie on the ice as compared to the other? If anyone could set me on the right path for how a calculation like this might unfold, that would be great. Thank you.


title: How to analyze healing time when the follow up times are not equal and we are looking at a specific time?

text: I would really appreciate if you could help me with your expertise: I have a study in which I evaluate patients for healing after discharge. So I have two group, healed and non-healed patients. The question I want to ask is how to analyze independent predictors of healing (or non healing) within a defined time interval like within two months of the insult. The time for follow exams is variable but most follow ups were in the period of two months of the insult with some extending up to 6 months of insult. I thought of using a cox model, where the censored data will be healed as non- healed patients which were examined after the period of interest were non- healed before that for sure but this is not the case for the healed ones. So wouldn't the Cox be reflective for my period of interest since most of the data lies within this period? Thank you very much.


title: Extract knowledge from the (possible unrelated) combination of multiple data sets

text: I am a beginner in this area so please bear with me. I have 2 data sets:   1. information about the cars which were sold in particular area (cost, model, number of seats, manufacturer, etc.)   2. information about social status of people in particular area (occupation position, average salary, education, health insurance, etc.) Is there some math that can show me some information e.g. how is the model of car related to occupation position of person (I am assuming that richer people buy more expensive cars)? Or even more crazy idea of how is the social status of person related to number of family members (I am assuming that people who buy cars with more seats (van etc.) have larger families)? Is this even possible with statistics/machine learning/data mining? In database world I would simply use selects and joins and get desired information, but what do do in this area? Can you point me to some resources with examples?


title: R - How to smooth a conversion of weekly sales into daily sales?

text: I have weekly sales figures, and would like to convert them into daily sales figures, making a simple hypothesis that there are 7 days with equal sales "power". Let's imagine that I have:               Week, Sales     1,11     2,15     3,9      Now I can easily convert this into daily figures with "flat" weekly patterns:               Week,Day,Sales     1,1,11/7     1,2,11/7     1,3,11/7     1,4,11/7     1,5,11/7     1,6,11/7     1,7,11/7     2,1,15/7     2,2,15/7     2,3,15/7     2,4,15/7     etc      But the main improvement I am looking for is to have a linear trend inside the week, so that I smooth the jump from week 1 to week 2 (from 11 to 15). Is there an R function that can do this work for me?


title: Does one lose information or precision when calculating effect sizes from statistical values rather than from the raw data?

text: Let's say that I want to compute the effect size of a certain intervention, for which I have both the raw data and some statistical value that can be transformed to an effect size. For example, if one has the t score and the df value of a t-test, one can transform it into Cohen's d by using the equation $$d=\frac{t*2}{\sqrt{df}}$$ This can also be done for other statistical values, like for example f statistics and Chi square. Do I lose any kind of information or precision if I chose to make the conversation from the statistical value rather than from the raw data? That is, will these transformations result in exactly the same effect size as if I would make the computations with the raw data instead, or are some of them merely estimates?


title: Pseudoinverse for time series filtering

text: Suppose I have a few time series, that are strongly correlated to each other, and I consider this to be noisy data. How do I use the pseudoinverse technique from image processing to filter this data?


title: Value Iteration exact convergence

text: I'm running Gauss Seidel value iteration on a weapon-target assignment problem (I can explain this more later if necessary, but i don't think it is). My VI is converging _exactly_ in 2 iterations. I'm not sure if it is because my problem has ~1/2 absorbing states (i.e. if you expend all the weapons or if you destroy the primary target you no longer transition, nor do you obtain a further reward). Being VI, I don't think it should be converging exactly (meaning that the inf-norm(V_2 - V_1) = 0). I've hand-computed a couple of the states and things seem to be working properly, but I can't figure out theoretically why it would happen. Matlab code posted below.               lambda = 0.95;     epsilon = 0.001;     tol = epsilon*(1-lambda)/2*lambda;     MaxDelta = 1000;          %%%%%%% Gauss-Seidel Value iteration %%%%%%%%          %Initialization     VS = zeros(cardS,1);     VS_up = zeros(cardS,2);     argmax = VS;     count = 1;     Vsa_T = cell(cardS,1);          while MaxDelta > tol          for sindex = Sindex         for aindex = 1:Aindex(sindex)             Vsa_T{sindex}(aindex) = r{sindex,aindex} + lambda*(Pij{sindex,aindex}'*VS);         end         [VS(sindex), argmax(sindex)] = max(Vsa_T{sindex});     end          VS_up(:,2) = VS;     MaxDelta = norm(VS_up(:,1) - VS_up(:,2),inf);     MD(count) = MaxDelta;     VS_up(:,1) = VS_up(:,2);     count=count+1     end


title: Rule of Thumb for minimum length of time series for AR(1) estimation

text: I have a data set of 350 points, I want to estimate the lag 1 auto correlation for different sub-sets of the data. More precisely I want to take non overlapping windows of length 1,2,3....n and calculate the auto correlations of these sub-sets. So for n equal to 2 I take every second observation similarly for n=3 I take every third observation and so on. Just wondered roughly what the maximum value of n I can realistically work with? Even with n=10 I'm down to 35 observations. Kind Regards Baz


title: Data transformation

text: I was writing with a question regarding a time-varying state space model of the form: \begin{align} y(t) &= \mu_1(t) + A(t)x(t) + v(t); &v(t) &\sim (0, R(t)) \\\ x(t) &= \mu_2(t) + \Phi(t)\\!\times\\!x(t-1) + s(t); &s(t) &\sim (0, Q(t)) \end{align} $x(t)$ is of the form: $x(t) = \log[1-X(t)]$ obviously, $$ |X(t)| \leq 1 \\\ \Rightarrow x(t) \text{ lies in } (-\infty, \log(2) ] $$ Could anyone suggest a transformation that I could use, so that $x(t)$ becomes unbounded?


title: Paired Samples T Test Suitability

text: I have data for 110 companies re. the price they have paid for a service annually between 2004 and 2011. Now I would like to find out if there is a statistically significant difference between the price they are charged pre and post 2008 . . . how should I deal with this, given that I have 8 observations per company (Can I pool them for pre and post 2008 and use the paired samples t test?)


title: R gls() vs. SAS proc mixed with interaction: Why does R complain about a singular matrix when SAS does not?

text: I like to keep analyses all in SAS or all in R when I can help it and lately have been using R more and more, but there's one analysis that I do somewhat routinely that has given me trouble in R. I have repeated measures data where I would like to fit the following model: $$Delta = Day + Group + Day\times Group$$ where $Delta$ is the change from baseline, $Day$ is the number of days from the beginning of the study, and $Group$ is the experimental group. I fit a variance-covariance matrix to account for repeated measures (for this example I'm using compound symmetry, but the difference is the same using other variance-covariance matrices). I have the data at the end of the post. If I don't include the interaction, I can get the analysis to run as I want it to in both SAS and R. In SAS:               proc mixed data=df;       class group day id;       model delta = day group;       repeated day / subject=id type=cs;       lsmeans group / diff=all;     run;      In R:               library(nlme)     library(lsmeans)     fit.cs <- gls(Delta~Day+Group,                   data=df,                   corr=corCompSymm(,form=~1|ID))     anova(fit.cs,type="marginal")     lsmeans(fit.cs,pairwise~Group)      Obviously the results differ in terms of denominator DF, but I don't intend to start that discussion (unless that difference is causing the problem). When I add interaction in SAS, everything is fine:               proc mixed data=df;       class group day id;       model delta = day | group;       repeated day / subject=id type=cs;     run;      But when I do the same from R...               fit.cs <- gls(Delta~Day*Group,                   data=df,                   corr=corCompSymm(,form=~1|ID))     # Error in glsEstimate(object, control = control) :      #   computed "gls" fit is singular, rank 19      Why does R complain about the fit being singular but SAS doesn't? Here are some fake data that are representative of data that I work with (from R):               df <- structure(list(Delta = c(-1.27, -0.34, 1.92, 0.45, 1.21, 0.43, -0.41, 0.16, -0.35,     1.49, -0.85, -0.86, 1.04, 0.49, 2.32, 0.13, -0.32, 0.5, 0.48, 1.21, -0.82, 0.93,     -0.58, 2.3, -0.9, 0.21, -0.72, 0.11, -0.28, -0.33, -0.7, -1.16, -0.23, -0.88, 0.97,     0.25, 0.8, 0.16, 0.63, -0.49, -0.63, -0.9, 1.1, -1.45, 0.38, -0.93, 0.4, 0.45, 0.48,     0.14, 1.02, -0.01, -1.98, 2.19, -1.53, -0.49, -1.57, -1.02, 1.09, 1.74, 0.54, -1.57,     -1.5, -0.48, 0.26, 0.2, -0.36, -1.05, -1.73, -0.77, -0.65, -1.07, -0.45, -0.14,     -0.56, 0.84, -2.66, -0.52, 1.44, 0.45, 0.24, -0.92), Day = structure(c(1L, 2L, 3L,     1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 6L, 7L, 1L, 2L, 3L, 6L, 7L, 1L, 2L, 3L, 6L, 7L,     1L, 2L, 3L, 6L, 7L, 1L, 2L, 3L, 6L, 1L, 2L, 3L, 6L, 7L, 1L, 2L, 3L, 6L, 7L, 1L, 2L,     3L, 6L, 7L, 1L, 2L, 3L, 6L, 7L, 1L, 2L, 3L, 6L, 7L, 1L, 2L, 3L, 6L, 1L, 2L, 3L, 6L,     7L, 1L, 2L, 3L, 6L, 7L, 1L, 2L, 3L, 6L, 7L, 1L, 2L, 3L, 6L, 7L), .Label = c("4",     "7", "10", "12", "14", "16", "28"), class = "factor"), Group = structure(c(1L, 1L, 1L,     1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,     2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,     3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,     4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L), .Label = c("1", "2",     "3", "4"), class = "factor"), ID = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L,     4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 8L,     8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 11L,     12L, 12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L, 13L, 14L, 14L, 14L, 14L, 15L, 15L, 15L,     15L, 15L, 16L, 16L, 16L, 16L, 16L, 17L, 17L, 17L, 17L, 17L, 18L, 18L, 18L, 18L, 18L),     .Label = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14",     "15", "16", "17", "18"), class = "factor")), .Names = c("Delta", "Day", "Group", "ID"),     class = "data.frame", row.names = c(NA, -82L))      Here's the same data for SAS:               DATA  df;     INPUT Delta Day Group $ ID $ @@;     CARDS;     -1.27 4 1 1 -0.34 7 1 1 1.92 10 1 1 0.45 4 1 2 1.21 7 1 2 0.43 10 1 2     -0.41 4 1 3 0.16 7 1 3 -0.35 10 1 3 1.49 4 2 4 -0.85 7 2 4 -0.86 10 2     4 1.04 16 2 4 0.49 28 2 4 2.32 4 2 5 0.13 7 2 5 -0.32 10 2 5 0.5 16 2 5     0.48 28 2 5 1.21 4 2 6 -0.82 7 2 6 0.93 10 2 6 -0.58 16 2 6 2.3 28 2 6     -0.9 4 2 7 0.21 7 2 7 -0.72 10 2 7 0.11 16 2 7 -0.28 28 2 7 -0.33 4 2 8     -0.7 7 2 8 -1.16 10 2 8 -0.23 16 2 8 -0.88 4 3 9 0.97 7 3 9 0.25 10 3 9     0.8 16 3 9 0.16 28 3 9 0.63 4 3 10 -0.49 7 3 10 -0.63 10 3 10 -0.9 16 3 10     1.1 28 3 10 -1.45 4 3 11 0.38 7 3 11 -0.93 10 3 11 0.4 16 3 11 0.45 28 3 11     0.48 4 3 12 0.14 7 3 12 1.02 10 3 12 -0.01 16 3 12 -1.98 28 3 12 2.19 4 3 13     -1.53 7 3 13 -0.49 10 3 13 -1.57 16 3 13 -1.02 28 3 13 1.09 4 4 14 1.74 7 4 14     0.54 10 4 14 -1.57 16 4 14 -1.5 4 4 15 -0.48 7 4 15 0.26 10 4 15 0.2 16 4 15     -0.36 28 4 15 -1.05 4 4 16 -1.73 7 4 16 -0.77 10 4 16 -0.65 16 4 16 -1.07 28 4 16     -0.45 4 4 17 -0.14 7 4 17 -0.56 10 4 17 0.84 16 4 17 -2.66 28 4 17 -0.52 4 4 18     1.44 7 4 18 0.45 10 4 18 0.24 16 4 18 -0.92 28 4 18     ;     RUN;


title: Transforming the data before the use of a maximum likelihood estimation

text: This might sound dumb but if I have $d_i$, where $i=1 \dots n$ observations and I assume they are exponentially distributed, before I use the MLE, should I transform my data to follow an exponential distribution? The reason I ask is because I often see others who normalize their dataset prior to the MLE if they assume their data is normally distributed.


title: lasso with missing values and categorical variables

text: I have a dataset with a lot of missing values and mix of continues and categorical variables. I want to use something like group lasso to do features selection. Probably the output is binary 0,1 and so grouped lasso logistic regression seems to be the more sensible choice. My problem is the very large number of missing values. Deleting non complete rows is not an option. Is there any R implementation that can be used similarly to the lasso and that can handle missing values and categorical variables at the same time? A possible solution has been proposed here but it does not refer to any R package.


title: Computation of variance based on the conditional variance

text: Can anyone give some advice on solving the following problem? $X$ is distributed as $N(\mu_1, \sigma_1^2)$; $Y \big| X > c$ is distributed as $N(\mu_2, \sigma_2^2)$; $Z \big|\min(X, Y) >c$ is distributed as $N(g(X,Y), \sigma_3^2)$. Here $\mu_1$, $\sigma_1$, $\mu_2$, $\sigma_2$, $\sigma_3$ are constants. How to find the variance of $Z$? Thank you very much.


title: When do you stratify an analysis versus including an interaction term?

text: I’m not very familiar with when and why you would stratify on a variable or set of variables in a regression analysis generally and would like to know what the issues are particularly in contrast to including the variable (by itself or as an interaction term) in the model without stratifying. I’m somewhat familiar with the example in Cox regression where you would stratify when the baseline hazard functions differ across levels of a variable, the variable is discrete and you don’t care if you have an estimate of it. In this case stratifying may be preferable to including an interaction with time term. But I’ve never really considered it in other regression contexts. What am I missing and what do I need to pay attention to? Also, I’m only talking about the analysis side of things, not in the design phase or sampling. I am only working with secondary analysis on observational data that hasn’t been taken from a stratified sample or design if that helps narrow the focus. And a specific question within this – would I really check each predictor in my model on whether to stratify? I understand with interactions that you wouldn’t necessarily check for interactions if you didn’t have an a priori reason to and they are not of primary interest. Thanks.