QK Box 10.11

Van der Geest et al (2020) set up a field experiment to examine the mutualistic relationship between seagrasses (Zostera noltei) and lucinid bivalves (Loripes orbiculatus), in particular the extent to which the bivalves and their endosymbiotic bacteria could oxidize harmful sulfides and prevent their accumulation by seagrasses. In a coastal lagoon in the western Mediterranean, they established two fixed factors: presence or absence of bivalves (added vs control/background levels) and organic matter (added to increase pore water sulfide production vs control), with a single experimental plot for each combination replicated across six randomly chosen locations (blocks). The response variable we will analyze was the percentage of sulfur in plant tissue originating from the sediment (F_Sulfide).

You can find an image of the seagrass Zostera noltei here, but if you just would like to see what a typical Zostera looks like, here’s Zostera nigricaulis from southern Australia.

Mick Keough, CC BY 4.0

Lucinid bivalve, in this case Loripes aberrans. Muséum national d’histoire naturelle, CC BY 4.0

Preliminaries

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

Import vandergeest data file (vandergeest.csv)

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

NA

set contrasts from afex

set_sum_contrasts()

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

Fit OLS anova model to check residuals

vandergeest.aov <- aov(fsulfide~om+bivalve+block+om*bivalve+om*block+bivalve*block, data=vandergeest)

Check residuals - variance increases with mean

plot(vandergeest.aov)

residuals look fine

Fit models to get type III SS using car

options(digits = 3)
vandergeest.lm <- lm(aov(fsulfide~om+bivalve+block+om*bivalve+om*block+bivalve*block, data=vandergeest))
Anova(vandergeest.lm, type='III')

Anova Table (Type III tests)

Response: fsulfide
              Sum Sq Df F value  Pr(>F)    
(Intercept)    11350  1 1202.12 3.8e-07 ***
om               101  1   10.67   0.022 *  
bivalve           25  1    2.66   0.164    
block             19  5    0.39   0.835    
om:bivalve         5  1    0.54   0.497    
om:block          24  5    0.51   0.763    
bivalve:block     23  5    0.50   0.770    
Residuals         47  5                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get correct F-ratio and P value for om: F = 100.711/4.785. These are the mean squares for om and om*block

1-pf(100.711/4.785, 1, 5, lower.tail = TRUE, log.p = FALSE)

[1] 0.00591

Get correct F-ratio and P value for bivalve 25.122/4.681

1-pf(25.122/4.681, 1, 5, lower.tail = TRUE, log.p = FALSE)

[1] 0.0683

Get main effect means

emmeans(aov(fsulfide~om+bivalve+block+om*bivalve+om*block+bivalve*block, data=vandergeest), "om")

NOTE: Results may be misleading due to involvement in interactions

 om emmean    SE df lower.CL upper.CL
 C    19.7 0.887  5     17.4     22.0
 OM   23.8 0.887  5     21.5     26.1

Results are averaged over the levels of: bivalve, block 
Confidence level used: 0.95 

emmeans(aov(fsulfide~om+bivalve+block+om*bivalve+om*block+bivalve*block, data=vandergeest), "bivalve")

NOTE: Results may be misleading due to involvement in interactions

 bivalve emmean    SE df lower.CL upper.CL
 Added     20.7 0.887  5     18.4       23
 C         22.8 0.887  5     20.5       25

Results are averaged over the levels of: om, block 
Confidence level used: 0.95 

Using VCA package to get anova var comps

vandergeest.vca <- anovaMM(fsulfide~om+bivalve+om*bivalve + (block) + (om*block) + (bivalve*block) , vandergeest)

Convert variable block from "character" to "factor"!
Convert variable bivalve from "character" to "factor"!
Convert variable om from "character" to "factor"!

VCAinference(vandergeest.vca, alpha=0.05, VarVC=TRUE, excludeNeg=FALSE, constrainCI=FALSE)

Inference from Mixed Model Fit
------------------------------

> VCA Result:
-------------

    [Fixed Effects]

              int               omC              omOM      bivalveAdded          bivalveC  omC:bivalveAdded omOM:bivalveAdded      omC:bivalveC     omOM:bivalveC 
            25.28             -5.02              0.00             -2.97              0.00              1.84              0.00              0.00              0.00 


    [Variance Components]

  Name          DF      SS      MS     VC      %Total  SD     CV[%]   Var(VC)
1 total         19.5352                10.3647 100     3.2194 14.8042        
2 block         5       18.5806 3.7161 0.923   8.9052  0.9607 4.4178  3.694  
3 om:block      5       23.9258 4.7852 0*      0*      0*     0*      11.2043
4 bivalve:block 5       23.4033 4.6807 0*      0*      0*     0*      11.1054
5 error         5       47.2085 9.4417 9.4417  91.0948 3.0727 14.1297 35.6582

Mean: 21.7 (N = 24) 

Experimental Design: balanced  |  Method: ANOVA | * VC set to 0 | adapted MS used for total DF


> VC:
-----
              Estimate CI LCL CI UCL  One-Sided LCL One-Sided UCL
total         10.3647  6.0336 21.838  6.5691        19.2672      
block         0.923    0*     4.69    0*            4.0844       
om:block      0        0*     6.5606  0*            5.5058       
bivalve:block 0        0*     6.5315  0*            5.4814       
error         9.4417   3.6788 56.7948 4.2643        41.213       

> SD:
-----
              Estimate CI LCL CI UCL One-Sided LCL One-Sided UCL
total         3.2194   2.4563 4.6731 2.563         4.3894       
block         0.9607   0*     2.1656 0*            2.021        
om:block      0        0*     2.5614 0*            2.3464       
bivalve:block 0        0*     2.5557 0*            2.3412       
error         3.0727   1.918  7.5362 2.065         6.4197       

> CV[%]:
--------
              Estimate CI LCL  CI UCL  One-Sided LCL One-Sided UCL
total         14.8042  11.2952 21.4889 11.7858       20.1844      
block         4.4178   0*      9.9585  0*            9.2933       
om:block      0        0*      11.7782 0*            10.7899      
bivalve:block 0        0*      11.7521 0*            10.766       
error         14.1297  8.8198  34.6546 9.4958        29.5205      


95% Confidence Level  | * CI-limits constrained to be >= 0
SAS PROC MIXED method used for computing CIs 

Fit mixed effects model using REML

vandergeest.lmer1 <- lmer(fsulfide~om+bivalve+om*bivalve + (1|block) + (1|om:block)
                     + (1|bivalve:block) , REML=TRUE, vandergeest)

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

summary(vandergeest.lmer1, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: fsulfide ~ om + bivalve + om * bivalve + (1 | block) + (1 | om:block) +      (1 | bivalve:block)
   Data: vandergeest

REML criterion at convergence: 104

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.850 -0.567 -0.096  0.723  1.774 

Random effects:
 Groups        Name        Variance Std.Dev.
 bivalve:block (Intercept) 0.00     0.00    
 om:block      (Intercept) 0.00     0.00    
 block         (Intercept) 0.00     0.00    
 Residual                  5.66     2.38    
Number of obs: 24, groups:  bivalve:block, 12; om:block, 12; block, 6

Fixed effects:
             Estimate Std. Error     df t value Pr(>|t|)    
(Intercept)    21.747      0.485  5.000   44.80    1e-07 ***
om1            -2.048      0.485  5.000   -4.22   0.0083 ** 
bivalve1       -1.023      0.485  5.000   -2.11   0.0889 .  
om1:bivalve1    0.460      0.485  5.000    0.95   0.3872    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) om1   bivlv1
om1         0.000              
bivalve1    0.000  0.000       
om1:bivalv1 0.000  0.000 0.000 
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

anova(vandergeest.lmer1, ddf="Kenward-Roger")

Type III Analysis of Variance Table with Kenward-Roger's method
           Sum Sq Mean Sq NumDF DenDF F value Pr(>F)   
om          100.7   100.7     1     5   17.81 0.0083 **
bivalve      25.1    25.1     1     5    4.44 0.0889 . 
om:bivalve    5.1     5.1     1     5    0.90 0.3872   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Refit without 2-way interactions to solve singularity - unsuccessful

vandergeest.lmer2 <- lmer(fsulfide~om+bivalve+om*bivalve + (1|block), REML=TRUE, vandergeest)

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

summary(vandergeest.lmer2, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: fsulfide ~ om + bivalve + om * bivalve + (1 | block)
   Data: vandergeest

REML criterion at convergence: 104

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.850 -0.567 -0.096  0.723  1.774 

Random effects:
 Groups   Name        Variance Std.Dev.
 block    (Intercept) 0.00     0.00    
 Residual             5.66     2.38    
Number of obs: 24, groups:  block, 6

Fixed effects:
             Estimate Std. Error     df t value Pr(>|t|)    
(Intercept)    21.747      0.485  5.000   44.80    1e-07 ***
om1            -2.048      0.485 15.000   -4.22  0.00074 ***
bivalve1       -1.023      0.485 15.000   -2.11  0.05231 .  
om1:bivalve1    0.460      0.485 15.000    0.95  0.35867    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) om1   bivlv1
om1         0.000              
bivalve1    0.000  0.000       
om1:bivalv1 0.000  0.000 0.000 
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

anova(vandergeest.lmer2, ddf="Kenward-Roger")

Type III Analysis of Variance Table with Kenward-Roger's method
           Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)    
om          100.7   100.7     1    15   17.81 0.00074 ***
bivalve      25.1    25.1     1    15    4.44 0.05231 .  
om:bivalve    5.1     5.1     1    15    0.90 0.35867    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1