This script is to instruct people to analyze data from split-plot design and hierachical design
First create a data set about wheat and herbicide, learning the effects of wheat variety(A0,A1),and rates of a herbicide(B0,B1,B2,B3)
design as follows: 3 blocks(replicates), in each block there are two plots (whole-plot) to grow a varietey of wheat; Besides, the whole-plot is divided into four sub-plots to apply different rates of a herbicide.
rm(list=ls()) # clean up memory
Dat <- data.frame(list(
observation = 1:24,
replicate = rep(c("I","II","III"),each =8),
A = rep(rep(c("A0","A1"),each = 4),3),
B = rep(c("B0","B1","B2","B3"),times = 6 ),
yield = c(13.8,15.5,21.0,18.9,19.3,22.2,25.3,25.9,
13.5,15.0,22.7,18.3,18.0,24.2,24.8,26.7,
13.2,15.2,22.3,19.6,20.5,25.4,28.4,27.6)
))Without knowledge of the design, the analysis below is not accurate.
aov.bad <- aov(yield~replicate+ A*B, data =Dat)
summary(aov.bad)## Df Sum Sq Mean Sq F value Pr(>F)
## replicate 2 7.87 3.93 4.486 0.0312 *
## A 1 262.02 262.02 298.862 7.68e-11 ***
## B 3 215.26 71.75 81.843 4.08e-09 ***
## A:B 3 18.70 6.23 7.109 0.0039 **
## Residuals 14 12.27 0.88
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Library
library(agricolae)
## if we think replicates, no blocking effect
aov.good<- aov(yield~ A*B+Error(replicate:A), data =Dat)## Warning in aov(yield ~ A * B + Error(replicate:A), data = Dat): Error()
## model is singularsummary(aov.good)##
## Error: replicate:A
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 262.0 262.02 81.24 0.000839 ***
## Residuals 4 12.9 3.23
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## B 3 215.26 71.75 118.96 3.43e-09 ***
## A:B 3 18.70 6.23 10.33 0.00121 **
## Residuals 12 7.24 0.60
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## if we think we need consider block effect. put replicate in the model
aov.good2<- aov(yield~ replicate + A*B + Error(replicate:A), data =Dat)## Warning in aov(yield ~ replicate + A * B + Error(replicate:A), data = Dat):
## Error() model is singularsummary(aov.good2)##
## Error: replicate:A
## Df Sum Sq Mean Sq F value Pr(>F)
## replicate 2 7.87 3.93 1.562 0.39032
## A 1 262.02 262.02 104.062 0.00947 **
## Residuals 2 5.04 2.52
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## B 3 215.26 71.75 118.96 3.43e-09 ***
## A:B 3 18.70 6.23 10.33 0.00121 **
## Residuals 12 7.24 0.60
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Do some diagnosis to check the assumptions, mainly on the residuals
#diagnosis
aov.good3<- aov(yield~ replicate + A+replicate:A + B + A:B, data =Dat)
summary(aov.good3)## Df Sum Sq Mean Sq F value Pr(>F)
## replicate 2 7.87 3.93 6.520 0.01211 *
## A 1 262.02 262.02 434.388 8.61e-11 ***
## B 3 215.26 71.75 118.956 3.43e-09 ***
## replicate:A 2 5.04 2.52 4.174 0.04206 *
## A:B 3 18.70 6.23 10.333 0.00121 **
## Residuals 12 7.24 0.60
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1layout(matrix(1:4,2,byrow = T))
plot(aov.good3, 1)
plot(aov.good3, 2)
plot(aov.good3, 5)
library(car)
boxCox(aov.good3)
layout(1)Part 2: R code to analyze hierarchical design. The table is from table 9.8 on page 350.
The hierarchical model is : \(Y_{ijk} = \mu +\alpha_i + \beta_{(i)j}+\epsilon_{ijk}\)
where: \(\mu\): constant; overall mean
\(\alpha_i\): constant for ith treatment group; deviation from mean of i
\(\beta_ij\): a random effect due to the ith group nested witin the jth experimental unit
\(\epsilon_{ijk}\): andom deviation associated with each observation
pigment <- data.frame(list(
batch = rep(1:15,each=4),
sample = rep(1:30,each =2),
test = 1:60,
y = c( 40,39,30,30,26,28,25,26,29,28,14,15,30,31,24,24,
19,20,17,17,33,32,26,24,23,24,32,33,34,34,29,29,
27,27,31,31,13,16,27,24,25,23,25,27,29,29,31,32,
19,20,29,30,23,23,25,25,39,37,26,28)
))
pigment$batch <- as.factor(pigment$batch)
pigment$sample <- as.factor(pigment$sample)
pigment$test <- as.factor(pigment$test)
aov.pigment <- aov(y~batch+sample+test,data=pigment)
summary(aov.pigment)## Df Sum Sq Mean Sq
## batch 14 1216.2 86.87
## sample 15 871.5 58.10
## test 30 27.0 0.90#another method:
anova(lm(y~batch/sample,data=pigment)) #test part is just random noise## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## batch 14 1216.2 86.874 96.526 < 2.2e-16 ***
## batch:sample 15 871.5 58.100 64.556 < 2.2e-16 ***
## Residuals 30 27.0 0.900
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1