[R-sig-ME] Fitting a ``trivial'' model with lmer().
Rolf Turner
r.turner at auckland.ac.nz
Wed Jul 30 00:48:20 CEST 2008
I sent the message, to be found below, to Doug Bates just now.
Since he is flat out preparing a number of workshops, he suggested
that I post to this SIG to see if I could get some enlightenment.
Ergo, here goes:
I wrote (to the general R-help list; 28 July 2008):
<snip>
>> How can I (is there any way that I can)
>> tell lmer() to fit the most general possible covariance structure?
>>
and Doug Bates responded:
> It sounds like you want a model formula of
>
> lmer(y ~ tstnum + (0 + tstnum|stdnt), data=schooldat)
>
> but that model will have 21 variance-covariance terms to estimate (22
> if you count the residual variance but that one gets profiled out of
> the optimization). I would not be surprised if the estimated
> variance-covariance matrix for the random effects turns out to be
> singular.
>
I tried the above formulation and it (as I posted yesterday) gives
a warning about a false convergence message from nlminb. I think
there are indeed singularity problems, or rather problems with a
certain indeterminateness.
I think that the forgoing formulation translates into ``mathematical
expressions'' as
y_ij = mu + alpha_i + s_j + e_ij
(with treatment ``contrasts'' giving the constraint that alpha_1 = 0)
where alpha_i corresponds to the i-th time, s_j corresponds to the j-th
student, and the e_ij are (Gaussian) errors with mean 0 and a covariance
structure
Cov(e_ij, e_{i',j'}) = 0 if j != j'
= gamma_{i,i'} if j = j'
Translated into my ``trivial model'' thinking, this can be expressed as
Y_j = MU + s_j + E_j
where the Y_j and E_j are 6-dimensional vectors and the addition in the
forgoing follows the S/R convention whereby the scalar s_j is added to
each component of the vectors. And of course MU is the 6 dimensional
vector of mean values.
The vectors Y_j are independent N(MU,Sigma) where Sigma = Gamma +
sigma^2
where Gamma = [gamma_ij] and sigma^2 is the variance of the
(independent,
mean 0) s_j. I believe that sigma^2 is the residual variance in the
lmer() formulation of the model.
The covariance matrix Sigma is well-determined/estimable but Gamma
and sigma^2
are not --- we can ``sensibly'' have Gamma = Sigma - sigma^2 for any
value
of sigma^2 such that Sigma - sigma^2 is positive definite. In general
there
will be a wide range of such values of sigma^2.
It seems to me that the s_j are redundant in the model; we just need:
y_ij = mu + alpha_i + e_ij
Is there a way to formulate the foregoing model in lmer()? Effectively,
from one p.o.v., one wants to constrain the sigma^2 variance
component to
be 0.
Of course one doesn't need lmer() to fit such a model. I just want to
get some understanding of the modelling syntax in lmer() and perhaps
be able to fit more complicated models which are related to the
foregoing
``trivial'' one.
Thanks.
cheers,
Rolf Turner
######################################################################
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}
More information about the R-sig-mixed-models
mailing list