example1: EXAMPLE 1: Split-plot design with one qualitative and one...
In RNED/agriTutorial: Tutorial Analysis of Some Agricultural Experiments
Description
Details
References
Examples
Description
Gomez & Gomez (1984, p. 143) report a rice experiment with three management practices
(minimum, optimum, intensive), five different amounts of nitrogen (N) fertilizer
(0, 50, 80, 110, 140 kg/ha), and three varieties (V1, V2, V3). The experiment
involved variety and management as qualitative treatment factors and nitrogen
fertilizer as a quantitative treatment factor. Overall, there were 45 treatments.
The fertilizer treatments were applied to main plots, the management practices to
split-plots and the varieties to split-split plots.
Details
Section 1
Section 1 examines treatment effects by fitting qualitative factorial models and the first analysis
calculates a full analysis of variance (Table 1) for main plots (nitrogen),
split-plots (management) and split-split plots (variety). Each type of experimental unit (or "stratum")
requires a separate error term in the fitted analysis.
The second analysis (Table 2) uses a REML mixed model analysis to find treatment means and SE's for each marginal
treatment classification averaged over all the other treatment factors, together with estimates of
pairwise contrasts of treatment means and the SE's of the pairwise treatment comparisons. This analysis fits
the full set of nitrogen by variety interaction effects assuming additive managment effects and the fit of the model
is tested by a graphical plot of the model residuals.
Residual plots provide an important check on model assumptions but many more options for model testing are
available and further methods for diagnostic testing will be examined in subsequent examples.
The third analysis (Table 3) shows a mixed model analysis of the full factorial model fitted by REML using the
lmer function of the lme4 package. Generally with mixed models, determination of the denominator degrees of
freedom for Wald-type F- and t-statistics becomes an issue, and here we use the method proposed by
Kenward & Roger (1997).
Section 2
Section 2 examines treatment effects by fitting polynomial models and the first step calculates
a full set of four raw polynomials for the 5-levels of N using the poly() function.
The N rates are re-scaled by division by 100 to improve numerical stability.
The second step fits a mixed model polynomial analysis of nitrogen effects assuming additive management effects (Table 7).
In this analysis, most of the nitrogen treatment effect can be explained by linear and quadratic trend effects.
but it is important to note that there is a non-negligible Variety x Cubic N interaction effect. This suggests
that not all the varieties responded in a similar way to the N treatments and that some further analysis of the
data may be required (see also the N plots of individual varieties and replicates in Fig 1).
The third step fits the required model for the actual fitted model coefficients (Table 8).
When estimating model effects, only effects that are significant for the fitted model or that are marginal to
those effects (functional marginality) should be included in the model therefore only linear and quadratic nitrogen effects
are included in this model. The fitted model for the nitrogen effects fits the actual actual nitrogen levels used
in the experiment therefore this model provides the required coefficients for the actual applied nitrogen levels.
Finally, individual plots of the nitrogen response for each variety in each block averaged over all management effects are
shown (Fig 1). These plots show the nitrogen response of each variety in each replicate block and
show some evidence of anomalous behaviour by variety 1 in blocks 1 and 2 compared with block 3. In practice,
this anomaly would need to be investigated by further discussion with the research workers.
agriTutorial : back to home page
References
Gomez, K.A., & Gomez, A.A. (1984). Statistical procedures for agricultural
research, 2nd edn. New York: Wiley.
Kenward, M.G., & Roger, J.H. (1997). Small sample inference for fixed effects from restricted maximum likelihood. Biometrics, 53, 983–997.
Piepho, H. P. & Edmondson R. N. (accepted). A tutorial on the statistical analysis of factorial experiments with qualitative and quantitative
treatment factor levels.Journal of Agronomy and Crop Science. Accepted
Examples
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
## Copy and paste the following code into a R console or GUI to run examples
## Packages lmerTest, lsmeans and pbkrtest MUST be installed
## Not run:
## *************************************************************************************
## Preliminaries
##**************************************************************************************
## sink("F:\\tutorial2\\OutputsR\\outExample1.txt") #sink file for outputs
## pdf("F:\\tutorial2\\OutputsR\\outExample1_Fig_S1.pdf") #opens a graphical pdf output file
options(contrasts=c('contr.treatment','contr.poly'))
require(lmerTest)
require(lsmeans)
require(pbkrtest)
data(rice)
## write.table(rice, "c:/rice.txt", sep="\t") # export data to a text file
## write.xlsx(rice, "c:/rice.xlsx") # export data to a spread sheet
## *************************************************************************************
## Section 1: Qualitative analysis of factorial treatment effects
##**************************************************************************************
## Table 1 Full analysis of rice data assuming qualitative nitrogen effects
rice.aov1 = aov(yield ~ Replicate + management*variety*nitrogen +
Error(Replicate/Main/Sub),rice)
summary(rice.aov1, ddf="Kenward-Roger",type = 1)
## Table 2 REML means and se's for additive management and qualitative nitrogen effects
rice.means= lmer(yield ~ Replicate + management+ nitrogen*variety +
(1|Replicate:Main)+ (1|Replicate:Main:Sub), data=rice)
anova(rice.means, ddf="Kenward-Roger",type = 1)
plot(rice.means,sub.caption=NA,ylab="Residuals",xlab="Fitted",
main="Full analysis with full nitrogen effects")
## dev.off()# closes graphical device
lsmeans::lsmeans(rice.means,~nitrogen)
lsmeans::lsmeans(rice.means,~variety)
lsmeans::lsmeans(rice.means,~nitrogen*variety)
## REML contrasts and sed's for additive management and qualitative nitrogen effects
n.v= lsmeans::lsmeans(rice.means,~nitrogen|variety)
contrast(n.v, alpha=0.05, method="pairwise")
v.n= lsmeans::lsmeans(rice.means,~variety|nitrogen)
contrast(v.n, alpha=0.05, method="pairwise")
## Table 3 Mixed model effects for rice data with significance tests
rice.lmer= lmer(yield ~ Replicate + nitrogen*management*variety + (1|Replicate:Main) +
(1|Replicate:Main:Sub), data=rice)
anova(rice.lmer, ddf="Kenward-Roger",type = 1)
## *************************************************************************************
## Section 2: Quantitative analysis of factorial treatment effects
##**************************************************************************************
## adds raw N polynomials to data frame: note that the nrate is re-scaled
N=poly((rice$nrate/100),4,raw=TRUE)
colnames(N)=c("Linear_N","Quadratic_N","Cubic_N","Quartic_N")
rice=cbind(rice,N)
## Table 7: Mixed model fitting raw polynomials for nitrogen effects
rice.fullN= lmer(yield ~ Replicate + management + variety*(Linear_N + Quadratic_N +
Cubic_N + Quartic_N) + (1|Replicate:Main)+ (1|Replicate:Main:Sub) , data=rice)
anova(rice.fullN, ddf="Kenward-Roger",type = 1)
## Table 8 Coefficients for separate linear and common quadratic N with additive management
rice.quadN= lmer(yield ~ Replicate + management + variety*Linear_N + Quadratic_N +
(1|Replicate:Main) + (1|Replicate:Main:Sub), data=rice)
summary(rice.quadN, ddf="Kenward-Roger")
## Fig 1 Nitrogen response per variety per plot showing anomalous behaviour of Variety 1
Rice=aggregate(rice$yield, by=list(rice$Replicate, rice$nitrogen,rice$variety),
FUN=mean, na.rm=TRUE)
colnames(Rice)=c("Reps","Nlev","Vars","Yield")
wideRice <- reshape(Rice, timevar="Nlev",idvar=c("Vars","Reps"),direction="wide")
wideRice=wideRice[,-c(1,2)]
N=c(0,50,80,110,140)
par(mfrow=c(3,3),oma=c(0,0,2,0))
for (i in 1:3)
for (j in 1:3)
plot(N,wideRice[(i-1)*3+j,],type = "l",ylab="yield",main=paste("Variety",i,"Block",j),
ylim=c(0, max(wideRice)))
title(main="Fig S1. Variety response to nitrogen in individual replicate blocks",
outer=TRUE)
## *************************************************************************************
## Closure
##**************************************************************************************
## dev.off()# closes graphical device
## sink() #closes sink file
## End(Not run)
RNED/agriTutorial documentation built on May 28, 2019, 2:26 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Details References Examples
Description
Gomez & Gomez (1984, p. 143) report a rice experiment with three management practices (minimum, optimum, intensive), five different amounts of nitrogen (N) fertilizer (0, 50, 80, 110, 140 kg/ha), and three varieties (V1, V2, V3). The experiment involved variety and management as qualitative treatment factors and nitrogen fertilizer as a quantitative treatment factor. Overall, there were 45 treatments. The fertilizer treatments were applied to main plots, the management practices to split-plots and the varieties to split-split plots.
Details
Section 1
Section 1 examines treatment effects by fitting qualitative factorial models and the first analysis calculates a full analysis of variance (Table 1) for main plots (nitrogen), split-plots (management) and split-split plots (variety). Each type of experimental unit (or "stratum") requires a separate error term in the fitted analysis.
The second analysis (Table 2) uses a REML mixed model analysis to find treatment means and SE's for each marginal treatment classification averaged over all the other treatment factors, together with estimates of pairwise contrasts of treatment means and the SE's of the pairwise treatment comparisons. This analysis fits the full set of nitrogen by variety interaction effects assuming additive managment effects and the fit of the model is tested by a graphical plot of the model residuals. Residual plots provide an important check on model assumptions but many more options for model testing are available and further methods for diagnostic testing will be examined in subsequent examples.
The third analysis (Table 3) shows a mixed model analysis of the full factorial model fitted by REML using the lmer function of the lme4 package. Generally with mixed models, determination of the denominator degrees of freedom for Wald-type F- and t-statistics becomes an issue, and here we use the method proposed by Kenward & Roger (1997).
Section 2
Section 2 examines treatment effects by fitting polynomial models and the first step calculates a full set of four raw polynomials for the 5-levels of N using the poly() function. The N rates are re-scaled by division by 100 to improve numerical stability.
The second step fits a mixed model polynomial analysis of nitrogen effects assuming additive management effects (Table 7). In this analysis, most of the nitrogen treatment effect can be explained by linear and quadratic trend effects. but it is important to note that there is a non-negligible Variety x Cubic N interaction effect. This suggests that not all the varieties responded in a similar way to the N treatments and that some further analysis of the data may be required (see also the N plots of individual varieties and replicates in Fig 1).
The third step fits the required model for the actual fitted model coefficients (Table 8). When estimating model effects, only effects that are significant for the fitted model or that are marginal to those effects (functional marginality) should be included in the model therefore only linear and quadratic nitrogen effects are included in this model. The fitted model for the nitrogen effects fits the actual actual nitrogen levels used in the experiment therefore this model provides the required coefficients for the actual applied nitrogen levels.
Finally, individual plots of the nitrogen response for each variety in each block averaged over all management effects are shown (Fig 1). These plots show the nitrogen response of each variety in each replicate block and show some evidence of anomalous behaviour by variety 1 in blocks 1 and 2 compared with block 3. In practice, this anomaly would need to be investigated by further discussion with the research workers.
agriTutorial : back to home page
References
Gomez, K.A., & Gomez, A.A. (1984). Statistical procedures for agricultural research, 2nd edn. New York: Wiley.
Kenward, M.G., & Roger, J.H. (1997). Small sample inference for fixed effects from restricted maximum likelihood. Biometrics, 53, 983–997.
Piepho, H. P. & Edmondson R. N. (accepted). A tutorial on the statistical analysis of factorial experiments with qualitative and quantitative treatment factor levels.Journal of Agronomy and Crop Science. Accepted
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 | ## Copy and paste the following code into a R console or GUI to run examples
## Packages lmerTest, lsmeans and pbkrtest MUST be installed
## Not run:
## *************************************************************************************
## Preliminaries
##**************************************************************************************
## sink("F:\\tutorial2\\OutputsR\\outExample1.txt") #sink file for outputs
## pdf("F:\\tutorial2\\OutputsR\\outExample1_Fig_S1.pdf") #opens a graphical pdf output file
options(contrasts=c('contr.treatment','contr.poly'))
require(lmerTest)
require(lsmeans)
require(pbkrtest)
data(rice)
## write.table(rice, "c:/rice.txt", sep="\t") # export data to a text file
## write.xlsx(rice, "c:/rice.xlsx") # export data to a spread sheet
## *************************************************************************************
## Section 1: Qualitative analysis of factorial treatment effects
##**************************************************************************************
## Table 1 Full analysis of rice data assuming qualitative nitrogen effects
rice.aov1 = aov(yield ~ Replicate + management*variety*nitrogen +
Error(Replicate/Main/Sub),rice)
summary(rice.aov1, ddf="Kenward-Roger",type = 1)
## Table 2 REML means and se's for additive management and qualitative nitrogen effects
rice.means= lmer(yield ~ Replicate + management+ nitrogen*variety +
(1|Replicate:Main)+ (1|Replicate:Main:Sub), data=rice)
anova(rice.means, ddf="Kenward-Roger",type = 1)
plot(rice.means,sub.caption=NA,ylab="Residuals",xlab="Fitted",
main="Full analysis with full nitrogen effects")
## dev.off()# closes graphical device
lsmeans::lsmeans(rice.means,~nitrogen)
lsmeans::lsmeans(rice.means,~variety)
lsmeans::lsmeans(rice.means,~nitrogen*variety)
## REML contrasts and sed's for additive management and qualitative nitrogen effects
n.v= lsmeans::lsmeans(rice.means,~nitrogen|variety)
contrast(n.v, alpha=0.05, method="pairwise")
v.n= lsmeans::lsmeans(rice.means,~variety|nitrogen)
contrast(v.n, alpha=0.05, method="pairwise")
## Table 3 Mixed model effects for rice data with significance tests
rice.lmer= lmer(yield ~ Replicate + nitrogen*management*variety + (1|Replicate:Main) +
(1|Replicate:Main:Sub), data=rice)
anova(rice.lmer, ddf="Kenward-Roger",type = 1)
## *************************************************************************************
## Section 2: Quantitative analysis of factorial treatment effects
##**************************************************************************************
## adds raw N polynomials to data frame: note that the nrate is re-scaled
N=poly((rice$nrate/100),4,raw=TRUE)
colnames(N)=c("Linear_N","Quadratic_N","Cubic_N","Quartic_N")
rice=cbind(rice,N)
## Table 7: Mixed model fitting raw polynomials for nitrogen effects
rice.fullN= lmer(yield ~ Replicate + management + variety*(Linear_N + Quadratic_N +
Cubic_N + Quartic_N) + (1|Replicate:Main)+ (1|Replicate:Main:Sub) , data=rice)
anova(rice.fullN, ddf="Kenward-Roger",type = 1)
## Table 8 Coefficients for separate linear and common quadratic N with additive management
rice.quadN= lmer(yield ~ Replicate + management + variety*Linear_N + Quadratic_N +
(1|Replicate:Main) + (1|Replicate:Main:Sub), data=rice)
summary(rice.quadN, ddf="Kenward-Roger")
## Fig 1 Nitrogen response per variety per plot showing anomalous behaviour of Variety 1
Rice=aggregate(rice$yield, by=list(rice$Replicate, rice$nitrogen,rice$variety),
FUN=mean, na.rm=TRUE)
colnames(Rice)=c("Reps","Nlev","Vars","Yield")
wideRice <- reshape(Rice, timevar="Nlev",idvar=c("Vars","Reps"),direction="wide")
wideRice=wideRice[,-c(1,2)]
N=c(0,50,80,110,140)
par(mfrow=c(3,3),oma=c(0,0,2,0))
for (i in 1:3)
for (j in 1:3)
plot(N,wideRice[(i-1)*3+j,],type = "l",ylab="yield",main=paste("Variety",i,"Block",j),
ylim=c(0, max(wideRice)))
title(main="Fig S1. Variety response to nitrogen in individual replicate blocks",
outer=TRUE)
## *************************************************************************************
## Closure
##**************************************************************************************
## dev.off()# closes graphical device
## sink() #closes sink file
## End(Not run)
|
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.