[R-sig-teaching] Analyze Nested and Split-Plot Designs with R
Ramon Diaz-Uriarte
rdiaz02 at gmail.com
Fri Jan 9 13:46:25 CET 2015
Dear Rich and Steven,
I think Steven's comment about differences refer to the F statistics and
p-values for the Suppliers term. In particular,
summary(MontEx14.1.aov1)
does not give the same F statistic for Suppliers as that shown in Table
14.4 of Montgomery (F = 0.97 in Table 14.4, F = 2.85 in the output from
R). But that is because the R output is computing the F statistic as the
ratio of Suppliers Mean Squares over the Error Mean Squares. In contrast,
in Table 14-4 of Montgomery[1] the F-ratio is from the ratio of Suppliers
to Batches Mean Squares. This can be seen from the column labeled "Expected
Mean Squares" (where we see we need this ratio to test the \tau), and this
in turn corresponds to the middle column of Table 14.1, where A is Fixed
and B is Fixed.
So the F-ratio and p-value differences for Suppliers are not an artifact of
rounding, but of using different denominator Mean Squares (the rationale
for which is explained in the previous pages in Montgomery, those where he
goes over the entries in Table 14.1). Thus even if, as Steven says,
"(Pvalue 0.05)", they are larger than 0.05 in very different ways :-)
Of course, the third model that Rich provided:
MontEx14.1.aov3 <- aov(Purity ~ Supplier + Error(Supplier:Batch),
data=MontEx14.1)
summary(MontEx14.1.aov3)
does give identical F statistics and p-values for Supplier (F = 0.97, p =
0.42) as those in Table 14.4 and table 14.6 of Montgomery.
In addition, and to try to make this analysis more clear (though this is
not in Montgomery), I often find that in cases such as this it is
instructive to average the batches: this way, we are forced to assess
variation among Suppliers relative to the among-Batch-within-Supplier
variation (i.e., given the perfect balance, etc, etc, here by decree we use
the right denominator in the ratio to find the F statistic for Supplier)
MontAverage <- aggregate(Purity ~ Supplier + Batch,
data = MontEx14.1, mean)
## F is same as Table 14.4 for Supplier. One can also check that the Sums of
## Squares here are 1/3 of those in Table 14.4 for Supplier, and that the
## Residuals SS here are 1/3 of the SS for Batches in Table 14.4.
summary(aov(Purity ~ Supplier, data = MontAverage))
Best,
R.
[1] I only have access right now to the 5th edition, so I am looking at
Table 13-4 in p. 561, but I think this is exactly the same as Table 14-4 in
the 8th edition.
On Thu, 08-01-2015, at 21:15, Richard M. Heiberger <rmh at temple.edu> wrote:
> Nothing is wrong. The "Warning" makes sure you know that there is no
> numerator term in
> the "Within" stratum. Since that is intended by this design, all is well.
>
> Any differences in numerical values are most likely an artifact of rounding.
>
> Please be specific about what you typed and what you received when
> asking questions on this
> list. In this case, there are three models and three tables in
> Montgomery. Your question doesn't
> say which you are asking about, nor what you differences you are
> seeing. The general comment
> from the R-help list applies here:
>
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
> Rich
>
> On Thu, Jan 8, 2015 at 3:04 PM, Steven Stoline <sstoline at gmail.com> wrote:
>> Dear Rich:
>>
>> First thank you very much. it helped a lot.
>>
>> The F value and its p-values corresponding to the factor "Supplier"are
>> different from the ones in the textbook, Montgomery 8th ed. But conclusion
>> stay same (Pvalue > 0.05)
>>
>>
>>
>> In the last model (the third one) used for Table 14.6 , I am not sure what
>> is wrong.
>>
>>
>>> MontEx14.1.aov3 <- aov(Purity ~ Supplier + Error(Supplier:Batch),
>>> data=MontEx14.1)
>>
>> Warning message:
>>
>> In aov(Purity ~ Supplier + Error(Supplier:Batch), data = MontEx14.1) :
>> Error() model is singular
>>
>>
>>
>>
>> many thanks
>> Steven
>>
>> On Thu, Jan 8, 2015 at 1:58 PM, Richard M. Heiberger <rmh at temple.edu> wrote:
>>>
>>> Steven,
>>>
>>> I assume you mean Montgomery 8th edition (he changed chapter numbers
>>> recently).
>>>
>>> Please state what you expect as output.
>>>
>>> For your first attempt, you have the case wrong (purity instead of
>>> Purity).
>>>
>>> Are you reading Supplier and Batch as factors. It can't be determined
>>> from the
>>> printed table in your email. Use dump (or dput) next time. dput is
>>> designed for R.
>>>
>>> With the above corrections, your first model formula gives Table 14.4
>>> and your second
>>> formula gives Table 14.5.
>>>
>>> To get Table 14.6 you need to use the Error() function in the model
>>> formula.
>>> Here are statements for all three tables.
>>>
>>> Rich
>>>
>>> ## dump("MontEx14.1","")
>>> MontEx14.1 <-
>>> structure(list(Supplier = structure(c(1L, 1L, 1L, 1L, 1L, 1L,
>>> 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
>>> 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("1",
>>> "2", "3"), class = "factor"), Batch = structure(c(1L, 1L, 1L,
>>> 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 1L, 1L, 1L, 2L, 2L, 2L, 3L,
>>> 3L, 3L, 4L, 4L, 4L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L,
>>> 4L), .Label = c("1", "2", "3", "4"), class = "factor"), Purity = c(1L,
>>> -1L, 0L, -2L, -3L, -4L, -2L, 0L, 1L, 1L, 4L, 0L, 1L, -2L, -3L,
>>> 0L, 4L, 2L, -1L, 0L, -2L, 0L, 3L, 2L, 2L, 4L, 0L, -2L, 0L, 2L,
>>> 1L, -1L, 2L, 3L, 2L, 1L)), .Names = c("Supplier", "Batch", "Purity"
>>> ), row.names = c(NA, -36L), class = "data.frame")
>>>
>>>
>>>
>>> MontEx14.1$Supplier <- factor(MontEx14.1$Supplier)
>>> MontEx14.1$Batch <- factor(MontEx14.1$Batch)
>>>
>>> MontEx14.1.aov1 <- aov(Purity ~ Supplier/Batch, data=MontEx14.1)
>>> summary(MontEx14.1.aov1)
>>>
>>> MontEx14.1.aov2 <- aov(Purity ~ Supplier*Batch, data=MontEx14.1)
>>> summary(MontEx14.1.aov2)
>>>
>>> MontEx14.1.aov3 <- aov(Purity ~ Supplier + Error(Supplier:Batch),
>>> data=MontEx14.1)
>>> summary(MontEx14.1.aov3)
>>>
>>>
>>> On Thu, Jan 8, 2015 at 1:06 PM, Steven Stoline <sstoline at gmail.com> wrote:
>>> > Dear All:
>>> >
>>> > example 14.1, Montgomery, chapter 14. *Supplier* is a *fixed factor*,
>>> > *Batches* is a *random factor* nested within the fixed factor Supplier.
>>> >
>>> > How to analyze these data in R in two ways:
>>> >
>>> > 1- Nested Design
>>> >
>>> > fit <- aov(purity~Supplier/Batch)
>>> >
>>> >
>>> > it did not give me the expected output.
>>> >
>>> >
>>> > 2- as a factorial (suppliers Fixed, Batches Random)
>>> >
>>> > fit.out <- aov(Purity~Supplier*Batch, data=have)
>>> >
>>> > it did not give me the expected output.
>>> >
>>> >
>>> > Here is the data set:
>>> > ===============
>>> >
>>> >> data
>>> > Supplier Batch Purity
>>> > 1 1 1
>>> > 1 1 -1
>>> > 1 1 0
>>> > 1 2 -2
>>> > 1 2 -3
>>> > 1 2 -4
>>> > 1 3 -2
>>> > 1 3 0
>>> > 1 3 1
>>> > 1 4 1
>>> > 1 4 4
>>> > 1 4 0
>>> > 2 1 1
>>> > 2 1 -2
>>> > 2 1 -3
>>> > 2 2 0
>>> > 2 2 4
>>> > 2 2 2
>>> > 2 3 -1
>>> > 2 3 0
>>> > 2 3 -2
>>> > 2 4 0
>>> > 2 4 3
>>> > 2 4 2
>>> > 3 1 2
>>> > 3 1 4
>>> > 3 1 0
>>> > 3 2 -2
>>> > 3 2 0
>>> > 3 2 2
>>> > 3 3 1
>>> > 3 3 -1
>>> > 3 3 2
>>> > 3 4 3
>>> > 3 4 2
>>> > 3 4 1
>>> >
>>> >
>>> > many thanks
>>> > Steven
>>> >
>>> > --
>>> > Steven M. Stoline
>>> > 1123 Forest Avenue
>>> > Portland, ME 04112
>>> > sstoline at gmail.com
>>> >
>>> > [[alternative HTML version deleted]]
>>> >
>>> > _______________________________________________
>>> > R-sig-teaching at r-project.org mailing list
>>> > https://stat.ethz.ch/mailman/listinfo/r-sig-teaching
>>
>>
>>
>>
>> --
>> Steven M. Stoline
>> 1123 Forest Avenue
>> Portland, ME 04112
>> sstoline at gmail.com
>
> _______________________________________________
> R-sig-teaching at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-teaching
--
Ramon Diaz-Uriarte
Department of Biochemistry, Lab B-25
Facultad de Medicina
Universidad Autónoma de Madrid
Arzobispo Morcillo, 4
28029 Madrid
Spain
Phone: +34-91-497-2412
Email: rdiaz02 at gmail.com
ramon.diaz at iib.uam.es
http://ligarto.org/rdiaz
More information about the R-sig-teaching
mailing list