R help - May 2004 - GLMMs & LMEs: dispersion parameters, fixed variances, design matrices

Three related questions on LMEs and GLMMs in R:

(1) Is there a way to fix the dispersion parameter (at 1) in either glmmPQL
(MASS) or GLMM (lme4)?

Note: lme does not let you fix any variances in advance (presumably because it
wants to "profile out" an overall sigma^2 parameter) and glmmPQL
repeatedly calls lme, so I couldn't see how glmmPQL would be able to fix the
dispersion parameter. The section on glmmPQL in V&R4 says that the default
is to estimate the dispersion parameter, but didn't seem to say how to
change the default.

(2) Is there a way to tell lme (either in nlme or lme4) to just use a specified
design matrix Z for the random effects, rather than constructing one itself from
factors? Sometimes I would really like to use my own funny-looking Z matrix
(e.g. with non-integer coefficients), and even with contrasts() I haven't
managed to do this.

(3) Are there any plans to allow some variances to be fixed in lme? It would be
useful e.g. for meta-analysis (and indeed for glmms with fixed dispersion).

Note: it has occurred to me that lme can possibly be tricked into fixing the
measurement error variance (i.e. var[y|b] where b is the random effects and y
the observed data) at some specified value e.g. 1 by adding two
pseudo-observations at +/-1, with all zeros in the corresponding rows of the X
and Z matrices, and with huge weights. Then sum( w*(y-E[y|b,params])^2) / sum(w)
will be approximately 1, and any attempt to change the estimate of sigma^2 away
from 1 will be "deterred" by a large penalty. Similar tricks might be
possible for fixing other variances. However this approach is not nice and
perhaps might cause computational problems-- and I haven't actually tried it
yet.

Thanks for any help

Mark

*******************************

Mark Bravington
CSIRO (CMIS)
PO Box 1538
Castray Esplanade
Hobart
TAS 7001

phone (61) 3 6232 5118
fax (61) 3 6232 5012
Mark.Bravington at csiro.au

On Thu, 13 May 2004 Mark.Bravington at csiro.au wrote:
> Three related questions on LMEs and GLMMs in R:
> 
> (1) Is there a way to fix the dispersion parameter (at 1) in either
> glmmPQL (MASS) or GLMM (lme4)?
not glmmPQL in R (can be done in S-PLUS).
> Note: lme does not let you fix any variances in advance (presumably
> because it wants to "profile out" an overall sigma^2 parameter)
and
> glmmPQL repeatedly calls lme, so I couldn't see how glmmPQL would be
> able to fix the dispersion parameter. The section on glmmPQL in V&R4
> says that the default is to estimate the dispersion parameter, but
> didn't seem to say how to change the default.
It's done in the same way as for lme via the control parameter (that is, 
not at all in R).
> (2) Is there a way to tell lme (either in nlme or lme4) to just use a
> specified design matrix Z for the random effects, rather than
> constructing one itself from factors? Sometimes I would really like to
> use my own funny-looking Z matrix (e.g. with non-integer coefficients),
> and even with contrasts() I haven't managed to do this.
> 
> (3) Are there any plans to allow some variances to be fixed in lme? It
> would be useful e.g. for meta-analysis (and indeed for glmms with fixed
> dispersion).
It has been possible for a while in S-PLUS.
> Note: it has occurred to me that lme can possibly be tricked into fixing
> the measurement error variance (i.e. var[y|b] where b is the random
> effects and y the observed data) at some specified value e.g. 1 by
> adding two pseudo-observations at +/-1, with all zeros in the
> corresponding rows of the X and Z matrices, and with huge weights. Then
> sum( w*(y-E[y|b,params])^2) / sum(w) will be approximately 1, and any
> attempt to change the estimate of sigma^2 away from 1 will be
"deterred"
> by a large penalty. Similar tricks might be possible for fixing other
> variances. However this approach is not nice and perhaps might cause
> computational problems-- and I haven't actually tried it yet.
-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

Hi Mark,
> (2) Is there a way to tell lme (either in nlme or lme4) to just use a
> specified design matrix Z for the random effects, rather than
> constructing one itself from factors? Sometimes I would really like to
> use my own funny-looking Z matrix (e.g. with non-integer coefficients),
> and even with contrasts() I haven't managed to do this.
Here's a simple example showing a `standard' call to lme, and a
less standard one doing the same thing with a specified Z....

library(nlme)
data(Rail)
lme(travel~1,Rail,~1|Rail) # simple r.e. model example

## now repeat fit with home-made Z...
Rail$all<-factor(rep(1,18)) # all one group
## create a simple Z ...
Rail$Z<-model.matrix(~factor(Rail$Rail,ordered=FALSE)-1)
## use pdIdent to force r.e.s to be independent ...
lme(travel~1,Rail,list(all=pdIdent(~Z)))


best,
Simon

<Mark.Bravington at csiro.au> writes:
> Three related questions on LMEs and GLMMs in R:
> 
> (1) Is there a way to fix the dispersion parameter (at 1) in either
> glmmPQL (MASS) or GLMM (lme4)?
> 
> Note: lme does not let you  fix any variances in advance (presumably
> because it wants to "profile out"  an overall sigma^2 parameter)
and
> glmmPQL repeatedly calls lme, so I couldn't see how glmmPQL would be
> able to fix the dispersion parameter. The section on glmmPQL in V&R4
> says that the  default is to estimate  the dispersion parameter, but
> didn't seem to say how to change the default.
At the core of the lme calculations is the solution of a penalized
least squares problem defined by the relative dispersion matrix of the
random effects and the model matrices for the random effects and the
fixed effects.  In versions 0.6-1 and later of the lme4 package (the
first release candidate is available from my web site
http://www.stat.wisc.edu/~bates/) the components of the log-likelihood
or the REML criterion are available as the devComp slot of the S4
object that represents the model and that is used to solve the
penalized least squares problem.  If, using these components, you can
write the log-likelihood for the model that you wish to fit then you
can give it to an optimizer such as optim or nlm to fit.

In the notation of Bates and DebRoy (2004), "Linear mixed models and
penalized least squares" (to appear in J. of Multivariate Analysis,
available in preprint form from my web site), the components are 

log(|Z'Z + \Omega|), log(|\Omega|), log(|R_{XX}|^2), and log(r_{yy}^2)

The C code that uses these to evaluate the deviance form of the
profiled log-likelihood criterion or the profiled REML criterion from
these components is in src/ssclme.c from the Matrix package.

Modifying the criteria for a fixed dispersion parameter may be trivial
or it may not.
> (3) Are there any plans to allow some variances to be fixed in lme? 
> It would be useful e.g. for meta-analysis (and indeed for glmms with
> fixed dispersion).
The method = 'Laplacian' version of the GLMM function fixes the
dispersion parameter in those families where it should be fixed.  As
we continue to develop lme4 we will provide a further enhancement
using an adaptive Gauss-Hermite evalution of the log-likelihood for
GLMMs that will also have this property.