Re: [R] Conservative "ANOVA tables" in lmer

From: Douglas Bates <bates_at_stat.wisc.edu>
Date: Wed 13 Sep 2006 - 15:02:43 GMT

On 9/13/06, Dimitris Rizopoulos <dimitris.rizopoulos@med.kuleuven.be> wrote:
>

> ----- Original Message -----
> From: "Manuel Morales" <Manuel.A.Morales@williams.edu>
> To: <A.Robinson@ms.unimelb.edu.au>
> Cc: "Douglas Bates" <bates@stat.wisc.edu>; "Manuel Morales"
> <Manuel.A.Morales@williams.edu>; <r-help@stat.math.ethz.ch>
> Sent: Wednesday, September 13, 2006 1:04 PM
> Subject: Re: [R] Conservative "ANOVA tables" in lmer
>
>
> > On Wed, 2006-09-13 at 08:04 +1000, Andrew Robinson wrote:
> >> On Tue, September 12, 2006 7:34 am, Manuel Morales wrote:
> >> > On Mon, 2006-09-11 at 11:43 -0500, Douglas Bates wrote:
> >> >> Having made that offer I think I will now withdraw it. Peter's
> >> >> example has convinced me that this is the wrong thing to do.
> >> >>
> >> >> I am encouraged by the fact that the results from mcmcsamp
> >> >> correspond
> >> >> closely to the correct theoretical results in the case that
> >> >> Peter
> >> >> described. I appreciate that some users will find it difficult
> >> >> to
> >> >> work with a MCMC sample (or to convince editors to accept
> >> >> results
> >> >> based on such a sample) but I think that these results indicate
> >> >> that
> >> >> it is better to go after the marginal distribution of the fixed
> >> >> effects estimates (which is what is being approximated by the
> >> >> MCMC
> >> >> sample - up to Bayesian/frequentist philosophical differences)
> >> >> than to
> >> >> use the conditional distribution and somehow try to adjust the
> >> >> reference distribution.
> >> >
> >> > Am I right that the MCMC sample can not be used, however, to
> >> > evaluate
> >> > the significance of parameter groups. For example, to assess the
> >> > significance of a three-level factor? Are there better
> >> > alternatives than
> >> > simply adjusting the CI for the number of factor levels
> >> > (1-alpha/levels).
> >>
> >> I wonder whether the likelihood ratio test would be suitable here?
> >> That
> >> seems to be supported. It just takes a little longer.
> >>
> >> > require(lme4)
> >> > data(sleepstudy)
> >> > fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
> >> > fm2 <- lmer(Reaction ~ Days + I(Days^2) + (Days|Subject),
> >> > sleepstudy)
> >> > anova(fm1, fm2)
> >>
> >> So, a brief overview of the popular inferential needs and solutions
> >> would
> >> then be:
> >>
> >> 1) Test the statistical significance of one or more fixed or random
> >> effects - fit a model with and a model without the terms, and use
> >> the LRT.
> >
> > I believe that the LRT is anti-conservative for fixed effects, as
> > described in Pinheiro and Bates companion book to NLME.
> >
>
> You have this effect if you're using REML, for ML I don't think there
> is any problem to use LRT between nested models with different
> fixed-effects structure.

There are two issues here: how the test statistic is evaluated and what reference distribution is used for the test statistic (i.e. how do you convert a value of the test statistic to a p-value). Manuel is addressing the latter issue. If you compare the difference of -2*logLik for the models to a chi-square with degrees of freedom determined by the difference in the number of parameters the test will be anti-conservative when the number of observations is small. The use of the chi-square as a reference distribution is based on asymptotic properties.

The other question is how does one evaluate the likelihood-ratio test statistic and that is the issue that Dimitris is addressing. The REML criterion is a modified likelihood and it is inappropriate to look at differences in the REML criterion when the models being compared have different fixed-effects specifications, or even a different parameterization of the fixed effects. However, the anova method for an lmer object does not use the REML criterion even when the model has been estimated by REML. It uses the profiled log-likelihood evaluated at the REML estimates of the relative variances of the random effects.  That's a complicated statement so let me break it down.

The optimization in lmer is done with respect to the elements of the variance-covariance matrix of the random effects relative to $\sigma^2$. Given these values the conditional estimates of the fixed-effects parameters and of $\sigma^2$ can be evaluated directly with some linear algebra. In the summary or show output of an lmer model there are two quantities called the MLdeviance and the REMLdeviance. Those are based on the same relative variances but different conditional estimates of $\sigma^2$ (and hence different estimates of the elements of the variance-covariance of the random effects). It turns out that there is very little difference in the value of the profiled log-likelihood at the ML estimates and at the REML estimates. This is not to say that the log-likelihood is similar at the two (complete) sets of estimates - it is the profiled log-likelihoods that are similar and these are what are used to create the likelihood ratio test statistic, even when the models have been fit by REML.

As I said, this is complicated and until I went to reply to this message I hadn't really sorted it out for myself. I knew the empirical result, which I had checked before releasing the code, but I hadn't worked out the details of why it occurred. This discussion has been very helpful to me at least.

If the explanation is too complicated, I suggest trying an example. Run the code that Andrew sent and look at the logLik that is quoted in the anova test. It is a log-likelihood not a REML criterion.

Now try

> logLik(fm1ML <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy, method = "ML"))
'log Lik.' -875.9697 (df=5)
> logLik(fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy), REML = FALSE)
'log Lik.' -875.993 (df=5)

and

> logLik(fm2ML <- lmer(Reaction ~ Days + I(Days^2) + (Days|Subject), sleepstudy, method = "ML"))
'log Lik.' -875.1408 (df=6)
> logLik(fm2 <- lmer(Reaction ~ Days + I(Days^2) + (Days|Subject), sleepstudy), REML = FALSE)
'log Lik.' -875.1588 (df=6)

Thanks to everyone for there comments on this issue.

> >> 2) Obtain confidence intervals for one or more fixed or random
> >> effects -
> >> use mcmcsamp
> >>
> >> Did I miss anything important? - What else would people like to do?
> >>
> >> Cheers
> >>
> >> Andrew
> >>
> >> Andrew Robinson
> >> Senior Lecturer in Statistics Tel:
> >> +61-3-8344-9763
> >> Department of Mathematics and Statistics Fax: +61-3-8344
> >> 4599
> >> University of Melbourne, VIC 3010 Australia
> >> Email: a.robinson@ms.unimelb.edu.au Website:
> >> http://www.ms.unimelb.edu.au
> >>
> >> ______________________________________________
> >> R-help@stat.math.ethz.ch mailing list
> >> https://stat.ethz.ch/mailman/listinfo/r-help
> >> PLEASE do read the posting guide
> >> http://www.R-project.org/posting-guide.html
> >> and provide commented, minimal, self-contained, reproducible code.
> >
> > ______________________________________________
> > R-help@stat.math.ethz.ch mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> > http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> >
>
>
> Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm
>
>



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Thu Sep 14 18:29:24 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Thu 14 Sep 2006 - 09:30:05 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.