title: Required conditions for using a t-test

text: The conditions that I have learned are as follows:   1. If the sample size less than 15 a t-test is permissible if the sample is roughly symmetric, single peak, and has no outliers.   2. If the sample size at least 15 a t-test can be used omitting presence of outliers or strong skewness.   3. With a larger sample the t-test can be use even if skewed distribution if the sample is greater than 30, but less than 10% of the population.  Why can't you use a the t-test when the sample size is larger than 10% of the population size? What happens then? Do you use the z-test?


title: Finding the effects of certain levels of a factor predictor

text: I have fitted a binary classification `gbm` model, and one of the predictor variables, `Affiliate` has 50 different levels. Given the following readout from `gbm.perf`:               var                             Importance     Affiliate                       52.939994     ProgramName                     25.765384     Distributor                     17.502216     MarketingCategoryName           2.212374     ProductName                     1.580031      I am now wondering which levels within `Affiliate` are associated with more `1` classifications and which are associated with more `0` classifications. I was thinking about fitting a `glm(..., binomial("probit")` model by each predictor to give me a general sense, but I think this is quite a rough way to go about it. Where can I find some writings on how to get a more granular understanding of the effects of the different variables?


title: minimum abundance of a bacterial species

text: I conducted a study in which I used sequencing of marker genes to study bacterial diversity and I would like to calculate the minimum abundance for a species to still be detected. The **total volume** of my experimental unit is **2000 ml**. My **sample size** is **150 ml**. The **total abundance** (all species confounded) is **10^6 cells/ml**. I can identify **4000 individuals in each sample** (4000 reads). As I understand it, I sampled **4000 individuals out of 150 ml** which is about **27 (= 4000/150) ind./ml**. With a **total abundance of 10^6 ind./ml** any given species needs to have a **concentration of at least 10^6/27 = 3.7*10^4 ind./ml** to be represented (in average) with at least 1 ind. in my sample. Is that logic correct or do I need to take into account my total experimental volume, as having sampled **4000 ind. out of 2000 ml** , hence only **2 ind./ml** with an according detection limit of **5*10^5 ind./ml**. sorry for that really simple question, but I get very confused and can't make up my mind. Thanks for your help! Fabian


How to transform a normal random variable such that I can simulate normal samples between the range of 1 and 45? Do I need to do a jacobian transformation ?


title: Inter-rater agreement with multiple raters

text: **Compare inter-rater reliability between 8 human experts and a computer** I am running a study where I use automated analysis to analyse an EEG signal (brain activity) and quantify certain properties thereof. The computerized method quantifies each property and places it into one of three categories, for example: normal, moderate or abnormal (for some properties, two categories exist). My goal is to compare the automated method with expert reviewers (i.e., humans), and for this I have sent the same set of recordings to 8 neurologists and asked them to categorize each of the properties for all the recordings. Now my question: I want to compare the reliability of automated analysis with human reviewing, but because there is no gold-standard, I want to do this by comparing inter-rater variability between expert reviewers using the kappa statistic, and then between the system and reviewers by using the kappa statistic again. How exactly should I do this? My first idea was to calculate the inter-rater agreements between each of reviewers individually (i.e., A vs B, A vs C, etc.) and find the mean thereof, and then do the same for computer analysis by taking the mean of its agreement with each of the reviewers. But this seems wrong to me. Will this make a fair comparison between human experts and a computer algorithm? I found somewhere on Google that similar measures are needed for comparing algorithms with humans in communication studies such as interpreting texts, and although it seems like they also use kappa values, its not clear to me how fair comparisons are made between the computer and the human.


title: How to use chi-squared (or Bayesian) to detect non-chance variations in matrix of frequencies?

text: I've got a database showing signatures added to a petition, with each signature linked to a location (city, state), a date, and a sequence in the date (e.g. the seventh signature on 3/26 came from Rochester, NY). I'd like to detect geographic patterns in the data - e.g. when one signature from Rochester leads to further signatures from Rochester, and also when a signature from Rochester leads to a signature from San Francisco or Fairbanks). To simplify for the moment, I can look just a signatures by day, ignoring sequence w/i day and construct a matrix - one row for each location, one column for each location. The value of each cell is the observed frequency of one or more signatures from the column location on a day, when there was at least one signature from the row location on the previous day. Clearly, if there is no geographic linkage, the values for each cell in a given column will be identical. Clearly, also, chance variations will lead to slight variations in cell values down the column, even if there is no geographic linkage. What statistic should I use to determine whether, in general, there is some geographic dependency - that it helps to know which locations provided signatures yesterday in order to predict which locations will provide signatures today? (I.e. my null hypothesis is that knowing which locations generated signatures yesterday doesn't help to predict which locations will generate signatures today.) What statistic should I use to detect specific dependencies? In particular, I'm guessing that some locations provide signatures at random, while other locations are driven either by activity within that location (e.g. a bunch of activists in Fairbanks, AK work every day to round up signatures), or even between locations (San Francisco activists press for lots of SF signatures on Monday, then contact their colleagues in LA on Tuesday, etc.). In this case, I have a null hypothesis for each location X: Knowing which locations (incl. X) generated signatures yesterday doesn't help to predict whether X in particular will generate signatures today. From my long-ago stats course, I'm thinking chi-squared, but I'd be grateful for more detailed pointers.


title: Problem with negative values when simulating multivariate data using rmvnorm

text: I'm simulating a multivariate dataset using a modified form of the `mvrnorm()` function (from the MASS package in R). The problem is that I'm getting some negative values after the eigenvalue transformation (I assume) because I have many negative correlations and a number of small means. Is there a special way to deal with this phenomenon? My idea is to just add some factor to all the datapoints (such as the 1st quartile) because I don't care so much what the exact means are, as long as the correlation structure remains intact. A simple example:               require(MASS)     # For mvrnorm()     require(Matrix) # For nearPD()          corr <- diag(5)     corr[5,1] <- .5     corr[1,5] <- .5     corr[4,2] <- -.5     corr[2,4] <- -.5          set.seed(1000)     mm <- mvrnorm(n=10, mu=rep(1,5), Sigma=nearPD(corr, corr=TRUE)$mat,               empirical=TRUE)       As you can see, mm has nonpositive examples. Since I want to model physical measurements, this makes no sense. edit2: multivariate sampling from log-normal distribution still results in negative examples:               mvrlnorm <- function (n = 1, mu, Sigma, tol = 1e-06, empirical = TRUE) {         require(Matrix)         p <- length(mu)         if (!all(dim(Sigma) == c(p, p)))             stop("incompatible arguments")         eS <- eigen(Sigma, symmetric = TRUE, EISPACK = TRUE)         ev <- eS$values         if (!all(ev >= -tol * abs(ev[1L])))              stop("'Sigma' is not positive definite")           # HERE BE log-normal distribution         X <- matrix(rlnorm(p * n), nrow=n)         if (empirical) {             X <- scale(X, TRUE , FALSE)             X <- X %*% svd(X, nu = nrow(X), nv = ncol(X))$v             X <- scale(X, FALSE, TRUE)         }         retval <- eS$vectors %*% diag(sqrt(ev), length(ev)) %*% t(eS$vectors)         retval <- X %*% retval         retval <- sweep(retval, 2, mu, "+")         X <- retval                   nm <- names(mu)         if (is.null(nm) && !is.null(dn <- dimnames(Sigma)))              nm <- dn[[1L]]         dimnames(X) <- list(nm, NULL)         if (n == 1)              drop(X)         else t(X)     }


title: Can anybody help me in calculating the relative importance of predictors in a model using the t-value?

text: I have a model with three variables(no intercept).The coeeficients are -2.7398,-.0190,.3044 respectively.The relative importance calculated from t value are 16%,3%,3%.Can anyone help me in finding out these numbers?


title: Chi-squared fitting with uncertainties estimated from the data

text: I have a question about how to compute the goodness of a fit. Suppose I have a bunch of data $(x_i, y_i, \sigma_i)$ and that I want to fit a non-linear model $f$ (a stretched exponential) to them by minimizing the $\chi^2$: $$\chi^2 = \sum_{i=1}^N \frac{\left(y_i - f\left(x_i, \beta\right)\right)^2}{\sigma_i ^ 2}$$ At the minimum, I know that the above follows the $\chi^2$ distribution with $N - \nu$ degrees of freedom, where $\nu$ is the number of parameters of the model. Thus I could compute a $p$-value based on that and evaluate the goodness of my fit. This requires however that the uncertainties $\sigma_i$ are known. Now suppose that my observations of $y_i$ are actually obtained by averaging over samples of different sizes. I mean that for each $i$ I have $n_i$ replicates $y_{ij}$, $j = 1, ..., n_i$, and that, instead of $y_i$, I have $$\overline{y}_i = \frac{1}{N}\sum_{j=1}^{n_i} y_{ij}$$ and instead of the true $\sigma_i$ I have $$\hat\sigma_i = \frac{s_i}{\sqrt{n_i}}$$ where $s$ is the unbiased estimator of the standard deviation of the $i$-th sample. Would using $\hat\sigma^2$ in place of $\sigma^2$ violate the assumptions behind the usage of $\chi^2$ for assessing the GoF?


title: 2x2 contingency table - maximum likelihood estimate of the odds ratio and exact confidence intervals

text: I have two questions regarding the following 2x2 contingency table: ![enter image description here](http://i.stack.imgur.com/Sy8fk.png) 1, How can I derive the maximum likelihood estimate of the odds ratio ($OR = (a*d)/(c*b)$) 2, How can I derive exact confidence intervals (not wald normal approximations) It is my understanding that a distributional assumption must be made, i would have thought this would be binomial, but am not sure where to go from there.


title: Markov Chain Queuing Problem

text: Suppose the times between successive arrivals of customers at a single-server station are independent random variables having common distribution $F$. Suppose when a customer arrives, he or she either immediately enters service if the server is free or else joins the end of the waiting line if the server is busy with another customer. When the server completes work on a customer that customer leaves the system and the next waiting customer, if there are any, enters service. Let $X_{n}$ denote the number of customers in the system immediately before the $n^{th}$ arrival and $Y_{n}$ denote the number of customers that remain in the system when the $n^{th}$ customer departs. The successive service times of customers are independent random variables (which are also independent of the inter-arrival times) having common distribution G.   1. If $F$ is the exponential distribution with rate $\lambda$, which, if any, of the processes $\\{X_{n}\\}, \\{Y_{n}\\}$ is a Markov chain?   2. If $G$ is the exponential distribution with rate $\lambda$, which, if any, of the processes $\\{X_{n}\\}, \\{Y_{n}\\}$ is a Markov chain?   3. Give the transition probabilities of any Markov chains in parts (a) and (b) I have solutions to the problem but I have read over them multiple times and still don't understand the working. Any help would be greatly appreciated. Cheers!


title: Variance explained by random effects using lme4

text: I am using the `glmer()` function from the lme4 package to run a GLMM using the poisson distribution. In all the examples that I see, the random effects part of the output has a residual part that has been estimated from the data (surrounded by 2 asterisks on either side in the example below). This information can then be used in interpreting the amount of variation explained by the random effect. Here is an example:               > summary(M1)     Linear mixed model fit by REML      Formula: Richness ~ NAP * fExp + (1 | fBeach)      Data: RIKZ      AIC   BIC    logLik   deviance REMLdev     236.5 247.3 -112.2    230.3    224.5     Random effects:     Groups   Name          Variance Std.Dev.     fBeach   (Intercept)   3.3072   1.8186       **Residual             8.6605   2.9429**     Number of obs: 45, groups: fBeach, 9          Fixed effects:                   Estimate Std. Error t value     (Intercept)   8.8611     1.0208   8.681     NAP          -3.4637     0.6279  -5.517     fExp11       -5.2556     1.5451  -3.401     NAP:fExp11    2.0005     0.9461   2.114          Correlation of Fixed Effects:                (Intr) NAP    fExp11     NAP        -0.181                   fExp11     -0.661  0.120            NAP:fExp11  0.120 -0.664 -0.221      However, when I use my own data, I get output that does not include this information, and I am not sure why. I want to know how much variation is explained by my random effects, but can't figure out how to access the information necessary to answer the question. Any clues? Is this a data/statistics issue or is this a knowing how to access the information issue? I apologize if I'm asking in the wrong place. The output I get looks similar to the following output:               Generalized linear mixed model fit by the Laplace approximation      Formula: y ~ z.score(x1) + z.score(x2) + z.score(x3) + z.score(x4) + z.score(x5) +      z.score(x6) + (1 | RE)      Data: p      AIC   BIC logLik deviance     419.5 454.7 -201.8    403.5     Random effects:     Groups Name        Variance Std.Dev.     RE     (Intercept) 0.021605 0.14699      Number of obs: 600, groups: RE, 40          Fixed effects:                       Estimate   Std. Error z value Pr(>|z|)         (Intercept)       1.70591    0.02911    58.60   < 2e-16 ***     z.score(x1)       0.19087    0.03595    5.31    1.10e-07 ***     z.score(x2)      -0.14302    0.04083   -3.50    0.000460 ***     z.score(x3)      -0.16562    0.04020   -4.12    3.79e-05 ***     z.score(x4)       0.13229    0.03425    3.86    0.000112 ***     z.score(x5)      -0.10588    0.03985   -2.66    0.007885 **      z.score(x6)       0.17600    0.05798    3.04    0.002401 **      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Correlation of Fixed Effects:                 (Intr) z.(x1) z.(x2 z.s(x3) z.(x4 z.(x5     z.scr(x1)  -0.051                                        z.s(x2)     0.038  0.259                                 z.scr(x3)   0.045  0.156  0.113                          z.(x4      -0.040  0.144 -0.052  0.044                   z.(x5       0.026 -0.368 -0.339 -0.072 -0.073            z.scor(x6) -0.031 -0.020  0.002 -0.143 -0.004  0.004      Here is some sample data, to be fit with `glmer(y ~ x1 + (1|RE), data=d, family=poisson)`.               d <- data.frame(       y  =  c(3, 5, 2, 6, 3, 7, 2, 3, 0, 4, 0,10, 1, 4, 0, 4, 2, 3, 0, 6,                3, 4, 2, 3, 2, 3, 3, 4, 0, 5, 6, 5, 4, 4, 0, 3, 1, 6, 0, 3, 2,                2, 1, 6, 2, 7, 0, 2, 0, 4, 0, 6, 4, 5, 1, 5, 1, 4, 1, 2, 3, 6,                6, 7, 0, 5, 0, 9, 1, 4, 5, 6, 1, 7, 1, 4, 1, 4, 0, 4, 1, 6, 1,                4, 0, 7, 1, 4, 0, 6, 0, 7, 2, 6, 0, 6, 1, 5, 0, 4, 1, 7, 2, 4,                1, 5, 1, 7, 2, 5, 0, 4, 3, 5, 1, 4, 0, 3, 0, 6, 0, 8, 3, 9, 0,                2, 3, 8, 0, 1, 0, 3, 0, 5, 0, 4, 4, 5, 0, 5, 1, 5, 3, 5, 1, 4,                3, 4, 4, 4, 4, 4, 4, 7, 1, 8, 1, 4, 0, 2, 2, 5, 1, 4, 1, 5, 1,                4, 2, 4, 2, 4, 0, 6, 1, 6, 0, 6, 1, 2, 1, 3, 1, 8, 1, 6, 1, 6,                0, 6, 1, 6, 2, 6, 2, 4, 0, 1, 1, 1, 1, 6, 5, 5, 1, 5, 2, 4, 2,                6, 1, 7, 1, 8, 2, 8, 1, 8, 2, 4, 1, 7, 3, 6, 4, 7, 3, 7, 1, 6,                3, 5, 1,10, 1, 7, 2, 5, 1, 5, 0, 6, 1, 8, 4, 7, 1, 6, 1, 9,                0, 9, 1, 3, 2, 5, 2, 9, 3, 5, 0, 2, 2, 3, 0, 5, 0, 5, 0, 4, 3,                6, 1,10, 2, 8, 0, 6, 0, 4, 2, 6, 2, 4, 2, 6, 1, 4, 0, 5, 2,                6, 1, 5, 2, 5, 1, 5, 1, 5),       x1 = rep(c(0.1008, 0.0511, 0.1792, 1.0345), c(80, 80, 80, 60)),       RE = rep(c(37, 88, 139, 190, 241, 292, 343, 394, 91, 142, 193, 244, 295,                   346, 397, 40, 94, 145, 43, 196, 247, 298, 349, 400, 301, 352,                   403, 250, 148, 199, 46, 97, 355, 406, 253, 304, 49, 100, 151,                   202, 37, 88, 139, 190, 241, 292, 343, 394, 91, 142, 193, 244,                   295, 346, 397, 40, 43, 94, 145, 247, 298, 349, 196, 400, 199,                   250, 301, 352, 406, 46, 97, 148, 403, 49, 100, 151, 202, 253,                   304, 355, 37, 88, 139, 190, 241, 292, 343, 394, 193, 244, 346,                   397, 295, 40, 91, 142, 43, 94, 145, 196, 46, 97, 148, 151, 247,                   400, 298, 349, 352, 199, 250, 301, 403, 253, 304, 355, 202, 406,                   49, 100, 37, 88, 139, 190, 241, 292, 343, 394, 346, 397, 193,                   244, 295, 40, 91, 142, 43, 94, 145, 196, 247, 298, 349, 400,                   97, 148, 46, 199, 250, 301), each=2)     )


title: Rectifier function of random variable

text: Let $X$ be a random variable with distribution $F_X$ and density $f_X$. Define $$g(x) = \left\\{ \begin{array}{lr} x & : x \ge 0\\\ 0 & : x < 0 \end{array} \right\\}$$ and let $Y=g(X)$. What is the PDF of $Y$? **My approach:** via CDF $$F_Y(y) = \Pr(Y\le y) = \Pr(g(X)\le y) = \Pr(X\le y) = F_X(y), \ y\ge 0$$ After taking the derivative, $$f_Y(y) = f_X(y), \ y \ge 0$$ This seems to be incorrect because $f_Y(y)$ is not a PDF (integral is not equal to one). I know that it probably has to do something with the fact that all $x\lt 0$ are mapped onto $y=0$. I just don't know how I can incorporate this in a constructive way.


title: Descriptive statistics show more than 100%

text: I have race as a variable in my model and my measurement reflects the executives' race in the company's management (score 1 if any specific race is in the management). The problem is that managements of companies have more than one race in their executives. So when I want to statistically summarize race which are four, the total percentage is more than 100% What is the best way to demonstrate the results? Your suggestions are highly appreciated


title: Logistic regression power analysis with moderation between categorical and continuous variable

text:
I've been reading Simulation of Logistic Regression Power Analysis - Designed Experiments, http://sas-and-r.blogspot.com/2009/06/example-72-simulate-data- from-logistic.html, and Power analysis for ordinal logistic regression and I'm still a little lost on how to do a power analysis for my data. I want to be able to determine what N I should have if I have an interaction between a categorical variable (with 3 levels) and a continuous variable. Simulation of Logistic Regression Power Analysis - Designed Experiments provides some information, but I can't figure out how to simulate the relationship between the categorical and continuous variables and outcome. > set.seed(1)   > > > repetitions = 1000   >  N = 10000   >  n = N/8   >  var1 = c( .03, .03, .03, .03, .06, .06, .09, .09)   >  var2 = c( 0, 0, 0, 1, 0, 1, 0, 1)   >  rates = c(0.0025, 0.0025, 0.0025, 0.00395, 0.003, 0.0042, 0.0035, 0.002)   > > > var1 = rep(var1, times=n)   >  var2 = rep(var2, times=n)   >  var12 = var1**2   >  var1x2 = var1 *var2   >  var12x2 = var12*var2   > > > significant = matrix(nrow=repetitions, ncol=7)   > > > startT = proc.time()[3]   >  for(i in 1:repetitions){   >  responses = rbinom(n=N, size=1, prob=rates)   >  model = glm(responses~var1+var2+var12+var1x2+var12x2,   >  family=binomial(link="logit"))   >  significant[i,1:5] = (summary(model)$coefficients[2:6,4]<.05)
 > > significant[i,6] = sum(significant[i,1:5])
 > modelDev = > model$null.deviance-model$deviance   >  significant[i,7] = (1-pchisq(modelDev, 5))<.05   >  }   >  endT = proc.time()[3]   >  endT-startT   > > > sum(significant[,1])/repetitions # pre-specified effect power for var1   >  [1] 0.042   >  sum(significant[,2])/repetitions # pre-specified effect power for var2   >  [1] 0.017   >  sum(significant[,3])/repetitions # pre-specified effect power for var12   >  [1] 0.035   >  sum(significant[,4])/repetitions # pre-specified effect power for var1X2   >  [1] 0.019   >  sum(significant[,5])/repetitions # pre-specified effect power for var12X2   >  [1] 0.022   >  sum(significant[,7])/repetitions # power for likelihood ratio test of model   >  [1] 0.168   >  sum(significant[,6]==5)/repetitions # all effects power   >  [1] 0.001   >  sum(significant[,6]>0)/repetitions # any effect power   >  [1] 0.065   >  sum(significant[,4]&significant[,5])/repetitions # power for interaction > terms   >  [1] 0.017   > I feel like I should be able to adapt the code from Power analysis for ordinal logistic regression and that this would be a better, more succinct option > library(rms) > > tmpfun <\- function(n, beta0, beta1, beta2) {   >  x <\- runif(n, 0, 10)   >  eta1 <\- beta0 + beta1*x   >  eta2 <\- eta1 + beta2   >  p1 <\- exp(eta1)/(1+exp(eta1))   >  p2 <\- exp(eta2)/(1+exp(eta2))   >  tmp <\- runif(n)   >  y <\- (tmp < p1) + (tmp < p2)   >  fit <\- lrm(y~x)   >  fit$stats[5]   >  }   > > > out <\- replicate(1000, tmpfun(100, -1/2, 1/4, 1/4))   >  mean( out < 0.05 )   > but I'm not completely sure how to do so. I'm assuming tmpfun(100,-1/2, 1/4,1/4) is specifying the N and betas that you want, but how do I adjust tmpfun to another (categorical) variable and include an interaction term? Ultimately the equation should include 6 betas: the intercept, the beta for x, the beta for z1, the beta for z2, the interaction term between x and z1, and the interaction term between x and z2. Finally, I can't find any reliable sources on what sorts of "effect sizes" I should be using as small or medium. Let me know if I can provide more information!


title: Ok to use different levels of aggregation for different prediction periods?

text: Using regression. For simplicity, let's assume stationary non-seasonal data. Example: If we wish to predict sales volume for specific months, we aggregate daily data to monthly data and fit our model, etc. If we also want to predict by year, would it be valid to then aggregate that data into years, fit a model and predict? Or should we stick with the monthly aggregates/model? Then we'd have to predict for 12 months at a time. We do not wish to compare or relate monthly and yearly predictions.


title: Multi-level model with varying intercept vs. fixed effect regression

text: I'm working through Gelman and Hill, Data Analysis and Regression using Multilevel/Hierarchical Models (2007), using the `arm` package, and trying to relate multilevel models to the econometric framework I'm more familiar with. I expected a multilevel model with a non-varying slope coefficient and a varying intercept coefficient to provide identical results to a fixed effect regression with no constant. I expected the following R and Stata code to produce the same results. They do not - can you tell me why? R code:               M1 <- lmer(y ~ x1 + x2 + (1 | county))      Stata code:               reg y x1 x2 i.county, noconstant      The coefficients produced by these two approaches are quite different. The Stata code regresses y on x1, x2 and K additional indicator variables for each county. What is R doing that is different? Is there an OLS regression analog?


title: Correlation between variables and outcome

text: It is a retrospective study consisting of 47 horses with a specific disease. I have no data on the control group, just knowing that during this period of time the clinic had 382 horses where 47 of those were diagnosed with the current disease. I was asked to look at correlation between three variables and the diseas (breed, time of year - Jan, Feb, etc and an oral stereotypic behaviour known as cribbing in the horse). But how can I see if the variables correlates with the disease when I don't have any control group?


title: Robust monotonic regression in R

text: I have the following table in `R`               df <- structure(list(x = structure(c(12458, 12633, 12692, 12830, 13369,      13455, 13458, 13515), class = "Date"), y = c(6080, 6949, 7076,      7818, 0, 0, 10765, 11153)), .Names = c("x", "y"), row.names = c("1",      "2", "3", "4", "5", "6", "8", "9"), class = "data.frame")          > df                x     y     1 2004-02-10  6080     2 2004-08-03  6949     3 2004-10-01  7076     4 2005-02-16  7818     5 2006-08-09     0     6 2006-11-03     0     8 2006-11-06 10765     9 2007-01-02 11153      I can plot the points and a Tukey's linear fitting (`line` function in `R`) via               plot(data=df,  y ~ x)     lines(df$x, line(df$x, df$y)$fitted.values)      which produces: ![enter image description here](http://i.stack.imgur.com/n5TGd.png) All fine. The above plot shows energy consumption values, expected only to increase, so I'm happy with the fit not passing through those two points (which will be subsequently flagged as outliers). However, "just" removing the last point and replot again               df <- df[-nrow(df),]     plot(data=df,  y ~ x)     lines(df$x, line(df$x, df$     )$fitted.values)      The result is completely different. ![enter image description here](http://i.stack.imgur.com/4Ddau.png) My need is to have ideally the same result in both scenarios above. R doesn't seem to have ready to use function for monotonic regression, besides `isoreg` which however is piecewise constant. EDIT: As @Glen_b pointed out the outliers-to-sample size ratio is too big (~28%) for the regression technique used above. However, I believe there might be something else to consider. If I add the points at the beginning of the table:               df <- rbind(data.frame(x=c(as.Date("2003-10-01"), as.Date("2003-12-01")), y=c(5253,5853)), df)      and recalculate again like above `plot(data=df, y ~ x); lines(df$x, line(df$x,df$y)$fitted.values)` I get the same result, with a ration of ~22% ![enter image description here](http://i.stack.imgur.com/hwIZp.png)


title: Are the results of OLS linear regression correct?

text: How could we know if the coefficients and R-Squared in the following table are reliable and are correctly computed or not? Is there any way to validate the correctness of the results, considering that we do not have access to the data set? ![enter image description here](http://i.stack.imgur.com/6BfCN.jpg) The descriptive analysis of the original variables is as follows, but in regression they transformed variables! ![enter image description here](http://i.stack.imgur.com/O6X1v.jpg)


title: Why is there a dramatic difference between aov and lmer?

text: I have a mixed model and the data looks like this:               > head(pce.ddply)       subject Condition errorType     errors     1    j202         G         O 0.00000000     2    j202         G         P 0.00000000     3    j203         G         O 0.08333333     4    j203         G         P 0.00000000     5    j205         G         O 0.16666667     6    j205         G         P 0.00000000      Each subject provides two datapoints for errorType (O or P) and each subject is in either Condition G (N=30) or N (N=33). errorType is a repeated variable and Condition is a between variable. I'm interested in both main effects and the interactions. So, first an anova:               > summary(aov(errors ~ Condition * errorType + Error(subject/(errorType)),                      data = pce.ddply))          Error: subject               Df  Sum Sq  Mean Sq F value Pr(>F)     Condition  1 0.00507 0.005065   2.465  0.122     Residuals 61 0.12534 0.002055                         Error: subject:errorType                         Df  Sum Sq Mean Sq F value   Pr(>F)         errorType            1 0.03199 0.03199   10.52 0.001919 **      Condition:errorType  1 0.04010 0.04010   13.19 0.000579 ***     Residuals           61 0.18552 0.00304                           Condition is not significant, but errorType is, as well as the interaction. However, when I use lmer, I get a totally different set of results:               > lmer(errors ~ Condition * errorType + (1 | subject),                         data = pce.ddply)     Linear mixed model fit by REML      Formula: errors ~ Condition * errorType + (1 | subject)         Data: pce.ddply          AIC    BIC logLik deviance REMLdev      -356.6 -339.6  184.3     -399  -368.6     Random effects:      Groups   Name        Variance Std.Dev.      subject  (Intercept) 0.000000 0.000000      Residual             0.002548 0.050477     Number of obs: 126, groups: subject, 63          Fixed effects:                            Estimate Std. Error t value     (Intercept)            0.028030   0.009216   3.042     ConditionN             0.048416   0.012734   3.802     errorTypeP             0.005556   0.013033   0.426     ConditionN:errorTypeP -0.071442   0.018008  -3.967          Correlation of Fixed Effects:                 (Intr) CndtnN errrTP     ConditionN  -0.724                   errorTypeP  -0.707  0.512            CndtnN:rrTP  0.512 -0.707 -0.724      So for lmer, Condition and the interaction are significant, but errorType is not. Also, the lmer result is exactly the same as a glm result, leading me to believe something is wrong. Can someone please help me understand why they are so different? I suspect I am using lmer incorrectly (though I've tried many other versions like (errorType | subject) with similar results. (I have seen researchers use both approaches in the literature with similar data.)


title: Box-Muller algorithm and Monte Carlo Integration

text: ![Question](http://i.stack.imgur.com/Z2g98.png) Can someone help me out here? Using Box-Muller would generate N(0,1) however X ~ N(-1, 4). How do I transform the variables to the distribution of X? And would the pseudocode include the usage of the rejection algorithm too? Appreciate your help!


title: Interpreting the Pearson Chi-Square test of a constant Poisson event rate over independent counting experiments

text: ## The Experiment Consider a set of $N$ independent time segments.   1. Measure the amount of time in given segment: $t_i$   2. Count the number of events that occur in a given segment: $m_i$  An estimate of the overall trigger rate is $R=M/T$ where $T$ and $M$ are the total time and total number of events in this analysis respectively. The estimated number of events that will in the $i^{th}$ segment, given computed rate above, is $e_i=t_iR$. One can then use Pearson $\chi^2$ test $$ X^2 = \sum_{i=1}^N \frac{(m_i-e_i)^2}{e_i} $$ to compute a test statistic for the data and from that calculate a p-value for $\chi^2 \geq X^2$ given $N-1$ degrees of freedom. * * * ## A subset of the data               t = numpy.array([  41964.,   65105.,   16178.,   16018.,   27868.,    3293.,           7153.,   10808.,    7228.,   37028.,   83448.,   86558.,          25638.,   15268.,   14753.,   25773.,   69998.,    2298.,           1858.,   28058.,    3788.,   20653.,    5483.,   22238.,           2933.,    4168.,   31823.,    6358.,   26908.,    7388.,           4558.,    8068.,   32718.,   34648.,   13588.,   22298.,          12588.,    2588.,   19972.,    7826.,   26798.,    6998.,          39253.,   62278.,   71018.,   51338.,   17838.,  106228.,          19228.,   50038.])     m = numpy.array([ 185.,  322.,   86.,   82.,  132.,   16.,   34.,   66.,   44.,         176.,  387.,  439.,  138.,   73.,   62.,  115.,  331.,    9.,           9.,  146.,   16.,  101.,   25.,  111.,   14.,   21.,  160.,          29.,  111.,   39.,   21.,   34.,  148.,  168.,   64.,  107.,          58.,   19.,  100.,   29.,  124.,   30.,  173.,  298.,  328.,         233.,   84.,  493.,   95.,  229.])      * * * ## Code & the results               import numpy     from scipy import stats     R = numpy.sum(m)/numpy.sum(t)     e = t*R     X2, p_rs = stats.chisquare(m,e)        * $R = 0.004755$   * $X^2 = 30.71$   * $p_{rs} = 0.9811$ It is unclear to me how to interpret the right-sided (one tail) p-value being so large. If I include all 115 segments in the calculation rather than just the first 50,   * $R = 0.004757$    * $X^2 = 66.59$    * $p_{rs} = 0.9999$ which only worsens my confusion. I cannot see how it is possible that I have such good agreement between the estimated number of events and the measured values. What simple (or not so simple) idea am I missing that makes sense of this?


title: Correct expression of confidence level

text: A very basic question: when noting the statistical significance of fitting some model, does one report that results are significant at the 95 percent confidence level, or at 0.05 significance?


title: Compare two sets of data

text: I have three sets of data, measured by three different devices: `A,B` and `C` of air balloon whose fall is influenced by wind. Each data sheet looks like:               A:          Longitude    Latitude    Altitude    Weight    Rotation(about main axis)   Time         ...       ...           ...       ...       ...                         ...          B:          Longitude    Latitude    Altitude    Weight    Rotation(about main axis)   Time         ...       ...           ...       ...       ...                         ...      and similarly for `C`. Are there standard techniques to compare the error measurement of `B` and `C` with respect to `A` (considered standard)?. Note that except time, all other parameter can both increase and decrease, and longitude/latitude and rotation are the only parameters that can be positive and negative. The problem with regression is that I do not have independent/dependent variables. The method of error analysis I generally use is to calculate ${\rm Err}_{x}$, ${\rm Err}_{y}$, and ${\rm Err}_{z}$, and combine them suitably when I have a function explaining something. I do not have 'the function'. (I mean I do not have a function to use `'propagation of error'`.) Precisely, I am trying to measure accuracy of devices `B` and `C` Side note: I could not find an appropriate tag.


title: Probability of getting a specific straight in poker (texas hold'em)

text: For those unfamiliar with texas hold'em, a straight is simply 5 cards in a numbered sequence - for example, `{2,3,4,5,6}`. We get those cards by choosing the best 5 of 7 cards, two being in your starting hand, and 5 in the table, drawn in succession. My goal is to calculate the probability of getting a straight, given a specific starting hand (pre-flop). For example, let's say we get the starting hand `{2,3}`. How can we calculate the probability of getting `{2,3,4,5,6}` using the community cards? We know we need the cards `{4,5,6}` to come up at least once on the community cards. In order to (hopefully) simplify the calculations, I've decided to consider each card draw independent - let's pretend to be playing with a shuffled stack of a infinite number of standard 52-card decks. We can get the exact same card twice, including those we already have in our hand. On such a game, having 5 random draws, and considering only the card number, we have a total of `13^5` possible results. My question is two-fold:   * How to efficiently compute the number of times specific cards (`{4,5,6}`) come up in the `13^5` possible outcomes. I'm wondering if - since the draws are independent - the probability of getting `{4,5,6}` is simply `(1-(12/13)^5)^3`? I'm struggling to come up with something more plausible.    * Is the "infinite number of decks" a reasonable approximation, and is there a way to (still efficiently) compute the probability without it? \--EDIT-- Can this, by chance, be calculated using a `binomial(5,1/13)` distribution, where the value of the "reverse CDF" is 1 (that is, where the sum of the PDF from 1 to infinity == 1)? A simulation seems to come up with the number `0.021`


title: How do I prove the mutual independence of a set of random variables given another set of mutually independent random variables?

text: Suppose $X_{1},\,\ldots,\, X_{n}$ and $T_{1},\,\ldots,\, T_{r}$ and two sets of random variables, with each $X_{i}$ being categorical, and each $T_{j}$ being continuous taking on all values in the range $\left[0,\,\infty\right)$. It is known a priori that the T 's are a set of observations that give rise to the distributions of the X 's. How would I prove $\overset{n}{\underset{i=1}{\bot}}X_{i}\mid\overset{r}{\underset{j=1}{\bot}}T_{j}$?


title: Multiple covariates in an ANCOVA design

text: I have found a group difference (3 groups) for my dependent variable and now I am supposed to find out if possible covariates influence this finding. What is the best way to "correct" for 5 different covariates (age, duration of illness, medication, symptom score, gender)? Shall I run 5 separate ANCOVAs (1 for each covariate) and see if the group difference remains significant or should I better run one single ANCOVA including all covariates?


title: The usage of Kolomogorov-Smirnov statistic in predictive modelling

text: Kolmogorov-Smirnov statistic... I read in wiki that it is used to compare a sample distribution with reference. However, I've interned in a company which was supposed to do predictive modeling and they keep talking about _KS for a model_. What were they talking about?


title: How the optimal separating Hyperplane can be constructed for a Support Vector Machine

text: i refer to http://www.markowetz.org/florian/FlorianMarkowetz_DiplMath_thesis.pdf on page 31-35 he tries to reformulate Vapniks proof that the primal an dual problem for the maximal margin hyperplane (optimal hyperplane) are equivalent. Unfortunately the stated version is wrong in my opinion because on page 35 the chain of equation right before (7) is not correct. At least i do not get why second last ">" should hold. substituting the supposition in the equation results in a "<" as far as i can see. Does anybody know a alternative way to prove the desired result or can tell me what i am getting wrong? I had a look in Vapniks statistical learning Theory (1998) on page 403 but this kind of notation really confused me. Thanks a lot!


title: Showing $ \frac{1}{n} \sum_{i=1}^n (X_i - \bar{X})^q \to E((X-EX)^q)$ in probability

text: I am trying to show for $X_n$ iid st. $E|X|^q < \infty$ that $$ \frac{1}{n} \sum^n (X_i - \bar{X})^q \to E((X-EX)^q) $$ in probability. I used the following approach: We note   * $ \frac{1}{n} \sum_{i=1}^n X_i \to EX$ in prob. by the WLLN, i.e., $\bar{X}_n \to EX$ in probability   * $ Y_i := (X_i - EX)^q $ are iid i.e. we use Slutsky on $\frac{1}{n} \sum^n (X_i - \bar{X})^q $ to obtain iid random variables to then use WLLN as per below: $$ \frac{1}{n} \sum_{i=1}^n (X_i - \bar{X})^q \quad \xrightarrow{\text{Slutsky}} \quad \frac{1}{n} \sum^n (X_i - EX)^q \quad\xrightarrow{\text{WLLN}} \quad E(X-EX)^q \>. $$ Now, I realised that, of course, once I applied Slutsky I have taken $ n \to \infty $ so I cant really take the last step via the WLLN. Can somebody help me fix this proof ?


title: FIML assumptions and violations

text: From what I know about full information maximum likelihood (FIML), the assumptions are missing at random (MAR), which can't be tested, and multivariate normality. The advantages of FIML are asymptotic efficiency, asymptotic normality and unbiased (or consistent?) parameter estimates. Disadvantages: standard errors are biased downward. In a regression model with binary (dummy) variables, e.g. gender, the assumption of multivariate normality is violated. What are the consequences of this violation, and what is the best way to deal with these?


title: Simple descriptive statistics in SPSS or Excel

text: Apologies if this is a stupid question. I have a data-set with 3109 people They were asked to say whether they thought Reason 1, Reason 2, Reason 3 and Reason 4 were good reasons for having a genetic test. The percentages are as follows: Reason 1: 37% yes; 28% no; 25% not sure Reason 2: 45% yes, 31% no; 24% not sure Reason 3: 45% yes; 30% no; 25% not sure Reason 4: 46% yes; 29% no; 25% not sure How can I tell if it was the same people answering 'yes', 'no' and 'not sure' for each reason? I've had a long think and tried to find the answer by reading around and can't figure it out! The data are in Excel but can be put into SPSS easily Thanks!


title: Bayes Theorem - Is this Marketing Example Correct?

text: I'm new to whole concept of Bayes Theorem and its applications to marketing. I've been trying to learn this on my own but unsure if I'm making dumb mistakes or if I'm applying the formula correctly - hopefully you can tell me! I want to get the probability that a certain age group (say 18 to 25) was the group that converted given a transaction. This is the data I used (numbers slightly changed):               group 1: age 18-24, 92 transactions, 0.65% conversion rate     group 2: 25-34, 458, 0.87%     group 3: 35-44, 480, 1.10%     group 4: 45-54, 499, 1.36%     group 5: 55-64, 582, 1.38%     group 6: 65+. 382, 1.43%      Is the following correct?   * I formulated the question as: Given a transaction, what is the probability that it came from group X?   * I took the transaction proportion and multiplied it by the conversion rate and then divided the whole thing by the sum of all the conditionals   * With this data (rounded to the nearest percent), I got 2%, 13%, 18%, 23%, 27% and 18%. So for group 1, the probability that a transaction came from 18-24 is only 2%, is this correct usage? This feels like it's right but feeling that it's right and it actually be used correctly are two different things!


title: Covariates in Repeated Measures ANOVA

text: I have made a repeated measures ANOVA which I have two independent variables, three covariates which I will calculate the effects individually for my thesis. I basically conducted two tasks with the same sample which is the within subjects part. The factor is statistifically significant. There is a mean difference for the task one and task two. The indepent variables can predict it. However, when I put a covariate,the income level, the factor becomes not significant. But the interaction between factor and the covariate is also not significant. The fact is my sample have usually higher scores for that covariate. I mean 80% of them high income, only 20% have low income. Can it be caused by this? If not, what does it mean? Your response will be greatly preciated.


title: Squared random variable $X^2$ vs $X\times X$

text: As I understand a random variable represents all possible outcomes of an experiment with their associated probabilities. Why $X^2$ is understood as squaring outcomes of experiments instead of as multiplying the results of two identical experiments? Let's X be a throw of a die. I've seen examples where expected value of $X^2$ is calculated as: $$ 1^2\dfrac{1}{6} + 2^2\dfrac{1}{6} + 3^2\dfrac{1}{6} + 4^2\dfrac{1}{6} + 5^2\dfrac{1}{6} + 6^2\dfrac{1}{6} = 15.16 $$ So here we just squaring the outcome of a single experiment. But in math square can be replaced by multiplication. Like: $X\times X$. And for multiplication of two random variable we probably need to construct a product distribution. As I understand it should be something like throwing 2 dice, multiplying the results, and then calculating the average of these values. So as both X are independent in this case expected value can be found as $E[XY]=E[X]E[Y]$. I.e. $3.5 *3.5 = 12.25$


title: Data Mining technique for preferred vendor problem?

text: Currently I use open auction method for selecting a transporter and want to switch to preferred transporter where I would select a transporter based on my preferred transporter list.   The project includes identification of preferred transporter based on historical analysis. The only data/ parameters I've is : 1) Cost: How many times I've allocated work to this transporter? I've always allocated work to transporter based on lowest price through open auction method. 2) Quality & Reliability: Difference in Detention Cost charged by transporter and actual detention cost calculation. 3) Lead Time: Which transporter helps me transporting material in lowest time. Using the 3 parameters I decided to develop a C5.0 Tree in SPSS Modeler but results are not satisfactory. What according to all would be a good data mining technique of solving this problem and any literature you would suggest for modelling this problem ?


title: Confidence intervals for predictors in multivariate logistic regression

text: I've got a question. I am dealing with medical data which contain 5 predictors and 1 binary outcome. When I try to classify the data using all 5 predictors I get 0.84 area under roc-curve which is quite good for me. However when I consider confidence intervals for these 5 parameters, I see that just 1 of them is significantly different from 0, its confidence interval doesn't include value "1". For all other predictors confidence interval include value 1 (from 0.78 to 1.2, for example). When I try to classify the dataset based on just this 1 predictor (the one with good confidence interval), I still have good confidence interval for it but the area under roc-curve is 0.74 and combinations of specificity and sensitivity are really worse than in first example. The question is: how should I treat the situation when not all of the predictors are significantly different from 0? Can I consider such deciding rule as a good one and use it in further analysis? Thanks a lot!


title: Good text on Clinical Trials?

text: I'm an undergraduate statistics student looking for a good treatment of clinical trials analysis. The text should cover the fundamentals of experimental design, blocking, power analysis, latin squares design, and cluster randomization designs, among other topics. I have an undergraduate knowledge of mathematical statistics and real analysis, but if there's a fantastic text that requires a bit higher level of statistics or analysis, I can work up to it.


title: What is the difference between PCA and asymptotic PCA?

text: In two papers in 1986 and 1988, Connor and Korajczyk proposed an approach to modeling asset returns. Since these time series have usually more assets than time period observations, they proposed to perform a PCA on cross-sectional covariances of asset returns. They call this method Asymptotic Principal Component Analysis (APCA, which is rather confusing, since the audience thinks immediately of asymptotic properties of PCA). I have worked out the equations, and the two approaches seem numerically equivalent. The asymptotics of course differ, since convergence is proved for $N \rightarrow \infty$ rather than $T \rightarrow \infty$. My question is: has anyone used APCA and compared to PCA? Are there concrete differences? If so, which ones?


title: Modelling longitudinal data where the effect of time varies in functional form between individuals

text: **Context** : Imagine you had a longitudinal study which measured a dependent variable (DV) once a week for 20 weeks on 200 participants. Although I'm interested in general, typical DVs that I'm thinking of include job performance following hire or various well-being measures following a clinical psychology intervention. I know that multilevel modelling can be used to model the relationship between time and the DV. You can also allow coefficients (e.g., intercepts, slopes, etc.) to vary between individuals and estimate the particular values for participants. But what if when visually inspecting the data you find that the relationship between time and the DV is any one of the following:   * different in functional form (perhaps some are linear and others are exponential or some have a discontinuity)   * different in error variance (some individuals are more volatile from one time point to the next) **Questions** :   * What would be a good way to approach modelling data like this?    * Specifically, what approaches are good at identifying different types of relationships, and categorising individuals with regards to their type?   * What implementations exist in R for such analyses?   * Are there any references on how to do this: textbook or actual application?


title: Forecast R package producing flat predictions

text: I've just started playing with the R forecast package and found I must be doing something wrong because I can't get a decent prediction for a simple sinus.               weightData <- data.frame(weight = sin(seq(1:100)), week=1:100)     weight <- as.numeric(weightData$weight)     predicted <- forecast(weight,h=3,level=95)           # see the predicted values by forecast predicted          myplot <- forecast(weight,h=10,level=95)     plot(myplot)      And I get a flat prediction. I understand the generic forecast methods selects the best method for my data. Isn't that true? Am I missing something? Thanks in advance!


title: Using the sde package in R to simulate a SV model with leverage

text: Using the sde package in R, I would like to simulate the following model for stock prices $p_t$: $\mathrm{d}\sigma^2_t = (\theta_1 - \theta_2\sigma^2_t)\mathrm{d}t + \theta_3\sigma_t\mathrm{d}W_{\sigma,t}$ (CIR model used for stochastic volatility) $\mathrm{d}\log{p_t} = \sigma_t\mathrm{d}W_{p,t}$ where correlation between $W_{p,t}$ and $W_{\sigma,t}$ is possibly non-zero. Q: Is that possible? Is it even possible to simulate multivariate SDEs in the "sde" package the way one can in S+FinMetrics using the `gensim` functions? Is there another package in R that might do this? I'm able to simulate the variance process simply enough with the following code:               library(sde)     sig2 <- sde.sim(X0=0.04,  theta=c(0.3141, 8.0369, 0.43),  model="CIR")      I can then simulate the price series (starting with initial price = 100) using:               logr <- rnorm(n=length(sig2),sd=sqrt(sig2))     logp <- cumsum(c(log(100),logr))     p <- exp(logp)      But this approach seems unnecessarily clunky and can't capture non-zero correlation between the Brownian Motions.


title: What makes the formula for fitting logistic regression models in Hastie et al "maximum likelihood"?

text: I am learning logistic regression from The elements of statistical learning: data mining, inference, and prediction, by Trevor Hastie, Robert Tibshirani, Jerome H. Friedman. Suppose $G$ is a random variable representing response taking class label in $\\{1,\ldots,K \\}$, and $X$ is a random vector representing predictor taking values in $\mathbb{R}^n$. After modeling $P(G = k|X = x; \theta), k=1,\ldots,K, x \in \mathbb{R}^n$ by logistic model in section 4.4 Logistic Regression, the parameter $\theta$ is estimated from samples $\\{(x_i,k_i), i=1,\ldots,N \\}$ in section 4.4.1 Fitting Logistic Regression Models: > Logistic regression models are usually ﬁt by maximum likelihood, using the > conditional likelihood of $G$ given $X$. Since $P(G|X)$ completely specifies > the conditional distribution, the multinomial distribution is appropriate. > The log likelihood for $N$ observations is $$ℓ(θ) = \sum_{i=1}^N \log P(G = > k_i|X = x_i; \theta) $$ **My question** is how the estimation method is maximum likelihood. Generally speaking, MLE is to estimate parameter $\theta$ of a **single** distribution $p(z;\theta)$ from its **i.i.d. samples** $\\{z_i, i=1,\ldots,N \\}$ by maximizing $\log p(x_1,\ldots, x_N; \theta) = \sum_{i=1}^N \log p(x_i; \theta) $, the log of joint distribution of the samples. The estimation for logistic regression as in the quote seems not using MLE, because:   1. $P(G = k|X = x_i; \theta), i=1,\ldots,N $ are $N$ different distributions, instead of just one distribution as for MLE;   2. in $ℓ(θ) = \sum_{i=1}^N \log P(G = k|X = x_i; \theta) = \log \prod_{i=1}^N P(G = k|X = x_i; \theta)$, I wonder how to explain that "${ (k_i|x_i), i=1,\ldots,N }$" are independent and therefore $\prod_{i=1}^N P(G = k_i|X = x_i; \theta)$ is their joint distribution, as I don't know how to understand $(k_i|x_i)$ which seems not make sense. **My guess answer** is that $\\{(x_i,k_i), i=1,\ldots,N \\}$ are assumed to be i.i.d. samples of the single distribution for $(X,G)$. So MLE for the distribution of $(X,G)$ is to maximize $$ \log P((x_1,k_1),\ldots,(x_N,k_N)) = \log \prod_{i=1}^N P(x_i,k_i) $$ $$ = \log \prod_{i=1}^N P(k_i | x_i) P(x_i) = \log \prod_{i=1}^N P(k_i | x_i) + \log \prod_{i=1}^N P(x_i) $$ $P(k_i | x_i)$ is parameterized by $\theta$ as $P(k_i | x_i;\theta$, while $P(x_i)$ is assumed to be independent of $\theta$, so to solve for $\theta$ by MLE for the distribution of $(X,G)$ is equivalent to maximize $\log \prod_{i=1}^N P(k_i | x_i; \theta)$. I wonder if this is a reasonable explanation, and what yours are? Thanks and regards!


title: How to interpret the margin of error in a poll?

text: Recently the media reported on a political poll that stated that "46% of Republican voters in Mississippi think that interracial marriage should be illegal". One example story (of many around the 'net) is this one from the NY Daily News: Interracial marriage should be illegal, say 46% of Mississippi Republicans in new poll. From that article: > A survey conducted last month shows 46% of GOP voters in the state believe > interracial marriage should be illegal - a plurality of the people > questioned Given that I know many people in interracial marriages (including in my own family) I was motivated to track down the full story rather than fall to the obvious hysteria. I located the press release on the poll from Public Policy Polling here (PDF) which has: > 46% of these hardcore Republican voters believe interracial marriage should > be illegal, while 40% think it should be legal. With Barbour included, > Huckabee gets more support (22%) from the former than the latter (15%), as > does Palin (13-6). The support for Bachmann (10-2), Gingrich (13-8), and > Pawlenty (4-1) works the opposite way. And at the end of the first page they state: > PPP surveyed 400 usual Mississippi Republican primary voters from March 24th > to 27th. The survey’s margin of error is +/-4.9%. Now what I don't understand is how to interpret the results especially the margin of error value. From the raw numbers, 46% (~50%) of 400 voters is 200 voters who responded who think that interracial marriages should be illegal. But from stats of Mississippi in the 2008 Federal election, McCain polled about 725,000 Republican votes (to Obamas 556,000). My gut feeling is that you can't extrapolate that 46% number to the other 724,600 Republican voters and still retain a 5% margin of error. So is the media misrepresenting the numbers (really? the media would do that?!?!?!) Or are there statistics at work that I don't understand? Thanks for your help!


title: Bound on variance of sum of variables

text: Suppose I have two finite sets of data `A` and `B`, with equal length `n`. What's the best upper and lower bound I can make on `var(A+B)`, in terms of `var(A)` and `var(B)`?


title: Singular value decomposition of a three-dimensional array

text: The _**Singular Value Decomposition (SVD)**_ of a matrix is $$A_{m\times n} = U_{m\times m}\Lambda_{m\times n} V_{n\times n}'$$ where $U$ and $V$ are orthogonal matrices and $\Lambda$ has ( _i_ , _i_ ) entry $\lambda_i \geq 0$ for $i = 1, 2, \cdots , min(m, n)$ and the other entries are zero. Then the _left singular vectors_ $U$ for rows of matrix and _right singular vectors_ $V$ for columns of matrix can be plotted on the same graph called _**bi- plot**_. I'm wondering how to do the _**SVD**_ of a three dimensional array and plot the singular vectors on the same graph like _**bi-plot**_. Thanks


title: Does GSVD implement all linear multivariate techniques?

text: I came across the article by Hervé Abdi about generalized SVD. The author mentioned: > The generalized SVD (GSVD) decomposes a rectangular matrix and takes into > account constraints imposed on the rows and the columns of the matrix. The > GSVD gives a weighted generalized least square estimate of a given matrix by > a lower rank matrix and therefore, with an adequate choice of the > constraints, the GSVD implements all linear multivariate techniques (e.g., > canonical correlation, linear discriminant analysis, correspondence > analysis, PLS-regression). I'm wondering how the GSVD is related to all linear multivariate techniques (e.g., canonical correlation, linear discriminant analysis, correspondence analysis, PLS-regression).


title: Gradient decay in neural networks

text: I read that in traditional feed-forward neural nets the gradients in the early layers decay very quickly and that this is 'a bad thing'. But I don't understand why. Can someone please explain what does that mean? or at least direct me to places to read about it and understand.


title: How can I even out a random distribution while minimising how far each data point is moved?

text: I have a software application which uses a queue and multiple processors to process those jobs. Jobs get re-run on a daily basis for customers, but we also have new customers signing up regularly. The problem is that customers generally sign up during office hours, and the daily re-running of the jobs simply schedules the job to run at the same time every day. This means that we have built up quite large spikes of work during the day, with very few jobs running overnight. I would like to even out the spikes in this job queue, but I also want to minimise how far each individual job is moved from it's original time slot. In other words, I would rather move 100 jobs by 5 minutes each than move 1 job by an hour to achieve an even distribution of jobs over the whole 24 hour period.


title: How to perform text mining, sentiment mining, and business category identification, and where to obtain a categorization library

text: Currently I am working with Text Mining which includes sentiment identification and assigning corresponding business categories using open source tool R. I found these two documents which helped me to some extent:   * http://www.jstatsoft.org/v25/i05/   * http://epub.wu.ac.at/1923/1/document.pdf My approach is to tokenize the text and then lookup for sentiment and business category. To do this I require positive and negative libraries for sentiment mining and category file which contains word and category. I was able to get positive and negative words but was unable to get a categorization library.   * Where can I get a categorization library?    * Is the above approach appropriate? Is there a better way to do this?


title: Comparing two populations over time (longitudinal population study)

text: I have two sets of data for men and women, measuring the turnover of the most popular points of one characteristic of the population, for each year since 1912. This is popular data from an insurance company, and there is no sample size information available, just the compiled rankings data. Intuitively, it looks like turnover is getting faster for women, but I would like to know how to test this for statistical significance. I have not used longitudinal data before in this way, so I would appreciate any help! Thanks very much in advance.


title: Partial $\eta^2$ in multivariate analysis

text: I wonder why the partial $\eta^2$ values are so small in my SPSS GLM Multivariate analysis. I have 3 DVs and 5 IVs and the sample size is large ($n=5000$). The main effects and significant but the effect sizes are above $0.01$ (partial $\eta^2$). When I do the same analysis with Univariate (one DV at a time) the significance remains but the effect sizes increase up to $0.60$. Why? The factors correlate with each other ($0.03-0.06$) so, I guess I need to do Multivariate although the effect sizes are so small. Can you explain why the effect sizes change so much doing Multivariate versus Univariate?


title: Regression and point-specific p-values (using R for explanation)

text: Consider the following (in R):               library(MASS)     plot(stack.loss~Air.Flow,data=stackloss)     regression <- rlm(stack.loss~Air.Flow,data=stackloss)     abline(regression)      For each of the points I would like to to test against the null hypothesis that it really is located on the line (rather than where it really is) - in essence quantifying the degree to which these are outliers with respect to the regression model. Possibly using hopeless search term combinations, I have been unable to identify appropriate methodology. Would it be valid to use for each residual a one-sample t-test (ar equivalent) against the residual population to establish such a measure? Edit: Further reading seems to indicate that "Tolerance Interval" is what I might be looking for - the "tolerance" R package provides calculations of that. I am however puzzled by the apparent contradictory nature of such intervals for regressions as defined by different publications. In the R package case the tolerance bands seem to track the regression line parallely (see Fig. 12 and below example), while alternative sources such as this operate with bands that are reminiscent of confidence intervals in shape (see Fig. 4). The latter seems more intuitive, but other's opinions would be appreciated. For each point in a sample I am actually looking for "which x/95% tolerance band is running through this point" ... Here's how usage of "tolerance" might look like (adding to above code) for 95%/95% 2-sided nonparametric regression tolerance bounds:               library(tolerance)     tol.bounds <- npregtol.int(x=stackloss$Air.Flow,                 y=stackloss$stack.loss, y.hat=regression$fit, side=2, alpha=0.05, P=0.95,          method="WILKS")     lines(tol.bounds$x,y=tol.bounds$"2-sided.lower",col="red")     lines(tol.bounds$x,y=tol.bounds$"2-sided.upper",col="red")


title: Underdispersed count data and factor analysis

text: I am working with count data from a psychological experiment. Data are from several subtests, i.e. scores on these tests. The way these scores are generated is clearly not Gaussian (not only because they are discrete, but also since the scores are basically sums of correctly solved tasks). As a baseline I worked with models using the Poisson distribution. For example the Rasch (1960) Poisson counts model (which is a type of GLM), for which I also estimated the dispersion from the data (1.01, bootstrapped CI for the dispersion parameter includes 1). To go one step further I would like to fit factor analytic models. I am aware of work of Wedel et al (2003, Journal of Multivariate Analysis), who deal with the situation with and without overdispersion. Also MPlus is able to handle the overdispersed situation. However my data shows only little dispersion (in fact the raw scores on some subtests are clearly underdispersed). I went ahead anyway and tried to use MPlus for (ML) factor analysis with Poisson errors and log-link. I used MPlus mainly for convenience, since it is the only package that handles count data and factor analysis (feel free to correct me here). Using MPlus did not work well with my real data, but did OK with simulated data, at least when the latent factors had unit variance. Factor loadings are harder to recover when the variance is reduced (say to 0.25^2) and is not working (at least with MPlus) when the variance is very low (say 0.1^2). I suspsect that factor analytic models for count data need at least some variance in the underlying (latent) factors so that they can be estimated, since otherwise the observed variability is attributed to the Poisson (neg- Bin, zero-inflated Poisson, ...) distribution rather than to the latent factor. I could not find a reference for this claim and my own simulations are rather limited. 1) Is there any reference for this phenomenon? Is it maybe trivial from another point of view? 2) Clearly I could treat the data as continuous, use standard exploratory and confirmatory factor analysis, but I feel this is neither correct nor satisfying. What are alternative strategies?


title: Percentage change vs. first difference

text: I want to regress (using an OLS regression) survey results on stock returns. How do I decide whether I use the first difference or the percentage change of the survey results?


title: Algorithm for generating hyperexponential distribution

text: I need to generate a hyperexponential distribution for my project. I have already implemented a poisson generating algorithm given by Donald Knuth, but I couldn't find an algorithm for generating a hyper exponential random variable. I am provided with the mean and variance required of the distribution and I need an algorithm which can generate a random variable from this distribution when I execute it.


title: 2x2 design for ANOVA

text: I am beginner with statistics and I need some help for my data analysis. In my research I am trying to find an impact of colors on websites and types of webstores on customer trust, purchasing intention, preference etc. For example websites selling insurance service will be trusted more if the color of websites is blue than red. Design of my research is 2x2 /color: red vs. blue $\times$ type of website: service vs. outlet. I am not sure if I should use one-way anova or two-way anova and how the model should look like. Or is there a better solution? Could you help me please? I am using SAS.


title: Is it possible to calculate mutual information by moments generating functions?

text: I went to listen to a workshop and some audience asked the presenter how the moments can improve the mutual information. I am learning the MI(Mutual Information) and moments so don't have enough knowledge to understand what it means. As I know so far, the moments generating functions find moments mathematical presentation for a random variable (Please correct me if I am wrong). Then, I have done some research but still have some confusions. I am wondering if anyone has more knowledge about this can clarify for me. Here are my questions.   * Mutual information is usually calculated by bin functions to estimate the probability of two random variable which can be a case of two vectors X and Y. Is the moment generating function another way to estimate probability ?   * If moment generating functions can present the probability of X and Y, how do we calculate it ?   * Does a MI have a moments generating function ?   * If MI has a moments generating function, how can we present a MI of X and Y by its moments functions ? Thank you very much,


title: Reject null hypothesis or not?

text: Tony conducts a hypothesis test on whether the proportion of all students in his highschool who bike to school (denoted as $p$) equals $30\%$. Specifically, Tony has $H_0:p=0.3$ versus $H_A:p≠0.3$. He obtains a P-value of $0.01$. On the other hand, John would like to test if there is sufficient evidence to support that $p$ is greater than $0.3$ at the $10\%$ significance level. Based on Tony's result, will the null hypothesis of John's test be rejected? Here is what I have done. Tony: $H_0:p=0.3$ , $H_A:p≠0.3$ with P-value of $0.01$ John: $H_0:p=0.3$ , $H_A:p>0.3$ with $\alpha = 0.1$ From Tony's result, p-value = $0.01 $, which is less then $\alpha = 0.1$. Hence the null hypothesis of John's test should be rejected. Do I have the right procedures?


title: R-squared in linear model verses deviance in generalized linear model?

text: Here's my context for this question: From what I can tell, we cannot run an ordinary least squares regression in R when using weighted data and the `survey` package. Here, we have to use `svyglm()`, which instead runs a generalized linear model (which may be the same thing? I am fuzzy here in terms of what is different). In OLS and through the `lm()` function, it calculates an R-squared value, the interpretation of which I do understand. However, `svyglm()` does not seem to calculate this and instead gives me a Deviance, which my brief trip around the internet tells me is a goodness-of-fit measure that is interpreted differently than an R-squared. So I guess I essentially have two questions on which I was hoping to get some direction:   1. Why can we not run OLS in the `survey` package, while it seems that this is possible to do with weighted data in Stata?   2. What is the difference in interpretation between the deviance of a generalized linear model and an r-squared value? Thanks.


title: Simulating data for a multiple linear regression

text: How can I simulate a data for multiple linear regression (one IV and more then 2 IV) when we know the skewness and kurtosis of the data? Another question is, how can we determine the type of distribution (ex: gamma dist) based on the skewness and kurtosis? We know that for normal distribution, the skewness and kurtosis are (0,0)... Thanks


title: From joint cdf to joint pdf

text: We can get the joint pdf by differentiating the joint cdf, $\Pr(X\le x, Y\le y)$ with respect to x and y. However, sometimes it's easier to find $\Pr(X\ge x, Y\ge y)$. Notice that taking the complement doesn't give the joint CDF, so we can't just differentiate and flip signs. Is there still some simple rule to differentiate this and get the pdf?


title: exponential density & bernoulli distribution

text: I had asked a question on Maximum Likelihood earlier. Now I have two questions that are related to this question which had: Let $x$ have an exponential density $p(x|\theta) = \theta e^{-\theta x} \text{ if }x \ge 0;\quad 0 \text{ otherwise.}$ Now suppose that n samples $x_1, ... , x_n$ are drawn independently according to $p(x|\theta)$ Show that the maximum-likelihood estimate for $\theta$ is given by $$\hat{\theta} = \frac{1}{\frac{1}{n} \,\, \sum^n_{k = 1} x_k}$$ The second question is similar Let $x$ be a $d$-dimensional binary (0 or 1) vector with a multivariate Bernoulli distribution $$p(x|\theta) = \prod^d_{i=1} \theta_i^{x_i}(1-\theta _i)^{1-x_i}$$ where $ \theta = ( \theta_i, ..., \theta_d)^t $ is an unknown parameter vector, $\theta_i$ being the probability that $x_i = 1$. Show that the maximum-likelihood estimate for $\theta$ is $$ \hat{\theta} = \frac{1}{n} \sum^n_{k = 1} x_k $$ how would I show that the maximum-likelihood estimate is given for either of them. I've been reading tutorials and stuff and I still don't quite understand it.


title: Which value is the p-value?

text: It's a right-tailed Population Proportion test. The critical value is 1.645 and the calculated test statistic is -1.14. The area from minus infinity to -1.14 is 0.1271. **The question is, since it's a right tailed test is the value 0.8729 (1-0.1271) or is it just 0.1271?** Any help would be much appreciated!


title: Mouse models - 'replicates' and analysis

text: I'm hoping to hear from someone who has worked on mouse models or similar biological analyses where there is a tendency to run 'replicates' of an experiment. I know multiple testing is a sizeable kettle of fish which is definitely relevant to this discussion. I have some applications for projects where they talk about running 3 replicates of an experiment, where each experiment has n = 3 to 7. However there seems to be no mention of what they will do with the multiple sets of results, as in how they will handle success, failure, success vs failure, failure, success etc. It seems like this 'replicates' approach is quite common practise in this field. What are your thoughts / experiences with this situation? I know there are different types of replication, technical vs biological, however I've found little useful reading on this issue.


title: Cause-specific mortality (or survival) analysis

text: I am looking for the statistical test that shows the significant contribution of one factor to the overall rate (cause and death rate). Our study is about the number of people who died from year 2008 to 2012 and the cause of death. What test am I going to use to know the individual contribution of the causes to the total number of death?


title: Interplay of Categorical Data

text: I apologize if my phrasing is un-educated, I don't have much of a stats background. I'm struggling to find info answering my question because I lack the terminology to accurately describe it. Let's say that I was reviewing data on mp3 player sales, and wanted to see if certain combinations of categorical data resulted in returns. Things like color, memory size, brand. I don't simply want to see the relationship between returns and the color red, but rather which combinations indicate a likely return. For example, imagine that a certain manufacturer uses a different factory to produce each color mp3 player - and the factory that produces blue ones has a defect that shows up 3 months later and is causing returns. I'd be able to see a relationship between that color and brand with returns in that example. Does that make sense? What test could I leverage for this? My jumping off point is a chi-square, except that I can't see the interplay between the different variables with that. I haven't been able to figure out where to go from there.


title: When is there no point in using regression?

text: Under what circumstances would using regression with two given variables not increase accuracy of prediction?


title: How to compute margin of error with a given confidence interval?

text: I was given the following question: A survey found that 89% of a random sample of 1024 American adults approved of cloning endangered animals. Find the margin of error for this survey if we want 90% confidence in our estimate of the percent of American adults who approve of cloning endangered animals. I know that for 90% Confidence, $\text{ME}\sim 0.82/\sqrt{n}$. I attempted using this formula with n equal to both 1024 and (.89)1024. I got 0.025625 and 0.02716, respectively. The answer given for the problem is 1.61%. I do not understand where I went wrong. Perhaps I am using the Margin of Error formula incorrectly? Thanks. :)


title: Sample size determination for Bayesian logistic regression

text: Suppose we have three different treatments A, B, and C and we are modeling response(0,1). I want to know how to do sample size determination for each group(A,B,C) given say normal priors on each group. I need to show how we can (and how much we can reduce) the sample size and get similar testing of differences in main effects if we have prior info. the most important thing is being able to show how we can reduce the sample size if we incorporate the prior info.


title: How to assess whether weighted samples come from the same distribution?

text: Suppose I have two samples, and I'd like to compare them to see if they come from the same underlying distribution. Normally, I would probably do this using a KS test. However, suppose also that I have some 'weight' for each observation in each of the two samples indicating how 'sure' I am that the sampled value is 'accurate'. In other words, some observations in the samples are more representative of the sampled distribution than others (note that these weights will differ between samples). With this in mind, my question is: Are there any tests (like a KS test) that I can use to test if these two samples come from the same underlying distribution, taking into account the weights associated with each sample? Hopefully I'm not being too vague here... but any help and/or comments are greatly appreciated.


title: Observation symbols for training a set of HMMs

text: If we are to classify 2 separate classes/actions using HMMs, we design 2 separate HMMs (one for each class). Do they share **same set or a different set** of observations-symbols for each of the HMM? If they share different set of symbols, during the time of extracting a symbol (using vector quantization) how do we know which set of symbols to generate in order to classify them?


title: Chi-square independence test for A/B split testing

text: I have the results of an A/B split test of a web page.               Page  Impressions   Clicks    CTR            A     56,000        10        0.018%         B     78,000        21        0.027%          I'd like to test the ctr for statistically significant difference. Null hypothesis - there is no difference in CTR between pages A & B. How would I go about doing this on a more theoretical level? As in what should be rows, column, setting expected values?


title: Count data distribution

text: I am currently dealing with count data with $n=360$. The data fits neither the Poisson nor the negative binomial distribution. Additionally, I don't think my data has excess zeroes, so there is no point of considering zero-inflated model. Now I am considering an alternative model (COM Poisson). What else do you suggest?


title: Using R and plm to estimate fixed-effects models that include interactions with time

text: I am using `plm()` to estimate fixed-effects models of the form               y ~ x + time + time:fixed_trait      where `fixed_trait` is a variable that varies across individuals but is constant within individuals. The point of interacting `time` with `fixed_trait` is to permit the effect of `fixed_trait` to vary across time. (I am working here from Paul Allison's recent booklet on fixed effects. Citation appended.) `plm()` has no trouble estimating coefficients and standard errors for such models. But `summary.plm()` can't calculate R^2 for these models. This is the problem that I would like to fix. Here is a minimal example:               library(plm)     tmp <- data.frame(ID=rep(1:3, 2), year=rep(0:1, each=3),                        y=rnorm(6), const=rep(1:3, 2))     fe1 <- plm(y ~ year,              index=c('ID', 'year'), data=tmp)     fe2 <- plm(y ~ year + year:const, index=c('ID', 'year'), data=tmp)     summary(fe1)  # works fine     summary(fe2)  # Error in crossprod(t(X), beta) : non-conformable arguments      Delving into `plm:::summary.plm` makes the problem clearer. To calculate R^2, `plm` tries to do this:               beta <- coef(fe2)     X <- model.matrix(fe2)     crossprod(t(X), beta)      This doesn't work because `beta` only includes estimates for `year1` and `year0:const`, while `X` also includes a column for `year1:const`. In other words, `X` includes columns for both `year0:const` and `year1:const`, and it's impossible to estimate both of those coefficients. A workaround is to create the interaction term "by hand" before entering it into the formula:               tmp$yearXconst <- tmp$year*tmp$const     fe3 <- plm(y ~ year + yearXconst, index=c('ID', 'year'), data=tmp)     summary(fe3)  # works fine      But this is cumbersome. Short of this, is there anything that I can do to make `summary.plm` work with such models? === Allison, Paul D. 2009. Fixed Effects Regression Models. Los Angeles, CA: Sage. See especially pages 19-21.


title: Is correlation equivalent to association?

text: My statistics professor claims that the word "correlation" applies strictly to linear relationships between variates, whereas the word "association" applies broadly to any type of relationship. In other words, he claims the term "non- linear correlation" is an oxymoron. From what I can make of this section in the Wikipedia article on " _Correlation and dependence_ ", the Pearson correlation coefficient describes the degree of "linearity" in the relationship between two variates. This suggests that the term "correlation" does in fact apply exclusively to linear relationships. On the other hand, a quick Google search for "non-linear correlation" turns up a number of published papers that use the term. Is my professor correct, or is "correlation" simply a synonym of "association"?


title: Does it make sense to "cluster" when you use a regression discontinuity?

text: One of the breakthroughs of econometrics over the past two decades has been to employ "clustering" to take into account the correlation of error terms across observations. For instance, if you're evaluating the effect of an educational intervention where you have data on individual students but you suspect that teachers implemented the intervention differently, it is common to analyze the data in a way that recognizes that there are common effects at the "class" level. A common correction is to use clustering. When you run a regression discontinuity, do you similarly have to take into account that your observations may be clustered? If so, how is the estimator implemented differently?


title: Similarity metric not based on statistical correlation

text: I'm building a simple recommender algorithm and as a prerequisite to calculating a prediction, I need to assign a similarity metric between each two pair of users. The next step being calculating a weighted average between all users. My problem with Pearson's or cosine is this:   A, B are users whose vectors represent some ranking they have made on 6 items in a scale of 5.   Let **A = [1,1,1,0,1,1]**   Let **B = [5,5,5,0,5,5]**   Then the similarity s is:               s(pearson)=1       s(cosine)=1      So basically there's a full correlation between A,B, although user A ranks low and user B ranks high, and the users aren't literally “similar”. I would like to use a metric that relates how 'close' A,B are to each other.   I used the following algorithm:   **1.** for each item `i` both users A,B have assigned a value `r` to, I used   $${{q}_{i}}=\frac{\left| {{r}_{A,i}}-{{r}_{B,i}} \right|}{scale}$$ **2.**   $$c(A,B)=1-\frac{\sum\limits_{i}{{{q}_{i}}}}{\sum\limits_{i}{{}}}$$   **3.** Then from a scale of 0 to 1 (1 being 'alike') `c(A,B)` is the similarity metric (c for closeness). **EDIT:**   **I just noticed that what I did is just a variation of a normalized Minkowski Distance with a dimension of`r=1`** I know it's a bit general, but is the metric `c` makes any sense in respect to the function I want it to serve? More to the point, will I see any improvement in the prediction if instead of doing a weighted average on let's say 10 users with the highest `s` to the user I try to accommodate, I'll use the ones with the highest `c`?


title: Estimation of the signal to noise ratio of an image

text: I originally had a multispectral satellite image consisting of 5 bands. I applied principal components to these bands. So an example of the process would go something like this: reality -> satellite image -> first principal component Is there any way to estimate the signal to noise ratio here? Of this first component? Of the satellite image in itself? I know very little of signal processing so please bear with me if I dont have a good grasp of these concepts.


title: Error count limit for known residual failure rate and desired false-positive rate?

text: We perform a known fixed number $K$ of independent experiments, each of which has known/assumed residual odds $p$ to fail. We count the number $N$ of failures. How can I choose $M$ so that $N>M$ has odds at most $q$, for some low desired overall false-positive rate $q$? I'd prefer an approximation based on $q$ and $E = K \cdot p$ only, with a quantitative rule on when this approximation is valid; this should include $10^{-7}\le q\le 10^{-4}$, $0.5 \le E \le 10$, $K\ge 10^4$. Pointer to an authoritative source is prefered. My real-life application: I'm designing devices with a self-test (of a physical random number generator) that has a known residual failure rate $p$ when there is no defect. A device self-destructs if the self-test fails more than $M$ times in its life. I want to choose (one of) the lowest $M$ such that failure rate after performing $K$ self-tests is no more than $q$ for devices with no defect. Please retag as approriate!


title: method and needed tests for assumption in analyzing a paired sample Likert measure?

text: I am writing an exploratory paper considering an argument that a certain reaction that can arise in the viewing of art can serve as an empirical indicator that another particular cognitive process has occurred within their experience (in this case that one has changed some component of their own self-conception). (for example: say that when people report that they "felt like sneezing" in an art encounter, we would also expect that they should have changed some component of their self image. Those who did not feel like sneezing should show no change.) To measure this, we gave a pre-test with twenty 7-point bipolar adjectival scales (e.g., good...bad) assessing people's conception of self image, then had them view art, and then administered a second test with the same questions as well as a question assessing (Y/N) whether they felt like sneezing. The population size was 21 individuals taken from a museum population, with about 40% feeling like sneezing. I am wondering what would be the best test to compare changes in self conception between the sneeze and no-sneeze populations. I originally did a paired sample t-test on both populations, and showed a significant change in certain measures for the sneezing group. However have been told that this method is not appropriate. Any suggestions would be appreciated.


title: Ljung-Box test: Basic Questions

text: I'm interested into the Portmanteau tests for fitting ARIMA models and ended up to Ljung-Box test and their implementations for R-statistics. I've read already some relative question on the subject and came to the conclusion that the safest way to reject the null hypothesis is when p<0.05. However I was browsing this site, where the author gets the following results for the test and h=20: X-squared = 13.5844, df = 20, p-value = 0.851 and states that: " _and the p-value for the Ljung-Box test is 0.9, we can conclude that there is very little evidence for non-zero autocorrelations in the forecast errors at lags 1-20_ " what am I missing here? thank you in advance for your patience.


title: Count data and binomial distribution

text: My problem described in buckets and balls (this is stats after all)... I have 1 million buckets.   In each bucket I have some number of balls. I want to model this system so that I have the expected number of balls for a bucket. (I expect my distrubtion of balls per bucket to be "gaussian" shaped) Since my counts are integers I'm immediately pointed to the binomial distribution. However, I can't seem to wrap my head around the notion of successes and failures in my system since every observation is a success. Can anyone point out the vital piece of information that would help me figure out what distribution to use? **Edit** : I want to fit a distribution to my system to model the buckets and balls because these represent the noise in another problem. If I can parameterize this distribution (i.e., if it is binomial), then I can use that to pull my signal from the noise in another problem.


title: LMNN: Simple Algorithm Description

text: B"H Hello, I understand the principles of LMNN (Large-Margin Nearest Neighbor) classification very well, but not all that well Weinberger's pseudo-code for its implementation. Can anyone simply (but thoroughly) describe a simple implementation of the LMNN algorithm (i.e. calculation of the Mahalanobis matrix - **M** ), or direct me to such a document? Thanks much.


title: Significant effect in lme4 mixed model

text: I use lme4 in R to fit the mixed model               lmer(value~status+(1|experiment)))      where value is continuous, status and experiment are factors, and I get               Linear mixed model fit by REML      Formula: value ~ status + (1 | experiment)        AIC   BIC logLik deviance REMLdev      29.1 46.98 -9.548    5.911    19.1     Random effects:      Groups     Name        Variance Std.Dev.      experiment (Intercept) 0.065526 0.25598       Residual               0.053029 0.23028      Number of obs: 264, groups: experiment, 10          Fixed effects:                 Estimate Std. Error t value     (Intercept)  2.78004    0.08448   32.91     statusD      0.20493    0.03389    6.05     statusR      0.88690    0.03583   24.76          Correlation of Fixed Effects:             (Intr) statsD     statusD -0.204            statusR -0.193  0.476      How can I know that the effect of status is significant? R reports only t-values and not p-values Thanks a lot


title: Survey weights in logistic regression models

text: We are interested in fitting a multiple logistic regression model using data obtained from a national survey of children with special health care needs. The data has an accompanying weight variable intended to standardize children to the national population in which we intend to make inference. This weight variable does not sum to 1 nor are the weights integral (they may take values such as 23.2). This model is being fit using SAS v9.2. In consulting the documentation for the logistic procedure, I notice in the syntax description the following statement: > Caution: PROC LOGISTIC does not compute the proper variance estimators if > you are analyzing survey data and specifying the sampling weights through > the WEIGHT statement. The SURVEYLOGISTIC procedure is designed to perform > the necessary, and correct, computations. I don't understand why this should be an issue. If model based standard errors are being computed, then the weighted maximum likelihood estimator should give standard errors which are correct for the population of interest. Is this correct? What likelihood function is SAS's logistic regression solver optimizing if the above statement is correct?


title: Combination of normal and gamma distribution in R

text: I have to solve this problem: _A machine comes with a 2 year warranty. As a statistician you are told that (i) The total number of hours for which a machine is used over the course of 2 years follows a N(10, 5) distribution. (ii) The number of hours for which a machine motor runs until the first fault occurs has a Gamma(0.006, 7) distribution. Moreover, the number of hours of operation before the first fault is independent of how frequently the machine is used. Run a simulation study to estimate the proportion of machines whose motors will develop a fault during their 2 year warranty._ Actually, I don't understand what exactly the exercise asks me to do. Could you please explain it to me? Thank you


title: Why does HoltWinters() in R behave differently when parameters are set as 0 vs set as False?

text: Why is it that the `HoltWinters` function in R gives _vastly_ different parameter estimates, say of `alpha` and `beta`, when `gamma` is set as `0` versus set as `False`? It would seem from the structure of the model the 2 should be specifying equivalent models?


title: Best optimization package for employee scheduling problem?

text: I am looking to solve the optimization problem described below. Which optimization software package would be best suited for this, considering the requirements specified below? **Requirements:** 1) Can solve the complex problem described below. 2) The package must be free of charge and ideally open source. 3) Ideally no heavy reading on the model syntax required / fast learning curve. 4) Ideally the model can be run using Python or Excel as front ends. 5) Ideally a visual front end is available such that the end users (with zero optimization knowledge) can easily modify the optimization constraints. **Optimization problem:** The optimizer should find a monthly work schedule for the 50 employees of a hospital which satisfies all of the constraints below. For each calendar day of the relevant calendar month, the optimizer determines which of the 50 employees are scheduled to work (=1) and which ones are not (=0). So there are 30 days * 50 employees = 1500 variables to solve for, each of which can either be 0 or 1. It is not required to maximize or minimize a certain variable; finding a solution which satisfies all the constraints is sufficient. **Constraints:** 1) Each employee has a prespecified min. and max. number of work days per month. 2) No more than 4 consecutive work days for each employee. 3) At least 2 free days in between work days for each employee. 4) Some employees can only work on some of the days of the month (for example only Monday to Wednesday). These constraints are prespecified. 5) Each employee is either a doctor or a nurse. On each calendar day, a prespecified minimum number of doctors and nurses must be present. 6) At least 2 weekends (=Sat and Sun) per month must be free for each employee.


title: unbiased sampling of expecation over maximization operator

text: My problem setting is as follows, I have a set of M random variable               X = \{ {X_1},{X_2},...{X_M}\}      where each variable X_i is estimated via stochastic sequence u^t_i, aka, u^t_i is the unbiased estimator for X_i, {X_i} = {E_t}[u{_{}^t_i}] And my task is to estimate {\max _i}E_t[{X_i}]. If at time t, I compare all the u_i over i, this is actually doing ${\max _i}u_{_i}^t$ at round t, however, this is not the unbiased estimator of {\max _i}E_t[{X_i}], since {E_t}[{\max _i}u_ {_i}^t]! = {\max _i}E[{X_i}] We can see that the sequence of the expectation operator and max operator is flipped. BTW, how to post equations on CrossValidated？ I previously use metaoptimize and not familiar with CV. Thanks!


title: Estimating the parameters of a beta distribution with zeroes and ones in the sample

text: I have a list of values in [0,1] that I want to fit to a Beta distribution in order to get the corresponding alpha parameter. I can't use a beta fitting function because my values might be 0's and 1's. On wikipedia I found this formula for the method-of-moments estimation of alpha given the sample mean and sample variance: ![enter image description here](http://i.stack.imgur.com/dw4pi.png) I don't need any exact fitting, just a good approximation. By using sample mean and variance, I should get rid of all the problems arising with 0s and 1s in my sample, right? Will this still be acceptable?


title: Treating numeric as categorical variable in regression

text: I need a little bit of help and confirmation that I have the right idea. I have some fake data of 8 tribes; within each tribe members work hard to gain food for their own tribe. No one can speak to these tribes, but people suspect that each tribe has one of the two strategies presented below for gaining food:   1. Members of a tribe who travel farther away from the tribe's main location are given more food, so they face less of a chance of starving before coming back.   2. Members of a tribe who travel far are given less food; that way if they are lost, there is less of a food loss to the tribe as a whole. The variables (columns) I am working with are `Tribe number`, `Distance from the tribe location` when sample was taken (10, 20, 30, or 40 miles), `Weight of each member` that we are studying (related to the amount of food is given), `height` (taller people use energy more efficiently, and there is a strong positive correlation between weight and height in arbitrary units and inches), finally I have each observation categorized by height (group `1`: 56–62in, group `2`: 62–64...). I want to find out if the tribes use different strategies, and also if there is a difference among the classes `pf` height. In addition I want to find out the strategies that are in use. I am having a hard time with how to classify each tribe as using either strategy `1` or `2`. I was thinking of doing a one- way ANOVA to check if there is a difference in `mean` within each group based on `distance`. (In a particular tribe is there a difference in the mean of weight between those who were 10, 20 , 30, or 40 miles?) I don't know how to figure out if each colony uses a different strategy. Finally, I want to build a linear model of mass on colony, distance, and height. I know how to build a model and run diagnostics. My concern here is, can I use distance as a categorical variable since its values are 10, 20, 30, or 40 miles?


title: Regression Discontinuity Group Level Program Evaluation

text: I'd like to use a regression discontinuity design to evaluate a program where the discontinuity/assignment to treatment occurs at the group level. However, I'd like to measure the outcome at the individual level (as opposed to measuring the outcome as an average at the group level). I have three questions (feel free to answer any):   1. What are the criteria/assumptions for doing an RD with assignment at the group level?   2. Would using HLM be appropriate?   3. Do you know of any studies that do this? For example, consider a very large firm where workers make widgets independently, but are divided into teams where they can consult with each other or ask for advice a few times a week. The firm calculates how long it takes each worker to make a widget on average, and then calculates an average widgets per hour for each team. They provide additional training for groups where the group average is below a certain threshold and no additional training for groups above the threshold, hence the discontinuity. The firm then wants to evaluate the effect of additional training on an individual.


title: Random forests vs boosting

text: I thought it would be interesting to talk about two of the best ensemble methods off-the-shelf: **Random Forests** and **Boosting**.   * When would you apply one method rather than the other one?


title: Irregularly spaced time-series in finance/economics research

text: In financial econometrics research, it is very common to investigate relationships between financial time series that take the form of **daily data**. The variable will often be made $I(0)$ by taking the log difference, for example; $\ln(P_t)-\ln(P_{t-1})$. However, daily data means that there's $5$ data points each week, and Saturday and Sunday are missing. This seems to get no mention in the applied literature that I'm aware of. Here's some closely related questions that I have that come from this observation:   * Does this qualify as irregularly spaced data, even though financial markets are closed over the weekend?   * If so, what are the consequences for the validity of extant empirical results garnered thus far in the gigantic number of papers that ignore this issue?


title: What is the income distribution in the U.S., in 1% increments?

text: Where can I find a roughly complete CDF of income in the U.S. (not simply the Gini coefficient or a point summary like "X% of the population make Y% of the income?") _Note: I think this may be off-topic here and have proposed aseparate QA site for these sorts of questions, but someone else suggested that these types of questions are on-topic here._