title: GLMM model specification help gender effects + an effect that is nested only within female

text: Main question: "What are the contributing variates to daily movement distances?" Specifically my question today relates to: "What is the contribution related to gender, and then within female what are the effects related to those females having young?" I have multiple readings per animal (hundreds per individual), with 13 individuals (5 females and 8 males). The females sometimes have young with them, and I know this contributes to the distance they move. I have several contributing factors in a GLMM; I am using the nlme::lme function. The current form is:               lm1 <- lme(movedistance ~ Gender+YoungPresent+x3+x4+x5+x6, random  = ~1|AnimalID/Month, data = df1)      There is no significant gender effect in this current model; however I know that this is mis-specified, because males never have young; but I don't know how to fix it. "YoungPresent" is a binary term, it is always 0 for males, and 0 or 1 for females, 1 when they have young. What I want is to somehow remove the attribution of variation by "YoungPresent" to males in the "Gender" term, but not from females. Please let me know what is the correct term for what I am looking for (Crossed? Nested?), and how I can correctly specify this structure in lme. (EDITED) After suggestions from the first response, the code now looks like this               > str(df1$YvsNY     num [1:6308] 0 0 0 0 0 0 0 0 0 0 ...          > str(df1$MvsF)     num [1:6308] -0.667 -0.667 -0.667 -0.667 -0.667 ...     > dmd <- lme(dist~Age+MvsF+TempMax+MeanRain+herb1_dens+herb2_dens+YvsNY+herb3_dens, random  = ~1|ANIMALID/Month, data = df1          > summary(lm1)     Linear mixed-effects model fit by REML     Data: df1            AIC      BIC    logLik      57915.15 57995.23 -28945.58          Random effects:     Formula: ~1 | ANIMALID        (Intercept)     StdDev:    7.923558     Formula: ~1 | Month %in% ANIMALID        (Intercept) Residual     StdDev:    7.150394  33.4111          Fixed effects: dist ~ Age + MvsF + TempMax + MeanRain + herb1_dens + herb2_dens +          YvsNY + herb3_dens                  Value Std.Error   DF   t-value p-value     (Intercept)  86.08050 10.468338 5639  8.222939  0.0000     Age           1.47128  0.967371   10  1.520906  0.1593     MvsF        -10.80126  5.214172   10 -2.071520  0.0651     TempMax      -0.58513  0.136191 5639 -4.296398  0.0000     MeanRain     -0.08233  0.020589  197 -3.998523  0.0001     herb1_dens    0.53651  0.327763  197  1.636886  0.1033     herb2_dens   -0.04928  0.032569  197 -1.513059  0.1319     YvsNY        13.07835  4.435959  197  2.948257  0.0036     herb3_dens    3.51159  1.797992  197  1.953061  0.0522      However, when I compare this to the original model formulation, I suspect that the dummy variables are not properly addressed in the lme(). I coded "MvsF" as -0.66667 for males and 0.3333 for all females, and yet the estimate, s.e. and probaliity is the same as using the original "gender" variate.               > str(df1$Gender)     Factor w/ 2 levels "Female","Male": 2 2 2 2 2 2 2 2 2 2 ...          > dmd2 <- lme(dist~Age+Gender+TempMax+MeanRain+herb1_dens+herb2_dens+YvsNY+herb3_dens, random  = ~1|ANIMALID/Month, data = df1                    Linear mixed-effects model fit by REML      Data: df1        AIC      BIC    logLik       57915.15 57995.23 -28945.58          Random effects:      Formula: ~1 | ANIMALID         (Intercept)     StdDev:    7.924217           Formula: ~1 | Month %in% ANIMALID         (Intercept) Residual     StdDev:    7.150576 33.41108          Fixed effects: dist ~ Age + Gender + TempMax + MeanRain + herb1_dens +      herb2_dens + YvsNY + herb3_dens                 Value Std.Error   DF   t-value p-value     (Intercept) 82.47974 10.290996 5639  8.014748  0.0000     Age          1.47131  0.967439   10  1.520825  0.1593     GenderMale  10.80111  5.214561   10  2.071337  0.0651     TempMax     -0.58512  0.136191 5639 -4.296341  0.0000     MeanRain    -0.08233  0.020590  197 -3.998450  0.0001     herb1_den    0.53652  0.327768  197  1.636892  0.1032     herb2_des   -0.04928  0.032570  197 -1.513020  0.1319     YvsNY       13.07854  4.436096  197  2.948209  0.0036     herb3_dens   3.51152  1.798025  197  1.952988  0.0522      I suspect that I need to call the MvsF using some kind of signal to lme() to let it know the values of MvsF are important, similarly to the way I might use factor(variate) or I(variate) inline to denote the way lme() should handle each variate. factor(MvsF clearly has no effect (basically what I have shown above), lme() does not treat the variables MvsF and Gender any differently. If what my @JakeWestfall suggests has been used correctly, then the only thing new that I have added to the model is the YvsNY where 'Males' are coded differently to 'Females with No Young', where in the original variate they were coded the same, 0. This has changed the model for sure, and looks more like its on the right path, but why did I code MvsF to THOSE values, if it changes nothing? I could easily have ONLY changed YoungPresent (0/1) to YvsNY (0, -.5, .5).... One of the problems as I see it, is that males are _still_ included in the YvsN variate - the parameter YvsN estimates a line that goes through three points on the x axis: the three levels of that variate - (-.5,0,.5 = Young, Male, No Young), and therefore Males are still contributing to the estimate of this variate - when I think they should not. I believe what I may need is similar to a grouping structure (in the random term?) where YvsN is nested within Gender (or MvsF, I think it doesn't matter) such that Males do not contribute to estimation of the YvsN parameter. \--------NEW Values of NvsNY and MvsF------------               > table(df1$YvsNY)          -0.5    0  0.5      1180 3172 1581      > table(df1$MvsF)          -0.666666667  0.333333333          3172         2761                           #Finally to check if this worked, I added a value of 2 to all the male response variates:       > df2 <- df1       > df2 <- df2[df2$Gender=="Male",]$dist <- df2[df2$Gender=="Male",]$dist +2                 # and checked that only the males' data was affected:            > tapply(df1$dist, df1$Gender, mean)        Female     Male        81.01595 92.07785        > tapply(df2$dist, df$Gender, mean)       Female     Male        81.01595 94.07785             > dmd1 <- nlme(dist~Age+MvsF+TempMax+MeanRain+herb1_dens+herb2_dens+YvsNY+herb3_dens, random  = ~1|AnimalID/Month, data = df1)       > dmd2 <- nlme(dist~Age+MvsF+TempMax+MeanRain+herb1_dens+herb2_dens+YvsNY+herb3_dens, random  = ~1|AnimalID/Month, data = df2)       > summary(dmd1)            #(truncated)                           Value Std.Error   DF   t-value p-value       (Intercept)  86.08050 10.468338 5639  8.222939  0.0000       Age           1.47128  0.967371   10  1.520906  0.1593       MvsF        -10.80126  5.214172   10 -2.071520  0.0651       TempMax      -0.58513  0.136191 5639 -4.296398  0.0000       MeanRain     -0.08233  0.020589  197 -3.998523  0.0001       herb1_dens    0.53651  0.327763  197  1.636886  0.1033       herb2_dens   -0.04928  0.032569  197 -1.513059  0.1319       YvsNY        13.07835  4.435959  197  2.948257  0.0036       herb3_dens    3.51159  1.797992  197  1.953061  0.0522            > summary (dmd2)                  #truncated                       Value Std.Error   DF   t-value p-value       (Intercept)  86.74714 10.468406 5639  8.286567  0.0000       Age           1.47128  0.967379   10  1.520896  0.1593       MvsF        -12.80125  5.214219   10 -2.455065  0.0340       TempMax      -0.58513  0.136191 5639 -4.296397  0.0000       MeanRain     -0.08233  0.020589  197 -3.998520  0.0001       herb1_dens    0.53651  0.327763  197  1.636889  0.1033       herb2_dens   -0.04928  0.032569  197 -1.513057  0.1319       YvsNY        13.07837  4.435970  197  2.948254  0.0036       herb3_dens    3.51158  1.797993  197  1.953054  0.0522            #VERY close, but a miniscule difference in coefficient had me a little worried, so I multiplied the response my 2:            > df3 <- df1       > df3 <- df1[df1$Gender=="Male",]$dist <- df1[df1$Gender=="Male",]$dist *2            # and checked that only the males' data was affected:       > tapply (df3$dist, df3$Gender, mean)          Female      Male          81.01595 184.15570                       > dmd3 <- nlme(dist~Age+MvsF+TempMax+MeanRain+herb1_dens+herb2_dens+YvsNY+herb3_dens, random  = ~1|AnimalID/Month, data = df3)       > summary(dmd3)            #(truncated)                  Value Std.Error   DF    t-value p-value          (Intercept)  121.22306 17.048079 5639   7.110658  0.0000          Age            2.75032  1.550867   10   1.773407  0.1066          MvsF        -101.41686  8.296464   10 -12.224107  0.0000          TempMax       -1.14168  0.232908 5639  -4.901840  0.0000          MeanRain      -0.13735  0.035870  197  -3.829013  0.0002          herb1_dens     0.62596  0.570363  197   1.097478  0.2738          herb2_dens    -0.12191  0.056353  197  -2.163386  0.0317          YvsNY         14.71697  7.506579  197   1.960543  0.0513          herb3_dens     6.15790  3.110133  197   1.979948  0.0491      The effect seems small, but it still seems possible to push around the YvsNY estimate, by changing males response values. This is what worries me.


title: Can someone help me interpret these findings?

text: So here is what I did, I input the command, `davis.mod <- lm(weight ~ repwt, data=Davis)`. It's a bi-variate regression that is using reported weight to estimate actual weight. When i summarized the data, this is what I get.               Call:     lm(formula = weight ~ repwt, data = Davis)          Residuals:         Min      1Q  Median      3Q     Max       -7.048  -1.868  -0.728   0.601 108.705           Coefficients:                 Estimate Std. Error t value Pr(>|t|)     (Intercept)   5.3363     3.0369   1.757   0.0806     repwt         0.9278     0.0453  20.484   <2e-16          Residual standard error: 8.419 on 181 degrees of freedom       (17 observations deleted due to missingness)     Multiple R-squared:  0.6986,    Adjusted R-squared:  0.697      F-statistic: 419.6 on 1 and 181 DF,  p-value: < 2.2e-16          > confint(davis.mod)                      2.5 %    97.5 %     (Intercept) -0.6560394 11.328560     repwt        0.8384665  1.017219      I do not know how to interpret the findings at all... Can someone explain what the intercepts and the confidence intervals mean?


title: Is this classifier random?

text: I trained a classifier and got the following stats on some test data:                   Correctly Classified Instances        1059               95.1482 %         Incorrectly Classified Instances        54                4.8518 %         Kappa statistic                          0              Mean absolute error                      0.0485         Root mean squared error                  0.2203         Total Number of Instances             1113                   === Detailed Accuracy By Class ===                     TP Rate   FP Rate   Precision   Recall  F-Measure   ROC Area  Class                  1         1          0.951     1         0.975      0.5      0                  0         0          0         0         0          0.5      1     Weight Avg.  0.951     0.951      0.905     0.951     0.928      0.5      Since my `kappa` and `ROC Area` (for both classes) are $0$ and $0.5$, respectively, have I learned a random classifier? **EDIT**               === Confusion Matrix ===              a    b   <-- classified as      1059    0 |    a = 0        54    0 |    b = 1


title: How can I 'dodge' the position of geom_point in ggplot2?

text: I am using ggplot2 in R to make plots like the following ones: ![enter image description here](http://i.stack.imgur.com/mwRYh.png) The errorbars overlap with each other which look really messy. How can I separate the errorbars for different indices? I have used position="dodge" but it seems to be not working. Here is the main part of my code:               plot =  ggplot(data,aes(x=ntrunc,y=beta_best,group=ntrunc,colour=INDEX))            +geom_point(aes(shape=detectable),na.rm=TRUE,position="dodge")              +geom_errorbar(aes(x=ntrunc,ymax=beta_high,ymin=beta_low),na.rm=TRUE,position="dodge")


title: Fuzzy Neural Network Implementation

text: Is it possible to implement the Takagi-Sugeno (T-S) fuzzy model using PyBrain ?


title: Percentile bootstrap confidence interval

text: I have seen that a number of Statistics texts state that in order to get the $1-a$ confidence interval, we have to arrange our $B$ bootstap sample estimators in order statistics: $$ \theta_{(1)} \leq \theta_{(2)} \leq \ldots \leq \theta_{(B)} $$ Then choose $m= \left[ \left( a/2 \right) B \right]$ where $[.]$ denotes the greatest integer value and $a$ is the desired significance level. Then the interval: $$ \left( \theta_{(m)}, \theta_{(B+1-m)} \right) $$ gives the desired CI. I understand the lower value, what I do not seem to understand is the upper limit of the interval, especially why we need to add 1. I would think it was like $B(1-a/2)=B-m$ but that 1 there confuses me. Thanks for your help.


title: Is there an alternative representation for $E[\max\{X,Y\}]$?

text: Suppose we have $X,Y$ i.i.d. Is there a simplified form of $E[\max\\{X,Y\\}]$? Is it just $\max\\{E[X],E[Y]\\}$? That doesn't seem right, because the latter would just be $E[X]$, and it seems like taking a max should increase the expectation. Here's what I have so far: $$\begin{align} E[\max\\{X,Y\\}] &= \int\int \max\\{x,y\\}p(x)p(y)dxdy \\\ &= \int_{-\infty}^{\infty}p(x)\left[x\int_{-\infty}^xp(y)dy + \int_x^{\infty}yp(y)dy \right]dx \end{align}$$ The term on the left is kind of $xP(x>y)$ and on the right it's sort of $E[Y|Y>X]$ but I don't really know how to simplify this. If it's relevant, I ask because I'm trying to understand this paper.


title: Forecasting differentiated timeseries with ARMA?

text: If I apply ARMA on a stationary differentiated time series and want to make forecasts with this model, the forecast will be on the differentiated values. I need the values to be non-differentiated, is it possible?


title: Which of these two ways is better to report comparison of paired data?

text: Suppose I am comparing the computational time of two algorithms **A** and **B**. They are solving randomized problem **P**. Let's say I solved the problem for seeds $ 1, 2, \dots, n $ with both A and B and got the times: $ a_1 , a_2, \dots , a_n $ for problem **A** and $ b_1 , b_2, \dots , b_n $ for problem **B**. I can use the Wilcoxon signed rank test to see which algorithm performed better. But suppose that I want to provide the reader with a sense of _magnitude_. That is I want to show just _how much_ one is better than the other. For example: > A is **12 times** faster than B on the average. How should I compute this **average**? Should I use: $$ avg = \frac{1}{n} \sum_{k=1}^{n} \frac{b_k}{a_k} $$ or: $$ avg = \frac{\sum_{k=1}^{n} b_k }{\sum_{k=1}^{n} a_k } $$ ?


title: An example: Lasso regression using glmnet for binary outcome

text: I am having confusion and difficulties using `glmnet` with lasso where my outcome of interest is dichotomous. I have created a small bogus data frame below:                age <- c(4,8,7,12,6,9,10,14,7)       gender <- c(1,0,1,1,1,0,1,0,0)      bmi_p <- c(0.86,0.45,0.99,0.84,0.85,0.67,0.91,0.29,0.88)      m_edu <- c(0,1,1,2,2,3,2,0,1)      p_edu <- c(0,2,2,2,2,3,2,0,0)      f_color <- c("blue", "blue", "yellow", "red", "red", "yellow", "yellow", "red", "yellow")      asthma <- c(1,1,0,1,0,0,0,1,1)      df <- data.frame(age, gender, bmi_p, m_edu, p_edu, f_color, asthma) # df is a data frame for further use!      The columns (variables) in the above dataset are as follows:                age (age of child in years) - continuous      gender - binary (1 = male; 0 = female)      bmi_p (BMI percentile) - continuous      m_edu (mother highest education level) - ordinal (0 = less than high school; 1 = high school diploma; 2 = bachelors degree; 3 = post-baccalaureate degree)      p_edu (father highest education level) - ordinal (same as m_edu)      f_color (favorite primary color) - nominal ("blue", "red", or "yellow")      asthma (child asthma status) - binary (1 = asthma; 0 = no asthma)      The goal of this example is to make use of lasso to create a model predicting child asthma status from the list of 6 potential predictor variables (age, gender, bmi_p, m_edu, p_edu, and f_color). Obviously the sample size is an issue here, but I am hoping to gain more insight into how to handle the different types of variables (ie. continuous, ordinal, nominal, and binary) within the `glmnet` framework when the outcome is binary (1 = asthma; 0 = no asthma). As such, would anyone being willing to provide a sample `R` script along with explanations for this silly example using lasso with the above data to predict asthma status? Although very basic, I know myself, and likely many others on CV would greatly appreciate this! Thanks!


title: How can I model optimal time at which a treatment becomes effective?

text: **In a nutshell:** Two treatments are thought to result in comparatively different outcomes depending on when they are administered, with one possibly being better for patients who are treated soon after symptom onset, and the other being better for patients who are treated later. How can I determine what the cut-off time is ? * * * **More detail:** Interest is centred on the effect of a new treatment for an emergency medical condition. Of primary concern is the time elapsed between the onset of a patient's symptoms, and the time at which the treatment is given. The outcome is binary - death. The time to treatment is measured in minutes and its log is approximately normally distributed. Data are available over five years, on all patients admitted who received either treatment. The treatment status is binary coded - 0 for the old treatment, 1 for the new treatment. In year 1, approximately 10% of patients received the new treatment. In year 5, approximately 70% of patients received the new treatment and there was an approximately linear change in the proportion during the study period. There are approximately 15,000 cases per year over the 5 year period. Preliminary analysis using a GLM with death as the outcome, time-to-treat and treatment status as the exposures, shows statistically significant estimates for both, also controlling for several confounders.                             Est    SE     Treatment    -0.72  0.04     TimeToTreat   0.13  0.03      It has been suggested that I split the dataset according to the time to treatment, for example into 2 subsets (short and long times to treatment), and look for opposite signs in the treatment effect estimate for each subset, then narrow it down by splitting at different time points until, magically, some "ideal" is reached. Apparently this would be appealing to clinicians. This does not seem like a great way forward to me - for one thing, it could be quite time consuming, and for another it seems it would give "average" effects in the short and long time groups, rather than instantaneous effects. So instead I thought about including the interaction between time-to-treat and the treatment status.                             Est    SE     Treatment    -0.06  0.35     TimeToTreat   0.19  0.05     Interaction  -0.14  0.07      `0.19` is the slope of the regression line for the `Treatment=0` group. `0.19-0.14=0.05` is the slope of the regression line for the `Treatment=1` group. Since the `Interaction` term has a negative estimate, it means that the regression line for the `Treatment=1` group is below that for the `Treatment=0` group, and when the time-to-treat variable is 0, the difference in the estimates is 0.14. Given this, it means that the regression lines for each group diverge above zero for increasing time to treat. Furthermore, since the `Treatment` effect is not significant, there is no significant difference between the estimates for the two groups, when time-to-treat is zero. Therefore the new treatment should be preferred to the old treatment for all time-to-treat. Does this interpretation make sense ? Are there better ways I could analyse these data ?


title: MLE of odd Poisson distribution

text: I have two series of $m$ observations each $n_i,X_i$ where $ i=1,....m$, and a related probability $p$. $n_i$ is large and $p$ is small, so we can assume that $X_i \sim Poi(n_ip)$. What is the maximum likelihood estimator of p in this situation and what is the expectation and variance of the estimator? Thanks in advance for any help!


title: How to interpret the significance code?

text: I'm studying for some time and I´m trying to do a logistic regression (using GLM in R) and now it´s extremely difficult to know what to do. I have a binary dependant variable and 15 independent variables. As result of running GLM I've got:               glm(formula = y ~ ., family = binomial(link = "logit"), data = crs$dataset[crs$train,      c(crs$input, crs$target)])          Deviance Residuals:          Min       1Q   Median       3Q      Max       -2.9601  -0.5093   0.3074   0.4809   2.9348          Coefficients:                    Estimate  Std. Error z value Pr(>|z|)         (Intercept) -0.52793176  0.19178032  -2.753  0.00591 **      x1          -0.03566232  0.01605379  -2.221  0.02632 *       x2           0.00497392  0.00110514   4.501 6.77e-06 ***     x3          -0.00002352  0.00000944  -2.491  0.01272 *       x4           0.01004249  0.01174335   0.855  0.39246         x5           0.10133956  0.01674787   6.051 1.44e-09 ***     x6           0.11445741  0.01819984   6.289 3.20e-10 ***     x7           0.06258882  0.01386824   4.513 6.39e-06 ***     x8           0.02266609  0.00103133  21.978  < 2e-16 ***     x9           1.05134745  0.11131339   9.445  < 2e-16 ***     x10         -0.46848579  0.09126661  -5.133 2.85e-07 ***     x11         -0.63923543  0.09490545  -6.735 1.63e-11 ***     x12         -0.36519602  0.08772172  -4.163 3.14e-05 ***     x13          0.02825176  0.00712733   3.964 7.37e-05 ***     x14          0.05050850  0.01282973   3.937 8.26e-05 ***     x15         -0.05083407  0.03117515  -1.631  0.10298         ---     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1           (Dispersion parameter for binomial family taken to be 1)              Null deviance: 31488  on 26813  degrees of freedom     Residual deviance: 19489  on 26798  degrees of freedom     AIC: 19521          Number of Fisher Scoring iterations: 5          Log likelihood: -9744.441 (16 df)     Null/Residual deviance difference: 11999.293 (15 df)     Chi-square p-value: 0.00000000     Pseudo R-Square (optimistic): 0.66283672      I´d like to ask for help to see at this result to know if its a acceptable model. I evaluated this model using ROC, and it gave me the result of 0,8843. Does it means that the prediction is 88,43% right? Can I use as reference these significance codes `(***,**,*,., )`? Which one should I choose to keep in my model?


title: Error in variables problem and maximum likelihood estimation in R

text: I am estimating cross-sectional regressions - fragment: > lm(rate~liqamih.log+cap.log+F1+F2, data=x) of the R code listed below. F1 and F2 are the coefficients estimates of time series model. In such case we need to deal with so called "error in variables problem". In literature (links: gendocs.ru/docs/23/22031/conv_1/file1.pdf (page 1091) and papers.ssrn.com/sol3/papers.cfm?abstract_id=6992 (whole article)) the MLE method is one of the effective solution of this problem. I would like to implement method introduced by KIM (reference to literature- links above) in my R code below. HOW TO DO THAT ?               data<-read.table("reg5-dane.csv", head=T, sep=";", dec=",")       data$indx <- as.numeric(gl(123*334,334,123*334))     lst1 <- split(data[,-53],data[,53]) #max 53 variables in "lst2" regression     any(sapply(lst1,nrow)!=123)     #[1] FALSE     lst2 <- lapply(lst1[sapply(lst1,function(x) !(all(rowSums(is.na(x))>0)))],                    function(x) summary(lm(rate~liqamih.log+cap.log+F1+F2, data=x)) )     capture.output(lst2,file="nooldor_regr_summ.txt")     # - f2 -f6 /+f2 -f5 +F6 - f6 - f2 - liq + f6 - f6 +f2 - f4 - f4 +f3          capture.output(lst2,file="nooldor_regr_summ.csv")     # HERE - above- IS THE CORE / most important part of question ... for now you can skip what is below      f4 <- function(meanmod, dta, varmod) {       assign(".dta", dta, envir=.GlobalEnv)       assign(".meanmod", meanmod, envir=.GlobalEnv)       m1 <- lm(.meanmod, .dta)       ans <- ncvTest(m1, varmod)       remove(".dta", envir=.GlobalEnv)       remove(".meanmod", envir=.GlobalEnv)       ans     }     library(car)     lst3 <- lapply(lst1[sapply(lst1,function(x) !(all(rowSums(is.na(x))>0)))],                    function(x) f4(rate~cap.log, x))     lst4 <- lapply(lst1[sapply(lst1,function(x) !(all(rowSums(is.na(x))>0)))],                    function(x) durbinWatsonTest(lm(rate~cap.log, x)))     lst5 <- lapply(lst1[sapply(lst1,function(x) !(all(rowSums(is.na(x))>0)))],                    function(x) vif(lm(rate~cap.log, x)))          res5 <- do.call(rbind,lst5)     res4 <- do.call(rbind,lapply(lst4,function(x) unlist(x[-4])))     library(tseries)          res6 <- do.call(rbind,lapply(lst1[sapply(lst1,function(x) !(all(rowSums(is.na(x))>0)))],function(x) {resid <- residuals(lm(rate~.,data=x)); unlist(jarque.bera.test(resid)[1:3])}) )          capture.output(lst2,file="nooldor_regr_summ.txt")     capture.output(lst3,file="nooldor_arch_test.txt")     capture.output(lst4,file="nooldor_durbin.txt")     capture.output(lst5,file="nooldor_vif.txt")      There is 123 time observations on 334 subjects observations. For each of 123 time points I am running one regressions with 334 subjects (so I am repeating it 123 times for each point of time). I post this topic also in stackoverflow.com - because I need help from one who have strong background in statistics/econometrics and also in R programming. I would appreciate your valuable help. Thank you.


title: Parameter estimation of gaussian function kernel using cross-validation

text: I need to estimate (using cross-validation), the parameters $\sigma$ and $\lambda$ of the Gaussian kernel: $K_G(x,y) = \sigma^2 \exp{(-\frac{1}{2\lambda^2}\sum_{i,j}(x_{ij}-y_{ij})^2})$ where $x$ and $y$ are the adjacent matrices of a Markov chain. How can I do that? Thank you.


title: When to use weighted Euclidean distance and how to determine the weights to use?

text: I have a set of data where each data consist of $n$ different measures. For each measure, I have a benchmark value. I would like to know how close each data is to the benchmark value. I thought of using the Weighted Euclidean Distance like this: $\hspace{0.5in} d_{x,b}=\left( \sum_{i=1}^{n}w_i(x_i-b_i)^2)\right)^{1/2} $ where $\hspace{0.5in}x_i$ is the value of the i-th measure for the particular data $\hspace{0.5in}b_i$ is the corresponding benchmark value for that measure. $\hspace{0.5in} w_i$ is the value of the weight between I will attach to the i-th measure subject to the following: $\hspace{1in}0<w_i<1$ and $\sum_{i=1}^{n}1$ However, base on this document, I found out that the weight to use is the reciprocal of i-th measure's variance. I don't think this sort of weighting will account for the importance that I will attach to each measure. Therefore:   1. Are there methods to come up with a set of weights that reflects the observer's relative importance of a measure or can the observer assign any arbitrary values for the weights?   2. Is it appropriate to use the Weighted Euclidean Distance to solve this problem? Thank you.


title: Comparing Averages of Pearson Correlation Coefficients for Significance

text: How can I test if two averages of Pearson's correlation coefficients are significantly different? For instance, I have data pertaining to participants' estimates of certain percentages (e.g. how many percent of US-citizens are taller than 6 feet) in a number of different conditions. For instance, I want to know whether the average degree of correlation between participants is significantly different across conditions. I have data from 20 participants who each answer the same 10 questions in two different conditions. Within each condition it is straightforward how to estimate the average between participant correlation; just calculate a Pearson's for the answers of each pair of participants and then divide the sum of these correlations with the number of comparisons. But how do I determine whether or not these averages are significantly distinct? At first I thought that I could use a test for the difference between non- independent r's (p.281 in Howell's "Statistical Methods for Psychology"). But occasionally I get results where the equation requires me to take the square root of a negative number which I take as evidence that using the equation is inappropriate. So I'm at a loss of how to do this. EDIT: Here is some background information: I'm interested in the relation between the average validity (understood as the correlation between a participant's answers and the correct answers) of participants in a group and the validity of the mean answer for that group. The reason I use Pearson is that there is an equation (found in e.g. Hogarth 'A Note on Aggregating opinions' in Organizational Behavior and Human Performance 21, 1978) that states the validity of the mean answer in a group as a function of the average degree of correlation between group members and group members' average validity.


title: How to interpret $R^2$ above 1?

text: > **Possible Duplicate:**   >  R^2, ESS, TSS - how well my model predicts outcome I have already posted longer version of my question, however I have obtained no answer, so I will try to make it simpler. I have calculated ESS and TSS for many cycles and my R^2 values where from 0.99 to 1.05. How can I interpret this? Do I have an error in my calculations? I am operating on values within range from 0 to 1 if it matters.   The exact way of how I calculate ESS and TSS and what I am working on is available in the question, to which I have already given link at the beginning.   Of course I do not believe that my model is so good - what I am afraid that I use inapriopriate way of verifying how well my model reflects real values.


title: Robust outlier detection in financial timeseries

text: I'm looking for some robust techniques to remove outliers and errors (whatever the cause) from financial time-series data (i.e. tickdata). Tick-by-tick financial time-series data is very messy. It contains huge (time) gaps when the exchange is closed, and make huge jumps when the exchange opens again. When the exchange is open, all kinds of factors introduce trades at price levels that are wrong (they did not occur) and/or not representative of the market (a spike because of an incorrectly entered bid or ask price for example). This paper by tickdata.com (PDF) does a good job of outlining the problem, but offers few concrete solutions. Most papers I can find online that mention this problem either ignore it (the tickdata is assumed filtered) or include the filtering as part of some huge trading model which hides any useful filtering steps. Is anybody aware of more in-depth work in this area? **Update:** this questions seems similar on the surface but:   * Financial time series is (at least at the tick level) non-periodic.   * The opening effect is a big issue because you can't simply use the last day's data as initialisation even though you'd really like to (because otherwise you have nothing). External events might cause the new day's opening to differ dramatically both in absolute level, and in volatility from the previous day.   * Wildly irregular frequency of incoming data. Near open and close of the day the amount of datapoints/second can be 10 times higher than the average during the day. The other question deals with regularly sampled data.   * The "outliers" in financial data exhibit some specific patterns that could be detected with specific techniques not applicable in other domains and I'm -in part- looking for those specific techniques.   * In more extreme cases (e.g. the flash crash) the outliers might amount to more than 75% of the data over longer intervals (> 10 minutes). In addition, the (high) frequency of incoming data contains some information about the outlier aspect of the situation.


title: Hierarchical clustering with mixed type data - what distance/similarity to use?

text: In my dataset we have both continuous and naturally discrete variables. I want to know whether we can do hierarchical clustering using both type of variables. And if yes, what distance measure is appropriate?


title: High Censoring Rate in Survival Analysis; Much higher survival time among censored patients

text: I am trying to understand censoring in survival analysis and wondering about how to tell when standard use of censoring breaks down. In one case, the number of censored patients is fairly high (low death rate), yet the median or mean "survival" time (times of last confirmed observation of the censored patient) among these censored patients with death unconfirmed is nearly twice the equivalent survival times found among the uncensored patients with confirmed dates of death. In the limiting case of a "wonder-drug" where nearly all of the patients survived beyond the trial close-out date, nearly all patients would be censored and using only the remaining uncensored patients (who may have died in the short time before the drug "kicked-in") might lead to a very poor looking survival curve that belies the activity of the drug. I understand that censoring is usually a useful method, but it seems a bit counter-intuitive to assess the effectiveness of a drug based on patients for whom the drug is least effective.


title: Optimizing parameters via simple visualisation technique

text: I hope I am at the right place, I was pointed here for further feedback from stackoverflow.   I would like to ask you, whether it makes sense to you what I am trying to do in order to optimize my model. Here is my problem I am building ranking, based on Elo approach. I have extended it with 2 additional parameters (basic form has only one - K) and I would like to check, which values for those parameters make the best predictions. In order to achieve this, I decided to make some visualisations. I am calculating (predicted score - real score)^2 for each game, in which all players had already played 20 games. Predicted score is the probability that the player will win (in a long run), and the real score is 1 for win, 0.5 for draw and 0 for lose. I decided to test several values, that I guess will work best. I have processed my whole dataset 75 (3 different values for parameter A, 5 different B and 5 different C = 3 * 5 * 5 = 75 ) times and I stored the following results in the .csv file (the file name contains information about values for parameters that were used):                ID, (predicted score - real score)^2, number of observation      ID is just the number of the cycle - in one cycle I process several thousands rows. I use the second and third column to determine what was the average error in predictions - I divide 2nd column by 3rd one. It is the metric, that I would like to use to compare different values of parameters. To sum up my questions:   1. Does my way of determining whether my model is reasonable make sense to you from statistical point of view (namely calculating (predicted score - real score)^2/(number of observation)?    2. What are your suggestions to compare obtained results in order to determine the best values for parameters? Should I make few plots in 1 chart and color them on the grounds of parameters or is there a better way?


title: What is the pdf of sample mean and variance, given the population mean and variance?

text: Suppose I have drawn n samples from a population of known mean and variance ( for example, a normal distribution with mean zero and variance 1.0 ). I then calculate the mean and standard deviation of the sample. How do I calculate the pdf of these sample values, given that I know the population values?


title: Maximum number of significance tests

text: I have a data set comprising of responses from approximately 385 people on several variables such as their shopping habits, mode of transport to the shopping centre etc. I want to do Pearson correlation on the variables (by age, gender, income level etc. of the respondents). **Is there any rule that says what would be the maximum number of tests that can be reliably performed on this data set?** I am unsure if "reliable" is the correct term; I mean "valid" or "sound" in this case. My alpha is 0.05 and I have corrected it using Sidak's method, so the experiment wise error is 0.05 for the 80 tests that I want to perform on this data set.


title: Chi-square / Fisher test variable independence for an individual

text: I have a data set of $\approx$100 individuals showing pulse rate before exercise, pulse rate after exercise, sex, smoker status, weight etc. for each individual. I've been asked to select an individual at random to test whether **their** sex and smoker status are independent, with a hint that I might want to think about chi-square and Fisher tests. If I were looking at the whole dataset, I could construct a $2\times2$ table with sex and smoker status and then run a Fisher test on that. However, I have no idea how to do it for a single individual. I can't find anything helpful in my textbooks or online either. Can anyone help? Thanks.


title: What are the mean and variance of the ratio of two lognormal variables?

text: Ratios of two normally distributed variables (e.g X/Y) have no moments (e.g. means and variances) because Y can have zero values. However, lognormal variables have no zero values. How can I calculate the mean and variance of the ratio of two lognormal variables?


title: Linearity of logistic regression

text: Given an unknown classifier, could someone make assumptions for the linearity of this model from the fact that $$y=\left\\{ \begin{aligned} 1 &\;\mathrm{if}\;& p(y=1\,|\,x) \geq 0.5 \\\ -1 &\;\mathrm{otherwise} \end{aligned} \right.$$ ? Is that the reason we conclude logistic regression is linear?


title: Testing for significance of difference between two populations

text: I am trying to test the significance of the difference between two populations. Here is a link to the population 1 data: http://pastebin.com/914mCRqA and population 2 data: http://pastebin.com/PmDgLLAb Could someone tell me whether it is best to use median or mean and hence what test to use? I was thinking of using the sign test for medians, but I can't find any tables with large n and I'm not sure if median is the most suitable measure. I've done a Z-test but not sure if that is suitable as well. Thanks.


title: statistical distribution that would generate combination of numbers without repitition

text: I have set of 'n' numbers in an array. Is there a statistical distribution that would generate randomly picked 'n' numbers from the same array such that there are no duplication. For example, if I have an array {1,2,3,4,5}, I should be able to generate 5! ways of combinations of the same numbers like {5,4,3,2,1} {2,4,3,5,1} {1,5,2,3,4} etc.. Now is it possible to generate combinations of values from a statistical distribution? apologies for any mistakes that I assumed/made. correct me if i am wrong.


title: Fitting parameters of a function

text: Imagine a graph where the x-axis is time in minutes (continuous) from 0 to 360. And the y-axis represents 'z' values. Of this continuous time course graph, I have 8 time points w/ their corresponding real z-values. Here is the data frame 'real_df' with that real data. real_df:               timepoints     real_z              0  1.0000000             20  0.7202642             40  0.8000305             60  1.7430986            120 16.5172242            180 25.6521268            240 33.9140056            360 34.5735984      Now, I have 'z' values that I'm making artificially for the same time points of the real 'z' values. I use 'a' and 'b' values as shown below, to find the modeled z-values. And this 'mod_df' contains hundreds of combinations of 'a' and 'b' values, but here is an example of 2 combos of 'a' and 'b': mod_df                 a   b   t     model_z             30 120   0    0.000000       30 120  20    0.000000       30 120  40    2.881133       30 120  60    8.643400       30 120 120   25.930199       30 120 180   34.573598       30 120 240   34.573598       30 120 360   34.573598       30 121   0    0.000000       30 121  20    0.000000       30 121  40    2.857322       30 121  60    8.571967       30 121 120   25.715900       30 121 180   34.573598       30 121 240   34.573598      30 121 360   34.573598       What I want to do is find a single 'a' and 'b' combination which results in the best fit of all 8 modeled z-values to the 8 real z-values (i.e., fit the modeled z-values to real z-values at corresponding time points). I've been thinking a least squares fit would be just fine, but I'm curious if there is a better way to fit. So far, I've been making 'mod_df' by running all combos of 'a' and 'b' that I put into an equation I've created. But I hear that there may be a way to get a better estimate for my variables 'a' and 'b'. Because now, obviously, they're limited to the values of 'a' and 'b' I put into the equation. Maybe what I'm saying isn't true could someone shed some light on this? Least squares should work fine for me now, but could I get some help w/ writing a least squares fit based on the 8 modeled z-values vs. 8 real z-values? Here is some code that I have, but it actually squares the difference of the sum of all 8 values. Which will not give me anything very precise.               mod_df[seq(from=1, to=nrow(mod_df), by=8),1:2][which.min(colSums((matrix(mod_df$model_z,                 nrow=8) - real_df$real_z)**2)),]      There's probably a better way of doing this, could I get some help?


title: New to R: computing the correlation of the residuals with the fitted values of a data set

text: I am new to R and trying to practice with some exercises. Given a data set with 40 observations and 5 variables. Spending is the the response and there are 4 predictors. I started with a linear model Residuals:                   Min      1Q  Median      3Q     Max      -51.082 -11.320  -1.451   9.452  94.252           Coefficients:                  Estimate Std. Error t value Pr(>|t|)         (Intercept)  22.55565   17.19680   1.312   0.1968         sex         -22.11833    8.21111  -2.694   0.0101 *       status        0.05223    0.28111   0.186   0.8535         income        4.96198    1.02539   4.839 1.79e-05 ***     verbal       -2.95949    2.17215  -1.362   0.1803         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 22.69 on 42 degrees of freedom     Multiple R-squared: 0.5267, Adjusted R-squared: 0.4816      F-statistic: 11.69 on 4 and 42 DF,  p-value: 1.815e-06       First, is this what they mean by fit regression model and Secondly, how do I compute the correlation of the residuals with the fitted values?


title: Clarification of an equality involving conditioning

text: In this wikipedia section, the first block equation claims that $P(n_b \leq n^* | s+b)=P(n\leq n^* | b)$ Some context (also found in that linked section): $n_b$ follows $Pois(b)$ and $n_s$ independently follows $Pois(s)$, $b$ is known and $s$ is to be estimated, and only $n=n_b + n_s$ can be measured. Wikipedia claims that the equality above is "easily seen", but I don't see it even after banging my head for an hour. How is that true?


title: Classification More Robust Than Regression

text: Are neural networks using classification more robust / reliable than using regression to produce a single value? The only reason I would think so is that it would be easier for the network to adjust the weights to accommodate a fixed number of outputs as opposed to manipulating the weights in order to output any number (within a reasonable range). I am wondering if it is worth while to reword a problem into classification? Also... if so, does using a catch all output, like an else statement, bad practice?


title: Is binning a skewed Likert scale variable justifiable?

text: I have recently advised some colleagues on the malpractice of binning a continuous variable, which was used in order to put it as a covariate in a regression model and retained as a significant predictor. They responded that when summed responses to several Likert item are not normal (in this case extremely skewed) it is questionable whether these can be used as continuous variables. For this reason, they treated this variable as ordinal data with dummy categories. Is it ok to do so? Isn't it better to transform the variable in order "normalize" it (with methods like square root or loglinear transformations, or aggregating the less frequently observed responses)? Doesn't the binned variables lead to several shortcomings (biased estimates etc.)? Thanks in advance for your answers


title: Predictor variables in Multinomial Logistic Regression

text: Is it always better to have categorical predictor variables when performing a multinomial logistic regression analysis? or can it be done using continuous predictor variables? If a predictor variable should always be categorical, what is the reason behind it?


title: Chi Square test for survey data

text: I am working with a product that turns survey data into useful statistics. Reviewing their code, has made me somewhat nervous, and I'm not a statistician, so I hope I can ask for clarity of the following problem: Out of a survey S, for a product P. Respondents where asked if they   1. liked the product   2. were indifferent   3. hated the product The group of respondents were separated into men and women. The chart supplied by the software when crunching some survey data, says that "Men are significantly more likely to be likers." OR "Men and Women..." Or "Women..." For me this already raises issues:   1. Men are significantly more likely to be likers than what?   2. Men and Women are more signicantly likely to be likers than what?    3. How are these things measured?   4. What test is being used?... etc. When I had a look at the code, I noticed they were using a chi-test(!). I had to ask what exactly the null hypothesis was, because this was making less and less sense. Apparently the null hypothesis is that " _the chance that men and women are likers is the same_ " ...ok, fine. But wait. So, we have the following table:                               Men       Women      Total      likers         54        46         100           indifferent    23        26         49           non-likers     22        31         53           Total          99        103        202      We can populate the expected distributions for all three rows:                               Men          Women      Total      likers         54-49        46-50      100           indifferent    23-24        26-24      49           non-likers     22-25        31-27      53           Total          99           103        202      The code then populates a matrix with chi values based on the above. The programmer decided that the degrees of freedom when doing these calculations was (m-1)(n-1) = 2, which at this point made me think the null hypothesis was rather that _if you are a liker, indifferent or a non-liker, there is equal probability that you are a man or a woman_. We're using a 90% confidence level, so all I imagined we needed to do was to sum over all the 6 chi values, and compare that with a critical value given by the degrees of freedom and the confidence interval. From that point we could say with 90% certainty that men and women were equally as likely to be a liker, etc... or reject the N.H. This is what the code does instead:   1. It uses 1 degree of freedom instead of 2 (still at 90%), so we have a new critical value 2.706   2. For each row (liker, etc...) of the chi value matrix, if an element is greater than the critical value reject the null hypothesis, and add the element to a 'significance' list.  To illustrate, it looks at [likers;men] > cv i.e. chi_value[0][0] > cv, if that is true, reject N.H., and add 'men' to the list. On the chart this result is reflected as: _men are more likely to be likers_. For me this single evaluation of men and women for each row seems wrong. It doesn't make sense to make pronouncements about two variables when you're only looking at one... I am not nearly as smart as my boss, but I feel like something has gone wrong here and I would appreciate it if someone could help clarify this. Lastly the client has asked to know the _% more likely men are to be likers than women_ \-- I think this is an erroneous request, as a chi square test does not address questions of which is greater or smaller, but only serves to confirm that a set of variables are independent. Am I right? I just want to add, that I used the following statement to guide my thinking: _Cautionary Note It is important to keep in mind that the chi-square test only tests whether two variables are independent. It cannot address questions of which is greater or less. Using the chi-square test, we cannot evaluate directly the hypothesis that men are likers more than girls; rather, the test (strictly speaking) can only test whether the two variables Like and Gender, are independent or not._


title: How to interpret a control chart containing a majority of zero values?

text: I am using a control chart to try to work on some infection data, and will raise an alert if the infection is considered "out of control". Problems arrive when I come to a set of data where most of the time points have zero infection, with only a few occasions of one to two infections, but these already exceed the control limit of the chart, and raise an alert. How should I work on the control chart if the data set is having very few positive infection counts? Thanks!


title: Hierarchical clustering: is it possible to combine single-linkage clustering and average linkage clustering?

text: A "seismic section" shows amplitude for m discrete x values along its horizontal axis times n discrete time values along its vertical axis: ![Seismic section](http://i.stack.imgur.com/L6zwW.jpg) Peaks in amplitude (black) are centered on horizons; interfaces between different geological layers. The horizons are mostly continous but not necessarily horizontal nor flat lines. In the image above somebody has "tracked" the green and yellow horizons (using a crayon on a piece of paper?) that form the top and bottom of the "Ordovian" geological layer. These two lines are continous everywhere except for where they cross one of the vertical black lines: faults. It is easy to implement a computer program that "draws" the colored lines by following the peaks, except for at the faults. The problem here is that the peak "jump" in time and by how much differs from case to case. One possible way to automate this is to assume that the shape of the signal around the peak is almost constant along a horizon. We can then reconstruct the amplitude around each peak sample by e.g. a Cosine series expansion with say 5 terms: ![Cosine approx](http://i.stack.imgur.com/1XZNe.png) The result is a feature vector for each peak sample. In the case above the feature vector would be:               F1 = [94.2496  192.7706 -211.4520  -82.8782   11.2105]      The feature vectors can then be used as input to a classification/clustering algorithm that figures out which peak samples are on the same horizon. The problem with the feature vector above is that it does not force the peak samples to be linked together; in fact it disregards their positions in space altogether. This would be fine if each horizon had a very different shape around the peak. Unfortunately this is often not the case. So we want to include the x postion into the feature vector. However many clustering algorithms are not able to deal with this since they would require that all observations within a cluster have almost the same x position. **Perhaps hierarchical clustering may be a solution to this problem?** As I understand "single-linkage clustering" would only require that each peak had a peak within the same cluster that was its neighbour (adjacent x positions). This is almost what we want. Ideally I would want single-linkage clustering on the x-position and average linkage clustering on the Cosine coefficients. **Is this possible?** Thanks in advance for any answers!


title: Weighting first and second level

text: I have a two level nested structure data and two weights: for the first and second level. I normalized each weight so that the sums of each weight is equal to the total sample size at the first level. Is this correct or should the second level weight only be equal to the level 2 sample size?


title: Is the AR(1) process always Gaussian given Gaussian innovations?

text: I found that $AR(1)$ process $x_t=\phi x_{t-1}+\epsilon_t$ is not always Gaussian given Gaussian innovations $\epsilon_t$. This only happens when the $AR(1)$ model coefficient is very large. This goes against with the theory (If the white noise $\epsilon_t$ is a Gaussian process then $x_t$ is also a Gaussian process). How to explain this then? The following R codes gives that only $363$ out of the 1000 generations of $AR(1$) processes are gaussian process by the `Shapiro-Wilk test`.               decision=c()     for (k in 1:1000){ # 1000 runs        phi=0.9 # AR(1) coefficient       x=0.5 # value of x_0       X=c() # generated AR(1) process       for (i in 1:1000){  # length of the process is 1000         x=phi*x+rnorm(1) # Gaussian innovations         X=c(X,x)        }         decision=c(decision,shapiro.test(X)$p.value) # Gaussian test     }     sum(decision>0.05)/1000 # percentage of Gaussian processes


title: Is this process an AR(1)?

text: The following process $X$ is generated as the sum of a structure $y$ and an $AR(1)$ process with $\phi=0.4$ and $\sigma_{\epsilon}^2=1$:               N=200 # length of time series     t=1:N     A=20 # amplitude of the structure     y = rep(0,N)         a = 30     b = 170     y[a:b] = A * (t[a:b] - a)/ (b-a) # the structure (see the plot below)          X=y+arima.sim(list(order = c(1,0,0),ar=0.4),n=N,sd=1) # the process          model=ar.ols(X) # fit AR process          Box.test(model$resid,lag=round(log(N)),type="Ljung-Box",fitdf=1) # white noise test for the model residuals          par(mfrow=c(1,3))     plot(y,type='l',ylab='y',xlab='')     plot(X,type='l',ylab='X',xlab='')     plot(model$resid,xlab='')      ![enter image description here](http://i.stack.imgur.com/pSgqP.png) In the end, the test tells the model residuals are independent with $p=0.8327$ (please click here), which means $X$ can be modelled as $AR(1)$ process?


title: Algorithm to measure statistical probability of proximity

text: I am trying to find a simple algorithm to measure the probability of the proximity distance between location of randomly placed actual target and subject's guess of its location on an 8 x 8 grid (64 squares). I.e. if the guess is a MISS, what is the probability of the guessed location being located on any given square?


title: Limitations of Cohen's kappa for sparse data?

text: Are there any limitations for using Cohen's kappa with sparse data? I need inter-rater agreement between 2 raters for ~15 items, and the data in the contingency table is quite sparse (0 in some cells, lots of values < 5).


title: What is a recommended Python framework for Bayes Nets

text: What is a good Python framework to use for statistical analysis using Bayes Nets. The following statistical frameworks in one way or other do not support expected features. Python Scikit does not have an implementation of Bayes Net. Looked into other Python implementation such as: Libpgm - not clear if it has all functionality of predictions Pebl - no parameter inference and not active Python Bayes Network Toolbox - not clear about its licensing BayesPy - not sure if it fully implements structure learning Orange - not much documentation PyMC/OpenBUGS - There are good implementations in non-python community such as "Bayes Net toolbox" and Weka. I was wondering if there is a reason for the gap; especially given Python and its foot prints in machine learning, I feel I am missing a link I need. Please share your experience.


title: How do I calculate a weighted standard deviation? In Excel?

text: So, I have a data set of percentages like so:               100   /   10000   = 1% (0.01)     2     /     5     = 40% (0.4)     4     /     3     = 133% (1.3)      1000  /   2000    = 50% (0.5)      I want to find the standard deviation of the percentages, but weighted for their data volume. ie, the first and last data points should dominate the calculation. How do I do that? And is there a simple way to do it in Excel?


title: Testing the parameter of a sub-sample

text: I have some survival data from which I compute the annual failure rates. The parameter in which I am interested is the mean failure rate. Now I select certain subjects that have a given characteristic, so the sub- sample is not random. For example, I have very few subjects having the characteristic in the beginning and I have more and more subjects with time. I want to "prove" that the change in the mean failure is not statistically significant. To do so I use bootstrapping. I select randomly sub-samples from the data having the same number of subjects per year as the original sub-sample and compute the mean failure rate. Then, I can compute the a p-value using the bootstrap distribution and eventually reject or not the null hypothesis. I have some questions regarding this procedure:   * is it reasonable?   * what would be the alternatives?   * what is in this case (bootstrapping) the difference between hypothesis testing and confidence intervals? ## **EDIT:** Let me give an example of type data I have. A number of N subjects are followed over the years. The population is not constant, each year some new subjects can be added or some subjects can be removed. The subjects can be in 3 states {"healthy","ill","dead"}. Obviously, "dead" is an absorbing state, but recovery is possible from "ill"ness. I am interested in the mean failure rate per state ("healthy","ill"). I estimate annual failure rates per state with a cohort approach: I compute the number of subjects that are in states "healthy" and "ill" at the beginning of each year and look which ones have died at the end of the year. The estimator is imposed, I do not want to change it for now. But advices are welcome! Now, assume that I have 2 types of subjects: presence of a certain gene or not. I am interested in the effect of this gene on the mean failure rates. Due to selection of subjects, it happens that there are few subjects with the gene at the beginning of the study and more and more subjects with the years. I want to test if the mean failure rates of the subjects with the gene are the same as the global ones. I have used bootstrapping to compute a p-value. Assuming that I have $h_t^w$ healthy and $i_t^w$ ill subjects with the gene at year t, I sample with replacement from $h_t$ (healthy) and $i_t$ (ill) subjects (with or without the gene) $h_t^w$ and $i_t^w$ subjects respectively and look how many have died.


title: Can I use Kolmogorov-Smirnov test and estimate distribution parameters?

text: I've read that Kolmogorov-Smirnov test should not be used to test the goodness of fit of a distribution whose parameters have been estimated from the sample. Does make sense to split my sample in two and use the first half for parameter estimation and the second one for KS-test? Thanks in advance


title: Bayesian, MDL or ML interpretation of cross-validation?

text: Is there any known Bayesian, ML or MDL interpretation of cross-validation? Can I interpret cross validation as performing the right update on a specifically crafted prior?


title: Collinearity diagnostics disagree - VIF, condition index, and correlation matrix

text: I'm working with a large dataset consisting of just over 1 million cases. The data are longitudinal covering 14 years and hierarchical with about 500 of the level 2 units. Each case is a criminal sentence outcome (prison vs. no prison) and the level 2 units are the judges handling the cases. I'm trying to quantify the effects of experience. For example, I want to know what happens when a judge handles more murder cases - does he become more lenient with the subsequent murder offenders or more punitive? I see two ways to operationalize experience. 1\. include a running count of the number of xyz type cases the judge has handled. 2\. calculate a running proportion of the number of each judge's cases that are of each type. In both cases, a cross level interaction term between the crime type variable and count/proportion variable is necessary. There are 9 different crime categories, so I include 8 dummy variables. Each crime type dummy has a corresponding count (or proportion) variable that goes along with it. However, one of the 8 dummies is a catchall category - 'other' - so I didn't include a count for that because it doesn't really make sense. The problem I'm having is that collinearity diagnostics suggest that the 7 count/proportion variables are problematic. The proportion variables look ok on some diagnostics and bad on others while the count variables look others and bad on some. I don't know why the diagnostics disagree and which ones I should trust. In models with either set of variables (count or proportion) the variables are generally statistically significant and the standard errors are not large but that's probably attributable to the large number of cases. Here's how the diagnostics shake out. Counts:   Very high bivariate correlations of .7, .8, and even .9. The condition index calculated from uncentered, unstandardized variables with a constant (Stata's 'coldiag2) is 24.5 - which is high but acceptable. The VIFs on the other hand, calcuated using Stata's 'collin', are impossibly high - 18, 11, 17, 10, and a few 'low' ones under 5 Proportions: The correlations are all modest at about .2 or less, although one is as high as .5. The VIFs are also great, the mean VIF is just 1.44. But the condition number is even worse than the count variables - 32.5. So the diagnostics disagree. I know the correlation matrix is not the best tool, but the VIFs and condition number wildly disagree and I'm not sure what to make of that. * * * Putting aside all that, I 'want' to use the count variables instead of the percentages. Read on if you want to see why, but these issues are secondary to the collinearity issue outlined above. The proportion variables take on extreme values (i.e. 100% or 1) for the first case and then regress toward the mean as the judge handles more cases. Thus, the change in the percentage variables really just reflects regression toward the mean. Every additional case is an additional 'unit' of experience and yet additional cases _only_ change the percentage variables if they change the _distribution_ of the judge's caseload. That's something that's hard to do when the judge has heard a thousand or more cases. I think these limitations of the percentage variables can be at least partly addressed by the inclusion of a running count of the total number of cases the judge has handled but that's not a panacea and it will make the collinearity problem worse. On the other hand, the count variables are overdispersed as you would expect any count variable to be. * * * So, to summarize, the problems are twofold - difficulty operationalizing "experience" and disagreeing collinearity diagnostics. Perhaps there's a better way to operationalize experience than either counts or proportions that would not generate collinearity problems. I would welcome any advice about alternative model specification OR correct interpretation of the collinearity diagnostics.


title: Interpreting 2D correspondence analysis plots

text: I've been searching the internet far and wide... I have yet to find a really good overview of how to interpret 2D correspondence analysis plots. Could someone offer some advice on interpreting the distances between points? Perhaps an example would help, here is a plot that's found on many of the websites I've seen that discuss correspondence analysis. The red triangles represent eye colour and the black dots represent hair colour. ![alt text](http://i.stack.imgur.com/nVOUx.jpg) Looking at the graph above, could you make a few statements about what you see in these data. Points of interest about the different dimensions and relationships between triangles and dots? An explanation of row-points verses column points, and the use of the word "profile" with particular focus on the example would be instrumental.


title: What is the intuition behind the variation of information (VI) metric for cluster validation?

text: For non-statisticians like me, it is very difficult to capture the idea of `VI` metric (variation of information) even after reading the relevant paper by Marina Melia "Comparing clusterings - An information based distance" (Journal of Multivariate Analysis, 2007). In fact, I am not familiar with many of the clusterings' terms out there. Below is a MWE and I would like to know what does the output mean in the different metrics used. I have these two clusters in R and in the same order of id:               > dput(a)     structure(c(4L, 3L, 4L, 4L, 4L, 4L, 4L, 3L, 4L, 4L, 3L, 3L, 3L,      4L, 4L, 4L, 4L, 4L, 3L, 4L, 4L, 4L, 3L, 4L, 4L, 3L, 4L, 4L, 4L,      1L, 1L, 4L, 4L, 4L, 4L, 4L, 4L, 3L, 4L, 4L, 3L, 4L, 4L, 2L, 2L,      4L, 3L, 3L, 2L, 2L, 2L, 4L, 3L, 4L, 4L, 3L, 1L, 4L, 3L, 4L, 4L,      4L, 3L, 4L, 4L, 4L, 4L, 2L, 2L, 2L, 4L, 3L, 4L, 4L, 4L, 4L, 4L,      4L, 4L, 4L, 4L, 3L, 4L, 4L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 3L, 4L,      4L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 3L, 4L, 4L, 4L, 4L, 2L, 2L, 4L     ), .Label = c("1", "2", "3", "4"), class = "factor")     > dput(b)     structure(c(4L, 3L, 4L, 4L, 4L, 4L, 4L, 3L, 4L, 4L, 3L, 3L, 3L,      4L, 4L, 3L, 4L, 4L, 3L, 4L, 4L, 4L, 3L, 4L, 4L, 3L, 4L, 4L, 4L,      1L, 1L, 4L, 4L, 4L, 4L, 4L, 4L, 3L, 4L, 4L, 3L, 4L, 4L, 2L, 2L,      4L, 3L, 3L, 2L, 2L, 2L, 4L, 3L, 4L, 4L, 3L, 1L, 4L, 3L, 4L, 4L,      3L, 3L, 4L, 4L, 4L, 4L, 2L, 2L, 2L, 4L, 3L, 3L, 4L, 4L, 4L, 4L,      4L, 4L, 4L, 4L, 3L, 4L, 4L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 3L, 4L,      4L, 3L, 4L, 4L, 4L, 4L, 4L, 3L, 3L, 4L, 4L, 4L, 4L, 2L, 2L, 4L     ), .Label = c("1", "2", "3", "4"), class = "factor")      Now doing comparisons based on the `VI` as well as other metrics / indices and in chronological order of their appearance in literature.               library(igraph)       # Normalized Mutual Information (NMI) measure 2005:     compare(a, b, method = c("nmi"))      [1] 0.8673525       # Variation of Information (VI) metric 2003:     compare(a, b, method = c("vi"))      [1] 0.2451685       # Jaccard Index 2002:     clusteval::cluster_similarity(a, b, similarity = c("jaccard"), method = "independence")      [1] 0.8800522       # van Dongen S metric 2000:     compare(a, b, method = c("split.join"))      [1] 8       # Adjusted Rand Index 1985:     compare(a, b, method = c("adjusted.rand"))      [1] 0.8750403       # Rand Index 1971:     compare(a, b, method = c("rand"))      [1] 0.9374788      As you can see, the `VI` value was different from all the others.   * What does this value tell (and how is it related to the figure below)?    * What are the guidelines for considering this value low or high?    * Are there any guidelines defined?  Maybe experts in the field can provide some sensible descriptions for laymen like me when trying to report such results. I would really appreciate if someone would provide also guidelines for other metrics as well (when to consider the value is large or small, i.e., in relation to a similarity between two clusters). I have read related CV threads here and here, but still couldn't grasp the intuition behind `VI`. Can someone explain this in plain English? The below figure is figure 2 from the above mentioned paper about `VI`. ![enter image description here](http://i.stack.imgur.com/OnCzM.png)


title: What statistical test would be appropriate to analyse these proportions?

text: I have this matrix:               habitats <- matrix(c(0.4,0.1,0.5, 0.2,0.6, 0.2), nrow = 3, ncol = 2,                    dimnames = list(c("grassland","intertidal","saltmarsh"),                                    c("site1","site2")))          habitats                     site1 site2     grassland    0.4   0.2     intertidal   0.1   0.6     saltmarsh    0.5   0.2      I want to know if the habitat in site1 differs from the habitat in site2. Each site is split into three habitat types and the matrix shows the proportion of each habitat in each of site1 and site2. Initially, I thought a chi-squared would be appropriate but I have been made aware chi-squared should only be used if data are frequencies. What test could I use to tell if the habitat differs between site1 and site2?


title: Is this a 2-ways repeated measures ANOVA?

text: Imagine a long, horizontal, transparent acrylic tube which is marked with duct tape every 50 cm in such a way that you can distinguish 10 "sections" of the tube. I want to study the movements of ants in these plastic tubes, which are under four combinations of two different conditions, namely each tube has different texture of the inner wall and used a different (transparent) paint as repellent. Therefore I consider TEXTURE and PAINT my two factors, each one with two levels: Texture1-Texture2 and Paint1-Paint2. We choose 100 random ants and put them in each tube, evenly distributed along the tube. We wait 5 hours and count how many ants we find in each one of the 10 sections of the tube. So we have 4 tubes: Texture1 + Paint1 Texture1 + Paint2 Texture2 + Paint1 Texture2 + Paint1 In each tube we have 10 values (one number of ants per section). Eg. in a given tube I have Section1 with 10 ants, Section2 with 34 ants, etc. until Section10. In total, all 10 sections of a tube sum 100 ants. In this way we can see if the horizontal distribution of the ants (eg. where is the maximum number of ants located inside the tube? in which section?) is related to the texture and/or paint used. So my initial thought was a 3-ways ANOVA with Texture, Paint and Section as factors. My dependent variable is simply ANTS: the number of ants in a given section of a given tube. But after thinking for a while I realize that the number of ants in each section is actually related, because the initial number of ants is fixed and therefore if I have 50 ants in one section then I have 50 less ants for the rest of the sections, right? Then I thought maybe this would be for a 2-Ways ANOVA, but with repeated measures. In this case the factors are Texture and Paint, and the repeated measures are Ants in each Section. So for each combination of Texture and Paint I would have 10 repeated measures. However I am not sure because repeated measures is almost always related to time in the textbooks and examples (eg. measure something "before" and "after").. Is my reasoning correct? Can I use repeated measurements with the justification that the measurements are not independent even when time is not a factor or variable here?


title: How do I find co-occurring labels?

text: I have a matrix that is of the shape `n_samples` by `n_features`, where `n_samples` ~ 7000 and `n_features` ~ 800. To clarify, `n_features` are categorical labels, each representing an amino acid position. I'd like to know which columns have amino acid positions which are highly correlated with each other. For example, one might have a protein where amino acid position 255 = H when 288 = Y ("HY" combination), and 255 = N when 288 = P ("NP" combination). This is a completely hypothetical scenario, btw. In the data set, I might find, for example, 3500 instances of HY, 100 instances of HP, 3400 instances of NP. In this case, we would say that position 255 and 288 are correlated with each other. The problem is, I'm not sure how to compute this. Using scikit-learn, I have transformed the data into numerical categorical labels, where an amino acid position is now represented by numbers (i.e. from 0 to 20). So now, it's `rows=each sequence`, and `columns=each position`, and the `numbers in each cell=0 to 20`, each letter represented by a number. What would your advice be for going from this data matrix to finding columns that "covary" together?


title: SEM using composites or indicators?

text: I am new to SEM. I am analyzing a dataset with five constructs and 21 indicators. I found there could be two ways to run the structural model in AMOS. One is to include all the measured indicators in the model and the other one is to use the composites imputed by AMOS. Model fit of the model with indicators was poor, while the one with the composites were great.   * What is the difference between the two ways?   * Which one is more appropriate?


title: How to compare a sample against some baseline data?

text: First let me describe the situation I'm dealing with: I'm looking at performance data for a software system. I have data for many versions of the software, including ongoing. For each version I have a series of values for the run-times of tests run against it, usually about 20 (of the same test). These generally look normally distributed (for each version). What I want to know is, given some data for runs against a new version, is the new distribution _different_ in a way that merits investigation. The comparison could be either against the previous version, or against a set of previous versions that have been selected as a having "stable" performance: a kind of baseline. Any kind of change could be relevant: a change in the mean, variance or shape could be significant. Now, the different versions are, well, different, but for those in the baseline I think I can assume that they're effectively samples from the same distribution. So I've got $X_1, ... X_n, X'$, for $n \geq 1$, and I want to test whether $X'$ is relevantly similar to the $X_i$s in an automated fashion. From what I've seen on the internet I've come up with a few options: a) Kolmogorov-Smirnov test: either of $X'$ vs $X_n$, or of $X'$ vs $\bigcup X_i$. b) T-test: similarly. c) Mann-Whitney/Wilcoxon test? Firstly, I'm not clear which would be better for my situation, as they both test for different kinds of "similarity", or whether I should use both and report some combination of the results. Secondly, looking at the data, it looks like the $X_i$, while normally distributed, tend to move around a bit: so their means vary, but e.g. the variance is similar. If I just lump them into one big sample, then this information is lost; for example, it will look like the typical variance is much bigger than it actually is. For that reason I wondered whether I might be able to instead look at the distribution of the means and variances of the $X_i$s, and compare the mean and variance of $X'$ against that, somehow. However, I'm unsure how to do that: I think the means should follow a t-distribution, so I ought to be able to estimate the probability of getting $\bar{X'}$ given that the sample came from the same distribution, but that's the wrong conditional probability! (although that's kind of just what a p-value is...) I can't do a proper t-test as that requires either more than one value to compare against, or the assumption that the variance is shared, which I'm not sure I have. Finally, my inner Bayesian feels like I ought to be able to do better than producing p-values for rejecting the null hypotheses: surely I ought to be able to calculate a posterior probability that $X'$ is, say, drawn from the same distribution as the $X_i$s? Apologies for the huge question; I hope this gives a reasonable idea of where I'm coming from! I'm mathematically trained, but I'm pretty unfamiliar with statistics, so I can cope with some maths. Edit: I'm also familiar with R; I'm going to be using it to do the calculations... once I figure out what to calculate!


title: How to setup a laboratory experiment in Ecological Research under high natural variability

text: I am about to do a laboratory experiment in the scientific field of soil ecology and hydrology. Beforehand I want to make sure not to make any crucial mistakes, and therefore I would appreciate any hints and comments from your side. The main issue is how to deal with high natural variability (~25%) and a rather small sample size (total max=20). **Short Description of the Experiment:** We will put soil cores in cylinders and keep them under different moisture scenarios. Three different kind of Carbon forms will be measured for about one month.   We want to know if the variables change within groups and between groups. The experimental design that was proposed is as follows:   There will be 2 different soil types, 2 different treatments + one control, and each treatment should be replicated three times. The total number of cylinders is thus 18 cylinders. The variability in field measurements can be as high as 25% within one group. Due to pracitcal reasons there can not be more than 20 cylincers in total. My Questions:   1. Would it make more sense to have only one treatment and one control, but each one is replicated five times?    2. Will I be able to draw any reliable conclusions from this kind of experiment under these conditions?    3. How should I set the parameters to make some power calculations using e.g. G*Power 3? Which test should I choose? What should I set the effect size to and what should the numbers for df be?    4. How should I analyze the data after the experiment is done? Should I use ANOVA? Can I use a mixed effect model?


title: Non-response bias

text: I have the following question: > In a survey, a simple random sample of 1,000 households was drawn to > determine the distribution of household size in a city. Interviewers were > required to visit the house to carry out the survey. Interviewers could only > reach 653 households in the sample. The rest of households was not reachable > due to no one was at home to answer the house visit. That is, there was a > non-response rate of 34.7%. > > The researcher involved in this survey decided to draw a second batch of > simple random sample of 1,000 households, and used the first 347 completed > interviews in the second batch to bring up the sample up to its original > planned size of 1,000 households. The survey result showed a total of 4,087 > people in these 1,000 households, and the average household size in the city > was estimated to be about 4.1 persons. > > Choose one of the following options which is the most appropriate. > > (a) The above described method of drawing sample would be a good way to fix > the problem of non-response bias. > > (b) The 347 households used in the second batch would be quite different > from those unreachable 347 households in the first batch. > > (c) The 653 households from the first batch would be a good representative > of the population. > > (d) The final sample of 1,000 household used in the survey can be considered > as a random sample, as both batches of 1,000 households were drawn by simple > random sampling method. > > (e) The estimation of 4.1 persons per household in the city is likely to be > about right. It seems (b) is the most plausible answer, but I don't know how to prove the other statements wrong. Any help is appreciated, thanks!


title: ANCOVA: How does a covariate give insight when comparing the mean to the adjusted mean?

text: When comparing the direction and magnitude of the difference between the original means and the adjusted means what does the change say about the covariate? For example: IV = Grade, DV = Test Score, Covariate = IQ Score Am I correct in assuming that if a group's **adjusted mean** is higher than its original mean that the effect of the covariate (higher IQ) is that it decreases the DV (Test Score) thus resulting in the higher adjusted mean? Thanks.


title: Is it possible to use factor analysis on census data?

text: I have some census data of 20 urban regions. I want to use them to combine 5 criteria to determine which region is most deprived. Is it acceptable to use factor analysis?


title: Express product as convolution? Or, given $f(s)$, find $g$ satisfying $f(s)=\mathbb{E}[g(X)]$ where $X\sim \mathcal{N}(0,s^2)$

text: Given a function $f(\mu)$ (satisfying certain properties), it is possible to find a function $g(x)$ such that $f(\mu)=\int_{-\infty}^{\infty} g(x)\phi(x-\mu) dx $, where $\phi$ is the standard normal density. In fact, $g=\widehat{\widehat{f}/\widehat{\phi}}$, where hat denotes the Fourier transform (might be missing some constants there, depending on convention). I'm trying to figure out if something similar can be derived for the standard deviation instead of the mean, i.e. given some $f(\sigma)$, find $g$ s.t. $f(\sigma)=\int_{-\infty}^{\infty} g(x)\frac{1}{\sigma}\phi\left(\frac{x}{\sigma}\right)dx=\mathbb{E}g(\sigma X)$ where $X\sim\mathcal{N}(0,1)$. If we pretend that $X$ is positive for a second, then writing $g(\sigma X)=g(\exp(\log(\sigma)+\log(X)))$ we have, defining $g_0(z)=g(\exp(z))$, that $\mathbb{E}g(\sigma X)=\mathbb{E}g_0(\log(\sigma)+\log(X))$ which (after finding the distribution of $\log(X)$) would look like a convolution, and the same approach as was used above for $\mu$ might work. Now, $X$ isn't positive, but it might be reasonable to say that $g$ should be symmetric anyway -- no reason for it not to be. So we can look for $g$ such that $f(\sigma)=\mathbb{E}g(\sigma |X|)$, and then the above approach will require computing the Fourier transform of the density of the log of the "folded" normal. But before I try diving into that -- does this sound familiar to anyone? Am I missing something that would make this either much easier or impossible? (Should I be asking this over at mathoverflow instead?) Bonus question -- multivariate version of the same problem, i.e. with $f(\Sigma)$. This is really the version I need to solve, and since the absolute value approach doesn't seem like it would work here, I'm trying to avoid it...


title: How much can searching for the best parameters improve SVM performance? (linear,polynomial,rbf)

text: Is there a rule of thumb? How much improvement in accuracy can I expect with extensive grid search? Can it be a HUGE difference? Does this differ with different kernels? ( linear, polynomial, RBF? ) I am wondering whether it's something I should focus on ( I am using SVMLight so it's really hard to do - as there doesn't appear to be any solution to do so )


title: Non-linear model fitting in many dimensions

text: I am interested in comparing a non-linear model with up to 12 parameters to many datasets. However, each instance of the model takes a significant amount of time to compute (~1 hour), so I am pre-computing instances of the model for various parameter values and then comparing these to all the different datasets. There are various ways to sample parameter space - so far, I've come across regular grids (impossible here), sparse grids + interpolation, Monte-Carlo random sampling, and there are probably others. Which approach would be optimal for a fixed amount of computing resources, and therefore a fixed number of model instances?


title: Ways to integrate user input into clustering algorithm

text: I want to integrate user input into clustering algorithm. So that users can control which docs should be in the same cluster. What are some ways to achieve this? Can metric learning do this?


title: What is a good analogy to illustrate the strengths of Hierarchical Bayesian Models?

text: I'm relatively new to bayesian statistics and have been using JAGS recently to build hierarchical bayesian models on different datasets. While I'm very satisfied of the results (compared to standard glm models), I need to explain to non-statisticians what the difference with standard statistical models is. Especially, I would like to illustrate why and when HBMs perform better than simpler models. An analogy would be useful, especially one that illustrates some key elements:   * the multiple levels of heterogeneity   * the need for more computations to fit the model   * the ability to extract more "signal" from the same data Note that the answer should really be an **analogy** enlightening to non-stats people, not an easy and nice-to-follow example.


title: How many of the biggest terms in $\sum_{i=1}^N |X_i|$ add up to half the total?

text: Consider $\sum_{i=1}^N |X_i|$ where $X_1, \ldots, X_N$ are i.i.d. and the CLT holds.   How many of the biggest terms add up to half the total sum ?   For example, 10 + 9 + 8 $\approx$ (10 + 9 + 8 $\dots$ + 1) / 2: 30 % of the terms reach about half the total. Define   $ \qquad\text{sumbiggest( j}; X_1 \dots X_N ) \equiv \text{sum of the j biggest of } |X_1| \dots |X_N| $   $ \qquad\text{halfsum}( N ) \equiv \text{the smallest j such that sumbiggest( j )} \approx \text{sumbiggest}( N ) / 2 . $ Is there a general asymptotic result for halfsum( $N, \mu, \sigma$ ) ?   A simple, intuitive derivation would be nice. (A little Monte Carlo suggests that sometimes halfsum( $N$ ) $\approx N$/4 or so;   that is, the biggest 1/4 of the $X_i$ add up to 1/2 the total.   I get 0.24 $N$ for halfnormal, 0.19 $N$ for exponential, for $N$ = 20, 50, 100.)


title: Is it ever appropriate to fit a regression on 3 points?

text: Lets assume we only have 2 variables, `x` and `y`, and our model (in R code) is `lm(y~x)`. Given that any 2 points can be used to form a line, I know it doesn't make much sense to try to fit this model if I only have 2 observations. What about the case where I have 3 observations? Would this model pass the laugh test?


title: Why do these power functions for difference in proportions give different answers?

text: I am comparing two power functions in R, a base function `stats::power.prop.test`, and the function `pwr::power.2p.test` in the `pwr` package. I would think that they would give the same answer, but they are slightly different. And the equations being solved in the source code are different. Why are they different, and in what circumstances should I use one over the other? If I understand the source code right, here is the power formula when the alternative hypothesis is that $p_1 < p_2$ for `pwr`: $$ \Phi \left( \Phi^{-1}\left( \alpha \right) - 2 \left( \arcsin \sqrt{p_1} - \arcsin \sqrt{p_2} \right) \sqrt{\frac{n}{2}}\right) $$ and for `stats`: $$ \Phi \left( \frac{ \sqrt{n} \left| p_1 - p_2 \right| + \Phi^{-1} \left( \alpha \right) \sqrt{ \left( p_1 + p_2 \right) \left( 1 - \left( p_1 + p_2 \right) \right)} } { \sqrt{p_1 \left( 1 - p_1 \right) + p_2 \left( 1 - p_2 \right)} }\right) $$ Below is some code comparing the output for the two functions.               require(pwr)     stats::power.prop.test(p1 = .50, p2 = .75, power = .90,        sig.level = 0.05, alternative = "two.sided")     ## gives n = 76.71     pwr::pwr.2p.test(h = ES.h(0.50,0.75), power = .90,        sig.level = 0.05, alternative = "two.sided")     ## gives n = 76.65


title: Tolerance interval for linear regression

text: I have two positively correlated variables, Q95 (independent, 95-percentile flow in cumecs from the previous summer) and LIFE (dependent, a score describing the health of some invertebrates). I would like to be able to predict, to some degree of confidence, what value of Q95 is required to attain a given LIFE. From what I have read online, what I need is a linear regression (which I'm quite happy with) with a tolerance interval (which is new to me). However, I can't find a good description (i.e., one I can understand) of how to do this. Ultimately, I'd like to be able to say, "A Q95 of x is required to achieve the target LIFE of y 90% of the time, with 95% certainty". Ideally I'd like to be able to perform these calculations in Excel. Any help would be greatly appreciated. Edit 14/12/11: I've since discovered how to do this in R, but it doesn't look like it's possible to do with Excel, due to a lack of functions for noncentral t and chisq distributions.               reg <- lm(y~x)     lower <- regtol.int(reg, side = 1, alpha = 0.05, P = 0.90)$'1-sided.lower'


title: How to test whether percent values differ according to a time condition?

text: I'm currently stuck on this data, trying to figure out what test is needed. I started thinking it needed an ANOVA test but on reading, figured it was a Chi Squared test. On more reading, that doesn't seem right either. The data I have is in a table as such:               Damage Group    0 minutes     2 minutes     11              40.17         28.38     12              10.08         9.47     13              9.30          12.12      and so on. The data in the minute columns represent percentages of the results over the damage group. (Damage Group has 9 levels, I haven't shown them all and there's 15 columns of minute data) I want to see if the different time distributions are significantly different to each other. What test is needed for this?


title: How best to analyze length of stay data in a hospital-based RCT?

text: I am interested in knowing whether or not there is a consensus about the optimal way to analyze hospital length of stay (LOS) data from a RCT. This is typically a very right-skewed distribution, whereby most patients are discharged within a few days to a week, but the rest of the patients have quite unpredictable (and sometimes quite lengthy) stays, which form the right tail of the distribution. Options for analysis include:   * t test (assumes normality which is not likely present)   * Mann Whitney U test   * logrank test   * Cox proportional hazards model conditioning on group allocation Do any of these methods have demonstrably higher power? Thanks.


title: Assumptions for robust linear regression, specifically ANCOVA type models

text: I have started looking into robust methods of data analysis. In particular, robust regression methods for ANCOVA type models using the _"lmRob"_ function from the _"robust"_ library in R. I'm having a hard time finding what assumptions I need to satisfy in order to use this. I've read that they can be useful in situations where there are outliers, heterogeneity of variance and for non-normal distributions, i.e., http://psych.colorado.edu/~willcutt/res_meth/Wilcox_1998.pdf (more recently http://www.springerlink.com/content/x761787633512440/). However, I'm unsure of how much departure robust measures will allow for and what other _traps_ might be awaiting me.   References and/or examples will be much appreciated. Data description:   sample size ~ 1000, differs substantially between groups   1st factor has 3 levels; 2nd factor has 3 levels   3 covariates, 1 of which is related to both factors, 1 related to only one factor The model I'm especially interested in examining is something like this:               lmr1 <- lmRob(dv ~ cov1 + cov2 + cov3 + fac1 * fac2, data)


title: How to detect a significant change in time series data due to a "policy" change?

text: I hope this is the right place to post this, I considered posting it on skeptics, but I figure they'd just say the study was statistically wrong. I'm curious about the flip side of the question which is how to do it right. On the website Quantified Self, the author posted the results of an experiment of some metric of output measured on himself over time and compared before and after abruptly stopping drinking coffee. The results were evaluated subjectively and the author believed that he had evidence that there was a change in the time series and it was related to the change in the policy (drinking coffee) What this reminds me of is models of the economy. We only have one economy (that we care about at the moment), so economists are often doing essentially n=1 experiments. The data is almost certainly autocorrelated over time because of this. The economists generally are watching, say the Fed, as it initiates a policy and trying to decide if the time series changed, potentially on account of the policy. What is the appropriate test to determine if the time series has increase or decreased based on the data? How much data would I need? What tools exist? My initial googling suggest Markov Switching Time Series Models, but not my googling skills are failing me at helping do anything with just the name of the technique.


title: Designing an inference engine for unknown items with known probabilities

text: I have done some AI projects involving inference, neural networks, etc., but I'm not sure what the right algorithm would be for this problem.   * Say there is a set of pills, where each type is visually unique. So, for example, a pain reliever might be a red pill and an allergy medication might be a blue pill. One is seeing all of these pills for the first time and they are unlabeled, so you don't know what the red pill is or the blue pill is for, other than that you know know all red pills are the same and that they are all different from blue pills.    * The way you are given new pills is that some assortment is put into a bag according to type (e.g., 2 allergy pills, 3 pain pills, 5 sleeping pills) and you take one from the bag at random. You are told what the assortment is, so when you pull a red pill from the bag you know (independent from anything else) that there was a 20% chance for it to have been allergy, 30% pain, etc.    * This process repeats itself with new independent bags and different percentages, so that you can make inferences about the function of the pill, such as knowing that you pulled a red from an {allergy (20%), pain(30%), sleeping(50%)} bag and another red from an {allergy(90%), pain (5%), antibiotic(5%)} bag means that it's either allergy or pain with allergy being more likely. Testing a pill to rule out possibilities of labels is possible, e.g., to completely rule out the possibility that a green pill is a pain pill. So, how would I set up some sort of inference engine to track this sort of thing? The output would be, for each of the colors, a percentage value for how likely it is to have a given label based on all of the information known at that point. The system would need to handle being given new random draw events with the new result and underlying percentages, along with events where labels are ruled out or confirmed outside of the random drawings. I could probably figure out how to do the individual calculations using conditional probability, but I want to make an efficient algorithm to do it automatically. If it matters, I would be planning on implementing this in Python.


title: When to drop a term from a regression model?

text: Could anyone advise if the following makes sense: I am dealing with an ordinary linear model with 4 predictors. I am in two minds whether to drop the least significant term. It's $p$-value is a little over 0.05. I have argued in favor of dropping it along these lines: Multiplying the estimate of this term by (for example) the interquartile range of the sample data for this variable, gives some meaning to the clinical effect that keeping this term has on the overall model. Since this number is very low, approximately equal to typical intra-day range of values that the variable can take when measuring it in a clinical setting, I see it as not clinically significant and could therefore be dropped to give a more parsimonious model, even though dropping it reduces the adjusted $R^2$ a little.


title: Subsample of a random sample: random sample?

text: Let's say you have a large random sample of soccer players in Europe but you are only interested in what happens in Spain. Could you reduce your sample to players in Spain and still call it a random sample (but of a different population)? If not, how would you call that subsample and which specific precautions should you take to be able to make inferences on the population of spanish soccer players? My feeling is that using that subsample would be fine as long as it is large enough, but maybe I am missing something.


title: How to describe a test of differences between correlation of the same two variables across two groups?

text: If I have two (Pearson) correlation coefficients between Variables A and B (for example .657 for males and .876 for females) and I want to test if there is a difference based on gender (as in this case), which of the following phrases should I use to describe this process?   * Test for equality of the correlation coefficients   * Test for the difference in the correlation coefficients I want to use one of the above as a row heading in a table that describes the above results, so I am looking for the one that is more theoretically grounded in statistical terminology.


title: Assessing statistical significance of a rare binary event in time series

text: I have a univariate time series which is an index of susceptibility to failure, and a binary variable which indicates whether a failure actually occurred in a given time window or not. I want to carry out a statistical test that would quantify the performance of the susceptibility index (e.g., compared to a null model where there is no correlation between the time series and the binary variable, or to another competing index). I am looking for methods to achieve this goal. I understand that if the failure is a very rare event obtaining a statistical significance may be hard, but even before that, my main problem is that the data points I sample from the time series are not independent (e.g., if the index is very high at a given time period, it's very likely high also in the next period). Therefore the length of the time window that I employ should be important. Any ideas?


title: Reference for the sum and difference of highly correlated variables being almost uncorrelated

text: In a paper I've written I model the random variables $X+Y$ and $X-Y$ rather than $X$ and $Y$ to effectively remove the problems that arise when $X$ and $Y$ are highly correlated and have equal variance (as they are in my application). The referees want me to give a reference. I could easily prove it, but being an application journal they prefer a reference to a simple mathematical derivation. Does anyone have any suggestions for a suitable reference? I thought there was something in Tukey's EDA book (1977) on sums and differences but I can't find it.


title: What is "reject inferencing" and how can it be used to increase the accuracy of a model?

text: Can anybody explain in detail:   1. What does reject inferencing mean?   2. How can it be used to increase accuracy of my model? I do have idea of reject inferencing in credit card application but struggling with the thought of using it to increase the accuracy of my model.