[R-sig-ME] random effects specification

[Reinhold Kliegl]

Looking at your data, the lmer tells you very nicely the pattern of
means shown in your plot:

> (n.lmerA)
Linear mixed model fit by REML
Formula: n ~ treatment + (1 | id)
   Data: dta
 Subset: community == "A"
    AIC    BIC logLik deviance REMLdev
 -181.7 -164.9  96.83   -218.5  -193.7
Random effects:
 Groups   Name        Variance   Std.Dev.
          (Intercept) 5.5211e-13 7.4304e-07
 Residual             9.8079e-03 9.9035e-02
Number of obs: 120, groups: id, 30

Fixed effects:
            Estimate Std. Error t value
(Intercept) 18.78058    0.01808  1038.7
treatmentb   0.33300    0.02557    13.0
treatmentc   0.01150    0.02557     0.4
treatmentd   0.32939    0.02557    12.9

Treatment b is significanlty higher, c is at the same level, and d is
again significantly higher than treatment a. (With t-values of 13 we
may say that even without HPD intervals.)

The same is true for community B, except that here treatment b and d
are significantly lower than a and that what the estimates tell you.
If you learn how to combine estimates, you can derive the pattern of
means directly from your first analyses.

Reinhold


On Tue, Apr 8, 2008 at 6:22 PM, Sebastian P. Luque <spluque at gmail.com> wrote:
> Hi,
>
>  I think I'm misunderstanding something because I expected a different
>  result below (building up a reproducible example).
>
>  ---<---------------cut here---------------start-------------->---
>  ## Simulating data
>  library(lattice)
>  library(lme4)
>  set.seed(1000)
>  rCom <- rnorm(2, mean=5, sd=0.5)
>  rTr <- rep(rCom / 1.1, 2)
>  nbase <- rnorm(240, 10, 0.1)
>  dta <- within(expand.grid(community=LETTERS[1:2], treatment=letters[1:4],
>                           id=factor(1:30)), {
>                               n <- rCom[as.numeric(community)] +
>                                   rTr[as.numeric(treatment)] + nbase
>                           })
>  dta <- dta[order(dta$community, dta$treatment), ]
>  ## Simulate an interaction
>  dta$n[dta$community == "A"] <- rev(dta$n[dta$community == "A"])
>  ## Have a look
>  xyplot(n ~ treatment | community, data=dta, groups=id,
>        type="b", pch=19, cex=0.3)
>
>  ## We fit LME with community and treatment fixed effects, their
>  ## interactions, and use random effects for subject.
>  n.lmer1 <- lmer(n ~ community * treatment + (1 | id), dta)
>  ## Remove the interaction, for comparison.
>  n.lmer2 <- lmer(n ~ community + treatment + (1 | id), dta)
>  ## Compare these models.
>  anova(n.lmer1, n.lmer2)                 # interaction term needed
>
>  ## So let's test the treatment effect at each community separately.
>  n.lmerA <- lmer(n ~ treatment + (1 | id), dta,
>                 subset=community == "A")
>  ## Here I expected some terms to be significantly different
>  ## from zero, but that's not the case:
>  summary(n.lmerA)
>  ---<---------------cut here---------------end---------------->---
>
>  Looking at the data, there seems to be an obvious treatment effect, yet
>  this doesn't show in the lmer summary.  What am I misunderstanding here?
>  Thanks.
>
>
>
>
>  Cheers,
>
>  --
>  Seb
>
>  _______________________________________________
>  R-sig-mixed-models at r-project.org mailing list
>  https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>