title: Combining summary statistics

text: Currently I add Var1+Var2=VarSUM12 and then perform a linear regression on VarSUM12~x to obtain test statistics for x and get insight into weak associations that are present across Var1 and Var2 with x. I would like to combine summary statistics (coefficient and standard error) from the individual regressions (Var1~x and Var2~x) to obtain the coefficient/se of VarSUM12~x. I guess the coefficient is easy, but is it possible to combine the standard errors?


title: Choosing a statistical test based on the outcome of another (e.g. normality)

text: So I have heard it said that it is not a good idea to choose one statistical test based on the outcome of another. This seems strange to me though. For example, people often choose to use a non parametric test when some other test suggests that the residuals are not normally distributed. This approach seems pretty widely accepted but does not seem to agree with the first sentence in this paragraph. I was just hoping to get clarification on this issue.


title: Statistical inference about degree of a node in a genetic network

text: I am working on Gene-Gene interaction networks. I build a graph by adding edges between genes (nodes) representing statistical interaction in prediction of a quantitative parameter value (say, brain volume) in a multiple regression model. For creating a well-connected graph, I have lowered the p-value threshold of the linear model interaction term in order to include as many edges as possible, so it is likely that a proportion of these links are false-positives (actually I have omitted multiple comparison correction and used a nominal p<0.05). My question is that, if some nodes show a very high degree centrality which is very unlikely to happen in a random Erdős–Rényi model (aka Poisson distribution), can I make a statistical inference that this node is biologically relevant to the brain size parameter because it has "received" more than random edges, even if the whole graph does not show a Erdős–Rényi degree distribution pattern? What is the best test for statistical evaluation of this hypothesis, and if there exists one, is it also necessary to perform a multiple-comparison?


title: How to measure the "average" autocorrelation of a time-series signal with itself

text: I have a short time series signal (say around 30 samples), and I would like to check whether or not it's oscillating. One approach I came up with was to measure the autocorrelation of the signal with itself. Given a random input signal, I would expect an output of 0. Given an oscillating signal I would expect some value consistently greater than 0. I have run the autocorrelation across all possible offsets of the signal (by "wrapping around" the index and assuming it's periodic), but when I take the mean of the autocorrelations across all possible offsets, I get (almost exactly) 0 no matter what the input signal looks like. This is likely "correct" (although not what I expected), but I can't seem to wrap my head around why that would be. Is there a better or more correct way of examining the autocorrelations to measure whether or not a signal tends to "agree with a shifted version of itself"? What am I missing?


title: Fitting a zero-inflated negative binomial regression with R

text: In this thread, I laid out a problem involving fitting a model that attempts to use minor league baseball statistics to predict success at the major league level (explained in full in the thread). After doing further research outside of the thread, I have come to the conclusion that a zero-inflated negative binomial model is likely the best fit given that I believe there are two processes generating the data. The first process determines whether a player will make the majors and once the player reaches the majors, a second process governs their success (as measured by WAR - also explained in the linked thread). I ran the model and the ZINB model appears to be a reasonable fit given the following diagnostic plot of fitted values vs residuals (broken into two plots to make it easier to inspect visually). ![image3](http://i.imgur.com/GO28KA3.jpg) ![image1](http://i.imgur.com/6J1QGOj.jpg) ![image2](http://i.imgur.com/Qe3vmGx.jpg) **EDIT 2: Here is a plot of the Pearson residuals.** ![image3](http://i.imgur.com/eyRBF7C.jpg/) **EDIT 3: Here is a plot of the Pearson residuals vs the fitted values.** ![image4](http://i.imgur.com/nNsgFPl.jpg) Although this is a big improvement over my previous models, there is clearly a bias for the model to underestimate a player's career WAR i.e., a majority of the residuals are greater than zero. **EDIT: It turns out that the plot is misleading due to the high number of overlapping residuals. In reality, only 10% of the residuals are greater than zero.** I am guessing that this may be improved by either a) a better specification of the functional form of the covariates or b) using a different yet similar model e.g., ZIP. Given that this type of regression goes well beyond what I have worked with before, I would appreciate any suggestions on further diagnostics I can use to both test this model and compare it to others and how to improve the functional form of the covariates given that the R function I have used, zeroinfl (from the pscl package), does not appear to be that flexible. Thank you!


title: Is the square root of a positive semi-definite matrix a unique result?

text: I am trying to decompose a time series of $n$ observations $\bf{\mathrm{v_c}}$ into the $n \times n$ variance-covariance structure $\sum$ and a random series $\bf{\mathrm{v}}$. So, I can derive the variance-covariance matrix $\sum$ from the autocorrelation function of $\bf{\mathrm{v_c}}$. This will be a Toeplitz matrix, which is positive semidefinite. Therefore, I am able to compute a suitable matrix $\sum^{-\frac{1}{2}}$ to transform my correlated series into a random signal. $\bf{\mathrm{v}} = \sum^{-\frac{1}{2}}\bf{\mathrm{v_c}}$ I am able to do this using the _sqrt(m)_ function in MATLAB, but can also find a Cholesky factorisation of the variance-covariance matrix and use this to induce the correlations. However, I get different (but somewhat similar) results for the random series using the sqrtm and Cholesky methods. I have read through several texts to determine how I might ascertain the square root of various matrices, and have looked at eigenvalue decomposition methods and so on. I see there are only unique solutions under certain prescribed conditions - but I presume that these unique solutions are still only one of many roots? My question is this: is there any way to argue that one particular square root is preferable over another. If not, is there a way to extract all possible solutions, such that all possible random functions may be obtained?


title: MA on a non-stationary time series

text: All, I have some data I would like to do some simple forecasting on. Its is non- stationary, looking at the time plot & from ADF & KPSS tests. After differencing I now have a stationary series. I want to use a moving average to forecast some future values. Is it recommended to do this on the stationary series and then use that forecast to apply back to my original one? Or can I go ahead and use on the original? Also, good reason to believe that more recent terms have more influence, should I use EWMA? many thanks!


title: Mixed model repeated measures, nested

text: I have a study in which I have 3 groups of crews (4 crews in each group) in which there are two people (L/R). Each crew experienced 4 lighting treatments (D1-4), once in one position, and again in another position (F/M). Each one of the crews had three things to view (N1-3). Just to complicate matters, it is probably more appropriate to say that the Ns were nested inside lighting treatments, as they were slightly different depending on the lighting treatment. D is within each crew. N is between crews (groups of). After each exposure, they are given a questionnaire with 5 items. My assumption is that the 4 crews within the same conditions should be considered a replication. So i have 4 reps of each Lighting/N/Position combination. Or, can i just consider each crew separately? **My first thought was to do**               L - within - fixed (4 levels)     N - between - fixed (3 levels)     Reps - between - random (4 levels)     Position - between - fixed (2 levels)      **Then, I thought perhaps i should consider that i have just 12 unique (L(4) x N(3)) treatments, rather than thinking they are two factors. But wasn't sure I could do that in a simple anova since one is between and the other is within. So it seemed a nested model was required**               N - between - fixed        N(L)-within- fixed       Rep - random       Position - fixed       N(L)*Rep        N(L)*Pos       N(L)*Rep*Pos      Is that right? And even if it is, ... I'm not sure I'd understand how to do that with repeated measures in spss. What would I lose if i just did it in GLM? My thought is to take the averages across all crews for the same N-group and position; so I would have essentially 6 'subjects' (defined by which N they got (3) x which position they are in (2)) that gave repeated measures of 5(questions)*4(L levels) = 20 columns of dependent variables; where each cell is an average of the 4crews' data. But that seems to mean you're missing some of the variability in the analysis that you'd be able to test in the Rep factor above.. yes? help .. help please ... thank you...


title: Evaluation of LMT (Logistic Model Tree) classifier results

text: I have been using the LMT Logistic Model Trees algorithm in some classification experiments. However, after reading the reference/documentation regarding the algorithm (Niels Landwehr, Mark Hall, Eibe Frank (2005). Logistic Model Trees. Machine Learning. 95(1-2):161-205.) and implementing it myself, I would still require some more clarification on the best way to evaluate the trees that are produced, especially in the manner of observing/evaluating how each leaf can contribute to a particular classification. Here is an example of a result: > ## Logistic model tree > >> : LM_1:12/12 (453) _// <\- here what would the 1:12/12 ? (as the 453 is number of instances)_ > > Number of Leaves : 1 > > Size of the Tree : 1 > > LM_1: Class 0 : 0.45 + > > [[cue_1]] * -72.14 + > > [[cue_2]] * -2.71 + > > [[cue_3]] * -3.88 + > > [[cue_4]] * 54.6 + > > [[cue_5]] * -12.28 + > > [[cue_12]] * -0.57 + > > [[cue_18]] * -7.24 > > Class 1 : -0.45 + > > [[cue_1]] * 72.14 + > > [[cue_2]] * 2.71 + > > [[cue_3]] * 3.88 + > > [[cue_4]] * -54.6 + > > [[cue_5]] * 12.28 + > > [[cue_12]] * 0.57 + > > [[cue_18]] * 7.24 In the leaf above, these are the cues that were taken into consideration from a list of 20. Thus, I suppose that these were the cues that had viable information, correct? From what I understand, each leaf is representative of a logistic regression model. In this sense, do the pluses and minuses that occur signify positive or negative weights toward that particular class? Does it signify the amount of weight that a particular cue (either positive or negative) provides toward classification toward one class or another? Any information regarding the topic would be extremely appreciated. Thank you, in advance, for any assistance you would be able to lend me.


title: Statistically fuzzy version of a checksum for file text signature

text: **Background:** Often I end up downloading the same pdf article twice since I do not remember I've already downloaded it. One way around is to maintain an index of cheksums (say md5 etc.) based on the plaintext conversion of a pdf and if a match is found do not re-download. Problem is online downloaded papers often have a timestamp, IP, or other embedded information in the footers or an extra beggining page from the downloading service. This precludes an exact checksum. Further if the pdf consists of scanned images an OCR pre-processing step is needed and the output might vary a tiny bit every time. Even the pdf->plaintext conversion utilities are not error free and might add some noise while converting formatting, equations, tables etc. **The Statistics:** What might be a good way to signature an article's text that is robust to such small noise? One option I can think of is, say, the frequency spectrum of the top 20 most common words in any given paper. But that is prone to get diluted by "the, a" and all the common words. What else might be good ideas for a fuzzy variant of a cheksum? By design for md5sum etc. the output changes drastically for even tiny changes in input. I don't want this property.


title: PCs scores from Correlation and Covariance matrices through matrix computations and prcomp

text: I'm want to get PCs scores through matrix approach. My calculated PCs scores for correlation matrix matches with `prcomp` results but the PCs scores for covariance matrix do not match with the results of `prcomp`. Could you point out what am I missing? Thanks **PCA on Correlation matrix**               # PCA on Correlation matrix     X <- USArrests     Cor <- cor(X)          EigenCor <- eigen(Cor)     ECor <- EigenCor$vectors     head(t(t(ECor) %*% t(scale(X))))                                [,1]       [,2]        [,3]         [,4]     Alabama    -0.9756604  1.1220012 -0.43980366  0.154696581     Alaska     -1.9305379  1.0624269  2.01950027 -0.434175454     Arizona    -1.7454429 -0.7384595  0.05423025 -0.826264240     Arkansas    0.1399989  1.1085423  0.11342217 -0.180973554     California -2.4986128 -1.5274267  0.59254100 -0.338559240     Colorado   -1.4993407 -0.9776297  1.08400162  0.001450164          PCACor <- prcomp(x = X, retx = TRUE, center = TRUE, scale. = TRUE)     summary(PCACor)               Importance of components:                               PC1    PC2     PC3     PC4     Standard deviation     1.5749 0.9949 0.59713 0.41645     Proportion of Variance 0.6201 0.2474 0.08914 0.04336     Cumulative Proportion  0.6201 0.8675 0.95664 1.00000          head(PCACor$x)                       PC1        PC2         PC3          PC4     Alabama    -0.9756604  1.1220012 -0.43980366  0.154696581     Alaska     -1.9305379  1.0624269  2.01950027 -0.434175454     Arizona    -1.7454429 -0.7384595  0.05423025 -0.826264240     Arkansas    0.1399989  1.1085423  0.11342217 -0.180973554     California -2.4986128 -1.5274267  0.59254100 -0.338559240     Colorado   -1.4993407 -0.9776297  1.08400162  0.001450164      **PCA on Covariance matrix**               # PCA on Covariance matrix     Cov <- var(X)     EigenCov <- eigen(Cov)     ECov <- EigenCov$vectors     head(t(t(ECov) %*% t(X)))                     [,1]      [,2]      [,3]       [,4]     Alabama    -239.7035 -46.45394 -5.873077  5.7840485     Alaska     -267.7288 -39.91901 16.748431 -0.7178995     Arizona    -298.9695 -66.73235 -5.065592 -0.9775376     Arkansas   -193.2414 -41.19804 -3.167955  2.8551540     California -282.3243 -80.42202  3.367729  0.5643217     Colorado   -209.8773 -71.62154  8.901219  1.6546839                    PCACov <- prcomp(x = X, retx = TRUE, center = TRUE, scale. = FALSE)     summary(PCACov)                    Importance of components:                                PC1      PC2    PC3     PC4     Standard deviation     83.7324 14.21240 6.4894 2.48279     Proportion of Variance  0.9655  0.02782 0.0058 0.00085     Cumulative Proportion   0.9655  0.99335 0.9991 1.00000          head(PCACov$x)                      PC1        PC2        PC3        PC4     Alabama     64.80216 -11.448007 -2.4949328 -2.4079009     Alaska      92.82745 -17.982943 20.1265749  4.0940470     Arizona    124.06822   8.830403 -1.6874484  4.3536852     Arkansas    18.34004 -16.703911  0.2101894  0.5209936     California 107.42295  22.520070  6.7458730  2.8118259     Colorado    34.97599  13.719584 12.2793628  1.7214637


title: Is it ever appropriate to put dependent variables in a model when they calculate the independent variable directly?

text: Is it ever appropriate to fit a multivariable regression model with dependent variables that directly calculate the independent variable? For example, I know that fitting a model to predict BMI by using height and weight and other variables in the model is a No-No since BMI=weight/height^2, however... what if I have a longitudinal data set where individuals have BMI measurements reported over time and I wanted to use a higher order term using height and weight (like time*height and time*weight and time*height*weight)? Is it still inappropriate to include the lower order term of height and weight in the model? Is it even appropriate to use the higher order terms that include height and weight? I appreciate the feedback.


title: What is conditioning in spatial statistics?

text: Could someone explain to me: the concept of " **conditioning** " in spatial statistics in a fairly advanced context?   Here is an example to clarify the question:   **Step 1)** generate a 2D point process, here 6 realization are shown: ![enter image description here](http://i.stack.imgur.com/4t2Pr.png) **Step 2)** choose a region from one realization ![enter image description here](http://i.stack.imgur.com/vSfBe.png) **Step 3)** try to fill the rest of same size region, here I used the same PP: ![enter image description here](http://i.stack.imgur.com/x0eZW.png) **Step 4)** remove the overlapped region and merge with given part: ![enter image description here](http://i.stack.imgur.com/LHJyn.png) **Step 5)** assume you did right! ![enter image description here](http://i.stack.imgur.com/ZX6Yd.png) **Step 6)** ask question now: **Is the result is a conditioned PP on a set of given points?** **Step 7)** wait keeping hope that somebody will answer your question thoroughly:)


title: Inter-feature ratio explicitly or implicitly?

text: Let's say I have two numerical features were is suspect that the ratio between them is the most meaningful way of looking at them. I have a NN learner. Should I add the ratio as a third feature or is this just idiotic? In other words: Would it make sense to do a feature selection with all three of them and see if one or two of them suffice? Or is this just a waste of time because it doesn't matter what feature of the three of them I remove?


title: Distribution of volatility

text: Is there a test/way to determine the distribution of stochastic volatility? For example, I have a random walk where the increments are non-normally distributed and heteroskedastic. I would like to try a model which includes stochastic volatility but have no reason to believe it is normal - and don't most models assume this?


title: Generative modeling of a mix of continous and discrete variables

text: I'm trying to build a generative model to run a Monte Carlo simulation. The existing data consist of a combination of discrete and continuous variables. Suppose I have a number of people...               Age Sex Non-white     21  1   1     35  1   1      I can easily use a mixture of Gaussians to model this data set and just use EM to estimate the mixture coefficients. But I'm not sure if the underlying assumptions are sound when the underlying components of the joint aren't necessarily (or at all) Gaussian. I know a couple of more complicated methods, such as using Bayesian networks or things coming out of machine learning (e.g. Boltzmann Machines), but I doubt that they'll be useful for my rather small data set. I am wondering if there's a compact way to build generative models of multivariate, correlated categorical variables. Or is mixture of Gaussians generally sufficient?


title: Simple introduction to MCMC with Dirichlet process prior?

text: I'm looking for a simple and easy to read introduction to using MCMC with a dirichlet process prior. Or perhaps using MCMC in any machine learning scenario, eg Gaussian Process. I've been circling around various papers and tutorials (Neal, Teh, Sahu, Ghahramani, Ferguson, Escobar and West, ... ), but I can't find anything that explains in straightforward procedural terms, specifying explicitly each equation I would need to implement MCMC for DP or GP. Edit: I got this far: my mcmc dp progress


title: Can statistics help in inventory management decisions?

text: I'm in the process of calculating inventory levels for a bunch of SKUs and I'm wondering if you all-mighty statisticians had any card up your sleeves that I could borrow. Here's the deal: service levels are inevitably going to be affected by the amount of material I decide to put in stock. The more I stockpile, the higher the fill-rate will be (but of course inventory doesn't come cheap). The less I keep on-hand, the higher the chance of having a stock-out. I have to strike a balance between holding costs, and service level. Based on past sales, I was able to plot the expected fill-rate against the volume of on-hand materials. Check it out: ![enter image description here](http://i.stack.imgur.com/wVMjE.png) As you can all see, the service level increases very rapidly at first, but then the growth rate starts decreasing quite abruptly. Going from 90% to 100% requires 9 times the amount of material it took to go from 0% to 90%. Now, getting to the point, is there any scientific way to tell where I should stop? How do I know when the marginal increase in fill-rate is too little for that much of an increase in stocked volume? I know choosing the right service level is mostly a "political" decision, but I really hoped that you guys could come up with some kind of suggestion to help me in the process. I'm kinda out of my element here. Thank you guys for your time!


title: Time-series data evaluation by PCA / EOF

text: I have got time series data of a crop attribute, each time step with 500 data points on a spatial grid. 2 measurements were conducted in the first year, 3 in the consecutive year. I want to evaluate the variability of the measured crop attribute by topographic and soil parameters. This is simple for each time step, but to achieve an overall pattern, I think that PCA/EOF is the right choice (right?). Basically I am following the process stated on: Performing a principal component analysis in R First I manually scaled the data by `w_ph_s <- (w_ph - 0.87) / 0.11`, where 0.87 is the **mean of all data of one year** and 0.11 the SD of one year. Then, I used `pca_w <-princomp (w_ph_s)`, which revealed 3 components, where the first reflects 82% of the variability. This example is for the second year with 3 measurements. Should I use the scores of the first component, which reflects most of the variability as the input data for correlation with the stated parameters? Do the loadings than reflect the time fluctuation at the specific point? Is the scaling of the data right? Can I do the same for both years with adjusted values?


title: Bernoulli Confidence Intervals for p very close to 0

text: Let's say I have the following observations from many Bernoulli distributions with different p (p1, p2, ..): Observations from Distribution 1: 10 successes, 100,000 trials, p_hat = 0.0001 Observations from Distribution 2: 0 successes, 100 trials, p_hat = 0 Observations from Distribution 3: 4 successes, 60,000 trials, p_hat = 0.00007 I want to order these distributions by their true probabilities of success and get rid of the ones that have low probability of success. However, because of the inherent nature of these distributions, the probability of success is so low, that if I use a standard Wald and Wilson confidence interval for Bernoulli distributions, the results don't make too much sense. Is there a standard statistical way to deal with these types of problems? Or do I have to resort to some self defined heuristics to remove distributions with low probability of success?


title: Multinomial confidence interval

text: I am trying to help out a quality auditing department and figure out sample sizes and confidence intervals and such. The sample size and confidence level (95%) is pretty well documented on the web and I am muddling through it. I even figured out how to do the exact confidence level of a binomial sample which I would use for a pass/fail condition but I wanted to figure out how to do a multinomial of that same sample so that I could generate some type of information about a situation with three possibilities: pass / fail(critical) / fail(non-critical). Can someone point me in the right direction here? Ultimately I am going to want to program this into an application so a formula would get me going, a link to a way to program it would be a home run.


title: Forecast accuracy measures for different forecast horizon h in R

text: I have a yearly time series data, from 1980 to 2005. The data is splitted into an training sample and a out of sample; the out-of sample consists of the 6 most recent observations and the rest is considered for training sample. I need to fit a ETS model and compare different accuracy measures for different forecast step aheads h=1,2,3,4,5 and 6. Something like this:                   h=1  h=2  h=3  h=4  h=5   h=6      MSE .. .. .. .. .. .... .. .. .. .. .... .. .. .. .. .. MASE .. .. .. .. .. ...... .. .. .. .. .... .. .. .. .. .. The following code gives me the accuracy measures for h=6:               trainx<- window(x,end=1999.99)     testx<- window(x,start=2000)     fit<- ets(trainx)     accuracy(forecast(fit,h=6),testx)      The questions are:   1. How can I calculate the accuracy measures for h=1,2,3,4,5 ? For instance, when h=2, I fit a model to training data and I produce the forecast that correspond to 2000 and 2001.    2. Now, how should I produce the forecast for 2002 and 2003, etc?   Should I suppose that the observations for the year 2000 an 2001 are known and then fit a new model (this time I need to add the observations of 2000 and 2001 to the training set), then, to produce the forecast for 2002 and 2003?


title: Extremal serial dependence

text: As part of my analysis of heavy-tailed time series of company returns, I would like to check whether extreme returns exhibit serial dependence, i.e. if extreme events are followed by extreme events. For now, I am only interested in testing uni-variate time series, i.e. I just want to check if an extreme return of company A has an impact on the returns of company A in the following days. Since the autocorrelation function (ACF) of a time is looking at the whole return series (not only the extremes), it has limited value in assessing dependence dependence between extreme events. In my search so far, I have found recent papers about the "extremogram" (http://arxiv.org/abs/1001.1821 and http://www.uis.no/getfile.php/Konferanser/paper_RichardDavis.pdf), something they view as the extreme-value analog of the ACF of a stationary process. They also write that the extremogram is similar to the tail dependence coefficient of a uni-variate time series. While looking for more information about the tail dependence coefficient, I think tail dependence coefficients are mostly used in the analysis of co-movements in the tails of two distributions - which I am not interested in (for now). Now finally my question:   1. Are there any other methods available for estimating extremal serial dependence?   2. Are there any implementations of the extremogram available in R/MATLAB? (Note: I have also contacted the authors and asked if they can share their code with me.)   3. If there is no implementation available for the extremogram, are there any other applications available? (e.g. how could I calculate the tail dependence coefficient of a univariate time series in R/MATLAB). Edit: Author of the extremogram papers replid to me and told me there is no "official" implementation of the extremogram in R/MATLAB yet, but he assured me that it is easy to implement. Will be working on my own implementation for the extremogram then.


title: Having trouble with nls function in R

text: I've tough luck with the use of nls() in _R_ for the following model $$N_e = N_o\\{1-exp[\frac{(d+bN_o)(T_h N_e - T)}{(1+c N_o)}]\\}$$ where $b>0$, $c\geq 0$, $T_h>0$, and $T=72$. This code               T <- 72     NLS.Fit3 <- nls(Ne~No*(1-exp((d+b*No)*(Th*Ne-T)/(1+c*No))), data = Data,                 start = list(d = 0.01, b = 0.01, Th = 0.01, c = 0.01),                 control = nls.control(maxiter=50, tol=1e-05, minFactor=1/1024))      gives the following error message: Error in nls(Ne ~ No * (1 - exp((d + b * No) * (Th * Ne - T)/(1 + c * : singular gradient And the following               NLS.Fit31 <- nls(Ne~No*(1-exp((d+b*No)*(Th*Ne-T)/(1+c*No))), data = Data,                 start = list(d = 0.01, b = 0.01, Th = 0.01, c = 0.01),                 control = nls.control(maxiter=50, tol=1e-05, minFactor=1/1024),                 algorithm = "port", lower=c(0, 0, 0, 0))     summary(NLS.Fit31)      code converges but provides the wrong results (drastically different from PROC NLIN)               Formula: Ne ~ No * (1 - exp((d + b * No) * (Th * Ne - T)/(1 + c * No)))          Parameters:     Estimate Std. Error t value Pr(>|t|)       d  0.008325   0.003488   2.387   0.0192 *     b  0.000000   0.000064   0.000   1.0000       Th 0.000000   0.614220   0.000   1.0000       c  0.020670   0.034439   0.600   0.5500       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 4.631 on 85 degrees of freedom          Algorithm "port", convergence message: relative convergence (4)      I'd prefer to do this in _R_ rather than _SAS_ and how the constrains can be placed on only few paramters. Any help in this regard will be highly appreciated. Thanks Data is here:               No  Ne     5   0     5   1     5   1     5   2     5   2     5   2     5   2     5   3     7   0     7   0     7   1     7   1     7   2     7   2     7   2     7   3     10  1     10  1     10  2     10  2     10  3     10  3     10  3     10  4     10  7     15  1     15  1     15  3     15  3     15  4     15  5     15  5     15  5     20  3     20  4     20  7     20  7     20  8     20  8     20  9     20  11     25  4     25  5     25  6     25  7     25  9     25  9     25  13     25  14     30  5     30  8     30  10     30  11     30  11     30  12     30  14     30  20     45  4     45  7     45  8     45  10     45  11     45  14     45  15     45  19     60  9     60  14     60  14     60  16     60  18     60  21     60  24     60  26     80  7     80  11     80  12     80  15     80  17     80  12     80  21     80  23     100 7     100 8     100 10     100 11     100 15     100 24     100 26     100 33


title: How to perform logistic regression with lasso using GLMSELECT?

text: I am trying to perform logistic regression with lasso. For the logistic regression part I am using `PROC LOGISTIC` but I am not sure how to do lasso with `PROC LOGISTIC`. I searched online and found that `PROC GLMSELECT` allows us to do lasso. But I am not sure how to do logistic regression with lasso using `PROC GLMSELECT`. **Note:** I posted this question in the SAS Discussion Forum.


title: Formula for recommended sample size for multivariate testing

text: Context: To recommend a minimum sample size when performing multivariate testing of a web page. The sample size would vary based on the number of factors being tested (e.g. A heading, and an image) and the number of variations of a factor (e.g. Two different headings, and three different images). The goal could be to see which combination caused the most purchases of a product. Based on a recommended minimum sample size of 100 for a single factor with perhaps two variations, I'm trying to work out a formula that recommends a sample size with multiple factors and variations. The formula I first came up with is shown below where n is the number of variations for that factor. $ samplesize = 100*\left((n_{f1}-1)^{n_{f1}-1}*(n_{f2}-1)^{n_{f2}-1}*\ldots\right)$ Does this seem reasonable and/or is there a simpler formula that would be as reasonable? The intended audience are online business owners who aren't necessarily strong at maths. Thanks for reading!


title: How to plot an ellipse from eigenvalues and eigenvectors in R?

text: Could someone come up with _R_ code to plot an ellipse from the eigenvalues and the eigenvectors of the following matrix $\begin{equation*} \mathbf{A} = \left( \begin{array} {cc} 2.2 & 0.4\\\ 0.4 & 2.8 \end{array} \right) \end{equation*}$ Thanks


title: Different non-parametric methods for estimating the probability distribution of data

text: I have some data and was trying to fit a smooth curve to it. However, I do not want to enforce too many prior beliefs or too strong pre-conceptions (except the ones implied by the rest of my question) on it, or any specific distributions. I just wanted to fit it with some smooth curve (or have a good estimate of the probability distribution that it might have come from). The only method that I know for doing this is kernel density estimation (KDE). I was wondering, if people knew of other methods for estimating such a thing. I just wanted a list of them and from that I can do my own research to find out which ones I want to use. Giving any links or good references (or intuitions on which ones are good) are always welcome (and are encouraged)!


title: Estimator as a Vector

text: Suppose we are performing a linear regression of $Y$ on $X$. That is: $$E[Y|X] = \beta_{0} + \beta_{1}X$$ Is it correct to say that an estimator is $\mathbf{\hat{e}} = (\hat{\beta_0}, \hat{\beta_1})$?


title: Detecting time-shifted time series

text: Assume there are two time series of real values. How to test hypothesis, that these series are quite the same, but there is a time shift between them? Sorry, if the question is too basic. Thanks in advance


title: How to support choice for hold out/forecast period/window size in forecasting stock returns

text: Are there any literature concerning what in-sample or out-sample period i should choose and how long the window should be when i consider a "moving window forecast"


title: Root finding for stochastic function

text: Suppose we have a function $f(x)$ that we can only observe through some noise. We can not compute $f(x)$ directly, only $f(x) + \eta$ where $\eta$ is some random noise. (In practice: I compute $f(x)$ using some Monte Carlo method.) What methods are available for finding roots of $f$, i.e. computing $x$ so that $f(x) = 0$? I am looking for methods which minimize the number of evaluations needed for $f(x)+\eta$, as this is computationally expensive. I am particularly interested in methods that generalize to multiple dimensions (i.e. solve $f(x,y) = 0, g(x,y) = 0$). I'm also interested in methods that can make use of some information about the variance of $\eta$, as an estimate of this may be available when computing $f(x)$ using MCMC.


title: Can I partition a frequency table and use Fisher's exact test instead of the chi-squared test?

text: I am working with very small cell counts, and am wondering if Fisher's exact test can be used in place of the chi-squared test when partitioning IxJ tables. I am interested if the same process for partitioning the tables is acceptable.


title: Is anything wrong with this approach to PCA?

text: I'm working on an implementation of PCA that works on very large data sets. Based on my understanding of the algorithm, the first step is to do an SVD of the input `m x n` matrix, `X`. This SVD looks like X = WΣVT. The "interesting" output `Y` of this process -- from Wikipedia, "The PCA transformation that preserves dimensionality (that is, gives the same number of principal components as original variables)" -- is given by the following equation: ![equation](http://i.imgur.com/orkKvMo.png) Based on my reading, if I can compute the `W` component of the SVD, then I can compute Y as: ![equation](http://i.imgur.com/itBvuNY.png) The upshot here is that I'd only have to compute `W`. In terms of computational and memory complexity, this approach is significantly more efficient because the only matrix above and beyond the initial data set I'd have to load has size `m x m`, which (at least in my case) is much smaller than `V`, which would be `n x n`. Is there some reason why this derivation won't work that I'm not seeing?


title: Different approaches of linear discriminant analysis

text: What are the differences in three different approaches of classification of LDA for two groups? The three are: 1\. Fisher's approach 2\. regression approach 3\. Bayes's approach They will be exactly the same when the two groups have same size. But when they are of different size, how three approaches differ?


title: Uncertainty and sensitivity analysis

text: I have the following problem: > Given inputs $x$ ($n$-dimensional vector) of scalars, ordered integers and > unordered integers (i.e., labels) and one or several outputs $y$, I would > like to estimate: > >   1. Which inputs explain the best the outputs. >   2. To what extent variations of one input imply variations of the outputs. > This is supposed to be related to uncertainty and sensitivity analysis, which are quite broad fields. Do you know of any methods/resources with an approach related to my problem?


title: Good summaries (reviews, books) on various applications of Markov chain Monte Carlo (MCMC)?

text: Are there any good summaries (reviews, books) on various applications of Markov chain Monte Carlo (MCMC)? I've seen Markov Chain Monte Carlo in Practice, but this books seems a bit old. Are there more update books on various applications of MCMC in areas such as machine leearning, computer vision, and computational biology?


title: How to do a multiple regression with ARIMA using R?

text: I am analyzing some tree physiology data (transpiration) in relation to a number of environmental variables (many of which are predictors such as temperature, PAR and vapour pressure deficit). I have fine-scale (30 min intervals) data of these various measurements, and there are two objectives I am trying to achieve:   1. Use the various predictors (glm?) to see which among these explain the most amount of variation in transpiration. However, since there is clear autocorrelation at this scale (i.e., trans at time $t$ is highly correlated with trans at $t+1$ etc.), I am looking to use ARIMA models with regressors.    2. I would like to construct a final predictive ARIMA model that explains the highest variation in trans, from all the different candidate models. So far, I have noticed that ccf plots show -ve lags between trans and a number of variables (rightly so, e.g., as you expect temp at time $t$ to influence transpiration at $t+1$). My questions are:   1. How do you perform an ARIMA with transpiration as the response variable and several regressors?    2. How do you know which one of the regressors to leave out? Does this have to be done manually in R (as in, add each regressor to the model, and inspect the resulting AIC)?    3. Is `auto.arima` the best way to determine the differencing term (etc.)? (E.g., `auto.arima(trans, xreg=temp+vpd+......)`.)   4. How do you account for the lag between response variable at time $t$ and predictors at $t-1$?


title: What is the difference between a Normal and a Gaussian Distribution

text: Is there a deep difference between a Normal and a Gaussian distribution, I've seen many papers using them without distinction, and I usually also refer to them as the same thing. However, my PI recently told me that a normal is the specific case of the Gaussian with mean=0 and std=1, which I also heard some time ago in another outlet, what is the consensus on this? According to Wikipedia, what they call the normal, is the standard normal distribution, while the Normal is a synonym for the Gaussian, but then again, I'm not sure about Wikipedia either. Thanks


title: What are some good resources for the history of time series analysis?

text: I have checked out the answer to this question on stats.stackexchange: What are good resources providing a history of statistics? Indeed, the Stigler book "Statistics on the Table" looks excellent and I'm looking forward to reading it. But I'm more interested in the development of modern ARIMA models. I think I remember hearing that much progress was stimulated in trying to predict random inaccuracies with artillery guns around WWII. Also, of course, to some extent astronomers throughout the second half of the millennium were using some sort of time series for understanding the motion of heavenly bodies. However, I can't remember where I heard about the artillery application of time series and I have a background in physics, and I don't really know what statistical methods the astronomers were working with. So, I would love to hear what you think are the most pertinent historical influences on the development of time series methods, e.g., were they stimulated mainly by finance, defense, geology/geophysics, or some combination of all these and more? Are there informative books or websites on the history of ARIMA?


title: Find the minimum of an expensive-to-sample noisy paraboloid

text: I have a function on the unit square that's expensive to evaluate, and want to estimate where it reaches its minimum in at most $N$ evaluations, $N$ ~ 10 say. Assuming a model of the form paraboloid + noise,   $ \qquad \text{f}( x; x_{min}, a, c, \sigma ) \approx a (x - x_{min})^2 + c + \mathcal{N}( 0, \sigma^2 ) $ what's an algorithm to generate points $x_0 \dots x_{N-1}$ to sample f() at, for a best estimate of $x_{min}$ ? (My real problem is 2d, but it may be helpful to do 1d first, so $ x\ x_{min}\ a\ c\ \sigma $ are 1d or 2d as appropriate.) After sampling f() at the 4 corners and the middle, one can alternate:   * fit a paraboloid to the data so far $\rightarrow x_{min}$   * generate the next sample point $x_i$ ... how ? Heuristics come to mind, e.g. "near $x_{min}$ but not too near", but this is surely a well-known problem ? (I imagine the cases $\sigma$ known and $\sigma$ unknown will be different.) * * * Added 11 Feb: In grid search, one typically evaluates f() 5 times at each point (5-fold cross-validation), on a say 10 x 10 grid, for 500 evaluations in all. The 5 values would give an estimate of $\sigma$ at each point, but then what -- how would a statistician combine the 5 x 100 values ? This question on SVM grid search shows 12 _different_ bumpy not-so- paraboloids, but does not show the spread at each point. Also, search grid+search on stats.stackexchange.


title: Where is the derivative for the transfer function in the delta rule?

text: First off: I understand derivatives and the chain rule. I'm not great with math, but I have an understanding. Numerous tutorials on backpropagation (let's use this and this) using gradient descent state that we use the derivative of the transfer function (sigmoid or tanh) to calculate the gradient and hence which way to travel next. In these tutorials I see $(t-o)(1-o)(o)$ as the formula for calculating the error for the output neurons, which seems to be the derivative of the error calculation $\frac{1}{2}(t-o)^2 \times (1-o)$. (t = target, o = output_actual) Why do I not see the derivative of the transfer function (assuming sigmoid): $\dfrac{e^x}{(e^x + 1)^2}$ anywhere? Or for when tanh is used as the transfer function: $\operatorname{sech}^2(x)$ ... where x = weighted input? Also, some tutorials use (target - actual), (target - actual)2 [ Sum of Squares - useful for negative outputs] or the squared error function: (1/2)(target - actual)2. Where is the derivative of the transfer function and which is the correct error formula to use?


title: Clustering of accelerometer events

text: Refer to the images below. I have an accelerometer attached to a door that logs events everytime someone opens and closes a door. I'm attempting to predict the individual who opened or closed the door based on how they exert force on the door during the event. In the images I've given examples of 2 separate people opening the door and you can see minor variations in their behaviour. I'm at this stage exploring a cluster/classification solution to identifying the individual related to each event. Any help on finding an approximate solution would be great! ![enter image description here](http://i.stack.imgur.com/KpBbD.png) ![enter image description here](http://i.stack.imgur.com/Eghwg.png) ![enter image description here](http://i.stack.imgur.com/BIMPD.png)


title: Difference in tf-idf values in R

text: I am playing around in R to find the `tf-idf` values. I have a set of `documents` like:               D1 = "The sky is blue."     D2 = "The sun is bright."     D3 = "The sun in the sky is bright."      I want to create a matrix like this:                  Docs      blue    bright       sky       sun        D1 tf-idf 0.0000000 tf-idf 0.0000000        D2 0.0000000 tf-idf 0.0000000 tf-idf        D3 0.0000000 tf-idf tf-idf tf-idf      So, my code in `R`:               library(tm)     docs <- c(D1 = "The sky is blue.", D2 = "The sun is bright.", D3 = "The sun in the sky is bright.")          dd <- Corpus(VectorSource(docs)) #Make a corpus object from a text vector     #Clean the text     dd <- tm_map(dd, stripWhitespace)     dd <- tm_map(dd, tolower)     dd <- tm_map(dd, removePunctuation)     dd <- tm_map(dd, removeWords, stopwords("english"))     dd <- tm_map(dd, stemDocument)     dd <- tm_map(dd, removeNumbers)     dtm <- DocumentTermMatrix(dd, control = list(weighting = weightTfIdf))     as.matrix(dtm)           Terms                 Docs      blue    bright       sky       sun                 D1 0.7924813 0.0000000 0.2924813 0.0000000                 D2 0.0000000 0.2924813 0.0000000 0.2924813                 D3 0.0000000 0.1949875 0.1949875 0.1949875      If I do a hand calculation then the matrix should be:                           Docs  blue      bright       sky       sun                 D1    0.2385    0.0000000  0.0880    0.0000000                 D2    0.0000000 0.0880     0.0000000 0.0880                 D3    0.0000000 0.058      0.058     0.058       I am calculating like say `blue` as `tf` = `1/2 = 0.5` and `idf` as `log(3/1) = 0.477121255`. Therefore `tf-idf = tf*idf = 0.5*0.477 = 0.2385`. In this way, I am calculating the other `tf-idf` values. Now, I am wondering, why I am getting different results in the matrix of hand calculation and in the matrix of R? Which gives the correct results? Am I doing something wrong in hand calculation or is there something wrong in my R code?


title: How to use posterior function within Matlab

text: I would like to create a GMM within Matlab and then input an observation x to get the probability of this x. I know the equation is the following: $$p(x) = \sum_{k}p(C_{k})p(x|C_{k})$$ I would like to get the value of $p(x|C_{k})$ first then calculate $p(x)$   I would like to ask that if the output of function "posterior" in matlab is equal to $p(x|C_{k})$?


title: How likely are various outcomes in a lottery with multiple prizes

text: This question is (I think) simpler than most posed here, but it's beyond my ability to solve. I'm trying to calculate the probability of various outcomes for a charter school lottery. There are two classes of entrants: students otherwise zoned for under-performing schools, and all other students. Let x=the number of students zoned for under-performing schools Let y=the number of students not zoned for under-performing schools Let z=the number of spaces available in the school It's easy enough to determine the expected number of students admitted from each class, but how can I determine the probability of various other outcomes. For instance, how likely is it that the number of students admitted from a class is half the expected number? For what it's worth, I'm a math teacher who has never had a stat, probablity, or discrete math course, although I've taught myself the basics. If you have any suggestions for further self-guided study in probability I'd appreciate them, but right now I just need to solve this practical problem.


title: Where can I download sample data sets for mobile app downloads?

text: For statistical testing, I am in need for android and apple mobile app download data. I have checked on AWS (Amazon) and Infochimps. No luck so far. Can anyone suggest other resources or is there a code that can be used in Visual Studio C# to scrape the data out Google API? e.g.               AppID, AppName, Number of Downloadds, Type of phone etc


title: apply unscented kalman filter

text: I have a time series of electricity spot price over a 2 years period. I want to create a model in order to predict the future values of the time series using unscented Kalman filter. I found out that I could use dual estimation algorithms. My questions are: 1\. Is there another (simpler) approach to solve the modeling problem? 2\. Can anyone provide me a detailed numerical example of applying UKF? Thank you for your time.


title: Creating a Multivariate time series or matrix containing the bottom level series

text: Am using the hts and gts for hierarchical forecasting in R. I have a 2 level structure which comprises with Y_t(My multivariate time series for the total energy loss) which in turn can be either be due to UCLF or OCLF and the total of 3 stations(Ankerlig,Anort and Camden). I want to construct Y_t= **S** *Y_k,t where **S** is my summing matrix.


title: MLE for simulated case of binomial p with constant labelling rate

text: I simulated the following example: 2000 binomial runs (p(Heads)=0.6). Each run has a sample size ranging from 500 to 300 0. (But that is not important: we have the same result if the sample size is the same). We imagine that there is some constant counting error among the tails, the heads are over-reported. For example, if in a run, there are N=2000 coins, the actual number of heads is 1205, but the reported number of heads is 1205 + (2000-1205)*0.10, 0.10 being the constant over-counting error. My question is : can MLE recover the p and the error rate (0.6 and 0.1) in this case? I have tried doing this, but the global MLE (by trial and error) is at 0.66 (with error rate=0). Can I set up my MLE, or use any other way to recover the 0.6 and 0.1 ?


title: Way to Measure Groupings Using Distances Between Individuals?

text: I am working on a problem that requires me to measure groupings of people. I have the location of every individual in my sample at every point in time. It's therefore trivial to calculate the distances between each individual for a certain point in time. Is there a way I can think about groupings using these pairs? Currently I am defining a group as two individuals $\le$ a certain distance away from one another. This is the way I'd like to define my groupings. Basically, I'm trying to find a way to measure "groupishness" of all these characters over a set time span. I know there are _n_ choose 2 pairings, and therefore I could calculate for every possible pair the average distance between them, but I'm trying to get a measure for the entire group. Any ideas would be helpful. Thanks.


title: Calculating whether a disease is probable using Bayes rule?

text: I want to compute whether it is more probable that a patient has a disease or the contrary. If I am given the following information: $$P(\mathrm{disease})= 0.008$$ $$P(+|\mathrm{disease})= 0.98$$ $$P(-|¬\mathrm{disease})= 0.97$$ To find the answer would it be correct to do the following: $$\frac{0.008 \times 0.98}{(0.008 \times 0.98)+(0.992 \times 0.03)}=0.2085$$ and then calculate $1-0.2085$ to see which is larger?


title: selecting the bins for extremely skewed data

text: I have a data that exhibits nearly a power law distribution, and I want to know a good binning technique to summarize the statistics. For example consider the following data: $$ \begin{array}{rr} \textbf{Object} & \textbf{Count}\\\ A & 82952\\\ B & 8832\\\ C & 1801\\\ D & 555\\\ E & 231\\\ F & 123\\\ G & 60\\\ H & 43\\\ I & 28\\\ J & 20\\\ K & 8\\\ L & 13\\\ M & 4\\\ N & 4\\\ O & 4\\\ P & 3 \end{array} $$ The count is the number of times "Object X" appears. As you can see from "Object J" and onwards, the frequencies of "Object K–P" are too low. So I want to combine "Object K-P" into one bin. The regular techniques like "equal width" and "equal frequency" binning do not work. I can manually write a program to do this, but is there any statistical discretization technique that can effectively bin such patterns?


title: MLE for distribution where some data not available

text: say we have a random sample of a random variable, but we are told that no information is available for values less than or equal to some number, perhaps k. That is, we don't have the full sample. We only get the values that are greater than k. So, how do we account for this in the likelihood function? Do we multiply by (1 - F(k)) or something like that? Any help is greatly appreciated!


title: What are the canonical data sets used for testing robust linear fitting?

text: The UCI database (link) is one of the repositories for canonical data. It has ~295 data sets for use. There are many others. (link) While data can be useful not all data is relevant for all problems. The Fischer Iris data could be considered canonical for a classification test, but is not necessarily good for line fitting. It could be hacked, subsetted, and abused to make a fit, but its best use is the classification/clustering. So what are the textbook datasets used to benchmark performance of robust linear estimators? If I wanted to (again) demonstrate to myself the improvement of bootstrapping the Theil-Sen estimator over its "textbook" implementation, which data-sets show exactly what the nature of the improvement is, and where the "cliff" in magnitude of improvement is located? I prefer open data, not closed. If you want to point to a paper as a reference and say "they used x data and it was good" that is fine as long as the data was open data. Thank you in advance for your help.


title: How to get an R-squared for a loess fit?

text: How to calculate R-squared ($r^2$) statistic in R for `loess` and/or `predict` function output? For example for this data:               cars.lo <- loess(dist ~ speed, cars)     cars.lp <- predict(cars.lo, data.frame(speed = seq(5, 30, 1)), se = TRUE)      `cars.lp` has two arrays `fit` for model and `se.fit` for standard error.


title: Life after the Box-Cox transformation

text: Suppose, we have a set of measurements of some quantity in some units of measurement. We also have a nice model that heavily relies on the properties of the Gaussian distribution. The model is tailored for data in some specific units of measurement with some physical meaning behind (like watt, ohm, etc.). It turns out that the distribution of the data does not exactly follow the normal distribution and has some undesired features (like skewness). We apply the popular Box-Cox transformation and obtain a more or less normally distributed data set. The problem now is that we have logarithms, powers, etc. of the original measurements, which contradicts with our nice model. The question is, what can one do in such a situation? I need to change the model such that it can handle the new data? And in general, if I got everything correctly, why do people what to study transformed data that have lost their physical meaning? Because, at the end of the day, one will, probably, have to return back to the original units of measurement. Thank you.


title: Results on Monte Carlo estimates produced by importance sampling

text: I have been working on importance sampling fairly closely for the past year and have a few open-ended questions that I was hoping to get some help with. My practical experience with importance sampling schemes has been that they can occasionally produce fantastic low-variance and low-bias estimates. More frequently, however, they tend to produce high-error estimates that have low sample variance but very high bias. I am wondering whether anyone can explain exactly what kinds of factors affect the validity of importance sampling estimates? In particular, I am wondering: 1) Are importance sampling estimates guaranteed to converge to the correct result when the biasing distribution has the same support as the original distribution? If so, why does this appear to take so long in practice? 2) Is there a quantifiable relationship between the the error in an estimate produced through importance sampling and the "quality" of biasing distribution (i.e. how much it matches the zero-variance distribution) 3) Partially based on 1) and 2) - is there a way to quantify 'how much' you have to know about a distribution before you were better off using an importance sampling design than a simple Monte Carlo method.


title: Choosing the right test when comparing two groups: Mann-Whitney, Kolmogorov-Smirnov or other?

text: I need help, since my skills and knowledge in statistics are really limited. Using Pearson correlation I found that there is some relationship between a binary one and a second variable taking values between 1 and 4. So now I want to understand if there are significant differences between the two groups defined by the binary variable. The 2 variables are the following:               X = (0 0 0 0 0 0 0 0 0 0 1 1 1 1)     Y = (4 3 4 3 3 3 3 3 3 3 4 4 4 4)      To determine if there's a significant difference between $Y_{X=0}$ and $Y_{X=1}$, I used a Mann-Whitney test. Was this the correct choice? Would it have been preferable to use a K-S test? And why not a t-test? I am actively looking for an answer to these questions.


title: What's the chance for one to die in a car accident in their life-time?

text: Considering the high number of car accidents, it must be around 5%, right?


title: Classification table for ordinal logistic regression in R

text: I want to create a classification table regarding an ordinal response variable with three levels but I don't know how to do it. Searching on the site I fell on the question posted by Brandon Bertelsen that covers only the case of the binary logistic regression (link at the end of the post).Does anyone knows how I can create such that table in my case? I don't know if it is important but I used the `rms` package to run the olr and using the `predict(fit,type="fitted.ind")` command I get the next table with probability for each case                     grade=1    grade=2   grade=3     1  0.08042197 0.28380601 0.6357720     2  0.08086877 0.28475584 0.6343754     3  0.41472656 0.40802584 0.1772476     4  0.39680650 0.41484517 0.1883483     5  0.25402385 0.43644283 0.3095333     6  0.13539881 0.37098177 0.4936194     7  0.12591996 0.35959459 0.5144855     8  0.50489952 0.36489760 0.1302029     9          NA         NA        NA     10 0.34757283 0.42969971 0.2227275     11 0.24690054 0.43539812 0.3177013     12 0.17325212 0.40529586 0.4214520     13 0.45795712 0.38900855 0.1530343     14 0.03594015 0.16033637 0.8037235     15         NA         NA        NA     16 0.50188652 0.36653955 0.1315739     17 0.48710163 0.37441720 0.1384812     18 0.38094725 0.42028884 0.1987639     19 0.04134659 0.17894428 0.7797091     20 0.12844729 0.36275605 0.5087967     21 0.23991274 0.43410413 0.3259831     22 0.20506362 0.42316514 0.3717712     23 0.45457929 0.39061326 0.1548075     24         NA         NA        NA     25 0.31269786 0.43606610 0.2512360     26 0.20905830 0.42483513 0.3661066     27 0.05240710 0.21353381 0.7340591     28 0.26569967 0.43759072 0.2967096     29 0.21258621 0.42621415 0.3611996     30 0.11407246 0.34347156 0.5424560     31 0.34656138 0.42993750 0.2235011     32 0.01813256 0.08978609 0.8920813     33 0.44034224 0.39716470 0.1624931     34 0.12213714 0.35468488 0.5231780     35 0.40888783 0.41032190 0.1807903     36 0.33901842 0.43161582 0.2293658     37 0.13275554 0.36793345 0.4993110     38 0.32091057 0.43492411 0.2441653     39 0.45984161 0.38810515 0.1520532     40 0.55550665 0.33564053 0.1088528     41 0.02812293 0.13122652 0.8406505     42 0.46250424 0.38681892 0.1506768     43 0.07352751 0.26852580 0.6579467     44 0.04330967 0.18541327 0.7712771     45 0.45457929 0.39061326 0.1548075      Logistic Regression: Classification Tables a la SPSS in R


title: Time series analysis on login data to forecast CPU demand using R

text: Motivation: I was hired as an intern a few weeks ago to figure out if my company needed to buy new machines six months in advance. Database machines take up to 4 months to install and there is a 2 month grace period. I signed an NDA, so I don't think I can give any actual data. The only reliable information I have now, is information on the number of logins and registrations for an education company from 2002 to 2011. I think I can get more recent information on registrations, and people are working on getting login information. We stopped logging login information in 2011 so there will be a gap of no data when I try to forecast :( The information is collected daily. I've created a time series forecast of the data using R. I used this tutorial http://a-little-book-of-r-for-time- series.readthedocs.org/en/latest/src/timeseries.html#arima-models To make a holt winters exponential model with daily frequency (frequency = 365). I've removed February 29 from the data. Unfortunately the gap in login data means I will have to try a more specific ARIMA right? Will I be able to use arima if there are long gaps in the data? Also, the arima function in R doesn't allow for frequencies greater than 350, and it runs out of memory quickly, so I'd have to use a monthly model (freq = 12). I have tried using fourier but the predictions didn't look right intuitively. Since I want to know what the peak usages are though, I think I might want to be more specific. Is it ok to use a weekly frequency (freq = 52) and just remove Dec 31? Is daily frequency allowable? Like can I use exponential smoothing with daily frequency even though Sept 7, 2012 might fall on a Sunday, whereas Sept 7, 2011 and 2010 and 2009 might all be weekdays. There is a daily, weekly, and yearly seasonality in demand and number of logins. Eg. 6pm, and Monday, and September are more loaded in general than 4am, and Saturday, and May. There is a yearly seasonality in number of registrations. I've been having some issues with the login predictions The problem is that variability increases too much before 6 months have even passed. At the 80% confidence interval. The projection line extends into 2012 and the orange area is the 80% confidence interval. Logging and using additive exponential smoothing gave me much more variability than multiplicative exponential smoothing. It's not useful to the company to say that "well you might have 8 jillion logins sometime in the next 6 months and you might have 20% more than you had last year." How do I reduce the variance in the projection? http://img836.imageshack.us/img836/8460/holtwintersloginmultipl.png Finally, I was thinking that after I got accurate projections, I'd put logins and registrations in a neural network, and I'd put something like average wait time on a few machines as the ouput variable, and I'd forecast peak projected processing power demand in 6 months. There are other variables to consider, like software releases that change cpu demand per user, but I'm hoping the neural network will learn when these happen, or that they are easy to detect and account for. I don't have any good data on average wait time yet, but assuming I find some, is this a good plan?


title: Sample size calculation for truncated normal distribution

text: So, I am new here. I need to perform a sample size calculation for a clinical trial. The study sample will be select according to criteria of person's height. Persons within the particular height range (female, 1.6m to 1.7m) will be invited to participate in trial. We know the expected sample standard deviation from previous trial. But my concern is that sample is not from a normal distribution. Usual power/sample size calculation need the assumption of normal distribution of test statistic under $H_0$ and $H_1$, but here I believe we have truncated normal distribution. So how may I modify `power.t.test`, or make some other calculation in R, to accommodate this? My colleague says to just rely on central limit theory and assume normal with 1.65 mean and known standard deviation, but I believe this is wrong due to the truncation. Any advice would be appreciated.


title: Statistical tests for low sample numbers in both indendent groups (n=3 and n=1)

text: I've been asked to analyze some data for an experiment that's already been done. 15 rats were raised in 5 groups of 3   1. No treatment   2. Saline treatment (control drug is delivered via saline)   3. Drug treatment   4. Saline shock (high dose)   5. Drug shock (high dose) Due to budget constraints, groups 2~5 were physically pooled (i.e. the three samples in each group were mixed and analyzed as a single sample) so as to fit onto a single iTRAQ 8-plex, whereas group 1 was not pooled so I have 8 channels   1. No treatment   2. No treatment   3. No treatment   4. Saline treatment (control drug is delivered via saline)   5. Drug treatment   6. Saline shock (high dose)   7. Drug shock (high dose)   8.  For those of you who don't know what iTRAQ is, it's a mass spectrometry protein identification and quantification method, which generates a list identified proteins and an 'absolute' quantification value. I say absolute, but it is meaningless in an absolute sense, and only useful when compared to another quantification value in the same experimental run. So for a given protein X I can calculate the ratio of X in channel 5 to the ratio of X in channel 1, or even the ratio of X in channel 5 to X in group 1 (average of channels 1,2 and 3) $$X_5/(\sum_{n=1}^{3}X_n)/3$$ But what I want to know is if there are any statistical methods I can use to determine whether a given protein ($X_5$) is significantly different (in quantity) from ($X_{1,2,3}$). Having only three samples with no treatment obviously means establishing any variance is probably meaningless, so I was thinking that some non-parametric rank based test might be the way to go? What is the best test for this setup?


title: Likelihood of (bivariate?) Bernoulli distribution

text: In 50 individuals we measure the presence of 2 SNP's ("single nucleotide polymorphisms") and 1 trait. We assume that both SNP's and the trait are binary. We assume the model $Y_{ijk} = \text{Bernoulli}(p_{ij})$, where $$Y_{ijk} = 1\text{ if person }k\text{ with SNP 1 }= i \text{ and SNP 2 }= j\quad (i,j\in \\{0,1\\})\text{ has trait 1}$$ and $Y_{ijk} = 0$ otherwise. I want to write down the likelihood of the data but I'm not quite sure what it looks like. Is it true that $P(Y_{ijk} = y) = p_{ij}^y(1-p_{ij})^y$, $y =0,1$, or is this wrong? On the internet I find only one dimensional distributions so I just want to check if the density is correct. Thanks in advance and sorry for the silly question!


title: Removal of statistically significant intercept term boosts $R^2$ in linear model

text: In a simple linear model with a single explanatory variable, $\alpha_i = \beta_0 + \beta_1 \delta_i + \epsilon_i$ I find that removing the intercept term improves the fit greatly (value of $R^2$ goes from 0.3 to 0.9). However, the intercept term appears to be statistically significant. With intercept: > >     Call: >     lm(formula = alpha ~ delta, data = cf) >   >     Residuals: >          Min       1Q   Median       3Q      Max >     -0.72138 -0.15619 -0.03744  0.14189  0.70305 >   >     Coefficients: >                 Estimate Std. Error t value Pr(>|t|)   >     (Intercept)  0.48408    0.05397    8.97   <2e-16 *** >     delta        0.46112    0.04595   10.04   <2e-16 *** >     --- >     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 >   >     Residual standard error: 0.2435 on 218 degrees of freedom >     Multiple R-squared: 0.316,    Adjusted R-squared: 0.3129 >     F-statistic: 100.7 on 1 and 218 DF,  p-value: < 2.2e-16 >   Without intercept: > >     Call: >     lm(formula = alpha ~ 0 + delta, data = cf) >   >     Residuals: >          Min       1Q   Median       3Q      Max >     -0.92474 -0.15021  0.05114  0.21078  0.85480 >   >     Coefficients: >           Estimate Std. Error t value Pr(>|t|)   >     delta  0.85374    0.01632   52.33   <2e-16 *** >     --- >     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 >   >     Residual standard error: 0.2842 on 219 degrees of freedom >     Multiple R-squared: 0.9259,   Adjusted R-squared: 0.9256 >     F-statistic:  2738 on 1 and 219 DF,  p-value: < 2.2e-16 >   How would you interpret these results? Should an intercept term be included in the model or not? **Edit** Here's the residual sums of squares:               RSS(with intercept) = 12.92305     RSS(without intercept) = 17.69277


title: Beyond observable and unobservable data, is there any term "semi-observable" defined?

text: When dealing with data in fields such as Natural Language Processing(NLP) or Speech Recognition(ASR) and trying to model the data using Hidden Markov Model(HMM) one should first make it clear that if the available data is observable or unobservable. When applying Markov property to the data, that is the current data probability at time T is only dependent on the data at time T-1 (or more). If all the variables in probability function are observed in the training data, it is called observed data, thus could be trained simply using Maximum Likelihood criterion. In cases, that there are variables in the probability function, which are not observed in the training data it is called a "hidden variable" and so using HMM and Expectaion-Maximization(EM) training criterion, by just defining the number of hidden states, it can be trained. I'm trying to address a special case, in which there are many hidden states (e.g. 100 000 hidden states) which makes it impossible for EM to train the model. By some assumption on the model, by observing each training instance and using the assumptions, the model can restrict the model to a limited number of states among all those hidden states, e.g. every instance is limited to 10-20 states among 100 000 hidden states. This special case has been worked in some problems such as Hidden Vector States proposed by Yulan He. Is this called "semi-observable" case or it's just a hidden variable with some simple assumptions made it possible to be trained? How is it converted to a standard HMM problem? EDIT: In case there is some ambiguities about the terms 'hidden' or 'observable', please read this tutorial on Hidden Markov Model by Rabiner describing all these terms a great manner. EDIT2: This cool tutorial on Expectation-Maximization has a good and brief explanation on 'hidden variables' too.


title: How Conditional Random Fields and Logistic Regression could be the same?

text: If you have read this tutorial about `CRF`, on page 4 under the section `Classification`, it wants to relate `CRF` to `Logistic Regression` (or `Maximum Entropy`, as it is known by this name in NLP community). The formulation of `Logistic Regression` is said to be $ p(y|x) = \frac{1}{Z(X)}exp \\{ \lambda_y + \sum_{j=1}^{K}\lambda_{y,j}x_{j}\\} $. On the other hand, on the Wikipedia page for logistic regression, the formulation is $ f(t) = \frac{1}{1+e^{-t}} $ and the classifier's formulation is the explanation of this formulation which is derived from Logistic Function. My question is, first, how Logistic Regression and Maximum Entropy Classifier are related to each other and are identical? Second, the CRF tutorial says Logistic Regression (or Maximum Entropy Classifier) is simple CRF model, how?


title: why to do structure learning for Bayesian networks?

text: Given a very-large data set, if our goal is to do probabilistic inference, what are the main advantages of learning a Bayesian network from data and then, use the Bayesian network to compute conditional probabilities? I see that we could also approximate these probabilities directly from the data set by counting. Further, if the data set is large enough, one could also try to use the CLT/WLLN to compute some confidence intervals. Why would it be better to construct a Bayesian network (hard optimization problem) and do inference (message passing algorithms)? Some kind of over-fitting argument? Thanks!


title: From a statistical perspective, can one infer causality using propensity scores with an observational study?

text: **Question:** From the standpoint of statistician (or a practitioner), can one infer causality using propensity scores with an observational study ( **not an experiment** )? Please, do not want to start a flame war or a fanatical debate. **Background:** Within our stat PhD program, we've only touched on causal inference through working groups and a few topic sessions. However, there are some very prominent researchers in other departments (e.g. HDFS, Sociology) who are actively using them. I've already witnessed some pretty heated debate on this issue. It is not my intention to start one here. That said, what references have you encountered? What viewpoints do you have? For example, one argument I've heard against propensity scores as a causal inference technique is that one can never infer causality due omitted variable bias -- if you leave out something important, you break the causal chain. Is this an unresolvable problem? **Disclaimer:** This question may not have a correct answer -- completely cool with clicking cw, but I'm personally very interested in the responses and would be happy with a few good references which include real-world examples.


title: How to determine error bars given number of flips

text: I have a biased coin; it's going to return either heads or tails at some percentage. If I run a test and get back say 28 heads out of 40 flips, how can i best add error bars to indicate i don't have a lot of samples, and the real results could be higher or lower?


title: Smoothing time series data

text: I am building an android application that records accelerometer data during sleep, so as to analyze sleep trends and optionally wake the user near a desired time during light sleep. I have already built the component that collects and stores data, as well as the alarm. **I still need to tackle the beast of displaying and saving sleep data in a really meaningful and clear way, one that preferably also lends itself to analysis.** A couple of pictures say two thousand words: (I can only post one link due to low rep) Here's the unfiltered data, the sum of movement, collected at 30 second intervals ![](http://electricsleep.googlecode.com/files/es-ss2.png) And the same data, smoothed by my own manifestation of moving average smoothing ![](http://electricsleep.googlecode.com/files/es-ss1.png) _edit) both charts reflect calibration- there is a minimum 'noise' filter and maximum cutoff filter, as well as a alarm trigger level (the white line)_ Unfortunately, neither of these are optimal solutions- the first is a little hard to understand for the average user, and the second, which is easier to understand, hides a lot of what is really going on. In particular the averaging removes the detail of spikes in movement- and I think those can be meaningful. So why are these charts so important? These time-series are displayed throughout the night as feedback to the user, and will be stored for reviewing/analysis later. The smoothing will ideally lower memory cost (both RAM and storage), and make rendering faster on these resource-starved phones/devices. Clearly there is a better way to smooth the data- I have some vague ideas, such as using linear regression to figure out 'sharp' changes in movement and modifying my moving average smoothing according. I really need some more guidance and input before I dive headfirst into something that could be solved more optimally. Thanks!


title: Creating plots for String type columns in R

text: I have imported a datafile in R. It has different columns. There is one column that has the name of Operating System belonging to that row information. I wanted to get the **percentage** share of each unique OS (windows, linux, ios etc) from this column and plot it. I am very new to R and would like to know if there is any inbuilt way to do this in R.


title: Comparing and visualizing rates of subsystems

text: I am trying to track the performance of a set of homogeneous subsystems vs the system as a whole. The performance metric is measured by a mean value based on a set of tasks performed by a the subsystem, for example, the mean wait time before scheduling task. Each subsystem may be responsible for varying numbers of tasks. I can calculate the mean wait time for each subsystem and for the system as a whole. My question is how to measure and visualize how each subsystem is performing compared to the whole system average (or maybe there's a better aggregate measure for the whole system?) Obviously being below the mean shows the subsystem is performing relatively better and above the mean vice versa, and the distance from the whole system mean represents that, but I struggle to get real meaning or provide a good visualization from that. Is there a better way? Or if I'm on the right track, how do I assign meaning / visualize the metric?


title: degree to become a statistician

text: what is the minimum entry degree to be employed as a statistician


title: Relative influences of terms in the law of total cumulance

text: I posted this to mathoverflow and no one's saying anything. $\newcommand{\cum}{\operatorname{cum}}$ Denote the joint cumulant of several random variables by $\cum(A,B,C,\ldots)$ (more precisely, the cumulant of the joint probability distribution). In particular, what is called the $n$th cumulant of (the probability distribution) of a random variable $A$ is $\cum(\ \underbrace{A,\ldots,A}_{\text{$n$}} \ )$. Then Brillinger's law of total cumulance (Brillinger didn't use that term) says (in the concrete case of exactly four random variables) the following. There is one term for each partition of the set of random variables (so in the case of just four of them there are $15$ terms): $$ \begin{align} \cum(A,B,C,D) & = \cum(\cum(A,B,C,D\mid X)) \\\\[8pt] & \phantom{={}} + \underbrace{\cum(\cum(A,B,C\mid X),\cum(D\mid X)) + \cdots}_{\text{4 terms}} \\\ & \phantom{={}} + \overbrace{\cum(\cum(A,B\mid X),\cum(C,D\mid X)) + \cdots}^{\text{3 terms}} \\\\[8pt] & \phantom{={}} + \underbrace{\cum(\cum(A,B\mid X),\cum(C\mid X),\cum(D\mid X))+\cdots}_{\text{6 terms}} \\\\[8pt] & \phantom{={}} + \cum(\cum(A\mid X),\cum(B\mid X),\cum(C\mid X),\cum(D\mid X)), \end{align} $$ and similarly for numbers of random variables other than $4$, so for example, for six random variables, there are $203$ terms since there are that many partitions of the set of six. Now:   * If $A,B,C,D$ are independent of $X$, then all terms vanish except the first, and the right and left sides of the equality become identical; and   * If $A,B,C,D$ are completely determined by $X$, then all terms vanish except the last, which again simplifies to be identical to the expression to the left of "$=$". So one can **GUESS** that as one moves from one extreme to the other, with varying degrees of dependence, the first term goes from being predominant to contributing little, and the last from contributing very little to predominating. Statisticians have written innumerable volumes about the special case $\cum(A,A)$, where that is the same as $\operatorname{var}(A)$, so you're just talking about the "explained" and "unexplained" components of the variance, the the proportion that is explained is just the thing usually denoted $R^2$. The first of the $15$ terms above would be the component that's wholly unexplained by $X$ and the last would be the component that's wholly explained. Is there some way to make that guess precise? Should one just regard the quotients from dividing each term divided by the total as a sort of summary measure of dependence? In the case $\cum(A,A)$, every number in $[0,1]$ can occur as the quotient $$\frac{\text{explained component of variance}}{\text{total variance}},$$ and no numbers $>1$ or $<0$ can occur. Which $15$-tuples of numbers summing to $1$ can occur as the tuple whose components are the terms on the right side divided by the one term on the left side? And more generally, for other numbers of random variables than $4$? Is that a hard problem comparable to moment problems?


title: Slack values in Data Envelopment Analysis

text: We're using Excel Solver to model a DEA problem. We've worked with a few resources and feel confident that the model is returning correct results but I need some help interpreting them. When solver completes its analysis, it provides an answer report that gives the slack found in some of the constraints -- which are binding and which are not. Here's a sample result for a DMU that was found to be inefficient when Solving to Minimize the Inputs:                             WeightedSum             Outputs = DMU                of Reference Set          Inputs = DMU * E   Slack     Output1        2,892,230     >=        2,892,230        binding     Output2        66,229,077    >=        66,229,077       binding     Output3        9,724,273     >=        5,372,010        4,352,263     Output4        3,123,300     >=        2,511,864        611,436     Output5        21,609        >=        16,039           5,570     Output6        120,527       >=        33,234           87,292     Output7        8,548         >=        2,769            5,779     Output8        17,314        >=        12,605           4,709     Output9        748           >=        616              132      Input1        19,775        <=        23,314           3,539      Input2        9,132         <=        9,132            binding      Input3        20,028        <=        143,604          123,576      Input4        74,293        <=        74,293           binding      Input5        15,558        <=        15,558           binding      Is it correct to assert that those inputs found to be binding could not be further optimized but those that contained slack are where the inefficiencies are found? Can the same be said for Outputs? or does that work the opposite? Or is slack just a bi-product of the analysis and not relevant? Thanks for your help!


title: Convergence in Probability

text: In the following question, I was able to compute part (i), and part (ii), but I am finding hard time to compute with part (iii). ![enter image description here](http://i.stack.imgur.com/19FdI.jpg) So far I have computed the MLE of $\theta$, which is $\frac{\sum_{i=1}^n Y_i -n}{\sum_{i=1}^n Y_i - M} $ I was trying to use the following concept for convergence ![following concept for convergence](http://i.stack.imgur.com/GdCIv.jpg) but I am finding it hard to proceed. Any help would be highly appreciated.


title: Which statistical method to use: Kruskal-Wallis, one-way ANOVA, or another?

text: I have a csv file with time data values (elapsed times) of different lengths for three groups X, Y and Z, knowing that:   * the samples are independent   * normality checking via Shapiro-Wilk fails for Y and Z groups (it does not fail for X but according to X's Q-Q plot X should not be normally distributed)   * the stardard variations are (sd(x) in R): sd(X)=0.003577565 , sd(Y)=0.013948358 , sd(Z)=0.006929825 (Does these tell the presence of unequal variances?) I have generated some plots (boxplot, and Q-Q plots for X, Y and Z) that you can download from here. What statistical method should I use to determine statistical significant differences between the three groups? Must I use Kruskal-Wallis instead of one-way ANOVA? Thanks


title: Short online videos that could assist a lecture on statistics

text: When lecturing statistics, it can be useful to incorporate the occasional short video. My first thoughts included:   1. Animations and visualisations of statistical concepts   2. Stories regarding the application of a particular technique   3. Humorous videos that have relevance to a statistical idea    4. Interviews with a statistician or a researcher who uses statistics Are there any videos that you use when teaching statistics or that you think would be useful? Please provide:   * Link to the video online   * Description of content   * To what statistical topic the video might relate


title: Huge ΔAIC between GAM and GAMM models

text: I'm working with spatial fisheries catch data and environmental variables, and I'm correlating the abundance in the catch to some oceanographic parameters. I'm using a Generalised Additive Model (GAM) and a Generalised Additive Mixed Model (GAMM) with one and two random effects (`mgcv` package in `R`), in particular:               m1 <- gam(kg ~ s(var1, k=10) + s(var2, k=10) + ... + offset(time),                family = Gamma(link=log), method=REML)          m2 <- gamm(kg ~ s(var1, k=10) + s(var2, k=10) + ... + offset(time),                 family = Gamma(link=log), method=REML, random= list(vessel = ~ 1))          m3 <- gamm(kg ~ s(var1, k=10) + s(var2, k=10) + ... + offset(time),                family = Gamma(link=log), method=REML,                random= list(vessel = ~ 1, dayOfTheYear = ~ 1))      I get quite reasonable results in terms of residuals and covariates partial effects, but I get a huge difference in AIC between m1 and m2/m3, i.e. `AIC(m1) = 58998`, `AIC(m2$lme) = 9410`, and `AIC(m3$lme) = 8751`. From what I understood, I should be able to compare GAM and GAMM from `mgcv` using `AIC`... Am I doing something wrong? Thanks!


title: Expression for conditional expectation

text: Suppose we have two random variables: $X$ is a continuous r.v.; $Y$ is a discrete r.v. taking values $0$ and $1$. Is the following expression true? $E[(E[X|Y])^{2}]= [(E[X|Y=1])^{2}]\times P(Y=1)+ [(E[X|Y=0])^{2}]\times P(Y=0)$ ?


title: HLM vs linear mixed models vs random effects model

text: would anyone please tell me the differences among the three methods? also i have a data in which the dependent variable is a count variable. is it classified as a level 1 variable? one of my classmates tried using HLM, and he clustered the DV into groups. Im not really familiar with hierarchical modeling. so i would really appreciate if someone could help me with this one?


title: How can I (numerically) approximate values for a beta distribution with large alpha & beta

text: Is there a numerically stable way to calculate values of a beta distribution for large integer alpha, beta (e.g. alpha,beta > 1000000)? Actually, I only need a 99% confidence interval around the mode, if that somehow makes the problem easier. **Add** : I'm sorry, my question wasn't as clearly stated as I thought it was. What I want to do is this: I have a machine that inspects products on a conveyor belt. Some fraction of these products is rejected by the machine. Now if the machine operator changes some inspection setting, I want to show him/her the estimated reject rate and some hint about how reliable the current estimate is. So I thought I treat the actual reject rate as a random variable X, and calculate the probability distribution for that random variable based on the number of rejected objects N and accepted objects M. If I assume a uniform prior distribution for X, this is a beta distribution depending on N and M. I can either display this distribution to the user directly or find an interval [l,r] so that the actual reject rate is in this interval with p >= 0.99 (using shabbychef's terminology) and display this interval. For small M, N (i.e. immediately after the parameter change), I can calculate the distribution directly and approximate the interval [l,r]. But for large M,N, this naive approach leads to underflow errors, because x^N*(1-x)^M is to small to be represented as a double precision float. I guess my best bet is to use my naive beta-distribution for small M,N and switch to a normal distribution with the same mean and variance as soon as M,N exceed some threshold. Does that make sense?


title: Follow up: In a mixed within-between ANOVA plot estimated SEs or actual SEs?

text: I am currently finishing a paper and stumbled upon this question from yesterday which led me to pose the same question to myself. Is it better to provide my graph with the actual standard error from the data or the one estimated from my ANOVA?   As the question from yesterday was rather unspecific and mine is pretty specific I thought it would be appropriate to pose this follow-up question. **Details:**   I have run an experiment in some cognitive psychology domain (conditional reasoning) comparing two groups (inductive and deductive instructions, i.e., a between-subjects manipulation) with two within-subjects manipulations (problem type and content of the problem, each with two factor levels). The results look like this (left panel with SE-estimates from the ANOVA Output, right panel with SEs estimated from the data): ![alt text](http://i.stack.imgur.com/fHmrK.png)   Note that the different lines represent the two different groups (i.e., the between-subjects manipulation) and the within-subjects manipulations are plotted on the x-axis (i.e., the 2x2 factor levels). In the text I provide the respective results of the ANOVA and even planned comparisons for the critical cross-over interaction in the middle. **The SEs are there to give the reader some hint about the variability of the data.** I prefer SEs over standard deviations and confidence intervals as it is not common to plot SDs and there are severe problems when comparing within- and between-subjects CIs (as the same surely applys for SEs, it is not so common to falsely infer significant differences from them). To repeat my question: **Is it better to plot the SEs estimated from the ANOVA or should I plot the SEs estimated from the raw data?** **Update:**   I think I should be a little bit clearer in what the estimated SEs are. The ANOVA Output in SPSS gives me `estimated marginal means` with corresponding SEs and CIs. This is what is plotted in the left graph. As far as I understand this, they should be the SDs of the residuals. But, when saving the residuals their SDs are not somehow near the estimated SEs. So a secondary (potentially SPSS specific) question would be:   **What are these SEs?** * * * **UPDATE 2:** I finally managed to write a R-function that should be able to make a plot as I finally liked it (see my accepted answer) on its own. If anybody has time, I would really appreciate if you could have a look at it. Here it is.


title: How many observations are enough to perform linear regression with fixed effects

text: I am new to econometrics. I am studying earnings management in banks during the financial subprime crisis. I manage to collect data from 23 banks from 2005 to 2010. Few years of them are missing but I end up with 126 observations. Is that enough to perform a regression analysis with fixed effects (I do not know if that last actually has an impact in the number of minimum required observations).


title: Stepwise introduction of predictors to mixed-effects models

text: As the title says, what I'd like to do is stepwise introduction of predictor variables to a mixed-effects model. I'm going to first say what I'd be doing if it were stepwise linear regression, just to make sure I've got that part right, and then describe the full model to which I want to apply an analogous approach. I have a student population who took a pretest, then a tutorial, then a posttest. The tutorial involved doing problems from several categories with feedback, and the users could control which category the next problem would come from and when to stop the tutorial. I want to create a model that will account for posttest performance using pretest score and some measures of behavior during the tutorial, including total number of problems done, accuracy, and probability of switching category. The last of these is of greatest theoretical interest. There are other variables I'm not mentioning for simplicity. For the linear regression approach, I first did a simple regression using posttest score as the DV and including the main effects (only) of pretest score, tutorial accuracy, and number of problems as predictors. Then, I added probability of switching as an additional predictor, and compared the resulting model to the previous one to see if it had significantly better explanatory power (it did). The R code I used is below.               lm1 <- lm( posttestScore ~ pretestScore + practiceAccuracy + practiceNumTrials, data=subj.data )     lm2 <- lm( posttestScore ~ pretestScore + practiceAccuracy + practiceNumTrials + probCategorySame, data=subj.data )     anova( lm1, lm2 )      So far so good? OK, next, I switched to a mixed model in order to include a binary within-subjects factor, 'test question type'. Both pretest and posttest have values for each level of this factor for every subject. (It's unrelated to the 'problem category' I mentioned for the tutorial.) The other predictors, however, only have one value for each participant. My models then became:               library( nlme )     lm1 <- lme( posttestScore ~ pretestScore + questionType + practiceAccuracy + practiceNumTrials, random=~1|sid, method="REML", data=D )     lm2 <- lme( posttestScore ~ pretestScore + questionType + practiceAccuracy + practiceNumTrials + probCategorySame, random=~1|sid, method="REML", data=D )      However, I don't know how to test whether the second model resulted in a significant improvement over the first model. Is that the right question I should be asking and, if so, how should I do it?


title: How to test a yes/no outcome with different inputs?

text: I have what I think is a very simple beginner question, but I do not have any formal training or knowledge of statistics or design of experiments. Let's say I have a yes,no (0,1) outcome of an experiment, say whether or not a customer purchased a product given a coupon. There are 5 kinds of coupons, each with different colors, or text, i.e. we want to find the most effective coupon. Say we collect this data:               coupon sample-size yes-count           1        2589       153          2        1789       212          3        3097       298          4        2654       281          5        2103       211      What is the right / best way to go about telling which is the best coupon? As a follow up what if I have two dimensional inputs, say coupon and store-type where store-type can be one of 3 different values.               coupon  store-type  sample-size yes-count          1           1          387        48          1           2          675       117          1           3          467        88      ...      Then what would be the best / right way to say which combination of coupon and store-type is most effective?


title: $\chi^2$ tests to compare the fit of large samples logistic models

text: Does anyone know of any $\chi^2$ tests to compare the fit of logistic models which factor out the sample size? I'm dealing with a very large sample and I fear the significant $\chi^2$ test I get when adding a single variable to the model is simply the result of the sample size (>200,000 cases). I'm doing what is known as differential item functioning analysis with logistic regression. Basically it's as if I'm checking whether giving the right answer to a question (dependent variable) depends on your ethnicity when controlling for the total exam score. Model 1 `Q1~TotalexamScore` Model 2 `Q1~TotalexamScore+ Group` I'm basically using a chi-squared test to compare model1 to model2. The coefficient significance is not that important but $\chi^2$ and sometimes $R^2$ are generally recommended to check differential item functioning. My problem is that my sample is very large. In theory (for the question I'm considering) there should be no real difference across groups, so I suspect it's simply the sensitivity of the $\chi^2$ to sample size. I'd rather use the whole dataset instead of taking (small) random samples as it is highly skewed. I've seen things like Phi and Cramer's V for crosstabs but I'm not sure whether they have been used before on logistic regression, if there are better ones and if there are any packages (I generally use Spss, Mplus, Stata, R).


title: Generate data from a bivariate power-law distribution in R

text: I need to generate data from a random vector that follows a bivariate power- law: $$ f_{X,Y}(x,y) = \frac{C}{XY} \left(\frac{X}{X_0} \right)^{-\alpha} \left(\frac{Y}{Y_0} \right)^{-\beta} , $$ where $X_0,Y_0,C>0$ are constants. I can generate data from a univariate power-law by using a uniform random variable and transform it (I am using $\texttt{R}$), but I have no idea how to do it for a random vector. Do you have any suggestions ?


title: Interpreting error when trying to import dta file in R

text: I'm trying to solve a couple of exercises in Gelman and Hill 2007. The datasets they use are all on a website, however, every time I try to import the files I get an error: > pollution=read.table("pollution.dta") Error in scan(file, what, nmax, sep, > dec, quote, skip, nlines, na.strings, : line 1 did not have 2 elements Here's the dataset: http://www.stat.columbia.edu/~gelman/arm/examples/pollution/ I guess there's a problem with the datasets, but I'm wondering if there's a modification that could enable me to import it?