F-values in nested designs

Hi:

INB4: if I have a nested design with treatment A and treatment B
within A, F-values are MSA/MSA(B) and MSA(B)/MSE, correct? How can I
make R give these values directly, without further coding?

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)))

a <- lmer(y ~ areas + (1 | areas:sites), data = dn)
a

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()
)

2) why I don't have an F-value for the nested effect? I realize that R
call it as Residuals in the first part of the summary, but there is a
way to make R consider it s another factor?

To get the fixed effects part of Minitab's ANOVA table with lmer(),

anova(a)

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 - pf(0.1579, 2, 9)

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

Hello all

I'm trying to learn how to fit a nested model in R. I found a toy
example on internet where a dataset that have?3 areas and 4 sites
within these areas. When I use Minitab to fit a nested model to this
data, this is the ANOVA table that I got:

Nested ANOVA: y versus areas, sites

Analysis of Variance for y
Source ?DF ? ? ? ?SS ? ? ? MS ? ? ?F ? ? ?P
areas ? ?2 ? ?4.5000 ? 2.2500 ?0.158 ?0.856
sites ? ?9 ?128.2500 ?14.2500 ?3.167 ?0.012
Error ? 24 ?108.0000 ? 4.5000
Total ? 35 ?240.7500

When I use R, this is the ANOVA table that I got:

summary(aov(y ~ areas + Error(areas%in%sites)))

Error: areas:sites
? ? ? ? ?Df Sum Sq Mean Sq F value Pr(>F)
areas ? ? ?2 ? 4.50 ? ?2.25 ?0.1579 0.8563
Residuals ?9 128.25 14.25

Error: Within
? ? ? ? ?Df Sum Sq Mean Sq F value Pr(>F)
Residuals 24 ? ?108 ? ? 4.5
Warning message:
In aov(y ~ areas + Error(areas %in% sites)) : Error() model is singular

The results are the same, except for one F-value and I don't
understand why. Hence, these are my questions:

1) I searched google and I can't find a reason to have this warning in
my code. Why is this happening?

2) why I don't have an F-value for the nested effect? I realize that R
call it as Residuals in the first part of the summary, but there is a
way to make R consider it s another factor?

INB4: if I have a nested design with treatment A and treatment B
within A, F-values are MSA/MSA(B) and MSA(B)/MSE, correct? How can I
make R give these values directly, without further coding?

Thanks for your help.

Below is my code and information about my system.
----------------------
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 = as.factor(rep(c("m1", "m2", "m3"), each=12))
#sites = as.factor(c(rep(c(1, 2, 3, 4), 3), rep(c(5, 6, 7, 8), 3),
rep(c(9, 10, 11, 12), 3)))
sites = as.factor(c(rep(c(1, 2, 3, 4), 9)))
repl ?= as.factor(rep(c(1, 2, 3), each=4, 3))

summary(aov(y ~ areas + Error(areas%in%sites)))

summary(aov(y ~ areas + Error(areas%in%sites)))
Error: areas:sites
? ? ? ? ? Df Sum Sq Mean Sq F value Pr(>F)
areas ? ? ?2 ? 4.50 ? ?2.25 ?0.1579 0.8563
Residuals ?9 128.25 ? 14.25
Error: Within
? ? ? ? ? Df Sum Sq Mean Sq F value Pr(>F)
Residuals 24 ? ?108 ? ? 4.5
Warning message:
In aov(y ~ areas + Error(areas %in% sites)) : Error() model is singular



sessionInfo()
R version 2.13.1 Patched (2011-08-25 r56798)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] splines ? stats ? ? graphics ?grDevices utils ? ? datasets ?methods
[8] base

other attached packages:
[1] car_2.0-11 ? ? ? ? survival_2.36-9 ? ?nnet_7.3-1
[4] MASS_7.3-14 ? ? ? ?lme4_0.999375-40 ? Matrix_0.999375-50
[7] lattice_0.19-33 ? ?nlme_3.1-102

loaded via a namespace (and not attached):
[1] grid_2.13.1 ? stats4_2.13.1 tools_2.13.1
--
Marcus Nunes
marcus.nunes at gmail.com