title: Unable to understand derivation of Expectation Maximizaton

text: In Paper,               System Identification using Symbolic Chaotic     Sequence, Authored by A. Kurian and H. Leung      download link under section II B, can somebody please explain how to go from Step (8) to the equations (10)? The paper explains how to estimate states and parameters using EM-UKS (Expectation Maximization with Unscented Kalman Smoother) for Moving Average (MA) model driven by deterministic input. The Authors assume that               "dynamics of underlying chaotic system that     generates the input sequence"      is known. So, the set of unknown parameters are the model coefficients $h$, the initial condition of the deterministic signal $c_0$ and the variance of the noise $\sigma_w^2$. Assuming the model in Eq(1): $y(n) = h^T z(n) + w(n)$ to be corrupted with Gaussian White noise $w(n)$, the MA model regresses over the output of a deterministic system $z_n$. So, the input to the MA model is $z(n)$ which is obtained from a chaotic deterministic system. Now, $z(n)$ is converted into binary sequence $s_n$ having values only in [-1,1] by Eq(3). The Authors are estimating the $\theta$ which consists of : states of $z(n)$, the unknown parameters/coefficients $h$ of the MA model, initial condition $c_0^s$ used to generate $z(n)$ and the variance of the measurement noise. They say that since $z(n)$ is obtained from a deterministic system, estimating the initial condition would regenerate the entire time series. So, the unknown parameter set $\theta = [h, c_0^s,\sigma^2_w]$ Q1: Why is the assumption that the dynamics should be known required? Does knowing the dynamics mean that the functional form of the deterministic model in Eq(11) is known and why this assumption is required? Q2: My problem is that I am unable to understand how the estimates are obtained after the Expectation Step in the Maximization Step. The results of the Maximization Step are substituted in the State Equation of the Unscented Kalman Filter. Any help will be very useful.


title: Finding the value of the pressure to maximize mean

text: The peel strength ($Y$) of photo resist in the lamination process of PCBs depends on the pressure ($x$) and temperature ($z$) by the relation: $$Y = 50 + 5x + 10z + (20xz)-(.01(x^2)\cdot (z^2))$$ The temperature is a noise variable and has some distribution with mean 100 and variance 25. Find the value of the pressure to maximize the mean peel strength (Hint: maximize E(Y), where the expectation is taken with respect to $z$) So basic info:   The mean of $z$ is 100 and the var of $z$ is 25. But where do I go from here?


title: How to model a multiplicative effect from this hypothesis?

text: Suppose we have outcome variable Y and two predictors X1 and X2, where X1 and X2 sort of come from a type of category. For example, X1 is the relationship score with friends and X2 relationship score with family. Not the same, but they reflect some thing together (social relationship, in this case). I want to set up a model to test the hypothesis that the diversity/balance of predictors can influence the outcome. Some thing like this: Suppose there are two persons who both have scores of 10 on X1 and X2 total. I suspect that the one with score of 5 on X1 and 5 on X2 will have better outcome than the one with 10 on X1 and 0 on X2. I think a normal model (Y = X1 + X2) or an interaction term (Y = X1 + X2 + X1*X2) seem not to do the job. Maybe there needs to create new variable(s) but I'm not sure which to create.


title: How do I relate the variance of a Normally distributed variable to the variance of the exponentiated transform

text: If I have a Normally distributed variable with a Mean of zero and a given variance, how will this variance relate to the variance of the exponentiated form?


title: 10-fold CV performance test: need to use same folds for both evaluated algorithms?

text: I want to compare the performance of 2 algorithms with 10-fold cross validation. What is the correct procedure and why?   * Procedure A:      1. generate the 10 folds     2. test algo1 with these 10 folds and register results     3. generate 10 new folds     4. test algo2 with these 10 new folds and register results     5. compare performance   * Procedure B:      1. generate 10 folds     2. test algo1 on these 10 folds and register results     3. test algo2 on these 10 folds and register results     4. compare performance Does the fact that in procedure A both algorithms are tested with different folds have an impact on our ability to compare results?


title: Simulation of covariates in regression

text: I am doing a simulation study on a simple linear regression model to see how well OLS estimator performs. I am planning to use coverage probability to assess the estimation of slope parameter using 1000 replications with sample size of 100 Suppose the model is simple as $y_i=0.5+2*x_i+\epsilon_i$ I am not sure which of the following two options to choose in simulation.   1. simulate all 100 $x_i$ values, say from a normal distribution. Then in each of the 1000 replications, use the same 100 $x_i$ values for generation of $y_i$.   2. In each of 1000 replications, first simulate 100 $x_i$ values, and based on these covariate values, generate $y_i$. Any help is appreciated.


title: (Quantile regression) Which standard error for heteroscedasticity & serial correlation

text: I have heteroscedastic and autocorrelated residuals in my multivariate quantile regression model. What's the quantile regression standard error estimator that's robust to this? Something hopefully like HAC Newey West but for quantile regression, or perhaps a bootstrap.


title: How does a Poisson distribution work when modeling continuous data and does it result in information loss?

text: A co-worker is analyzing some biological data for her dissertation with some nasty Heteroscedasticity (figure below). She's analyzing it with a mixed model but is still having trouble with the residuals. Log-transforming the response variables cleans things up and based on feedback to this question this seems to be an appropriate approach. Originally, however, we had thought there were issues in using transformed variables with mixed models. It turns out that we had been misinterpreting a statement in Littell & Milliken's (2006) _SAS for Mixed Models_ that was was pointing out why it is inappropriate to transform _count_ data and then analyze it with a normal _linear_ mixed model (full quote is below). An approach that also improved the residuals was to use a generalized linear model with a Poisson distribution. I've read that the Poisson distribution can be used for modeling continuous data (eg, as discussed in this post), and stats packages allow it, but I don't understand what is going when the model is fit. For the purpose of understanding how the underlying calculations are being made, my questions are: When you fit a Poisson distribution to continuous data, **1)** does it the data get rounded to the nearest integer **2)** does this result in the loss of information and **3)** When, if ever, is it appropriate to use a Poisson model for continuous data? Littel & Milliken 2006, pg 529 "transforming the [count] data may be counterproductive. For example, a transformation can distort the distribution of the random model effects or the linearity of the model. More importantly, transforming the data still leaves open the possibility of negative predicted counts. Consequently, inference from a mixed model using transformed data is highly suspect." ![enter image description here](http://i.stack.imgur.com/iHBoe.png)


title: Generalization of the Signal-Noise ratio for non-Gaussian processes

text: The signal to noise ratio is simple, and is usually defined in the context of simple Gaussian local-level models. In the cause of non-gaussian signal or noise models, do people do things more complicated then the ratio of the variances of the two distributions?


title: Java library for Bayesian Networks

text: Hello fellow Number crunchers As the headline suggest, I am looking for a java library for learning and inference of Bayesian Networks. I have already found some, but I am hoping for a recommendation. Requirements in a quick overview:   * written in Java (my overlord tells me that this is no point of discussion)   * configuration is possible via code (and not solely via a GUI).   * source code available   * project is still maintained   * the more powerful, the better Which one do you recommend ? PS: This should be cw.


title: Use coefficients of thin plate regression splines in a clustering method

text: There are numerous procedures for functional data clustering based on orthonormal basis functions. I have a series of models built with the GAMM models, using the `gamm()` from the mgcv package in R. For fitting a long-term trend, I use a thin plate regression spline. Next to that, I introduce a CAR1 model in the random component to correct for autocorrelation. For more info, see eg the paper of Simon Wood on thin plate regression splines or his book on GAM models. Now I'm a bit puzzled in how I get the correct coefficients out of the models. And I'm even less confident that the coefficients I can extract, are the ones I should use to cluster different models. A simple example, using:               #runnable code     require(mgcv)     require(nlme)     library(RLRsim)     library(RColorBrewer)          x1 <- 1:1000     x2 <- runif(1000,10,500)          fx1 <- -4*sin(x1/50)     fx2 <- -10*(x2)^(1/4)     y <- 60+ fx1 + fx2 + rnorm(1000,0,5)          test <- gamm(y~s(x1)+s(x2))     # end runnable code      Then I can construct the original basis using smoothCon :               #runnable code     um <- smoothCon(s(x1),data=data.frame(x1=x1),              knots=NULL,absorb.cons=FALSE)     #end runnable code      Now,when I look at the basis functions I can extract using               # runnable code     X <- extract.lmeDesign(test$lme)$X     Z <- extract.lmeDesign(test$lme)$Z          op <- par(mfrow=c(2,5),mar=c(4,4,1,1))     plot(x1,X[,1],ylab="Basis function",xlab="X",type="l",lwd=2)     plot(x1,X[,2],ylab="Basis function",xlab="X",type="l",lwd=2)     plot(x1,Z[,8],ylab="Basis function",xlab="X",type="l",lwd=2)     plot(x1,Z[,7],ylab="Basis function",xlab="X",type="l",lwd=2)     plot(x1,Z[,6],ylab="Basis function",xlab="X",type="l",lwd=2)     plot(x1,Z[,5],ylab="Basis function",xlab="X",type="l",lwd=2)     plot(x1,Z[,4],ylab="Basis function",xlab="X",type="l",lwd=2)     plot(x1,Z[,3],ylab="Basis function",xlab="X",type="l",lwd=2)     plot(x1,Z[,2],ylab="Basis function",xlab="X",type="l",lwd=2)     plot(x1,Z[,1],ylab="Basis function",xlab="X",type="l",lwd=2)     par(op)     # end runnable code      they look already quite different. I can get the final coefficients used to build the smoother by               #runnable code     Fcoef <- test$lme$coef$fixed     Rcoef <- unlist(test$lme$coef$random)     #end runnable code      but I'm far from sure these are the coefficients I look for. I fear I can't just use those coefficients as data in a clustering procedure. I would really like to know which coefficients are used to transform the basis functions from the ones I get with `smoothCon()` to the ones I extract from the lme-part of the gamm-object. And if possible, where I can find them. I've read the related articles, but somehow I fail to figure it out myself. All help is appreciated.


title: Which test to find out best concentration (the one having maximum effect)?

text: Hi my apologies first, I'm a biologist and not so good in statistics. In my study I'm studying the effect of concentration of feed on growth of a certain specimen. I have with me the different concentrations and the growth data for each concentration. What do you suggest is the best statistical test to find out the best concentration (i.e the one having maximum growth). Thanks in advance.


title: Interpreting correlation from two linear mixed-effect models

text: I have a linear mixed-effect model which I hope will answer the question of whether an increase in the frequency of use of one word leads to an increase of the frequency of use of that word by another person in a conversation, factoring out random effects of subject and topic of conversation. The basic model I've come up with looks like this:               Linear mixed model fit by REML      Formula: relative.sents.A ~ relative.sents.B + (1 | id.A) + (1 | abstract)         Data: nonzero          AIC    BIC logLik deviance REMLdev      -698.8 -683.9  354.4   -722.6  -708.8     Random effects:      Groups   Name        Variance   Std.Dev.       id.A     (Intercept) 1.0790e-04 0.0103877      abstract (Intercept) 3.0966e-05 0.0055647      Residual             2.9675e-04 0.0172263     Number of obs: 146, groups: id.A, 97; abstract, 52          Fixed effects:                      Estimate Std. Error t value     (Intercept)      0.017260   0.003046   5.667     relative.sents.B 0.428808   0.080050   5.357          Correlation of Fixed Effects:                 (Intr)     rltv.snts.B -0.742      The "dependent" variable is relative frequency of use by one person, and the fixed variable is relative frequency of use by another. I decided to see what the R^2 would be:               > cor(nonzero$relative.sents.A, fitted(nonzero.lmer))^2     [1] 0.6705905      To see what proportion of this is due to the fixed effect, I made a new model with only the random effects:               > summary(r.only.lmer)     Linear mixed model fit by REML      Formula: relative.sents.A ~ 1 + (1 | id.A) + (1 | abstract)         Data: nonzero          AIC    BIC logLik deviance REMLdev      -678.2 -666.3  343.1   -696.7  -686.2     Random effects:      Groups   Name        Variance   Std.Dev.       id.A     (Intercept) 1.2868e-04 0.0113435      abstract (Intercept) 7.8525e-06 0.0028022      Residual             3.7643e-04 0.0194017     Number of obs: 146, groups: id.A, 97; abstract, 52          Fixed effects:                 Estimate Std. Error t value     (Intercept) 0.029149   0.002088   13.96      ...and tried the same thing:               > cor(nonzero$relative.sents.A, fitted(r.only.lmer))^2     [1] 0.6882534      To my surprise, without that fixed effect, R^2 seems to increase! Does this mean my model is useless? If so, any suggestions on what might be wrong? Or am I somehow misinterpreting these results?


title: Bootstrapping interaction effects with boot() and the effects package - environment problems?

text: UPDATE AT BOTTOM - workaround and exciting new roadblock! UPDATE #2: SOLVED. Details at bottom. This is a combination stats/programming question, so I apologize if it should go somewhere else. I'm using the effects package and ggplot2 to visualize interaction effects from a multivariate linear model. I'm using effect() to extract the values to plot in ggplot. My response variables are very non-multivariate-normal, so I've used bootstrapping to estimate confidence intervals for all the regression coefficients so far. I'd therefore also like to bootstrap the confidence intervals for plotting the interaction effects rather than using the CIs supplied by effect(). I'm having trouble with this part - I keep getting an error that the input to my boot function is not found:                Error in is.data.frame(data) : object 'd' not found       I suspect this is caused by an issue with different environments and namespaces and, um... stuff. I have only the most rudimentary understanding of that part. I've looked around for solutions, and tried adapting what I've found: attaching the data frame that has my predictors, assigning values and functions to different environments and/or using do.call within the function. But frankly I'm just thrashing about, and haven't had any luck. I'd be grateful for any guidance. I'm open to other approaches (i.e. that don't use the effect package) toward the goal of boostrapped CIs for interaction plots. I'm not clear on the best way to holding the other variables at 'reasonable' values, or on whether those values should shift along the x axis of the plot, etc. Here is a reproducible example of what I've been doing so far:               # Some response variables:     y1 <- rlnorm(n=100, mean=0, sd=2)     y2 <- rlnorm(n=100, mean=20, sd=2)     y3 <- rnorm(n=100, mean=20, sd=3)     resp <- cbind(y1, y2, y3)          # Some predictor variables:     x1 <- sample(0:1, size=100, replace=TRUE)     x2 <- rnorm(n=100, mean=0, sd=3)     x3 <- factor(sample(c("Thanks", "For", "Your", "Help"), size=100, replace=TRUE))     pred <- data.frame(cbind(x1, x2, x3))          require(effects)     require(boot)          # fitting the model     mod <- lm(resp ~ .^2, data=pred)           # custom function so effect() will return a data frame instead of an eff object     effectdf <- function(...) {       suppressWarnings(as.data.frame(effect(...))) }          # make sure everything is working so far     effectdf(c("x1:x2"), mod, response="y1")     effectdf(c("x1:x3"), mod, response="y2")     # effects seems to work well, and you can see that column 3, fit,      # is the one I'd like to bootstrap and obtain CIs          # make the boot function     bf <- function(data, indices, formula, interaction, response) {       d <- data[indices,]       fit <- lm(formula, data=d)       effects <- effectdf(interaction, fit, response)       eff.fit <- effects[, 3]     }          x12_y1 <- boot(formula=resp ~ .^2,                      data=pred,                     R=99,                     statistic=bf,                     interaction="x1:x2",                     response="y1")          # Error in is.data.frame(data) : object 'd' not found       Thanks so much for your time and attention UPDATE: I haven't increased my effective understanding of the problem, really, but I have found a workaround. By offloading the parameters from the call to boot(), and placing them directly into the custom boot function, I found a way around whatever namespace problem was happening. I also have to run it as multiple univariate model - which is how effect() does it anyway. I _think_ the problem has to do with the way effect() looks for the original data to refit the multivariate response as univariate when it calculates effects. So now it looks like this:               # make the boot function         bf <- function(data, indices) {           d <- data[indices,]           fit <- lm(y1 ~ .^2, data=d)           effects <- effect("x1:x2", fit)           eff.fit <- effects[, 3]         }              x12_y1 <- boot(data=pred,  R=99, statistic=bf)      It doesn't throw the same error anymore. Now it throws the dreaded > Error in boot(data = pred, R = 999, statistic = bf) : number of items to > replace is not a multiple of replacement length Which apparently I can only solve by removing factors from my regression model - which is not an option. Why can I boot CIs for regression coefficients but not for these???? Any insight is greatly appreciated. UPDATE #2: As Diane Cook would later suggest (thank you for chiming in!), I used debug() and figured out the error was occurring when model.frame() was called. Googling my error message and model.frame led me to believe it was a problem with NAs generated by the resampling of unbalanced factors. I changed the na.action in lm() in the custom boot function from na.omit to na.exclude and that fixed it! I've seen many posts about the error message for boot, but none mentioning this solution. So let it be said plain and simple: if boot is throwing the error "number of items to replace is not a multiple of replacement length," it may help to change the way the functions in your custom boot function are dealing with NAs.


title: Linear Regression and Spatial-Autocorrelation

text: I want to predict Tree Heights in a certain area using some variables obtained through remote sensing. Like approximate Biomass, etc. I want to first use a linear regression (I know it's not the best idea but it is a must step for my project). I wanted to know how badly can spatial autocorrelation affect it and what is the easiest way to correct this if it is even possible. I'm doing everything in R by the way. Thank you in advanced. J.


title: How to interpret parameter estimates correlated with the intercept parameter estimate?

text: I have performed a regression in SAS and extracted the estimated correlation matrix of the parameters, which includes the intercept. One of my variable parameters has a strong correlation to the intercept parameter. I have not seen this before and want to be sure I understand what it means (if anything). Can anyone provide an intuitive explanation of what it means for a variable parameter estimate to be correlated to the intercept parameter estimate?


title: How to visualize changes in time series variance

text: I am exploring a time series with the aim of detecting obvious changes in variance over time (I have little experience in time series analysis). The hypothesis is that the series has become more stable over time (stable=decreased dispersion). I am wondering whether there are any strategies to visually assess how a measure of dispersion (does not have to be variance) changes as the series progresses. I understand there are tests for the constant variance assumption. My outcome is a daily measure captured over 10 years, with changes between months being the relevant measure. As a first step, I have plotted the monthly MAD (median absolute deviations) over time. However, I'm not sure that's a sensible approach.


title: Modelling samples with complex unsure distribution origin using PyMC

text: I've been learning PyMC for the last couple of days and have been trying to model the following problem: I have got a number of effectors that regulate targets with known values. For some targets the effector is known, but for others they are only partially known. The total number of effectors and targets is large, so I can't simply enumerate them in my code (or use a Bernoulli distribution with uniform prior). You can think of the targets to be samples from **one** effector (Normal) distribution. As an example, let's assume we have got the effectors E1 and E2 and the targets T1, T2, and T3. We know that:               E1 -> T1, T2     E2 -> T2, T3      which means that E1 is regulating T1 and E2 is regulating T3 (i.e., T1 is sampled from E1 and T3 is sampled from E2). For T2, we don't know if it is regulated by/sampled from E1 or E2, so I think it makes sense to use a categorical distribution with a Dirichlet prior. One idea was to model the data set analogous to this question and use one observed `Normal` with a categorical combination of `mu` and `tau` for each effector. In this case:   * for T1, it would be `1*mu(E1) + 0*mu(E2)`   * for T3, it would be `0*mu(E1) + 1*mu(E2)`   * for T2, it would be a Dirichlet with `tau=(1,1)` But how can I combine this into one framework? An obvious solution would be to use `tau=(1,0)` and `tau=(0,1)` for T1 and T3, but I can't assign zero probabilities in a Dirichlet. I guess I have to modify the `Categorical` to return `0` for not regulating effectors and a Dirichlet of all possible effectors (or modify the Dirichlet?) but I'm not sure how to do that. Another possibility would be to model each `E` as a `Normal` and assign uncertain origin to each value of `T`, but I don't know if that's possible. Any ideas?


title: Three level mixed-effects (crossed and nested) model in stata

text: I am trying to fit a 3 level random and crossed effects model to a continuous outcome to estimate the respective variance components of each level of the 3 (for doctors)   * Level 3: Physicians ID   * Level 2: Patient's ID   * Level 1: stage ( 1 and 2) For this study, each doctors have predicted the mortality rate for all the patients in two times. Thus clinician and patient is crossed effect and patient and stage is nested ( I think). I presumed that Physicians id is crossed with patient and patient is nested with stage. I tried the crossed random effects model for this data Model1:               xtmixed  mortality_est || doctorsid:R. patient|| doctorsid:R.stage      Model2:               xtmixed  mortality_est || doctorsid:R. patient|| doctorsid:R. stage|| patient:          xtmixed  mortality_est ||stage:||doctorsid:R.patient||z:, variance      but did not converge Any suggestion will be very appreciated.


title: Multivariate densities for whose we know the marginal HPD intervals

text: I am investigating the performance of a pipeline to compute marginal HPD intervals of multivariate densities (typically posterior densities). In this context, I am looking for some densities for whose these quantities are known in order to compare what the pipeline gets to the ground truth values. I would like the densities to exceed dimension 10. Does anybody have an idea?


title: Use predicted values with or without random part to plot Residuals with binnedplot of a logistic regression in glmer (lme4 package) in R?

text: Which binnedplot of the glmer should I use to check the model? The residuals against the predicted values without random part(REform=NA) or residuals against the predicted values with random part(REform=NULL)? I have one binary response variable (y.10) derived from one continuous variable with around 50 to 75% of zeros. I want to model the probability to exceed the limit of 10. For this example I used only one predictor "fragments" which is transformed by taking the logarithm to get a normal distribution an later a better fit . All variables are measured in tree regions (region). Within this regions are different plots (plot) and a set of samples were taken from some objects (object). To inspect the residuals I used binnedplot like discribed in the answer of the question: Unexpected residuals plot of mixed linear model using lmer (lme4 package) in R. To save calculation time with very complex models I modeled at first with glm {stats} and based on this results the model with less variables with glmer{lme4}. Doing this I could observe a big difference in the binnedplot of residuals. To examination the differences I created this example with only one variable. Like you can see in the picture bellow the models of glm and glmer without the random part show a very similar behavior. At the end the random part is not for interest. I need the random part only during model selection. Which binnedplot of the glmer should I use to check the model? The residuals against the predicted values without random part(REform=NA) or residuals against the predicted values with random part(REform=NULL)? The code and the resulting picture is given here:               fit.glm=glm(y.10 ~ x.t , data=data, family="binomial")     fit.glmer=glmer(y.10 ~ x.t + (1|region) + (1|plot) + (1|object),                data=data, family="binomial")          y.glm=predict(fit.glm, type ="response")     y.glmer=predict(fit.glmer,REform=NA,type ="response")      y.glmer.ran=predict(fit.glmer,REform=NULL,type ="response")          par(mfrow=c(2,2))     plot(y.10~x, data=data, type="n", main="Models")     points(y.glmer.ran~data$x, col="green", pch=4, cex=0.5)         lines(y.glm~data$x, col="red")     lines(y.glmer~data$x, col="darkgreen")          binnedplot(fitted(fit.glm),resid(fit.glm), main="Binned residual plot glm")     binnedplot(y.glmer.ran,resid(fit.glmer), main="Binned residual plot glmer(REform=NULL)")     binnedplot(y.glmer,resid(fit.glmer), main="Binned residual plot glmer(REform=NA)")      ![Shows in the first first Plot the glm fitted model \(red\) and glmer fitted models with  \(green crosses\) and without \(dark green line\) random part. The other plots show the corresponding residuals.](http://i.stack.imgur.com/HoW04.jpg)


title: Simplifying numerous dummy variables

text: I'm looking into a regression analysis to compare the time it takes to award public contracts across different countries but holding some other variables constant is proving challenging. The set of variables I'd like to hold constant relates to the type of tender and this is coded in the dataset according to a hierarchy of business activities along the following lines:   * Computing Services (code 1)      * Software development (code 1-1) \-- web development (code 1-1-1) \-- ...     * Support (code 1) \-- ...   * Telecoms (code 2)      * .. (2-1) \-- .. (2-1-1)      * .. (2-2) (actually the Common Procurement Vocabulary) I could create a dummy variable for all but one of the top-level categories (there are about 30). This would be quite interpretable, but I can see that a lot of explanatory information may be lost if I do this. Within each category there are some activities that lean towards major projects and other that lean towards smaller chunks of work. Is there a recommended technique to help with this problem? Thanks,


title: Why does the model change when using relevel?

text: When computing regression models with R, I regularly use the relevel function to get my model to give me results for the other level, too. I noticed that sometimes, but not often, this changed the model in the sense that levels of other factors that were significant before the relevelling are not any more. Is this inherent to relevelling or exceptional and maybe due to some problem with my data? Does it show that my data likely does not meet one of the prerequisites of linear models? Related to that, is it alright if I use relevel, recompute my model, and then report significance values from both models in my article? If significance differs between the two models for a certain factor, I suppose I should then go with one that is less optimistic? I suppose my question betrays that I don't know enough about lm to grasp the need for a base level. I thought I understood it pretty well ;) Somehow none of the introductions I read explained that point, or I was too daft to grasp it. So if someone could direct me to a site where the point of having base levels in lm is explained or explain it themselves, that would be great, too! Edit: Here's a minimal example:               library(datasets)     sprays<-OrchardSprays     model<-lm(decrease~treatment+rowpos+colpos,data=sprays)     summary(model)      Part of the summary says               treatmentC    20.625      9.731   2.120  0.03866 *      So if treatment == C this has significant positive influence on 'decrease'. Now I relevel 'treatment' to B to find out what influence treatment == A has:               sprays$treatment<-relevel(sprays$treatment,"B")     summary(model)      And now treatment == C is not significant in this new model:               treatmentC    17.625      9.731   1.811  0.07567 .      Sorry for posting in the wrong place! Can I move my question to stats statexchange or should I open a new one there?


title: How can I compare my model to a technically invalid model?

text: I've created nice little nonlinear model relating survival probability to length in salmon. I fit it assuming binomial errors and minimizing the negative log likelihood. I've been asked to compare it to someone else's model, where they binned the data and fit a straight line to it. However, the lowest bin includes the long left tail of the length distribution, and would predict 0 (or negative) chance of survival for those fish, were they not lumped into a bin with higher length average---but some of those fish do survive. That said, for some data sets, the linear model does quite well on the binned data. I'd like to compare these models, but I can't use AIC because the linear model's invalidity makes its AIC explode. I could truncate the data--it is a very small proportion of the data, or I could bin the data and calculate an AIC for my model assuming normal errors, but I don't really feel great about either of those. Are there other options, or are these choices not so bad?


title: Details regarding the delete-a-group jackknife

text: I was reading a paper by Phillip S. Kott on DAGJK: The delete-a-group jackknife. Journal of Official Statistics, 17 (4):521-526. (full text is freely available) I don't have much of a survey/sampling background so I'm having a bit of trouble understanding some of the paper.   1. What exactly is a primary sample unit (PSU)? An example would be incredibly helpful.   2. In the second section, Kott assumes that $t=\sum_{h=1}^{H}\sum_{j=1}^{n_h}t_{hj}$ and he defines $q_{hj}=t_{hj}-t_{h+}$, where $t_{h+}=\sum t_{hg}/n_h$ with the summation going over all PSUs in $h$. How did he arrive at the variance formula, $\textrm{Var}(t_{h+})=\frac{n_h}{n_h-1}\sum_{j=1}^{n_h}q_{hj}^2$?


title: Using a PCA to reduce response variables or multivariate multiple regression?

text: Does it make sense to use a PCA (principal component analysis) on a set of response Y variables and then conduct a multiple regression, or carry out a multivariate multiple regression all response variables as they are? I have 4 response variables and 2 predictor variables, which are all continuous.


title: How to combine states of viterbi path?

text: I am using HMM in my project. After Viterbi, I can get a sequence, such as 0000000000011111011111100000000000111111111111....... The 0 in subsequence 111110111111 may be wrongly decoded. Because in real situation, the same state should occur together. Is there any algorithm to solve this problem? In addition, how can calculate the probability of obeserving a subsequence. Thanks


title: Guidelines for splitting positively-skewed DV for Linear Modeling?

text: I am creating a multiple Linear Regression model using a Dependent variable which is extremely positively-skewed: My skewness statistic is 27.610, std. deviation is 2832.139, range of values from 0 - 100,000 in a population of 234,270 (over 95% of my cases have values at or below 500.00). Is splitting the file by the DV in order to reduce skewness and then modeling separately a reasonable course of action? Are there rules or guidelines for splitting the DV I should follow? For example, taking a certain number or percentage by Std. Deviation to inform where I should make the split? Note: I'm new at this - be gentle!


title: Computing the steady state probability vector of a random walk on $\{0, 1, \dots, n\}$

text: Suppose we have a random walk on $\\{0, \dots, n\\}$ with transition probabilities $$P(x, x + 1) = p \\\ P(x, x - 1) = 1 - p$$ for $1 \le x \le n -1$, $P(0, 1) = a$, $P(0, 0) = 1 - a$, $P(n, n - 1) = b$, and $P(n, n) = 1- b$ I need to compute the steady state vector $\pi$ for this random walk. First I wrote down the five difference equations $$\pi(x + 1)(1-p) + \pi(x - 1)p = \pi(x) \ \ \ \ \ \ \ 2 \le x \le n -2$$ $$\pi(0)(1 - a) + \pi(1)(1 - p) = \pi(0)$$ $$\pi(0)a + \pi(2)(1 - p) = \pi(1)$$ $$ \pi(n-1)p + \pi(n)(1 - b) = \pi(n)$$ $$\pi(n - 2)p + \pi(n)b = \pi(n-1)$$ and we also have the requirement $$\sum_{k=0}^n \pi(x) = 1$$ The general solution to the first equation is $$\pi(x) = k_1 + k_2\bigg(\frac{p}{p-1}\bigg)^x , \ \ \, p \ne 1/2 $$ $$\pi(x) = k_1 + k_2x, \ \ \ \ p = 1/2$$ I am unsure what do to next. Any help is appreciated.


title: Selecting the best subset of variables for parsimonious binary logistic regression models

text: In addition to PROC VARCLUS, randomForest, glmnet, and assessing multicollinearity among potential predictor variables (without regards to the outcome of interest), I am seeking other methods of variable selection in lieu of using stepwise methods for building more parsimonious binary logistic regression models (containing 8 to 12 variables to predict outcomes such as loan payment/default or current/late payment history) from a wide array of potential predictor variables (500+ variables, 200k+ records). Below I have included an R script using FSelector to select the 8 highest "ranked" variables:               library(FSelector)     fit <- information.gain(outcome ~ ., dataset)     fit2 <- cutoff.k(fit,8)     reducedmodel <- as.simple.formula(fit2,"outcome")     print(reducedmodel)      I have two questions regarding this script and the `FSelector` algorithm in general:   1. Is the `information.gain` criteria in the above script synonymous with Kullback-Leibler divergence? If so, can someone explain this in more layman terms than Wikipedia as I am relatively new to this area of statistics and would like to start off with the right idea of this concept as I may likely use this approach a great deal in the future?   2. Is this a valid approach, if there is such a thing as a valid approach, to select a desired number of variables for a binary logistic regression model (e.g., selecting the 8 highest "ranked" variables for use in a parsimonious model)? If not, can you provide an alternative approach to do so? Any insight or references regarding this topic and/or these questions will be greatly appreciated!


title: Are there default functions for discrete uniform distributions in R?

text: Most standard distributions in R have a family of commands - pdf/pmf, cdf/cmf, quantile, random deviates (for example- dnorm, pnorm, qnorm, rnorm). I know it's easy enough to make use of some standard commands to reproduce these functions for the discrete uniform distributions, but is there already a preferred built-in family of functions for modeling discrete uniform distributions in R that I'm unaware of?


title: Pairwise Mahalanobis distance in R

text: I'm trying to calculate a Mahalanobis-type pairwise distance matrix in R. I have 33 individuals, each with 10 variables. The idea is to get a distance matrix D, where $$D_{i,j}=(\mathbf{X}_i-\mathbf{X}_j)W^{-1}(\mathbf{X}_i-\mathbf{X}_j)^T$$ However I haven't been able build proper code for it.


title: How to interpret F- and p-value in ANOVA?

text: I am new to statistics and I currently deal with ANOVA. I carry out an ANOVA test in R using               aov(dependendVar ~ IndependendVar)      I get – among others – an F-value and a p-value. My null hypothesis ($H_0$) is that all group means are equal. There is a lot of information available on how F is calculated, but I don't know how to read an F-statistic and how F and p are connected. So, my questions are:   1. How do I determine the critical F-value for rejecting $H_0$?   2. Does each F have a corresponding p-value, so they both mean basically the same? (e.g., if $p<0.05$, then $H_0$ is rejected)


title: Interpret F-Value in ANOVA

text: When we conduct an ANOVA, we get F-Value and P-Value. If P-value if smaller than our alpha level of .05 for example, we reject our null hypothesis. What is the importance of F-Value that is obtained in the table as output? How it should be interpreted especially when conducting a within-subject ANOVA test? Thanks


title: Algorithms for predicting a couple points in the future

text: I'm familiar with supervised learning algorithms like regression and neural networks which look at a bunch of input points and learn a function which outputs a value (the value varying depending on whether the algo is a classifier, logistic regression, or standard regression)..but I am facing a problem now where each observation I have is a couple of data points (TAQ data), and then after a fixed amount of these chronologically sorted events there is a spike and what appears to follow an exponential decay mean reversion. Every such instance described here is one of my observations and my goal is to predict/model this exponential decay mean reversion. What type of algorithms would do something like this? time series models a-la ARCH/GARCH followed by a prediction look-ahead of N steps after training them or what else? Thank you, any type of suggestions or advice/references are greatly appreciated.


title: Fit a robust regression line using an MM-estimator in R

text: **Context.** I'd like to fit a regression line to study to relation between some response variable $y$ and some continuous covariate $x$. Because of the presence of bad leverage points, I have opted for an MM-estimator instead of the usual LS-estimator. **Methodology.** Basically, MM-estimation is M-estimation initialised by an S-estimator. Hence, two loss functions have to be picked. I have chosen the widely used Tukey Biweight's loss function $\rho ( u ) = \left\\{ \begin{array}{ll} 1 - \left[ 1 - \left( \tfrac{u}{k} \right)^{2} \right]^{3} & \textrm{if } | u | \leq k \\\ 1 & \textrm{if } | u | > k, \end{array} \right.$ with $k = 1.548$ at the preliminary S-estimator (which gives a breakdown point equal to $50 \%$), and with $k = 2.697$ at the M-estimation step (to guarantee $70\%$ Gaussian efficiency). I'd like to use R to fit my robust regression line. **Question.**               library(MASS)     rlm(y~x,          method="MM",         k0=1.548, c=2.697,         maxit=50)        * Is my code consistent with the previous paragraph?    * Would you use other optional arguments? **EDIT.** Following my discussion with @Jason Morgan, I realise that my previous code is wrong. (@Jason Morgan: Thank you very much for this!) However, I am still not convinced by his proposal. Instead, here is what I propose now:               library(robustbase)     lmrob(y~x,            tuning.chi=1.548, tuning.psi=2.697)      I think it sticks to the methodology now. **Do you agree?** Thanks!


title: Is it possible to use the caret package in R with non-numerical data?

text:
I was wondering if it is possible to use the caret package with non numerical data. I know, for example, if I want to use a simple linear regression `lm` I could have a factor variable for classification. However, caret blows up if I attempt this. I'm also following the step outlined here The caret Package: A Unifed Interface for Predictive Models for illustration I'm attempting to run `stepDurationlm <- train (x= trainDescr, y=trainClass, method="lm")` on               str(trainDescr)     'data.frame':   589235 obs. of  2 variables:      $ Anon.Student.Id    : Factor w/ 574 levels "02i5jCrfQK","02ZjVTxC34",..: 7 7 7 7 7 7 7 7 7 7 ...
      $ Step.Duration..sec.: num  5 5 5 4 5 5 5 5 4 4 ...      I get               Error in train.default(x = trainDescr, y = trainClass, method = "lm") :      All predictors must be numeric for this model. Use the formula interface: train(formula, data)      alternatively, could anyone explain how to have a test set for model performance in R? That's what's motivating me to get caret working.


title: Sum of subpopulation estimates not equal to estimate of total. Is this correct?

text: I am studying a population that consists of several sub-populations. I need models for the entire population and each of the sub-populations. Strangely, the total for the estimated variable of interest across the sub-populations differs from the estimated variable of interest for the population as a whole. **Some details.** The variable of interest (VOI) is the number of people who have done something. The observations are for several years and for several countries -- this is panel data. I defined the dependent variable (DV) as the logit of the ratio (VOI / population of the country). This transformation makes the most sense to me and produces the best fit for VOI. I calculate the estimated VOI by estimating DV using the model and then undoing the transformation. For each sub-population, I found the models using stepwise regression with fixed effects. I need to predict what would happen to VOI if the independent variables had certain values. The predicted VOI is the estimated VOI with the independent variables set to those values. **The problem** is with the predicted VOI. Logically, the sum of VOI across the sub-populations should equal to the VOI for the population as a whole, since these are numbers of people. In practice, the sum of the predicted VOI for the sub-populations is **very** different from the overall predicted VOI. I don't think that this is wrong per se. In each case, I am getting the best prediction conditional on my knowledge (or lack of it) of the sub-population. Is this correct? But even if it's not wrong, it looks bad. People who are reading the paper will complain. Ideas?


title: sampling for estimation - using random numbers - homework help

text: I have a figure like so :![enter image description here](http://i.stack.imgur.com/Xs7a0.png) Now the question asks to 1) generate 100 samples of iid - 2D uniform random variables in the unit- square. 2) count how many samples generated fall within the quarter unit- circle centered at the origin. What I've done so far: **approach 1** generate two vectors of random samples using `runif` in R.               samples_x(100)     samples_y(100)      and then I'm finding out if the area of the rectangle `(samples_x * samples_y)` is less than the area of the quarter circle **approach 2** generate two vectors of random samples using runif in R and then `samples_x ^ 2 + samples_y ^ 2` is less than or equal to radius square (`1`). Is either train of thought correct? 3) the question asks later,   * use these samples generated to estimate the area of the inscribed quarter circle.   * use the estimated area to estimate the value of pi. Which I have no idea as to how to solve. I've been looking at estimation of population but it's mostly count the samples in a small area and divide by the total area. By that method, would that mean I should estimate using samples falling within the circle and divide by 100 to estimate area? any pointers, help ,suggestions, opinions? much appreciated.


title: What is an "Unpaired Bland-Altman plot"?

text: I found this quote: > Unpaired Bland-Altman plots were used to check systematic or bias errors in > the SPECT myocardial uptake coefficients with respect to the reference > standard PET data. in the following paper: ALHASSEN, Fares, et al. Myocardial blood flow measurement with a conventional dual-head SPECT/CT with spatiotemporal iterative reconstructions-a clinical feasibility study. _American journal of nuclear medicine and molecular imaging_ , 2014, 4.1: 53. My knowledge of "Bland-Altman plot" is the following: I have $n$ distinct objects and two instruments, A and B; the instruments measure a property $w$ of the objects. $\mathbf{w}_A\in{\mathbb{R}}^n$, $\mathbf{w}_A=[w_{1A},w_{2A},\ldots,w_{nA}]$ and $w_{iA}$ is the measurement of the $i$-th object made with instrument A. $\mathbf{w}_B\in{\mathbb{R}}^n$, $\mathbf{w}_B=[w_{1B},w_{2B},\ldots,w_{nB}]$ and $w_{iB}$ is the measurement of the $i$-th object made with instrument B. The plot is made by the $n$ points $p_i(x_i,y_i)$ where $\mathbf{x}\in{\mathbb{R}}^n$, $\mathbf{x}=\frac{\mathbf{w}_A+\mathbf{w}_B}{2}=[x_1,x_2,\ldots,x_n]$ $\mathbf{y}\in{\mathbb{R}}^n$, $\mathbf{y}=\mathbf{w}_A-\mathbf{w}_B=[y_1,y_2,\ldots,y_n]$ So, what is an "Unpaired Bland-Altman plot" and how it's made?


title: Estimation of logit-normal parameters (II)

text: This question is related to the below question. Estimation of logit-normal parameters In the previous question, I showed only simple part of the whole problem. The log-likelihood function that I would like to maximize is $ n_1 log \int_{-\infty}^{\infty} f_1(y,\mu,\sigma) dy + n_2 log \int_{-\infty}^{\infty} f_2(y,\mu,\sigma) dy + n_3 log \int_{-\infty}^{\infty} f_3(y,\mu,\sigma) dy $ $n_1,n_2,n_3$ are known. $f_1,f_2,f_3$ are as follows. $ f_1 = \frac{exp(2y)}{(1+exp(y))^2} \frac{1}{\sqrt{2\pi}\sigma} exp(-\frac{(y-\mu)^2}{2\sigma^2} )$ $ f_2 = \frac{exp(y)}{(1+exp(y))^2} \frac{1}{\sqrt{2\pi}\sigma} exp(-\frac{(y-\mu)^2}{2\sigma^2}) $ $ f_3 = \frac{1}{(1+exp(y))^2} \frac{1}{\sqrt{2\pi}\sigma} exp(-\frac{(y-\mu)^2}{2\sigma^2}) $ How can I find $\mu$ and $\sigma$ that maximize the above equation?


title: Clustering time-shifted sales time-series

text: I need to perform clustering and classification of time series of weekly sales of different products. My data are weekly sales of different products in differest areas (stores). The challenges on this problem are:   \- Time-series are usually short: 10-52points(weeks).   \- Time-series may have a lot of zeros - sparce data. Products do not sell every week.   \- Not all products start to sell on the same date. This can result in time- shifted time-series. Even the same typical lifecycle of a product can be time- shifted in calendar along different stores.   \- Sales may have noise such as random events, promotions etc. A sample of data is like this: 20140105,prod1,store1,5   20140119,prod1,store1,10   20140126,prod1,store1,2   ....   20140105,prod1,store2,2   20140112,prod1,store2,3   ....   20140112,prod2,store3,4   20140126,prod2,store3,7   Can somebody share any insight on how to do this? Is it good to use a method such as DTW to compare time-series?If so, how am I going to handle the timeshifts? As for the implementation R seems a good way to go. Which packages would you recommend?


title: Identification of useless questions from a questionnaire

text: I'm developing a questionnaire. To improve its reliability and validity I want to use statistical methods. I want to eliminate questions whose answers are always the same. This means that nearly all participants gave the same answers on those questions. Now my questions are:   1. What's the technical term for such useless questions whose answers are always the same, independently from the context of use?   2. What are methods to identify such questions?


title: Coding categorical variables for regression

text: I'm not sure of the best way to code my categorical predictor variable for use in a hierarchical regression in order to test my specific hypothesis. This categorical variable has 3 levels representing 3 groups. I want to compare group 1 to group 2, group 1 to group 3 and group 2 to group 3. I know that for dummy coding I create k-1 variables, so 2 dummy variables in my case and code these variables with 0s and 1s while choosing one level of the categorical variable to be a reference category. However, I'm not sure this is the best way of making the comparisons I wish to make as it appears I could only compare each group to the reference category, am I correct? So if group 3 was the reference category I could compare group 1 to group 3 and group 2 to group 3 but I could not compare group 1 to group 2. What alternative method of coding should I use to make these comparisons? My regression model will also contain continuous variables. I'm an undergrad psychology student and statistics are not my strong point simple answers would be best for me. I use SPSS. Thank you!


title: Modelling non-stationary random walks

text: I am trying to model some data using a random walk, but after the standard increment and log (financial data) transformations for stationarity have found that, over long time frames, there is still not constant variance. What is the best way to now proceed? Or, are there further transformations I could consider to address the non-stationarity the standard transforms didnt address.


title: Compare event count in two different time periods

text: We have over 10 years of data on a number of hospitals that we regularly transport patients to. This past February, a new hospital (1830) opened that has substantially altered our transport patterns, especially taking away from the previous primary destination (1247).               Destination   Dec     Jan     Feb     Mar     Apr     1247          117     140     131      65      77     1017           75      69      54      45      31     1818           56      58      43      46      39     1830            -       -      30     125      86           I want to evaluate the overall impact (improved in-service times, reduced mileage/vehicle costs, etc...) of the new facility. What is the appropriate manner to evaluate the new facility? Should I compare equal periods prior to and after it opened? Or, can I evaluate the prior year and the period since it opened? Any help is appreciated. EDIT: I'm not sure that I have too much seasonal/cyclical variance in my long- term data. Here's a chart that has the previous 5 year history (each point is a 2-week period). I've highlighted total patients transported (with 26 period moving average), and the destination hospital in question. ![Patient Transport History](http://i.stack.imgur.com/mEd0A.png)


title: How to correct for generated regressor bias?

text: Dear Stack Exchange heroes, For my thesis I am writing a paper on the financial crisis. In my model, I use two regressions, which look like this: $$CONF = α + β_1 DEF_t + β_2 DIV_t + β_3 INF_t + β_4 IP_t + β_5 CONS_t + β_6 TBILL_t + β_7 URATE_t + β_8 NBER_t + β_9 DEF_{t-1} + β_{10} DIV_{t-1} + β_{11} INF_{t-1} + β_{12} IP_{t-1} + β_{13} CONS_{t-1} + β_{14} TBILL_{t-1} + β_{15} URATE_{t-1} + β_{16} NBER_{ t-1} + ε_t $$ Where $CONF$ is confidence index data and the rest are macroeconomic variables. $$ΔS_t = α_1 + λΔε_{t-1} + η_t$$ Where $S$ is a stock price index and $Δε$ uses the residuals from the previous equation. I use Eviews to run these equations and it works fine, but apparently I need to correct for the fact that $Δε$ is a generated regressor. The Pagan (1984) paper is often referred to for this, but after attempting to read it I'm still clueless on what to adjust in my covariance matrix to fix for this bias. Is there any program that can do this for me, or someone who can explain this in easy steps? Any help is appreciated :)


title: Loglikelihood and unusual transformations with binary response models

text: I was just doing some thought on the transformation for likelihoods to log likelihoods for binary response models and I realized what I am sure many people have realized before that the transformation is non-standard. For instance, given a response pattern $x \in 0,1$ and a parameter set $\theta$ with $P$ defined as the probability $x=1$, the likelihood of a binary response is given as $$l(\theta;x)=xP+(1-x)(1-P)$$ It is fairly straightforward to break the above separately as $l(\theta;x=1)=P$ and $l(\theta;x=0)=(1-P)$. The next step people take for various reasons is to take the log likelihood $$ll(\theta;x)=x\log(P)+(1-x)\log(1-P)$$ Which we can easily see is equivalent to taking the log of $l(\theta;x)$ observing that x can only be $1$ or $0$. However, this is a non-standard transformation ($\log(ab + cd) \ne a\log(b) + c\log(d)$ generally) and only achieved through reasoning (which I do not have a problem with). I am just wondering if there is some law, rule, or lemma in which the above transformation falls into?


title: Different Likert scales used in statistical anlysis

text: I am using different Likert Scales in questionnaire. Three dependent variables have Likert scales of 4 points, 2 points and 2 points. Three independent variables have 3 points scales. One moderating variable has 3 points Likert Scale whereas the other moderating variable has three questions with 4 points, 3 points and 3 points scales. If I use average of each variable in statistical analysis to find correlation, regression and moderated regression, is it allowed and results are considered valid or not?


title: Using factors scores in a multiple regression

text: Can factor scores be used as dependent variables in regression? Also, I want to see if demos (e.g. age, SES etc) predict performance in five different cognitive ability tasks. For each task I have the accuracy and speed of response. Would it be worthwhile attempting to run a factor analysis on these ten variables. So I could use these factor scores as my dependent variable and the personality variables as my predictors in multiple regressions. Or would it be more logical to run a multiple regression for each composite score (one for every cognitive task)? Any help would be appreciated!


title: Regression analysis with factor scores as explanatory variable

text: I intend to use factor scores as derived from exploratory factor analysis in subsequent multivariate regression analysis, as an explanatory variable. I've read in multiple books/papers that Structural Equation Modeling is the appropriate method to model relationships AMONG LATENT VARIABLES (using factor scores being incorrect). I haven't, so far, found a negative answer as to whether it is a problem to use factor scores solely as an explanatory variable, my dependent variable being observed. Does anyone have a clear answer regarding this issue?


title: Advice/literature on combining items with different response scales into composite scales?

text: Let's say I have some self-report items measured on a 5-point Likert scale (Strongly Disagree to Strongly Agree) and other items measured on a 4-point Likert scale (Never, Rarely, Sometimes, Often). Can anyone point me to literature (or practical advice) on combining these items into a composite scale? Let's assume for the sake of argument that we have some empirical evidence that items should be combined. Some ideas: **1\. Sum raw scores** Pro: Easy Con: Max responses on the 5-point scale (originally 4's) drive up the total scale score more than max responses on the 4-point scale (originally 3's). **2\. Rescale and sum** Put all of the items on a 0-1 scale and sum. So the 4-point items (0,1,2,3) would be multiplied by (4/3)/4 and the 5-point items (0,1,2,3,4) would be multiplied by 1/4, resulting in possible values of (0,.33,.66,1) and (0,.25,.50,.75,1), respectively. This way max responses on the 5-point scale (originally 4's) would not drive up the total scale score more than max responses on the 4-point scale (originally 3's). Pro: Items would have equal weighting. (could be a con, depending on your perspective). Con: Ignores differences in variability between items on different metrics? **3\. Standardize and sum** A related approach would be to standardize all of the items (z score) and then sum. Pro: Addresses differences in variability between items on different metrics Con: Total scale score becomes less interpretable and sample-specific. The latter makes it hard to benchmark as a measure to be used in other settings/other samples. **4\. PCA or other data reduction** 4a. EFA to get factor loadings. Multiply scaled items by factor loadings. 4b. PCA to get score of first principal component. Pro: Items weighted by influence. Con: Same as #2. EFA-derived scores could vary a lot depending on rotation/extraction choices. Some would not advise on ordinal data. * * * Overall: I like #2 because it seems easier to compare results across different samples. Thoughts? Alternative ideas or concerns about the ideas presented?


title: How can I replace this condition by a probability?

text: I want to see if a datapoint x should (or not) be assigned to a nearest component y using the following condition: if ($d > T$) then {do not assign x to y}. With $d = distance(x,y)$ and $T = \bar{d}+\sigma$ where $\bar{d}$ is the mean distance of y to datapoints that were already assigned to it in the past, and $\sigma$ the associated stdev. I want to replace this condition (i.e. ($d > T$)) by an equivalent expression which is expressed as a probability $p$. How can I do that ? I've tried for exemple $p = exp(\frac{-d^2}{2\sigma`^2}) / (\sigma`\sqrt{2\pi})$ but it does not behave like the first condition unless we manually find a good value for $\sigma`$ Note: it does not really matter if it is not a true probability (that sum to one).


title: What is the accepted concept behind measuring the contribution of subindexes to aggregate level indexes?

text: I am trying to find an additive decomposition of a Tornqvist index - to measure the contribution of individual industries to sector-level productivity indexes - and I can find formulae for this, but not the explanation for what these are based on.   My initial intuition for the contribution of a particular industry to the overall sector's index is that we would calculate the index with and without the industry. If I do this I can derive a relatively simple decomposition. From discussions with someone I work with, though, it appears that the more accepted way of thinking about contributions to an index is leaving the industry (or product, or whatever the subindex is) in the index, but with a 0 movement. Can anyone point me to any literature on why/if this is the concept that should be applied?


title: How do betting houses determine betting odds for sports?

text: Let's take football (soccer) for example. There are 3 possible outcomes, home win, draw, away win. I took a random game from bet365                 Turkey vs Ukraine     hwin, draw, awin     2.20  3.40  3.20      So for investment of 100\$ on given result, you either loose 100\$ or win: 220\$, 340\$ or 320\$. Their probability assessment doesn't add up to 100%, they take extra 5%-12%, but how did they come to these numbers (2.20, 3.40, 3.20)? Is it the betting patterns of people that bet, for example, if 90% of people put money on Turkey, `hwin` coefficient would be lower, or is it some kind of a calculation?   The problem with calculations is that the sample is very poor, national teams play very few games in long periods of time, between whole range of teams of different strengths, many outside parameters contribute, like injuries, current form and motivation of individual players, etc.   Is their strategy for national championships any different, you could find more regularity as the games are more frequently played, though 4 national league games per month isn't that much (and are also played home/away, which are two very different things).   So basically, the question is on what do they rely the most, how do they come up to these numbers, is it the calculation, betting patterns of other players, combinations etc?      One sub question, if other gamblers have a strong influence how coefficients are put, it seems to me that such assessments would be with significant error. I don't know if you can tell the difference between 65% and 70% for a given outcome, but that difference to me is indistinguishable. To be clear, I believe that Turkey in given example is a favorite, mostly because they play at home, but are their chances for a win 45% or 55% is too abstract, if they played against Monaco national team, then I'd give you a probability for the win with much more confidence.


title: generating covariance matrices from multiple priors

text: In many optimisation problems, one typically uses many forms of regularisations over the parameters that is being estimated. For example, a typical cost function (to maximise) may look like this: $$ C = f(w) - \lambda_1 r_1(w) - \lambda_2 r_2(w) $$ Here $f(w)$ can be interpreted as some likelihood term and $r_1$ and $r_2$ enforce some sort of structure on parameters $w$. $\lambda_1$ and $\lambda_2$ are the respective weights for these penalty/regularisation terms. In the probabilistic framework, one can imagine modelling these with two zero mean multivariate normal probability distributions: $P(w|\lambda_1) \sim \mathcal{N}(0, (\lambda_1\Sigma_1))$ and $P(w|\lambda_2) \sim \mathcal{N}(0, (\lambda_2\Sigma_2))$ Now, my questions are: So, applying Bayes law, assuming the observed data is denoted by y, can I write something like: $$ P(w|y) \propto P(y|w) \times P(w|\lambda_1) \times P(w|\lambda_2) $$ Is it possible to combine these two distributions $P(w|\lambda_1)$ and $P(w|\lambda_2)$ into one prior as I have defined them (with the $\lambda$ scale factors for the covariance matrices)? Doing some search online, there is some literature on Gaussian process priors. Is that something that might be relevant here? I just need some starting intuitive explanation of whether it is something I should look into as I need to use many priors like this over my parameters (each prior controlling some specific aspect of $w$).


title: Is it correct to use 'Ln' instead of 'ln' for natural logarithm?

text: In some research papers, authors use 'Ln' for natural logarithm instead of 'ln'. Is it correct?


title: Unknown p-value calculation

text: I was recently debugging an R script and I found something very weird, the author defined their own p-value function               pval <- function(x, y){         if (x+y<20) { # x + y is small, requires R.basic             p1<- nChooseK(x+y,x) * 2^-(x+y+1);             p2<- nChooseK(x+y,y) * 2^-(x+y+1);             pvalue = max(p1, p2)         }         else { # if x+y is large, use approximation             log_p1 <- (x+y)*log(x+y) - x*log(x) - y*log(y) - (x+y+1)*log(2);             pvalue<-exp(log_p1);         }         return(pvalue)     }      Where X and Y are values positive values greater than 0. The <20 case seems to be a calculation for some kind of hypergeometric distribution (something similar to Fisher test?) and does anyone know what the other calculation is? As a sidenote, I am trying to optimize this code so trying to figure out the proper R function to call and replace this with. Edit: Paper detailing formula for p-value calculation can be found here (need to click pdf to see the formulas) Methods start on page 8 of the pdf and the formula in question can be found on page 9 under (1). The distribution they assume is a Poisson.


title: Testing equality of two binomial proportions proportion (one near 100 %)

text: I am testing equality of two binomial proportions 87/88 and 48/60. I use this online calculator and it noted that the Standard Normal approximation is not valid in this case. It seem quite strange to me. Is it some problem in testing equality of two binomial proportions if one of the proportions is very near to 100 %?


title: What are the limitations of graphical displays with a small sample size (n=18)

text: I have a small sample size (n=18) and am performing univariate analysis. I have graphical displays (stem and leaf plot, dotplots, boxplots and a cumulative distribution plot). I am trying to understand the limitations of the graphical displays for the small data set. Can anyone enlighten me?


title: How does Gower distance work with free text?

text: The Gower distance measure is a good measure for mixed-type data (i.e., data attributes can be qualitative, categorical, ordinal or binary). But can data attributes be free-text (e.g., names of people)? I understand that string edit distance is a better measure for comparing free-text, but I have a table of records, where every column (attribute) in the table has a different type. I want to compute a similarity matrix of the table (pairwise record comparisons) for clustering. Gower distance seems suitable, except I'm not sure if it can handle free-text strings. I understand that I can just perform Gower on all the attributes except the free-text string columns, and perhaps add the string edit distance as a normalized component to the overall distance measure. But I want to know if Gower can handle free-text strings. Or any other suggestions for mixed-type distance measures are appreciated.


title: Interval regression with sample selection in Stata

text: I'm wondering if anyone knows how to correct for sample selection when using the interval regression, preferably with Stata. I have firm-level bribery data of which 70% did not pay a bribe. For the rest of the firms which paid a bribe, I have information about the frequency of the bribery, namely: once, 2-5 times, 6-10 times, and more than 10 times. Due to the nature of the data, I don't think I can simply use 'intreg' command. Is there any way that I can control for selection with interval regression model please? Thank you in advance.


title: What type of test do I use?

text: I have created 39 question Likert scale survey and want to cross-reference results with 5 demographic questions. I plan on assigning values for each answer 1=strongly disagree, 2=disagree, 3=neutral, etc. I am required to put info into SPSS. I am rather confused on which test I should use: ANOVA, ANCOVA etc. Any thoughts?


title: How can I tell if modelling/testing with Binomial/Beta is correct?

text: My standard disclaimer: I'm a programmer, not a stats guy, mostly wanting to understand statistical tests for optimizing web apps. Let's say I want to measure monetization rate of users. "How many users monetize?". This is a yes/no question on each user. After reading Wikipedia, I might say, "This yes/no question is a Bernoulli variable. I can model groups of users as instances of the Binomial distribution, and use Binomial and Beta distributions to estimate _p_ and compare the groups of users." But, what if _p_ is not the same for every user in the group? Some users will always monetize, some users will never monetize... Some will monetize after their next payday. Even though the question is a yes/no question, is it sometimes wrong (or rather, how safe is it) to use Binomial/Beta distributions to model a collection this way? These distributions assume a constant _p_ , but when _p_ might vary for each user am I getting incorrect results?


title: How to evaluate results of linear regression

text: I have a linear regression problem. In short, I have a dataset, I divided it into two subsets. One subset is used to find the linear regression (training subset), another is used to evaluate it (evaluation subset). My question is how to evaluate the result of this linear regression after applying it to the evaluation subset of data? Here are the details: In the training subset, I do linear regression: $y = ax + b$, where $y$ is groundtruth (also known as target), $x$ is an independent variable. Then I found $a$ and $b$. ($x$ and $y$ are given in the training subset). Now, using $a$ and $b$ found above from the training subset, apply them to the evaluation subset, I found $y' = ax' + b$. In other words, these $y'$ are found from linear regression with $x'$. Now, in addition to $y'$, I also have $y$ from the evaluation set. How do I evaluate my result (how much $y'$ differ from $y$)? Any general mathematical model to do that? It needs to be some sort of math model/formula. I can think of different ways to do it, but they are all kinda ad-hoc or simple, but this is for a scientific work, so things that sound ad-hoc cannot be used here, unfortunately. Any idea?


title: The R code for showing histogram/plot of conditional distribution of Y|X

text: Would you please write the command/code for drawing conditional distribution of Y|X in R programming? I need this code to see if Y and X relationship is skewed or not, and see if I could have linearity assumption for them. Thank you


title: How to make glmnet give the same results as glm?

text: What parameterization to `glmnet` will give the same results as `glm`? (I'm mainly interested in logistic and linear regressions, if that matters.)


title: Minimal number of samples/conversions for statistical validity

text: We are measuring conversion rates (% of visitors who bought) on an e-commerce site. The test apply to a segment of visitors who meet specific criteria (for example people from a certain country). The people from the segment are divided into 2 groups. Part of them see a banner and the other don't (control group). Usually the control group is 30% of visitors. The test begin after the banner is shown to all users in the segment for a while so the data of exposed people extends much longer than the data of the control group. So at a given time we have for example X people exposed and Xb of them converted; likewise for Y people who were not exposed, Yb of them were converted. Y and Yb are much smaller than X and Xb. The conversion Rate for X is Xb/X and for Y is Yb/Y. My **first question** is how to determine statistical validity. We used Chi-square test to do it, and got results similar to those from this on-line calculator (implemented similar to table1 here). However, sometimes it looks like the number of purchases is extremely small yet the chi-square test says its valid (>95%). Here is a real life example:               X=189 Xb=1 Y=93 Yb=3     Conversion X= 0.5% Conversion Y=3.2%      Statistical confidence (chi-square based) 92.8%      Altough the confidence is below 95%, it seems too close for determining its confidence while there was only 1 conversion for X. My **second question** is then: Should there also be a requirement for a minimal number of conversions for the confidence to be valid? If so, how do we calculate it? Thanks!


title: What is meant by 'weak learner'?

text: Can anyone tell me what is meant by the phrase 'weak learner'? Is it supposed to be a weak hypothesis? I am confused about the relationship between a weak learner and a weak classifier. Are both the same or is there some difference? In the adaboost algorithm, `T=10`. What is meant by that? Why do we select `T=10`?


title: How to interpret R stl() output

text: I want to have a logical interpretation of the results of my stl analysis.                summary(BR.stl)              Call:         stl(x = ts, s.window = "per")               Time.series components:             seasonal              trend            remainder                  Min.   :-2.9127520   Min.   : 2.87902   Min.   :-32.97469            1st Qu.:-1.7160211   1st Qu.:10.31844   1st Qu.: -2.80863            Median : 0.0410429   Median :13.86636   Median : -1.08351            Mean   :-0.0257600   Mean   :15.90906   Mean   : -0.05876            3rd Qu.: 1.0959295   3rd Qu.:18.57809   3rd Qu.:  2.05679            Max.   : 2.1946197   Max.   :59.88744   Max.   : 18.82706            IQR:          STL.seasonal STL.trend STL.remainder data           2.812        8.260     4.865         9.875          %  28.5         83.6      49.3         100.0          Weights: all == 1               Other components: List of 5          $ win  : Named num [1:3] 1821 19 13              $ deg  : Named int [1:3] 0 1 1          $ jump : Named num [1:3] 183 2 2              $ inner: int 2          $ outer: int 0      What is the correct interpretation of the percentages of variation. They add to >100%, does this mean they can only have a relative interpretation similar to the bars on an stlplot? My perfect result would be to have a proportion of variation for each component of stl. Any help much appreciated.


title: Comparing correlation coefficients

text: I have two sets of data where I have at ~250.000 values for 78 and 35 samples. Some of the samples are members of a family and this may have an effect of the data. I have calculated pairwise correlation and it varies between 0.7 and 0.95 but I would like to know if there is significant difference in correlation coefficients intra vs inter family? What is the best way to do this? Thanks


title: Hypothesis testing on values that are functions of multiple distributions

text: I have a hypothetical experiment where I am comparing scores on some measure at 1 year to scores on the measure at 4 years. I use a non-paired t-test to see if there is a significant difference between the two. Then, I do:               (4 year mean) - (1 year mean)      Now, I want to take that value, a function of two different distributions that have different n and SD, and compare it to another value that was obtained in the same way from two similar distributions (1st & 4th year but under a different independent variable). How can I do this? And is it valid do 4yrMean - 1yrMean to obtain the value, or should I do p-value4 - p-value1 or something?


title: Fitting a GARCH(1,1) model

text: How much data is needed to properly fit a GARCH(1,1) model?


title: How to estimate time-per-product in a factory?

text: Workers in a factory are assembling items built from several parts. For each item they fetch the parts from the warehouse and then assemble the items. I suppose they need a certain time per item for assembly plus a certain time per part for fetching them. $D_{order} = n_{items} * D_{fetch} + n_{parts} * D_{assemble}$ I have a list of past orders which tells me   * $n_{items}$: how many items there are in that order   * $n_{parts}$: how many parts there are in that order and   * $D_{order}$: how long it took the worker to complete the order. I want to estimate, based on the number of items and parts, how long a future order will take. So I need to find out $D_{fetch}$ and $D_{assemble}$. Which is the best way to plot the data to estimate those two factors? Here is some sample data (order completion time on the Y-axis): This chart shows the time to complete an order vs. the overall number of parts in it: ![time_vs_parts](http://i.stack.imgur.com/CA5wM.jpg) This chart shows the time to complete an order vs. the number of assembled items in it: ![time_vs_items](http://i.stack.imgur.com/3qPbG.jpg) What I didn't tell before: Usually several of the items are the same product, so maybe there is a overhead per product for getting/reading the appropriate blueprints. So this chart shows the time to complete an order vs. the number of unique products in it: ![time_vs_unique_products](http://i.stack.imgur.com/cgQFf.jpg)


title: Proper way of using recurrent neural network for time series analysis

text: Recurrent neural networks differ from "regular" ones by the fact that they have a "memory" layer. Due to this layer, recurrent NN's are supposed to be useful in time series modelling. However, I'm not sure I understand correctly how to use them. Let's say I have the following time series (from left to right): `[0, 1, 2, 3, 4, 5, 6, 7]`, my goal is to predict `i`-th point using points `i-1` and `i-2` as an input (for each `i>2`). In a "regular", non-recurring ANN I would do process the data as follows: > >      target| input >           2| 1 0 >           3| 2 1 >           4| 3 2 >           5| 4 3 >           6| 5 4 >           7| 6 5 >   I would then create a net with two input and one output node and train it with the data above. How does one need to alter this process (if at all) in the case of recurrent networks?


title: ANOVA and multiple testing correction in gene screening

text: I want to discover which genes are expressed in only one of five treatments. This is my pipeline:   1. ANOVA between the five treatments    2. Holm multiple testing correction    3. Tukey for significant genes discovered in step 2 My question is: Should I also correct Tukey p-values for example multiplying the p-value by the number of significant ANOVA p-values (Bonferroni correction on the number of tests performed), or should I only correct at the level of ANOVA? Thanks a lot in advance, Rossella


title: Predicting user interaction with interface

text: If I have a GUI with a set of buttons and I want to predict what button the user will select next based on a set of features including the previous n button clicks and other meta data (e.g. user's name), what area of ML is this? Does this count as classification (i.e. classifying the next move) with the x classes being the x buttons and the features as mentioned above?


title: What kind of statistic analysis fits in my data?

text: I am not mathematician, but I want to do a small statistic analysis in my data. I have data from basketball matches in the following format:               Points for Team A / Points for Team B / Winner (A,B or tie) / Theoretical percentage for Team A winning / Theoretical percentage for Team B winning / Theoretical percentage for tie / Theoretical points for Team A / Theoretical points for Team B      What I want to check is if the theoretical percentages follow the real results and somehow check the same for the points and theoretical points for each team. I would like also to check the following: If theoretical points for Team A is 50 and for team B is 45 (Team A win) and the real results are Team A:100 and Team B:20, even if Team A win, it is far away of what the theoretical was.


title: Model validation for multilevel logistic regrssion?

text: I have designed a multilevel logistic model using PROC GLIMMIX in SAS 9.3 for hierachical data based on pupil attainment where level-two is the school the pupil attends. I'm quite sure that my model is a good fit, I have a generalized chi-square/DF value of 1.07 but have no other way of validating the model. Are there any other methods of checking that the model is a good fit? I have read a lot about possibly extending the ROC curve to multilevel models but believe this to be invalid. Is there any multilevel equivalent of the Hosmer & Lemeshow test statistic used for single level logistic regression?


title: Multivariate regression with weighted least squares in python?

text: I have a multivariate regression problem that I need to solve using the weighted least squares method. In particular, I have a dataset `X` which is a 2D array. It consists of a number of observations, `n`, and each observation is represented by one row. Each observation also consists of a number of features, `m`. So that means each row has `m` columns. Therefore my dataset `X` is a `n`×`m` array. Given a test data observation, multivariate regression should produce a function that predicts the response vector `y`, which is a 2D array as well. This function will consist of `m` coefficients, i.e. one coefficient/parameter for each of the `m` features of the test input. This solution is already available here: Statsmodels' WLS, except that they don't support a 2D response vector yet. In other words, when I fit the data, I have to provide my dataset `X`, but can only provide a 1D array as the response `y`. In addition, I also need a 2D weights vector, similar in dimension to the response vector `y`. Is there a Python implementation of WLS multivariate regression where `y` and the `weights` can be 2D vectors? Or if not a direct implementation, can any of the existing packages be used as an implementation somehow, by a small amount of adjustment? ## Edit To make my question clearer, these are the parameters I would give in and the results I would need to get out: Input:   * `X` : a 2D dataset, like 10x3, which is 10 observations with 3 features each.   * `y` : which is also a 2D vector, in this case 10x2. In other words, a 2-value response vector for each observation. (I'm doing classification and there are two possible classes).   * `weights`: a 2D response vector which is also 10x2, like `y`. The `10` above is an arbitrary number of rows. Ultimately, however many observations I have, that's just what the number of rows is going to be, for all vectors above. Output I need:   * the coefficients of the regression. Given that my response and weight vectors are 2D, I believe the coefficients would also be a 2D array, probably 3x2 or 2x3.


title: Hypothesis test for Bernoulli

text: Let $X_{1},..,X_{36}$ be a sample from a Bernoulli distribution with parameter $p$. The sample proportion is $\frac{1}{3}$. Consider a normal approximation test $H_{0}:p=0.5$ vs $H_{1}:p\neq 0.5$, with confidence level $0.95$. How will the criterion change if $H_{1}$ is replaced by $p=1/3$?


title: Example of invariance of Wald Test

text: In order to make us understand that the Wald doesn't work with a transformation of your parameter at H0, we have to test the following $H_0$'s: A) $H_0: exp( \mu)=0.62$ B) $H_0: \mu=ln(0.62)$ I can write the Wald test statistic for case B, which is $W_B= (\mu - ln(0.62))^2 / Var( \mu)$ which matches the result we get out of Stata, But what is the formula for $W_A$? I think you have to apply the delta method, but I can't figure it out unfortunately. Thank you in advance!


title: From standard deviation to linear regression error

text: I am trying to understand how to obtain linear regression error for a football shaped scatter plot. Suppose I have two distributions correlated by coefficient $r$. Moreover these two distributions have mean $\mu_1$ and $\mu_2$ and standard deviation $\sigma_1$ and $\sigma_2$. How to obtain error for linear regression starting from this data?


title: ARIMA analysis - decomposition

text: Decomposition can be used to generate de-seasonalised data for input into an ARIMA analysis. But why is this done? Does it have a large impact on the analysis and therefore resulting forecast??


title: Calculating Hellinger Divergence from Results of Kernel Density Estimates in Matlab

text: Using the ksdensity function in matlab returns a density estimation in the form of 2 vectors _f_ and _xi_. Where _f_ are the density values and xi the corresponding points for the density values. How do I calculate the hellinger distance between 2 density estimations based on their _f_ and _xi_ ?


title: Does PCA followed by LDA make sense?

text: This is a question about classification. I am a neuroscience student with little experience of classification methods and I'd be grateful for any advice about the best way to implement a linear classifier (LDA) on this data. I have a magnetoencephalography dataset, recorded from people as they perform a cognitive task. This has the following properties:   * 306 channels of data, but a preprocessing step has reduced the dimensionality to 64, and then reprojecting the data onto the sensors.   * The data is fine-grained in the time domain (1000Hz)   * The data is chopped up into short segments termed 'trials'. These correspond to the cognitive task the subjects were performing. The trials can differ, for example, in trial type A (of which I will have ~50) the subject may have been asked to pay attention to stimuli on the left hand side of space, and in trial type B (of which I have also 50) to pay attention to stimuli on the right hand side of space). I want to classify my data at particular timepoints within trials (e.g. 0.5 seconds after subjects are told to attend left/right), training the classifier to discriminate between trials of type A and type B. The feature vector for each observation (i.e. trial) is the vector of instantaneous activity at each sensor at that timepoint. I have more features than trials (306 sensors) so I either need to do feature selection or use regularized LDA (or use something like Hastie's sparse discriminant analysis which seems to me to do both). This I could do, BUT it seems like it throws away the information about the statistical structure of the data contained in the datapoints NOT from the exact time I am trying to classify. Also, I know the data has a lower dimensionality than 306 - less than 50 components will likely capture the vast majority of the variance in the data. I was thinking therefore about using a dimensionality reduction step, probably PCA, before passing the reduced-dimensionality data to an unregularized LDA or naive bayes classifier. The idea is that the dimensionality reduction step exploits the fact that I have a large amount of data sampled over time. And that's where I got confused. PCA projects the data onto orthogonal dimensions, so does it makes sense to then do LDA (which estimates a covariance matrix in order to use information about the correlation between the features)? Or should I do naive bayes following the dimensionalty reduction? In which case, I could probably have just done naive bayes to start with, given that it's not sensitive to the number of features in the way that LDA is. If anyone can advise on a good approach here, given the particular structure of this data, I'd be very grateful. The key question is: does PCA followed by LDA make sense? Thanks! George


title: Cross-validation accuracy interpretation -( accuracy of 100%

text: Here the setup: I have 90% of the data for training and other 10% for testing. I am doing stratified cross-validation on the 90% Tranining. It is a 10-class dataset. I am using LibSVM for that. When doing 10 fold cross-validation for tuning the Hypermarameters (C in the C-SVM) I get accuracies of 100%. Basically something like this:                       Training with 0.03125   -  Cross Validation Accuracy = 68.5097%             Training with 0.12500   -  Cross Validation Accuracy = 98.3%             Training with 0.50000   -  Cross Validation Accuracy = 100%             Training with 2.00000   -  Cross Validation Accuracy = 100%             Training with 8.00000   -  Cross Validation Accuracy = 100%      It is ok to have 100% accuracy in the cross-validation on the TRAINING data? In this case, should I chose C = 0.5 as the best hyper-paramter? or instead should I move away from parsmeters that ge me 100% in the cross- validation? and why? if I don't take those with 100%, should I take what 98%? 90%? Thanks,


title: Quality assurance and quality control (QA/QC) guidelines for a database

text: ## Background I am overseeing the input of data from primary literature into a database. The data entry process is error prone, particularly because users must interpret experimental design, extract data from graphics and tables, and transform results to standardized units. Data are input into a MySQL database through a web interface. Over 10k data points from > 20 variables, > 100 species, and > 500 citations have been included so far. I need to run checks of the quality of not only the variable data, but also the data contained in lookup tables, such as the species associated with each data point, the location of the study, etc. Data entry is ongoing, so QA/QC will need to be run intermittently. The data have not yet been publicly released, but we are planning to release them in the next few months. Currently, my QA/QC involves three steps:   1. a second user checks each data point,   2. visually inspect histogram each variable for outliers    3. users report questionable data after spurious results are obtained ## Questions   1. Are there guidelines that I can use for developing a robust QA/QC procedure for this database?   2. The first step is the most time consuming; is there anything that I can do to make this more efficient?


title: Sequential clustering algorithm

text: I want to cluster elements in array. The crucial difference from a normal clustering algorithm is that the order of elements is significant. For instance if we look at a simple sequence of numbers like this:               1.1, 1.2, 1.0, 3.3, 3.3, 2.9, 1.0, 1.1, 3.0, 2.8, 3.2      It is obvious that there are two clusters in there (1.1, 1.2, 1.0, 1.0, 1.1) and (3.3, 3.3, 2.9, 3.0, 2.8, 3.2). What I want is to find sequential groups of similar elements               (1.1, 1.2, 1.0), (3.3, 3.3, 2.9), (1.0, 1.1), (3.0, 2.8, 3.2)      4 in this case. Of course I can run some variant of a normal clustering algorithm and then split clusters according elements' indices, but there's probably a simpler way to do this. Is there any algorithm that I can use for this?


title: how to complement the results of cluster analysis with known groups

text: I have some prior knowledge of grouping, but this may be incorrect or is not sufficient as I need larger number of groups (i.e. subgroups). For example in the following data I have 3 groups in addition to two variables. I would like to use the group information (as prior knowledge) (here 3 groups) to create meaningful groups (here 9 groups/clusters). Is there a correct way to perform such analysis.               # Dummy data      group <- rep(1:3, each =3000)     X <- c(rnorm(1000, 0.1, 0.04), rnorm(1000,0.2, 0.04), rnorm(1000, 0.4, 0.02),            rnorm(1000, 0.4, 0.04), rnorm(1000,0.5, 0.08), rnorm(1000, 0.6, 0.12),             rnorm(1000, 0.7, 0.08), rnorm(1000,0.8, 0.1), rnorm(1000, 0.9, 0.06)     )          Y <-  c(rnorm(1000, 0.5, 0.04), rnorm(1000,0.6, 0.04), rnorm(1000, 0.7, 0.04),            rnorm(1000, 0.35, 0.12), rnorm(1000,0.45, 0.04), rnorm(1000, 0.3, 0.02),             rnorm(1000, 0.55, 0.09), rnorm(1000,0.65, 0.12), rnorm(1000, 0.65, 0.04)     )      Prior information of 3 clusters:               col = c("red", "cyan", "green")     plot(cbind(X,Y), col = col[group], pch = ".")      ![enter image description here](http://i.stack.imgur.com/v5xUI.jpg) Clustering analysis assuming 9 clusters.               cl <- kmeans(cbind(X,Y), 9)          colrs <- c("red","purple", "yellow", "tan", "pink", "cyan", "blue", "green", "black")     plot(cbind(X,Y), col = colrs[cl$cluster], pch = ".")      ![enter image description here](http://i.stack.imgur.com/tecHn.jpg)


title: logistic regression with related predictor variables

text: I have two simulations for physical ability testing that add up to a total test score, coded as dichotomous pass or fail 0, 1. Can I treat the two simulations continuous times as independent predictors in binary logistic regression even though they add to total score.


title: Alpha inflation in the context of usability studies

text: In a hypothetical usability experiment, two interface changes are tested against a baseline. The interface changes are cumulative, i.e. variant B is an extension of A. We measure usability according to its definition by measures for efficiency, effectivity and satisfaction. To make things worse, the interface variants are measured by having human participants process three different sets of tasks. If the hypothesis is formulated as "variant A improves the usability over baseline", the metric to test would be an index of the usability measures over all three sets of tasks. Since all three measures are subject to noise, a signal in one of the measures might be missed. As an alternative, we could examine each usability measure separately. Is it correct to assume that in the latter case we have to correct for alpha- inflation? If so, would it be better to state only few hypotheses a priori (e.g. "B improves efficiency over baseline"), report the results and examine other hypotheses in an exploratory analysis? Also, how about the three different sets of tasks? Is it correct to assume that examining the measures for each set of tasks separately would result in even higher alpha inflation? **tl;dr:** In a 3x3 design (variant, task set) measuring three variables, do we have to adjust for a family-wise error rate of $1-(1-p)^{3\cdot3\cdot3}$? If so, can we get around it by stating fewer, more clearly-define hypotheses and use the remaining hypotheses as a basis for an exploratory analysis?


title: Generating random samples from a given distribution

text: What are the different methods by which we can generate random samples from the existing distribution which are computationally efficient? One of the common methods given in statistics text is by applying a transformation called the inversion method. We start with a Uniform(0,1) distribution from which we can pick an arbitrary number, and then use the inverse of the CDF of the distribution (its quantile function) to find the member of the desired distribution corresponding to the randomly chosen number between 0 and 1. What other methods exist for doing the same?


title: Probability of what is being predicted in glm

text: In documentation to `glm` I read: " _For binomial and quasibinomial families the response can also be specified as a factor (when the first level denotes failure and all others success)_ " Does it mean that probability of failure or success is being modeled? I'm trying to apply simple logistic model to "german credit scoring" dataset where there are levels "good" and "bad". To get correct results (higher probability means higher likelihood of being good) I have to assume that `Failure=Good` and `Success=Bad`. This works, but it is really counterintuitive. I interpret this as - this will model probability of Failure (failed to be bad).               require(ggplot2)          german_data <- read.csv(file="http://archive.ics.uci.edu/ml/machine-learning-databases/statlog/german/german.data",                   sep=" ", header=FALSE)          names(german_data) <- c('ca_status','mob','credit_history','purpose','credit_amount','savings',     'present_employment_since','status_sex','installment_rate_income','other_debtors',     'present_residence_since','property','age','other_installment','housing','existing_credits',     'job','liable_maintenance_people','telephone','foreign_worker','gb')          str(german_data)          german_data$gb <- factor(german_data$gb,levels=c(2,1),labels=c("bad","good"))          levels(german_data$gb)[1]           table(german_data$gb)          model <- glm(data=german_data,formula=gb~.,family=binomial(link="logit"))          german_data$prob <- predict(model,newdata=german_data, type="response")          ggplot(data=german_data) + geom_boxplot(aes(y=prob,x=gb))  + coord_flip()      ![enter image description here](http://i.stack.imgur.com/kcbao.png)


title: Regression in projective space?

text: Is there a method for (nonlinear? kernelized?) regression of functions with output in projective space? That is, given a series of examples $x_i\in\mathbb{R}^n$ (or $x_i\in\mathbb{P}^n$) and $y_i\in\mathbb{P}^m$ (so, $\forall c\in\mathbb{R}, y_i=c\cdot y_i$) such that $f(x_i)=y_i$, can I write down some function that reasonably approximates $f$? Thanks!


title: PyMC: how can I define a function of two stochastic variables, with no closed-form distribution?

text: I'm learning PyMC and basically I have a random variable $Z = X + Y$ where (say) $X \sim \mathrm{Normal}(\theta_X)$ and $Y \sim \mathrm{Lognormal}(\theta_Y)$ and $Z$ has no simple closed-form distribution. Now I have observations $z_i,\,i=1...N$ of $Z$ and I want to infer $\theta_X$ and $\theta_Y$. What's the most straight-forward way of doing this with PyMC? If I had the distribution of $Z$ available, then I think I could do:               Z = DistZ('Z', param_x=theta_x, param_y=theta_y, value=z, observed=True)      and then do inference, but I don't know `DistZ`. It's also easy to define the sum as:               @pymc.deterministic     def z_sum(x=Y, Y=y):         return x + y      but then I don't think I can define an observed deterministic function. I _think_ I could do something like:               @pymc.stochastic(observed=True)     def z_sum(value=z, x=X, y=Y):         def logp(z, x):             # return log-likelihood      but I'm not clear on the details. I do know the joint likelihood $\mathcal{L}(z, x)$, but I was hoping it wouldn't be needed. I was able to do this with a custom Gibbs sampler (using the joint likelihood), but I'm looking for a more "elegant" solution with PyMC. * * * **EDIT** : found a similar question in the BUGS FAQ that says functions of random variables aren't supported. Not sure if that applies to PyMC, and what the standard approach is.


title: How to analyse this data obtained from a simple physics experiment on attractive forces?

text: I did a simple physics experiment that measures the attractive force a plate experiences towards the other plate as a function of the applied voltage and distance between the plates. Now I have to know whether the gathered data confirms the relation $F \propto V^2/d^2$. Suppose the following is what I gathered.               voltage\distance 8.000 10.000 12.000 14.000 16.000                  4.0   3.3    2.0    1.7    1.2    0.8                  6.0   8.1    4.8    3.8    2.5    1.9                  8.0  13.4    9.2    6.0    4.4    3.5                 10.0  22.9   14.2    9.4    7.0    5.6                 12.0  32.7   20.1   13.5   10.5    7.9      The uncertainties are $\pm 0.05$ for the voltage, $\pm 0.005$ for the distance, and $\pm 0.1$ for the force. As far as I know, there are a few ways to analyze the data with R.   1. `lm(log(force) ~ log(voltage) + log(distance))`   2. `lm(I(force * distance^2) ~ -1 + I(voltage^2))`   3. `aov(I(force * voltage^-2 * distance^2) ~ voltage * distance)` Which one should I use?


title: If x = y*y, and you know var(y), var(z), and cov(y,z), do I know cov(x,z)?

text: If I know that x = y*y, and I know a whole of statistics pertaining to y, such as the variance and its covariance with other variables, can I analytically solve for the variance and covariance of x? Also, if I know a bunch of statistics regarding x, what can I know about y analytically? Thanks for the help!