title: Basic idea of zero inflated two part models(hurdel) and application to big data (machine learning)

text: I'm currently working on the data which has 90% 0s in response variable. Based on my research, it seems zero inflated models could be a solution to this. However, while I was reading related documents, it seems the basic idea of two part model( hurdel) is $ E[Y|X] = P(Y>0|X) E[Y|X,Y>0]$ The question is will it be okay if I use any model to get both `P(Y>0|X)` and `E[Y|X,Y>0]` ,and multiply the result of those two? More specifically, in R code would be like               #for first model     FirstModel=randomforest(Y,X) #random forest      First_pred=predict(FristModel,newdata,type='prob')          #for second model     SecondModel=glmnet(Y~X) #glm with L1 L2 penalty     Second_pred=predict(SecondModel,newdata)               #for E[Y|X] = P(Y>0|X) * E[Y|X,Y>0]     First_pred * Second_pred      Will this make sense in view of machine learning modeling strategy? In other wards, will this improve accuracy of prediction without considering interpretation?


title: Assistance with hazard ratios

text: I want to compare car accident deaths in two groups of people aged 18-24 (lets use this as the reference) and 24-30 after a certain law was changed using a Cox hazard model using SPSS version 22. I want unadjusted results, but then want to enter covariates such as whether or not they are on medications, gender, and speed in the adjusted model. How do I input these into SPSS? I know to go to analyze --> survival analysis --> cox regression but keep getting stuck after. The variables I have are age, whether or not there is a death, time to death after the law change in addition to the ones above. Thank you!


title: Number of ngrams for vocabulary

text: I read the great book Foundations of Statistical Natural Language Processing. by Christopher D. Manning and Hinrich Schütze. My question is regarding counting the number of ngrams for a given vocabulary. For example, in the book, the given vocabulary has a size of 20000 words and the number of ngrams are as following bigram model: $20, 000 × 19, 999 = 400 $ million trigram model: $20, 000^2 × 19, 999 = 8 $ trillion four-gram model: $20, 000^3 × 19, 999 = 1.6 × 10^{17}$ But why it actually so, why the last element is one less than the previous, I thought is should be simpler for bigrams just square of vocabulary.


title: Maximum number of independent variables in Logistic Regression

text: Is there a measure in logistic regression that maybe penalizes you for having too many independent variables like in multiple regression with the adjusted R squared? That is, does having too many independent variables in a logistic regression hurt the model? What about dummy variables? Can you have too many of those to the point of unpredictability?


title: Mixed Effects Model with Nesting

text: I have data collected from an experiment organized as follows: Two sites, each with 30 trees. 15 are treated, 15 are control at each site. From each tree, we sample three pieces of the stem, and three pieces of the roots, so 6 level 1 samples per tree which is represented by one of two factor levels (root, stem). Then, from those stem / root samples, we take two samples by dissecting different tissues within the sample, which is represented by one of two factor levels for tissue type (tissue type A, tissue type B). These samples are measured as a continuous variable. Total number of observations is 720; 2 sites * 30 trees * (three stem samples + three root samples) * (one tissue A sample + one tissue B sample). Data looks like this...                       ï..Site Tree Treatment Organ Sample Tissue Total_Length         1        L  LT1         T     R      1 Phloem           30         2        L  LT1         T     R      1  Xylem           28         3        L  LT1         T     R      2 Phloem           46         4        L  LT1         T     R      2  Xylem           38         5        L  LT1         T     R      3 Phloem          103         6        L  LT1         T     R      3  Xylem           53         7        L  LT1         T     S      1 Phloem           29         8        L  LT1         T     S      1  Xylem           21         9        L  LT1         T     S      2 Phloem           56         10       L  LT1         T     S      2  Xylem           49         11       L  LT1         T     S      3 Phloem           41         12       L  LT1         T     S      3  Xylem           30      I am attempting to fit a mixed effects model using R and lme4, but am new to mixed models. I'd like to model the response as the Treatment + Level 1 Factor (stem, root) + Level 2 Factor (tissue A, tissue B), with random effects for the specific samples nested within the two levels. In R, I am doing this using lmer, as follows               fit <- lmer(Response ~ Treatment + Organ + Tissue + (1|Tree/Organ/Sample))       From my understanding (...which is not certain, and why I am posting!) the term:               (1|Tree/Organ/Sample)      Specifies that 'Sample' is nested within the organ samples, which is nested within the tree. Is this sort of nesting relevant / valid? Sorry if this question is not clear, if so, please specify where I can elaborate. Many thanks to any help you can offer. ~Erik


title: How to design questions to implement parametric test?

text: I am planning to do research on micro-finance institutes in India. Due to my limited knowledge on questionnaires I am facing challenges as to what kind of questions should be asked (Dichotomous, Likert scale) so that T- test, Anova , chi square or such parametric test can be implemented. Any help or sample questions shall be appreciated. Here are the research questions. RQ1: What are the categories of people that take loans from micro-finance institutions based on their income? RQ2: Has microfinance evolved in India? RQ3: What is the impact of microfinance on Indian economy?


title: How to get predictions in terms of survival time from a Cox PH model?

text: I want to develop a prediction model (Cox PH) for all-cause mortality in a dataset of participants of whom (almost) all have died at the end of follow-up (e.g. 1-year). Instead of predicting the absolute risk of dying at a certain timepoint, I would like to predict the survival time (in months) for each individual. Is it possible to obtain such predictions in R (from e.g. a coxph-object) and, if yes, how can I do that? Many thanks in advance!


title: Factored Joint Distribution of Tree Augmented Naive Bayes Algorithm

text: I asked this question at Mathoverflow but they recommend me to ask here. I need to find factored joint distribution of Tree Augmented Naive Bayes algorithm. I read the paper but I couldn't figure out the answer. Any help or pointers appreciated.


title: What is the precise definition of a "Heywood Case"?

text: I had been using the term "Heywood Case" somewhat informally to refer to situations where an online, 'finite response' iteratively updated estimate of the variance became negative due to numerical precision issues. (I am using a variant of Welford's method to add data and to remove older data.) I was under the impression that it applied to any situation where a variance estimate became negative, either due to numerical error or modeling error, but a colleague was confused by my usage of the term. A google search doesn't turn up much, other than that it is used in Factor Analysis, and seems to refer to the consequences of a negative variance estimate. What is the precise definition? And who was the original Heywood?


title: How to take many samples of 10 from a large list, without replacement overall

text: I've got a large set of data (20,000 data points), from which I want to take repeated samples of 10 data points. However, once I've picked those 10 data points, I want them to not be picked again. I've tried using the `sample` function, but it doesn't seem to have an option to sample without replacement over multiple calls of the function. Is there a simple way to do this?


title: How to make forecasts for a time series?

text: I am not that familiar with the analysis of time series data. However, I have what I think is a simple prediction task to address. I have about five years of data from a common generating process. Each year represents a monotonically increasing function with a non-linear component. I have counts for each week over a 40 week cycle for each year. The process starts, the function begins at zero, increases rather quickly over the first half of the function, slowing over the second half before leveling during the last five weeks. The process is consistent across years with small differences in rate of change and volume across the segments from year to year. y1={0, Nt1, Nt2, ... Nt39, Nt40) ... ... y5={0, Nt1, Nt2, ... Nt39, Nt40) Where Ntx equal the count at time x. The goal is to take the N at tx (or better t0 to tx, or the slope to that point) and predict the N at t40. For example, if Nt10 is 5000 what is the expected value of Nt40 for that year. So, the question is, how would you model such data? It's easy enough to summarize and visualize. But I'd like a model to facilitate predictions and incorporate a measure of error.


title: Graphical data overview (summary) function in R

text: I'm sure I've come across a function like this in an R package before, but after extensive Googling I can't seem to find it anywhere. The function I'm thinking of produced a graphical summary for a variable given to it, producing output with some graphs (a histogram and perhaps a box and whisker plot) and some text giving details like mean, SD, etc. I'm pretty sure this function wasn't included in base R, but I can't seem to find the package I used. Does anyone know of a function like this, and if so, what package it is in?


title: Variance of sample correlation coefficient

text: On wiki page about Fisher transformation I read that variance of sample correlation coefficient becomes smaller as population correlation coefficient (in absolute value) gets closer to 1. Could anybody, please explain why this happens? Thanks a lot for any relevant formulas or intuitive explanations!


title: What is the reasonable methodology to study machine learning techniques efficiently?

text: Nowaday, there are lots of online courses related to machine learning, and there are lots of programming languages which can be used to do data mining, such as R, Python or Matlab. For example there are lots of data science related courses on coursera, some use R and some use Matlab. Lots of the time I just download the videos to my disk, and never open it again, other times I spend all my time to finish the courses and just forget the real problem I want to study. When I am using machine learning tools, I find it's better to know some underlying mathematics, then I would pick up some book such as convex optimiztion, but this time I find it's too slow to master such materials. Seems there is no plan at all when I study machine learning techniques. So anybody here who can share something about your methodology when you try to do machine learning or data mining, do you have systemic and effective plans?


title: Probability of a sale in a given month based on past activity including order frequency

text: The sales history for customer A looks like this:               Jan : $500 (4 sales)          Feb : $200 (2 sales)          Mar : $0 (0 sales)      etc.. I want to know:   1. The probability of customer A having a sale occurring in a month.   2. Can we do anything interesting as far as taking into account the frequency of purchases in a month?   3. Can we do anything interesting as far as taking into account the value of the purchases in a month? The goal is to find, out of a large set of customers, those that have a high probability of purchasing or a high sale potential based on frequency/value.


title: R GMM - Error in solve.default(x$v, gb) : system is computationally singular: reciprocal condition number

text: I'm having the following problem estimating something in GMM in R. I have created a "Hello World" below. In principle, I would not need GMM to estimate the parameters, but I want to use it to obtain standard errors as in my application there is likely time series dependence (there is not in the example below). Thank you in advance. The code is below:               #data     set.seed(12345)     r1 <- 0.5     x1 <- rnorm(n=5000,mean=0,sd=0.05)     x2 <- rnorm(n=5000,mean=0,sd=0.05)     x3 <- 0.5 * x1 + sqrt(1-0.5**2) * x2          #gmm estimating function     g <- function(theta,x) {         m1 <- theta[1] - x[,1] # mean 1         m2 <- theta[2] - x[,2] # mean 2         m3 <- theta[3] - (x[,1]-theta[1]) * (x[,2]-theta[2]) # covariance         m4 <- theta[4] - x[,1] * x[,2]         m5 <- theta[5] - (x[,1]-theta[1]) * (x[,2]-theta[2]) / (theta[1]*theta[2]+theta[3])         f <- cbind(m1,m2,m3,m4,m5)         return(f)     }          #call the gmm routine     library(gmm)     iniVal <- c(mean(x1),mean(x3),cov(x1,x3),mean(x1*x3),cov(x1,x3)/mean(x1*x3))     dat <- as.matrix(cbind(x1,x3))     gmmres <- gmm(g,dat,iniVal,wmatrix="ident")     summary(gmmres) # singularity problem


title: Instrumental Variables for Logistic Regression

text: I am trying to regress a Ratio variable,Y, on an independent variable,X; Variable X is endogenous, so I have to use an Instrumental Variable, Z; both X and Z are continues variables. How can I run this model in SAS. I know proc SYSLIN easily does it if the dependent variable is continues, but what if the dependent variable is a ratio or dichotomous? PROC QLIM only allows the left hand variable to be a binary outcome and endogenous! again, when the right hand side variable is endogenous, QLIM is not useful. If there is no built in procedure, is there any alternative way to run this model in SAS? A general recommendation is also appreciated , I do not need the answers to be specifically for SAS. Thanks for your help in advance.


title: Do kernel methods "scale" with the amount of data?

text: I've been reading about kernel methods, where you map original $N$ data points to a feature spaces, compute the kernel or gram matrix and plug that matrix into a standard, linear algorithm. This all sounds good when the feature space is infinite dimensional or otherwise very high-dimensional (much much larger than $N$), BUT the kernel matrix itself is also pretty large at $N \times N$, meaning if you double the number of points you quadruple the amount of memory required. Does this mean kernel methods do not scale well to larger data sets? Or is it not necessary to compute the entire kernel matrix and hold the whole thing in memory for most algorithms?


title: Random variable as a function

text: I'm reading Schaum's outline of probability, random var. and random processes. In the second chapter they make it clear that a random variable $X$ is not a variable in the traditional sense but rather a function $X(\zeta)$ where $\zeta$ is a sample point in sample space $S$. So $S$ is the domain of r.v. $X$ and collection of all numbers $X(\zeta)$ is the range of r.v. $X$ I'm trying to relate this to distribution of features used on audio. Take MFCCs for example. Let's assume we have two different audio signatures each having a range of MFCC values, so we could have a MFCC distribution for each of the two classes. If we define the numerical value of r.v. as the values taken by MFCCs, it's clean that range of r.v. $X$ is from $-\infty$ to $+\infty$ (in N dimensions, N = dimensions of MFCC), but what is the sample space here? And what is the $\zeta$? i.e. If $X(\zeta)$ = MFCC_n, what is this $\zeta$ we are referring to?


title: Interpretation of two residual QQ plots from poisson regression

text: I have what I believe are an interesting QQ plot of the residuals of some sports data fitted using Poisson regression. My model is actually a live model to predict the number of remaining events in the game (basketball, soccer etc) from a set of explanatory variables. It looks also very much like the data is Poisson distributed, with equal mean and variance. Residuals are as we know $r_i=Y_i-\mu_i$, where $Y_i$ are the observations, and $\mu_i$ the fitted mean of $Y_i$ ($Y_n$ may have different explanatory variables than, for example, $Y_{n-1}$). Here is the first residual QQ plot of fitting my model to soccer match data (whole 90 minutes): ![Residual QQ plot number 1](http://i.stack.imgur.com/nDi2u.png) Interpreting QQ plots if I want to check for normality of some data is not too hard. It is easy to see that the residuals have too thick tails. But assuming the data is without a doubt Poisson distributed, something must therefore be wrong with $\mu_i$. Looking at the lower tail, I cannot reach any other conclusion than that the values of $\mu_i$ generally are too big - such that $r_i$ becomes too small. But also, looking at the upper tail, I cannot reach any other conclusion than that $\mu_i$ is too small - such that $r_i$ becomes too big. Is this a contradiction, or merely just a reflection that $\mu_i$ is not fitted good enough for extreme cases?


title: Books to refer enhancing my knowledge in statistics

text: I am a new entry in statistical field with lots of interest in this area. I would like to know any statistical book which can enhance my statistical knowledge meaning; wish to know the name of the book which can help me in model prediction....for eg. detailed explanation and steps of regression (incl multiple, logistic, poisson...) , Anova, etc.... The book should contain each and every concept of stats so that I can practically handle projects in my company.


title: time series is obviously periodic, but seasonal decomposition is not working in R

text: my time series is obviously periodic, but the seasonal decomposition using stl() is not working in R:               a <- c(6.7, 20.3, 23.5,  7.9,  3.3,  2.0,  2.5,  2.9,  2.3,  5.0, 15.0, 20.1, 27.0, 28.2, 18.3,  7.8,  1.6,  0.8,  1.3,  1.2,  0.6,     1.6,  4.9, 24.2, 28.8, 23.6, 18.6,  5.3,  1.8,  0.4,  0.5,  0.2,  0.1,  0.3,  3.5, 17.6, 26.1, 22.7, 18.2,  7.2,  2.1,  1.0,     1.1,  1.4,  0.7,  3.2,  7.6, 22.3, 30.4, 26.5, 23.2,  8.6,  4.3,  1.5,  1.7,  2.4,  0.9,  3.0,  8.0, 21.9, 25.9, 25.6, 25.2,     13.1,  5.4,  2.2,  2.8,  2.7,  1.9,  5.1, 12.4, 28.4, 27.8, 22.2, 21.1, 13.6,  6.1,  4.0,  3.7,  4.5,  3.6,  6.3, 11.1, 25.0,     31.1, 24.5, 22.1, 10.9,  5.2,  2.6,  4.4,  4.1,  3.1,  7.0, 17.3, 25.5, 30.1)     a.ts <- as.ts(a, frequency=12, begin=c(1980,1))     stl(a.ts,s.window = "periodic")      Output:               > stl(a.ts,s.window = "periodic")     Error in stl(a.ts, s.window = "periodic") :        series is not periodic or has less than two periods      ![Plot](http://i.stack.imgur.com/GmejJ.png) My main goal is to remove the seasonality of the time series. Thanks for any help!


title: Dice probability - amount of "hits" in a pool

text: What I want is to calculate the probability of a minimum number of rolls with a minimum value for at given number of dice, e.g. "if I roll 20 dice, what is the probability of getting at least 3 dice with a value of 5 or more?" I have a brute force setup, that creates the Cartesian product and just enumerates through that, but it completely explodes at 9 dice. So what I'm really looking for is some kind of function that does not require me to expand all possibilities, so that I can calculate for 20-30 dice.


title: What is it meant by the "rejection region" and "power" of a likelihood ratio test?

text: Suppose that $X_1,...,X_n$ are i.i.d. data from a $N(\mu, 100)$ distribution. I am trying to find the rejection region for the likelihood ratio test for level $\alpha= 0.10$ of the test: $H_0: \mu = 0$ versus $H_1: \mu= 1.5$ Finally, if $n = 25$, I would like to find the power of the test. What confuses me about this question is that it is sort of non-standard. My approach is to first write out the likelihood ratio: $$\newcommand{\EXP}[1]{\exp\left[#1\right]} \frac{L(0)}{L(1.5)} = \EXP{\frac{\sum(X_i-.1.5)^2-\sum{X_i^2}}{2(100^2)}}$$. Then, would my rejection region be: $$\EXP{\frac{\sum(X_i-.1.5)^2-\sum{X_i^2}}{2(100^2)}} \leq k\text{ ?}$$ And finally, I am very confused about what the power of the test is, because it seems that I have already computed it in the process of writing out: $$P\left(\EXP{\frac{\sum(X_i-.1.5)^2-\sum{X_i^2}}{2(100^2)}} \leq k\right)\leq \alpha = 0.10$$ and then, I would rearrange the above until I get $\sum{X_i}$ by itself, then apply the central limit theorem, then work backwards to find $k$. However, I still dont know what they mean by the power of a test. I feel like it is staring right at me, but I am not sure what it is. Any help would be greatly appreciated!


title: SUTSE DLM on daily mean water & air temperature TS

text: I have two time series: (1) daily mean water temperature from 1988 to 2014 and (2) daily mean air temperature from 1968 to 2010. The water temperature time series has missing data, occurring on individual days and for periods of consecutive days (e.g., data for all of 1993 is missing). The air temperature data has data missing from 1990-02-01 - 1991-09-30. The data are plotted below: ![Water temperatures. Red bands indicate missing data. Smooth is annual GAM.](http://i.stack.imgur.com/7c2bF.jpg) ![Airtemperatures. Red bands indicate missing data. Smooth is annual GAM.](http://i.stack.imgur.com/12Wur.jpg) I want to use the air temperature series to inform missing values in the water temperature series. From my understanding, this could be achieved using a state-space or Dynamic Linear model (DLM). I’m imagining a Seemingly Unrelated Time Series Equations (SUTSE) DLM in which both temperature series follow a stochastic linear trend and annual seasonal component. So far (I’m only learning DLM), I have managed to fit a simple 2nd order stochastic linear trend SUTSE DLM to the water and air time series (both reduced to the same period 1986-08-19 to 2010-04-15) as follows (I’m afraid a data-sharing agreement forbids me to provide the data):               # the 'base' univariate model     uni <- dlmModPoly()          # to get the matrices of the correct dim     temps_mod <- uni %+% uni           ## now redefine matrices to keep levels together and slopes together     FF(temps_mod) <- FF(uni) %x% diag(2)     GG(temps_mod) <- GG(uni) %x% diag(2)     W(temps_mod)[] <- 0 # 'clean' the system variance          ## define a build function for MLE     buildSUTSE <- function(psi) {       U <- matrix(0, nrow = 2, ncol = 2)       U[upper.tri(U)] <- psi[1:2]       diag(U) <- exp(0.5 * psi[3:4])       W(temps_mod)[3:4, 3:4] <- crossprod(U)       diag(V(temps_mod)) <- exp(0.5 * psi[5:6])       temps_mod     }          ## estimation     temps_est <- dlmMLE(temps, rep(-2, 9), buildSUTSE,                      control = list(maxit = 500)) ## few iters for test          ## check convergence     print(temps_est$conv)          ## set up fitted model     temps_mod <- buildSUTSE(temps_est$par)          ## look at the estimated variance / covariance matrices     print(W(temps_mod)[2:3, 2:3])     print(cov2cor(W(temps_mod)[2:3, 2:3])) # slopes     print(sqrt(diag(V(temps_mod)))) # observation standard deviations          ## smooths     tempsSmooth <- dlmSmooth(temps, temps_mod)      which results in the following fits: ![SUTSE DLM linear stochastic trend fits](http://i.stack.imgur.com/e4O2j.jpg) While it is encouraging to see that the fitting has worked, the 1993 water temperature estimates look poor / too low relative to the surrounding values. Here are my questions:   * Am I taking the right approach, or is there a better / simpler way to estimate missing water temperature data from air temperature?   * How do I extend the above code to include a seasonal trend, with (presumably) a frequency of 365 (even though leap years are present)? Should I be considering trigonometric functions, e.g., dlmModTrig?   * Has anyone seen this attempted before, and if so, where can I read about it (preferably with an R example)? Thanks in advance for any help.


title: Is this hypothesis test somehow "optimal"?

text: I have a system whose output I am analyzing. I've used a very simple hypothesis test with great results and now am curious _why_ did I get such great results -- did I perhaps stumble on some kind of an "optimal" test without much formal statistics training? My system can be in two states. When it's in state 0 the output can is independently drawn from standard normal distribution $\mathcal{N}(0,1)$ (well, technically, very, _very_ close approximation to it-- but it's normal enough for my purposes). When it's in state 1, the output is still independently normally-distributed with variance 1 but the means are perturbed. Thus, the output $i$ in state 1 has a mean-shifted normal distribution $\mathcal{N}(\mu_i,1)$. For large number $n$ of observations of the system in state 1, the average $\frac{1}{n}\sum_{i=1}^n\mu_i=0$. I also know that quantity $\frac{1}{n}\sum_{i=1}^n\mu_i^2=M<\infty$ since this is a real-world system (though I do not know the exact value of $M$, just the fact that it's finite). Besides those two facts, there is no structure to means in state 1. I devised the following hypothesis test to determine whether the system is in state 0 or 1, where the null hypothesis $H_0$ corresponds to the system being in state 0, and alternate $H_1$ to state 1. First, I collect a sequence of $n$ observations $\{x_i\}_{i=1}^n$. Then I compute the following test statistic: $$S_n=\left(\frac{1}{n}\sum x_i^2\right)-1$$ I then pick a threshold $t>0$ and accept null hypothesis if $S_n<t$, rejecting it if $S_n\geq t$. This test just made sense to me and is surprisingly easy to analyze (at least surprising to me without much formal background in mathematical statistics). The mean of $S_n$ when the system is in state 0 is 0, and the variance is $2/n$. In state 1 the mean of $S_n$ is $M$ and the variance is $(4M+2)/n$. Thus, I can upper-bound both the false-positive error probability $\alpha$ and the probability that I will accept the null hypothesis in error $\beta$ for a given $n$ and $M$ using the Chebyshev's Inequality: $$\begin{array}{rcl}\alpha&\leq&\frac{2}{nt^2}\\ \beta&\leq&\frac{4M+2}{n(M-t)^2} \end{array}$$ Using the two inequalities I can pick an appropriate threshold $t$ and the number of observations $n$ I need to classify the state of the system for some particular value of $M$ given some tolerance for errors $\alpha$ and $\beta$. The math is simple and completely makes sense: for some fixed $\alpha$, as I increase $n$, my threshold $t$ decreases, and so does $\beta$. Is the reason for this beauty and simplicity that the test I devised is in some way "optimal"? If so, can the community suggest a way prove this optimality? I did find "Uniformly Most Powerful tests" in a stats book, but I am not sure if I have that (the book was quite confusing)...


title: LASSO or other regularized regression with censored (missing) data

text: Here is my problem. I am looking at various time series curves. Let's call them total spend aggregated over all customers on various products versus time. At any given time, I want to predict the spend on any product for which I don't have data, based on the spends from the products for which I do have data. The number of products is much larger than the number of time data- points so this requires a highly regularized solution. I am starting with using LASSO to do this. The trouble is that the products all have a finite lifetime. Outside of this lifetime they show zero spend or are just missing. Sometimes within their lifetime, they also have zero spend but in that case, these zeros really mean zero. Outside of their lifetime, the zeros just indicate that the product was not available and so the customer spend likely went to another product that was available with similar functionality. So the question is: what kind of algorithm can fit such a regularized linear model while treating these periods outside of the lifetime as missing data. Assume that, I know where to draw the boundaries on the lifetimes. I'd like implementations in either R or Python. Perhaps the flare package in R is the place to start. In there it is mentioned that the Dantzig selector and CLIME are both tolerant of missing values in the design matrix. Can anyone make a suggestion? Note that I am trying to make a predictor for all variables so all the better if I can do this all at once by estimating the precision matrix of a Gaussian graphical model directly.


title: Individual and overall RMSE for multivariate data

text: I have a dataset which contains missing values, and I'm using imputation packages (`R`s `mi` and `mice`) to fill the missing values. I'd like to measure their performance on my data set, which may look like this (my actual data has more rows, but this example should serve just fine):               "var1"  "var2"  "var3"     0.183689466952776   0.388415623304919   -1.3390868493301     NA  -2.0495669969489    NA     NA  NA  1.1107054143715     1.29820089212697    0.736777347408364   -1.19623852541909     -1.17191167149872   -0.744790411450254  -1.96820040415179     -0.686058069998857  NA  NA     0.96219165971458    -1.26927815931595   1.13102353621198     NA  -0.181582994079309  -1.88246768436578     0.133837989978951   -0.298476696697043  0.887731971394049     0.42775517098228    -1.91391435336026   NA     NA  NA  0.0027473853587295     0.605986709105715   0.297153545105678   -1.03855048360928     NA  -1.18987831904712   1.0500895435177     NA  -0.219325915775778  1.54228872681253     NA  NA  -0.976306655339306     NA  NA  0.440861027292491     -1.92738847897133   -0.779770748497074  0.403377851347805     NA  -0.8839601961621    0.0382354592857369     -1.79066885776893   0.723084216521015   0.287610507512217     NA  -2.70392018097682   0.744853382274342      Note the different ratios of missing data (50%, 25%, 10%) and that (in this case, by construction) the pattern of missing data is random. In order to measure the imputation error, I replace some of the non-NA values in each variable by `NA`. To keep the structure of the data the same, I opted to replace 10% of the non-`NA` entries in each variable by `NA`, i.e. 1, 2 and 2 in `var1`, `var2` and `var3`, respectively. Other possible ways to create a test set of values would be to pick 10% of all non-missing values (regardless of value), or to pick the same absolute number in every variabe, i.e., say, 1 in every colum. I measure the quality of the imputation by calculating the RMSE of every column separately. Is there a way to calculate an overall rmse which takes into account that that my test set was created in a stratified way? Is it OK to calculate it as $\sqrt{\frac{RMSE_1^2 +RMSE_2^2 + RMSE_3^2 }{3}}$, or does it need to be weighted to reflect the different sample sizes used to calculate the individual RMSEs? Calculating the root of the mean of the MSEs is something I found in What is the RMSE of k-Fold Cross Validation?, but the answer is quite short why this is the right way to calculate the overall RMSE, and I'm also not sure if this formula from CV is applicable to my case.


title: MRF MAP inference for non-submodular pairwise terms

text: I have a multilabel MRF MAP inference problem (a labeling problem). The graph has relatively few nodes, about a thousand or so. The pairwise term is (very) not submodular (it does not satisfy the triangle inequality). GraphCut with truncation does not work as about 40% of the entries are not submodular. What is the state of the art for this type of problem? Loopy Belief Propagation (BP), Tree Reweighted Message Passing (TRW-S), Quadratic Pseudo-Boolean Optimization (QPBO), Iterated Conditional Modes (ICM)?


title: How to perform Gaussian process regression when  function being approximated changes over time?

text: What are good strategies for performing Gaussian process regression when the function I am trying to approximate changes over time? The naive approach that springs to my mind is to only use the N most recent data points to perform the regression. What are better strategies?


title: Large performance variance when using back-propagation in neural networks for feature extraction

text: I'm using a multi-layered neural network for feature extraction (similar to deep belief network). I test the performance of my model with cross-validation. When I'm using back-propagation to train my network, I get a very large variance in the performance. Seems like it is really crucial which samples (and maybe in what order?) are used for back-propagation. Is it this a known issue with NN? I tried using methods from here (I have unequal classes) http://www.springerlink.com/content/2wmmd867cyvbbk3h/ which helped performance but not the variance. UPDATE: I'm using 256 features and I have about 2000 training samples. I'm training my network with layer-wise pretraining and then back-propagation. I'm using weight decay (around 0.02) and a learning rate of about 0.01 . Thanks.


title: PCA, LDA, CCA, and PLS

text: How are PCA, LDA, CCA, and PLS related? They all seem "spectral" and linear algebraic and very well understood (say 50+ years of theory built around them). They are used for very different things PCA: dimensionality reduction, LDA for classification, PLS for regression but still they feel very closely related.


title: Average case analysis of learning algorithms

text: Typically analysis of learning algorithms is in the worst-case setting, for example regret bounds in online learning, or generalisation error bounds in classification. Whilst worst case performance is clearly an important question, it seems that average case analysis would be much more indicative of an algorithm's real-world performance. Is such analysis even possible? What assumptions would needed to make it work?


title: Continuous Multistate Ada-boost method?

text: I recently read this simple introduction to Adaboost as a review after learning about it a few years ago. This is in preparation for actually using it to solve a problem I am working on which is slightly different than anything addressed in this paper or that I have found online. Now in the problem I am trying to solve there is a continuous scale, rather than the binary problem addressed in the paper (and almost everywhere online) or even the distinct multistate (kind of like an enum) problems they mention at the end of the paper. Imagine experts scoring some performance on a continuous scale between 1-10 where scores like 8.4 mean twice as good as 4.2. The way that they did the multistate version in the paper was to break it into a series of binary tasks but this completely defeats the point since it throws out a lot of the information about what the scale means and just makes it a named class. I was wondering what methods people know of for this continuous problem and if an adaboost idea can still be used. thanks.


title: Bootstrap confidence bands for time series data

text: I have time series data for 30 people. The data is knee joint angle versus percent of gait cycle, with data at every 1% of the cycle. I can easily find the mean joint angle at every 1% of the gait cycle. I would like to use bootstrapping to find the 90% confidence bands over the entire gait cycle. I have read in the literature that confidence bands formed by bootstrapping will include close to the desired percentage of the population, whereas gaussian estimation of the confidence bands will include much less than the desired percentage of the population (90% in this case). I can calculate the point-by-point (i.e. every 1% of the cycle) confidence bands using bootstrap, but I do not think that is appropriate. Can anyone explain in plain language what are the steps for calculating the confidence bands, given that I have 30 subjects? I have read a bit on block bootstrapping, but could not figure out how that would work with my 30 samples. I prefer to write the code myself, but if there is any canned software that will calculate the confidence bands for a set of 30 samples such as I have, I would be willing to give it a try. Any advice you can provide is appreciated.


title: Low intraclass correlation although between variance is very high in random effects model

text: I don't understand why a random effects model can give you a very low `rho` (intraclass correlation coefficient), if the between variance (between subjects) in your dataset is very high on all independent variables compared to the within variance (over time). Maybe my interpretation of `rho` is incorrect and somebody can help me understand it correctly. As you can see from this table, between variane is much higher than within variance for almost all independent variables:               . xtsum IV_1-IV_10          Variable         |      Mean   Std. Dev.       Min        Max |    Observations     -----------------+--------------------------------------------+----------------     IV_1     overall |  9703.536   92324.06          0    3551455 |     N =  128238              between |               105455          0    3509721 |     n =    2497              within  |             2094.317  -226060.1   131299.5 | T-bar = 51.3568                      |                                            |     IV_2     overall |  725.3808   3339.206          0   106827.5 |     N =  128238              between |             3769.248          0     106817 |     n =    2497              within  |             71.34145  -2631.032   13295.87 | T-bar = 51.3568                      |                                            |     IV_3     overall |  960.2351   574.6654   .0028935   2656.622 |     N =  128238              between |             614.7887   .4450117   2635.694 |     n =    2497              within  |             69.42291   714.5346    1170.56 | T-bar = 51.3568                      |                                            |     IV_4     overall |  .0446572   .2065295          0          1 |     N =  128238              between |             .2076612          0          1 |     n =    2497              within  |             .0183641  -.9040607   .9840511 | T-bar = 51.3568                      |                                            |     IV_5     overall |  .5832515   .4930224          0          1 |     N =  128238              between |              .498173          0          1 |     n =    2497              within  |                    0   .5832515   .5832515 | T-bar = 51.3568                      |                                            |     IV_6     overall |  .8897909   .3131515          0          1 |     N =  128238              between |             .3459811          0          1 |     n =    2497              within  |                    0   .8897909   .8897909 | T-bar = 51.3568                      |                                            |     IV_7     overall |    .70823   .4545788          0          1 |     N =  128238              between |             .4657856          0          1 |     n =    2497              within  |                    0     .70823     .70823 | T-bar = 51.3568                      |                                            |     IV_8     overall |  .1469527   .3092312          0          1 |     N =  128238              between |             .2002412          0          1 |     n =    2497              within  |             .2398179  -.8360982   1.137862 | T-bar = 51.3568                      |                                            |     IV_9     overall |  5.341053   15.03976          1        453 |     N =  128238              between |             17.59767          1   308.2857 |     n =    2497              within  |             7.586072  -242.5589   329.5824 | T-bar = 51.3568                      |                                            |     IV_10    overall |  13.00411   35.60811          0       2230 |     N =  128238              between |             27.50115          0        279 |     n =    2497              within  |             30.47548  -265.9959   2098.004 | T-bar = 51.3568                 etc      For my in dependent variable, I get:               . xtsum DV          Variable         |      Mean   Std. Dev.       Min        Max |    Observations     -----------------+--------------------------------------------+----------------     DV       overall |  .3408194   .4739866          0          1 |     N =  128238              between |             .1320563   .0191388   .8740157 |     n =    2497              within  |             .4584589  -.5331963   1.321681 | T-bar = 51.3568      My model output however states a `rho` of 0.0407. Please correct me if that's wrong but this means that the share of between-variance on the total variance in this model is ~4% .               .xtlogit DV IV_1-IV_10, re              Random-effects logistic regression              Number of obs      =    128238     Group variable: nameID                          Number of groups   =      2497          Random effects u_i ~ Gaussian                   Obs per group: min =         2                                                                    avg =      51.4                                                                    max =       209                                                          Wald chi2(21)      =   9381.54     Log likelihood  =  -75336.99                    Prob > chi2        =    0.0000          --------------------------------------------------------------------------------------                  DV      |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]     ---------------------+----------------------------------------------------------------             IV_1         |  -9.35e-08   1.15e-07    -0.81   0.418    -3.20e-07    1.33e-07             IV_2         |   2.33e-07   2.86e-06     0.08   0.935    -5.37e-06    5.84e-06             IV_3         |  -.0000156   .0000191    -0.82   0.413     -.000053    .0000218             IV_4         |   .0249436   .0546815     0.46   0.648    -.0822301    .1321173             IV_5         |  -.0384232   .0237309    -1.62   0.105     -.084935    .0080885             IV_6         |    .089663    .036282     2.47   0.013     .0185516    .1607743             IV_7         |  -.0131013   .0255441    -0.51   0.608    -.0631668    .0369641             IV_8         |   -.231409   .0253735    -9.12   0.000    -.2811401   -.1816779             IV_9         |   .0073497   .0006319    11.63   0.000     .0061112    .0085882            IV_10         |  -.0000409   .0002392    -0.17   0.864    -.0005096    .0004279                    _cons |  -1.924186   .0481483   -39.96   0.000    -2.018555   -1.829817     ---------------------+----------------------------------------------------------------                 /lnsig2u |   -1.96814    .054337                     -2.074639   -1.861641     ---------------------+----------------------------------------------------------------                  sigma_u |   .3737867   .0101552                      .3544034      .39423                      rho |   .0407386   .0021234                      .0367744    .0451101     --------------------------------------------------------------------------------------     Likelihood-ratio test of rho=0: chibar2(01) =  1534.08 Prob >= chibar2 = 0.000


title: What is "estimated unbiased variance of the error term"?

text: Disclosure: _This is a homework question._ I have fit a multiple linear regression model in eviews, and I am asked to calculate "estimated unbiased variance of the error term, i.e., $\hat\sigma^2$". What confused me was "estimated unbiased" part. I have done a Google search for exact term "estimated unbiased variance of the error term". However, I didn't get an exact result. My approach was to calculate variance of residuals through `genr varresid = @var(resid)` (eviews command). However, I am not sure if this is the number I am trying to get. So, my question is, **what is the formula for the estimated unbiased variance of the error term?**


title: Clustering hetrogeneous data types: ordinal, interval

text: Say I have data that I'd like to cluster that has different dimensions that are of different data types. For example:   1. ordinal: You mood today: Very happy, happy, neutral, sad, very sad   2. Ratio: Age: ...   3. Interval: Weight: 100-120, 121-180, 181-200, 201-600 Is there an accepted step wise procedure to normalizing the data before I can cluster it. Thanks.


title: Multiple regression or partial correlation?

text: I don't even know if this question makes sense, but what is the difference between multiple regression and partial correlation (apart from the obvious differences between correlation and regression, which is not what I am aiming at)? I want to figure out the following:   I have two independent variables ($x_1$, $x_2$) and one dependent variable ($y$). Now individually the independent variables are not correlated with the dependent variable. But for a given $x_1$ $y$ decreases when $x_2$ decreases. So do I analyze that **by means of multiple regression or partial correlation**? _edit to hopefully improve my question:_ I am trying to understand the difference between multiple regression and partial correlation. So, when $y$ decreases for a given $x_1$ when $x_2$ decreases, is that due to the combined effect of $x_1$ and $x_2$ on $y$ (multiple regression) or is it due to removing the effect of $x_1$ (partial correlation)?


title: R's stats::aov and nlme::lme produce conflicting results

text: I want to use both anova and linear models to test the assumption that my some of my categories have different means than the rest. I am using stats::aov for anova and nlme::lme for linear modelling. The full code is available in this notebook. Basically I end up getting:               Error: ID               Df  Sum Sq  Mean Sq F value Pr(>F)     Residuals  6 0.01198 0.001997                         Error: Within               Df  Sum Sq  Mean Sq F value Pr(>F)       COI        7 0.01926 0.002752    1.92 0.0903 .     Residuals 42 0.06020 0.001433       and               Fixed effects: ER ~ COI                            Value  Std.Error DF   t-value p-value     (Intercept)      0.00200000 0.01465712 42 0.1364525  0.8921     COIemotion-hard  0.05600000 0.02023699 42 2.7672098  0.0084     COIscrambling-06 0.04218045 0.02023699 42 2.0843242  0.0433     COIscrambling-10 0.02094737 0.02023699 42 1.0351029  0.3065     COIscrambling-14 0.00685714 0.02023699 42 0.3388420  0.7364     COIscrambling-18 0.02085714 0.02023699 42 1.0306445  0.3086     COIscrambling-22 0.04885714 0.02023699 42 2.4142493  0.0202     COIscrambling-26 0.02742857 0.02023699 42 1.3553681  0.1825      Which means that lme is telling me 3 groups may significantly lie outside the range of the others, while aov is telling me the probability of ANY GROUP AT ALL being different from the others is not significant (using a p=0.05 cutoff). What am I to make of this? Do I understand these results correctly?


title: Conditional expected lifetime in survival analysis

text: I want to do survival analysis with the Colon Cancer data in survival package:               library(survival)     Colon<-subset(colon,etype==2 & rx!='Obs')     fit<-survreg(Surv(time,status)~rx+age+sex,data=Colon, dist="weibull")     Now how can I predict conditional expected lifetime given time, rx, age and sex?      For example, if a 50 years old, male patient under Levamisole treatment has already survived for 2000 days, what would be his expected life time? I've tried several parameters in predict function, none of them seem to give what I want directly, and I don't know which to choose so I can make some transformation.


title: FE shrinkage estimates

text: In what circumstances can Fixed effects models produce BLUE/shrinkage estimates? I am referring to models discussed in Tekwe et al.(2004)...Journal of Educational and Behavioral statistics. I would appreciate any help on this. Thanks!


title: Real examples of Correlation confused with Causation

text: I'm looking for specific, real cases in which a causal relationship was inappropriately inferred from evidence of a correlation. Specifically, I'm interested in examples that meet the following criteria:   * Existence of the causal relationship was accepted as fact widely enough to have **notable effects** (on public policy, discourse, individual decisions, etc.).   * The link was **inferred solely on the basis of correlative evidence** (perhaps along with the existence of a coherent but unproven causal mechanism).    * Causality has been **objectively falsified** or at least called into serious doubt. The two examples that came to mind for me aren't quite ideal:   1. **Sodium intake and blood pressure:** As I understand it, it has since been determined that salt intake only increases blood pressure in sodium-sensitive individuals. The existence of a valid causal relationship (although not quite the one that was originally accepted) make this example less compelling.   2. **Vaccines and autism:** I may have the background wrong, but I believe this link was surmised on the basis of both correlations and (fraudulent) experimental evidence. This example is weakened by the fact that (fake) direct evidence existed. **Note:** I've seen this similar question: Examples for teaching: Correlation does not mean causation My question differs primarily in that it focuses on **notable, real-world examples** and not on examples in which a causal link is clearly absent (e.g., weight and musical skill).


title: How are confidence intervals for proportions computed under SPSS?

text: I wonder whether any of you can help me. SPSS puts confidence intervals in graphs of frequencies and proportions, but it clearly is not using a normal approximation, i.e. it is not using the formula $\text{CI} = m \pm 1.96 \sqrt{(p(1-p)/n}$. apart from anything else the SPSS interval is not symmetrical around the value. So what is it doing? I could not track down the formula anywhere, not even searching on the web (which usually solves 99% of problems in life, plus or minus 1%). :)


title: Hotelling T^2 test derivation question

text: I am reading about the Hotelling $T^2$ test ( _A primer of multivariate statistic s by Richard J. Harris_ ). It says here that the test can be seen as creating a linear combination of your variables and then obtaining the set of coefficients that maximizes the 'combined t radio'. So I start on the math and I see how the t test is defined as: $t(a) = \frac{a^T(\overline{X}-\mu_0)}{\sqrt{a^TSa/N}}$ When taking the square of this to get rid of the square root it looks like this: $t^2(a) = \frac{N[a^T(\overline{X}-\mu_0)]^2}{a^TSa}$ So far so good. Now we formulate an optimization function like this: $h(a)= N[a^T(\overline{X}-\mu_0)]^2 - \lambda(a^TSa-1) $ So basically the idea here is to maximize $h^2(a)$ subjected to the condition that $a^TSa = 1$. _The lambda term appears as later on the idea is to show that the characteristic root of $NS^{-1}(\overline{X}-\mu_0)(\overline{X}-\mu_0)^T$ is the solution of the optimization problem._ **Now here comes my question/problem:** The next thing that is done in this book is to derive $h(a)/da$. When I try to work this out knowing that for a symmetric matrix $A$: $\frac{\partial x^TAx}{\partial x} = 2x^TA$ I obtain: $h(a)/da = 2Na^T(\overline{X}-\mu_0)(\overline{X}-\mu_0)^T - 2\lambda a^TS$ However, the book and pretty much everywhere else shows this derivative as: $h(a)/da = 2N(\overline{X}-\mu_0)(\overline{X}-\mu_0)^Ta - 2\lambda Sa$ What am I doing wrong? I know that $(\overline{X}-\mu_0)(\overline{X}-\mu_0)^T$ is symmetric (product of a matrix by its transpose always is), and also $S$ is symmetric because it is the estimated covariance matrix.


title: White noise for level, log and log differences data sets

text: I am using eviews 7 and I have 3 data sets for DAX stock market index: level (dax), log (ldax), and log differences (dldax). I need to check whether the error terms of these data sets are white noise or not. **Attempt:** I am estimating an equation $y_t=c+\beta y_{t-1} +v_t$ for each data set. Then, I perform the residuals diagnostics on each. For level data set (dax) looking at the graph of the residuals is enough to see that the error term is not white noise: the variance is not constant. However, I am having difficulties with determining whether the error terms of log (ldax) and log differences (dldax) data sets are white noise or not. The log (ldax) data set is not stationary and the auto correlation is very high while log differences (dldax) data set is stationary and has no auto correlation. I also checked the stationarity of these sets by Augmented Dickey Fuller test. So when I estimate the equation I listed above and check the residuals it looks like the error terms of both data sets are white noise. Here are the graphs for the log: ![log](http://i.stack.imgur.com/GAgM3.jpg)   and for the log difference:   ![log differences](http://i.stack.imgur.com/z0ghM.jpg)   Also the correlogram of the residuals shows no auto correlation for both. So the error terms of both log and log differences are white noise, but this does not make sense because log data set is not stationary. Am I missing something? Thank you for your help. If something is confusing or you need more info let me know. I will make necessary edits.


title: Frequency of time series in R

text: I have 10 minutes-intervaled wait-times data for a coffee shop for 4 weeks between 9am-5pm. I use R's ts for my analysis. What should be the frequency parameter? Is it (48=# of intervals per day) or just 1?


title: For which distributions does uncorrelatedness imply independence?

text: A time-honored reminder in statistics is "uncorrelatedness does _not_ imply independence". Usually this reminder is supplemented with the psychologically soothing (and scientifically correct) statement "when, nevertheless the two variables are **jointly normally distributed** , then uncorrelatedness does imply independence". I can increase the count of happy exceptions from one to two: when two variables are **Bernoulli-distributed** , then again, uncorrelatedness implies independence. If $X$ and $Y$ are two Bermoulli rv's, $X \sim B(q_x),\; Y \sim B(q_y)$, for which we have $P(X=1) = E(X) = q_x$, and analogously for $Y$, their covariance is $$\operatorname{Cov}(X,Y)= E(XY) - E(X)E(Y) = \sum_{S_{XY}}p(x,y)xy - q_xq_y $$ $$ = P(X=1,Y=1) - q_xq_y = P(X=1\mid Y=1)P(Y=1)-q_xq_y$$ $$= \Big(P(X=1\mid Y=1)-q_x\Big)q_y $$ For uncorrelatedness we require the covariance to be zero so $$\operatorname{Cov}(X,Y) = 0 \Rightarrow P(X=1\mid Y=1) = P(X=1)$$ $$\Rightarrow P(X=1,Y=1) = P(X=1)P(Y=1) $$ which is the condition that is also needed for the variables to be independent. So my question is: **Do you know of any other distributions (continuous or discrete) for which uncorrelatedness implies independence?** **Meaning:** Assume two random variables $X,Y$ which have _marginal_ distributions that belong to the same distribution (perhaps with different values for the distribution parameters involved), but let's say with the same support eg. two exponentials, two triangulars, etc. Does all solutions to the equation $\operatorname{Cov}(X,Y) = 0$ are such that they also imply independence, by virtue of the form/properties of the distribution functions involved? This is the case with the Normal marginals (given also that they have a bivariate normal distribution), as well as with the Bernoulli marginals -are there any other cases? The motivation here is that it is usually easier to check whether covariance is zero, compared to check whether independence holds. So if, given the theoretical distribution, by checking covariance you are also checking independence (as is the case with the Bernoulli or normal case), then this would be a useful thing to know.   If we are given two samples from two r.v's that have normal marginals, we know that if we can statistically conclude from the samples that their covariance is zero, we can also say that they are independent (but only because they have normal marginals). It would be useful to know whether we could conclude likewise in cases where the two rv's had marginals that belonged to some other distribution.


title: Survival regression and prediction using median

text: I'm using artificially generated hazard curves (that is, I know the true hazard curve) and Aalen's additive model to fit the covariates. For example, below is an individual's hazard curve and my estimate of it: ![enter image description here](http://i.stack.imgur.com/LC6nJ.png) It's a good fit to my eye, and here is the individual's estimated survival curve (exp of the negative of the above ) ![enter image description here](http://i.stack.imgur.com/3112g.png) Again, this looks fine to me. What happens next though is odd. To make predictions of individuals' lifetimes, I use the median of the survival curve (the intersection of the curve and 0.5). I then plot this prediction against the observed value, and what I was expecting to see was a high correlation between the two, instead I see this: ![enter image description here](http://i.stack.imgur.com/Gk1Di.png) This seems very wrong to me. A few questions I have:   * Is using the median correct? Intuitively I see nothing wrong with it, but my results suggest otherwise.   * Is predicting lifetimes just too noisy of a challenge and this is indeed that best I can do? Has anyone experienced something like this before?


title: Covariance with conditional expectation

text: > Suppose $X$ and $Y$ are random variables, $E(Y^2) < \infty$ and $\varepsilon > = Y - E(Y|X)$ so that $Y = E(Y|X) + \varepsilon$. > > Given that $E(\varepsilon | X) = E(\varepsilon) = 0$, show that > $Cov(\varepsilon , E(Y|X)) = 0$. This question has multiple parts so $E(Y^2) < \infty$ may or may not be applicable in this case. Here's what I tried so far. I used the fact that $Cov(X,Y) = E(XY) - E(X)E(Y)$ and $Cov(X,Y) = E[(X - E(X))(Y-E(Y))]$ and concluded that $Cov(\varepsilon , E(Y|X)) = E(\varepsilon E(Y|X))$ or in other words, $E(\varepsilon E(Y)) = 0$. From there, I guess the only thing I have to show is that: $E(\varepsilon E(Y|X)) = 0$, but I'm having trouble doing this.   Am I going in the right track or is this completely the wrong approach to tackling this problem?


title: Poisson process, time and probabilities

text: An average of **n** birds fly through an area, in an hour, following a Poisson process. (I think this means the hours don't matter; there's no influence in the number of birds that fly through, at different parts of the day, hypothetically. Correct me if I'm wrong.) P1 is the probability that exactly **m** birds fly through between 12:00-14:00 (2 continuous hours). P2 is the probability that exactly **m** birds fly through between 15:00-16:00 and 17:00-18:00 (a total of 2 hours still, but not together). Please notice m is the same for both situations. Are the probabilities of P1 and P2 the same or I'm assuming something that's wrong? Thanks! (Please tell me if you need more information, I can make it up)


title: Sampling and confidence calculations

text: I have 100,000 metal parts on a production line a week. Parts can either be fine, or defective. Random parts are automatically pushed onto a conveyor belt to be checked by a human. We then do a simplistic calculation like: After checking 1000 parts, 5 are defective, 5/1000 = 0.5% defective, so potentially 0.5% * 100,000 = 500 are defective in all. We assume the human checking is flawless. Can anyone give any pointers on the type of statistics required to calculate "how many parts need to be checked to be 90% confident of the defect rate"?


title: How should I compare average store sales change across time?

text: I have what seems like a fairly common business statistics scenario: I need to compare one group of stores to another group of stores and be able to say if their difference in sales is statistically different. For example: Group A ($n_A$ = 30 stores) participated in a promotion and saw an avg sales increase for this month compared to the same month last year of $\bar{x}_A$ and a standard deviation of $s_A$. Group B ($n_B$ = 50 stores) did not participate in the promotion and had avg sales increase of $\bar{x}_B$ and corresponding $s_B$. I realize there are a number of other variables but, in theory, I should be able to say with some certainty that there is or isn't any difference between stores that took or did not participate in a promotion, right? Can I do a standard comparison of means test? Does it make a difference that these stores comprise of the entire population? Or is it not the entire population and I should be looking at average increase for multiple months and multiple years?


title: Limiting distribution of a squared sum of random variables

text: Let $S_n = \frac{1}{n}\sum_{i=1}^n X_i$, and $T_n = \frac{1}{n}\sum_{j=1}^nY_i$, where The $X_i$ are iid, the $Y_i$ are iid (with a different law) $X_i$, and $Y_i$ are dependent For $i\neq j$, $X_i$ and $Y_j$ are independent. Is there a central limit type result for $S_n^2 - T_n^2$?


title: What is the closed form solution for the inverse CDF for Epanechnikov

text: Is there a closed form solution for this inverse CDF?


title: Topic and subject classification

text: I have a set of documents that are OCR-ed and represented as a text file. I want to find out what are the documents that are talking about the same subject and maybe about the same person. I started to look if this problem has been addressed in webpages. For example, I'd like to search for all the webpages that talk about the subject "terrorism." It doesn't have to mention the word "terrorism" but yet it should be implicit. Do you have any idea on where I can start my review?


title: coefficients and p-value in logistic regression

text: In logistic regression, I have a variable with larger coefficient and larger p-value and another variable with smaller coefficient and smaller p-value. If use p-value then the latter one is more significant, but if calculate odds ratio, the first one is more influential. How should I interpret this? Which one is more important in the model? Thanks! High: coef:-0.0153869 p-value:0.0000257266774040 Low : coef:-0.0052793 p-value: < 0.0000000000000002


title: Normalized Data vs. Non Normalized Data in SVM Training

text: Im training an SVM on a dataset containing 17000 events with 30 features. as I see it, this amount of features against the dataset size should suffice by itself for avoiding overfitting. Im training an RBF kernel and I get somewhat strange results. Im pretty sure my training flow is correct, so it is not the main issue. The main issue is the difference between normalized to non normalized data and it's effect on the training error. Here is what I do: Im training a classifier with 'rbf' kernel, C is 16 and gamma is 8. Im checking the training error by plotting the ROC on the training set. Im doing this once for the original data, and once for the normalized data (I normalize by removing the average and dividing by (max-min) for each column) I get the following graphs: For normed data: ![Normalized Data](http://i.stack.imgur.com/fXfHJ.png) For the original Data: ![None Normed Data](http://i.stack.imgur.com/576C6.png) This is quit different. Any ideas why it is like that?


title: plot effects on glmer

text: I done a GLMM with the function glmer form lme4, the variable paternity is binomial (yes or no)(model<-glmer(paternity~factor1, factor 2 ....,family=binomial). My model run correctly, but after I want plot the effect. So I used the package effects. When I plot my effect, the Y-axis is numeric but I do not understand what it mean because this axis is binomial (paternity yes or no)... How I must interpreted this ? thank you for your help


title: Fitting multiplicative model in a Frequentist way in R

text: Say, I have sensory panel data for 30 products, which are scored by 11 panelist in 2 replicates for different traits. To fit this data,I am using model of the form: yijr = alpha_j + beta_j*theta_i + error_ijr i= panelists, j= products, r= replicates In this model, I am considering multiplicative interaction for beta_j(for panelists) and theta_j(for products). Error term is heterogeneous and panelist-dependent. I am trying to fit this model in R, so I can get separate coefficients for alpha_j, beta_j, theta_i and error_ijr. Does anyone can suggest me is it possible to fit this model in frequentist way and what would be appropriate R package or function to do these. Thanks


title: Stationary and ADF for unit root

text: If I calculate the returns in my time series as follows, $\ln(P_t / P_{t-1})$, where $P_t$ is the price at time _t_ and $P_{t-1}$ is the price at $t-1$ then I do the ADF test, my $H_0$ is rejected. So my conclusion is that I do not have a unit root, so my time series have stationarity properties; is this conclusion correct? If yes, I will use a mean-reverting process, like: $d(\log \bar{P}_{0,j}(t))=dX_i(t)=k_i (\delta -X_i (t))dt+ \sigma dW_i^1 (t)$ I am not quite sure the conclusions make sense. So first, taking returns, and rejected the ADF $H_0$, can I conclude that my time series has stationarity properties? Second, because of the stationarity properties, my time series will in the end go back to the long-term mean, therefore a mean-reverting process to describe the dynamics of the returns is the right process. Is that right? So, are the steps that I take (from returns $\rightarrow$ reject ADF test $\rightarrow$ conclude stationarity properties $\rightarrow$ model stationarity properties using a mean-reverting model) correct and logical?


title: Logistic Regression Performance on training data set V/s AIC

text: I am fitting a logistic Regression on data set having 700 variables (after Chisquare test) and 15000 rows. For that I did best subset analysis using glmulti package in R on first 70 variables and got a subset based on min AIC, after that I added manually variables to my Model, this I repeated for another 30 variables. During the latter analysis sometimes AIC values of model decreases but the performance on test data remains same. So should I keep on like this or should I stop here only ? Is it okay to select variable like this ?


title: Using PyMC to obtain the posterior for the parameters of a parametric model?

text: I am a newbie with **pyMC** and I want to program an **MCMC** sampling for a complicated problem. I have a function given by in the following: $$g(r,c,\theta_1,\theta_2)=\frac{\delta_c}{\Sigma x^2}\log{\frac{x}{2}}$$ The parameter $\Sigma$ and $\delta_c$ and $x$ are given by $$\Sigma(z)=\frac{const*z+\sqrt z}{(1+z)^2}$$ $$\delta_c=\frac{c}{\log(1+c)-c}$$ $$x=\frac{c\sqrt(\theta_1^2+\theta_2^2)}{r}$$ The free parameters which I would like to find their posterior probabilities are $r$, $c$ and $\theta_1$ and $\theta_2$, while for $r$ the prior is **exponential** , $c$ has **lognormal** prior and $\theta_1$ and $\theta_2$ have **uniform** priors. My **main problem** is how to define the model, because for each of point in data set I have a probability distribution for parameter $z$ which I would like to **marginalize** over and I don't know how it should be done. I have written so far the following piece of code, I have written a **class** for $g$ :               import pymc as pm     import g_class     import numpy as np     import math     import random     from scipy.stats import expon     data=np.loadtxt('test.cat') #my data     for i in range(len(data[:,0])):     if (i<1):         theta=np.array([data[i,0],data[i,1]])         observed_g=np.array([data[i,2],data[i,3]])     else:         theta=np.vstack((gal_pos, np.array([data[i,0],data[i,1]])))          observed_g=np.vstack((observed_g, np.array([data[i,2],data[i,3]])))           z_pdf=np.loadtxt('pdf.z.dat') #probability distribution of z in the given range     z=np.arange(0,1.5,0.001)          @pm.stochastic(dtype=np.float, observed=False, trace=True)     def thata1(value=235,t1_l=210,t1_h=270):         if ((value>t1_h) or (value<t1_l)):            return -np.inf         else:            return -np.log(t1_h-t1_l+1)                  @pm.stochastic(dtype=np.float, observed=False, trace=True)     def theta2(value=340,t2_l=320,t2_h=360):         if ((value>t2_h) or (value<t2_l)):            return -np.inf         else:            return -np.log(t2_h-t2_l+1)                r=10*pm.Exponential('r', beta=10)               @pm.deterministic  #     def sigma(value = 1, r=r):         if r < 10**15:            return .09        else:            return .06          cExpected = (r/10)**(-.1) # parameter c and r related via this formula, I use this to compute exoected value for log(c)     conc = pm.Lognormal("concentration", cExpected, sigma)     #model     @pm.deterministic( name='g', observed=True, trace = True )     def g(r=r,c=conc,z=zh,pos=[thata1,thata2],theta):         nfw = g_class.g(r=r,c=conc,z=zh,pos=[thata1,thata2])         value1=0;value2=0         n=theta[0]         ID=[i for i,l in enumerate(data) if n in l]         for i in range(len(z)):             g1,g2,mag=g_class.getG( theta, z[i])             value1+=g1*z_pdf[ID,i]/sum(z_pdf)#first component              value2+=g2*z_pdf[ID,i]/sum(z_pdf)         return np.array([value1,value2])               std_dev = pm.Normal('tau', mu=0.0 , 0.3)     #likelihood     obs = pm.Normal("obs", mu=g, tau=1./std_dev**2, value=observed_g, observed=True)      The problem is that when I define my model as **Deterministic** and use observation to construct the model, I get the error message that `theta` is not defined properly as a **parent**. I am wondering, how could I include `theta`?


title: Comparing to vectors of labels that contain different labels, but the same number, and the same places

text: Problem: After doing two cluster analysis, I get two arrays $V_1$ and $V_2$ of length $N$. I have a groundtruth vector as well $GT$. I want to make the comparison, which is label invariant. Say first cluster analysis constrains two labels blue and yellow, and the second cluster analysis come up with 1 and 0. I want a comparison that takes these values, and looks for if the element white is contained the same places and 1 is, and so forth. My idea: Make the cross-product $ A = V_1 \times GT$, which will consist of a matrix, where each element is $$a_{i,j} = (v_i, gt_i) \qquad v_i \in V, gt_j \in GT $$ for $ i = 1,..., N$. And then run some compare function of two variables over each element in this matrix. I cannot make the cross product in Matlab, furthermore I am unable to think of how to yield any information from this pairwise comparison matrix. Anyone would like to assist?


title: necessary condition for the existence of the posterior mean

text: I wonder if there exists results relating properties of the prior and the likelihood to the existence of the posterior mean (or more generally to the posterior moments). Any hint or sufficiant condition will be also appreciate.


title: Where to start in estimating the days until payment of an invoice?

text: I am in the process of modeling the random variable X as the days from issuing an invoice until its payment. This variable is dependent on the days of credit for the invoice, so there are distinct X's depending on the days of credit, let's call them X_c. I'm starting with the most common, X_30, which should be centered around 30 with a very heavy tail and a very rapid ramp up of the density getting to 30. A basic histogram looks like this: ![alt text](http://i.stack.imgur.com/p10ia.png) I have lots of data to try and fit but would like to get some pointer into which distributions might model this better from a conceptual point of view.


title: Sample selection and variation in the variable of interest when using panel data

text: I wonder how to reason when selecting samples. I am doing a panel data regression analysis about how the euro membership correlates with the budget deficit in the member countriues. Using R that is:               BUDGETDEFICIT ~ EURODUMMY + some_controll_variables      EURODUMMY is the variable of interest. Using panel data I have anual observations for many countriies. My question is how selecting what sample to use effects the regression analysis. Should I only use countries that has introduced the Euro? Is adding other European countries or western like countries benefitial as they add more datapoints and/or acts like a reference gropup? Or is id only bad for the data to add countries that has not introduced the euro as they have no variation in the variable of interest (EURODUMMY)? Other comments regarding this problem is also welcome.


title: Are Auto-Associative Regression Trees Distinct from Auto-Regressive Trees?

text: After some reading in the field I was confused as to whether these two models are distinct or really the same. I'm just looking for a simple yes/no with a brief explanation. Note that _Auto-Associative_ is a term used mostly in the context of Neural Networks, only rarely in the context of decision or regression trees.


title: Tobit model with negative values

text: I have a dependent variable which I would like to estimate using a tobit model. The dependent variable measures change between two points in time, and thus can take on either positive or negative values. Some observations choose not to participate in the process generating the DV at all, so they are observed with a value of 0. The same independent variables determine whether to participate (not 0) as well as change in the DV. How does the tobit model handle negative values? An alternative could be to model the DV as the end value, and include the starting value as a control variable.


title: Modeling user behavior in a system

text: Consider a system, where each user **enters** the system, performs a series of predefined **actions** , and then **exits**. For instance, consider a system with 5 predefined action. The action log of some user is detailed below: > **enter** , **action 3** , **action 1** , **action 5** , **action 1** , > **action 4** , **exit**. Notice that the order of actions matters. Assume that we have the action logs for thousands of users, and we want to graphically model their behavior.   1. What is the proper model, and    2. How can we translate the action logs into this model? * * * **More info.** The Markov chain does not seem to be right, as Markov chains are memoryless, while in our system, the user remembers his/her actions, and decides on the next action based on the complete history after entering the system. I guess Bayesian networks might be the appropriate choice, but I'm not sure. Specifically, I'm not sure if the acyclicity of Bayesian networks imposes a restriction on the desired model.


title: Introductory materials in statistical analysis and data visualisation

text: To get started on statistical analysis -- where would you recommend I go to learn?


title: Mixed model, repeated measures, nested?, nominal data

text: I am not sure how to approach the statistical analysis of the following study. Any insights would be greatly appreciated. Thanks! So, I have two independent groups, group A and group B. I am collecting data for 3 different tasks (Task 1, Task 2, Task 3) at 3 different time points, Pre, During and Post. There are 3 trials for each of the tasks, so I have 9 observations per subject per task (3 pre, 3 during and 3 post per task), 27 observations in total per subject. The problem is that the dependent variable, which is Pattern P, is nominal (each pattern coded with numbers from 0-5). Group A Task 1: Pre1, Pre2, Pre3, During1, During2, During3, Post1, Post2, Post3 Task2: Pre1, Pre2, Pre3, During1, During2, During3, Post1, Post2, Post3 Task 3: Pre1, Pre2, Pre3, During1, During2, During3, Post1, Post2, Post3 Group B Task 1: Pre1, Pre2, Pre3, During1, During2, During3, Post1, Post2, Post3 Task2: Pre1, Pre2, Pre3, During1, During2, During3, Post1, Post2, Post3 Task 3: Pre1, Pre2, Pre3, During1, During2, During3, Post1, Post2, Post3 I am interested in pattern differences between the 3 time conditions, between the 3 different tasks and between the 2 groups. I would normally do a repeated measures MANOVA in this case, with 2 factors (Task with 3 levels and Time with 3 levels) Within-Subjects and SubjectGroup as a Between-Subjects factor, but in this case I have 3 observations and not 1. Since the measure is nominal, I can't take the mean of the 3 trials. Perhaps, treating each trial as nested would solve the problem? If yes, can I do this with repeated measures in SPSS? Or another type of analysis would be preferable? Any suggestions? Many thanks!


title: Multiple regressions in one proc statement

text: I am attempting to run a mixed model regression multiple times on various outcomes. As opposed to copying and pasting a new proc mixed statement for each outcome how can I write a statement to cycle through it all. I took a first shot at it, but this fails to work.               data vitd.level;set vitd.level;     array outcome A--OOO;     do i=1 to dim(outcome) by 1;           PROC MIXED DATA = vitd.level; *method = ml covtest ic;                  CLASS  Protocol Time;                        Model outcome{1} = Protocol Time Protocol*Time;                       REPEATED Time/SUBJECT = Patient type=arh(1);                       LSMEANS PROTOCOL*TIME;      end; run; Any one know a way to do this? -K


title: What is a vector autoregressive model?

text: I'm looking to understand this from a managerial perspective. For example if I was explaining linear regression I would say it is a line of best fit through some data points and it can be used to predict a "y" value for some given value of "x". Is there an analogous explanation for VAR? I do not have a strong background in statistics.


title: What level to use when comparing subjects in a hierarchical Bayesian analysis?

text: Say that I have an experiment where I test the reaction time of a number of subjects where each subject makes many reaction time trials. In a Bayesian framework the reaction times ($y$) could be modeled by a hierarchical model with prior distribution both on the subject level and for the whole group of subjects. A diagram of the model, Kruschke style, could be: ![Model Diagram](http://i.stack.imgur.com/xb3Ks.png) ... and the corresponding BUGS/JAGS code would be:               for(i in 1:length(y)) {       y[i] ~ dnorm(mu[subj[i]], tau[subj[i]])     }          for(j in 1:nbr_of_subjects)       mu[subj[i]] ~ dnorm(M_mu, P_mu)       tau[subj[i]] ~ dgamma(S_tau, R_tau)     }          M_mu ~ dnorm(M_M, P_M)     P_mu ~ dgamma(S_P, R_P)          S_tau <- pow(m , 2) / pow(sd, 2)     R_tau <- m / pow(sd, 2)     m ~ dgamma(S_m, R_m)     sd ~ dgamma(S_sd, R_sd)      If I wanted to compare the reaction time of two subjects I would then compare the their respective $\mu$ distributions. If the reaction time trials were split up into four blocks I could also model that by adding an extra block level with priors between the subject level and trial level in the diagram (as it might be the case that the subjects reaction time differs slightly between blocks for some reason). My question is now, if I would want to compare two subjects what distributions should I compare? I could compare the distribution of the means on the subject level (which now partly defines the prior for the mean on the block level) but I could also compare the distribution of the means on the block level which corresponds to $\mu$ in the old model. In one way it seems more logical to compare the subjects on the subject level, but does it make any difference? And if there are very few blocks, say two, wouldn't the distribution of the means on the subject level be very "wide"?


title: Using QAICc with Poisson, or AIC with Poisson lognormal, in information theoretic approach?

text: I am trying to use an IT approach to analyse some ecological data. I have a mixed model with nested random effects (I'm using glmer in package lme4 in R). I initially fit the model with a Poisson error distribution, but it was overdispersed so I added an observation-level random effect to deal with the overdispersion. Model check plots are then very good. However, my next step in the IT approach is to generate a set of candidate models (using function 'dredge'). I have two choices now, and I don't know which is correct / best: Choice 1. Start with the Poisson lognormal model (i.e. I have already dealt with the overdispersion) and generate a model set using AICc. Choice 2. Start with the Poisson model, which is overdispersed, and generate a model set using QAICc and specifying the dispersion parameter. Is it acceptable to do either of these, or should I stick to Choice 2, which is what the example code does in the MuMIn package help documentation? (However, if I go for choice 2, then the model check plots of the global model aren't great, and yet there's no way of checking the averaged model for goodness of fit - so how do I know the model fit is OK?). Thank you very much for any help you can give.


title: Randomly given a coin from an unfair set of 6 coins, how can I determine which coin I was given?

text: Suppose I have 6 different coins that have wildly different increasing chances of flipping a heads, and I know the chances of flipping a heads for each. (ie. Coin 1: 1/280 chance of heads ... Coin 6: 1/120 chance of heads) Without looking, if I choose one of these coins at random and flip it some limited number of times, how can I determine the probability it is each of the 6 coins in relation to each other from the results alone? (ie. A result like [Coin1: 5%, Coin2: 7%, Coin3: 10%, Coin4: 14%, Coin5: 25%, Coin6: 39%]) My statistics knowledge is limited to what I learned in first year university, so I originally thought I could solve this problem using the chi-square test with the expected results of each of the 6 coins. The problem is this gives me the probability that the results were simply due to chance vs it not actually being that coin (possibly not even one of the original 6 coins, when I know it must be one of the 6 coins!). So for example if I flipped a coin 10,000 times and (by some extreme chance) every result was a heads, I should be getting a result reporting the confidence of coin 6 as being something extremely high with a very tiny % chance of being coin 1 to 5. Instead the result I get is that there's a 0.000001% or so chance of it being any of the 6 coins at all. Am I using the wrong approach to solving this problem?


title: I want to show $M(t)=e^{\frac{t}{2}}\sin(B(t))$ is a martingale by using Ito s formula

text: Show that $M(t)=e^{\frac{t}{2}}\sin(B(t))$ is a martingale by using Ito s formula. $B(t)$ is brownian-motion. i must show $s \le t $ then $E(M(t)|\mathcal F_{s})=M(s)$ but i dont know how use ito formula   thanks for help


title: Expectation Propagation and multivariate priors

text: I have been struggling with this for weeks now. All EP examples that I have found on the net seem to deal with univariate priors and I am really at a loss as to how to make it work with a multivariate prior. The scenario is as follows: I am trying to estimate a bunch of parameters where the prior over these parameters form a multivariate Gaussian with off diagonal covariance terms. The likelihood can be written down as a product over individual likelihood for each of these parameters. So, the equation is as follows: $$ P(w|y) = \frac{1}{Z}P(w)\prod_iP(y_i|w_i) $$ where $y$ is the observed data and $w$ are the parameters I am trying to estimate. So, when I want to perform inference on this using EP, what should be done with $P(w)$. $P(w)$ is a multivariate Gaussian with 0 mean and a dense covariance matrix. So, at some point I will have to multiply terms between $P(w)$ and $P(y_i|w_i)$ or some similar approximates of it (but the forms are very different). So effectively to compute the functional form of the posterior, I am multiplying a multivariate Gaussian ($P(w)$) with a univariate Gaussian function ($\prod_i{q(y_i|w_i)}$) where $q(y_i|w_i)$ are the approximating factors. How can this be achieved? Does $P(w)$ need to be broken down as well? I am not sure how to make this work.


title: Boxplots as tables

text: We are preparing a manuscript and the editor asked us to convert a figure with a boxplot into a table "because of the more exact data content". While I think boxplots are pretty decent in revealing something about the data, I was wondering what you guys think about this? Do you (often) opt for tables over boxplots for presenting your data?


title: Binomial Response Difference Between Samples

text: I have a data set of 500 companies which have been clustered into groups based on their answers to survey questions. I'm trying to gauge whether the profitability between groups is statistically different, and if so what the correct test is to use. Profitability is indicated as a binary variable (1,0).   1. Sample      * n = 500     * profitable = 203   2. Cluster 1      * n = 198     * profitable = 78   3. Cluster 2      * n = 159     * profitable = 55   4. Cluster 3      * n = 143     * profitable = 70 So far I have used the binomial propotions test to determine whether the sample mean statistically different than the cluster mean. Is this the correct route?


title: Interaction effect in stepwise regression

text: I am trying to creat a multiple regression model with a forward stepwise procedure. Predictors are air temperature, soil temperature, PAR and snow depth. I also want to see if there are some interaction effect between snow depth and other factors that can also be included in the model. Should I include the interaction as factors (e.g. Ta:Snow depth) when going through the stepwise procedure (main effect may not be included) or make another model to illustrate the interaction effect separately (include both main and interaction effect)? PS:I am using R to do this. Thanks a lot in advance!


title: How to test hypothesis  that correlation is equal to given value using R?

text: Is there a function to test the hypothesis that the correlation of two vectors is equal to a given number, say 0.75? Using cor.test I can test cor=0 and I can see whether 0.75 is inside the confidence interval. But is there a function to compute the p-value for cor=0.75?               x <- rnorm(10)     y <- x+rnorm(10)     cor.test(x, y)


title: Can I use a difference score as my outcome variable (pre-post change) in a moderated multiple regression equation?

text: I am wanting to examine whether religiosity moderates intervention effects on stigmatized attitudes. Here are my variables: X = group (1 = experimental; 0 = control) Z = religiosity (14-item scale - using total score) Y = either post- test scores or pre-post change scores - Social Distance; Genderism/Transphobia I realize that I need to center my moderator variable and will also need to center my outcome variable (if I use post-test scores). My question is whether to use post-test scores or to use pre-post test change scores (and if I use the latter would I still center)


title: Unit root tests for panel data in R

text: I have the `plm` package and would like to run unit root tests on some variables. I get the following error:               > purtest(data$tot.emp)     Error in data.frame(baldwin = c(59870, 61259, 60397, 58919, 57856, 57227,  :        arguments imply differing number of rows: 14, 19, 11, 12, 1, 20, 18, 10, 13      I assume that I'm getting this error because my panel is unbalanced. Two questions:   * Can you use panel unit root tests (Levin, Lin and Chu (2002), Im, Pesaran and Shin (2003), or others) for unbalanced panels?   * If so, is it implemented in R?


title: Estimating the influence of different features on the outcome

text: I been trying, with no luck, to find the correct algorithm for the following 2 scenarios and I can't seem to get it right. **First scenario** Every day I get data like the following:               +---------+----------+----------+----------+----------+---------+     | day     | keyword1 | keyword2 | keyword3 | keyword4 | success |     |         | clicks   | clicks   | clicks   | clicks   |         |     +---------+----------+----------+----------+----------+---------+     | day1    | 10       | 3        | 5        | 9        | 76      |     | ...     | ...      | ...      | ...      | ...      | ...     |     +---------+----------+----------+----------+----------+---------+      Success is a value which measures how well that day's 'clicks' performed on different values, conversions, etc. Is there any algorithm I can use to assign a weight to each keyword in relation to the success across various days? I though of association rules and decision trees but I can't seem to get the idea right on how would those help me. **Second scenario** ... which is pretty similar were data has the following structure               +------------+---------+---------+---------+---------+     | importance | value 1 | value 2 | value 3 | value 4 |     +------------+---------+---------+---------+---------+     |          1 |      18 |      21 |      35 |      25 |     |          2 |      93 |      36 |      11 |      56 |     |          3 |      34 |      26 |      47 |      47 |     |          1 |      19 |      20 |      10 |      23 |     |          1 |      17 |      20 |       3 |      25 |     +------------+---------+---------+---------+---------+      In this case what I am trying to do is understand how different values affect the `importance` value. In the previous case you can easily see that for `importance=1` value 1, value 2, and value 4 are "close" to one another while value 3 is not.


title: Why is the Fisher exact test not always used in the cases where a chi squared is appropriate?

text: According to Wikipedia: > Although in practice [Fisher's exact test] is employed when sample sizes are > small, it is valid for all sample sizes. If Fisher's exact test can provide the exact p-values in lieu of the approximation that is given by the chi squared, what is the reason to ever use the $\chi^2$ test?


title: What is Shannon's source entropy?

text: Suppose that ${X_n; Y_n}$ is a random process with a discrete alphabet, that is, taking on values in a discrete set for $n$ data length. They correspond to the input and output of a communication process. Assuming Y to be the discretized output, what is meant by Shannon's source entropy for a given discretization bin $k$ : $H_{source} = lim_{k \rightarrow \infty} \frac{1}{k} H_k$ where $H_k$ stands for Shannons entropy. I found this in Stojanovski et al's paper: Applications of Symbolic Dynamics in Chaos Synchronization. Now, Jun Chen also mentions about source entropy but the formula is very different! Its the same as Shannon's entropy. So, what exactly is _source entropy_? Is it the entropy of $X$ or $Y$ or is it the same as Shannon's entropy?


title: ANOVA for non-normal heterogenous unequally replicated data

text: I need some advice on how to proceed with my data analysis. I have 3 groups (Archaea, Bacteria and Eukarya). Each group has unequal number of individual species (70, 651, and 244 respectively). Each specie was given 7 different treatments (A, B, C, D, E, F, and G) and I counted the number of domains (you can call it 'y'). Now the problems are: Counts for treatment A are extremely high compared to counts of other treatments. Treatment G has the lowest counts. Other treatments have more or less similar counts. Data is non-normal, does not have homogenous variances, unequally replicated and any help would be appreciated.   * Which test should I use?    * Any transformation?


title: Eigenfunctions of an adjacency matrix of a time series?

text: Consider a simple time series:               > tp <- seq_len(10)     > tp      [1]  1  2  3  4  5  6  7  8  9 10      we can compute an adjacency matrix for this time series representing the temporal links between samples. In computing this matrix we add an imaginary site at time 0 and the link between this observation and the first actual observation at time 1 is known as link 0. Between time 1 and time 2, the link is link 1 and so on. Because time is a directional process, sites are connected to (affected by) links that are "upstream" of the site. Hence every site is connected to link 0, but link 9 is only connected to site 10; it occurs temporally after each site except site 10. The adjacency matrix thus defined is created as follows:               > adjmat <- matrix(0, ncol = length(tp), nrow = length(tp))     > adjmat[lower.tri(adjmat, diag = TRUE)] <- 1     > rownames(adjmat) <- paste("Site", seq_along(tp))     > colnames(adjmat) <- paste("Link", seq_along(tp)-1)     > adjmat             Link 0 Link 1 Link 2 Link 3 Link 4 Link 5 Link 6 Link 7     Site 1       1      0      0      0      0      0      0      0     Site 2       1      1      0      0      0      0      0      0     Site 3       1      1      1      0      0      0      0      0     Site 4       1      1      1      1      0      0      0      0     Site 5       1      1      1      1      1      0      0      0     Site 6       1      1      1      1      1      1      0      0     Site 7       1      1      1      1      1      1      1      0     Site 8       1      1      1      1      1      1      1      1     Site 9       1      1      1      1      1      1      1      1     Site 10      1      1      1      1      1      1      1      1             Link 8 Link 9     Site 1       0      0     Site 2       0      0     Site 3       0      0     Site 4       0      0     Site 5       0      0     Site 6       0      0     Site 7       0      0     Site 8       0      0     Site 9       1      0     Site 10      1      1      The SVD provides a decomposition of this matrix into Eigenfunctions of variation as different temporal scales. The figure below shows the extracted functions (from `SVD$u`)               > SVD <- svd(adjmat, nu = length(tp), nv = 0)      ![Eigenfunctions](http://i.stack.imgur.com/MGUOx.png) The eigenfunctions are periodic components at various temporal scales. Trying `tp <\- seq_len(25)` (or longer) shows this better than the shorter example I showed above. _Does this sort of analysis have a proper name in statistics?_ It sounds similar to Singular Spectrum Analysis but that is a decomposition of an embedded time series (a matrix whose columns are lagged versions of the time series). **Background:** I came up with this idea by modifying an idea from spatial ecology called Asymmetric Eigenvector Maps (AEM) which considers a spatial process with known direction and forming an adjacency matrix between a spatial array of samples that contains 1s where a sample can be connected to a link and a 0 where it can't, under the constraint that links can only be connected "downstream" - hence the asymmetric nature of the analysis. What I described above is a one-dimensional version of the AEM method. A reprint of the AEM method can be found here if you are interested. The figure was produced with:               layout(matrix(1:12, ncol = 3, nrow = 4))     op <- par(mar = c(3,4,1,1))     apply(SVD$u, 2, function(x, t) plot(t, x, type = "l", xlab = "", ylab = ""),           t = tp)     par(op)     layout(1)


title: Prediction interval = credible interval?

text: I am wondering if prediction interval and credible interval evaluate the same thing. For instance with a linear regression, when you estimate the prediction interval of a fitted values, you estimate the $(1-\alpha)\%$ limits of the interval in which you expect your value to fall. Conversely to a confidence interval, you do not focus on a distribution parameter such as the mean value, but on the value that your explained variable could take for a given X value (supposing that $\ Y = a + b.X$). When you estimate the fitted value for a given $X$ value within a Bayesian framework, from the posterior probability distribution, you can estimate a credible interval. Does this interval give you the same information on the fitted value or not?


title: Are models identified by auto.arima() parsimonious?

text: I have been trying to learn and apply ARIMA models. I have been reading an excellent text on ARIMA by Pankratz - _Forecasting with Univariate Box - Jenkins Models: Concepts and Cases_. In the text the author especially emphasizes the priciple of parsimony in choosing ARIMA models. I started playing with `auto.arima()` function in **R** package **forecast**. Here is what I did, I simulated ARIMA and then applied `auto.arima()`. Below are 2 examples. As you can see in both example `auto.arima()` clearly identified a model that many would consider non-parsimonious. Especially in example 2, where `auto.arima()` identified ARIMA(3,0,3) when actually ARIMA(1,0,1) would be sufficient and parsimonious. Below are my questions. I would appreciate any suggestions and recommendations.   1. Are there any guidance on when to use/modify the models identified using automatic algorithms such as `auto.arima()`?   2. Are there any pit falls in just using AIC (which is what I think `auto.arima()` uses) to identify models?   3. Can an automatic algorithm built that is parsimonious? By the way I used `auto.arima()` just as an example. This would apply to any automatic algorithm. **Below is Example #1:**               set.seed(182)     y <- arima.sim(n=500,list(ar=0.2,ma=0.6),mean = 10)          auto.arima(y)          qa <- arima(y,order=c(1,0,1))     qa      Below are the results from `auto.arima()`. Please note that all the coefficients are insignificant. i.e., $t$ value < 2.               ARIMA(1,0,2) with non-zero mean           Coefficients:              ar1     ma1      ma2  intercept           0.5395  0.2109  -0.3385    19.9850     s.e.  0.4062  0.4160   0.3049     0.0878          sigma^2 estimated as 1.076:  log likelihood=-728.14     AIC=1466.28   AICc=1466.41   BIC=1487.36      Below are the results from running regular `arima()` with order ARIMA(1,0,1)               Series: y      ARIMA(1,0,1) with non-zero mean           Coefficients:              ar1     ma1  intercept           0.2398  0.6478    20.0323     s.e.  0.0531  0.0376     0.1002          sigma^2 estimated as 1.071:  log likelihood=-727.1     AIC=1462.2   AICc=1462.28   BIC=1479.06      **Example 2:**               set.seed(453)     y <- arima.sim(n=500,list(ar=0.2,ma=0.6),mean = 10)          auto.arima(y)          qa <- arima(y,order=c(1,0,1))     qa      Below are the results from `auto.arima()`:               ARIMA(3,0,3) with non-zero mean           Coefficients:              ar1      ar2     ar3     ma1     ma2     ma3  intercept           0.7541  -1.0606  0.2072  0.1391  0.5912  0.5491    20.0326     s.e.  0.0811   0.0666  0.0647  0.0725  0.0598  0.0636     0.0939          sigma^2 estimated as 1.027:  log likelihood=-716.84     AIC=1449.67   AICc=1449.97   BIC=1483.39      Below are the results running regular `arima()` with order ARIMA(1,0,1)               Series: y      ARIMA(1,0,1) with non-zero mean           Coefficients:              ar1     ma1  intercept           0.2398  0.6478    20.0323     s.e.  0.0531  0.0376     0.1002          sigma^2 estimated as 1.071:  log likelihood=-727.1     AIC=1462.2   AICc=1462.28   BIC=1479.06