Skip to content
Prev 273367 / 398506 Next

F-values in nested designs

Dennis Murphy 
Hi:
This is how to get an equivalent model in lme4, but it probably isn't
what you expect (particularly the 'without further coding' part).
Using your example,

library('lme4')
dn <- data.frame( y = c(10, 12, 8, 13, 14, 8, 10, 12,
                        9, 10, 12, 11, 11, 13, 9, 10, 14,
                        11, 10, 9, 8, 9, 8, 8, 13, 14, 7,
                        10, 10, 13, 9, 7, 16, 12, 5, 4),
                 areas = factor(rep(c("m1", "m2", "m3"), each=12)),
                 sites = factor(rep(1:4, 9)))
Linear mixed model fit by REML
Formula: y ~ areas + (1 | areas:sites)
   Data: dn
   AIC BIC logLik deviance REMLdev
 171.1 179 -80.56      167   161.1
Random effects:
 Groups      Name        Variance Std.Dev.
 areas:sites (Intercept) 3.25     1.8028
 Residual                4.50     2.1213
Number of obs: 36, groups: areas:sites, 12

Fixed effects:
            Estimate Std. Error t value
(Intercept)   10.750      1.090   9.865
areasm2       -0.750      1.541  -0.487
areasm3       -0.750      1.541  -0.487

Correlation of Fixed Effects:
        (Intr) aresm2
areasm2 -0.707
areasm3 -0.707  0.500
##------

lme4 reports the estimated variance components and their square roots,
the standard deviation components (Std.Dev). The estimated residual
variance component is 4.5, which is the same as the residual MSE from
the Minitab output. The estimated variance component associated with
sites nested within areas (areas:sites) is 3.25. Since the design is
balanced, the expected mean square of this term (assuming the model
assumptions are correct) is $\sigma_e^2 + 3 \sigma_s^2$, which is
estimated by 4.5 + 3(3.25) = 14.25, the observed mean square for sites
within areas, again coinciding with the Minitab output. However,
lmer() does not report the result of an F-test for the 'significance'
of the sites variance component, because the null hypothesis
$\sigma_s^2 = 0$ is on the boundary of the parameter space and there
are questions about the reliability of p-values for such tests. See
http://rwiki.sciviews.org/doku.php?id=guides:lmer-tests   In other
words, don't accept the reported p-value re the sites variance from
Minitab on faith. This answers (even in multi-stratum models in aov()
)
To get the fixed effects part of Minitab's ANOVA table with lmer(),
Analysis of Variance Table
      Df Sum Sq Mean Sq F value
areas  2 1.4208 0.71039  0.1579

Once again, the p-value is not reported (by design). Assuming that the
specified normal-theory based model is correct, the conventional F
test for testing the null hypothesis of equal area means would be the
mean square ratio of areas to sites, which would have an F(2, 9)
distribution under the null hypothesis. The p-value of that test would
be
[1] 0.85625

Apart from the needless test of the sites within areas variance
component, the lmer() output corresponds to that of the Minitab table.
The output from lmer() gives you the capacity to do much more, but it
helps to understand some of the theory behind mixed models first.

The transition from fixed effects ANOVA to random effects and mixed
models is not a smooth one - multiple sources of random variation
complicate both testing and confidence/prediction interval procedures,
with several messages on R-sig-mixed-models (including the one cited
above) discussing such issues at length.

As I said, this is probably not what you expected.

Dennis
On Tue, Oct 4, 2011 at 11:17 AM, Marcus Nunes <marcus.nunes at gmail.com> wrote:

Thread (3 messages)

Marcus Nunes F-values in nested designs Oct 4 Dennis Murphy F-values in nested designs Oct 4 Marcus Nunes F-values in nested designs Oct 5