title: Correlated historical time-series with missing data

text: I am trying to propose a method of determining whether health quality changed between 1700-1820 (roughly) in the United States. Some available data/time- series include: mortality rates, height, life expectancy, and weight - which are all correlated to varying degrees. There are gaps in the time series, and some only cover specific regions. This is for a paper for a third year economic history class. The purpose is to test the hypothesis of a paper that claims there was little to no economic growth during the period, and that food consumption and productivity was near constant. I propose that these claims could be tested by determining whether or not there was a change in overall health quality (if health greatly improved nutrition probably had some role, and productivity would have increased, for example). My first thought was some sort of yearly health index that aggregates the data. The problem is that there may be a ten year period with no mortality rates and only height data, for example (but I might know the relationship between height and mortality, or can estimate it). What sort of statistical methods might be useful in painting a larger picture of health quality? I would just like to propose some ways researchers could overcome the problem of inconsistent data in order to decide if individuals were healthier in 1820 than in 1700 (discover the long term trend). I do not have to do any actual computation.


title: Should plots for time series use a logarithmic scale for duration/time?

text: I was reading this tutorial on `ggplot2` and was struck by a bit where it said that > durations are typically best displayed along a logarithmic scale Is this an established principle? The author states the claim as if it were, but I've never seen them plotted like this (at least not in phonetics), and to me, a linear scale makes more sense (time moving at a constant rate and all). If it is an established principle, I would appreciate some references and maybe an explanation of why this is the case.


title: Literature on nonparametric density estimation

text: I am about to write my bachelor thesis about non-parametric density estimation, especially kernel density estimators and their application in classification. As I am quite new to looking for academic literature, I am having a hard time finding the most important and modern papers, or other resources, and would be glad if someone could give me a hint. Right now I am mainly working with older works (especially by Silverman and Devroye).


title: Significant result with one-way anova test but not with Scheffe test?

text: I conducted one-way anova test and get a significant result, but when I conducted Scheffe test there was no significant result. Is this possible? How should I proceed or report the results?


title: Time-varying Coefficients

text: I have time series data on fish catches from 1950-2011. I wish to implement a regression model with varying coefficients. I'm aware that cox models etc. exist and implementation via the `survival` package in R. My data is not survival data, it's just several variables with fish catches and year. Is there a way in R to implement such models? I've yet to come across this but I don't think it's unreasonable to want to model such data without it being survival data. I want to model `inlandfao` from `marinefao`. Here is my data and some plots:               fishdata <- read.csv("http://dl.dropbox.com/s/4w0utkqdhqribl4/fishdata.csv", header=T)          require(reshape2)     require(ggplot2)     theme_set(theme_bw())     require(scales)          df2 <- data.frame(cbind(year,totalmarinefao, totalinlandfao))     df2     dd <- melt(df2, id.vars = "year")     dd     pp <- ggplot(dd, aes(year, value, colour=variable)) + geom_point() + geom_line(size=1)     pp_final <- pp +  xlab("Year") + ylab("Catches (Tons)") + ggtitle("Time Series of Variables (1950-2011)")      pp_final     pp_final2 <- pp_final +  scale_colour_discrete(name = "Variable - Catches (FAO)", breaks = c("totalmarinefao", "totalinlandfao"),                                                    labels=c("Marine", "Inland")) +        scale_shape_discrete(name = "Variable (FAO)", breaks = c("totalmarinefao", "totalinlandfao"), labels=c("Marine", "Inland")) +        scale_x_continuous(breaks=seq(1950,2011,10)) + scale_y_continuous(labels=comma)          pp_final2     pp_3 <- pp_final2 + theme(axis.text.x  = element_text(vjust=1, size=16)) + theme(axis.title.x = element_text(size=20))     pp_4 <- pp_3 + theme(axis.text.y = element_text(vjust=0, size=16)) + theme(axis.title.y = element_text(size=20, vjust=0.2))     pp_5 <- pp_4 + theme(plot.title = element_text(lineheight=.8, face="bold", size=20))     pp_5          qplot(marinefao, inlandfao, data=fishdata, main="Scatterplot of the Marine & \n Inland      fish Catches (Tons)", xlab="Marine Catches", ylab="Inland Catches") +      scale_x_continuous(labels = comma) + scale_y_continuous(labels = comma)      From these plots, a linear model isn't appropriate. I have fitted GAMs etc. to these data. Let me more if you require details.


title: Getting the next observation from a HMM gaussian mixture distribution

text: I have a continuous univariate xts object of length 1000, which I have converted into a data.frame called `x` to be used by the package `RHmm`. I have already chosen that there are going to be 5 states and 4 gaussian distributions in the mixed distribution. **What I'm after is the expected mean value for the next observation. How do I go about getting that?** So what I have so far is:   1. a transition matrix from running the `HMMFit()` function   2. a set of means and variances for each of the gaussian distributions in the mixture, along with their respective proportions, all of which was also generated form the `HMMFit()` function   3. a list of past hidden states relating to the input data when using the output of the HMMFit function and putting it into the `viterbi` function How would I go about getting the next hidden state (i.e. the 1001st value) from what I've got, and then using it to get the weighted mean from the gaussian distributions. I think I'm pretty close just not too sure what the next part is...The last state is state 5, do I use the 5th row in the transition matrix somehow to get the next state? All I'm after is the weighted mean for what is to be expect in the next observation, so the next hidden state isn't even necessary. Do I multiply the probabilities in row 5 by each of the means, weighted to their proportion for each state? and then sum it all together? here is the code I used.               # have used 2000 iterations to ensure convergence     a <- HMMFit(x, nStates=5, nMixt=4, dis="MIXTURE", control=list(iter=2000)     v <- viterbi(a,x)     a     v      As always any help would be greatly appreciated! p.s. this is a copy from stack overflow http://stackoverflow.com/questions/11644522/getting-the-next-observation-from- a-hmm-gaussian-mixture-distribution


title: How to calculate with tiny probabilities and large samples?

text: Is it possible to calculate or approximate the probability of something extremely unlikely happening once over a large sample, i.e., in situations where the probability is smaller than machine error? For example, I was trying to calculate the approximate likelihood of someone sharing my genome. Apparently an individual genome can be losslessly compressed to about 4MB (2^25 bits). So the chance of one of the ~7 billion humans on the planet sharing my genome is about: $(1-\frac{1}{2^{2^{25}}})^{(7\times10^9)}$. Or using the Birthday problem method, the likelihood of there being at least two people with identical genomes (let's ignore twins and so forth) is something like: $\frac{(7\times10^9)!\cdot{{2^{2^{25}}} \choose 7\times10^9}}{({2^{2^{25}})}^{(7\times10^9)}}$ The problem here is that the numbers are so small or so large that it's impossible to guess roughly where they'd lie. So, is there any way of approximating these or similar calculations? I realise that in some cases the assumptions behind the problems might be out by multiple orders of magnitude, but even being able to approximate to "more likely than not" would be interesting.


title: Increase sample size for significant correlation

text: At the moment I have $5$ paired samples for correlation. Spearman's R is $0.2$ and the $p = 0.78$. How do I calculate the number to extra samples I would need to get a more significant p-value?


title: Does a canonical link function always exist for a Generalized Linear Model (GLM)?

text: In GLM, assuming a scalar $Y$ and $\theta$ for the underlying distribution with p.d.f. $$f_Y(y | \theta, \tau) = h(y,\tau) \exp{\left(\frac{\theta y - A(\theta)}{d(\tau)} \right)}$$ It can be shown that $ \mu = \operatorname{E}(Y) = A'(\theta)$. If the link function $g(\cdot)$ satisfies the following, $$g(\mu)=\theta = X'\beta $$ where $X'\beta$ is the linear predictor, then $g(\cdot)$ is called the canonical link function for this model. My question is, does a canonical link function always exist for a GLM? In other words, can $A'(\theta)$ always be inverted? What are the necessary conditions for a canonical link function to exist?


title: What is the hardest statistical concept to grasp?

text: This is a similar question to the one here, but different enough I think to be worthwhile asking. I thought I'd put as a starter, what I think one of the hardest to grasp is. Mine is the difference between _probability_ and _frequency_. One is at the level of "knowledge of reality" (probability), while the other is at the level "reality itself" (frequency). This almost always makes me confused if I think about it too much. Edwin Jaynes Coined a term called the "mind projection fallacy" to describe getting these things mixed up. Any thoughts on any other tough concepts to grasp?


title: Correcting for multiple comparisons when running a bivariate correlation in SPSS

text: I am running a bivariate correlation analysis in SPSS, and I am performing multiple comparisons (there are 8 variables in total). I want to correct for multiple comparisons because I am aware that any 'significant' results could simply be flukes. However, the Bonferroni correction is not appropriate in this case (it is too strict). Does anyone know how to correct for multiple comparisons in SPSS? Consider some sample data as follows. There are 4 'independent variables' and 4 'dependent variables'. Independent variables:   1. Blood flow through middle cerebral artery   2. Blood flow through anterior cerebral artery   3. Blood flow through posterior cerebral artery   4. Blood flow through anterior communicating artery.  Dependent variables   1. Performance on cognitive test #1   2. Performance on cognitive test #2   3. Performance on cognitive test #3   4. Performance on cognitive test #4 The 4 'independent' variables are **not** uncorrelated to each other. The 4 'dependent' variables are also **not** going to be uncorrelated to each other (i.e., if a person does well on one test, chances are they will also do well on another test). (I realize that it is wrong to call these variables 'independent' and 'dependent', since correlation does not prove causality, but this is the way that I have framed them in my mind) I suppose that there are 2 problems here:   1. how to correct for multiple comparisons (from a statistics point of view) and   2. how to actually implement this in SPSS (a practical problem). Any help would be much appreciated (especially for problem #2).


title: MLE for bivariate distribution

text: I am trying to estimate parameters defining the rate of a Poisson bivariate distributions. (...following Maher's 1982 paper "Modelling association football scores"). Basically I have random variables $Z_{i,j} = X_{i,j}\cdot Y_{i,j}$ (here the $X_{i,j}$ and the $Y_{i,j}$ are independent and follow Poisson distributions with rates $\alpha_i\cdot \beta_j$ and $\alpha_j\cdot \beta_i\cdot\rho$ respectively, and $\alpha_i,\beta_j$ and $\rho$ are positive numbers for $1\le i,j \le N$. I have $M$ samples $(x,y)$. Depending on the case, $N$ is under one hundred, and $M$ in the thousands. So what is the best approach here to estimate the parameters? I haven't found a solver for the bivariate case. I'd need to optimize functions of the sort $$ {\mathcal L}(x_i,y_j; \alpha_i,\beta_j:1\le i,j\le N) =\sum_{k=1}^M -(\alpha_{i_k}\beta_{j_k}) + x_k\log(\alpha_{i_k}\beta_{j_k}) - (\alpha_{j_k}\beta_{i_k}\rho)+y_k\log(\alpha_{i_k}\beta_{j_k}) $$ where the index $k$ is an ordering of the sample and, e.g., the $k$-th sample is drawn from the distribution $P(\alpha_{i_k}\beta_{j_k})P(\alpha_{j_k}\beta_{i_k}\rho)$. Alternatively, I considered this as a univariate problem, first estimating the $\alpha$s, then the $\beta$s. I have successfully done this using scipy.optimize: basically   1. I start with the initial guesses $\hat\alpha_i=\frac{1}{\sqrt{\sum_{1\le i,j\le N,i\neq j} x_{i,j}}} \sum_{j\neq i}x_{i,j}$ and $\hat\beta_j=\frac{1}{\sqrt{\sum_{1\le i,j\le N,i\neq j} x_{i,j}}} \sum_{j\neq i}x_{i,j}$ suggested by Maher,   2. I estimate the $\alpha_i$s using Nelder-Mead optimization method from scipy.optimize with function $$ {\mathcal L}_\alpha(\alpha_i;x_k,\beta_j:1\le i,j\le N,1\le k\le M) =\sum_{k=1}^M -(\alpha_{i_k}\beta_{j_k}) + x_k\log(\alpha_{i_k}\beta_{j_k}) $$ I update the initial guess $\hat \alpha_i$s with the newfound estimation, and next I estimate the $\beta_j$s.   3. I loop over the above procedure.   4. I do the same with the $\gamma_i$ and $\delta_j$. As I write this, I am running experiments checking the prediction's accuracy from historical data. Results tend to be a little over %50.


title: Correction for non-linear trend seen in residuals plot when predictor is categorical

text: I'm running a linear regression analysis in R. One variable is a continuous outcome variable (`score2`) and the other is a categorical variable for treatment group (`group`). I also have a covariate which is the time 1 measure of the outcome variable (`score1`). So the model in R looks like this:               lm(score2~group + score1, data=dataSet)      Here's the question: When I plot the residuals, the plots shows evidence of non-linearity (a U-shaped plot). If both of my variables were continuous rather than categorical, I'd be sure that I need to correct for non-linearity (using Box-Cox for example). Does having a categorical predictor change this or do I still need to make the correction?


title: How to do a repeated measures multinomial logistic regression using SPSS?

text: ### Question:   * Can you do a repeated measures multinomial logistic regression using SPSS? ### Context: I need to do a regression on data at two points of time and I think this maybe the only way to go(?). To elaborate: I work for a national health service supporting individuals with psychosis. I want to investigate whether any factors (age, gender etc) predict the vocational outcome of a client group (3 categories: partial vocation, full vocation or no vocation) at entry of service and at 18 months into the service. I want to consider whether there is a difference at 18 months compared to entry to the service in regards to vocational status.


title: Quantile extrapolation?

text: Suppose you wanted to estimate the $q$ quantile of a distribution by observing $n$ independent draws from that distribution, but with $q < \frac{1}{n}$. What methods are available, and for what ranges of the 'undersampling coefficient' $q n$ are they expected to work? I can imagine some parametric approaches: a fully parametric approach would assume the form of the distribution, then fit _e.g._ the scale and location, and estimate the quantile from the fitted distribution. A slightly less parametric form would try to estimate the form ( _e.g._ fit to a power law) of the left tail, and use that to extrapolate. Are either of these sensible? Do they work? What are some other approaches?


title: AR(2) simulation problem

text: Take covariances $Cov[X_{t-2},X_{t}]$, $Cov[X_{t-1},X_{t}]$ and $Cov[X_t,X_t]=Var[X_t]$ and calculate the parameters for the AR(2) process ($a_1$, $a_2$ and $\sigma^2$ (the variance of the error term)) using the Yule- Walker equations. Once I have the parameters which define the AR(2) process I want to simulate a time series. For some covariances the parameter $a_2$ turned out to be very close to zero. The problematic cases arise when $Cov[X_{t-2},X_{t}]>Cov[X_{t-1},X_{t}]$. For example: $Cov[X_{t-2},X_{t}]=0.5$ $Cov[X_{t-1},X_{t}]=0.3$ $Var[X_t]=1$ results in the following AR(2) process: $$X_t=0.467X_{t-1}+0.0067X_{t-2} + \epsilon_t$$ with $\epsilon_t \sim N(0,0.864)$. Can this be regarded as an AR(2) process, or is it closer to an AR(1) process given the small parameter $a_2$? Can I pick any covariances and calculate parameters for the AR process, or are there constraints such as $Cov[X_{t-1},X_{t}]>Cov[X_{t-2},X_{t}]$? Thanks in advance. Cheers!


title: Attitude scaling

text: I wanted to know how to analyze the data to reflect the respondent's attitude. I 've got several questions or items about the trust in judicial system each of which has the responses 1. strongly agree 2. Agree 3. Not sure 4. Disagree 5. Strongly disagree. All the items have similar responses and is unidimensional to measure the attitude of trust towards judicial system. What is the best possible way to analyze it?


title: Data separability

text: I have a problem that can assume two "almost" equal inputs pointing to different targets (in two different situations). Like,               a) 141 442 53 15 151 023 => 25     ....     b) 152 442 53 15 151 023 => 74     ....      1) Is there a way to use a single (not different bases), to resolve this problem ? I thougth that if I use another input (that can differenciate the two sets) would be nice. But, it didn´t work well. Like:               a) X 141 442 53 15 151 023 => 25     ....     b) Y 152 442 53 15 151 023 => 74     ....      2) In this concept, can I use some kind of a Priority, that an algorithm can see some feature and goes that way ? 3) For example, I have 50 kind of "almost" equal data, that would take me to generate 50 different NNs to solve this. Thanks, Andre


title: Maximum likelihood estimation (MLE) for Markov Chain initial distribution?

text: I am working on using MLE to estimate a Markov Chain, I have successfully estimated the transition matrix $A$, using the method provided in http://www.stat.cmu.edu/~cshalizi/462/lectures/06/markov-mle.pdf Which is simply counting $A_{i,j} = N_{i,j}/\sum N_{i,j}$ But is there a way to estimate the initial distribution $\pi$ or named $\mu^{0}_{i}$s ? It seems like MLE can do it but I dont know how.


title: Which measure could tell me about what is the best predictor in survival analysis?

text: I have data of cancer survival, like that:               > my.data     survival stage my_class     27       3     2     221      2     1     43       3     3     ...      The survival is the time of survival since diagnosis in months (for example, in the first line the patient died 27 months before the diagnosis) and stage and my_class are two factors which I use to construct two different survival curves (considered to me as two models: one by stage, and other by my_class). I want to know which measure should I use to determine which is the best predictor the "stage" or "my_class". In R fitted the models, like that:               > fit.stage <- coxph(Surv(survival) ~ stage, data = my.data)     > fit.my_class <- coxph(Surv(survival) ~ my_class, data = my.data)      The summary of these objects include Rsquare and concordace, and I still can extract AIC (with extractAIC from the stats package) and AUC from SurvivalROC:               > extractAIC(fit.stage)          > extractAIC(fit.my_class)          > survivalROC(my.data$survival,         replicate(length(my.data$survival), 1),         as.numeric(my.data$stage), NULL, 251, method = "KM")          > survivalROC(my.data$survival,         replicate(length(my.data$survival), 1),         as.numeric(my.data$my_classifier), NULL, 251, method = "KM")      I am newer in the field survival analysis models evaluation, and I am in doubt about which of these measures I should use to select the best model between these two for example. Another predictos or combinations will be involved in the study too. Thanks a lot!


title: How do I test for a symmetric distribution?

text: I collect numbers from generators that yield different ranges of whole numbers with an unknown distribution. I want to estimate the mean of the numbers outputted by this generator. I'm convinced the distributions are symmetric (more specifically, uniform) though, so I can just average the min and max after I have a decent sample size.   What would be the best way to test that the distribution of numbers is symmetric or uniform? Things I've tried:   -Checked that the mean, median, and midpoint converge to the same values, but I can't quantify how close the values would have to be.   -Histogram, but my sample sizes (n<100 because I must manually collect) are too low for me to tell   -Chi-square test with the frequencies of each number in the range and testing them against the expected value for a uniform distribution, but I've read there are problems with this in other posts   -Pearson's skewness coefficient and comparing that with the standard error of skewness, but I think there are limits to this kind of approach -Comparing kurtosis of data with that of a uniform distribution I've also read about the Kolmogorov-Smirnov tests and Shapiro-Wilk tests but these seem too complicated for such a simple seeming task.. Thanks


title: Biased coin toss simulation -- which random generator is most appropriate?

text: I have a doubt whether the uniformly distributed random generator is the most appropriate to use for simulating biased coin toss. http://stackoverflow.com/questions/477237/how-do-i-simulate-flip-of-biased- coin-in-python Let's say a chance of Democrat winning in a state is 0.6. Ok, so I thought at first it makes sense to use uniform..., in infinite number of trials, in 60% occasions the value will be below 0.6 (democrat winning). However I was also wondering, should I use normal distribution instead? So let's assign the value (x) = 0 for democrat, and value 1 for republican. In the chart, we can draw a bar whose height is 0.6 on x(0), and another bar 0.4 on x(1). That makes the expected value (mu) = 0x0.6 + 0x0.4 = 0.4 The std deviation would be: square root of pow(0 - 0.4,2)*0.6 + pow(1-0.4,2)*0.4 = sort(24/100) = 0.49 Now, I execute random.gauss on python:               import random     random.gauss(0.4, 0.49)      I got: 0.9193140340408493 What do I make of that value? Should I interpret it as a win for Republican (under the rule: anything to the right of mean value goes to republican, the rest to democrat) ? If that's not true, please help me pointing out the flaw in my reasoning, and a pointer to a topic in stats book to fix that. Thanks in advance, Raka ADDITION: background: I want to implement spinner similar to the one I see in NYT site: http://www.nytimes.com/newsgraphics/2014/senate-model/ . The reason I was wondering if we should use normal random is because when I saw the source code of that page, it also uses normal random. Here's a snippet of the source code from that page (it's the click event handler; function named 'click'):               // Compute the national bias.       var nationalErrorScale = data.parameters.nationalErrorScale,           nationalBias = normalRandom() * nationalErrorScale,           localErrorScale = Math.sqrt(1 - nationalErrorScale * nationalErrorScale);            // For each race,       race.each(function(race) {              // Pick one of the potential matchups based on probabilities of each.         var matchup = race.matchups[bisectCumulative(race.matchups, Math.random())];              // Compute the local bias.         var localBias = normalRandom() * localErrorScale;              race.spinnerResult = matchup.mean + matchup.scale * (localBias + nationalBias) > 0 ? "dem" : "rep";       });


title: How to interpret overlapping interaction plots?

text: I have an interaction plot that seems to be overlapping on top of one another, what does this signify? ![enter image description here](http://i.stack.imgur.com/Lk7ak.png)


title: Long Run VAR alpha and beta significance levels

text: I am using a VAR with 2 variables and 4 lags. I am combining the coefficients of these variables to get an overall alpha and beta value for in the form $Y = \alpha + \beta X$. In order to get the long run alpha, I am using the following formula: $\alpha = \frac C {1-\alpha_1-\alpha_2-\alpha_3-\alpha_4}$, where the alphas are the coefficients for lagged $Y,\ C$ is the constant.![enter image description here](http://i.stack.imgur.com/CwDdL.png). and the following formula for the beta calculation: $\beta = \frac{\beta_1+\beta_2+\beta_3+\beta_4}{1-\alpha_1-\alpha_2-\alpha_3-\alpha_4}$, where $\beta$s are coefficients for lagged $X$ and $\alpha$ are coefficients for lagged $Y$. Am I correct in saying that these provide the LONG RUN alpha and beta values, or are these equations incorrect? **Also, and most importantly, what equation would I use to find the significance of these headline alpha and beta values? Are they derived from the standard errors/t-stats?** Attached is an example of the above, with the headline alpha and beta coefficients derived from the above formulas. Many thanks in advance.


title: How to compare the price of two sets of business transactions?

text: A business transaction can be defined as a function of 6 independent categorical variables (representing the customer and the product). Sales representatives are allowed to charge whatever price their customers are willing to pay. Alice and Bob are two sales representatives. They both conduct 100 business transactions in the same time period. For practical reasons, we cannot expect that the 100 transactions from Alice are exactly identical over all 6 variables from the 100 transactions of Bob.   * What is the statistical approach of determining whether Alice was able to charge a higher price than Bob?   * Is the number of observations needed to compare the two populations dependent on the number of variables?   * What is the statistical way of determining whether the "mix" of transactions from Alice are comparable to the ones of Bob?


title: How to guess a curve distribution from count data

text: I have a sample composed by 2500 count data values. I've plotted in R the corresponding histogram and ecdf. I've run the One-Sample Kolmogorov-Smirnov test to check if the distribution is either exponential or poisson, but I got negative results. Is it possible that the One-Sample Kolmogorov-Smirnov test, it's not the correct test to work with this count data. What is the best way to identify a distribution for this kind of data? ![ecdf and histogram](http://i.stack.imgur.com/DY1tJ.png)


title: modeling rates with machine learning tools (svm, gbm, nnet)

text: I have a numeric integer variable that is knowly proportional to an exposure measure plus other continuous / categorical covariates. If I were to use classical log-linear glms i would model ln(number)=offset(log(exposure))+explatatory factors. But, If I would like to use machine learning tools like gbm or SVM or Neural Network is it possible to handle this?


title: Interpreting results of a factor analysis

text: I performed factor analysis on R using factanal. Following advice I found on this tutorial, I chose the number of factors as being the number of principal components that capture 90% of the variability. I got the following table of results:                              Factor1 Factor2 Factor3 Factor4 Factor5 Factor6 Factor7 Factor8 Factor9     SS loadings     29.677   9.530   6.340   5.638   4.188   3.636   1.484   0.804   0.651     Proportion Var   0.309   0.099   0.066   0.059   0.044   0.038   0.015   0.008   0.007     Cumulative Var   0.309   0.408   0.474   0.533   0.577   0.615   0.630   0.638   0.645          Test of the hypothesis that 9 factors are sufficient.     The chi square statistic is 30148.02 on 3732 degrees of freedom.     The p-value is 0       I tried with one factor and I got a p-value of 0 as well...   1. Why is the cumulative var smaller than with 9 PCs for a pca?   2. How can I interpret the above table to choose the correct number of factors? (since the p-value is the same for 1 of 9 factors)


title: Determining Significance of Proportions

text: I've done a test where from a cohort of 112 people I've determined whether they respond or do not respond to the treatment. Now, the definition of them responding or not responding is whether they've lost a set value in that parameter. Now with that parameter, I've divided into groups to see at what value does a loss in that range have more responders than non-responders. However, how do I determine the significance of this? My table looks like this: ![Data requiring stats assistance](http://i.stack.imgur.com/Svxd4.jpg)


title: Validation: Data splitting into training vs. test datasets

text: I was naively validating my binomial logit models by testing on a test dataset. I had randomly divided the available data (~2000 rows) into training (~1500) and validation (~500) datasets. I now read a post in another thread ( Frank Harrell) that causes me to question my approach: > Data splitting is not very reliable unless you have more than 15,000 > observations. In other words, if you split the data again, accuracy indexes > will vary too much from what you obtained with the first split. How serious is this worry and what are ways around it? The OP speaks of "resampling" but not sure how that works here for validation. Edit: Adding context as per @Bernhard's comment below: Comparing logistic regression models


title: Moment of random variable on a integral form

text: Let $X_1,\dots,X_n$ are i.i.d. samples from uniform distribution on $(0,1)$. Let $\hat F_n$ be their _modified_ empirical distribution function defined by $$ \hat F_n(x)=\frac1{n+2}\left(1+\sum_{i=1}^n1_{\\{X_i\le x\\}}\right) $$ for every $x$, where $1_A$ is the characteristic function of the event $A$. Define $T_n$ by $$ T_n=\int_0^1\left\\{\Phi^{-1}(\hat F_n(u))-u\right\\}du $$ where $\Phi$ is the cummulative standard normal distribution. I wish to compute $E(T_n^k)$ for $k\ge 1$. For $k=1$, $E(T_n)$ can be computed as follows: Let $X_{(1)}\le \dots\le X_{(n)}\le X_{(n+1)}=1$ be the ordered statistics of our observation. Then \begin{align} E(T_n)&=\sum_{i=0}^n\int_0^1\left\\{\Phi^{-1}\left(\frac{i+1}{n+2}\right)-u\right\\}P\left(X_{(i)}\le u\le X_{(i+1)}\right)du\\\ &=\sum_{i=0}^n\int_0^1\left\\{\Phi^{-1}\left(\frac{i+1}{n+2}\right)-u\right\\}{n\choose i}u^i(1-u)^{n-i}du. \end{align} But then, how do I compute the second moment of $T_n$, i.e, $E(T_n^2)$? Could anyone help me? Any help will be vastly apprecaiated.


title: Expected prediction error - derivation

text: I am struggling to understand the derivation of the expected prediction error per below (ESL), especially on the derivation of 2.11 and 2.12 (conditioning, the step towards pointwise minimum). Any pointers or links much appreciated ![enter image description here](http://i.stack.imgur.com/rD05O.png)


title: Difference between linear regression prediction intervals and logistic regression targets

text: I am trying to understand the difference between logistic regression probabilities and linear regression prediction intervals. For example, let's say we have a database of student test scores in the range of 1 to 100 and some predictors. The goal of this study is to build a model to predict if other students will reach _at least a score of 60 with an 80% confidence_. To simplify, we are assuming that all the linear modeling requirements in the data are verified. The first method would be to run a linear regression on the observed data, then calculate the 80% prediction intervals and finally determine whether a student will reach a score of 60 or higher based on the lower end of the prediction interval. The other approach is to categorize our data and run a logistic regression on each student score < 60 (0) or >= 60 (1) observation. Is there any benefit in using the logistic regression approach in this case? Or does linear regression will result in the same level of accuracy when using prediction intervals?


title: Questions about Nolan's stable.exe program

text: I'm using J. Nolan's free stable.exe program to fit the pareto-levy stable model to some income data. This is option 3 "fit a bivariate sample" in the program. In the instruction text file it is said that after selecting option 3 and uploading the data (data = a text file containing a single column of numbers-- the income levels at each centile): "Some descriptive statistics will be displayed and you will be prompted for what kind of fit you want. The options are maximum likelihood fit, quantile fit of McCulloch, sample characteristic method of Koutrovelis-Kogon-Williams, user supplied parameters, and normal fit." The instructions also mention that I will be prompted for a number of other things, like the search method and the bounds of the four parameters. None of this happens when I do it. After uploading my data What I get is a line saying the number of values read from the data file ("14 data vectors read" for eg.), and then a line that says "Projection direction: 0-360 degrees" My first question is: What the heck does "Projection direction: 0-360 degrees" mean? No mention of this in the instructions. Secondly, I enter an angle, say 90 degrees, and then get a table of "Projection data" and values for the four parameters (alpha, beta, gamma, delta) in three categories: S0 parameters, S parameters, and S* parameters. Are these the fitted values of the parameters? And what do these three different categories mean? There are some other things bewildering me but I'll start with this. Thanks


title: What's a good approach to teaching R in a computer lab?

text: There have been several good questions and sets of answers on introductory books or approaches to learning R eg here and here. But I have a slightly different problem - the best way to run an hour long session (or several such sessions) in a computer lab that will get people started in R, familiar with its basic approach etc. My current plan would be to effectively work through the introductory chapter/s of something like Verzani's SimpleR and then introduce a familiar dataset, but is there any other approach that people have found useful? For example, is it good to introduce real data straight away, or address issues in more abstract way? Should I exhaustively go through how to use the square brackets, or excite people with examples of lattice graphics? My target audience are familiar with statistics (although not experts) and competent SPSS users; not familiar with programming languages beyond the sort of macro and scripting you'd get in SPSS and similar things. Any tips or references to lessons plans would be appreciated. However, I don't want to duplicate the many good lists of on-line material introducing R - strictly references to the face to face instructional question.


title: Baseline Classifier Evaluation

text: if i want to learn the baseline score of a accuracy precission recall and f-measure of a classifier...which classifies instances into 2 classes Class 1 has 100 instances Class 2 has 200 instances the baseline accuracy and precision is 200/(100 +200)= 0.66 right?? the baseline Recall is 200/(200+0) = 100 % right; the baseline F-measure 2* [(0.66*1)/(1+0.66) ] thank u


title: Confidence intervals for transition probabilities between known discrete states

text: I am trying to calculate confidence intervals for transition probabilities between a set of discrete states over time (involving multiple time steps). Here is a toy example to get the transition probabilities:               # create a toy dataset     states<- c("A","B","C")     t1 <- sample(states,30,c(2/3,1/6,1/6),replace=T)     t2 <- sample(states,30,c(2/3,1/6,1/6),replace=T)     t3 <- sample(states,30,c(2/3,1/6,1/6),replace=T)     df.t <- data.frame(t1=t1,t2=t2,t3=t3)              transition.matrix.t1 <- table(t1,t2)     transition.matrix.t2 <- table(t2,t3)              prob.trans.t1 <- t(apply(transition.matrix.t1,1,function(x) x/sum(x)))     prob.trans.t2 <- t(apply(transition.matrix.t2,1,function(x) x/sum(x)))      Each row of _df.t_ is a sample observation over 3 time periods. _prob.trans.t1_ contains the transition probabilities for the time step t1->t2. In reality, I don't have any previous knowledge about the underlying probabilities and I only have a randomly picked sample. I would like to calculate a 95% confidence interval for the above transition probabilities _transition.matrix.t1_ , _transition.matrix.t2_.


title: binarization of variable - experimental threshold choice. Is it good approach?

text: I have some ratings averages values from 1 to 5(users were rating on 1,2,3,4,5 scale). I would like to split them into two classes: `credible`, `not- credible`. I know that respondents have a bias toward evaluating objects as credible to willingly.               distribution:     1   1%     2   2%     3   5%     4   77%     5   14%      I will then use this binary variable **as explained variable in model** , where I use about 20 features to train a classifier to predict whether an object is credible or not. **QUESTION** : Is it acceptable to experimentally choose a threshold by comparing how well the classifier(e.g k-NN) performs (F-measure) with different cut offs? Or is such an approach a data-leakage crime?:) HOW DATA LOOKS LIKE: (sorted) ![enter image description here](http://i.stack.imgur.com/bV2hD.png) (unsorted) ![enter image description here](http://i.stack.imgur.com/yMZe2.png) I also consider aggregation with use of min: (sorted min) ![enter image description here](http://i.stack.imgur.com/HMWED.png) (unsorted min) ![enter image description here](http://i.stack.imgur.com/O6I7J.png)


title: Why is the MLE of N of the discrete uniform distribution the value you choose?

text: If $X$ is taken from a discrete uniform distribution, {1,2,.....,N}, why is the MLE of N=X? I can give an intuitive argument in that since you are only viewing one observation, it is the largest value and hence is the maximum. But, is there a nice mathematical derivation via argmax or anything like that? Thanks!


title: Coefficients by group

text: I have a data set with 8 different treatments and there is unequal number of observation within each group. I'd like to calculate regression coefficients for each group but I can not do it in other way than by subsetting data. I assume that I need some loop to run that but I have no idea how to do that. Can you help me somehow? Best regards Mateusz


title: Why does eigenvalue decomposition of a correlation matrix maximizes possible variance?

text: I was reading up on principal component analysis and I was wondering how does eigenvalue decomposition of the correlation matrix maxmimizes the possible variance that is captured? Can someone refer to a proof or give a simple intuitive proof? Thanks!


title: How to combine effect size and uncertainty in one measure?

text: The p-value and posterior probability measure the level of certainty in a particular outcome such as "the web traffic to variant A is greater than that to variant B." Of course these measures tell you nothing about the effect size, such as, "how much bigger is the web traffic to A than B?" Ideally you should measure both the size of the effect and the certainty. Suppose I want to compute a score for each of n variants of a website that incorporates information about both the size of the differences among the variants, and the uncertainty thereof. I have considered simply multiplying the posterior probability of an effect of interest by the expectation of the effect. Is this valid? Are there better ways to construct a score that incorporates both certainty and effect size? Thank you.


title: Estimation of mean of ratio of empirical distributions

text: I am trying to estimate the mean of the ratio of two empirical count distributions $f(x)$ and $b(x)$ defined on the limited range $[-1, 1]$ $$ \left<\frac{f(x)}{b(x)}\right> $$ $f(x)$ and $b(x)$ are measured in a counting experiment by filling unweighted histograms which then have Poisson uncertainties in each bin. To naively calculate the mean I can choose a binning for the histograms which is coarse enough that bins have enough counts to estimate the ratio and its uncertainty (e.g. more than zero counts in $b$) and calculate the ratio in each bin in $x$ and then average over the full $x$ range in the resulting ratio histogram. At the same time using a coarser binning means that I lose resolution when extracting the mean. Are there other possible ways to estimate this mean which allow me to avoid rebinning just to keep $b$ from ever being 0? Computing power isn't an issue so iterative approaches are perfectly fine.


title: Variances not homogeneous problem in ANOVA mixed between-within subjects, too bad?

text: I'm conducting a 2x3 mixed between-within subjects ANOVA (group x time). I have 20 participants, 10 for each group, who undergo all the 3 measures over time. Levene's test of equality of variances is not significant for both time1 and time3. Unfortunately, on time 2, group1 (M = 2,869, SD = 1,571) and group2 (M = 6,706, SD = 3,521) have variances significantly not homogeneous (p = 0.016). Is this too bad? Can I trust the ANOVA results? Thank you!


title: Which statistical method of anlysis should I use?

text: I am conducting an experiment where I will measure the strength of a material used for roads. Ideally, I would like to use a large surface area to measure this but this takes too long and is really expensive. I therefore will scale it down to one area which is $50cm^2$ and one which is $12cm^2$. I then want to compare the data I get from both these samples and see how (if) they differ. For the smaller sample, I'll also be doing less observations, i.e I'll get around 30 values for the $50cm^2$ one and only like 3 or 4 for the $12cm^2$ due to how much it costs. I was going to use the Normal Distribution, but as I won't have a lot of data, I can't justify using the Normal Distribution through the Central Limit Theorem. Would I be able to use a non-parametric test like the Sign test? Or, would it actually be something like the Wilcoxon signed rank test as the data is paired isn't it? How would I compare the data from there? Thank you


title: Proper way to deal with several sources of variability

text: What is the proper way of combining various sources of variability? For example, if we know that under the normal conditions, a machine produces parts with normally distributed diameter with some $\mu$ and $\sigma$, and that measuring this diameter has its own normally distributed "noise" with $\sigma_{meas}$. We than sample $n$ parts and measure each exactly once. What would be the mean and standard deviation of the measurements under the assumption that the machinery works fine? I am solving a homework assignment. I have stripped down the assignment to what I think is its pivotal point, from which I will be able to work my way out.


title: Estimating ROC/AUC on large data sets?

text: Plotting an ROC curve of a classifier compared to cases requires that the data set be sorted first on the classifier score. I am in a position where I need to calculate ROC on a large data set very quickly and sorting is the bottleneck (even using quicksort in C or F90). If instead of calculating ROC by thresholding at each case in the data set I instead threshold at every 100 cases then my execution time decreases by orders of magnitude based upon how I can write the code. The result is an ROC curve with let's say 10,000 points instead of 1,000,000. My tests show that the area under these two curves are the same out to > 5 decimal places. I would like to use this method but have not ran into anyone trying to speed up calculation in this way. Most of the lit. is on uses of ROC analysis where the data sets are relatively small and execution time is not an issue, so I have not found anyone else using this method or another to speed calculation by "thinning" out the points on the curve. Has anyone ran into a reference/study that has used or evaluated this or another method for speeding up ROC analysis? If so, or if you have other thoughts, please share.


title: Explanatory power of variables in multiple regression

text: Suppose that I have a regression model with two explanatory variables: $Y = b_0 + b_1X_1 + b_2X_2 + \epsilon$ The determination coefficient $R^2$ measures the amount of variability in the response variable $Y$ that can be explained by the explanatory variable. However, how can I find out the explanatory power that is exclusively attributable to each of the two explanatory variables as well as their common explanatory power? I can think of the following approach. Orthogonalize $X_1$ relative to $X_2$ and call it $X_1^*$. In the following new model $Y = b_0 + b_1X_1^* + b_2X_2 + \epsilon$ I obtain the determination coefficient $R^2$ for $X_1^*$, which would be the unique explanatory power for $X_1$. Do the same for $X_2$. Is this a reasonable approach? If yes, how can I obtain the common explanatory power for both $X_1$ and $X_2$? If not, do you have any solutions?


title: Sample size required for bivariate normal distribution

text: I am new here. Would like to ask a question on the sample size requirement for hypothesis testing following. If I am drawing samples with 2 non correlated variables (say x,y) from a bivariate normal distribution(assumed) with unknown mean and variance. I would need to determine the CEP or 50% of the population would fall into X < x and Y < y. How much sample size should I be taking? **What I have thought of:** I have tried to think of it as a binomial distribution problem as well (either sample is in or not in a certain radius r = sqrt(x^2 +y ^2) , and setting probability to be p =0.5 variance as p(1-p). Matching against the table for 50% confidence. I have calculated the sample size required to be ~480. While examining some literature, I have found sample size for Gaussian distribution to be only around 20-30. Am I making some wrong assumption? What is the more appropriate way to follow? THanks


title: Internal reliability for an ordinal scale

text: My teacher is asking whether it is possible to look at the Cronbachs's alpha when looking at the internal reliability of an **ordinal** scale. She thinks, because you use the mean of it, it is not possible, but I've seen it in previous research as well. I think it is possible, but I cannot give a technical explanation for this. Who can help me?


title: How to test for parameter stationarity?

text: Suppose I have time series observations from distributions drawn from some population. That is, I observe $X_{t,i}$ for $t=1,2,...,T,$ and $i=1,2,...,n$, where I believe that $X_{t,i}$ have pdf $f(\theta_i)$. (I have some idea about the distribution of the $\theta_i$, but that may not be important here.) I have some sample statistic which is a good estimator of $\theta_i$ given some observations. However, there is the suspicion that, in fact, the $\theta_i$ are not stationary, rather the observations come from $f(\theta_{t,i})$, where the $\theta_{t,i}$ are changing slowly over time. How can I test this, either by a formal hypothesis test or an 'eyeball' test? The amount of data available in the time domain is not so great ( _i.e._ $T$ is not so large), thus partitioning the time domain and computing the sample estimate on each partition would only be advisable for a small number (say 5) of partitions (because otherwise the standard error of the estimate is too great). However, the number of series, $n$, is largeish, say 10,000. I realize there are a number of gaps in this question, _e.g._ how the $\theta_{t,i}$ might be varying with time, the standard error of the parameter estimator, etc. However, any hints would be appreciated. To be concrete, one could think of the $X_{t,i}$ as being normally distributed with mean $\theta$ and standard deviation $1$, and the sample statistic is the sample mean.


title: How to do weighted pair hierarchical clustering in R?

text: Here is an example of hierarchical clustering of genes in the microarray data using the **weighted pair gene method** in `Spotfire`. I am not sure how to do this in `R`. In the `hclust` function, I see `ward", "single", "complete", "average", "mcquitty", "median" or "centroid"` as the methods. Also, lets say I have performed hierarchical clustering and found groups of genes using `cuttree` method. I wanted to plot the expression of genes in a group across columns (which may represent treatment, time, etc.). And I want to do this for all the groups separately. In a way similar to the Mfuzz package's way of showing clusters. Can any one please help me? TIA for any pointers.


title: Bayesian model comparison for randomly sampled sets of models

text: I have two models $M_1$ and $M_2$ that I am using to try and compare to observed data $D$. $M_1$ is an $n_1$-dimensional model, and $M_2$ is an $n_2$-dimensional problem. The Bayes factor $K$ to compare the models can be calculated using: $K = P(D|M_1)/P(D|M_2) $ assuming no prior preference for either model. The numerator and denominator can be written as $P(D|M_i) = \int P(D|\mathbf{w},M_i) P(\textbf{w}|M_i) d\mathbf{w}$ where $\mathbf{w}$ is the parameter vector, so the integral is over parameter space. Now say that due to e.g. computational constraints, one can only compute $M_1$ and $M_2$ for a finite number of random samples of the parameter vector $\mathbf{w}$, where the number of samples is given by $s_1$ and $s_2$. Would it be acceptable to then say that the integral above becomes a summation over the random samples, and assuming the random samples are uniformly distributed through parameter space, $P(\textbf{w}|M_i)$ becomes $1/s_i$, so that: $P(D|M_i) = \sum_{j=1}^{s_i} P(D|\mathbf{w_j},M_i) / s_i$ and so what is being compared in the Bayes factor $K$ is the ratio of the average probability over all the samples for each model?


title: Nonparametric Bayesian analysis in R

text: I am looking for a good tutorial on clustering data in `R` using hierarchical dirichlet process (HDP) (one of the recent and popular nonparametric Bayesian methods). There is `DPpackage` (IMHO, the most comprehensive of all the available ones) in `R` for nonparametric Bayesian analysis. But I am unable to understand the examples provided in `R News` or in the package reference manual well enough to code HDP. Any help or pointer is appreciated. A C++ implementation of HDP for topic modeling is available here (please look at the bottom for C++ code)


title: Sampling covariance matrix using Gibbs sampling

text: I am sampling covariance matrix from a Inverse Wishart distribution. In one dimensional case, after doing sufficient iterations I am taking the mode value for variance (after removing the burn-in values). How to do the same in a multivariate case?


title: How can I estimate theta for the inverse hyperbolic sine transformation?

text: I would like help with R code to estimate theta for the Inverse Hyperbolic Sine Transformation. This transformation is useful to transform skewed data that contain negative values or zeros. There are a couple of related posts that discuss the IHS and suggest there is a maximum likelihood approach to estimating theta but I can't work out how to apply it. These related posts are:   * How should I transform non-negative data including zeros?   * Inverse hyperbolic sine transformation: estimation of theta Please see example code below. I've made skewed data by using the inverse of the IHS. Now, given the skewed data, and no prior knowledge of theta, how can I work out what theta should be? I would be most grateful for R code to undertake this analysis.               # Define the IHS transformation and its inverse     IHS <- function(x, theta){  # Inverse IHS transformation       (1/theta)*asinh(theta * x)     }          Inv.IHS <- function(x, theta){  # IHS transformation       (1/theta)*sinh(theta * x)     }          set.seed(1)     # generate some normal data     x <- rnorm(1000)     hist(x, breaks='FD')          # skew it by applying the Inverse of the IHS transformation     xt <- Inv.IHS(x, theta=2)     hist(xt, breaks='FD') # yep this is skewed.  How could we estimate theta?      ![usng inverse of the IHS to skew data](http://i.stack.imgur.com/VusX2.png)


title: Ridging a singular covariance matrix (SEM)

text: Structural equation models sometimes cannot be fit due to a singular sample covariance matrix. Now some authors suggest to apply a "gentle ridging" to the diagonal which helps (this is for example done in Mplus).   1. How does ridging exactly work? and   2. If your preferred SEM software does not support ridging, how can you do it "by hand"?


title: Confidence interval for gam model

text: Reading mgcv::gam's help page: "confidence/credible intervals are readily available for any quantity predicted using a fitted model" However I can't figure a way to actually get one. I thought `predict.gam` would have a `type=confidence` and a `level` parameter but it doesn't. Can you help me on how to create it ?


title: Pairwise comparisons in accordance with regression fit

text: I have a unconventional question arising from my current master thesis on regression modeling. Suppose we have fitted a (linear) relationship between a dependent variable $y$ and an independent variable $x$. Now we choose two points on the $x$-axis, i.e. according to percentiles $x_{10}$ and $x_{90}$. These two points are chosen to select the data points for two groups in order to perform pairwise comparisons. In the first run, we choose 20 data pairs around the $x$-values of $x_{10}$ and $x_{90}$ and run a statistical test in order to infere if their $y$-values differ significantly. In a subsequent step, we choose, say 30 data points around $x_{90}$ and test against $x_{10}$; in the next step 40 data points around $x_{90}$ and test against $x_{10}$ and so on. The intention for this is to determine the group size and consequently the corresponding lowest $x$-value where the comparison turns out to be significant. How would you approach this problem. I thought of a many-to-one procedure like the Dunnett test where multiple groups are tested against the same control (which could be $x_{10}$ in our example). However, these multiple groups would contain partly the same subjects. Or is some sort of adaptive design the right choice. Or something completely different? I'd be very happy for any kind of help since I'm completely stuck with this question and have no idea how to solve it. * * * Edit: The idea behind this procedure would be to determine the lowest x-value where the pairwise comparison would turns out to be significant. Suppose the first and smallest group around $x_{90}$ is tested against the control $x_{10}$ and is not significant. So you take a look at the next group comprising the data points from group1 plus additional 10 and so on. If the pairwise comparison is significant, one determines the smallest x-value of the corresponding group as some sort of "change point". How could such a problem be approached? Thanks a lot! Andres


title: Quantify volatility on a time series data

text: I am now working on a set of time series data and I would like to quantify the series' volatility. I have think of an approach which is first model it as some auto-regressive model like ARCH and compute the sum of squared residuals between the model and the real data. Does it make any sense?


title: Sufficient, Complete Sufficient, UMVUE, Rao-Blackwell, Admissible. What are ties between these?

text: I am taking stat inference course. I have some trouble understanding some these terms:   1. Sufficient Statistics: a stat that does not depend on the parameter, say $\Sigma X$ for normal distribution    2. Complete Sufficient Statistics: do not understand this. Why complete? Is sufficient not enough? What do we use it for?    3. UMVUE: I understand this is the best unbiased estimator. It is the result of Rao-Blackwell Theorem?    4. Rao-Blackwell Theorem: we use this to show an unbiased estimator has the smallest variance, thus, the best.   5. Admissible: estimator $\delta(X)$ has a lower risk than any estimator. They both seem to imply best estimator. Could anyone tell me if I understood them correctly? Or what am I missing?   Is there any working examples that anyone could point to me, so I can better understand it. I have an exam coming up. Thanks a lot!


title: Dimensionality Reduction Algorithm for Large Dataset?

text: I have a reasonably large (5k variables x 120k cases) that I'd like to run a dimensionality reduction algorithm on. I tried doing a simple Factor Analysis on it in SPSS, but it (predictably) barfed on a 3GB machine. I got an answer by truncating the data set to 2.5k variables and 25k cases, but I started wondering if there is an algorithm other than Factor Analysis/PCA that would handle the data set better. My background in software development is much stronger than my background in statistics. It's my hope that there's an algorithm that could handle the original, unsampled data set more easily in 3GB of RAM. Does anyone know of such an algorithm? Forgive me if this is a FAQ; I tried finding the answer before I posted. **EDIT** : Ideally I'd like an existing implementation in something like R or SPSS, but since I'm a dev, even a software-based solution (like "try this feature of numpy") would be extremely helpful.


title: Zero-inflated two-part models for semi-continuous data

text: I am trying to study predictors of companies' pollution output of some specific chemicals. The data I am using have many 0's (i.e., the company did not pollute at all with those chemicals) and then are continuous with a long right tail. I have seen others model this data by logging the dependent variable after adding 1. My sense is that this is wrong, but I don't understand why. Could someone explain? This approach is much simpler than what I think I should be doing - using zero-inflated two-part models for semi- continuous data - so I'd be thrilled if it turned out simply adding 1 and logging is right. Second, I have found a Stata ado file to run zero-inflated two-part models for semi-continuous data. Is there a way to incorporate fixed effects into this type of model?


title: Trouble finding good model fit for count data with mixed effects - ZINB or something else?

text: I have a very small data set on solitary bee abundance that I am having trouble analysing. It’s count data, and almost all the counts are in one treatment with most of the zeroes in the other treatment. There are also a couple of very high values (one each at two of the six sites), so the distribution of the counts has an extremely long tail. I’m working in R. I have used two different packages: lme4 and glmmADMB. Poisson mixed models didn’t fit: models were very overdispersed when random effects were not fitted (glm model), and underdispersed when random effects were fitted (glmer model). I don't understand why this is. The experimental design calls for nested random effects so I need to include them. A Poisson lognormal error distribution did not improve the fit. I tried negative binomial error distribution using glmer.nb and couldn’t get it to fit – iteration limit reached, even when changed the tolerance using glmerControl(tolPwrss=1e-3). Because a lot of the zeroes will be due to the fact that I simply didn’t see the bees (they are often tiny black things), I next tried a zero-inflated model. The ZIP didn’t fit well. The ZINB was the best model fit so far, but I am still not too happy with the model fit. I am at a loss as to what to try next. I did try a hurdle model but couldn’t fit a truncated distribution to the non-zero outcomes– I think because so many of the zeroes are in the control treatment (the error message was “Error in model.frame.default(formula = s.bee ~ tmt + lu + : variable lengths differ (found for 'treatment')”). In addition, I think that the interaction I have included is doing something strange to my data as the coefficients are unrealistically small – although the model containing the interaction was best when I compared models using AICctab in package bbmle. I am including some R script that will pretty much reproduce my data set. Variables are as follows: d=Julian date, df=Julian date (as factor) , d.sq=df squared (number of bees increases then falls throughout the summer), st=site, s.bee=count of bees, tmt=treatment, lu=type of land use, hab=percentage of semi natural habitat in surrounding landscape, ba=boundary area round fields. Any suggestions as to how I can obtain a good model fit (alternative error distributions, different types of model etc) would be very gratefully received! Thank you.               d <- c(80,  80,  121, 121, 180, 180, 86,  86,  116, 116, 144, 144, 74,  74, 143, 143, 163, 163, 71, 71,106, 106, 135, 135, 162, 162, 185, 185, 83,  83,  111, 111, 133, 133, 175, 175, 85,  85,  112, 112,137, 137, 168, 168, 186, 186, 64,  64,  95,  95,  127, 127, 156, 156, 175, 175, 91,  91, 119, 119,120, 120, 148, 148, 56, 56)     df <- as.factor(d)     d.sq <- d^2     st <- factor(rep(c("A", "B", "C", "D", "E", "F"), c(6,12,18,10,14,6)))     s.bee <- c(1,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,4,0,0,0,0,1,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,3,0,0,0,0,5,0,0,2,0,50,0,10,0,4,0,47,3)     tmt <- factor(c("AF","C","C","AF","AF","C","AF","C","AF","C","C","AF","AF","C","AF","C","AF","C","AF","C",     "C","AF","AF","C","AF","C","C","AF","AF","C","AF","C","AF","C","AF","C","AF","C","AF","C",     "C","AF","AF","C","AF","C","AF","C","AF","C","C","AF","C","AF","C","AF","AF","C","AF","C",     "AF","C","AF","C","AF","C"))     lu <- factor(rep(c("p","a","p","a","p"), c(6,12,28,14,6)))     hab <- rep(c(13,14,13,14,3,4,3,4,3,4,3,4,3,4,15,35,37,35,37,35,37,35,37,0,2,1,2,1,2,1),              c(1,2,2,1,1,1,1,2,2,1,1,1,1,1,18,1,1,1,2,2,1,1,1,14,1,1,1,1,1,1))     ba <-  c(480,6520,6520,480,480,6520,855,1603,855,1603,1603,855,855,12526,855,5100,855,5100,2670,7679,7679,2670,     2670,7679,2670,7679,7679,2670,2670,7679,2670,7679,2670,7679,2670,7679,1595,3000,1595,3000,3000,1595,1595,3000,1595     ,3000,4860,5460,4860,5460,5460,4860,5460,4860,5460,4860,4840,5460,4840,5460,3000,1410,3000,1410,3000,1410)     data <- data.frame(st,df,d.sq,tmt,lu,hab,ba,s.bee)     with(data, table(s.bee, tmt) )          # below is a much abbreviated summary of attempted models:          library(MASS)     library(lme4)     library(glmmADMB)     library(coefplot2)          ###     ### POISSON MIXED MODEL              m1 <- glmer(s.bee ~ tmt + lu + hab + (1|st/df), family=poisson)         summary(m1)              resdev<-sum(resid(m1)^2)         mdf<-length(fixef(m1))          rdf<-nrow(data)-mdf         resdev/rdf     # 0.2439303     # underdispersed. ???                    ###     ### NEGATIVE BINOMIAL MIXED MODEL              m2 <- glmer.nb(s.bee ~ tmt + lu + hab + d.sq + (1|st/df))     # iteration limit reached. Can't make a model work.                    ###     ### ZERO-INFLATED POISSON MIXED MODEL              fit_zipoiss <- glmmadmb(s.bee~tmt + lu + hab + ba + d.sq +                         tmt:lu +                         (1|st/df), data=data,                         zeroInflation=TRUE,                         family="poisson")     # has to have lots of variables to fit     # anyway Poisson is not a good fit                    ###     ### ZERO-INFLATED NEGATIVE BINOMIAL MIXED MODELS          ## BEST FITTING MODEL SO FAR:              fit_zinb <- glmmadmb(s.bee~tmt + lu + hab +                         tmt:lu +                         (1|st/df),data=data,                         zeroInflation=TRUE,                         family="nbinom")         summary(fit_zinb)     # coefficients are tiny, something odd going on with the interaction term     # but this was best model in AICctab comparison          # model check plots         qqnorm(resid(fit_zinb))         qqline(resid(fit_zinb))              coefplot2(fit_zinb)              resid_zinb <- resid(fit_zinb , type = "pearson")         hist(resid_zinb)              fitted_zinb <- fitted (fit_zinb)         plot(resid_zinb ~ fitted_zinb)               ## MODEL WITHOUT INTERACTION TERM - the coefficients are more realistic:              fit_zinb2 <- glmmadmb(s.bee~tmt + lu + hab +                         (1|st/df),data=data,                         zeroInflation=TRUE,                         family="nbinom")          # model check plots         qqnorm(resid(fit_zinb2))         qqline(resid(fit_zinb2))              coefplot2(fit_zinb2)              resid_zinb2 <- resid(fit_zinb2 , type = "pearson")         hist(resid_zinb2)              fitted_zinb2 <- fitted (fit_zinb2)         plot(resid_zinb2 ~ fitted_zinb2)                    # ZINB models are best so far     # but I'm not happy with the model check plots


title: Finding similar users

text: I am working on a problem in the online advertising space. I am trying to identify consumers similar to the set of consumers who have bought a product in the past (have 'converted'). If I can identify 'similar' users, I can deliver them targeted messaging with the goal of turning them into paying customers. I have approached this as a predictive modeling problem. I have been using KNN and logistic regression to determine major factors that influence the target variable ('being a converter'). I feel that I am approaching this wrong, as in my set of 'non-converters' I have a set of people who are quite similar in every way to 'converters'. So if I built a perfect model it would not be helpful in detecting 'non-converters' who should be 'converters'. My ultimate goal is to detect the non-converters who behave similarly to converters. Should my goal with the model be to predict converters with high accuracy, and argue that users classified as converters who were actually non- converters fall into the set of consumers I should be targeting? Can someone help push me in the right direction?


title: Adjusting time series for methodological change during the time period

text: The issue is that the national bureau of statistics changed permanently the variable calculation methodology during the overall period. Now I have two official time series for the same variable. One of them comprises the period from 1985 to 2004 using the old methodological approach, and the other one comprises the period from 1995 to 2011 using the new methodological approach. This troubles immediate comparisons between those two periods, so I'd Like to know if someone is aware of any well supported statistical method that: Using the data for each series in the 'time intersection period' (1995-2004) as _a priori_ information, allows me to calibrate the 1995-2011 and backdate (or ascribe values) to the 1985-1995 period, so I could compare if the original 1985-1995 series behaves like the backdated 1985-1995 series. References would be welcome. Best regards.


title: What should this poor guy be told?

text: > **Possible Duplicate:**   >  What are the chances my wife has lupus? A day or two ago a man posted a question on math.stackexchenge about whether his wife likely has lupus, given the diagnosis. First, http://math.stackexchange.com/questions/256079/what-are-the-chances-my-wife- has-lupus It appears I will be able to post a screenshot of the question and of my answer, which i deleted after the OP expressed a lack of satisfaction. However, the screenshots take up quite a lot of room, so I will try to make a question here. On the odd chance that this man ever looks at his MSE question again, is there anything else of value he could be told? Note that there are two answers after mine, explaining rather better than i did that he has not given enough information. =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= ![enter image description here](http://i.stack.imgur.com/ePN2r.jpg) =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= ![enter image description here](http://i.stack.imgur.com/D14gk.jpg) =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=


title: Independence of errors assumption in paired & independent samples t-test

text: first of all, this is a great site and I've already learned a lot from threads here wile searching for answers to my stats questions on the web. Yet here's another one: I understand that the independence of errors assumption in the t-test means that you shouldn't be able to predict the value of the outcome variable in any observation from the values in the other observations.   1. This refers to observations regardless of group membership (as in group 1 and group 2 of the t-test), right? That is, you should neither be able to predict the outcome variable based on other observations of the outcome variable WITHIN one group, nor BETWEEN the groups?   2. And is paired samples t-test necessary because of the specific violation of this assumption in the sense that you could predict observations in one group from those in the other?   3. And finally, then, the paired t-test should still require independence of observations within the 2 groups, right? If I understand correctly, the first answer in the thread below basically confirms my presumptions, but since I'm not quite sure, I figured I'd ask about it explicitly. Question about independence assumption for ANOVA, t-test, and non-parametric tests Thanks and cheers, Hans


title: What type of non-linear regression would you use on this data

text: I am trying to model a wind turbine power curve, or in other words, given the windspeed (X), I would like predict the produced power (Y). I want to use historical data to estimate this model. You can read more about the power curve here: http://www.windpowerwiki.dk/index.php?title=Power_curves. On the given rescaled power curve below, the power can never go below 0 and never above 4, and the windspeed can never go below 0. **Question 1** : What type of model would you use to predict Y with X (the red line on the plot is a LOESS fit)? **Question 2** : Is it possible, in your suggested model, to use lagged variables of Y? The dots on the graph are measurements over time, so if one point is above the red line, it is very likely that the next measurement in time is also above the red line (this is because another important predictor, air pressure, is not included in the model, so I assume there is important autocorrelation). How can I take this into account? Thank you! ![enter image description here](http://i.stack.imgur.com/JegJv.jpg)


title: Deducing a single variance from a covariance matrix

text: I use linear discriminant analysis and I want to compare different assumptions about the covariance matrix. I want to explore what happens if I use a covariance matrix proportional to the identity matrix although my data does not come from such a covariance matrix. I get the following sample covariance matrix: $$ \begin{bmatrix} \sigma_1^2 & b \\\ b & \sigma_2^2 \\\ \end{bmatrix} $$ However, I do not want to use this matrix. I want to use something like: $$ \begin{bmatrix} \sigma^2 & 0 \\\ 0 & \sigma^2 \\\ \end{bmatrix} $$ So how can I estimate $\sigma^2$ from the first matrix. Is this correct? $$ \sigma^2 = (\sigma_1^2 + \sigma_2^2)/2 $$


title: Logistic regression in R resulted in Hauck Donner phenomenon. Now what?

text: I'm trying to predict a binary outcome using 50 continuous explanatory variables (the range of most of the variables is -∞ to ∞). My data set has almost 24,000 rows. When I run glm in R, I get: Warning messages:   1: glm.fit: algorithm did not converge   2: glm.fit: fitted probabilities numerically 0 or 1 occurred I've read the other responses that suggest perfect separation might be occurring, but I'm confident that isn't the case in my data (though quasi- complete separation could exist; how can I test to see if that's the case?). If I remove some variables, the "did not converge" error might go away. But that's not always what happens. I tried using the same variables in a bayesglm function and got the same errors. What steps would you take to figure out exactly what's going on here? How do you figure out which variables are causing the problems?


title: Poisson fixed and random effect

text: I have claim frequency panel data for five years consisting of age, CC, make of car, age and gender of insured person and geographical area. I have run a Poisson regression to model the claim frequency in R. Now I want to test for a Poisson fixed and random effect. How can I proceed?


title: Using Student or Wilcoxon test on small sample without information on the normality of distribution

text: I have three subjects in a qualitative study. I am comparing the features of their writing before and after treatment. I have no idea about the normality of the distribution since they were selected purposefully. I have counted the number of everything (like the complex sentences, correct sentences, etc.) and now I want to compare each case in the subjects before and after treatment. When you look at the numbers, the difference is too high (like for example around 800 before the treatment but more than 2000 after). T-test results showed a p-value of 0.017 while wilcoxon showed 0.109. Hence t-test showed a significant difference while wilcoxon showed a nonsignificant difference. Now my question is, due to the small sample size (only three) and not being sure of the normality of distribution, is t-test a correct procedure? I have also run chi square and it has shown a significant difference.


title: At what point does LDA (Latent Dirichlet Allocation) not make sense to use?

text: I was wondering if anyone could provide some intuition about when the "documents" in LDA are too small for it to provide any benefit. For instance, I have seen papers where people have you used context windows of 5 words as 'documents" which I find to not make too much sense as the model is only taking words into account (no structure). Clearly at some point clustering that takes more than word counts must be used instead of LDA. By this I mean taking word similarity into account or structural info. Is there any intuition about the nature of this relationship?


title: Effect of dismissing data points in regression

text: I would like to know what is the effect of dismissing data points in a regression analysis. For example, say originally a sample of 100 data points is collected (say they are all useful, equally reliable observations, no outliers) but one decides to fit a model to only 50 or 60 points from this sample. Let us say also, that the dismissed points are not chosen randomly. I am looking at a case where the parameters of a normal distribution $\Phi$, $\sigma$ and $\mu$, are estimated using simple linear regression: $$ X=\sigma \Phi^{-1}(p(X)) + \mu $$ The points that are systematically dismissed in this example are the points which have probabilities $p(X)$ of 0 or 1, because the standard normal distribution does not converge there. From reading Green (1991), I understand that dismissing data points reduces the _power_ of the analysis. 1/ Firstly, in the context of regression analysis I am not 100% sure what this means (increases the probability of type II erros? On what?). 2/ Any other consequences of dismissing data points apart from this? Ref. Green, S.B. (1991). How many subjects does it take to do a regression analysis? Multivariate Behavioral Research, 26(3), 499-510.


title: Distance from bivariate Gaussian mean in terms of variance

text: Not sure if my question is a valid one but I will just put it out here. Consider a bivariate data set $(x_i, y_i)$ $[i=1,...,n]$ to which a bivariate Gaussian Distribution is fitted. Now, consider a data point $z=(x_p, y_p)$. I wish to express the Euclidean distance of $z$ from the bivariate Gaussian mean $(\mu_1, \mu_2)$ in terms of $\sigma$, the variance-covariance matrix here. Can this be done? My motivation for this question comes from the univariate case wherein the aforementioned distance can be obtained easily. Can this be done for the bivariate case? If yes, any assistance in terms of theory or R code would be helpful as I am clueless about both. Thanks!


title: Comparing two items with various parameters of different units

text: Not a big stats guy, but I dabble in programming and frequent stackoverflow and figured this would be a good place to turn. Here's what I'm looking for. I have a set of food items. Each has a number of different parameters (right word? I'm not sure) such as calories, grams of fat, mg of sodium etc; some have different units. I'd like to be able to rank these items in terms of their similarity in nutrient composition to a another food item. My naive approach was creating a similarity score for each item as the sum of the percent differences of each parameter, and then ordering that vector. Then I had an idea where for each parameter I'd rank the parameter in terms of its percent difference, generating a rank matrix. But, then my mind ran out of ink and I couldn't quite think of what to do with the rank matrix - create rank products? sums? I dunno. Again, not much stats background just spitballing in my head. If someone could just suggest a method I should look into, I'd really appreciate it.


title: Measuring treatment effect on top-ranked subjects selected at point in time from longitudinal data

text: I am trying to measure the effect of treatment after some time on a group of subjects who were selected for treatment because they were the highest (top 1200) ranked subjects during a single period marking the start of the experiment, the hope for the treatment being that they will have lower ranks in the future than they otherwise would. (Note: I did not design the experiment.) No controls were set aside from this selection, making this essentially an observational study rather than a randomized experiment for all intents and purposes. The only candidates for controls are those subjects who were not in the top 1200 that period, meaning all other subjects in previous time periods where no experiment was being conducted. However, there are obviously a few confounding factors:   * Because they were selected based on top rank, there are going to be regression-to-the-mean effects that must be accounted for;   * Any trend in the metric of interest might have a trend that would make matching on raw values problematic. The initial plan of action was as follows:   * All subjects of similar rank from prior time periods (who are also test subjects) are candidates for control. This is a minimal requirement to deal with regression effects.   * Any trend in the metric should be found and the series should be detrended, i.e. work with the residuals, rather than the raw data. Use only data from periods prior to the start of the experiment.   * Match test subjects to control candidates based on similarity over some selection of points that correspond to a "pre" treatment period, e.g. t-9, t-6, t-3, t0, using a caliper of, say, 0.25 SD. This is to establish that they are sufficiently similar prior to the start of treatment.   * Use a Wilcoxon test to compare treatment subjects with the averaged matched controls for each subject on some "post" period change, e.g. the change from t0 to t+6. Problems I'm running into are that the trend is not a simple linear trend, first of all. I'm also at a loss over which data should be included in my trend determination. Also, since the treatment group is expected to have some impact on the trend post-intervention, and since there are few examples of people of similar rank during that time frame, it's difficult to determine if the trend selected is holding well as it's projected forward. On clue to the fact that things aren't working out is that I tested the protocol using only older data and a dummy test group. This experiment should result in something not significantly different than zero, but I get results that are have a significant difference. This casts doubt on all my actual results. There also might be holes in the protocol I've described that need improving. I could be going about things all wrong and I'm beginning to think that's the case, that it's impossible to approach this way. Does anyone have any suggestions on how to better go about tackling this problem? I'm open to changing how I'm doing anything here, even totally different approaches to the problem. At the end of the day I need a statistically defensible method of saying there is or isn't a treatment effect and some way of talking about its magnitude. **To summarize: I have a longitudinal data set tracking an irregularly trending metric. At some point along the way, some number of subjects were selected for being the current leaders in the metric and treated. Then more data was collected on everyone for some more time. How can I best determine the presence and magnitude of any treatment effect (90-95% confidence interval with p-value ideally).** FWIW, I'm doing the analysis in R, if anyone wants to suggest helpful methods or packages. General advice is also welcome.


title: Mathematical Explanations behind ANOVA

text: I was wondering: could someone could direct me to papers (preferably) or textbooks that give strong explanations of the mathematical background behind ANOVA?


title: Does it make sense to plot train and test results on a common ROC graph?

text: I am building a model using SPSS with 80% of my data. I can plot a ROC curve for both the training set and separately for the test set. I can´t find a way to plot both on the same graph.   Since both sensitivity and sensibility vary between 0 and 1, this seems posible to do. My question is: does it make sense to directly compare both curves, in other words, is it statistically valid?


title: Crossed vs. nested & fixed vs. random factors

text: A barn has 4 sections of animals. Within each section are 4 goats. Each goat is given one of four types of food. One of each goat's kidneys is randomly selected to be inspected after 3 days for level of carbohydrates. To measure the level, the goats' kidneys are divided into three pieces, and three measurement techniques are used. I understand two treatments are applied: $\text{Food}$ (4 fixed levels) and $\text{Measurement Technique}$ (3 fixed levels). I am wondering a couple of things:   * Are the kidneys "factors"?      * If so, are they nested or fixed?  My original thinking is that they are fixed...We have a left and right kidney for each animal. Some of my classmates are arguing that the kidneys are nested because each kidney differs for every goat. I don't understand this logic, because EVERY experimental unit differs (even if they are homogeneous, they are different still). So, I'm thinking:   1. $\text{Kidney}$ is a crossed factor with 2 random levels.   2. $\frac 1 3$ of $\text{Kidney}$ is a split-plot experimental unit, that is, a factor that is crossed with $\text{Kidney}$ and has 3 fixed levels. My classmates:   1. $\text{Kidney}$ is nested (2 levels) with the goat.   2. $\frac 1 3$ of $\text{Kidney}$ is nested, with the $\text{Kidney}$ (3 levels). What are your thoughts?


title: ANOVA repeated measurements

text: I am a master's student, and I need your help. I did experiment on one group of rats, the group was tested during 8 days, this is the first stage. The same group was also tested again but this time the test phase was 7 days. I tried to do my analysis using ANOVA repeated measures but I had a problem with the days.. I want to know the difference in the 2 tests but the days are not the same number. I don't know how to handle it. If there is any way to do the analysis using repeated measures ANOVA. Thank you. I tried to use repeated measure anova within subject effect My first factor will be test>>> with 2 levels (1 & 2 test) second factor is day with the level of ?? Here is the problem! if I write 8 there will be extra day in test 2 and the analysis cant carry!! Please help me. Thank you.


title: Reasons not to use fuzzy numbers instead of pds to represent uncertainty

text: Can someone explain why (if at all) it would be a bad idea to use fuzzy numbers to represent uncertainties in model parameters instead of probability distributions. I'll try to explain the motivation behind my question. Assume a decision model where the parameters are provided as pds. The information to support the decision is computed by running a MC simulation to evaluate that evaluates the model at each step and the result is, of course, a probability distribution. From the resulting pd, one can compute the expected value, risk, and other statistical measures. Could one simply use fuzzy numbers in place of pds and apply fuzzy arithmetic instead of MC?


title: How to measure the consistency of clustering results

text: I'm clustering data on a daily basis and would like to measure the consistency of the clustering method. Let's say following clusters result in method A:               On day 1: {a,b,c} {d,f} {g}     On day 2: {a,b,c} {d,f} {g}     On day 3: {a,b,c} {d,f} {g}      With method B:               On day 1: {a,b,c} {d,f} {g}     On day 2: {a,b} {c,d,f} {g}     On day 3: {a,b} {c,d} {f} {g}      With method C:               On day 1: {a,b,c} {d,f} {g}     On day 2: {a} {b} {c,g} {f}     On day 3: {a,g,d} {b} {c} {e,f}      The amount of variables stays the same, but the cluster sizes and count varies. Obviously the grouping is less consistent in the latter examples than in the first one. Ideally I'd like to have a measure that assigns a value of 1.0 to a completely consistent method and 0.0 where the clustering seems random. I'm struggling to find any literature and pointers to how this can be achieved. How could it be done? ### Background In my case the clusters are on correlation matrices of financial instruments.


title: Hypothesis testing with very high z value

text: Bags of concrete mix labeled as containing 100 lb have a population mean weight of 100 lb and a population standard deviation of 2.551 lb. What is the probability that the mean weight of a random sample of 50 bags is less than 95 lb? In solving this, I tried to use the formula $z = \dfrac {\bar x - \mu} {\frac {\sigma} {\sqrt n}}$, but I end up getting $z = \dfrac {95 - 100} {\frac {2.551} {\sqrt 50}} = -13.859$. This seems like a ridiculously large z score, can it possibly be right? Or am I missing something?


title: How well can multiple regression really "control for" covariates?

text: We’re all familiar with observational studies that attempt to establish a causal link between a nonrandomized predictor X and an outcome by including every imaginable potential confounder in a multiple regression model. By thus “controlling for” all the confounders, the argument goes, we isolate the effect of the predictor of interest. I’m developing a growing discomfort with this idea, based mostly on off-hand remarks made by various professors of my statistics classes. They fall into a few main categories: **1\. You can only control for covariates that you think of and measure.**   This is obvious, but I wonder if it is actually the most pernicious and insurmountable of all. **2\. The approach has led to ugly mistakes in the past.** For example, Petitti & Freedman (2005) discuss how decades’ worth of statistically adjusted observational studies came to disastrously incorrect conclusions on the effect of hormone replacement therapy on heart disease risk. Later RCTs found nearly opposite effects. **3\. The predictor-outcome relationship can behave strangely when you control for covariates.** Yu-Kang Tu, Gunnell, & Gilthorpe (2008) discuss some different manifestations, including Lord’s Paradox, Simpson’s Paradox, and suppressor variables. **4\. It is difficult for a single model (multiple regression) to adequately adjust for covariates and simultaneously model the predictor-outcome relationship.** I’ve heard this given as a reason for the superiority of methods like propensity scores and stratification on confounders, but I'm not sure I really understand it. **5\. The ANCOVA model requires the covariate and predictor of interest to be independent.** Of course, we adjust for confounders precisely BECAUSE they're correlated with the predictor of interest, so, it seems, the model will be unsuccessful in the exact instances when we want it the most. The argument goes that adjustment is only appropriate for noise-reduction in randomized trials. Miller & Chapman, 2001 give a great review. So my questions are:   1. How serious are these problems and others I might not know of?   2. How afraid should I be when I see a study that "controls for everything"? (I hope this question isn't venturing too far into discussion territory and happily invite any suggestions for improving it.) _EDIT_ : I added point 5 after finding a new reference.


title: Contrasts in repeated mesures ANOVA in R

text: I am a refugee from SPSS in the process of re-learning how to do everything in R. Mostly it's been fun, as R is great for a lot of things but I've run into a snag that I can't seem to find a solution for; any help appreciated. I want to run repeated-measures ANOVA with contrasts. I know linear mixed models are generally the better way to go but I work in a community where that is only just starting to catch on and my students need to understand what they're reading when they encounter older techniques. SPSS makes contrasts on repeated measures easy, but R does weird things I don't understand. Here's a sample:               df.wide = data.frame(Subj = 1:8, Day1 = c(6,5,5,6,7,4,4,5),                           Day2 = c(5,5,6,5,3,2,4,7), Day3 = c(2,4,3,4,3,1,1,2))     attach(df.wide)          ## contrasts as t-tests     Linear = -Day1 + 0*Day2 + Day3     Day1vs23 = -2*Day1 + Day2 + Day3     t.test(Linear)     t.test(Day1vs23)          ## repeated-measures with long-format data     library(reshape2)     df.long = melt(df.wide, id.vars="Subj", variable.name="Trial", value.name="DV")     contrasts(df.long$Trial) = contr.poly(3)     df.long$Subj = factor(df.long$Subj)     Anv = aov(DV ~ Trial + Error(Subj/Trial), data=df.long)     summary(Anv, split=list(Trial=list("Linear"=1, "Quad"=2)))      The t-tests show t-values of 7.51 for Linear and 1.27 for quadratic, but this does NOT match the `aov()` output, which provides F-values of 25.28 and 2.51, and of course substantially different p-values. SPSS contrasts on the repeated measures yield t-values and F-values that match up. It looks like R isn't partitioning the error term like SPSS does for contrasts, and both contrasts are being tested against a 14 d.f. error term, which is why it doesn't match the t-test. That seems wrong to me. So, what am I doing wrong in R and how do I fix it?


title: What statistical methods are best suited for missing data problems?

text: The data I'm analyzing has been collected and measured from about 2'000 animals. For every animals, I have a number of _discrete_ values, like species, gender, colour etc., as well as some _continuous_ values, like age, size, weight, and some medicinal values like the level of some specific chemicals in the blood. I want to build a statistical model which can predict some missing values. For example if the data of an additional animal is missing the `weight` value, I want the model to predict it based on the other existing data. Coming from a software development background, I assume that **Neural Networks** would be one possible solution for the problem. But there must be others. **What's the name of this kind of problem? What statistical methods are used for it?**


title: Division by zero in the Neyman–Pearson lemma

text: I understand that in theory the Neyman-Pearson lemma gives the most powerful test for a prescribed level, when we have a simple null hypothesis against a simple alternative. However, when I try to apply it in the following situation, the problem is that I get a division by zero, so I'm doubly-stuck: (1) I don't know if the lemma applies at all in this case, I would say yes, as there are no such restrictions in the statement; (2) still I can't see how to formulate the test, i.e. how to define the rejection region. Here is my situation: _We draw n times from the uniform distribution on $(0,\theta)$._ $H_0$: $\theta=\theta_0$ $H_a$: $\theta=\theta_a$ _We need to find the most powerful test at level $\alpha$. (Suppose it's also known that $\theta_a <\theta_0$)_ If I apply the lemma directly, I see that for a given sample I can get that either the denominator or both likelihoods simultaneously become zero. Any ideas?


title: How do you find weight of edges between individuals based on co-participation in the same groups? (SNA)

text: I have a data set with 20,000 discussion forum threads, many only one or two posts, some up to 400-posts. I have 5,000 individuals who participated in these threads. I want to calculate the strength of the relationship between two people based on how many threads they have participated in for export to Gephi and clustering to see if I can find any clear groups - are there any pre-existing algorithms to do this? I hacked up my own function in R, which looks at the length of the thread (smaller thread means stronger connection to other participants - too large threads are excluded), number of times posting (two people posting twice each in a thread are more strongly related), etc. It kind of works, but I'm sure I'm reinventing the wheel, there must be some existing algorithms/approaches etc? (I'm using R, but even a general algorithm is welcome).


title: How to prove that this process is strict stationnary

text: $W_t$ is strict stationary process and $X_t=g(W_{t-1},W_{t-2},...,W_{t-q})$ where $g$ is a measurable function. Prove that the process $X_t$ is strictly stationary. Remark : to prove that a process $X_t$ is strictly stationary we have to prove that: $P(X_1 \leq C_1,X_2\leq C_2,...,X_n\leq C_n)= P(X_{1+h}\leq C_1,X_{2+h}\leq C_2,...,X_{n+h}\leq C_n)$ for any $C_1,C_2,..C_n$ and any $h$


title: Best test of change between aggregate scores for two sampling periods

text: I have a numerator and denominator for a population at two different measurement periods. The % score for each time period represents the percent of the population that got a specified service (e.g., percent of members in population that got a mammogram). The members in the population change between the two measurement periods, with some members dropping out and new members joining between the two measurement periods. What is the best test of statistically significant change between the two scores?


title: What does a "closed-form solution" mean?

text: I have come across the term "closed-form solution" quite often. What does a closed-form solution mean? How does one determine if a close-form solution exists for a given problem? Searching online, I found some information, but nothing in the context of developing a statistical or probabilistic model / solution. I understand regression very well, so if any one can explain the concept with reference to regression or model-fitting, it will be easy to consume. :)


title: Winbugs and other MCMC without information for prior distribution

text: What happens when you don't have an idea of the parameters distribution? What approach should we use? Most of the time we aim to undersatnd if a certain variable has any influence over the presence/absence of a certain species, and the variable is accepted or not according to the variable importance. This means that most of the times we are not thinking on the expetcted distribution a parameter should have. Is it correct to assume that all parameters follow a normal distribution, when all i know is that b1,b2,b3 and b4 should vary between -2 and 2, and b0 can vary between -5 and 5 ?               model {         # N observations         for (i in 1:N) {             species[i] ~ dbern(p[i])             logit(p[i]) <- b0 + b1*var1[i] + b2*var2[i] +                  b3*var3[i] + b4*var4[i]         }         # Priors         b0     ~ dnorm(0,10)         b1   ~ dnorm(0,10)         b2 ~ dnorm(0,10)         b3  ~ dnorm(0,10)         b4  ~ dnorm(0,10)     }


title: Which test to use in a single group pre- and post-treatment analysis

text: I am attempting to analyse pre-existing data (using SPSS). I have one group of participants, 4 quantitative, non-parametric outcome measures, and am looking to see if there is a change in outcomes post intervention. I also have qualitative outcomes, but am not worrying about them yet! I understand that t-tests (and the non-parametric equivalents) can be used for comparing several groups (i.e. if I had a control group as well), but don't know what to do with my single group of participants. Any help much appreciated!