[R-sig-ME] mcmcpvalue and contrasts

[Ken Beath]

On 26/02/2008, at 9:42 AM, Hank Stevens wrote:

> Hi Folks,
> I wanted to double check that my intuition makes sense.
>
> Examples of mcmcpvalue that I have seen use treatment "contrast"  
> coding.
> However, in more complex designs, testing overall effects of a factor
> might be better done with other contrasts, such as sum or Helmert
> contrasts.
>
> My Contention:
> Different contrasts test different hypothesis, and therefore result in
> different P-values. This consequence of contrasts differs from
> analysis of variance, as in anova( lm(Y ~ X1*X2) ).
>
> *** This is right, isn't it? ***
>

The main problem is testing for a main effect in the presence of  
interaction. While it looks like it gives sensible results in some  
cases like balanced ANOVA, they really aren't sensible and the effect  
of parameterisation in other cases makes that clear.

The difference for the interaction is probably just sampling  
variation, increasing samples fixes this.

Ken



> I provide a self-contained example.
>
> library(lme4) # v. 0.99875-9
> library(coda)
> library(SASmixed)
>
> mcmcpvalue <- function(samp) {
>    ##From Bates pers comm. September 14, 2006 2:59:23 PM EDT
>    ## elementary version that creates an empirical p-value for the
>    ## hypothesis that the columns of samp have mean zero versus a
>    ## general multivariate distribution with elliptical contours.
> # samp <- rnorm(10000, m=3)
>    ## differences from the mean standardized by the observed
>    ## variance-covariance factor
>    std <- backsolve(chol(var(samp)),
>                     cbind(0, t(samp)) - colMeans(samp),
>                     transpose = TRUE)
>    sqdist <- colSums(std * std)
>    sum(sqdist[-1] > sqdist[1])/nrow(samp)
> }
>
>
> options(contrasts=c("contr.treatment","contr.poly"))
> print(fm1 <- lmer(resistance ~ ET * position + (1|Grp),  
> Semiconductor),
>       corr=F)
> c1 <- mcmcsamp(fm1, 50000)
>
> options(contrasts=c("contr.helmert","contr.poly"))
> print(fm2 <- lmer(resistance ~ ET * position + (1|Grp),  
> Semiconductor),
>       corr=F)
> c2 <- mcmcsamp(fm2, 50000)
>
>  mcmcpvalue(c1[, 2:4 ])
> [1] 0.32458
>> mcmcpvalue(c2[, 2:4 ])
> [1] 0.0249
>> mcmcpvalue(c1[, 5:7 ])
> [1] 0.45376
>> mcmcpvalue(c2[, 5:7 ])
> [1] 0.16302
>> mcmcpvalue(c1[, 8:16 ])
> [1] 0.6373
>> mcmcpvalue(c2[, 8:16 ])
> [1] 0.88808
>>
>> sessionInfo()
> R version 2.6.2 (2008-02-08)
> i386-apple-darwin8.10.1
>
> locale:
> C
>
> attached base packages:
> [1] stats     graphics  utils     datasets  grDevices methods
> [7] base
>
> other attached packages:
> [1] SASmixed_0.4-2    xtable_1.5-2      reshape_0.8.0
> [4] coda_0.13-1       lme4_0.99875-9    Matrix_0.999375-4
> [7] lattice_0.17-4    CarbonEL_0.1-4
>
> loaded via a namespace (and not attached):
> [1] grid_2.6.2
>>
>
> Dr. Hank Stevens, Assoicate Professor
> 338 Pearson Hall
> Botany Department
> Miami University
> Oxford, OH 45056
>
> Office: (513) 529-4206
> Lab: (513) 529-4262
> FAX: (513) 529-4243
> http://www.cas.muohio.edu/~stevenmh/
> http://www.cas.muohio.edu/ecology
> http://www.muohio.edu/botany/
>
> "If the stars should appear one night in a thousand years, how would  
> men
> believe and adore." -Ralph Waldo Emerson, writer and philosopher
> (1803-1882)
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>