QK Box 11.2

Legault et al. (2018) experimentally examined the effects of nutrient loading on the outcome of competitive interactions between two species of saltmarsh plants, the native Spartina alterniflora and the exotic Phragmites australis, in northeast USA. The between-plots factor was nutrient addition (fixed, with three groups: none, low and high) with five 50-gallon bins containing half-strength seawater as the plots. The within-plots factor was competition (fixed with two groups: no competition, with competition), with sub-plots being small pots containing one (no competition) or both (with competition) species within each bin. This design was replicated at the sub-plot level with two pots for each competition group within each bin. We will analyze biomass of P. australis at ambient temperature as the response variable.

Phragmites australis. Mick Keough, CC BY 4.0

The paper is here

Legault, R., Zogg, G. P. & Travis, S. E. (2018). Competitive interactions between native Spartina alterniflora and non-native Phragmites australis depend on nutrient loading and temperature. PLoS One, 13, e0192234.

Preliminaries

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

Import legault data file (legault.csv)

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

Set contrasts from afex

Re-arrange nutrient levels so that no nutrients is reference

set_sum_contrasts()

setting contr.sum globally: options(contrasts=c('contr.sum', 'contr.poly'))

legault$nutrient <- factor(legault$nutrient, levels=c("Null (0 g N/m2/year)","Low (30 g N/m2/year)","High (120 g N/m2/year)"))
legault$comp<-factor(legault$comp)

Fit aov crossed model to check approx residuals

legault1.aov <- aov(biom~nutrient+comp+nutrient*comp, legault)
plot(legault1.aov)

Some evidence for unequal residuals but use untransformed biomass to match original paper

Fit model using aov (balanced design so all SS types the same)

legault2.aov <- aov(biom~nutrient+Error(bin/comp)+comp+nutrient*comp, legault)
summary(legault2.aov)

Error: bin
          Df Sum Sq Mean Sq F value   Pr(>F)    
nutrient   2  15071    7535   46.45 2.24e-06 ***
Residuals 12   1947     162                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: bin:comp
              Df Sum Sq Mean Sq F value  Pr(>F)   
comp           1 1400.4  1400.4   9.726 0.00888 **
nutrient:comp  2  285.4   142.7   0.991 0.39966   
Residuals     12 1727.9   144.0                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: Within
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals 30   2314   77.12               

print(summary(legault2.aov))

Error: bin
          Df Sum Sq Mean Sq F value   Pr(>F)    
nutrient   2  15071    7535   46.45 2.24e-06 ***
Residuals 12   1947     162                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: bin:comp
              Df Sum Sq Mean Sq F value  Pr(>F)   
comp           1 1400.4  1400.4   9.726 0.00888 **
nutrient:comp  2  285.4   142.7   0.991 0.39966   
Residuals     12 1727.9   144.0                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: Within
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals 30   2314   77.12               

Fit split-plot anova model using ez for comparison of output

This uses mean biomass for each competition group within each bin - same answers as above due to balanced design

legault.ez <- ezANOVA(dv=biom, wid=bin, within=comp, between=nutrient, type=1, legault)

Warning: Converting "bin" to factor for ANOVA.Warning: Collapsing data to cell means. *IF* the requested effects are a subset of the full design, you must use the "within_full" argument, else results may be inaccurate.

print(legault.ez)

$ANOVA
         Effect DFn DFd          F            p p<.05        ges
1      nutrient   2  12 46.4539769 2.239943e-06     * 0.80398143
2          comp   1  12  9.7258295 8.877441e-03     * 0.27595248
3 nutrient:comp   2  12  0.9909378 3.996639e-01       0.07206656

Rerun as classic split-plot design based on sub-plot means - matches above

This example illustrates a short-cut that produces a simpler analysis for a paper, etc.

legault_means <- summarySE(legault, measurevar= 'biom', groupvars= c('nutrient','bin','comp'))
legault_means

legault3.aov <- aov(biom~nutrient*comp+Error(bin), legault_means)
summary(legault3.aov)

Error: bin
          Df Sum Sq Mean Sq F value   Pr(>F)    
nutrient   2   7535    3768   46.45 2.24e-06 ***
Residuals 12    973      81                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: Within
              Df Sum Sq Mean Sq F value  Pr(>F)   
comp           1  700.2   700.2   9.726 0.00888 **
nutrient:comp  2  142.7    71.3   0.991 0.39966   
Residuals     12  863.9    72.0                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Fit mixed model using lme4

legault1.lmer <- lmer(biom~nutrient+comp+nutrient*comp+(1|bin)+(1|bin:comp), REML=TRUE, legault)

Get residuals from lmer model fit - some evidence on variance heterogeneity but again use untransformed biom to match original paper

plot(legault1.lmer)

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

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: biom ~ nutrient + comp + nutrient * comp + (1 | bin) + (1 | bin:comp)
   Data: legault

REML criterion at convergence: 426.7

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.79289 -0.47345 -0.04004  0.51095  1.96250 

Random effects:
 Groups   Name        Variance Std.Dev.
 bin:comp (Intercept) 33.432   5.782   
 bin      (Intercept)  4.556   2.134   
 Residual             77.124   8.782   
Number of obs: 60, groups:  bin:comp, 30; bin, 15

Fixed effects:
                Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)       30.267      1.644  12.000  18.408 3.67e-10 ***
nutrient1        -21.426      2.325  12.000  -9.214 8.61e-07 ***
nutrient2          5.014      2.325  12.000   2.156  0.05208 .  
comp1              4.831      1.549  12.000   3.119  0.00888 ** 
nutrient1:comp1   -3.071      2.191  12.000  -1.402  0.18629    
nutrient2:comp1    1.781      2.191  12.000   0.813  0.43213    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) ntrnt1 ntrnt2 comp1  ntr1:1
nutrient1    0.000                            
nutrient2    0.000 -0.500                     
comp1        0.000  0.000  0.000              
ntrnt1:cmp1  0.000  0.000  0.000  0.000       
ntrnt2:cmp1  0.000  0.000  0.000  0.000 -0.500

anova(legault1.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)    
nutrient      7165.4  3582.7     2    12 46.4540 2.24e-06 ***
comp           750.1   750.1     1    12  9.7258 0.008877 ** 
nutrient:comp  152.8    76.4     2    12  0.9909 0.399664    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Contrast no nutrients with low and high separately

legault1.emm <- emmeans(legault1.lmer, ~nutrient)

NOTE: Results may be misleading due to involvement in interactions

legault1.con <- contrast(legault1.emm, "trt.vs.ctrl", ref=1)
summary(legault1.con, adjust="none")

 contrast                                          estimate   SE df t.ratio p.value
 (Low (30 g N/m2/year)) - (Null (0 g N/m2/year))       26.4 4.03 12   6.565  <.0001
 (High (120 g N/m2/year)) - (Null (0 g N/m2/year))     37.8 4.03 12   9.395  <.0001

Results are averaged over the levels of: comp 
Degrees-of-freedom method: kenward-roger 

Get variance components with CIs

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

Computing profile confidence intervals ...

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

                              2.5 %       97.5 %
sd_(Intercept)|bin:comp   0.0000000   69.8785882
sd_(Intercept)|bin        0.0000000   46.7902117
sigma                    48.3597992  132.4605475
(Intercept)             739.2584171 1111.8018654
nutrient1               664.4696090  291.5109541
nutrient2                 0.4378363   87.7089963
comp1                     3.8408980   59.3286228
nutrient1:comp1          50.8634666    0.9791549
nutrient2:comp1           5.1977743   34.1233859

Evaluate whether to omit bins(nutrients) x competition interaction

legault1a.lmer <- lmer(biom~nutrient+comp+nutrient*comp+(1|bin)+(1|bin:comp), REML=FALSE, legault)
legault2.lmer <- lmer(biom~nutrient+comp+nutrient*comp+(1|bin), REML=FALSE, legault)
anova(legault1a.lmer,legault2.lmer)

Data: legault
Models:
legault2.lmer: biom ~ nutrient + comp + nutrient * comp + (1 | bin)
legault1a.lmer: biom ~ nutrient + comp + nutrient * comp + (1 | bin) + (1 | bin:comp)
               npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
legault2.lmer     8 461.66 478.41 -222.83   445.66                     
legault1a.lmer    9 462.82 481.67 -222.41   444.82 0.8365  1     0.3604

AICc(legault1a.lmer,legault2.lmer)

Get results from model omiting bins(nutrients) x competition interaction

legault2a.lmer <- lmer(biom~nutrient+comp+nutrient*comp+(1|bin), REML=TRUE, legault)
summary(legault2a.lmer)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: biom ~ nutrient + comp + nutrient * comp + (1 | bin)
   Data: legault

REML criterion at convergence: 428.5

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.87288 -0.46690 -0.01945  0.52482  2.17844 

Random effects:
 Groups   Name        Variance Std.Dev.
 bin      (Intercept) 16.50    4.062   
 Residual             96.23    9.810   
Number of obs: 60, groups:  bin, 15

Fixed effects:
                Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)       30.267      1.644  12.000  18.408 3.67e-10 ***
nutrient1        -21.425      2.325  12.000  -9.214 8.61e-07 ***
nutrient2          5.014      2.325  12.000   2.156 0.052077 .  
comp1              4.831      1.266  42.000   3.815 0.000441 ***
nutrient1:comp1   -3.071      1.791  42.000  -1.715 0.093753 .  
nutrient2:comp1    1.781      1.791  42.000   0.994 0.325754    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) ntrnt1 ntrnt2 comp1  ntr1:1
nutrient1    0.000                            
nutrient2    0.000 -0.500                     
comp1        0.000  0.000  0.000              
ntrnt1:cmp1  0.000  0.000  0.000  0.000       
ntrnt2:cmp1  0.000  0.000  0.000  0.000 -0.500

anova(legault2a.lmer, type="3", ddf="Kenward-Roger")

Type III Analysis of Variance Table with Kenward-Roger's method
              Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
nutrient      8940.4  4470.2     2    12 46.4540  2.24e-06 ***
comp          1400.4  1400.4     1    42 14.5530 0.0004407 ***
nutrient:comp  285.4   142.7     2    42  1.4828 0.2386491    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Fit mixed model on means (classic split-plot design) using lmer

legault3.lmer <- lmer(biom~nutrient+comp+nutrient*comp+(1|bin), REML=TRUE, legault_means)
summary(legault3.lmer)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: biom ~ nutrient + comp + nutrient * comp + (1 | bin)
   Data: legault_means

REML criterion at convergence: 190.4

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.85116 -0.31715 -0.06402  0.09750  2.64333 

Random effects:
 Groups   Name        Variance Std.Dev.
 bin      (Intercept)  4.556   2.134   
 Residual             71.994   8.485   
Number of obs: 30, groups:  bin, 15

Fixed effects:
                Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)       30.267      1.644  12.000  18.408 3.67e-10 ***
nutrient1        -21.425      2.325  12.000  -9.214 8.61e-07 ***
nutrient2          5.013      2.325  12.000   2.156  0.05208 .  
comp1              4.831      1.549  12.000   3.119  0.00888 ** 
nutrient1:comp1   -3.071      2.191  12.000  -1.402  0.18629    
nutrient2:comp1    1.781      2.191  12.000   0.813  0.43213    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) ntrnt1 ntrnt2 comp1  ntr1:1
nutrient1    0.000                            
nutrient2    0.000 -0.500                     
comp1        0.000  0.000  0.000              
ntrnt1:cmp1  0.000  0.000  0.000  0.000       
ntrnt2:cmp1  0.000  0.000  0.000  0.000 -0.500

anova(legault3.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)    
nutrient      6688.9  3344.4     2    12 46.4540 2.24e-06 ***
comp           700.2   700.2     1    12  9.7258 0.008877 ** 
nutrient:comp  142.7    71.3     2    12  0.9909 0.399664    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get variance components with CIs

legault3.ci <- confint.merMod(legault3.lmer)

Computing profile confidence intervals ...

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

                     2.5 %       97.5 %
.sig01            0.000000   46.7886429
.sigma           30.382596  105.9913106
(Intercept)     739.258090 1111.8022773
nutrient1       664.470048  291.5106558
nutrient2         0.437825   87.7091600
comp1             3.841433   59.3264999
nutrient1:comp1  50.860716    0.9787693
nutrient2:comp1   5.196895   34.1211091