[R] Fitting mixed-effects models with lme with fixed error term variances

Douglas Bates bates at stat.wisc.edu
Tue Nov 21 21:06:33 CET 2006


On 11/21/06, Virgilio Gómez-Rubio <v.gomezrubio at imperial.ac.uk> wrote:
> Dear R users,

> I am writing to you because I have a few question on how to fix
> the error term variances in lme in the hope that you could help me. To
> my knowledge, the closest possibility is to fix the var-cov structure,
> but not the whole var-cov matrix. I found an old thread (a few years
> ago) about this, and it seems that the only alternative is to write the
> likelihood down and use optim or a similar function to get the MLE of
> the parameters of the model.

> I was also trying to compute Small Area Estimators using area level
> models, so that I have one observation per area. I tried to fit the
> following model:

> Yhat_i = X beta + u_i + e_i

> where Yhat is a direct estimator of the target variable, X, the area-
> level  covariates, u_i the random effects independent and distributed as
> N(0, sigma2_u) and e_i the 'error' terms, which are distributed as N(0,
> sigma2_e/n_i), where n_i is the sample size in area i.

> This model should be unidentifiable because we have one observation per
> group and we have to estimate two area level terms: u_i and e_i. But, to
> my surprise, I have been able to fit that model using lme and get
> sensible results.

> Am I missing something here? I have checked Pinheiro and Bates' book but
> I have not been able to find an answer to my questions.

It is more like we missed something there.  We should have checked
that the number of levels in each grouping factor is less than the
number of observations.  That check is done in the lmer function in
the lme4 package.

Because lme and lmer both optimize a profiled likelihood that
substitutes the conditional mle of the residual variance in the
expression for the likelihood, you can't easily make that value fixed.

Can you be more specific about which parameters you want to fix and
which you want to vary in the optimization?



More information about the R-help mailing list