# [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
>

```