Chapter 45 Split-Split Plots

45.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)

45.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)
Residual plot for weed biomass data

Figure 45.1: Residual plot for weed biomass data

> qqnorm(residuals(WB.lmer))
NPP of residuals for weed biomass data

Figure 45.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)
Side by side plot of effects in weed biomass data

Figure 45.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))
Nitrogen by weed interaction plot for weed biomass data

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