title: Show that $\hat\theta=\frac{2 \bar Y- 1}{1- \bar Y}$ is an unbiased estimator for $\theta$

text: Let $Y_1,Y_2,...,Y_n$ denote a random sample from the probability density function $$f(y| \theta)= \begin{cases} ( \theta +1)y^{ \theta}, & 0 < y<1 , \theta> -1 \\ 0, & \mbox{elsewhere}, \end{cases}$$ Find an Estimator for $\theta$ by using the method of moments and show that it is consistent. I have found the estimator but unsure how to show that it is consistent estimator. E(Y)=$ \frac{ \theta +1}{ \theta +2}$ and $m_1'(u)= \frac{1}{n} \sum_{i=1}^{n}Y_i= \bar Y$ Now, E(Y)=$ \frac{ \theta +1}{ \theta +2}$=$m_1'(u)= \frac{1}{n} \sum_{i=1}^{n}Y_i= \bar Y$ So $\bar Y=\frac{ \theta +1}{ \theta +2} \to \hat\theta=\frac{2 \bar Y- 1}{1- \bar Y} $ Now I am unsure how to show that $\hat\theta=\frac{2 \bar Y- 1}{1- \bar Y}$ is an consistent estimator for $\theta$ Can someone please help me?


title: Is it possible to specify the indices for K-fold cross-validation with the tune function in R?

text: I would like to use the `tune()` function in R to specify the hyperparameters in an SVM model. The particular SVM implementation that I am using is the one that is provided in the `e1071` package. By default, `tune()` runs a K-fold cross-validation, and chooses the testing and training sets for each fold by itself. I am wondering whether I can override this behavior and specify the indices of the testing and training set for each fold myself. If not, does anyone know if there is a package or combination of packages which allows me to run a K-fold CV on SVM **and** specify the indices of the training and testing set.


title: Heteroskedasticity - residual plot interpretation

text: I am plotting a residual plot to test for heteroskedasticity. The Breusch- Pagan test is significant and therefore I am suspecting there is evidence on heteroskedasticity. The question is: (a) How can I interpret such graph? I know some points seem to be on top of each other etc. is it as simple as that? (b) Also the data has a lot of binary dummy variables. Could this be the reason for heteroskedasticity? ![Residual plot](http://i.stack.imgur.com/owiAH.png)


title: If any parametric test does not reject null, does its nonparametric alternative do the same?

text: If nonparametric tests are assumed to have less power than their parametric alternatives, does this imply that if any parametric test does not reject null, then its nonparametric alternative does not reject null too? How can this change if assumptions of parametric test are not met and the test is used anyway?


title: Decision boundary equation of the perceptron

text: As I know the standard linear equation has the following form in $R^2$: $w_1 x_1 + w_2 x_2 = b$ where $b$ is the intercept with $x_2$ Also I know that a decision boundary in $R^2$ for a perceptron is a line. However I read everywhere that the decision boundary equation is the following: $w_1 x_1 + w_2 x_2+ b = 0$ Why is that? Shouldn't it be $w_1 x_1 + w_2 x_2 - b = 0$!?


title: Test for markov-property in a time-series

text: Given an (observed) time-series $X_t$ with $X_t\in\{1,...,n\}$, is there a statistical test for testing the null-hypothesis that $P(X_t|X_{t-1},X_{t-2},...,X_1)=P(X_t|X_{t-1})$ (i.e. the markov-property)?


title: CLT and stable distributions

text: I have a few questions about generalizations of the CLT and stable distributions. I'm trying to correct my understanding and make it precise. Please forgive my naivete, I am not a professional statistician :-) If I take the sum of a large enough sequence of independent R.V.'s, do they always converge to a stable distribution? (I've heard about generalizations of the CLT, but I'm looking for more precision). When working with real data, what would be a hint that I need to model with a stable distribution? Is it possible to perform max likelihood with stable distributions?


title: How to evaluate dimension reduction from n-space to d-space?

text: I'm performing dimension reduction on some data sets and would like to evaluate how has a particular dimension reduction algorithm performed in terms of how much data is lost. If we are given 1000 dimensions, and we reduce it to 2, then how effective is it? I'm trying to figure till what should you do DR such that your results don't go bad? Is there a metric which does this? I'm using PCA. **Edit:** Can I use some distance metric to do the evaulation?


title: Is there a law that says if you do enough trials, rare things happen?

text: I'm trying to make a video about loaded dice, and at one point in the video we roll about 200 dice, take all the sixes, roll those again, and take all the sixes and roll those a third time. We had one die that came up 6 three times in a row, which is obviously not unusual because there should be a 1/216 chance of that happening and we had about 200 dice. So how do I explain that it's not unusual? It doesn't quite seem like the Law of Large Numbers. I want to say something like "If you do enough tests, even unlikely things are bound to happen" but my partner said people might take issue with the "bound to" terminology. Is there a standard way to state this concept?


title: Distribution with $n$th cumulant given by $\frac 1 n$?

text: Is there any information out there about the distribution whose $n$th cumulant is given by $\frac 1 n$? The cumulant-generating function is of the form $$ \kappa(t) = \int_0 ^ 1 \frac{e^{tx} - 1}{x} \ dx. $$ I've run across it as the limiting distribution of some random variables but I haven't been able to find any information on it.


title: Unpaired repeated measures data

text: I realise this is one of those "Don't start from here" type questions, and I did indeed tell them not to "start from here" but it would be a shame to junk the data. Basically I have some interval normal data, Likert scale type stuff, from two groups, at time one and time two, and I know that some subjects in time two are the same ones as from time one, so really the data should be paired. However, the data was collected anonymously so we just have two groups to compare. My best thought would be that I could report a very conservative estimate of difference and then say "for all the problems with data collection, this is the smallest the effect is likely to be". Edit: added- If significance is possible too, i.e. "and I am reasonably confident this effect is nonzero" then, great. I don't even know where to begin looking up such a stupid way of analysing data. Would bootstrap be useful at all? Any suggestions gratefully received, as I say highly conservative ones I think would be fine.


title: Can you do a paired t test in SPSS or GraphPad with only means, errors and N?

text: I have to verify some results. I have only the mean and its st. error, not the raw data. In a research paper they did a paired t test. Is it possible to check whether their p value is right by only using the means, st. errors and N, in SPSS or GraphPad?


title: Perturbations of rational numbers with random variables

text: For $a_1,a_2$ and $b_1,b_2$, $\in\Re^+$, if $a_1<b_1$ , then for any perturbation by a random variable $\epsilon\in \Re^+$, $$r_1=\frac{a_1+\epsilon}{b_1+\epsilon}>\frac{a_1}{b_1} $$ and if $a_2>b_2$, $$r_2=\frac{a_2+\epsilon}{b_2+\epsilon}<\frac{a_2}{b_2} $$ If $\epsilon$ is generated from a density as, $\epsilon_1={f_\epsilon(t=x)}$ and $\epsilon_2={f_\epsilon(t=y)}$ for $r_1$ and $r_2$ respectively or as a special case where $\epsilon=\epsilon_1=\epsilon_2$, when $x=y$, What would be a way to characterize $r_1-r_2$ in terms of $f_\epsilon(t)?$ _Note_ : The parameters of the density function $f_\epsilon(.)$ are generic here. Feel free to use multi-parameter densities, if you would like to solve the problem with a chosen distributional assumption for $f_\epsilon$


title: How can I calculate central moments of a joint pdf?

text: Let's say I have two signals $x_1$ and $x_2$, each having $N$ samples, i.e.: $$ x_1 = \{ x_{11}, x_{12}, ..., x_{1N} \} $$ $$ x_2 = \{ x_{21}, x_{22}, ..., x_{2N} \} $$ The signals are both zero-mean. (Each signal's mean has already been subtracted from all of that signal's elements.) Since the signals are both zero-mean, their moments and central moments are equal. How can I calculate central moments of the joint probability density functions of $x_1$ and $x_2$? In other words: How can I calculate $E[x_{1}^p x_{2}^q]$ for all positive integer values of $p$ and $q$? I know how to compute $E[x_{1}^p x_{2}^q]$ when $p = q = 1$, but that's about it. I don't know how to do the calculation for larger values of $p$ and/or $q$. The signals $x_1$ and $x_2$ are NOT statistically independent. If this were the case, it would be easy to compute as, for statistically independent signals: $E[x_{1}^p x_{2}^q] = E[x_{1}^p] E[x_{2}^q]$. This is not what I want. Thank you very much for your help. PS: Even the formulae which can be used to calculate/approximate these central moments (for finite samples) would be helpful, as I haven't been able to find them yet. It doesn't matter if there aren't any ready-to-use out-of-the-box Matlab functions, I'll then take care of coding the formulae myself. **UPDATE** When talking about statistical independence, does $E[x_{1}^p x_{2}^q]$ mean the first (central) moment of $x_1$ to the power of $p$, times $x_2$ to the power of $q$??? If this is the case, an example would be: $$ x_1 = \sin(x)$$ $$ x_2 = \cos(x)$$ $E[x^2 y^3]$ is the first central moment of $\sin(x)^2 \cos(x)^3$, which is the mean of $\sin(x)^2 \cos(x)^3$ [for finite samples] Can someone please verify this? Thank you.


title: Best PCA Algorithm For Huge Number of Features?

text: I previously asked this on StackOverflow, but it seems like it might be more appropriate here, given that it didn't get any answers on SO. It's kind of at the intersection between statistics and programming. I need to write some code to do PCA (Principal Component Analysis). I've browsed through the well-known algorithms and implemented this one, which as far as I can tell is equivalent to the NIPALS algorithm. It works well for finding the first 2-3 principal components, but then seems to become very slow to converge (on the order of hundreds to thousands of iterations). Here are the details of what I need:   1. The algorithm must be efficient when dealing with huge numbers of features (order 10,000 to 20,000) and sample sizes on the order of a few hundred.   2. It must be reasonably implementable without a decent linear algebra/matrix library, as the target language is D, which doesn't have one yet, and even if it did, I would prefer not to add it as a dependency to the project in question. As a side note, on the same dataset R seems to find all principal components very fast, but it uses singular value decomposition, which is not something I want to code myself.


title: Aspect data in linear regression

text: I have a dataset of various ecological variables on which I want to run linear regression. The variables are continuous, but also include aspect data (sun exposure), in grades. My problem is that the aspect values ranges from 0 to 400, with 0=400=north. How can I include "cyclic" values in the regression? Should I keep them as they are, or cos/sin-transform them? Does anyone of you have experience with these problem? I have had difficulties to find any reference papers... All ideas and references are welcome and appreciated! I hope I was not too unclear, if you need any further information, I'd be happy to provide, as far as I can.


title: Find matching/common questions within a set of surveys

text: I'm currently working on cleaning/preprocessing a bunch of survey data from a collection of similar but distinct surveys. In order to combine the survey results this involves, among other things, identifying common (as in shared) questions (which might, obviously, not be identically phrased). While the data set is small enough[1] to do this manually, it got me curious about what models/methods might be appropriate for automating this. Some quick searches revealed nothing, though this might just be do to not knowing the appropriate terminology. I was thinking about applying a clustering method on the collection of all questions, using some similarity measure taking into account both the survey answers (which can be of a number of different types, though in this case primarily open ended text, multiple choice, and 1 to 5 scale) and the question text. The obvious problem with this is that each survey contains a large number of similar questions (e.g. a slew of "How much do you like ___ on a scale 1 to 5"). My question is therefore whether this approach can be salvaged (for example by first identifying question types and then matching them between surveys) and/or whether there are other, more appropriate, solutions. EDIT/CLARIFICATION: My problem is one of finding similar/same questions among the surveys. For example, while multiple surveys contain a question asking for year of birth, these questions might be worded differently (e.g. "Born:" vs. "Year of birth:" vs. "When were you born?"). EDIT/ADDITIONAL BACKGROUND: The surveys are simple questionnaires, all sent to the same set of schools. The exact wording of the questions have changed over time (the survey spanning multiple years) and might differ depending on the exact audience (students, teachers, management e.t.c.). All this being said, I will most probably do this manually, and my reason for asking is mostly curiosity, so feel free to relax/add assumptions. [1] About 1400 "questions" in 16 surveys, though most of these are "artifacts" of the way the survey software exports multiple choice questions, the real value being something like 10 to 20 questions per survey.


title: Using PCA to merge and grade correlated items

text: I have a real estates' condos sold dataset with the following fields   1. DOM: Date on the market   2. sellPct: Percentage difference between the original and final price.   3. other fields such as Exposure( S,E,N,W), Parking (Yes,No), number of rooms,... I am trying to combibe DOM and sellPct into one score measure (between 0-100). For example a good condo is assumed to have small DOM and sellPCt, so it should have a high score like 85 for example.   My approach was to use PCA to combine/scale DOM,sellPct into one measure, I then use kmeans to divide the PCA components into 4 clusters (Very Good, Good, Fair, Bad). From here on, i am not sure how to map these clustered pca scores to a human readable score (something between 0-100 for example)?   To ultimate goal is to use these scores in a regression model with other fields as predictors.


title: Is there LASSO type model in which only some of the regressors are regulated?

text: I am a beginner with LASSO. Is there any way (paper/code) to perform a LASSO type model in which only some of the explanatory variables are regularized? Or imposing different regulation parameters to several explanatory variables groups?


title: R package random forest: can corr.bias be harmful?

text: There exists an option corr.bias for making regression analyses with the R package randomForest. The manual warns “Experimental. Use at your own risk.” What does corr.bias precisely do and why can it be dangerous to use it? As far as I know the bias correction is based on a linear regression between the out- of-bag-prediction and the observed values which is used for correction of all predictions made by the model, but this does not appear so harmful to me.


title: The "Risk" game dice problem

text: This definitely sounds like a homework, but I assure you that it's not. You're probably familiar with the Risk game. Now, friend of mine rolled 3 aces in one hand. I reckon that probability of such event is $$C(n,k) = \frac{(n-1+k)!}{(n-1)!k!}$$ so that's $8!/(5!*3!) = 56$, so the probability is $1/56$. Am I correct? Problem starts here: he rolled 3 aces in his 2nd attack, so he asked me: now, tell me 'bout the odds, you do statistics! And I must admit that I'm still stunned by his tremendous attacks (he lost 3 tanks both times). Is the probability of such two consecutive events $1/56^{2}$?


title: What would be a better alternative to two balanced latin squares (human subject experiment)?

text: Instead of running a certain balanced latin square twice (if you have _m_ experimental conditions but 2 _m_ subjects), what would be a better approach (to cover twice as many of the possible orderings and get the maximum protection against ordering effects etc.)? ...or in general terms, what is a better alternative to running a certain balanced latin square _n_ times in a psychological research experiment with _n_ _m_ subjects (where _n_ is an whole number)? Is there anything better than seeding the first columns (or rows) of other squares randomly?


title: Cox model with missing covariates

text: I was wondering if any of you have experience in fitting a Cox model to data with missing covariates. Do you know of any reference addressing the issues associated with Cox regression with missing covariates? I know of the approach using the Expectation Maximization (EM) algorithm, but I want to know if there is a publication comparing the different methods out there.


title: What do additional arguments indicate in  the simple.sim function in R?

text: In R suppose that we have the following code:                require("UsingR") # required for "simple.sim" function below      f= function(n=100,mu=0,sigma=1) {          nos = rnorm(n,mu,sigma)          (mean(nos)-mu)/(sigma/sqrt(n))         }              simulations = simple.sim(100,f,100,5,5)         hist(simulations,breaks=10,prob=TRUE)      Does this basically simulate 100 normal random variables with mean 5 and standard deviation 5? In the simulations argument, what do the 5's indicate?


title: Which Analysis to use? Likert Scale

text: I have 3 independent variables (1-5 Likert Scale) questions and I want to check how well these three can predict/explain my DV (1-5 Likert scale) The three independent variables are: 1. Quality of information 2. Accessibility of staff 3. Quality of technical advice My DV is: Overall evaluation of service center All variables are ordinal (1 = low... 5 = high) Which analysis would be appropriate to run here? I would prefer an easy approach and I think Ordinal Logistic Regression is way too complicated. Can I use a Linear Regression? Basically, I want to be able to say that (for example) "quality of technical advice" is better at predicting "overall evaluation" than "Accessibility of staff" Also, I have a 0 value on all variables ("No opinion", so in fact all are measured on 0-5 Likert scale)). How should I treat this variable? Can I replace the 0s with the mean of the observations? Many thanks! Fredrik


title: Finding comparison metric/methods

text: I want to compare the performance of two algorithms. Please suggest me some metrics/methods to compare between the algorithms. Below is given a short description of the output of the algorithms. Suppose there are 100 products (product1, product2,...,product100) each of which has a true rating (true_rate1,true_rate2,...,true_rate100) associated with them. Both the algorithms output a list of best commenter/rater. There are 10,000 raters/reviewers. Not all of them rate all the products. For example, **User1** may give _product1,product3_ , and _product98_ ratings _user1rate1, user1rate3, user1rate98_ respectively. **User2** may give _product5,product7, product9,_ and _product10_ ratings _user2rate5, user2rate7, user2rate9,and user2rate10_ respectively. algorithm 1 returns a list of 100 top raters: alg1_best_rater1,...,alg1_best_rater100, algorithm 2 returns another 100 top raters: alg2_best_rater1,...,alg2_best_rater100. Now we want to decide which algorithm is better. A solution is to calculate and compare the total and avg sum of error (or squared error). the error by user1 is (abs(true_rate1- user1rate1)+ abs(true_rate3- user1rate3)+ abs(true_rate98-user1rate98))/3. Another possible solution could be to use the measure of dispersion of the errors: do those errors wildly fluctuate or do the tend to cluster around the mean. Is there other possible approaches to measure the performances of the algorithms?


title: basic anova for all columns

text: I would like to perform an anova for all the columns in my data set. I tried this:               aov.ex2 <- aov(base~1*2*3,data=data1)      summary(aov.ex2)      results in               8   1   4.83E+08    482880306   29.802  7.90E-08    ***     2   1   7.34E+07    73427001    4.532   0.03381 *     7   1   5.39E+07    53881269    3.325   0.068875    .     5   1   1.53E+07    15340263    0.947   0.331062         0   1   1.39E+07    13937464    0.86    0.354179         3   1   7.05E+06    7053679 0.435   0.509716         6   1   5.99E+06    5990139 0.37    0.543473         3   1   1.53E+06    1526782 0.094   0.759008         1   1   5.26E+05    526334  0.032   0.85705       but the results do not match:               aov.ex2 <- aov(base~1,data=data1)      summary(aov.ex2)          aov.ex2 <- aov(base~2,data=data1)      summary(aov.ex2)          aov.ex2 <- aov(base~3,data=data1)      summary(aov.ex2)      results in:               0   1   1.39E+07    13937464    0.536   0.465                Residuals   709 1.85E+10    26021810                         1   1   9.70E+06    9701265 0.373   0.542                Residuals   709 1.85E+10    26027785                         2   1   1.57E+07    15712139    0.604   0.437                Residuals   709 1.85E+10    26019307                                        aov.ex2 <- aov(base~col1+col2+col3,data=data1)       results in               0        1 1.394e+07  13937464   0.549  0.4591         1            1 5.263e+05    526334   0.021  0.8856         2           1 7.343e+07  73427001   2.890  0.0896 .       3           1 1.527e+06   1526782   0.060  0.8064         4          1 7.054e+06   7053679   0.278  0.5984         5             1 1.534e+07  15340263   0.604  0.4374         6            1 5.990e+06   5990139   0.236  0.6274         7             1 5.388e+07  53881269   2.121  0.1457         8            1 4.829e+08 482880306  19.007 1.5e-05 ***     Residuals   701 1.781e+10  25404903                          Could someone explain the difference to me? I got the code from this webpage: http://personality-project.org/r/r.anova.html I am trying to test the my ability to characterize base using the variables in col 1-n I thought Anova gave a significance to the variable ability to correlate with the base. If I can find the significance of each variable working together that would be great too but I don't fully understand the output at the moment. Here is some sample data from the file It is space delimited to be read in with read.table.               data1 <- read.table(file=args[1], sep='\t', header=T, row.names=1)         base    0   1   2   3   4   5   6   7   8     707 1   1389    430 493.6   0.3096  0.3554  525.47  0.3783  170.29  0.1226     708 4   1392    440 472.8   0.3161  0.3397  669.32  0.4808  197.07  0.1416     709 9   1403    437 487.4   0.3115  0.3474  581.1   0.4142  174.45  0.1243     710 49  1469    476 658.1   0.324   0.448   378.88  0.2579  123.59  0.0841     711 4   1428    435 513.5   0.3046  0.3596  503.33  0.3525  151.95  0.1064


title: Time Series Cross Sectional Analysis and Forecasting With R

text: cty time tl Argentina 2009_Q4 3 Argentina 2010_Q1 2 Argentina 2010_Q2 7 Argentina 2010_Q3 7 Argentina 2010_Q4 10 Argentina 2011_Q1 7 Argentina 2011_Q2 7 Argentina 2011_Q3 1 Argentina 2011_Q4 7 Argentina 2012_Q1 5 The data set has around 40 countries with each country having quarterly data for 5 years. As I am new to R, can someone assist me with how to perform a cross sectional time series ARIMA/ARMA models etc. and forecast in R? Some basic codes to prepare the data, to perform analysis and forecast for this series as an example would be helpful.


title: projecting future survival rate

text: I'm working on a customer retention project that predicts the probability a customer is still subscribing to our services at time T. Unfortunately, we only have the most recent two years of customer data. Estimating the probability of survival at end of year 1 and 2 using Kaplan-Meier is straight- forward. But now we need to project the probability of survival at end of year 3, 4, and 5. Is it possible to do such a task using survival analysis or some other technique (preferably using R)?


title: Regularization $L_1$ norm and $L_2$ norm empirical study

text: There are many methods to perform regularization -- $L_0$, $L_1$, and $L_2$ norm based regularization for example. According to Friedman Hastie & Tibsharani, the best regularizer depends on the problem: namely the nature of the true target function, the particular basis used, signal-to-noise ratio, and sample size. Is there any empirical research comparing methods and performance of various regularization methods?


title: Error metrics to assess a stochastic model with Monte Carlo simulations

text: Suppose we have constructed a model of some stochastic system; we are also able to perform Monte Carlo simulations of this system. Now, we have two sets of samples: one from our model and one from the MC-based approach, and we would like to assess the accuracy of our technique. The first two accuracy measures, which come to my mind, are the differences of expectations and variances. Then, one can try to construct the empirical PDFs/CDFs and compare them at certain points, e.g., find the RMSE. What are the most preferable and exhaustive ways to draw this comparison? Which metrics are the first ones to look at, which are a must?


title: Explanation of a censored regression model

text: I'm wondering how to explain the `logSigma` and its p.value in a censored regression model:               require(cenReg)     data( "Affairs", package = "AER" )     estResult <- censReg( affairs ~ age + yearsmarried + religiousness +                                occupation + rating, data = Affairs )     summary(estResult)          Call:     censReg(formula = affairs ~ age + yearsmarried + religiousness +          occupation + rating, data = Affairs)          Observations:              Total  Left-censored     Uncensored Right-censored                 601            451            150              0           Coefficients:                   Estimate Std. error t value  Pr(> t)         (Intercept)    8.17420    2.74145   2.982  0.00287 **      age           -0.17933    0.07909  -2.267  0.02337 *       yearsmarried   0.55414    0.13452   4.119 3.80e-05 ***     religiousness -1.68622    0.40375  -4.176 2.96e-05 ***     occupation     0.32605    0.25442   1.282  0.20001         rating        -2.28497    0.40783  -5.603 2.11e-08 ***     logSigma       2.10986    0.06710  31.444  < 2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Newton-Raphson maximisation, 7 iterations     Return code 1: gradient close to zero     Log-likelihood: -705.6 on 7 Df      My question is how to explain the logSigma and its significant p.value in the above model.


title: Adjust for everything you have in propensity score?

text: I have a methodological question, and therefore no sample dataset is attached. I'm planning to do a propensity score adjusted Cox regression that aims to examine whether a certain drug will reduce the risk of an outcome. The study is observational, comprising 10,000 individuals. The data set contains 60 variables. I judge that 25 of these might affect treatment allocation. I would never adjust for all 25 of these in a Cox regression, but I've heard that you can include that many variables as predictors in a propensity score and then only include the propensity score subclass and treatment variable in the Cox regression. (covariates that will not be equal after prop score adjustment would of course have to be included in the Cox regression). Bottom line, is it really smart to include that many predictors in the prop score?


title: Closed form expression for the quantiles of $\alpha_1\sin(x)+\alpha_2\cos(x)$

text: I have two random variables, $\alpha_i\sim \text{iid }U(0,1),\;\;i=1,2$ where $U(0,1)$ is the uniform 0-1 distribution. Then, these yield a process, say: $$P(x)=\alpha_1\sin(x)+\alpha_2\cos(x), \;\;\;x\in (0,2\pi)$$ Now, i was wondering whether there is a closed form expression for $F^{-1}(P(x);0.75)$ the theoretical 75 percent quantile of $P(x)$ for a given $x\in(0,2\pi)$ --i suppose i can do it with with a computer and many realizations of $P(x)$, but i'd prefer closed form--.


title: Use regressors in an STL forecast

text: I'm using a STL decomposition to make forecasts in R (using the `forecast` package), but I'm not sure how to incorporate my regressors into the model. I'm using the forecast function:               f = forecast(stl(my.data, s.window = 'periodic'), h = 12, method = 'arima')      The function above accepts the `xreg` parameter, but how do I specify the model regressors and the forecast regressors? Is there a function that does this or do I have to do it 'by hand'?


title: Is there a way to estimate the distribution from only a small number of outliers?

text: For example, I was looking at this list of the 93 people who have broken the "10-second-barrier", after reading that sprinter Christophe Lemaitre was the first person of purely European decent to break the barrier, which got me to wondering what the difference between the mean sprinting times for whites vs blacks. Unfortunately, that number is probably not known since it would require making thousands of average people sprint, and even then it wouldn't necessarily reflect the "true genetic" difference, since the people sprinting were not training for sprinting. So if you wanted to measure the genetic component of the difference, it might be more accurate to measure only the fastest people in the world who are equally motivated and have been training for years, and therefore have eliminated the non-genetic disadvantages, at least that would be my theory. So if you could get a list of say the top 1000 fastest times in the 100 meter dash, and say 20 people on that list are white, could you use that data to give some estimation on what the full distributions look like and/or find what the mean of those distributions are? How? * * * QUESTION is above^^ this is just some rambling: I would guess that if you were trying to find the difference in means between blacks and whites 100 meter sprinting times AFTER everyone in each population had trained for years, lost excess weight ect, i.e. you only want to measure the genetic difference in maximum potential, then measuring average people will not be the way to go, since none of them will have trained to reach their maximum potential, thus there will be many non-genetic factors causing differences, and the other problem is that the difference between trained and untrained, may not be the same, so if you want to measure difference in max potential, it would be better to look at the 1000 fastest, rather than 1000 average people. Also, the data for the 1000 fastest is very high quality, since it was done with laser timing and under official supervision, whereas data gathered from some fitness survey done at a few high schools would probably be of low quality.


title: Checking whether or not a variable has impact

text: For a statistics assignment, I've been given a data set (regarding drug prevention) and a few questions. One of the questions is to check whether or not the choice of treatment, treatment A or treatment B, impacts the odds of staying drug-free. We're told to be cautious about possible confounds, like a having history of intravenous drugs use. Now, my problem is, I'm not really sure how to check if this has an impact or not, and even less sure how to keep the confounds in mind. The given variables in the dataset are: Age, beck depression score, history of IV drug use (never, in the past or recently), number of previous treatments, race, duration of treatment (long/short), which treatment (A/B), and whether or not they were still drug-free after 12 months.


title: Test data for differences

text: I am trying to figure out the best way to test my data for differences in R. The data I have look like this:               domain               variable    lifecycle     162                      CSP1         5     162               old Species        10     162                        BD        16     162                       DBH        14     162                    Height        16     162           Height comments        16     162            Other comments         4     162                      Year         5     123                   species        12     123                  den bran        14     123                  den core        14     191                      CSP2         4     191                   C kg ha         5     185                  location         5     185                   biomass         6                  82                   CSP num         6       82                  depth lb         8      82                       C t         7     190                     Layer         8          190                Dry weight         9     204                 item plot         5     204               item volume         7     204 within central plot small         7     205                      plot         3     205        successional stage         3     205            rarefy hundred         3      Now I want to use the "variable" column as factorial group and the "lifecycle" as the characteristiom to compare. Because the data is not normally distributed and I have more than two groups to compare I chose a Kruskal- Wallis test, which shows no significant differences.               Kruskal-Wallis rank sum test            data:  lifecycle by variable     Kruskal-Wallis chi-squared = 25, df = 25, p-value = 0.4624      My question is now if I can use the Kruskal-Wallis in this case to find differences in "lifecycle" for the "variable" groups, or is the sample size to small with only one value per group? And if it is not usable here in this case what would be a better solution? EDIT: I did the Kruskal-Wallis now with the domain as grouping factor. Which gives me a significant difference. The interpretation would be that there are domains with significant more or less impact than other. After that I performed a post-hoc test to see between which of the domains the differences occur. But supprisingly the paired wilcoxon test with bonferroni correction showed up no significant grup differences. How could this be?


title: What model can be used when the constant variance assumption is violated?

text: Since we can not fit ARIMA model when the constant variance assumption is violated, what model can be used to fit univariate time series?


title: No minimum in cross validation error

text: I am evaluating both non-linear LASSO and SVM for a high-dimensional non- linear dataset. In order to do that I am applying 10-fold cross validation. The training error decreases as usual as the model gets more complex. However, the validation/error test always decreases as well, and no minimum can be found. My question is: is it possible that no minimum for the test error is found during cross-validation? How to interpret the results in that case?


title: Quantiles of a compound gamma/negative binomial distribution

text: Are there any formulae for the quantiles of a compound gamma/negative binomial distribution? That is, suppose we have $$N \sim \text{NegBin}(\alpha, \lambda)$$ and conditional on $N$, $$Y = \begin{cases} 0 & \text{if } N = 0 \\ \sum_{i=1}^N X_i & \text{otherwise} \end{cases}$$ where each of the $X_i$ is iid $\text{Gamma}(\mu, \sigma)$. Given values of $\alpha, \lambda, \mu$ and $\sigma$, I'd like to calculate the quantiles of $Y$. Background: I'm modelling ATM withdrawals per day. The number of withdrawals I assume has a negative binomial distribution, and the amount of each withdrawal has a gamma distribution. I can estimate all the parameters from the data, and now I want to estimate how much cash to keep in the ATM. Keeping too little cash risks not having enough to meet demand, but too much is also bad. So, I want to estimate how much cash is needed to meet demand for a given level of cash-out risk. That is, I want an upper quantile of the compound distribution, say at the 99%, 99.9% or 99.99% level. Right now I'm using simulation to get an estimate, but this takes a long time especially for very high quantiles. Eg if I want to get the 99.99% quantile, I need at least 10,000 replications, and ideally much more than that; I also need to repeat this for each ATM. So some kind of analytical approximation would be great. (I realise that the actual distribution is probably more complicated than a simple negbin/gamma combination. For example the observed withdrawals would be right censored since if an ATM runs out of cash then there's no more withdrawals for the day. I'm keeping things simple to start with.)


title: Control for previous day values when doing across day analysis with generalized estimating equations?

text: I am a graduate student and have a question regarding generalized estimating equation analysis. When conducting across day analysis with GEE, with variable a on day $x$ predicting variable $b$ on day $x+1$, is it necessary to control for previous day's values for variable $b$, or is this correlation already accounted for with the GEE technique? Your comments and references would be much appreciated.


title: Non-significant group effect, but LSD post hoc is significant. Why?

text: I am performing repeated ANOVAs in SPSS and came across something I don't understand. Basically, the main effect of group (I have three groups) is not significant (p = .062), but the LSD _post hoc_ test indicates that one of the pair-wise comparisons is significant (p = .032). How is that possible?


title: How to handle sets with less than $k$ elements when using a single hash function to minhash?

text: A minhash implementation with multiple hash functions can easily handle comparisons between sets with a vastly different number of elements because the denominator of the unbiased estimator $k$ is based on the number of hash functions, which is always fixed. However, when calculating a minhash using a single hash function, you need to choose a value $k$ which represents the minimum $k$ values when all the items in $S$ are run through a single hash function. How can you process a set $S$ that has less than $k$ elements and get an unbiased estimate? I see two possible approaches. The first is to simply void any sets that don't have more than $k$ elements for comparison. This is undesirable because even for sets with a small number of values, I'd still like to compare it to other sets (which I can do with the multiple hashing variant). The second is to get the intersection of $h_k(S_1)$ and $h_{k}(S_2)$ where $|S_{1 and/or 2}| < k$ as you would normally. This seems like it would work, but for small sets that are very similar (if not outright duplicates of each other) it seems it would provide erroneous values; you'd possibly miss matches, but at least you won't get false positives for reasons other than your hash function not being a perfect hash.


title: How do I identify this simultaneous equations model?

text: I have the following model (which is in this form is not identifiable if the $y$'s are indeed endogenous): (1) $y_1 = a_0 + a_1y_2 + a_2y_3 + \boldsymbol{Xa} + \boldsymbol{u}$ (2) $y_2 = b_0 + b_1y_1 + b_2y_3 + \boldsymbol{Xb} + \boldsymbol{v}$ (3) $y_3 = c_0 + c_1y_1 + c_2y_2 + \boldsymbol{Xc} + \boldsymbol{w}$ where $\boldsymbol{X}$ is a $n\times k$ matrix of exogenous variables $\boldsymbol{a, b, c}$ which are each $k\times 1$ vectors and the $\boldsymbol{u,v,w}$ are the random noise. Since every endogenous variable (i.e., $y_i$ for $i=1,2,3$) is a function of the other $y$'s and each equation contains the whole set of exogenous variables (i.e., $X$) one cannot identify the equations given above. So I checked my empirical data and the theory I have in mind, and I cannot relax the dependency on which of the endogenous variables are related to each other. My next guess would be that I could check the exogenous variables to find a way to exclude at least 2 of them in each equation and each of these 2 should in the best case not be excluded in any of the other equations. This would be just amazing, but I won't get out 6 different exogenous variables anyway. How should I deal with this problem?


title: R software implementation of combining mixed treatment comparisons and meta-regression

text: I am currently conducting a meta-analysis in which I need to use a mixed treatment comparison method. As I understand it, this method works in the following way: Say you have a group of studies that make the following set of treatment comparisons:   * Intervention 1   * Intervention 2   * Intervention 3   * Control You are interested in all possible comparisons between these treatments. So, not only are you interested in intervention 1 versus control, intervention 2 versus control, and intervention 3 versus control, but also intervention 1 versus intervention 2, intervention 1 versus intervention 3, etc. The problem occurs in that not all of the studies in your meta analysis include each intervention type. So, while study 1 may have tested intervention 1, intervention 2, and a control group, study 2 tested intervention 2 and intervention 3 versus a control group. And so on. Mixed treatment comparisons (Caldwell, Ades, & Higgins, 2005; Lu & Ades, 2004; Mills et al., 2011) arose as a way of using the indirect information from your sample of studies to estimate the magnitude of the missing comparisons. For my study, I am interested in how several different moderators affect the magnitude of the various treatment comparisons. I stumbled across a paper (Nixon, Bansback, & Brennan, 2007) that combines the mixed treatment comparison method with meta-regression. My problem is finding a good software implementation for this method (preferably an implementation in R, since I'm most familiar with R). As far as I can tell, the metafor package isn't able to handle mixed treatment comparisons. Does anybody know whether there's a package out there that's able to handle both mixed treatment comparisons and meta-regression? Thanks in advance! Patrick


title: GPS coordinates rotated stats calculation

text: I have a set of GPS coordinates along with their heading angle. These coordinates are within 15 m of each other and so I can use cartesian math to determine the center, std deviation of coordinates etc. I am then trying to find a subset of GPS coordinates that are "above" or "below" the center coordinate. This works well if the road is aligned with the X and Y axis (see fig a below). If the road is at an angle to the X-Y axis (fig b below) how can I determine the center, SD etc for the rotated road? ![Two roads with GPS coordinates](http://i.stack.imgur.com/pQF0t.jpg)


title: How can I estimate the shape of a curve where the predictor variable is right censored interval variable?

text: Are there methods available to conduct interval-regression style analysis where a predictor variable is the interval variable and the outcome variable is a dichotomous variable? To give a bit of background, I would like to plot the relationship between two variables in a dataset where the x axis is a continuous measure has been collected in an ordinal format. The category widths are non-uniform (and increasing). To give a concrete example, the categories are about expenditure on a category of goods and were (in dollars): 0-50; 51-100; 101-250; 251-500; 501-1000; > 1000. The y-axis is a yes/no style variable. I am particularly interested in determining whether the relationship is likely to be linear or exponential. Are there any methods for doing this, other than arbitrarily picking a point in each x-axis category? Obviously, the final category (> 1000) is the most problematic. Bonus points (metaphorically speaking) if you can point me in the direction of an R package. If it helps, I am happy to look for an answer that requires a strong distributional assumption (e.g. that the data are zero-inflated gamma or log- normal distributed).


title: Clustering with dependent features

text: I have $n$ observations which I want to cluster based on $p$ features. Those $p$ features are dependent to each other (or correlated). If I use Euclidean distance, the weights of all features would be the same, independent from the correlations between them. But what I want is to decrease the weight of a feature in clustering process if this feature is correlated with many other features. So, I want to arrange the weight of a feature so that it is somewhat inversely proportional to the amount of correlation of this feature with other features. I know that Mahalanobis distance can be used to take care of the correlations between the features, but I am not sure what exactly Mahalanobis does with those correlations. So, I wonder if it does what I need to do (arranging the weight of a feature so that it is somewhat inversely proportional to the amount of correlation of this feature with other features). Thanks!


title: Predicting and calculating test mse using cforest R

text: I'm new to the `cforest` package and am trying to create a cforest model to predict a new test set and calculate the model test MSE. My data is split into d.train and d.test. If I use randomForest, the code for predicting and determining the MSE is as follows:               set.seed(1)     rfp.data <- randomForest(d.train$Y.Weight ~ . -d.train$Y.Weight, data = d.train, mtry = sqrt(p), importance = TRUE)     yhat.rfp <- predict(rfp.data, newdata = d.test)     mean((yhat.rfp - d.test$Y.Weight)^2)     #MSE  = 9.648027      My code for cforest is similar but it doesn't look like cforest is returning a model like randomForest.               cf.data <- cforest(d.train$Y.Weight ~. -d.train$Y.Weight, data = d.train, controls =   cforest_unbiased(mtry = p/3, ntree = 500))     yhat.cf <- predict(cf.data, newdata = d.test)      > Error in model.frame.default(formula = ~(Wrist.Diam + Wrist.Girth + > Forearm.Girth + : variable lengths differ (found for 'data$Y.Weight') I also tried using the entire dataset and then using subset for the training set and received the same error.               cf.data <- cforest(d.train$Y.Weight ~. -d.train$Y.Weight, data = d.train, controls = cforest_unbiased(mtry = p/3, ntree = 500))     yhat.cf <- predict(cf.data, newdata = d.test)      > Error in model.frame.default(formula = ~(Wrist.Diam + Wrist.Girth + > Forearm.Girth + : variable lengths differ (found for 'd.train$Y.Weight') I haven't been able to find too many examples of cforest other than the cran party pdf. Any pointers to an example would be greatly appreciated as well.


title: Equivalence of simple random sampling with- and without- replacement

text: We know that in simple random sampling without replacement the probability of drawing n units from a population of N units(where ordering is considered) is 1/(N)n or 1/(N n) (where ordering is not considered). How is this equivalent to drawing units independently one by one and with replacement where each unit has probability of selection 1/N at each draw, until n distinct units are drawn and then discarding the repeated units?


title: How do I show that the final point estimates of a variable satisfy the invertibility conditions?

text: I cannot seem to figure out how to `"identify the final point estimates of theta1 and theta2 and verify that they satisfy the invertibility conditions... [and] identify the final point estimate of delta"`. All of these are obviously values in the model $Y_t=\text{ ...}$ I am given a MINITAB output of the estimation, diagnostic checking, and forecasting of daily readings using the model. I want to be able to figure it out myself, so I wont post the actual tables I am given to make sure I understand it when I have to do it again later. I know how to find point estimates...but how does the 'final' point estimate differ? And how is it calculated? Lastly, how do I show that they satisfy the invertibility conditions?


title: NNLS in R correlation of variables

text: I am trying to fit some variables which are high collinear using nnls(non negative least squared routine) in R. I have about 1 million observations and 200 highly collinear variables for the coefficients to be found. My problem is that in the final result, the algorithm seems to select variables that have pretty low correlation with the independent variable, it seems to throw the variables that have a high correlation in favor for the ones with a low correlation. I dont seem to understand why. Has anyone run into something similar, any way I can modify to pick variables in order of correlation, I am using the standard NNLS package in R.


title: In R, what does anova(f1, f2, f3) compare?

text: I'm trying to understand the reason why anova(f1, f2, f3) and anova(f1, f2) gives me different result while anova(f1, f2, f3) and anova(f2, f3) gives me the same. Here's the code:               > data(swiss)     > fit1 <- lm(Fertility ~ Agriculture, data=swiss)     > fit3 <- lm(Fertility ~ Agriculture + Examination + Education, data=swiss)     > fit5 <- lm(Fertility ~ Agriculture + Examination + Education + Catholic + Infant.Mortality, data=swiss)          > class(fit1)     [1] "lm"     > class(fit3)     [1] "lm"     > class(fit5)     [1] "lm"          > anova(fit1, fit3, fit5)     Analysis of Variance Table          Model 1: Fertility ~ Agriculture     Model 2: Fertility ~ Agriculture + Examination + Education     Model 3: Fertility ~ Agriculture + Examination + Education + Catholic +          Infant.Mortality       Res.Df    RSS Df Sum of Sq      F    Pr(>F)         1     45 6283.1                                       2     43 3180.9  2    3102.2 30.211 8.638e-09 ***     3     41 2105.0  2    1075.9 10.477 0.0002111 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          > anova(fit1, fit3)     Analysis of Variance Table          Model 1: Fertility ~ Agriculture     Model 2: Fertility ~ Agriculture + Examination + Education       Res.Df    RSS Df Sum of Sq      F    Pr(>F)         1     45 6283.1                                       2     43 3180.9  2    3102.2 20.968 4.407e-07 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          > anova(fit3, fit5)     Analysis of Variance Table          Model 1: Fertility ~ Agriculture + Examination + Education     Model 2: Fertility ~ Agriculture + Examination + Education + Catholic +          Infant.Mortality       Res.Df    RSS Df Sum of Sq      F    Pr(>F)         1     43 3180.9                                       2     41 2105.0  2    1075.9 10.477 0.0002111 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      To me, this is read like this:   * `anova(fit1, fit3, fit5)` shows me P value 8.638e-09 for the comparison of f1 and f3.    * `anova(fit1, fit3)` shows me 4.407e-7 for the same comparison. And this doesn't make sense!   * `anova(fit1, fit3, fit5)` and `anova(fit3, fit5)` shows the same P value for the comparison of fit3 and fit5 and it's 0.0002111. Probably I missed something. What is it?


title: Why are the following probabilities equal to each other?

text: The following is a Bayesian model, there are variables $Y_1,Y_2 \in \{ 0,1 \}$ and variable $X \in \{ 0, 1\}$ and $\theta$ between $0$ and $1$. The model is:   1. draw $\theta$ from Beta distribution with parameters $a$ and $b$   2. draw $Y_1$ from Bernoulli $\theta$   3. draw $Y_2$ from Bernoulli $\theta$   4. draw $X$ as following: with prob. $\theta$, $X = Y_1\,\, XOR\,\, Y_2$ and with probability $1-\theta$ $X = 1 - (Y_1\,\, XOR\,\, Y_2)$. The joint distribution is: $$p(\theta,Y_1,Y_2,X) = p(\theta |a,b)p(Y_1 | \theta)p(Y_2 | \theta)p(X | Y_1,Y_2,\theta)$$ I calculated the joint $p(X,Y_1,Y_2,\theta)$. The interesting part which I don't have an intuitive explanation for is why $$p(Y_1 = 0 | Y_2 = 0, X = 0) = P(Y_1 = 1 | Y_2 = 0, X = 0) = 1/2$$ which is what I get from my calcs. It seems counter-intuitive, because I would expect $a$ and $b$ to bias this posterior. The probabilities do not equal each 1/2, for example, if we compare $$p(Y_1 = 0 | Y_2 = 0, X = 1)$$ and $$p(Y_1 = 1 | Y_2 = 0, X= 1)$$. Does anyone have an intuitive explanation why the above probabilities, each equals 1/2?


title: Uncertainty consistency?

text: 1. The problem statement, all variables and given/known data I am given a set of x and y values x: (1,2,3,etc.) y: (1.2,2.2,3.1,etc.) with a given uncertainty and am asked a) find the best fit b) at what prob can you rule out a 5% higher slope c) is the stated uncertainty consistent with the data? I can find the best fit relatively easily by minimizing chi-sqd and setting the derivatives to 0. I am confident in my result as it matches with the graph given by excel. For part 2 I tried putting a higher value for the given slope into the chi-sqd equation and checking a chart but that didn't give me a reasonable answer (not completely sure about the degree of freedom) I am not sure how to approach the third part. I have computed the errors in the fitting coefficients but that doesn't seem to play into it. Thanks,


title: Is there a data decay method that decays data based on how much new data it receives?

text: So here's the situation: There are n balls. Each ball has been voted on a scale of 1-N by an arbitrary number of people. Every day (or every week, or more), some/most/all of the previous voters and some new voters are able to vote on the balls again. The thing is, I don't want all previous votes to all disappear because the votes that were valid yesterday are still valid today. The re-votes/new votes will just be weighted "more". There can be more OR fewer people who vote during the new voting session. What is a good method of data decay for this? I have thought of taking the square root of "yesterday"'s votes every day/week, but that doesn't seem like it would work well if I have lots of voters in one day because the square root of the sum of the votes doesn't decay that "much" and if I have few voters, the square root of the sum of the votes can decay less and so the new, initial votes, will skew the result set. I could be wrong; I'm not sure.


title: Is there an intuitive interpretation of a negative variogram "nugget" value?

text: A variogram plots the variance of the difference between sample pairs on a field (any dimensionality) against spatial separation (the "lag") of those samples. The extrapolation from observed small-lag variances to a zero-lag variogram value is generally termed the "nugget" value. Variograms I've come across before (geostatistics) always had positive "nuggets", and these values were generally interpreted as arising from (and being some sort of guide to the magnitude of) a combination of measurement error and sub-observation-scale variation in the observed field. However, I'm just now looking at variograms of some imaging data and any reasonable extrapolation of the low-lag variogram clearly results in a **negative** "nugget" value. Of course it's always possible to "force" a non- negative nugget by, say, fitting an exponential model... but this doesn't look very convincing when plotted. Is there any simple interpretation of negative nugget values (or a simple explanation of how they arise) as there is for positive ones ?


title: Interpretation of p-values & reduction of model deviance for factors in anova() vs summary()

text: My data has 3 major inputs: `BLDDAY` (a factor), `BLDMNT` (a factor), and `D_BLD_SER` (days as an integer variable). The output is whether input variable has any impact on failure. My model is: `model = glm(FAILED~BLDDAY+BLDMNT+D_BLD_SER, family="binomial", data=data_list)`. (I used UCLA's statistics help site's guide to logistic regression in R to build this model.) Output:               Deviance Residuals:          Min       1Q   Median       3Q      Max       -1.3282  -0.9123  -0.8128   1.4056   2.1053            Coefficients:                       Estimate Std. Error z value Pr(>|z|)         (Intercept)     -0.7672583  0.1317132  -5.825 5.70e-09 ***     BLDDAYMonday    -0.1545646  0.0839380  -1.841  0.06556 .       BLDDAYSaturday  -0.1257976  0.2028259  -0.620  0.53511         BLDDAYSunday    -0.1183008  0.1868713  -0.633  0.52669         BLDDAYThursday  -0.2007452  0.0772653  -2.598  0.00937 **      BLDDAYTuesday    0.0480453  0.0758603   0.633  0.52651         BLDDAYWednesday -0.0358585  0.0760027  -0.472  0.63707         BLDMNTAug        0.3009445  0.1405545   2.141  0.03226 *       BLDMNTDec        0.5562170  0.1338467   4.156 3.24e-05 ***     BLDMNTFeb        0.3334978  0.2133475   1.563  0.11801         BLDMNTJan        0.4076504  0.2277978   1.790  0.07353 .       BLDMNTJul        0.1306585  0.1415302   0.923  0.35591         BLDMNTJun       -0.0357361  0.1428105  -0.250  0.80241         BLDMNTMar        0.4570491  0.1949815   2.344  0.01907 *       BLDMNTMay       -0.2292620  0.1614577  -1.420  0.15562         BLDMNTNov        0.3060012  0.1334034   2.294  0.02180 *       BLDMNTOct        0.2390501  0.1341877   1.781  0.07484 .       BLDMNTSep        0.2481405  0.1384901   1.792  0.07317 .       D_BLD_SER       -0.0020960  0.0003367  -6.225 4.82e-10 ***          (Dispersion parameter for binomial family taken to be 1)         Null deviance: 10288  on 8182  degrees of freedom     Residual deviance: 10154  on 8164  degrees of freedom     AIC: 10192     Number of Fisher Scoring iterations: 4      The ANOVA table is the following:               anova(model, test="Chisq")     Analysis of Deviance Table     Model: binomial, link: logit     Response: FAILED     Terms added sequentially (first to last)                    Df Deviance Resid. Df Resid. Dev  Pr(>Chi)         NULL                       8182      10288                   BLDDAY     6   20.392      8176      10268  0.002357 **      BLDMNT    11   70.662      8165      10197 9.142e-11 ***     D_BLD_SER  1   43.797      8164      10154 3.642e-11 ***      My questions are:   1. Although the p-values for all three components are less than 0.05, which can be considered as significant, the deviance reduced due to each component is less than 1% of the total deviance. **Normally the interpretation of output like this is input parameter affects output and it's better to use this parameter then using noting.** But does it really make sense of taking this parameter as significant input?   2. The p-values for `BLDDAY` and `BLDMNT` given by `anova()` is the overall p-value, which is significant, but `summary()` gives detailed impact of each factor level. If I consider the p-values for each factor overall `BLDDAY` is significant but individually only `BLDDAYThursday` is significant. I am bit confused not as whether to consider `BLDDAY` as significant input, or Thursday only, or Thursday & Friday both.


title: Testing dispersion of values between condition on a group

text: I may have a tricky question. I have a within design where every participant is tested in 2 conditions. I want to compare the dispersion of the measure, for the two condition, with the hypotheses that one condition has a larger dispersion than the other one. I though about taking the sd to get a measure of dispersion by condition and participant. After that I don't know if I have the right to test for a statistical difference by using a t test on the group (just a paired t test). I'd be grateful if anyone can help me. Thank you.


title: Expectation of the product of two log normal variables

text: I am struggling with a proof, and I am wondering if anyone can help or point me to the right direction. Suppose that we have two variables, $X$ and $Y$, and they follow a multivariate normal distribution with covariance matrix $Q$. Now, suppose that we have two more variables, $Z_1$ and $Z_2$ that follow a normal distribution with covariance $\Omega$, whose elements are, $\Omega_{11} = exp(X)$, $\Omega_{22} = exp(Y)$, $\Omega_{21} = \Omega_{12}=\exp (X/2)\exp(Y/2)$, And the question is, what is the expectation of $\Omega$? For $\Omega_{11}$ and $\Omega_{22}$ it is clear that the expectation is $\exp(0.5Q_{11})$ and $\exp(0.5Q_{22})$ using the properties of log normal variables, but I am not sure how to compute it for the covariance term. Thanks!


title: K-fold vs. Monte Carlo cross-validation

text: I am trying to learn various cross validation methods, primarily with intention to apply to supervised multivariate analysis techniques. Two I have come across are K-fold and Monte Carlo cross-validation techniques. I have read that K-fold is a variation on Monte Carlo but I'm not sure I fully understand what makes up the definition of Monte Carlo. Could someone please explain the distinction between these two methods?


title: Issues when using interaction term with a lagged variable in R

text: Today I tried to estimate models using both `plm` and `pgmm` functions in the plm R package, with an interaction between `X1` and `lag(X2, 1)`. And I notice two issues. Let $Y=b_1 X_1 + b_2 X_2 + b_3 X_1 X_2 + e$ be our model.   1. When using `plm`, I got different results when I coded the interaction term with `I(X1 * lag(X2, 1))` and when I just saved this multiplication `X1 * lag(X2, 1)` in a different variable of the dataset and then used it in the regression.   2. With `pgmm` it is not even possible to run a formula which contains `I(X1 * lag(X2, 1))`. How can I pass such interaction?


title: What's the null hypothesis in a one-sided Kolmogorov-Smirnov test?

text: Suppose I use K-S to figure out if the CDF of $X$ is greater than the CDF of $Y$. I get the statistic $D^+ = \max_u\{C_x(u) - C_y(u)\}$ where $C_x$ is the ECDF for $X$ and similarly for $C_y$. There are various tables which translate this into a p-value. Intuitively, if I'm testing if $C_x < C_y$ it seems like the null hypothesis should be $C_x \geq C_y$. (If it were $C_x = C_y$ then I don't understand how it would differ from a two-sided test.) Yet for my data set I find low p-values for both the proposition that $C_x < C_y$ and $C_x > C_y$, leading me to believe that the null hypothesis is not the complement of the alternative. Am I completely misunderstanding something?


title: R test a regression with only statistically significant variables

text: I'm performing a multiple logistic regression with a large amount of columns. I want to set all the variables over a certain p value to have 0 coefficients, then test the model on the test data. Since there's so many variables I'm using, it's impractical to manually go in and set the coefficients. So I'm looking for something like the following.               log_results <- glm(formula, data, family);     log_results_sig <- get_sig_only(log_results, p value threshold);      This way I can use the results with the predict() function, and export the results easily for use in another program. Additionally, if anyone knows a way to automatically extract the significant variables and refit with them, I would appreciate that too. Thanks.


title: What happens to the constant in Least Squares

text: **Preliminaries** For simplicity's sake assume we are dealing with a 2-dimensional dataset of examples $(x_i, y_i) \in \mathbb{R}^2$ which are split into a training (objects and their label known) and test set (only objects known). Prediction using _least squares_ can be represented using the following matrix procedure:   1. From our training set find the optimal regressor constant $w = Y'X(X'X)^{-1}$   2. Apply the above to successive unlabeled examples from the test set to find their label $\hat{y} = w'x = Y'X(X'X)^{-1}x$ **The Question** After you've finished part 1. of the above procedure how do you find the optimal linear least squares line? If a linear line can be represented using $y = w'x + c$ how do you find the constant $c$?


title: goodness of fit? how to test if a set of simulated data distributions are well-fitted to an empirical data distribution in R or Python?

text: I have a data set, whose elements (say 1,000 elements) are numeric with values between 0 and 1. It is safe to think the elements are independent from each other, i.e, the order of these elements does not matter at all. Also, from the histogram of the data set, I can tell it does not follow any parametric distribution (normal or others). I used a modeling tool (not in R) to simulate 100 realizations of this data set given a set of parameters. In each realization, the number of elements is not necessarily equal to 1,000. My goal is simple: Does the model fit the data well? I plotted the histograms of the simulated data against the empirical data one (using both counts and probability density) and was sort of able to tell the quality of the fit, but would like to have a solid statistical way to quantitatively reach my goal (say showing the p-value of this fit is > 10% or so). Is there a good way in R or Python I can do so? I looked up the standard R fitting methods, such as lm() and glm() in R, but they are for fitting a line to the data, not doing something I would like to have. I believe that this is simply because I am not familiar with R and statistics given this question does not seem to be weird and difficult. Thanks, Ben


title: Hypothesis testing on Poisson Binomial distribution

text: Imagine $N$ baskets each filled with a different number of products chosen from a set $P$. It is possible to have repeated products in the same basket. For instance:               basket 1 = [A,B,B,C]     basket 2 = [A,B,C,D,F,J]     basket 3 = [C]     ...     basket N = [E,E,E,F,F]      Now we ask someone to chose one item from every basket. For instance:               basket 1 = [A,B,B,C] - A     basket 2 = [A,B,C,D,F,J] - A     basket 3 = [C] - C     ...     basket N = [A,E,E,E,F,F] - E      How can we test the preference of this person for a given product? **EDIT** : The problem can be simplified into this one: We have $N$ independent trials where every trial $i$ follows a Bernoulli distribution with probability $p_i$. That is, a Poisson Binomial distribution. Imagine that we observed $k$ successful trials. Is there a way to compute the deviation of the observation respect to the distribution, and how significant is that deviation? (I'm thinking of p-values and z-scores in a normal distribution)


title: Change in output depending on mirt version

text: I have a matrix which I'm trying to run through the `mirt` function of the `mirt` package:               resp.freq <- data.frame(matrix(c(11, 46, 12, 31, 13, 8, 21, 20, 22, 68, 23, 12,                                      31,  1, 32, 12, 33, 11), nrow = 9, ncol = 2,                                    byrow = T))     dados3 <- matrix(rep(resp.freq$X1, resp.freq$X2), ncol = 1)     dados3 <- data.frame(cbind(as.numeric(substr(dados3, 1, 1)),                                as.numeric(substr(dados3, 2, 2))))          require(mirt)     m1 <- mirt(dados3, 1, D = 1, SE = TRUE)     coef(m1)      I used to run this through **mirt 0.5.0** and would get the following end results:               $X1             a1    d1     d2     pars 1.913 0.633 -3.152     SE   0.228 0.181  0.322          $X2             a1    d1     d2     pars 1.659 1.121 -2.488     SE   0.200 0.190  0.261      However, my workstation has been updated and on **mirt 1.2.1** and now I get the following output upon running `m1 <\- mirt(...`:               Iteration: 1000, Log-Lik: -387.670, Max-Change: 0.00007     EM iterations terminated after 1000 iterations.               Calculating information matrix...     Warning message:     In loadESTIMATEinfo(info = info, ESTIMATE = ESTIMATE, constrain = constrain) :       Negative SEs set to NaN.      And `coef(m1)` gives me:               $X1                 a1    d1     d2     par     11.515 2.732 -13.98     CI_2.5     NaN 0.685    NaN     CI_97.5    NaN 4.778    NaN          $X2                a1    d1     d2     par     1.099 0.903 -2.139     CI_2.5  0.712 0.561 -2.611     CI_97.5 1.485 1.245 -1.668          $GroupPars             MEAN_1 COV_11     par          0      1     CI_2.5      NA     NA     CI_97.5     NA     NA      I've read the changelog, but couldn't find a reason for the behavior change. I've tried desperately changing the parameters of `mirt()`, but couldn't reach even a similar level of convergence. What gives? The data doesn't look like it would be inappropriate for this kind of prodecude, am I missing anything?


title: Observed variables in R's lavaan

text: I am running a SEM on lavaan in R. The model consists of latent variables [N1 to N3], two dummy choice variables [du1 and du2] and observed satisfaction variable [SV] (beside the observed variables). Without the dummy choice variables, this is visually what the model should look like (without intercept and residuals):               N1 -> observed vars.     N2 -> observed vars.     N3 -> observed vars.     SV <- beta1 * N1 + beta2 * N2 + beta3 * N3      Utilizing lavaan, this would look something like:               N1 =~ SV + observed vars.     N2 =~ SV + observed vars.     N3 =~ SV + observed vars.      I want to add the observed choice dummies to the model so that it visually becomes:               N1 -> observed vars.     N2 -> observed vars.     N3 -> observed vars.     SV <- beta1 * N1 + beta2 * N2 + beta3 * N3 + beta4 * du1 + beta5 * du2      How should I adjust the lavaan model in order to achieve that? Thank you, Deuterium


title: Interaction term as a dependant variable in LMM with R

text: In a longitudinal study, two groups of subjects have been measured over a period of two years at 6 months intervals. During these measurements subjects have been assessed with a series of $k$ measures ($M_1$ to $M_k$). From other sources we expect to see a difference between the groups and a decline in the performance of only one group. The question is: which of the $M_1$ to $M_k$ measures can pick up the decline of performance over the time. If we had only one measure M, I suppose we could simply get an answer using:               fit.0 <- lme(M ~ Group * Time, random=~1|Subject)      Then we could simply check the significance of the interaction effect. However, I guess looking at $M_1$ to $M_k$ measures one by one might not be ideal. Although differences between the two groups are assumed to be real according other sources, still some form of multiple comparison might be going on. Another way to look at the problem might be to try to address our research problem directly, i.e. trying something like this:               fit.1 <- lme(Group:Time ~ Group + Time + M0 + M1 + ... Mk, random~1|Subject)      This way, I could directly see which $M_1$ to $M_k$ terms is a significant predictor of the interaction term (or I could add the $M_1$ to $M_k$ terms one by one or use some other model selection method). However the main problem is that it does not work! It seems R does not like to have an interaction term as the dependent variable. So here is what I want to know: is the basic idea behind the second approach sound? If yes, is there a way to make `fit.1` work? Comment 1: The problem of finding the measures that can predict decline in longitudinal studies is not that rare, at least not in my corner of the world. For instance consider the case of neuro-degenrative diseases like Parkinson's or MS. From all kind of sources we know that these diseaes progress constantly and slowly. Thus it is interesting to find the measures that reflect the progress of the disease.


title: Rough estimates for training time of deep belief networks

text: I'm still learning about deep learning. However I'm currently interested to know if deep learning architectures scale well or not. Suppose I have a dataset with 1 million training examples, can you give me some rough estimates how much time would it take to train a deep architecture on such a dataset? I know there are lot of hyper-parameters that affect the training time (#layers...etc), but I just want some rough mean estimates since I'm still at the beginning and not aware of all parameters and so on. I just want to know how well they scale to large datasets.


title: Are there techniques to merge two cumulative distribution functions?

text: I'm trying to estimate quantiles from cumulative distribution functions. Given N CDFs are there techniques that can be used to merge them to form another CDF from which a quantile can be estimated ? Thanks everyone for looking at the question. Here is a more concrete example: Imagine a web service being served by N server instances. Each server constructs an empirical CDF of the response times of a request. Now I would like to aggregate or merge the CDFs to find say the 95th percentile response time of the web service. The only data available is the empirical CDF and number of samples used to construct the CDF.


title: Discretizing all attributes vs. discretizing only the class label

text: I have a numeric dataset and I would like to apply classification algorithms on it. So I am using Weka's Discretize filter to convert the numeric values to discrete. However I wonder what can be the implications of converting all the attributes (including the class label) to discrete vs. converting only the class label to discrete. Is there anyone who experienced such a scenario?


title: Relationship between two time series: ARIMA

text: Given the following two time series ( **x** , **y** ; see below), what is the best method to model the relationship between the long-term trends in this data? Both time series have significant Durbin-Watson tests when modelled as a function of time and neither are stationary (as I understand the term, or does this mean it only needs to be stationary in the residuals?). I have been told that this means I should take a first-order difference (at least, maybe even 2nd order) of each time series before I can model one as a function of the other, essentially utilising an arima(1,1,0), arima(1,2,0) etc. I don't understand why you need to detrend before you can model them. I understand the need to model the auto-correlation, but I don't understand why there needs to be differencing. To me, it appears as though detrending by differencing is removing the primary signals (in this case the long-term trends) in the data that we are interested in and leaving the higher-frequency "noise" (using the term noise loosely). Indeed, in simulations where I create an almost perfect relationship between one time series and another, with no autocorrelation, differencing the time series gives me results that are counterintuitive for relationship detection purposes, e.g.,               a = 1:50 + rnorm(50, sd = 0.01)     b = a + rnorm(50, sd = 1)     da = diff(a); db = diff(b)     summary(lmx <- lm(db ~ da))      In this case, **b** is related strongly with **a** , but **b** has more noise. To me this shows that differencing _doesn't work_ in an ideal case for detecting relationships between low frequency signals. I understand that differencing is commonly used for time-series analysis, but it appears to be more useful for determining relationships between high-frequency signals. What am I missing? **Example Data**               df1 <- structure(list(     x = c(315.97, 316.91, 317.64, 318.45, 318.99, 319.62, 320.04, 321.38, 322.16, 323.04, 324.62, 325.68, 326.32, 327.45, 329.68, 330.18, 331.08, 332.05, 333.78, 335.41, 336.78, 338.68, 340.1, 341.44, 343.03, 344.58, 346.04, 347.39, 349.16, 351.56, 353.07, 354.35, 355.57, 356.38, 357.07, 358.82, 360.8, 362.59, 363.71, 366.65, 368.33, 369.52, 371.13, 373.22, 375.77, 377.49, 379.8, 381.9, 383.76, 385.59, 387.38, 389.78),      y = c(0.0192, -0.0748, 0.0459, 0.0324, 0.0234, -0.3019, -0.2328, -0.1455, -0.0984, -0.2144, -0.1301, -0.0606, -0.2004, -0.2411, 0.1414, -0.2861, -0.0585, -0.3563, 0.0864, -0.0531, 0.0404, 0.1376, 0.3219, -0.0043, 0.3318, -0.0469, -0.0293, 0.1188, 0.2504, 0.3737, 0.2484, 0.4909, 0.3983, 0.0914, 0.1794, 0.3451, 0.5944, 0.2226, 0.5222, 0.8181, 0.5535, 0.4732, 0.6645, 0.7716, 0.7514, 0.6639, 0.8704, 0.8102, 0.9005, 0.6849, 0.7256, 0.878),     ti = 1:52),      .Names = c("x", "y", "ti"), class = "data.frame", row.names = 110:161)          ddf<- data.frame(dy = diff(df1$y), dx = diff(df1$x))     ddf2<- data.frame(ddy = diff(ddf$dy), ddx = diff(ddf$dx))     ddf$ti<-1:length(ddf$dx); ddf2$year<-1:length(ddf2$ddx)     summary(lm0<-lm(y~x, data=df1))      #t = 15.0     summary(lm1<-lm(dy~dx, data=ddf))    #t = 2.6     summary(lm2<-lm(ddy~ddx, data=ddf2)) #t = 2.6


title: An unclear regularity condition for the multivariate case of the Cramér-Rao inequality

text: I'm trying to make sense of the multivariate case of the Cramér-Rao inequality. One of its conditions of regularity states: > $$ \frac{\partial}{\partial\theta}\left[\int > T(x)f\left(x;\theta\right)\space dx\right]=\int > T(x)\left[\frac{\partial}{\partial\theta}f\left(x;\theta\right)\right]\space > dx\hspace{10mm}(*) $$ here $\theta\in\mathbb{R}^d$ and $T(x)=\left(T_1(x),\dots,T_n(x)\right)\in\mathbb{R}^n$. I don't understand equality $(*)$:   1. What kind of multiplication is performed inside the integrals?    2. What does it mean to integrate a function that does not take its values in $\mathbb{R}$, as do, presumably, the integrands on both sides?    3. What does the functional $\frac{\partial}{\partial\theta}$ mean on either side of the equality?  Any help will be appreciated.


title: Cointegration structure

text: I have two time series that I am investigating, acc and amb, the time frequency is daily data. They are both non stationary, as evidenced by the follows:               adf.test(df$acc)              Augmented Dickey-Fuller Test          data:  df$acc     Dickey-Fuller = -2.7741, Lag order = 5, p-value = 0.2519     alternative hypothesis: stationary          > adf.test(df$amb)              Augmented Dickey-Fuller Test          data:  df$amb     Dickey-Fuller = -1.9339, Lag order = 5, p-value = 0.6038     alternative hypothesis: stationary      I am looking to test for cointegration between the two time series but the problem I am running into is that the cointegrating vector seems to change in time. 1) First 200 points **Johansen-Procedure** Test type: maximal eigenvalue statistic (lambda max) , with linear trend               Eigenvalues (lambda):       [1] 0.0501585398 0.0003129906      Values of teststatistic and critical values of test:                         test 10pct  5pct  1pct     r <= 1 |  0.06  6.50  8.18 11.65     r = 0  | 10.19 12.91 14.90 19.19      Eigenvectors, normalised to first column: (These are the cointegration relations)                       acc.l2    amb.l2     acc.l2  1.0000000  1.000000     amb.l2 -0.9610573 -2.237141      Weights W: (This is the loading matrix)                      acc.l2       amb.l2     acc.d -0.03332428 -0.002576070     amb.d  0.03986111 -0.001591227      2) First 1000 points **Johansen-Procedure** Test type: maximal eigenvalue statistic (lambda max) , with linear trend               Eigenvalues (lambda):     [1] 0.019211132 0.001959403      Values of teststatistic and critical values of test:                         test 10pct  5pct  1pct     r <= 1 |  1.96  6.50  8.18 11.65     r = 0  | 19.36 12.91 14.90 19.19      Eigenvectors, normalised to first column: (These are the cointegration relations)                       acc.l2   amb.l2     acc.l2  1.0000000  1.00000     amb.l2 -0.8611314 15.76683      Weights W: (This is the loading matrix)                       acc.l2        amb.l2     acc.d -0.008993595 -0.0002419353     amb.d  0.027935684 -0.0002067523      3) Whole History **Johansen-Procedure** Test type: maximal eigenvalue statistic (lambda max) , with linear trend               Eigenvalues (lambda):     [1] 0.0144066813 0.0008146258      Values of teststatistic and critical values of test:                         test 10pct  5pct  1pct     r <= 1 |  1.16  6.50  8.18 11.65     r = 0  | 20.64 12.91 14.90 19.19      Eigenvectors, normalised to first column: (These are the cointegration relations)                       acc.l2    amb.l2     acc.l2  1.0000000   1.00000     amb.l2 -0.8051537 -25.42806      Weights W: (This is the loading matrix)                      acc.l2       amb.l2     acc.d -0.01003068 7.009487e-05     amb.d  0.02128464 6.980209e-05      You can see the marginal change the coefficient values, from -0.96 to -0.86 to -0.80. My question is how to interpret this, what is the optimal look back period, what is the true relationship I should use for future prediction?


title: Why does k-NN perform better than SVR and linear regression?

text: I have a data set used in a regression with 30 attributes and 30K instances. I am trying out a bunch of algorithms (SMO regression, Linear Regression and K-NN) but it was quite surprising to see that K-NN performs much better than SMO and Linear Regression. Since K-NN is a very naive algorithm compared with at least SMO, I was expecting to get better results for SMO algorithm. All results are with 10-fold cross validation:               Algorithm | Correlation coef. | Root mean square error     SMO       | 0.6861            | 2225232        K-NN      | 0.9133            | 1192372     Linear    | 0.7621            | 1853795      What can be the reason for this? Can there be a problem in my dataset or K-NN just happened to be the best algorithm that can explain my data? On a side note: when I depict the actual SMO regression model, I can see that some attributes are rated negatively. Does it shows that some of my attributes are not independent from each other? So that SMO build a bad model because of this? Here is my model for SMO:               weights (not support vectors):      +       0.3724 * (normalized) a1      +       0.0062 * (normalized) a2      -       1.6342 * (normalized) a3      +       0.0909 * (normalized) a4      -       0.2963 * (normalized) a5      +       0.8106 * (normalized) a6      -       0.1123 * (normalized) a7      -       0.4488 * (normalized) a8      +       0.0406 * (normalized) a9      -       1.2222 * (normalized) a10      +       1.7119 * (normalized) a11      -       0.1321 * (normalized) a12      +       0.0292 * (normalized) a13      +       0.54   * (normalized) a14      -       0.3266 * (normalized) a15      -       0.3521 * (normalized) a16      -       0.1737 * (normalized) a17      +       0.0829 * (normalized) a18      -       0.1706 * (normalized) a19      -       0.0909 * (normalized) a20      -       0.4802 * (normalized) a21      +       0.0077 * (normalized) a22      +       0.029  * (normalized) a23      -       0.03   * (normalized) a24      -       0.1366 * (normalized) a25      +       0.1647 * (normalized) a26      -       0.0312 * (normalized) a27      +       0.1972 * (normalized) a28      -       0.104  * (normalized) a29      +       0      * (normalized) a30      -       0.0575 * (normalized) a31      +       0.0215      A second sidenote: Also the distribution of my attributes are pretty skewed. Can this also be the reason why K-NN beats others? (because of the majority vote problem in K-NN algorithm?)


title: Interpretation of interaction effect in negative binomial regression

text: I am trying to interpret my interaction effects, which are all negative. One example: Experience (variable A) x absolute size of the acquired knowledge base (variable B): B= -0.002, exp(B)= 0.998. Can I interpret this interaction in the following way: As variable B decreases, the effect of A increases. How can this interaction be interpreted in terms of exp(B)?


title: Using econometrics, how do I solve out the endogeneity problem?

text: I had a project from before I graduated from college that was effectively me acting as a research aid to a politician. The project served as my senior thesis, as well. The portion of the report that I was working on focused on determining the effect of new fuel economy legislation for heavy-duty vehicles, such as 18-wheelers or large lorries and their drivers. I managed to get the data cleaned up and put together, but when I ran a regression (with dummy variables for year, region, etc) I realized that I was seeing the police-and-crime problem in action: wages affected employment which affected wages etc. Indeed, I checked the regression by running it once with wages and employment against one another and had a _very_ significant p-value spit back out at me, with a 99% R-value -- making me realize that I had a problem. I know that I need to find some variable that only affects one (wages or employment) but not the other. (In the Levitt paper, he uses elections, which affect police presence and not crime.) What type of variable would this be in my wage/employment model? Are there other ways of working around the endogeneity problem? I can post data if it is requested, though my question is more conceptual than data-focused.


title: Do variations of ELO system exist for non-mirror games?

text: The ELO rating system and later improvements like Glicko are specifically designed, as far as I know, for mirror games. That is, we assume that no player has a systematic advantage over another. In chess, this is typically the case because a player will be black/white a roughly equal number of times. But you can imagine that if in some variant of chess, the player who was white/black was not randomly determined that some players could gain a systematic advantage and this would be reflected in their ELO rating. What about games like for example Starcraft, a computer game with 3 races (Terran, Zerg, Protoss)? The races are unlikely to be perfectly balanced. One is probably better than the others. Would it be possible to calculate an ELO rating that took into account racial advantages and thus more closely measured the 'true skill' of a player? How would you do so and what is the best way? (References to existing, well- documented CS/Stat/Math/etc techniques are preferred in the response, but if you have your own unique way of tackling this problem, I'd like to know too.) Note that you cannot assume that the player pools have the same 'true skill' distribution in real life. For example, in Starcraft 2, there are probably more pro-level players in the Terran pool than the rest, because many of the professionals grew up watching Starcraft 1, in which the most admired players were mostly Terran.


title: Sample size for unequal groups in logistic regression

text: I am conducting a study, in which I will be using Logistic Regression for analysis. I will be collecting data in an education setting (about cancer prevention), then I will follow subjects to collect data on who took action after their education (cancer prevention screening) among those who were educated. The two groups (DV) will be "screened" and "not screened." Among the surveys that are administered, I am sure the number of subjects in a group that were screened will be much higher than those who did not. Without unequal groups considerations, my sample size, using Green's formula, i.e. $N\geq 104+m$, is 114. Do the number of subjects in each group have to be equal? How much of a difference is considered not equal? In the case when groups will be unequal, how should it be dealt with for sample size calculation or analysis?


title: Relation between autocorrelation function and periodogram in time series analysis

text: I was wondering if anyone could give me some insight on the relation between the ACF and the periodogram of a time series. I have a bunch of timeseries and their ACF's and periodograms are typically much like the examples below. For my analysis, I'm mostly interested in periodicity at lag 8 and 16 (for theoretical reasons) The frequencies 'B' and 'HB' correspond to lag 16 and lag 8 respectively. The time series actually concerns interresponse intervals in musical performance of a piece that consists solely of eighth notes (16 of them in a 4:4 bar so 'B' stand for bar and 'HB' for half bar). The thing I actually wanted to ask: in my periodograms, I consistently get very large peaks at frequency 0.25 (which corresponds to lag 4). However, the ACF peak at lag 4 is much smaller than those at lag 8 or 16. I was wondering how to interpret this finding. A lot of time series variance can be explained at this frequency even though the lag 4 autocorrelation is quite low? I hope I was sufficiently clear in my question. If not, don't hesitate to ask me. ![enter image description here](http://i.stack.imgur.com/K6bsc.jpg) ![enter image description here](http://i.stack.imgur.com/RKYMl.jpg)


title: How can I get the prior of a random variable that's a function of a random variable in Bayesian data analysis?

text: I have a model which includes the following priors: $\lambda_C \rightarrow \dfrac{1}{\sigma_C^2}$ and $\sigma \sim \text{uniform}(0,500)$ Where $\sigma$ is the standard deviation and $\lambda_C$ is the precision of a normal distribution. Being a beginner, I have problems conceptualizing the prior of $\lambda_C$. How would I even calculate the single values for the prior of $\lambda_C$?


title: When doing systematic sampling, what should be done if the sampling interval (i.e. the skip) is not an integer?

text: Let: population size $=N$; sample size $=n$; sampling interval $=\frac{N}{n} = k$, which can be non-integer; and $r=$ random starting point, which can be non-integer, $0 < r < k$. http://en.wikipedia.org/wiki/Systematic_sampling says we round up $r + mk$ (where $m$ is an integer between $0$ and $n-1$, both inclusive) although the values given (11, 26, 41, 56, 71, 86, 101, and 116) show some rounded-down values. ocw.jhsph.edu/courses/statmethodsforsamplesurveys/PDFs/Lecture2.pdf says: 1) "if $k=5$ is considered, stop the selection of samples when $n=175$ achieved." But this means the last few members would not have any chance of being picked. 2) "if $k=6$ is considered, treat the sampling frame as a circular list and continue the selection of samples from the beginning of the list after exhausting the list during the first cycle." This doesn't give equal chance to each member of being picked, does it? 3) "An alternative procedure is to keep $k$ non-integer and continue the sample selection as follows: Let us consider, $k=5.71$, and $r=4$. So, the first sample is 4th in the list. The second $=(4+5.71) =9.71$ ~ 9th in the list, the third $=(4+2\times5.71) =15.42$ ~ 15th in the list, and so on. (The last sample is: $4+5.71\times(175-1) = 997.54$ ~ 997th in the list)." This uses rounding down of $r + mk$ (different from the rounding up suggested by the Wikipedia page mentioned above). Shouldn't we be rounding off instead to give equal chance to each member of being picked? An even better way is to let random starting point be $R$, randomly selected from the integers 1 to $N$, both inclusive, and use $r + mk$, rounded off and modulo $N$?


title: Leave-out data approach to get intervals for predictive performance of a regression model

text: When measuring predictive performance of a regression model, I am thinking about using repeated data splitting (or leave-out data at random) instead of using bootstrapping. By "repeated data splitting" (not sure if that's even term) I understand:               FOR N=10000     Sample 70% of the data points and fit the model     Predict for the remaining 30%     Store performance     END FOR     Calculate quantiles from stored performance data      I am currently using a GLM in order to rank my data according to a ratio scaled response variable Y (which my model is supposed to predict). Later I am evaluating the quality of the ranking using a specific metric. I would like to calculate the prediction (or at least) confidence intervals for this metric -- that's the thing I care about (and not, say, the variance of the $\beta_3$ estimate) **Would the repeated data-splitting provide useful intervals for the performance of my ranking? i.e. Can I consider them as prediction intervals?** I would argue that this approach makes more sense than non-parametric bootstrap, as long as the data set is large enough, since we are using only actual data (that really occurred at one point in history, as opposed to a bootstrapped point the quality of which completely depends on the quality of the empirical approximation of the underlying distribution). Since my data set is fairly large ($>20,000$) but also extremely skewed and long tailed, I would prefer using the data-splitting but I am not sure if this is a plain heuristic, or whether the intervals are actual prediction intervals.


title: How to calculate spatial correlation between two variables?

text: I have a dataset of point coordinates of individuals and different variables of these individuals. I want to calculate if the spatial distribution of a certain variable is correlated to the spatial distribution of a different variable. For example: do patches of large plants have higher predation rates, even small plants within a patch of large plants. Is it possible to calculate a correlation or a level of significance for these two variables and what software should I use?


title: Updating classification probability in logistic regression through time

text: I am building a predictive model that forecasts a student's probability of success at the end of a term. I’m specifically interested in whether the student succeeds or fails, where success is usually defined as completing the course and achieving 70% or more points out of the total points possible. When I deploy the model, the estimation of success probability needs to be updated through time as more information becomes available -- ideally immediately after something occurs, like when a student submits an assignment or gets a grade on one. This updating sounds sort of Bayesian to me, but given my training in educational statistics, that is a little outside my comfort zone. I have so far been using logistic regression (actually lasso) with a historical data set containing week-based snapshots. This data set has correlated observations, since each student has $TermLength/7$ observations; observations for one student are correlated. I am not specifically modeling the correlation within a particular student’s weekly observations. I believe that I would only need to consider that in an inferential setting since standard errors would be too small. I think--but not sure on this--that the only problem arising from the correlated observations is that I need to be careful when I cross-validate to keep clustered observations in one subset of the data, so that I don’t get artificially low out-of-sample error rates based on making predictions about a person the model has already seen. I am using R’s glmnet package to do a lasso with a logistic model to generate a probability of success/failure and to automatically pick predictors for a particular course. I have been using the week variable as a factor, interacted with all other predictors. I don't think this differs in general from just estimating individual week-based models except that it gives some idea of how there may be some common model that holds throughout the term that is adjusted via various risk adjustment factors at different weeks. My **main question** is this: Is there a better way of updating classification probabilities over time rather than just dividing up the data set into weekly (or other interval-based) snapshots, introducing a time-period factor variable interacted with every other feature, and using cumulative features (cumulative points earned, cumulative days in class, etc)? My **second question** is: am I missing something critical here about predictive modeling with correlated observations? My **third question** is: how can I generalize this to a real-time updating, given I'm doing weekly snapshots? I'm planning on just plugging in variables for the current weekly interval, but this seems kludgey to me. FYI, I'm trained in applied educational stats but do have a background in mathematical stats from a long time ago. I can do something more sophisticated if it makes sense but I need it explained in relatively accessible terms.


title: Fixed-effects or mixed-effects models to population based data

text: There are various definitions of fixed and random effects in linear mixed- effects models in the statistics literature, and even in the fixed-effects models from econometrics. The data I have now represent the entire consumption of some products in a country. I aggregated them to annual level for 6 years by each state. A longitudinal analysis was planned. The primary question here is for the population data by every state in the country, whether fixed-effects models (N states fixed-parameters) or mixed- effects models (1 fixed-effect parameter, N states random-effects) should be used? Personally, I prefer the definition by Fitzmaurice et al. (Applied Longitudinal Analysis), where a fixed-effect represents the population- average, while random-effects represent subject/state-specific effects. Therefore even the data is based on a fixed number of states in a country, each state will have its own growing trajectories over time, which is regarded as random-effects. Some authors argue that for individual patient or small unit data, where an individual or unit was drawn from a large population, mixed-effects models or random-effects models are appropriate. But state-level data for a country can hardly be considered as a subset, then fixed-effects models is more appropriated. Please correct me if any of the statements above is wrong. I wish to hear more opinions on this.