title: Linear mixed model when interest is ONLY in population

text: Lets say we have 1 to 6 repeted measures over time for some (many) subjects, and the outcome is values between 0 and 23 (so we assume normal distribution). For a normal distribution the fixed effects has interpretation of population average in a mixed model. So we can either use a mixed model or a marginal model (gee). If interest is only in population we could stick with a marginal model, but is it better than a mixed model when we are only interested in population? I'm thinking that we would use up alot of unnecessary degree of freedom with a mixed model? What should I say to someone that wants to use a mixed model (normal) when only interested in population parameters?


title: What do the coefficients of the crossproduct of regression mean?

text: How can I interpret the coefficients of the crossproduct of each of the following codes? What do they mean? How can I deduce that they correspond to our expectation?               fit.model=lm(formula = CO ~ ., data = cigarettes)     Call:     lm(formula = CO ~ ., data = cigarettes)          Coefficients:     (Intercept)          Tar     Nicotine       Weight            3.2022       0.9626      -2.6317      -0.1305       X=model.matrix(fit.model)     > X                      (Intercept)  Tar Nicotine Weight     Alpine                     1 14.1     0.86 0.9853     Benson&Hedges              1 16.0     1.06 1.0938     BullDurham                 1 29.8     2.03 1.1650     CamelLights                1  8.0     0.67 0.9280     Carlton                    1  4.1     0.40 0.9462     Chesterfield               1 15.0     1.04 0.8885     GoldenLights               1  8.8     0.76 1.0267     Kent                       1 12.4     0.95 0.9225     Kool                       1 16.6     1.12 0.9372     L&M                        1 14.9     1.02 0.8858     LarkLights                 1 13.7     1.01 0.9643     Marlboro                   1 15.1     0.90 0.9316     Merit                      1  7.8     0.57 0.9705     MultiFilter                1 11.4     0.78 1.1240     NewportLights              1  9.0     0.74 0.8517     Now                        1  1.0     0.13 0.7851     OldGold                    1 17.0     1.26 0.9186     PallMallLight              1 12.8     1.08 1.0395     Raleigh                    1 15.8     0.96 0.9573     SalemUltra                 1  4.5     0.42 0.9106     Tareyton                   1 14.5     1.01 1.0070     True                       1  7.3     0.61 0.9806     ViceroyRichLight           1  8.6     0.69 0.9693     VirginiaSlims              1 15.2     1.02 0.9496     WinstonLights              1 12.0     0.82 1.1184     attr(,"assign")     [1] 0 1 2 3     (1)result=t(X) %*% X                      (Intercept)       Tar  Nicotine    Weight     (Intercept)     25.0000  305.4000  21.91000  24.25710     Tar            305.4000 4501.2000 314.67100 302.17874     Nicotine        21.9100  314.6710  22.21050  21.63176     Weight          24.2571  302.1787  21.63176  23.72096     (2)XbyX <- crossprod(X)                 (Intercept)       Tar  Nicotine    Weight     (Intercept)     25.0000  305.4000  21.91000  24.25710     Tar            305.4000 4501.2000 314.67100 302.17874     Nicotine        21.9100  314.6710  22.21050  21.63176     Weight          24.2571  302.1787  21.63176  23.72096      Also which crossproduct is correct? (1) or (2)? What's the difference? Here's the data:               > dput(cigarettes)     structure(list(Tar = c(14.1, 16, 29.8, 8, 4.1, 15, 8.8, 12.4,      16.6, 14.9, 13.7, 15.1, 7.8, 11.4, 9, 1, 17, 12.8, 15.8, 4.5,      14.5, 7.3, 8.6, 15.2, 12), Nicotine = c(0.86, 1.06, 2.03, 0.67,      0.4, 1.04, 0.76, 0.95, 1.12, 1.02, 1.01, 0.9, 0.57, 0.78, 0.74,      0.13, 1.26, 1.08, 0.96, 0.42, 1.01, 0.61, 0.69, 1.02, 0.82),          Weight = c(0.9853, 1.0938, 1.165, 0.928, 0.9462, 0.8885,          1.0267, 0.9225, 0.9372, 0.8858, 0.9643, 0.9316, 0.9705, 1.124,          0.8517, 0.7851, 0.9186, 1.0395, 0.9573, 0.9106, 1.007, 0.9806,          0.9693, 0.9496, 1.1184), CO = c(13.6, 16.6, 23.5, 10.2, 5.4,          15, 9, 12.3, 16.3, 15.4, 13, 14.4, 10, 10.2, 9.5, 1.5, 18.5,          12.6, 17.5, 4.9, 15.9, 8.5, 10.6, 13.9, 14.9)), .Names = c("Tar",      "Nicotine", "Weight", "CO"), class = "data.frame", row.names = c("Alpine",      "Benson&Hedges", "BullDurham", "CamelLights", "Carlton", "Chesterfield",      "GoldenLights", "Kent", "Kool", "L&M", "LarkLights", "Marlboro",      "Merit", "MultiFilter", "NewportLights", "Now", "OldGold", "PallMallLight",      "Raleigh", "SalemUltra", "Tareyton", "True", "ViceroyRichLight",      "VirginiaSlims", "WinstonLights"))


title: Binomial proportion confidence interval

text: I have a sample and I've classified each member into 3 groups A,B,C. To get uncertainties on the percentages of these groups, I've been told I should use the wilson score interval. As someone who is not well versed in this area, I'm not really sure how to implement this in my code! I've had a look at this, but is an acceptable method for my case? (I'm using python, so any hints in that area are very welcome!)


title: How to fit log-linear poisson autoregressive mixed model?

text: I have time-series count data $N_{i,j}$ (population sizes in site $i$ and year $j$) and I want to correlate year-to-year changes with the environmental conditions $x_{i,j}$. For this, I want to fit this model: $$\begin{eqnarray} \mbox{log} ( \mu_{i,j+1} ) &=& \mbox{log} ( \mu_{i,j} ) + \alpha + \beta x_{i,j} + \gamma_j \\\ \\\ N_{i,j} &\sim& \mbox{Poiss} ( \mu_{i,j} ) \\\ \gamma_{j} &\sim& \mbox{Norm} (0, \sigma ) \end{eqnarray} $$ So I'm interested in parameter $\beta$, the slope of the relationship. $\gamma_{j}$ is the random effect for year (as the residuals within single year were correlated). This is an autoregressive model. **How can I fit such a model?** I tried to look at the PESTS project (PESTS R code here) but I am not able to find out if and how to fit my model using it. (Note that I am trying to avoid Bayesian tools because of computation time - I have thousands of such models).


title: Linear regression - are my models nested or non-nested?

text: I’ve got two main groups, which I am comparing, a) rats and b) rabbits. Each group consists of two subgroups – breed 1 and breed 2. Using linear regression, I would need to find out, does a model taking into account the subgroups account for more variance than a model with only the main groups. So I have created two models: \begin{align} Y &= a_1+b1_1(\text{rab})+e_1 \tag{1} \\\ Y &= a_2+b1_2(\text{rat2})+b2_2(\text{rab1})+b3_2(\text{rab2})+e_2\ \tag{2} \end{align} Where all variables are dummy variables:   * $\text{rab}\ \ = 1$ for a rabbit and zero otherwise (e.g., rats)   * $\text{rab1} = 1$ for rabbit breed 1   * $\text{rab2} = 1$ for rabbit breed 2   * $\text{rat2}\, = 1$ for rat breed 2 In order to compare the models statistically, I would need to know, whether my models can be considered nested or not. If the models are considered nested, I would also need advice in how to compare them statistically - preferably with SPSS. With non-nested models I can handle the rest.


title: z-test in paper false? can Williams Hotelling test be used instead?

text: I am looking at the statistics of this paper: http://emotional.intelligence.uma.es/documentos/6-Extremera2011Emotional.pdf I am wondering how you can use the z-test to compare the two correlation coefficients? But here, we talk about a significant difference in correlation between X,Y and X,Z while controlling for 5 other variables at two different time points. Actually, you would need to use the Williams-Hotelling correlation coefficient comparison and not the z-test, right? However, W-H correlation significance testing does not control for any other variables, thus, i have no idea whether this is really the most appropriate thing to do? any comments on that would strongly be appreciated. Moimeme


title: Estimating the probability of causation based on finding a correlation, including experimental details

text: Say there is a hypothesis that A causes B (A -> B), and some likelihood that the hypothesis is correct (AB1%). Now, an experiment is run that claims to find a correlation between A and B. What I would like to know is whether my new probability (AB2%) is greater than, or less than, my previous probability (AB1%). Common sense tells me that AB2% **must** be greater than AB1% -- that there's no way that finding a correlation would make it less likely that A -> B. Is my common sense correct? I would like to account for the possibilities of sampling bias, huge experimental error, reporting bias, etc. -- everything except for falsified data. For example, knowing that things like these can happen:   * extremely biased experiment in which any correlation is much more likely to be due to biased sampling than to A -> B.   * very noisy experimental design, where random correlations between noise is much more likely than observing any true correlation   * totally random data, cherry-picked to find correlations how would I incorporate my knowledge of those (an any other) factors into an estimate for AB2%? * * * tl,dr; Is this statement true: a correlation is found, therefore it's more likely that there's a real causal link between two variables than it was before the correlation was found? * * * I'm not a statistician -- please be gentle! :)


title: Using Lorenz curve / Gini coefficient for (non-ecomoical) distribution data

text: I have distributional data which I represent as a density. The data represents frequencies of user activities on a computer screen (e.g. amount of clicks on the y or x-axis of that screen but also other activities that can be related to coordinates and can therefore be binned by those coordinates (e.g. 5 pixels bins)). I would like to compare two kinds of that behavior and find out how compatible their distributions are. Very general. No assumptions exist. I can't assume parametric conditions such as linearity or normality. I read about Lorenz curves and the Gini coefficient to be very much like what I need to compare distributions but also know that those methods find application primarily for economic and sociological problems and are usually not applied for general distributions. Am I applying the wrong tool for the job? What is your opinion about this? What alternatives do you recommend in order to find out how similar two distributions are?


title: Are there associated variants in these columns?

text: I have data that look like this in a 365*804 matrix:                Sample     Variant1      Variant2        Variant3      Person1      1/0            0/0             1/1      Person2      0/1            1/1             0/0      Person3      0/0            0/0             0/0      I want to see if any Variants are significantly associated with each other. By that I mean, if you see one how likely are you to see another? The possible values are either (0/0, 1/0, 0/1 or 1/1) with 1/0 and 0/1 being the same thing. If weighting is possible, 1/0 and 0/1 should be worth half what 1/1 is. I would prefer to do this in R.


title: Are log-log models the same as lognormal models?

text: I have a dataset that I want to fit according to $$\log(y) = a + b_1\log(x_1) + b_2\log(x_2) +\cdots + b_k\log(x_k).$$ My statistical package has options to do a linear regression and lognormal. I am not sure which one I should choose.


title: Can I do a t-test Power Analysis for Unequal Size Groups which Produces 2 Different Minimum n's?

text: It's usually straightforward to do a `Power Analysis` to compute the `minimum sample size`, especially in R which is my preferred statistical computing environment. However, I am being asked to conduct a Power Analysis that's a little different than anything I've done or which I can find reference to online. I'm wondering if what I'm being asked for is even possible/valid. The project basically has two `unequal groups` of states and the hypothesis is that these two groups are significantly different in terms of an outcome variable (which is the duration of phone calls to customers). The "control" group consists of 40 states and produced about 2,500 observations. The "test" group has about 10 states and 500 observations. Initially, I found group `means` \+ `pooled standard deviation`, which I used to calculate an `Effect Size`. Then I used a package called `pwr` in `R` and found that I needed a minimum sample size of about 135 observations per group, given .05 significance and .8 power. However, they rejected my answer because they want one group to be much bigger than the other like it is now, and they are expecting either two different minimum numbers of observations per group or a minimum % of the population in terms of numbers of states or observations that have to go into their "test" group. I see Power Analyses for two sample t-tests (i.e. the R function `pwr.t2n.test`), but I'd have to specify at least one of the sample sizes whereas they want me to tell them the minimal sample size for both groups (either as numbers or percentages) and this function doesn't reflect the differences in standard deviations for the two groups. Is this possible or do I just tell them that's not how it works (i.e. the best I can do is tell them that given one of the sample sizes and a pooled standard deviation the second group has to be at least a certain size)?


title: Need help selecting appropriate statistical tests

text: I have several data related questions I am hoping to get answers to. I have a 4-week intervention (no control group) and all participants are given and pre and post test assessment consisting of two highly reliable, valid, and correlated measures (BAI and BDI). Here are the questions I am trying to answer:   1. Is the observed difference in BAI / BDI scores between pre and post tests significant?   2. Is the observed difference in BAI / BDI scores between pre and post test due to chance or my 4-week intervention?    3. Does the number of sessions attended have an impact on the change between pre and post test scores (e.g. is there a magic # of weeks to attend, do only those who complete all 4 get significant change, etc.) _(Note: I'm not sure this one is answerable from data but thought it wouldn't hurt to ask on here)_   4. Is there a way to tell if EXTENDING the intervention from 4 to 6 weeks would have an impact of improving post test scores?


title: Can I use non-linear PCA (CATPCA?) on my multivariate dataset that contains nominal AND ordinal data?

text: In my multivariate dataset I have over 100 objects/cases that have been coded on 20 different variables. Several variables are ordinal and several variables are nominal. Is it possible to use non-linear principal components analysis as a dimension reduction technique on this dataset? I am not certain how to approach optimal scaling with the two types of variables simultaneously.


title: scatterplot smoothing in r with big dataset: different methods

text: I have a large dataset (>300,000 rows) with two variables. y is binary and x is continuous & numeric. I'd like to plot y and add smooth curve against x. I understand that loess(y~x) is a solution, but since I have such a big dataset, it takes too long to run, even if I set the 'cell' parameter to 500. Using scatter.smooth, it runs much faster and I think it also uses loess. but I have trouble understanding the parameter 'evaluation = 50'. Does this mean that it only uses 1/50 of data to produce the smooth curve? I also tried using geom_smooth, it would automatically switch to 'method=gam' since I have more than 1000 data points. but the curve looks different from the one I got using scatter.smooth (I guess that's normal as they are different models). My goal was just to see the pattern of the data. Which smoothing method should I use? Can I trust scatter.smooth? what's the difference between using loess and gam? below is the plot from scatter.smooth. It looks good, but it runs so much faster than the regular loess(). I'm not sure how it works... ![enter image description here](http://i.stack.imgur.com/n28Q6.png) Using the method whuber provided: ![enter image description here](http://i.stack.imgur.com/pEFV4.png) any help would be highly appreciated! Thanks


title: create bins for lognormal data for cluster analysis

text: I have a series of dollar amounts that are highly right skewed, but are roughly log-normal. i want to put this grouped dollar amount as a predictor variable into a latent class cluster analysis. In this case I'm not sure whether transforming into log-normal makes sense for a predictor variable, as opposed to just binning the raw, skewed data. Raw distribution: ![Frequency Distribution of Raw Dollar Amounts](http://i.stack.imgur.com/lxJyk.png) ![QQ plot of Raw Dollar Amounts](http://i.stack.imgur.com/9G3Jl.png) Log-normal distribution: ![Frequency Distribution of Logged Dollar Amounts](http://i.stack.imgur.com/iUohS.png) ![QQ plot of Logged Dollar Amounts](http://i.stack.imgur.com/4vlN6.png)


title: What are the Mild Regularity Conditions in the context of GEEs?

text: I was reading the Wikipedia Article on Generalized Estimating Equations, where I came across the following sentence: > Parameter estimates from the GEE are consistent even when the covariance > structure is misspecified, under mild regularity conditions. In this context, what is meant by _mild regularity conditions_?


title: How to derive the first order autocorrelation coefficient of an AR(1) process?

text: I want to derive the first-order autocorrelation coefficient of an AR(1) process, where the AR(1) process can be expressed as: $$y_t = \Theta y_{t-1} + u_t$$ with $|\Theta|<1$. Please give me a hint on how to start with this problem?


title: Combining information from multiple studies to estimate the mean and variance of normally distributed data - Bayesian vs meta-analytic approaches

text: I have reviewed a set of papers, each reporting the observed mean and SD of a measurement of $X$ in its respective sample of known size, $n$. I want to make the best possible guess about the likely distribution of the same measure in a new study that I am designing, and how much uncertainty is in that guess. I am happy to assume $X \sim N(\mu, \sigma^2$). My first thought was meta-analysis, but the models typically employed focus on point estimates and corresponding confidence intervals. However, I want to say something about the full distribution of $X$, which in this case would also including making a guess about the variance, $\sigma^2$. I have been reading about possible Bayeisan approaches to estimating the complete set of parameters of a given distribution in light of prior knowledge. This generally makes more sense to me, but I have zero experience with Bayesian analysis. This also seems like a straightforward, relatively simple problem to cut my teeth on. 1) Given my problem, which approach makes the most sense and why? Meta- analysis or a Bayesian approach? 2) If you think the Bayesian approach is best, can you point me to a way to implement this (preferably in R)? Related question EDITS: I have been trying to work this out in what I think is a 'simple' Bayesian manner. As I stated above, I am not just interested in the estimated mean, $\mu$, but also the variance,$\sigma^2$, in light of prior information, i.e. $P(\mu, \sigma^2|Y)$ Again, I know nothing about Bayeianism in practice, but it didn't take long to find that the posterior of a normal distribution with unknown mean and variance has a closed form solution via conjugacy, with the normal-inverse- gamma distribution. The problem is reformulated as $P(\mu, \sigma^2|Y) = P(\mu|\sigma^2, Y)P(\sigma^2|Y)$. $P(\mu|\sigma^2, Y)$ is estimated with a normal distribution; $P(\sigma^2|Y)$ with an inverse-gamma distribution. It took me a while to get my head around it, but from these links(1, 2) I was able, I think, to sort how to do this in R. I started with a data frame made up from a row for each of 33 studies/samples, and columns for the mean, variance, and sample size. I used the mean, variance, and sample size from the first study, in row 1, as my prior information. I then updated this with the information from the next study, calculated the relevant parameters, and sampled from the normal-inverse-gamma to get the distribution of $\mu$ and $\sigma^2$. This gets repeated until all 33 studies have been included.               # Loop start values values            i <- 2       k <- 1          # Results go here            muL      <- list()  # mean of the estimated mean distribution       varL     <- list()  # variance of the estimated mean distribution       nL       <- list()  # sample size       eVarL    <- list()  # mean of the estimated variance distribution       distL    <- list()  # sampling 10k times from the mean and variance distributions          # Priors, taken from the study in row 1 of the data frame            muPrior  <- bayesDf[1, 14]    # Starting mean       nPrior   <- bayesDf[1, 10]    # Starting sample size       varPrior <- bayesDf[1, 16]^2  # Starting variance            for (i in 2:nrow(bayesDf)){          # "New" Data, Sufficient Statistics needed for parameter estimation              muSamp    <- bayesDf[i, 14]          # mean         nSamp     <- bayesDf[i, 10]          # sample size         sumSqSamp <- bayesDf[i, 16]^2*(nSamp-1)  # sum of squares (variance * (n-1))          # Posteriors              nPost   <- nPrior + nSamp         muPost  <- (nPrior * muPrior + nSamp * muSamp) / (nPost)           sPost   <- (nPrior * varPrior) +                      sumSqSamp +                     ((nPrior * nSamp) / (nPost)) * ((muSamp - muPrior)^2)         varPost <- sPost/nPost         bPost   <- (nPrior * varPrior) +                      sumSqSamp +                     (nPrior * nSamp /  (nPost)) * ((muPrior - muSamp)^2)     # Update               muPrior   <- muPost         nPrior    <- nPost         varPrior  <- varPost          # Store              muL[[i]]   <-  muPost         varL[[i]]  <-  varPost         nL[[i]]    <-  nPost         eVarL[[i]] <- (bPost/2) / ((nPost/2) - 1)          # Sample              muDistL  <- list()           varDistL <- list()              for (j in 1:10000){           varDistL[[j]] <- 1/rgamma(1, nPost/2, bPost/2)           v             <- 1/rgamma(1, nPost/2, bPost/2)           muDistL[[j]]  <- rnorm(1, muPost, v/nPost)         }          # Store               varDist    <- do.call(rbind, varDistL)         muDist     <- do.call(rbind, muDistL)         dist       <- as.data.frame(cbind(varDist, muDist))         distL[[k]] <- dist          # Advance              k <- k+1          i <- i+1            }            var     <- do.call(rbind, varL)       mu      <- do.call(rbind, muL)       n       <- do.call(rbind, nL)       eVar    <- do.call(rbind, eVarL)       normsDf <- as.data.frame(cbind(mu, var, eVar, n))        colnames(seDf) <- c("mu", "var", "evar", "n")       normsDf$order <- c(1:33)      Here is a path diagram showing how the $E(\mu)$ and $E(\sigma^2)$ change as each new sample is added. ![enter image description here](http://i.stack.imgur.com/18Jkv.png) Here are the desnities based on sampling from the estimated distributions for the mean and variance at each update. ![enter image description here](http://i.stack.imgur.com/tNNRL.png) ![enter image description here](http://i.stack.imgur.com/dWqG3.png) I just wanted to add this in case it is helpful for someone else, and so that people in-the-know can tell me whether this was sensible, flawed, etc.


title: Distributions whose domain are themshelves a distribution

text: Dirichlet process, Pitman Yor processes and their variants are special types of distributions whose domain are themshelves a distribution. Is there a name to refer to these distributions in general? What other distributions have the same property?


title: Best approach to correlate discrete and continuous variables using machine-learning technique

text: I am using a scale which consists of discrete values: 0 (normal), 1 (mild), 2 (moderate), 3 (severe) and these values are for 4 components. I have measured 200 patients on this scale. Each patient has four discrete numbers. I am going to find the correlations between this scale and some continuous variables; these continuous variables are also measured for 200 patients. Goal: To generate these discrete values with a machine-learning algorithm based on the continuous variables.   1. What is best method for feature selection (for the continuous variables), Spearman correlation?  If the continuous variables were highly correlated with discrete values; does it mean I should consider them as features?   2. Which classification method would be the best?   3. Do you suggest I use Weka?


title: adf.test returning p-value > 0.99 even when data is stationary

text: I have a pair of time series that I differenced and now take on values {-1, -0.5, 0, 0.5, 1}. My goal is to test them for Granger causality and discover any lead/lag relationships. Since this requires the time series to be stationary, I ran adf.test on them and was surprised to find that the p-value for both were 0.99. I also got a warning saying p-value greater than printed p-value. I was hoping to learn what I'm doing wrong here, please let me know what further information I can provide.


title: Characterizing uncertainty in empirical PMF of unknown discrete random distribution

text: Suppose a discrete random distribution $X$ which takes integer values in some small (but not binary) range, e.g. $[0, 255]$ and whose true pmf is unknown (however, the null hypothesis is that it is uniform). If I draw a large sample $\\{x\\}$ from this distribution, I can estimate the pmf by counting the number of $x_i$ that take each value and dividing by the total. But if I draw another large sample, I will most likely get a different estimate. I'd like some way of characterizing the uncertainty in each point in the pmf, ideally one that can be updated as I (that is, the computer) continually draws new samples. "Confidence intervals for empirical CDF" looks closely related to what I'm trying to do but AFAICT only applies to _continuous_ distributions.


title: How hard is to convert a library from S-PLUS 8.0 to R?

text: I have a S-Plus library which I'd like to convert to R. I am a programmer, but I don't know anything about S-Plus or R. From my research it seems that they are highly compatibile. Is that true? The code I want to convert only uses core S-Plus libraries. I have attached a picture of the library as seen in the S-Plus 8.0 Object Explorer. Besides function source files, there are a few entries which I'm not sure how to transpose to R. For example the last 5 (oneDay, ...), which seem to be some sort of global variables, and have specific values assigned to them. What would be the equivalent of them in R? ![alt text](http://i.stack.imgur.com/102FK.png)


title: Statistical Tests That Incorporate Measurement Uncertainty

text: Suppose I am given two groups of mass measurements (in mg), which are referred to as y1 and y2. I want to do a test to determine if the two samples are drawn from populations with different means. Something like this for example (in R):               y1 <- c(10.5,2.9,2.0,4.4,2.8,5.9,4.2,2.7,4.7,6.6)     y2 <- c(3.8,4.3,2.8,5.0,9.3,6.0,7.6,3.8,6.8,7.9)     t.test(y1,y2)      I get a p-value of 0.3234, and at a significance level of 0.05 don't reject the null hypothesis that the two groups are drawn from populations with the same mean. Now I am given uncertainties for each measurement:               u1 <- c(2.3,1.7,1.7,1.7,2.0,2.2,2.1,1.7,2.3,2.2)     u2 <- c(2.4,1.8,1.6,2.3,2.5,1.8,1.9,1.5,2.3,2.3)      where u1[1] is the combined standard uncertainty in measurement y1[1] (and so on). How do I incorporate these uncertainties into the statistical test?


title: How to replace the missing data from AMELIA results

text: I have run a AMELIA imputation for a data set including missing data. I need to replace the missing point by the result of amelia(). But it content 5 group of imputed values. How can i choose the best one to replace the missing values (to plot a graph of data set after imputing)


title: Is $f_{X,Y}(X,Y)$ the same as $f_{X,Y|X}(X,Y|X)$

text: This might be a very simple question to answer, but for some reason I have just been banging my head against the wall for a little while now. To me, it makes some intuitive sense that $f_{X,Y}(X,Y)$ is same as $f_{X,Y|X}(X,Y|X)$, but I cannot for the life of me figure out how to show it. However at the same time, if those two joint distributions were the same, wouldn't that imply that $X$ and $Y|X$ are independent? Which definitely seems off to me. Thanks in advance!


title: Statistical similarity of time series

text: Supposing one has a time series from which one can take various measurements such as period, maximum, minimum, average etc. and then use these to create a model sine wave with the same attributes, are there any statistical approaches one can use that could quantify how closely the actual data fit the assumed model? The number of data points in the series would range between 10 and 50 points. A very simplistic first thought of mine was to ascribe a value to the directional movement of the sine wave, i.e. +1 +1 +1 +1 -1 -1 -1 -1 -1 -1 -1 -1 +1 +1 +1 +1, do the same to the actual data, and then somehow quantify the degree of similarity of directional movement. Edit: Having given more thought to what I really want to do with my data, and in light of responses to my original question, what I need is a decision making algorithm to choose between competing assumptions: namely that my data is basically linear (or trending) with noise that could possibly have cyclic elements; my data is basically cyclic with no directional trend to speak of; the data is essentially just noise; or it is transitioning between any of these states. My thoughts now are to maybe combine some form of Bayesian analysis and Euclidean/LMS metric. The steps in this approach would be Create the assumed sine wave from data measurements Fit a LMS straight line to the data Derive an Euclidean or LMS metric for departures from the original data for each of the above Create a Bayesian prior for each based on this metric i.e. 60 % of the combined departures attach to one, 40 % to the other, hence favour the 40 % slide a window one data point along the data and repeat the above to obtain new % metrics for this slightly changed data set - this is the new evidence - do the Bayesian analysis to create a posterior and change the probabilities that favour each assumption repeat along the whole data set (3000+ data points) with this sliding window (window length 10-50 data points). The hope/intent is to identify the predominant/favoured assumption at any point in the data set and how this changes with time Any comments on this potential methodology would be welcome, particularly on how I could actually implement the Bayesian analysis part.


title: Constructing a Gaussian from the posterior distribution

text: I have a Linear Dynamical System as in the following graphical model: ![enter image description here](http://i.stack.imgur.com/qUDHu.png) Every random variable is a scalar; we know $B$ and $x_1$ a priori. $A$ comes from a Gaussian, the variances $R$ and $Q$ come from Inverse Gamma distributions. There a $N$ hidden states; $x_{1:N}$ and $N$ observations; $y_{1:N}$. The state transition and emission densities are Gaussians with: $$ x_t \sim N(x_t|Ax_{t-1},R) $$ $$ y_t \sim N(y_t|Bx_{t},Q) $$ What I want is to calculate the posterior distribution of $x_t$ given all other variables. I brought the calculations up to the following: $$ P(x_t|x_{1:t-1},x_{t+1:N},y_{1:N},A,R,B,Q) = \frac{N(x_t|Ax_{t-1},R)N(y_t|Bx_t,Q)N(x_{t+1}|Ax_{t},R)}{\int_{-\infty}^{\infty}N(x_t|Ax_{t-1},R)N(y_t|Bx_t,Q)N(x_{t+1}|Ax_{t},R)dx_t} $$ $$ \propto \exp\left(\frac{-(A^2Q+Q+RB^2)}{2RQ}\left(x_t - \frac{AQx_{t-1}+AQx_{t+1}+RBy_t}{A^2Q+Q+RB^2}\right)^2\right) $$ If I made all calculations correctly, this looks like the exponential part of a Gaussian distribution. Now what I am wondering is whether sampling from this Gaussian is equivalent to sampling from the full posterior $P(x_t|x_{1:t-1},x_{t+1:N},y_{1:N},A,R,B,Q)$, since this final Gaussian is only equal to the posterior up to a normalizing constant. If the answer is yes, how come it can be possible since the distributions are not exactly equal and only proportional? Thanks in advance.


title: Likelihood vs. noise kernel hyperparameter in GPML Toolbox

text: I'm using GPML toolbox by C.E.Rasmussen to solve the basic GP regression problem (presented in the book) with noisy observations. That is to say, estimate the underlying function $f$ of a static noisy mapping $$y = f(\mathbf{x}) + e, \qquad e \sim \mathcal{N}(0, \sigma^2)$$ from a set of training examples $\\{ (\mathbf{x}_i, y_i) \\}_{i=1}^{n}$. As far as I understand it, I should respect the noisiness of observations by choosing the kernel as the sum $$ k(\mathbf{x}_i, \mathbf{x}_j) = k_f(\mathbf{x}_i, \mathbf{x}_j) + \sigma^2_{e}\delta_{ij}$$ where the final term in the sum is the kernel of the white noise (that is, noise of the observations). When using GPML toolbox, for those who are familiar, you have to specify a likelihood. In my case I chose Gaussian likelihood, which has one hyperparameter - in the code documentation this correponds to formal parameter $s_n$. So all together, when I perform optimization, I have one hyperparameter for the noise kernel ($\sigma_e$), one for the likelihood ($s_n$) and $d$ (say) hyperparameters for the $k_f$. I am confused about the meaning of the hyperparameters $\sigma_e$ and $s_n$. Which one of the hyperparameters ($\sigma_e$ or $s_n$) represents the variance of the noise in the observations? If the Gaussian likelihood is the measurement model, then $s_n$ should be the variance of the observations $y_i$, but then for what reason do we add the noise kernel (with additional hyperparameter ($\sigma_e$), which I think is redundant, at this point since we already have $s_n$ to do the job)? Perhaps they're one and the same and should be tied together during optimization. I'm confused. GPML code for exact inference:               [n, D] = size(x);     K = feval(cov{:}, hyp.cov, x);  % evaluate covariance matrix     m = feval(mean{:}, hyp.mean, x); % evaluate mean vector     sn2 = exp(2*hyp.lik); % noise variance of likGauss     if sn2<1e-6           % very tiny sn2 can lead to numerical trouble     L = chol(K+sn2*eye(n)); sl =   1; % Cholesky factor of covariance with noise     pL = -solve_chol(L,eye(n));  % L = -inv(K+inv(sW^2))     else     L = chol(K/sn2 + eye(n)); sl = sn2; % Cholesky factor of B     pL = L; % L = chol(eye(n)+sW*sW'.*K)     end     alpha = solve_chol(L,y-m)/sl;      `sn2` is the likelihood parameter, `hyp.cov` contains the kernel hyperparameters (including the noise kernel hyperparameter $\sigma_e$)


title: Finding the posterior pdf

text: > Suppose $X$ has probability density function $$f(x, \theta) = \theta > e^{-\theta x}$$ when $x > 0$ and $\theta > 0$, and $0$ otherwise; given > $\Theta = \theta$. Suppose the prior probability density function of > $\Theta$ is $$h(\theta) = 1$$ when $0 < \theta < 1$, and $0$ otherwise. Find > the posterior probability density function of $\Theta$ given $X = x$ (for $x > > 0$). Let $k( \theta | x)$ denote the posterior pdf. We have $k( \theta | x) = \frac{L(x | \theta)h(\theta)}{f_1(x)}$, where $f_1(x)$ is the joint pdf of **$X$**. We have $L(x | \theta)h(\theta) = (\theta^n e^{ - \theta \sum^{n}_{i=1} x_i})(1) = \theta^n e^{ - \theta \sum^{n}_{i=1} x_i}$. Let $Y= \sum^{n}_{i=1}$ We also have $f_1(x) = \int^{\infty}_0 \theta^n e^{-\theta Y} d \theta$. But since $\Gamma(n+1) = \int^{\infty}_0 \frac{\theta^{(n+1)-1} e^{-\theta Y}}{(1/Y)^{n+1}} d\theta$, we have $f_1(x) = \int^{\infty}_0 \theta^n e^{-\theta Y} d \theta = \frac{\Gamma(n+1)}{Y^{n+1}}$. So $k(\theta | x) = \frac{(\sum x_i)^{n+1} \theta^n e^{-\theta \sum x_i}}{\Gamma(n+1)}$ for all $x_i > 0$ and $0< \theta < 1$, and $0$ otherwise. Do you think my answer is correct? Thanks in advance


title: Dependent (Paired) T-Test Sample Variance

text: When conducting a repeated measures test where one measures the before/after values for several pairs, how does one account for variance within the pairs themselves? For example lets look at this example of 3 student scores SAT math section as the result of a prep-course. ![example](http://i.stack.imgur.com/AdDwo.png) At this level, it seems that the test was significant in improving their overall score. However, what if I now add in the information that each student has taken the exam 3 times prior to the class and their scores varied +/- by 50 points each time. Now the 30 point increases no appear to be relevant. How does one account for these types of scenarios when running paired t-tests?


title: AIC doesn't agree with model checking

text: I have two glm, one with a gaussian distribution and identity link and one with gamma family and log link. The predictors are the same, the only thing that change is the response that in the gaussian glm is log-transformed and in the gamma glm is not. The deviance is lower in the second model, and when plotting the residuals the second model looks much better than the first. However, when looking and the AIC, the first model has half the AIC of the second! How is that possible? Is it wrong to compare AIC of glm with different distribution families?               summary(mod1)          Call:     glm(formula = log_RS ~ DIET + log_Disp + log_II + log_LS + log_SM +          log_AS + log_LONG, family = gaussian)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -5.1050  -0.2672   0.1063   0.3945   1.3012            Coefficients:                  Estimate Std. Error t value Pr(>|t|)         (Intercept)  1.061888   0.008574 123.856  < 2e-16 ***     DIETOH      -0.077913   0.012762  -6.105 1.07e-09 ***     log_Disp     1.529989   0.006960 219.822  < 2e-16 ***     log_II      -0.422808   0.006786 -62.310  < 2e-16 ***     log_LS       0.325270   0.006093  53.382  < 2e-16 ***     log_SM      -0.268292   0.008372 -32.048  < 2e-16 ***     log_AS       0.185077   0.010166  18.205  < 2e-16 ***     log_LONG     0.067347   0.010494   6.418 1.45e-10 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for gaussian family taken to be 0.3246263)              Null deviance: 24060.5  on 9345  degrees of freedom     Residual deviance:  3031.4  on 9338  degrees of freedom     AIC: 16018          Number of Fisher Scoring iterations: 2          > summary(mod2)          Call:     glm(formula = RS ~ DIET + log_Disp + log_II + log_LS + log_SM +          log_AS + log_LONG, family = Gamma(link = log))          Deviance Residuals:           Min        1Q    Median        3Q       Max       -2.94338  -0.35264  -0.00239   0.27669   1.40915            Coefficients:                  Estimate Std. Error t value Pr(>|t|)         (Intercept)  1.186181   0.006435 184.332  < 2e-16 ***     DIETOH      -0.075002   0.009579  -7.830 5.41e-15 ***     log_Disp     1.536406   0.005224 294.104  < 2e-16 ***     log_II      -0.364580   0.005093 -71.585  < 2e-16 ***     log_LS       0.306213   0.004573  66.955  < 2e-16 ***     log_SM      -0.233847   0.006283 -37.216  < 2e-16 ***     log_AS       0.156994   0.007630  20.575  < 2e-16 ***     log_LONG     0.053129   0.007876   6.745 1.62e-11 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for Gamma family taken to be 0.1828775)              Null deviance: 23484.5  on 9345  degrees of freedom     Residual deviance:  2349.8  on 9338  degrees of freedom     AIC: 33167          Number of Fisher Scoring iterations: 8          >


title: Probability of being included in a Ranked Set Sample

text: I am doing a research work which needs a probability calculation. I'll be grateful if the community cross validates the formula that I have derived using the idea of binary strings and following a great answer by whuber to this question. But I really apologize to you for such a long thread! In my design, I am planning to allocate a patient to either treatment policy B1 or to treatment policy B2 upon their remission times measured in days. Say, there are $mn^2$ observations in total. For the following picture $m=2$ and $n=4$. We divide the $mn^2$ patients randomly into $m$ groups of size $n^2$ and then randomize the $n^2$ patients again into $n$ blocks each having $n$ observations. This is all the same as dividing the $mn^2$ patients into $mn$ blocks of size $n$. Then we allocate the patients with the smallest remission time in block-1, 2nd smallest in block-2, 3rd smallest in block-3 up to the largest of block-n (here n=4) to treatment B1. Likewise we allocate patients to treatment B2. ![enter image description here](http://i.stack.imgur.com/hVyTN.jpg) Now the original data may contain tied observations. So I wanted to extend whuber's answer for the case where $x_1<=x_2<=\ldots<=x_n$ can be possible. Here $x_i$ are the remission times that may or may not have tied observations. This time I want to calculate the probability that $x_i$ will be allocated to treatment B1. **So, here is my generalization:** After we sort the data in `R`, let us define the following notations: $r1$=Number of tied observations that come before $x_i$ in the sorted data $r2$=Number of tied observations that come after $x_i$ in the sorted data $p$=Number of ones in the binary string (please see whuber's answer that was mentioned earlier) that are placed in the $r1$ places before $x_i$. $q$=Number of ones in the binary string that are placed in the $r2$ places after $x_i$. $k$ denotes the order of $x_i$ in the RSS block. So when, $r1=r2=0$ then $p=q=0$ and $(p+1+q)$=The number of tied observations along with $x_i$ in the RSS block. With these notations according to my calculation, the probability that $x_i$ will get B1 is given by- $prob(i)=\sum_{k=1}^{n}\frac{\sum_{p=0}^{min(r1,n-1)}\sum_{q=0}^{min(r2,n-1-p)} \dbinom{i-r1-1}{k-1-p} \dbinom{r1}{p} \dbinom{r2}{q} \dbinom{mn^{2}-i-r2}{n-k-q}}{\dbinom{mn^2}{n}}$ > Since each RSS block may contain tied observations to $x_i$ and then all of > the tied observations will have the same probability of being selected as > the k-th smallest of the RSS block (with probability $1/(p+1+q)$, I guess), > so do I need any adjustment? Although the the sum of the probabilities (over > all $i$) of being the k-th smallest of a particular RSS block seems to be 1 > always, which is the total law of probability. **For further generazation...** Can anyone suggest me if there were censoring and I had $C_i$, the censoring time for each observation, then what adjustment to the formula should be made?


title: How to differentiate hype from word of mouth in Google Trends?

text: I'm new to data analysis. I was wondering how I could recognize hyped (fake) growth and word of mouth (true) growth in Google Trends specifically (but I guess it applies even to stock markets and the rest). For example, the image below is based on the Google Trends plot when you search for 'iPhone': ![iPhone trend](http://i.stack.imgur.com/pYHRS.gif) _(Orange, magenta, and green annotation were added to the original.)_


title: Kaplan-Meier p-values

text: I have breast cancer trial data where, for subgroups 1 and 2, two treatment arms are being compared. The results are:   * For **subgroup 1** , `N=33`, median survival in `trt1=10mo`, median survival in `trt2=5.7mo`, p-value for difference`=0.18`.    * For **subgroup 2** , `N=16`, median survival in `trt1=10mo`, median survival in `trt2=5.5mo`, p-value for difference`=0.003`. How can the p-value be significant for subgroup 2 with smaller N and median OS being similar in both subgroups?


title: One-sample median test:: Signrank seems unsuitable

text: I have some data on successive bets made by customers. I want to see whether there is a statistically significant change in bet stake with each subsequent bet. The data is skewed, so I have equalized all initial bets to 1. Thereafter, I calculate the proportional change in bet stake. So if someone had a first bet of 15 and a third bet of 30, the value for the stake in the third bet shows 2.00. For a series of groups, I'm showing the mean standardized stake at the second, third and fourth round of betting. Although the data is quite skewed, I've used a t-test to test whether there is a statistically significant difference between this mean and 1 (i.e. no change in bet size since the initial bet). Given the data is skewed, I also want to perform a non-parametric test, but I'm having an issue with the median and a one-sample median test. Although for many groups the median is 1, the median test shows a statistically significant difference from 1. For example, one groups has approximately 4000 observations, of which 1800 are `ties/zero` or equal to exactly 1. When performing a `signrank` test, 1800 observations are showing us as zero, while of the remaining 2200 observations, there are 200 more positive counts than negative, so the result of the one-sample median test is highly significant. It doesn't seem right. Should I be using the signrank test for this at all? Would another test be more appropriate?


title: Can likelihood be changed when the prior changes?

text: I have a data which follows gamma distribution and want to know the uncertainty of the parameters of this data.   * $\text{Data} \sim \text{Gamma} (\alpha, \beta)$   * Parameters $\alpha \sim \text{Gamma} (k_\alpha, \theta_\alpha)$ $\beta \sim \text{Gamma}(k_\beta, \theta_\beta)$ I used Winbugs (code below).               model{       for (i in 1:N){          Y[i] ~ dgamma(k, theta)       }            k ~ dgamma(0.1, 0.1)       theta ~ dgamma(0.1, 0.1)     }        1. To plot the likelihood I used a uniform prior then divided the posterior by the prior which makes the posterior same as the likelihood. (Figure1) ![enter image description here](http://i.stack.imgur.com/GCkaq.png)   2. Next I changed the prior several times and looked what happens, but the problem is when I plotted the likelihood by dividing the posterior by the prior, likelihood changes which shouldn’t, whenever I change the prior. (Figure 2, 3) ![enter image description here](http://i.stack.imgur.com/uiB9R.png) ![enter image description here](http://i.stack.imgur.com/1cAQo.png) Question Is it possible that the likelihood changes when the prior changes? Can somebody help me what the problem is? If the prior is too narrow, is there a possibility that the posterior might be wrong?


title: Common or regional time dummies as controls: Does it make a difference?

text: I am analyzing policy diffusion using panel data, and the inclusion of temporal controls substantially influences my results. I'm wondering if it can theoretically make a difference for the estimates of my **other** variables if I include common time dummies for all units of observation, or if, e.g. I use separate time dummies for EU and non-EU countries. Shouldn't the common time dummies capture EU specific shocks as well?


title: Non-normally distributed data -- what can I do with the z-scores?

text: layman here... I have some sets of several hundred scores that are not distributed normally. I want to normalize these scores onto a percentile scale. The mean-median-mode for these look like this:               3-2-0,     1-1-0,     9-3-0,     1-5-5,      I have calculated the variance/std dev and z scores for these, but now I'm not sure what to do with this info. What good are z scores when the distribution isn't normal? What's the best way to get a normalized percentile score for these?


title: What is the meaning of the PRESS statistic in the case of an exact interpolation?

text: I'm trying to teach myself about cross-validation – I'm new at this – and I'm having a conceptual problem. I'm reading the book "Classical and Modern Regression" by Myers, chapter four, which is on cross-validation for model selection, specifically, the section on the PRESS statistic. The intuitive definition of the PRESS statistic – which is defined as the sum of the square residuals that result when you drop a data point, recalculate the regression with the data point removed, and then calculate the residual between the new regression and the dropped data point – makes perfect sense to me. But the PRESS statistic is not usually calculated this way; it's normally calculated by taking advantage of the Sherman–Woodbury–Morrison theorem, which allows you to compute the PRESS residuals in terms of the ordinary residual, which looks like: $(y_i - \hat{y}_i)/(1 - x_i^t (X^t X)^{-1} x_i) $ Where the numerator is just the residual of the original regression at data point _i_. But what if the original regression exactly interpolates the data point, so that the numerator is zero? The intuitive definition of the PRESS seems to indicate that the data point gets dropped, and the new regression would be unlikely to still interpolate the point, so the PRESS residual should be nonzero. But the equation given above indicates that it is either exactly zero or, if the denominator is zero, it's undefined. Clearly I'm missing something. But what?


title: R package to handle gene expression and proportional hazards model

text: I have survival data with RNA expression. A lot of the data is censored. Is there an R package that does a proportional hazards model where the number of covariates is much larger than the number of samples? The only one I could find is glcoxph. Are there any others?


title: What is the best way to simulate data for a linear regression model?

text: I am concerned with simulating data for a linear regression model. I need to control the means, variances, and correlations (covariances) between the predictors and the criterion variable. In addition, I need to be able to vary the explained variances ($R^2$). It is obvious to me that the latter must be a function of the earlier, so at least one correlation (covariance) in $\Sigma$ is perhaps dependent on the choice of $R^2$, where $\Sigma=E((Y,X)(Y,X)^T)$, for centered $X$ and $Y$, is the variance-covariance matrix of all variables. My plan is thus as follows:   1. Specify $\Sigma$, means, and $R^2$   2. Simulate data with these sufficient statistics, e.g. by sampling from the multivariate normal   3. Check estimated $\beta$ (regression coefficient) vector against population (theoretical) coefficients and use the model for unrelated tests/science. Hence, my approach does not suggest specifying $\beta$ but letting the coefficients be a function of population $\Sigma$, means, and $R^2$. The reason I need to do this is to attribute some realisitc scale to $X$ and $Y$ (e.g., let $Y$ assume an 'income' scale and give $X$ a realistic scale for years of eduction). Therefore, I specify the sufficient statistics instead of regression coefficients. But maybe there is a better way. Moreover, I have two specific questions:   1. Given the population variance-covariance matrix $\Sigma$ of one criterion variable $Y$ and a series of predictors (covariates) $X$ , I would like to calculate the vector of true population regression coefficients. Of course, I could simulate data $X$ and $Y$ and use the OLS estimator, but there should be a direct way to use $\Sigma$ in the estimation of population $\beta$?   2. Which options are there to specify covariances (correlations) in $\Sigma$ given I need a fixed $R^2$ of a linear regression of $Y$ on $X$? This, to systematically vary the explanatory power of the regression model.


title: Time series with same variables but different respondents

text: I would like to measure the impact of an intervention using a form of time series or some related technique (like segmented regression). I would like to use the Schools and Staffing Survey, which assesses my variable of interest over time, however, the respondents are different at each time point. I assume different respondents would present a problem for this kind of analysis. Can anyone suggest another technique that would allow me to gauge the impact of the intervention?


title: Predictions on statistical series

text: I need to summarize possible approaches (as well as their pros and cons) to the problem of statistical series prediction. It's pure theoretical question, an example below was given only for convenience and specificity. First of all, sorry for the term 'statistical series'. Question originally asked in another language, and I had hard time translating it. But what kind of series it is should be obvious from example. An example represents series of dice throwing experiment. We have 20 observations and the observed frequency distribution is described by the table below. $\begin{matrix} \text{dice face} & 1 & 2 & 3 & 4 & 5 & 6 \\\ \text{sum} & 2 & 6 & 8 & 2 & 1 & 1\\\ \text{frequency} & 0.1 & 0.3 & 0.4 & 0.1 & 0.05 & 0.05\end{matrix}$ Let's return to our task now. We have to derive some predictive model now.   I was able to think of only 1 approach:   * We calculate standard error and mean, then able to make and prove any hypothesis about future observations with use of t-test. A negative aspect of this approach is the involvement of uncertainty in the answer (confidence intervals, etc.)  When I gave it further thought, I may not understand very good process of prediction here. If continue sampling, our histogram will get more and more precise, but what would be the value to predict at first place? density function, population mean? Due to last paragraph, I don't see data mining machinery applicable here (regressions, neural nets, decision trees). I would appreciate insights and clarifications.


title: Generating functional-form PDF from Max Likelihood Estimation

text: For the purpose of this question, please consider me a stats newbie. I'm working on a (very fun!) research project which involves estimating a pdf of "personal values" -- i.e. how much a certain person values a certain good they can buy. I wanted to find the functional form of the PDF so that I could generate data according to the distribution I saw in the empirical data I've been working with. But I received a suggestion to follow Maximum Likelihood Estimation, which isn't something I've done in the past. I've been experimenting but my data doesn't seem to likely fit any typical distribution... it's more of a "bactrian camel" distribution. **So what I'm interested in knowing is:**   1. Is my current approach correct? (see below)   2. How might I approach this with Maximum Likelihood Estimation? **Here's what I've done so far:**   1. Constructed a CDF of the data I have (taking them as "valuations")   2. Literally... found the OLS estimator of the CDF, in functional form, so now I have a polynomial estimation of the CDF   3. Took the derivative of that function (which should equal the PDF)   4. Generated data according to that derivative. Overall, this approach has performed fairly well. I don't have any metrics to compare my generated data with the actual data -- I don't know any! -- but from comparing the two datasets, they're fairly close in mean and variance. Much appreciated! \-- An enthusiastic college kid Edit: Adding details: Here's what my density function looks like. I'd like to approximate a functional form of this, because I need to do two things:   1. Infer the distribution of data above 150; according to theory, there exist data that is above 150, but it's not in this distribution.   2. Generate random data (say a sample size of 5,000 -> 10,000) ![](http://i.stack.imgur.com/pA3vd.png) The distribution I approximated through my take-derivative-of-CDF-estimation looks like this (ignore the improper scaling of y-axis -- this is scaled up) ![enter image description here](http://i.stack.imgur.com/Od6tK.png)


title: Sample size needed for test of binomial using Fishers exact test

text: I am running two different advertisements served up randomly. I monitor how many people respond to each advertisement. I am trying to determine if one ad is better than the other. If I run a Fisher's significance test on the results I can find if the null hypothesis can be rejected at any point in time. I think it can be confidently said that the two ads will have different response rates no matter how small the difference, but how can I test that the results are no more than 5% different in terms of response rate?


title: How can you summarize growth in the "production" rate per person?

text: I have been reading the book, _Good to Great_ by Jim Collins. It says that you should focus on information that "cannot be ignored," and that you should focus on the right "denominator." I thought it would be good to measure "Average Growth Rate per Producing Entity." If we are providing value and support, we would expect to be helping them "produce more." I have production data in an excel spreadsheet that I can summarize. I primarily use "pivot tables" for my data analysis. I know that you can do "percentage change over previous time period" but I end up with distortions for the first year, when people don't have a full period to produce. I also am not sure how to summarize it. Is there a better program that I should use for this? Is there a good way to do this in excel? I was thinking you could take the geometric mean (what growth rate would you need on the initial production number in order to end up with the final number -- or so I believe it is explained) and besides that I'm not sure where to start. I would love to use pivot tables or something similar where I can just hit refresh and get the updated version of the numbers. Any ideas out there?


title: Specify number of linear discriminants in R MASS lda function

text: I have a matrix with 1024 features and 10000 samples with a label vector of three different classes. I use R MASS library `lda` function to compute the model and get the coefficients for linear discriminants:               > library (MASS)     >      > ldamodel = lda(data$X, data$y)     > head (ldamodel$scaling)                  LD1         LD2     V5  -0.053074978  0.14565211     V6  -0.009618016 -0.11198306     V7  -0.003863230  0.28189459     V8   0.063191889 -0.26726050     V9  -0.029950632  0.16121364     V10  0.017988965 -0.01584389     ...      The model has two linear discriminants, which is ok. But I'd like to know, if it is possible make the `lda` function create more than two linear discriminants and how to do it in R? I'm not yet very faimiliar with LDA, so I'm sorry, if this question is a statistical nonsense.


title: How to understand kernel functions and how to choose a suitable kernel?

text: I am trying to describe my understand of kernels in the Support Vector Machine(SVM) and why some of them are more popular, but I am not sure if I misunderstand these concepts: 1) There are a large number of kernel functions in SVM, such as RBF, polynomial, linear,etc. Generally, they can be separated into two groups: linear and non-linear. For linear kernels, they can separate data linearly in original space. Non-linear kernel functions separate data using non-linear boundaries in the original space (but after mapping, fitting the Mercer's theorem, the data are still separated linearly in the higher dimensional feature space). This makes non-linear kernel more flexible if data cannot be separated linearly in original space. 2) Compared to other non-linear kernels, why the RBF kernel is more popular in SVM? Can I assume the reason is that this kernel is infinitely differentiable (Gaussian), and this means that the data can be mapped to an infinite dimensional feature space. Other non-linear kernels can also do this, for instance, by adding parameters, a polynomial kernel can also map data to a very high dimensional feature space, but this is not convenient. If my understand is correct? 3) In Gaussian Processes (GP), we use kernel functions to measure the similarity between different data points. What characteristics of the squared exponential (SE) kernel make this kernel popular? (I know it is Gaussian and infinitely differentiable, any other reasons?) By the way, if I want to make dynamic data analysis (the data are not periodic), is this kernel still suitable and how to select a suitable kernel function? Thank you very much for your help and patience in advance.


title: Sampling with replacement

text: Suppose you have 40 different books (20 math books, 15 history books, and 5 geography books). Let M = math books, H = history books, G = geography books You pick 5 books at random, with replacement, one at a time. What is the probability that you've picked books from at most two disciplines? I know this is basic, but... The way I have approached this problem is: The only way you can have at most two disciplines is to have   1. { H H M M G }   2. { H H G G M }   3. { G G M M H } so the answer I got is: $$ \left(\left(\frac{15}{40}\right)^2\cdot\left(\frac{20}{40}\right)^2\cdot\left(\frac{5}{40}\right)\right) + \left(\left(\frac{15}{40}\right)^2\cdot\left(\frac{5}{40}\right)^2\cdot\left(\frac{20}{40}\right)\right) + \left(\left(\frac{5}{40}\right)^2\cdot\left(\frac{20}{40}\right)^2\cdot\left(\frac{15}{40}\right)\right) $$ I have a feeling I am doing this problem wrong. Is there an easier way to do this? I tried to do it with combinations, but the "with replacement" aspect threw me off a bit. I feel as though order should not matter here, but I do not know whether that is the case. Any help would be appreciated! Thanks.


title: Dynamically generating regression coefficients and tracking their change over time

text: I'm running a generalized linear model (quasi-poisson regression) as a cron job in R that trains on data from an SQL query. The SQL query pulls data from the last 30 days. Depending on the sample of data from the last 30 days, the regression coefficients of course change. As a result, I'm then writing the regression coefficients generated in R back to an SQL table, with year/month/day as the the primary key. Thus, I can see the daily change in regression coefficients based on the most recent 30 days of data. My question is this: how do I best interpret the change in regression coefficients over time? If I can see a broad overview of techniques, I can hopefully find the one most appropriate given my data. My objective is to reduce the RMSE of a predictive model. The response variable is the number of events in the 30 day interval from the current time. I have multiple predictors, including interaction terms. Edit: it is not a requirement to only query the last 30 days of data. I simply wish to weight data closer in time more heavily. A discrete-time (days, in this case) weighted regression on all data would probably be ideal.


title: Problem plotting GLM data of binomial proportional data 2

text: I am trying to plot the predicting line for a GLM with proportional data. After having troubles with the example from "Crawley" I found this link: Problems plotting GLM data of binomial proportional data My data consists out of 3 columns: Altitude, number of beetles and number of other insects. I want to investigate the change of beetle proportion throughout the different altitudes. After importing the table I calculated the proportion as shown in the topic linked above. Then               y=cbind(mean$beetles,mean$rest)        model=glm(y~altitude,family=binomial(logit),data=data)        plot(p~data$altitude, ylab="Beetles (proportional)", xlab="Altitude")     p is the propotion i calculated before     lines(mean$altitude,predict(model,type="response"))      I then get this plot: ![plot](http://s12.postimg.org/5osjjbjyl/Rplot.png) this is obviously not correct but i cannot figure out what is wrong here. I would be very pleased if you helped me. this is the little data set:                altitude,beetles,rest      304,88,35      842,7,18      877,21,27      1505,2,25      1118,11,19      702,25,24      806,16,17      960,16,32      975,34,61      1112,11,17      1348,11,32      502,58,24      1428,2,13      1521,0,35      1040,9,6      1412,8,19      1231,5,19      980,14,16      1053,2,10      Cheers Peter


title: Standard deviation of standard deviation

text: What is an estimator of standard deviation of standard deviation if normality of data can be assumed?


title: Incorporating boolean data into analysis

text: I have a data set of about 3,000 field observations. The data collected is divided into 20 variables (real numbers), 30 boolean variables, and 10 or so look up variables and one "answer" variable We have about 20,000 objects in the field, and i'm trying to produce an "answer" for the 20,000 objects based on the 3,000 observations. What are some of the available methods that incorporate booleans and look up tables? any suggestions on how i should proceed? **EDIT** the answer variable is a boolean as well **EDIT 2** a sample of the variable data:   * Age of specimen   * length, area, volume   * time since last inspection   * height   * design life Lookup table   * material type    * coating type   * design standard   * design effectiveness a sample of the boolean   * is it inspected?   * is it in bad shape   * does it need repairs soon the answer variable which is my f(x) is:   * is it useable


title: Alpha adjustment for main and interaction effects in multiple factorial ANOVAs

text: let's say that I performed 10 2x2 factorial ANOVAs. In total I get 30 effects (two main effects and one interaction for each ANOVA). When adjusting the alpha level for multiple comparisons, should I consider 30 hypothesis tests or only 10? Formally they are not 30 comparisons, since they are "nested" into triplets in each ANOVA.


title: Looking for 2D artificial data to demonstrate properties of clustering algorithms

text: I am looking for datasets of 2 dimensional datapoints (each datapoint is a vector of two values (x,y)) following different distributions and forms. Code to generate such data would also be helpful. I want to use them to plot / visualise how some clustering algorithms perform. Here are some examples:   * star like cloud data   * four clusters, one easy seperable   * a spiral (no cluster)   * a ring   * two barely seperated clouds   * two parallel clusters forming a spiral   * ... etc


title: How are piecewise cubic spline bases constructed?

text: There are words from the _The Elements of Statistical Learning_ on page 119: > It is not hard to show that the following basis represents a cubic spline > with knots at $\xi_1$ and $\xi_2$:   * $h_1(X)=1$    * $h_2(X)=X$   * $h_3(X)=X^2$   * $h_4(X)=X^3$   * $h_5(X)=(X-\xi_1)_+^3$   * $h_6(X)=(X-\xi_2)_+^3$ Then my question is that, according to the numeric analysis, piecewise cubic polynominals tries to fit the function in region $[\xi_1,\xi_2]$ with polynomials at order 3, so I think the basis should be   * $h_1(X)=f(\xi_1)$   * $h_2(X)=(X-\xi_1)$   * $h_3(X)=(X-\xi_1)^2$   * $h_4(X)=(X-\xi_1)^3$. There are at most 4 basis, so where am I wrong?


title: Adjusting data for age and bmi

text: I am doing a study on the effects of bad blood sugar control on patients. BMI, age, type of diabetes (Type 1 and 2 or gestational diabetes) are all variables. Outcomes include C section, early labour etc. A reviewer told me > Comparisons should be adjusted to maternal age and baseline BMI. How does one do this with an online calculator? I have been doing my calculations only with online calculators and have never ventured to use SPSS. Can anyone advise me?


title: normalizing dataset for extrapolation - sample or population mean and standard deviation?

text: I am currently fitting models that are intended to be used for extrapolating from a limited sample to a large population. For a specific example, one model is predicting water temperature in rivers based on attributes of the river. The sample includes data from approximately 1,000 river segments (a river segment is a unique length of a river), which will be used to extrapolate to over 100,000 river segments. This is a neural network model, but I have the same question for a different model using hierarchical logistic regression models, so the question is not statistical method specific. I have read and have found through experience that normalizing my modeling dataset can improve model fit. However, I am unsure in this instance because the goal is extrapolating to the population. Should the predictors be normalized [(y - mean)/stdev] based on their mean and sdev in the sample or in the population of river segments? I have taken the opinion that I should normalize based on the population to help ensure that the population variability is represented adequately in the modeling dataset. I apologize if there is a cross posting, as I have searched for this answer thoroughly and have not seen a similar question or discussion anywhere. So I am asking this here. Please point me to an answer if one exists elsewhere.


title: Transforming a uniform-on-sphere random vector

text: Consider the 3-D real random vector $(X_1,X_2,X_3)$ which is uniformly distributed on the surface of a unit sphere. What can be told about the distribution of $(aX_1,bX_2,cX_3)$, where $a,b,c,$ are non-zero and non- identical real constants? Is it true to say that they are distributed uniformly on the surface of an ellipsoid with the corresponding parameters $a,b,c$?


title: How to create the class range given 1 as the middle class, 3 as the highest value and 0.65 as the lowest value?

text: Good day. I know getting the class interval given 3 as the highest value and 0.65 as the lowest value is easy. Here's the catch, the distribution of the interval starts at 1 which is considered as the base line or the acceptable value. So first I tried calculating for the range using the usual ways of doing way of doing things.. Number of desired classes: 5 Class Interval = highest value - lowest value / number of classes ... = 3 - 0.65 / 5 ... = 0.47 The class interval becomes 0.47. The problem is this, since the class interval is to be added or subtracted to one, here is how it looks like... calculating class range above 1: 1 + 0.47 = 1.47 (first class above 1) 1.47 + 0.47 = 1.94 1.94 + 2.41 = 2.88 calculating class range below 1: 1 - 0.47 = 0.53 (first class below 1) 0.53 - 0.47 = 0.06 the range: .. 2.88 - 3.58 .. 2.41 - 2.87 .. 1.94 - 2.40 .. 1.47 - 1.93 .. 1.00 - 1.46 .. 0.53 - 0.99 .. 0.06 - 0.52 .. 0.00 - 0.05 this looks very wrong so I asked for an assistance from my adviser and come up with another method. instead of using 0.65 as the lowest value, we used 0. So... Number of desired classes: 5 Class Interval = highest value - lowest value / number of classes ... = 3 - 0 / 5 ... = 0.60 and again we tried to create the ranges still having 1 as the acceptable value or should be the lower limit of the class in the middle. calculating class range above 1: 1 + 0.60 = 1.60 (first class above 1) 1.60 + 0.60 = 2.20 2.20 + 0.60 = 2.80 calculating class range below 1: 1 - 0.60 = 0.40 (first class below 1) outcome :: .. 2.80 - 3.00 .. 2.20 - 2.79 .. 1.60 - 2.19 .. 1.00 - 1.59 .. 0.40 - 0.99 .. 0.00 - 0.39 again, we wanted 5 class ranges but we got six.. and another problem is that the number of classes above 1 should be equal to the number of classes below 1. obviously in the given outcome, it is not.. again, we tried asking another adviser, and it was confirmed that the distribution of the scale should be equal above and below 1. so the next suggestion was to use standard deviation and use it as the class interval.. standard deviation needs to use the mean for all the data sets. there are 800 files in our data sets so i wont bother posting it.. mean for all the data = 1.85 a. ) Summation of the Square of the Difference Between the Calculated Mean and Data = 158.39 b.) Number of Data Sets Minus One (1) = 799 Division of a and b = 0.20 Square Root of 0.20 = 0.45 Creating the ranges... calculating class range above 1: 1 + 0.45 = 1.45 (first class above 1) 1.45 + 0.45 = 1.90 1.90 + 0.45 = 2.35 2.35 + 0.45 = 2.80 2.80 + 0.45 = 3.25 calculating class range below 1: 1 - 0.45 = 0.55 (first class below 1) 0.55 - 0.45 = 0.10 so... .. 2.80 - 3.25 .. 2.35 - 2.79 .. 1.90 - 2.34 .. 1.45 - 1.89 .. 1.00 - 1.44 .. 0.45 - 0.99 .. 0.10 - 0.44 .. 0.00 - 0.09 problems with what i'm doing so far The number of classes should only be 5, but my calculation always exceed There should be equal distribution with the classes above and below 1, but the the number of ranges above 1 is always greater. All interval must be equal (according to my readings) and must not be negative (according to our study) i'm really getting confused as to what to do and what method or statistical tool to use.. can anyone help me with this problem asap.. it would be very helpful..


title: Idea of the Nyblom-Hansen test?

text: The Nyblom-Hansen test gives information about the stability of the estimated parameters in a model. As far as I understand this test, it looks at the score of the ML at evaluates, how near to zero it is. But could anyone give me a more intuitive explanation for this? Especially I am not getting one thing: How does this test evaluate the stability of the parameters by just looking at the score? I imagine that the test has to include different time points of the data, so it needs to take into account the information of time, to evaluate the changing of the parameters over time, but I cannot see how this is connected to the test statistics. So how does the test take the time into account?


title: Confusion related to feature engineering

text: I was reading this tutorial where they mentioned               Make predictions on the development set.     Look at contingency table; where are the errors?     What do you miss? Error analysis!     What feature would the classifier need to get this right? What features are confusing the classifier?     ... If it never appears in the development set, it isn’t useful     ... If it doesn’t appear often, it isn’t useful     Make a contingency table for that feature (should give you good information gain)     Throw it into your classifier (accuracy should improve)      I didn't get how to see the contingency table and see which features are needed for the classifier to get it right? Like lets say I defined some features myself, I can tell that ok these features are not seen most of the time so throw them away. But how do I know which further features it need. Further what is meant by contingency table of a feature


title: clustering gene expression data

text: I have a question about clustering. I' m managing gene expression microarray data and I would like to cluster them in classes. I searched around to find the best clustering algorithm for my data, but since as in all gene expression microarray experiments the number of genes exceed greatly the number of samples being analyzed, the majority of this algorithms assume independency of the variables (so the genes). But in reality we all know that the activity of one gene is not independent from the activity of another gene. My question is: is there an algorithm or a technique able to take into account with simplifications such a dependencies between genes? Thanks a lot! Best e.


title: different results using R or Libre Office calc

text: I have this set of data               x       y       0.5 306.3     1   622.3     2   1230.1     3   2017.7     4   2589     5   3175.4     6   3751.9     9   5585.7     12  7388.8     15  9218.8     20  12042.8     25  15056.5     30  18081.2     45  27095.3     60  35609.2     80  47263.1      When I plot x vs y with R I obtain an almost linear plot, when I plot them with Libre Office Calc I obtain a parabolic plot. How is that possible? Which one is the right one?


title: Logistic regression: can weights be used as a predictor variable?

text: I counted the number of birds in a flock, which gave counts like these:               set.seed(1)     number.birds.in.flock <- sample(51:500, 10)     number.birds.in.flock     # [1] 170 218 307 456 140 450 470 343 329  78      I recorded the species that formed each flock:               species <- c(rep("species.x", 5), rep("species.y", 5))      I then counted the number in the flock that were feeding:               set.seed(1)     number.birds.feeding <- number.birds.in.flock - sample(1:50, 10)      And the number of birds in the flock that were vigilant:               number.birds.vigilant <- number.birds.in.flock - number.birds.feeding      I then went on to model the number of the flock that were feeding. First I converted `number.birds.feeding` to a proportion:               proportion.feeding <- number.birds.feeding/number.birds.in.flock      I used a binomial generalised linear model, with `number.birds.in.flock` as a weight:               glm(proportion.feeding ~ species, weights = number.birds.in.flock, family = "binomial")      Here's my question, can I also `number.birds.in.flock` as both a weight and a explanatory variable like this:               glm(proportion.feeding ~ species + number.birds.in.flock, weights = number.birds.in.flock, family = "binomial")


title: Regression model for periodic data

text: Moved from Stack Overflow Could someone help me to find adequate regression model for my data? I tried to find one by changing the model and initial approximation (ln 15-16) in this simple Python program:               # -*- coding: utf-8 -*-          import matplotlib.pyplot as plt     import numpy as np     import scipy.linalg as la          from numpy import random          from scipy.optimize import minimize               def main():         a = np.loadtxt('group_all_tweets.dat', dtype=np.float32, delimiter='\t')              sin_model = lambda p, x: (p[0] + p[1] * np.exp(np.sin((np.pi / p[2]) * x + p[3]))**3)         x0 = np.array([1.58e+04, -1.72e+03, 24.0, 7.59])              res = minimize(lambda p: la.norm(sin_model(p, a[:, 0]) - a[:, 1]),                        x0=x0,                        method='Powell')         print res              sin_params = res['x']              plt.plot(a[:, 0], a[:, 1])         plt.plot(a[:, 0], sin_model(sin_params, a[:, 0]))              plt.figure()              rss = a[:, 1] - sin_model(sin_params, a[:, 0])              pol_model = lambda p, x: sum([p[i] * x**i for i in xrange(4)])         x0 = np.array([0.0 for i in xrange(4)])         res = minimize(lambda p: la.norm(pol_model(p, a[:, 0]) - rss),                        x0=x0,                        method='Powell')         print res              pol_params = res['x']              rss = a[:, 1] - sin_model(sin_params, a[:, 0]) - pol_model(pol_params, a[:, 0])         plt.plot(a[:, 0], rss)         plt.figure()              plt.plot(a[:, 0], a[:, 1])         plt.plot(a[:, 0], sin_model(sin_params, a[:, 0]) + pol_model(pol_params, a[:, 0]))              plt.show()               if __name__ == '__main__':         main()      Here, firstly I try to find a periodic pattern `((p[0] + p[1] * np.exp(np.sin((np.pi / p[2]) * x + p[3]))**3))` while `p[i]` are varying parameters, then approximate the remains of the first regression with second, polynomial regression. The best result that I managed to get with the method described is shown in in the graph below: ![graph](http://i.stack.imgur.com/bE5A1.png) I'm pleased with how the fit is approaching the bottom part of the graph, but the top parts I just do not like. Has anyone here an experience of finding of regression models? I would be grateful for any help. Thank you. The datafile is here. I need to find a dependence of the second column from the first. I think I want to build a model which contains a periodical component, with "top-trend" and "bottom-trend" components, last two are independent.


title: Good machine learning algo for partial derivatives?

text: Does anyone know of good robust algos to estimate partial derivatives of a regression model? I am talking about a general regression model like this: $\mathbb{E}(y|x_1, x_2, ... x_n) = f(x_1, x_2, ... x_n)$, where I want estimates of $\frac{\partial f}{\partial x_k}$. These are of great importance in medicine, economics, social sciences, etc. Now a linear model would give an approximation of these partial derivatives (its parameters) around the mean point. This approximation can be made valid further around the mean by the inclusion of polynomial terms, but the exponential increase in numbers of parameters as the polynomial degree rises makes this infeasible with a linear model (except perhaps with Lasso?). For better local estimates over the whole input space, natural candidates would be kernel regression and feedforward neural networks, which have smooth forms (unlike tree-based methods) and from which partial derivatives are easy to recover. But I don't know about the robustness of these estimates, as the main aim of these models is to reduce prediction error, not find the right derivatives... I know gradient boosting (and really all regression models) give a "partial dependance", but they are not partial derivatives: they are a function $f_k(x_k)$ that gives the mean value of the predicted regression function $\widehat{f}$ for a given value of $x_k$, averaging over the realisations of the other variables $(x_1, x_2 ... x_n)$ in the training dataset. Any help would be greatly appreciated!


title: Is it possible to determine the maximum predictive accuracy for a given data set with a linear regression? e.g. max adjusted R^2

text: For example say there are N independent variables, and you make a fit with three of them that has a decent Adjusted R2, how do you know when to stop? This is a theoretical question.


title: Fitting ARMA(p,q) by conditional maximum likelihood

text: I have found a lot of information on fitting ARMA(p,q) using exact maximum likelihood, but very little for conditional maximum likelihood. Can anyone direct me to an R function that fits ARMA(p,q) by conditional maximum likelihood as well as some documentation that explains the concept?


title: Durbin-Watson test for first order autocorrelation when interaction effects are present

text: I have a multiple linear regression model (OLS) that has a constant, 5 variables, and an interaction term between two of the dummy variables. I need to perform the Durbin-Watson test for first-order correlation, and in order to find the critical values, I need to know the number of explanatory variables (not taking into account the constant). Is the interaction term an explanatory variable in itself for these purposes (so I have 6 explanatory vars), or do I find the critical values for 5 explanatory vars? Thanks!


title: Applying an RBF kernel first and then train using a Linear Classifier

text: I will start off by saying that I don't have a concrete understanding of whats under the hood of a SVM classifier. I am interested in using an SVM with the RBF kernel to train a two class classifier. I however find that the training (and even prediction) takes a lot of time when working with the RBF kernel (have been implementing on matlab using libsvm and python using sklearn). My question is that is it possible to project my data into the higher dimension using an RBF kernel on its own and then apply a linear SVM to this transformed data, i.e. is that going to yield the same results as using an RBF SVM, as long as I use the same C and Gamma? I am not too sure how kernels are applied so I hope this part makes sense. If that is true, that I could pre-process the data into the higher dimensional feature space using the RBF kernel and then the training and prediction much fast by using a simple linear classifier.


title: Relationship between the kernel and the value of C in SVM's

text: How exactly does the value of C relate across different kernels that we can use for SVM's? As in, how does it vary when changing the polynomial degree of a kernel or while using a Gaussian kernel?


title: Reducing data size to cross correlate with another data set?

text: I have three matrices A, B, C. A is a matrix of 200 X 32. There are 2000 such different A matrices which make up the B matrix. B is a matrix of 2000 x A. That is there are 2000 x 200 rows in matrix B with 32 columns. C is a matrix of 2000 x 14. Each of the A matrix corresponds to a row in C matrix. Let's say matrix A1 of 200 X 32 corresponds to first row in C matrix. I want to see the correlation among the matrix A1 and the row of C. One way to do it is to have a one to one mapping for each of the A matrices to C's rows. But I'm not sure how to shrink A to a vector which can be compared one to one. Is this feasible and possible to do? If yes how? Are there some other approaches to dealing with this kind of data? Corelation Details: The row of C has 14 attributes, these 14 attribute values is because of the data in one of the A matrices. Think of it like this, there is a process which takes the input matrix A1 of size 200 x 32 and produces one row of C with 14 attributes. These 14 attributes might have large values depending on the input and low due to input. I don't know what output I might get. It might be a direct relation or it might be inverse.


title: Simple OLS with two samples

text: I want to obtain an unbiased estimator for $b_1$ in a simple regression like that: $Y_i = B_0 + B_1X_i + u_i$ when I have two samples, always same size for Y and X, but once the sample size is l and once the sample size is m. The respective sample means $\bar{Y_l},\bar{X_l}$ and $\bar{Y_m},\bar{X_m}$ are given. Now I wonder how I cans tart to get an unbiased estimator? My idea was to use the 'normal/one-sample' formula and just put weights (correcting for different sample size between the two independent sets of data) in front. An estimator for $b_1$ would be: (X'X)$^{-1}$X'Y without matrices: $\frac{\sum X_iY_i - N \bar{Y}\bar{X}}{\sum X_i^2 -N \bar{X}^2}$ which I wanted to modify to $\frac{l}{m+l} \frac{\sum X_iY_i - L \bar{Y_l}\bar{X_l}}{\sum X_i^2 - L \bar{X_l}^2} + \frac{m}{m+l} \frac{\sum X_iY_i - M \bar{Y_m}\bar{X_m}}{\sum X_i^2 - M \bar{X_m}^2}$ The capital M and L denoting the respective sample size. Now I am not sure if my result is right, as I cannot show if it is unbiased, to be honest. Is it unbiased in probabilistic terms? Or is it just a wrong estimator?


title: Chi Squared Results in R and Python

text: Based on this answer, Python requires expected values in a chi square test to be absolute frequencies. Consider the following in Python:               import numpy     import scipy.stats     # chisquare function requires (observed, expected)     scipy.stats.chisquare(numpy.array([0,0,23,0]), numpy.array([1,1,1794,1]))     (1751.2948717948718, 0.0)      results in a p-value of 0 (whatever that means). The same calculation in R, which requires that the expected values be proprotions:               chisq.test(c(0, 0, 23, 0), p=c(1/1797,1/1797,1794/1797, 1/1797))                  Chi-squared test for given probabilities          data:  c(0, 0, 23, 0)     X-squared = 0.0385, df = 3, p-value = 0.998      resulting in a p-value of 0.998. Which is correct?


title: How to interpret negative ACF (autocorrelation function)?

text: ![Returns](http://i.stack.imgur.com/vvOCA.jpg) So I plotted the ACF/PACF of oil returns and was expecting to see some positive autocorrelation but to my surprise I only get negative significant autocorrelation. How should I interpret the above graph? They seem to indicate that there is a tendency for oil returns to increase when it decreased previously and vice-versa, thus the oscillating behaviour. Please correct me if I'm wrong.


title: Test to compare large and small datasets

text: I send Alice and Bob out to record people's eye colour (blue, brown, green, other). Alice does a great job and writes down the eye colour of 2000 people. Bob only records the eye colour of 20 people. Alice found all 2000 people had brown eyes. Bob found 19 people with brown eyes and one with green eyes. How important is the one green result? How can I compare the results? Is the green result significant? **Update** I now send Alice and Bob out to a different area on another day. Alice is thorough and finds 2000 people again, while Bob is still lazy and only finds 20 people. Their results would look something like this:               | Day | Dataset | Blue | Brown | Green | Other |     |  1  | Alice   | 0    | 2000  | 0     | 0     |     |  1  | Bob     | 0    | 19    | 1     | 0     |     |  2  | Alice   | 10   | 1900  | 45    | 45    |     |  2  | Bob     | 2    | 10    | 6     | 6     |      I repeat this process, sending them out to different areas on different days. Alice always finds 2000 people and Bob always finds 20. Considering each day separately, how can I determine if Alice or Bob's results are more representative of the population they sampled on a given day? How can I compare their results for a given day?


title: Do I need to aggregate data when using linear mixed effects model?

text: I have a 2 by 2 design with 12 subjects. The two factors are within-subjects variables. For each cell, 24 responses were collected from each subject. My question is whether I need to aggregate every 24 responses into a mean and then use the mean in lmer as the dependent variable. This webpage http://talklab.psy.gla.ac.uk/simgen/faq.html seems to suggest that aggregation is not needed if use mixed effects model, but I doubt that it is true. Also, if the responses are binary responses, can I just average them and use generalized lmer with binomial family? Thanks.


title: Probability that multiple confidence intervals contain the true population mean

text: Suppose you take 3 samples and report a 99% confidence interval for each dataset. How would you go about calculating the probability that all the datasets contain the true population mean? No other statistics are provided.


title: How to define weights when testing an exponential trend?

text: In a linear model including ANOVA one can test a trend (e.g., linearity, quadratic effect, etc.) among the ordered effects (regression coefficients or factor levels) through assigning proper weights. For example, suppose that we have a time series regression model $Y=\alpha+\beta_1X_1+\beta_2X_2+\beta_3X_3+\beta_4X_4+\epsilon$ where $Y$ is a vector of length $n$, $X_i$ is an explanatory variable of length $n$ ($i=1,2,3,4$), and $\epsilon$ is the residual vector. The linear and quadratic trends among four equally-spaced (and ordered) effects $\beta_1, \beta_2, \beta_3, \beta_4$ can be tested with weights of -3, -1, 1, 3, and 1, -1, -1, 1 respectively. When those effects are not equally-spaced, such weights can be obtained through the coefficients of orthogonal polynomials. My understanding is that typically the exponential trend among some sequential effects can be tested through log-transforming the response variable. However, the situation I'm dealing is not directly about the response variable, but the effect estimates (or regression coefficients, e.g., 30 ordered effects similar to a linear model above). So the challenge I'm facing here is, how to find the weights when testing for an exponential trend, when the effects are equally- or unequally-spaced? P.S. There are two time scales involved here that may cause some confusions. At the local time level, each explanatory variable $X_i$ in the time series regression model codes for a causal event that lasts for a short period of time. At the higher level of time, there is the sequence of those causal events and their estimated effects $\beta_1$, $\beta_2$, .... It is the exponential trend or saturation effect among those estimated effects that is of interest to me. For a linear trend, I could perform a general linear hypothesis testing with appropriate weights for those effects. However, I'm not so sure whether it's feasible to do the same for an exponential trend (or saturation effect).


title: Goodness of fit test for exponential distribution

text: The given data below shows the time in intervals in seconds between successive white cars in flowing traffic in an open road.Can these be modeled by an exponential distribution.               Time        0-     20-      40-     60-      90-     120-180     Frequency   41     19       16      13        9        2      My question is when I calculate Expected frequencies though the total should add up to 100 it does not do so because of **rounding off** errors right?   I get a value of 98.89.   Here should I need to add another category myself as the interval being **180-infinity** and get its expected value as 100-98.89=1.11. Is it necessary to add this category as this is an **exponential distribution** which goes for infinity?   Of course as the Expected value of the created category is less than 5, in this case I have to add this up with the upper category.   But if I do not consider this new category then the degrees of freedom changes.Is it necessary to add this new category as 180-infinity


title: Regression for outcomes truncated from both sides

text: I have an outcome which is truncated from both left and right sides. I would like to know which regression method may account this kind of analysis. Here is an example:               set.seed(123)     z<-rnorm(300)     y<-z[z>=-0.8 & z<=1.2] #truncated     x1<-sample(1:3, length(y), T)     x2<-rbinom(length(y), 1, 0.3)      What's the best regression method for modeling `y~x1+x2`? The R package `truncreg` seems just have arguments for either left or right truncation rather than both sides like this example.


title: Help with the Ljung-Box test for time independence of residuals

text: I fit a simple linear model $y = bX$ to a data set today, and that produced 24 residuals (I have 24 data points, one for each year from 1984-2007). I would like to test the time-independence of the residuals of my model, and I was recommended by my supervisor to use the Ljung-Box test. The `Box.test` function in R takes 4 arguments:   * `x`: a numeric vector or univariate time series.   * `lag`: the statistic will be based on lag autocorrelation coefficients.   * `type`: test to be performed: partial matching is used.   * `fitdf`: number of degrees of freedom to be subtracted if x is a series of residuals. What does `lag` mean, and what value would you guys recommend I use for the test? Also, what does `fitdf` represent, and what would the value for that parameter be in my case? Finally, the value of `x` is a vector of my 24 residuals, correct?


title: How to simulate Signal-Noise Ratio?

text: I have a vector of simulated data:               data = c(0.47, 0.45, 0.30, 1.15, 0.82, 0.38, 0.51, 1.36, 1.72, 0.36)      I've been adding noise to this by generating random numbers centered at 0 with different standard deviations:               noise = rnorm(10, mean = 0, sd = 0.1)     data_wNoise = data + noise      I've been setting the standard deviation arbitrarily (between 0.001 and 1.5). Is there a better way to simulate this by setting a specific signal-to-noise ratio? I don't know anything about the power of the signal in the data.


title: Clustering 5 dimensions into 4 based on mean?

text: In my project, one variable named Organizational Culture has 5 dimensions namely Employee Development, Harmony, Customer Orientation, Social Responsibility and Innovation. In SPSS, I need to configure these 5 dimensions into 4 types so that:   1. The highly integrative culture has high scores on all 5 dimensions.   2. The market oriented culture describes the organizations that emphasize customer orientation more than any other cultural dimension.   3. The moderately integrative culture shows average scores on all ?ve culture dimensions.   4. Lastly, hierarchy culture show low score on all culture dimensions. In short, here is what I need: I need to cluster the firms from which I got data into 4 types (highly integrative, moderately integrative and so on) on the basis of the score assigned by their (firms') managers to the the 5 culture dimensions. As a result I will have four types of firms which will have four different culture types. I do see some hope in the method you suggest in your above post. Thank you very much for your time and I appreciate any help :)


title: Standard error for the mean of a sample of binomial random variables

text: Suppose I'm running an experiment that can have 2 outcomes, and I'm assuming that the underlying "true" distribution of the 2 outcomes is a binomial distribution with parameters $n$ and $p$: ${\rm Binomial}(n, p)$. I can compute the standard error, $SE_X = \frac{\sigma_X}{\sqrt{n}}$, from the form of the variance of ${\rm Binomial}(n, p)$: $$ \sigma^{2}_{X} = npq$$ where $q = 1-p$. So, $\sigma_X=\sqrt{npq}$. For the standard error I get: $SE_X=\sqrt{pq}$, but I've seen somewhere that $SE_X = \sqrt{\frac{pq}{n}}$. What did I do wrong?


title: Bayesian and frequentist optimization and intervals

text: I realize the methodology pursued by the Frequentist and Bayesian camps generally differ. However, one method of estimation that they do share is optimization of a certain function:   * Frequentists maximize the _likelihood_ function, giving the Max. Likelihood (ML) estimator.   * Bayesians maximize the _posterior_ function, giving the Max A-Posteriori (MAP) estimator. Both functions will typically have been constructed using Baye's rule/theorem, which is universally agreed upon, and which might have been applied once (in "batch mode") or multiple times iteratively. Similarly, both Frequentists and Bayesians will deduce their interval (confidence/credibility) from this function. So if the prior is _uninformative_ (assuming we can formulate such a prior), there should be no distinction between the "results" obtained by Bayesians and Frequentists, even though the interpretation of said results will be different. If this is right, then the only practical difference between Bayesians and Frequentists is the prior. Is this true? * * * Edit: Actually, the optimization bit of my question is a bit misleading, as it is only a specific example of differences between Bayesian and Frequentist thinking. My question could be posed simply as the difference simply between the likelihood function and the posterior. For example, would frequentists ever use MCMC to calculate the likelihood function?