Many thanks to Reinhold Kliegl, my sole respondent.

I have thought a bit more about this since my original post. It seems 
that factors and interactions are a bit harder to handle than continuous 
covariates in hierarchical situations and it is easy to overdo the 
number of random effects. Also it may have some advantages to take 
direct control over the model matrices that you are setting up (as 
Reinhold has done) so that you can specify which particular contrasts 
you want to be random and which fixed.

I also experienced the convergence failures that Reinhold speaks of, 
except when I simulated an implausibly large number of sites.

Reinhold's suggestion

mod = lmer(y ~ year/season + season/month + (1 | site) + (0 + year/season | site) )

does come close to answering my original question but the large number of parameters in this model means that it is not suitable for real data. I suppose for many purposes oversimplifying the model to

y ~ year/season + season/month + (1 | site) 


would be adequate. The intercept is an obvious thing to make random but after that I guess that a bit of subject area knowledge is needed to pick a small number of contrasts to be made random.

Murray Jorgensen


Reinhold Kliegl wrote:
> Two more things:
> (1) I accidentally left out the main effects for "year" in the
> fixed-effect part of the lmer model. So the call should actually be:
>   
>> mod.2 = lmer(y ~ year + ys + sm + (1 | site) + (0 + ys | site) )
>>     
> or for direct estimates of (almost all) cell means
>   
>> mod.3 = lmer(y ~ 0 + ys + sm + (1 | site) + (0 + ys | site) )
>>     
> In real data, year should probably be specified with trends, e.g., as
> poly(year, 2).
>
> (2)  I think , in principle, replacing
>   
>> mod = lmer(y ~ year/season + season/month + (1 + year/season | site) )
>>     
> with
>   
>> mod = lmer(y ~ year/season + season/month + (1 | site) + (0 + year/season | site) )
>>     
> should also work for the sigmas in your script. I did run into a
> convergence problem with this specification.
> Also, the fixed effects were not very transparent to me in this specification.
>
> Best
> Reinhold Kliegl
>
>
> On Thu, May 14, 2009 at 11:38 PM, Reinhold Kliegl
> <reinhold.kliegl@gmail.com> wrote:
>   
>> The following parameterization recovers both the fixed effects and
>> also sigma_mu, sigma_beta, and sigma. The correlation does not look
>> bad either. There are probably also ways of forcing some of the
>> irrelevant correlations to zero.
>>
>> Reinhold Kliegl
>>
>> ############################################
>>
>> ....
>>
>> P = rmvnorm(nsite,mean=theta,sigma=V, method="chol")
>> ...
>>
>> # specifying seasons nested within year
>> ys  <- factor(paste(year, season, sep="_"))
>> cmat.ys <- matrix(c(          -1,  1, rep(0,8),
>>                            rep(0,2), -1,  1, rep(0,6),
>>                            rep(0,4), -1,  1, rep(0,4),
>>                            rep(0,6), -1,  1, rep(0,2),
>>                            rep(0,8), -1,  1             ),  10,  5)
>> colnames(cmat.ys) <- c(".s2-s1|y1", ".s2-s1|y2", ".s2-s1|y3",
>> ".s2-s1|y4",".s2-s1|y5")
>> contrasts(ys, 5) <- cmat.ys
>>
>> # specifying months nested within season
>> # contrasting months 2 vs 1 and 3 vs 2 within season
>> sm  <- factor(paste(season, month, sep="_"))
>> cmat.sm <- matrix(c(-2/3,  1/3,  1/3,    0,     0,     0,
>>                                  -1/3, -1/3,  2/3,    0,     0,     0,
>>                                      0,     0,     0, -2/3,  1/3,  1/3,
>>                                      0,     0,     0, -1/3, -1/3,  2/3  ),  6,  4)
>> colnames(cmat.sm) <- c(".m2-m1|s1", ".m3-m2|s1", ".m2-m1|s2", ".m3-m2|s2")
>> contrasts(sm, 4) <- cmat.sm
>>
>> mod.2 = lmer(y ~ ys + sm + (1 | site) + (0 + ys | site) )
>> print(summary(mod.2), cor=FALSE)
>>
>> Linear mixed model fit by REML
>> Formula: y ~ ys + sm + (1 | site) + (0 + ys | site)
>>  AIC  BIC logLik deviance REMLdev
>>  1377 1745 -621.4     1200    1243
>> Random effects:
>>  Groups   Name  Variance  Std.Dev. Corr
>>  site   (Intercept) 7.894963 2.80980
>>  site   ys1_1        7.401678 2.72060
>>          ys1_2        4.987042 2.23317   0.691
>>          ys2_1        3.228583 1.79683   0.126  0.256
>>          ys2_2        5.913934 2.43186  -0.177 -0.205  0.419
>>          ys3_1      11.414706 3.37857  -0.239 -0.242  0.018  0.651
>>          ys3_2      10.336715 3.21508   0.109  0.070  0.244  0.426
>> 0.661
>>          ys4_1        4.113868 2.02827   0.392  0.180 -0.115 -0.159
>> -0.015  0.521
>>          ys4_2        4.801364 2.19120   0.391  0.108 -0.268 -0.215
>> -0.061  0.217  0.503
>>          ys5_1        2.867584 1.69339   0.186 -0.073 -0.356 -0.190
>> 0.107  0.373  0.097  0.268
>>          ys5_2        3.745493 1.93533  -0.002 -0.107 -0.049  0.075
>> 0.174  0.224 -0.069  0.046  0.241
>>  Residual             0.010422 0.10209
>> Number of obs: 1800, groups: site, 60
>>
>> Fixed effects:
>>             Estimate Std. Error t value
>> (Intercept) 13.892119   0.375663   36.98
>> ys.s2-s1|y1 -4.447045   0.122316  -36.36
>> ys.s2-s1|y2 -4.164777   0.120216  -34.64
>> ys.s2-s1|y3 -3.161133   0.155749  -20.30
>> ys.s2-s1|y4 -3.047036   0.129606  -23.51
>> ys.s2-s1|y5 -3.863464   0.136543  -28.29
>> sm.m2-m1|s1  1.988407   0.008335  238.55
>> sm.m3-m2|s1 -0.981789   0.008335 -117.79
>> sm.m2-m1|s2 -1.990586   0.008335 -238.81
>> sm.m3-m2|s2 -0.008253   0.008335   -0.99
>>
>> ###########################################
>>
>>
>> On Wed, May 13, 2009 at 6:19 AM, Murray Jorgensen
>> <maj@stats.waikato.ac.nz> wrote:
>>     
>>> I am trying to help a friend towards an appropriate lmer() analysis of her
>>> data and help myself understand how to use lmer syntax in the right way.
>>>
>>> She says "...here is what I think is a very standard situation in
>>> ecology....
>>>
>>> A normally distributed measure (say bird density) is collected for:
>>>
>>>  5 years of data (fixed factor),
>>>  2 seasons within each year (fixed factor),
>>>  3 months within each season (fixed factor)..."
>>>
>>> The response is measured at the same six sites at every visit. What follows
>>> is my interpretation of what she might want, which may not be what she
>>> really wants, but this may not matter for the purpose of creating an
>>> example.
>>>
>>> I will simulate some data from a known model, then use lmer() to fit the
>>> model and compare the estimated parameters with those used in the
>>> simulation.
>>>
>>> I would be interested in comments about
>>>
>>> a)  whether I have simulated the data from my model correctly
>>>
>>> and
>>>
>>> b)  whether my lmer syntax is appropriate for my model.
>>>
>>>
>>> The systematic part of my model is:
>>>
>>>  \mu + \beta_{is} + \gamma_{sm}
>>>
>>> where i varies over years, s over seasons, and m over months.
>>>
>>> (I assume that months are nested within seasons but crossed with years.)
>>>
>>> Now let me decide how these coefficients should vary with site. I am
>>> going to assume that the constant term and the year-season term vary
>>> with site but that month-within-season does not.
>>>
>>> The random part of the model needs to be generated from a 1 + 5*2 =
>>> 11-dimensional multivariate normal. I will assume 1 correlation
>>> parameter, that between two seasons in the same year. All other
>>> correlations are zero.
>>>
>>> Values for fixed effects:
>>>
>>> Suppose  \mu = 10
>>>
>>>   the beta_{i1} are   7, 8, 9, 6, 7.5
>>>   the beta_{i2} are   1, 3, 4, 2, 2.5
>>>
>>>   the gamma_{1m} are  0, 2, 1
>>>   the gamma_{2m} are  0, -2, -2
>>>
>>> The s.d. parameters are  sigma_mu = 3, sigma_beta = 2, sigma = 0.1
>>> The correlation parameter is 0.4.
>>>
>>> Instead of 6 sites I will assume nsite = 60.
>>>
>>> I have written the following to generate data and fit the model:
>>>
>>> library(mvtnorm)
>>> library(lme4)
>>> V = diag(11)
>>> for ( i in 2:10) {
>>>  V[i,i+1] = 0.4
>>>  V[i+1,i] = 0.4
>>>  }
>>> sig = diag(c(3,rep(2,10)))
>>> V = sig%*%V%*%sig
>>> nsite = 60           #Why not?
>>> sigma = 0.1
>>> theta = c(10,7, 1, 8, 3, 9, 4, 6,2, 7.5,2.5)
>>> P = rmvnorm(nsite,mean=theta,sigma=V)
>>> a =  rep(1:10,rep(3,10))
>>> TH =  cbind(rep(1,30),model.matrix(~factor(a)-1))
>>> gg = rep(c(0,2,1,0,-2,-2),5)
>>> y = as.vector(TH%*%t(P)) + rep(gg,nsite) + rnorm(30*nsite,0,sigma)
>>> year = gl(5,2*3,30*nsite)
>>> season = gl(2,3,30*nsite)
>>> month = gl(3,1,30*nsite)
>>> site = gl(nsite,30,30*nsite)
>>> mod = lmer(y ~ year/season + season/month + (1 + year/season | site) )
>>> summary(mod)
>>>
>>>
>>> The output I get is
>>>
>>>       
>>>> library(mvtnorm)
>>>> library(lme4)
>>>>         
>>>       
>>>> V = diag(11)
>>>> for ( i in 2:10) {
>>>>         
>>> +    V[i,i+1] = 0.4
>>> +    V[i+1,i] = 0.4
>>> +    }
>>>       
>>>> sig = diag(c(3,rep(2,10)))
>>>> V = sig%*%V%*%sig
>>>>         
>>>> nsite = 60
>>>> sigma = 0.1
>>>> theta = c(10,7, 1, 8, 3, 9, 4, 6,2, 7.5,2.5)
>>>> P = rmvnorm(nsite,mean=theta,sigma=V)
>>>>
>>>> a =  rep(1:10,rep(3,10))
>>>> TH =  cbind(rep(1,30),model.matrix(~factor(a)-1))
>>>> gg = rep(c(0,2,1,0,-2,-2),5)
>>>> y = as.vector(TH%*%t(P)) + rep(gg,nsite) + rnorm(30*nsite,0,sigma)
>>>> year = gl(5,2*3,30*nsite)
>>>> season = gl(2,3,30*nsite)
>>>> month = gl(3,1,30*nsite)
>>>> site = gl(nsite,30,30*nsite)
>>>> mod = lmer(y ~ year/season + season/month + (1 + year/season | site) )
>>>> summary(mod)
>>>>         
>>> Linear mixed model fit by REML
>>> Formula: y ~ year/season + season/month + (1 + year/season | site)
>>>  AIC  BIC logLik deviance REMLdev
>>> 1332 1717 -596.1     1153    1192
>>> Random effects:
>>> Groups   Name          Variance  Std.Dev.
>>> Corr
>>> site     (Intercept)   14.874181
>>> 3.85671
>>>         year2          8.955943 2.99265
>>> -0.234
>>>         year3          8.751238 2.95825  -0.343
>>> 0.504
>>>         year4          8.671103 2.94467  -0.208  0.450
>>> 0.422
>>>         year5          7.388752 2.71823  -0.324  0.402  0.516
>>> 0.462
>>>         year1:season2  4.283068 2.06956  -0.156  0.638  0.322  0.517
>>> 0.358
>>>         year2:season2  5.046195 2.24637  -0.149 -0.469  0.223 -0.053
>>> 0.015 -0.308
>>>         year3:season2  5.490231 2.34312   0.078 -0.139 -0.418  0.243
>>> -0.007  0.019 -0.294
>>>         year4:season2  5.353655 2.31380   0.135 -0.204 -0.056 -0.583
>>> 0.086 -0.254  0.009 -0.318
>>>         year5:season2  5.400134 2.32382  -0.053  0.047  0.068 -0.044
>>> -0.472 -0.068  0.074 -0.018
>>> Residual                0.010340
>>> 0.10169
>>>
>>> -0.301
>>>
>>> Number of obs: 1800, groups: site, 60
>>>
>>> Fixed effects:
>>>               Estimate Std. Error t value
>>> (Intercept)    17.473505   0.497980   35.09
>>> year2           0.731252   0.386499    1.89
>>> year3           1.965976   0.382058    5.15
>>> year4          -1.014649   0.380309   -2.67
>>> year5           0.688820   0.351081    1.96
>>> season2        -6.255910   0.267479  -23.39
>>> year2:season2   1.889409   0.451063    4.19
>>> year3:season2   1.117068   0.400037    2.79
>>> year4:season2   2.564384   0.448766    5.71
>>> year5:season2   1.055936   0.415366    2.54
>>> season1:month2  1.993234   0.008303  240.07
>>> season2:month2 -2.029299   0.008303 -244.41
>>> season1:month3  1.001748   0.008303  120.65
>>> season2:month3 -2.027852   0.008303 -244.24
>>>
>>> Correlation of Fixed Effects:
>>>           (Intr) year2  year3  year4  year5  seasn2 yr2:s2 yr3:s2
>>> yr4:s2 yr5:s2 ssn1:2 ssn2:2 ssn1:3
>>> year2
>>> -0.235
>>>
>>> year3       -0.343
>>> 0.504
>>>
>>> year4       -0.208  0.450
>>> 0.422
>>> year5       -0.324  0.402  0.516
>>> 0.462
>>> season2     -0.157  0.637  0.322  0.516
>>> 0.358
>>> year2:sesn2 -0.003 -0.680 -0.047 -0.341 -0.203
>>> -0.790
>>> year3:sesn2  0.163 -0.531 -0.531 -0.162 -0.245 -0.654
>>> 0.377
>>> year4:sesn2  0.183 -0.515 -0.229 -0.696 -0.156 -0.765  0.575
>>> 0.343
>>> year5:sesn2  0.062 -0.376 -0.159 -0.364 -0.572 -0.693  0.572  0.444
>>> 0.377
>>> sesn1:mnth2 -0.008  0.000  0.000  0.000  0.000  0.016  0.000  0.000
>>> 0.000  0.000
>>> sesn2:mnth2  0.000  0.000  0.000  0.000  0.000 -0.016  0.000  0.000
>>> 0.000  0.000  0.000
>>> sesn1:mnth3 -0.008  0.000  0.000  0.000  0.000  0.016  0.000  0.000
>>> 0.000  0.000  0.500  0.000
>>> sesn2:mnth3  0.000  0.000  0.000  0.000  0.000 -0.016  0.000  0.000
>>> 0.000  0.000  0.000  0.500  0.000
>>>       
>>>
>>> The variance estimates for the random effects look large, but maybe
>>> because of the way I have parameterized the model I am estimating some
>>> combination of sigma_mu and sigma_beta?
>>>
>>>
>>> Cheers,
>>>
>>>
>>>
>>> Murray Jorgensen
>>>
>>>
>>> --
>>> Dr Murray Jorgensen      http://www.stats.waikato.ac.nz/Staff/maj.html
>>> Department of Statistics, University of Waikato, Hamilton, New Zealand
>>> Email: maj@waikato.ac.nz                                Fax 7 838 4155
>>> Phone  +64 7 838 4773 wk    Home +64 7 825 0441   Mobile 021 0200 8350
>>>
>>> _______________________________________________
>>> R-sig-mixed-models@r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>
>>>       

-- 
Dr Murray Jorgensen      http://www.stats.waikato.ac.nz/Staff/maj.html
Department of Statistics, University of Waikato, Hamilton, New Zealand
Email: maj@waikato.ac.nz                                Fax 7 838 4155
Phone  +64 7 838 4773 wk    Home +64 7 825 0441   Mobile 021 0200 8350


	[[alternative HTML version deleted]]


