QK Box 11.5

Hey et al. (2020) studied the effect of artificial light at night (ALAN) on the growth of a wildflower species (Asclepias syriaca) in an outdoor field experiment. The plant, a species of milkweed, is an important host plant for caterpillars of the monarch butterfly.

They used a split-plot design for their experiment. The between-plots factor was ALAN, with five plots receiving artificial light and five control plots with the same set-up but only receiving ambient light. The sub-plots were 11 L small pots and there were two within-plot factors arranged in a crossed design; planting density (one or three plants per pot) and soil moisture (weekly addition of water vs no watering). There were four pots in each plot, so each combination of density and moisture had a single pot; the design is unreplicated at the sub-plot level. The response variable was total biomass (root and shoot) per pant in each pot, averaged across the plants within the higher density pots.

Asclepias syriaca. Cbaile19, CC0, via Wikimedia Commons

The paper is here

Hey, M. H., DiBiase, E., Roach, D. A., Carr, D. E. & Haynes, K. J. (2020). Interactions between artificial light at night, soil moisture, and plant density affect the growth of a perennial wildflower. Oecologia, 193, 503-10.

Preliminaries

First, load the required packages (afex, car, lattice, lme4, lmerTest, nlme, VCA, ez, emmeans, Rmisc, MuMIn)

Import heymean data file (heymean.csv)

heymean <- read.csv("../data/heymean.csv")
heymean

Set contrasts from afex

Make plot a factor

set_sum_contrasts()

heymean$plot <- factor(heymean$plot)

Check residuals by leaving out error term

hey1.aov <- aov(dwtotal~light*density*water, heymean)
plot(hey1.aov)

See if transformation improves the residuals

hey2.aov <- aov(log10(dwtotal)~light*density*water, heymean)
plot(hey2.aov)

Some improvement in plot; outliers still present but transform response

Note Biomass is always positive, so we don’t really need to add anything before transforming. Hey et al. used a log(x+1) transformation, so we’ll do the same here for consistency. As an exercise, you could modify the next line of code to remove the +1 and rerun the model-fitting to see how much difference there is.

heymean$ldwtotal <- log10(heymean$dwtotal+1)

Fit factorial split-plot with separate within-plot error terms

hey3.aov <- aov(ldwtotal~light*density*water+Error(plot/(density*water)), heymean)
summary(hey3.aov)

Error: plot
          Df Sum Sq Mean Sq F value Pr(>F)
light      1 0.0121  0.0121    0.33   0.58
Residuals  8 0.2925  0.0366               

Error: plot:density
              Df Sum Sq Mean Sq F value Pr(>F)
density        1 0.0021 0.00212    0.21   0.66
light:density  1 0.0017 0.00172    0.17   0.69
Residuals      8 0.0820 0.01025               

Error: plot:water
            Df Sum Sq Mean Sq F value Pr(>F)  
water        1 0.0386  0.0386    4.38   0.07 .
light:water  1 0.0258  0.0258    2.93   0.13  
Residuals    8 0.0704  0.0088                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: plot:density:water
                    Df Sum Sq Mean Sq F value Pr(>F)
density:water        1 0.0029 0.00289    0.29   0.60
light:density:water  1 0.0211 0.02112    2.12   0.18
Residuals            8 0.0796 0.00995               

Fit equivalent mixed effects model

hey1.lmer <- lmer(ldwtotal~light+density+water+light*density+light*water+density*water+light*density*water+(1|plot)+(1|plot:density)+(1|plot:water), REML=TRUE, heymean)

boundary (singular) fit: see help('isSingular')

plot(hey1.lmer)

residuals OK

summary(hey1.lmer, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: ldwtotal ~ light + density + water + light * density + light *      water + density * water + light * density * water + (1 |      plot) + (1 | plot:density) + (1 | plot:water)
   Data: heymean

REML criterion at convergence: -17.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.9570 -0.5011  0.0056  0.2921  2.4266 

Random effects:
 Groups       Name        Variance Std.Dev.
 plot:water   (Intercept) 0.000000 0.0000  
 plot:density (Intercept) 0.000436 0.0209  
 plot         (Intercept) 0.006579 0.0811  
 Residual                 0.009376 0.0968  
Number of obs: 40, groups:  plot:water, 20; plot:density, 20; plot, 10

Fixed effects:
                       Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)             0.25082    0.03023  8.00000    8.30  3.4e-05 ***
light1                 -0.01736    0.03023  8.00000   -0.57    0.582    
density1                0.00728    0.01601  8.00000    0.45    0.661    
water1                 -0.03105    0.01531  8.00000   -2.03    0.077 .  
light1:density1         0.00655    0.01601  8.00000    0.41    0.693    
light1:water1          -0.02539    0.01531  8.00000   -1.66    0.136    
density1:water1         0.00850    0.01531  8.00000    0.56    0.594    
light1:density1:water1  0.02298    0.01531  8.00000    1.50    0.172    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) light1 dnsty1 water1 lght1:d1 lght1:w1 dns1:1
light1      0.000                                               
density1    0.000  0.000                                        
water1      0.000  0.000  0.000                                 
lght1:dnst1 0.000  0.000  0.000  0.000                          
light1:wtr1 0.000  0.000  0.000  0.000  0.000                   
dnsty1:wtr1 0.000  0.000  0.000  0.000  0.000    0.000          
lght1:dn1:1 0.000  0.000  0.000  0.000  0.000    0.000    0.000 
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

anova(hey1.lmer, ddf="Kenward-Roger")

Type III Analysis of Variance Table with Kenward-Roger's method
                    Sum Sq Mean Sq NumDF DenDF F value Pr(>F)  
light               0.0031  0.0031     1     8    0.33  0.582  
density             0.0019  0.0019     1     8    0.21  0.661  
water               0.0386  0.0386     1     8    4.11  0.077 .
light:density       0.0016  0.0016     1     8    0.17  0.693  
light:water         0.0258  0.0258     1     8    2.75  0.136  
density:water       0.0029  0.0029     1     8    0.31  0.594  
light:density:water 0.0211  0.0211     1     8    2.25  0.172  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get variance components

hey1.ci <- confint.merMod(hey1.lmer, oldNames=FALSE)

Computing profile confidence intervals ...
Warning: Last two rows have identical or NA .zeta values: using minstepWarning: Last two rows have identical or NA .zeta values: using minstepWarning: non-monotonic profile for sd_(Intercept)|plot:waterWarning: bad spline fit for sd_(Intercept)|plot:water: falling back to linear interpolationWarning: collapsing to unique 'x' values

hey1.vc <- (hey1.ci)^2
print(hey1.vc)

                               2.5 %   97.5 %
sd_(Intercept)|plot:water   0.00e+00 6.34e-03
sd_(Intercept)|plot:density 0.00e+00 6.91e-03
sd_(Intercept)|plot         3.08e-04 1.84e-02
sigma                       4.37e-03 1.33e-02
(Intercept)                 3.71e-02 9.56e-02
light1                      5.73e-03 1.68e-03
density1                    5.62e-04 1.46e-03
water1                      3.49e-03 9.04e-06
light1:density1             5.97e-04 1.41e-03
light1:water1               2.86e-03 7.03e-06
density1:water1             3.82e-04 1.34e-03
light1:density1:water1      2.57e-05 2.60e-03

Compare this model to model with no plot by density and moisture interactions

hey1a.lmer <- lmer(ldwtotal~light+density+water+light*density+light*water+density*water+light*density*water+(1|plot)+(1|plot:density)+(1|plot:water), REML=FALSE, heymean)

boundary (singular) fit: see help('isSingular')

hey2.lmer <- lmer(ldwtotal~light*density*water+(1|plot), REML=FALSE, heymean)
anova(hey1a.lmer, hey2.lmer)

Data: heymean
Models:
hey2.lmer: ldwtotal ~ light * density * water + (1 | plot)
hey1a.lmer: ldwtotal ~ light + density + water + light * density + light * water + density * water + light * density * water + (1 | plot) + (1 | plot:density) + (1 | plot:water)
           npar   AIC   BIC logLik deviance Chisq Df Pr(>Chisq)
hey2.lmer    10 -47.7 -30.8   33.8    -67.7                    
hey1a.lmer   12 -43.7 -23.4   33.8    -67.7  0.03  2       0.99

AICc(hey1a.lmer, hey2.lmer)

Simpler model fits equally well

hey2a.lmer <- lmer(ldwtotal~light*density*water+(1|plot), REML=TRUE, heymean)
summary(hey2a.lmer, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: ldwtotal ~ light * density * water + (1 | plot)
   Data: heymean

REML criterion at convergence: -17.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.0042 -0.4504 -0.0212  0.2962  2.4499 

Random effects:
 Groups   Name        Variance Std.Dev.
 plot     (Intercept) 0.00672  0.0820  
 Residual             0.00967  0.0983  
Number of obs: 40, groups:  plot, 10

Fixed effects:
                       Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)             0.25082    0.03023  8.00000    8.30  3.4e-05 ***
light1                 -0.01736    0.03023  8.00000   -0.57    0.582    
density1                0.00728    0.01555 24.00000    0.47    0.644    
water1                 -0.03105    0.01555 24.00000   -2.00    0.057 .  
light1:density1         0.00655    0.01555 24.00000    0.42    0.677    
light1:water1          -0.02539    0.01555 24.00000   -1.63    0.115    
density1:water1         0.00850    0.01555 24.00000    0.55    0.590    
light1:density1:water1  0.02298    0.01555 24.00000    1.48    0.152    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) light1 dnsty1 water1 lght1:d1 lght1:w1 dns1:1
light1      0.000                                               
density1    0.000  0.000                                        
water1      0.000  0.000  0.000                                 
lght1:dnst1 0.000  0.000  0.000  0.000                          
light1:wtr1 0.000  0.000  0.000  0.000  0.000                   
dnsty1:wtr1 0.000  0.000  0.000  0.000  0.000    0.000          
lght1:dn1:1 0.000  0.000  0.000  0.000  0.000    0.000    0.000 

anova(hey2a.lmer, ddf="Kenward-Roger")

Type III Analysis of Variance Table with Kenward-Roger's method
                    Sum Sq Mean Sq NumDF DenDF F value Pr(>F)  
light               0.0032  0.0032     1     8    0.33  0.582  
density             0.0021  0.0021     1    24    0.22  0.644  
water               0.0386  0.0386     1    24    3.99  0.057 .
light:density       0.0017  0.0017     1    24    0.18  0.677  
light:water         0.0258  0.0258     1    24    2.67  0.115  
density:water       0.0029  0.0029     1    24    0.30  0.590  
light:density:water 0.0211  0.0211     1    24    2.18  0.152  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Compare to factorial OLS split-plot pooling within-plot error terms

Same as above for fixed effects

hey4.aov <- aov(ldwtotal~light*density*water+Error(plot), heymean)
summary(hey4.aov)

Error: plot
          Df Sum Sq Mean Sq F value Pr(>F)
light      1 0.0121  0.0121    0.33   0.58
Residuals  8 0.2925  0.0366               

Error: Within
                    Df Sum Sq Mean Sq F value Pr(>F)  
density              1 0.0021  0.0021    0.22  0.644  
water                1 0.0386  0.0386    3.99  0.057 .
light:density        1 0.0017  0.0017    0.18  0.677  
light:water          1 0.0258  0.0258    2.67  0.115  
density:water        1 0.0029  0.0029    0.30  0.590  
light:density:water  1 0.0211  0.0211    2.18  0.152  
Residuals           24 0.2320  0.0097                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1