R help - Nov 2006 - Fitting mixed-effects models with lme with fixed error term variances

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.

Many thanks for your help. Best regards,

Virgilio

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?

Douglas Bates <bates <at> stat.wisc.edu> writes:
...> 
> Can you be more specific about which parameters you want to fix and
> which you want to vary in the optimization?
It would be nice to have the ability to fix all variances i.e. variances of
random effects.

Gregor

Hello again,

And excuse me for not replying in real time. I am subscribed in digest
mode and the data that generated my question are in a computer not
connected to the Internet.

Douglas Bates replied to my original mail and he requested the
following:

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


My model is

Y_i = X beta + u_i + e_i

where u_i ~ N(0, sigma2_u I) are the random effects and 
e_i ~ N(0, SIGMA) the error terms. What I would like to fix 
is matrix SIGMA, which, in addition, is diagonal. Hence, I only want to  
estimate beta and sigma2_u. The reason why I want to fix SIGMA is because,
otherwise, the model is unidentifiable (I have one observation per group)
and I can get an estimate of SIGMA by other means.

Is that possible? If not, could I use any trick to do so?

Many thanks.

Best regards,

Virgilio