Gustaf Granath (PhD)
Post doc
McMaster University
> Dear R-sig-mixed-models list,
> I am new to using lmer and was hoping someone could comment on whether the model I am using is properly coded to address my questions and also to answer a specific question about random nested interactions.
> Main questions: 1) Do crops differ in height compared to wild species and 2) is the effect of domestication consistent across different domestication events.
> Data:
> Six individual plants were measured per plant species. Half of the 62 species are crops and the other half are wild plants.
> Each crop species is paired with a close wild relative (PAIR) and thus I have 31 independent domestication events.
>
> Setup:
> - DOMESTICATION as a fixed effect (crop or wild)
> - PAIR as a random effect, I want to know the general consistency of domestication and not the patterns in each pair specifically
> - I need a DOMESTICATION by PAIR interaction term (which would be random) to see if domestication effects change with PAIR
> - I need to avoid pseudo replication and thus I nested SPECIES within PAIR
>
> Mathematical model:
> Height = global mean+ domestication + pair[random] + domestication*pair[random] + species(pair)[random] + residual error
>
> Working LMER Model??
> m1<- lmer( HEIGHT~ DOMESTICATION + (1 | PAIR : SPECIES) + (DOMESTICATION | PAIR), data=X1, na.action=na.omit)
>
> I included the term (1 | PAIR : SPECIES) to avoid the problem I was having with a previous model (see m2 below).
> Does this model make sense given my design and questions? Have I properly avoided pseudo replication?
>
> Linear mixed model fit by REML
> Formula: HEIGHT ~ DOMESTICATION + (1 | PAIR:SPECIES) + (DOMESTICATION | PAIR)
> Data: X1
> AIC BIC logLik deviance REMLdev
> 2481 2508 -1233 2472 2467
> Random effects:
> Groups Name Variance Std.Dev. Corr
> PAIR:SPECIES (Intercept) 28.769 5.3637
> PAIR (Intercept) 86.789 9.3161
> DOMESTICATIONWild 13.586 3.6859 0.317
> Residual 26.240 5.1225
> Number of obs: 374, groups: PAIR:SPECIES, 62; PAIR, 31
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 20.351 1.967 10.346
> DOMESTICATIONWild 1.388 1.605 0.865
>
> Correlation of Fixed Effects:
> (Intr)
> DOMESTICATI -0.228
>
>
> Problematic LMER model :
> m2<- lmer( HEIGHT~ DOMESTICATION + (DOMESTICATION | PAIR / SPECIES), data=X1, na.action=na.omit)
>
> This is the first model I tried. It has all the effects I want but includes a random effect of SPECIES (nested within PAIR) on the slope (DOMESTICATION) (underlined in output below). I do not understand how this can work given that the data for any given SPECIES is only a single level of DOMESTICATION (all crop or all wild). Am I missing something? Is the m2 approach better than the m1? If so why?
>
> Linear mixed model fit by REML
> Formula: HEIGHT ~ DOMESTICATION + (DOMESTICATION | PAIR/SPECIES)
> Data: X1
> AIC BIC logLik deviance REMLdev
> 2485 2520 -1233 2472 2467
> Random effects:
> Groups Name Variance Std.Dev. Corr
> SPECIES:PAIR (Intercept) 29.923 5.4702
> DOMESTICATIONWild 19.646 4.4323 -0.381
> PAIR (Intercept) 85.633 9.2538
> DOMESTICATIONWild 10.092 3.1768 0.410
> Residual 26.240 5.1225
> Number of obs: 374, groups: SPECIES:PAIR, 62; PAIR, 31
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 20.351 1.967 10.346
> DOMESTICATIONWild 1.388 1.605 0.865
>
> Correlation of Fixed Effects:
> (Intr)
> DOMESTICATI -0.228
>
> thanks for your help
>
> Mart
>
>
> Martin Turcotte, Ph. D.
> Dept. of Biology
> University of Toronto at Mississauga
> 3359 Mississauga Road North, Mississauga,
> Ontario, Canada, L5L 1C6
>
> http://individual.utoronto.ca/martinturcotte