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

[Emmanuel Charpentier]

Le dimanche 10 mai 2009 à 09:17 -0500, Douglas Bates a écrit :
> Sorry to be so late in responding to this thread.  It has been a busy time.

Thank you for your answer !

> 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.

I think that mine is but one ...

> 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?

Profiling  comes to mind ... Wald tests ? A rough estimator for score
tests ? I'm no big fan of hypothesis testing, but (some) journal editors
live and die by it... and won't take a bootstrap for an answer, alas !

Again, thank you !

					Emmanuel Charpentier


> 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
> >>
> >
> > _______________________________________________
> > R-sig-mixed-models at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
>