[R-sig-ME] mcmcpvalue and contrasts
Ken Beath
kjbeath at kagi.com
Tue Feb 26 01:05:02 CET 2008
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
>
More information about the R-sig-mixed-models
mailing list