Message-ID: <CA+QGQvsNpXa365FyWioNgFeNQqwRit7191OXOXiu=h5obVmuSQ@mail.gmail.com>
Date: 2011-10-05T14:51:19Z
From: Marcus Nunes
Subject: F-values in nested designs
In-Reply-To: <CADv2QyH5FdAwe3YWQezbthqTXRvEQb9dZvrNZdypRzywo26rJA@mail.gmail.com>
Dennis,
thanks for your help. I've read your email and the references you gave
and things are more clear to me.
Best,
Marcus
On Tue, Oct 4, 2011 at 19:28, Dennis Murphy <djmuser at gmail.com> wrote:
>
> 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
> >
> > ______________________________________________
> > R-help at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> >
--
Marcus Nunes
marcus.nunes at gmail.com