title: Moderated Mediation with categorical moderator

text: Hi I'm trying to make the following moderated mediation but I'm missing the final step on testing differences between groups. The code below tests the indirect effects for each level of the categorical variable, but I have no clue how to test if the indirect effect differs between categories. I want to test the indirect effect of "Spatial" on "MCUnd" through "TotCorr" in which "Condition" moderates the effect of "Spatial" on "TotCorr". "Spatial", "TotCorr" and "MCUnd" are all continuous and "Condition" is categorical with four categories. This is what I have:               # Now it's time to do the actual moderated mediation. First     # I created a model/object representing the treatment variable      # (spatial ability) to the mediator variable (dynamic). As it     # is a moderated mediation, SA is interacting with condition in this model          medmod.fit <- lm(TotCorr ~ Spatial * Condition, data=Masters)          # I then create a model/object for the effect from mediator     # to Y (Multiple Choice Understanding scores), again, as it     # is a modmed, variables are crossed with condition.          outmod.fit <- lm(MCUnd ~ TotCorr + Spatial * Condition +                              TotCorr * Condition, data=Masters)          # At this point I run one overall mediation analysis and four     # mediation analyses that split the data by condition. This is     # consistent with Tingley, Yamamoto, Koole, & Imai (2012)          medmod <- mediate(medmod.fit, outmod.fit, treat = "Spatial",                       mediator = "TotCorr")          medmod.cond1 <- mediate(medmod.fit, outmod.fit, treat = "Spatial",                             mediator = "TotCorr", covariates = list(Condition = 1))          medmod.cond2 <- mediate(medmod.fit, outmod.fit, treat = "Spatial",                             mediator = "TotCorr", covariates = list(Condition = 2))          medmod.cond3 <- mediate(medmod.fit, outmod.fit, treat = "Spatial",                             mediator = "TotCorr", covariates = list(Condition = 3))          medmod.cond4 <- mediate(medmod.fit, outmod.fit, treat = "Spatial",                             mediator = "TotCorr", covariates = list(Condition = 4))          #This is just an overall summary     summary(medmod)          # This is a summary of the mediation effects for each of the four     # conditions. If there is moderated mediation, the mediation effects     # should be different across conditions. I do find that there     # mediation effects in conditions 2 and 3(moderate effect), but no     # mediation effect in conditions 1 and 4.          summary(medmod.cond1)     summary(medmod.cond2)     summary(medmod.cond3)     summary(medmod.cond4)      Thank you!


title: Estimate the parameters of an ellipse in the presence of large outliers

text: I have a number of points $x_1,\ldots,x_m\in\mathbb{R}^n$ with weights $w_1,\ldots,w_m$ between 0 and 1. There is an ellipse which contains a very high concentration of points with weights close to one, but there are also many points scattered outside the ellipse with large weights, and not all the points in the ellipse have large weights. Still, from inspection, it's very obvious where the boundaries of this ellipse are. I'm trying to estimate $\mu\in\mathbb{R}^n$ and $C\succ0$ such that the ellipse is defined by all points $\left\{x\middle|(x-\mu)^T C (x-\mu) < 1\right\}$. I am trying to do this by solving: $\min_{\mu,C\succeq 0} \sum_{i=1}^m w_i L((x_i-\mu)^T C (x_i-\mu))$ where $L(x) = 0$ if $|x|<1$, $L(x) = |x|-1$ if $x\ge 1$. It might also be worth adding a penalty on the log determinant of $C$, which is equivalent to penalizing the volume of the ellipse. and I have two questions 1) Does this make sense as a robust way of finding a bounding ellipse? Is there a more natural objective given the problem I described? I know there are solvers out there for finding the minimum volume enclosing ellipse, but those don't include weights or outliers. 2) How would I actually go about implementing this optimization? I haven't seen much out there about combining semidefinite programming with hinge-type losses.


title: Problem with BootCV for coxph in pec after feature selection with glmnet (lasso)

text: I am attempting to evaluate the prediction error of a coxph model that was built after feature selection with glmnet. In the preprocessing stage I used na.omit (dataset) to remove NAs. I reconstructed all my factor variables into binary variables with dummies (using model.matrix) I then used glmnet lasso to fit a cox model and select the best performing features. Then I fit a coxph model using only these feature. When I try to evaluate the model using pec and a bootstrap I get an error that the prediction matrix has wrong dimensions. In the original dataset there are 313 variables. In the coxph model (glmnet.cox in my code) there are 313 df. However when I run pec the prediction matrix is noted as having 318 variables and therefore doesn't fit the testset which has 313 variables (+the status and time variables). It seems that pec does something to the variables while retraining the coxph model on the bootstrap samples. I tried setting singular.ok=FALSE in the coxph model without avail. Here is a summary of key variables (prior to restructuring for glmnet) Followed by the code I ran and the errors:               > summary (trainSet)          time              status         Age_at_Dx       CREATININE         Min.   :  0.1429   Min.   :0.0000   Min.   :14.81   Min.   :0.400        1st Qu.: 22.0714   1st Qu.:1.0000   1st Qu.:52.51   1st Qu.:0.800        Median : 52.9286   Median :1.0000   Median :63.79   Median :0.900        Mean   : 96.5415   Mean   :0.7826   Mean   :61.01   Mean   :1.042        3rd Qu.:154.6071   3rd Qu.:1.0000   3rd Qu.:73.01   3rd Qu.:1.125        Max.   :437.7143   Max.   :1.0000   Max.   :87.23   Max.   :6.000        Performance_Status    ALBUMIN          Cyto             WBC              Min.   :0.0000     Min.   :0.70   Min.   : 1.000   Min.   :  0.60        1st Qu.:1.0000     1st Qu.:3.00   1st Qu.: 4.000   1st Qu.:  3.20            Median :1.0000     Median :3.60   Median : 4.000   Median :  9.20        Mean   :0.9728     Mean   :3.48   Mean   : 6.802   Mean   : 28.35        3rd Qu.:1.0000     3rd Qu.:4.00   3rd Qu.:11.000   3rd Qu.: 33.10        Max.   :4.0000     Max.   :5.20   Max.   :17.000   Max.   :373.00         PRKCD_pT507             HGB            Maxblast          CD19             Min.   :-2.429605   Min.   : 5.400   Min.   : 0.00   Min.   : 0.000         1st Qu.:-0.627005   1st Qu.: 8.500   1st Qu.:24.00   1st Qu.: 1.000        Median :-0.013117   Median : 9.550   Median :40.00   Median : 2.000        Mean   : 0.006432   Mean   : 9.689   Mean   :43.74   Mean   : 6.312        3rd Qu.: 0.559782   3rd Qu.:10.800   3rd Qu.:65.25   3rd Qu.: 5.000        Max.   : 2.963658   Max.   :14.300   Max.   :98.00   Max.   :98.000            GAPDH                CD74              TP53               Fli1                Min.   :-2.391021   Min.   :-1.7902   Min.   :-1.64244   Min.   :-3.723143        1st Qu.:-0.662164   1st Qu.:-0.5502   1st Qu.:-0.53685   1st Qu.:-0.559783        Median : 0.003236   Median :-0.1546   Median :-0.10928   Median :-0.002154        Mean   : 0.015759   Mean   : 0.0235   Mean   :-0.02285   Mean   :-0.016555        3rd Qu.: 0.632472   3rd Qu.: 0.4759   3rd Qu.: 0.29869   3rd Qu.: 0.554521        Max.   : 3.512741   Max.   : 3.7226   Max.   : 5.39242   Max.   : 4.415324             LDH              SQSTM0            BM_BLAST         CCND3               Min.   :    8.4   Min.   :-1.56639   Min.   : 0.00   Min.   :-2.44772        1st Qu.:  562.8   1st Qu.:-0.48762   1st Qu.:31.00   1st Qu.:-0.59042        Median :  952.5   Median :-0.05746   Median :46.50   Median :-0.08412        Mean   : 1617.9   Mean   :-0.02272   Mean   :50.64   Mean   :-0.02506        3rd Qu.: 1596.8   3rd Qu.: 0.33718   3rd Qu.:72.00   3rd Qu.: 0.48353        Max.   :36708.0   Max.   : 3.59456   Max.   :98.00   Max.   : 3.58761             GAB2                BAX               ABS_BLST           PRKAA1_2             Min.   :-3.201564   Min.   :-3.230737   Min.   :     0.0   Min.   :-1.90671        1st Qu.:-0.640279   1st Qu.:-0.593174   1st Qu.:    99.8   1st Qu.:-0.66896        Median : 0.007018   Median :-0.106097   Median :  1836.5   Median :-0.10586        Mean   :-0.013284   Mean   :-0.007051   Mean   : 15024.1   Mean   :-0.03531         3rd Qu.: 0.635609   3rd Qu.: 0.588098   3rd Qu.: 11178.0   3rd Qu.: 0.44906        Max.   : 2.396419   Max.   : 3.439942   Max.   :358080.0   Max.   : 3.60037           RPS6KB1              SPP1               MAPT                ARC                Min.   :-2.05490   Min.   :-1.39728   Min.   :-1.826101   Min.   :-2.67081        1st Qu.:-0.65736   1st Qu.:-0.55249   1st Qu.:-0.573470   1st Qu.:-0.63908        Median :-0.12007   Median :-0.15156   Median :-0.046306   Median :-0.07845        Mean   :-0.06115   Mean   : 0.02759   Mean   : 0.009604   Mean   :-0.01082        3rd Qu.: 0.49055   3rd Qu.: 0.35195   3rd Qu.: 0.495406   3rd Qu.: 0.63638         Max.   : 3.54709   Max.   : 4.44919   Max.   : 4.063568   Max.   : 2.52810         FOXO1_pT24_FOXO3_pT32      ATG7               SMAD5           RPS6_pS235_236            Min.   :-2.2506       Min.   :-2.916018   Min.   :-2.321250   Min.   :-3.229642         1st Qu.:-0.5590       1st Qu.:-0.621408   1st Qu.:-0.480924   1st Qu.:-0.633143         Median :-0.1281       Median : 0.003795   Median : 0.003186   Median :-0.003345         Mean   :-0.0121       Mean   :-0.029199   Mean   : 0.010663   Mean   :-0.015028         3rd Qu.: 0.4415       3rd Qu.: 0.590721   3rd Qu.: 0.502824   3rd Qu.: 0.660224         Max.   : 3.1239       Max.   : 2.757525   Max.   : 2.238310   Max.   : 2.816105             LSD1             HDAC2               SFN                  Min.   :-4.5434   Min.   :-2.40682   Min.   :-2.082592         1st Qu.:-0.5654   1st Qu.:-0.49165   1st Qu.:-0.554862         Median : 0.0477   Median : 0.02828   Median :-0.116248         Mean   :-0.0165   Mean   : 0.04641   Mean   : 0.006255         3rd Qu.: 0.6044   3rd Qu.: 0.47826   3rd Qu.: 0.506527         Max.   : 1.8736   Max.   : 4.22724   Max.   : 3.733445                      ##-----------------------------------------------------------------------------------         -------------------        ##                                    FEATURE SELECTION BY GLMNET-LASSO        ##----------------------------------------------------------------------------------     ---------------------        ##         ## This function takes a predictor matrix and a Surv obj and uses the glmnet lasso      regularization         ## method to select the most predictive features and create a coxph object        ## predict_matrix is the 'dataset' reformated to replace categorical variables to     binary with dummy          #creat Y (survival matrix) for glmnet       surv_obj <- Surv(trainSet$time,trainSet$status)                 ## tranform categorical variables into binary variables with dummy for trainSet       predict_matrix <- model.matrix(~ ., data=trainSet,                                       contrasts.arg = lapply (trainSet[,sapply(trainSet,     is.factor)], contrasts, contrasts=FALSE))            ## remove the status/time variables from the predictor matrix (x) for glmnet       predict_matrix <- subset (predict_matrix, select=c(-time,-status))            ## create a glmnet cox object using lasso regularization       glmnet.obj <- glmnet (predict_matrix, surv_obj, family="cox")            # find lambda for which dev.ratio is max        max.dev.index <- which.max(glmnet.obj$dev.ratio)            optimal.lambda <- glmnet.obj$lambda[max.dev.index]             # take beta for optimal lambda        optimal.beta  <- glmnet.obj$beta[,max.dev.index]             # find non zero beta coef        nonzero.coef <- abs(optimal.beta)>0        selectedBeta <- optimal.beta[nonzero.coef]             # take only covariates for which beta is not zero        selectedVar   <- predict_matrix[,nonzero.coef]            # create a dataframe for trainSet with time, status and selected variables in binary      representation for evaluation in pec      reformat_dataSet <- as.data.frame(cbind(surv_obj,selectedVar))            # create coxph object with pre-defined coefficients         glmnet.cox <- coxph(surv_obj ~ selectedVar,data = reformat_dataSet,model=TRUE,      singular.ok=FALSE,init=selectedBeta,iter=0)                ##------------------------------------------------------------------------------------   -------------------       ##                                    MODEL PERFORMANCE       ##-----------------------------------------------------------------------------------    --------------------       ##        set.seed(17743)       pec.f <- as.formula(Hist(time,status) ~ 1)             ## Brier score for the cox-glmnet model        brierGlmnet <- pec(list(glmnet.cox), data= reformat_dataSet , splitMethod="BootCV",     B=50, formula = pec.f)      Here is the result:               Prediction matrix has wrong dimensions:     368 rows and 318 columns.           --->(but there are only 313 predicting vars     in the dataset)      But requested are predicted probabilities for      138 subjects (rows) in newdata and 315 time points (columns) --- (315=313 predicting     vars+status+time)     This may happen when some covariate values are missing in newdata!?     Warning in FUN(1:2[[2L]], ...) :      Oddly, when I call print on the glmnet.cox object I get that it has 313 variables and not as pec outputs (318)               print(glmnet.cox)     Likelihood ratio test=0  on 313 df, p=1  n= 368, number of events= 288


title: Interpreting Coefficients from Logistic Mixed Effects Model

text: I am trying to use glmer to create logistic mixed effects models (LMEMs) and then run a model comparison between them. Respondents can click "Yes" (coded as 1) or "No" (0), which is used as the DV. They can make three responses to each stimulus across time. The fixed effects are viewing duration of a stimulus at the point of each rating (continuous, the main variable of interest), the rating number (1st, 2nd, 3rd, which will be highly correlated with cumulative viewing duration) and whether the true answer was "Yes" or "No" (which I'll call 'veracity'). I have specified my models as: **complex model** -- response = 1 + rating number (1-3) + cumulative viewing duration + veracity + (1 + cumulative duration | participant number) + (1 + cumulative viewing duration | stimulus item number) **simple model** -- response = 1 + rating number (1-3) + veracity + (1 + cumulative duration | participant number) + (1 + cumulative viewing duration | stimulus item number) The simple model is missing the cumulative viewing duration main effect. asking for anova(simple, complex) gives the following:               Data: data     Models:     simpleML: resp ~ 1 + responseno + veracity + (1 + cumdur + veracity | pno) +      simpleML:     (1 + cumdur | itemno)     complexML: resp ~ 1 + responseno + cumdur + veracity + (1 + cumdur + veracity |      complexML:     pno) + (1 + cumdur | itemno)               Df    AIC    BIC  logLik  Chisq Chi Df Pr(>Chisq)        simpleML  12 2248.1 2313.8 -1112.0                                 complexML 13 2241.1 2312.3 -1107.5 8.9558      1   0.002766 **     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       I understand this to mean that the addition of the cumulative viewing duration variable significantly enhances the model. I am interested in knowing how the odds of responding "Yes" changes with each additional unit of cumulative viewing duration (measured in seconds). So to do that, I figured I should look at the complex model in more detail: the beta will surely be reported there, and then I can use the exponent link function to give me a log odds (hopefully I'm right so far). So here's the complex model:               Generalized linear mixed model fit by the Laplace approximation      Formula: resp ~ 1 + responseno + cumdur + veracity + (1 + cumdur + veracity | pno) + (1 + cumdur | itemno)      Data: data        AIC  BIC logLik deviance      2241 2312  -1108     2215     Random effects:      Groups    Name        Variance   Std.Dev.  Corr                pno       (Intercept) 3.1359e-01 0.5599916                               cumdur      3.3268e-05 0.0057679 -0.014                        veracity    6.8532e-01 0.8278393 -0.900  0.063       itemno (Intercept) 3.5810e-01 0.5984182                               cumdur      7.7778e-05 0.0088192 -0.052             Number of obs: 1764, groups: pno, 49; itemno, 12          Fixed effects:                  Estimate Std. Error z value Pr(>|z|)     (Intercept)  0.326981   0.314613   1.039    0.299     responseno  -0.224734   0.148566  -1.513    0.130     cumdur      -0.001513   0.005310  -0.285    0.776 <<-------- **     veracity     0.387681   0.432779   0.896    0.370          Correlation of Fixed Effects:                (Intr) rspnsn cumdur     responseno -0.014                   cumdur     -0.144 -0.759            veracity   -0.702 -0.262  0.217      From this, I first notice the p-value is non-significant. I would have assumed the p-value here is a partial value: the effect of cumulative duration when the other variables are held at zero. So this worries me. Second, the beta is incredibly small, and the exponential of -0.001513 = 0.9985. In real terms, this means for each 1 second of cumulative viewing time, the odds of making a 'yes' judgment is lower by (1 - 0.9985 = ) 0.002%. Or for each 100 seconds of cumulative viewing time (which is probably about half the length of a stimulus item on average), the odds of making a 'yes' judgment is lower by 0.2%. This seems so incredibly small that the p=.776 makes sense to me. But then why does the Likelihood Ratio Test (i.e. "anova(complex, simple)") find a strong significant effect? I have tried googling for various answers, reading Bates' documentation and even asking work colleagues, but to no avail. Any help at all would be greatly appreciated! * * * ## EDIT - Possible Answer? Dave Vinson has been kind enough to help me out, and I now have further information. This information is entirely credited to him (and any mistakes entirely credited to my misunderstanding). I have for now ignored the random effects, and regressed response (yes/no, y-axis) on the rating number (1, 2 or 3, x-axis). This gives y = 0.61 - 0.05x, where x is the rating number and y the predicted response: *[apparently I cannot post images without sufficient reputation. Annoying. Hopefully the model above gives you some idea of the effect.] Having done that, I then calculated the residuals (y minus response). These residuals, the unexplained variance, were then regressed on the cumulative viewing duration, shown below. Residual is plotted on the Y-axis, cumulative duration on the X-axis. As can be seen, the slope is remarkably shallow. *[sorry, not able to post images.] Perhaps the structure of the mixed effects models was incorrect. Specifying a slope for cumulative duration in the simple model could be wrong, resulting in the variance that would normally be explained by the fixed effect of cumulative duration to start being accounted for in the random effect. Plausible? And how might I check for this?


title: Asking Mallet Language Tool Kit

text: I am working to Text Document. I try to filter the document using Mallet Language Tool Kit. My questions are. 1. I use macintosh os. should i install ant ? 2. How can i use the packages in this directory class/cc/mallet/classify on Java eclipse? Thank you.


title: Practical PCA tutorial with data

text: Searching the internet for PCA tutorial gives thousands of results (even videos). Many of the tutorials are very good. But I am not able to find any practical example where PCA is explained using some data-sets which I can use for demonstration. I need a tutorial which provides some small data set which is easy to plot (not 10000s of lines of data with 100s of dimension) before and after the PCA analysis and the difference/results can be clearly demonstrated. (I think a step-by-step practical example with a data having about some 100 lines and 3 dimension will be great). Do you have any suggestion??


title: Sample size for Cluster sample

text: I'veen asked to claculate a necessary sample size and I tought the best way to reduce costs and time is using cluster samplig protocol but I have few information and I don't have much experience with cluster sampling. The situation is as follows. I have a population of 83,201 individuals. I need to calculate a sample to have 90% confidence level. My population is divided in 5 clusters and the size of each cluster is: Cluster 1: 22248 Cluster 2: 14692 cluster 3: 7912 Cluster 4: 5086 Cluster 5: 33263 I want to know how many clusters i need to choose and how many individuals of each cluster to have a sample with 90% of confidence level but I don't have more data. I dont have variance of the population... Can someone help me to solve this?? Thank's a lot


title: Comparing two histograms

text: I would like to know what are the common techniques to compare two histograms? I have histogram of two images and I want to see are they similar or not meaning that is there any correlation between them or not. The histograms are for two different parts of tissue. ![enter image description here](http://i.stack.imgur.com/octv6.jpg)


title: How to find pdf of a joint distribution in R?

text: $F(x,y) =\frac{1}{6}(x^2\, y+x\, y^2)\,,\quad 0\leq x\leq 2,\, 0\leq y\leq 1$ Above is the joint distribution given,   1. how to find out cumulative distribution function of y?   2. how to obtain joint probability density function of x and y? I am a beginner in R, I know basic commands. thanks for the help


title: Bootstrapping results for p-values of trading rule profit significance

text: I want to obtain p-values for the significance of trading rule profits as described by LeBaron and Brock (1992). The author's use what is developed by Efron (1973) called Bootstrapping. It goes like this:   1. Run a trading rule (i.e. SMA crossovers) on a price series    2. Calculate the mean return generated by this rule and call it "ret.mean.original"   3. Sample from the vector of original returns with replacement and add these returns up, starting from the first value of the original series.   4. Run the trading rule on the "new price" series and calculate the mean return. Call it "ret.mean.new"   5. If ret.mean.new > ret.mean.original a "1" enters a vector called "Loop". If ret.mean.new < ret.mean.original a "0" enters Loop at the corresponding position.   6. Repeat steps 3-5 10.000 times   7. Sum up all the ones from "Loop" and divide by 10.000 This value is then the approximate p-value for the profitability of the trading rule under consideration. The H0 is of course "Zero-Profitability". So here is my code:               require(quantmod)               # Loads price data from Yahoo     require(TTR)                    # for the trading rules     require(PerformanceAnalytics)   # for charts      first we run a backtest on the real market price series to obtain "ret.mean.original"               # Load S&P 500 Index data from Yahoo! Finance     getSymbols("SPY", from="2000-01-01", to="2014-05-01")     prices = SPY$SPY.Adjusted          # Obtain log returns from prices     ret     = ROC(prices,type="continuous")     ret[1]  = 0          # define strategies     # mean reversion     rsi2  = RSI(prices,2)          # trend following     sma.5   = SMA(prices,5)     sma.200 = SMA(prices,200)          # Create buy/sell signals     # mean reversion     sig.mr2   = ifelse(rsi2 < 50, 1, -1)     sig.mr2   = lag(sig.mr2) # lag the signal so todays signal is applied tomorrow          # trend following     sig.tf.200.5  = ifelse(sma.200 < sma.5, 1, -1)     sig.tf.200.5   = lag(sig.tf.200.5)               # replace missing values with zeros     sig.mr2[is.na(sig.mr2)]   = 0     sig.tf.200.5[is.na(sig.tf.200.5)]  = 0               #strategy returns     ret.mr2   = (ret*sig.mr2)     ret.tf.200.5  = (ret*sig.tf.200.5)          # cut the series      ret.tf.200.5  = ret.tf.200.5[-c(1:200), ]          # chart the return series to see how it worked     returns = cbind(ret,                       ret.tf.200.5,                       ret.mr2)                       ret.tf.50.15)     colnames(returns) = c("Buy and Hold","TF.200.5","MR2")          chart.CumReturns(returns,geometric=F,legend.loc="top left")      until here everything works fine Now to the tricky part where I don't get whats going wrong:               # Start Loop here     N     = 200 # Number of simulations     Loop  = mat.or.vec(N,2)     for (i in 1:N){            # sample with replacement from return distribution of index        ret.new = (sample(ret, length(ret), replace = T, prob = NULL))       # demeaning returns       ret.new = ret.new-mean(ret)       # new price series starting at same value as original series       prices.new = xts(prices[[1]]*exp(cumsum(ret.new)))            # define strategies       # mean reversion       rsi2.new  = (RSI(prices.new,2))             # trend following        sma.5.new   = (SMA(prices.new,5))        sma.200.new = (SMA(prices.new,200))             # Create buy/sell signals        # mean reversion        sig.mr2.new   = (ifelse(rsi2.new < 50, 1, -1))        sig.mr2.new   = lag(sig.mr2.new,na.pad = TRUE)                   sig.tf.200.5.new  = (ifelse(sma.200.new < sma.5.new, 1, -1))        sig.tf.200.5.new   = lag(sig.tf.200.5.new,na.pad = TRUE)                  # replace missing values with zeros        sig.mr2.new[is.na(sig.mr2.new)]   = 0        sig.tf.200.5.new[is.na(sig.tf.200.5.new)]  = 0             #strategy returns        ret.mr2.new   = (ret.new*sig.mr2.new)        ret.tf.200.5.new  = (ret.new*sig.tf.200.5.new)             # cut the series        ret.tf.200.5.new  = ret.tf.200.5.new[-c(1:200), ]             Loop[i,1] = if (mean(ret.tf.200.5.new)  > mean(ret.tf.200.5)) {1}else{0}        Loop[i,2] = if (mean(ret.mr2.new)       > mean(ret.mr2)) {1}else{0}     }     #Loop           p.value.ret.tf.200.5  = sum(Loop[,1]/N)      p.value.ret.mr2  = sum(Loop[,2]/N)                # Results      p.value.ret.tf.200.5       p.value.ret.mr2           # plots simulated series      returns.new = cbind(ret.new,                     ret.tf.200.5.new,                     ret.mr2.new)           chart.CumReturns(returns.mr.new,geometric=F)      chart.CumReturns(returns.tf.new,geometric=F)      So there are two things that I don't understand. The first is: The p-value for the SMA crossover profit is always around 0.4 which is very bad for the significance of the original profits. In the original paper the p-values are around 0.01 always. But the time series they used was from 1920 to 1990 or so. The high p-value might thus be the result of volatile markets. The second thing concerns the returns of the RSI indicator. The return series is completely nonsense (look at the return plots) It looks like there is not a single day without a gain. I just don't see where the mistake is. Funny, if I use a GARCH(1,1) model as null model the RSI indicator works fine...a riddle to me I know that this might be a long question and reading through this takes a while, but I think bootstrapping with this random walk (i.e.) as a null model is pretty enlightening as it tells you whether actual profits are just the result of chance or if the profits are statistically significant.


title: Quantiles for non-normal cdf

text: When I compute the five-number summary on my sample, I obtain quantiles that differ from the quantiles I got from the empirical cdf, since they are not normally distributed data. Can you help me in the interpretation of this difference? For instance, with a randomly-generated Poisson dataset x               x=rpois(50, 2)     summary(x)     Min. 1st Qu.  Median    Mean 3rd Qu.    Max.     0.00    1.00    1.50    1.82    2.75    6.00      y=ecdf(x)     summary(y)     Empirical CDF:    7 unique values with summary     Min. 1st Qu.  Median    Mean 3rd Qu.    Max.     0.0     1.5     3.0     3.0     4.5     6.0      What does it mean that the 3rd quantile of the sample is `2.75` while it is `4.5` for the ecdf?


title: Box-Cox transformation for mixed models

text: Does there exists a Box-Cox method for linear Gaussian models with random effects ?


title: One sentence explanation of the AIC for non-technical types

text: I need a one-sentence explanation of the use of the AIC in model-building. So far I have "Simply put, the AIC is a relative measure of the amount of observed variation accounted for by different models and allows correction for the complexity of the model." Any advice much appreciated. R


title: Unknown number of colours Bernoulli Urn

text: Okay, so, in the traditional Bernoulli Urn problem, we have an urn with a number N, possibly infinite, of coloured balls, and there are k possible colours. That one I grok. However, what if I don't actually _know_ what k is? That is, what if I have an urn with N balls and an unknown but finite and strictly positive number of possible colours? The main question is, in fact, what my priors should be. What's the prior that there is exactly one colour? Exactly two? At least two? How do I update on the relative frequencies of each colour? Is this problem even solvable? My first lines of thinking are to have a vector of parameters $\vec \theta \in \mathbb R^\infty$ such that the first parameter is the number of colours in the urn (let's call it $\alpha$) and the remaining are the relative frequencies of each colour. If $P(A=n|\vec\theta)$ is the probability that the colour of the next draw will be n given the knowledge contained by $\vec\theta$, we'd have:   * $\vec\theta = (\alpha, p_1, p_2, p_3, ...)$   * $\alpha \in \mathbb N^*$   * $\left(\sum\limits_{n=1}^\infty P(\alpha = n) \right)= 1$   * $\left(\sum\limits_{n=1}^\infty p_n\right) = 1$   * $\forall n > \alpha : p_n = 0$   * $\forall n \in \mathbb N^* : P(A=n|\vec\theta) = p_n$ However, this is just wild speculation on my part. I'm mostly curious about whether this is even in principle solvable. What I'd want to know is a way to compute both the prior (objective/uninformative) and posterior distributions of $P(\vec\theta)$ or, in other words, the pdfs $P(\alpha)$, $P(p_1)$, $P(p_2)$, etc. How to start with them and how to update on them.


title: Logistic regression - cutoff

text: If you are running logistic regression on data, and the prior probability of the event is 50%, and there is no preference for sensitivity or specificity... When would you still want to use a cutoff? Can you please give an example to help me, a beginner, understand?


title: Statistical significance of conditional probabilties

text: I used conditional probability to find the probability that                         A    not A         B    27        6     not B    12                    C    not C         B    29        4     not B    50                    C    not C         A    36        3     not A    43      \begin{align} P(A|B) &=.81 \\ ~\\ P(B|A) &=.69 \\ ~\\ P(B|C) &=.82 \\ ~\\ P(C|B) &=.36 \\ ~\\ P(A|C) &=.92 \\ ~\\ P(C|A) &=.45 \end{align} These numbers seem too high to be not statistically significant. Is there any inference I can draw from this information that would be judged to be significant? These numbers are drawn from the interaction of a ligand with a protein. Each number in my tables represents a unique contact between the instance and the protein. There are non unique duplicates of the connections sometimes as many as 23 per atom. I do not consider them because they are confounding. B is the standard for which 33 atoms are interact with the protein. B is a substandard and has 39 contact points. A is another standard to a lesser degree it is sometimes rejected. C made from A to excrete A and C is always rejected. With these probabilities I plan to identify novel ligands which significantly overlap with B and A and the like rejecting those who are similar to C.


title: The X quantile of the difference between two RV's

text: To first provide some context, I have two boxplots (displaying median, 25th, and 75th percentiles), and I'm wondering what a boxplot of the difference between these two boxplots would look like. The random variables (RV's) aren't necessarily distributed normally. Specifically, the "random variables" are samples from posterior distributions, but I think my question is ultimately pretty generic. Since I have many samples from these 2 RV's, one way to do this would be to compute the differences between randomized pairings, then compute the quantile of interest from these differences. Because my statistical terminology might be flawed, let me say this in the R language:               X <- rnorm(10, 1000) # The first RV, a numeric vector of length 1000     Y <- rnorm(20, 1000) # The second RV     Qs <- quantile(X-Y, probs=c(0.25, 0.5, 0.75) #The quantiles of their difference      However, I was wondering if there could be a way to compute the quantile of the difference between these RV's without going back to the original data. For example, from some combination of summary statistics (means, medians, variances, other quantiles, etc) computed from X and Y. An example of something that does not give the correct answer, but which might shed light onto what I'm looking for:               wrong_Qs <- quantile(X, probs=c(0.25, 0.5, 0.75)) - quantile(Y, probs=c(0.25, 0.5, 0.75))      I am looking for something akin to the formulas for the propagation of uncertainty (variance) , but I haven't found anything yet. The more I think about this, the more I'm thinking that there isn't a way to do what I'm asking. Any thoughts?


title: Is binary logistic regression the right choice?

text: Apologies for the rudimentary question. I'm taking on a project at work that's a bit out of my wheelhouse and I want to bounce my ideas off of those more experienced than myself. We use Salesforce.com at the software company where I work, and I want to identify which lead behaviors (whitepaper downloads, demo views, webinar attendances, etc.) are predictive of those leads turning into qualified sales opportunities. The idea is that we can use this data to create a model, on which we'll base a scoring model going forward. I've identified binary logistic regression, using stepwise selection, as the best choice, based on my research. Essentially, my thinking is that the dependent variable (opportunity status) is binary (Opportunity = 0, Not an Opportunity = 1), which would indicate that logistic regression would be the best approach. Also, I'm not sure which behaviors and data points will ultimately be predictive of the lead becoming an opportunity, so stepwise selection seems like a good approach. Can anyone think of a more appropriate analysis technique, or am I on the right track?


title: How to cope with missing data in logistic regression?

text: I'm investigating optimal bidding in auctions, and am using logistic regression to predict the probability of winning an auction given a set of explanatory variables (e.g. the price I bid, number of competing bids etc). One explanatory variable I want to use is the _second_ highest price that was paid. However, by the design of the auction, I only observe the second highest price paid when I am the highest bidder (i.e. when I win the auction). This missing data is a major issue as my dataset indicates that there is a winning bid only ~20% of the time, hence I don't know the second highest price paid 80% of the time. Yet intuitively, I don't want to drop this variable as it seems to me knowledge of the second highest bid is extremely valuable in determining my chances of being the winning bid. Thus are there any standard methods to cope with this kind of missing data for logistic regression?


title: sum and product rules of probability

text: I am reading Bishop's Pattern Recognition and Machine Learning. In page 73, chapter 2.1. I can't understand the formula 2.19 : $$p(x=1|\mathcal{D})=\int_0^1 p(x=1|\mu)p(\mu|\mathcal{D})\text{d}\mu $$ The author say, this is obtained by sum and product rules. The sum rule is: $$p(X) = \sum_Y p(X,Y)$$ and the product rule is: $$p(X,Y)=p(Y|X)p(X)$$ But from this, I can't deduce the formula. Could you help me please.. thanks very much.


title: Factor analysis for categorical target variable

text: I'm doing some research into factor analysis and I've hit that barrier where I don't know what search terms to use. I'm trying to see if something is possible. Basically I have a data set with ~100 numerical variables. Each record also has a target variable which is a category. My ultimate goal is to perform knn on this data set. To reduce the dimensionality I ran PCA which gave a minor performance boost. My understanding is that PCA reduces the data set to variables which capture the most variance in the multi-dimensional data set. But this method doesn't even consider the category that each record has. Is there any method to reduce the number of dimensions, or identify key factors that specifically contribute to the target variable? This may be poorly worded (or just an outright stupid question...), so let me know if it doesn't make sense and I'll try to clear it up. Thanks!


title: Partial derivative of bivariate normal cdf and pdf

text: How do I take the partial derivative of bivariate normal cdf and bivariate normal pdf with its arguments (i.e. $x_{1}$ ,$x_{2}$ , and $\rho$ in the following equations)? \begin{equation} y=\Phi(x_{1},x_{2},\rho) \end{equation} \begin{equation} z=\phi(x_{1},x_{2},\rho) \end{equation} where $x_{1}$ is normally distributed with mean 0 and variance 1 and $x_{2}$ is normally distributed with mean 0 and variance 1. $\rho$ is the correlation between $x_{1}$ and $x_{2}$.


title: Can we estimate value of missing correlation in a given correlation matrix?

text: This question has no practical importance for me, but anyway: I am interested whether we can predict somehow one missing value in a correlation matrix. Obviously, classical regression modeling is not appropriate here. Is there any other way, how we can say something about missing correlation? In addition, does the answer depend on the type of correlation (Pearson, Spearman, Kendall, MIC)? Here is an example (real data, Pearson correlation). ![enter image description here](http://i.stack.imgur.com/6I7Cu.jpg)               matrix <- structure(c(1, NA, 0.126, -0.163, 0.203, 0.172, -0.101, 0.068,      -0.125, 0.066, 0.195, NA, 1, 0.019, -0.169, 0.211, 0.265, -0.023,      0.088, 0.016, 0.14, 0.275, 0.126, 0.019, 1, 0.218, 0.094, -0.053,      0.189, 0.179, -0.07, -0.021, 0.172, -0.163, -0.169, 0.218, 1,      -0.066, -0.182, 0.302, 0.193, -0.116, -0.08, -0.151, 0.203, 0.211,      0.094, -0.066, 1, 0.251, 0.064, 0.122, 0.093, 0.214, 0.378, 0.172,      0.265, -0.053, -0.182, 0.251, 1, -0.101, 0.02, 0.192, 0.353,      0.304, -0.101, -0.023, 0.189, 0.302, 0.064, -0.101, 1, 0.652,      0.126, 0.023, -0.016, 0.068, 0.088, 0.179, 0.193, 0.122, 0.02,      0.652, 1, 0.127, 0.01, 0.083, -0.125, 0.016, -0.07, -0.116, 0.093,      0.192, 0.126, 0.127, 1, 0.295, 0.023, 0.066, 0.14, -0.021, -0.08,      0.214, 0.353, 0.023, 0.01, 0.295, 1, 0.299, 0.195, 0.275, 0.172,      -0.151, 0.378, 0.304, -0.016, 0.083, 0.023, 0.299, 1), .Dim = c(11L,      11L), .Dimnames = list(NULL, NULL))


title: What is fitted in a GARCH: residual or log-return?

text: Given a time-series of log-return of SP500, then to obtain the volatility process what should we do? Some people say that we need using the ARMA model to withdraw the residual series, then plug this residual series into the GARCH model to obtain the conditional variance process? Or directly plug the log- return plug the log-return process of SP500 into the GARCH model to obtain the conditional variance? I saw that in the book Introduction to Time Series in R, the author fits simulated residuals in GARCH model then he fits SP500 log returns in the GARCH model. I got confused...


title: Can explanatory variables be non-Gaussian?

text: I am conducting a logistic regression. However, the Shapiro-Wilk test has determined that one of the X variables is non-Gaussian. Do the explanatory (X) variables need to be Gaussian? If so what are the usual transformations that can be done?


title: Feature selection using correlation

text: I am trying to do some feature selection using correlation. However, I found that my features are not that correlated. The highest correlation was 0.08. So I am not sure if this is a useful thing to do. I tried to use lasso however, I have huge number of features and when I tried to use lasso, it took like ages. The same with other feature selection methods like wrapper based methods. What should I do under such circumstances?


title: Confusion related to multitarget regression

text: I have a bunch of features/xs and a bunch of targets ys. I want to run a linear regression on this. I don't know if I should take every y to be independent and learn separate regression models or just a single model with multitargets and xs. What would be the best? If I have to do with multiple ys how can I do it?


title: How to test relationship between categorical and numerical variables

text: I have a list of sites and a list of survival probabilities associated with those sites. The data looks like this:               Site     Probability     A        0.8     B        0.4     C        0.2     ...     (14 sites in total)      Is there any way to test if survival probability varies significantly between sites? What I'm finding difficult is that the factor (sites) has so many categories.


title: How to find/estimate probability density function from density function in R

text: Suppose that I have a variable like `X` with unknown distribution. In Mathematica, by using `SmoothKernelDensity` function we can have an estimated density function.This estimated density function can be used alongside with `PDF` function to calculate probability density function of a value like `X` in the form of `PDF[density,X]` assuming that "density" is the result of `SmoothKernelDensity`. It would be good if there is such feature in R.This is how it works in Mathematica http://reference.wolfram.com/mathematica/ref/SmoothKernelDistribution.html As an example (based on Mathematica functions):               data = RandomVariate[NormalDistribution[], 100]; #generates 100 values from N(0,1)          density= SmoothKernelDistribution[data]; #estimated density          PDF[density, 2.345] returns 0.0588784       Here you can find more information about PDF: http://reference.wolfram.com/mathematica/ref/PDF.html I know that I can plot its density function using `density(X)` in R and by using `ecdf(X)` I can obtain its empirical cumulative distribution function.Is it possible to do the same thing in R based on what I described about Mathematica? Any help and idea is appreciated.


title: Multiple imputation to part of a dataset

text: I am working with a survey dataset which contains hundreds of variables. The item-missing data rate ranges from 0.2% to 10%. In order to retain study units with missing values and to maintain a reasonable statistical power for my analyses, I attempted to do multiple imputation using sequential regression (chained equations) in Stata. However, with all the several variables with missing items in the equations, the models could never attain convergence. I spent days re-specifying the equations and the options within the equations but all in vain. In the end, I simplified the equations by selecting only variables with high item-missing data rates and variables which apply to all study units in the dataset. This time I succeeded in multiply imputing the data. However, I have two questions.   1. Since I did multiple imputation using only some of the variables in the dataset, there are still variables with missing values (though small). Is it acceptable to still have variables with missing values while analyzing multiply-imputed data?   2. I am also planning to do sub-population analysis. However, due to the convergence problem I stated earlier, variables pertaining only to sub-populations are excluded from the multiple imputation I did. Is it acceptable to do a separate multiple imputation to the sub-population of interest?


title: Why does using residuals as predictor not alter coefficient?

text: I am working with a data set describing a set of US Senators, their votes on a climate bill, and several characteristics of the state they represent. One of the variables, 'imptot', describes the combined importance for their state's economy of **four industries**. These industries cause greenhouse gas emissions and would be strongly affected by climate legislation. I also have an indicator of the **general greenhouse gas intensity** ('GHGintensity') of each state's economy, meaning its greenhouse gas emissions divided by its GDP. These two variables are positively correlated. I would like to know what's the influence of greenhouse gas intensity on a Senator's vote, **after controlling** for the total importance of these four industries. In other words, I would like to find out how much variance in voting behavior can be explained exclusively by greenhouse gas intensity, not taking into account the variance that can be explained by the presence of these four specific industries. I reasoned that if I regressed GHGintensity (dep. var.) on imptot (indep. var.), the **residuals** of this regression (a new variable called 'GHGresiduals') would represent the variance in GHGintensity that could not be explained by imptot. However, when I then ran regressions of these variables on voting behavior ('vote'), the following outcomes resulted for the fixed effects portion of the (multilevel logit) model:                              Estimate Std. Error z value Pr(>|z|)     (Intercept)      -0.045      0.607  -0.074    0.941     imptot            0.000      0.022   0.003    0.997     impbentot         0.085      0.530   0.160    0.873     GHGintensity     -0.176      0.076  -2.303    0.021     member            5.399      1.152   4.689    0.000     pubop             1.991      0.233   8.535    0.000                         Estimate Std. Error z value Pr(>|z|)     (Intercept)      -0.839      0.387  -2.169    0.030     imptot           -0.034      0.021  -1.622    0.105     impbentot         0.085      0.530   0.160    0.873     GHGresiduals     -0.176      0.076  -2.303    0.021     member            5.399      1.152   4.689    0.000     pubop             1.991      0.233   8.535    0.000      As you can see, GHGintensity in the first specification had the same coefficient as GHGresiduals in the latter. The coefficients of imptot, however, were different. Because GHGintensity and GHGresiduals are different variables, I did not expect this result. I had expected the effect of greenhouse gas intensity to change, because my intuition was that I had now removed the component of this variable that was correlated with imptot. Was this intuition incorrect? Should I do a different analysis to answer the questions I am asking? How should I interpret the results I obtained, and why did coefficients not change the way I expected? EDIT: I am unfortunately not in a position to share the data, but here's the associated R code:               dataSenate$GHGresiduals <- residuals(lm(GHGintensity~imptot,data=dataSenate))          Smodelgeo <- formula(vote ~       imptot +       impbentot +        GHGresiduals  +       member +       pubop +       (1 | state) + (1 | yearbill))          Sfitgeo <- glmer(as.formula(Smodelgeo), dataSenate, family=binomial,na.action = na.omit)


title: Autoregressive model with input variables in proc arima procedure

text: I am currently working on the time series analysis for series Y but I have to use other two variable A and B as an input variable in SAS proc arima procedure. But I am unable to interpret the cross correlation chart and backshift parameters. BackShift 1 for variable A and back shift 2 for variable B. Could anyone help me how to interpret the input variables outputs.


title: Cutoff and precision values of a binary classifier

text: Let's say I have fitted a binary classifier to some data and I'm varying the cutoff value, effectively producing a ROC-curve. Knowing the true proportions of positives and negatives, I can calculate the proportion of data classified as positive, i.e. $TP + FP$ for any given cutoff, and also its corresponding precision value. Am I correct in interpreting these precision values in relation to $TP + FP$ in the following way? Suppose there is a population of size $N$ that is assumed to be identically distributed to the fitted data. Suppose also that I have a constrain of "selecting" (i.e. classifying as positive) exactly $n$ individuals from this population. Would the corresponding precision value of this constrained classification be the very same as that corresponding to the cutoff for which $\frac{TP + FP}{TP+FN+FP+TN} = \frac{n}{N}$? In other words, can I interpret varying the cutoff value in a ROC-curve as making the best classification possible working under some constraint of classifying exactly $n$ individuals as positive? For instance, if I have to classify every individual as positive, this makes $\frac{n}{N} = 1$, which corresponds to the cutoff value $0$. If I have to classify half of the population as positive, this too corresponds to some cutoff value, precision etc. It makes sense to me but I can't find this kind of interpretation of a ROC- curve in any literature, which makes me a bit hesitant.


title: How can I do a test for proportions when the data has been collected weekly from the same general population?

text: So I have a series of proportions taken weekly for 7 weeks. This is a study of which way spiders were facing. So for example, the first week, 20 spiders were facing out, and 10 were facing in. Now I want to test whether the proportion of spiders facing out was consistent across the 7 weeks. The population of spiders sampled were the same every week, and so they are not really independent. I know you can do a binomial test for say one set of proportions, but is there anyway to test across the whole period? I am sort of familiar with R. Thanks a ton!


title: Determine if many coins are fair

text: Imagine I have many coins ($>100$), and for each coin I have the data: $(n_{\rm heads}, n_{\rm tails})$ where $\sum(n_{\rm heads}, n_{\rm tails})$ is decreasing for each next coin. Any one of the coins can be unfair ($P({\rm head})\ne 0.5$). Some coins might be an exact copy of any other coin, some might not be. How do I determine all coins are fair?


title: Missing values in SPSS for Chi-Square Tests

text: I'm a postgraduate student, doing a accent perception study, this is my first time using SPSS! I'm conducting the above tests on my data, where my participants used Likert scales to analyse the personality traits of 2 guises produced by speakers from Devon. Because of this, some of my responses are 'unknown' or sometimes the participants have put 'don't know' as a response. My question is, do responses such as these count as missing values in my data? Also, is there a specific command in SPSS to ensure that you 'exclude cases pairwise' or is it an automatic thing? Also, In order to set up my data set on SPSS, I have ordered my data in Excel, before I copy it across to SPSS. For each dependent variable, e.g. 'perceived class' I have 4 speakers producing 2 separate guises and I thought the easiest way to ensure that my scores were above 5 for the chi-square for independence test would be to simply add all the scores per participants for each guise, for example, all 4 of participant 1's 'perceived class' scores for the 2 guises etc... The Psych textbook I'm using suggests a score of 99 or 999 for missing values, so if I put 99 for a missing value in excel and I am adding those scores together, is it going to affect my scores in Excel?


title: plotting hurdle regression output

text: I am attempting to plot hurdle regression output with an interaction term in R, but I have some trouble doing so with the count portion of the model.               Model <- hurdle(Suicide. ~ Age + gender + bullying*support, dist = "negbin", link = "logit")      Since I am unaware of any packages that would allow me to plot the hurdle model in its entirety, I estimated each component separately (binomial logit and negative binomial count), using the `MASS` package, and I am attempting to plot the results. So far, I have had the best luck using the `visreg` package for plotting, but I would like to hear other suggestions. I have been able to reproduce and successfully plot the original logistic output from the hurdle model in `MASS`, but not the negative binomial count data (i.e., the parameter estimates from `MASS` are not the same as they are in the hurdle regression output). I would greatly appreciate any insight regarding how others have plotted hurdle regression results in the past, or how I might be able to reproduce negative binomial coefficients originally obtained from the hurdle model using `glm.nb()` in `MASS`. Here is what I am using to plot my data:               ##Logistic           logistic<-glm(SuicideBinary ~ Age + gender + bullying*support, data = sx, family="binomial")          data("sx", package = "MASS")     ##Linear scale     visreg(logistic, "bullying", by="support", xlab = "bullying", ylab = "Log odds (Suicide yes/no)")          ##Logistic/probability scale     visreg(logistic, "bullying", by="support", scale = "response", xlab = "bullying", ylab = "P(Initial Attempt)")          ##Count model          NegBin<-glm.nb(Suicide. ~ Age + gender + bullying*support, data = sx)          data("sx", package = "MASS")     ##Linear scale     visreg(NegBin, "bullying", by="support", xlab = "bullying", ylab = "Count model (number of suicide attempts)")          ##Logistic/probability scale     visreg(NegBin, "bullying", by="support", scale = "response", xlab = "bullying", ylab = "P(Subsequent Attempts)")


title: Looking for an adaptMCMC (R package for adaptive MCMC sampling) tutorial

text: I need to fit a very complex Bayesian hierarchical model and I want to use the adaptMCMC package of Andreas Scheidegger to do the posterior sampling. I can only find the adaptMCMC package reference manual, which doesn't provide a very helpful example of the core sampling function. Does anyone know if there is a more helpful tutorial out there? Thanks.


title: How to test for wide-sense stationarity with only one sample path of the process?

text: I have a univariate time series consisting of 70,000 observations (power consumption of a building) over equal time increments (15 minutes). How do I check whether this realization is wide-sense stationary? Note that I don't know the formula of the underlying stochastic process, so I cannot calculate the mean function E[x(t)] or the autocovariance function E[(x(t1)-m(t1))(x(t2)-m(t2))] (I'm aware of the definition of wide-sense/weak/covariance stationarity and I've read several threads on the topic here at Cross Validated, the closest one being this one.) I'm arriving to the conclusion that checking a sample path for weak stationarity comes down to some heuristic tests. The Engineering Statistics Handbook (which I like for its practical approach) states: > Stationarity can be defined in precise mathematical terms, but for our > purpose we mean a flat looking series, without trend, constant variance over > time, a constant autocorrelation structure over time and no periodic > fluctuations. (Also see the reply in this MATLAB Central thread.) If you agree with this practical approach, then my comments/questions are:   * We can easily do linear fitting to tell if there's a trend. We're good here.   * But how do we test for a constant autocorrelation structure over time?   * And how do we check for constant variance over time? Do we using moving windows, or non-overlapping windows? And what should the size of the window be? I'm looking for some good rules of thumb to test with. Thanks in advance for the help, and apologies for the long post. (First post here.)


title: Strange results in Holt forecast

text: I am trying to understand what could be causing these strange values to appear on applying a Holt model to a vector. The data represents actual sales of an item.               library(forecast)          sdata<-c(1955651,1691857,1617358,1509591,1452025,1387790,1630114,3522984,3322908,2315689,     4044634,11553435,9747305,5289976,3905946,1646080,1356731,1248135,1150756,848804,847815,     940673,1129751)          > holt(sdata)     Point Forecast         Lo 80         Hi 80          Lo 95         Hi 95     24  6.144027e+122 6.144014e+122 6.144041e+122  6.144007e+122 6.144048e+122     25  6.144033e+122 4.502199e+122 7.785866e+122  3.633064e+122 8.655001e+122     26  6.144038e+122 3.822133e+122 8.465942e+122  2.592991e+122 9.695084e+122     27  6.144043e+122 3.300300e+122 8.987785e+122  1.794914e+122 1.049317e+123      This involves an exercise where a very large number of forecasts has to be generated. This is merely an attempt to find errors before the forecast is run across all the data. A number of records generated such strange records. I am sure there is a good explanation behind why we see such data. Note that the data was coerced to a ts object prior to the model being applied on the same (and produced similar results) Thanks in advance.


title: Non-parametric version of paired t-test (Mann–Whitney U test)

text: I have data from 100 individuals doing a 2AFC under a control and a manipulated setting. I want to compare the means. I initially used a paired t-test. My supervisor noticed the distributions were not normal and recommended me using the Mann–Whitney U test instead of the t-test. However, as far as I can tell, the Mann–Whitney U test does not allow for paired data. Or does it? What non-parametric test should I use? What test generalizes the paired t-test? In other words, I have a non-normal distribution and I want to test how different its mean is from zero.


title: R nnet package Softmax/Entropy option

text: Can some one tell me when and why to use `softmax=true` option and `entropy=T` option? The model i am trying to build has binary outcome and `softmax=true` works only when the two outcomes are classified as two variables on its own. Otherwise `entropy=T` works. I could not get much help on this in the documentation. The question has been asked in many forums and no reply yet.


title: Uncertainty from Box–Cox estimation

text: Consider the model $$ Y_i^{(\lambda)} = \alpha+\beta x_i + \varepsilon_i,\qquad \varepsilon_i,\ (i=1,\ldots,n) \sim \mathrm{i.i.d.}\ N(0,\sigma^2) $$ where $$ \begin{align} y_i^{(\lambda)} & = \begin{cases} \dfrac{y_i^\lambda-1}{\lambda(\operatorname{GM}(y))^{\lambda -1}} , &\text{if } \lambda \neq 0 \\[12pt] \operatorname{GM}(y)\log{y_i} , &\text{if } \lambda = 0 \end{cases} \\[12pt] \text{and} & \operatorname{GM}(y) = (y_1\cdots y_n)^{1/n}\text{ is the geometric mean.} \end{align} $$ As my estimate of $\lambda$, I use the value that minimizes the sum of square of residuals when $\alpha$ and $\beta$ have likewise been estimated by least squares. **In finding a confidence region for $\alpha$ and $\beta$, how should one take into account the uncertainty in the estimate of the Box–Cox parameter $\lambda$?**


title: Priors for log-normal models

text: I am trying to determine what the most appropriate non-informative priors are for the two parameters of a log-normal distribution (for an accelerated failure time model). I had been working with a normal prior on the mean of the logged data (or log-normal prior on the median), and an independent normal prior on the log of the sd of logged data (or log-log normal prior on the dispersion). It seems to me that the parameterization of the second parameter is quite informative on the natural scale of the parameter but I am pretty new to this. Any ideas on better ways to approach this?


title: Using AUC to compare logistic lasso and elastic net

text: I've seen this question answered here but I do not understand the answer. Harrell recommends using deviance based measures. David Hand (referenced in the thread) says that that the AUC is inappropriate to compare different models, because it uses different misclassifications costs. I don't see how this would be the case when comparing lasso and elastic net, given that they are both trained on the same set of predictors.


title: Why is expectation the same as the arithmetic mean?

text: I am a beginner in statistics and probability study. Today I came across a new topic called the Mathematical Expectation. The book I am following says expectation is the arithmetic mean of any random variable coming from any probability distribution. But it defines expectation as the product of some data and the probability of it. How can these two be same? How can the probability times the data be the average of whole distribution?? Thanks in advance..


title: Univariate clustering of time series

text: I just want to know if its possible to cluster an univariate time series, in order , say, to detect anomalies? and do you have any online version for denstream code, in Matlab? * * * here is the time series > > +0003;+0004;+0004;+0005;+0004;+0004;+0004;+0003;+0003;+0004;+0003;+0002;+0004;+0002;+0002;+0002;+0002;+0002;+0002;+0005;+0004;+0004;+0004;+0004;+0004;+0004;+0004;+0004;+0004;+0004;+0004;+0005;+0004;+0004;+0004;+0004;+0003;+0003;+0000;+0002;+0004;+0002;+0004;+0004;+0004;+0004;+0003;+0003;+0004;+0003;+0002;+0000;+0003;+0002;+0003;+0002;+0002;+0002;+0002;+0002;+0003;+0002;+0002;+0003;+0000;+0003;+0003;+0002;+0003;+0002;+0003;+0002;+0004;+0004;+0004;+0003;+0004;+0004;+0004;+0004;+0004;+0005;+0004;+0002;+0003;+0002;+0005;+0004;+0003;+0003;+0001;+0005;+0003;+0003;+0005;+0005;+0004;+0005;+0002;+0005;+0005;+0006;+0005;+0005;+0004;+0003;+0003;+0005;+0002;+0002;+0001;+0003;+0004;+0002;+0002;+0002;+0001;+0001;+0003;+0004;+0004;+0004;+0005;+0004;+0005;+0002;+0002;+0005;+0003;+0002;+0003;+0002;+0001;+0002;+0003;+0005;+0004;+0005;+0004;+0004;+0003;+0003;+0004;+0005;+0003;+0004;+0005;+0004;+0005;+0003;+0003;+0005;+0005;+0005;+0003;+0002;+0004;+0004;+0004;+0005;+0005;+0004;+0005;+0005;+0004;+0004;+0003;+0003;+0005;+0005;+0005;+0005;+0004 * * * First you should now @IrishStat, that i should use a datamining technique for clustering, an exemple of times serie dataset is presented in this : ![picture](http://i.stack.imgur.com/xhK1u.jpg) What i'm using now , is a density base clustering, let's said that I have adapte the Dbscan for streaming Time series points, and its Should detect anomalies, and changes eacht time, after your introdudction , i'll trying to explain to u by pictures what i'm doing, thanks u , time data series sample :                    [,1]     V1    3.0     V2    4.0     V3    4.5     V4    5.0     V5    4.0     V6    4.0     V7    4.4     V8    3.0     V9    3.0     V10   4.0     V11   3.0     V12   2.0     V13   4.0     V14   2.0     V15   2.0     V16   2.0     V17   2.0     V18   2.0     V19   2.0     ..... so on


title: Ways to calculate the moments of a function of sample moments?

text: 1. Are there ways to analytically derive the moments of a function of sample moments? For example, my recent question here hasn't been addressed satisfactorily yet here: link   2. Do the moments of a function of sample moments depend on the sample distribution only through the moments of the sample distribution? This is the case in the linked question. Thanks.


title: Choose best time window to count events in order to produce an indicator

text: We want to create indicators for event based clinical conditions, like migraine or epilepsy. This conditions are characterized by events which can happen with various frequencies and we would like to classify patients according to the frequency of events. The problem is that the frequency of the events is connected with the length of the observation period. That is short periods produce high variability of results, (eg one week you can have 2 seizures, the week after none). Long periods instead are less feasible. So I was wondering whether there is a way to estabilish a proper time window, based on data and literature. PS: I know it's pretty vague now, but I'm sure your initial answers and comments will help me better define the problem in my mind.


title: What is the null hypothesis in the definition of Sensitivity and Specificity?

text: I have hard time to understand that what is the null hypothesis in the definition of Sensitivity and Specificity here? I'm student of Statistics and from my little knowledge in this field I can see that the here null hypothesis should be: $\textrm{H}_{0}:\textrm{ Person is healthy}$ Is this true? Thanks in advance for your help.


title: prediction of polls

text: Just as an example Scotland has poll to decide whether they need to be independent from UK or not. Here is BBC's summary of different polls: http://www.bbc.com/news/events/scotland-decides/poll-tracker ![enter image description here](http://i.stack.imgur.com/QUbq1.jpg) ![enter image description here](http://i.stack.imgur.com/m2Xjo.jpg) It looks like as this progresses the undecided fraction is reducing and yes component is increasing. Provided time and information we have is there a simple way we can predict the result with certainty ?


title: How to calculate 95% confidence interval of two group means and perform hypothesis test comparing group means?

text: Below is the information about the weight of coffee (in thousand kg) sold at two locations (location1 and location 2)               26.3, 26.5, 16.2, 16.2, 26.4, 16.8, 16.1, 16.9 (8 shops)     44.4, 44.7, 34.9, 54.1, 44.8, 22.6, 36.4 (7 shops)         1. Find the 95% confidence interval for the difference between the mean weight for the coffee sold at different locations.   2. where do you prefer to own a shop? Why? Explain   3. Now you need to prove whether preferred area has really got higher mean of weight of coffee sold over the week in question. carry out appropriate tests at 5% and 1% levels.


title: Should you use a regression line predicting drop from new car sale price from car age to set the price of a used car?

text: A regression showed a strong linear relationship between a car's age (measured in months) and a drop in its price from new. Your friend sees these figures and thinks he can sell his car at a price obtained from this regression analysis. What would you say to him?


title: How to evaluate the integral e^{(c ln(1/x))^s} dx?

text: Can anyone help me to evaluate $\int_{\alpha}^1 e^{-(c \cdot ln(1/x))^s} dx$, where $0 \leq \alpha \leq 1$ and $s$ is real. I tried changing variables, integration by parts etc. and got no clue how to handle this.


title: How to perform paired sample t-test when the standard error of the difference is zero?

text: I ran a paired T-test to determine the mean difference between pre and post intervention of one of the variables in my study. However the test was not run as the SPSS output showed "the correlation and t cannot be computed because the standard error of the difference is 0". What should I do to get the p-value of the mean difference for that variable?


title: What to say to a client that thinks confidence intervals are too wide to be useful?

text: Suppose I am a consultant and I want to explain to my client the usefulness of confidence interval. The client says to me that my intervals are too wide to be useful and he would prefer to use ones half as wide. How should I respond?


title: The formula for finding the mode of grouped data

text: Here is a link that describes the formula to find the mode of grouped data. Here is a link that gives a graphical method to finding the mode of grouped data. Question: Can someone please explain how the formula corresponds with the graphical method? The formula is an interpolation but I am not able to see how it expresses the idea of the graphical method. If it helps, I understand the formula for calculating the mean of grouped data. Thank you.


title: How can Gröbner bases used to describe discrete probability?

text: [Cross post from here, figured this community may be more relevant] I am working in the field of machine learning, and I have come across a few papers that show relationships between Gröbner bases and discrete probability. So I come here for help. Can you please explain how can Gröbner bases used to describe discrete probability? I have looked at Gröbner bases and I understand the general concepts (and used Maple to calculate a few examples). So it is the link that is missing for me.


title: Using lmer for prediction

text: Hello I have two problems that sound like natural candidates for multilevel/mixed models, which I have never used. The simpler, and one that I hope to try as an introduction, is as follows: The data looks like many rows of the form `x y innergroup outergroup` where x is a numeric covariate upon which I want to regress y (another numeric variable), each y belongs to an innergroup, and each innergroup is nested in an outergroup (i.e, all the y in a given innergroup belong to the same outergroup). Unfortunately, innergroup has a lot of levels (many thousands), and each level has relatively few observations of y, so I thought this sort of model might be appropriate. My questions are   1. How do I write this sort of multilevel formula?   2. Once **lmer** fits the model, how does one go about predicting from it? I have fit some simpler toy examples, but have not found a predict() function. Most people seem more interested in inference than prediction with this sort of technique. I have several million rows, so the computations might be an issue, but I can always cut it down as appropriate.  I won't need to do the second for some time, but I might as well begin thinking about it and playing around with it. I have similar data as before, but without x, and y is now a binomial variable of the form $(n,n-k)$. y also exhibits a lot of overdispersion, even within innergroups. Most of the $n$ are no more than 2 or 3 (or less), so to derive estimates of the success rates of each $y_i$ I have been using the beta-binomial shrinkage estimator $(\alpha+k_i)/(\alpha+\beta+n_i)$, where $\alpha$ and $\beta$ are estimated by MLE for each innergroup separately. This is has been somewhat adequate, but data sparsity still plagues me, so I would like to use all the data available. From one perspective, this problem is easier since there is no covariate, but from the other the binomial nature makes it more difficult. Does anyone have any high (or low!) level guidance?


title: Test if Two Categorical Samples Have the Same Joint

text: It seems to me that the standard way to deal with the problem of testing wether two multivariate samples of categorical data have the same joint probability distribution is to employ a goodness-of-fit test, like the chi- squared one: assuming that the frequencies of one sample (normalised counts per combination of variables) constitute the "theoretical" distribution $E$ and the others the "observed" $O$, one can compute $\chi^2=\sum\limits_{i=1}^{N} \frac{(O_i - E_i)^2}{E_i}$. Then, given the p-value obtained from the chi-squared probability distribution, one is able to say if the null hypothesis that the two samples have the same joint distribution can be rejected. However, I was asking myself if there are other approaches to this aim. In particular, I would be interested in some method that were also able to highlight which categorical variables are the most "distant" in the two samples. Do you have some advice? Thanks!


title: Varied experimental conditions in repeated measure designs

text: I'm running an experiment with an e-learning system using a repeated measures with crossover design. Subjects are shown two videos on a web page. In one session, a feature (a table of contents) is enabled, while the other session is a control:               C=Control F=Feature     subject     pretest1   video 1   posttest1    pretest2   video 2    posttest2     S1                        C                                  F     S2                        C                                  F     S3                        C                                  F     S4                        F                                  C     S5                        F                                  C     S6                        F                                  C      The amount learned is measured by subtracting pre-test result from post-test result. My issue is that I can't swap around the order of the two videos for the crossover because Video 2 requires knowledge of Video 1. If I had chosen videos on separate topics I would be introducing variability in aptitude, interest and prior knowledge within the subjects. It is also unlikely that both videos have equal difficulty, so the learning gains in one can't be directly compared to the learning gains in the other. I know that the crossover will cancel out this effect and allows me to get a mean difference, but this seems to prevent me from doing basic EDA like scatterplots, and I'm not sure how to proceed with statistical analysis from there. Is there anything that I can do to improve this experiment design, and do you have any advice or resources on how I should proceed with analysis in this scenario?


title: 2-way interaction with only one main effect in the model - is it a "real" interaction?

text: I'm estimating a GLM with a continuous variable $X_1$, a dummy variable $X_2$ and the interaction term $X_1*X_2$. After my understanding, when interpreting the coefficient of $X_1$ you should assume the value of zero for $X_2$, i.e. also for the interaction term. Now a colleague suggested that if I skipped one of the main effects, then the interaction would not be a real interaction, i.e. I could simply interpret the coefficients of both - the remaining main effect and interaction term as if there was no interaction term in the model. I could not find any literature that supported this view. Is it correct?


title: Probabilistic Output from AVA (Multi-class) Classifier

text: What is the typical method for creating a probabilistic output from an All-vs- All or One-vs-All multi-class classifier? For example, if my problem has 3 classes (Class1, Class2, Class3) and I build an AVA classifier, then I have a classifier for:   * Class1 vs Class2 ($f_{12}(x)$)   * Class1 vs Class3 ($f_{13}(x)$)   * Class2 vs Class3 ($f_{23}(x)$) Then, the classification is simply: $$ f(x)=\text{argmax}_{i} \sum_j f_{ij}(x) $$ However, we have no _probability_ output of the class label (as with binary class logistic regression). In other words, how does one typically compute $p(y=Class1|x)$, $p(y=Class2|x)$, $p(y=Class3|x)$.


title: How to compare the positive predictive value and negative predictive value of two diagnostic tests?

text: Im a surgeon who's trying to compare two diagnostic tests that are used to diagnose appendicitis. Two diagnostic tests were applied on 150 patients and the results were compared to a gold standard. I understand that the sensitivities and specificities can be compared by using a McNemars chi square test. But what about the positive and negative predictive values... how do you compare these? Any help would be appreciated! eg. ![enter image description here](http://i.stack.imgur.com/2ZrYE.jpg)


title: How to choose an error metric when evaluating a classifier?

text: I've seen different error metrics used in the Kaggle competitions: RMS, mean- square, AUC, amongst others. What's the general rule of thumb on choosing an error metric, i.e. how do you know which error metric to use for a given problem? Are there any guidelines?


title: Need clarity on alpha, beta, gamma optimization in Triple Exponential Smoothing Forecast

text: I asked a variation of this question, but I want to be more direct. Take the exact same Triple Exponential Smoothing Model (Holt-Winters with a moving level, trend, and seasonal component)--- Would or would not --- the smoothing parameters alpha, beta, and gamma (determining the weight of recent data on level, trend, and seasonality, respectively) ---- be optimized differently for different forecast horizons? ... aka t+1 periods vs. t+2 periods, t+6 periods, t+12 periods, etc.? My initial hunch ... is yes. The nature of Holt-Winters is that to forecast 6 periods out (weeks, months, whatever) -- you take the current level, add 6 times the calculated trend, and multiply (or add) by the seasonality of the last corresponding period. It just seems like common sense that a 'recent data' chasing model would, in many cases, be more accurate for a horizon of t+1, and a poorer for a long term, t+6 model. BUT, I definitely can be wrong on that. I just don't see this discussed anywhere, and perhaps I have some gaps in knowledge on this parameter optimization process, and perhaps I'm overfitting models. When specialized software typically optimizes these values --- I assume against test data or what-have-you, minimizing MSE or MAPE ... does this software minimize the MAPE against t+1 forecasts, or t+6 forecasts, or both, or every horizon? This is what confuses me.


title: Linear Regression With Groups vs. Points - Issue with Influence of Zeros

text: I have a data set with two columns, the first of which is to be used as a response variable and the second of which is to be used as a predictor variable. The predictor variable, however, is populated by integer values, many of which are zero. When I place data points for the predictor and response variables into groups corresponding groups, no group has a value of zero anymore. Performing a linear regression on the data when it is in this form produces a significantly higher $R^2$ than in the ungrouped regression. My question is: Is something seriously wrong with my methodology if the change in $R^2$ is massive between the grouped and ungrouped models? Can this be explained by the large number of zeros? To provide a slightly better idea of what I mean by grouping consider the following data:               A  B          0  .2     0  .3     0  .3     2  .6     4  .7     6  .8     7  .7     4  .4      A would be grouped as follows:               0-3     4-7      with the new B values being the group averages.


title: What's the proper way to do automatic iterative outlier rejection?

text: I know this is a touchy subject, so I'll proceed with caution... I am building regression models for some hundreds of thousands of data files containing only two columns. The sensor collecting these data is not great and sometimes things interfere with the measurement, so the data are extremely noisy; a typical file might contain 300 inliers and 50 outliers. The way I'm proceeding currently is as follows:   1. Build regression model for the inliers of the data (initially all of the data)   2. Obtain residuals by subtracting the model from all of the data   3. Calculate the scale of the residuals (e.g. using STD or MAD)   4. Reject (sigma-clip) all outliers exceeding some sigma * scale   5. Repeat the process until the same outliers are rejected two iterations in a row Sometimes the initial predictor is bad, so some of the flagged outliers are really inliers. When I check for outliers in the second iteration, I have it examine all of the data including the previously rejected outliers, and I ``unflag'' outliers if they now agree with the updated model. The problem is that this sometimes leads to an infinite loop. If I build a model, reject a point, rebuild the model, unflag that point and reject a different point, rebuild the model, unflag that point and reject the first point again, etc., this process never ends. I could limit it so that it only iterates some predetermined number of times, but this seems like a hack to me. So what's the right way to do it?


title: Non-parametric regression models and data transformations

text: I am trying to explain datasets gathered from simulation with various regression models. Linear regression seems not applicable because my data after transformations still violate the assumptions for linearity, normality of errors etc. I have built several non-parametric models (splines, regression trees, SVR and Gaussian processes), however I am not quite sure if I should transform the data or not for the analysis. From the theory I have so far understood that non-parametric models make no / very few assumptions about the structure of underlying data, however it is not clear to me if this means that I can apply non-parametrics to my original data without transformations of any kind. What does the term "few" assumptions mean? Unfortunately I haven't managed to find something more concrete on this. Another question is that splines for example derive the model with the least squares method. How is it possible that in this case no assumptions on the data are necessary as opposed to linear regression?


title: online semi-supervised learning with incomplete input data

text: Let's assume that one has the following inputs and corresponding outputs:               x1 = (a,b,c) with corresponding output y1, which is a number.     x2 = (d,e,?) with unknown output     x3 = (?,?,g) with known output y2, which is a number      and so on. So, by "incomplete input data" I don't just mean that some input data is unlabeled but some data-points from certain data-vectors are missing and are sometimes labeled and sometimes unlabeled. Is there at least one known online semi-supervised learning algorithm for this kind of "input data with gaps"? So far I only found algorithms for online semi-supervised learning for "input data without any gaps".


title: One-way Chi-Square goodness of fit and Cramér's V

text: Would it make sense to apply a measure of association test like Cramér's V on a 1xk table? I am a bit puzzled because in the literature phi and Cramér V are normally used on 2x2 or bigger contingency tables, but while I was looking for some more info, I found the following statement on Wikipedia: "Cramér's V may also be applied to goodness of fit chi-squared models when there is a 1×k table (e.g.: r=1). In this case k is taken as the number of optional outcomes and it functions as a measure of tendency towards a single outcome." No reference is provided though. What do you think about it? Could you please give me some bibliographic reference?


title: permutation test with a distribution of correct values

text: I want to test my classifier performance using a permutation test. I know how to test one actual value against a null distribution of values generated by randomly reassigning labels to my two classes. However, in this case, I don't have one value but also a distribution of performances on correctly labelled data, because I perform the classification several times (to use slightly different, but overlapping, data and to split the data up in different ways for training/testing and cross-validation).   1. Is it enough to have one random relabelling for each classification that I do, so that I have an actual and a null distribution with 100 (or 1000) entries each? Or should I, for some reason, make a null distribution for each actual decoding performance?   2. How do I test whether my decoding is significantly better than the null distribution in this case? Thanks!


title: Non-linear least absolute deviation regression with multiple global minima

text: I am fitting a single exponential decay formula with three parameters (a,b,c): y ~ $a \exp(-xb) + c$ using the LAD cost function: $ \min \sum |(y - f(x))| $. $x$ is in units of time (as is $b$), and $a$, $c$ are unitless. Measurements are taken far enough out in time that the last $y$ is effectively only noise from measurement process. Sometimes we have data points that are corrupted by sporadic errors that are most definitely not normally distributed when they occur at all, so this method seems a very reasonable approach to use. I include an example using simulated data below: ![Example using simulated data](http://i.stack.imgur.com/KkYDT.png) Ignoring parameter $c$, I plot the two-parameter cost function graphically, and it is easy to see that there is a plane of points in the parameter space that minimizes the cost function (see figure). I have confirmed this for the three-parameter model by brute force searching the parameter space. ![Parameter space plot with minimia scattered about](http://i.stack.imgur.com/gocL4.png) My questions is: Is there any way, given the data, the model and assuming I can find all of the reasonable minima, to choose one of those optimal parameter sets as the representative parameter set? Right now I am choosing the centroid of the hyperplane (if that's an appropriate term) of the minimizing parameters, and that's getting me extremely close to `R -> quantreg -> nlrq's` solutions with `tau=0.5` in about 90% of the cases.


title: Does the precision/recall metrics on train set must be better than on test set?

text: I ran a big Neural Network with the initialized weight matrix learned by RBM (Restricted Boltzmann Machine), and the result shows something like this. On train dataset   The best precision: 0.79564. And the best recall: 0.75685. The best f1: 0.775760   On test dataset   The best precision: 0.7356. And the best recall: 0.6568. The best f1: 0.69939   I know it's very suspicious. So, I doubt that if I implement the code wrong(all the RBM and NN code was implement by myself with Matlab, and I took the `DeepLearnToolbox` library as reference). Is it absolutely(or in majority such as 90% probability) correct that the precision/recall/f1 metrics should be better on train set than on test set? Or does it have the theory basis to prove that the evaluating metrics should/must be better on train set? Or it just seems like the code I have implemented has some mistake in it? Or could you just answer the question which whether you have seen the situation that the code is correct but the metrics is like the trend of my metrics?


title: Should I use the mean-squared-prediction-error from LOOCV for prediction intervals?

text: I have a question about which prediction variance to use to calculate prediction intervals from a fitted `lm` object in R. For a certain multiple linear regression model I have obtained an error variance with leave-one-out-cross-validation (LOOCV) by taking the mean of the squared difference between observed and predicted values (i.e., mean squared prediction error). I am aware of some of the drawbacks of LOOCV (e.g., When are Shao's results on leave-one-out cross-validation applicable?), but for my specific application this was the easiest (and probably the only realistically) implementable CV method. The final fitted linear model (`fitted_lm`) is fitted with all observations and with this model I would like to make predictions for new observations (`new_observations`). For this I am using the `predict.lm` function in R.               predict(fitted_lm, new_observations, interval = "prediction", pred.var = ???)      My questions are:   * What value do I use for `pred.var` (i.e., “the variance(s) for future observations to be assumed for prediction intervals”) in order to obtain realistic prediction intervals for my new_observations?    * Do I use the error variance obtained from the LOOCV, or do I use the function’s default (i.e., “the default is to assume that future observations have the same error variance as those used for fitting”)?    * Is the mean squared prediction error not appropriate in this case? Following up on Michael Chernick's answer hereunder, I had a look in the Draper & Smith (1998) book (“Applied regression analysis. 3rd Edition”). In this book _s 2_ is defined as “variance about the regression” (p 32). This is, I presume, what we describe below as the model estimate of residual variance. Furthermore, this book mentions: > “Since the actual observed value of _Y_ varies about the true mean value _σ > 2_ [independent of the _V(Ŷ)_ ], a predicted value of an individual > observation will still be given with _Ŷ_ but will have variance > > ![formula](http://i.stack.imgur.com/uUPXs.jpg) > > With corresponding estimated value obtained by inserting _s 2_ for _σ 2_” > (pp 82-81). Thus, as far as I understand, in the D & S book they only use the model estimate of residual variance to calculate confidence intervals. This would be the default setting in the `predict` function (function help: “the default is to assume that future observations have the same error variance as those used for fitting”). However, as fosgen states below, “although LOOCV mean squared prediction error is not equal to the real mean squared prediction error, it is much more close to real than error variance of fitted model”. To make this more concrete; in my dataset I get a model estimate of residual variance of `0.005998` and a LOOCV mean squared prediction error of `0.007293`. What should I then fill in as `pred.var` in the `predict.lm` function:   * Nothing (i.e. use the default, which would equal to the model estimate of residual variance)   * `0.007293` (i.e. the LOOCV mean squared prediction error)    * `0.005998 + 0.007293` (Michael Chernick: “The model estimate of residual variance gets added to the error variance due to estimating the parameters to get the prediction error variance for a new observation”).


title: Estimating the ratio of cell means in ANOVA under lognormal assumption

text: I am conducting a two-sample test (1-way ANOVA with 2 treatments), and the goal is to estimate the ratio of cell means assuming that the data are lognormal. A simple approach is to log the response and fit a model $\log Y = b_0 + b_1 * X$ and then estimate the ratio as $R = e^{b_1}$ However, that gives the ratio of geometric cell means rather than arithmetic cell means. I assumed that if I fit a "proper" lognormal model using either `gamlss` in R or `PROC GLIMMIX` in SAS, I will get the ratio of arithmetic means, but for some reason both procedures generate the same slope as the $\log Y$ regression. This is odd because when I use this approach with Poisson or Negative Binomial regression, I do get the ratio of arithmetic means. What am I missing? Thanks, James P.S. I think I identified the source of confusion, but I don't have an explanation for it. A lognormal setup with the identity link function is: $\log Y_1 \sim N(b_0, \sigma^2)$ $\log Y_2 \sim N(b_0 + b_1, \sigma^2)$ which implies $\frac{E[Y_2]}{E[Y_1]} = \frac{e^{b_0 + b_1 +\sigma^2/2}}{e^{b_0 + \sigma^2/2}} = e^{b_1}$ To me, it means that $e^{b_1}$ should have a point estimate equal to the ratio of arithmetic means for the original response. On the other hand, $E[\log Y_1] = b_0$ $E[\log Y_2] = b_0 + b_1$ $b_0$ is estimated as arithmetic mean of $\log Y_1$, $b_0 + b_1$ is estimated as arithmetic mean of $\log Y_2$. Hence, $e^{b_1}$ should have a point estimate equal to the ratio of geometric means for the original response, and it does, given the output from those two packages. Where did I make a mistake?


title: Multiple measurements of a single datapoint: what to do with variability?

text: I am interested in a general question about research design, but I will explain it on a specific example. Suppose I want to determine if one group of people has higher heart rate than the other group. Let's consider two methods to do that, two scenarios: Scenario 1: we get **one measurement** of heart rate from 5 people in each of the two groups. Then for each group we can calculate the mean and standard deviation, and run a t-test with that data to see if groups are different. Scenario 2: we get **several, let's say 3, measurements from each** of 5 people in each of the two groups, then average the 3 measurements for each person and use it a single datapoint for this individual. Then proceed as in scenario 1. I want to understand if the second scenario is valid. Such a design is commonly used and it seems to reduce the variability compared to the first scenario. But I am very concerned that in scenario 2, we are simply discarding, or ignoring, the intra-individual variability. So is it ok to do that? What is the cost of this? Can this lead to incorrect conclusions (false-positive, false-negative)? Or maybe the conclusion has to formulated in a very specific way? What exactly are we characterizing with this method? Average heart rate of a person, as opposed to a random heart rate reading of a person? Is it what we should be interested from a biological standpoint? Sorry, if the question is not totally clear, but it is not clear to me either, and may be you can point me in the right direction. Any thoughts are welcome. Update: I just realized that there is a third scenario: start as second, but use each measurement (not each person) as a datapoint, which increases the number of observations 3-fold. What are the advantages or limitations of doing that?


title: Is this mediation or a simple path?

text: ![model](http://i.stack.imgur.com/Lw0BS.png) This is my model adapted from a study. I want to know whether I can only study it as a path analysis without studying mediation effect (1 $\longrightarrow$ 5 direct effect, as well as indirect effect) by just checking effect of one variable effect on another? I am using amos for this purpose. All paths are significant except for the one between variables 3 and 5. I have to submit it for a thesis and mediation seems to be quite difficult at this time.


title: Interaction between percentage and population size in a regression

text: Assume I have a regression type model. Are their any statistical and/or substantial **reasons NOT to include a two-way interaction which consist of a percentage and the corresponding population size**. An example would be to include an interaction between the percentage of households which own a TV and the total number of households. The interaction term would be de facto the absolute number of households owning a TV. I appreciate any comments and references.


title: Specifying variable levels in multilevel repeated measures in R using lme4

text: I am trying to analyze a dataset in which there are three measures on patients within areal units, however I am having trouble in how I am thinking about random/fixed effects and including covariates at different levels of the model. **Dataset** Dependent variable:   * Presence of disease (1=Yes, 0=No) Independent variables:   * Disease indicator (1=Disease 1, 2=Disease 2, 3=Disease 3)   * Patient age   * Patient race/ethnicity   * Patient sex   * Areal unit exposure Example of my data:               ArealID | PatientID | DV | DV_Indicator | Patient_Age | Patient_Race | Patient_Sex | ArealExposure           A |      0001 |  1 |            1 |          57 |        White |        Male | 5           A |      0001 |  0 |            2 |          57 |        White |        Male | 5           A |      0001 |  0 |            3 |          57 |        White |        Male | 5           A |      0002 |  1 |            1 |          43 |        White |      Female | 5           A |      0002 |  1 |            2 |          43 |        White |      Female | 5           A |      0002 |  0 |            3 |          43 |        White |      Female | 5           A |      0003 |  0 |            1 |          60 |        Black |        Male | 5           A |      0003 |  0 |            2 |          60 |        Black |        Male | 5           A |      0003 |  0 |            3 |          60 |        Black |        Male | 5         ... |       ... | ...|          ... |         ... |          ... |         ...           Z |      5678 |  1 |            1 |          77 |        Black |      Female | 12           Z |      5678 |  1 |            2 |          77 |        Black |      Female | 12           Z |      5678 |  1 |            3 |          77 |        Black |      Female | 12           Z |      5679 |  1 |            1 |          70 |        White |      Female | 12           Z |      5679 |  0 |            2 |          70 |        White |      Female | 12           Z |      5679 |  1 |            3 |          70 |        White |      Female | 12           Z |      5680 |  0 |            1 |          64 |     Hispanic |        Male | 12           Z |      5680 |  1 |            2 |          64 |     Hispanic |        Male | 12           Z |      5680 |  0 |            3 |          64 |     Hispanic |        Male | 12      Note that the areal exposure is the same within each areal unit and that patient age, race, and sex are the same within each patient. **Model specification** I am trying to specify my model using the `glmer` function in the `lme4` package in R. Since I have multiple observations per patient and multiple patients in each areal unit, I am specifying my model as:               m1 <- glmer(DV ~ DV_Indicator + Patient_Age + Patient_Race + Patient_Sex + ArealExposure + (1 | ArealID/PatientID), data=mydata, family=binomial)      If I understand the syntax correctly, this model should have random intercepts for both patients and areal units, but I am not sure that I have included the variables correctly since they all seem to be included at the disease level (rather than patient or areal unit). The way I thought to indicate which level each variable should be included at was like this:               m2 <- glmer(DV ~ DV_Indicator + (1 + Patient_Age + Patient_Race + Patient_Sex | PatientID) + (1 + Areal_Exposure | ArealID), data=mydata, family=binomial)      but I think this is actually including these terms as random slopes(?). I have had trouble finding clear examples of this since it seems many people mix using other packages with `lme4`