title: P-value for ML estimate of error-autoregressive parameters in R using GLS with ARMA(p=1) correlation

text: I would like to know if there is a way to get p-values when using the GLS function (part of the nlme package) for the ML estimate of the error- autoregressive parameters. I looked through the package information on nlme and have not found a clear example of how to do this. Here is an example of some code that illustrates my question:               x <- rnorm(100)     y <- rnorm(100)     z <- rnorm(100)     df <- data.frame(x,y,z)          test.reg <- gls(x ~ y+z,data=df,correlation=corARMA(p=1),method="ML")          summary(test.reg)          Generalized least squares fit by maximum likelihood       Model: x ~ y + z        Data: df         AIC      BIC    logLik       289.3276 302.3534 -139.6638          Correlation Structure: AR(1)      Formula: ~1       Parameter estimate(s):        Phi      0.01887445           Coefficients:                      Value  Std.Error   t-value p-value     (Intercept) 0.08418063 0.10214246 0.8241492  0.4119     y           0.00236772 0.09410912 0.0251593  0.9800     z           0.07689263 0.09480738 0.8110405  0.4193           Correlation:        (Intr) y          y  0.135            z -0.013  0.074          Standardized residuals:             Min          Q1         Med          Q3         Max      -1.94366502 -0.66310726  0.09546043  0.64881478  2.78384648           Residual standard error: 0.9781187      Degrees of freedom: 100 total; 97 residual      What I would like to know more about is the Parameter estimate Phi. I have searched through the nlme package to get more information about p-values for Phi but have come up short. I am unsure if the gls function computes the p-values for Phi. I do know this gives the coefficients of phi:               coeff <- test.reg$modelStruct$corStruct          Correlation structure of class corAR1 representing        Phi      0.01887445       But that is all I have really found. Any help would be great! Thanks so much, John


title: How odd is a cluster of plane accidents?

text: Original question (7/25/14): Does this quotation from the news media make sense, or is there a better statistical way of viewing the spate of recent plane accidents? > However, Barnett also draws attention to the theory of Poisson distribution, > which implies that short intervals between crashes are actually more > probable than long ones. > > "Suppose that there is an average of one fatal accident per year, meaning > that the chance of a crash on any given day is one in 365," says Barnett. > "If there is a crash on 1 August, the chance that the next crash occurs one > day later on 2 August is 1/365. But the chance the next crash is on 3 August > is (364/365) x (1/365), because the next crash occurs on 3 August only if > there is no crash on 2 August." > > "It seems counterintuitive, but the conclusion follows relentlessly from the > laws of probability," Barnett says. Source: http://www.bbc.com/news/magazine-28481060 Clarification (7/27/14): What is counter intuitive (to me) is saying that rare events tend to occur close in time. Intuitively, I would think that rare events would not occur close in time. Can anyone point me to a theoretical or empirical expected distribution of the time between events under the assumptions of a Poisson distribution? (That is, a histogram where the y-axis is frequency or probability and the x-axis is time between 2 consecutive occurrences grouped into days, weeks, months, or years, or the like.) Thanks. Clarification (7/28/14): The headline implies it is more likely to have clusters of accidents than widely spaced accidents. Lets operationalize that. Let's say that a cluster is 3 airplane accidents, and a short period of time is 3 months and a long period of time is 3 years. It seems illogical to think that there is a higher probability that 3 accidents will occur within a period of 3 months than within a period of 3 years. Even if we take the first accident as a given, it is illogical to think that 2 more accidents will occur within the next 3 months as compared to within the next 3 years. If that is true, then the news media headline is misleading and incorrect. Am I missing something?


title: how to compare a user-based collaborative filtering with item-based

text: Apart from prediction accuracy, what are good evaluation metrics for comparing user-based and item-based collaborative filtering


title: MLE or instrumental variables

text: I'm trying to estimate a model in which one of the explanatory variables is correlated with the error term. As I see it there are two alternatives, specify the likelihood function and maximize it to get the estimates of the coefficients, or use instrumental variables. There are pros and cons to each approach, so I'd like to get some other opinions on the matter. the IV method is quicker, since it works with closed form equations, but the drawback is that I'll get different estimates for different instruments, which I don't really like. Any advice would be welcome.


title: Pegasos prediction time

text: I have been implementing the (kernelized) Pegasos algorithm, but am running into problems in terms of scalability. I will use notations as in the original manuscript. Here's a typical time measurement:               percent         m   t1      t2     0.05%           32  701     961     0.1%            64  676     1890     0.2%            128 648     2777     0.41%           256 725     3851     0.83%           512 861     5267     1.67%           1024 1293    7635     3.34%           2048 1928    11404     6.68%           4096 3256    18210     13.36%          8192 6312    31750     26.72%          16384 11420   56923     53.44%          32768 18838   96755      Where $m$ is the number of data points in the training set, t1 the training time (includes cross validation over many folds) and t2 the testing time (just 1 prediction over a larger, but fixed part of the data). As you can see, the train time scales reasonably well, but the prediction time scales very badly. I think this is due to the fact that the Pegasos algorithm requires one to compute the (kernel) product of every test-point with a large number of training inputs, that increases as the total training set ($m$) increases. Are there known ways to cope with this problem? I thought the L1-norm would take care of this problem due to sparsification of the model (#support vectors), but I don't think it does. I was thinking of simply dropping terms with low 'counts' ($\alpha$ value in Fig. 3 in the original manuscript), but that seems terribly unelegant ...


title: How can I sample from a distribution with incomputable CDF?

text: Semi-computer science simulation related problem here. I have a distribution where P(x) = $\frac{(e^b-1) e^{b (n-x)}}{e^{b n+b}-1}$ for some constants b and n, and x is an integer such that $0\leq x \leq n$. Now, I need to sample from this distribution. It has an invertible CDF, so it's possible to do this directly in theory. The problem is that the numbers involved are LARGE. So large in fact, that they both overflow conventionally formatted variables, and take at least minutes (at some point I gave up...) to compute using arbitrary precision formats. Basically, the inverse CDF still involves are term of $e^{b(n+1)}$, for $ 350 < n < 3500$. Despite this, the output numbers will still be in the range $0-n$, so it seems like there should be a way of doing this. What I'm looking for is a way of approximately sampling from this distribution that _is_ computable. Are there alternative methods of sampling? What are they?


title: Estimation After Selection on Non-central F Random Variables

text: Suppose that you observe $F_1,F_2,\ldots,F_k$ each _independently._ drawn from non-central F distributions with common, known, d.f. $\nu_1, \nu_2$, and with (unknown) non-centrality parameters $\lambda_1,\lambda_2,\ldots,\lambda_k$. Suppose that the sample is ordered by $F_{(1)} \le F_{(2)} \le \ldots \le F_{(k)}$. I would like to somehow estimate $\lambda_{(k)}$. That is, estimate the non-centrality parameter associated with the largest of the $F_i$. (Actually, procedures for estimating arbitrary $\lambda_{(i)}$ would be nice too.) There is a large body of literature on estimation after selection that typically assumes the RVs are normally distributed (Gupta and Miescke _inter alios_ ), or exponential, uniform, etc. I looked into the normalizing transforms of the non-central F, but these require you to know the $\lambda_i$ (they seem to exist to construct tables), and don't work well with estimates $\hat{\lambda_i}.$ Milton & Rizvi (1989) is a related reference.


title: Question on the appropriate way to calculate the CDF of a volatile measure

text: If I track the the volatility of a stock on a sample period and at every point in time where I record the volatility I would like to also record the CDF of the volatility. Should I take the entire volatility sample and compute its CDF or should I compute and update the CDF at every point in time? In MATLAB I did the following `[dens, xi]=ksdensity(volatility,volatility,'function','CDF')` and I get the CDF at every volatility. Is this the correct way to do it? What are the alternatives and the pros and cons?


title: What is better to transform doing linear regression, response or explanatory variable(s)?

text: Maybe it is a very basic question and already answered, but I could not find a clear answer. My plots of response vs. predictors show "curved" relationship, and log- transformation can help to achieve linearity. But it helps either when I log- transform response variable or when I log-transform explanatory variables, as well. The distributions of them all are more or less close to normal. So, what is better to transform, response or predictors, from statistical point of view? Or it doesn't matter for regression and I only need to seek to normal distributions? (I see here that when predictors are not transformed then R2 is related to the variance of the residuals and can be trusted. Is it the only reason?)


title: Does Poisson Regression have an error term?

text: I was just wondering if Poisson regression has an error term? Can a Poisson regression have random effects and an error term? I am confused about this point. In logistic regression, there is no error term because your outcome variable is binary. Is that the only glm model that doesn't have a residual term?


title: Finding the correct data mining approach

text: (I apologise for being a newb, but I'm a researcher introducing myself to data mining---any help or insight would be greatly appreciated. Also, this isn't technically a homework question, but I've marked it as such to distinguish that this is a broad question from a non-specialist). I'm having trouble finding the correct approach to my (fairly simple) example. Let's say I have months of data for log-in times of a certain website. The data has been selected and cleaned such that I have a list of Date_Time for each log-in. Now, suppose I wanted to predict the log-ins for the next two weeks by day and hour, based on these past trends. I imagine I would cluster the data by day (assuming beforehand that there will be different trends with respect to Monday vs. Friday) and make some regression analysis to predict the next two (say) Mondays. Similarly, I could cluster by the hour and do a regression analysis to extrapolate the trend of log-ins. Anyone know of a resource which tells you how to do this in Python? I want to keep this example fairly straightforward, but I'm open to any more ideas on how to model this behavior more efficiently.


title: Granger causality as a means of validation

text: I am busy validating a model for "run-off" using river flow information. I have approximately 5 years of daily data and have been looking at G-causes as an additional statistic to help explain the difference between observed and simulated run-off within a catchment. Using R and running the granger.test() from Library(MSBVAR) I have the following results > granger.test(Daily.sea, p=10) F-statistic p-value Observation -> Simulation > 3.363867 0.0002219447 Simulation -> Observation 44.726567 0.0000000000 Is it fair to say the both Observation g-causes Simulation and vice versa? This indicates that there is indeed a relationship between the two variables. I ask because my results using hydroGOF are not that good (as in weak relationships between the two).


title: Calculate chance of big drop in a timeline

text: I am not much of a mathematics person, but I have been asked to create a program that will calculate a chance for a big drop on a timeline. I'm pretty sure even the way I try to explain it is not really good, so I will give an example of the data I have, and what I want to know. Lets say I have the following data:               01/01/2012 - 10000         02/01/2012 - 10425         03/01/2012 - 9000         04/01/2012 - 11000         05/01/2012 - 200         06/01/2012 - 600         07/01/2012 - 654         08/01/2012 - 950      and so on. Basically what I want to know is, at what time range in the future the value might drop (big drop). For example, 03/01/2012 is a small drop (which would be nice if I could calculate it also, but thats not the main thing I need) and 05/01/2012 is a big drop (which is what I want to know). Now, I know this it is not possible to get the exact time and date the big drop will happen, but how can I calculate a small time range (few days?) that it might happen based on previous data?


title: How to read Cook's distance plots?

text: Does anyone know how to work out whether points 7, 16 and 29 are influential points or not? I read somewhere that because Cook's distance is lower than 1, they are not. Am, I right? ![enter image description here](http://i.stack.imgur.com/4NMSQ.png)


title: Good way to rank/sort data that falls in the lower range

text: I have a bunch of data like this:               P1 [1, 2, 3, 4, 5, 6, 7, 8, 10]     P2 [5, 8, 10, 12, 20]     P3 [10, 201, 440]     P4 [1, 2, 10]     P5 [1, 2]      Right now it's ranked by the size of its data set (as shown above), so values that are heavily represented on the chart are on top e.g., P1 has the most) I'm trying to rank this data so that it favors those that have lower values in data set. So for instance P5 (with low values) would be pushed higher than it's now. Also, P2 could be ranked higher because although it has less values than P1, the values are mostly in the higher range. Basically I want the data to be sorted by the size of it and its overall tendency to have lower values. I can average it all out but I don't wanna give weight to data with small set (e.g., 1 or 2 values) So something like this would be nice (guessing):               P1 [1, 2, 3, 4, 5, 6, 7, 8, 10]     P5 [1, 2]     P4 [1, 2, 10]     P2 [1, 5, 8, 10, 12, 20]     P3 [10, 201, 440, 500]      I'm looking for an R solution.


title: What book would you recommend for non-statistician scientists?

text: What book would you recommend for scientists who are not statisticians? Clear delivery is most appreciated. As well as the explanation of the appropriate techniques and methods for typical tasks: time series analysis, presentation and aggregation of large data sets.


title: Concise book/tutorial for Statistics (Data Analysis)

text: > **Possible Duplicate:**   >  What book would you recommend for non-statistician scientists? Could anyone suggest a concise book/tutorial on statistics with examples. I am a biologist and need basic concepts with examples. I have planned 10 days for this learning exercise so please keep that in mind while suggesting. Thank you


title: Use install.packages() to build and install all available architectures by default?

text: Is there a way to use `install.packages()` so that it will build and install both 32-bit and 64-bit versions of a package from **source** by default? I have multiple computers, all running Mac OS X on 32-bit and 64-bit Intel CPUs. I keep installed packages in my home directory, and I like the convenience of syncing my home directory between my computers. Right now I have to launch R twice—once for each architecture. It’s not a big deal, just annoying to do each and every time.


title: Predicting a semi-deterministic process

text: Say I have a process that gives me 3 outputs: $O^1$, $O^2$ and $O^3$. The outputs are generated from a semi-deterministic process, i.e. there is a deterministic component in the outputs, along with a random component. In particular, having $n$ measurements over time, the outputs $O_j \quad j=1,2,...,n$ are -at least in part- dependent on the previous outputs. So $O_j = f(O_{j-1}, O_{j-2}, O_{j-3}) + \epsilon$ (I'm not interested in going farther away than 2 or 3 measurements, $\epsilon$ is the random component). So now I have a set of ~150 consecutive measurements, how can I predict what are the likely outputs in the future? I can easily calculate the distribution of values following a certain output, for instance I could say that if $O^1_j$ is between 50 and 60 I have a certain probability of $O^1_{j+1}$ of being between 30 and 40, by looking at the measurements that I took in the past. I did construct some pdf for the distribution of these probabilities, but now I'm a bit stuck, especially because there is probably an interaction between the three outputs (so for instance, updating my previous statement $O^1_j = f(O^1_{j-l}, O^2_{j-l}, O^3_{j-l}) + \epsilon \quad\quad l=1,2,3$) I've been reading about Bayesian predictors and I tought they could be applied here, but I don't know enough about the topic to determine if this is a good choice or if there is something easier/more appropriate. I will appreciate any suggestion!


title: Correct way to calibrate means

text: I'm working on benchmarking the speed of various JavaScript methods. Part of the benchmark process requires repeating a test for a minimum time _(to reduce the percentage uncertainty to at or less than 1%)_. There is some overhead on each test _(the cost of the loop, incrementing a counter variable and so on)_. I currently benchmark an empty test to get the overhead cost. The result has it's own mean and margin of error. How should I apply this calibration? Do I take the calibration mean and subtract it from every measured value that compose other benchmark samples, or should I subtract it from the other benchmark's mean value even with differences in margin of error, or is their some formula to follow?


title: Sampling uniformly from n-sphere using spherical coordinates

text: It has been explained here why sampling from n-sphere is not achievable with naive parametrization. And it explains how to correct it for 3 dimensions. Can somebody please guide me what is the correct method for uniformly sampling in n dimensions using the spherical coordinates? i.e. \begin{align} &x_1 = r \cos(\phi_1) \,\\\ &x_2 = r \sin(\phi_1) \cos(\phi_2) \,\\\ &x_3 = r \sin(\phi_1) \sin(\phi_2) \cos(\phi_3) \,\\\ &{}\,\,\,\vdots\\\ &x_{n-1} = r \sin(\phi_1) \cdots \sin(\phi_{n-2}) \cos(\phi_{n-1}) \,\\\ &x_n = r \sin(\phi_1) \cdots \sin(\phi_{n-2}) \sin(\phi_{n-1}) \,. \end{align}


title: Contemporaneous effect in VAR / VEC models

text: Other than Granger causality, I would also like to test contemporaneous effects. So with a VAR (vector autoregression) model I use a simple t-test to test if the cross correlation of the level residuals is significant as it is a SVAR (structural vector autoregresion). But how can I test the contemporaneous effects in a VEC (vector error correction) model?


title: Must the SD of a 'bell shaped' distribution be less than the SD of a skewed distribution with the same range?

text: > 10) Which of the following distributions will have the smallest standard > deviation, assuming that none contain outliers?   >  a) A uniform distribution of integers with a mean of 5 and a range of 10.   >  b) A bell-shaped distribution of integers with a mean of 5 and a range of > 10.   >  c) A right-skewed distribution of integers with a mean of 4 and a range of > 10.   >  d) All three distributions have equal standard deviation.   >  e) The answer cannot be determined from the information given. The answer is given as (b) but I really can't figure out a simple reason why.


title: How many days are required for neural network prediction of dax future?

text: I am predicting the daily change of the dax future. Input data I am using:   * Dax, vdax, euro bund spy500 close prices as returns   * 8 series from investor series having values between 5 and -5 This is trained using Matlab time series Neural Network with:   * Feed recurrent forwards=10   * Delays=10   * Hiddenlayer(neurons?)=10 **Questions:**   1. How many points are recommended for full set (train, val, test)? $2000$?, $200$?, $500$? If Neural Network using more records is it more likely for overfitting?   2. Is it bad to add more series? Like ratios from above series? Will it be more noise and overfitting?


title: Can I make any tests statistically significant?

text: I am using a simple $\chi^2$ test to see if there is any impact of an independent variable (two charts) on the resulting choices (two choices). I found that as the number of sample increase, the $p$-value is decreasing (type I error is decreasing). What about type II error? Is this decreasing as well because we collect more samples? If this is true, does this mean that if I can collect more and more data, can I make any $\chi^2$ test statistically significant? If this is true, I feel like that I am cheating. (Actually, we collect data through Amazon Mechanical Turk, so collecting more data is not that difficult).


title: Testing coefficients separately and simultaneously

text: From the summary output by `lm`, I have the following from R:               Call:       lm(formula = DJIA[2:262] ~ DJIA[1:261])     Residuals:          Min       1Q   Median       3Q      Max     -176.878  -22.397   -0.641   26.478  125.139     Coefficients:                  Estimate Std. Error t value Pr(>|t|)         (Intercept) 36.898384  42.989100   0.858    0.392         DJIA[1:261]  0.994459   0.007477 133.002   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1     Residual standard error: 42.37 on 259 degrees of freedom     Multiple R-squared: 0.9856,     Adjusted R-squared: 0.9855     F-statistic: 1.769e+04 on 1 and 259 DF,  p-value: < 2.2e-16**      How do I test the coefficients separately and simultaneously in this case? What if I assume that the model has an intercept of 0 and slope of 1? Does it translate to B0 = 0 and B1 = 1?


title: How do you determine the effect of a simple predictor variable after a PLS analysis?

text: So, I am running PLS on a genetic dataset with phenotypic and genotypic information. I have about 1000 binary predictors (X), representing molecular markers, for each individual. My indicator variables (Y) are yield in pounds per acre for each individual. I am predicting yield (Y) using molecular markers (X) with about 3 latent variables. I have made the predictions and am satisfied with the model's ability to predict based on genotype. What I want to know is: How do I determine the effect that each marker is having on the Y prediction within each latent variable? Preferably in units of Y (lbs per acre). I should add that I know that there are more appropriate models to model the effects of individual markers, but I am primarily interested in comparing this model with other models, as well as knowing a few predicted marker effects for explanatory purposes.


title: Looking for a OLS-Equation if one Regressor is correlated with the error

text: How can I express a OLS-Estimator if I know about the correlation i.e. I know that $E(x_i u_i)=\rho$ (I'm not looking for IV or 2SLS). I'll explain my problem with an example: In a simple problem $\ \ \ y_i = \beta_0 + \beta_1 x_i + u_i \ \ \ $ i assume $E(u_i)=0$ $E(x_i u_i) = \rho$ $ \ \ \ i=1,...,n$ With the MM, I can derive the following equations: $\frac{1}{n}\sum_{i=1}^n(y_i - \hat{\beta}_0 - \hat{\beta}_1 x_i) = 0 \ \ \ \ \ \ \leftrightarrow \ \ \hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}$ $\frac{1}{n}\sum_{i=1}^n x_i(y_i - \hat{\beta}_0 - \hat{\beta}_1 x_i) = \rho \ \ \ \leftrightarrow \ \ \hat{\beta}_1 = \frac{\sum_{i=1}^n x_i(y_i - \bar{y})}{\sum_{i=1}^n x_i(x_i - \bar{x})} - \frac{\rho}{\sum_{i=1}^n x_i(x_i - \bar{x})}$ In terms of the ols-estimator and some factor, I can now write the following: $\hat{\beta}_1 = {\hat{\beta}_1}_{ols} - \frac{\rho}{m_{xx}}$ $\hat{\beta}_0 = {\hat{\beta}_0}_{ols} + \bar{y}\frac{\rho}{m_{xx}}$ We can now see by what particular factor the $\boldsymbol{\beta}_{ols}$ estimator is biased. For this very simple case, I could derive the estimators very easily but for the more common case with k regressors, I can't find a good approach to express the new estimator. For instance, I thought I could use $\boldsymbol{X^T(y-X\beta) = \rho} \ \ $ instead of $ \ \ \boldsymbol{X^T(y-X\beta) = 0}$ with $ \ \boldsymbol{\rho} \ $ a $\ \ k\times 1$-vector containing the information about the correlation between the regressors and the error-term. But the solution $\boldsymbol{\hat{\beta} = (X^TX)^{-1}X^Ty - (X^TX)^{-1}\rho}$ does not quite do the job.


title: ARCH, GARCH Forecasting in R

text: I tried to fit `auto.arima()` with a `ts` data. But it is not giving the right forecast. For many it is coming as `arima(0,1,0)` model which is not good at all. Can I fit a GARCH model to the original series in this case? How do you get fitted and forecasted values of original data using `garch(1,1,)` or some other model? I tried to use code for GARCH but it is not giving the fitted and forecast of original values.


title: How to fix the constant variance assumption?

text: We have a project where we have to find the best model using a large set of data. In our current model there are 10 variables, some quantitative, and a few that are qualitative. When we first do regression in Minitab we see that the constant variance assumption is violated because the residual graph is fanning out. Furthermore, the probability plot has an S shape which is not the straight line that we want. We therefore decided to do a natural log transformation on our dependent variable and this straightened our probability plot significantly, but now our residual graph looks like a big clump of points concentrated in one area. Any suggestions on how to fix it? We have tried transforming some independent variables as well, and it didn't fix it. Is there something that we should be looking for? What could be causing this?


title: Reading the output from a GLMM run in R

text: I am a complete novice and dummy when it comes to statistics so I apologise in advance... I have been asked to report the results of my GLMMs (I ran two) in a table. This table must state: effect, standard error, test statistic, and P value, for all fixed effects. Unfortunately I am struggling to read my output. The out put is as follows, if anyone would be kind enough to help I would be very grateful and will know for future reference which bit equates to what (also I have been told my degrees of freedom are different for both the tests, could someone explain why this is?).               GLMM 1-run for predictors of step length.      Response variable = step length.      fixed effects = depth and direction threshold.      random factor = individual          Models:     m2: step ~ (1 | ind)     m1: step ~ Depth * threshold + (1 | ind)     Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)      m2 3 373235 373259 -186615 373229      m1 8 373225 373290 -186605 373209 19.767 5 0.001382 **     ---     Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      .               GLMM 2 -run to investigate potential predictors of PDBA.     response variables = depth and step length.      fixed effect = direction threshold.     random factor = Individual          Models:     m3: PDBA ~ Depth + (1 | ind) + thresholdepth     m2: PDBA ~ step * threshold + Depth * threshold + (1 | ind)     Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)      m3 6 -48205 -48157 24109 -48217      m2 11 -48430 -48341 24226 -48452 235.1 5 < 2.2e-16 ***     ---     Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Models:     m4: PDBA ~ step + (1 | ind) + step:threshold     m2: PDBA ~ step * threshold + Depth * threshold + (1 | ind)     Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)      m4 6 -48206 -48158 24109 -48218      m2 11 -48430 -48341 24226 -48452 233.81 5 < 2.2e-16 ***     ---     Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      Hi, I think the package I used was was lme4? I have run a summary for the first GLMM and this is what I got, I have no idea which parts are relevant though, I assume it doesn't all go in a table?!               > summary(m1)     Linear mixed model fit by REML ['lmerMod']     Formula: step ~ Depth * threshold + (1 | ind)           REML criterion at convergence: 373184           Random effects:      Groups   Name        Variance Std.Dev.      ind      (Intercept) 196519   443.3         Residual             469370   685.1        Number of obs: 23473, groups: ind, 11          Fixed effects:                       Estimate Std. Error t value     (Intercept)      160.95895  134.80279   1.194     Depth              0.06438    0.44777   0.144     threshold2        51.18065   17.62222   2.904     threshold3         1.47733   21.43879   0.069     Depth:threshold2  -1.23654    0.60029  -2.060     Depth:threshold3  -0.09587    0.65088  -0.147          Correlation of Fixed Effects:                 (Intr) Depth  thrsh2 thrsh3 Dpth:2     Depth       -0.094                                 threshold2  -0.090  0.712                          threshold3  -0.075  0.588  0.567                   Dpth:thrsh2  0.071 -0.737 -0.745 -0.435            Dpth:thrsh3  0.064 -0.674 -0.490 -0.857  0.502      For GLMM 1 I ran this code -               m1<-lmer(step~Depth*threshold+(1|ind))     m2<-lmer(step~(1|ind))     anova(m1,m2)      For GLMM 2 I ran this code -               m2<-lmer(PDBA~step*threshold+Depth*threshold+(1|ind))     m3<-update(m2,~.-step*threshold)     anova(m2,m3)      and this one:               m2<-lmer(PDBA~step*threshold+Depth*threshold+(1|ind))     m4<-update(m2,~.-Depth*threshold)     anova(m2,m4)      The output from the Anova only gives me one p value for each GLMM and I think I need a p value for each of the fixed effects within the models? Does anyone know what code I can run to get this? Thank you


title: How to estimate the effect of changes happened in employee designation/title on salary?

text: Is it possible to estimate the beta changes corresponding to employee's departmental change or title change or Both. Can i follow linear regression directly to calculate slope ? I am using R statistical tool to compute beta's. ![enter image description here](http://i.stack.imgur.com/xMIR3.png)


title: correlation with logarithmic transformation

text: I have a dataset of 2.000.000 projects. Each project is defined by its size and the number of active developers. After applying a logarithmic transformation on the project size, I've plotted the data. It seems to me that the number of active developers decreases when the project size goes up. When I try to apply Spearman's correlation rank, I find only little evidence of it (0.09520886) I tried to use linear regression: `lm(size ~ active developers)` (in R), but I got negative results:               Coefficients:             Estimate Std. Error t value Pr(>|t|)         (Intercept)   41.672      2.871  14.513   <2e-16 ***     colls         18.381     11.034   1.666   0.0957     Residual standard error: 4089 on 2126091 degrees of freedom     Multiple R-squared:  1.305e-06, Adjusted R-squared:  8.351e-07      F-statistic: 2.775 on 1 and 2126091 DF,  p-value: 0.09572      Then I tried to use a negative binomial regression `glm.nb(commits ~ colls)`, but I got negative results again:               Coefficients:             Estimate Std. Error z value Pr(>|z|)            (Intercept)  3.73200    0.04862   76.76   <2e-16 ***        colls        0.32318    0.18683    1.73   0.0837 .           Null deviance: 2168.2  on 2126092  degrees of freedom     Residual deviance: 2164.6  on 2126091  degrees of freedom     AIC: 41924708          Number of Fisher Scoring iterations: 1                         Theta:  0.0002086134            Std. Err.:  7.43e-08           2 x log-likelihood:  -41924701.8870000020      qchisq(0.95, 2126091) =  2129484      Which method I should use to get more insights about the relation between these two variables? ![enter image description here](http://i.stack.imgur.com/dXGuP.png) I add the plot of mean log against file[,4] ![enter image description here](http://i.stack.imgur.com/zcpb8.png) I add the plot of mean reciprocal against number of active developers ![enter image description here](http://i.stack.imgur.com/lKQJF.png) ![enter image description here](http://i.stack.imgur.com/AwMMr.png)


title: sequential testing and hierarchical tesing

text: With respect to the hypothesis testing, I once heard sth such as sequential testing, hierarchical testing, and multi-level testing. But I could not find a good reference discussing these topics. Any recommendations are highly appreciated.


title: Does the degrees of freedom depend on whether the variable is continous or dummy variable?

text: In my dataset there is a lot of multicollinearity between the independent variables, therefore it is not possible for me to include them all in the same regression (especially because I have a small sample size). However, my supervisor said that if I categorize some of the variables I might be able to use more independent variables in the same regression. Why is this? Does it affect the degrees of freedom? Or what is the reason?


title: ANOVA with 3 groups- does adding a group close to the mean reduce power?

text: Sorry for that horrible title. I have run a psychology experiment where participants received one of three sets of instructions (between-subjects). Either they were told to approach the experiment a certain way, told to approach it a different way or it was not specified how they should approach the task (it turned out that some in this group did it one way and some the other, about 60/40). My ANOVA on performance in the experiment is nonsignificant, but my comparison between the 2 groups which received specified instructions were significant either as a planned comparison or in an independent t-test (the effect size between the two is pretty good, Cohen's d = ~0.45). I've been taught not to interpret ANOVA comparisons unless the overall effect is significant, but it seems to me that including this middle group (the means of which tend to be between the two groups) is underpowering the ANOVA. Is it "cheating" to just compare the two groups with a t-test? It certainly doesn't seem right to conclude there is no difference between groups when the ANOVA only breaks with inclusion of this middle group.


title: How can I change the title of a legend in ggplot2?

text: I have a plot I'm making in ggplot2 to summarize data that are from a 2 x 4 x 3 celled dataset. I have been able to make panels for the 2-leveled variable using `facet_grid(. ~ Age)` and to set the x and y axes using `aes(x=4leveledVariable, y=DV)`. I used `aes(group=3leveledvariable, lty=3leveledvariable)` to produce the plot so far. This gives me a visualization that is paneled by the 2-leveled variable, with the X axis representing the 4 leveled variable and different lines plotted within the panels for the 3-leveled variable. But the key for the 3-leveled variable is titled with the 3-leveled variable's name and I want it to be a title that has a character space in it. How can I rename the title of the legend? Things I've tried that don't seem to work (where `abp` is my ggplot2 object):                abp <- abp + opts(legend.title="Town Name")      abp <- abp + scale_fill_continuous("Town Name")      abp <- abp + opts(group="Town Name")      abp <- abp + opts(legend.title="Town Name")      Example data:               ex.data <- data.frame(DV=rnorm(2*4*3), V2=rep(1:2,each=4*3), V4=rep(1:4,each=3), V3=1:3)


title: Unexpected under-dispersion in Perl simulations of Poisson RV

text: The following piece of Perl code randomly maps a set of ranges onto a circumference of a circle. In the example, the circumference is of length 1000 and legal ranges are e.g. (0,8)=0,1,2,...,8 and (995,2)=995,996,...,999,0,1,2 (i.e. zero-based coordinates; both start and end are inclusive). I take some arbitrary position on the circumference (e.g. 36) and count how many ranges cover it in each simulation. finally, I calculate the mean and variance of this statistic.               use strict;     use warnings;          use Statistics::Descriptive;          my $n_simulations   = 1000;     my $circumference   = 1000;     my @lengths_distrib = (100) x 100;    # distibution of range lengths     my $some_pos        = 36;             # arbitrary position          my $stat = Statistics::Descriptive::Full->new();          foreach my $sim ( 1 .. $n_simulations ) {           # randomly map ranges onto circumference           my @random_ranges =       map { my $start = int( rand($circumference) ); [ $start, ( $start + $_ -1 ) % $circumference ] }       @lengths_distrib;           # count how many range contain $some_pos      my $num_covering_ranges = scalar(       grep { ( $_->[0] <= $some_pos and $_->[1] >= $some_pos ) or ( $_->[1] < $_->[0] and $_->[1] > $some_pos ) }        @random_ranges      );      $stat->add_data($num_covering_ranges);     }     print $stat->mean, ' ', $stat->variance, "\n";      To the best of my knowledge, this kind of random variable should follow Poisson distribution (law of rare events and so on). Hence, the mean and variance should be equal. However, the variance seems to systematically be a bit lower than the mean. What am I missing? **UPDATE** Following whuber response, a couple of notes and additions:   1. _the Poisson distribution is a good approximation of the binomial distribution if n is at least 20 and p is smaller than or equal to 0.05, and an excellent approximation if n ≥ 100 and np ≤ 10._ source: wikipedia. The example I gave here is a toy one. I usually use much larger circumferences (length ~ 3M), many more ranges (around 25k) an relatively small ranges (size ~ 2.5k), so I think in general Poisson distribution should be appropriate. I will give it a look again on some real data.   2. In my real data, I have e.g. some 25k ranges of different sizes. Currently, I'm simulating random mapping of these clones on the circumference, count the coverage of some point of interest in each simulation (this could be an arbitrary fixed point since all point are equal), then get the mean of this coverage over some 100 or 1000 simulations. Now, I use a CDF of a Poisson distribution with this mean to get a P-value for the actual observation (number of ranges covering a specific point of interest). I take lower tail since I'm interested in under-covered positions. How would you suggest going around this? Should I use a Binomial distribution? How can I evaluate its parameters? Note this is not exactly a Binomial RV, since while I do have n independent trials, the chances for each trial to succeed is not equal (it depends on the range size). Perhaps I should simply take the average success probability (average length) and not even have to use simulations, right? In other words, I guess what I'm asking is how does the sum of independent Bernoulli RV with **different p's** look like? I should also note that while this example only deals with the number of ranges covering a single position, I'm also interested in the number of ranges covering some range. For now I use the simulations scheme exactly the same way described before. Thanks!


title: What is the Unscented Kalman Filter?

text: What is the Unscented Kalman Filter and when is it used in preference to other types of filters? edit: I find the Wikipedia explanation a bit too technical to be readily understood.


title: Clustering structured data: Assessing the similarity of documents that appear in tree structure

text: Usually when performing text document clustering, similarities across documents are assessed based on the lexical content of documents. But, in my problem, I wish to consider both the lexical content and the ontology of documents while assessing similarities. I'll explain this with an example. I have documents arranged in a tree structure as in the figure below. ![Document Tree](http://i.stack.imgur.com/8mJQ4.png). Here, lets say, documents _admod_ and _customevent_ (leaf nodes) have a lexical similarity X (X is a real number). It could be noticed that these documents lie close to each other in the tree. Lets say, their distance in the tree is Y. Hence, while assessing similarity between these 2 documents, I wish to have > Similarity (admod, customevent) = f(X,Y), where f is an appropriate function > (mapping) Extending this example further, say the documents _admob_ and _Task Manager_ (root node) bear the same lexical similarity X. It could be observed that these 2 documents lie far apart in the tree. Say, they are Z units apart (where Z > Y), then, I require the similarity function to satisfy the following constraint, > Similarity (admod, customevent) > Similarity (admod, Task Manager) as admob & customevent lie closer in the tree. Please suggest on how to combine (include) the tree structure while assessing the similarities. NOTE: The tree structure is inferred through domain expertise.


title: Who to follow on github to learn about best practice in data analysis?

text: It is helpful to study the data analysis code of experts. I've recently been perusing github and there are a number of people sharing data analysis code there. This includes a few R Packages (which of course are available directly from CRAN), but also several examples of reproducible research, particularly using R (see this R list on github).   * Who are good people to follow on github to learn about best practice in data analysis?   * Optionally, what kind of code do they share and why is this useful?


title: Formula that R "predict" function uses to calculate intervals exponential smoothing

text: I have been trying to figure out the exact formula that the R "predict" function uses to calculate prediction intervals for simple exponential smoothing. The prediction interval formula seems to vary according to the software used (Gretl is different from Minitab is different from SAS). I have tried to look at the "predict" code by simply typing it in without parentheses but R won't display it. I spent the better part of a day trying to figure out what R's formula. Does anyone know the exact formula that R uses? Also, I'd appreciate a link showing the derivation of the formula, as there seems to be some disagreement in the literature about what the "right" formula to use is. Thank you.


title: Finding correlation, triggers and patterns between two (or more datasets)

text: I am trying to find correlations between multiple data sets, for this example i am using two pieces of data, the closing stock price for Coca-Cola and Pepsi for the last year. I've been pointed in the direction of using the `=CORREL()` function in excel which has given me an indication that there is a correlation as it returns the value `0.869783523` but does not give me any clue as to triggers or patterns in the data which is what i am actually interested in. Is there function or a workflow that i can implement in excel to give me a more detailed read out giving me an idea of triggers and/or patterns?


title: Small sample linear regression: Where to start

text: FULL DISCLOSURE: This is homework. I have been provided with a small data set (n=21) the data are messy, looking at it in a scatterplot matrix provides me with little to no insight. I've been provided with 8 variables that are metrics created from a longditudinal study (BI, CONS, CL, CR, ..., VOBI). The other measurements are of mutual fund sales, returns, asset levels, market share, share of sales, and proportion of sales to assets ![alt text](http://i.stack.imgur.com/hNTBd.jpg) Correlations, are everywhere.                              BI       CONS           CL          CR         QT        COM        CONV        VOBI          s            r           a          ms         ss       share      share2     BI      1.0000000  0.7620445  0.639830594  0.70384322  0.7741463  0.8451500  0.84704440  0.85003686  0.2106773 -0.238431047  0.36184548  0.40007830  0.4076563  0.31643802 -0.28283564     CONS    0.7620445  1.0000000  0.933595967  0.96979599  0.9892533  0.9069803  0.96781703  0.93416972  0.2316209 -0.074351798  0.31952292  0.40259511  0.4442877  0.24783884 -0.14788906     CL      0.6398306  0.9335960  1.000000000  0.88297431  0.8993748  0.8133169  0.89922684  0.81132166  0.1200420 -0.001107093  0.22132116  0.26729067  0.3033221  0.07650924 -0.25595278     CR      0.7038432  0.9697960  0.882974312  1.00000000  0.9788150  0.8965754  0.92335363  0.90848199  0.2934774 -0.119340914  0.35973640  0.46409570  0.5012178  0.32832247 -0.09005985     QT      0.7741463  0.9892533  0.899374782  0.97881497  1.0000000  0.9216887  0.95458369  0.94848419  0.2826278 -0.108430256  0.35520090  0.43290221  0.4823314  0.31761015 -0.12903075     COM     0.8451500  0.9069803  0.813316918  0.89657544  0.9216887  1.0000000  0.90302002  0.89682825  0.4305866 -0.255581594  0.50724121  0.55718441  0.5773171  0.40378679 -0.12085524     CONV    0.8470444  0.9678170  0.899226843  0.92335363  0.9545837  0.9030200  1.00000000  0.96097892  0.1993837 -0.065237725  0.32010735  0.41843335  0.4531298  0.28873934 -0.19668858     VOBI    0.8500369  0.9341697  0.811321664  0.90848199  0.9484842  0.8968283  0.96097892  1.00000000  0.2424889 -0.087126942  0.30390489  0.40390750  0.4845432  0.36588655 -0.07137107     s       0.2106773  0.2316209  0.120041993  0.29347742  0.2826278  0.4305866  0.19938371  0.24248894  1.0000000 -0.173034217  0.91766914  0.84673519  0.8596887  0.61299987  0.32072790     r      -0.2384310 -0.0743518 -0.001107093 -0.11934091 -0.1084303 -0.2555816 -0.06523773 -0.08712694 -0.1730342  1.000000000 -0.22512978 -0.18337773 -0.1030943 -0.17650579  0.51768144     a       0.3618455  0.3195229  0.221321163  0.35973640  0.3552009  0.5072412  0.32010735  0.30390489  0.9176691 -0.225129778  1.00000000  0.92445370  0.8656139  0.63049461  0.03876774     ms      0.4000783  0.4025951  0.267290668  0.46409570  0.4329022  0.5571844  0.41843335  0.40390750  0.8467352 -0.183377734  0.92445370  1.00000000  0.9572730  0.77582501  0.08435813     ss      0.4076563  0.4442877  0.303322147  0.50121775  0.4823314  0.5773171  0.45312978  0.48454322  0.8596887 -0.103094325  0.86561394  0.95727301  1.0000000  0.83931302  0.24371447     share   0.3164380  0.2478388  0.076509240  0.32832247  0.3176102  0.4037868  0.28873934  0.36588655  0.6129999 -0.176505786  0.63049461  0.77582501  0.8393130  1.00000000  0.20313930     share2 -0.2828356 -0.1478891 -0.255952782 -0.09005985 -0.1290307 -0.1208552 -0.19668858 -0.07137107  0.3207279  0.517681444  0.03876774  0.08435813  0.2437145  0.20313930  1.00000000       Now, I've tried running a number of "tests", for example:               summary.lm(share2 ~ BI + ...)       However, none of them provide any reasonable result (mostly negative adjusted R^2). I'm wondering, if you had data where it looked like there was no relationships (linear at least). What would your next steps be? P.S: I did try a number of model formulas that contained interaction effects and received much better results (R^2 Ra^2 > 80% and significant f-tests) but not all the interaction effects where significant.


title: Mixed effect linear regression model output interpretation

text: just fitted mixed effects model: Linear mixed model fit by maximum likelihood ['lmerMod'] Formula: price ~ variable + (1 | product) Data: podzbior                     AIC       BIC    logLik  deviance  df.resid      130840.14 130868.85 -65416.07 130832.14      9674       Scaled residuals:                   Min      1Q  Median      3Q     Max      -6.2824 -0.3099 -0.0547  0.2201 12.4291       Random effects:                Groups                    Name     Variance Std.Dev.      product                  (Intercept) 427375   653.7         Residual                              25930   161.0                  Number of obs: 9678, groups: product, 1222      Fixed effects:                                             Estimate Std.    Error   t value      (Intercept)                    9.362e+02  1.899e+01     49.29       variable                     -7.521e-04  1.171e-04     -6.42      Correlation of Fixed Effects:                                    (Intr)              variable    -0.050      that was output from summary(lmerModel), before that durring lmer() I got this warning: Warning:                 In checkScaleX(X, ctrl = control) :       Some predictor variables are on very different scales: consider rescaling      Q1 Predictor variable is numeric from 0 to something like 100k, how It should be scaled? Random effects with confidence intervals chart for this model looks like this, is it ok?: ![](http://ramza.tarchomin.pl/randomEffect.png) Im pretty sure residuals, aren't ok, what should I do in this case? ![](http://ramza.tarchomin.pl/randomDisturbance.png) How can I go deeper with this model diagnostic, besides checking p-values? appreciate any replys


title: How to accomplish unsupervised separation of subpopulations?

text: I have a dataset drawn from a social network that looks Bimodal on logarithmic scales for all attributes (I'll demonstrate only one here): ![Bimodal Gaussian Distribution](http://i.stack.imgur.com/f4R4N.png) I know the variable that would give me a clean separation for the two subpopulations (e.g., gender): ![enter image description here](http://i.stack.imgur.com/UzOC6.png)![enter image description here](http://i.stack.imgur.com/CLFs7.png) However, I need to come up with the same separation of the two clusters without using **that** binary variable (I could use internal knowledge of the two sub-populations statistics such as averages, stds etc .. ) in an unsupervised manner. I tried using Gaussian Mixture Models with EM, however, it didn't perform well on this dataset. Clearly what seems normal on logscale it isn't in reality :P I am looking for a package in python, R or Matlab that would be able to assign individual points to those two clusters.


title: Odd result for intercept in ANCOVA

text: I try to perform an ANCOVA for data with two factors and one covariate. The first factor has 2 levels and the second has 4.                   Coefficients:             Estimate Std. Error t value Pr(>|t|)           (Intercept)   7.9342    3.1123   2.704   0.0311 *         k[T.2]      -4.3127     1.8957  -1.248   0.2122           k[T.3]       -2.453     1.6786  -1.557   0.1432           hg[T.2]      2.2411     2.4353   0.891   0.3219           hg[T.3]      4.3456     2.2342   1.668   0.1263           hg[T.4]     -3.5233     2.2318  -0.682   0.5108        I seem to be getting 2 coefs for the first factor and 4 for the second. I would like to note that they are of type "factor" so that R would not treat them as numeric. Why is this happening and is it normal?


title: How can I find the coefficient of correlation between ordered residuals and their expected values under normality?

text: As the question says I have to find the coefficient of correlation between ordered residuals and expected values under normality in Minitab. How can I do so? I have just started using Minitab and am a novice at it.


title: What is this horizontal point & confidence graph called?

text: From page 166 of the Open Intro to Stats, what is this type of graph called? It represents 25 point estimates of a mean, together with their confidence intervals, centered around the actual population mean. The red line is a point estimate whose 95% confidence interval fails to capture the true population mean. Further, is there a default R command to plot these graphs? Most graphs in the book seem to have been generated with R... ![line graph](http://i.stack.imgur.com/vaZ3L.png)


title: Binary and multinomial logistic regression

text: Is it correct to say that binary logistic regression is a special case of multinomial logistic regression when the outcome has 2 levels?


title: Is this Poisson regression specified correctly?

text: I want to ask some confusing in my mind regarding the data analysis: I have conducting the experimental trial in repeated measure by hours as below:               Design: Completely Randomized Design--(in green houses)-there are 6 box unit     ========================================================     - 1 Treatment (2 levels): level A and B     - Replicate: 3 replication     - Response was recorded by hours (repeated measure) from 24 hours (poisson distribution)      This follows as linear mixed model. and I used $R$ to process my data by using this syntax:               > Random intercepts for 3 rep but time and treat effects are fixed      > model1 <-lmer(respon~ treat * time + rep + (1|rep), data=data, family = poisson)        OUTPUT: > >     model1 >       Generalized linear mixed model fit by the Laplace approximation >       Formula: respon ~ treat * time + rep + (1 | rep) >       Data: data >       AIC BIC logLik deviance >       720 740   -354      708 >       Random effects: >       Groups Name        Variance  Std.Dev. >       rep    (Intercept) 0.0078266 0.088468 >       Number of obs: 207, groups: rep, 3 >                   Fixed effects                  Estimate Std. Error z value Pr(>|z|)          (Intercept)  0.843199   0.311791  2.7044  0.00684 **       treat       -0.246484   0.124762 -1.9756  0.04820 *        time         0.010138   0.016759  0.6049  0.54523          rep          0.063730   0.082557  0.7719  0.44015          treat:time   0.009520   0.008028  1.1858  0.23569          ---       Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       **MY QUESTION IS: is my syntax above TRUE?** Thanks in advanced


title: Machine learning - Weighting the regularisation function

text: I am working on a machine-learning algorithm for estimating the performance of a student on various tasks. I'm using logistic regression with a L2 regularisation parameter. One parameter we used was a exponentially weighted moving average of the last 20 problem with exponent 0.3. We then tried some combinations of multiple weighted moving averages and found that using 4 or 2 different weighted averages as parameters actually made the performance worse. I suspect that this is because the other factors are better predictors so adding in extra weighted moving averages reduces their impact. We have a reasonably large data set, so I don't believe that it is related to overfitting. Now, I would expect that it should be possible to add in these extra factors in a way so that it should improve performance (please note that I don't actually intent to use these factors - I'm just asking this question out of interest). Like if we doubled the cost of the weights for the moving averages, then I would expect that this should stop them crowding out more useful effects and in most cases lead to the performance being no worse than a single feature. My question is how come no-one seems to do this? Every paper I see seems to just use standard L1 or L2 regression. This seems like a pretty straightforward technique for getting better performance from features that contain less information.


title: Conditional probabilities/expectations

text: A coin minting machine randomly produces unbalanced coins so that the probability of getting a head in tossing a coin is a random variably $Y$. Supposed $Y$ has a pdf $f(y) = 2y$ for $0 <= y <= 1$ and $0$ otherwise. Randomly take one coin.   1. Toss this coin, and let $X$ be 1 if the outcome is a head, and 0 if a tail. Find the probability $P(X=1)$.   2. Toss the coin $n$ times. Find the probability of getting $k$ heads in the $n$ tosses, where $n$ and $k$ are positive integers, and $k <= n$.   3. Toss this coin twice. If the first toss results in a tail, what is the conditional probability that the next toss is also a tail? Hint for (3): use the Beta function Edit: I've having trouble conceptualizing this question. For (1) Isn't P(X=1) equivalent to Y? What are the strategies behind solving (2) and (3)? For 3) isn't the second coin toss independent of the first toss? Why is there a conditional probability?


title: What is a robust way to find the max of $n$ independent, non-identical random variates?

text: Suppose I observe $n$ random variates along with their variance (but not mean) and I'd like to select the one with the largest mean as frequently as possible. The procedure must be memoryless--you cannot keep track of how often you selected each as the max nor can you keep track of any other statistics. What is the best way to do this? For example, suppose I draw from $n=5$ Normal distributions. Each has $\mu_i=0$ and $\sigma^2=\tfrac{1}{1+i}$. In this case I would like to have a scoring rule such that every variate has an equal chance of being chosen as the max--in short, I want to account for the risk present in some of the variates. If we were to always just select the max, then the lowest variance variable, ie, $i=5$ would be selected the least frequently. If $n=2$ we wouldn't have to worry about this problem. More generally it should have the property that variables with a higher mean get selected more often. There is obviously going to be some tradeoff between mean and variance, but I don't know what that should be. It seems like the answer to this problem has to do with the relative variances as well as the number of observations per ranking. One simple idea (which ignores $n$) is just use $x_i - \sigma_i^2$ as a score for ranking. Is there something better? How would I prove optimality? I've tried playing around with LINEX and Bayesian Decision Theory. I feel the answer lies in the asymmetric loss present in over-prediction (since $n>2$ implies there are more ways to get it wrong then right). I've also played around with the idea of using (in the above example) a Normal-Normal posterior. This has the effect of lowering the variance (making it more fair) at the expense of adding some bias (making it less fair). Is there something better out there? Is this a reasonable approach? What if I used LINEX loss under this posterior to derive an optimal rule? Note that this question is not just: what is the pdf of non-identically distributed random variates (or some variation). It is more about how to minimize the risk of choosing an outlier when we are looking for the "max" (which is clearly not robust). I feel this problem must come up in ranking whenever the score is calculated from a model trained with a univariate loss function (such as ranking using logistic regression).


title: Fast algorithm for variable selection

text: The (training) data contains 1280 observations with 1415 features. The test set has additional 380 observations. The data is sparse, that is, many of the variables has many zeros and few positive values. a typical variable\feature distributes with zeros in about 90% of the observations and small tail of few positive values The response is dichotomous, i.e zero-one classification. I have tried to fit:   * L1,L2 and mixtures of both regularized binomial(logistic) models with `cv.glmnet` from `glmnet` package.   * A simple `glm()` model with 100 first selected variables after selecting them with stepwise forward procedure `step()`. surprisingly the second fit (hardthresholded model) preformed better on the training and test sets, measuring to goodness of fit with auc (area under the curve in the ROC curve) and misclassification error rate. It even preformed better than a model with first 100 variables selected by L1 regulation.(any idea why?) For this reason and because of some regulatory conditions in the industry I would like to continue preform variable selection similar to stepwise algorithm but faster. Any suggestions?


title: Mixed model runs well in R whereas a random effect has only one level

text: How do you explain that ? There's only one operator but the mixed model returns an estimate for the `operator` random effect. Furthermore the `sample` effect is confounded with the interaction `sample:operator`. Below is the R code.               > dd        sample operator         y     9      10      SCF 0.9153188     10     10      SCF 0.9884982     19    100      SCF 2.0798781     20    100      SCF 2.0464027     29   1000      SCF 3.0401590     30   1000      SCF 3.0114448     39  10000      SCF 4.1348324     40  10000      SCF 4.0840063     49  1e+05      SCF 5.1235795     50  1e+05      SCF 5.1106381     59  1e+06      SCF 6.0803404     60  1e+06      SCF 6.2353263     > str(dd)     'data.frame':   12 obs. of  3 variables:      $ sample  : Factor w/ 6 levels "10","100","1000",..: 1 1 2 2 3 3 4 4 5 5 ...          $ operator: Factor w/ 1 level "SCF": 1 1 1 1 1 1 1 1 1 1 ...      $ y       : num  0.915 0.988 2.08 2.046 3.04 ...     > lmer(y ~ (1|sample)+(1|operator)+(1|sample:operator), data=dd)      Linear mixed model fit by REML      Formula: y ~ (1 | sample) + (1 | operator) + (1 | sample:operator)         Data: dd        AIC   BIC logLik deviance REMLdev      18.6 21.03 -4.302    9.932   8.605     Random effects:      Groups          Name        Variance   Std.Dev.      sample:operator (Intercept) 1.87954740 1.370966      sample          (Intercept) 1.87954925 1.370967      operator        (Intercept) 0.00063096 0.025119      Residual                    0.00283931 0.053285     Number of obs: 12, groups: sample:operator, 6; sample, 6; operator, 1          Fixed effects:                 Estimate Std. Error t value     (Intercept)   3.5709     0.7921   4.508      For those who are more familiar with SAS the corresponding code is:               PROC MIXED DATA=dd;     CLASS sample operator;     MODEL y=;     RANDOM sample operator sample*operator;     RUN;      This is nothing but the crossed 2-way ANOVA with random effects.


title: Linear contrast OLS - Simultaneous

text: I fit the following linear model using OLS, where $X_{5}$ and $X_{6}$ are both dummy variables and the rest are continuous: $Y=\beta_{0} + \beta_{1}X_{1}+ \beta_{2}X_{2} + \beta_{3}X_{3} + \beta_{4}X_{4} + \beta_{5}X_{5} + \beta_{6}X_{6} + u$ If I want to test a simple contrast, say that $\beta_{5}$ = 0 I can do the following (or look at the print out from any regression program, but just for illustration sake): Null Hypothesis: Ho: $\beta_{6}(1) - \beta_{6}(0)$ =0 Using the Car package in R:               cVec<-c(0,0,0,0,0,1,0)     car::linearHypothesis(mod,cVec, verbose=TRUE)      If I want to test that both $\beta_{5}$ = 0 and $\beta_{6}$ = 0 and set up the vector to multiply against the coefficient vector as follows: `cVec<-c(0,0,0,0,0,1,1)` the hypothesis test is not correct and one actually needs to use a matrix for this simultaneous test:               cMatrix<-t(cbind(c(0,0,0,0,0,1,0),c(0,0,0,0,0,0,1)))      car::linearHypothesis(mod,cMatrix, verbose=TRUE)      Why and what is `cVec<-c(0,0,0,0,0,1,1)` actually testing if not if there is a difference between $\beta_{0} + \beta_{1}X_{1}+ \beta_{2}X_{2} + \beta_{3}X_{3} + \beta_{4}X_{4} + \beta_{5}(1) + \beta_{6}(1)$ AND $\beta_{0} + \beta_{1}X_{1}+ \beta_{2}X_{2} + \beta_{3}X_{3} + \beta_{4}X_{4} + \beta_{5}(0) + \beta_{6}(0)$


title: Software quality control references

text: I'm looking for references on statistical quality control of software. My searches have turned up a lot of writing on how to predict the time that software projects will take, and also software which can be used to in quality control, but I haven't really found anything specifically for evaluating software quality. Some basic questions I'd want to ask are things like:   * Does department X create more bugs than department Y?   * Does "crunch time" cause an increase in the defect rate? In investigating these, I'm finding that data is not normally distributed, which makes me concerned that my limited stats knowledge is up to the task without help. Also, some of my problems (e.g. the non-normal distribution) apply to quality control in general, but I'm worried about "unknown unknowns."


title: How to cluster standard errors (by country) using SPSS?

text: I'm trying to obtain robust standard errors as my regression residuals seem to be correlated... My sample data has 8 countries and I would like to cluster the standard errors by these 8 countries using SPSS... I have read that you can use STATA to do it but I only have SPSS installed in my laptop.. help.. does anyone know any short cuts or macros written for this? Thank you!!!


title: Expectation notations

text: In Statistical Decision Theory, one often studies the following two measures (from "The Bayesian Choice"): **Average loss** (aka the _frequentist_ risk): $R\left(\theta,\delta\right) = \mathrm{E}_\theta\left[L\left(\theta,\delta(x)\right)\right] = \int_X L\left(\theta, \delta(x)\right)f(x|\theta)dx$ **Posterior expected loss** : $\rho(\pi, d|x) = \mathrm{E}^\pi\left[L\left(\theta,d\right)| x\right] = \int_\Theta L\left(\theta,d\right)\pi\left(\theta|x\right)d\theta$ I am confused by this notation: What do **subindices** (e.g. $\theta$ in $\mathrm{E}_\theta$), **superindices** (e.g. $\pi$ in $\mathrm{E}^\pi$) and **conditions** (e.g.$x$ in $L\left[.|x\right]$ ) represent when defining **expectations**? For reference, in the formulas above:   * $\delta(x)$ is known as the _decision rule_ (i.e. the allocation of a decision to each outcome x $\sim$ $f(x|\theta)$   * The value $\delta(x)$ is also known as the estimate of $\theta$   * The function $\delta$ is known as the _estimator_   * $\pi$ is the posterior distribution of $\theta$ given $x$


title: Reject or retain the null hypotheses with insufficient data

text: I stated a number of hypotheses and conducted a questionnaire. Only after the questionnaire it appeared I skrewed up with the questions. These questions don't actually allow me do reject or retain the set hypotheses. So my question is what do I do? Do I reject or retain?


title: Difference of Prediction between Graph Representation and Data Matrix Representation

text: In **data mining** , data can be usually represented in different forms such as records of a **matrix** , **graphs** or **ordered data**. While we find in research different papers addressing methods or solutions for these different representations, there is no clear description of the advantages of every representation compared to the other ( **i.e when different representations can be applied to solve the problem, under which conditions, a particular representation would have an advantage?** ). Here, I am interested in knowing what is, in particular, the advantage of the graph representation over the data matrix representation and vice versa. I realize that different problems would have an intuitive representation as one of the two ways. For example, a social network, is intuitively, represented as a graph while patients records, is intuitively, represented as a data matrix. However, I want to know how these representation compare when there is a prediction task and both representation can be used to solve the task. An example that may illustrate my interest is chemical-protein interaction network. In this network, chemicals that may have an effect on a specific protein target, will have an active relationship. This active relationship can be represented either as an edge of weight 1 between a chemical and a protein in a graph or as a positive label for a set of features describing the compound in a record. Another example would be the author-paper network. To predict the author of the paper, we may extract features from the papers and build our data matrix. Another way, would be building a graph in which a new paper is linked with the most similar paper and then, we try to predict who may be the authors based on traversing the graph. **One answer I would have once thinking about these two representations is the different levels of describing the data.** In a data matrix, there is the advantage of having many variables describing a given case or sample. In graphs, on the other hand, it is only one variable representing the similarity between the samples. Nevertheless, a graph topology may highlight important nodes in the network. What else ? **In summary, I am interested to know the expert advise on when to use a graph representation or a data matrix representation and why? If you are someone who likes graphs and prefer mining them, tell me why?**


title: Discriminant analysis for the validation of cluster analysis

text: I did a discriminant analysis for the validation of my cluster analysis. The cluster analysis is based on a PCA, so I used the components as the independant variables in the discriminant analysis. My question is: is there a restriction considering the number of the independant variables (components) compared to the number of my (cluster) groups? Thanks a lot!


title: Covariance and correlation matrix comparison

text: I am aware that this question may be too broad and that answers are scattered in various posts, but i need concise and organized answer. My dataset consists of linear measurements of cranial dimensions on 600 individual roe deer (50 distinct measurements with a dial caliper). I divide this dataset to unequal groups (corresponding to population membership), and caluclate correlation or covariance matrix so that every population is represented by 50x50 matrix. My question is, what is the best way to compare those matrices, both for equality and pattern (excluding Mantel test)? Problematic part may be the fact that those matrices are rarely of full rank since many measured characters are significantly correlated. Also, that comparison should include some kind of confidence intervals. **Edit:** In the meantime I have found one possible solution, I just need to implement it in R code. This paper suggests possible distances based comparisons that I really need. The _Euclidean distance_ as a simple method like this (Si is a sample covariance matrix): $d_e(S_1,S_2) = \sqrt {tr((S_1-S_2)^t(S_1-S_2))}$ which I implemented in R code like this (although unsure):               covDif <- sqrt(t(cov(malesMab)-cov(malesMbm))*(cov(malesMab)-cov(malesMbm)))     sqrt(sum(diag(cov(covDif))))      But this distance is not good for comparison and of all the distances suggested in the mentioned paper the Cholesky decomposition is the best but I don`t know how to program it in R. This is its form: $d_e(S_1,S_2) = chol(S_1)-chol(S_2)$ which I tried like this (just substituting in the upper equation)               covDif <- sqrt(t(chol(cov(malesMab))-chol(cov(malesMbm)))*(chol(cov(malesMab))-chol(cov(malesMbm))))     sqrt(sum(diag(cov(covDif))))      which works, but encounters rank deficiency problems which I hoped to avoid by using Cholesky decomposition. Any suggestions?


title: Frequency Distribution

text: I have a question that is very important to me related to the book _Basic Statistics for Business and Economics_ for organizing data into a frequency distribution: > Step 1: Decide on the number of classes. The goal is to use just enough > groupings or classes to reveal the shape of the distribution. Some judgment > is needed here. A useful recipe to determine the number of classes ($k$) is > the "2 to the $k$ rule". This guide suggests you select the smallest number > ($k$) for the number of classes such that $2^k$ is greater than the number > of observations ($n$): [$n \le 2^k$] I want to know, how can I prove this formula?


title: What kind of clustering algorithm for clustering samples according to OTUs

text: I have data in the following form:               OTU Sample1 Sample2 Sample3 ...     1   1       0       2     2   3       5       0     3   0       5       1     .   .       .       .     .   .       .       .      So samples are sequencing data. Then all the sequences are pooled together and OTUs (operational taxonimic units) are calculated by an algorithm. Now I want to see which samples are close to each other by clustering using OTU information as a "fingerprint". So in the example above Sample1 has 1 sequence in OTU1, Sample2 has five sequences in OTU2 and so on. What kind of clustering algorithm would you suggest (so far I understand hierarchical clustering and k-means, but I'd like to learn others if needed). I am not sure what fits here the best I would like to work with R. Some small side question: Does anyone know an algorithm called "Jacks knife" or so? If yes, where do I find info regarding this algo?


Can I compare AIC values of a linear function with a non-linear function? Because I get totally different results. One is 4000 other the 6000000. Estimation is done on the same data setvariables.


title: Modelling probabilties within friend sets

text: Recently, I was wondering about calculating the probability of a given individual in a given population "knowing" (let's say, present in individual's friend set) at least one person with a given trait A and at least one person with another given trait, B; where it is possible that any number of people in the given individual's friend set can possess both traits. For example, using genetic traits, in a given population, how could one calculate the probability that a given individual in a given population "knows" at least one person with grey eye colour and at least one person who is greater than 200cm tall; where, naturally, it is possible that any number of people in the friend set can possess grey eye colour _and_ be greater than 200cm tall. I have developed a sort of model, but it may not be correctly specified; it is as follows: **Assumptions and Qualifications:**   * First of all, for simplicity, let's assume that we define "knowing" as mutually connected friends on an online social network.   * Secondly, the frequencies of genetic traits are (a) going to be determined by ethnicity of a given population as well as environmental factors (nutrition, healthcare) and (b) unlikely to be independently distributed across a given individual's friend set (e.g. family members will have greater genetic similarity); however, for this problem, let's adopt a simple model where both of the above conditions are violated.   * Thirdly, assume that the individual's friend set provides a microcosmic representation of society; this facilitates a frequentist generation of probabilities from the instance rate in the population.   * Finally, instances of genetic traits in the population are fabricated, but are used to generate probabilities for the examples. **Model Formulation:** I have reasoned that a binomial random distribution can be applied to the probability of each genetic trait, where a "success" is defined as an individual in the friend set possessing that genetic trait. Thus for trait A, we have: $$P(A=k) = \binom{N}{k}p_{a}^{k}(1-p_{a})^{N-k}$$ and, for trait B: $$P(B=k) = \binom{N}{k}p_{b}^{k}(1-p_{b})^{N-k}$$ where $N$ is number of friends in the friend set, $k$ is number of people containing the genetic trait and $p_{x}$ is the probability associated with possessing genetic trait $X$. Because the scenario is concerned with the probability of _at least one_ person possessing trait A and _at least one_ person possessing trait B, it is easier to find the complement of no people in a friend set containing the trait; for both traits: $$P(A >= 1) = 1 - P(A = 0)$$ and $$P(B >= 1) = 1 - P(B = 0)$$ Furthermore, we know that probability of an intersection of events is given by: $$\mathbb{P}(A \cap B) = \mathbb{P}(A|B)\mathbb{P}(B) = \mathbb{P}(B|A)\mathbb{P}(A)$$ However, because we are assuming independence between genetic traits, $\mathbb{P}(A)$ and $\mathbb{P}(B)$ are independent, thus: $$\mathbb{P}(A \cap B) = \mathbb{P}(A)\mathbb{P}(B)$$ Combining the information above, we get the following model for the above scenario: $$\mathbb{P}(A >=1, B >=1) = \left(1-\binom{N}{k}p_{a}^{k}(1-p_{a})^{N-k}\right)\left(1-\binom{N}{k}p_{b}^{k}(1-p_{b})^{N-k}\right)$$ So, for the original example above, assuming a friend set of size $N = 300$, the instance of trait A in the population is $\frac{1}{800}$ and the instance of trait B in the population is $\frac{1}{5000}$; according to the model, we get the final probability: $$\mathbb{P}(A >=1, B >=1)$$ $$ = \left(1-\binom{300}{0}\left(\frac{1}{800}\right)^{0}\left(\frac{799}{800}\right)^{300}\right)\left(1-\binom{300}{0}\left(\frac{1}{5000}\right)^{0}\left(\frac{4999}{5000}\right)^{300}\right)$$ $$\approx 0.018 = 1.8\%$$ Does this model seem reasonable given the assumptions? Assuming this model is _not_ correctly specified, maybe somebody could provide a more accurate representation.


title: Using the rejection sampling with the method of inversion

text: I am hoping to write some rejection algorithm code in R to approximate a $\text{Gamma}(k,\lambda)$ distribution. The problem is more for educational purposes than real-world implementation. Given an $\text{Exponential}(\lambda)$ distribution with PDF: $$f_{\lambda}(x) = \lambda\text{exp}(-\lambda x)$$ the CDF can be expressed as: $$F_{\lambda}(x) = 1 - \text{exp}(-\lambda x)$$ where $x>0$ and $\lambda>0$. And the inverse of the CDF can be expressed as: $$F^{-1}(U) = -\text{log}(1-U)/\lambda$$ where $U$ is a random number generated such that: $U \sim \text{U}[0,1]$. Assuming that: $$X_{1}, X_{2}, \dots, X_{k} \stackrel{\text{ iid }}{ \sim }\text{Exponential}(\lambda)$$ and: $$Y = X_{1} + X_{2} + \dots + X_{k}$$ then: $$Y \sim \text{Gamma}(k, \lambda)$$ Now, the above approach for simulating a $\text{Gamma}(k,\lambda)$ from a sum of $k$ $\text{Exponential}(\lambda)$ random variables will, naturally, only work where $k \in {\bf N}$, but let's assume the objective is to simulate $\text{Gamma}(k,\lambda)$ where $k \notin {\bf N}$. It has been suggested to use an accept-reject approach with an envelope function being the $\text{Gamma}(\lfloor{k}\rfloor,\lambda-1)$ density, where $\lfloor{k}\rfloor$ is the function `floor(k)`. **EDIT:** In the above example, it's clear that target distribution is given by: $$f(x) = \frac{\lambda^{k}}{\Gamma(k)}x^{k-1}e^{-\lambda x}$$ But in order to solidify these concepts, could somebody provide the form of the simulating distribution ($h(x)$), the envelope function ($g(x) = M*h(x)$) and, of course, the optimal value of $M$? Solidifying these ideas would really help me out. ### Closely related Rejection sampling for a modified Gamma distribution How to quickly sample X if exp(X) ~ Gamma? (one answer provides a clever `R` implementation).


title: Understanding comparisons of clustering results

text: I'm experimenting with classifying data into groups. I'm quite new to this topic, and trying to understand the output of some of the analysis. Using examples from Quick-R, several `R` packages are suggested. I have tried using two of these packages (`fpc` using the `kmeans` function, and `mclust`). One aspect of this analysis that I do not understand is the comparison of the results.               # comparing 2 cluster solutions     library(fpc)     cluster.stats(d, fit1$cluster, fit2$cluster)      I've read through the relevant parts of the `fpc` manual and am still not clear on what I should be aiming for. For example, this is the output of comparing two different clustering approaches:               $n     [1] 521          $cluster.number     [1] 4          $cluster.size     [1] 250 119  78  74          $diameter     [1]  5.278162  9.773658 16.460074  7.328020          $average.distance     [1] 1.632656 2.106422 3.461598 2.622574          $median.distance     [1] 1.562625 1.788113 2.763217 2.463826          $separation     [1] 0.2797048 0.3754188 0.2797048 0.3557264          $average.toother     [1] 3.442575 3.929158 4.068230 4.425910          $separation.matrix               [,1]      [,2]      [,3]      [,4]     [1,] 0.0000000 0.3754188 0.2797048 0.3557264     [2,] 0.3754188 0.0000000 0.6299734 2.9020383     [3,] 0.2797048 0.6299734 0.0000000 0.6803704     [4,] 0.3557264 2.9020383 0.6803704 0.0000000          $average.between     [1] 3.865142          $average.within     [1] 1.894740          $n.between     [1] 91610          $n.within     [1] 43850          $within.cluster.ss     [1] 1785.935          $clus.avg.silwidths              1          2          3          4      0.42072895 0.31672350 0.01810699 0.23728253           $avg.silwidth     [1] 0.3106403          $g2     NULL          $g3     NULL          $pearsongamma     [1] 0.4869491          $dunn     [1] 0.01699292          $entropy     [1] 1.251134          $wb.ratio     [1] 0.4902123          $ch     [1] 178.9074          $corrected.rand     [1] 0.2046704          $vi     [1] 1.56189      * * * My primary question here is to better understand how to interpret the results of this cluster comparison. * * * Previously, I had asked more about the effect of scaling data, and calculating a distance matrix. However that was answered clearly by mariana soffer, and I'm just reorganizing my question to emphasize that I am interested in the intrepretation of my output which is a comparison of two different clustering algorithms. _Previous part of question_ : If I am doing any type of clustering, should I always scale data? For example, I am using the function `dist()` on my scaled dataset as input to the `cluster.stats()` function, however I don't fully understand what is going on. I read about `dist()` here and it states that: > this function computes and returns the distance matrix computed by using the > specified distance measure to compute the distances between the rows of a > data matrix.


title: Propagation of large errors

text: Everything I've read about combining errors in quadrature when multiplying or dividing quantities with associated errors says that this works for "small error". What about large error? Say I want to compute A/B where A is +/- 1% and B is +/- 50%, can I still reasonably add the errors in quadrature?


title: How to decide between a logistic regression or conditional logistic regression?

text: I have a case-control study in which the cases are firms with health insurance and the controls are firms with no health insurance. I am studying the factors affecting enrolment in health insurance and was therefore using a logistic regression, which includes several covariates on firm characteristics that were measured in a survey. I have randomly sampled the firms from a database that includes two strata: insured and uninsured firms. I selected 65 from each group. However, within the group I also sampled from four strata that correspond to industry. I am therefore wondering if I need to use conditional logistic regression, as opposed to unconditional logistic regression. However, I was under the impression that conditional logistic regression was for matched case-control studies or panel studies. In other feedback I've received I've been told that because I sampled on the outcome, I need to use the conditional model. Could someone please help me figure out which m? Any references would also be much appreciated. Thank you.


title: When/why should "rabbit ears hunting" density estimation "work"?

text: I came by Bernard Silverman concept of "rabbit ears hunting" interactive procedure for density estimation. Does anybody know:   1. A recommended reference for this procedure?   2. An explanation/justification for when it should "work"?   3. Any other suggested ad-hoc interactive criteria for density estimation?


title: Resources for beginners - how to determine probability of user action based on certain criteria?

text: Sorry if this is similar to another question, but I'm still trying to learn data analysis and I don't know what to search for. What I'd like to do is take some information that a potential customer enters into a web form and figure out the probability that they'll become a customer... giving our sales team an idea of which leads are most valuable. Specifically, being a travel company, there are criteria like origin and destination, date of travel vs. date they fill the form, and even how much/little they fill out the 3 subjective questions (about 100-150 characters seems ideal). I'm not even sure what the terminology is, so if anyone can provide some direction about where I can go to learn about this, I'd be very appreciative. Thanks!


In general common rejection region of t-test, there are three possible alternative hypotheses and rejection regions for the one-sample t-test:   1. For alternative hypothesis $H_1: μ1 ≠ μ2$, the rejection region is $|T| > t_{α/2},ν$   2. For alternative hypothesis $H_1: μ1 > μ2$, the rejection region is $|T| > t_{1-α},ν$   3. For alternative hypothesis $H_1: μ1 < μ2$, the rejection region is $|T| > t_{α},ν$ my question is, WHY for alternative hypothesis $(2) H_1: μ1 > μ2$ the criterion $t$ for rejection is using **$t_{1-α}$** than **$t_α$**???? thanks in advance!


title: Creat a confidence interval for non normal data

text: I have a data set based on percent time a fleet of machines were available each month, as a result the majority of the observations are close 100 with a few being lower and one or two being 0 this creates a J distribution which whose effect can be seen in the attached graph. ![enter image description here](http://i.stack.imgur.com/frU56.png) My question is therefore: As the distribution is not normal, what methods can I use to state with 95% confidence that the aggregate availability of the whole fleet in any given year will be X% or higher, given that the data is non normal? One solution might be to use bootstrapping. e.g. For a 10 machine fleet   1. Take 10*12 random samples find the mean for the fleets ‘machine   2. months’ repeat lots of times    3. find the confidence interval of the mean. Would this give a result that was meaningful? The assumptions that I am making are   * The machines are all identical.    * There are no significant downtimem differences between machines.    * There are no significant downtime differences between months.


title: Coefficient sign changes in fixed effect and first-difference estimation

text: I have a large empirical panel, where I basically want to regress the standard deviation of (equity) returns ($y_{it}$) of firm i at time t, on leverage (equity/debt) of firm i at time t ($x_{it}$). Further I would like to control for the equity size of the firm ($z_{it}$). Unfortunately the $\beta$ coefficient sign changes (significant in both specifications) if I run the following two models: Fixed Effects Model: $$ y_{it} = a_{i} + \beta * x_{it} + c * z_{it} + u_{it}$$ First-Difference Model: $$ \Delta y_{it} = \beta * \Delta x_{it} + c * \Delta z_{it} + \Delta u_{it}$$ Since I only have basic knowledge about econometrics and my colleagues are also puzzled, I was hoping that some statistic guys might help me. I identified some potentially problems (but there might be more):   * $y_{it}$ might be autocorrelated   * Since leverage is defined as Debt/Equity, controlling for equity brings some potential problems The true question now is: Which model (if any) is econometrically correct?


title: On the tractability of posterior distributions

text: I am trying to understand what makes estimating the posterior distribution such a hard problem. So, imagine I need to estimate the posterior distribution over a set of parameters given the data y, so a quantity $P(\theta|y)$ and $\theta$ is generally high dimensional. The prior over $\theta$ is a multivariate Gaussian i.e. $P(\theta) \sim N(\theta; 0, \Sigma)$ The likelihood i.e. $P(y|\theta)$ can be written down as product over Gaussian likelihoods. Now, it seems to be that the posterior distribution will also be Gaussian. Is that correct? Secondly, going through Bishop's book, it seems that the conditional posterior distributions and the marginal distributions will be Gaussian as well (assuming that the joint distribution over the parameters and data is Gaussian) and should have a closed form solution. If that is the case, why is this problem intractable? If I need to find the parameters of this posterior distribution, can this not be set as an optimisation problem where I estimate the mean and covariance of the posterior Gaussian? I am basically having trouble visualising why this problem is complicated?


title: Normality of dependent variable = normality of residuals?

text: This issue seems to rear its ugly head all the time, and I'm trying to decapitate it for my own understanding of statistics (and sanity!). The assumptions of general linear models (t-test, ANOVA, regression etc.) include the "assumption of normality", but I have found this is rarely described clearly. I often come across statistics textbooks / manuals / etc. simply stating that the "assumption of normality" applies to each **group** (i.e., categorical X variables), and we should we examining departures from normality for **each group**. _**Questions_** :   1. does the assumption refer to the **values** of Y or the **residuals** of Y?    2. for a particular **group** , is it possible to have a strongly non-normal distribution of Y **values** (e.g., skewed) BUT an approximately (or at least more normal) distribution of **residuals** of Y?  Other sources describe that the assumption pertains to the residuals of the **model** (in cases where there are groups, e.g. t-tests / ANOVA), and we should be examining departures of normality of these residuals (i.e., only one Q-Q plot/test to run).   3. does normality of residuals for the **model** imply normality of residuals for the **groups**? In other words, should we just examine the model residuals (contrary to instructions in many texts)? To put this in a context, consider this hypothetical example:     * I want to compare tree height (Y) between two populations (X).      * In one population the distribution of Y is strongly right-skewed (i.e., most trees short, very few tall), while the other is virtually normal      * Height is higher overall in the normally distributed population (suggesting there may be a 'real' difference).     * Transformation of the data does not substantially improve the distribution of the first population.   4. Firstly, is it valid to compare the groups given the radically different height distributions?   5. How do I approach the "assumption of normality" here? Recall height in one population is not normally distributed. Do I examine **residuals** for both **populations** separately OR residuals for the **model** (t-test)? * * * Please refer to questions by number in replies, experience has shown me people get lost or sidetracked easily (especially me!). Keep in mind I am not a statistician; though I have a reasonably **conceptual** (i.e., not technical!) understanding of statistics. P.S., I have searched the archives and read the following threads which have not cemented my understanding:   * ANOVA assumption normality/normal distribution of residuals   * Normality of residuals vs sample data; what about t-tests?   * Is normality testing 'essentially useless'?   * Testing normality   * Assessing normality of distribution   * What tests do I use to confirm that residuals are normally distributed?   * normality distribution   * Normality of residuals vs sample data; what about t-tests?


title: Relation between reduced chi square and Akaike criterion

text: In choosing the best fit between different curves, supposing the errors are distributed normally, I can do the following things:   * Calculate the reduced $\chi_r^2=\frac{\chi^2}{N_b-N_p}$, where $N_b$ is the number of bins and $N_p$ is the number of parameters, for all the fits and choose the smallest.   * Calculate the Akaike Information Criterion $AIC=\chi^2+\frac{2N_p(N_p+1)}{N_b-N_p-1}$, for all the fits and choose the smallest. What is the relation (if any) between the two methods?


title: how to code and interpret categorical time varying variable in Cox PH model

text: I'm working on a project looking at the relationship between exposure to several different drugs on the risk of preterm delivery in a cohort of pregnant women with a particular disease. I'm confident that I need to use a time-dependent variable in a Cox PH model for the drug exposures as women started the drugs at various points in pregnancy, and some were not exposed at all. So far, I've coded the drug variable in two ways: 1) as a categorical variable listing the drug types (they were never taken concurrently, but sometimes sequentially) with no drug exposure as the reference group, and 2) as a binomial variable for each drug with a 1 for exposed time intervals and a 0 otherwise. My problem lies in interpreting the HRs. For the case where the variable is categorical, the HR are relative to the reference group. But for the second case, are they relative to all other exposure types (including no exposure), or are they still relative to no exposure as the presence of all of the other drug variables in the model deals with the effects of those separately? In other words, is the interpretation of the HR more or less the same between the two types of models, or is the second model different, and less interpretable? Please forgive me if this is a basic question! I've scoured the intertubes and my reference library and can't seem to find quite this same problem addressed. If you know of a good source I should read, please do point it out. I'll say though that all of the time-dependent Cox PH resources I've found seem to deal with binary exposures... Thanks!


title: Combining correlated measurements - negative weights

text: I'm trying to get the analysis of this (I hope) simple case correct in my head. Lets say I measure the same thing two different ways, on the same set of data. Call these measurements $y_1$ and $y_2$. They have Gaussian errors $s_1$ and $s_2$. But since they are measured on common data, they have a correlation coefficient $\rho$. Say $\rho$ is measured with a toy monte carlo or similar, though I don't think that matters. Just to make sure we are on the same page, this will generate a co-variance matrix that has $s_1^2$ and $\rho s_1s_2$ on the top row, and $\rho s_1 s_2$ and $s_2^2$ along the bottom row. It is easy to calculate a weighted average, $ \bar y$. Under certain circumstances it is possible to find that $y$ is outside the range $y_1, y_2$. This happens, in particular, when $\rho >\frac {s_1}{s_2}$. Since these measurements are on the same data, and are measuring the same thing, this is just not possible. In short, "opps." I must have messed up my measurements of $s_1, s_2$, or $\rho$. I feel comfortable with this. :-) Now, lets say I do exactly as above, but now I take it to the next level and measure a number of systematic errors associated with each method. Say, for the sake of argument, that the systematic errors for method 1 are not correlated with any of the systematic errors for method 2. After thinking about this for a while... well, I'm not entirely sure what to think. Do I have such a neat test to see if I've made a mistake anymore? I actually think I must, as you might imagine that I have systematic errors but they are vanishingly small. So, I can do something like decrease $\rho$ - because they now have uncorrelated error (the systematic error) and then re- run the test. This is pretty straight forward, I think: I use $\rho$ to split the statistical error into a correlated and uncorrelated part, add the uncorrelated parts in quad with the systematic errors, and then re-calculate $\rho$. But does the test still make sense? I fear I'm missing a fundamental concept here or I'd be able to determine what was going on. Things in my current problem domain get arbitrarily complex. I have to combine measurements of the same thing in uncorrelated data sets, but with substantial correlated systematic errors. The procedure I describe just above for adjusting $\rho$ applies here just as well (the statistical error is uncorrelated, but some of the systematic errors will be), but here it seems incorrect. Another clue to my missing fundamental concept. :-) I'd appreciate some help getting my thinking straight.


title: dose response and lethal dose 50 analysis

text: I am going to perform an experiment to test for the pathogenicity of several bacterial strains. For this I will infect several animals (e.g. 5 per dose per bacterial strain) with increasing doses of bacteria (e.g. 10^5, 10^6, etc) and then obtain a binary response variable, which basically is whether the animal survived or died after exposure to a given bacterial strain/dose combination. I will then determine the dose response curves and their respective lethal dose 50 (LD50). Before I start doing this experiment I have been trying to understand what is the best method to analyse this data. I have been thinking that the best approach for me would be to use a binomial GLM, however I have started to read a bit more of the ecotoxicology literature and found the R package drc (by Christian Ritz & Jens C. Streibig) that implements several models for non- linear regression. Given that my knowledge of nonlinear regression is low, I am now wondering which of the two methods would be better to analyse the data (drc or glm), so any help on making a decision would be greatly appreciated. Cheers, Ricardo


title: How can I pool posterior means and credible intervals after multiple imputation?

text: I have used multiple imputation to obtain a number of completed datasets. I have used Bayesian methods on each of the completed datasets to obtain posterior distributions for a parameter (a random effect). How can I combine / pool the results for this parameter ? * * * More context: My model is hierarchical in the sense of individual pupils (one observation per pupil) clustered in schools. I have done multiple imputations (using `MICE`in R) on my data where I included `school` as one of the predictors for the missing data - to try to incorporate the data hierarchy into the imputations. I have fitted a simple random slope model to each of the completed datasets (using `MCMCglmm` in R). The outcome is binary. I have found that the posterior densities of the random slope variance are "well behaved" in the sense that they look something like this: ![enter image description here](http://i.stack.imgur.com/J522M.png) How can I combine/pool the posterior means and credible intervals from each imputed dataset, for this random effect ? * * * **Update1** : From what I understand so far, I could apply Rubin's rules to the posterior mean, to give a multiply imputed posterior mean - are there any problems with doing this ? But I have no idea how I can combine the 95% credible intervals. Also, since I have an actual posterior density sample for each imputation - could I somehow combine these ? * * * **Update2** : As per @cyan's suggestion in the comments, I very much like the idea of simply combining the samples from the posterior distributions obtained from each complete dataset from multiple imputation. However, I should like to know the theoretical justification for doing this.