[R-sig-ME] lmer: LRT and mcmcpvalue for fixed effects

[Julie Marsh]

Dear lmer expeRts,

I wondered if anybody could help me reconcile these two test results ?

I have been fitting models of the form:

BASIC MODEL
fl.basic <- lmer(y.fl~(I(TIME-200) + I((TIME-200)^2))*SEX + V65_A +
                       (I(TIME-200)|STUDYNO ),
            data=fl.data[!is.na(fl.data$TIME),], method="ML",  
na.action=na.omit,
            control=list(maxIter=1000, niterEM=0, gradient=FALSE,  
msVerbose = 1))


BASIC MODEL PLUS BASIC GENETIC EFFECT
fl.gen1 <- lmer(y.fl ~(I(TIME-200) + I((TIME-200)^2))*SEX + SNP_dom + V65_A +
                       (I(TIME-200)|STUDYNO),
            data=fl.data[!is.na(fl.data$TIME),], method="ML",  
na.action=na.omit,
            control=list(maxIter=2000, niterEM=0, gradient=FALSE,  
msVerbose = 1))


I have compared the models using both (1) LRT

anova(fl.gen1, fl.basic)
0.5*(1-pchisq(2*(logLik(fl.gen1)-logLik(fl.basic)),0)) +  
0.5*(1-pchisq(2*(logLik(fl.gen1)-logLik(fl.basic)),1))

OUTPUT (1)
Models:
fl.basic: y.fl ~ (I(TIME - 200) + I((TIME - 200)^2)) * SEX + V65_A + (I(TIME -
                     200) | STUDYNO)
fl.gen  : y.fl ~ (I(TIME - 200) + I((TIME - 200)^2)) * SEX + SNP_dom +
fl.gen  :           V65_A + (I(TIME - 200) | STUDYNO)
              Df     AIC     BIC  logLik  Chisq Chi Df Pr(>Chisq)
fl.basic     10  4335.4  4396.1 -2157.7
fl.gen1      11  4291.4  4358.2 -2134.7 45.929      1  1.226e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 0.5*(1-pchisq(2*(logLik(fl.gen1)-logLik(fl.basic)),0)) +  
> 0.5*(1-pchisq(2*(logLik(fl.gen1)-logLik(fl.basic)),1))
[1] 6.131984e-12
attr(,"nobs")
[1] 3189


and (2) using Prof Bates' wonderful mcmcpvalue function

mcmc = mcmcsamp(fl.gen1, n = 100000)
mcmcpvalue(as.matrix(mcmc[ ,5]))

OUTPUT (2): main effect SNP_dom: p=0.0263


Given that I am using 2 different tests for two different hypotheses I  
still would have expected these p-values to be more similar. I am so  
sorry that I can't post the data and the lmer output but I am bound by  
confidentiality. <big sigh>  I understand completely if it is not  
possible to provide any help given this lack of further information.

I have eagerly read and re-read the rwiki help page  .........

http://wiki.r-project.org/rwiki/doku.php?id=guides:lmer-tests

....... but still am unable to explain why the results should be so  
different. Much as I would love to argue against the reliance on  
p-values I'm afraid I am a resigned pragmatist when it comes to trying  
to get anything published.  <sorry!>  Needless to say I will swamp the  
article with far more informative plots and CI's.

Any help would be very much appreciated.

kindest regards,  julie marsh.