Chapter 44 Split-Split Plots
44.0.1 Weed Biomass Data
This is Example 16.7 of the text. Tables are blocks, nitrogen is the whole plot treatment factor; wetlands are split plots, and weed seeding is the split plot treatment factor; sections of wetlands (trays) are split split plots, and clipping is the split split plot treatment factor.
> library(cfcdae)
> library(lme4)
> library(car)
> data(WeedBiomass)
44.0.2 Modeling
Here is the basic split-split plot analysis. We can use the nesting notation in the random part because wetlands are nested in trays. We can do blocks as fixed or random.
> WB.lmer <- lmer(biomass ~ (1|table) + nitrogen*weed*clip + (1|tray/wetland),data=WeedBiomass)
boundary (singular) fit: see help('isSingular')
Constant variance looks good, and normality is not bad.
> plot(WB.lmer)

Figure 44.1: Residual plot for weed biomass data
> qqnorm(residuals(WB.lmer))

Figure 44.2: NPP of residuals for weed biomass data
We can see a big tray random effect, and a wetland effect that is larger than the residual noise. This also shows up in the different sizes of the SEs for estimates at different levels. On the other hand, the table (block) effect is estimated to be 0.
Notice that non-weed biomass decreases as nitrogen increases (that is, the weeds suck up the nitrogen and … grow like weeds).
> summary(WB.lmer)
Linear mixed model fit by REML ['lmerMod']
Formula: biomass ~ (1 | table) + nitrogen * weed * clip + (1 | tray/wetland)
Data: WeedBiomass
REML criterion at convergence: 168.3
Scaled residuals:
Min 1Q Median 3Q Max
-1.1489 -0.3925 0.0000 0.3925 1.1489
Random effects:
Groups Name Variance Std.Dev.
wetland:tray (Intercept) 2.617 1.618
tray (Intercept) 11.169 3.342
table (Intercept) 0.000 0.000
Residual 1.068 1.034
Number of obs: 48, groups: wetland:tray, 24; tray, 8; table, 2
Fixed effects:
Estimate Std. Error t value
(Intercept) 64.23333 1.23587 51.974
nitrogen1 10.98333 2.14058 5.131
nitrogen2 3.35833 2.14058 1.569
nitrogen3 -3.19167 2.14058 -1.491
weed1 17.04792 0.51244 33.268
weed2 -9.42708 0.51244 -18.396
clip1 -1.61667 0.14919 -10.836
nitrogen1:weed1 -8.58958 0.88757 -9.678
nitrogen2:weed1 -2.31458 0.88757 -2.608
nitrogen3:weed1 2.68542 0.88757 3.026
nitrogen1:weed2 4.08542 0.88757 4.603
nitrogen2:weed2 1.76042 0.88757 1.983
nitrogen3:weed2 -1.16458 0.88757 -1.312
nitrogen1:clip1 0.08333 0.25840 0.322
nitrogen2:clip1 -0.12500 0.25840 -0.484
nitrogen3:clip1 0.15833 0.25840 0.613
weed1:clip1 0.09792 0.21098 0.464
weed2:clip1 -0.02708 0.21098 -0.128
nitrogen1:weed1:clip1 0.46042 0.36543 1.260
nitrogen2:weed1:clip1 -0.53125 0.36543 -1.454
nitrogen3:weed1:clip1 0.18542 0.36543 0.507
nitrogen1:weed2:clip1 -0.56458 0.36543 -1.545
nitrogen2:weed2:clip1 0.24375 0.36543 0.667
nitrogen3:weed2:clip1 0.18542 0.36543 0.507
Correlation matrix not shown by default, as p = 24 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need it
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
weed
is very significant (no surprise: add more weed seeds
get more weeds). clip
is highly significant (although not
that large an absolute effect). Nitrogen
is only marginally significant, but it is involved in a
very significant interaction.
We can see the relative sizes of the terms in the side by side plot.
> Anova(WB.lmer,test="F")
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: biomass
F Df Df.res Pr(>F)
nitrogen 14.5360 3 3 0.02718 *
weed 555.4545 2 8 2.613e-09 ***
clip 117.4291 1 12 1.494e-07 ***
nitrogen:weed 24.5815 6 8 9.664e-05 ***
nitrogen:clip 0.2293 3 12 0.87419
weed:clip 0.1149 2 12 0.89246
nitrogen:weed:clip 0.7514 6 12 0.62033
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> sidebyside(WB.lmer)

Figure 44.3: Side by side plot of effects in weed biomass data
Let’s take a look at that interaction. What we see is that non-weed biomass decreases as you add nitrogen, but it decreases much more steeply when you have added weed seeds.
> with(WeedBiomass,interactplot(nitrogen,weed,biomass))

Figure 44.4: Nitrogen by weed interaction plot for weed biomass data