How to compute deviance of a (g)lmer model under different values of the parameters ?

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

Answering to myself, at least for archive users' sake :

I have a partial approximate solution, which seems to work (= gives "not
seemingly unreasonable answers") for *linear* mixed effect models.

For these models, the deviance (= -2 logLik(themodel)) is easily
computed as the sum on all observations of the squared residuals. This
residual r is the difference between the observed value y and the
estimate by the model y.hat.

If, as usual, we denote the fixed effects model matrix by X, beta the
unknown fixed effects coefficients, with estimator beta.hat, Z the
random effects model matrix and b the unknown random effects
coefficients with estimator b.hat, we all know that

y.hat = X%*%beta.hat + Z%*%b.hat (1)

Therefore, r = y-y.hat = y-y.hat = y-(X%*%beta.hat + Z%*%b.hat). (2)

The problem is that Z%*%b.hat might be *quite* complicated to rebuilt
from ranef(model) and Z.

Since we seek a solution allowing to compute the LL of the model under a
new estimator of beta (say beta.hat.new) *computed while ignoring the
random effects beyond what is reflected in beta.hat variances*, it seems
reasonable to ignore what modifications to Z%*%b.hat an adjustment of
the model to the new fixed effect estimators beta.hat.new would entail,
and use (2) with the original values to compute the new "residuals". To
be precise :

a) Zbh=y.hat-X%*%beta.hat (computed from the original model)
b) r.new=y-(X%*%beta.hat.new+Zbh) (in the supposed new model)

which can be "onelined", of course... at the expense of clarity.
Computing deviance an logLik is trivial afterwards...

Two questions :
? ? ? ?a) What do you think ? Am I crazy like a fox ? (I wouldn't mind being
crazy like a Fox(*) ... :)
? ? ? ?b) Could such a trick be applied to *generalized* linear models ?
Computing the deviance from the residuals doesn't seem as easy as for
linear models...

Sincerely,

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?Emmanuel Charpentier

(*) with apologies to John F, who has probably read|heard it a lot of
times...

Le mardi 05 mai 2009 ? 09:37 +0200, Emmanuel Charpentier a ?crit :

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

Sorry to be so late in responding to this thread.  It has been a busy time.

This thread and others like it have convinced me that I should make
available an evaluation of the deviance of a fitted or an intermediate
lmer model at arbitrary values of the parameters.  For computational
efficiency the deviance and the REML criterion of a linear mixed model
is always optimized as the profiled deviance or the profiled REML
criterion.  This is an important characteristic of the lme4 package
because it provides much faster and robust model fits.  Nevertheless,
there are circumstances in which it is desirable to evaluate the
deviance or the REML criterion at arbitrary values of the parameters.

In section 2 of Bates and DebRoy (2004, J Multivariate Analysis, 1-17)
we derive general expressions for the deviance (eqn. 8) and for the
REML criterion (eqn. 18) as functions of all the parameters in the
model.  (You do need to be careful about the REML criterion in that it
does not depend on the values of the fixed-effects parameters.  For
clarity it would be best to stick to the likelihood-based criteria
like the deviance.)  For parameter estimation these are collapsed to
the profiled deviance (eqn. 10) and the profiled REML criterion (eqn.
15) but it would be possible to allow for arbitrary values of the
fixed-effects parameters and the residual variance to be included in
the evaluation.  It turns out that arbitrary fixed-effects parameters
are relatively easy.  Arbitrary residual variance is a bit more
complicated because of the way that the variances and covariances of
the random effects are expressed in the model.  They are expressed
relative to the residual variance and in a parameterization with
simple non-negativity constraints.  Going from arbitrary
variance-covariance parameters to this particular parameterization
would be somewhat challenging but not insurmountable.

Is the evaluation of the deviance for arbitrary values of the
fixed-effects parameters but with the conditional estimate of the
residual variance of use by itself?

On Tue, May 5, 2009 at 4:02 PM, Emmanuel Charpentier
<charpent at bacbuc.dyndns.org> wrote:

Answering to myself, at least for archive users' sake :

I have a partial approximate solution, which seems to work (= gives "not
seemingly unreasonable answers") for *linear* mixed effect models.

For these models, the deviance (= -2 logLik(themodel)) is easily
computed as the sum on all observations of the squared residuals. This
residual r is the difference between the observed value y and the
estimate by the model y.hat.

If, as usual, we denote the fixed effects model matrix by X, beta the
unknown fixed effects coefficients, with estimator beta.hat, Z the
random effects model matrix and b the unknown random effects
coefficients with estimator b.hat, we all know that

y.hat = X%*%beta.hat + Z%*%b.hat (1)

Therefore, r = y-y.hat = y-y.hat = y-(X%*%beta.hat + Z%*%b.hat). (2)

The problem is that Z%*%b.hat might be *quite* complicated to rebuilt
from ranef(model) and Z.

Since we seek a solution allowing to compute the LL of the model under a
new estimator of beta (say beta.hat.new) *computed while ignoring the
random effects beyond what is reflected in beta.hat variances*, it seems
reasonable to ignore what modifications to Z%*%b.hat an adjustment of
the model to the new fixed effect estimators beta.hat.new would entail,
and use (2) with the original values to compute the new "residuals". To
be precise :

a) Zbh=y.hat-X%*%beta.hat (computed from the original model)
b) r.new=y-(X%*%beta.hat.new+Zbh) (in the supposed new model)

which can be "onelined", of course... at the expense of clarity.
Computing deviance an logLik is trivial afterwards...

Two questions :
       a) What do you think ? Am I crazy like a fox ? (I wouldn't mind being
crazy like a Fox(*) ... :)
       b) Could such a trick be applied to *generalized* linear models ?
Computing the deviance from the residuals doesn't seem as easy as for
linear models...

Sincerely,

                                       Emmanuel Charpentier

(*) with apologies to John F, who has probably read|heard it a lot of
times...

Le mardi 05 mai 2009 ? 09:37 +0200, Emmanuel Charpentier a ?crit :

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