Skip to content

lmer: LRT and mcmcpvalue for fixed effects

4 messages · Julie Marsh, Simon Blomberg, Rune Haubo +1 more

Original
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
[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.
Simon Blomberg
#
On Tue, 2008-07-15 at 13:43 +0800, Julie Marsh wrote:

            

      
      
      
    

        
Well, as Pinheiro and Bates say in their book (worth reading!), the LRT
for mixed effects models is anti-conservative. So your LRT p-value is
almost certainly too small. The posterior p-value might be more
accurate, if you accept the usual caveats re: priors and convergence
etc. Also, when calculating p-values by hand using pchisq, you should
probably use pchisq(..., lower.tail=FALSE) instead of 1-pchisq(...),
which is inaccurate. The log.p option might also be useful if you really
need to compare small probabilities. And why were you using pchisq with
0 df (which always == 1)? I don't understand that at all.

            

      
      
      
    

        

      
      
    

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Rune Haubo
    

    

      
      #
    

  


  

      
2008/7/15 Simon Blomberg <s.blomberg1 at uq.edu.au>:

            

      
      
      
    

        
Regarding the mixture of a chi-square distribution, this is more
appropriate when testing a single variance component. The ordinary
test with one df is conservative, since the test is on the boundary of
the parameter space. But Julie is not testing a variance component, so
the mixture is not appropriate here.

I can think of three reasons, that the p-values Julie obtains are
different in the likelihood ratio test and the posterior sampling.

1) the extensive use of control parameters indicates that convergence
may be an issue. If one or both the models have not reached
convergence, obviously the likelihood ratio test is based on wrong
likelihoods and will be misleading. I suppose the MCMC sampling will
also be inappropriate in this situation. The need for control
parameters could also indicate that some problems are related to the
data structure such as severe unbalance. Perhaps there is too little
information on some parameters and convergence is hard to achieve?

If the models converged nicely and problems with data or models are
unlikely, I would be inclined to trust the likelihood ratio test. It
is a test on one df with 3000+ observations, so as far as I know,
there should be no problems.

2) the MCMC sampling could be influenced by the priors or the chain
could get stuck in a specific region.

3) which version of lmer are you using? Douglas Bates recently posted
a message to this list reporting problems with MCMCsamp in the most
recent version of lme4.

/Rune

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Douglas Bates
    

    

      
      #
    

  


  

      
On Tue, Jul 15, 2008 at 3:33 AM, Rune Haubo <rhbc at imm.dtu.dk> wrote:

            

      
      
      
    

        
I think that the use of those control parameters is based on examples
in various postings to email lists or in papers.  As you mention
below, it is important to take note of the version of the lme4 package
being used.  Most of those control parameters are, to quote Ron
Ziegler, President Nixon's press secretary, "no longer operative".  In
a way, it's a pity.  Some of the best mathematics I have ever done was
the derivation of the analytic gradient and the ECME updates for the
profiled deviance and then I found that the use of these elegant
formulas actually slowed things down on large problems.  (There's a
lesson there for those who think that statistical computing research
consists of deriving a formula and letting others worry about the
trivial implementation details.)

I guess my question would be if it really is necessary to set the
maximum number of iterations so high.  (By the way, changing the
maximum number of iterations doesn't have an effect in the current
version of the lme4 package.  I will reimplement it in the next
release.)

            

      
      
      
    

        
I agree.  If convergence is smooth then the likelihood ratio statistic
of 46 on 1 degree of freedom is definitive.  I imagine that the
t-statistic for that coefficient is also large in magnitude (close to
+/- 7).  In those realms, assigning a numeric p-value is less
important than the information that "it is very, very small".

            

      
      
      
    

        
Yes, I would definitely want to look at the xyplot of the mcmcsamp
output, just to see if the chain managed to get itself stuck
somewhere.

            

      
      
      
    

        
Indeed - although I am more concerned about the sample of the variance
component parameters than about the sample from the distribution of
the fixed-effects parameters.  If the sample from the variance
component parameters does not have the correct properties that will,
of course, affect the sample from the distribution of the
fixed-effects parameters but I am more confident of the conditional
sampling of the fixed effects (it is just a multivariate normal - I
could still manage to botch that but it is less likely than is my
getting the more complicated distributions wrong) than of the
conditional sampling of the variance component parameters.

However, given the description of the problem the mcmcsamp would need
to have been created using lme4 version 0.99875-x (because versions
0.999375-y can't yet handle vector-valued random effects).

I think the bottom line is that the likelihood ratio test and the
t-statistic are definitive and there is no need to follow up with the
mcmcsamp approach.