title: Appropriate tests for finding a diffference between the scores of two groups of questions

text: I am testing whether there is difference between the scores of 2 groups of 12 questions on a Likert-scale questionnaire. Each question from one group has a corresponding question in the other group (12 pairs altogether) and my hypothesis is that the first group of questions has a higher score than the second one(one-tail). I've already entered my data into SPSS and I need help figuring out what would be the best tests to run (I am also curious about the difference in the score of individual pairs etc.). What would be most appropriate and what else I might want to look at? For what is worth I am thinking of a t-test or mann-whitney but I need help.


title: Test heteroskedasticity in Tobit model

text: I'm examing the relationship of payout policy and cash-flow uncertainty. My research' s data is cross section data and I use Tobit model to test the impact of explanatory variables (cash-flow uncertainty) to dependent variables (payout policy), and Logit model to test the probability of paying dividends. So when I use Tobit and Logit model with cross section data, I have to test the heteroskedasticity or not? It's said that there were two basic assumptions underlying the tobit model. And one of two basic assumptions is Tobit model' s data can be heteroskedasticity. That's right? Please hep me, thanks a lot


title: Analogous measure of AIC which uses the posterior distribution for model selection?

text: Suppose the following problem: I have $n$ models, $M_k$, each with parameters $\mathbf{\theta}_k$ for a data set $D$. There where previous observations of a subset of the parameters which are common to every model $M_k$ (i.e., I have well defined priors for a subset of the parameters $\theta_k$), so I performed an MCMC algorithm in order to obtain the posterior distribution of each model using that prior information, i.e., I have $p(\theta_k|D,M_k)$, and have to decide which of those models is the 'correct' one. I was thinking in defining what do I mean by 'the correct' one, and came up with the idea that I have to decide which of the posterior distributions is closer to the 'real' posterior distribution that generated the data (which may or may not be in my set of posterior distributions). I was thinking of using bayes factors, but I keep thinking that I need something like the AIC which, instead of using the likelihood and the corresponding MLE estimates, uses the posterior distributions and the corresponding maximum-a-posteriori estimates. My idea is to obtain an unbiased (or nearly unbiased) estimator of the KL divergence between the real posterior and my posteriors (understanding that the AIC is an estimator of the KL divergence between the 'real' likelihood and the likelihood of my models). Is there something like this in the statistical literature? I'm just kind of crazy of thinking the problem like this?


title: Specifying multiple (separate) random effects in lme

text: I was working in R packages **nlme** and **lme4** , trying to specify the models with multiple random effects. I found, that only **nlme** allows to specify the heterogeneous structure of the variance. Therefore, I got a model, where temperature (Y) depends on time (in hours), intercept varies by date and year, and variance also varies by year:               fit1 <- lme(Y ~ time, random=~1|year/date, data=X, weights=varIdent(form=~1|year))      However, if I need to add another random term (time varying by date), and specify the model like this:               fit2 <- lme(Y ~ time, random=list(~1|year, ~time-1|date,  ~1|date), data=X,                  weights=varIdent(form=~1|year))      the random effects become nested in each other: date in year; and then date in date and in year. I also tried               one  <- rep(1, length(Y))     fit3 <- lme(Y ~ time, random=list(one=pdBlocked(list(pdSymm(~1|year/date),                  pdSymm(~time-1|year)))), data=X, weights=varIdent(form=~1|year))      but it gives an error:               Error in pdConstruct.pdBlocked(object, form = form, nam = nam, data = data,  :       cannot have duplicated column names in a "pdMat" object      I understand that there have been already many questions related to the similar problem, but I really did not find the answer for my case. Could you help me with the right specification of the model?


title: show the variance of the estimator

text: I should show that the variance of the least squares estimator in the simple regression model $y_t = b x_t+u_t$ with the variables in deviation form is $\operatorname{Var}(\hat{b}) = \operatorname{E}[(\hat{b} - b)^2] = \sigma^2 / \sum(x_t^2)$ ($t$ is the time based indicator) I know that $\hat{b} - b = \sum(x_t*u_t) / \sum (x_t^2)$* but I cannot solve it to get this result. Can somebody help me? * _(The formula originally had $x_t+u_t$, i.e., the plus was changed to multiply)_


title: Variance-covariance matrix of the errors in linear regression

text: How is the var/cov error matrix calculated by statistical analysis packages in practice? This idea is clear to me in theory. But not in practice. I mean, if I have a vector of random variables $\textbf{X}=(X_{1}, X_{2}, \ldots, X_{n})^\top$, I understand that the variance/covariance matrix $\Sigma$ will be given the external product of the deviance-from-the-mean vectors: $\Sigma=\mathrm{E}\left[(\textbf{X}-\mathrm{E}(\textbf{X}))(\textbf{X}-\mathrm{E}(\textbf{X}))^\top\right]$. But when I have a sample, the errors of my observations are not random variables. Or better, they are, but only if I take a number of identical samples from the same population. Otherwise, they're given. **So, again my question is: how can a statistical package produce a var/cov matrix starting from a list of observations (i.e. a sample) supplied by the researcher?**


title: What is the opposite to "adaptive"?

text: In statistics, there are methods that are named with "adaptive", such as "adaptive randomization" in experimental design, and "adaptive learning". What is the word that is opposite to "adaptive", and used in the opposite methods? Thanks!


title: Stata: Test over time --> Which method?

text: I am working on a project regarding the influence of temperature and time variables (weekday, weekend, month, season etc.) on the sales of a local bakery. The research goal is to be able to better predict bakery sales (in my city) and to save food in this way. What data do I have? I do have the data of all variables for more than 3 years (>1000 days). The sales data is normalized (ratio scaled) and the independent variables have different scales of measurement with nominal as the lowest one (day of the week etc.). What have I done so far? So far I prepared the data and ran a few multiple linear regression tests e.g. regress sales temperature weekend_TrueFalse But I just saw in another feed that this may not be the right method whereas time series methods are. What would you say which method I should use (and if possible could you give me the stata commands)?


title: Finding a narrower confidence interval for a given CI, sample mean and size

text: I'm trying to understand confidence intervals but having some trouble. I've been doing some exercises I found online and I'm stuck on this question: I have been given a 95% conﬁdence interval for a population proportion: (0.35, 0.40), a sample size of 200, and I need to find a 99% confidence interval. The methods I would go to first involve using standard deviation, which I don't have. How can I approach this question without knowing variance? The whole quiz is about normal distributions, to give it some context. This isn't really homework but I'm tagging it this way because I'm looking for a similar outcome. The quiz has answers at the end; what I want is to be able to solve it myself.


title: Mediation model with linear regression

text: In my master thesis I have drawn a few hypotheses. I have answered them all with linear regression. In these linear regressions, I took control variables into account. My question is: do I have to run a mediation analysis? Or is it also possible to report the regressions of all relations separately (for example: X -> Y, X -> M, M -> Y and A -> Y)? Here is how my model look like: ![enter image description here](http://i.stack.imgur.com/gxHfX.jpg) My main hypothesis is about the relation between X and Y. I hope my question is clear. Thank you in advance!


title: How to test whether the majority of vertices belong to one giant component?

text: I have a sample of 200 independent networks and I want to test the hypothesis that majority of vertices belong to one giant component. I wonder what is the appropriate approach to do that. More formally, suppose we have two variables, $A$ and $B$, where $A$ denotes number of vertices belonging to a giant component in a network $i$ and $B$ denotes number of vertices not belonging to giant component in network $i$. The task is to test the hypothesis $\mu_A > \mu_B$. Because $A$ and $B$ are clearly dependent here, I doubt that classical statistical tests (e.g., Z-test) are appropriate here. Any suggestion would be greatly appreciated.


title: Exponential upper bound for $\sum_{k=1}^{N-1} \alpha^k \beta^{\frac{1}{N-k}}$

text: I will be too happy help me find the exact value or a very tight exponential upper bound for: $$\sum_{k=1}^{N-1} \alpha^k \beta^{\frac{1}{N-k}}$$ where $0 \leq \alpha < 1$, $\beta= \exp(-N^2\zeta)$ and $\zeta$ is a positive small value. Thanks a lot in advance.


title: Interpretation of Residual Covariance

text: I have recently begun using the Structural Equations Modeling (SEM) Method of confirmatory factor analysis for a research endeavor in educational science. My question is, suppose I have two latent variables each with many of its own manifest variables. If I observe a covariance between the two latent variables (I'm not sure if specifically what program I used to reach this point matters but if so I use the lavaan package in R), what does it mean? I believe it means that it is a residual covariance indicating the presence of a common factor not shown by their predictors but I am unsure. If this is the case what would the level of residual covariance mean? i.e. is there a cut-off point where it could be considered statistically insignificant? Thank you!


title: Standardizing standardized data?

text: Let's assume I standardize my data in SPSS through converting them to z-values. Then I save my dataset and load it into AMOS to do some regressions. Now, the data is already standardized, but AMOS will give me two sets of regression coefficients, labeled as standardized and unstandardized. Does this mean that the unstandardized values in the output already are standardized? What is the meaning of standardized regression coefficients calculated based on already standardized data?


title: Randomized trace technique

text: I've met the following _randomized trace technique_ in M. Seeger, “Low rank updates for the Cholesky decomposition,” University of California at Berkeley, Tech. Rep, 2007. $$\operatorname{tr}(\mathbf{A}) = {E[\mathbf{x}^T \mathbf{A} \mathbf{x}]}$$ where $\mathbf{x} \sim N(\mathbf{0},\mathbf{I})$. As a person without deep mathematics background, I wonder how this equality can be achieved. Moreover, how can we interpret $\mathbf{x}^T \mathbf{A} \mathbf{x}$, for example geometrically? Where should I look in order to understand the meaning of taking the inner product of a vector and its range value? Why is the mean equal to the sum of the eigenvalues? Besides theoretical property, what is its practical importance? I've written a MATLAB code snippet to see whether it works               #% tr(A) == E[x'Ax], x ~ N(0,I)          N = 100000;     n = 3;     x = randn([n N]); % samples     A = magic(n); % any n by n matrix A          y = zeros(1, N);     for i = 1:N         y(i) = x(:,i)' * A * x(:,i);     end     mean(y)     trace(A)      The trace is 15 where the approximation is 14.9696.


title: Seasonality and work days in time series

text: I have a time series which is daily data for workdays only. There is an almost obvious seasonality comparing the plot of the raw data that appears at the end of each month. However, there are not the same number of workdays for every month. Some times a month has 20 or 21 workdays. Do you have any idea how to treat the data? Do you suggest to manipulate the data in order to create a series with a constant frequency? The problem is that in R I have to set the frequency before running the ARIMA.


title: How to test whether a sample of data fits the family of Gamma distribution?

text: I have a sample of data which was generated from a continuous random variable X. And from the histogram I draw using R, I guess that maybe the distribution of X obeys a certain Gamma distribution. But I do not know the exact parameters of this Gamma distribution. My question is how to test whether the distribution of X belongs to a family of Gamma distribution? There exists some goodness of fit tests such as Kolmogorov-Smirnov test, Anderson-Darling test, and so on, but one of the restriction when using these tests is that the parameters of the theoretical distribution should be known in advance. Would anyone please tell me how to solve this problem?


title: The (statistically unsound) exploratory study (no particular hypothesis)

text: It seems like a proposed study gets penalised if it doesn't also predict the direction of an effect. For instance, if I want to see the effect of a certain medication on happiness ratings, but do not predict whether I think they'll be decreasd or increased by that manipulation. Why is it considered unsound to run an exploratory study, i.e. one that doesn’t have a specific prediction (hypothesis) but simply wants to see what effect a certain manipulation has on another dependent variable (i.e. you have a hunch there's an effect there but don’t quite know in which direction it'll be)? In other words, what exactly is the formal statistical principle that is violated if you run a study without formulating an apriori hypothesis? It seems that if you're happy to pay the 'penalties', e.g. use 2-tailed instead of 1-tailed tests and correct for all the multitude of tests, any exploratory result should be just as believable as those of a hypothesis-driven study!


title: Extracting meaningful information from a logistic regression model

text: Possible warning: basic question ahead. Let's say that I model whether I wear red shoes depending on the weather. Red shoes, which is my dependent variable, is a dichotomous variable as I either wear them or don't. Weather is a variable with five 'levels' and I'm trying to find the probability that I will wear read shoes. Let's say I model out this relationship with a logistic regression model in R:               mod = glm(shoes ~ weather, data=mydat, family=binomial(link="logit"))      Now, what I am interested in is finding what are the probabilities for wearing red shoes for each of the five 'levels' in weather. So perhaps I might have a 20% probability of wearing red shoes when cold, 30% when mild, 40% when warm, and so forth. I'm wondering if modeling is a requirement for finding this information? If so, how does one go from somewhat meaningful regression coefficients to meaningful probabilities in R?


title: What is the probability that the $k$th element falls in a specific interval?

text: The question I'm referring to comes from Stack Overflow: http://stackoverflow.com/questions/8723652/estimating-number-of-results-in- google-app-engine-query _In short:_ **With $N$ ordered samples of a uniform distribution, how to better estimated $N$ with $k$th sample info?**   1. How to compute the probability that the $k$th fall in a specify interval? For example, Is there a c.d.f for the order statistics?   2. Is there a better way to estimate the $N$? For example, use three samples 1000th, 2000th, 3000th together? **[EDIT1]** @whuber: I have read the link you provide, but I still have some question about how to use it correctly. > The probability U(k) falling in the interval [u,u+du] is equaled to ![enter > image description here](http://i.stack.imgur.com/l8Oez.png) (from http://en.wikipedia.org/wiki/Order_statistic#The_order_statistics_of_the_uniform_distribution)   1. Is this function a p.d.f or a c.d.f? If it is a p.d.f, could I integrate it to get a c.d.f?   2. While N is big (for example 60000), is there any approximation I can use?   3. Try to explain this issue more clearly: > We did N(unknown number) random samples for an uniform distribution [0,1]. > These samples are collected and sorted (but we still don't know the exactly > total number of them). How to guess the N based on the known kth value? For > example, if the 10th sample is 0.1, we may guess the N is 10. If the 10th > samples is 0.01, we may guess N is 100. But how accurate this guess will be? > If I can have more info such as the 10th sample is 0.1 and 20th sample is > 0.2, will that help to get better results?


title: What does "$R^2$ estimates the combined dispersion against the single dispersion of the observed and predicted series" mean?

text: In a paper I came across the description of $R^2$ as "it estimates the combined dispersion against the single dispersion of the observed and predicted series". I am not able to understand this statement. I understand that $R^2$ is the square of the correlation coefficient of predicted and observed values. Please help! Thanks in advance. :)


title: Finding matrix eigenvectors using QR decomposition

text: **First** , a general linear algebra question: Can a matrix have more than one set of (unit size) eigenvectors? From a different angle: Is it possible that different decomposition methods/algorithms (QR, NIPALS, SVD, Householder etc.) give different sets of eigenvectors for the _same_ matrix? **Second** , regarding QR decomposition: Are the columns of the Q matrix the eigenvectors? How can their eigenvalues be easily found (post the QR decomposition)?


title: Conversion rates over time

text: I'm looking at some data conversion rates for an ad over time.                         view      clicks         Day 1      100       10            Day 2      150       13     Day 3      90        9     Day 4      130       20     Day 5      150       21      Given that there is quite a bit of variation in the conversion rate of the ad from day to day, I want to know the confidence interval for the conversion rate across the 5 days. For each day I could develop confidence intervals for the conversion rate, and then compare across the days. However, I'm just not sure how to know with 90% confidence that the confidence interval for all the days was between so and so. Can anyone help! In R, I can calculate the confidence intervals for two seperate ads using the following function:               abtestfunc <- function(ad1, ad2){           sterror1 = sqrt( ad1[1] * (1-ad1[1]) / ad1[2] )           sterror2 = sqrt( ad2[1] * (1-ad2[1]) / ad2[2] )           minmax1 = c((ad1[1] - 1.96*sterror1) * 100, (ad1[1] + 1.96*sterror1) * 100)           minmax2 = c((ad2[1] - 1.96*sterror2) * 100, (ad2[1] + 1.96*sterror2) * 100)           print( round(minmax1,2) )           print( round(minmax2,2) )     }      Or is it right to calculate the means conversion rate for the conversion rates and then develop conf intervals. Then compare whether the majority of days fall within the conversion rate.


title: Generating data with a pre-specified odds ratio

text: I have to simulate (in javascript) the likelihood that members of a certain population would have a disease, given a contextual value. All information I have comes from papers that show odds ratios with $95\%$ confidence interval. I found Box-Muller method and other approaches to normal distribution, but I cannot realize the way to generate random numbers with OR, CI percentage and the interval.


title: Simple "yes / no" question. Uncertain how to analyze the results

text: I am in the midst of conducting a study which compares content in a printed catalog of records to content in a much larger online catalog of similar records. I have randomly selected 50 entries from the printed catalog and searched the online catalog to see if there are exact matches for these records. The results are 44 matches ("Yesses"), 6 non-matches ("Nos"). My hypothesis is that a significant number of records from the printed catalog will not be included in the online catalog. I have very, very little experience with statistical analysis and am unsure how to proceed. I appreciate any help or suggestions anybody can give me.


title: How to run use `lmer` in R to assess effect of condition, task order, and type of rating?

text: I would like to estimate a multivariate model with the following properties:   * Fixed effects:      1. Intercept      2. Main effect of condition     3. Main effect of Task order     4. Task order X Condition interaction    * Random effects:      1. Main effect of TypeOfRating      2. TypeOfRating X Condition interaction     3. TypeOfRating X Task order interaction     4. TypeOfRating X Task order X Condition interaction  Can this be done using `lmer`?


title: Support vector machines and regression

text: There's already been an excellent discussion on how support vector machines handle classification, but I'm very confused about how support vector machines generalize to regression. Anyone care to enlighten me?


title: R not testing certain variables in forward stepwise regression?

text: `FullModel<\- (lm(Fubar~.-Foo-Bar,data=BarFoo)) NullModel<-(lm(Fubar~1)) step(NullModel,scope=formula(FullModel),direction="forward",k=log(nrow(BarFoo)))` When doing the above forward stepwise regression, the forward steps halt before certain variables are added in. does this mean that they do not improve the AIC score or does it mean that the inputs are in error?


title: "Strength" of hypothesis tests

text: I was referred to this site from mathoverflow. After some brief searching, I haven't seen another thread on my particular question, but if it's already been answered elsewhere then I hope you'll forgive the spam. Has anyone ever "ranked" the various commonly-used parametric and non- parametric hypothesis tests in terms of the "strength" of the conclusions that one can ascertain from them? That is to say, if I have two populations that are wildly different from each other (and about which I have absolutely no assumptions whatsoever on the prior distributions), and I'd like to make the strongest possible case that they are different, what would be the right hypothesis test to use? An example would be if I have one population of samples that all lie between $0$ and $1$ and another population of samples that all lie between $10$ and $11$. Running a $t$-test would likely show that my sample means are different, but I would think that I could make stronger conclusions than that (based on the fact that largest element in one set is smaller than the smallest set in the other, for example). (as an aside, in the particular case I'm dealing with, one population has lots of samples and the other is quite small, and both have rather small standard deviations)


title: How to determine quantiles (isolines?) in a multivariate normal distribution

text: ![enter image description here](http://i.stack.imgur.com/ZoPHf.png) I'm interested in how one can calculate a quantile on a multivariate distribution. In the above picture, I have drawn the 5% and 95% quantiles on a given univariate normal distribution (left). For the right multivariate normal distribution, I am imagining that an analog would be an isoline that encircles the base of the density function. Below is an example of my attempt to calculate this using the package `mvtnorm` - but to no success. I suppose this could be done by calculating a contour on the results of the multivariate density function, but i was wondering if there is another alternative (e.g. analog of `qnorm`). Thanks for your help. ### Example:               mu <- 5     sigma <- 2      vals <- seq(-2,12,,100)     ds <- dnorm(vals, mean=mu, sd=sigma)          plot(vals, ds, t="l")     qs <- qnorm(c(0.05, 0.95), mean=mu, sd=sigma)     abline(v=qs, col=2, lty=2)               #install.packages("mvtnorm")     require(mvtnorm)     n <- 2     mmu <- rep(mu, n)     msigma <- rep(sigma, n)     mcov <- diag(msigma^2)     mvals <- expand.grid(seq(-2,12,,100), seq(-2,12,,100))     mvds <- dmvnorm(x=mvals, mean=mmu, sigma=mcov)          persp(matrix(mvds,100,100), axes=FALSE)     mvqs <- qmvnorm(0.95, mean=mmu, sigma=mcov, tail = "both") #?          #ex. plot        png("tmp.png", width=8, height=4, units="in", res=400)     par(mfcol=c(1,2))          #univariate     plot(vals, ds, t="l")     qs <- qnorm(c(0.05, 0.95), mean=mu, sd=sigma)     abline(v=qs, col=2, lty=2)          #multivariate     pmat <- persp(seq(-2,12,,100), seq(-2,12,,100), matrix(mvds,100,100), axes=FALSE, shade=TRUE, lty=0)     cont <- contourLines(seq(-2,12,,100), seq(-2,12,,100), matrix(mvds,100,100), levels=0.05^2)     lines(trans3d(cont[[1]]$x, cont[[1]]$y, cont[[1]]$level, pmat), col=2, lty=2)          dev.off()


title: Variance of log-odds probabilities

text: I am wondering how I could express the variance of log-odds into understandable terms. For example the variance in the log-odds of crime being reported to the police between neighbourhoods is 0.07 (0.01), and the mean log-odds is -0.65 (0.03). Just transforming this variance to probability given me 0.51. So if the probabilities have a variation of 0.51, the standard deviation is then 0.7, which would make the 95% coverage bound wider than 1, which is impossible for probabilities. Another idea was to calculate the 95% coverage bounds of my log-odds: $(-0.65) +/- [1.96\sqrt{0.07}] = [-1.19, -0.15]$ If I transform these to probabilities: $\exp(x)/[1+\exp(x)] = [0.23,0.46]$ So the chance of a crime being reported to the police can vary from 23% to 46%. So the 95% range is actually only 0.23 or 23%. Could someone explain me why this first calculation was wrong, and whether the second is correct? Would you have any other suggestions on how to make this variance more graspable?


title: Multilevel calibration curve in jags

text: # The problem I am trying to fit a calibration curve from known x and y with jags. I would then like to predict new xes, based on the measurement of new ys. For this, I am following Hamada et al. (2003), J Quality Technology 35:194-205. Since I have subject effect, I wanted to incorporate that in my calibration curve, but prediction has to be done for unknown subjects. I tried to follow Gelman et al. (2003): Bayesian Data Analysis, p.137, to expand the jags model accordingly. I am new to jags, and I wanted to know whether I am doing this properly. Any pointers to improve my analyses and predictions are very welcome! # R code with simulated data               library("rjags")     library("lme4")          # create data set     sim.data <- data.frame(x = c(rep(100, 10), rep(200, 10), rep(300, 10), rep(400, 10)))     sim.data$indid <- c(rep(LETTERS[1:10], 4))         sim.data$ranintr <- c(rep(rnorm(10, 0, 10), 4))     sim.data$resid <- rnorm(40, 0, 30)          a <- 200     b <- 1          sim.data$y <- a + b*sim.data$x + sim.data$ranintr + sim.data$resid          # prepare simulated data for jags     N1 <- length(data.sim$y)         x <- sim.data$x     y <- sim.data$y         indid <- factor(sim.data$indid)     K <- length(unique(indid))          # new data for prediction     N2 <- 2     x2 <- c(NA, NA)     y2 <- c(250, 350)          jags <- jags.model('simModelCalibLinearHierarchical_predict.bug',                        data = list('x' = x, 'y' = y,                            'N1' = N1, 'K' = K, 'indid' = indid,                            'N2' = N2, 'x2' = x2, 'y2' = y2),                        n.chains = 2,                        n.adapt = 100)          update(jags, 300000)          samples <- coda.samples(jags, c('alpha0', 'beta',                                      'sigma2', 'sigma2.a.indid', 'ICC', 'x2'),                              100000, thin = 100)          fit.1 <- lmer(sim.data$y ~ sim.data$x + (sim.data$x | sim.data$indid))          summary(samples[,c('alpha0', 'beta', 'x2[1]', 'x2[2]')])     summary(fit.1)      # BUGS code               model{            # calibration       for(i in 1:N1){       # the calibration observations         y[i] ~ dnorm(mu[i], tau)         mu[i] <- alpha0 + alpha[indid[i]] + beta * x[i]       }       # K should be the number of unique subjects       for(i in 1:K){         alpha[i] ~ dnorm(0, tau.a.indid)       }            # prediction       for(i in 1:N2){       # the number of predictions         y2[i] ~ dnorm(mu2[i], tau)         mu2[i] <- alpha0 + alpha2[i] + beta * x2[i]         x2[i] ~ dnorm(mux2, taux2)         alpha2[i] ~ dnorm(mean.alpha, tau.a.indid)       }            # priors for calibration       tau ~ dgamma(0.01, 0.01)       alpha0 ~ dnorm(0, 0.0001)       beta ~ dnorm(0, 0.0001)        tau.a.indid ~ dgamma(0.01, 0.01)            # priors for predictions       mux2 ~ dnorm(300, 0.0001)  # prior chosen so that (0, 600) contains 99.97%                                   # of its distribution; this reflects our knowledge                                   # that values must be positive        taux2 ~ dgamma(0.01, 0.01)            # calculations       mean.alpha <- mean(alpha)       sigma2 <- 1 / tau       sigma2.a.indid <- 1/tau.a.indid       ICC <- sigma2.a.indid / (sigma2 + sigma2.a.indid)     }


title: How to interpret logarithmically transformed coefficients in linear regression?

text: I have 1 continuous dependent and 1 continuous predictor variable that I've logarithmically transformed to normalise their residuals for simple linear regression. Can someone give me some pointers on how to relate these transformed variables to their original context? My situation: I want to use linear regression to predict the number of days of school that pupils missed in 2011 based on the number of days they missed in 2010. Because most pupils miss 0 days or just a few days the data is positively skewed to the left hence the need for transformation to use linear regression. I've used log10(var+1) for both variables (I used +1 for pupils who had missed 0 days school). I'm looking to use regression because I want to add in categorical factors - gender/ethnicity etc too. The audience I want to feed back to wouldn't understand log10(y) = log(constant) + log(var2)x (and frankly neither do I). Are there better ways of interpreting transformed variables in regression? I.e. for ever 1 day missed in 2010 they will miss 2 days in 2011 as opposed to for ever 1 log unit change in 2010 there will be x log units change in 2011? I'm using SPSS v20. * * * Actually I suppose my specific question is this. Looking at the following website http://www.ats.ucla.edu/stat/sas/output/sas_negbin_output.htm it quotes that: "This is the negative binomial regression estimate for a one unit increase in math standardized test score, given the other variables are held constant in the model. If a student were to increase her mathnce test score by one point, the difference in the logs of expected counts would be expected to decrease by 0.0016 unit, while holding the other variables in the model constant." Is this saying that for every one unit increase in the score of the UNTRANSFORMED variable math leads to a 0.0016 decrease from the constant (a)? So, if UNTRANSFORMED maths score goes up by two points, I subtract 0.0016*2 from the constant a? Therefore, to model the effect, I get the geometric means by using exponential(a)) and exponential(a+beta*2) and calculate the percentage difference between these two to say what effect the predictor variable(s) has/have on the dependent variable? Or have I got that totally wrong? Sorry for framing this in a long question.


title: Relationship between the characteristics of training data set and built decision tree

text: I have a training data set for a binary classification problem. There exist two possible scenarios, one is that all of the training data set are labeled as positive; another one is that the training data set includes labeled positive ones and labeled negative ones. Assume that I use this training data set to train a decision tree. How do these two different scenarios affect the trained tree? Moreover, if the ratio of positive and negative ones changes, how does this change affect the built tree model?


title: Checking the distribution of classes in a sample of a dataset

text: I am new here, so I might have missed a similar question. I am trying to do things good in my research and using proper statistical approaches, but in computer science, we had quite a weak training on statistical methods applied to our research. So, what I am asking is probably trivial. Here is the problem: We have a dataset of del.icio.us (the social bookmarking website) and want to work on a subset of this dataset for a manual study of the usage of tags. In this dataset, we have urls, users and tags that are linked by triplets defining an association of a tag by a user to a url. We have 401 970 328 such triplets in our dataset, so we can't look at all of them and want to choose a sample of it. Naively, at the beginning, we started with a purely random sampling to choose 500 user-bookmark pairs and all their associated tags. Because there are some more popular tags on del.icio.us, our sample has a long-tail distribution of tag usage, similar to the one from the whole dataset. What I would like to check is that the actual distribution of each tag is similar in the sample as in the original dataset. That is, if "java" appears 10 times more often than "indonesia" in the dataset, it should be a similar distribution in the sample. How would I go at doing this? one of the issue is that in the sample, some classes (tag) present in the full dataset will never appear. The second question is also: how to decide of a good size for the sample?


title: How to determine optimal number of clusters?

text: For a multilabel dataset, i would like to find the number of clusters involved in it. The below example gives more details about the problem:               Label_A: feature values     Label_B: feature values     Label_A, Label_C: feature values      Label_C: feature values ... etc      We have say $n$ datarecord. Label field may have single label/multilabel(as in the case of record 3). I would like to determine the number of cluster involved in the data. Assuming number of label as the number of cluster results in bad accuracy. This is because there may be case where single label can have multiple cluster. In this case, if we can find more cluster and assign two or more cluster to same label, we can increase the accuracy. Hence, how do you find the number of cluster present in the multilabel data?


title: Determining whether a group of spatially extensive data is centered around one point in space

text: Consider a single x axis representing points along a line in space. I've got a set of data (in this case the receptive fields of hippocampal place cells recorded as a rat runs along a linear track), which are themselves spatially extensive, ie. if displayed graphically against the axis they would look like a series of horizontal bars, distributed at different points along the axis. * * *                 ____         ____                ____ fields      _**___ ___ __ ___ __ ___ __ ___ __ ___ ____**_axis What I want to do is to define a point of interest on the line, and determine whether there are significantly more fields overlapping this point than elsewhere. I've achieved a simpler version of this by finding the centre of the fields and then measuring distances from this point and another arbitrarily defined point, but this approach loses all of the information that comes from the fact that they are spatially extensive. Can anyone point me in the right general direction? Even knowing the broad class of statistical test I should use for this kind of situation would be a big help! Thanks!


title: What does plotting residuals from one regression against the residuals from another regression give us?

text: I am working with the dataset of some heights and weights at different ages. My professor wants me to plot the residuals from regression of `soma.WT9` against the residuals from regression of `HT9.WT9` (don't mind the notation it's just two columns where soma is being regressed on `WT9` and `HT9` being regressed on `WT9`). What is the purpose of this plot.


title: Accounting for biases in the data (normalizing? detrending?)

text: What's it called when I'm trying to remove the effect of a variable on another variable? Am I in the right ballpark using the terms "detrending" or "normalizing" or "controlling for Y on X"? What's the general process for doing this? The problem I'm working with currently is as follows: I work at a mine that loads trains that are 160 wagons in length. After the trains are filled, they leave to get unloaded and then later return to the mine. Sometimes the wagons end up in different positions in the train. Our aim is to maximise tonnes in every wagon. We've noticed two things   1. some wagons IDs load consistently low in tonnes   2. wagons in poisitions near the start or end of the train load consistently low. i.e. tonnes go down as distance from middle position goes up I've got a linear model fitted in R like                   > lm(df$tonnes ~ df$dist_fr_middle)         Coefficients:             (Intercept)    df$dist_fr_middle               113.92001     -0.03915        How do I go about creating a new column in my dataframe for distance-from- middle-adjusted tonnes? Thanks for your help!


title: How different random number generators can be more similar than identical ones?

text: I have 6 random number generators. They are "black boxes", i.e. I do not know if they are the same or different. For example I do not know if they provide the same arithmetic averages and/or root mean square deviations. My goal is to check if they are the same or different. The problem is that I can generate a limited set of numbers (every random number generator gives me only 20 numbers). The procedure that I use is as following.   1. I choose a pair of generators.   2. I use each of them to generate 20 numbers.   3. For every generator I use these numbers to calculate the mean.   4. Then I calculate the difference between the means (mean1 - mean2). I will call this difference as "reference difference". Now I want to know if this difference is "real". In other words I want to know if the observed difference in means is caused by the fact that I do have two different random number generators or I have identical generators and the mean is different just because of chance (small number of numbers). To answer the last question I use the following procedure.   1. I combine the two sets of numbers (20 numbers + 20 numbers).   2. Then I randomly split this set into two subsets of equal size.   3. For each subset I calculate means and, as before, I calculate the difference between the two means (I will call this difference as "test difference").   4. I compare this test difference with the reference difference.   5. I repeat this procedure many times and to see in home many cases the test difference is smaller then the reference difference. I thought that if the reference difference is larger than the test difference in most of the cases than it is conditioned not by chance but by the fact that the compared random number generators are different. In contrast, if the generators are identical than the test difference has 50% of chances to be larger (or smaller) then the reference one. However, I see that for the most of the considered pairs of the generators the reference difference between the averages is smaller than the test differences. How could it be? **It looks like the the generators are more similar than identical generators?**


title: Standardizing unbalanced repeated measures data

text: I have a dataset of height and weight measurements for adult patients in hospital clinic and want to calculate a standardized measure of boy mass index (BMI). For a given pair of measurements, it is easy to calculate BMI: $BMI = weight (kg) / \left( height (m) \right)^2$ The problem I have is that some patients will have multiple visits to the clinic and therefore will have multiple height and weight measurements. So the number of BMI measures per patient can range from 1 to $i$. Example dataset:               Patient | BMI     --------|------        0001 | 24.0        0001 | 24.7        0001 | 25.3        0002 | 21.2        0003 | 29.9        0003 | 29.8        0004 | 25.5      My goal is to standardize the BMI measures, taking into consideration that BMI measures dependent on the patient and that there are different numbers of measures for each patient. Does anyone have experience with this type of problem?


title: Indirect measurement of the distributions

text: Suppose that there is a field where cabbages grow. We have weighed all the cabbages in spring and got Gaussian distribution of weights. The same we did in autumn and got another Gaussian distribution. The question is how to get distribution of growth rates (times over the summer) if cabbages are not numbered? Will it be Gaussian or not?


title: Not getting correct lag values in output

text: I am using R to calculate lagged correlation. Instead of getting lag values -20 to 20, I am getting the following values:                        [,1]      [1,] -0.83333333      [2,] -0.75000000      [3,] -0.66666667      [4,] -0.58333333      [5,] -0.50000000      [6,] -0.41666667      [7,] -0.33333333      [8,] -0.25000000      [9,] -0.16666667     [10,] -0.08333333     [11,]  0.00000000     [12,]  0.08333333     [13,]  0.16666667     [14,]  0.25000000     [15,]  0.33333333     [16,]  0.41666667     [17,]  0.50000000     [18,]  0.58333333     [19,]  0.66666667     [20,]  0.75000000     [21,]  0.83333333      I am using the following R code:               indiv_product_data_ts<-ts(indiv_product_data, frequency=12)     x<-indiv_product_data_ts[,2]     y<-indiv_product_data_ts[,3]     x_model<-auto.arima(x)     correlation<-prewhiten(x,y,x.model=x_model,plot=FALSE,lag.max=define_lag)     lag=data.frame(correlation$ccf$lag)     ccf=data.frame(correlation$ccf$acf)     len<-nrow(ccf)     product_list<-rep(product_list[prod_ID],len)     Output<-cbind(product_list,lag,ccf)     colnames(Output)<-c("BaseProduct","lag","ccf")      What could be the reason? I would appreciate any kind of help.


title: Evaluating Proportional Hazards with time-dependent variables

text: I'm doing survival analysis and using the survival package from R to do it. I have 7 variables, 5 of which are time-dependent (all except sexo and peso). Here is my model:               Call:     coxph(formula = Surv(tempo1, tempo2, status) ~ sexo + clima +          CODE_06 + TP_media7 + ndvi + peso + epoca, data = newftable)                                       coef exp(coef) se(coef)       z      p     sexoM                -0.1098     0.896   0.5472 -0.2006 0.8400     clima                -0.2368     0.789   0.0731 -3.2408 0.0012     CODE_06Regadio       -0.9961     0.369   0.6532 -1.5249 0.1300     CODE_06Sequeiro      -0.0711     0.931   0.5009 -0.1420 0.8900     TP_media7(12,22]     -0.5616     0.570   0.8064 -0.6964 0.4900     TP_media7(22,32]     -1.0945     0.335   0.8680 -1.2609 0.2100     TP_media7(32,41]      0.0533     1.055   0.9557  0.0558 0.9600     ndvi(3e+03,6e+03]    -0.1502     0.861   0.8386 -0.1791 0.8600     ndvi(6e+03,8.05e+03] -0.7635     0.466   1.0205 -0.7482 0.4500     peso[845,1005]       -0.3257     0.722   0.4961 -0.6565 0.5100     epocamid_inv_rep      0.4414     1.555   0.7077  0.6237 0.5300     epocamid_pos_inv     -0.0722     0.930   0.7601 -0.0950 0.9200     epocamid_rep_pos      0.4424     1.556   0.7770  0.5693 0.5700     epocapos_repr         0.3200     1.377   0.8114  0.3943 0.6900     epocarepr             0.1755     1.192   0.7742  0.2267 0.8200          Likelihood ratio test=28.4  on 15 df, p=0.019  n= 2267, number of events= 30         (94 observations deleted due to missingness)      I want to test the PH assumption (using cox.zph) only for sexo and peso, since time-dependent variables fail to satisfy this assumption (Kleinbaum & Klein, 2005). I did:               >cox.zph(teste)                           rho    chisq     p     sexoM                -0.14944 7.14e-01 0.398     clima                 0.19168 9.60e-01 0.327     CODE_06Regadio        0.06300 1.13e-01 0.737     CODE_06Sequeiro       0.24125 2.01e+00 0.157     TP_media7(12,22]      0.21771 1.29e+00 0.257     TP_media7(22,32]      0.26118 1.84e+00 0.175     TP_media7(32,41]      0.20317 9.91e-01 0.319     ndvi(3e+03,6e+03]     0.04029 5.40e-02 0.816     ndvi(6e+03,8.05e+03] -0.04816 7.40e-02 0.786     peso[845,1005]        0.08051 2.82e-01 0.595     epocamid_inv_rep      0.00243 2.26e-04 0.988     epocamid_pos_inv     -0.17620 8.17e-01 0.366     epocamid_rep_pos      0.10899 2.69e-01 0.604     epocapos_repr         0.27081 1.95e+00 0.162     epocarepr             0.09812 2.45e-01 0.621     GLOBAL                     NA 1.12e+01 0.737      And it seems that every variable satisfies the assumption (Although I know it's not true). When I have a model only with the 2 time-independent variables, the result is the following:                                  rho  chisq     p     sexoM          -0.1797 1.6287 0.202     peso[845,1005]  0.0181 0.0166 0.897     GLOBAL              NA 2.8486 0.241      My question is, what is the right procedure?


title: Signature plots definition

text: While reading some articles about High frequency statistics in financial markets, more precisely wile reading about microstructure noise, I was confronted by some signature plots simulated graphics. It seems me to be a basic concept in statistics or, at least, just a well known term. However for a lack of culture on general statistics, I don't understand what it means (so neither what it illustrates) and I could not find this definition anywhere on the web.


title: How to interpret negative factor loadings and residual variances in SEM?

text: I've been trying to fit a SEM model, using misconceptions and the correct conception as latent variables. My instrument is a set of items that ask a binary yes or no question. (Is something in a particular category?) The correct conception would have all items that would be correct if yes was chosen loaded onto it and the misconception would have all items we would expected to be answered yes if the subject had the given misconception. In this case, there would be some overlap between the misconception and correct conception factor loadings. i.e. cross- loading. Is cross-loading in SEM acceptable? and if so, how would I interpret a negative cross loading and negative residual variance for the items?


title: How does predictive model for the Eurovision Song Contest work?

text: I've encountered interesting prediction of Eurovision Song Contest http://mewo2.com/nerdery/2013/05/12/eurovision-2013-first-predictions/ it based on some kind of Bayesian model I assume but I don't understand how it works (the description is too sketchy for me). Along the lines it says that it uses Markov Chain Monte Carlo sampler but what exactly does it sample?


title: Bootstrapping power estimates for a bootstrap test

text: Assume I want to use a (nonparametric) bootstrap test for a hypothesis with a sample size of $n_1$ and I already have $n_0$ actual samples on which to base my power estimates. Usually, we would also have $n_1 > n_0$. Is it a valid procedure to estimate the power of my test by using nested bootstrapping? Basically, I would repeatedly sample $n_1$ samples with replacement from my $n_0$ samples and apply my non-parametric bootstrap test to each of these samples. Finally, I would look at the percentage of the bootstrap tests which was significant at my $\alpha$ level. Are there pitfalls I have to look out for? I already did some googling without much success, perhaps because I do not know the proper search terms. Therefore it would be also nice to know how the procedure (if it exists) is called, so I can find references.


title: Subject specific vs population average predictions

text: I am in doubt whether in my thesis I should report on the subject specific predictions of the probability to respond with an 'I don't know' answer, or the population average. Consider for example the following graphs (I have not gotten them smaller yet): ![enter image description here](http://i.stack.imgur.com/Df5cw.png) ![enter image description here](http://i.stack.imgur.com/UjCpH.png) They both have the same conclusion that the probability doubles for females. However, in terms of frequencies this means on average in the population I will have 3.3% missingness for females and 2.5 for males. For the average female individual though (in the median clusters= interviewed by the median interviewer, in the median country), I will have 1.9% missingness and 1.5% for this average male individual. I think the population average is more informative for my cause: evaluate whether missingness is related to gender and whether the rates are problematic. However, reading the following passage (Allison, P D (2009) Fixed Effects Regression Models, Quantitative Applications in the Social Sciences Series, Sage, Thousand Oaks, California): > “if you are a doctor and you want an estimate of how much a statin drug will > lower your patient’s odds of getting a heart attack, the subject-specific > coefficient is the clear choice. On the other hand, if you are a state > health official and you want to know how the number of people who die of > heart attacks would change if everyone in the at-risk population took the > stain drug, you would probably want to use the population–averaged > coefficients. Even in the public health case, it could be argued that the > subject-specific coefficient is more fundamental" _I can't seem to translate this to my specific situation. Anyone have any advice on my situation?_


title: Can you compare AIC values as long as the models are based on the same dataset?

text: I am doing some forecasting in R using Rob Hyndman's _forecast_ package. The paper belonging to the package can be found here. In the paper, after explaining the automatic forecasting algorithms, the authors implement the algorithms on the same data set. However, after estimating both an exponential smoothing and ARIMA model they make a statement I do no understand (on page 17): > _Note that the information criteria are not comparable._ I thought that **an advantage of using AIC for model selection is that we can compare AIC values from different models, as long as they are estimated using the same data set. Is this incorrect?** This matter is of particular interest to me, as I was planning on combining forecasts from different model classes (e.g. exponential smoothing and ARIMA) using so called Akaike weights (see Burnham and Anderson, 2002, for discussion on Akaike weights) ### References   * Burnham, K. P., & Anderson, D. R. (2002). Model selection and multi-model inference: a practical information-theoretic approach. Springer Verlag.


title: PCA on binary data (0's & 1's) -> what does it mean when a PC is correlated with the nr of 1's per subject?

text: I conducted a PCA on dichotomous variables (0's and 1's). The dataset consists of human subjects and a few thousand genetic variants, where the presence of a genetic variant is indicated with 0's and 1's. My first PC correlates >.9 with the nr of 1's in a subject. Is this expected? Could this be an consequence of the fact that PCAs are actually not meant for binary data? Or does this simply mean that subjects with more 1's resemble other subjects with more 1's (i.e., the more genetic variants are present, the more likely it is that those are the same variants as in another individuals with an approximately equal amount of genetic variants). Or could there be an alternative explanation? I hope the problem is well specified, otherwise please let me know! Many thanks!


title: Jensen-Shannon divergence calculation for 3 prob distributions: Is this ok?

text: I would like to calculate the jensen-shannon divergence for he following 3 distributions. Is the calculation below correct? (I followed the JSD formula from wikipedia):               P1  a:1/2  b:1/2    c:0     P2  a:0    b:1/10   c:9/10     P3  a:1/3  b:1/3    c:1/3     All distributions have equal weights, ie 1/3.          JSD(P1, P2, P3) = H[(1/6, 1/6, 0) + (0, 1/30, 9/30) + (1/9,1/9,1/9)] -                       [1/3*H[(1/2,1/2,0)] + 1/3*H[(0,1/10,9/10)] + 1/3*H[(1/3,1/3,1/3)]]          JSD(P1, P2, P3) = H[(1/6, 1/5, 9/30)] - [0 + 1/3*0.693 + 0] = 1.098-0.693 = 0.867      Thanks in advance... _**EDIT_** Here's some simple dirty Python code that calculates this as well:                   def entropy(prob_dist, base=math.e):             return -sum([p * math.log(p,base) for p in prob_dist if p != 0])              def jsd(prob_dists, base=math.e):             weight = 1/len(prob_dists) #all same weight             js_left = [0,0,0]             js_right = 0                 for pd in prob_dists:                 js_left[0] += pd[0]*weight                 js_left[1] += pd[1]*weight                 js_left[2] += pd[2]*weight                 js_right += weight*entropy(pd,base)             return entropy(js_left)-js_right          usage: jsd([[1/2,1/2,0],[0,1/10,9/10],[1/3,1/3,1/3]])


title: State Space models with Short Time Series

text: My problem is that I have a state space model that I estimate using the Berndt–Hall–Hall–Hausman (BHHH) algorithm. The state space model is relatively simple in that the hidden part follows a pure AR(1) process. From tests on simulated data I know I need a time series of 500 observations to get good convergence. However my empirical data set is only 375 observations, worse I need to be able to test the model in and out of sample so I have more like 188 observations!! Just wondered what my best options are? 1.) Use a different optimization method and hope that it is more efficient (seems labor intensive to try and find the right optimizer) 2.) Use bootstrapping similar to what they have done in the question below: Calculating confidence intervals via bootstrap on dependent observations Will the bootstrapping allow me to get better convergence with smaller data sets? If so how much shorter data set can I use? If 500 are required for good convergence will bootstrapping give me similar convergence with 350 observations or even 188 observations? Are there any other options I can try? Thanks Baz


title: What is a market basket problem?

text: Can somebody tell me or give a reference about the "Market Basket Problem" in data mining. What is this problem about? Any scientific reference with application will be greatly appreciated.


title: What are the differences between the linear regression and mixed models?

text: What are the key differences between the following two models?               lmefit = lmer(MathAch ~ SES + (1 |School) , MathScores)           lmfit = lm(MathAch ~ SES + factor(School) -1 , MathScores)       To me, they seem to be the same, except that lmefit takes less parameters (because it used Normal distribution to model the levels at the group level...) Am I right? * * * And what's the difference between these two models?               M0 <- lmer (y ~ 1 + (1 | county))     M1 <- lmer (y ~ -1 + (1 | county))


title: Regression: testing two predicted means on the same slope?

text: I'm surprised I haven't run into this question yet, but I was wondering if there is a sound way of testing the difference in predicted means between two points on the _same_ slope. For example: Let's say my predictor is the age of a car, and my predicted variable, its price. Let's assume that the downward slope is significant. What I would like to test is if a car 5 years old is significantly cheaper than one that is a year old. Can I just use the two predicted means, and divide them by some error term as a t-test? And if so, which error term? Thanks!


title: Resources and R code for analysing advanced experimental designs?

text: I switched from `SAS` to `R`. I was using "SAS Code for Some Advanced Experimental Designs" as a reference whenever I had data from experimental design studies. With advanced experimental designs I always have difficult to come up with right `R` code. Specifically I have a tough time with mixed effects analysis. I couldn't find any comprehensive material for analyzing advance experimental designs in `R`. I'd highly appreciate if someone point out any such material. It'd be nice if someone translate these SAS Codes for Some Advanced Experimental Designs to R Code for Some Advanced Experimental Designs.


title: Which robust correlation methods are actually used?

text: I plan to do a simulation study where I compare the performance of several robust correlation techniques with different distributions (skewed, with outliers, etc.). With _robust_ , I mean the ideal case of being robust against a) skewed distributions, b) outliers, and c) heavy tails. Along with the Pearson correlation as a baseline, I was thinking to include following more robust measures:   * Spearman's $\rho$   * Percentage bend correlation (Wilcox, 1994, [1])   * Minimum volume ellipsoid, minimum covariance determinant (`cov.mve`/ `cov.mcd` with the `cor=TRUE` option)   * Probably, the winsorized correlation Of course there are many more options (especially if you include robust regression techniques as well), but I want to restrict myself to the mostly used/ mostly promising approaches. **Now I have three questions (feel free to answer only single ones):**   1. **Are there other robust correlational methods I could/ should include?**   2. **Which robust correlation techniques are** _**actually_** **used in your field?** (Speaking for psychological research: Except Spearman's $\rho$, I have never seen any robust correlation technique outside of a technical paper. Bootstrapping is getting more and more popular, but other robust statistics are more or less non-existent so far).   3. **Are there already systematical comparisons of multiple correlation techniques that you know of?** Also feel free to comment the list of methods given above. * * * [1] Wilcox, R. R. (1994). The percentage bend correlation coefficient. _Psychometrika_ , 59, 601-616.


title: Deriving total (within class + between class) scatter matrix

text: I was fiddling with PCA and LDA methods and I am stuck at a point, I have a feeling that it is so simple that I can't see it. Within-class ($S_W$) and between-class ($S_B$) scatter matrices are defined as: $$ S_W = \sum_{i=1}^C\sum_{t=1}^N(x_t^i - \mu_i)(x_t^i - \mu_i)^T $$ $$ S_B = \sum_{i=1}^CN(\mu_i-\mu)(\mu_i-\mu)^T $$ Total scatter matrix $S_T$ is given as: $$ S_T = \sum_{i=1}^C\sum_{t=1}^N(x_t^i - \mu)(x_t^i - \mu)^T = S_W + S_B $$ where C is number of classes and N is number of samples $x$ are samples, $\mu_i$ is ith class mean, $\mu$ is overall mean. While trying to derive $S_T$ I came up to a point where I had: $$ (x-\mu_i)(\mu_i-\mu)^T + (\mu_i-\mu)(x-\mu_i)^T $$ as a term. This needs to be zero, but why?


I have a question about interaction terms in logistic regression: I find an interaction effect between two binary variables. One of the four categories (0*0) consists however of only 15 respondents (5% of the sample). Is this a problem?


title: Order of GARCH model

text: I am Ashenafi I would like to ask, Is it possible to know the maximum order of GARCH(p, q); value of p and q (For example, using R software)automatically? Example: In ArchTest, I got some where the result like the following                                  Chi-squared = 4.3044, df = 1, p-value = 0.0380                             Chi-squared = 4.3892, df = 2, p-value = 0.1114                             Chi-squared = 4.3927, df = 3, p-value = 0.2221                             Chi-squared = 4.4637, df = 12, p-value = 0.9735       Therefore, Is it possible to know the value of p and q? Sincerely Ashenafi,


title: Is the difference between conversion rates significant?

text: We ran advertising to a certain group of demographic profiles in one single television region. The advertising ran for 4 weeks. We measured the conversion rate for the period of the advertising and for the prior 4 weeks. We were then able to measure the difference in conversion rates. I would like to know if the difference in the conversion rate change from the demographic groups in the region we ran the advertising compared to the conversion rate from the same demographic groups in the other regions is significant. The data looks like Region Conversion Rate Change from Demographic Groups A (one in which we ran advertising) 11.0% B 6.4% C -12.7% D -23.4% E -18.6% F 0.1% My question is what is the appropriate statistical test to determine if the improvement in region A's conversion rate from the demographic profiles is significant?


title: What is the current research on missing not at random data?

text: I have a data set with N ~ 5000 and about 1/2 missing on at least one important variable. Comparing the people with missing data to those with complete data it is fairly clear that the data are not missing at random (nonignorable nonresponse). The missing data pattern is not monotone missing and at least some data is missing on about a dozen variables. The main analytic method will be Cox proportional hazards. I am using `SAS` which now offers an `MNAR` method for such data. If the pattern is monotone, then it offers options to estimate the missing data with either nearest neighbors or complete cases. If the pattern is non-monotone, then it offers only a method of 'adjusting' the parameters, which seems (from what I understand) like a pure sensitivity analysis. `SAS` also offers a MCMC method of creating a monotone missing pattern from a non-monotone pattern. My current plan is to first create a monotone missing pattern and then apply nearest neighbor, then analyze the multiply imputed data. However, I am not sure this is best, nor what determines which options to choose in a scenario like this. Advice welcome as are references to the literature.


title: Redundant feature generation

text: I'm evaluating a few supervised feature selection algorithms, with a focus on redundancy. As a base correlation statistic, I'm using Mutual Information - so discrete variables are also a focus for simplicity. Does anyone know any good methods of redundant feature generation for discrete data? Either from fully simulated data or, preferably, inserting new redundant features into a 'real-world' dataset. Some simple methods would be:   * Taking an existing variable and changing values with some probability distribution.   * Binarizing an existing variable by setting values to either 0 or 1   * Fusing two or more existing variables


title: Logistic Regression in R - Steps and Output

text: I am doing statistics for the first time in my life and I am not quite sure what to include and how to interpret the results. I am doing a logistic regression in R. Here is what I have so far:   1. `GLM` with family = binomial (dependent ~ indep1 + indep2 + ...+ indep7 +0) If I dont include the 0 I get NA for my last independent variable in the summary output..   2. `Update` the model (indep2 has a p-value > 0.05 and is left out)   3. I am applying anova              anova(original_model,updated_model, test="Chisq")             Resid.Df  Resid.Dev Df Deviance Pr(>Chi)     1     34067      18078                          2     34066      18075  1   2.4137   0.1203      Here I am not sure how to interpret it. What tells me if the simplification of the model is significant? the p-value is with 0.12 bigger than 0.05, does this mean that the simplification is not significant?   4. make a cross-table (compare predicted (probability >0.5) - observed)              fit           FALSE  TRUE       No  30572    68       yes  3407    31      I'd say that 31 values are predicted correctly (yes-true), resp 68 (no-true) but that most values are classified wrong, which means that the model is really bad?   5. then I make a wald test for each independent variable for the first independent variable it would look like this:              > wald.test(b = coef(model_updated), Sigma = vcov(model_updated), Terms     > = 1:1)      here I only look if the p-values are significant and if they are it means that all variables contribute significantly to the predictive ability of the model   6. I calculate the odds with their confidence intervals (this is basically exp(estimate)              oddsCI <- exp(cbind(OR = coef(model_updated), confint(model_updated)))      For all odds smaller than 1 i do 1/odd              Estimate        Odds Ratio      Inverse Odds     -0.000203       0.999801041     1.000198999      0.000332       1.000326571     odd bigger than 1     -0.000133       0.999846418     1.000153605     -3.48       0.008696665     114.9866056     -4.85       0.029747223     33.61658319     -2.37       0.000438382     2281.113996     -8.16       0.110348634     9.062187402     -2.93       0.062668509     15.95697759     -3.65       0.020156889     49.61083057     -5.45       0.033996464     29.41482359     -4.02       0.004837987     206.6975334      This O would interpret like that for the "odd bigger than 1" the case is over 1 times more likely to occur. (Is is incorrect to say that, or not?) Or for the last row you could say that t for every subtraction of a unit, the odds for the case to appear decreases by a factor of 206.   7. Then I look at               with(model_updated, null.deviance - deviance) #deviance     with(model_updated, df.null - df.residsual) #degrees of freedom      # pvalue     with(Amodel_updated, pchisq(null.deviance - deviance, df.null - df.residual,      lower.tail = FALSE))     logLik(model_updated)      But I don't really know what this tells me.   8. In a last step I do              stepAIC(model_updated, direction="both")      but also here I don't know how to interpret the outcome. I see that it looks at all interactions between my independent variables but I don't know what it tells me. After this, I can make a prediction by using the updated model and by separating it into training data and validation data I suppose?


title: Image similarity database

text: I'll preface this by stating I'm very new to computer vision. I took a Stanford Machine Learning course on Coursera, and that is the extent of my experience. I'm working on a free-time project to learn more about the subject. I fully acknowledge the steep learning curve. What I'm trying to do is to create a database of record covers that I can use to take a photo of a record and search the database to see if I can find the same album cover in it. The photo will not match what is in the database exactly. Expected stored images and search input images would look something like this; http://imgur.com/j14FBWM,1d2CxX1 where the first image is what is in the database, and the second is likely taken from an app on a phone. I figure first I'll need something that can simply recognize the square shape of the record in the search image and translate that to a proper square before doing a search. I'm looking at pHash right now, which sounds like it might be able to handle the second step. But I'd have to pull the cover out of the source first. I'm looking for any tools or libraries, algorithms, suggestions or any other input you are willing to give. If you think there's a better tool than pHash, I'd be interested in that too.


title: How to model occurrence counts in groups with latent membership?

text: I have a dataset that describes the daily history of occurrences of a certain phenomenon P among a certain population. (These are subdivided into certain classes or forms that the phenomenon can take: for every day, we have numbers of occurrences of P_1, P_2, ..., P_N.) The phenomenon is rather rare: more than 90% of days it is not observed at all. The dataset also contains information about attendance of certain members of that same population to recurring sessions of treatment. The sessions are organized approximately twice a week, following a pre-defined schedule. They are held in two groups, each group having its own schedule. Treatment is divided into three phases; at any given time, each group might be in a different phase from the other. The population changes as treatment finishes in phase 3 and new subjects are enrolled. We are interested in exploring relationships between the attendance and phase of the treatment and occurrences of the phenomenon. I have started by just computing Pearson correlation between number of attendees of the treatment sessions and number of occurrences of the phenomenon and found the correlation to be not significant. In more than 70% of days (I call them "zero days"), there is both zero attendance (treatment not taking place) and zero occurrences of the phenomenon. When I remove these and recompute the correlation on non-zero days only, the correlation comes out as high and statistically significant, but I think this is not a valid operation to perform as part of data cleaning. The problem is that there is no information on group membership of the subjects; only the number of people who attended treatment in each group and the phase of the treatment for every day is known. It is not known which members of population exhibited the phenomenon, whether they were subject to treatment, or if so which groups they were treated in. There are only counts referring to the population as a whole. I thought to use a zero-inflated Poisson model for this dataset to model the occurrences of the phenomenon, but I'm not sure how to take the above- mentioned issues into account, or indeed if they can be overcome. Can anybody offer a suggestion on what kind of models or statistical tests would be suitable for this kind of data?


title: Simulating outcome of 3 political parties

text: First I'm sorry I couldn't figure out the most accurate title for this question (suggestions welcome). Here's the case: I want to implement spinners like the ones on this page: http://www.nytimes.com/newsgraphics/2014/senate-model/ I more or less get the idea: for each state I will need a separate instance of a Gaussian random generator (each with its own mean and standard deviation). At each spin, I will need to execute every single one of them to obtain the outcome. I learned about such a function from this page: http://www.protonfish.com/random.shtml I will only need to know the probability of one party (P(Dem) for example). The other one (P(Rep)) is simply 1 - P(Dem). But then I figured it will only work for two-parties situation. What about 3-parties situation (let's say: Nationalist, Religious, and Communist)? Will the same model work (I have the gut feeling it wouldn't) ? Would you suggest which approach should I take? Or lead me directly to a chapter / topic in statistic that helps me understand better the problem?


title: Calculating conditional probabilities given a bivariate gaussian

text: This is a continuation of my previous question. * * * I have two classes, $C_1$ and $C_2$. $C_1$ is a bivariate Gaussian with mean $\mu = (0,0)$ and covariance $\Sigma = I$ $C_2$ is a bivariate Gaussian with mean $\mu = (1,3)$ and covariance $\Sigma = 2I$, where $I$ is the identity matrix. I am trying to calculate $P(x|C_1)$ and $P(x|C_2)$ so I can eventually calculate $P(C|x)$ To calculate $P(x|C_1)$ and $P(x|C_2)$ this I'm using the formula for a bivariate normal distribution found here. My covariance is zero, which makes this a little bit easier. **When I use this calculate $P(x|C_1)$ I have...** $z=x_1^2 - x_2^2$ $p=0$ $p(x_1,x_2) = (\frac{1}{2\pi})e^{-z/2}$ **When I use this to calculate $P(x|C_2)$ I have...** $z=\frac{1}{4}( (x_1-1)^2 + (x_2-3)^2 )$ $p=0$ $p(x_1,x_2) = (\frac{1}{8\pi})e^{z/8}$ Did I do this correctly? Also, I'm a bit confused as to whether what I'm doing even gives me $P(x|C_1)$ and $P(x|C_2)$. I'm a bit over my head in the class I'm in, so if I'm totally wrong there please correct me. Anyways, with these two values I'm supposed to calculate $P(C|x)$ using Bayes rule (I think). I have the priors of $C_1$ and $C_2$ (they are 0.4 and 0.6 respectively), but I'm lost exactly on how to calculate $P(C|x)$ with this. **Could somebody basically just check over some of my work and help me out with the process using Bayes rule to calculate $P(C|x)$?** EDIT: My end goal here is to calculate an optimal decision boundary between $C_1$ and $C_2$, if I'm going the complete wrong way here let me know, but from what I've gathered through time spent on this site, I think I'm headed the right way.


title: Logistic Regression: should I worry about coefficient correlation matrix?

text: With SPSS (v18), in the output of a logistic regression after the "Variables in the Equation" section, there is a Correlation Matrix which represents the correlation of the coefficients. What do high correlations imply? Should I try to avoid high values by removing variables? My AUC drops severely when I do this. I have found a similar Cross Validated question "interpretation of correlation matrix of regression coefficients in Cox regression?" which has no answer...


title: How do I optimize this scenario to achieve the lowest score?

text: Here is the scenario: I am iteratively seeking a goal state by selecting from one of a set of sets of options. The set of sets defines the possible routes to the goal state, and each inner set defines the probability I will succeed, and what the failure penalties and probabilities are. A failure costs me a variable amount of score, which is also a given. My objective is not only to achieve the goal state, but to do so with minimum score. For example, say I can select from two groups: In group 1, I have a 25% chance of success. However, I have a 50% chance of failing and adding 5 to my score, and a 25% chance of failing and adding 10 to my score. In group 2, I have a 50% chance of success. However, I have a 10% chance of failing and adding 5 to my score, a 30% chance of failing and adding 15 to my score, and a 10% chance of failing and adding 25 to my score. If I wish to succeed while minimizing my score, which of these two groups should I repeatedly select from? I'm looking for a general solution to this problem, where there may be an arbitrary number of groups, and each group has an arbitrary number of outcomes. As in the example, for each group, their relative probabilities of success and failure are known, along with the points added for failing. I _think_ I want to optimize for "mean score for success" in the group in order to minimize my overall score. However, I'm a bit confused because it seems to me that if the distribution of possible final scores is flat enough, I may be better off accepting a higher mean score from my selected group.


title: Generating time series with copula

text: I want to generate univariate time series in R using copulas. Let $X_1, X_2, \ldots X_N,\ldots$ denote random variables such that joint distribution of $X_t$ and $X_{t-1}$ is $C(F(X_t), F(X_{t-1}))$, where $F$ denotes a marginal distribution function, and $C$ a copula. I am a beginner in R, and at present I face some difficulties understanding the syntax. I think I should generate a first observation, then with the first observation use the conditional cdf and generate second observation, and so on. But I don't understand how to construct the conditional distribution (the function `ccopula` from the copula package is not working properly). I would be glad for any help.


title: Random forest cross validation for feature selection, imbalanced datasets

text: I have an 5297X26 imbalanced dataset, the class1 has 588 samples and class2 has 4709 samples. I used the following code to perform random forest:               rfp<-randomForest(label~.,data=data,importance=TRUE,proximity=TRUE,replace=TRUE,sampsize=c(588,588))          Thus I could solve the imbalanced problem by selecting 588 samples from each class in each iteration. But I also want to perform cross validation for feature selection. The function I plan to use is rfcv .I tried to add sampsize=c(588,588) to the function but it didn't work. How to perform the sampling in this function?


title: Compare Actual to Predicted

text: I would like to compare a series of actual values to predicted values in order validate the model used to generate the predictions. The model is a simple straight line and was derived from a regression analysis. The R squared value for the straight line is not brilliant at 0.62. What is the best method to use to determine if the actual values fall with the predicted taking into account the margin of error with the model? I have Excel and Minitab v16


title: Model uniformly distributed response in cumulative link (logit/probit) model

text: I have an ordered response variable (stated preference from strong dislike to strong like) that is virtually uniformly distributed. I am using package ordinal in R to estimate the probability of an observation matrix X to fit into a specific ordered category Y. Christensen (2013, p18), in his explanation of the probit model, states that the observed response distribution should not deviate too much from a bell-shaped curve, else it ruins the quality of the estimates. Given that I have about as many observations in each category, my response clearly does not meet this criteria. I'm not sure how to best deal with the data then. Any suggestion would be appreciated! At the moment, when using a probit or logit model, my errors are almost perfectly normally distributed but only 1/3 observations gets assigned to the right group, which is a rather poor data fit. This might obviously mean that my predictors are just no good but I think the assumption of a normal response decreases the likelihood of correct predictions on the extremes, which gives my error distribution pretty fat tails.


title: Questionnaire design - Agree/Disagree scale

text: I recently constructed a questionnaire to carry out a market survey on soap brands used. One of the questions which was asked was to rate few statements which influence the respondent to buy a certain brand. For each of the statements(Eg: price, advertisements etc) I gave four options strongly agree, agree, disagree, strongly disagree. At that time I decided not to add a "moderate" option. However, I am halfway through collecting the data an I am seriously worried if I had made a big mistake by going with a rating system with even number of responses. Is this okay? At the end of the survey I will have to justify why I did so hence I am in need for some good answers to justify my decision. I would really appreciate if someone could help me with this. Thanks in advance


title: Can cluster analysis cluster variables that both positively and negatively correlate with each other?

text: I received this question by email from a Neuroscience PhD student. > I would greatly appreciate if you could please let me know whether Factor > Analysis could load positively and inversely correlated variables onto the > same latent factor, whereas Cluster Analysis can only cluster into the same > factor either positively or inversely correlated variables (e.g., they would > be segregated into 2 factors).


title: When does the amount of skew or prevalence of outliers make the median preferable to the mean?

text: I know that "deviations in the data are devil", and when the distribution is highly skewed, it is better to consider median as average rather than mean, but how to decide these hard-limits. For example:   * CASE 1:     * Assume X = 10,20,30,40,50,60,70     * In this case, I think that it is better to use mean and that it will give very accurate results.   * CASE 2:     * Assume X = 10,20,30,40,50,60,70,7000     * In this case, I think that it is better to use median instead of using the mean.   * CASE 3:      * Assume X = 10,20,30,400,500,600,700     * In this case, I think it is better to use IQR (Inter Quartile Range) But I'm stuck with how to decide these hard-limits i.e. which to use in which condition, in general. I've found a tool working on subjected principle, which takes context-less sample-distribution as input and determines whether mean is close/moderate or against the null-hypothesis. Find References:-   * http://home.ubalt.edu/ntsbarsh/Business-stat/otherapplets/MeanTest.htm (For Mean Test)   * http://home.ubalt.edu/ntsbarsh/Business-stat/otherapplets/MediansTest.htm (For Median Test) What I'm really looking is a good answer which states how to derive these conclusions. Any help or learning pointers will be highly appreciated.


title: What is the meaning of rank in the context of change-detection?

text: In a technique that uses CUSUM for change-point detection in this paper, the first step is given below: > Let $x_1, x_2,..., x_n$ be the $n$ samples in an event-series. The samples > are ranked in increasing order and the rank $r_i$ for each sample is > calculated. In case of ties, we assign average rank to each sample. The > cumulative sums are computed as: $S_i = S_{i-1} + (r_i - \bar{r})$. If the > ranks are randomly distributed, then there is no change-point. However, if > there is indeed a change-point in the event-series, then higher ranks should > dominate in either the earlier or later part of the event-series. I am not quite sure I understand the meaning of rank in this context. For instance, if n=10 and my data points are: 5, 2, 4, 1, 9, 2, 9, 2, 10, 1 can someone please clarify what is really being done here?   * Sort in increasing order: 1, 1, 2, 2, 2, 4, 5, 9, 9, 10   * Assign Average Ranks to break ties: What does this step mean?   * Assign Final Ranks: What does this step mean?


title: Summarising wind vector statistics over time

text: My problem is to figure out an algorithm for summarising wind statistics over an area. The statistics are (24x) hourly but I want to end up with 2-3 'blocks' that summarise the conditions for the whole day. My data is in the form of separate magnitude and direction NumPy arrays but I could convert them to u,v or polar coordinates or whatever if it made the solution easier. There are other conditions like: each block must be at least 3 hours long, but the general problem is the trade off between accuracy and precision. And in a way, to maximise the similarity of the statistics inside each 'block' and the difference between each 'block'. Anyway I don't have a background in statistics and am wondering if there is an array-based approach I could use, maybe involving calculating minimum bounding box/convex hull, that kind of thing? I'm not looking for a complete algorithm but some hints as to what I should be searching for to do further research.


title: Resume buzz words

text: It has come to my attention that companies will rank resumes based on buzz words and only look at those that have high scores assuming enough people submit for the position. I don't like this system but still want a job, so instead of inserting these meaningless terms into my actual resume I've compiled a list of such terms and am going to incorporate the list into my resume in a text box with white text. The problem is that I don't know how those scoring engines work (is there an application most employers use?) and so I can't be sure if I should maximize the terms I put in, if I should repeat terms, or if I should maximize the effective terms put in because number of words matter. I have entered my resume on Rezscore with and without the buzz word metadata and found my resume score significantly lowered with buzz words because it was so verbose. Does anyone have more knowledge of these resume ranking systems?


title: Ways to understand 2-dimensional time-series data

text: I'm working on 2D time series data where two attributes are depth and temperature. When I plotted depth-vs-temp curve and saw its variation over time, the fluctuation occurs at few places only. I'm not saying temperature is dependent on depth. But given these data, what should I look into to establish relationships between depth, temperature, time? Ideas:   Look for depth regions where fluctuation in temperatures occurs over time. And develop time series prediction for these given regions. What models/papers should I study to get insights out of the data? I want to understand and study what are the various methods that could be tested and visualizations built over these data. It's a data exploration problem. Data shape: (2000,2,100) (samples,depth-temperature,time-sample) .i.e 2000 depth,temperature samples for each of 100 time stamps. Thanks.


title: Variance of the product of a random matrix and a random vector

text: If $X$ and $Y$ are independent random variables, then the variance of the product $XY$ is given by $ \mathbb{V}\left(XY\right)=\left\{ \mathbb{E}\left(X\right)\right\} ^{2}\mathbb{V}\left(Y\right)+\left\{ \mathbb{E}\left(Y\right)\right\} ^{2}\mathbb{V}\left(X\right)+\mathbb{V}\left(X\right)\mathbb{V}\left(Y\right) $ If $\mathbf{X}$ and $\mathbf{y}$ are independent matrix and vector of $m\times m$ and $m\times1$ dimension respectively, then what would be the variance of the product $\mathbf{X}\mathbf{y}$? **My Attempt** $ \mathbb{V}\left(\mathbf{X}\mathbf{y}\right)=\mathbb{E}\left(\mathbf{X}\right)\mathbb{V}\left(\mathbf{y}\right)\left\{ \mathbb{E}\left(\mathbf{X}\right)\right\} ^{\prime}+\left\{ \mathbb{E}\left(\mathbf{y}\right)\otimes\mathbf{I}_{m}\right\} ^{\prime}\mathbb{V}\left\{ \textrm{vec}\left(\mathbf{X}\right)\right\} \left\{ \mathbb{E}\left(\mathbf{y}\right)\otimes\mathbf{I}_{m}\right\} +\mathbb{V}\left\{ \textrm{vec}\left(\mathbf{X}\right)\right\} \left\{ \mathbb{V}\left(\mathbf{y}\right)\otimes\mathbf{I}_{m}\right\} $ I know this is not right, at least the last term is wrong. I'd highly appreciate if you give me the right identity or point out any reference. Thanks in advance for your help and time.


title: What is numerical overflow?

text: I came across an error of numerical overflow when running a maximum likelihood estimation on a log-linear specification. What does numerical overflow mean?


title: Techniques for Handling Incomplete/Missing Data

text: My question is directed to techniques to deal with _incomplete_ data during the classifier/model training/fitting. For instance, in a dataset w/ a few hundred rows, each row having let's say five dimensions and a class label as the last item, most data points will look like this: [0.74, 0.39, 0.14, 0.33, 0.34, 0] A few might look something like this: **[0.21, 0.68, ?, 0.82, 0.58, 1]** So it's those types of data points that are the focus of this Question. My initial reason for asking this question was a problem directly in front of me; however, before posting my Question, i thought it might be more useful if i re-phrased it so the answers would be useful to a larger portion of the Community. As a simple heuristic, let's divide these data-handling techniques based on when during the processing flow they are employed--before input to the classifier or during (i.e., the technique is inside the classifier). The best example i can think of for the latter is the clever 'three-way branching' technique used in Decision Trees. No doubt, the former category is far larger. The techniques i am aware of all fall into one of the groups below. While recently reviewing my personal notes on "missing data handling" i noticed that i had quite an impressive list of techniques. I just maintain these notes for general peace of mind and in case a junior colleague asks me how to deal with missing data. In actual practice, i don't actually use any of them, except for the last one.   1. **Imputation** : a broad rubric for a set of techniques which whose common denominator (i believe) is that the missing data is supplied directly by the same data set--substitution rather than estimation/prediction.   2. **Reconstruction** : estimate the missing data points using an auto-associative network (just a neural network in which the sizes of the input and output layers are equal--in other words, the output has the same dimension as the input); the idea here is to train this network on complete data, then feed it incomplete patterns, and read the missing values from the output nodes.   3. **Bootstrapping** : (no summary necessary i shouldn't think, given it's use elsewhere in statistical analysis).   4. **Denial** : quietly remove the data points with missing/corrupt elements from your training set and pretend they never existed.


title: Robust t-test for mean

text: I am trying to test the null $E[X] = 0$, against the local alternative $E[X] > 0$, for a random variable $X$, subject to mild to medium skew and kurtosis of the random variable. Following suggestions by Wilcox in 'Introduction to Robust Estimation and Hypothesis Testing', I have looked at tests based on the trimmed mean, the median, as well as the M-estimator of location (Wilcox' "one-step" procedure). These robust tests do outperform the standard t-test, in terms of power, when testing with a distribution that is non-skewed, but leptokurtotic. However, when testing with a distribution that is skewed, these one-sided tests are either far too liberal or far too conservative under the null hypothesis, depending on whether the distribution is left- or right-skewed, respectively. For example, with 1000 observations, the test based on the median will actually reject ~40% of the time, at the nominal 5% level. The reason for this is obvious: for skewed distributions, the median and the mean are rather different. However, in my application, I really need to test the mean, not the median, not the trimmed mean. **Is there a more robust version of the t-test that actually tests for the mean, but is impervious to skew and kurtosis?** Ideally the procedure would work well in the no-skew, high-kurtosis case as well. The 'one-step' test is almost good enough, with the 'bend' parameter set relatively high, but it is less powerful than the trimmed mean tests when there is no skew, and has some troubles maintaining the nominal level of rejects under skew. **background:** the reason I really care about the mean, and not the median, is that the test would be used in a financial application. For example, if you wanted to test whether a portfolio had positive expected log returns, the mean is actually appropriate because if you invest in the portfolio, you will experience all the returns (which is the mean times the number of samples), instead of $n$ duplicates of the median. That is, I really care about the sum of $n$ draws from the R.V. $X$.


title: Comparing multiple contingency tables, independent data

text: (migrating from math overflow, where no answers were posted) suppose I have $K$ different methods for forecasting a binary random variable, which I test on independent sets of data, resulting in $K$ contingency tables of values $n_{ijk}$ for $i,j=1,2$ and $k=1,2,...,K$. How can I compare these methods based on the contingency tables? The general case would be nice, but $K=2$ is also very interesting. I can think of a few approaches:   * compute _some_ statistic on each of the tables, and compare those random variables (I'm not sure if this is a standard problem or not),   * something like Goodman's improvement of Stouffer's method, but I cannot access this paper, and was hoping for something a little more recent (more likely to have the latest-greatest, plus computer simulations). any ideas?


title: Suggest a method for statistical data compression

text: There's a lot of work done in statistics, while state-of-art in lossless data compression is apparently this: http://mattmahoney.net/dc/dce.html#Section_4 Please suggest good methods/models applicable for data compression. To be specific:   1) How to estimate the probability of the next bit in a bit string?   2) How to integrate predictions of different models?   Update: > you should include a better description of what data you want to compress... Why, I'm talking about universal compression obviously. For data with known structure its really not a mathematical problem, so there's no sense to discuss it here. In other words, the first question is: given a string of bits, what do we do to determine the probability of the next bit, as precisely as possible? > otherwise we will have 10 different answers trying to summarize different > part of the huge theory of compression I'd written quite a few statistical compressors, and I'm not interested in that. I'm asking how a statistician would approach this task, detect correlations in given data, and compute a probability estimation for the next bit. > In addition, the two point you give to be more specific are not detailed > enough to be understood. What's not detailed in there? I'm even talking about bits, not some vague "symbols". I'd note though, that I'm talking about "probability of a bit" because computing a probability of bit==0 or bit==1 is a matter of convenience. Also, I'm obviously not talking about some "random data compression", or methods with infinite complexity, like "Kolmogorov compression". Again, I want to know how a good statistician would approach this problem, given a string of bits. Here's an example, if you need one: hxxp://encode.ru/threads/482-Bit-guessing- game