# [R-sig-ME] How to compute deviance of a (g)lmer model under different values of the parameters ?

Emmanuel Charpentier charpent at bacbuc.dyndns.org
Tue May 5 09:37:07 CEST 2009

```Dear list,

I'm trying to implement multiple imputations analysis of datasets with
missing values for mixed models. My interest goes mainly to the fixed
effects.

As you probably know, the principle of such an analysis is to impute
"reasonable" values to missing data a few times (with "reasonable"
variability), analyze these few completed datasets and pooling the
analyses.

Schafer(1999) states simple methods for pooling estimates of the model
and their variances. These methods can be applied to mixed models (with
a few snags : you need to *ignore* the variation of random effects
between imputed datasets beyond what is reflected in fixed effects' SE,
you need to guesstimate the residuals DoF, which is something Douglas
Bates has expressed grave doubts about). As it stands, these methods
give results which, at first look, does not seem unreasonable, and might
be used as a first approximation. More on this later.

My current goal is to build a function to build a "pooled test" for
likelihood ratio of two models. Such a test has been proposed by Meng &
Rubin (Biometrika, 1992), and, according to a presentation by RA
Medeiros, has been implemented in Stata under the name milrtest.

I'm trying to implement such a pooled test in R. My snag is that the
estimate is based on the estimation, for each imputed dataset, of the
(log)likelihood ratio of each of these datasets under the hypotheses of
the coefficients having precisely the values obtained in the coefficient
pooling step.

So what I need is a (if possible elegant and/or fast) way to compute
log-likelihood of a given dataset under a model (already built) under
this hypothesis, up to a quantity supposed constant between models.

The logLik(model) function will happily give me LL(model|parameters
giving makimum LL). What I need is something like logLik(model, coefs)
giving me (up to an additive constant) LL(model|parameters=coefs).

Do you have any suggestions ? I can always sum squares of residuals
recomputed under this alternative hypothesis, but I'm not quite sure
that's enough for my purposes, especially if I plan to include the
random effect estimates later...

Emmanuel Charpentier

```