example3: EXAMPLE 3: Polynomial regression model with two quantitative...
In RNED/agriTutorial: Tutorial Analysis of Some Agricultural Experiments

(Gomez & Gomez, 1984, p. 401) report a two-factor nitrogen uptake greenhouse experiment on rice involving duration of water stress (W) and level of nitrogen application (N) with four complete replicates of each treatment. The experiment had four water-stress levels (0, 10, 20 and 40 days) applied as main-plot treatments and four nitrogen rates (0, 90, 180 and 270 kg/ha) applied as sub-plot treatments. The four sub-plot treatments were randomized within main plots and the four main plot treatments were randomized within complete replicate blocks.

The first stage of the analysis is the calculation of polynomial powers of N and W using the poly() function. The N rates are re-scaled by division by 100 while the W rates are re-scaled by division by 10.

The second stage compares a Pearson residual plot of the untransformed N uptake data versus a Pearson residual plot of the log transformed N uptake data. Comparison of the two plots shows that the untransformed residuals increase as the fitted values increase, which shows that the variance of the untransformed data is non-constant. The log transformed N uptake residuals, however, are approximately constant over the full range of fitted values which shows that the log transformed N uptake data can be used to fit a simple unweighted analysis of variance for the effects of the treatment factors.

Sometimes the original scale of measurement is the proper scale of measurement for an analysis, e.g. an analysis of actual measured crop yields, and then it might be appropriate to fit a weighted analysis of variance in the original scale of measurement of the dependent variable (see Faraway 2002 Chapter 5). However, the log transformation model assumes a proportional rather than an additive model for treatment effects and, in this example, a proportional model for nitrogen uptake may well be a more natural physiological model than a simple additive model.

The next stage compares the fit of a first-order linear model (Table 9) versus a second-order quadratic model (Table 10). The first-order model shows significant lack-of-fit and is not adequate for the data. The second-order model is also not fully adequate for the data as there is a significant N lack of fit term indicating a significant cubic effect. However, the magnitude of the cubic effect is relatively small and it will be assumed here that a quadratic model is adequate for the data.

The final stage fits regression coefficients for the quadratic response surface model on the rescaled water stress and rescaled nitrogen rate treatments. The fitted coefficients are then used to plot the fitted quadratic log uptake curves versus the nitrogen rate treatments and the water stress treatments shown in Fig 4.

Note that in this analysis all the polynomial models are built by adding individual polynomial effects in acccordance with the requirements of functional marginality.

agriTutorial : back to home page

Faraway J (2002) Practical Regression and Anova using R. https://cran.r-project.org/doc/contrib/Faraway-PRA.pdf

Gomez, K.A., & Gomez, A.A. (1984). Statistical procedures for agricultural research, 2nd edn. New York: Wiley.

## Not run: 

## *************************************************************************************
##                            Preliminaries
##**************************************************************************************

## sink("F:\\tutorial2\\OutputsR\\outExample3.txt") #sink file for outputs
## pdf("F:\\tutorial2\\OutputsR\\outExample3_Fig_S2.pdf") #opens a graphical pdf output file
require(lmerTest)
require(lattice)
require(pbkrtest)
options(contrasts=c('contr.treatment','contr.poly'))
data(greenrice)
## write.table(greenrice, "c:/greenrice.txt", sep="\t") # export data to a text file
## write.xlsx(greenrice, "c:/greenrice.xlsx") # export data to a spread sheet

## *************************************************************************************
##            Section 1: Polynomial powers of N and W
##**************************************************************************************

greenrice$loguptake=log(greenrice$uptake)
greenrice$Nitrogen=factor(greenrice$N)
greenrice$Water=factor(greenrice$W)
PolW=poly((greenrice$W/10), degree=2, raw=TRUE)
colnames(PolW)=c("Linear_W","Quadratic_W")
PolN=poly((greenrice$N/100), degree=2, raw=TRUE)
colnames(PolN)=c("Linear_N","Quadratic_N")
greenrice=cbind(greenrice,PolW,PolN)

## residual plots of untransformed versus log transformed uptake data
greenrice.uptake = lmer(uptake ~ Replicate + factor(N) * factor(W)
  + (1|Replicate:Main), data=greenrice)
plot(greenrice.uptake,main="Pearson residual plot for untransformed N uptake",
 ylab="Residuals N uptake")

greenrice.loguptake= lmer(loguptake ~ Replicate +
factor(N) * factor(W) + (1|Replicate:Main), data=greenrice)
plot(greenrice.loguptake,main="Pearson residual plot for log transformed N uptake",
ylab="Residuals log N uptake")

## Tables 9 first-order model of log uptake and Wald tests
greenrice.lmer1 = lmer(loguptake ~ Linear_N + Linear_W +
Nitrogen*Water + (1|Replicate) + (1|Replicate:Main),data=greenrice)
anova(greenrice.lmer1, ddf="Kenward-Roger",type = 1)

## Tables 10 second-order model of log uptake and Wald tests
greenrice.lmer2 = lmer(loguptake  ~Linear_N * Linear_W + Quadratic_N +
Quadratic_W +  Nitrogen*Water + (1|Replicate) + (1|Replicate:Main),
data=greenrice)
anova(greenrice.lmer2, ddf="Kenward-Roger",type = 1)

## *************************************************************************************
##         Section 2 : Fitted regression models and quadratic log uptake curves
##**************************************************************************************

## Regression coefficients of quadratic response model of W and N
greenrice.lmer0 = lmer(loguptake  ~Linear_N * Linear_W + Quadratic_N +
Quadratic_W + (1|Replicate) + (1|Replicate:Main),data=greenrice)
summary(greenrice.lmer0, ddf="Kenward-Roger",type = 1)

panel.plot <- function(x, y) {
panel.xyplot(x, y) # show points
Water=c(0,1.0,2.0,4.0)[panel.number()]
panel.curve(-1.16+0.17603 * Water- 0.11599 * Water * Water +
0.68 * x - 0.0938 * x * x - 0.09072 * Water * x,
from = 0, to = 2.70, type="l", lwd=2)
}
xyplot(loguptake ~ Linear_N|factor(Linear_W),data=greenrice,
 scales=list(x=list(at=c(0,.9,1.8,2.7), labels=c(0,90,180,270))),
 main="Fig 4: Marginal model for nitrogen rate",
 xlab = " Nitrogen (kg/ha)", ylab = "Log nitrogen uptake (g/pot)",
 strip = strip.custom(strip.names = TRUE,
 factor.levels = c("0","10","20","40")),
panel = panel.plot)

panel.plot <- function(x, y) {
panel.xyplot(x, y) # show points
Nitrogen=c(0,.90,1.80,2.70)[panel.number()]
panel.curve(-1.16 + 0.17603 * x - 0.11599 * x * x + 0.68 * Nitrogen -
0.0938 * Nitrogen * Nitrogen - 0.09072 * x * Nitrogen,
from = 0, to = 4.0, type="l", lwd=2)
}
xyplot(loguptake ~ Linear_W|factor(Linear_N),data=greenrice,
 scales=list(x=list(at=c(0,1,2,4), labels=c(0,10,20,40))),
 main="Fig 4: Marginal model for water stress",
 xlab = " Water stress (days)", ylab = "Log nitrogen uptake (g/pot)",
 strip = strip.custom(strip.names = TRUE,
 factor.levels = c("0","90","180","270")),
panel = panel.plot)

## *************************************************************************************
##                                  Closure
##**************************************************************************************

## dev.off()
## sink() #closes sink file

## End(Not run)