[R] generalized linear mixed model by ML

[Douglas Bates]

On 12/15/05, Ben Bolker <bolker at ufl.edu> wrote:
> Abderrahim Oulhaj <abderrahim.oulhaj <at> pharmacology.oxford.ac.uk> writes:
>
> >
> > Dear All,
> >
> > I wonder if there is a way to fit a generalized linear mixed models (for
> repeated  binomial data)  via a direct
> > Maximum Likelihood Approach. The "glmm" in the "repeated" package (Lindsey),
> the "glmmPQL" in the
> > "MASS" package (Ripley) and "glmmGIBBS"  (Myle and Calyton) are not using the
> full maximum likelihood as I
> > understand. The "glmmML" of Brostrom uses the "full maximum likelihood" by
> approximating the integral
> > via  Gauss- Hermite  quadrature. However, glmmML is only valid for the random
> intercept model and the
> > binomial family must be represented only as  binary data. Does the lmer do the
> work?

The expression "full maximum likelihood" is confusing.  In the
multilevel modeling literature it is sometimes used in contrast to
restricted maximum likelihood or REML estimation but in that context
it represents what statisticians would simply term "maximum
likelihood".

I think you are referring to "exact maximum likelihood" where the
criterion being optimized is the exact log-likelihood associated with
the parameters.  I don't believe that can be done in the the
evaluation of the (marginal) log-likelihood associated with the
parameters to be estimated involves an intractable integral.  Various
methods for estimating parameters in this model involve approximations
to this integral.  Even the Bayesian methods such as MCMC methods are
using an approximation to this integral (in that case a stochastic
approximation).

So, as far as I know, there are no exact maximum likelihood estimation
methods for this model except in some special circumstances.
>
>   Hmmm.  I will be interested to hear what others have to say on
> this topic.
>
> * lmer() in the lme4 package (new version of nlme) can in
> fact do GLMMs with a choice of different
> integration methods (PQL is the default but not the only choice).
>
> * GLMMGibbs [sic] was actually
> using a full likelihood rather than an approximation, but was
> a Bayesian rather than a ML approach [GLMMGibbs is now in the
> "Devel" section of CRAN, apparently because of various unresolved
> compilation/installation problems.]
>
>   If you want to fit temporal correlation, as well as
> individual random effects, you may be out of luck: a GEE
> model will probably be your best bet in that case.  When I asked
> about the possibility of incorporating temporal/spatial correlation
> structures like those in nlme into lme4, Doug Bates said that he
> wanted to work first on getting the basic framework of the package
> really solid [can't blame him at all, and of course honor and
> glory to him for putting so much work into these tools in the first place]
>

Wow - honor and glory - thanks.

The lmer function in the currently released Matrix package (0.99-3)
does provide PQL, Laplace and Adaptive Gauss-Hermite Quadrature (AGQ)
methods for estimation of the parameters in a generalized linear mixed
model.  PQL is the fastest but least accurate method.  The Laplacian
approximation is actually a special case of the AGQ method where only
one quadrature point per random effect is used.  The AGQ method is the
most accurate and in many ways is the "gold standard".  It has a
tuning parameter which is the number of quadrature points per random
effect.  Usually this is an odd number (because you always evaluate at
the conditional modes of the random effects which, for odd numbers, is
one of the quadrature points).  The lmer function allows up to 11
quadrature points per random effect.  It would be easier to add the
capability for more but adding a whole lot more quadrature points
doesn't buy you that much.  It is definitely a case of diminishing
returns.

The released version of the Matrix package uses a partitioned
representation of the mixed model structures and sometimes encounters
errors with complex models based on multiple grouping factors that are
non-nested.  A couple of months ago I started developing an
alternative version (for those keeping count, yet this is the sixth
implementation - the one I said I wouldn't do) based on what is called
a supernodal Cholesky factorization for which Tim Davis recently
released code.  That is doing fine on the linear mixed models and PQL
for the generalized linear mixed model.  I'm adding Laplace now and
feel I am close to getting that working (although, as Dave Balsinger
said, "Programmers spend 95% of their time being `95% done'").  AGQ is
more difficult especially for multiple grouping factors.  At present
AGQ is only available for a univariate random effect associated with a
single grouping factor.  Once I get that implemented I will release a
new Matrix package based on the supernodal factorization.

>    good luck,
>     Ben Bolker
>
> ______________________________________________
> R-help at 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
>