title: How do I get a $p$-value from the Cochran-Armitage trend test?

text: So, I'm working with GWAS SNP data and want to perform several tests for association between genotype and phenotype. There are two phenotypes (case and control) and 2 or three genotypes. Most of them are Chi-squared tests with different contingency tables, $2 \times 2$ or $2 \times 3$, one of them is the Cochran-Armitage trend test (CATT) Once I have constructed the contingency table, I can easily get a $p$-value using the Apache commons math library for the Chi-squared tests. No problem. However, the explanation of the CATT on Wikipedia is not sufficient for me to implement it (my statistics knowledge is limited and I'm still learning). Like in the example, I suspect a linear trend, so my weights are $t = (0,1,2)$, which make the formula for $T$ to: $$ T \equiv (N_{12}R_2 - N_{22}R_1) + 2(N_{13}R2 - N_{23}R1) $$ and the one for the variance $$ Var(T) = {{R_1 R_2} \over N} ( N(C_2+4C_3) - (C_2 - 2C_3)^2) $$ I checked how the program PLINK does it, since it's already implemented there, but it differs slightly from the above formulas. The C++ source code there would correspond to this: $$ T = {(N_{12}R_2 - N_{22}R_1) + 2(N_{13}R2 - N_{23}R1)\over N} $$ and $$ Var(T) = {{R_1 R_2} \over N} {( N(C_2+4C_3) - (C_2 - 2C_3)^2) \over N^2} $$ Then it does calculates a chi-square value like this $$ \chi^2_{T} = {T^2 \over Var(T)} $$ and calculates the $p$-value like for any other chi-squared value with $df = 1$ I don't need to understand the theory completely, as long as my program calculates correctly, but understanding it would give me additional confidence. Is this correct or legitimite? Is this how I'll get the $p$-value?


title: What is the name of a plot of - log 10 of p values in multiple testing?

text: Under the null, the p-value is uniformly distributed between 0 and 1. Taking the negative of the log base 10 of many such p-values should follow an exponential distribution. You can show quantiles of their distribution under the null, but more importantly, logging them makes large differences very compelling. I can't remember the name of such a plot though: ![enter image description here](http://i.stack.imgur.com/eNhR2.jpg)               ## example R code     x <- replicate(100, rnorm(10), simplify=FALSE)     y <- replicate(100, rnorm(10), simplify=FALSE)          lm.sig <- function(x, y) {       coef(summary(lm(y ~ x)))[2, 4]     }          sigs <- mapply(lm.sig, x, y)          plot(-log(sigs, base=10), type='h', xlab='Comparison Index', ylab='-log p value')     abline(h=-log(0.05, base=10), lty=2)     legend('topleft', lty=2, 'Nominal 0.05 error rate')


title: K - means cluster always landing right on top of whole dataset mean

text: I have a so so sized data set - 30 000 observations. I would like to run K-means on them but to restrict the center(mean) of the data. This is, I would like to push the clusters away from this mean. As I have noticed that independantly of the # of clusters, one ends up landing right on top of the mean of all the variables, like a smaller version of the whole data set. Is there a way to restrict K-means to not behave this way? By the way I have two different ways to initialize cluster centers, random starting points and means of random samples of the data. I run each 500 times and the solution seems stable enough but maybe k-means++ would have a different outcome? I wouldnt think so... thank you in advanced reader.


title: Predicting Next Likely Outcome of Binary Time Series?

text: I'm trying to approach the following problem: Danny & Johnny are professional basketball players. Each day they meet, and play for a while. Whoever scores the most points is declared winner for the day, and wins a dollar. After a year of games, we are given a vector (including 365 elements) representing, say, Danny's winning days: [1, 0, 0, 1, 1, 0, 0, .......] Given the vector of prior 365 'basketball matches' between the two players, what methods can one use to try and predict the winner of the next (366th) game? Would greatly appreciate any insights / suggestions.


title: Box Cox transformed my data, now how to use it in my mixed model?

text: Sorry for the (most likely) simple question, but I have Box Cox transformed my data in SAS, but I am unsure how to use the transformed data in my mixed model. Do I output it to a new file, and use that in my model statement, or do I need to incorporate the proc transreg statement into my proc mixed? Details are appreciated as I am new to SAS. Thanks, Maggie


title: Cross-validated bootstrap samples

text: In one of the paper I am writing, I am looking at cross-validation on bootstrap samples. I wrote the following explanation.One of the reviewer wrote that he didn't understand where is the correlation. I don't know how to explain better. So I was wondering if someone has a reference I could cite about this issue when mixing cross-validation and bootstrap. Thanks! > > In the context of cross-validation, we have to be careful. Cross- > validation divides the data in two parts, the training set and the test set. > For the original data, every individual observation appears once in the > dataset (we assume all the (X, Y ) p + 1 tuples are distinct.) Hence there > is no possibility that an individual is used in the training and the test > set, a potential source of positive bias. This is not the case for a > bootstrapped dataset — there will almost certainly be ties. If we let some > of these tied values go to the training set, and some to the validation set, > we will have artificially created (additional) correlation. To avoid this, > we operate at the level of observation weights. In the original sample, each > observation has weight 1/n. In a bootstrapped sample, these weights will be > of the form k/n, for k = 0, 1, 2, . . .. When we cross-validate a bootstrap > sample, we randomly divide the original observations into the two groups, > but their bootstrap weights go with them (including the 0s).


title: sample size and the validity of prediction for each category in classification

text: This is a problem of classification. If I have 100 training data of category A and there are some features in A. Suppose feature 1 always shows up with category A. But there is only one training data for category B, with feature 2 in it. We also have many data for category C. Feature 1 and feature 2 never show up in C in the training data. Now we have a test data, with feature 1 in it and another test data with feature 2 in it. Could you please come up with a method to give the validity of the prediction for both test data? How to evaluate feature 1 and 2? Is there a way to calculate the size of the training data to guarantee a certain precision. Thank you!


title: $\arctan(-b/a)$ when $a$ and $b$ are i.i.d. normal random variables

text: If $a$ and $b$ are both independent, identically distributed normal random variables, what is the distribution of $\phi = \arctan{(-\frac{b}{a})}$? It is uniform, across $0$ to $2\pi$, but I am trying to see why. $a \sim \mathcal{N}(0,\sigma^2)$ and $b$ is identically distributed and independent to $a$. For example, setting $X = \arctan{(-\frac{b}{a})}$, my strategy is to first calculate the cumulative distribution function (CDF), and then compare the CDF to known CDFs. $$F_X(x) = P(X\leq x) = P\left(\arctan{\left(-\frac{b}{a}\right)}\leq x\right) = P\left(-\frac{b}{a}\leq\tan(x)\right)$$ Now, $a$ could be either positive or negative, since $a$ is a normal random variable with mean $0$ and non-zero variance. Then, $$\begin{aligned}P\left(-\frac{b}{a}\leq\tan(x)\right) & = P\left(-\frac{b}{a}\leq\tan(x)|a>0\right)P(a>0)\\ & +P\left(-\frac{b}{a}\leq\tan(x)|a\leq 0\right) P(a\leq 0) \end{aligned}$$ Each of the terms $P\left(-\frac{b}{a}\leq\tan(x)|a\leq 0\right) $ and $P\left(-\frac{b}{a}\leq\tan(x)|a > 0\right)$ have terms depending on the value of $x$. I am drawing a plot of $a$ on the x-axis and $b$ on the y-axis and trying to use polar coordinates to calculate the integral using the joint probability function of $a$ and $b$, which is simply the product of the two distributions. Can I have a hint on how to proceed? If I ignore the value of $a$ entirely (which is wrong) and ignore one of the probability terms (which is also wrong) and ignore the value of $x$, I am getting a result very similar to $\frac{x}{2\pi}$ (which is desirable). Thanks.


title: Procedure for variable selection + logistic regression when n is small, p is large, and data are unbalanced?

text: I have data that have been collected using case-control procedures, in which the population of positive cases is collected with a random sample of negative cases. This yields 62 positive cases and 179 controls. There are 58 possible predictor variables (mostly numeric, two factors). My goal is not classification per se--I am in the social sciences and significance testing is a big deal in my discipline, for better or worse. Rather, my collaborators want to make some inferences about the predictors that increase the probability of being a positive case. Obviously, the data are unbalanced, non-randomly sampled, and I have low power. I am trying to select variables for a logistic regression. As an exploratory option, I have run a random forest with various specifications to try and find important variables. My false positive rate is, as expected, only about .03, as there are many more positive cases than negative. My false negative rate is quite high, about .70. Adjusting the weights doesn't seem to help with this. What would be a good procedure for variable selection that could then be used to estimate rare events logistic regressions using the reduced set of variables? Thanks for any help.


title: Problems with a simulation study of the repeated experiments explanation of a 95% confidence interval - where am I going wrong?

text: I'm trying to write an R script to simulate the repeated experiments interpretation of a 95% confidence interval. I've found that it overestimates the proportion of times in which the true population value of a proportion is contained within the sample's 95% CI. Not a big difference - about 96% vs 95% but this interested me nonetheless. My function takes a sample `samp_n` from a Bernoulli distribution with probability `pop_p`, and then calculates a 95% confidence interval with `prop.test()` using continuity correction, or more exactly with `binom.test()`. It returns 1 if the true population proportion `pop_p` is contained within the 95% CI. I've written two functions, one which uses `prop.test()` and one which uses `binom.test()` and have had similar results with both:               in_conf_int_normal <- function(pop_p = 0.3, samp_n = 1000, correct = T){         ## uses normal approximation to calculate confidence interval         ## returns 1 if the CI contain the pop proportion         ## returns 0 otherwise         samp <- rbinom(samp_n, 1, pop_p)         pt_result <- prop.test(length(which(samp == 1)), samp_n)         lb <- pt_result$conf.int[1]             ub <- pt_result$conf.int[2]         if(pop_p < ub & pop_p > lb){             return(1)         } else {         return(0)         }     }     in_conf_int_binom <- function(pop_p = 0.3, samp_n = 1000, correct = T){         ## uses Clopper and Pearson method         ## returns 1 if the CI contain the pop proportion         ## returns 0 otherwise         samp <- rbinom(samp_n, 1, pop_p)         pt_result <- binom.test(length(which(samp == 1)), samp_n)         lb <- pt_result$conf.int[1]             ub <- pt_result$conf.int[2]          if(pop_p < ub & pop_p > lb){             return(1)         } else {         return(0)         }      }      I've found that when you repeat the experiment a few thousand times, the proportion of times when the `pop_p` is within the 95% CI of the sample is closer to 0.96 rather than 0.95.               set.seed(1234)     times = 10000     results <- replicate(times, in_conf_int_binom())     sum(results) / times     [1] 0.9562      My thoughts so far about why this may be the case are   * my code is wrong (but I've checked it a lot)   * I initially thought that this was due to the normal approximation issue, but then found `binom.test()` Any suggestions?


title: How to determine voter preference from a "small sample" poll using Bayesian statistics?

text: I created a poll on a sports forum asking people how they think their favorite NBA team did this off-season. The poll options were: Top 10, Middle 10, Bottom 10. I expect dozens of people to respond, which isn't a huge sample. The hypothesis I'm testing (unbeknownst to the voters) is that there will be a "homerism" bias which makes the "Top 10" vote more likely. In terms of Bayesian methodology, I would like to compare two models, one has preference weighted uniformly across all three poll choices, and the other has the weight heavily biased towards the "Top 10" choice. So, my question is how to create the prior and the likelihood functions for this type of problem. I don't need full details, but if you could point me in the right direction, I'd appreciate it.


title: How to write a poker player using Bayes networks

text: This is my first question on stackexchange and also my first time implementing a Bayesian network so I will apologize ahead of time for any novice mistakes I make. The goal of my project is to implement a Poker player that does Bayesian inference. There has been some work done on this by a group at Monash University in Australia led by Kevin Korb which I am using to reference off of. You can find their work here and here. The first reference, being a book, is more helpful and detailed (see Ch.5.5 & 11 for Poker). I am using the simplified version of Texas Hold'em called Leduc Hold'em to start. Leduc Hold’em is a two player poker game. The deck used in Leduc Hold’em contains six cards, two jacks, two queens and two kings, and is shuffled prior to playing a hand. At the beginning of a hand, each player pays a one chip ante to the pot and receives one private card. A round of betting then takes place starting with player one. After the round of betting, a single public card is revealed from the deck, which both players use to construct their hand. This card is called the flop. Another round of betting occurs after the flop, again starting with player one, and then a showdown takes place. At a showdown, if either player has paired their private card with the public card they win all the chips in the pot. In the event neither player pairs, the player with the higher card is declared the winner. The players split the money in the pot if they have the same private card. Each betting round follows the same format. The first player to act has the option to check or bet. When betting the player adds chips into the pot and action moves to the other player. When a player faces a bet, they have the option to fold, call or raise. When folding, a player forfeits the hand and all the money in the pot is awarded to the opposing player. When calling, a player places enough chips into the pot to match the bet faced and the betting round is concluded. When raising, the player must put more chips into the pot than the current bet faced and action moves to the opposing player. If the first player checks initially, the second player may check to conclude the betting round or bet. In Leduc Hold’em there is a limit of one bet and one raise per round. The bets and raises are of a fixed size. This size is two chips in the first betting round and four chips in the second. If you look at pg.185, Section 5.5.2.1 in Figure 5.14 there is a diagram for a Bayes Net for Poker. ![Bayesian network for a poker player](http://i.stack.imgur.com/f2xNV.png) This is essentially the same one I am using for my project, but granted that there are no up-cards in Leduc Hold'em the two corresponding nodes for up- cards are not applicable for it. I was able to compute the conditional porbability tables between the node pairs (BPP_Win, BPP_Fin, OPP_Fin), (OPP_Fin, BPP_Fin), (OPP_Fin, OPP_Curr), and (BPP_Fin,BPP_Curr) **but I am not sure how to compute the conditional probability of (OPP_Curr, OPP_Action)**.


title: Modeling a spline over time -- design matrix and survey of approaches

text: A response variable y is a nonlinear function of a number of predictor variables X (in my real data the response is binomially distributed, but here I'm using a normally-distributed value for simplicity). I can model the relationships between the predictors and the response using splines / smooths (for example, GAM models in the `mgcv` package in R). So far so good. However, each response is the result of processes that evolve over time. That is, the relationship between the predictors X and the response y changes over time. For _each_ response, I have data for the predictors over a number of time points around the response. That is, there is one response per group of time points (not that the _response_ evolves over time). Some illustrations are probably helpful at this point. Here is some data with known parameters (code below) and then plotted using ggplot2 (specifying the GAM method and an appropriate smoother), with time in the facets. For illustration, y is a quadratic function of x1, and the sign and magnitude of this relationship changes as a function of time. ![EDIT: since I'm new I can't post an image; run the code to see it.](http://i.stack.imgur.com/WzCDF.png) The relationship between x2 and y is circular, corresponding to an increase in y with a certain _direction_ of x2. The amplitude of this relationship modulates over time. (Modeled in ggplot using a gam specifying a "cc" circular cubic smoother). ![EDIT: since I'm new I can't post an image; run the code to see it.](http://i.stack.imgur.com/UV2TM.png) **I would like to model the (nonlinear) change in each predictor as a function of time using something like a two-dimensional spline.** I have considered using a two-dimensional smooth in the mgcv package (something like `te(x1,t)`), except that this would require the data in a long-form (i.e. a single column of time points). I think this is inappropriate, because _one_ response is associated with _all_ time points -- so arranging the data in long form (thus duplicating the same response across multiple rows of the design matrix) would violate the independence of observations. My data are currently arranged with columns `(y, x1.t1, x1.t2, x1.t3, ..., x2.t1, x2.t2, ...)` and I think this is the most appropriate format. I would like to know:   1. is there a better way to model this data   2. if so, what the design matrix / formula of the model would look like. Ultimately I would like to estimate model coefficients using Bayesian inference in an mcmc package like JAGS, so I would like to know how to _write_ a two-dimensional spline. R code to reproduce my example:               library(ggplot2)     library(mgcv)     #-------------------     # start by generating some data with known relationships between two variables,     # one periodic, over time.     set.seed(123)          nTimeBins <- 6     nSamples <- 500          # the relationship between x1, x2 and y are not linear.     # y = 0.4*x1^2 -1.2*x1 + 0.4*sin(x2) + 1.2*cos(x2)          # the relationship between x1, x2 and y evolve over time.      x1.timeMult <- cos(seq(-pi,pi,length=nTimeBins))     x2.timeMult <- cos(seq(-pi/2,pi/2,length=nTimeBins))          qplot(x=1:nTimeBins,y=x1.timeMult,geom="line") +            geom_line(aes(x=1:nTimeBins,y=x2.timeMult,colour="red")) +            guides(colour=FALSE) + ylab("multiplier")          df <- data.frame(setup=rep(NA,times=nSamples))     for (time in 1 : nTimeBins){       text <- paste('df$x1.t',time,' <- runif(nSamples,min=-3,max=3)',sep="")       eval(parse(text=text))            text <- paste('df$x2.t',time,' <- runif(nSamples,min=-pi,max=pi)',sep="")           eval(parse(text=text))         }         df$setup <- NULL          # each y is a function of x over time.     text <- 'y <- '     # replicated from above for reference:     # y = 0.4*x1^2 -1.2*x1 + 0.4*sin(x2) + 1.2*cos(x2)     for (time in 1 : nTimeBins){       text <- paste(text,'(0.4*x1.t',time,'^2-1.2*x1.t',time,') *                      x1.timeMult[',time,'] + (0.4*sin(x2.t',time,') +                      1.2*cos(x2.t',time,'))*x2.timeMult[',time,'] + ',sep="")     }          text <- paste(text,'rnorm(nSamples,sd=0.2)')     attach(df)     eval(parse(text=text))     df$y <- y          #-------------------     # transform into long form data for plotting:     df.long <- data.frame(y=rep(df$y,times=nTimeBins))         textX1 <- 'df.long$x1 <- c('     textX2 <- 'df.long$x2 <- c('          for (time in 1:nTimeBins){         textX1 <- paste(textX1,'x1.t',time,',',sep="")         textX2 <- paste(textX2,'x2.t',time,',',sep="")     }     textX1 <- paste(textX1,'NULL)',sep="")     textX2 <- paste(textX2,'NULL)',sep="")     eval(parse(text=textX1))     eval(parse(text=textX2))     # time stamp:     df.long$t <- factor(rep(1:nTimeBins,each=nSamples))          #-------------------     # plot relationships over time using GAM fits in ggplot:     p1 <- ggplot(df.long,aes(x=x1,y=y)) + geom_point() +                   stat_smooth(method="gam",formula=y ~ s(x,bs="cs",k=4)) +                   facet_wrap(~ t, ncol=3) + opts(title="x1 versus y over time")     p1          p2 <- ggplot(df.long,aes(x=x2,y=y)) + geom_point() +                   stat_smooth(method="gam",formula=y ~ s(x,bs="cc",k=5)) +                   facet_wrap(~ t, ncol=3) + opts(title="x2 versus y over time")     p2


title: Most appropriate statistical test of significance when comparing two groups of awareness data

text: I have 2 sets of brand preference data (in %) to measure impact of online advertising on target group F20-40 years i.e., 1.different time period (pre/post - not same consumers but same target group) & 2. Control (not seen online advertising) vs Exposed (exposed to advertising). Could you suggest what statistical significance test I should use and why?


title: Linear statistic for measuring association in the party R package

text: I had asked a general question about conditional inference trees via party a while back and gotten a great reply. I am revisiting this procedure and trying to make sense of the linear statistic that is being used (Hothorn et al., Unbiased Recursive Partitioning: A Conditional Inference Framework, Research Report Series, 2004, page 4, equation (1)). ![enter image description here](http://i.stack.imgur.com/g5uWU.jpg) I am not at all clear how this statistic is calculated. Can anyone help? Here is what it seems to me that if   * $g(\cdot)$ and $h(\cdot)$ are the identity functions   * and there is a predictor $x$ and a response $y$, both numeric, the statistic is simply a scalar. This is not what is shown, so I am wrong :) Example data:               x<-c(1,3,4,67,32,23,3,12,4)     y<-c(43,23,45,22,12,465,6,54,3)     w<-rep(1,9)          T<-0     for (i in 1:length(y))     {            T<-T+(w[i]*x[i]*y[i])            }     T #13523      I _THINK_ the result needs to be a vector of length 81.


title: Arima model giving high forecast values

text: I have some models built with the `auto.arima` function from the `forecast` package. I'm modeling a variable called 'natural efluent energy' (ena), which is how much energy you can extract from some Hydrography region. There are 2 regressor variables (rainfall precipitation from period $t$ and $t-1$.) Each region has it's own model - some series show positive trend, some shows negative trend, and some seems stationary. The problem is that some forecasts 'from `auto.arima`' are giving values higher/lower than usual (some forecasts give me negative values, which are not possible). My original call is below:                m1 = auto.arima(serie, xreg = regvars)      For the data on the link, I changed it to                m1 = auto.arima(serie, xreg = regvars, max.P = 0, max.Q = 0, stationary = TRUE)      Then I get good forecasts in this case. My question is, what these parameters(`max.P`, `max.Q`) actually control, and how they relate to the trend show by my model variable? Here is a link for the historic data: http://www.datafilehost.com/download-7718b3fc.html And here a link for the forecast regressors: http://www.datafilehost.com/download-ca44dfa4.html And here a link of mean historic values, the forecast must fall between these values: http://www.datafilehost.com/download-e1e265b7.html My data starts at 2001/Jun, so the serie is:                 y = ts(dframe$ena, freq = 12, start = c(2001, 6))


title: Exploratory factor analysis among items with different measurements

text: I intend to do an exploratory factor analysis among items that have different types of measurements: for example, some items have bounded scales (1-100%), some items take small and other large values (different units indeed, such as dollars or hours). I've been reading a lot and haven't found a negative answer, i.e. that it is not possible to do such analysis with such varied types of items; however, every example I came across is based on items that take the same values (e.g. 4-value likert scales from questionnaires). In sum, does anyone have a clear answer (and reference?) on this matter?


title: Issues with stochastic gradient descent

text: I am using stochastic gradient descent to learn a model. Here is the plot of the objective function for the iterations. I am trying to maximize the function value. ![enter image description here](http://i.stack.imgur.com/3Eq3T.jpg) Taking the average of 500 iterations, I have this next plot ![enter image description here](http://i.stack.imgur.com/vyo1j.jpg) As you can see the function value is increasing and we can say the algorithm is converging. However, looking at the first plot, we can see that the function never actually went past a certain threshold. Even though the minimum value it could reach went on converging with the iterations, the maximum value that it could attain remained almost the same. So can we say the algorithm is converging?


title: R caret difference between ROC curve and accuracy for classification

text: In case of _caret_ package test function metric option, one can use either accuracy or ROC as a metric that will be used to finalize values of tuning parameters. I felt that accuracy and ROC are the same Accuracy = total correct predictions/total samples ROC= looks at various probability cutoffs and gives probability cutoff at which accuracy will be the best. Am I correct?


title: How to work out Effect Size for a MANOVA using G*Power

text: I'm using G*Power to work out how many participants I need. I'll be using a MANOVA, I have two independent variables (Male and Female) and they will be answering 13 different questionnaires so 13 dependent variables. I put this data into G*Power and its telling me I only need 14 participants but I think that is too small?


title: Regression of proportions

text: I have to regress the mortality rate which is in number of deaths per thousand and percentage of vaccinated people. Have sample data for different districts. Can I use linear regression directly. The sample sizes are different for different districts. Please suggest the appropriate method.


title: Predicting Future Task Execution Times from Historical Data

text: In a manufacturing program I am collecting the estimated time required and the actual time taken for a series of tasks. Suppose I find that I often take twice as long, say, as I originally estimated for each task, I want to apply a typical difference between estimated and actual time taken to future planned tasks to improve the accuracy of my planning. Starting with no knowledge of stats, I'm not sure what to google or how to phrase the question - but so far I have thought a) determine a mean variance, positive and negative, for all past tasks between estimated and actual times and apply this value to future tasks to adjust the expected max and min times. b) determine the standard deviation of past task timing, and apply forward. Would appreciate any pointers- thanks.


title: Identifying a block design

text: In designing a email test, I have three binary factors (each have two levels, say low-high). Our customers are divided into 5 non-overlapping groups. Is it best to run a randomized block design, blocking out the customer groups (5 blocks, 2^3 factorial) or use the customer groups as an actual 5-level factor?


title: Survfit function in R to score a new dataset

text: I have built a cox proportional hazards model in the R survival package. I want to score new data set using this model. I thought the survfit function would doing this using survfit(original model, newdata), but it only seems to predict only the outcomes for the original data, not the original model using the new data. Any insight is appreciated!


title: Determining number of clusters with maximum jump or minimum variance

text: I want to determine a number of clusters, but I don't know how to do it. I want to use kmeans and select a number of cluster with minimum variance or with a maximum jump in variance. I plot variance; how can I do it with programming (Java)?


title: Including covariates increases spatial patterning rather than reducing it

text: I am fitting a series of Bayesian spatial models with an intrinsic conditional autoregressive (iCAR) structure. In these models, I begin by fitting an empty model (null) and then add successive sets of covariates based on a conceptual foundation. I then map quantiles of the mean posterior spatial random effects, highlighting those with 95% credible intervals that do not contain zero. Results show that including additional covariates _increases_ the number of areal units with 95% credible intervals that do not contain zero. I have not seen this before in spatial models – rather I typically see the number of areal units with significant exceedance probabilities _decrease_ when controlling for more variables. I am having difficulty explaining why controlling for more variables causes the spatial pattern to become more pronouced and am wondering if others have come across this in their own work.


title: Using percentiles as predictors - good idea?

text: I am thinking about a problem which is to predict log(spend) of a customer using linear regression. I am considering what features to use as input and wondering if it would be OK to use the percentile of a variable as inputs. For example I could use the companies revenue as a input. What I'm wondering is whether I could use the company revenue percentile instead. Another example would be a categorical industry classifier (NAICS) - if I were to look at median spend per NAICS code and then assign each NAICS code to a 'NAICS Percentile', would that be a valid explanatory variable I could use? Just wondering if there are any issues to be aware of when using percentiles? Is it in some ways equivalent to a type of feature scaling?


title: Simple linear regression on constrained variables

text: I have been asked whether there is any relationship between two variables $X$ and $Y$. I have done a simple regression and the confidence intervals suggest a small significant relationship. However, these variables are nonnegative and most of the points are "close" to (0,0) --- in the sense that in a graph of the points partitioned into a grid pattern most of the points will be within the square with the "south west" corner at (0,0). Does this violate the assumptions of regression? I believe "yes" because the errors in the regression might not be normally distributed (they would be bounded for a given $x\in X$). How can I best answer the question about the relationship between $X$ and $Y$? Should I just use simple regression or are there ways to get around this possible violation? My data looks a bit like                set.seed(5)      x <- abs(rcauchy(1000,0,1))      set.seed(5)      y <- rexp(1000,5)      plot(x,y)      though there is no reason to believe that the underlying distributions of either are cauchy or exponential.


title: How to test whether two means of categorisation generate significantly different results

text: If I have the following contingency table in which `Novel.Cat` is a new means of classifying some thing. `Ref.Cat` is the gold standard. How can I do some sort of test to see whether the two means of classification generate significantly different results? Is just a simple Fisher's test appropriate? Or is there some other, simple way?                         Novel.Cat      Ref.Cat  a  b  c  d             a 15 2  0  1             b 9  9  0  0             c 5  1  5  0             d 0  0  0  6      With a Fisher's test, all I am testing is the non-independence between the groups. However, what I want to specifically test is the concordance between `Novel.cat` categorisation and `Ref.cat` categorisation. So, the ideal (best concordance) contingency table would be:                         Novel.Cat      Ref.Cat  a   b   c   d             a 18  0   0   0             b 0   18  0   0             c 0   0   11  0             d 0   0   0   6      Fisher's however, would say that the following results are significant (i.e. `Novel.cat` and `Ref.cat` are non-independent)                         Novel.Cat      Ref.Cat  a   b   c   d             a 0   0   0   18             b 0   0  18   0             c 11  0   0   0             d 0   6   0   6      However, quite clearly, `Novel.cat` has misclassified everything if `Ref.cat` is to be taken as the gold standard.


title: How can you detect if a Gaussian process is over-fitting?

text: I am training a Gaussian process with an ARD kernel with lots of parameters by maximizing the marginal lielihood of the data, instead of cross-validation. I suspect that it is over-fitting. How can I test this suspicion in a Bayesian context?


title: High threshold coefficients

text: Consider the following code/output:                      library(ordinal)            clm(as.factor(y)~x, link="probit", data=data)                  0|1         1|2         2|3      x         59.5369     90.9923    112.1025    138.1162                 data          x    y          1    0          2    0          3    1          4    1          5    2          6    2          7    3          8    3          9    3         10    3      Why are the threhold coefficients so high? Shouldn't the threshold coefficients be approximately between $1$ and $10$? For example, $x$ is categorized as a $0$ if $0 \leq x^{*} \leq 2$ and a $1$ if $3 \leq x^{*} \leq 4$


title: Are probability graphic models useful for predictive modelling?

text: Are Probability Graphic Models (say specifically Bayesian Networks) useful for predictive modelling in terms of large data (100,000 - 1,000,000 rows) and many variables (hundreds)? Meaning, is this technique / methodology something that can be a competitor for, say random forests or some machine learning method for classification/regression?


title: Adding a quadratic term: should I use the squared original (and not the squared standardized?)

text: In a multiple logistic regression I need to standardize one of the variables because I need to add a quadratic term. Whether I add the quadratic term as the squared original or the squared standardized, I get very similar models, same AIC. Why? The estimates of the linear term change. When I square a standardized variable, I get positive numbers after squaring half of the data (which are negative) and that should completely mess up my data. Why does it not? Can someone please explain the mathematical reason behind it and tell me if I definitely should use the squared original and if I can use the squared standardized, why does it work? I notice this is a mathematics question rather than anything but I searched and could not find anything related to this.


title: How to select the best variables by RandomForest in R?

text: I have a table of `mRNA` levels of my target gene and it's transcription factors in many different condition. What I want to do is to select the most important conditions and most important transcriptional factors. I used `glm` in `R` but I am stuck at the point that I can't extract the important data from it now. Could anybody help me please? I would like to use also `Lasso`, `Elastic Net` and `Random Forest`. I've tried this:               data <- read.table("table.txt",header=T)      genenames=colnames(data)      attach(data_test)      f1=as.formula(paste(genenames[5],"~",genenames[1]))      glm(f1)      glm.linear <- glm(f1)      glm.linear                             V1  V2  V3      V4    V5     1   3.85  3.85  9.41  4.31  9.01      2   3.98  4.01  8.50  3.95  7.35     3   3.98  3.98  8.17  4.07  8.58     4   3.61  4.10  7.94  3.88  4.59     5   3.74  3.95  9.55  4.08  7.79     7   3.78  3.98  7.72  3.91  4.73     8   3.76  4.05  7.60  3.71  4.49     9   3.75  4.08  7.71  3.76  4.73     10  3.81  4.19  8.03  3.75  4.53


title: Conditioning on a variable

text: I am looking at some software code that performs conditioning on random variables. For example, one can have a set of random variables which have a multivariate normal distribution associated with them and then you can condition on a given variable to take on a certain value and then get the associated conditional distribution for the rest of the variables. Now, I notice that distributions which are univariate do not have this method implemented. But for the sake of completeness, for example for a univariate Gaussian when we actually observe that the variable takes on a certain value, shouldn't the conditional distribution be a point mass with zero variance?


title: Central limit theorem with scilab

text: I'm trying to illustrate the CLT with scilab but my results are weird. Did I make a mistake ?               function res = simul_ber(m,n,p)     // creates a vector simulating a binominal law with parameters n p     res = grand(1,m,"bin",n,p)     endfunction          function clt_binom(m,n,p)     clf()     v = ones([1:m])     histplot(floor(sqrt(m)),(cumsum(simul_ber(m,n,p)) - p*n*cumsum(v))./(sqrt(cumsum(v))*sqrt(n*p*(1-p))))     t = [-10:0.05:10]     plot2d(t,(1/sqrt(2*%pi))*exp((-t^2)/2))     endfunction      Here's the type of graph I obtain: ![enter image description here](http://img15.hostingpics.net/pics/465514Capture.png) There seems to be a problem of 'height'. I checked several times but I can't find my mistake.


title: Identifying subset of variables with most predictive power of the other variables

text: I have 31 numeric variables (e.g., A-AF) for which I am attempting to identify the smallest subset of those variables that will predict the values of the remainder of the variables with a CI of 90% or greater. For example, I would like end up with a model in which the values of variables A, B, C, D, and E, together, predict the values of F-AF, where variables A-E are the minimum number of variables needed to predict every other variable in the data. How might I approach this problem? I have the idea to create 31 different models to predict all 31 variables with 90% CI, then use some kind of "common variable" algorithm to identify which variables are common to all the models... But I don't know how to implement that. I'm using RapidMiner 6, but I could learn to use R or any other tool to accomplish this task.


title: Looking for beginner resources with effect sizes

text: First post here; first I would like to say that I have no background in stats whatsoever and not so much in math either (enriched high school math, and that's it). I started an MA and we have to do a meta-analysis; we've also been reading a lot of publications with effect sizes. I checked Wikipedia's page for effect sizes; some parts are helpful, and for others I feel the gap when I'm reading, and it almost looks like gibberish. I've found some lectures on YouTube which I think will be useful. I was wondering if some of you would have suggestions for a good point to start learning about this process? EDIT: I think I need to make my question more precise as I read the comments and when I skim the answers in the other questions, it still seems too high level sometimes. The effects sizes books seemed like it could help though, thanks for that! I do understand mostly what a meta-analysis is (I think!) but it's how people go from the data to the numbers in their effect sizes that I am wondering about. Like how do effect sizes, chi-square work etc. When I read the definitions on these formulas, realize I would also need a course on distribution, regression, multivariate as I have no idea what these mean and it doesn't ring much of a bell. What I think helped the most with the chi- square is when I read something about comparing product reviews on Amazon where one product had 5-stars but one review and another product had 4.5 stars but 350 reviews. If there are resources with worked cases for the type of formulas used in a meta-analysis, to understand what they're for and how they're used and with realistic examples so I can refer to something I know, I think that would greatly help. In fact, I need beginner resources in stats generally so I can build on that.


title: Varimax Rotated factors are very much NOT uncorrelated?

text: On the back of an earlier question i am having an issue. Software = SAS JMP Pro 11 Earlier Question My rotated factors (post Principal Components Analysis [PCA]) are not at all uncorrelated. I have used Variance Inflation Factor (VIF) scores to initially look at this. I expected some increase in the VIFs across the board but what i have is pretty dire with VIFS in the 10000s on occasion. The PCA factors unrotated are completely uncorrelated when looking at VIFs (all 1.5 or less). What problems can occur in the rotation step to cause this? The original variables load onto to many of the unrotated factors to too high a degree to be of much real world use. I have seen written there is no real difference between orthogonal or oblique rotation and I have also found this to be the case The only thing i can think is my non positive definitive matrix but I thought (I do not know) that if the unrotated factors are uncorrelated this meant the matrix was handled OK?


title: How can a probability distribution P not factorize over a graph H when P satisfies the independencies implied by H

text: How can a probability distribution P not factorize over a graph H when P satisfies the all the global independencies implied by H? Here's an example: Let $X_1, \dots X_4$ be 4 random variables that can take on 0 or 1. The graph $H$, is a circle: $X_1 \rightarrow X_2 \rightarrow X_3 \rightarrow X_4 \rightarrow X_1$. You can give a (valid) CPT that specifies the joint distribution of $X_1,X_2,X_3,X_4$ s.t. it satisfies $(X_1\bot X_3) \mid (X_2,X_4)$ and $(X_2 \bot X_4) \mid (X_1,X_3)$. However, you can specify an assignment to $X_1 \dots X_4$ s.t. such that the joint distribution does not equal the product of the clique factors. I'm trying to get some intuition for this, and I don't have any right now. (Note: I've been told that allowing zero probabilities in Markov models plays a role in understanding this.)


title: Bonferroni correction for gene functional category enrichment

text: I have performed a clustering on 20K genes, and I want to understand how coherent resulting clusters are by calculating the significance of the overlap of the clusters with a number of known GO gene categories. I use hypergeometric p-values for this purpose, and I also want to correct for multiple testing using Bonferroni correction. My question is about how to perform Bonferroni correction. I know I have to multiply each hypergeometric p-value with the number of experiments performed for calculating the p-values, but I am a little confused about what "the number of experiments" means. I intuitively think that in my case, the number I have to multiply each hypergeometric p-value is `#clusters * #GO-categories`; is that right? So, for example, if I am checking the enrichment of each of 100 clusters with each of 200 gene categories, I will multiply the p-value for each cluster-category pair by 20000. Could you please confirm that?? Thank you!


title: Poisson regression versus Poisson random effect regression models

text: i have a fitted a classical Poisson regression model to my 584 claim frequency data set for a period of 5 years (panel data)i have found out that predictors area,age and make of vehicle are siginificant (preferred model).Other predictors are age and gender of insured and cc of car. Now i want to fit a Poisson random effect and compare it to the preferred model fitted by classical Poisson regression. What can i include as my random effects?


title: What is a forward/inverse model?

text: The term forward model comes up a lot when reading about Bayesian modelling. I am yet to understand what exactly is the forward model? Is it the model that describes the output/observed variable and the inverse model is that given the observation, what is the underlying model that explains this observation?


title: How to convert daily variances to a monthly volatility and then annualize it?

text: I work out the conditional variance using a GARCH model based on daily returns as follows:                        conditional variance          1/1              0.024879711     2/1              0.02607681     3/1              0.026476405     4/1              0.025159491     5/1              0.023916521     .     .     .     31/1             0.022825424      How can I work out the monthly volatility and then annualize it based on these daily variances? I tried to sum the daily variance from 1/1 to 31/1 and then sqrt(252*sum(1/1 to 31/1)), but then the annualized volatility gets ridiculously large like 150%...


title: Which statistical analysis replaces MANOVA when there are too many dependent variables?

text: I want to understand if an experiment gives equal response at two different levels of an experimental parameter. What I did was to repeat the experiment 6 times at the first level of the parameter and 6 times at the second level. The problem is that the experiment measures 20 variables, so I have 12 experiments and 20 measured variables. I'm not sure, but I think that MANOVA will not work in this setting. Is this correct? What kind of alternative test should I use?


title: Bonferroni for outlier detection?

text: I am reading a book on time series analysis and I am having problems understanding the section about outlier detection. The authors say that when you want to know whether at a certain time $T$ there was an outlier, you should use a certain test statistic and a test with size less than $\alpha$. But when you don't know where an outlier could be and you have a time series of size $n$ then you should use the same test statistic for each point but you should use tests of size $\alpha/n$. They say that this is an application of the conservative Bonferroni correction. I just don't understand this. Doesn't this mean that there will be lots of outliers that you detect in short time series but don't detect in large ones? After all, spam filters don't have stronger spam criteria for people with more incoming email, right?


title: List of topics for a 'quantitative reasoning' course

text: I am planning to develop a course which I tentatively call as 'Quantitative Reasoning'. The goal of the course is to equip a typical undergraduate student with sound quantitative reasoning skills so that they can critically evaluate statistical and related quantitative claims they may encounter as part of their personal/professional lives. In contrast to the traditional 'Intro to Statistics' type of courses the focus of this course will not be on derivation of formulas, use of statistical tests etc but on interpretation and understanding of statistical concepts. Towards that end I have developed a rough list of topics that I would like to cover in this course:   * Logic of hypothesis tests   * Interpretation of p-values and clarify common misconceptions.   * Bayesian statistics and its relationship to frequency statistics (at an informal level)   * Conditions required for causality    * How to assess/test for causality?   * Advantages & disadvantages of experiments, surveys etc In light of the above, I have three questions:   1. What else should you think I should include/exclude from such a course?   2. Are there any textbooks that may be useful given the above goal?   3. Are you aware of any other course that attempts to accomplish the above? Links to syllabus would be very helpful. If it matters, the target student for the above course is an undergraduate student in the US possibly at the freshman or sophomore level.


title: Literature on generating "similar" synthetic time series from observed time series

text: The motivation for this question is from Finance. I have some market data (daily time series) for the price of some securities and I would like to generate synthetic versions of these which are statistically "similar" (in some sense) for testing trading strategies. Is there literature on this subject? I was hoping there would be a way of manipulating the market data that I have in a deterministic way (such as, say, taking the first difference between consecutive values and swapping these around) rather than extracting statistical information about the time series e.g. autocorrelation and then generating new random variables etc. to get a new time series. I am being deliberately vague about what I mean by "similar" as I don't know how realistic my question is and don't want to constrain it further.


title: Rebalancing discrete probability distribution

text: I'm working on a machine learning algorithm and have gotten stuck with how to rebalance a discrete probability distribution. I have a distribution represented as a simple array of $n$ numbers which are all bounded between $0$ to $1$ and that always sum to $1$. At a particular point I wish to update on of the numbers in the distribution by some random amount (it's still between $0$ and $1$) . At that particular point I must also update the remaining $n-1$ numbers so that the sum of the numbers in the distribution remain $1$. My problem is that any algorithm I come up with for rebalancing fails for certain cases. First algorithm was to uniformly distribute the difference to the remaining values. Obviously that didn't work (eg. can't subtract $0.1$ from $0.05$). Second approach was as follows: Let the distribution be represented by $\{p_{1},p_{2},...,p_{n}\}$ and let's assume I wish to update $p_{1}$ for the sake of argument. Let $p^{new}_{1}$ be the new value for $p_{1}$. Then I rebalance the remaining $n-1$ numbers by: $p_{i} = p_{i} + ( p_{1} - p^{new}_{1} ) \frac{p_{i}}{sum}$ Where $sum = \sum_{j=2}^{n}(p_{j})$ Which also fails. For example, if some values in the distribution are $0$ then they will not get increased in the case of a positive update. I tried to correct this by checking if the sum is $0$ but that only allows me to solve the case where all $n-1$ numbers are $0$. I feel like this should probably be very easy and that I'm just having a brain freeze at the moment...


title: PDF Manipulation for Bayesian analysis

text: This post pertains to Bayesian pdf manipulation.   * Firstly, assuming a prior probability specified as _Gamma_ distribution such that $\alpha = \mu_{0}^{2}/\sigma_{0}^{2}$ and $\beta = \mu_{0}/\sigma_{0}^{2}$; is it correct to assume, substituting into the Gamma distribution, that the prior has the following form: $$\pi(\mu) \propto \mu^{(\mu_{0}^{2}/\sigma_{0}^{2})-1}\text{exp}\{-(\mu_0/\sigma_{0}^{2})\mu\}1_{\{\mu>0\}}$$   * Secondly, as an aside, assuming the prior adhered not to a Gamma distribution, but, in fact, to a _Beta_ distribution, would it be correct to assume the following: Given that the pdf for a Beta distribution has the form: $$\frac{\mu^{\alpha-1}(1-\mu)^{\beta-1}}{B(\alpha, \beta)}$$ The prior distribution would have the following form: $$\pi(\mu) \propto \mu^{\alpha-1}(1-\mu)^{\beta-1}$$ And, given that the likelihood remains normally distributed, the posterior would have the form: $$\pi(\mu|\textbf{x}) \propto \mu^{\alpha-1}(1-\mu)^{\beta-1}\text{exp}\{-\frac{1}{2\sigma^{2}}\sum_{i=1}^{n}(x_{i}-\mu)^{2}\}$$ With $\frac{\alpha}{\alpha+\beta}=\mu_{0}$ and $\frac{\alpha\beta}{(\alpha+\beta)^{2}(\alpha+\beta+1)} \space{ } = \sigma_{0}^{2}$ Owing to: $E[\mu] = \frac{\alpha}{\alpha+\beta}$, $Var[\mu] = \frac{\alpha\beta}{(\alpha+\beta)^{2}(\alpha+\beta+1)}$   * Finally, is it correct to suggest that $\pi(\mu|\textbf{x})$ in the above example (Beta/Normal posterior) is a member of the exponential family given that it can be rewritten in the general form: $$h(\mu)^{n}\text{exp}\{\sum t(x_{i})\psi(\mu)\}$$


title: How should I analyze data on birth intervals of women?

text: I have data on birth intervals of women. Birth interval refers to the time interval from one child's birth date until the next child's birth date. Here, dependent variable is the birth interval, defined as the length of time between two successive births for different parities, which are measured in single months. Here four duration of birth are used. A variety of demographic and socio-economic characteristics influence women's spacing practices., therefore the above mentioned dependent variables depend on a number of demographic and socio-economic factors. We deals with many explanatory variables like: 1.survival of the previous child,   2.mother's age at the birth of the index child,   3.mother's education,   4.mothers job,   5.sex of the index child,   6.place of residence,   7.division etc.   Are these data correlated survival data? How can we analyze this? Is this multi-level frailty model reasonable?


title: How do I do multiscale HMM classification?

text: I'm using hidden Markov models to classify some accelerometer data. I take the Fourier transform of the raw data at a given window length, and then train an HMM for each class, and every test instance is classified as the class corresponding with the HMM giving the highest log likelihood. Now, I want to try considering multiple scales/window lengths at once. The approaches I tried so far:   1. Concatenate the observations for different window length in the time dimension.    2. Add the observations for different window length as additional features   3. Train an HMM for each window length and take the product of loglikelihoods across scales The results of the first two were extremely poor (chance level, basically, but I suspect there is a bug) and the third was mediocre (roughly equal to the worst single scale results). Any idea how I might improve this? Could I be doing something completely wrong? Which of the first two approaches is better? And is the product of loglikelihoods a sensible thing to use, or would I be better off with some sort of boosting, for example? Your advice would be much appreciated.


title: Statistical modelling

text: I am novice in Stats. I drew 5 numbers from 5 Odd and 4 Even numbers. After 10,000 runs, what will be statistical model look like; will be sure the ratio for odd to Even will end in ratio 3:2 ? If the same model is made for 5 numbers consisting 3 Odd and 2 Even and 3 numbers are drawn will it also end in 2:1 ratio? What model will look like for 1000 or smaller run? Thanks in advance, Rodney


title: Is R viable for production (deployed) code

text: I have read a number of articles that talk about companies such as Google, Facebook, and many others using R for research. The other scenario I have read about is companies using R to prototype an analytics solution and then re- implement it in another language. I am trying to find literature on companies using R for actual production analytics code. A use case might be a recommender system that user interacts with via a web page that gets a response from an R script executed on a remote server. The fact that I am having trouble finding such reports makes me wonder if it is unadvisable. If so why?


title: Regression Lines With Same Intercept

text: So, I struggle with Regression a lot. I just found out how to get 2 lines with the same slope, but I cannot manage to get 2 lines with the same intercept. I read about ANCOVA a lot (because I thought this was what I needed), but no one uses the same intercepts; just the same slope. Can someone help out with this?


title: How many random permutations to cover all possible permutations?

text: I have code that generates a random permutation. In my case, a permutation consists of N binary features, and each of the N features is set or unset randomly. How many times must I generate a random permutation in order to be reasonably assured that I have covered all possible permutations? I'm not sure how to define "reasonably assured" (50%, 95%?). For example, let's say that N = 10, so there are 1024 possible permutations, how many times should I randomly generate a permutation to be 50% confident that I have generated all 1024 permutations? It seems to me that this is related to the 36.8% duplicate hit rate when drawing a sample with resampling, but I'm not a statistician.


title: Gaussian clusters and original distributions

text: In Gaussian clustering (i.e. General Mixture Models) we model the data with some clusters. For example, in the below figure, we have two clusters $C_1, C_2$, each of which are modeled with a Gaussian (Normal) distribution: $p_{C_1} \sim \mathcal{N}(\mu_1, \Sigma_1)$ and $p_{C_2} \sim \mathcal{N}(\mu_2, \Sigma_2)$ with proportions $\pi_1$ and $\pi_2$. ![enter image description here](http://i.stack.imgur.com/A16KF.png) The probability that a given point $x$ belongs to $C_1$ (is generated from $C_1$) is: $$ p(C_1 | x ) = \frac { \pi_1 \mathcal{N}(\mu_1, \Sigma_1)} {\sum_{i=1}^2{ \pi_i \mathcal{N}(\mu_i, \Sigma_i) } } $$ The probability that $x$ belongs to $C_2$ can be similarly calculated. Now, assume that we have a group of points that are generated (sampled) from another Gaussian distribution $q(x) \sim \mathcal{N}(\mu_q, \Sigma_q)$. ![enter image description here](http://i.stack.imgur.com/4Z7um.png) Obviously, we can compute the probability of each point belonging to each cluster. However, we instead want to compute **the expected probability of a sample belonging to each cluster**. In other words, if we generate infinite samples from the $q(x)$, what is the **average (expected)** value of $p(C_1 | x )$ for all points? **P.S. ** I need the solution in closed-form.


title: How to draw a plot with mean and CI?

text: I want to draw a plot in Matlab or Stata. I have a mean and CI for some distribution. there should be a line with the upper end of the line representing the upper CI and the lower end the lower CI and the middle, the average. I am not sure how to do it or how to search for it. Can anyone kindly help me out ? Thanks. NB: is there a way to connect these CI ? (with straight lines ? 3 straight lines)


title: PDF for a function of random variables

text: If $g=f(x,y)$ is a function of independent random variables $x$ and $y$ then how do we arrive at the expression for the probability density function of $g$, $$f_G(g) = \iint f_X(x)f_Y(y)\delta(g-f(x,y))\mbox{d}x\mbox{d}y\ ?$$ I have been looking at some statistics books but cannot find it. Neither am I able to derive it myself. Any help would be appreciated.


title: Calculating likelihood from RMSE

text: I have a model for predicting a trajectory (x as a function of time) with several parameters. At the moment, I calculate the root mean square error (RMSE) between the predicted trajectory and the experimentally recorded trajectory. Currently, I minimise this difference (the RMSE) using simplex (fminsearch in matlab). While this method works to give good fits, I would like to compare several different models, so I think I need to compute the likelihood so that I can use maximum likelihood estimation rather than minimising the RMSE (and then compare the models using AIC or BIC). Is there any standard way of doing this?


title: Is sample mean the "best" estimation of distribution mean in some sense?

text: By (weak/strong) law of large numbers, given some iid sample points $\{x_i \in \mathbb{R}^n, i=1,\ldots,N\}$ of a distribution, their sample mean $f^*(\{x_i, i=1,\ldots,N\}):=\frac{1}{N} \sum_{i=1}^N x_i $ converges to the distribution mean both in probability and a.s., as sample size $N$ goes to infinity. When the sample size $N$ is fixed, I wonder if the LLN estimator $f^*$ is an estimator best in some sense? For example,   1. its expectation is the distribution mean, so it is an unbiased estimator. Its variance is $\frac{\sigma^2}{N}$ where $\sigma^2$ is the distribution variance. But is it UMVU?   2. is there some function $l_0: \mathbb{R}^n \times \mathbb{R}^n \rightarrow [0,\infty)$ such that $f^*(\{x_i, i=1,\ldots,N\})$ solve the minimization problem: $$ f^*(\{x_i, i=1,\ldots,N\}) = \operatorname{argmin}_{u \in \mathbb{R}^n} \quad \sum_{i=1}^N l_0(x_i, u)? $$ In other words, $f^*$ is the best wrt some contrast function $l_0$ in the minimum contrast framework (c.f. Section 2.1 "Basic Heuristics of Estimation" in "Mathematical statistics: basic ideas and selected topics, Volume 1" by Bickle and Doksum). For example, if the distribution is known/restricted to be from the family of Gaussian distributions, then sample mean will be the MLE estimator of distribution mean, and MLE belongs to the minimum contrast framework, and its contrast function $l_0$ is minus the log likelihood function.   3. is there some function $l: \mathbb{R}^n \times F \rightarrow [0,\infty)$ such that $f^*$ solve the minimization problem: $$ f^* = \operatorname{argmin}_{f} \quad \operatorname{E}_{\text{iid }\{x_i, i=1,\ldots,N\} \text{ each with distribution }P } \quad l(f(\{x_i, i=1,\ldots,N\}), P)? $$ for any distribution $P$ of $x_i$ within some family $F$ of distributions? In other words, $f^*$ is the best wrt some lost function $l$ and some family $F$ of distributions in the decision theoretic framework (c.f. Section 1.3 "The Decision Theoretic Framework" in "Mathematical statistics: basic ideas and selected topics, Volume 1" by Bickle and Doksum). Note that the above are three different interpretations for a "best" estimation that I have known so far. If you know about other possible interpretations that may apply to the LLN estimator, please don't hesitate to mention that as well. Thanks and regards!


title: What does truncated distribution mean?

text: In a research article about sensitivity analysis of an ordinary differential equation model of a dynamic system, the author provided the distribution of a model parameter as Normal distribution (mean=1e-4, std=3e-5) truncated to the range [0.5e-4 1.5e-4]. He then uses samples from this truncated distribution for simulations of the model. What does it mean to have a truncated distribution and sample from this truncated distribution? I could come up with two ways of doing this:   * Sample from a Normal distribution but ignore all the random values falling outside the specified range before simulations.   * Somehow obtain a special "Truncated Normal" distribution and get samples from it. Are these valid and equivalent approaches? I believe in the first case, if one were to plot the experimental cdf/pdf of the sample, it would not look like a Normal distribution because the curves do not extend to $\pm\infty$.


title: R: process data frame subgroups and merge results together

text: I need to process three data frames containing the same subgroups (name). So far I've been using this approach:               results = data.frame(name = factor("dummy"), col1 = 1, col = 2)      for( name in df1$name ) {       new.results = process(name, df1[df1$name == name, ], df2[df2$name == name,], df3[df3$name == name, ]       results = rbind(results, new.results)     }          return(results)      The problem with this code is that the process() function must return the same layout than the 'results' data frame. If I change the contents returned by process() I also have to change the outer function. Also the first row in the results data frame must be removed. Is there an easier way to do this?


In linear regression, is the $R^2$ value enough to assess whether the relationship between the independent and dependent variable is linear? It gives the amount of variability in the dependent variable explained by the independent variable. I know that you can plot residuals versus the x value or residuals versus the y value and see if there is a pattern (if there is a pattern then the relationship is not linear). But doesn't the correlation coefficient give enough information about linearity?


title: What does interaction depth mean in GBM?

text: I had a question on the interaction depth parameter in gbm in R. This may be a noob question, for which I apologize, but how does the parameter, which I believe denotes the number of terminal nodes in a tree, basically indicate X-way interaction among the predictors? Just trying to understand how that works. Additionally, I get pretty different models if I have a dataset with say two different factor variables versus the same dataset except those two factor variables are combined into a single factor (e.g. X levels in factor 1, Y levels in factor 2, combined variable has X * Y factors). The latter is significantly more predictive than the former. I had thought increasing interaction depth would pick this relationship up.


title: What is the maximum entropy probability density function for a positive continuous variable of given mean and standard deviation?

text: What is the maximum entropy distribution for a positive continuous variable, given its first and second moments? For example, a Gaussian distribution is the maximum entropy distribution for an unbounded variable, given its mean and standard deviation, and a Gamma distribution is the maximum entropy distribution for a positive variable, given its mean value and the mean value of its logarithm.


title: How to use gradient to obtain critical points of Time series?

text: How to determine the number of critical points for some time series (for example, using gradient)? As far as I understand I have to do the following steps:   1. Fit the time series with some _"curve"_ function (polynomial, sum of sin, etc.) _f(x)_.   2. Use this _f(x)_ as data input for gradient method to obtain critical points.


title: Software to plot a log graph

text: I have to plot a graph with log distribution on y axis. The values are: 10^-3..10^3. What software do you suggest me to use. My OS is Ubuntu, so I prefer software for Linux. Thanks.


title: How to model the sum of Bernoulli random variables for dependent data?

text: I have almost the same questions like this: How can I efficiently model the sum of Bernoulli random variables? But the setting is quite different:   1. $S=\sum_{i=1,N}{X_i}$, $P(X_{i}=1)=p_i$, $N$~20, $p_i$~0.1   2. We have the data for the outcomes of Bernoulli random variables: $X_{i,j}$ , $S_j=\sum_{i=1,N}{X_{i,j}}$   3. If we estimate the $p_i$ with maximum likelihood estimation (and get $\hat p^{MLE}_i$), it turns out that $\hat P\{S=3\} (\hat p^{MLE}_i)$ is much larger then expected by the other criteria: $\hat P\{S=3\} (\hat p^{MLE}_i) - \hat P^{expected} \{S=3\}\approx 0.05$    4. So, $X_{i}$ and $X_{j}$ $(j>k)$ cannot be treated as independent (they have small dependence).   5. There are some constrains like these: $p_{i+1} \ge p_i$ and $\sum_{s \le 2}\hat P\{S=s\}=A$ (known), which should help with the estimation of $P\{S\}$.  How could we try to model the sum of Bernoulli random variables in this case? What literature could be useful to solve the task? **UPDATED** There are some further ideas: (1) It's possible to assume that the unknown dependence between ${X_i}$ begins after 1 or more successes in series. So when $\sum_{i=1,K}{X_i} > 0$, $p_{K+1} \to p'_{K+1}$ and $p'_{K+1} < p_{K+1}$. (2) In order to use MLE we need the least questionable model. Here is an variant: $P\{X_1,...,X_k\}= (1-p_1) ... (1-p_k)$ if $\sum_{i=1,k}{X_i} = 0$ for any k $P\{X_1,...,X_k,X_{k+1},...,X_N\}= (1-p_1) ... p_k P'\{X_{k+1},...,X_N\}$ if $\sum_{i=1,k-1}{X_i} = 0$ and $X_k = 1$, and $P'\{X_{k+1}=1,X_{k+2}=1,...,X_N=1\} \le p_{k+1} p_{k+2} ... p_N$ for any k. (3) Since we interested only in $P\{S\}$ we can set $P'\{X_{k+1},...,X_N\} \approx P''\{\sum_{i=1,k}{X_i}=s' ; N-(k+1)+1=l\}$ (the probability of $\sum_{i=k+1,N}{X_i}$ successes for N-(k+1)+1 summands from the tail). And use parametrization $P''\{\sum_{i=k,N}{X_i}=s' ; N-k+1=l\}= p_{s',l}$ (4) Use MLE for model based on parameters $p_1,...,p_N$ and $p_{0,1}, p_{1,1}; p_{0,2}, p_{1,2}, p_{2,2};...$ with $p_{s',l}=0$ for $s' \ge 6$ (and any $l$) and some other native constrains. Is everything ok with this plan? **UPDATED 2** Some examples of empirical distribution $P\{S\}$ (red) compared with Poisson distribution (blue) (the poisson means are 2.22 and 2.45, sample sizes are 332 and 259): ![sample1](http://i.stack.imgur.com/y441B.png) ![sample2](http://i.stack.imgur.com/qAy56.png) For samples (A1, A2) with the poisson means 2.28 and 2.51 (sample sizes are 303 and 249): ![sample3](http://i.stack.imgur.com/jCNOE.png) ![sample4](http://i.stack.imgur.com/QWHiS.png) For joined samlpe A1 + A2 (the sample size is 552): ![sample 3 + sample 4](http://i.stack.imgur.com/VxVV6.png) Looks like some correction to Poisson should be the best model :).


title: Weighted community analysis in R

text: I am trying to do a weighted community analysis. I have association matrices calculated (in SOCPROG) year by year over a 30-year period for a population of 80-100 individuals. I would like to explore how community structure changes over time using a weighted analysis. I have been reading up on the best way to do this and I have found the fastgreedy.community algorithm in R. However, I am completely unfamiliar with R and am struggling to understand it. I was wondering if anyone could advice me on the best way to do this sort of analysis, or point me in the direction of any good tutorials please.


title: How to get a weighted-estimate for mean difference?

text: I wanted to estimate the difference of two means (say: $\mu_1$ and $\mu_2$). I have the estimate of means (say: $m_1$ and $m_2$) with variance estimates (say: $v_1$ and $v_2$). How do I calculate the weighted-estimate for difference of means ($\mu_2$ - $\mu_1$)? I am thinking these two should work (in the order of preference):   1. $\frac{v_1}{v_2} m_2 - m_1$   2. $\frac {v_1} {v_2 + v_1} m_2 - \frac {v_2} {v_2 + v_1} m_1$ This problem may seem like the Behrens-Fisher problem. However, in my case the data is NOT normally distributed. Any help is appreciated.


title: How to optimise multi-layer neural network architecture using the genetic algorithm in MATLAB

text: Can someone please provide me with a summary of how to optimise multi-layer feedforward neural network architecture using the genetic algorithm? I've read a few papers that say this is an effective method to select architecture, but i'm struggling to do it myself. By architecture I mean the number of layers, and the number of neurons in each layer. **My attempt so far** My neural network has 24 inputs, and 3 outputs, and I have 2000 examples that I will use to train the network. My network must have a minimum of 1 hidden layer, and a minimum of 1 neuron in each layer. It also must have no more than 3 hidden layers, and no more than 30 neurons in each layer. Therefore, the cost function needs to take in a 1 by 3 vector input, with each element in the vector being an integer between 1 and 30. My cost function takes a vector N with size 1 by 3, and then deletes an element of this matrix if it contains a zero. This is so that the neural network can select the number of layers. My fitness function is shown below.               Function [C] = fitness(N) %define the input and output arguments of the function          X=evalin('base','inputTrainingData'); %the data used to train the neural network     T=evalin('base','outputTargetData');          Xtest=evalin('base','inputTestingData'); %the data used to test the neural network     Ttest=evalin('base','outputTestingData');          N(N == 0) = [] %delete all zeros from the vector          net=feedforwardnet(N) %Create the network with architecture defined by N          net=train(net,X,T)          C = MSE(net(Xtest),Ttest)  %Average distance between net(Xtest) and Ttest      I then use the genetic Algorithm in the Matlab Global Optimisation Toolbox to minimise this fitness function. I have constrained the solutions so that each gene is an integer between 0 and 30. This gives pretty poor and inconsistent results, the average fitness of the population increases for the first 3 or 4 generations, and then ceases to increase for the next 100. I was hoping for better results. I once tried rounding $N$ at the start of the cost function, instead of constraining the solutions to be integers, and then just rounding the final optimal solution. This seems to always give better results but I do not know why. Is there anything I'm doing wrong, or anything I can change to get better results? Thanks!


title: Classification/Regression Tree with nonngeative response

text: I try to model the duration until a unit is inspected by a large number of possible explanatory variables. The duration is non-negative and the explanatory variables are factors and numerical variables. If I knew the best explanatory variables then I would use `glm` in R. I find family "Gamma" with link "log" appropriate for modelling something like a waiting time. Can I combine trees and `glm` in R? E.g. with `randomForest` I can use a formula but the regression is a simple linear model (like `lm`). Thanks!


title: Are there any drawbacks to using the normal equation in AR parameter estimation?

text: Is it ok to try and estimate the parameters $c$ and $p$ of a standard AR(1) model of the form $$ X_t = c+ pX_{t-1} +e_t $$ solving the normal equation $$ \Theta = (X^T X)^{-1}X^TY $$ where $\Theta$ is the parameter vector, $X$ a vector containing the $X_{t-1}$ values and $Y$ the $X_t$ vector? The $Y$ vector would obviously contain $t-1$ values (the first value removed) so that it would equal $X$ in length. We'd use something like               X          Y     NA         1     1          2     2          3     3          4      and of course ignore the first row. In other words, is it wrong to apply this "linear regression with multiple variables" method to a time series method of the AR (perhaps others, ARMA, GARCH, etc) sort? Since I haven't really seen the normal equation used for this sort of problem I reckon the answer to my question is Yes. A secondary question would then be "Why?" Why isn't this method OK?


title: Linear regression and non-invertibility

text: In linear regression there are two approaches for minimizing the cost function: The first one is using gradient descent. The second one is setting the derivative of the cost function to zero and solving the resulting equation. When the equation is solved, the parameter values which minimizes the cost function is given by the following well-known formula: $$ \beta = (X^TX)^{-1}X^TY $$ where $\beta$ is the parameter values, $X$ is the design matrix, and $Y$ is the response vector. Note that to have a solution $X^TX$ must be invertible. I think that even if $X^TX$ is non-invertible we can still minimize the cost function using the first approach (gradient descent). If this is true, what bothers me is which property of gradient descent makes it not vulnerable to this kind of problem. thanks


title: Poisson regression in R glm(): concept and practice

text: I have a dataset `data` about a group of persons developing events with Poisson distribution. These events are uniform and can be recurrent in some of the subjects. Their observation time is different. I tried to use R to calculate the incidence rate (in 1000 person years) of the event in the whole group and specific subgroups (e.g. male and female). `no.of.event` is the total number of observed events and `time.in.years` is the total observation time in years. As I do not understand Poisson regression well, I would like to ask:   1. If `no.of.event` equal to the total cumulative number of all events observed?   2. And `time.in.years` be the time from the start of observation till study end/ default or death? Or be the time from the start of observation till the development of event of interest?   3. In those persons with recurrent events what would be `time.in.years`? Say if someone has developed 3 events and did not default, is the time counted from the start of observation till the development of third events? Or from the start till the end of observation?   4. Is the following R script right?              model <- glm(no.of.event ~1, offset(time.in.years), family="poisson", data=data)        5. Should incidence (cases per 1000 person years) be calculated by this?              exp(model$coef[1])*1000        6. And the 95% confidence interval be like this?              exp(model$coef[1] + 1.96*sqrt(diag(vcov(model)))[1])     exp(model$coef[1] - 1.96*sqrt(diag(vcov(model)))[1])        7. If I want to calculate specific incidence, say among males and females and I have a column in the data frame called `sex` with 1 coded for males, should I write something like this?              model <- glm(no.of.event~1, offset(time.in.years), family="poisson",                   data=data[data$sex==1, ])


title: Kruskal Wallis test with high type 1 error

text: using this code               pvalue = rep(NA, 100000)     n = c(352, 198, 170)      group = rep(1:3, n)     for (i in 1:100000){     x=c(rnbinom(352, size=0.9563, mu=2.27), rnbinom(198, size=1.0468, mu=2.27), rnbinom(170, size=1.3264, mu=2.27))      kruskal = kruskal.test(x ~ group)     pvalue[i] = kruskal$p.value}     result = length(pvalue[pvalue<=0.05])/100000         result          I tried to examine the risk of type 1 error when comparing three neg. binomial distributions with same mu and different variances using the Kruskal Wallis Test. To my astonishment, I found out, that 8.3% of the 100,000 conducted tests had a p-value lower than 5%. I still remember from my statistics lessons, that non-parametric tests like the Kruskal Wallis are said to be more conservative and, if I remember correctly, less likely to lead to a type 1 error compared to parametric tests, when applied to data which violate certain assumptions like homogeneous variances or normal distribution. However, in this example, Kruskal Wallis apparently led to a higher type 1 error. I would be thankful, if someone could confirm to me, if the Kruskal Wallis Test could possibly be that anti-conservative or if my code is maybe wrong? However, the result seems logical, as different variances in the samples could lead to the conclusion, that they stem from different populations... BTW, I know there are better methods to examine the distributions at hand, I just tried out the Kruskal Wallis Test "for fun". Thank you very much in advance! added:   Doing the same with ANOVA instead of Kruskal Wallis, I get a lower result than 5%:               > pvalue = rep(NA, 100000)     > n = c(352, 198, 170)         > group = rep(1:3, n)       > for (i in 1:100000){       +   x=c(rnbinom(352, size=0.9563, mu=2.27), rnbinom(198, size=1.0468, mu=2.27), rnbinom(170, size=1.3264, mu=2.27))       +   result = aov(x ~ factor(group))       +   pvalue[i] = anova(result)[1,5]}       > robust = length(pvalue[pvalue<=0.05])/100000       > robust      [1] 0.0438          Probably there is something wrong with the ANOVA-code, but I can't figure out what. I would appreciate any help!


title: Plotting high dimensional data

text: I want to plot high dimensional data on x y plane. For that I know three methods: Principal component analysis (PCA), multidimensional scaling (MDS) and a method from spectral graph theory (using the second and third eigenvectors of the laplacian of the graph or w/e). What are the different aspects of each technique and which one would be better if I want to minimize distance between the original vectors and distance between their 2d representations? Are there other (better) techniques? (Preferably something with ready-made python code)


title: Under what conditions do Bayesian and frequentist point estimators coincide?

text: With a flat prior, the ML (frequentist -- maximum likelihood) and the MAP (Bayesian -- maximum a posteriori) estimators coincide. More generally, however, I'm talking about point estimators derived as the optimisers of some loss function. I.e. $$ \hat x(\,. ) = \text{argmin} \; \mathbb{E} \left( L(X-\hat x(y)) \; | \; y \right) \qquad \; \,\text{ (Bayesian) }$$ $$ \hat x(\,. ) = \text{argmin} \; \mathbb{E} \left( L(x-\hat x(Y)) \; | \; x \right) \qquad \text{(Frequentist)}$$ where $\mathbb{E}$ is the expectation operator, $L$ is the loss function (minimised at zero), $\hat x(y) $ is the estimator, given the data $y$, of the parameter $x$, and random variables are denoted with uppercase letters. Does anybody know any conditions on $L$, the pdf of $x$ and $y$, imposed linearity and/or unbiasedness, where the estimators will coincide?


title: What programming language do you recommend to prototype a machine learning problem?

text: Currently working in Octave, but due to the poor documentation progress is very slow. What language is easy to learn and use, and well documented to solve machine learning problems? I am looking to prototype on a small dataset (thousands of examples), so speed is not important. EDIT: I am developing a recommendation engine. So, I am interested in using Regularized Linear Regression, Neural Nets, SVN or Collaborative Filtering.


title: Can R geeglm handle proportion data?

text: I want to fit a regression model to see whether these is changes in the proportion of First-year students over years. I have count data for the total count of First-year students (FirstTimeStudents) and total count of all students (TotalStudents). I want to fit a GEE model with this data and wrote the following code. My output for coefficients are all NA               combFirstYear = geeglm(cbind(data$FirstTimeStudents, data$TotalStudents -     data$FirstTimeStudents) ~ Year, family=binomial(logit), data=data, id=DeptID,     corstr="independence")          summary(combFirstYear)      The reason I fit the model this way is because I read on this website about doing similar thing using glm and it works. But for some reasons, it does not work with geeglm. If anyone can think of different way of answering this question in terms of proportion, I'd really appreciate it too. It says: "R can handle this using glm with the binomial(link="logit") family, with a dependent variable that is actually a two-vector object, the first being the number of 'successes' and the second the number of 'failures'." > Here's an example by modeling hemlock cover with respect to total cover. The > example is a bit forced (as these are not actual samples of individuals) but > you get the idea. We start by generating the summed cover of all species for > each plot: > > `sumcover = tapply(dat`$`cover,dat`$plotID,sum)` > > `coverdat = data.frame(names(sumcover),sumcover)` > > Now we merge these data with our dat5 hemlock object, based on plotID: > > `dat6 = merge(dat5,coverdat,all.x=T,by.x=1,by.y=1)` > > And create our response variables of hemlock cover (successes) and total > cover minus hemlock cover (failures): > > `cover.y = cbind(dat6$cover,dat6$sumcover-dat6$cover)` > > We can now model these data as a binomial process of varying number of > trials per observation: `glm6 = > glm(cover.y~disturb*elev,data=dat6,family=binomial)` > > `summary(glm6)`


title: Correlation or independence on contingency table for large N

text: I have a dataset with about 35,000 individuals described by around 15 categorical variables. I'm trying to study the independence / correlation between these 15 categorical variables. My first idea was to, for each pair of variables, create a contingency table and calculate the $\chi^2$. Then, study the overall difference in the statistic. However, because the population is so large, $\chi^2$ is always significant. I'm having difficulty interpreting and comparing the results for each pair of variables. So, I can summarize my question as follows:   1. For large datasets, when I know $\chi^2$ will almost always be significant, is there an alternative test that will give more reasonable results? I have two ideas, as well   1. I was thinking of taking many bootstrap samples of say 1K individuals. On each sample calculate the correlation, then average over all the bootstrap samples. The average should be a good representation of the overall sample, but I _feel_ like I'm somehow cheating.   2. Can I simply compare the _magnitudes_ of the $\chi^2$ test between the different pairs of variables? The degrees of freedom are different (the categories are of different sizes), which leads me to think this won't make sense.


title: Fit power law for distributions with zeroes

text: I am pretty new to statistics and have some data that I think may follow a power-law distribution. However, it includes zeroes. I understand that mathematically zeroes can't work, but conceptually, would the point of a power law be violated if there are some zeroes in the data? If I were to scale the data by adding 1, would that be fine, or would that throw off everything? I am sure that this is a very newbie question, but that's what I am, and I would love some help in understanding how power laws work. If the distribution can't be power-law, how can I plot it (in R) to best show its characteristics? Would a log-log plot work? Thank you in advance.


title: Risk bounds for Gaussian regression

text: Does anyone know references for non-asymptotic risk bounds for Gaussian regression. Specifically, I am interested in Bayesian regression with a Gaussian prior on the estimator space, and a Gaussian likelihood. I would like to know any results that bound the difference between the risk of the Bayesian estimator and the risk of the "best" estimator, from among the space of estimators on which the prior is defined. I am new to the field, and the only results I have been able to find by googling are regret bounds for online Bayesian regression by Kakade and Ng. Is this a standard result, if so, can anyone recommend a textbook? Thanks much!


title: Fitting special variance structure in mixed model

text: I'm interested in fitting a linear mixed model with this special variance structure on the random effects $\mathbf{u}$: $\begin{eqnarray*} \mathbb{V}\left(\mathbf{u}\right) & = & \mathbf{A}\mathbf{G}\mathbf{A}^{\prime} \end{eqnarray*}$ This variance structure is very similar to Cholesky and Antedependence structures except few differences:   1. $\mathbf{A}$ is not a unit lower triangle or unit upper triangle matrix   2. $\mathbf{G}$ is not necessarily a diagonal matrix I can specify initial values of $\mathbf{A}$. I'd highly appreciate if someone give me some hints to fit this variance covariance structure in `R`. Thanks in advance for your help and time.


title: Moving from RCB to CR design, how much bias is introduced?

text: We designed a RCB experiment and assigned the factor levels to the experimental units randomly inside each block. Let's pretend we changed our mind and we would like to go for a Completely Randomized design. A new completely randomized assignment of factor levels to experimental units is not possible since the experiment has already begun. We have to stay with the RCB assigment (thus having 1 replicate in each former block).   The question is: how much bias would we introduce in the analysis of variance if we analysed the RCB configuration pretending it to be a CR configuration? I know that the RCB configuration can be considered as one of the $n$ equally probable random configurations, but is there a way to account for this bias in the further analysis? **Update**   The reason we would like to move from RCB to CR is that, from modelling studies, our block has no significant effect on the response varible, and, as gung pointed out, it would decrease the power of the test we will perform. Since Blocks do not significantly affect the response variable of the experimental units, in case we moved to a CR design, can the replicates (previously randomized within non-significant blocks) be considered as randomized within the whole population?   This question is relevant also for _post-hoc_ analysis, if a RCB ANOVA results in blocks not being significant, is a CR ANOVA appropriate on a configuration that was not completely randomized? Thank you.


title: Mixed model with repeated measures and both of two treatments with lme or lmer in R

text: I know just enough about `lme()` and `lmer()` to get myself by with simple models, and to get myself into trouble with more complex ones. I'm really confused about the proper way to model some data I'm currently working with, and would be grateful for some help. **The question** : do people give different estimates of a distance after being fatigued by exercise, and does ingesting sugar during the exercise change the relationship? **The setup** : participants are shown a sequence of 3 targets at different distances, and asked to provide a distance estimate for each. They then exercise for a while, and then make 3 more distance estimates to 3 new target distances. This whole procedure happens twice for each participant, on separate days, and one day they ingest sugar, while on the control day they do not. **The data** : DV: Distance estimates IVs:   * Actual distance: 6 possible distances   * Condition: 2 level factor: treatment/placebo   * Measurement occasion: 2 level factor: Pre/Post exercise So, for each subject, I have:   * 3 distance estimates before exercising w/ sugar   * 3 distance estimates after exercising with sugar   * 3 distance estimates before exercising w/out   * 3 distance estimates after exercising w/out. Ultimately, I'd like to know if the estimates people make post-exercise are different depending on the Condition, but those estimates should also be dependent on the starting point of Pre-exercise estimates. Since the actual target distances before and after are different, I can't do a direct comparison (the means of each group of 3 are the same, however). The thing that is tripping me up is that both factors are within subjects. I have repeated measures at different levels (repeated estimates in each state, repeated estimates separated by time, and repeated estimates separated by treatment) and I'm not sure what to treat as fixed effects and what to treat as random effects, whether I have a nested structure or not, and whether I should be looking at random intercepts only or also random slopes (and for which variables). Any help or suggestions are greatly appreciated. Thanks!


title: Variance estimation in Random Effects model

text: I'm studying panel data models in my introductory econometrics class, especially random effects models. Consider the model: $$y_{it}=x_{it}'\beta +c_i+u_{it}$$ with the assumptions $E[c_i]=0=E[u_{it}]$, $E[c_i^2]=\sigma^2_c$, $E[u_{it}]=\sigma_u^2$, $\text{Corr}(x_{it},c_i)=0$ and $\text{Corr}(x_{is},u_{it})=0$. Denote composite error term $c_i+u_{it}=v_{it}$ and $\Omega=E[v_iv_i']$ We've learned that we can estimate $\beta$ consistently using feasible generalised least squares where we estimate $\Omega$ from $\Omega=\sigma_u^2 I+\sigma_c^2 \iota \iota'$ where $I$ is identity and $\iota$ is appropriate sized vector of ones. Using fixed effect residuals $\hat{e}$, with number of cross section observations $N$ and time observations $T$, we get an estimate: $$\hat{\sigma}_u^2=\frac{1}{N\cdot(T-1)}\sum_{i=1}^N \sum_{t=1}^T \hat{e}_{it}^2$$ Furthermore, estimating the between model and estimating error variance, we get and estimate of variance of $\eta=c_i+\bar{u}_i$. Now we can estimate $\hat{\sigma}_c^2$ from: $$\hat{\sigma}_c^2=\hat{\sigma}_{\eta}^2-\frac{1}{T}\hat{\sigma}_u^2$$ (I hope that was enough details, if not please let me know). Now my question is, what happens if this is negative? Is it an indication of some error, or just something that can happen in finite samples? What should I do about it? Last, what is this estimation called, and are there other ways of estimating that variance? I hope my question is well formulated and on the right Stack Exchange—I am by no means proficient in statistics/econometrics.


title: Simulation vs analytic methods to forecast time series process

text: Just asking if someone knows why the prediction intervals are quite different when one uses a time series analytic method of estimation _versus_ when one simulates such time series. For example, I used the forecast package's `auto.arima` function to get the best fit to my data, say it was an ARIMA(1,1,1), and then, on the one hand, I simulated such process doing around 10 thousand simulations and then calculating the 95% percentile with "quantile" function, and on the other hand, I used R's `forecast` package to do it. So I realized that these different approaches gave prediction intervals with different width (actually, those related with simulation approach are closer than those obtained with forecast package). The way I simulated such time series process is simulating the parameters as random variables distributed normally with mean equal to its estimated value and standard deviation equal to its related standard error. The "white noise" variables related with the Moving Average (MA) part of the process were simulated as normally distributed with mean zero and variance equal to the variance of the residuals. Thanks in advance for your help.


title: How to choose train/test sample ratio, for machine learning?

text: I am building a real time machine learning module, which is not based on a huge** sample size, with hyper parameter grid search and cross validation process. I am looking for any insight/advice, as I`m considering one of these options:   1. Use cross validation grid search to look for the best hyper-parameter (HP) fit, and once I found the best HP combination, use it to re-train my classifier on the whole sample set.   2. Break in advance my training set to 2 sub-set, use cross validation by iteratively re-break my training set into train/test sets which searching for the best HP combination, and use the trained classifier with the best HP without going through the last re-train stage described in 1.    3. Do the same as 1, but keep the same random seed as I re-train the classifier. The trade off, as I see it, is between getting a bit more sample size for training, while at the same time loosing assurance that the performance will be good and I`m not overfilling. Again any though/insight into my dilemma are welcome. *Note that I`m using random-forest and extra trees classifiers during grid search. **My sample size is typically between a few hundreds and a few thousands for each class, and the number of features is between 70 and 1500, typically.


title: Why the CDF for the Normal Distribution can not be expressed as a closed form function?

text: I am working my way through Think Stats, where the author states that > "there is no closed form expression for the normal cumulative density > function" but does not provide any further details as to why this is the case, simply saying that the alternative is to write it in terms of the _error function_. Is there some way to intuit why the Normal Distribution can not be expressed as a closed form function?


title: What is the difference between a "nested" and a "non-nested" model?

text: In the literature on hierarchical/multilevel models I have often read about "nested models" and "non-nested models", but what does this mean? Could anyone maybe give me some examples or tell me about the mathematical implications of this phrasing?