[R-sig-ME] Fitting a ``trivial'' model with lmer().
Doran, Harold
HDoran at air.org
Wed Jul 30 16:45:32 CEST 2008
Can you post str(Rolf's_data)? If I understand correctly, y is a
temporally ordered vector of the measurements on each student, correct?
I further assume that tstnum is a factor that associates each element in
the vector y with a point in time. So, your data are something like:
Student y tstnum
1 y1 1
1 y2 2
1 y3 3
2 y1 1
2 y2 2
2 y3 3
...
If this is your data structure and your lmer syntax is:
lmer(y ~ tstnum + (0 + tstnum|stdnt), data=schooldat)
Then you will get the marginal effects for each time point as (mu +
alpha_i) and the variance/covariance matrix of the ranodm effects will
be very large and probably unstable. I don't know what you mean by the
"most general covariance matrix possible". As someone who models
longitudinal educational data on a daily basis, I would approach this
problem differently.
Again, I don't think I understand your data structure entirely correct,
but assume the data structure I pose above is what you have. I would
first start with a model where tstnum is an integer and not a factor in
the data frame. Because it seems you want the marginal effects for time,
I might do something like:
lmer(y ~ factor(tstnum) + (tstnum|stdnt), data=schooldat)
In this formulation, treating time in this manner accounts for the
serial correlation between scores and reduces the number of elements in
the variance/covariance matrix.
Now, I experimented quickly with the egsingle data set that is in
library(mlmRev). Here is how my thinking plays out:
> fm1 <- lmer(math ~ factor(year) + (year|schoolid), egsingle)
> summary(fm1)
Linear mixed-effects model fit by REML
Formula: math ~ factor(year) + (year | schoolid)
Data: egsingle
AIC BIC logLik MLdeviance REMLdeviance
20546 20608 -10264 20503 20528
Random effects:
Groups Name Variance Std.Dev. Corr
schoolid (Intercept) 0.1933133 0.439674
year 0.0098065 0.099028 0.555
Residual 0.9664253 0.983069
number of obs: 7230, groups: schoolid, 60
Fixed effects:
Estimate Std. Error t value
(Intercept) -3.20081 0.10074 -31.77
factor(year)-1.5 1.18359 0.09304 12.72
factor(year)-0.5 2.19556 0.09552 22.99
factor(year)0.5 2.88863 0.09991 28.91
factor(year)1.5 3.48153 0.10660 32.66
factor(year)2.5 4.29051 0.11451 37.47
Correlation of Fixed Effects:
(Intr) f()-1. f()-0. f()0.5 f()1.5
fctr(y)-1.5 -0.831
fctr(y)-0.5 -0.816 0.914
fctr(yr)0.5 -0.785 0.893 0.927
fctr(yr)1.5 -0.739 0.853 0.903 0.933
fctr(yr)2.5 -0.691 0.809 0.872 0.915 0.932
> fm2 <- lmer(math ~ factor(year) + (factor(year)|schoolid), egsingle)
Warning messages:
1: In .local(x, ..., value) :
Estimated variance-covariance for factor 'schoolid' is singular
2: In .local(x, ..., value) :
nlminb returned message false convergence (8)
The number of elements in the variance/covariance matrix in fm2 is large
and this model also has a singularity problem, similar to what you have
encountered. Let me stop here and let you add info or experiment with
what I pose.
> -----Original Message-----
> From: r-sig-mixed-models-bounces at r-project.org
> [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf
> Of Rolf Turner
> Sent: Tuesday, July 29, 2008 6:48 PM
> To: R-Sig-mixed-models
> Subject: [R-sig-ME] Fitting a ``trivial'' model with lmer().
>
>
> 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}}
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list