title: Larger ngrams vs nested ngrams?

text: Whenever I see people using ngrams - I often see them looking for higher chains like 4-grams, 5-grams and so on. However, I'm wondering why I never see any mention of "nested ngrams" (I'm not sure what the real name is). For example, suppose you tracked the bigram "to follow" followed by the bigram "after dark" instead of the 4-gram "to follow after dark". Is there any reason to look into checking multiple nested levels of smaller ngrams instead of continuing to grow larger and larger ngram chains? It seems that one benefit might be storage costs.


title: How to choose the number of nodes in layer recurrent networks?

text: When using NNs for estimation purpose, we have three layers: input, output and hidden layer, and we must choose the number of nodes for every layer. Is there a method or rule of thumb to choose the correct number of nodes?


title: How to get proper randomization in website A/B testing?

text: In the statistical practice of experimental design, you separate your tests by 'blocks' if you can control the factors (say routing to one page or another, or categorizing by browser) or through randomization of trials to counteract variables you don't control. When doing A/B testing, you obviously can't control who comes to your site when. How can good statistical practice be done?


title: Are there any R functions which support Reversible Jump MCMC for a GLM or SGLM?

text: I would like to perform reversible jump on a network model, but before arriving there, I'm wondering if there are any R packages which support reversible jump for a user specified generalized linear model or spatial-GLM? Something as simple as an RJMCMC procedure (in R) for the selection of predictors in a logistic regression would be a nice place for me to start? Does such a function exist? Through googling, I've only found RJaCGH which appears to be a bit more complicated (and application specific) than I was hoping for.


title: Random numbers and the multicore package

text: When programming in R, I've used the multicore package a few times. However, I've never seen a statement about how it handles it's random numbers. When I use openMP with C, I'm careful to use a proper parallel RNG, but with R I've assume that something sensible happens. Can anyone confirm that something sensible does happen? **Example** From the documentation, we have               x <- foreach(icount(1000), .combine = "+") %do% rnorm(4)      How are the `rnorm``s generated?


title: Unbalanced repeated measure design for the given data?

text: I have 3 factors factor1 (treatment) with 2 levels (control, stress) factor2 (Variates) with 12 Levels (Var1, Var2,....Var12) factor3 (Time) with 12 levels (Week1, Week2,..., Week12) The treatment has 3 replicates for control and 6 replicates for stress. Would it be unbalanced Design? What design would you suggest in this case? looking forward for your positive response! Thanks in advance! Regrads, Jacki


title: Inter-rater reliability for ordinal or interval data

text: Which inter-rater reliability methods are most appropriate for ordinal or interval data? I believe that "Joint probability of agreement" or "Kappa" are designed for nominal data. Whilst "Pearson" and "Spearman" can be used, they are mainly used for two raters (although they can be used for more than two raters). What other measures are suitable for ordinal or interval data, i.e. more than two raters?


title: How well can the Dirichlet process cluster really small datasets?

text: I have been debating between a model-based parametric clustering approach (e.g. HMMs), and a hierarchical Dirichlet/Pitman-Yor process for clustering sequential data. I understand the latter has been very popular in the recent past (no model selection, heavy tails, etc.), but I could not find any analysis of how it performs on really small datasets. I suspect that HMMs would be better, but I would really appreciate if someone could provide additional insight or point me to relevant texts on this topic.


title: How do I identify this trap error in WinBUGS?

text: I am currently working on my thesis and interested in estimating a Weibull conditional Hazard Frailty Model for recurrent event data using WinBUGS. I wrote the model but when am trying to use it with the following script it hangs unexpectedly. It returns a messages titled "undefined real result". What am I doing wrong? Also, I have tried different prior and initial values but I could not solve the problem!               model     {     for(i in 1:633){     ht[i]<-gama*pow((t[i]),(gama-1))*exp(age[i]*b[1]+sex[i]*b[2]+mar[i]*b[3]+back1[i]*b[4]+back2[i]*b[5]+form[i]*b[6]+u[subject[i]]+delta*enu[i])     }     for(i in 1:633){     beta[i]<-pow((1/gama),((gama*pow((t[i]),(gama-1)))/ht[i]))     }     for(i in 1:633){     t[i]~dweib(gama,beta[i]) I(cen[i],)     }     for (j in 1:6){     b[j]~dnorm(0,0.001)     }     for(i in 1:159){     u[i]~dnorm(0,tau)     }     delta~dnorm(0,0.01)     gama~dgamma(1,0.01)     tau~dunif(0,100)     sigma2.subject<-1/tau     }     list(b=c(0,0,0,0,0,0),gama=1,tau=1,delta=0)      I would appreciate it if you could help me for solving this problem.


title: Analysis of Customer satisfaction surveys

text: I have customer feedback data about 2-3 products from 100 customers. Number of questions are around 160. I have data in excel format. Header row contains the question and row below contains the responses from customers. Not all questions have been answered by each customer. I have to categorize the feedback as either positive or negative. The analysis needs to be done in R. This data contains boolean responses, numerical responses and text responses. I've collected the data for text responses in a separate sheet and trying to do sentiment analysis using sentiment analysis package and the general text mining package in R. How to do the sentiment analysis of boolean and numerical responses ? Also, Can I combine the score from all the above three analysis and can categorize the response as positive or negative overall ?


title: Linear Regression makes impossible predictions

text: I have created a multiple linear regression model to predict prices and about 17% of the predicted prices have come out to be negative. Is there a way to correct for this error, or does it mean that the independent variables are not good predictors? The R^2 coefficient is 0.865.


title: Odds ratios from logistic regression on subsets of data using R

text: I have 15 subjects each with 200 trials & I'd like to run separate regressions for each subject. If I just run the regression on the whole dataset I am able to generate odds ratios / confidence intervals from it:               acc.log = glm(response ~ phdot + sdot, data=fullerr, family="binomial")     summary(acc.log)          exp(acc.log$coefficients)     exp(confint(acc.log))      But then when I fit each subject separately (I've tried using `by`), I am an unsure how to then extract information like odds ratios (standard errors & p values are also no longer shown in the output).               split.acc <- by(fullerr, fullerr[,"subject"], function(fullerr){                        glm(response ~ phdot + sdot, data=fullerr, family="binomial")                     })


title: Should types of data (nominal/ordinal/interval/ratio) really be considered types of variables?

text: So for instance here are the definitions that I get from standard text books Variable - characteristic of population or sample. ex. Price of a stock or grade on a test Data - actual observed values So for a two column report [Name | Income] the column names would be the variables and the actual observed values {dave | 100K} , {jim | 200K} would be the data So if I say that the [Name] column is nominal data and that [income] is ratio data, wouldn't I be more accurate describing it as a type of variable instead of a type of data like most textbooks do? I understand that this might be semantics, and that's fine it that's all there is too it. But I fear that I might be missing something here.


title: How to see if one (or one set of commonly answered) question(s) impacts/statistically influences another question/set of questions?

text: I'd like to review the results of a recent survey I conducted. I'd like to see how one set of answers may/may not have influenced another set of answers. For example: If users who answered question #1 on a 1-8 Likert scale with answers from 6-8, does this correlate to how they answered questions 2-8? Or: Users who answered false for question 10 did/did not similarly answer questions 11-12.


title: Division of two poisson variables

text: Let's say I have four variables $x_1, x_2, x_3$ and $x_4$ all distributed as Poisson, with means $\lambda_i$, and standard deviations respectively. Now I want to calculate two further variables $y_1= \frac{x_1}{x_2}$ and $y_2 = \frac{x_3}{x_4}$. Finally I have a function $f(y_1, y_2)=(y_1-n)/(y_2-m)$. How can I calculate the standard deviation of this function? Is it a simple propagation of uncertainty?


title: Which SVM kernel to use for a binary classification problem?

text: I'm a beginner when it comes to support vector machines. Are there some guidelines that say which kernel (e.g. linear, polynomial) is best suited for a specific problem? In my case, I have to classify webpages according to whether they contain some specific information or not, i.e. I have a binary classification problem. Can you say in general which kernel is best suited for this task? Or do I have to try several of them on my specific dataset to find the best one? By the way, I'm using the Python library scikit-learn that makes use of the libSVM library.


title: What is the two-sample CDF of of $D^{+}$ and $D^{-}$ from the one-sided Kolmogorov-Smirnov Test?

text: I am trying to understand how to obtain $p$-values for the _one-sided_ Kolmogorov-Smirnov test, and am struggling to find CDFs for $D^{+}_{n_{1},n_{2}}$ and $D^{-}_{n_{1},n_{2}}$ in the two-sample case. The below is cited in a few places as the CDF for $D^{+}_{n}$ in a one-sample case: $$p^{+}_{n}\left(x\right) = \text{P}\left(D^{+}_{n} \ge x | \text{H}_{0}\right) = x\sum_{j=0}^{\lfloor n\left(1-x\right)\rfloor}{ \binom{n}{j} \left(\frac{j}{n}+x\right)^{j-1}\left(1 - x - \frac{j}{n}\right)^{n-j}}$$ Also, whuber sez there's a slightly different formulation of this one-sample CDF (I am substituting $x$ for $t$ in his quote for consistency with my notation here): > Using the probability integral transform, Donald Knuth derives their > (common) distribution on p. 57 and exercise 17 of TAoCP Volume 2. I quote: $$\left(D^{+}_{n}\le \frac{x}{\sqrt{n}}\right)=\frac{x}{n^{n}}\sum_{c\le k\le x}\binom{n}{k}\left(k-x\right)^{k}\left(x+n-k\right)^{n-k-1}$$ This would apply to one-sided hypotheses in the one-sample case, such as: H$_{0}\text{: }F(x)-F_{0} \le 0$, where $F(x)$ is the empirical CDF of $x$, and $F_{0}$ is some CDF. I _think_ the $x$ in this case is the value of $D^{+}_{n}$ in one's sample, and that $\lfloor n\left(1-x\right)\rfloor$ is the largest integer in $n-nx$. (Is that right?) But what is the CDF for $D^{+}_{n_{1},n_{2}}$ (or $D^{-}_{n_{1},n_{2}}$) when one has two samples? For example, when H$_{0}\text{: }F_{A}(x)-F_{B}(x) \le 0$ for the empirical CDFs of $A$ and $B$? How to obtain $p^{+}_{n_{1},n_{2}}$?


title: y-axis length for LOWESS curve in Stata

text: I have a problem considering the length of the _y_ -axis in a LOWESS plot in Stata. The _y_ variable ranges from 0 to 50. I would like to run the LOWESS on all the data but only display the graph for the range of _y_ from 30 to 50. I've searched throughout the net, but I just found that I would also have to restrict the data to $\ge$ 30, but that of course changes the curve. Do you have any advice? The graphs are simply unusable for that _y_ -axis interval. Thanks!


title: Distribution of norm of random matrix

text: I am curious to know whether central-limit theorem like considerations hold true for special functions like the norm of a matrix. Specifically, I'm interested in the Spectral norm of a matrix $$\|M^\prime\| = \text{Max eigenvalue of } (M^{\prime H}M^\prime),$$ where $(\bullet)^H$ represents the conjugate of the transpose of $(\bullet)$. Here, $M^\prime$ is close to matrix $M$ such that each element of ${M}^\prime$ has random error ${M}^\prime_{ij} = \mathcal{N}(M_{ij}, \sigma^2)$. Question: **Will the distribution of $\|\tilde{M}\|$, which is a function of random variables $\tilde{M}_{ij}$ be a normal distribution?** I have tried to figure this out numerically. My numerics fit the normal distribution rather well. Is there a way to prove this analytically?


title: Power analysis for Kruskal-Wallis or Mann-Whitney U test using R?

text: Is it possible to perform a power analysis for the Kruskal-Wallis and Mann- Whitney U test? If yes, are there any R packages/functions that perform it?


title: Random vector times random matrix

text: $$ Y_i = \sum_j B_{ij} X_j$$ $$ covar(X_i, X_j) = V_X = \delta_{ij} var(X_i)^2 $$ $$ covar(B_{ij},B_{kl}) \neq 0 $$ $$ V_Y = ? $$ I know, that if $B$ was fixed, it is straight forward, but I would like to generalize it. I managed only by inventing some notation myself, which a physicist should not do I guess. I will check my method later with Monte Carlo, but maybe there is a formal way to write it down? Maybe the answer should be generalized from variance-of-product- of-2-independent-random-vector? **Edit** We can also assume that $covar(X_i,B_{kl})=0$. **Edit2** Replaced Einstein Summation rule with a sum. I don't just skip the indexes as I want to to stress the dimensions and make the latter more readable.


title: Is two-way ANOVA appropriate?

text: This is the description of my study. I'm experimenting with three plants: A, B, and C. These plants are supposed to reduce blood glucose for diabetic patients. I want to determine which of these three plants has a longer effect on blood glucose reduction after a single administration into mice. This is done by measuring the blood glucose from mice at 7 time points (day 1, 2, 3, 5, 7, 10, and 14). So there are 4 groups (untreated, treated with A, treated with B, and treated with C). For each group 3 mice were used (n=3). My goals are:   1. To determine if the effect of each plant treatment is significant compared to untreated.   2. To compare the effect of the among the groups for each day.   3. To determine which treated group has the longest effect after 14 days. My solution is to use two-way ANOVA since there are more than 2 groups and I want to compare the groups in each day and finally the overall effect. Is the this the correct method? Will I be able to rank which will be the best plant followed by the 2nd and the 3rd? Or should I use time series analysis?


title: Likelihood Function for Complicated Transformations

text: Suppose that data _X_ have a Normal distribution with some mean $\mu$ and some variance $\sigma^2$. However, you don't get to see _X_. Instead, you see $Y = g(X)$ where $g$ is a known function. Assume that $g$ is complicated and is not invertible. How do you write the likelihood function for $Y$? Or how do you use some numerical method to approximate a likelihood in a way that allows inference about $\mu$ or $\sigma^2$? Seems like a basic question that has totally stumped me.


title: A randomized, double-blind crossover (AB/BA) observational study of new use of an FDA approved drug

text: I found a new use of a drug already in use. The results are as follows. Serail No. of Patient-- Attacks while on Placebo-- Attacks while on drug --Difference in the no. of attacks 1- 23- 2- 21 2- 17- 1- 16 3- 13- 3- 10 4- 26- 0- 26 5- 16- 0- 16 6- 11- 1- 10 7- 22- 3- 19 8- 15- 1- 14 9- 19- 2- 17 10- 16- 1- 15 11- 22- 1- 21 12- 27- 0- 27 13- 31- 3- 28 14- 24- 2- 22 15- 22- 1- 21 16- 19- 2- 17 17- 16- 1- 15 18- 25- 3- 22 19- 19- 1- 18 20- 21- 2- 19 21- 6- 6- 0 22- 28- 3- 25 23- 19- 1- 18 24- 21- 2- 19 25- 27- 1- 26 26- 20- 1- 19 27- 17- 1- 16 28- 25- 1- 24 29- 18- 1- 17 30- 19- 1- 18 31- 5- 5- 0 32- 17- 1- 16 33- 28- 1- 27 34- 4- 4- 0 I am sorry I am not able to format the table correctly. I have very basic knowledge of SPSS. Kindly advise how to go about.


title: ncvTest from R and interpretation

text: I have done a ncvTest, but I am not really sure how to interpret this properly. I look documentation and examples online but was not able to find anything that clearly explains about ncvTest. Here are my results from ncvTest:               Non-constant Variance Score Test      Variance formula: ~ fitted.values      Chisquare = 3.203987    Df = 1     p = 0.073459       ![enter image description here](http://i.stack.imgur.com/msyem.png)


title: Rule of thumb about how many data points are needed to train a GMM

text: I have about $10^6$ data vectors, each with 13 dimensions. I'm trying to train a GMM of 256 components using EM algorithm (using GMMBAYES). Is there any rule of thumb about how many data points are generally sufficient to train a N dimensional, M component GMM?


title: State of the art method(s) to find zero mean portions of a time series

text: I have noisy time series which I need to segment into those portions with a zero mean and those portions without a zero mean. Finding the boundaries as accurately as possible is important (clearly where the boundary precisely lies is a bit subjective). I think a cusum variant could be adapted to do this but as cusum is primarily about finding single changes that leaves the whole segmentation strategy completely unaddressed. I'm sure a bunch of research has been done on this problem but have not been able to find it. P.S. The amount of data in these time series is quite large, i.e. up to hundreds of millions of samples, and an individual sample can be a vector with a couple of hundred components, so a method that can be computed reasonably quickly is a significant factor. P.P.S There isn't a segmentation tag, hence the classification tag.


title: Statistical methods for analyzing symptoms correlation with nutrients

text: I am working on a web based project that will track a person's symptoms and drugs and nutrients they intake. For example thyroid patients symptoms of fatigue vs the dose of thyroid and also how much selenium, iodine, vitamin A, protein, etc they get from food or supplements. I've added a basic set of descriptive statistics functions, linear regression functions, but am not sure these are appropriate for this type of analysis. Most biological systems are very complex and multiple factors may influence a symptom. The math/statistics books I've skimmed seem to talk about science experiments where one variable is isolated which is not the case here. Could anybody recommend some reading materials or information that could help with this type of analysis? My background is in biology and programming, not statistics. Thanks, Dan


title: Limit on number of time varying covariates in extended Cox regression?

text: I'm working on an extended Cox proportional hazards regression model which is predicting a recurrent event, inpatient hospitalization, during a 2 year period and will include multiple time-varying covariates (at least 10). My understanding is that the counting process method incorporates multiple records per subject, with each record corresponding to an interval of time during which all variables (including time-varying variables) remain constant and each interval either ends in the event of interest or is right censored. Is there a limit on the number of time-varying covariates which can be reasonably included in a Cox regression model? I'm wary of sacrificing power by slicing up my time period into too many intervals.


title: Search algorithms for time series in machine learning

text: Apart from k nearest neighbour and its variants like fast K NN, what are the other search algorithms? I tried googling and reading but I did not get any effective answer. I am aware of binary search, heap search, simulated annealing, tabu search. But I do not know what are the time complexities and was looking for a minimum time complexity algorithm. In general what are the iterative search algorithms that can be applied to search a solution space say for example a parameter space that is computationally efficient and robust. For instance, application of simulated annealing in a solution space where the objective is to find a minimum nearest neighbor etc. Any resources will be highly useful.


title: Evaluating a logistic regression

text: Recently I’ve been working on a logistic model and I´m having some difficulties evaluating the results. I am still learning all this, so I apologise in advance for the mistakes. I would still appreciate any insight though. So, my model is a binomial logit. My explanatory variables are: 1 categorical variable with 15 categories, 1 dichotomous variable, and 2 continuous variables. My N is large +8000 I am trying to model the decision of Firms to invest. Dependent variable in Investment (yes/no), the 15 level variables are different obstacles for investments reported by managers. The rest of the variables are controls for sales, credits and used capacity. Below are my results, using the rms package in R.                 Model Likelihood     Discrimination    Rank Discrim.                                  Ratio Test            Indexes          Indexes            Obs          8035    LR chi2     399.83    R2       0.067    C       0.632          1           5306    d.f.            17    g        0.544    Dxy     0.264          2           2729    Pr(> chi2) <0.0001    gr       1.723    gamma   0.266         max |deriv| 6e-09                          gp       0.119    tau-a   0.118                                                    Brier    0.213                                         Coef    S.E.   Wald Z Pr(>|Z|)     Intercept -0.9501 0.1141 -8.33  <0.0001      x1=10     -0.4929 0.1000 -4.93  <0.0001      x1=11     -0.5735 0.1057 -5.43  <0.0001      x1=12     -0.0748 0.0806 -0.93  0.3536       x1=13     -0.3894 0.1318 -2.96  0.0031       x1=14     -0.2788 0.0953 -2.92  0.0035       x1=15     -0.7672 0.2302 -3.33  0.0009       x1=2      -0.5360 0.2668 -2.01  0.0446       x1=3      -0.3258 0.1548 -2.10  0.0353       x1=4      -0.4092 0.1319 -3.10  0.0019       x1=5      -0.5152 0.2304 -2.24  0.0254       x1=6      -0.2897 0.1538 -1.88  0.0596       x1=7      -0.6216 0.1768 -3.52  0.0004       x1=8      -0.5861 0.1202 -4.88  <0.0001      x1=9      -0.5522 0.1078 -5.13  <0.0001      d2         0.0000 0.0000 -0.64  0.5206       f1        -0.0088 0.0011 -8.19  <0.0001      k8         0.7348 0.0499 14.74  <0.0001       Basically I want to assess the regression in two ways, a) how well the model fits the data and b) how well is the model in predicting the outcome. for a) (goodness of fit): Deviance tests based on Chi-squared are not appropriate in this case because the number of unique covariates approximates N, so we cannot assume a X2 distribution. **Is this interpretation correct?** I can see the covariates using the epiR package. require(epiR) logit.cp <- epi.cp(logit.df[-1]))               id n x1   d2 f1 k8      1 1 13 2030 56  1      2 1 14  445 51  0      3 1 12 1359 51  1      4 1  1 1163 39  0      5 1  7  547 62  0      6 1  5 3721 62  1     ...     7446      I have also read that the Hosmer-Lemeshow gof test is outdated, as it divides the data by 10 in order to run the test, which is rather arbitrary. Instead I use the le Cessie–van Houwelingen–Copas–Hosmer test, implemented in the rms package. I not sure exactly how this test is performed, I have not read the papers about it yet. In any case, the results are:               Sum of squared errors     Expected value|H0                    SD                     Z               1711.6449914          1712.2031888             0.5670868            -0.9843245                          P                  0.3249560       P is large, so there is no sufficient evidence to say that my model doesn´t fit. Great! However.... b) When checking the predictive capacity of the model, I draw ROC curve and get the following results: AUC is 0.6320586 ![enter image description here](http://i.stack.imgur.com/EzYWC.png) That doesn’t look very good. So, to sum up my questions: 1- Are the tests I run appropriate to check my model? What other test could I consider? 2- Do you find the model useful at all, or would you dismiss it based on the relatively poor ROC analysis results?


title: Lognormal and Logit models in JAGS

text: I'm trying to run this two-part model in JAGS. I have two components, one logit and one lognormal regressions in there:               model {          # Lognormal regression model      for (i in 1:N) {         for (j in 1:M) {             y[(i-1)*8+j] ~ dlnorm(mu[(i-1)*8+j], prec.y)             mu[(i-1)*8+j] <- beta0 + delta1[i] + delta2[j] + inprod(x[(i-1)*8+j,], beta[])         }         delta1[i] ~ dnorm(0, tau1)     }          for (j in 1:M) {         delta2[j] ~ dnorm(mu.delt2[j], tau2)     }          #Prior distributions on regression parameters for lognormal model     prec.y ~ dgamma(0.001, 0.001)      sigmasq.y <- 1/prec.y      beta0 ~ dnorm(mbeta0, precbeta0)         for (l in 1:W) {         beta[l] ~ dnorm(m[l], prec[l])     }          # Logistic model     for (i in 1:N) {         for (j in 1:M) {             y0[(i-1)*8+j] ~ dbern(p[(i-1)*8+j])              logit(p[(i-1)*8+j]) <- alpha0 + delta3[i] + inprod(x[(i-1)*8+j,], alpha[])         }         delta3[i] ~ dnorm(mu.delt3[i], tau3)     }          #Prior distributions on regression parameters for logistic model     alpha0 ~ dnorm(malpha0, precalpha0)      for(l in 1:W) {         alpha[l] ~ dnorm(meana[l], preca[l])     }          # Trivariate Normal prior – correlation between random effects      mu.delt2 <- lambda1*delta1 + lambda2*delta3     mu.delt3 <- psi1*delta1     tau1 ~ dgamma(0.001, 0.001)     tau2 ~ dgamma(0.001, 0.001)      tau3 ~ dgamma(0.001, 0.001)           psi1 ~ dnorm(0, 1.0E-2)      lambda1 ~ dnorm(0, 1.0E-2)      lambda2 ~ dnorm(0, 1.0E-2)          sigmasq1 <- 1/prec.delt1      nusq2 <- 1/prec.delt2      nusq3 <- 1/prec.delt3     sigmasq3 <- nusq3 + sigmasq1*pow(psi1, 2)      sigmasq2 <- nusq2 + sigmasq1*pow(lambda1, 2) + sigmasq3*pow(lambda2, 2)      cov13 <- psi1*sigmasq1      cov12 <- (lambda1 + lambda2*psi1)*sigmasq1      cov23 <- psi1*(lambda1 + lambda2*psi1)*sigmasq1 + lambda2*nusq2      corr12 <- cov12/(sqrt(sigmasq1*sigmasq2))      corr13 <- cov13/(sqrt(sigmasq1*sigmasq3))      corr23 <- cov23/(sqrt(sigmasq2*sigmasq3))      }      In R I call               y <- br[,37]     y0 <- br[,50]     x <- as.matrix(br[,c(64,65,66,67,38,13,76,75,44)]) # gdp.cap, bf.cap, ifdm, lula.vote, year     x[is.na(x)] <- 0          N <- length(unique(br$munname))     M <- 8     W <- dim(x)[2]          data <- list(N=N, M=M, W=W, y=y, y0=y0, x=x, m=c(-2,2,4,3,1,0,0,0,0), prec=c(.20,.20,.20,.20,.20,.2,.2,.2,.2), mbeta0=1, precbeta0=.01, malpha0=2, precalpha0=.01, meana=c(-2,2,4,3,1,0,0,0,0), preca=c(.20,.20,.20,.20,.20,.2,.2,.2,.2), prec.delt1=.2, prec.delt2=.2, prec.delt3=.2)          param <- c("alpha", "alpha0", "delta3", "beta0", "beta", "delta0", "delta1")          sim <- jags.model(file="2partmodel.bug", data=data, n.chains=5, n.adapt=1000)      But then I get this error:               Compiling model graph        Resolving undeclared variables        Allocating nodes        Graph Size: 116332          Initializing model     Deleting model          Error in jags.model(file = "2partmodel.bug", data = data, n.chains = 5,  :        Error in node y0[258]     Observed node inconsistent with unobserved parents at initialization      What's wrong?


title: Numerical accuracy of multivariate normal distribution

text: In MATLAB, I've written two snippets of code that compute the PDF of a multivariate normal distribution. However there's a difference in the values these two methods produce and I can't figure out why. I've narrowed the problem down to something having to do with computing the inverse of the covariance matrix. Inaccurate code                 function p = mvnpdf_inacc(X, mu, sigma)         xc = bsxfun(@minus, X, mu);         [n, k] = size(xc);         twopic = (2 * pi) ^ (-k / 2);         sqrtdetsig = sqrt(det(sigma)) ^ -1;         c = twopic * sqrtdetsig;         p = zeros(n, 1);         for i = 1:n             xci = xc(i, :);             p(i) = c * exp(-0.5 * (xci / sigma * xci'));         end     end      Accurate code                 function p = mvnpdf_acc(X, mu, sigma)         [R, err] = cholcov(sigma, 0);              if err             error('%s', 'sigma is not both symmetric and positive definite');         end              X0 = bsxfun(@minus, X, mu) / R;         d = min(size(X));         slogdet = sum(log(diag(R)));         p = exp(-0.5 * sum(X0 .^ 2, 2) - slogdet - 0.5 * d * log(2 * pi));     end      Testing code                 function iseq_func = test_mvnpdf(n)         x = linspace(-2, 2, n);         y = x;         [X, Y] = meshgrid(x, y);         XY = [X(:), Y(:)];         mu = [0, 0];         sigma = [1.0, 0.5; 0.5, 0.4];         p_inacc = mvnpdf_inacc(XY, mu, sigma);         p_acc = mvnpdf_acc(XY, mu, sigma);         p_diff = abs(p_inacc - p_acc);         iseq_func = nnz(p_diff) == 0;     end      I get a value of false from running `iseq_func(25)`. What the heck is going on here? Thanks!


title: Comparison over time of categorical data and between two groups

text: I've got 2 treatment groups with 2 possible outcomes when tested (positive or negative). I also have different time points for the testing. Is there anyway I can analyze if the change from positive to negative is significant between treatment at a given period of time? Thanks!


title: How to display concentration data in space?

text: I have a data frame of chemical concentrations that are measurements taken from 12 locations all in the same vicinity. I have manually assigned x and y coordinates to each location (on a 0 to 20 scale) and each observation has a concentration value. a small sample is given below:               dat<-data.frame(x = rep(c(2, 5.5, 8, 14), each=2),         y = rep(c(3, 9, 14, 17), 2),         z = log(rnorm(8, 10, .75))      I want to display the concentration over space in an image by sort of averaging over the 2 dimensions so as to fill the whole plot space. I think what I want to do here is a contour plot, but I can't figure out how to make this work, so maybe this isn't the right kind of data for such a plot? I'm trying to implement this in R, so any R advice on how to do this would be awesome.


title: Using MGF for multivariate random variables

text: How do you use MGF for solving moment based questions for multivariate random variables? For the single variable case, we: find $E(e^{tX})$, find the interval in which it exists (around 0), differentiate it n times and plug in t as 0 to receive the n'th moment $E(X^n)$. With a general multivariate distribution like $f_{XY}(x,y)$, how do I use MGF to find stuff like $E(X), E(XY), E(Y)$ etc? I read wikipedia and some other pages but they keep losing me at some point. I understand the MGF will be $E(e^{t'X})$ where t and X are now vectors but where do I go from here?


title: Comparing Beta weights in HLM?

text: I want to compare the strengths of different predictors (all at Level 1) using HLM and am not sure what statistical test would be appropriate. More specifically, if I have 3 predictors at level 1 and I want to know which of them are statistically different from the others, what is the best way to do this?


title: Metric and Clustering Method

text: I need some suggestions regarding what kind of metric and clustering analysis I should use. I read a lot of posts but didn't get any hints about this type of data. I have a 3000*5000 matrix, where 3000 is signals from various different sources and 5000 is time point from 0 to 5000 seconds. for example: a signal intensity from a source looks like: `0 0 0 0 0 0 0 0 0 2.3 3.2 4.5 6.5 4.2 3.5 2.5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...` and so on indicating "0" for no signal and a positive rational number for a signal intensity at a particular second. I am interested in capturing similar patterns rather intensity values i.e. having spikes/intensity > 0 at a similar time point. I tried using sq. Euclidean metric and then performed MDS (multi-dimensional scaling) on the distance matrix. When I plotted the graph, I don't discern any similar patterns. I did K-means (k = 3-12 clusters) with the Euclidean metric and tried comparing the signals between the two sources but I didn't get similar patterns between them. What type of metric and clustering algorithm should I use? Would you recommend to me 1-cor metric?


title: How to generate iid samples from the linear congruent method?

text: Given a uniform random number generator (such as the linear congruent method), how shall I generate a sequence of i.i.d. random samples? Are samples generated in a successive sequence i.i.d., or are they just realization of a single random variable (i.i.d. is for several random variables, and therefore doesn't apply to samples of a single random variable?)? Given a seed, we can get a sequence of points. Can these points be taken as iid uniform samples? I don't think they are independent, because each point can be determined by the previous point as $$ X_{n+1} \equiv \left( a X_n + c \right)~~\pmod{m}. $$ If you also don't think they are independent, how do you generate iid samples from the linear congruent method then? If you think they can be good approximation to independence, why? Thanks and regards!


title: Difference between canonical correpondence analysis and canonical correlation analysis

text: I am no statistical geek and this might be a very naive question.....I am sorry before hand... please share your knowledge. I am bit confused between two terms **Canonical Correpondence Analysis** and **Canonical Correlation Analysis**. Are the two some how related or they are entirely different techniques? Do they point to the same outcome or different outcomes (if they are applied to same datasets)? Thank you


title: How to evaluate fit of a logistic regression

text: I have a set of data points, which exhibit a solid linear correlation $r\approx 0.9$. I am basically plotting population in certain areas against the number of occurrences of a certain phenomenon (so in other words, I think the number of people should predict the number of occurrences of this phenomenon). I don't like the linear regression though, because a linear model can take on negative or huge values. In my case the values need to be positive, and there is an upper bound (which I don't know; but logically, the values can only get so large). Unfortunately one of the inputs for which I need a predicted output is **far** larger than the inputs which I used to make a regression. This made me choose a logistic model instead, $$y=\frac{A}{1+Be^{-Cx}}$$ because I am familiar with its shape (always positive, and approaching a limiting value) from basic differential equations, and it seems to be what I want for this situation. Also the data has a slight "point of inflection", so the curvature of the logistic regression (visually to me) is even a little better fit for the data points than the line was. So finally my question(s):   1. I keep reading on the internet that logistic models are meant for probabilities, which take values between $0$ and $1$. I would like, in addition to having a high $r$, to try and reason why this type of model "makes sense" for my situation: am I going the completely wrong direction here?   2. How can I measure the "goodness of fit" of my curve to the data points? I have read Wikipedia but it is too vague for me to understand. An explicit formula - must be computable by hand calculator - would be helpful.  I have done my research but absolutely nothing I have found on the internet has been accessible to me, since I know virtually nothing about stats. An introductory/easy explanation would be so nice... and my math background is far stronger than my stats background, so use all the math you want but assume I'm a beginner at stats. I want it to be somewhat rigorous too - like if I could plug a number into a t-Test for example, that would be good to get an "objective" value.


title: CLT can be used for weighted sum of different Bernoulli variables?

text: Suppose $$ z_i \sim Bernoulli (p_i) $$ Can we use CLT for the following weighted sum? $$ S = \sum_i w_i z_i $$ i.e. can $S$ be approximated with a normal distribution? If yes, with which theorem? (I suppose the classical CLT holds only for average of iid variables)


title: Only binary predictors -- ANOVA, regression, or other?

text: I am trying to fit a model to predict a quantitative response variable, using several binary variables. In particular, I am interested in measuring the relationship between one of the binary variables and the outcome variable, holding the other binary variables constant. My inclination is that ANOVA is the appropriate method to use. The advantage of ANOVA over a series of t-tests should be avoiding Type 1 errors. But I seem to recall from stats that ANOVA is primarily designed for predictor variables with more than 2 levels, aka not binary. Which is preferable? Final question: is there any reason not to use multiple regression in this case? Thanks.


title: Using generalized method of moments (GMM) to calculate logistic regression parameter

text: I want to calculate coefficients to a regression that is very similar to logistic regression (Actually logistic regression with another coefficient: $$ \frac{A}{1 + e^{- (b_0 + b_1 x_1 + b_2 x_2 + \ldots)}},$$ when $A$ could be given). I thought of using GMM to calculate the coefficients, but I'm not sure what are the moment conditions I should use. Can anyone help me with that? Thanks!


title: Are significance tests for the assumption of constant variance too strict when sample size is large?

text: My question is regarding linear regression and non-constant variance. I've heard that even though a large data set fails a normality test, this does not necessarily mean the data is not normal. This is because the tests for normality can be 'overly' sensitive.' Thus, I've heard that, in practice, people overlook when normality fails. Are the tests for constant variance like that? Or should I be able to keep transforming the data until it passes tests for constant variance? Also, I read that the constant variance tests like the breusch pagan test are sensitive to normality. If my data set fails the normality test, then will BP also fail? My data set is over 2500 observations and has many predictor variables. When I take a look at the scatter plots, they look reasonable.


title: Is a biased or unbiased estimator used for pooled SD in calculating Cohen's d?

text: When calculating Cohen's $d$ for independent samples, you must use a pooled $SD$. However, I have seen both of these: $$SD_{\text{pooled1}} = \sqrt{\frac{ (n_1 - 1)s_1^2 + (n_2 - 1)s_2^2}{n_1 + n_2}}$$ vs. $$SD_{\text{pooled2}} = \sqrt{\frac{ (n_1 - 1)s_1^2 + (n_2 - 1)s_2^2}{n_1 + n_2 -2}}$$ Supporting the use of $SD_{\text{pooled1}}$ , some website have $SD_{\text{pooled}}$ listed as $\sqrt{\frac{s_1^2 + s_2^2}{2}}$, which is the same as $SD_{\text{pooled1}}$ when sample sizes are equal. From online discussions, it sounds like $SD_{\text{pooled1}}$ is a "biased" estimator of $SD$, and $SD_{\text{pooled2}}$ is less biased, that is, I think it means $SD_{\text{pooled1}}$ underestimates $SD$. In addition, some sites suggest that some effect size metrics use $SD_{\text{pooled1}}$ (Cohen's $d$), and others $SD_{\text{pooled2}}$ (Hedges' $g$). Is this true? And if so, why would one effect size metric (Hedges' $g$) use the unbaised estimator, $SD_{\text{pooled2}}$ , while the other (Cohen's $d$) not?


title: Correlations between continuous and categorical (nominal) variables

text: I would like to find the correlation between a continuous (dependent variable) and a categorical (nominal: gender, independent variable) variable. Continuous data is not normally distributed. Before, I had computed it using the Spearman's $\rho$. However, I have been told that it is not right. While searching on the internet, I found that the boxplot can provide an idea about how much they are associated; however, I was looking for a quantified value such as Pearson's product moment coefficient or Spearman's $\rho$. Can you please help me on how to do this? Or, inform on which method would be appropriate? Would Point Biserial Coefficient be the right option?


title: Finding P(A|Ū) from known P(A|U), P(A) and P(U)

text: I am wondering I may get some advice. If I wanted to find the value of P(A|Ū), and have values for P(A|U), P(A) and P(U), is there any probability functions that I can use? Many thanks once again.


title: Dealing with correlating fixed effects in a linear mixed-effects analysis

text: My question is about the best way to estimate the effect of a predictor on a dependent variable, while accounting for several other predictors that may correlate with the predictor of interest. I'm using a linear mixed-effects model, using the `lmer` function from the R `lme4` package. (Warning: I'm fairly new at this, so their may be some misunderstandings woven through my question.) ## The problem To make things a bit more specific, I'll just explain the actual data that I'm working with. I have eye-movement data of participants freely viewing natural scenes. I want to determine whether pupil size predicts the 'visual saliency' (i.e. the conspicuity) of the locations in the image that participants are looking at. But there are many other things that correlate with pupil size, such as luminosity, and this makes the analysis tricky (or does it?). ## Option 1 (simple): Looking at fixed effects One option would be to simply create a linear mixed-effects model that has all predictors of saliency that I can think of, including the predictor of interest (`pupil_size`), as fixed effects and `subject` and `scene` as random effects. (To keep things manageable, I'm using a purely additive model, although I suppose that this is a whole topic in itself.)               my_lmer = lmer(saliency ~ brightness + (.. lots of predictors ...)         + pupil_size + (1|subject) + (1|scene))      This will give me a t-value for the fixed effect `pupil_size`. From what I understand, this fixed effect will already be partial, so it's the unique predictive power of pupil size, with any correlations between fixed effects already taken into account. Is my understanding correct? ## Option 2 (complex): Using model comparison An alternative approach, which I have from Baayen et al. (2008), is to compare a model without pupil size as fixed effect (`simple_model`) to a model with pupil size as fixed effect (`complex_model`).               simple_model = lmer(saliency ~ brightness + (.. lots of predictors ...)         + (1|subject) + (1|scene))     complex_model = lmer(saliency ~ brightness + (.. lots of predictors ...)         pupil_size + (1|subject) + (1|scene))      Now I can use the `anova` function to compare these two models (see Baayen's paper for an example). This will give me a `Chisq Chi` value, and I can use this to determine whether adding `pupil_size` as fixed effect is a justified addition to the model. Clearly, this model comparison approach is more complex than simply looking at the t-values for fixed effects in a single model. And it seems to me that if `pupil_size` is a significant predictor (per Option 1), then it must also be a significant addition to the model (per Option 2). In sum, my question is: _Is there any reason to do a model comparison (Option 2), or am I better off just creating a single linear mixed-effects model and seeing whether the t-value associated with`pupil_size` as fixed effect is sufficiently high (Option 1)?_


title: Four tricky time series questions with a "seasonal twist"

text: A ski-hotel has the most guests in the third quarter in every year (check the data below after the four questions).   Can you answer these four questions (every year has 4 values, the first is quarter one and the second is quarter two, the third is quarter three and the fourth is quarter four and it's from the 2007 - 2012):   1. Presume that there is no autocorrelation. what is the proper theoretical model for this type of data?   2. Validate the model assumptions.   3. Describe the estimation of the theoretical model in question number 1 above.   4. Make a prediction of the occupancy on the hotel for season 3, year 2014.               Year         2007:  250 260 370 300        2008:  270 280 375 300        2009:  280 285 390 320        2010:  290 320 395 320        2011:  295 325 430 340        2012:  320 360 450 340


title: Reporting Mann-Whitney U-Test without homogenity of variance

text: I am really stuck with this question: The Mann-Whitney test requires homogenity of variance if a median difference is suppossed to be statistically significant. In case homogenity of variance is not met, but the test is significant: Which aspects of the test can I report? Thanks a lot in advance.


title: How to calculate the MGF for a poisson random variable

text: $\phi (t)=\sum_{n=0}^{\infty}e^{tn}\cdot \frac{\lambda^{n}}{n!}e^{-\lambda}$ $=e^{-\lambda}\sum_{n=0}^{\infty}\frac{(e^{t}\lambda)^{n}}{n!}$ $=e^{-\lambda+\lambda e^{t}}$ I'm confused here as how $\sum_{n=0}^{\infty}\frac{(e^{t}\lambda)^{n}}{n!}$ simplifies to $e^{\lambda e^{t}}$, can someone please explain?


title: Is the goodness of fit test in JMP the Hosmer-Lemeshow goodness of fit test?

text: I'm working with an organization that is using JMP in their analysis, and I can't tell from the description in JMP's help files if the test for goodness of fit in their logistic regression is the Hosmer-Lemeshow test. If it makes a difference, my data set has only one predictor variable, so we aren't considering complex models. > The next questions that JMP addresses are whether there is enough > information using the variables in the current model or whether more complex > terms need to be added. The Lack of Fit test, sometimes called a Goodness of > Fit test, provides this information. It calculates a pure-error negative > log-likelihood by constructing categories for every combination of the > regressor values in the data (Saturated line in the Lack Of Fit table), and > it tests whether this log-likelihood is significantly better than the Fitted > model. > > The Saturated degrees of freedom is m–1, where m is the number of unique > populations. The Fitted degrees of freedom is the number of parameters not > including the intercept. For the Ingots example, these are 18 and 2 DF, > respectively. The Lack of Fit DF is the difference between the Saturated and > Fitted models, in this case 18–2=16. > > The Lack of Fit table lists the negative log-likelihood for error due to > Lack of Fit, error in a Saturated model (pure error), and the total error in > the Fitted model. Chi-square statistics test for lack of fit. (I'm an engineer, not a statistician, the saturated line parts are confusing me).


title: Complement of Conditional Probability

text: I am wondering I may get some advice if P(A|Ū) is equals to 1- P(A|U)? Many thanks once again.


title: determining the order of Box-Jenkins modeling process

text: I have this problem on what model class (AR,MA,ARMA,ARIMA,etc) will I use on my data I'm using Box-Jenkins process and what order( say 1,0,1) will I use. I already done many transformations on my data but the errors are so large and the correlation is somewhat small. My data are stationary using (ADF test and KPSS test) but not normally distributed (Anderson-Darling, Wilk-Shapiro and Kolmogorov- test). So I apply natural log and then test it again but is not normally distributed. So I differenced it by one and it is now stationary and normally distributed. I already satisfy the requirements of using a Box- Jenkins process then I use auto.arima in R to know what order to use and I also try SPSS using its expert modeler to cross check but my problem is I still get large errors and small R-squared. I need to know what to do for determining order? I also have problems in understanding ACFs and PACFs. this is my actual data               Harvest      > 60477 29323 51369 15800 58994 45496 17227 92103 138573 39181 51192 13132 400 > 18258 54553 7220 1418 6807 17915 89015 122154 122853 63398 27246 27013 36317 > 65735 94744 78763 39769 20422 27398 33552 10000 6500 5300 5700 4800 5300 > 6450 9300 5834 29200 39975 65000 45494 79000 7900 54758 70581 31505 45437 > 29691 110947 40498 71238 42170 38723 64813 122992 17929 11652 134137 110043 > 60153 7625 25967 38918 1621 14946 76610 84516 72223 40399 63482 34918 63098 > 105388 135809 31345 66880 160511 40238 35767 105560 119276 154348 86935 > 73728 167119 128709 97040 21780 9906 62213 99940 72626 117783 58037 68756 > 25721 19853 4943 2027 20251 114718 27801 80868 94761 18914 119632 187924 > 56950 52886 141456 141507 ![Harvest graph](http://i.stack.imgur.com/ASQah.png) This is the differenced data               d_Harvest      > -31154 22046 -35569 43194 -13498 -28269 74876 46470 -99392 12011 -38060 > -12732 17858 36295 -47333 -5802 5389 11108 71100 33139 699 -59455 -36152 > -233 9304 29418 29009 -15981 -38994 -19347 6976 6154 -23552 -3500 -1200 400 > -900 500 1150 2850 -3466 23366 10775 25025 -19506 33506 -71100 46858 15823 > -39076 13932 -15746 81256 -70449 30740 -29068 -3447 26090 58179 -105063 > -6277 122485 -24094 -49890 -52528 18342 12951 -37297 13325 61664 7906 -12293 > -31824 23083 -28564 28180 42290 30421 -104464 35535 93631 -120273 -4471 > 69793 13716 35072 -67413 -13207 93391 -38410 -31669 -75260 -11874 52307 > 37727 -27314 45157 -59746 10719 -43035 -5868 -14910 -2916 18224 94467 -86917 > 53067 13893 -75847 100718 68292 -130974 -4064 88570 51 ![enter image description here](http://i.stack.imgur.com/WzyDi.png)


title: Probability of Rolling a number that changes over time on One Die

text: I've seen the following question asked and answered many times on various forums: > Knowing that the odds of rolling a six with one die are one in six, what are > the odds of rolling a six in two, three, four (etc.) consecutive rolls? 1/6 > + 1/6 + ... wouldn't make sense as you'd have a 100% chance of rolling a six > in six rolls. The answer I've seen is that the odds of rolling a six after n rolls are 1 - (5/6)^n. Now my question is this - lets say that instead of rolling for a six, we pick a number out of a fixed sequence of six unique numbers after each roll. So for example, lets say our sequence is 1, 2, 3, 4, 5, 6. At the first roll we're looking to roll a one, the second roll we're looking for a two, then a three etc. until finally at the seventh roll we go back to one. What are the chances of matching any of the numbers in the sequence in this way after n rolls?


title: What is an "uninformative prior"? Can we ever have one with truly no information?

text: Inspired by a comment from this question: What do we consider "uninformative" in a prior - and what information is still contained in a supposedly uninformative prior? I generally see the prior in an analysis where it's either a frequentist-type analysis trying to borrow some nice parts from Bayesian analysis (be it some easier interpretation all the way to 'its the hot thing to do'), the specified prior is a uniform distribution across the bounds of the effect measure, centered on 0. But even that asserts _a_ shape to the prior - it just happens to be flat. Is there a better uninformative prior to use?


title: Zero trick in WinBUGS for truncation

text: I am using the Zero-trick to create a prior for correlation parameter. Data `y` is a bi-variate normal with covariance matrix with same diagonal elements and correlation phi. The covariance matrix is positive definite for `|phi| < 1`, the usual range for correlation. But instead of common prior `phi ~ uniform(-1,1)`, I am using               phi ~ dnorm(0,1)     O <- 0     O ~ dbern(cons)     cons <- step(phi*phi-1)      The zero trick should force `|phi| < 1` but the Gibbs sampler picks up illegitimate values for `phi` and the code crashes. I have a more complicated model than this where I can't use uniform priors. Here is the code if you want to check yourself.               model     {     for (i in 1:n) {         y[i,1:2] ~ dmnorm(mu[1:2], Sinv[,])        }           # Priors      for (i in 1:2) {mu[i] ~ dnorm(0, 0.01) }          # Get the precision matrix Sinv from the var-cov matrix     Sinv[1:2,1:2] <- inverse(Sig[1:2,1:2])           # Inverse Gamma prior for variance parameter sigma      tau ~ dgamma(1,1)      sigma <- 1/tau          # Build the var-cov matrix Sig from variance and correlation (phi)      Sig[1,1] <- sigma      Sig[2,2] <- sigma      Sig[1,2] <- sigma*phi      Sig[2,1] <- Sig[1,2]          # Constrained priors for phi using the "Zero-trick"     # phi ~ dunif(-1,1) # We could use this but want to see if zero-trick works          phi ~ dnorm(0,1)     O <- 0     O ~ dbern(cons)     cons <- step(phi*phi-1)          }          # Initial values     list(mu = c(0,0),     tau=1,phi=0)          # This dataset causes OpenBUGS code to crash     list(n=5,     y = structure(.Data = c(     5.8, 2.7,     1.3, 5.2,     9.1, 7.9,     8.5, 7.6,     2.7, 5.2),.Dim = c(5,2)))


title: Estimating and using reliability metric with the same data

text: I'd like to know if it is reasonable to use the same set of data to estimate and use a reliability metric from repeated measures? I am mostly interesting by the minimal detectable change (MDC). But I thought that the MDC was a property of a measuring instrument BEFORE using that instrument. Is that wrong? Is it cheating to estimate the MDC and then use it to determine if a change has really happened?


title: `bGAMM` and other `GMMboost` algorithms for large data sets

text: Regularized generalized linear mixed models and generalized additive mixed models are exactly what I need. I'm an R user, so it looks like `bGAMM` and maybe `glmmLASSO` are possibilities outside of doing my own implementation (which I don't have time to do for this project). My question is whether these packages can handle reasonably large numbers of observations (up to 150K) and random effects (over 500) along with, say, up to 1000 parameters. I'll be using a beefy Amazon instance for this use case.


title: Textbooks on the numerical solution of stochastic differential equations

text: Is there a canonical textbook or textbooks on the use of numerical methods in solving stochastic differential equations that you would recommend?


title: Sample distribution for Kullback-Leibler distance

text: For two $n$ dimensional multivariate normal distributions $X_{1}\sim N\left(\mu_{1},\Sigma_{1}\right)$ and $X_{2}\sim N\left(\mu_{2},\Sigma_{2}\right)$, the Kullback-Leibler distance is given by $$KL=\frac{1}{2}\left(tr\left(\Sigma_{2}^{-1}\Sigma_{1}\right)+\left(\mu_{2}-\mu_{1}\right)'\Sigma_{1}^{-1}\left(\mu_{2}-\mu_{1}\right)-ln\left(\frac{det\left(\Sigma_{1}\right)}{det\left(\Sigma_{2}\right)}\right)-n\right) $$ When $\Sigma_{1}=\Sigma_{2}$, the above will be similar to half the Mahalanobis distance squared. The Mahalanobis distance squared is chi-square distributed with $n-1$ degrees of freedom, so I imagine that is part of the answer. When the parameters of $KL$ are replaced with estimates from on a random sample, is it possible to derive the distribution in the more general case when $\Sigma_{1}$ is not necessarily equal to $\Sigma_{2}$?


title: Time series modeling the number of users of a mobile app

text: I want to model the number of users of an mobile app. This app has two kinds of users: _free_ and _paid_. I thought of this autoregressive model: $x_t = Ax_{t-1}$ with $x_t$ being a 4-dimensional vector, with quantities (all in period $t$) for the numbers of   * mobile users   * visitors of the app page   * free users   * paid users I imagined $A$ as follows: $A = \begin{bmatrix} a & 0 & b & c \\ d & e & 0 & 0 \\ 0 & 1-e & f & 0 \\ 0 & 0 & 1-(b+f) & g \end{bmatrix}$ Each parameter is the proportion of people in a group going to another group. Example: $f$ is the proportion of _free users_ that remain _free users_ in the next period; $d$ is the number of mobile users that visited the app page, and so on. Problems I see with this model:   1. $b$ and $c$ are the proportion of users that leave the app. Not sure if it's nice to put them in the group of _mobile users_ , or make another dimension for the group _former users_ , that don't register back in the app.   2. The number of mobile users grow in time, so the way I thought to cope with this is by not demanding that $a$ and $d$ sum is 1. Not sure if it's valid. So, the question is: how can I make a better model that I could feed the numbers after the launch of the app and get good estimates for these parameters, and good estimates for future numbers of users?


title: Prediction based on repeated measure with binary outcome

text: I'd like to make a prediction about a subject's likelihood that a certain outcome will occur based on a series of measurements taken over time. The easiest way to explain the problem is as follows. Imagine that I'm tracking an individual over time to generate a probability of ever getting lung cancer. The training data would look something like this with the outcome being "has cancer"               ID           age   cigarettes  asbestos   lung        has cancer                          per day   exposure   capacity     S_1           5       0           N          20             N     S_1           17      20          N          50             N     S_1           31      50          N          30             N     S_1           42      0           N          40             Y     S_2           15      0           N          55             N     S_2           37      0           N          45             N     S_3           5       3           Y          15             N     S_3           50      0           Y          30             Y     etc...      The test data would be taken before an individual has gotten cancer and would predict their probability of getting cancer at any time in the future. What's the best way to model this? I'm currently experimenting with mixed models, treating this like a repeated measure problem. However, the dependent variable (whether someone ever gets cancer) is not being repeatedly measured over time. Instead it is fixed, so it seemed this might not be the best approach. As you can see there are additional problems of the dataset being unbalanced, with measures taken at irregular periods. **EDIT 1:** Would it be better to think of this in terms of survival analysis? I assume the event in this case would be contracting cancer. If so how do I deal with deaths not related to lung cancer? **Edit 2:** I'm thinking now of using a joint model that combines a mixed effect model with a survival analysis model: http://rwiki.sciviews.org/doku.php?id=packages:cran:jm. Any thoughts as to whether this might work for my problem? It seemed that I could use the hazard function to get a prediction for likelihood of ever getting cancer.


title: Mixed effect model

text: I have a model: $y=a+bx+\varepsilon_1$ for each observation (tree) $b=c+dz+\varepsilon_2$ for each forest (group of trees) $y$, $x$, $z$ are continuous variables. $a$, $b$, $c$, $d$ are coefficients. $z$ is a property of forests (not individual trees). Error $\varepsilon_1$ is normally distributed and varies by tree. Error $\varepsilon_2$ varies by forest and not much is known about its distribution. Is this lme4 model specification correct?               lmer(y~x+x:z+(0+x|forest))      I will greatly appreciate your help!


title: Sum of autoregressive processes?

text: I am working on a research topic where I need to add together two AR processes and I was wondering if the distribution of these processes is of a recognizable form/structure. More formally, if $x_t$ is a AR(p) process with characteristic polynomial $\Phi_x(u)$ and $y_t$ is a AR(q) process with characteristic polynomial $\Phi_y(u)$, then what is the structure of $z_t=x_t+y_t$?


title: Zero-inflated negative binomial mixed-effects model in R

text: Is there such a package that provides for zero-inflated negative binomial mixed-effects model estimation in R? By that I mean:   * Zero-inflation where you can specify the binomial model for zero inflation, like in function zeroinfl in package pscl:               zeroinfl(y~X|Z, dist = "negbin") where Z is the formula for the zero inflation model;   * Negative binomial distribution for the count part of the model;   * Random effects specified similar to function lmer of package lme4. I understand glmmADMB can do all that, except the formula for zero inflation cannot be specified (it is just an intercept, i.e. Z is just 1). But are there any other packages that can do it all? I will be very thankful for your help!


title: Modification of "corrected repeated k-fold cv test” when also averaging Random Forest results across multiple sampling seeds?

text: I would be _very_ grateful for any ideas concerning the following problem. I would be even more grateful if someone could point me to a literature reference proposing a solution for a scenario like the one I describe below. In short, I am looking to adapt Bouackert and Frank's 2004 article on corrected repeated k-fold cv test (PDF version) for a scenario in which the classifier is also intrinsically random. I am a student in cheminformatics (i.e. would appreciate mathematical jargon kept to a minimum if possible). The scenario:   1. I am comparing two descriptor/attribute/feature sets (for the purpose of building classifiers using Machine Learning - the classes might be toxic vs. nontoxic). For example - encoding molecules (instances) in my dataset by a vector of numbers corresponding to a) a 2D vs. b) a 3D representation of molecular structure.    2. In both cases, I generate k-fold cross-validation results (non-overlapping folds) using Random Forest acros R repetitions of cross-validation. Furthermore, for each train/validation set pair (corresponding to a given value of R and k), I repeat model generation and validation Q times - each time, using a different seed to initialise the random number generator for Random Forest.  E.g. for a given value of k and R (I am working in the R programming language):              set.seed(CURRENT_SEED_VALUE)      rf <- randomForest(current_train_x,current_train_y)     predictions <- predict(rf,current_test_x)     figure_of_merit <- fom_function(predictions,current_test_y)        3. I can then take the arithmetic mean of all figure_of_merit values (for both sets of descriptors). I wish to get a p-value corresponding to the null-hypothesis: The mean difference in this figure of merit for both methods is zero. Supposing I was not considering multiple randomForest RNG seeds, I could obtain such a p-value using Bouckaert and Frank's 2004 “corrected repeated k-fold cv test”. But, how could I modify this to get a p-value for the overall mean difference??? (The only approach I can currently think of is getting one p-value per RNG seed value, then (after correcting for multiple testing using `p.adjust(...)`) seeing whether the p-value is statistically significant for all RNG seeds.)


title: Logistic regression performs better on validation data

text: Recently I've been building a model using logistic regression. To my suprisise LIFT chart looks better on the validation data than on the training data, the same is with ROC. All variables in the model are statistically significant. The question is: Is this really a serious problem? If so, what are the methods to detect what cause it? Does this problem have a name? I've tried to solve it by building all possible models without a single variable, but the problem still occured in each model. **Additional Info**   The size of the sample is about 20000, the proportions are: Training 70%, Validation 30%, I don't have results in front of my eyes, but as I recall the lift in second percentil in training data was about 8 and in validation was about 9, the a priori is 4.2%.


title: Correlation of two samples from single subject?

text: I have a small data set which contains two tissue samples from ~20 patients. Each patient has contributed a tissue sample from left and right hip joint. My aim is to study certain cell amounts in these samples. My aim is to investigate whether findings of the hips in each patients correlates, ie. that both hips have similar cell counts (=there is only a minimal within-subject variation). This can be simply done creating two variables, Var1 and Var2. Var1 is the cell count of right hip and Var2 left hip. If there is significant correlation, it is shown that cell counts in same patient does not vary in random. However I have two covariates. First one, Cov1, is a blood marker, ie. each patient contributes one value, ie. there is no within-subject variation. The other variable (Cov2) is hip related variable, (time interval between primary surgery operation and point where the tissue sample was retrived), ie. there is within-subject variation. How can i investigate the correlation between Var1 and Var2 and adjust with Cov 1 and/or Cov2? Thanks!


title: how to calculate precision, recall, and accuracy, if the label is not binary?

text: I am new in statistics, so maybe my question kind of stupid.   I know how to evaluate the prediction performance if the label is binary, such as the classification label is `sick(yes/no)`.   The situation is the label is `sick(level1/level2/level3)`,in this situation, could we still calculate the precision? and how?   I feel confused in how to calculate TN/TP/FN/FP. I hope anyone can help me. And I understand R code, if needed, you can write in R or use any R package. For a better descript of data, it is kind like               label  x1   x2     1      21   123     1      223  22     2      22   1     3      8    9


title: a question on multiplicative SMV kernel

text: I'm new to SVM and would like to use it to solve a problem formulated, for example, as follows: The patterns are four-dimensional vectors $(x_1,x_2,x_3,x_4)$, and the kernel is $K(\bf{x_i}, \bf{x_j})=k1([{x_i}_1, {x_i}_2], [{x_j}_1, {x_j}_2]) \times k2([{x_i}_3, {x_i}_4], [{x_j}_3, {x_j}_4)$ In other words, the kernel is the multiplication of two other kernels, each of which takes a pair of sub-patterns as input. I came across libSVM and would like to know if the library is able to solve this kind of problem. Thanks.


title: Are all values within a 95% confidence interval equally likely?

text: I have found discordant information on the question: " **If one constructs a 95% confidence interval (CI) of a difference in means or a difference in proportions, are all values within the CI equally likely?** Or, is the point estimate the most likely, with values near the "tails" of the CI less likely than those in the middle of the CI? For instance, if a randomized clinical trial report states that the relative risk of mortality with a particular treatment is 1.06 (95% CI 0.96 to 1.18), is the likelihood of 0.96 being the correct value the same as 1.06? I found many references to this concept online, but the following two examples reflect the uncertainty therein:   1. Lisa Sullivan's module about Confidence Intervals states: > The confidence intervals for the difference in means provide a range of > likely values for ($μ_1-μ_2$). It is important to note that all values in > the confidence interval are equally likely estimates of the true value of > ($μ_1-μ_2$).   2. This blogpost, titled Within the Margin of Error, states:  > What I have in mind is misunderstanding about “margin of error” that treats > all points within the confidence interval as equally likely, as if the > central limit theorem implied a bounded uniform distribution instead of a > _t_ distribution. [...]   >  The thing that talk about “margin of error” misses is that possibilities > that are close to the point estimate are much more likely than possibilities > that are at the edge of the margin". These seem contradictory, so which is correct?


title: Why is correlation between y and $\hat{y}$ in a model with and without intercept equal?

text: I consider a very simple example with two models like this;               set.seed(66)     n   <- 100     x   <- rnorm(n,0,1)     y   <- 10 + 2*x + rnorm(n,0,1)          mod1 <- lm(y ~ x)     mod2 <- lm(y ~ -1 + x)      I already checked the following: If I write down $r_{y\hat{y}}$ in Matrix notation with $P_j = j(j'j)^{-1}j'$ and $P_x = x(x'x)^{-1}x'$ for the model without intercept and $P_X = X(X'X)^{-1}X'$ for the model with intercept ill gain $r_{y\hat{y}}$ by $r_{y\hat{y}} = \frac{ (y'M_j\hat{y})^2 }{y'M_jy\hat{y}M_j\hat{y}}$ Here I just exchange $\hat{y}$ with $P_xy$ for the model without intercept and $P_Xy$ for the model containing the intercept. I already now that stuff like FWL does not apply here. Hence $P_x y \ne P_X y$. I already checked via simulation and found that the nominator and denominator for $r_{y\hat{y}}$ is different for both models. But apparently the difference is just a proportional amount which cancels out if I do the fraction. What amount is that? I don't get why? What am i missing? I tried already a lot of different approaches. For instance I tried to write $r_{y\hat{y}}$ for the model with intercept as a function of $r_{y\hat{y}}$ from the model without intercept but I could not figure it out quite right...


title: Sum of Squares reduction test: Convergence criteria met, but not all parameters of the model estimated

text: Background: I am using weighted non-linear regression to model the growth of plant organs, with dummy variables for different species. I am using a sum of squares reduction test (SSRT) to compare the parameter estimates of the model between the species in my study. For instance, for this logistic equation: w=(Wo*Wf)/(Wo + (Wf-Wo)*exp(-k*Time)) I want to test if the Wf parameter (asymptotic organ weight) for Species A = Species B, and so on for the rest of the parameters. In my study I have 5 different species, and for a model with 3 parameters, comparing 2 species at a time, I would need 30 SSRTs to make all the relevant comparisons (plus the full model). Here are the SSRT F-tests that I perform to make each comparison (as presented in p.5167 of the SAS 9.3 user's manual): Fr= ((SSEr-SSEf)/q)/(SSEf/(n-p)) Where SSEr and SSEf are the error sums of squares of the reduced and full models, respectively. q is the number of constraints on the full model to get the reduced model (in my case, setting the values of a parameter equal for two species), n is the number of observations and p is the number of parameters in the full model. I have the code in SAS to do all of this, but I find that in some cases while the convergence criteria is met, SAS cannot estimate all the parameters in the model. Hence, I get that instead of reducing the model by 1 parameter (q=1, since I am making a given parameter equal for two species), I get that the model is actually reduced by 2 or 3 parameters (q=2 or 3). In most (if not all) cases, this makes sense, since some of the species have quite different parameter estimates (e.g. some have very different asymptotic organ weights), and when I make a given parameter equal for two species, I also affect the estimation of the other parameters in the model. **Question:** Can I still use these SSRTs that converge, even though I cannot obtain adequate estimates for all the parameters in the reduced models? I am concerned that the SSEr in these situations may not be reliable (and thus, the Ftests), since SAS could not find estimates for all the parameters. Again, I am not interested in the parameter estimates from the SSRTs, I just want to know if the species are different for individual parameters or not. _As an aside, I should mention that I have tried everything from reparameterizing the model, transforming variables/data, different numerical methods, parameter starting values, bounds... and still nothing. I did not change the standard convergence criteria either_


title: ARIMA Forecast Restrictions

text: Is it possible to fit an ARIMA under a restriction that bounds the minimum to positive values? I'm using Minitab.


title: Classification optimal decisions considering a loss function

text: > Suppose we're given data from three different classes which are normally > distributed with the following means and variances: > > $C_1: \mu_1=(1,2)^T, \Sigma_1^{-1}=( \begin{array}{ccc}2 & 1 \\1 & 2 > \end{array})$   >  $C_2: \mu_2=(2,-2)^T, \Sigma_2^{-1}=( \begin{array}{ccc}1 & 0 \\0 & 2 > \end{array})$   >  $C_3: \mu_3=(1,-1)^T, \Sigma_3^{-1}=( \begin{array}{ccc}7 & 5 \\5 & 6 > \end{array})$ > > And the loss function of those three classes is $L=\left > (\begin{array}{ccc}0 & 1 &3 \\2 & 0 & 2\\ 4&3&0\end{array}\right)$ > > How should criterion be modeled to make optimal decisions for such a > problem? And to which class the point $(0,0)^T$ must be assigned? I think I can solve the problem without considering the loss function (from prior and posterior probabilities) but I have no idea what can be done about the loss function and consider it in making optimal decisions.


title: Real examples of multinomial distribution

text: I found multinomial distribution really powerful in theory. But I don't know if my real case can be treated as multinomial. There are two cases:   1. There are $n$ students, $30\%$ chose apple, $40\%$ chose banana, $30\%$ chose grape. Then can I say the choice of single student follow a $\mathrm{Multinomial}(n,0.3,0.4,0.3)$ distribution?   2. Similar, say I have $n$ books, $30\%$ will be distributed to group 1, $40\%$to group 2, the rest to group 3. Then can I assume the allocation of a single book follows a multinomial distribution? What if I was to allocate money instead of books? In this case, money is continuous not discrete.


title: Decision theory - reject option

text: In decision theory, we define a reject option ($\theta$) so that when making decision is difficult, the case will be ignored. Suppose $1/k \leq \theta \leq1$:   * If $\theta=1/k$ no cases will be rejected and   * If $\theta=1$ all cases will be rejected. Why this happes and how can I prove it for $\theta=1$ and $1/k$? Section 1.5.3, chapter 1 of Pattern recognition by Bishop: (reject option) > We have seen that classification errors arise from the regions of input > space where the largest of the posterior probabilities p(Ck|x) is > significantly less than unity, or equivalently where the joint distributions > p(x, Ck) have comparable values. These are the regions where we are > relatively uncertain about class membership. In some applications, it will > be appropriate to avoid making decisions on the difficult cases in > anticipation of a lower error rate on those examples for which a > classification decision is made. This is known as the reject option. For > example, in our hypothetical medical illustration, it may be appropriate to > use an automatic system to classify those X-ray images for which there is > little doubt as to the correct class, while leaving a human expert to > classify the more ambiguous cases. We can achieve this by introducing a > threshold θ and rejecting those inputs x for which the largest of the > posterior probabilities p(Ck|x) is less than or equal to θ. This is > illustrated for the case of two classes, and a single continuous input > variable x, in Figure 1.26. **Note that setting θ = 1 will ensure that all > examples are rejected, whereas if there are K classes then setting θ < 1/K > will ensure that no examples are rejected. Thus the fraction of examples > that get rejected is controlled by the value of θ.**


title: How to nest in a generalized linear model in R?

text: I have 108 counts taken from 36 field points. I took three counts from each point at different times, also, the 36 field points are divided in a altitude factor of three levels (12 points at each level). I have also six numerical independent variables measured at each count. I want to make a generalized linear model to see which of my independent variables contribute significantly to the counts observed (time is a factor of three levels). I know that if I had only factors I would use a split-plot design to see differences among levels, where I know how to nest, but I have six numerical independent variables. Transform data to do a lm is not working, that´s why I am trying to use a `glm()`. My question is, how can I nest my data in a `glm()` function?


title: Reducing size of PMML DataDictionary and MiningSchema elements

text: The DataDictionary and MiningSchema elements for a PMML model requires quite a bit of metadata for each field. With sparse, high dimensional data these could each be many times larger than either the training data or the trained model. Are there any conventional extensions (or evil non-standard kludges, depending on how you think about it) to the PMML syntax that, for instance, just says that all fields have the same metadata? Or, even better, do something like specify all fields whose names have the same prefix get the same metadata? Also, any typical way of having the MiningSchema just say "use everything in the DataDictionary that's not the "predicted" feature as an input feature"?


title: Why not perform meta-analysis on partially simulated data?

text: ### Background: A typical meta-analysis in psychology might seek to model the correlation between two variables X and Y. The analysis would typically involve obtaining a set of relevant correlations from the literature along with sample sizes. Formulas can then be applied to compute a weighted average correlation. Then, analyses can be performed to see whether correlations vary across studies by more than would be implied by mere effects of random sampling. Furthermore, analyses can be made a lot more complex. Estimates can be adjusted for reliability, range restriction, and more. Correlations can be used in combination to explore meta structural equation modelling or meta regression , and so on. However, all these analyses are performed using summary statistics (e.g., correlations, odds ratios, standardised mean differences) as the input data. This requires the use of special formulas and procedures which accept summary statistics. ### Alternative approach to meta-analysis Thus, I was thinking about an alternative approach to meta-analysis, where raw data is used as input. I.e., for a correlation the input data would be the raw data used to form the correlation. Obviously, in most meta-analyses several if not most of the actual raw data is not available. Thus, a **basic procedure might look like this:**   1. Contact all published authors seeking raw data, and if provided, **use actual raw data**.   2. For authors that do not provide raw data, **simulate raw data** so that it has identical summary statistics as those reported. Such simulations could also incorporate any knowledge gained from the raw data (e.g., if a variable is known to be skewed, etc.). It seems to me that such an approach might have several benefits:   * Statistical tools that use raw data as input could be used for analyses   * By at least obtaining some actual raw data, authors of meta-analyses would be forced to consider issues related to the actual data (e.g., outliers, distributions, etc.). ### Question   * Are there any problems with performing meta-analysis studies on a combination of true raw data and data simulated to have identical summary statistics to existing published studies?   * Would such an approach be superior to existing methods of performing meta-analyses on summary statistics?   * Is there any existing literature discussing, advocating, or critiquing this approach?


title: What is the best way to use a 2-class classifier for a multi-category case?

text: I have a 3 class sample labeled data set, which I have divided into 2 parts. I am using the first part to train the two-class perceptron classifier. One approach is to train $\binom{3}{2}$ two-class classifiers on training data for two classes at a time, and then use voting to classify the test samples for a multi-category case. The obvious problem with the above approach is the presence of ambiguous regions like the ones shown in the image below: ![Linear decision boundary for a 4-class problem](http://i.stack.imgur.com/dP67S.png) Pattern classification by Duda, Hart, Stork suggests training $c$ different linear discriminant functions, where $c$ is the number of unique classes, such that $$g_{i}(x) = W_{i}^tX + w_{i0} \quad \quad \quad i = 1,..., c$$ and assigning $X$ to $\omega_{i}$ if $g_{i}(X) > g_{j}(X)$. The resulting classifier is called in the text book as a _linear machine_. Following is an illustration from the book showing decision boundary produced by a linear machine for 3-class problem. ![decision boundary produced by a linear machine for 3-class problem](http://i.stack.imgur.com/JYlNC.png) My doubt is how is the training process of the $c$ linear discriminant functions $g_{i}(X)$ different from training of $\binom{c}{2}$ 2-class classifiers? **UPDATE:** The one vs rest approach mentioned in DavidDLewis's answer is not the one I am referring to. With the one vs rest approach, there are infact even more ambiguous regions. I am referring to linear machine classifier using which there are no ambiguous regions. See illustration of 1 vs rest approach below: ![1 vs rest](http://i.stack.imgur.com/qUP9Y.png)


title: How to understand the plotting of the cox.zph function in R?

text: I use a Cox proportional hazard model (`coxph()`) which gives me HR of about 2.9 for presence of factor B (factors can be A, B, C, A as baseline in the model), with 95% CI 1.8-4.8 , p<0.001. When checking the proportionality assumption there is significant evidence that there is a violation of that assumption. Result of `plot(cox.zph)` for the model with factor A is shown below. ![](http://img31.imageshack.us/img31/7213/coxzph.jpg) My question is how should I understand the smoothing line of the graph, and what is its relation (and the relation of the values on the y-axis) to the beta estimate the `coxph()` function gives me (2.9 for the above example)? If there was no violation and the line of the `cox.zph` plot was straight, would the y-value of the line be (in this example) log(2.9)=0.46? If there is no violation of the proportionality assumption, does the "intercept" of the line equal the log of the HR that `coxph()` outputs?


title: Working with correlation coefficients - Fisher's transformation

text: Further to my earlier question on Working with correlation coefficients, I am still confused. I used the online calculator mentioned in the above post and if I compare .8978 (sample size: 105) and .7865 (sample size: 120), I come up with z=2.94 and p (one tail) = .0016 and p (two tail) = .0016. How do I interpret this result? How do I use the above calculator to compare three or more coefficients? Many thanks for assisting me.


title: Multiple testing correction of already corrected values

text: I have a list of say 50,000 e-values (multiple testing corrected P values) for a statistical test I am doing that essentially matches patterns together. Essentially the e-value is the expectation that these two patterns would match purely by chance. Most people who do this will do another round of multiple testing correction on their list of e-values which are already corrected for the number of tests that generate that e-value. For example each E-value may be the result of 100,000 tests matching different patterns together. The theory is when you have 50,000 of these e-values you will still have a good number of type I error's associated with your list. Can someone give me a some theory on why it's a good idea to apply multiple testing correction to values that were already corrected for multiple testing? Is there a good paper I can read that will help me understand? Conceptually I can understand why you want to do this, because that list of 50,000 e-vales will have some false positives, but if there is a paper someone can point me to that would be great.


title: How to describe the failure of this linear modelling?

text: I have a time series $X_t$, which is shown in the first plot. In the second plot, I am doing a linear regression on $X_t\sim X_{t-1}$. The regression line is very close to $y=x$. But this is tricky since if if we look at the bottom left or the top right part of the data, they are almost random. From the diagnosis of residuals, it is not a good regression either. But the model passes all the $t$ tests and $F$ tests. How can I say it's not a good model then? Is there a statistic to describe (not visually) the failure of this modelling? Here are the `R` codes I used to generate the plots:               # Generating X_t     x=c(arima.sim(list(order = c(1,0,0),ar=0.1),n=64,sd=1),3,5,7,11,14,17,rep(20,64)+arima.sim(list(order = c(1,0,0),ar=0.1),n=64,sd=1))     # Regression X_t~X_{t-1}     reg=lm(x[2:length(x)]~x[1:(length(x)-1)])     # Plotting     par(mfrow=c(3,2))     plot(x,xlab='',ylab=expression(X[t]),ty='l')     plot(x[1:(length(x)-1)],x[2:length(x)],xlab=expression(x[t-1]),ylab=expression(x[t]) ,main=paste('coeff= ',round(reg$coefficients[2],2)))     # Plotting the regression line     abline(reg,col=2)     # Plotting the residual diagnose     plot(reg)      ![enter image description here](http://i.stack.imgur.com/2xJ65.png)


title: What does ItemAverage Recommender do in Mahout

text: The description sounds to me as if it makes a MostPopular recommendation. But the MostPopular recommendation I did myself got much better results. So what does this recommender really return? It is a Boolean data set. And I just count the occurrences of an item and sort them with the counts.


title: Correlation Method

text: What statistical equation should I use to correlate two sets of segmentation data per the attached example. This would be similar to a dating site, where I have two people looking at one another’s profiles. Factoring in the following:   1. In some cases a variable might be non-existent as in the case of “Interests” as that field might not been completed by the user.    2. Also, many of the attributes are weighted differently by the opposite party i.e. Person A weights the importance of Person B’s and vice versa. * * *               ROW ATTRIBUTES     Age     Job Title     Company Size     Industry     Business Function     Seniority     Years of Experience     Interests     Revenue     Company Type     Skills          COLUMN - PERSON A – ATTRIBUTE VALUES     23     125     10     15     3     5     10     0     100     3     0          COLUMN - PERSON B – ATTRIBUTE VALUES     35     34     10     7     4     3     2     4     100     4     3          COLUMN - PERSON A WEIGHTED IMPORTANCE OF PERSONS B’S ATTRIBUTES     1        4        3        5        4        2        0        5        3        3        5          COLUMN - PERSON B WEIGHTED IMPORTANCE OF PERSONS BAS ATTRIBUTES     1     2     5     5     3     4     2     2     3     3     3


title: Neural network for multivariate time series

text: I am currently analysing multivariate time series and have worked on VAR models in R. I need to know if there is any way to analyse it using neural networks in R. PS: I am aware of the `nnet` package and `tsDyn` package, but I guess they are applicable for univariate time series only.


title: How does dimensionality of network topology affect self organising maps?

text: I'm pretty new to self organising maps (SOMs). From reading around a bit, it seems like most implemetations use a 2D grid for the SOM network structure (or less often a 1D grid). This is nice, cause it's easy to visualise, but is there any particular reason why barely anyone uses 3D grids, etc? In particular, how does a higher dimensional network topology affect computational difficulty, results, and interpretation of results?


title: visualization of interaction effects: using predicted values?

text: I am currently working on a study in which I am using beta regression (in Stata). My depend variable is bounded between 0 and 1. To interpret the effects of my independent variables I am using the margins command. However, the problem is that interaction-terms cannot be interpreted using the margins command. My model includes a quadratic term, i.e.: y=a+bX+cX^2. I want to visualize the non-linear connection between the variables by using graphs. My question might be a bit simple, but I am not sure which values I should choose to visualize the effect. Should I use the original depend variable or the predict values? Thanks in advance for the help.


title: Why is that the mean forecasting error is minimized when Conditional Expectation is used

text: As the title suggests, why is that the mean forecasting error is minimised when conditional expectation is used for forecasting in ARMA model? Need some guidance on this..


title: What does it mean to 'align data frequencies'?

text: If I have two variables that have different sampling frequencies, one of the first steps is to align the frequencies; could someone explain the intuition of this alignment without being too technical. I am reading some theory on doing MIDAS regression and can't grasp this concept. What does it mean to align data frequencies? for my case; I have a lower frequency variable Yt which is sampled monthly, and one higher frequency variable that is sampled daily. ![dayum](http://i.stack.imgur.com/yreaQ.png) For my purposes I need only one explanatory variable, so my mapping should only require the first row of the matrix. A simple example would help me a lot! Thank you


title: Constant RMSE margin between training and teseting set

text: I have a large number of independent datasets of varying size but same feature meaning. Features and outcome are both binary. I am trying to fit logistic regression to the data. I estimate regularization parameters with 10-fold random data partitioning (with classes proportion same as in training/testing) and select the best model based on accuracy. When I look at the results, I see that the margin between train and test accuracy decreases with the increase of the training set size, while the SD of accuracy in different folds decreases - as expected. ![](http://i.stack.imgur.com/C4UZy.png) ![](http://i.stack.imgur.com/mViET.png) However, when I look at the same plots for RMSE, the margin keeps constant regardless of the dataset size. Why could that be happening? ![](http://i.stack.imgur.com/6ahdr.png) ![](http://i.stack.imgur.com/T9nGL.png) In the histograms below we can also see that accuracy difference is centered around 0, while rmse difference is not. ![](http://i.stack.imgur.com/dROtk.png) ![](http://i.stack.imgur.com/AY7Zs.png) Thanks in advance!


title: Hat matrix,$H=X(X^{\prime}X )^{-1}X^{\prime}$

text: What is the importance of hat matrix, $H=X(X^{\prime}X )^{-1}X^{\prime}$ in regression analysis? Is it only for easier calculation ?


title: Finding improvement of a test population

text: I am gathering data across two variants of my website (akin to a split test). There is an action that I log each time a user invokes it along with which variant they are included in. I also log all of the users included in my test (even if they don't perform any action). Given a set of data for each variant now how can I tell a) What the average percentage increase/decrease is b) Whether the average percentage increase/decrease is statistically significant? Note: This is not a binomial distribution because it isn't just about people successfully performing an action. Instead I care how many times the users in each variation perform the action. Could I use something like eye-balling a box plot (and its IQR) to see if the two samples overlap any if at all and then just report on the percentage difference between the means?


title: What is the difference between the forward-backward and Viterbi algorithms?

text: I want to know - what are the differences between the forward-backward algorithm and the Viterbi algorithm for inference in hidden Markov models (HMM)?