example4: EXAMPLE 4: One qualitative treatment factor with repeated...
In RNED/agriTutorial: Tutorial Analysis of Some Agricultural Experiments
Description
Details
References
Examples
Description
Milliken & Johnson (1992, p. 429) describe an experiment with four sorghum varieties in which the leaf area
index was assessed in five consecutive weeks starting two weeks after emergence. The experiment was laid out
in five randomized complete blocks and had one qualitative factor (variety) and one quantitative factor (week).
The week factor was a repeated measurement taken on each plot on five consecutive occasions over time,
which means that successive measurements on the same plot are likely to be serially correlated. For a valid
analysis, the serial correlations between the repeated measures must be modelled by assuming a
suitable correlation structure for the repeated observations on the individual experimental units.
Details
The first stage of the analysis is the calculation of raw polynomials for weeks and orthogonal
polynomials for blocks using the poly() function. We use orthogonal block polynomials as these give
contrasts which are orthogonal to the overall block mean.
The second stage fits and compares five different correlation structures for the repeated measures
using the gls() function of the nlme package. The goodness of fit of the models is compared by AIC statistics
where the smaller the AIC the better the fit. Here, the AR(1)+nugget model fitted by the corExp() function gave
the best fitting model.
The third stage fits a full regression model over weeks (Table A1) to test for possible interactions between
variety and week effects. The full regression model is then decomposed into individual polynomial contrasts over
weeks (Table A2) to find the most parsimonious model adequate for the data. The analysis into single degree
of freedom polynomial contrasts shows that the variety-by-weeks interaction is mainly due to the linear and
quadratic interaction effects although there is some evidence of higher-degree polynomial intraction effects.
The fourth stage estimates the fitted model coefficients assuming a quadratic regression model
for variety-by-week interaction effects (Table 15).This model uses orthogonal polynomial contrasts to fit block
and block-by-weeks interaction contrasts before fitting the quadratic variety-by-week interaction effects.
Finally, studentized residuals from the quadratic regression model are plotted to test the model assumptions.
The residual plot shows some evidence that the smallest fitted values have the largest
positive residuals and this suggests that some further investigation of the adequacy of the fitted model
would be valuable.
A problem with the gls() function is that it must contain the same polynomial terms for weeks in the
blocks regression model as in the treatments regression model,
which is why we have used raw polynomials for the blocks regression model.
However for a long series of repeated measures, raw polynomials can become numerically unstable
and will eventually fail. The final generalization shows how higher-degree orthogonal
polynomials CAN be used for the blocks regression model PROVIDED that raw polynomials are used
for the blocks regression model terms that match the raw polynomials used for the treatments regression model.
Then the block and the treatments model will both contain the same set of raw polynomials allowing
the gls() algorithm to fit the same set of raw polynomials for both the blocks and the treatments model.
This formulation allows numerically stable orthogonal polynomials to be used for higher-degree blocks model
effects while still allowing raw polynomials to be used for the treatments model.
agriTutorial : back to home page
References
Milliken, G.A., & Johnson, D.E. (1992). Analysis of messy data. Volume I: Designed experiments. Boca Raton: CRC Press.
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
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
## Not run:
## *************************************************************************************
## Preliminaries
##**************************************************************************************
## sink("F:\\tutorial2\\OutputsR\\outExample4.txt") #sink file for outputs
## pdf("F:\\tutorial2\\OutputsR\\outExample4_Fig_S1.pdf") #opens a graphical pdf output file
options(contrasts=c('contr.treatment','contr.poly'))
require(nlme)
data(sorghum)
## write.table(sorghum, "c:/sorghum.txt", sep="\t") # export data to a text file
## write.xlsx(sorghum, "c:/sorghum.xlsx") # export data to a spread sheet
## *************************************************************************************
## Section 1: Polynomials for weeks and blocks contrasts
##**************************************************************************************
sorghum$factblock=factor(sorghum$varblock)
PolWeek=poly(sorghum$varweek, degree=4, raw=TRUE)
colnames(PolWeek)=c("linWeek","quadWeek","cubWeek","quartWeek")
sorghum=cbind(sorghum,PolWeek)
PolBlocks=poly(sorghum$varblock, degree=4, raw=FALSE)
colnames(PolBlocks)=c("linBlock","quadBlock","cubBlock","quartBlock")
sorghum=cbind(sorghum,PolBlocks)
## *************************************************************************************
## Section 2: Compares correlation stuctures for repeated measures
##**************************************************************************************
AIC=NULL
logLik=NULL
Model=c("ID","CS","AR(1)","AR(1)+nugget","UN")
## independent uncorrelated random plots
full_indy = gls(y ~ factweek * (Replicate+variety),sorghum)
anova(full_indy)
AIC=c(AIC, AIC(full_indy))
logLik=c(logLik,logLik(full_indy))
## corCompSymm compound symmetry
corCompSymm = gls(y~factweek*(Replicate+variety),corr=corCompSymm(form=~varweek|factplot),sorghum)
anova(corCompSymm)
AIC=c(AIC, AIC(corCompSymm))
logLik=c(logLik,logLik(corCompSymm))
## corExp without nugget
corExp=gls(y~factweek*(Replicate+variety),corr=corExp(form =~ varweek|factplot),sorghum)
anova(corExp)
AIC=c(AIC, AIC(corExp))
logLik=c(logLik,logLik(corExp))
## corExp with nugget
corExp_nugget=gls(y~factweek*(Replicate+variety),corr=corExp(form=~varweek|factplot,nugget=TRUE),
sorghum)
anova(corExp_nugget)
AIC=c(AIC, AIC(corExp_nugget))
logLik=c(logLik,logLik(corExp_nugget))
## corSymm unstructured
corSymm = gls(y~factweek*(Replicate+variety),corr=corSymm(form = ~1|factplot),
weights=varIdent(form = ~ 1|varweek), sorghum)
anova(corSymm)
AIC=c(AIC, AIC(corSymm))
logLik=c(logLik,logLik(corSymm))
## Table 11 Comparison of log Likelihood and AIC statistics for different correlation structures
dAIC=AIC-AIC[4]
logLik=-logLik*2
dlogLik=logLik-logLik[4]
AICtable=data.frame(Model,round(logLik,2),round(dlogLik,2),round(AIC,2),round(dAIC,2) )
colnames(AICtable)=c("Covar_Model","-2logLr","-diff2logLr","AIC","diffAIC")
AICtable
## *************************************************************************************
## Section 3: Tests for interactions between variety and and week effects
##**************************************************************************************
## Table A1 Sequential Wald tests for full model sorghum data
full_Wald = gls(y ~ (factblock + variety) * factweek ,
corr = corExp(form = ~ varweek | factplot, nugget=TRUE), sorghum)
anova(full_Wald)
## Table A2 (cf Table 14) Sequential Wald tests for full model sorghum data
pol_Wald = gls(y ~ (factblock + variety) *(linWeek + quadWeek + cubWeek + quartWeek),
corr = corExp(form = ~ varweek | factplot, nugget=TRUE), sorghum)
anova(pol_Wald)
## *************************************************************************************
## Section 4: Fitted quadratic model for variety-by-week interaction effects
##**************************************************************************************
## Table 15 quadratic model coefficients
quad_Wald = gls(y ~ PolBlocks + PolBlocks:(linWeek + quadWeek + cubWeek + quartWeek) +
variety * (linWeek + quadWeek),corr = corExp(form = ~ varweek | factplot, nugget=TRUE), sorghum)
anova(quad_Wald)
summary(quad_Wald)$tTable
plot(quad_Wald,sub.caption=NA,main="Residuals from quadratic model")
## *************************************************************************************
## Section 5: Adding orthogonal polynomials for high-degree block effects
##**************************************************************************************
## Generalization: fitting orthogonal polynomials for a long series of repeated measures:
orthoPolWeek=poly(sorghum$varweek, degree=4, raw=FALSE)
quad_orthog_Wald = gls(y ~ PolBlocks + PolBlocks:(linWeek + quadWeek + orthoPolWeek[,3:4]) +
variety * (linWeek + quadWeek),corr = corExp(form = ~ varweek | factplot, nugget=TRUE), sorghum)
anova(quad_orthog_Wald)
summary(quad_orthog_Wald)$tTable
## *************************************************************************************
## 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
Milliken & Johnson (1992, p. 429) describe an experiment with four sorghum varieties in which the leaf area index was assessed in five consecutive weeks starting two weeks after emergence. The experiment was laid out in five randomized complete blocks and had one qualitative factor (variety) and one quantitative factor (week). The week factor was a repeated measurement taken on each plot on five consecutive occasions over time, which means that successive measurements on the same plot are likely to be serially correlated. For a valid analysis, the serial correlations between the repeated measures must be modelled by assuming a suitable correlation structure for the repeated observations on the individual experimental units.
Details
The first stage of the analysis is the calculation of raw polynomials for weeks and orthogonal polynomials for blocks using the poly() function. We use orthogonal block polynomials as these give contrasts which are orthogonal to the overall block mean.
The second stage fits and compares five different correlation structures for the repeated measures using the gls() function of the nlme package. The goodness of fit of the models is compared by AIC statistics where the smaller the AIC the better the fit. Here, the AR(1)+nugget model fitted by the corExp() function gave the best fitting model.
The third stage fits a full regression model over weeks (Table A1) to test for possible interactions between variety and week effects. The full regression model is then decomposed into individual polynomial contrasts over weeks (Table A2) to find the most parsimonious model adequate for the data. The analysis into single degree of freedom polynomial contrasts shows that the variety-by-weeks interaction is mainly due to the linear and quadratic interaction effects although there is some evidence of higher-degree polynomial intraction effects.
The fourth stage estimates the fitted model coefficients assuming a quadratic regression model for variety-by-week interaction effects (Table 15).This model uses orthogonal polynomial contrasts to fit block and block-by-weeks interaction contrasts before fitting the quadratic variety-by-week interaction effects.
Finally, studentized residuals from the quadratic regression model are plotted to test the model assumptions. The residual plot shows some evidence that the smallest fitted values have the largest positive residuals and this suggests that some further investigation of the adequacy of the fitted model would be valuable.
A problem with the gls() function is that it must contain the same polynomial terms for weeks in the blocks regression model as in the treatments regression model, which is why we have used raw polynomials for the blocks regression model. However for a long series of repeated measures, raw polynomials can become numerically unstable and will eventually fail. The final generalization shows how higher-degree orthogonal polynomials CAN be used for the blocks regression model PROVIDED that raw polynomials are used for the blocks regression model terms that match the raw polynomials used for the treatments regression model. Then the block and the treatments model will both contain the same set of raw polynomials allowing the gls() algorithm to fit the same set of raw polynomials for both the blocks and the treatments model. This formulation allows numerically stable orthogonal polynomials to be used for higher-degree blocks model effects while still allowing raw polynomials to be used for the treatments model.
agriTutorial : back to home page
References
Milliken, G.A., & Johnson, D.E. (1992). Analysis of messy data. Volume I: Designed experiments. Boca Raton: CRC Press.
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 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 | ## Not run:
## *************************************************************************************
## Preliminaries
##**************************************************************************************
## sink("F:\\tutorial2\\OutputsR\\outExample4.txt") #sink file for outputs
## pdf("F:\\tutorial2\\OutputsR\\outExample4_Fig_S1.pdf") #opens a graphical pdf output file
options(contrasts=c('contr.treatment','contr.poly'))
require(nlme)
data(sorghum)
## write.table(sorghum, "c:/sorghum.txt", sep="\t") # export data to a text file
## write.xlsx(sorghum, "c:/sorghum.xlsx") # export data to a spread sheet
## *************************************************************************************
## Section 1: Polynomials for weeks and blocks contrasts
##**************************************************************************************
sorghum$factblock=factor(sorghum$varblock)
PolWeek=poly(sorghum$varweek, degree=4, raw=TRUE)
colnames(PolWeek)=c("linWeek","quadWeek","cubWeek","quartWeek")
sorghum=cbind(sorghum,PolWeek)
PolBlocks=poly(sorghum$varblock, degree=4, raw=FALSE)
colnames(PolBlocks)=c("linBlock","quadBlock","cubBlock","quartBlock")
sorghum=cbind(sorghum,PolBlocks)
## *************************************************************************************
## Section 2: Compares correlation stuctures for repeated measures
##**************************************************************************************
AIC=NULL
logLik=NULL
Model=c("ID","CS","AR(1)","AR(1)+nugget","UN")
## independent uncorrelated random plots
full_indy = gls(y ~ factweek * (Replicate+variety),sorghum)
anova(full_indy)
AIC=c(AIC, AIC(full_indy))
logLik=c(logLik,logLik(full_indy))
## corCompSymm compound symmetry
corCompSymm = gls(y~factweek*(Replicate+variety),corr=corCompSymm(form=~varweek|factplot),sorghum)
anova(corCompSymm)
AIC=c(AIC, AIC(corCompSymm))
logLik=c(logLik,logLik(corCompSymm))
## corExp without nugget
corExp=gls(y~factweek*(Replicate+variety),corr=corExp(form =~ varweek|factplot),sorghum)
anova(corExp)
AIC=c(AIC, AIC(corExp))
logLik=c(logLik,logLik(corExp))
## corExp with nugget
corExp_nugget=gls(y~factweek*(Replicate+variety),corr=corExp(form=~varweek|factplot,nugget=TRUE),
sorghum)
anova(corExp_nugget)
AIC=c(AIC, AIC(corExp_nugget))
logLik=c(logLik,logLik(corExp_nugget))
## corSymm unstructured
corSymm = gls(y~factweek*(Replicate+variety),corr=corSymm(form = ~1|factplot),
weights=varIdent(form = ~ 1|varweek), sorghum)
anova(corSymm)
AIC=c(AIC, AIC(corSymm))
logLik=c(logLik,logLik(corSymm))
## Table 11 Comparison of log Likelihood and AIC statistics for different correlation structures
dAIC=AIC-AIC[4]
logLik=-logLik*2
dlogLik=logLik-logLik[4]
AICtable=data.frame(Model,round(logLik,2),round(dlogLik,2),round(AIC,2),round(dAIC,2) )
colnames(AICtable)=c("Covar_Model","-2logLr","-diff2logLr","AIC","diffAIC")
AICtable
## *************************************************************************************
## Section 3: Tests for interactions between variety and and week effects
##**************************************************************************************
## Table A1 Sequential Wald tests for full model sorghum data
full_Wald = gls(y ~ (factblock + variety) * factweek ,
corr = corExp(form = ~ varweek | factplot, nugget=TRUE), sorghum)
anova(full_Wald)
## Table A2 (cf Table 14) Sequential Wald tests for full model sorghum data
pol_Wald = gls(y ~ (factblock + variety) *(linWeek + quadWeek + cubWeek + quartWeek),
corr = corExp(form = ~ varweek | factplot, nugget=TRUE), sorghum)
anova(pol_Wald)
## *************************************************************************************
## Section 4: Fitted quadratic model for variety-by-week interaction effects
##**************************************************************************************
## Table 15 quadratic model coefficients
quad_Wald = gls(y ~ PolBlocks + PolBlocks:(linWeek + quadWeek + cubWeek + quartWeek) +
variety * (linWeek + quadWeek),corr = corExp(form = ~ varweek | factplot, nugget=TRUE), sorghum)
anova(quad_Wald)
summary(quad_Wald)$tTable
plot(quad_Wald,sub.caption=NA,main="Residuals from quadratic model")
## *************************************************************************************
## Section 5: Adding orthogonal polynomials for high-degree block effects
##**************************************************************************************
## Generalization: fitting orthogonal polynomials for a long series of repeated measures:
orthoPolWeek=poly(sorghum$varweek, degree=4, raw=FALSE)
quad_orthog_Wald = gls(y ~ PolBlocks + PolBlocks:(linWeek + quadWeek + orthoPolWeek[,3:4]) +
variety * (linWeek + quadWeek),corr = corExp(form = ~ varweek | factplot, nugget=TRUE), sorghum)
anova(quad_orthog_Wald)
summary(quad_orthog_Wald)$tTable
## *************************************************************************************
## 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.