title: Tricky question about MLE

text: > $X_1, \ldots, X_n$ iid ~ Pois($\lambda$). Suppose, you don't know the value > of each $X_i$, but you know if $X_i = 0$ or not for every i. > > Find MLE for $\lambda$. Does MLE always exist? I know how to find MLE. My problem is that I don't understand how to deal with the limited knowledge of the sample. I used to think that we don't care about what sample we have when creating an (not only ML, but any) estimator. So, my answer is MLE would be the same. Am I right?


title: Does a big difference in sample sizes matter for an independent t-test?

text: There is a very confusing question in my mind. I have data, and would like to compare numeric scores between men and women. There is a big difference in those two groups: the number of men is 34, while the number of women is 310, and variances are not equal. As far as I know, when the variances are not equal, I can use the Welch–Satterthwaite equation (independent t-test assuming unequal variances). My question is: can I still use this equation although there is a really big difference in sample size between my two samples? Or is there a certain limit for the difference in sample size between two samples?


title: package R.basic

text: I would like to use the fuction **weighted.median()** in package **R.basic** (and some other ones too, although this one is most important).               install.packages("R.basic")       gives me:               package ‘R.basic’ is not available      and               library(R.basic)       says               Error in library("R.basic") : there is no package called 'R.basic'      There is a downloadable version here. But, typing:               install.packages("~/Documents/R.classes_0.62.tar.gz", CRAN=NULL)      Gives               Error in do_install(pkg) :        this seems to be a bundle -- and they are defunct     Warning message:     In install.packages("~/Documents/R.classes_0.62.tar.gz",  :       installation of package '~/Documents/R.classes_0.62.tar.gz' had non-zero exit status      How to get thru ?


title: Prediction Intervals for General Linear Model

text: How do I derive prediction intervals for a general linear model? My general linear model written in matrix form is, $$ \mathbf{Y} = \mathbf{X} \mathbf{B} + \mathbf{R}$$ with each of the rows of the residual matrix $\mathbf{R}$, distributed as $\mathbf{r}_i \sim \mathcal{N}(\mathbf{0},\Sigma_{R})$. Alternatively, $row(\mathbf{R}) \sim \mathcal{N}(\mathbf{0}, \mathbf{I} \otimes \Sigma_{R})$. I have a least squares estimator as $$\hat{\mathbf{B}} = (\mathbf{X}^{T} \mathbf{X})^{-1} \mathbf{X}^{T} \,, $$ from affine transformations on $\mathbf{R}$, this estimator is distributed as, $$ row(\hat{\mathbf{B}}) \sim \mathcal{N}(\mathbf{B}, \Sigma_{\hat{\mathbf{B}}}) $$ What I wish to do is form a prediction interval on a new instance of $\mathbf{Y}$ that is within the range (an interpolation, not extrapolation). This new instance is $\tilde{\mathbf{y}}_{i}$. What I wish to do is find a prediction interval (or confidence interval?) on this new instance. On my own, I have used affine transformations again on $\hat{\mathbf{B}}$ to find that $$\tilde{\mathbf{y}}_{i} \sim \mathcal{N}(\tilde{\mathbf{x}} \hat{\mathbf{B}}, \mathbf{X}^{T} (\mathbf{X}^{T} \mathbf{X})^{-1} \mathbf{X} \otimes \hat{\Sigma}_{R} ) $$ However I am uncertain in my solution. My textbooks on regular multiple linear regression have formulas for prediction intervals. But if the LS estimator has a distribution, why not just use this result? From it I can make $1\sigma$ or $3 \sigma$ intervals etc. In short, $\textbf{How does one derive prediction intervals for the general linear model? }$


title: Is it a problem when scale normalization yields lot of zeros?

text: I use scale normalization for the features of my neural network. As many of the values occur very infrequently, i.e. the minimum count within all samples is often also the count of the feature within many samples, I end up with a lot of zeros in my sample. E.g., out of 15 features 8 to 10 might be zero. My question is: Is this is still ok, or should I add some smoothing parameter, in order to let the infrequent items also contribute?


title: How to explore/model this problem?

text: I have data extracted from a state machine. Each row of data looks like this:               label1,label2,label3,impedance,did_state_change      The column data types are described as follows:               label1,label2,label3: These are nominal (categorical) values     impedance: This is a ratio value (float number)     did_state_change: This is a categorical value (boolean)      What I'm trying to do is to come up with a model that will help to **estimate the probability of a state change, given a 4-tuple of (label1,label2,label3,impedance)**. I am not sure which statistical model to use and will appreciate some help in choosing the correct analysis/model.


title: Real values in Multinomial Naive Bayes classifier

text: > **Possible Duplicate:**   >  Naive Bayes for two continuous features I've seen that Multinomial Naive Bayes classifier in Sklearn accepts real valued feature vectors. I previously assumed that Multinomial Naive Bayes classifier only handles discrete values If I were to design a multinomial Naive Bayes classifier of my own, how would I handle real valued feature vectors?


title: How to report accuracy of predicted classifications with CI when using k-fold CV?

text: I modelled a binary response against several explanatory variables using logistic regression. There were various combinations of the explanatory variables that produced reasonable fits. I would like to validate and compare the models by using 10-fold CV by comparing the mean positive prediction rates. When doing such a comparison, is it enough to report the mean prediction rate and provide the standard deviation of the accuracy across the ten folds? Or is it possible to report confidence intervals?


title: Separate ANOVAs on subsets of data

text: I often read that, after obtaining complex (e.g., four-way) interactions in a factorial ANOVA, researchers decide to split their data by a factor (e.g., gender) and run separate ANOVAs for these two groups. I realize that this can be very helpful in understanding and explaining complex interactions. However, I have the feeling that this is statistically not completely sound. Isn't it better to perform the "big" ANOVA, including contrasts of interest to help understand the effects? Also, how do you justify the way you split your data? For example, can you only include variables into the separate, "small" ANOVAs that were part of the significant interaction in the "big" ANOVA, or can you chuck in the whole bunch of independent variables again? This seems like effect-fishing to me. I cannot find any helpful literature on this at all so I am unsure whether my gut-instinct is simply wrong. I'd be very grateful for any insight someone may have into the matter!


title: Significant output in Levene's test for equality of variances in MANOVA; what to do?

text: I want to perform a MANOVA in SPSS as follows: i have a independent variable that consists of three groups (- of which the group sizes are not equal), and six dependent variables. When checking the Levene's test for equality of error variances, there are two dependent variables in which the test is significant; so, the error variances are not equal. Because of the significant Levene's test output, my thoughts are that I can't perform this analysis, and have to alter my type of analysis. I know that in a univariate ANOVA, we have the Brown-Forsythe's F and Welch's F which adjust the F-test for inequality of variance. But, in the MANOVA, there seems no such option. My question is what options I _do_ have. Should I instead perform several ANOVA's (and by doing so omitting the power of a MANOVA)? Or should I eliminate the two variables from the MANOVA and run it with only four dependent variables?


title: Cointegrated series

text: My problem is this: I have 2 series $y$ and $x$ that I want to verify are cointegrated. So first I verified them with ADF unit root test and revealed that one variable is integrated at first difference $I(1)$ and the other at second difference $I(2)$ using the test with an intercept! But when I run the ADF test on the variables not including an intercept (the last option – none), both series are integrated at first difference $I(1)$. So after this I want to test for cointegration. Instead of running a Johansen test, which has some difficult options to set up (I set up option 3 and 1 2 for the lagged endogenous variables and found one cointegration equation, but I don't know if it is the correct combination for me; this decision is very subjective, I think), I want to ask you if this is correct: Knowing that $y$ and $x$ are $I(1)$, I estimate the equation $y = a_0 + a_1 x + e$ and save the values for error variable $e$. Then I run ADF test for unit root in residual variable $e$ and find out that the error is $I(0)$ or stationary. Does this not prove that the series $y$ and $x$ are cointegrated?


title: Latent Dirichlet Allocation in PyMC

text: As an exercise to improve my skills in PyMC (Python's Markov chain Monte Carlo library), I am trying to implement Latent Dirichlet Allocation as described here: https://en.wikipedia.org/wiki/Latent_Dirichlet_allocation. The model can be compactly described as $$\boldsymbol\phi_{k=1 \dots K} \sim \operatorname{Dirichlet}_V(\boldsymbol\beta)\\\ \boldsymbol\theta_{d=1 \dots M} \sim \operatorname{Dirichlet}_K(\boldsymbol\alpha)\\\ z_{d=1 \dots M,w=1 \dots N_d} \sim \operatorname{Categorical}_K(\boldsymbol\theta_d) \\\ w_{d=1 \dots M,w=1 \dots N_d} \sim \operatorname{Categorical}_V(\boldsymbol\phi_{z_{dw}}) $$ I came up with the following toy code:               import numpy as np     import pymc as pm          K = 2 # number of topics     V = 4 # number of words     D = 3 # number of documents          data = np.array([[1, 1, 1, 1], [1, 1, 1, 1], [0, 0, 0, 0]])          alpha = np.ones(K)     beta = np.ones(V+1)          theta = pm.Container([pm.Dirichlet("theta_%s" % i, theta=alpha) for i in range(D)])     phi = pm.Container([pm.Dirichlet("phi_%s" % k, theta=beta) for k in range(K)])     Wd = [len(doc) for doc in data]          z = pm.Container([pm.Categorical('z_%i' % d,                                   p = theta[d],                                   size=Wd[d],                                  value=np.random.randint(K, size=Wd[d]),                                  verbose=1)                   for d in range(D)])               w = pm.Container([pm.Categorical("w_%i" % d,                                  p = pm.Lambda('phi_z_%i' % d, lambda z=z, phi=phi: [phi[z[d][i]] for i in range(Wd[d])]),                                  value=data[d],                                   observed=True,                                   verbose=1)                   for d in range(D)])          model = pm.Model([theta, phi, z, w])     mcmc = pm.MCMC(model)     mcmc.sample(100, burn=10)      The tricky part is within $w_{d=1 \dots M,w=1 \dots N_d} \sim \operatorname{Categorical}_V(\boldsymbol\phi_{z_{dw}})$. Given the output of the sampling, I must be doing something wrong, because the model does not converge and I get many warnings about the probabilities in categorical_like that do not sum to one. Is there a PyMC expert around who can shed some light on this all?


title: How to compute the Greenhouse-Geisser Epsilon for the interaction term of the 2 repeated factors in a Repeated Measures ANOVA

text: I’m trying to implement a program in C# to get the Greenhouse-Geisser Epsilon for a 2 factor ANOVA design with 3 levels each factor. Following the procedure on http://www.real-statistics.com/anova-repeated- measures/sphericity/additional-information-sphericity/ I was able to get the GG Epsilon for the single term but for the interaction, I have no idea how to compute the sample covariance matrix. I have tried to partition the data similar to what SPSS requires: Eg. For factors A and B with 3 levels each: A1B1 A1B2 A1B3 A2B1 A2B2 A2B3 A3B1 A3B2 A3B3 4 6 4 5 3 4 9 6 3 5 5 3 3 5 3 6 2 2 6 8 6 4 8 1 6 3 9 7 3 0 4 0 8 3 4 8 4 4 5 3 5 6 2 8 5 3 3 4 2 3 6 2 5 8 Then compute the covariance matrix for each group: 2.166667 0.7 0.333333 0.566667 -0.4 -1.66667 -1.46667 -0.06667 1.366666667 0.7 1.1 2.2 1.1 -0.2 0 -0.2 -0.6 -0.5 0.333333 2.2 7.866667 3.933333 1.8 -1.73333 1.066667 -4.53333 -4.066666667 0.566667 1.1 3.933333 3.766667 4 -1.26667 2.533333 -4.46667 -0.233333333 -0.4 -0.2 1.8 4 7.2 -1 4.8 -5.8 0 -1.66667 0 -1.73333 -1.26667 -1 4.666667 1.066667 2.466667 -0.466666667 -1.46667 -0.2 1.066667 2.533333 4.8 1.066667 4.266667 -3.53333 0.133333333 -0.06667 -0.6 -4.53333 -4.46667 -5.8 2.466667 -3.53333 6.266667 1.133333333 1.366667 -0.5 -4.06667 -0.23333 0 -0.46667 0.133333 1.133333 8.566666667 Then if I perform the procedure as if it was a one within subject factor ANOVA, I then get GGEpsilon = 0.3594. While SPSS Two Way repeated measures gives GGEpsilon = 0.603. This is obviously not the right way to do this. Can anyone point me to the right direction? Or direct me to some reference material where I can find the information required. Thanks much


title: Why does Ridge Regression work well in the presence of multicollinearity?

text: I am learning about ridge regression and know that ridge regression tends to work better in the presence of multicollinearity. I am wondering why this is true? Either an intuitive answer or a mathematical one would be satisfying (both types of answers would be even more satisfying). Also, I know that that $\hat{\beta}$ can always be obtained, but how well does ridge regression work in the presence of exact collinearity (one independent variable is a linear function of another)?


title: Linear regession problem

text: I am doing a simple linear regression on SPSS and have run Spearman and Pearson tests which came out at .40's Now on the linear regression my B and bets have come out negative I am not sure what to infer from this? Amy help would be appreciated, I am only doing a short course in this so my background knowledge is not much


title: How to specify random effects in lme?

text: I have searched for this online for hours but none of online posts is what I am looking for. My question is very easy to implement in SAS Proc mixed procedure but I am not sure how to do it in lme and/or lmer packages. Assume, I have a model, $y = \mu + \alpha + \beta +\alpha\beta + e$, where $\alpha$ is fixed but $\beta$ and $\alpha\beta$ are random. My R code is                f1 = lme(y ~ factor(a), data = mydata,          random = list(factor(b) = ~ 1, factor(a):factor(b) = ~ 1))      Error: unexpected `=` in:                f1 = lme(y ~ factor(a), data = mydata,               random = list(factor(a) =        Could someone please tell me how to specify these random effects in lme? many thanks in advance


title: Sine link with binary regression

text: I have used the SIN link to estimate probabilities, mostly with `Program MARK`. However, I am not sure how the SIN link works. I know the SIN link enables parameter estimation on boundaries (when the probability is 0 or 1), while the logit link seems to have problems in those situations. Given that the SIN link works in `Program MARK`, it seems to me that it should also work with logistic regression in glm in `Program R` since both are estimating proportions of successes out of `n` trials. However, there does not appear to be a built-in SIN link in `R` for glm. I asked a question about this on `Stack Overflow` when I thought it would be a simple programming issue and learned `R` does allow user-specified links in glm. However, nobody posted the `R` code for the SIN link and I have encountered problems with the inverse of the link when trying to do it myself. Namely, there may be multiple solutions given the periodic nature of the SIN function. I imagine the search space must be constrained somehow to avoid multiple solutions. The issue now seems more like a statistical problem than a programming problem and I am posting here hoping someone will either post the `R` code for the SIN link in glm or provide a description of how the SIN link works sufficient for me to codify it myself. The other option I considered was to post a bounty on my `Stack Overflow` question and I might still do that later. Specifically, I am stuck on how to create the inverse of the link and whether the search space must be constrained to avoid multiple solutions resulting from the periodic nature of the SIN function. I am aware of the `family` statement in glm and the `make.link` statement and that user-specified links are possible in `R`, for example with the `logexp` link (`?family` in `R`). Here is a link to my earlier `Stack Overflow` question: http://stackoverflow.com/questions/25431923/using-the-sin-link-with-logistic- regression The SIN link is:               (sin(X * Beta) + 1) / 2      where:               X = design matrix     Beta = parameter estimates      From my earlier post, note that `sin(z)` is bounded between -1 and +1 inclusive. Therefore, `(sin(z) + 1)` is bounded between 0 and 2 inclusive. Therefore, `(sin(z) + 1)/2` is bounded between 0 and 1 inclusive. This inclusivity is probably what enables the link to return estimates of proportions even when probability is 0 or 1. My `Stack Overflow` post also contains `R` code for an 'exhaustive search' further suggesting the link probably does work for estimation of proportions if the search space is somehow constrained. Thank you for any help. Hopefully this is an acceptable use of cross-posting given the statistical focus of the present question and the programming focus of the earlier question. I should point out that in `Program MARK` the SIN link seems to be reserved for models that contain only one entry per row in the design matrix, i.e., without covariates. **EDIT** The inverse of the SIN link, as near as I can tell, is: `lin.pred <- asin((2*y) - 1)` This returns the correct linear predictor as long as the linear predictor, `z`, falls in the range `-1 <= z <= 1`. If covariates are not allowed, as alluded to above, then perhaps this inverse function will always be acceptable. Of course, I am still going to assume for now that dummy variables are allowed in the linear predictor (at least if there is no intercept). The above proposed inverse link might solve the immediate problem of the inverse and constraints on the search space. The next step is to incorporate this inverse link into a user-specified link function in glm in R along with any additional required information to estimate proportions.


title: Interaction terms in regression-theory

text: I am learning categorical variables and interaction terms in the regression course and there few things I don't understand. Suppose a salary data set with predictors   experience $X$ in years   Education coded 1 for diploma , 2 for degree, 3 for higher degree   Management $M$ coded as 1 for a person with management responsibility and 0 otherwise.   Using dummy variable   $E_{i1}=1 $ if ith person has diploma ,0 otherwise.   $E_{i2}=1 $ if ith person has degree ,0 otherwise. Then the regression line would be.   $Y=\beta_0+\beta_1X+\beta_2E_1+\beta_3E_2+\beta_4M+\epsilon$ Then in the book , 6 regression equations are considered as when   * $E=1,M=0 $   * $E=1,M=1$   * $E=2,M=0$   * $E=2,M=1$   * $E=3,M=0$   * $E=3,M=1$  Doesn't these 6 combinations look like the case where there's a interaction between E and M predictors. Why 6 equations?Shouldn't it be 2 equations for regardless of the education level if M=0 and if M=1.   At this stage of the book interaction is not considered.   What is the difference of these 6 combinations and when the interaction term $[E=1*M=0]$ and $[E=2*M=0]$ My second question is,   For this regression line if coefficient of $E1=-3000$ then is it interpreted as   **For fixed level of experience, regardless of the management position, a higher degree is worth 3000 than a diploma**   is this a correct interpretation.


title: Principles of Time Series Clustering

text: I would like to understand complexity of time series clustering. Clustering is similarity based, so as a basic step we evaluate distance between to points in a multidimensional space. In time series, each point is expanded in time - a function of time. Instead of finding distance between two points we are finding distance between two lines? Then a distance component in one dimension is a difference, but for a line that would be something like sum of square differences (like in regression)? Let's consider a simple example of classifying human beings. A person has a number of descriptors at a given fixed moment of time (age, occupation, income, residence, marital status, gender, medical conditions, etc.) We could do clustering on these as dimensions. Now, let's introduce a time component as nearly all of these characteristics change in a lifetime. Usually, we consider a one-dimensional time series. But above is a practical example. Some dimensions will change in time, some will not. They will change differently in different people. For some people some descriptors remain constant, for others - not. There could be also a non-linear higher-order dependencies between the descriptors. So, how do we deal with these complex problems? Clustering is a very computationally intensive task. But time component makes it a lot more complicated. We could even generalize this problem for supervised and unsupervised learning. What if one or more descriptors (attributes, features, independent variables) we have are functions of other variables. How do we model such problems then? Which multi-model solutions exist?


title: How can I constrain the coefficients from linear regression to sum to a given value (in R)?

text: I'm trying to do a linear regression in R, however I have the added constraint the coefficients from the linear regression need to sum to a user given value between 0 and 1. (I understand forcing the coefficients like this will make the fit rather poor, but it's a needed constraint) I've been though the documentation for lm() and the glmc package, as well as some similar questions, but none seem to tackle how to get the coefficients to sum to a specific value. For data T, N, and target sum for coefficients p:               Regression <- function(T, N, p) {              fit <- lm(T[,1] ~ N)          coef <- coef(fit)     ...     }      Ideally I want sum(coef[-1]) = p (ignoring the intercept) Sorry I can't provide more code, but I don't think I can do what I would like with lm(). Does anyone know of a way to do this, or of a package that will let me do this?


title: Methods to analyze a cross-over interaction between a factor and a continuous variable?

text: I'm wondering what would be the best method to analyze a "cross-over" interaction between a factor and a continuous variable. Here's my experimental set-up and hypotheses in a nutshell: 58 participants were randomly assigned to one of two conditions (each n=29) with different manipulations. It was hypothesized that a continuous individual difference variable would make participants more susceptible to the manipulations. The different manipulations were intended to have the opposite effects on the DV. Thus, it was hypothesized that in Condition 1 the DV and the IDV would have a negative correlation, and in Condition 2 they would have a positive correlation. Indeed, correlation analyses (between the DV and the IDV) revealed that r = .-37, p < .05 in Condition 1, and r = .26, p = .18 in Condition 2. I've tried the General and _Generalized_ Linear Models on SPSS to properly test my hypothesis. In the models, I included main effects of Condition and the IDV (which I entered as a covariate), and an interaction effect, but the results don't appear to be very robust (i.e., even the smallest changes in the models seem to have a big effect on the results). I've used dummy coding (0s and 1s) for Condition, and the DV and the IDV appear to be normally distributed. Does this seem like an appropriate way of testing my hypothesis, or can someone suggest a better approach?


title: R algorithm or function to train using a LogLoss error

text: Is there a way to train a RandomForest, GBM or other classification model using a LogLoss error measure? So far I have trained it, but I get an OOB error rate. From there I calculate the sigmoid and use that as my probability in order to minimize the LogLoss... but I have a hard time to find literature or examples on how to optimize an algorithm training using log loss. Here is the LogLoss function in R:               llfun <- function(actual, prediction) {         epsilon <- .000000000000001         yhat <- pmin(pmax(prediction, epsilon), 1-epsilon)         logloss <- -mean(actual*log(yhat)                          + (1-actual)*log(1 - yhat))         return(logloss)     }


title: Clustering longitudinal (trajectory) data

text: I am hoping to implement an unsupervised technique that identifies distinct clusters of individuals based on longitudinal data: 100 continuous or categorical variables measured at different ages. A lot of the functionality provided by R packages seems to have been developed for simpler cases (eg with just one variable measured at different time points), so I was wondering what the best way to approach such a problem with R might be and which techniques (eg Latent Class Modelling) are considered to perform best.


title: Fixed vs Random Effects

text: I have very recently started learning about Generalised Linear Mixed Models and was using R to explore what difference it makes to treat group membership as either fixed or random effect. In particular, I am looking at the example dataset discussed here: http://www.ats.ucla.edu/stat/mult_pkg/glmm.htm http://www.ats.ucla.edu/stat/r/dae/melogit.htm As outlined in this tutorial, the effect of Doctor ID is appreciable and I was expecting the mixed model with a random intercept to give better results. However, comparing AIC values for the two methods suggest that this model is worse:               > require(lme4) ; hdp = read.csv("http://www.ats.ucla.edu/stat/data/hdp.csv")     > hdp$DID = factor(hdp$DID) ; hdp$Married = factor(hdp$Married)     > GLM = glm(remission~Age+Married+IL6+DID,data=hdp,family=binomial);summary(GLM)          Call:     glm(formula = remission ~ Age + Married + IL6 + DID, family = binomial,      data = hdp)          Deviance Residuals:      Min       1Q   Median       3Q      Max       -2.5265  -0.6278  -0.2272   0.5492   2.7329            Coefficients:                   Estimate Std. Error z value Pr(>|z|)         (Intercept) -1.560e+01  1.219e+03  -0.013    0.990         Age         -5.869e-02  5.272e-03 -11.133  < 2e-16 ***     Married1     2.688e-01  6.646e-02   4.044 5.26e-05 ***     IL6         -5.550e-02  1.153e-02  -4.815 1.47e-06 ***     DID2         1.805e+01  1.219e+03   0.015    0.988         DID3         1.932e+01  1.219e+03   0.016    0.987             [...]          DID405       1.566e+01  1.219e+03   0.013    0.990         DID405       1.566e+01  1.219e+03   0.013    0.990         DID406      -2.885e-01  3.929e+03   0.000    1.000         DID407       2.012e+01  1.219e+03   0.017    0.987         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for binomial family taken to be 1)          Null deviance: 10353  on 8524  degrees of freedom     Residual deviance:  6436  on 8115  degrees of freedom     AIC: 7256          Number of Fisher Scoring iterations: 17               > GLMM = glmer(remission~Age+Married+IL6+(1|DID),data=hdp,family=binomial) ; m          Generalized linear mixed model fit by the Laplace approximation      Formula: remission ~ Age + Married + IL6 + (1 | DID)      Data: hdp      AIC  BIC logLik deviance     7743 7778  -3867     7733     Random effects:     Groups Name        Variance Std.Dev.     DID    (Intercept) 3.8401   1.9596       Number of obs: 8525, groups: DID, 407          Fixed effects:                  Estimate Std. Error z value Pr(>|z|)         (Intercept)  1.461438   0.272709   5.359 8.37e-08 ***     Age         -0.055969   0.005038 -11.109  < 2e-16 ***     Married1     0.260065   0.063736   4.080 4.50e-05 ***     IL6         -0.053288   0.011058  -4.819 1.44e-06 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Correlation of Fixed Effects:              (Intr) Age    Marrd1     Age      -0.898                   Married1  0.070 -0.224            IL6      -0.162  0.012 -0.033               > extractAIC(GLM) ; extractAIC(GLMM)          [1]  410.000 7255.962     [1]    5.000 7743.188      Thus, my questions are: (1) Is it appropriate to compare the AIC values provided by the two functions? If so, why does the fixed effect model do better? (2) What is the best way to identify if fixed or random effects are more important (ie to quantify that the variability due to the doctor is more important than patient characteristics?


title: Estimation of mean, variance and mean squared error of a histogram that poorly models real distribution

text: I'm having a basic statistics problem. I'm running an experiment using a photon counting detector to measure a diffraction pattern. I want to calculate the centroid of each diffraction spot and the error in the centroid. In 1D my data can be considered to be a histogram of counts within equally spaced bins. I'm calculating the centroid and error as the intensity weighted mean and standard error of the distribution. My problem occurs when I have a spot whose distribution is narrower than the pixel size (i.e. in 1-D all the counts fall within a single histogram bin) and I get a variance, and therefore an error, of zero in the determination of the mean. My first thought was that since the position of each count has an error of $\pm 0.5$ pixels, I could just propagate that and add it to the standard error; however, the fact that no-one else seems to be doing that worries me somewhat. I know that the $\text{MSE} = \text{variance} + \text{bias}^2$, but I'm equally unclear about how to calculate the bias of my data, although I'm aware that in histogram statistics the bias varies with the bin size. I'm pretty sure I'm missing something simple; if anyone has any ideas that would be great.


title: Measuring the possibility of occurrence of tree according to a transition matrix

text: I have a $n\times n$ transition matrix where each value represents a node (of the tree) to node transition probability. I want to measure of possibility of occurrence of a tree according to the transition matrix. How can I do that?


title: Advice needed on auction system simulation

text: I am trying to simulate an auction system in which a number of competitors, $N$, independently offer a discount from a reference price previously published by the buyer. The order is awarded to the a certain competitor based on a complex formula involving the max discount not exceeding a certain threshold, the median of the bids, and other parameters. I would like to generate a certain set of discounts by drawing them from a positive distribution, centered around a prescribed mean, and always between a given min and max discount.   * The probability of having a value around min or max should be approximately $1/N$.   * The normal distribution is not probably the best solution as it can give negative values. I could artificially “shrink” it between min and max by imposing $\text{std.dev}=(\text{max}-\text{min})/6$, but this would give a very low probability of occurrence for the extremes. **Do you have any suggestion for a distribution with the above properties?** Ideally it should be easy to handle in MS Excel, where the simulation will be done.


title: Explanation of density rewriting?

text: Can somebody please explain the math behind this statement to me? I am not sure how they represent the left hand side by that integral and finally how it is proportional to that.   \begin{align} p(S_{t+1}|z_{t+1}) &=\int p(S_{t+1}|S_t,z_{t+1})d\mathbb{P}(S_t|z^{t+1}) \\\ &\propto \int p(S_{t+1}|S_t,z_{t+1})p(z_{t+1}|S_t)d\mathbb{P}(S_t|z^t) \end{align} where $z^{t+1}=(z_1,...,z_t)$ Also, are these two equivalent? $$\int p(S_{t+1}|S_t,z_{t+1})d\mathbb{P}(S_t|z^{t+1}) \stackrel{?}{=} \int p(S_{t+1}|S_t,z_{t+1})p(S_t|z^{t+1})dS_t$$ as in, is that the correct way to interpret integrating with respect to a probability measure?


title: Checking autocorrelation in non-parametric methods

text: I would like to forecast short-term electric load by using Artificial Neural Network and Support Vector Regression. However, there's one question that sticks in my mind. In such forecasting with non-parametric methods, can one check the autocorrelation of the errors, which is actual values subtracted from predicted values? I have checked papers about load forecasting in the literature, but I have run into one and only paper by "Espinoza et al" that checks the autocorrelation of the NARX errors. ftp://ftp.esat.kuleuven.be/sista/ida/reports/06-84.pdf Other papers do not check the distribution of the errors. I believe it's the correct move because you are already using non-parametric methods. I am, nevertheless, curious about your thoughts on this matter.


title: Drawing from Dirichlet distribution

text: Let's say we have a Dirichlet distribution with $K$-dimensional vector parameter $\vec\alpha = [\alpha_1, \alpha_2,...,\alpha_K]$. My question is, how can I draw a sample (a $K$-dimensional vector) from this distribution? I will be very greateful for possibly simple explantation.


title: Quasi Maximum Likelihood estimator for GARCH with jump (Compoud poisson process with normal distribution)

text: I have to modify some R code to take jumps into account in my modelisation. First I build a simulation of a Garch 1,1: $r_t=\mu+\sqrt{h_t}\epsilon_t$ $h_t=w_0 +\alpha(r_{t-1}-\mu)^2+\beta h_{t-1} $ where $\epsilon$ follows a normal distribution. I then try to estimate the parameters, so I build and optimise the log likelihood that way:                   garch_loglik<-function(para,x,mu){     # Parameters     omega0=para[1]     alpha=para[2]     beta=para[3]     # Volatility and loglik initialisation     loglik=0     h=var(x)     # Start of the loop     vol=c()     for (i in 2:length(x)){     h=omega0+alpha*(x[i-1]-mu)^2+beta*h     loglik=loglik+dnorm(x[i],mu,sqrt(h),log=TRUE)     }      Now I have a modified GARCH: $r_t=\mu+\sqrt{h_t}(\epsilon_t+\sum_{i=0}^{N_{t}} x_{i,t})$ $h_t=w_0 +\alpha(r_{t-1}-\mu)^2+\beta h_{t-1}$ Where epsilon follows a normal distribution, N a poisson(0.1) distribution, $x_{i,t}$ a normal distribution(-0.1,0.3). I have modified the simulation of my GARCH process consequently.               pois=rpois(n,0.1)     eps=rnorm(n,0,1)+ sum(rnorm(pois[i],-0.1,0.3))      But I strugle to modify the loglikelihood. From what I understand I have to modify the line with dnorm in the loglikelihood, but I don't know how to find the distribution of a normal law + a sum of normal depending on a poisson process. Is it possible to compute directly ? (a function to 'mix' distributions) Do I have to build it manually ? (simulating multiple variables to find the distribution.)


title: Easy book to understanding basic concepts

text: I have a medium-strong background on programming and logic, however I'm trying to start using R, and other tools to make machine learning based studies of some problems. I did take probability and statistics at university, but it was a very poor taught one and never realized until now. I found two books on Amazon I liked:   * Naked Statistics: Stripping the Dread from the Data   * The Complete Idiot's Guide to Statistics Both seem like good starting points, at least to quickly understand key concepts and then use different automatic tool techniques for my purposes, but has anyone read those? Which one is better? P.S. I would like to learn: distributions (normal, poisson, etc), descriptive statistics, inferential statistics, hypothesis, $R^2$, probability, Bayes, Markov chains, and all related concepts.


title: Should it be A|B or B|A?

text: I came across a question in condition probability which left me wondering. A departmental store reports that 30% of sales are transacted via cash, 60% are credit card and 10% with a debit card. 20% of the above cash purchases, 90% of the above credit card purchases and 80% of the above debit card purchases are made for purchases more than $200.               So Let $C = {cash}$, D = {credit card}, E = {debit card} and M = {purchases > $200}          P(C) = 0.3     P(D) = 0.6     P(E) = 0.1      What I am not really sure is if it should be $P(M|C) = 0.2$ or is it $P(C|M) = 0.2$. Was wondering if I may get some advice on this please.


title: Comparing effects of slightly different sets of independent variables on multiple dependent variables in logistic regression

text: I have a problem with logistic regression. I want to find out if there is any effect of social capital variables (civic participation, generalized trust...) on different dependent variables (namely self-reported health, depression, cognitive functioning etc.). Similar problem had been answered here: Comparing logistic coefficients on models with different dependent variables?, but my problem is slightly different. I had created range of different models for this dependent variables, but to achieve best fit of my models I want to add for example depression measure to my self-rated health model. I presume that being clinicaly depressed can affect self-rated health. So my questions are: Is effect of social capital variables on this different dependent varibales comparable with slightly different independent variables lists? Is there something I can do to make my models comparable? Thank You very much!


title: coefficients extracted by effects package in R

text: Consider this logistic regression:               mtcars$vs <- as.factor(mtcars$vs)     log_reg_mtcars <- glm(am ~ vs*wt +vs*mpg, family = "binomial", mtcars)      I tried using the effects package to extract some coefficients from the model:               as.data.frame(effect("vs", log_reg_mtcars))           NOTE: vs is not a high-order term in the model       vs       fit       se       lower     upper     1  0 0.3957869 2.916877 0.002150238 0.9950031     2  1 0.0336822 2.181009 0.000484832 0.7146704      Could anyone explain why the coefficients/standard error for `vs` given by this code are different. Or in other words, what is the effects package doing?               as.data.frame(summary(log_reg_mtcars)$coef)[2,1:2]               Estimate Std. Error     vs1 -40.70047   40.53197


title: How to compute the binomial confidence interval for non integer parameters?

text: I am trying to implement C4.5 and I use Quinlan's book called also "C4.5: Programs for Machine Learning". The pruning method proposed by C4.5 is to use a pessimistic CI-based estimate of the error from a node. So, given a confidence level ($\alpha$), the total weight of the cases considered $N$, and the total weight of the mis-classified cases $E$, the proposed estimator is the upper limit of the confidence interval for binomial distribution with $N$ trials and $E$ errors. I fully understand that we are not talking here about a sample, that the estimated value is more or less a heuristic, and so on. I am not interested in this, I try to met all the original specifications. The problem which I do not understand is the following one: due to how the missing values are handled by C4.5, it is usual to have quantities $N$ and $E$ as non-integer values. If this is the case, how can one compute CI for binomial? Note: Now I am currently studying this article: http://www.sigmazone.com/binomial_confidence_interval.htm to give me some light on this topic. Still, I believe that this will not solve my problem.


title: Item Response Theory alternatives

text: What are other approaches than Item Response Theory to model learning of students in standardised tests?


title: How to simulate data with values of Mardia's Kurtosis?

text: I have two (related) questions regarding the simulation of correlated, non- normal data with a user-specified value of Mardia’s 1970 multivariate kurtosis. (1) It is very common in my field to use the method developed by Vale and Maurelli (1983) of doing a polynomial transformation of standard, normal random variables ($Y=a + bX + cX^{2} + dX^{3}$ where $X \sim N(0,1)$) where the coefficients c and d control the (univariate) skewness and kurtosis estimates. The non-normality is then achieved through the 1-dimensional marginal distributions but there is no control over the multivariate skewness/kurtosis estimates (as defined in Mardia (1970). So here comes my first question: How can I obtain the population values of multivariate kurtosis if I have only the population values of all the univariate marginals? Here's my attempt at a solution. It's using the semTools package in R (calculates Mardia's measure) and the lavaan package (implements the Vale & Maurelli method)                   library(semTools)         library(lavaan)               ### three variables all correlated at 0.5     model <- 'x1 ~~ 0.5*x2               x2 ~~ 0.5*x3               x3 ~~ 0.5*x1'               ### vector to store Mardia's kurtosis values     mardia_values <- double(100)          for (i in 1:100){          ### generate simulated data with N=100,000 population skewness of 2 and kurtosis of 7     dataz <- simulateData(model, sample.nobs=100000, skewness=c(2,2,2), kurtosis=c(7,7,7))          mardia_values[i] <- as.numeric(mardiaKurtosis(datz)[1])     }          mean(mardia_values)      So essentially I'm approximating it via a simulation (at a large sample size) Would this method be acceptable? I'm very open to suggestions (particularly if it can be derived analytically (2) Now, on a related note... does anyone know of any method to simulate data where the value of Mardia Kurtosis can be specified by the user? So far, the only thing I have been able to come up with is, once again, using the Vale and Maurelli (1983) method and sort of just go by trail and error until I get the value that I want. Since Mardia defined these measures in the population, I assume it should be possible to have some method to generate data where the value can be specified in advance, but I do not know how. Thank you! References Mardia, K. V. (1970). Measures of multivariate skewness and kurtosis with applications. Biometrika, 57(3), 519-530. Vale, C. D., & Maurelli, V. A. (1983). Simulating multivariate nonnormal distributions. Psychometrika, 48(3), 465-471.


title: Is there enough information provided? If there are serial correlations, build a mean equation for the log returns

text: `If there are serial correlations in the data, build a mean equation for the data.` As shown in the below figure, the PACF and ACF are basically identical. How should I know if I'm supposed to be building an AR(p) model or MA(q) model? Does serial correlation imply that I should be sticking to an AR(q) model? The next problem in the set asks us to try to fit an ARMA model, which would imply that I shouldn't try to build an ARMA model for this question (even though I see that it fits the data much better than an MA or AR model alone). ![enter image description here](http://i.stack.imgur.com/emFll.png)


title: Using ACF or PACF to build an AR(p) model?

text: In my text, we're told to use PACF to find the order of MA(q) models and ACF to find the order of AR(p) models. But in a homework problem, the professor specifies the order of an AR(p) model using the PACF. This is pretty confusing. Can anyone shed light on the apparent contradiction?


title: What statistical analysis should I use to compare shooting percentages in basketball?

text: If I want to compare shooting percentages (basketball) in two conditions, what type of statistical analysis would I use? If I enter data in columns condition 1 makes   condition 1 attempts   condition 2 makes   condition 2 attempts Is that a bad way to do it? Should I assign a 1, 0 to make and miss and compare by groups? Any advice? I want to look at shooting percentages on lay-ups off of one foot vs lay-ups off of two feet in girls varsity basketball.


title: Dynamic decision networks

text: I was wondering if anybody had some good references on dynamic decision networks, is there a nice recent survey on this topic? Thanks in advance.


title: Comparing questionnaires with actual observations: is it worth it?

text: First of 2 part query: I have observed a 100 Drs and recorded the number of times they touch surfaces while in a room: $$ \begin{array}{l|c|c|c|c} Dr\text{ number}&Bed&Table&Chair&Door\\\ Dr 1&3&2&2&1\\\ Dr 2&4&1&0&1\\\ \vdots&\vdots &\vdots &\vdots &\vdots\\\ Dr 100& 2&3&1&1 \end{array} $$ I would like to calculate the probability of him touching the Table for example. Would it be:? $$ \begin{array}{l|c|c|c|c|c} &Bed&Table&Chair&Door&\text{Total}\\\ \text{Average contacts}&3&2&1&1&\text{7}\\\ P(surface)&3/7&2/7&1/7&1/7&1 \end{array} $$ This is confusing me but I think a fresh pair of eyes will spot it instantly. 2nd Part: I gave questionnaires out to some doctors who ticked boxes instead of giving contact numbers ie: $$ \begin{array}{l|c|c|c|c} Dr\text{ number}&Bed&Table&Chair&Door\\\ Dr &\checkmark&\checkmark&\checkmark&\checkmark\\\ Dr&\checkmark&\checkmark&\times&\checkmark\\\ \vdots&\vdots &\vdots &\vdots &\vdots\\\ Dr 100& \checkmark&\checkmark&\checkmark&\checkmark \end{array} $$ What can I deduce from this table of ticks to compare with the table of exact values? Or is it not worth anything? Please ask for any clarification needed, Regards


title: What test should I perform on a structural equation model?

text: ![model](http://i.stack.imgur.com/9wV0U.png) This model is a replication of an earlier study. After factor analysis in SPSS, I need to test my hypothesis (as you can see there are 7) all should give significant positive relation with each other. For example, `Innovativeness` has a significant positive effect on `Ease of use` etc. My supervisor insists that this is mediation analysis and you have to see whether it's partial or full or indirect by AMOS but the research paper from which I've taken it did not mention any of those things. They just used structural equation modeling and mentioned each variable's effect on another one by one. I want to know what test I should perform? is it really mediation or just SEM? Should I use SPSS or AMOS? My Prime objective is to see if these variables affect the dependent variable " **Future shopping intention** ". Especially, how does `Innovativeness` directly effect " **Future shopping intention** " and it's indirect effect via `Ease of use`. But `Usefulness` and the other effects on `Future shopping intention` are also important. I do not have much time left for my thesis and I have already mentioned couple of times that I am quite poor in SPSS and even worst in AMOS. Suggest and help me out with the right test. I need step by step guide. If you want to see what original research paper has used for testing then please check it out on emerald insights 10.1108/14684520810914025


title: SPSS Amos Binary Outcome - Model Fit

text: I was just wondering if anybody could help me, I'm conducting a path analysis in SPSS Amos with a binary outcome. I've specified the variable as binary using Tools>Data Recode, I've then calculated the variable estimates using MCMC. My question is: how do I know calculate overall model fit? The model is as follows, with all variables except psychosis (binary) being continuous.   * Esteem <\--- Bullying   * Intrusion <\--- Bullying   * Intrusion <\--- Esteem   * PsychosisRisk <\--- Intrusion   * PsychosisRisk <\--- Bullying Many thanks.


title: State-space model with strict positive values

text: I've got a random process of positive observations. I've built a state space model (following an ARMA(1,3) structure), found the parameters that fit the process observations through log-likelihood maximization, and after simulating the obtained model I've seen that negative values occur. I really don't have any idea on how to correct these results. Is it that some type of lower bound must be formulated at the stage of parameter estimation?


title: Generalized Binomial Models

text: With the standard Binomial probability distribution we consider n trials each with a probability of success p. This can be somewhat generalized to the Beta- Binomial Distribution which is effectively the same, but this time the probability of success is a random variable with a Beta Distribution. My question is: Could I in theory use a p which samples from some arbitrary probability distribution (continuous or discrete) so long as the support is on [0,1], or is there something special about the Beta distribution which allows for such a distribution?


title: Effect of Violating Spearman_Brown Prophecy Formula Assumptions

text: The Spearman-Brown Prophecy formula is sometimes used to predict the reliability of the average assessment based on N raters using the interrater reliability calculated from the N sets of individual ratings in a given setting. My question concerns application, say in an interview setting where there are many candidates and N assessors. In real life, assessors influence each other. Even if they do not speak between the end of the interview and the rating, the assessors have communicated by the questions they asked during the interview and by their body language (sometimes, even to the point of sighing or rolling their eyes). So, in this setting we calculate interrater reliability based on ratings that are not independent and then apply the SB formula to estimate the reliability of the mean rating. But if the interrater reliability is inflated, perhaps the SB approach does not make sense. Perhaps the mean of the N ratings is not well predicted using the interrater reliabilities that are not based on independent ratings. I think interrater reliability is taken as a (very indirect) indication of validity. (At the least, validity can be no higher than the square root of reliability.) Is this a critique that has been raised before? Any thoughts or can you suggest places to find material on this aspect of the use of the SB approach, specifically on the extent of the danger in a quantitative sense?


title: building a decision tree after association rules

text: I built association rules using R Arules package. Then I filtered those rules and kept rules that have a specific variable on RHS. That variable is my Y variable. In entire dataset, my variable occurs only 30% of times. 99% of rules that I got had a lift <1.5 Does this mean that my X variables don't have a good predictive power? In such case does it make sense to build a decision tree on that data? would it provide us any more information? Or does it mean that there is a class imbalance and I need to somehow increase frequency of Y variable? (over sampling)


title: What is a correct way to test these data for significance?

text: I have data on 1548 disease cases and 1137 controls. Specifically, for each person, I have a certain score $n$. I am interested in the following:   * First, for either group, I want to know if the difference $D$ between the sample variance (I'll call it $SV$ here) of the $n$ variable and a special additive variance of the $n$ variable ($AV$) is significantly less than $0$. In other words:  $$ D = SV - AV \\\ ~ \\\ H_0: D = 0 \\\ H_A: D < 0 $$   * Second, I am curious about the difference between the two groups. If I let $D_1$ be the D property above for cases and $D_2$ be likewise for controls: $$ H_O: D_1 - D_2 = 0 \\\ H_A: D_1 - D_2 \ne 0 $$ I'm not too acquainted with the fine details / considerations of selecting a test for this, but I suppose it would be similar to a test for the difference between 2 means--something one-sided in the first case, and two-sided in the second case. Can someone point me in the right direction?


title: How to handle news coverage ratio in a regression analysis?

text: I have collected data from Twitter, Google and Wikipedia to look for any correlations (and perhaps predictive value) against the stock market. I have collected daily Twitter, Google and Wikipedia data for 148 days. In regression analysis, I am using a dummy variable, "sector". All 100 companies from which I have collected data can be classified in one of 6 industry sectors. Now, I started thinking that a regression analysis may be much stronger for companies that have had higher news coverage. It is acceptable to think that if there's much news about a company, people tend to tweet more about that company. I have collected news headlines mentioning each of the 100 companies over the same 148 days. There are many days for many of the companies on which there is no news at all, so trying to correlate non-news days to number of tweets for instance, doesn't seem logical. I have now aggregated all the headlines, for instance: during the 148 days, there were 846 headlines about Apple, 137 about Microsoft and 0 about Abbott Laboratories. I'm now thinking to of adding a dummy variable to the dataset, something like `dNews10`, `dNews100`, meaning companies with less than 10 news articles, 10 to 100 articles, etc. Would this be an acceptable variable to add to the analysis? Usually, as I have been taught, dummies are properties to the research subject (man: hight, sex, eye-color, etc.). I don't know if news coverage would qualify.


title: Predicting outlier series in Kalman filter

text: I have built a Kalman Filter model for flu forecasting as shown below. Y - Target Variable X1 - Predictor1 X2 - Predictor2 While forecasting into the future, I will NOT have data for all three variables. So, I am predicting X1 and X2 using two Kalman filters. The code is below               x1.model <-  dlmModSeas(52) + dlmModPoly(1, dV=5, dW=10)     x2.model <-  dlmModSeas(52) + dlmModPoly(1, dV=10, dW=10)          x1.filter <- dlmFilter(c(train$x1, rep(NA, noofsteps)), x1.model)     x2.filter <- dlmFilter(c(train$x2, rep(NA, noofsteps)), x2.model)      Now, I am forecasting Y using the predicted X1 and X2 as below               pred  <- cbind(c(train$x1, x1.filter$f[260:312]),c(train$x2, x2.filter$f[260:312]))          Y.model   <- dlmModSeas(52) + dlmModReg(pred, dV=0.5, dW=c(10,0.05,0.1))     Y.filter  <- dlmFilter(c(train$Y, rep(NA, noofsteps)), Y.model)        I have two questions 1) The flu season usually starts in early September and it usually peaks during the last week of December or the first week of January. The following are the stages               Stage1 - September, October     Stage2 - November, December     Stage3 - January, February     Stage4 - March, April      The dlmModSeas captures this seasonality very well. But what do I do to forecast an off season wave that has the same stages as above. Say it starts in May and it peaks at end of August or early September. I am able to get the information of this off season wave from the user during the start, i.e. in May. How do I introduce this change in seasonality? 2) Is there a better way to do this forecasting? Please Advise


title: Overall significance test for the effect of an independent continuous variable on a categorical dependent variable

text: If I have a dependent variable having more than two categories (the categories can be ordered) and a few independent variables which are all continuous, then how can I see whether the independent variable affects the dependent variable significantly or not? I can see that if I run an ordinal logistic regression, then compared to the reference category, a beta coefficient and its significance level is given for each of the other categories of the dependent variable. But my purpose is to comment in an overall manner whether the independent variable affects the dependent variable. For example, whether income(continuous) will affect social status(categorical) significantly. I don't want to find out in comparison to any reference category, I just want to make an overall comment. If this is possible, then how can I find that overall significance?


title: How could I extract the distribution of one RV when given a set of sums of two RVs?

text: I'm sitting on a rather complicated task and I'm not sure which method to use here. For a university class, we are consulting a psychologist who has shown people images like this: ![enter image description here](http://i.stack.imgur.com/fsUne.png) The task was to   * a) Find which of the diamonds has a different color than the others. This takes an amount of time called "selection time"   * b) Find which side of the diamond is chipped off. This takes an amount of time called "discrimination time"   * c) Report the side of the missing edge by clicking the left or right mouse button. Now, we were only able to measure the total reaction time, i.e., selection _plus_ discrimination time, because we don't know how much of a person's total reaction time is spent searching or, respectively, discriminating. But we are actually interested in _only_ the selection time, part `a)` of the above list. This is why we additionally performed the following experiment: Before showing people the picture, we let a little square (called a 'cue') flash up in the spot where the target diamond will appear. The person will know at what location the target diamond appears. By this we hoped to eliminate the selection time, thus only measuring part `b)`, the discrimination time. But note that we are still interested in _selection_ time, not _discrimination_ time. So our data set consists of 16 persons, with reaction times of around 360 shown pictures per person. About half of them were presented without the little flashing 'cue' square (which means the total reaction time includes the selection time), and the other half had the flashing square (measuring only discrimination time). Our primary interest is the _distribution_ of the selection times. We heard a few keywords for which method we could use to extract that selection time. One was "deconvolution of distributions", then we heard about a "fast Fourier transform", which we're not sure how to use yet. Simple subtraction of reaction times without cue from the reaction times with cue might work, but we don't know how to pair those reaction times. I know this sounds confusing, but I hope I could make our task somewhat clear. If something is not explained properly, I can of course clarify. If you had any suggestions on which method might be appropriate here, we would be most grateful.


title: Question about correct ANOVA design

text: We have carried out an experiment that measures how many references it takes within a group to find the most expert person in that group on a specific topic (DV). We have 7 topics (4 of which are considered "visible" and 3 of which are considered "invisible"). The visible/invisible distinction is our treatment. We have 42 groups (school classes). We conducted the experiment with all 7 topics in each group so that groups and topics are crossed. We believe that the number of references required to find the most expert person in a class depends on the visibility/invisibility of the search topic. I am struggling to come up with the proper design. I believe that search topics are nested in the treatment (visible/invisible). I would also like to consider the fact that I have 7 measures for each class. I am not interested in assessing a direct effect of class but just want to control for the fact that the 7 measures within one class may not be independent from each other. Can this still be done with an ANOVA (and if so what type?) or should I rather use a GLM?


title: Beginner question: available space in a ware house

text: I have a year's worth of data for the available slots in a warehouse. There is a data point for every minute of the year.   * timestamp   * available slots The data seems to follow a certain pattern, namely, less space during the working hours, more space during the night and on the weekends. My question is, can probability theory help me make predictions here? How can I calculate the likelihood of the warehouse being empty, half full, full during a certain time_span / weekday combination. For instance,   * how likely is it, that the warehouse will be mostly full between 10am and 11am on a regular Monday?   * how likely is it, that the warehouse will be empty on 2pm on a regular Thursday


title: Interpreting the regression output from a mixed model when interactions between categorical variables are included

text: I have a question about my use of a mixed model/lmer. The basic model is this:               lmer(DV ~ group * condition + (1|pptid), data= df)      Group and condition are both factors: group has two levels (groupA, groupB) and condition has three levels (condition1, condition2, condition3). It's data from human subjects, so pptid is a random effect for each person. The model found the following with p value output:                                                  Estimate MCMCmean HPD95lower HPD95upper  pMCMC Pr(>|t|)     (Intercept)                          6.1372   6.1367     6.0418     6.2299 0.0005   0.0000     groupB                              -0.0614  -0.0602    -0.1941     0.0706 0.3820   0.3880     condition2                           0.1150   0.1151     0.0800     0.1497 0.0005   0.0000     condition3                           0.1000   0.1004     0.0633     0.1337 0.0005   0.0000     groupB:condition2                   -0.1055  -0.1058    -0.1583    -0.0610 0.0005   0.0000     groupB:condition3                   -0.0609  -0.0612    -0.1134    -0.0150 0.0170   0.0148      Now, I know that the rows listed compare each level of the factors to the reference level. For group, the reference is groupA and for condition, the reference is condition1. Would I be correct in interpreting this output in the following way:   * No overall differences between the groups (hence groupB having a p of >.05)   * Overall differences between condition 1 and condition 2, and between condition 1 and condition 3.   * Differences between groupA, condition 1 versus groupB, condition 2 and also between groupA, condition 1 versus group B, condition 3. Is that correct? I think I'm a little confused about how to interpret this with regards to interactions between levels of two different factors. I've read various questions on here and done some web searches, and managed to get contrasts set up with glht: would that be a better way to look at the differences between the groups and conditions? I figured that would be the case given the signs of interactions being present here.


title: Survival analysis/LTV analysis with few deaths

text: I'm tasked with performing a lifetime value analysis for the clients of a small to medium sized business, however I have an issue that client turnover is actually quite slow and there are not many 'deaths' so far to beging the process of estimating mean survival time. Are there any techniques I can apply or should be aware of for estimating the lifetime value of these clients, givent he relatively young age of the company and the low level of client drop off. (The company has approximately 450 clients and probably fewer than 50 drop offs so far, for reference).


title: Are nominal attributes strict classifications and equivalent to enumerations in programming languages?

text: I've been looking random forest algorithms for text classification and referencing the Mahout random forest decision tree description. In it, there is a reference to two types of variables, nominal attributes, and real-valued attributes. However, I'm not particularly sure what a nominal attribute is. Real-valued attributes seem to be just numeric attributes, things you can measure with a real numeric value. Looking around, I've found the "Level of Measurement" entry on Wikipedia, specifically the section on "Nominal Scale". It states: > At the nominal scale, i.e., for a nominal category, one uses labels; for > example, rocks can be generally categorized as igneous, sedimentary and > metamorphic. For this scale, some valid operations are equivalence and set > membership. Nominal measures offer names or labels for certain > characteristics. > > Variables assessed on a nominal scale are called categorical variables; see > also categorical data. Later on, it states (and I think this is the most crucial part): > The central tendency of a nominal attribute is given by its mode; neither > the mean nor the median can be defined. This sounds a great deal to me like enumerations in programming languages (or any key-value pair where the key is a distinct numeric value and the values are categories/labels); I would rarely perform mean, median, or any other arithmetical operations on the numeric values, but I'd certainly be able to determine the mode (how many things that classification is applied to a certain set). Does this sound about right?


title: Significance testing of three or more correlations using Fisher's transformation

text: Following on from my earlier posts, as far as I can understand, if I have three correlation coefficients, I will have to test them in pairs to see if there is a significant difference among them. This means that I would have to use Fishers transformation to work out the z score of r and then the p value of z (which the recommended calculators in the earlier posts do, thankfully) and then ascertain if the p value is higher or lower than my alpha value (0.05) for each pair. i.e. If 21 to 30 years olds is Age Group 1, 31 to 40 years is Age Group 2, and 41 to 50 years is Age Group 2, my comparison of the correlations between their shopping habits and weight loss would be:   * Group 1 vs. Group 2   * Group 1 vs. Group 3   * Group 2 vs. Group 3 Instead of doing three separate calculations, is there a way to do all these calculations in a single step?


title: Regression analysis for multiple related outcomes

text: This is my first post in CrossValidated :) I have the dependent variables: 1\. Coughing when being fed liquid (Yes/No) 2\. Coughing when being fed solid (Yes/No) and an independent variable: 1\. Age group I am thinking of using logistic regression method to answer the question "Is coughing related to age?". However, I am not sure how to accommodate both dependent variables in the model. I note the discussion here about multiple related outcome analysis but would like ask for your advice if it's relevant to my analysis. Thanks in advance.


title: Precision and recall for clustering?

text: Im confused about calculating precision and recall for clustering mentioned in this paper (last paragraph of column 1 on page 5). Suppose, if we are given with resultant cluster indicator matrix $C$ of size $n \times k$ where $n$ is the number of data points and $k$ is the number of cluster. A data point may be assigned to multiple cluster. Can any one provide me the exact details of calculating precision and recall for this resultant cluster indicator matrix $C$ as mentioned in the paper. Assume the true cluster indicator matrix is also available, say $True\\_C$


title: Profile variable in collaborative filtering

text: I'm trying to create a recommendation system based on purchases. I did some tests and I found that for some groups of customers, the recommender works very well, but not for others. How can I implement a clustering approach or a profile variable in a item-based collaborative filtering algorithm?


title: What is the best method for checking convergence in MCMC?

text: What is your preferred method of checking for convergence when using Markov chain Monte Carlo for Bayesian inference, and why?


title: Creating a translated sample by subtracting mode

text: > Consider a sample $x_1,...,x_n$ which we will call the original sample. To > create a translated sample $x_1+c,...,x_n+c$, we add a constant $c$ to each > sample point. > > Let $y_i=x_i+c, \;i=1,...,n$ > > Then $\bar y=\bar x +c$ > > From the table below, to compute the arithmetic mean of the time interval > between menstrual periods, it would be more convenient to work with numbers > that are closer to zero rather than $28$, which is the mode. > > ![enter image description here](http://i.stack.imgur.com/e951L.jpg) My questions are:   * Why is $28$ subtracted for the translated sample?    * Do you always subtract the mode to create a translated sample or would any value work? Why is this?    * Is it more accurate when you choose the mode rather than subtracting, say, 23?


title: Importing Live sound stream in R

text: This is the link for the live audio stream from a VLF receiver               http://67.207.143.181/vlf19.m3u      I need to know the method for importing live stream. If R is not suitable for such kind of signal processing, Please to suggest alternatives for processing and visualizing live data.


title: Sample size to tell if more than X% of the population can do <thing>

text: I want to perform a test to determine (with 95% confidence) whether at least 70% of a population can perform some task. The test involves sitting a randomly chosen person down and them attempting a task, which they either pass or fail. Equivalently, I flip a weighted coin. The tests are expensive, so we will only perform a dozen or so of them. The problem is to calculate pairs of (sample size, number of passes) that would give the required significance. I want to be able to say something of the form > "Ask 10 people. If 8 or more pass then you can say with 95 confidence that > the true ratio in the population is greater than 70%"


title: How do you know which predictor value to use?

text: I'm trying to **automate** linear regression with R, although I don't really have a concrete background in statistics. I was wondering: Are there numerical techniques in determining whether the predictor values is even worth trying to fit against a response value before attempting to do this in R:               lmfit <- lm(x ~ y)      My initial guess was covariance and correlation, but again what are the accepted values? Also, what if the relationship between response and predictor isn't linear. In most scenarios (as mentioned by commenters already) one would base predictors based on an a priori hypothesis. However; in knowledge discovery, there are times where you **_don't know_ what the hypothesis is**, hence my motivation towards formulating a set of rules whereby a quick numerical parameter would allow me to decide if a certain predictor variable should be added or removed in the linear regression model. Note that the discussion above is not data dependent, it should work with any general data that may or may not have linear relationships between variables.


title: Confusion regarding CRF

text: I was reading this paper related to conditional random fields http://www.inference.phy.cam.ac.uk/hmw26/papers/crf_intro.pdf. However, I have some confusion related to the section CRF probability as matrix computation. I couldn't actually grasp how it all fits in Can anyone please provide an example or two so that I can understand how it all fits in?


title: Help Simulating Random Effects in SAS

text: I'm trying to create a simulation of drug concentration based on the dose of a drug given. I have some preliminary data and I used a random effects model to analyze the relationship between log(dose), predicting log(drug concentration), modelling subject as a random effect. The results of that analysis are below. I want to take these results and simulate similar data in SAS, so I can look at the effect of changing doses on the resulting concentration of drug in the body. I know that when I simulate the data, I need to ensure the random slope is correlated with the random intercept, but I'm unsure exactly how to do that. Random effects:               Formula: ~LDOS | RANDID       Structure: General positive-definite, Log-Cholesky parametrization                  StdDev             Corr        (Intercept) 0.15915378        (Intr)       LDOS        0.01783609          0.735       Residual    0.05790635                Fixed effects: LCMX ~ LDOS                   Value   Std.Error  DF   t-value      p-value       (Intercept) 3.340712     0.04319325   16    77.34339        0       LDOS        1.000386      0.01034409   11    96.71090        0            Correlation:                    (Intr)       LDOS -0.047


title: How to do goodness-of-fit and residual diagnosis for rlm in R?

text: I am reading What are common statistical sins?, and especially @jebyrnes answer: > Failing to test the assumption that error is normally distributed and has > constant variance between treatments. These assumptions aren't always > tested, thus least-squares model fitting is probably often used when it is > actually inappropriate. and @DikranMarsupial's subsequent comment: > If the data are heteroscedastic you can end up with very innacurate out of > sample predictions because the regression model will try too hard to > minimise the error on samples in areas with high variance and not hard > enough on samples from areas of low variance. This means you can end up with > a very badly biased model. It also means that the error bars on the > predictions will be wrong. I am worried about my regression: I used "rlm" with weights. My questions are:   1. Please give me some pointers about how to do goodness-of-fit and residual diagnosis for rlm in R? Any good tutorials/examples, etc?   2. How to do the residual diagnosis etc. when regression weights are used?   3. My end goal was to get the predictions $\hat y$... and I am less concerned about the bias or efficiency of the betas... Do my predictions $\hat y$ get messed up if there are heteroskedasticity, non- normality, etc? i.e. what's the relationship between efficiency/unbiasedness and good out- sample prediction? Thank you!


title: Dimension reduction: preserving original characteristics of data or not?

text: I am new in dimension reduction techniques, I need to know what are the advantages and disadvantages of   * Preserving the original pattern of data after making data reduction    * Making reduction and affecting the original pattern of data And please, could you name for me these methods. Is it possible to give me some cases where it is preferable to preserve or not to preserve the pattern of original data.


title: Dummy Variables and Learning algorithms

text: Suppose a predictor variable $x$ is nominal/categorical with three levels: $1,2,3$. Thus we create two dummy variables $x_2$ and $x_3$ with level $1$ as the reference variable. Let $y$ be a binary response variable. Then a logistic regression model could be: $$ \text{logit}(p) = \beta_0 +\beta_1 x_{2} + \beta_{2}x_{3}$$ In SVMs and decision trees, is it necessary to create the dummy variables $x_2,x_3$ for $x$? Or can we just keep it as is?


title: AB testing vs testing the null hypothesis

text: I'm trying to understand the difference between   1. testing the null hypothesis (i.e. testing that the probability of a "goal" is the same across 2 different populations, similar to prop.test in R)    2. an A/B test using a bayesian formula such as described here: http://www.evanmiller.org/bayesian-ab-testing.html Is there a difference? Is one preferable? The problem I'm facing looks something like this: control group has 100,000 impressions and 100 reactions test group has 50,000 impressions and 55 reactions


title: J stat problem, GMM

text: I have recently performed a GMM estimations, my problem is that all the J-stats are 0.0000. It means that the IV are overrefined right or the model is not well specified. I used one-period lags of the independent variables as interments. Should I just report that the GMM model is misspecified or even if the J stat is 0.0000, i can make some inferences about the coefficients.?


How is covariation and homoscedasticity related? In my lecture notes for stats class we're told that an assumption pf Pearons correlation coefficient is that covariation must be linear, which I assume must be uniform? And to test for linear regression the data must be homoscedastic. How are these terms interrelated?


title: Null hypothesis of t-test and ANOVA

text: I've recently started to brush up my statistics, and have begun doing so with two of the most often used statistical tools in my field: the t-test and ANOVA. Reading up on popular information sources such as wikipedia (none of which provided citable references for the following statements - sadly) I have concluded that the null hypotheses tested by these tools are:   * **t-test** : _The t-test is a method used to test the null hypothesis that means of two compared groups are equal._   * **(one-way) ANOVA** : _Formally, ANOVA tests the null hypothesis that the means of multiple defined groups are equal._ My adviser told me both these statements are very wrong, and that _“inference- statistical methods do not test the equality of means”_ ; but offered no further clarification. I'm wondering whether I'm misrepresenting the information I have read. So, could you tell me what the null hypotheses of these tests are then? also, do you know what sources I could cite for this (since it's obvious I'm not stating my own independent conclusions)?


title: Reconciling two definitions of 'uncorrelatedness'

text: In this paper, the authors defined uncorrelatedness in the following way: Let $\mathbf{X}=(X_1,...,X_n)$, and $\mathbf{Y}=(Y_1,...,Y_n)$, where $X_i\sim X$ and $Y_i\sim Y$. $\mathbf{X},\mathbf{Y}$ are uncorrelated iff $(\mathbf{X}-\bar{X}\mathbf{1})'(\mathbf{Y}-\bar{Y}\mathbf{1})=0$, where $\bar{X}=1/n \sum X_i$ and $\mathbf{1}=(1,...,1)$. My question is what is the link between this definition of uncorrelatedness and the one usually used (I think), that of (for vectors r.v.) covariance matrix - a matrix where each entry $(i,j)$ is $\text{Corr}(X_i,X_j)$ - being an identity matrix? Any help would be appreciated.


title: Rule of thumb to rule out reverse causality in the OLS model

text: Let' say I have a regression model: $y=a+b*x+error$ Suppose $x$ is income and $y$ is consumption. The hypothesis is that higher income leads to higher consumption and hence, the coefficient on $x$ should be positive, other things remaining the same.Let's also say the estimated coefficient is 0.60. This model obviously suffers from omitted variable bias. Please ignore this issue. My question: a) Does this model suffer from reverse causality? In other words, is it the case that the relationship is because higher consumption is driving down income? My first guess is that this is not the case because coefficient is positive which means correlation between income and consumption is positive. See here. b) Given (a), can I use this as a rule-of-thumb to rule out the reverse causality in this case? Is this generalizable to other cases with two variables? Thanks. P.S. One can also just use correlation rather than running a simple regression as mentioned earlier.


title: Routing Table and Splitting Traffic in omnet++

text: I am new to Omnet++ 4.3, I have a WMN network and I want to send traffic through two shortest paths, from one of my client nodes to two ONU(wgRouter).   1. First I want to know how can I find the 2 shortest paths, is it possible to run the network and find the first shortest path through one of the wgRouter and save the routing table and then turn off that wgrouter and then run the model again and find the second shortest path through another wgrouter and again save the routing table ,then turn on the off wgRouter, use both routing tables in initialization part ( I dont know how can I use them both in one time, please let me know if you know) and run the model ?   2. let consider that we found the two shortest path through two wgRouters, Now how can I split the traffic equally between both of them?  any information would be very useful for me. Thanks.


title: Statitical test on count data

text: I have a questions regarding a statistical analysis of a biologic experiment. In the experiment I have collected imaging data on the duration of contacts between immune cells. We have compared two distinct types of immune cells and their contact time with a third cell type. The duration of contacts was measured in minutes. Now I want to determine if there are significant differences number of contacts that lasted longer than 30 minutes. This gives mit counts for the two cell types. Example: Cell type 1: 23 out of 140 contacts measured; cell type 2: 34 out of 126 contacts measured. Which statistical test do I need to use in this case?


title: Posterior covariance of Normal-Inverse-Wishart not converging properly

text: I am trying to implement a simple normal-inverse-Wishart conjugate prior distribution for a multivariate normal with unknown mean and covariance in numpy/scipy such that it can take a data vector and construct a posterior. I'm using the update equations specified by Wikipedia for a NIW: http://en.wikipedia.org/wiki/Conjugate_prior My distribution class is as follows:               import numpy as np     from scipy.stats import chi2          class NormalInverseWishartDistribution(object):         def __init__(self, mu, lmbda, nu, psi):             self.mu = mu             self.lmbda = float(lmbda)             self.nu = nu             self.psi = psi             self.inv_psi = np.linalg.inv(psi)              def sample(self):             sigma = np.linalg.inv(self.wishartrand())             return (np.random.multivariate_normal(self.mu, sigma / self.lmbda), sigma)              def wishartrand(self):             dim = self.inv_psi.shape[0]             chol = np.linalg.cholesky(self.inv_psi)             foo = np.zeros((dim,dim))                  for i in range(dim):                 for j in range(i+1):                     if i == j:                         foo[i,j] = np.sqrt(chi2.rvs(self.nu-(i+1)+1))                     else:                         foo[i,j]  = np.random.normal(0,1)             return np.dot(chol, np.dot(foo, np.dot(foo.T, chol.T)))              def posterior(self, data):             n = len(data)             mean_data = np.mean(data, axis=0)             sum_squares = np.sum([np.array(np.matrix(x - mean_data).T * np.matrix(x - mean_data)) for x in data], axis=0)             mu_n = (self.lmbda * self.mu + n * mean_data) / (self.lmbda + n)             lmbda_n = self.lmbda + n             nu_n = self.nu + n             psi_n = self.psi + sum_squares + self.lmbda * n / float(self.lmbda + n) * np.array(np.matrix(mean_data - self.mu).T * np.matrix(mean_data - self.mu))             return NormalInverseWishartDistribution(mu_n, lmbda_n, nu_n, psi_n)      I am running a simple sanity check to see if the posterior converges to the true distribution:               x = NormalInverseWishartDistribution(np.array([0,0])-3,1,3,np.eye(2))     samples = [x.sample() for _ in range(100)]     data = [np.random.multivariate_normal(mu,cov) for mu,cov in samples]     y = NormalInverseWishartDistribution(np.array([0,0]),1,3,np.eye(2))     z = y.posterior(data)          print 'mu_n: {0}'.format(z.mu)          print 'psi_n: {0}'.format(z.psi)      The mean is appropriately converging, but the scale matrix appears to be converging to incorrectly large values along the diagonal, rather than the true value of 1. As far as I can tell, I'm copying the update rule exactly. Am I implementing something inappropriately here? **Edit** : It looks like in fact the posterior _is_ converging, but the sample routine is returning biased samples. Am I doing something wrong in my sampling method? **Edit2** : I've confirmed that the same phenomenon happens in the MCMCpack riwish function in R:               > library(MCMCpack)     > samples <- replicate(100000, riwish(3, matrix(c(1,0,0,1),2,2)))     > mean(samples[1,1,])     [1] 4.889211      This leads me to believe I must be misunderstanding something. From the Wikipedia page (http://en.wikipedia.org/wiki/Inverse-Wishart_distribution), we have: $\newcommand{\E}{\mathrm{E}}$ $\E[{\Sigma}] = \frac{\Psi}{\nu - p - 1}$ However, in my test case, $\nu = p + 1$, so $\E[{\Sigma}] = \Psi$. Thus, if I sample $\Sigma$ a bunch of times, shouldn't I recover $\Psi$ on average? **Edit 3** : Realized my problem was simple algebraic oversight: I needed to set $\nu = p+2$. Problem solved


title: Code for conditional linear regression

text: I've just run a linear regression on an entire data set, but now I need to run the regression with data just from females within the data. Females are denoted under the `female` column of the data set by a `1`. Males are denoted by a `0` under the same column. I don't know how to remove the male data so I can run the regression on female data only.


title: Finding clusters to fit least squares and produce a piecewise equation

text: In the figure below, I've manually drawn an approximate solution to a least squares fit of the associated regions separated by black lines. The data appears to be bounded by two asymptotes (y=-18 and y = -75), which I then use for locating the regions. 1) What is a method for find these asymptotes and the associated three regions? 2) Is there a special name for this type of piecewise equation, which will be the result of LS fit on the data in each region? ![model](http://i.stack.imgur.com/5z7gA.png)


title: Mean comparisons following multiple imputation

text: I need to do some simple mean comparisons between groups (basic ANOVA F-tests) on data with missing values. I use the mice package in R for multiple imputation, but I can only pool results for the linear model coefficients, or the $R^2$. Does anyone know how to combine to pool multiple F-statistics from each linear model fit? Or, how can I compute the standard errors for the F-test?


title: Differences for PROC GENMOD, COUNTREG and FMM for count data model

text: I want to fit a negative binomial model for my count data. I am wondering what the differences for the computation of these three PROCs are in SAS.


title: Vowpal Wabbit: understanding prediction output

text: I've trained a classifier model using Vowpal Wabbit to decide if a person is Male or Female based on name alone. I assigned labels Male=0 and Female=1. When I ran Vowpal Wabbit in prediction mode, the output had values varying between 0.0 and 1.0. I interpreted this as p <= 0.5 ==> Male and p > 0.5 ==> Female. With this the accuracy of prediction was 85%. I'm now trying to train a classifier model for another problem where there are 200 labels instead of just two like above. The labels go from 0 to 199. The prediction output has values between 0.0 and 199.0. When I get a value such as 99.56 (for example), how do I interpret this? Does it map to label 100?


title: Assessing accuracy of randomForestSRC on a survival data

text: In the randomForestSRC in R, we have the example               data(veteran, package = "randomSurvivalForest")     pt.train <- sample(1:nrow(veteran), round(nrow(veteran)*0.80))     veteran.out <- rsf(Surv(time, status) ~ ., data = veteran[pt.train , ])     veteran.pred <- predict(veteran.out, veteran[-pt.train , ])      How can I tell if my model predictions are are good using the veteran.pred results?


title: Interpreting odds in a logistic regression

text: I am having some problems interpreting the odds. I run a logistic regression for an out come 'Yes' or 'no'. My reference category is 'No'. I have 2 variables and this are the log(odds) and the odds: Variable A -> It is a integer with values 2 to 80. logodds: -0.014078787; odds: 0.9860199 Variable B -> It is a integer with values 1 to 30. logodds: 0.214099984; odds: 1.2387465 What is a correct statement to interpret these odds? And how do I transform this in real probabilities? I would like to say something like: Every extra point of variable A causes an decrease of x in the probabilities of going from 'No' to 'Yes' If I do exp(coef(model)) /(1+exp(coef(model))) for Variable A I obtain 0.49. This is the probability right? Is it correct to say that each point of variable A makes the probability go down 0.01% I am really confussed about this so would thank any help!


title: Restriction on Cholesky decomposition

text: In my model I'm using a Cholesky decomposition to produce a certain covariance structure of a vector of normal variables: $\Omega=LL'$ Where L is the lower triangular Cholesky factor and $\Omega$ is a covariance matrix. Right now $\Omega$ is unrestricted, but I want to restrict $\Omega$ to have a specific structure. For instance let's assume that $\Omega$ is $2 \times 2$ and I want to set the variance of the first variable to one: $\Omega=\left(\begin{array}{cc} 1 & \sigma_{21} \\\ \sigma_{21} & \sigma_{22} \end{array} \right)$ It is easy to restrict the covariance matrix itself, but the problem is in my R program, I'm using the Cholesky factor to produce the covariance matrix (this is because among others that the Cholesky factor ensures a positive definite covariance matrix), so I need to figure out how to restrict the Cholesky factor in order to produce the correct restricted covariance matrix. In the very easy case above, I can make the restriction on the covariance matrix by setting the first element of the Cholesky factor equal to zero: $\Omega=\left(\begin{array}{cc} 1 & \sigma_{21} \\\ \sigma_{21} & \sigma_{22} \end{array} \right)=LL'=\left(\begin{array}{cc} 1 & 0 \\\ l_{21} & l_{22} \end{array} \right)\left(\begin{array}{cc} 1 & l_{21} \\\ 0 & l_{22} \end{array} \right) =\left(\begin{array}{cc} 1 & l_{21} \\\ l_{21} & l_{21}^2+l_{22}^2 \end{array} \right)$ So my question is basically: Is there a general approach to restricting the Cholesky factor in order to get the correct covariance matrix?


title: How do you calculate the amount of parameters needed to be estimated?

text: I don't quite understand this. A question was, pretend we have 4 predictors and all of them are binary - for the Naive Bayes method, how many parameters are there to estimate in the training step? The answer was 2*4. Then, another question was, if we draw an edge between predictors 1 and 2 to create a graphical model, how many parameters do we need to estimate in the training step? The answer was 1+1(2)+(2*7)=11. Finally how about if all of the predictors have edges in between? The answer was 1+2(15)=31. I don't get how these equations are being constructed. Could someone explain please?


title: How do you interpret the interquartile range?

text: My data is: Q1: 8.5 Q2: 13 Q3: 17.5 IQR:9. The math problem says to interpret the interquartile range. I have no idea what this means or how to do it.


title: Does preclustering help to build a better predictive model?

text: For the task of churn modelling I was considering:   1. Compute k clusters for the data   2. Build k models for each cluster individually. The rationale for that is,that there is nothing to prove, that the population of subsribers is homogenous, so its reasonable to assume that data-generating process may be diffrent for diffrent "groups" My question is, is it an appropriate method? Does it violate anything, or is it considered bad for some reason? If so, why? If not, would you share some best practices on that issue? And 2nd thing - is it generally better or worse to do preclustering than model tree (As defined in Witten,Frank - classification/regression tree with models at the leafs. Intuitively it seems that decision-tree stage is just another form of clustering, but idk if it has any advantages over "normal" clustering.).


title: What are some general rules or guidelines about choosing between color and 3 axes for 3D plotting?

text: For plotting in 3 dimensions, I would say that there are two top level approaches.   1. Plot with 3 axes: one dimension per axis.   2. Plot with 2 axes with the 3rd dimension represented in color. What are some general rules or guidelines about choosing between these two approaches?


title: How does one interpret SVM feature weights?

text: I am trying to interpret the variable weights given by fitting a linear SVM. (I'm using scikit-learn):               from sklearn import svm          svm = svm.SVC(kernel='linear')          svm.fit(features, labels)     svm.coef_      I cannot find anything in the documentation that specifically states how these weights are calculated or interpreted. Does the sign of the weight have anything to do with class?


title: Unwritten laws and dirty tricks to influence the outcome of a regression analysis

text: Background: Some years ago, I was working for a professor specialized in macroeconometrics. As a student research assistant, it was my task to replicate other papers and "play around" with the data. I learned quickly that many published results were not even robust to small changes in the modeling setup (e.g. if you omit one variable, change the data range, add more lags). In theory, the researcher hypothesizes, does one regression and reports the results. But in practice it seems to be an unwritten law to estimate hundreds of regressions and publish one. My question: What are unwritten laws and dirty tricks to influence the outcome of a regression analysis? Do you know some reference (academic/ semi-academic) where the author describes such an informal procedure to obtain preferable results?


title: Quality Control Chart for non-normal distributed data

text: I am interested in building quality control charts for time series data where every data point has different gradation within the range of 0 to 1. For example first data point **_X1_** may be 0.1, 0.2,..,1 and second data point **_X2_** may be 0, 1/3, 2/3, 1. Third data may have some another interval. The data is not seasonal and there is no pattern in how the intervals look like (e.g., 0.1( ** _X1_** ), 1/3( ** _X2_** ) etc.). In other words I don't know what is the gradation of future data points but know that the range is between 0 and 1. My questions are 1) Wether the limits that I will build based on some training data (UCL, LCL) will be effected as a result of having different gradation - the variance will be higher since the data is not continuous? 2) Is it a problem given the fact that this represents the nature of future data and that I have a training set of 30 samples? 3) The data points represent performance of student over different tests in consequent time, and were normalised to be between 0 and 1 (10 questions in test will lead to interval of 0.1, 3 questions will lead to interval of 1/3). Maybe there are some ways to normalise the data to avoid the problem of blowing up the variance? Many thanks for your patience and help!


title: consequences of lagged dependent variables in panel data and how to deal with it?

text: I have some elementary problems understanding the consequences of using/adding a lagged dependent variable in my predictive model. I’m trying to predict values $Y_{i,t+\tau}$ for $\tau=1-3$ with: $Y_{i,t+1}=a+bY_{i,t}+cX_{i,t}+e_{i,t+1}$ $Y_{i,t+2}=a+bY_{i,t}+cX_{i,t}+e_{i,t+2}$ $Y_{i,t+3}=a+bY_{i,t}+cX_{i,t}+e_{i,t+3}$ I already performed a pooled regression where you basically ignore individual firm effects and time-effects and treat every subject equally. As I am trying to forecast different levels (in USD) and my data appears to be extremely tailed as it covers a few extremely large subjects (with extremely high values) but also many small subjects the predictions of the model perform rather poor as the intercept $a$ that is equal for all subject seems largely responsible for this. A fixed model however with individual intercepts is not valid with Lagged dependent variables as the LDV is correlated with the within errors. To account for the heavy tailed errors I already estimated the pooled model with the rlm package (robust lm) that produced slightly better results but overall they appear still very unsatisfactory. I further read that adding LDV’s results in biased and inconsistent estimators as there is severe correlation between the predictor variables and the model errors and that regular procedures for autocorrelation are not valid anymore. One solution I came across is the use of Instrumental Variables with an Anderson-Hsiao Estimator (i.e using a lag -2 that is not correlated with the error term (with non-autocorrelation assumed but how can you assume no autocorrelation if you incorporate a lag?) Another one is the Arellano Bond GMM estimator, however applying GMM you have to set up moment conditions and I have no idea how to do that and I don’t know exactly how this methods work. What I care about is to obtain an unbiased estimator with valid coefficients and not about standard errors as I don’t do inference. Are there any other strategies to cope with LDV’s I am currently unaware of and what is the best/ideal/easiest way to deal with such matter? Do you best take care of some issues while you ignore others (e.g. autocorrelation)? I’m a little bit lost here.


title: Arithmetic vs Geometric Mean

text: There are verbal suggestions everywhere on when one should use a geometric average or when an arithmetic average should be preferred, but I can't find any formal statistical treatment of this question. Is it possible to formally test which one of these averages should be used for a particular sample?


title: Finding predictions from a linear model

text: I'm trying to educate myself about predictive analytics and am using `R` to generate a linear model with the following data.               age     <- c(23,   19,   25,   10,9, 12,   11, 8, 13)     steroid <- c(27.1, 22.1, 21.9, 10.7,  7.4, 18.8, 14.7, 5.7)     gpa     <- c( 2.1,  2.9,  2.8,  3.5,  3.2,  3.9,  2.8, 2.6)     sample  <- data.frame(age, steroid, gpa)     fit2    <- lm(steroid~age+gpa)     summary(fit2)     newdata <- data.frame(age=15, gpa=3.2)     predict(fit2, newdata, interval="predict")  # I want the fitted / predicted value      From the summary for the linear model, I have information on the coefficients associated with each predictor. However, I want to go further and find the predicted probability of information from this data. So for someone who is 13 and has a GPA of 3.3, what's the predicted probability of them ranking high on the steroid scale? What about for someone who is age 10? etc. I understand explanatory modeling and have no problem implementing that in `R`, I just have issues with using `R` to construct meaningful predictive analytics. If you have any insights on this, feel free to share.


title: Is there a minimum sample size required for the t-test to be valid?

text: I'm currently working on a quasi-experimental research paper. I only have a sample size of 15 due to low population within the chosen area and that only 15 fit my criteria. Is 15 the minimum sample size to compute for t-test and F-test? If so, where can I get an article or book to support this small sample size? This paper was already defended last Monday and one of the panel asked to have a supporting reference because my sample size is too low. He said it should've been at least 40 respondents.