title: Test to determine significance for non normal data with different variance and unequal n?

text: How can I find out if the differences in length and width of cut marks left in fabric are significantly different between groups (knife types)? **Additional Information** Groups:               KnifeType1 n=15     KnifeType2 n=15     KnifeType3 n=11     KnifeType4 n=15     KnifeType5 n=6     KnifeType6 n=4     KnifeType7 n=4     KnifeType8 n=5      I want to test length and width separately, and all measurements are in mm. The data is not normally distributed, and the groups have very different variances; do I have to rely on Welch's t test? I'm trying to determine if the lengths and widths of cut marks would help you determine what type of knife was used-possibly by seeing if the mean lengths and widths are sig different from each other. I'm using SPSS, and have a year of undergraduate stats. **Update** I read that Welch's t test was pretty robust agains nonnormality, would it be okay to use that one and state that the normality assumption was violated?


title: Factor analysis to group questions into blocks?

text: So we have a questionaire of 32 questions (answer 1-7 on each question) where we present the results in grouped in four areas (Question 1-8 in one, 9-16 etc), since we think that those questions belong together. We are now trying to find out of these grouping are correct or if maybe some should switch between groups (say: question 12 should be in group 4 instead of 2). And on the other hand: If 4 groups aren't correct and we should use 3 or 5 instead? Can we use factor analysis for that? I have a feeling that is more used to answer the question: Can we replace our 32 questions with less and still get a good answer if I have understood it correctly.


title: Determining Relative Weights

text: I am looking for some recommendations and more specifics about how to do the following: Objective: To determine the weights of a number of stock valuation metrics. I am looking at doing this across 10 different industry sectors as I expect that each sector will merit a different weighting scheme. The data I have are: Monthly average metrics for each sector over approximately 20 years. I also have return data for each sector on a monthly basis. I am looking to do this in R Studio. Thank you for your help on this.


title: Model assumptions and diagnostics for proportional hazard regression model with frailty in R

text: I am wondering are there any assumptions that must be met in the proportional hazard regression model with frailty? I remember that in regular proportional hazard model without frailty all variables need to significant in log-rank test which means that the survival curves for each variable differ significantly and are proportional. Moreover I was thinking about model diagnostics. Is Schoenfeld test appropriate for model with frailties? Also do scored Schoenfeld's residuals (or martingale residuals) are a good measure of model fit like in model without frailty? Maybe there are no assumptions and diagnostics for this kind of model? And quick question about model conclusions. Below is syntax of a model fit and I am curious about p-value next to frailty. It is higher than 0.05. Does it mean there is no significant frailty in the data? Also `Variance of random effect= 0.493`. Does it stand for frailty or opposite that the frailty model is inappropriate? Thanks for any help!               Call:     coxph(formula = Surv(time, inf) ~ age + sex + frailty(patient,          dist = "gaus") + as.factor(type), data = kidneyExt)                                    coef     se(coef) se2    Chisq DF   p           age                        0.00489 0.015    0.0106  0.11  1.0 0.74000     sex                       -1.69728 0.461    0.3617 13.55  1.0 0.00023     frailty(patient, dist = "                          17.89 12.1 0.12000     as.factor(type)1           0.21308 0.516    0.3828  0.17  1.0 0.68000     as.factor(type)2          -1.31617 0.781    0.5817  2.84  1.0 0.09200     as.factor(type)3          -0.17986 0.545    0.3927  0.11  1.0 0.74000          Iterations: 7 outer, 42 Newton-Raphson          Variance of random effect= 0.493      Degrees of freedom for terms=  0.5  0.6 12.1  1.7      Likelihood ratio test=47.5  on 14.9 df, p=2.82e-05  n= 76


title: Defining the Stochastic and Deterministic variables with pymc3

text: I am trying to use write my own `stochastic` and `deterministic` variables with `pymc3`, but old published recipe for `pymc2.3` explained how we can parametrize our variables no longer works. For instance I tried to use this direct approach and it failed:               def x_logp(value, x_l, x_h):         if ((value>x_h) or (value<x_l)):             return -np.inf         else:             return -np.log(x_h-x_l+1)     def x_rand(x_l,x_h):         return np.round((x_h-x_l)*np.random.random_sample())+x_l          Xpos=pm.stochastic(logp=x_logp,                        doc="X position of halo center ",                        observed=False,                         trace=True,                        name='Xpos',                        random=x_rand,                        value=25.32,                        parents={'x_l':0,'x_h'=500},                        dtype=float64,                        plot=True,                        verbose=0)      I got the following error message:               ERROR: AttributeError: 'module' object has no attribute 'Stochastic' [unknown]     Traceback (most recent call last):       File "<stdin>", line 1, in <module>     AttributeError: 'module' object has no attribute 'Stochastic'      I am wondering how could I define my own variables with `pymc3` without using the available `pymc` distributions?


title: Running an AR(4) model based on engle-granger cointegration results

text: I have the below results based on the engle-granger test for cointegration, where the dependent and independent variables themselves are not stationary. The results show that null (process is not cointegrated) is rejected for lags 0,1,2 for both single mean and trend. However, lags 3,4,8 fail to reject the null. Is it then appropriate to run an AR(1,4) model, if we cannot reject the null that the process is cointegrated at lag 4? Engle-Granger Cointegration Test ![enter image description here](http://i.stack.imgur.com/hqvX2.png)


title: Please help me plug values into an MA(q) model

text: yt = u + 0et + 0et-1 + 0et-2 01= .5 02 = .25 My series of existing data is : 1,0,2,0,1,0 So I estimate the first forecast value will be : yt+ 1 = .66 + .5(-.66) + .25 * .34 yt = 1 = .4125 I have had a go at this in an excel document but I am certain it is not correct. At this stage I am not estimating the first error term just trying to understand how this works by plugging numbers into an existing formula. I don't understand how 0et would possibly work how can you know the error of something you're trying to forecast ? I will work on the rest of it as I go along. Thanks


title: Which algorithm given this plot?

text: I computed a feature `x` that I use to predict `y` which is a **probability** of being a certain class. **Raw data in R format for (x,y) is pasted here:** http://tny.cz/a97b3fd0 (500 samples) plot(x,y) looks like this: ![enter image description here](http://i.stack.imgur.com/sKoJo.png) It seems similar to a previous question I asked, so I am tempted to try this advice, which I understand is the "Least Squared Fit"on log(y) (because y is a probability). Am I correct? (My objective is to minimize the MSE.)


title: Calculating discrete hazard rates problem

text: I am working on an assignment for a Stochastic Modeling class and am stuck on the following question: > Let $X$ have probability mass function $p_j = P \lbrace X = j \rbrace $ for > $j \geq 1$. Let > > $\lambda_n = P \lbrace X = n|X > n − 1 \rbrace = > \frac{p_n}{1-\sum^{n-1}_{j=1} p_j}$. > > The quantities $\lambda_n$, $n \geq 1$ are called the _discrete hazard > rates_. The following is called the discrete hazard rate method of > generating samples of $X$: > >   * Step 1: Set $X = 1$ >   * Step 2: Generate a uniform random number $U$ on $[0, 1)$ >   * Step 3: If $U < \lambda_X$ stop and output $X$. >   * Step 4: Else, set $X = X + 1$ and go to step 2. > > > (a) Show that the output of the above algorithm has the desired probability > mass function. > > (b) Assume X is a geometric random variable with parameter p. Determine λ n > and explain how the above algorithm operates in this case. I am trying to follow the method for generating samples of $X$ but get stuck when trying to evaluate the _discrete hazard rates_ $\lambda_n$. For $X = 1$, I get $\lambda_1 = P \lbrace X = 1|X > 0 \rbrace = \frac{p_1}{1-\sum^{0}_{j=1} p_j}$ Given that the summation has poorly defined limits, I increment to $X = X + 1 = 2$, which gives: $\lambda_2 = P \lbrace X = 2|X > 1 \rbrace = \frac{p_2}{1-\sum^{1}_{j=1} p_j}$ now I am not sure how to calculate $p_2$. The assigment gives the pmf as $p_j = P \lbrace X = j \rbrace $ for $j \geq 1$, which means $p_2 = P \lbrace X = 2 \rbrace $, but I feel like there is information missing; am I supposed to assume a distribution for $P \lbrace X = 2 \rbrace $?


title: Optimising test re-test reliability through randomly generating hypothetical repeated trials

text: I have had a look around the forum, and despite some similar questions on measurement error, no one appears to have asked this question specifically. I have developed a test where a small group of participants $(n = 40)$ have undertaken three trials each of the test in total. See the two attached figures for an example Bland Altman plot representation of score comparisons for both (a) Trial 1 (on the x axis) vs Trial 2 (change in performance on y axis) and (b) Trial 2 vs Trial 3, respectively. Whilst in the second figure, the Trial 3 scores are much more similar to the Trial 2 scores as opposed to those in the first figure, on the whole a learning effect still seems to exist, with the majority of participants producing higher scores in the subsequent testing trial. ![enter image description here](http://i.stack.imgur.com/L6Fuu.jpg) ![enter image description here](http://i.stack.imgur.com/ou57q.jpg) So, as these figures show (along with the mean difference and limit of agreement values—not shown here), there appears to be on the whole a learning effect for undertaking repeated trials of the test, that may not ‘stabilise’ even after having undertaken three trials in total (not totally unexpected, however, given the complex requirements of the test).   Unfortunately, collecting further data (i.e. trials 4, 5, 6 etc.) is not feasible in this instance. Therefore, my question relates to whether there is an approach that could be used to determine how many trials would be needed in order to reach ‘acceptable’:   * mean change (improvements or decrement) when compared to the preceding trial and   * reliability indices (i.e. 95% limits of agreement or ICC’s).  I understand the subjectivity in determining what ‘acceptable’ values means in this context; however, my question is around the methodology. My thoughts at this stage—could Monte Carlo simulations be used to generate random scores for a hypothetical Trial 4 (and perhaps trial 5) based on scores reported by the initial cohort for trials 1–3 (or the trends from these trials)? This would first allow identification of when the familiarisation/learning effect from undertaking repeated trials is reduced or plateaus/stabilises (i.e. trial 4, 5 or otherwise). Second, it may also allow identification of when and if ‘acceptable’ reliability can be achieved (i.e. ICC’s of $>$ 0.80, narrower 95% limits of agreement) based on the results discovered thus far.


title: Stata question: Combining rows

text: I have something that should be a super simple one-line command in Stata. I have a bunch of rows that are disaggregated simply because one variable is different. That one variable isn't important, so I want to combine all of the rows into one and just have NA or whatever for that one variable. One tricky part -- this is only for a subsection of the data. For other rows, that variable is important. If that wasn't descriptive enough, here's a simple example (this isn't actually my dataset... ) NAME --- BREED --- WEIGHT weenie --- weiner --- 50 lbs libby --- mutt --- 40 lbs spot --- terrier 1 --- 10 lbs spot --- terrier 2 --- 10 lbs spot -- terrier 3 --- 10 lbs I want to make all those different terrier type collapse into one row with just 'terrier' or whatever as a label for that var.


title: Residual plot with a slope of 1

text: Had a question on a quiz I chose a (I and III) but the teacher said d.(all the above) > Which of the following statements about residuals are true? > > I. The means of the residuals is always zero > > II. The regression line for a residual plot is a horizontal line. > > III. A definite pattern in the residual plot is an indication that a > nonlinear model will show a better fit to the data than the LSRL I don't think II is correct. If you have three residuals (-1,0,1) then you would have a line threw all three points. Could someone help me understand what my mistake is? The regression line minimizes the distance between the points and the line, thus with my example (-1,0,1) that fits on a straight line with r =1.


title: Why does nobody use the Bayesian multinomial Naive Bayes classifier?

text: So in (unsupervised) text modeling, Latent Dirichlet Allocation (LDA) is a Bayesian version of Probabilistic Latent Semantic Analysis (PLSA). Essentially, LDA = PLSA + Dirichlet prior over its parameters. My understanding is that LDA is now the reference algorithm and is implemented in various packages, while PLSA should not be used anymore. But in (supervised) text categorization, we could do exactly the same thing for the multinomial Naive Bayes classifier and put a Dirichlet prior over the parameters. But I don't think I have ever seen anyone do that, and the "point estimate" version of multinomial Naive Bayes seems to be the version implemented in most packages. Is there any reason for that?


title: Interpretation of the regression coefficient of a proportion type independent variable

text: If I want to use some proportion type independent variables in a logistic regression, then what will be the interpretation of the regression coefficients corresponding to those proportion type variables? Will that mean- "The change in log odds for per unit change in the proportions"? But what will be meant by "per unit change" in this case? As the proportions lie within [0,1], I am getting a little confused with what a "per unit change" will mean in this scale. Does it mean 0.01 or 1%? (I am sorry for my noob thoughts!) In my data the range of the proportions is 0 to 1, not multiplied by 100. Do I need to multiply them by 100? So that I can say "per unit change" means 1% change? I have seen that the coefficients do differ in scale if I multiply the proportions by 100. For example, a coefficient of -1.3 for proportions becomes -0.013 for percentages (when the proportions are multiplied by 100).


title: What principled arguments are made with data?

text: Consider the frequentist argument: "I will believe A since not-A makes the data improbable". Consider the Bayesian argument: "I will believe A since given my beliefs on A and the data generating process, it is more probable than not-A". Are there any other argument types and how would you construct an argument? (Fiducial? Fuzzy logic? Dempster-Shafer? Possibility theory? Bounds? ...)


title: Propensity score weighting in Cox PH analysis and covariate selection

text: Regarding propensity score weighting (IPTW) when doing Cox proportional hazard modeling of time-to-event survival data: I have prospective registry data where we're interested in looking at treatment effect of a medication that in most cases the patients were already taking at baseline. I'm therefore not sure how to best analyze the data. Potentially, some of the baseline variables are to a substantial degree influenced _by_ the treatment and not the other way around (e.g. certain biomarkers). I'm a bit lost as to which covariates I should include in the propensity score model for estimating weights and which covariates I should include as covariates in the `coxph` model (if any at all). Any hints in the right direction would be helpful! I haven't been able to find any literature regarding this in CoxPh modeling as of yet. I'm thinking that covariates that represent treatments instituted at baseline that (might) influence the outcome should be included as Cox PH covariates, but I'm not sure of this. **How do I determine which variables should be included as covariates in the Cox model instead of being used in calculating the propensity score weights?** * * * _Follow-up question:_ I understand the inherit problem of evaluating a treatment effect of a certain intervention that has already begun - i.e. is prevalent among the patients, prior to start of observation. Both in regards to introducing bias related to time-variation of the risk (e.g. adverse side effects more common the first year of therapy) and covariates being affected by treatment. If I'm not mistaken - this has been proposed as a cause of discrepancy between observational and randomized that in regards to cardiovascular endpoints and hormon replacement therapy. In my dataset on the other hand, we're interested in looking at an possible adverse effect of the treatment. **If I use propensity score adjustment to investigate the treatment effect among prevalent users, i.e. already using the medication before observation starts, in cohort data and we observe an adverse effect of a pharmacological therapy (and this is what we were looking for). Can I rule out the possibilty of overestimating the risk associated with the treatment? I.e. as long as the risk is significantly elevated, it's most "definitely" not protective?** I can't quite picture an example where this kind of bias can introduce an overestimation of risk of falsy risk association in this context.


title: Multivariate meta-analysis in R: how to investigate network of variables

text: I would like to conduct a meta-analysis to investigate the interaction of three variables: **hair color** (dark/light), **gender** (male/female) and **size** (continuous). ![enter image description here](http://i.stack.imgur.com/FJVU6.png) I have three studies reporting effect sizes for the difference in size between light and dark hair:               my_data_1 <- data.frame(studies = c("study_1","study_2","study_3"),     group_1 = c(rep("light_hair",3)),     group_2 = c(rep("dark_hair",3)),     outcome = c(rep("size", 3)),     smd = c(1.2,1.3,0.8),     var_smd = c(0.2,0.23,0.15))          > my_data_1            studies    group_1   group_2 outcome smd var_smd     1 study_1 light_hair dark_hair    size 1.2    0.20     2 study_2 light_hair dark_hair    size 1.3    0.23     3 study_3 light_hair dark_hair    size 0.8    0.15      Three studies reporting effect size for the difference in size between male and female:               my_data_2 <- data.frame(studies = c("study_4","study_5","study_6"),     group_1 = c(rep("male", 3)),     group_2 = c(rep("female", 3)),     outcome = c(rep("size", 3)),     smd = c(1.8,0.9,1.55),     var_smd = c(0.25,0.35,0.1))          > my_data_2             studies group_1 group_2 outcome  smd var_smd     1 study_4    male  female    size 1.80    0.25     2 study_5    male  female    size 0.90    0.35     3 study_6    male  female    size 1.55    0.10      Three studies reporting effect sizes for the difference in gender between light and dark hair:               my_data_3 <- data.frame(studies = c("study_7","study_8","study_9"),     light_male = c(10,10,20),     light_female = c(5,7,12),     dark_male = c(15,7,14),     dark_female = c(19,21,25))     my_data_3 <- escalc(ai=light_male, bi=light_female, ci=dark_male,     di=dark_female, data=my_data_3, measure="OR2D")          > my_data_3             studies light_male light_female dark_male dark_female     yi     1 study_7         10            5        15          19 0.5125     2 study_8         10            7         7          21 0.8023     3 study_9         20           12        14          25 0.6013           vi     1 0.1275     2 0.1317     3 0.0744      All studies represent independent samples! I would like to set up one multivariate meta-analytic model to investigate the mutual influence of the three factors hair color, gender and size. Could somebody point me to potential statistical approaches to this as well as potential implementations in R (metafor, gemtc, metaSEM)?


title: Eucledian Distance and Mahalanobis distance

text: If I have a compositional data set and I wish to compare them using a metric, which is better Eucledian Distance or Mahalanobis Distance ? Comparing them after using a centered log ratio transformation. Example of Data set:   * [ 0.016, 71.2, 0.123, 1.74, 14.0, 0.002, 2310, 0.064, 0.29, 0.32, 5.63, 96.5, 0.044, 4360, 1110, 585, 0.052, 62.9, 4.45, 1110, 1.50, 15.10, 783, 0.015,78.9, 5.61, 0.007, 0.022, 0.007, 0.53, 29.3 ]   * [ 0.073, 245.0, 0.299, 2.77, 17.4, 0.039, 2460, 0.145, 0.85, 0.99, 20.40, 359.0 0.062, 4040, 1530, 148, 0.113, 217.0, 18.10, 1310, 4.61, 4.56, 880, 0.069, 230.0, 12.20, 0.028, 0.025, 0.013, 9.92, 34.1]


title: Calculating the probability of a sum of Poisson-distributed random variables

text: In Java, I have a set of users, each of which has a Poisson-distributed demand with a known mean:               import org.apache.commons.math3.distribution.PoissonDistribution;          public class User {       private int mean;       private PoissonDistribution dist;       public void setMean(int mean) {           this.mean = mean;           this.poissonDistribution = new PoissonDistribution(mean);       }     }          ...     User u1 = new User(); u1.setMean(20);     User u2 = new User(); u2.setMean(30);     User u3 = new User(); u3.setMean(40);      Now I'd like to calculate the probability that all `User`s have a cumulated demand below a certain value:               double probabilityBelowX = calculateCumulatedProbability(50, u1, u2, u3); // <- what must this method look like?      I am stuck at the question, how to solve this problem in Java. Am I missing something in the math package? I know that the demands of the customers are independent, so according to my knowledge about the Poisson distribution I can simply add up the values. But I only have `cumulativeProbability(int x)` for every single one of the user demands, but not for several at once.


title: R - fit new data to Ets models

text: I need to fit new data to a previous Ets model. From the documentation (?ets into RStudio) i can read this parameter: **model:** It is also possible for the model to be equal to the output from a previous call to ets. In this case, the same model is fitted to y without re- estimating any parameters. Lets try an example:               tmp_x_1 [1] 22 50 37 19 25 31 23 12 10 34 23 19 22 20 11 21 20 6 10 16 9 [22] 23 20 34 16 24 31 62 27 11 36 21 30 36 25 18 86 17 39 24 14 54 [43] 18 36 28 16 44 17 11 26 28 44 19 8 51 28 14 19 63 52 41 53 38 [64] 28 47 88 19 74 56 38 66 36 68 93 62 27 32 37 41 37 40 71 20 17 [85] 30 18 17 87 32 40 12 47 72 11 34 18 25 19 22 39 24 32 32 15 12 [106] 27 19 11 37 17 12 17 37 35 19 11 31 20 10 13 30 11 11 20 15 7 [127] 31 13 24 6 19 9 3 8 6 1 3 26 8 7 6 12 23 9 9 18 5 [148] 1 3 9 9 5 22 22 11 28 14 10 24 15 24 33 14 37 14 84 145 91 [169] 52 21 61 10 10 32 29 20 28 43 24 22 28 40 20 14 13 41 24 15 31 [190] 5 22 33 37 2 34 31 32 2 22 30 54 18 28 17 25 25 37 23 15 11 [211] 21 28 20 39 24 33 27 30 31 9 27 31 12 25 10 29 1 6 6 37 1 [232] 27 52 48 33 27 55 68 53 69 33 101 39 34 55 41 22 41 45          ETS <- ets(ts(tmp_x_1, frequency=1),ic="aic")          summary(ETS)          ETS(M,A,N)           Call:     ets(y = ts(tmp_x_1, frequency = 1), ic = "aic")           Smoothing parameters:      alpha = 0.2432       beta  = 1e-04           Initial states:      l = 37.9625       b = 0.7106           sigma:  0.5908           AIC     AICc      BIC       2787.545 2787.709 2801.615           Training set error measures:                 ME     RMSE    MAE       MPE     MAPE      MASE       ACF1     Training set -2.654641 18.22777 13.399 -83.84249 104.7498 0.8067377 0.03035755ETS          ets(ts(tmp_x_1, frequency=1),ic="aic",model=ETS)          summary(ETS)          ETS(M,A,N)           Call:      ets(y = ts(tmp_x_1, frequency = 1), model = ETS, ic = "aic")           Smoothing parameters:      alpha = 0.2432       beta  = 1e-04           Initial states:      l = 24.4151       b = 0.7094           sigma:  0.592           AIC     AICc      BIC       2781.344 2781.393 2788.379           Training set error measures:                 ME     RMSE      MAE       MPE    MAPE      MASE       ACF1     Training set -2.447278 18.19555 13.31933 -83.05161 104.252 0.8019408 0.0316133      I don't understand why if the vector is the same of the previous execution the output is different. How can i get the same result?


title: Using KNN for 2 Dimentional data

text: I am dealing with sentiment results of web articles. Sentiment is represented with two int values, +ve and -ve. For given article "xyz" I am getting different sentiments while testing at different time period. Each result has +Ve and -ve value. I represent them using 2-D Graph on x-y axis. My purpose is to get sentiment result with higher density. Is KNN also is suitable for this kind of scenario? Or any othe approach if someone can suggest ..


title: Importance of removing correlated predictors in non linear regression and classification

text: Removing linear combinations and correlated predictors from a dataset is something I do when I have to fit a linear regression. I wonder though how this technique affects in general the outcome of e.g. SVM or Random Forest.


title: Binary Logistic Regression Multicollinearity Tests

text: I like Peter Flom's answer to an earlier question about multicollinearity in logistic regression, but David Garson's Logistic Binomial Regression states that there is no valid test for multicollinearity for binary-dependent logistic regression, even if the independent variables are ratio scale. Can anyone supply one or more references? My own experience is that OLS correlation matrices and VIF worked for me, as my logistic coefficients went haywire before removing entangled independent variables based on the OLS tests for multicollinearity. But I have to publish my results and methods, and would like a reputable way to cite the practice, if one or more exist.


title: Coefficients from a Dynamic Panel Data Model of Economic Growth

text: I am having difficulties with the interpretation of the regression results from estimating a growth regression in a dynamic panel data set-up, estimated using Stata. (I'm using difference GMM and system GMM estimators by Arellano- Bond and Arellano-Bover/Blundell-Bond respectively). For my thesis I am interested in the effect of a set of regressors $X$ on growth. My specification is a follows: $$\ln Y_{i,t} = a \ln Y_{i,t-1} + bX_{i,t} + e_{i,t}$$   * $Y$ measures income per capita   * $X$ is the independent variable of interest, which represents a percentage where 1% is written as 0.01 for example. My difficulty lies in the interpretation of coefficient $b$, therefore I would like to ask for your help.   * Is it correct that b represents the effect on the growth rate(as the lagged income per capita is also incorporated in the model)?   * Also if the outocome is for example 0.5 how do i calculate the actual effect?   * How do i calculate the rate of convergence in this model?


title: A case of survivorship bias?

text: With the recent FIFA world cup, I decided to have some fun and determine which months produced world cup football players. Turned out, most footballers in the 2010 world cup were born in the first half of the year. Someone pointed out, that children born in the first half of the year had a physical advantage over others and hence "survivorship bias" was involved in the equation. Is this an accurate observation? Can someone please explain why he says that? Also, when trying to understand the concept, I found most examples revolved around the financial sector. Are they any other everyday life examples explaining it? Thanks!


title: Simplex Random Walk

text: This link describes how to perform a random walk on the simplex using the Metropolis-Hastings algorithm: http://en.wikipedia.org/wiki/User:Skinnerd/Simplex_Point_Picking The description says: "The stationary distribution of its samples is the unit-exponential distribution." As well as "This procedure effectively samples x_new from a gamma random variable with mean of x_old and standard deviation of h*x_old." Finally the actual update step is: "x_new <- x_old * exp( Random_Normal(0,h) )" I'm trying to understand the algorithm, specifically how it calculates the hastings ratio given the general description of the algorithm here http://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm. My question is what is the proposal distribution g(x_old -> x_new)? Is the proposal distribution actually g(x_old -> x_new) = x_old * exp( Random_Normal(0,h) ), which would make the hastings ratio g(x_new -> x_old) / g(x_old -> x_new) = x_new / x_old? Thanks for any help.


title: Use of linear regression or logistic regression when testing conditioned group differences

text: ### Background I have a colleague interested in a particular disease and specifically if the continuous variable $X$ is different between controls and patients. Preliminary results suggest that patients have higher values of $X$ than controls. The straightforward approach would be a unpaired **t-test** to test if there is a difference in mean values between patients and controls. However, the literature suggests that **age** and **sex** are correlated with $X$ so my colleague wanted to control for this when testing. Hence, she decided to use regression to solve the problem. ### Problem My colleague's supervisor reasons that since they do not know if it is the disease that causes raised $X$ or if raised $X$ causes the disease they could use either linear regression or logistic regression to solve the problem. The supervisor also argued that logistic regression was the preferable approach. If $d$ is a categorical variable with levels {control, disease} then the first model could be written $$X = \beta_0 + \beta_1 d + \beta_2 \text{age} + \beta_2\text{sex} + \epsilon$$ Where the interpretation was that there is a correlation between disease and $X$ if $\beta_1$ was found significant. and the second model $$\log(\text{Odds}(d=\text{disease})) = \beta_0 + \beta_1 X + \beta_2\text{age} + \beta_2\text{sex} + \epsilon$$ where the interpretation was that there is a correlation between disease and $X$ if $\beta_1$ was found significant. What is the opinion among the experts on cross validated? Are there other methods? ### comments The patient and control groups were not perfectly matched and my colleague wants to make sure that she has controlled for both sex and age to avoid upsetting the reviewers. I do not know if there were any significant differences between the sex and age distributions between patients and controls. There were also indications that the variance in $X$ was different in the two groups and questions were asked if this would influence the regression models. My personal opinion is that because my colleague is interested in the conditional expected value $E(X|d,\text{age},\text{sex})$ she should use the ordinary regression model or just a t-test (if there are no age or sex differences between groups).


title: Conditional Logit for recommender systems?

text: Are conditional multinomial logits used for recommendation engines? Although they are commonly used in econometrics, I've never heard it used or discussed in the context of recommender systems. Economists use multinomial conditional logits to model which of several options a person would choose and how much they value each characteristic of the items being chosen. This is often referred to as the hedonic model. The classic multinomial logit deals with discrete items (predict whether a commuter would "walk," "take the bus," or "take the subway." The conditional multinomial logit uses data on ANY sets of observed choices and does not require that each person choose among the same set of things. It also puts values on various characteristics/variables. For example, you may see people decide which of several houses to buy. Each house is different -- square footage, number of rooms, price, etc. Based on observed choices, the model estimates the importance of various characteristics and you can derive a predicted "utility" score that each person has for each house. The model then predicts that the house with the highest score is chosen. Here is a description: http://data.princeton.edu/wws509/notes/c6s3.html


title: Model averaging effect sizes of Gamma family GLMs

text: I'm trying to get some model averaged effect sizes from a set of candidate models, all of them assuming a Gamma error distribution, according to the theory given by the book from Burnham and Anderson (2002). For the calculation I'm using the function `modavg` from the R package `AICcmodavg` by Marc J. Mazerolle. In the help pages of the `modavg` function it says: > if gamma GLM is used, the dispersion parameter should be specified here to > apply the same value to each model. This seems somewhat counterintuitive, as applying the same dispersion to all candidate models would result in a "non-MLE" situation. I did not find any literature addressing this issue. Can you confirm Mazerolles restriction to one equal dispersion estimate and/or suggest some literature on that? * * * Burnham, K. P., Anderson, D. R., 2002. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach, 2nd Edition. Springer- Verlag. * * * Edit: Or, to ask in a different way: When averaging parameter estimates and standard errors - Is there a reason that for GLMs with Gaussian family you should take the different estimated dispersions, whereas for Gamma family you should take only one single dispersion?


title: Understanding tail dependence coefficients

text: How can I analyze the $\lambda_U$ and $\lambda_L$ results (estimated by non- parametric method)? What does higher or lower coefficients mean? Does $\lambda_U = 0.5$ mean there's some kind of linear dependence between $X$ and $Y$? This is what I mean by $\lambda_U$ and $\lambda_L$: ![enter image description here](http://i.stack.imgur.com/i63q0.jpg)


title: A problem on Test of Hypothesis

text: This is same as http://math.stackexchange.com/questions/647713/a-problem-on- test-of-hypothesis Let $X$ & $Y$ be two distributions take values from 1 to 100. I have a machine $M$ which follows either $X$ or $Y$ and generates integers from 1 to 100. I want to know whether $M$ follows $X$ or $Y$ by Maximum likelihood ratio test. What is the approach? How many numbers do I need from $M$ to make the success probability more than 80%?


title: How would I assign the optimal weight matrix to combine the following moment conditions in MM?

text: The following information is given: Based on a sample x(n), with n = 1,..., 100 from the exponential distribution, we want to estimate lambda but the original data was lost; we only know that 38 observations are >= 1 and 11 observations are >= 2. Combining the information based on the 38 greater than or equal to 1 and the 11 greater than or equal to 2, what is the optimal weight matrix to combine the two moment conditions?


title: Non-homogeneity in 2-way ANOVA in SPSS

text: I am trying to determine differences in fish weights from 1st to 2nd sampling between a control group and 4 different treatment groups. My groups are unbalanced (control = 1 group, t1= 2 groups, with t2,t3,t4 = 3 groups). Originally I was told this shouldn't be a problem. I've run a 2-way ANOVA and the Levene's test in SPSS is significant (0.017) but ANOVA is (0.018) and the significant factor is (.001). Is there a way to adjust for this or do I need to run a different stats test? Otherwise I was thinking I may need to simplify and combine all the treatment groups and run a t-test, then run another t-test for the control to demonstrate there are differences between treated and non-treated fish. What do you think? Ok so there are a total of 21 tanks of fish i.e. n= 21. We didnt have room for 22 tanks so one treatment group only had to 2 tanks. We attempted to run everything in triplicate so 3 of the treatment groups are triplicated. 1 treatment group is duplicated, and 1 group is the control (single). The weights of fish are tanks averages. The largest difference is a total of 22 grams, the smallest is a difference of less than 1 between tanks


title: Sample weighting with continuous variables

text: Given two sets of data:   * Total usage of a few hundred products across the population as a whole   * Per-consumer usage of those products for a sample of consumers I'd like to choose weights for each sample member such that the sum product of their weights and product consumption matches the national consumption, and such that the variance among the weights is minimized. For example, my sample might include three consumers that used products at these rates:   * Consumer A: 1.0 pound of butter; 2.0 pounds of flour   * Consumer B: 3.0 pounds of butter; 3.0 pounds of flour   * Consumer C: 3.0 pounds of butter; 1.0 pounds of floor And I might know that in my population as a whole there were 40 pounds of butter and 50 pounds of flour consumed. This can be thought of as an indeterminate system of linear equations, so there are infinitely many sets of weights that could be chosen (e.g, A=10, B=10, C=0; or A=20, B=0, C=10), but I'm specifically interested in sets of weights with low variance. If the variables were binary I might use RIM weighting or Iterative Proportional Fitting to calculate the weights, but it's not obvious how to extend those methods to data sets with continuous variables. Is there a weight selection method that works well with continuous variables?


title: interpret the PLS regression

text: I'm giving you a bit of background before asking my question. I've done a univariate PLS regression where I came out with many models. My boss asked to interpret the PLS regression for the best model. As I've never done a PLS before, I don't know how to interpret that. What does he mean by interpretation? Thanks for any leads you will be able to give me. Cheers.


title: Text mining discrete paragraphs

text: I have a spreadsheet with Titles in the first column, descriptions in the second column and whether or not they were included in our screening in the third column. We are attempting to find out which individual words were related to inclusion versus exclusion so that we can use them in a much larger data set. From my brief understanding of text mining, I think this would require a Naive Bayes or Logistic Regression method. Is this correct or is there a more appropriate method? * * * Preferably, without knowing the words in advance, I would hope to identify which are related to inclusion vs. exclusion. An example: Title: **Excel** 2010 Advanced Description: This free **Excel** 2010 eBook should be used as a point of reference after following attendance of the advanced level **Excel** 2010 training course. It covers all the topics taught and aims to act as a support aid for any tasks carried out by the user after the course. Title: **Excel** 2003 Formulas Description: Everything you need to know about Mastering operators, error values, naming techniques, and absolute versus relative references; Debugging formulas and using the auditing tools; Importing and exporting XML files and mapping the data to specific cells; Using **Excel** 2003's rights management feature; Working magic with array formulas; Developing custom formulas to produce the results you need So if these two were considered to be included then for example the word "Excel" would be a common predictor. I am sure that my question might be naive but am trying to understand the scope of the project and whether or not I can handle it, and how much assistance I will require to accomplish my goals. Thanks in advance for all your help and insight. AMAS


title: How can I find the number of support vectors in an SVM, depending on kernel?

text: Here on slide 3 I see the claim that for linearly separable cases the number of support vectors will be d+1, where d is the dimension we are working on. What if my data is not linearly separable and I use polynomial or RBF kernel? Is there a way to know what will be the number of support vectors?


title: Inventory control - first order quantity of periodic order quantities model

text: I Have read " _Manufacturing Planning and Control for Supply Chain Management 6th Edition_ " by F.Robert Jacobs. I found the problem on the **Periodic Order Quantities (POQ)**. As I know **POQ** calculation is after calculate the **Economic Order Quantity (EOQ)**. In the book contained an example problem : ![Data](http://i.stack.imgur.com/qPSBf.jpg) To calculate EOQ from the above data using the following formula : \begin{align} EOQ= \sqrt{\frac{2C_P\bar{D}}{ C_H }} \end{align} And the results of **EOQ** is `166` For **POQ** using the formula : \begin{align} POQ= \frac{EOQ}{\bar{D}} \end{align} And the result of **POQ** is `1.8` From the above data it can be said that ordering interval is 2 weeks, and this is the data from 12 period (12 weeks) : ![POQ](http://i.stack.imgur.com/wI8MM.jpg) Question : Can someone tell me, Where the amount of the first order obtained ? From 1st week plus 2nd Weeks `(10+10)` ? If I have estimation (requirement) for 13th weeks is 70, how much to order ? so for 14th weeks I don't need to order again.


title: Is Xorshift RNG good enough for Monte Carlo approaches? If not what alternatives are there?

text: I recently stumbled across an article on pseudorandom numbers in Java which mention potential weaknesses in the default algorithm, called linear congruential generator (LCG), and gives some alternatives. Among those I find Xorshift generators interesting. So I put on my researcher's hat and dove in to find more information. It turns out that it's extremely simple to implement, and very fast at runtime, which sounds great. The question is whether or not the method is a significant improvement than the default LCG algorithm. To figure that out I read through the original paper where Marsaglia introduces the method, as well as a rather tough critique/analysis by Panneton and L'ecuyer. I cannot claim that I followed the math the whole way, so I am not sure as to whether the method is a step up from LCG or not, which brings me to my first question: **_Question 1:_** Is Xorshift a step up from LCG based default RNG that java utilizes? Panneton and L'ecuyer raise the point of parameter choice, and also point out that only 3-bitshifts may not be enough. I have recently dug up a bit on hash functions (here's a relevant question I asked on StackOverflow). I wonder if one could improve the 3-bitshift method by something like this? > `long lhash = prime * (hash1 ^ hash2);` then using `(int)((lhash >> 32) ^ > lhash)` As per @whuber's comment, it's tricky business to change nuts and bolts of a RNG so this might be tricky, but I would appreciate any ideas/leads on how to amend weaknesses in Xorshift RNG. **_Question 2:_** What can be done to improve the "overall performance" of Xorshift algorithm? (note that I do not mean computational performance here) As a final note, what I want to do is to generate pseudorandom Gussians, to simulate $e^{rand}$ type numbers, where `rand` is a Gaussian with $N(0,\sigma)$. I am not too keen on adding a dependency just for this reason, and I also think using SecureRandom is pretty much overkill at this point.


title: Confidence interval does not include 0, but Wilcoxon test is sgnificant

text: I've created an error plot with CI 95% to visualize the difference (post test - pre test) between achieved scores. ![error bars](http://i.imgur.com/C28LGA0.png) n = 64 Judging from the confidence intervals one could conclude that there's is no significant improvement between the pre and post test scores (the CI does include 0). When I use a paired t-test for the pre and post scores in combination, none of them return significant as expected. But since the differences of the scores are not normally distributed for Total Score Change (p < .001) and Sub Score Change I (p < .05) according to Shapiro- Wilk, i've use the paired Wilcoxon test, which turns out significant for Sub Score Change I (p < .05). So which test should I believe and why?


title: Confusing related to kernel method of Regression

text: In the Support Vector Machine (SVM), non-linear data are mapped to a higher- dimensional feature space to be separated linearly, a kernel is used to compute the inner product in the lower dimensional space, and this is called 'kernel trick'. In Gaussian Process Regression (GPR), there are two ways to think about kernels in GP: 1) Function view, which regard the kernel as a metric to measure the similarity between two data points and 2) weight view. I am not clear on this concept. It looks this is also a kind of **'kernel trick'** , and my question is: **why does GPR need to map the data to a higher dimensional space**?


title: $\phi$-divergence?

text: I am frustrated of looking for a simple explanation of this term $\phi$-divergence, but I cannot find any. Therefore I would be really grateful if somebody could introduce a reference or write a definition (even qualitivative for it). Here is the paper that I first noticed the term in: Jiménez, R., and J. E. Yukich. "Statistical distances based on Euclidean graphs." Recent Advances in Applied Probability. Springer US, 2005. 223-239. Thanks in advance.


title: Generate a random variable with a defined correlation to an existing variable

text: For a simulation study I have to generate random variables that show a prefined (population) correlation to an existing variable $Y$. I looked into the `R` packages `copula` and `CDVine` which can produce random multivariate distributions with a given dependency structure. It is, however, not possible to fix one of the resulting variables to an existing variable. Any ideas and links to existing functions are appreciated! * * * **Conclusion:** Two valid answers came up, with different solutions:   1. An `R` script by caracal, which calculates a random variable with an _exact_ (sample) correlation to a predefined variable   2. An `R` function I found myself, which calculates a random variable with a defined _population_ correlation to a predefined variable


title: Comparing classifier performance

text: I used a SVM to perform classifications on the same sample and I get the following results: Good classifications:                       baseline : 78%     transformation 1 : 79%     transformation 2 : 63%      I do not know which test to use to tell if the performance difference is significant (ie. is the transformation 1 really better than transformation 2?). Which is the adequate test? Note: transformations are transformations performed on datas (i.e. : different imputation method for missing values).


title: Maximum likelihood and sufficient statistics

text: fT(t;B,C) = exp(-t/C)-exp(-t/B) / C-B where our mean is C+B and t>0. so far i have found my log likelihood functions and differentiated them as follows: dl/dB = sum[t*exp(t/C) / (B^2(exp(t/c)-exp(t/B)))] +n/(C-B) = 0 i have also found a similar dl/dC. I have now been asked to comment what you can find in the way of sufficient statistics for estimating these parameters and why there is no simple way of using Maximum Likelihood for estimation in the problem. I am simply unsure as to what to comment upon. Any help would be appreciated. Thanks, Rachel


title: What test to use for p-value?

text: I have a 200 samples from a population, and I want the p-value to test for significance. What test should I use? Is t-test ok? How do you determine what test to use?


title: Simple joint probability distribution table calculations

text: Given this joint probability distribution table                    |   A   |   B   |   C       X    | 0.15  | 0.03  | 0.07     notX | 0.05  | 0.25  | 0.45      I've to calc the following probabilities               P(B or X) = ?     P(C or notX) = ?      So far I got this solution:               P(B or X) = 0.15 + 0.03 + 0.07 + 0.25      P(C or notX) = 0.07 + 0.45 + 0.05 + 0.25       but I'm very unsure...


title: Ensemble learning combined with item response theory

text: I am recently learning item response theory, which is always used to select some combinations of tests for students of different specific levels. Then I have an immature idea. **Could I apply item response theory to model selection for prediction?** In this case, model prediction ability would be Y and features maybe X in the test information curve. Thus we can select a more useful combination of models for prediction according to a specific observation. And I also compared it to the traditional ensemble learning method, such as random forest, boosting. In random forest, we generate a bunch of models from the subsets, do prediction with all of them and then make a majority vote. After using my expected model, the procedure might be   1. generate a bunch of models from subsets   2. for each model, make an item information curve   3. receive an observation to be predicted, do model selection with item response theory   4. use generated combination of models to predict. I am thinking it could be a new application of item response theory and could reduce the computation thanks to model selection. But I start learning machine learning and have not thought this application in detail. Hope to get suggestions about the idea and implementation of this idea from anyone. Thanks a lot.


title: Quantify the information lost given by the Kullback-Leibler divergence measure

text: Consider there are $N$ individuals and these measure a quantity $X\in \mathbb{R}^{N\times M}$ where $M$ is the number of measurements and let $P(X)$ denote a probability distribution over $X$. The objective is to select a subset of individuals with size $n$ $(n\leq N)$ such that the probability distribution of the measurable quantity, denoted by $P(\hat{X})$ where $\hat{X}\in\mathbb{R}^{n\times M}$, becomes almost equal to $P(X)$, for any $P(X)$. To conduct the above optimization task, one criterion that I am using is to measure the distance between $P(X)$ and $P(\hat{X})$. This can be achieved by using the Kullback-Leibler divergence measure which is given by \begin{equation} D_{KL}(p\| q)=\sum^{T}_{i=1} p(x_i) \log{\left\{\frac{p(x_i)}{q(x_i)}\right\}} \end{equation} where $x_i$ is an element of the discrete space $\mathbb{L}^{N}$ and the cardinality of $\mathbb{L}^{N}$ is $T$. In this situation, the observed distribution $p$ corresponds to $P(\hat{X})$ while the theoretical distribution $q$ corresponds to $P(X)$. $D_{KL}(p\| q)$ is then a measure of the information lost (or error) when q is used to approximate p. Also, note that when $q(x_i)$ is zero then $p(x_i)$ will be zero as well, since the distribution $p(x_i)$ derives from $q(x_i)$, that is, $\hat{X} \subseteq X$. Now lets say that there is in total 100 individuals and that the subset size is 90. This means that 10 individuals have been omitted and lets assume that the Kullback-Leibler divergence measure gives for this situation something like 0.2357. My question now relates to the quantification of the error that the measure gives. In that, what does 0.2357 actually means? Could I say something like by omitting 10 individuals then 23.57% of the dataset variance has been lost? If no, is there any measure capable of estimating how much variance gets lost (or left unaccounted) when the mentioned individuals are omitted? Could the mutual information measure be utilized?


title: Correlation between approximate classification and standard classification

text: Let's say I have a large collection of emails and I want to classify in two different class: _Non-Spam_ , _Spam_. Assume this classification process is very costly, but I have access to an approximate classification procedure that outputs a value $\kappa$ in the range $[0;1]$. In general the method behaves as follows:   * If the email $i$ is spammy $\rightarrow \kappa_i < 0.5$   * If the email $i$ is non-spam $\rightarrow \kappa_i \ge 0.5$ Now what I want is some statistics assurance that the approximate metric is actually a good approximation for the original classification problem. What kind of method should I use? Initially I was thinking to take a random sample from the initial dataset, consisting of 100 emails. The point is that the entire dataset is skewed towards spammy email:               0,1,0.508     1,0,0.580     2,1,0.684     3,0,0.717     4,1,0.575     5,0,0.427     6,1,0.791     7,0,0.694     8,1,0.752     9,1,0.538     10,1,0.174     11,1,0.021     12,1,0.795     13,1,0.412     14,0,0.668     15,1,0.560     16,1,0.714     17,0,0.425     18,1,0.367     19,0,0.538     20,1,0.104     21,1,0.555     22,0,0.512     23,0,0.722     24,0,0.300     25,0,0.803     26,1,0.517     27,1,0.649     28,1,0.125     29,0,0.220     30,1,0.711     31,1,0.253     32,1,0.229     33,1,0.516     34,1,0.460     35,1,0.046     36,1,0.453     37,1,0.343     38,0,0.453     39,1,0.549     40,1,0.654     41,1,0.516     42,0,0.725     43,0,0.735     44,1,0.606     45,0,0.795     46,1,0.444     47,1,0.141     48,1,0.724     49,0,0.635     50,1,0.614     51,1,0.506     52,1,0.572     53,1,0.591     54,1,0.169     55,0,0.521     56,1,0.909     57,1,0.489     58,0,0.694     59,1,0.121     60,1,0.648     61,1,0.415     62,1,0.446     63,1,0.509     64,0,0.501     65,0,0.187     66,1,0.537     67,1,0.196     68,0,0.565     69,1,0.373     70,1,0.153     71,1,0.490     72,1,0.350     73,1,0.317     74,1,0.772     75,1,0.726     76,1,0.374     77,1,0.483     78,1,0.348     79,1,0.743     80,0,0.897     81,1,0.463     82,1,0.015     83,1,0.798     84,1,0.539     85,1,0.456     86,1,0.507     87,1,0.028     88,1,0.397     89,1,0.565     90,1,0.245     91,1,0.501     92,1,0.748     93,1,0.475     94,0,0.764     95,1,0.679     96,1,0.347     97,1,0.646     98,1,0.206     99,1,0.406      Any pointer to page/tool with a practical example will be greatly appreciated since I am relatively new to this field. # Update I have tried all the methods you have suggested but I am **still confused**. I firstly used R to cross tabulate the dataset, and then used Chi-Square test. I found out that there is a relationship between the two variables, namely the approximate function and real function. I also tried to fit a logistic regression model (and a SVM linear classifier), registering an F1 score of 0.8. But **this is somehow wrong** , since the dataset is really skewed as I previously pointed out and the system always returns the label 1 independently from the input. Therefore I tried to use an Anomaly Detection algorithm (one class SVM with RFB kernel) that seems to be more applicable here, but I don't think this is a good way to prove that the approximate method resembles the original one with high accuracy. What do you think should I use to prove my point? In order to avoid the problem with skewed label distribution, I thought that may be it's better to create a new training dataset composed as follows:   * 10 emails with $\kappa$ between $[0;0.10]$   * 10 emails with $\kappa$ between $[0.10;0.20]$   * ... Do you think this approach can improve the situation?


title: Impute values with Amelia in R from factor variable

text: I've a dataset of individuals from which I would like to impute the missing values for 'Age'. Althought the set has several columns, I noticed the most relevant one in respect to Age is the column 'Title' (ie. 'Miss', 'Mrs', 'Mr', 'Dr' and so on). My imputation is quite trivial atm:               individual$Age[which(is.na(individual$Age) & individual$Title == "Mrs")] <-      mean(individual$Age[which(individual$Title == "Mrs")], na.rm=TRUE)      and this for every Title level. Basically I impute Age for a given title using the mean of Age of the observations with that title I would like to improve my imputation a bit, perhaps using Amelia 2 for R. The question is: how can I drive amelia to make a multiple-imputation, given the apriori knowledge of a very strong regressor (Title)?


title: MANOVA with unequal sample sizes

text: I have data for 352 response ratings that have been categorized as "Low," "Medium," and "High." I would like to compare the differences of five variables across the three categories to see which ones have significant differences. In the end I would like to be able to say something like, "variable 1 differed significantly between the response categories, suggesting that it is an important component to consider..." It seemed like a MANOVA was what I wanted to use for this, but the problem is that the sample sizes between the categories are quite different: "Low" = 66, "Medium"= 275, "High" = 11. Can I still use a MANOVA for this even with such a disparity? Any suggestions?


title: Which one to use Splines / Interaction or both?

text: I am modeling a binary event of whether a sale will happen or not for an online retailer. I have millions of clicks to refer to. Needless to say the rate of sales happening is very very low (<1%). Now, to model this event, I am taking the product price to be one of the predictor variables that i am passing in the logistic model.   1. Right now, I am using the linear piece wise models or linear splines. However, I read in a paper that Cubic splines are better if the distribution is highly non linear. What should be an objective way to gauge that cubic splines should be used or not?    2. Currently I am creating knots based on the business sense and intuition. Is there a better way to know the number of knots and their placements?   3. In one of the test runs, I found that interaction between price and sellers (again one of the predictor variables!) plays an important role. Do you think that it can also explain some non linearity?  Now, I know that interaction and splines both should help me with the non linearity in my independent variable. Which one should I choose or how to decide if I should go with a combination of both of them?


title: Type II anova with poly contrasts matrix

text: I need to conduct a type II anova. Therefore, I use the `Anova` function from the `car` package. In addition, for `Factor1` I want to use a polynomial contrast matrix to see linear and quadratic effects. However, using the `Anova` function from the `car` package does not display the linear and quadratic effects compared to the normal `summary.lm(lm(...))` functions. Before I apply the `Anova` function from the `car` package I test for non significant interaction as that is required for doing a type II anova. I was wondering how I can extract the linear and the quadratic effects from the `Anova` output. Here is a working example in R: http://pastebin.com/mfk5PEq1 Cheers


title: How to measure if a mean is stable

text: ### Background I am a physics major, however I am currently interning at a psychiatry/neuro- imaging laboratory. The primary area of research in my lab is diffusion tensor imaging (DTI). A lot of the studies that are conducted here are group comparisons between "normal" controls, and patients whose brains are to some degree "abnormal" (Schizophrenia, Traumatic Brain Injury, etc.). A common procedure is to segment the brains into several regions of interest (ROI) and then calculate some characteristic numbers for each ROI, most importantly the Fractional Anisotropy (FA). To put it very simple, every region of the brain is is assigned a number between 0 and 1. Since the FA of an ROI is calculated as the mean of all voxels in that ROI, it can be assumed that the FA value for a specific ROI is normally distributed among all healthy controls (Central Limit Theorem). What happens next is that for each ROI, the mean FA is calculated from all healthy controls, i.e. for every ROI in the brain we find a number that is the "standard value" for this ROI. It is then investigated how much the patients with "abnormal" brains differ from this "standard value". ### Question/Problem An important question here that has been given to me as a summer project is how many healthy controls it takes to make this mean (i.e. the "standard value") _stable_. "Stable" means that the mean of $N$ and $N+1$ controls does "not differ much", e.g. adding another control contributes only negligible new information. I know that this is a very vague formulation, and this is also the reason why I have come here: **_What do you think would be a suitable way of characterizing if the mean of $N$ and $N+1$ controls does "not differ much"?_** My supervisor has described to me what she would ideally like to have at the end of my internship here: A program into which she can put the data of all the controls she has so far, then the program does its magic, and in the end it tells her "You have enough controls, your mean is stable", or "Your mean is not stable, and you need at least $X$ more controls to make it stable." I realize this probably sounds pretty wishy-washy and not very well-defined to you, but it's the same for me :-/ I have already started to do some research myself, and tried a couple of things, but nothing so far seemed particularly promising. Therefore, I would be very grateful for any kind of advice how to tackle this problem. Algorithms, references to literature, wild ideas, ... anything that gives me a starting point for further research would be greatly appreciated! Thanks in advance for your efforts!


title: Cross validation and prediction for unknown data

text: How do we build a model, cross validate it and use it to predict for unknown data? Say I have a known dataset of 100 points. Steps for 10 fold cross-validation are-   1. Divide the data randomly into training and test datasets in a ratio of 90:10   2. Make a model on the training dataset (90 points). (I used libSVM grid.py to optimize `C` and `gamma`)   3. Test the optimized model on the test dataset (10 points) and calculate the error.   4. Repeat steps (1,2,3) 10 times for 10-fold cross validation. Average the error from each repeat to get the average error. Now, after repeating the steps 10 times, I will have 10 different optimized models. To predict for an unknown dataset (200 points), should I use the model which gave me minimum error OR should I do step 2 once again on the full data (run grid.py on full data) and use it as model for prediction of unknowns? Also I would like to know, is the procedure same for other machine-learning methods (like ANN, Random Forest, etc.)


title: Autocorrelation from multiple time series samples

text: I have multiple samples of a time series (for example, the time series might be minutely samples from 12am to 3pm, and I have that for ten different days) and I'd like to compute the autocorrelation function $\rho(k)$ together with confidence intervals. I can think of two "obvious" things to do:   1. Chain all of the samples together end-to-end and compute the autocorrelation function and confidence intervals for the combined sample.   2. Compute the autocorrelation function for each sample individually. Average the values pointwise to get the total autocorrelation function, and apply a square root rule to get the confidence intervals. Both of these have disadvantages. Option (1) will have artifacts from where I have joined together the time series, which will become more import as I compute $\rho(k)$ for large $k$. Option (2) seems too ad-hoc - I wouldn't know whether to believe my confidence intervals. Is there a canonical or correct way to do this?


title: Recover data after PCA

text: I have a linear regression problem with about 120 predictors and I tried to remove a number of predictors from it. First I tried to remove multi- collinearities by calculating the variance inflation factor. This left me with about 20 different (hopefully not collinear anymore) predictors. Then I used a PCA to reduce dimensionality even further. Because the predictors' variances are very different to one another I used the correlation matrix for this. I can get the 'final' data when I multiply the eigenvectors with the largest eigenvalues with my original data, right? In the end I want to find out which original predictors are left and how I can recover the 'new' original data. But for some reason I am not able to recover correct numbers.


title: Growth analysis

text: I am not a statistician but a Java/R programmer. So even the topic could be wrong. Please bear with me. I am collecting details from production servers. These details could be the number of active sessions, CPU utilization etc. I draw graphs. So for example, I have a R graph showing Time(x), Web Server hits(left y-axis) and MBytes transferred(right y-axis). Another one is the number of active sessions over time. I need to understand the statistical growth pattern of this data before and after an event. These are all distributions.   1. How do I go about measuring the growth ? I use R.   2. How do I understand what causes that growth ? I think this is about regression. The more complex question from my perspective.   1. What is the statistical process to predict future growth in these cases ? I already read Capacity Planning books. This is not about Capacity planning which probably is the next step. I have come across topics like 'Growth Analysis curves and Visualization' but couldn't access any material even after searching. Mohan


title: Standardization in SPSS

text: I have a survey based on a Likert-scale (1-5) with around 20 variables and 140 cases. Since some people only cross from 1 to 3 for all questions (variables) and others cross 1-5, I want to standardize the values. So I want to calculate the average of all variables of one person to use the result for the standardization (to have a standard deviation of one and a mean of zero). Does anybody know how this standardization is called and how I can do it in SPSS? Thanks a lot!


title: Partial Dependency plots and Gradient boosting (GBM package)

text: Is it possible to plot a partial dependency plot to display the class probability and estimate the effects of a predictor for a GBM model? Something similar to `partialPlot` from `randomForest` package. ![enter image description here](http://i.stack.imgur.com/ema8x.png) According to this article, a partial plot is doable with gbm. Thanks in advance for your help.


title: Generating artificial experimental data to test model fitting software

text: I wish to test some model fitting software and I would like to generate synthetic datasets for this test. The synthetic data is supposed to originate from a experiment that measures the change in concentration of a substance over time. I can generate error free data for such an experiment using a differential equation model. My question however is what kind of noise should I add so that the data resembles a real experimental dataset? My current two choices are to either add noise from a gaussian or exponential distribution, any preferences? (Eg adding gaussian noise could potentially generate negative concentrations) My other concern is should the degree of noise be equal for all data points? For example I would imagine there would be more uncertainty in small concentrations that larger ones, simply because it is more difficult to measure small concentrations.


title: Modelling for soccer scores

text: In Dixon, Coles (1997), they have used the maximum likelihood estimation for the two modified independent Poisson models in (4.3) to model the scores in soccer. I am trying to use R in order to "reproduce" the alpha and beta as well as the home effect parameters (pg. 274, Table 4) without using any packages (using the usual independent Poisson models are fine too). I have tried using `bivpois` package but I am not sure on how to modify its parameters. I would greatly appreciate it if anyone can help me with the R code to model the data - Scores from the home and away team for Season 2012/13 in English Premier League.


title: Negative Binomial Regression: is parameter theta (R) the reciprocal of parameter kappa (SAS)?

text: After some frantic googling I do believe the answer is yes, but more so I am frustrated that the relation between the two parameters seems to be nowhere described explicitely so I do it here. (I hope this isn't against the rules of stackexchange.) This very nice article states: we will denote the random variable Y having a negative binomial distribution as Y ~ NB($\mu, \kappa$) with a parameterization such that E(Y) = $\mu$, var(Y) = $\mu + \kappa \mu^2$. I take this latter equation as the definition of $\kappa$. Apparently this kappa is implemented in SAS. Now turning to R, the function `glm.nb` in the `MASS` package contains a parameter $\mu$ which is obviously the same $\mu$ as above and a parameter $\theta$. The question is how $\theta$ and $\kappa$ are related. The documentation for `glm.nb` only refers to it as an "additional parameter". The answers to this and this stackexchange questions directly imply that $\theta = 1/\kappa$, but this question seems to suggest that $\theta = \kappa$. The help page for negative binomial in R is nice and introduces a parameter called `size` that equals $1/\kappa$. Fitting `glm.nb` on random data generated by `rnbinom` for various choices of $\mu$ and `size` seems to support the thesis that $\theta = 1/\kappa$ (i.e. that $\theta$ = `size`) but also that for large values of size the estimation is poor. Summarizing: I do believe that $\theta = 1/\kappa$ but it would be nice if there were an easily googlable place on the internet stating this explicitly. Maybe one of the answers to this questions can serve as such a place?


title: MCMC direct comparison of difference of two parameters

text: Say I have run a Hierarchical Bayesian model in STAN (or JAGS or BUGS) and I have the posterior samples of two slope parameters that I want to compare: $\beta_1$ and $\beta_2$. The model appears to have converged properly. So I look at the density of the difference of the MCMC samples of the posteriors of the two slopes ($\beta_1-\beta_2$). The 95% HPD interval of this difference includes zero, and from what I've read this means something similar to a Frequentist failure-to-reject: zero is credible and hence we cannot say that $\beta_1 \ne \beta_2$. Yet, I can see that 80% of the differences are greater than zero. Which tempts me to say that, given my assumptions and model the odds are 4-1 that $\beta_1>\beta_2$. But I don't think I've ever read anything like that. EDIT2: OK, I found an answer to my question. The previous paragraph describes a decision rule, which is distinct from the model, MCMC, etc. This particular decision rule has two problems:   1. It doesn't take into account the density of the differences. That is, the 80% of the differences that are greater than zero might all be very close to zero, or might be quite distant from zero. This should make a difference, but the simplistic decision rule doesn't differentiate between these conditions.   2. The 4-1 odds are not all that strong to begin with, considering that the conventional 95% CI is essentially 19-1 odds. I'm not sure how to make a better decision rule, though.


title: Is clustering (kmeans) appropriate for partitioning a one-dimensional array?

text: I want to group the outcome of a function into 2 (or 3) categories. I have a function efficiency=f(weight,speed,#refueling_stops) that takes 3 input parameters and the output tells me how "efficient" a truck is. My goal is to take the most inefficient trucks off the road. However, I don't know which truck to keep and which to reject. In other words, I want to divide all possible output values of my function into the category "keep" or the category "reject" (or the category "between"). Furthermore I have no way to rate how suited my decision was and therefore the point where I draw the line(s) is more or less arbitrary. Nevertheless, I'm looking for a science-based approach to this problem. Is there a name for this kind of problem? So far I've stumbled upon clustering (kmeans and natural breaks / Jenks) which is completely new to me. Also I've read that my problem may be similar to converting a color image into black and white (and gray). But I couldn't find out what the current practice for this process is. Up to now, I've calculated all possible outcomes of my function. The histogram and PDF of the resulting one-dimensional array look like this: ![histogram and pdf](http://i.imgur.com/C84knXg.png) Then I partitioned them into 2 (or 3) categories via R:               library(classInt)     x <- read.table("all_possible_outcomes")     classIntervals(b, n=2, style = "kmeans")     classIntervals(b, n=3, style = "kmeans")      Now I'm curious if this approach to my problem is a current method, or if not, what is the best practice for this? I guess what I'm looking for is some kind of confirmation that it's appropriate to use clustering. If not, what alternatives can you think of?


title: Calculating variation of one variable within one experiment

text: I am conducting a preliminary study, where people assess (7-point likert scala) different pictures. Besides asking for my 2 Independent variables (which are the manipulations for my main study, an experiment) I included a third question, the covariate variable. I need to calculate / quantify if this covariate variable changes for the several pictures or if the scores are similar. All the analysis of variance (ANOVA, MANOVA, etc) I find compare always 2 or more different variables. In my case I don't want to compare Independent and Dependent variables. I just need to see if the rankings for the same variable differ significantly for the different pictures. Are descriptive statistics (variance and SD) enough to statistically support the variation of the variable? Thanks!


title: Intuitive explanation for periodicity in Markov chains

text: Can someone explain me in a intuitive way what the periodicity of a Markov chain is? It is defined as follows: For all states $i$ in $S$ $d_i$=gcd$\{n \in \mathbb{N} | p_{ii}^{(n)} > 0\} =1$ Thank you for your effort!


title: How to select tuning parameter for regularized regressions for interpretation?

text: I'm using linear regression to predict a continuous variable using a large number (~200) of binary indicator variables. I have around 2,500 data rows. There are a couple of issues here:   * When I run univariate regressions, most of the indicators have significant associations with the dependent variable (by design - I collected data I thought would be associated with the dependent variable)   * Some of the indicators have very low incidence (~1% or less) (but I try to exclude these variables from analysis)   * Some of the indicators co-occur I'm looking for a stable, parsimonious model. I've been using AIC or BIC selection. I'm aware this can be frowned upon, but I'm comfortable with it, primarily because I'm only feeding in variables I'd expect to have some association. I don't want to use a standard Enter regression because of IDV (independent variable) co-occurrence. Unfortunately, if I construct a new dataset using random draws of cases from the original dataset, and refit my BIC model, I find that after 20 such loops the majority of the terms each individual BIC model retains are not preserved across all 20 models. I'm thinking I should use some sort of regularized regression - say ridge (for its ability to deal with IDV correlations), lasso (for its ability to set parameters to 0) or elastic net (to combine these advantages). My understanding is that this is akin to using a softer thresholding criteria than stepwise AIC or BIC that can yield more repeatable results. I've never used these techniques before, but after doing some research I'm comfortable selecting a tuning parameter using cross validation to maximise performance in terms of prediction error (though if anyone feels compelled to provide a walkthrough here, I wouldn't object). **My main question:** I'm actually more interested in _interpreting_ my model than I am in using it to make predictions. I've read that selecting a tuning parameter for interpretation is harder than it is for prediction, but haven't found any more information. Can someone point me in the correct direction? I'm using R, if anyone is curious. **A secondary question:** R provides (I guess Wald) significance estimates for the models in a lasso model. Do I trust them? Can I, at a pinch, interpret them? I have an engineering background, so I'd particularly like referenced answers (web resources are fine) if possible, and math is welcome as long as its provided with an explanation/intuition. Thanks!


title: Recommendation for sequential least squares programming book

text: I would like to understand better how to solve least squares in a sequential/recursive way. How to use weights in the least squares if the following problem:               y = A*a + B*b      where the estimation of `a` should be recursively updated, and `b` should be estimated each time separately. Any recommendation for a book that can explain it well?


title: Introduction to recurrent neural networks?

text: I have two questions: 1- What are the applications of recurrent neural networks? 2- Can you recommend some good resources/papers/tutorials that introduce recurrent neural networks?


title: Model estimation using ACF and PACF

text: Can anyone help in model estimation ? The following are the ACF,PACF and the plot of the sample respectively. ![plot of data](http://i.stack.imgur.com/1RxDS.png) ![acf of data](http://i.stack.imgur.com/ft0Ls.png) ![pacf of model](http://i.stack.imgur.com/9QqEN.png)


title: Determine penalty parameter based on AIC values in lasso regression for linear mixed model

text: I have a dataset with 14 variables, of which one random effect (Patient) I fitted following code: "fit1_B2=lmmlasso(y=y,x=x.matrix,z=z,grp=grp,lambda=0.1,pdMat="pdIdent") summary(fit1_B2) fit1_B2$aic" I have, however, chosen the penalty parameter at random. The purpose is to obtain the penalty parameter that minimizes the AIC value. Does anyone have an idea how this might be done? Many thanks in advance. Pieter


title: Simulated annealing - why when encountering a lower energy state we move to it with probability 1?

text: I'm doing simulated annealing, following a recipe I have -   1. Compute the energy of current state and a neighbor state.   2. If the energy of the neighbor is lower - move to neighbor.   3. If not - move to neighbor with probability $e^{-\Delta/\tau}$ where $\Delta$ is the difference between new energy and old energy. This seems a bit arbitrary to me. Why not compute the exponent for both the current state and the neighbor, and then just move/stay with these probabilities? This means using a sigmoid, which is derived by dividing the exponent of one of the options with the sum exponents of the two: So the probability to move is $e^{-E1/\tau} / e^{-E1/\tau} + e^{-E2/\tau}$ which can also be written $1 / e^{-\Delta/\tau}$


title: How the hypergeometric distribution sums to 1?

text: The hypergeometric distribution is defined for $\max(0, n+K-N)\leq k\leq \min(K,n).$ But, when we use Vandermonde's identity to prove that probabilities sum to $1$, then we use the range of $0\leq k \leq n.$ I wonder how this is justified?


title: Should the final R glm include only significant levels of factors

text: I am running a glm in R on data with quite many predictors (~50), both initially continuous and factors. The response is binary and the volume of the data is OK (~100K rows), in order to model non-linear relationships, I convert the continuous variables to factors as well. In the end this results in some levels of some variables being insignificant. For example, variable "age" is binned into 20-30, 40-50, 50-60, 60+ and only the first two are significant. All this results in a reasonable model with AUC 0.86 and residuals looking random with a small bias. I understand that R converts each factor level to a binary variable and then runs the regression, and I'd still like to improve the model. Would you please help me understand if:   * converting all the variables' factor levels to binary myself,   * selecting only the variables that were significant in the initial model   * then re-fitting the model with the new (reduced) variable set sound like a good idea? Or should I continue to use predict() on the model with many insignificant coefficients for factor levels happily as before?


title: cointegration with time trend and binary variables

text: I have an equation composed of one explained variable and three explanatory variables, which are all non-stationary. Is it possible to add for the equation another two variables: a time trend and a binary variables and then check for cointegration (with proper results)?


title: Linear regression in R

text: I am confused about the usage of `lm` and `predict` in R. I have a predictor variable $x$ and an output variable $y$, and I want to perform a univariate regression analysis. I have $15$ predictor samples from $x$, and $15$ output samples from $y$ which I will use for training. Then, I want to give a number as a predictor to the trained regression model and get the output value from the model. Here is an example R code:               set.seed(23)     x <- rnorm(15)     y <- x + 1     predict(lm(y ~ x), .3)      But the result of this `predict` call is a vector of size $15$, which has actually the same values as in $y$, although I was expecting it to be just a single number, which is the output of the test predictor $.3$. What's wrong with that code, and how should I change it to get what I want to get? I know this is a very basic question, but I am new to that stuff. Thanks!


title: Equivalent of Kolmogorov-Smirnov test for integer data?

text: Is there an equivalent of the two-sample Kolmogorov-Smirnov test for integer data (not count data, as it can include negative integers)? The Kolmogorov-Smirnov test does not perform well in the presence of lots of ties, which are obviously common with integers.


title: In regression analysis what's the difference between data-generation process and model?

text: In regression analysis what's the difference between 'data-generation process' and 'model'?


title: How to determine statistically whether an effect saturates?

text: I have a behavioral measure obtained from N participants over 5 levels of a condition C -- 10, 30, 50, 70, an, 90 degrees. With a one-way repeated measures ANOVA, I find a significant effect of C. Now qualitatively it appears that the behavioral measure increases from 10 to 30, and slightly from 30 to 50, but then saturates at 50 degrees. Is there a way to show statistically that I'm finding an effect that increases at first and then saturates? Would it be sufficient to say that 10 & 30, and 30 & 50 are statistically different (through paired t-tests), but 50 & 70, and 70 & 90 are not (and that the resulting p-values for these two are very large)? Thanks!


title: What are the problems with small sample size based predictions when using time-series data?

text: I had to put the majority of the question in the textbox as Cross Validated didn't like the fact my question was too long....but here it is, in terms of what I wanted to ask: Due to a small sample size (i.e. large earthquakes are rare), does predicting 'post-large earthquake aftershocks' suffer from the same statistical problems of predicting the actual earthquake itself? or Does the consistency (assuming there is some) of 'post-large earthquake aftershock' effect make it valid to make 'strong' (with high probability of success) and accurate predictions about when/where and how large these aftershocks might be? And finally are the statistical techniques different for predicting 'post-(rare) event events' vs predicting the actual (rare) event itself? Any thoughts, links or pointers would be much appreciated. P.S. This question is partially related to my previous question Can survival analysis be used to predict earthquakes?


title: Differential expression gene analysis

text: I have gene expression data in two samples like this:               Gene           s1   s2     C9orf152    2.96295 2.99861     RPS11       11.2968 11.0775     ELMO2       8.16786 7.6039     CREB3L1     6.24832 6.69416     PNMA1       8.11559 7.90189     MMP2        8.32177 7.06863     TMEM216     5.48133 5.68662     C10orf90    2.91206 3.11626     ZHX3        5.47916 5.31362     ERCC5       5.94527 5.92997     GPR98       3.09796 3.24092     RXFP3       5.17484 5.14564     CTAGE10P    2.71189 2.60001     APBB2       7.16389 7.40401      Now I want to know if any test can figure out the differential expressed genes between `S1` and `S2`. T-test can only give one p-value between two samples, right? I want to figure out each p-value for each gene between `s1` and `s2`.


title: Variance component analysis with lme4

text: I am beginner with lme4 and I am seeking some advice on how to carry a variance component analysis. The data come from a real-world scenario, not from a designed experiment, resulting in a complex nesting structure:   * _"Students"_ are examined on 4 different practical _"Topics"_.    * For each topic the performance ( _"Scores"_ ) of the students is rated by a different _"Rater"_.   * Also, since the cohort of students is very large, this examination process is further divided multiple times: two _"Circuits"_ (first and second) are running in parallel twice a day ( _"Session 1"_ and _"Session 2"_ ). Finally the exams goes on for several _"Days"_ ("1" , ... , "D").  To sum up, here are the factors with their levels:   * Days : {1, ..., D}   * Sessions : {1, 2}   * Circuit : {First, Second}   * Topics : {1, 2, 3, 4}   * Raters : {1, ... , R}   * Students : {1, ... , S} Therefore, the data set (just for Day == 1) I am working on has the following structure:               Day    Session  Circuit      Topics      Raters      Students    Scores     1          1  First              1            1           1      0.78     1          1  First              1            1           2      0.72     1          1  First              1            1           3      0.60     1          1  First              2            2           1      0.91     1          1  First              2            2           2      0.74     1          1  First              2            2           3      0.16     1          1  First              3            3           1      0.95     1          1  First              3            3           2      0.89     1          1  First              3            3           3      0.85     1          1  First              4            4           1      0.51     1          1  First              4            4           2      0.96     1          1  First              4            4           3      0.40     1          1  Second             1            5           6      0.06     1          1  Second             1            5           7      0.98     1          1  Second             1            5           8      0.65     1          1  Second             2            6           6      0.51     1          1  Second             2            6           7      0.06     1          1  Second             2            6           8      0.14     1          1  Second             3            7           6      0.85     1          1  Second             3            7           7      0.66     1          1  Second             3            7           8      0.17     1          1  Second             4            8           6      0.58     1          1  Second             4            8           7      0.17     1          1  Second             4            8           8      0.63     1          2  First              1            9           9      0.41     1          2  First              1            9           10     0.59     1          2  First              1            9           11     0.33     1          2  First              2           10           9      0.39     1          2  First              2           10           10     0.33     1          2  First              2           10           11     0.87     1          2  First              3           11           9      0.50     1          2  First              3           11           10     0.09     1          2  First              3           11           11     0.97     1          2  First              4           12           9      0.74     1          2  First              4           12           10     0.59     1          2  First              4           12           11     0.47     1          2  Second             1           13           14     0.42     1          2  Second             1           13           15     0.77     1          2  Second             1           13           16     0.05     1          2  Second             2           14           14     0.48     1          2  Second             2           14           15     0.88     1          2  Second             2           14           16     0.53     1          2  Second             3           15           14     0.95     1          2  Second             3           15           15     0.85     1          2  Second             3           15           16     0.89     1          2  Second             4           16           14     0.41     1          2  Second             4           16           15     0.17     1          2  Second             4           16           16     0.53      What I would like to estimate is how much of the variation in the "Scores" can be associated with a given factor ("Day", "Session", "Raters", etc...). For this purpose I was trying to apply a linear-mixed effects model. What I find confusing is the nesting structure of my variables. Using the coding from the library lme4 I would be inclined to write:               lmer( Scores ~ 1 +  ( 1 | Days ) +                         ( 1 | Sessions ) +                         ( 1 | Circuit ) +                         ( 1 | Topics ) +                         ( 1 | Raters ) +                         ( 1 | Students ) +                                        ( 1 | Students : Topics ) + # Students are nested into Topics, Circuit, Sessions and Days                           ( 1 | Students : Circuit ) +                         ( 1 | Students : Sessions ) +                         ( 1 | Students : Days ) +                         ( 1 | Topics : Circuit ) + # Topics are nested into Circuit, Sessions and Days                          ( 1 | Topics : Sessions ) +                         ( 1 | Topics : Days ) +                         ( 1 | Circuit : Sessions ) + # Circuit are nested into Sessions and Days                         ( 1 | Circuit : Days ) +                         ( 1 | Sessions : Days ) + # Sessions are nested into Days                         ( 1 | Topics : Circuit : Sessions : Days ) +                         ( 1 | Topics : Circuit : Sessions ) )      Is this the proper way to extract the variances components associated with each factor? Any help would be very much appreciated many thanks, again


title: Using pooled OLS when running a model with panel data?

text: How bad is it to use pooled OLS instead of fixed effects when you have 7 years of panel data? From what I have understood, the risk is that the coefficients will be correlated with the error term, thus making the estimates biased. There will be some form of endogeneity. Would it help if I include year dummies in the pooled OLS regression? It still wouldn’t capture the effects of varying intercept in the individual dimension, right? One of my major explanatory variables is significant at the 5% level in FE regression. In the pooled OLS it is significant at the 0.001 level. Is this result negligible or could it still be used with the reservation that it is overestimated? I ask this because most of the estimated parameters are strongly significant in the pooled OLS regression. Also, two of my explanatory variables that are constant get dropped in the FE regression. Although they are of secondary interest they contribute by explaining quite a lot of the variation in the dependent variable. (The sample is btw not congruent with a random effects model). Is there some way to decide which model might be more suitable? If you know some things I should keep in mind when implementing the models I would be very grateful to hear them! (I asked this question at another forum. I'll update if I get an answer. http://www.talkstats.com/showthread.php/56320-Using-pooled-OLS-when-running-a- model-with-panel-data?p=159061#post159061 )


title: Evaluating relative accuracy/error for continuous value prediction (and assessing relative average difference/error)

text: I was trying to understand what are good (or standard) way of evaluating relative accuracy for continuous data. Say for example, say I have some statistical algorithm S that outputs some real number from a vector of data. Something like $S:\mathbb{R}^k \mapsto \mathbb{R}$ Intuitively what I am interested in finding is a measure of how often S gets its approximation right. Though, since the value is not discrete, its not easy to compare it to the true value (we can't just count the times its not equal to the true value because then the accuracy would always be zero with very very very high probability). However, we could make the algorithm S output say $n$ predictions with value $p_i$ and each has its own true value $t_i$ and then, calculate the difference from $t_i$ for each true value $t_i$. The problem is that since they all not have the same true value $t_i$, that would be maybe kind of a weird way of calculating the "average relative distance form the true value". Instead maybe calculating an average relative distance made more sense to me but was not sure if such a measure was something reasonable. Something like: $$\delta = \frac{1}{n}\sum^{n}_{i=1}\frac{|p_i - t_i|}{t_i}$$ If you think this is not a good metric, it would be useful to point out its weaknesses and what other metrics are better and why. Also, as n approaches infinity, does that metric approaches 0? Please address those points in your answer (in as much detail as possible). I was interested in doing something of that sort and was wondering what the statistics community thought about it. * * * On a side note, alternatively something that also seamed reasonable could be counting the number of times S is more than 10% (some percent) away from the true value.


title: Find variables selected for each subset using caret feature selection

text: I am doing feature selection using the command 'rfe' in the caret package (http://caret.r-forge.r-project.org/featureselection.html). This command uses a metric to find the optimal amount of variables and which variables that is. However, I would like to also see the other steps in the feature selection than simply the last one. For instance, I would like to know which variables were the optimal ones if I wanted exactly 10 variables. My code is the following:               ctrl <- rfeControl(functions = rfFuncs,                        method = "cv",                        verbose = FALSE)     subsets <- c(5,10,15,20,25)     lmProfile <- rfe(dat2_X, dat2_Y,                      sizes = subsets,                      rfeControl = ctrl)


title: Need a rigorous statistical framework for automating visualization

text: I am faced with a challenging problem. Suppose I have a large dataset with many attributes and I can filter the data using a set of attributes. The problem is in the event we have a large number of filters how can we use statistics to determine which data to show as visualization to the end user. The end user cannot have the time to look at each filter. If I have unlimited computation power I can potentially evaluate each filter using cluster analysis or Kohonen's self-organizing maps, but need further statistics to eliminate the ones that are not "interesting" in an automated manner.


title: Does sphericity require the covariances between difference scores to be the same?

text: Field's Discovering Statistics Using SPSS (2013) defines sphericity as follows: > Sphericity: a less restrictive form of compound symmetry which assumes that > the variances of the differences between data taken from the same > participant (or other entity being tested) are equal. This assumption is > most commonly found in repeated-measures ANOVA but applies only where there > are more than two points of data from the same participant. This definition doesn't mention anything about the covariances between difference scores needing to be equal. Is that because the definition is incomplete, or because the covariances don't need to be equal? Why is it necessary/not necessary for covariances to be equal? I'm thinking here about a within-measures ANOVA. EDIT: I now include a table from Field (2013) to illustrate more clearly what I mean by difference scores. ![enter image description here](http://i.stack.imgur.com/6Z022.png) I have calculated the covariances and they are as follows: AB -8.05 AC 7.65 BC 2.65 It turns out that the variances of the difference scores are a bit dissimilar (15.7, 10.3, 10.7). However, let's say those variances were exactly the same, but the three covariances remained very different from one another. Would sphericity thereby be violated? Or are the covariances between difference scores not relevant to sphericity? Why are they relevant or not relevant to sphericity?


title: pseudo Anova Table for a GLS regression

text: I fit a couple of GLS models, using the `nlme` pakage (`gls()` function), and I would like to obtain the pseudo ANOVA table as described by Griffis & Stedinger (2007). `nlme` doesn't calculate the pseudo Anova table. Despite I can understand the method described by Griffs and Stedinger, I am unable to program this. I am wondering if there is any R package that allows obtain this table. Any help would be appreciated.


title: Optimization of MLE for mixture problems

text: I have about 1000 data points from some thick tailed distribution that I would like to fit a parametrized distribution to. From my data, I've made some adjustments and constructed an empirical distribution (so I have percentiles). What is the best way to fit a mixture of parametrized distribution functions (pareto, lognormal, gamma, etc) to this empirical distribution? So far I have been using Excel to maximize the grouped MLE function; using solver to maximize the parameters subject to sum(weights = 1). I have R as well but am new to it. It is pretty obvious that excel is getting stuck in local maxima. How would you maximize a MLE function for a mixture distribution?