Q & K Box 10.4

Lozada-Misa et al. (2015) compared two different morphology types of the coral genus Porites for lesions caused by white syndrome disease. They collected ten random colonies of the branching P. cylindrica and ten colonies of the massive Porites spp., and five random lesions per colony were photographed and measured. Morphology type was a fixed factor, with colony a random factor nested within morphology type, and size of individual lesions as the observations.

Branching coral. Mick Keough, CC BY 4.0.

Massive coral. Emmet Keough, CC BY4.0.

Massive coral. Liam Keough, CC BY-SA 4.0

Coral surface showing lesions. Mick Keough CC BY 4.0

The paper is here

Lozada-Misa, P., Kerr, A. M. & Raymundo, L. (2015). Contrasting lesion dynamics of white syndrome among the scleractinian corals Porites spp. PLoS One, 10, e0129841.

Preliminaries

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

Import lesions data file (lesions.csv)

lesions <- read.csv("../data/lesions.csv")
head(lesions, 10)

Set contrasts from afex

Make colony a factor

set_sum_contrasts()

lesions$col <-factor(lesions$col)

Calculate summary stats by groups (to be used later)

lesions_sum <- summarySE(lesions,measurevar='size', groupvars= c('col','morph'))

Do preliminary checks

Run for all data and summary data with means for each colony

boxplot(size~morph, data=lesions)

boxplot(size~morph, data=lesions_sum)

boxplot(log10(size)~morph, data=lesions)

boxplot(log10(size)~morph, data=lesions_sum)

Size distributions look better with log10 transform

lesions$lsize <- log10(lesions$size)

Fit nested anova model

lesions.aov <- aov(lsize~morph+Error(col), lesions)

Plot residuals based on colonies

plot(resid(lesions.aov[[2]])~fitted(lesions.aov[[2]]))

Get anova summary

options(digits = 3)
tidy(lesions.aov)

To get a test for differences among colonies, construct an F-ratio using the colony variation (“Error:col Residuals) and the within-colony term. Use pf to get p value

f <- 0.3041/0.2137
f
1-pf(f, 18, 80, lower.tail = TRUE, log.p = FALSE)

Demonstrate that the nested analysis is the same as one-way on colony means

This equality holds as long as the design is balanced. When there are unequal sample sizes within the nested factor (lesions for each colony in this example), the results will be very similar, though not identical

lesions_sum <- summarySE(lesions,measurevar='lsize', groupvars= c('col','morph'))
lesions.aov1 <- aov(lsize~morph, lesions_sum)
tidy(lesions.aov1)

using VCA package to get anova variance components

This VCs are estimated from an OLS analysis, so it’s possible that lower CIs can be negative

lesions.vca <- anovaMM(lsize~morph/(col), lesions)

Convert variable morph from "character" to "factor"!

lesions.vca

ANOVA-Type Estimation of Mixed Model:
--------------------------------------

    [Fixed Effects]

           int morphBRANCHING   morphMASSIVE 
         1.232         -0.627          0.000 


    [Variance Components]

  Name      DF        SS        MS       VC       %Total    SD       CV[%]    
1 total     94.107336                    0.231766 100       0.481421 52.388834
2 morph:col 18        5.473944  0.304108 0.018085 7.803304  0.134482 14.634504
3 error     80        17.094469 0.213681 0.213681 92.196696 0.462256 50.303293

Mean: 0.919 (N = 100) 

Experimental Design: balanced  |  Method: ANOVA

VCAinference(lesions.vca, alpha=0.05, VarVC=TRUE, ci.method="satterthwaite")

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

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

    [Fixed Effects]

           int morphBRANCHING   morphMASSIVE 
         1.232         -0.627          0.000 


    [Variance Components]

  Name      DF      SS      MS     VC     %Total  SD     CV[%]   Var(VC)
1 total     94.1073                0.2318 100     0.4814 52.3888        
2 morph:col 18      5.4739  0.3041 0.0181 7.8033  0.1345 14.6345 5e-04  
3 error     80      17.0945 0.2137 0.2137 92.1967 0.4623 50.3033 0.0011 

Mean: 0.919 (N = 100) 

Experimental Design: balanced  |  Method: ANOVA


> VC:
-----
          Estimate    DF CI LCL CI UCL One-Sided LCL One-Sided UCL
total       0.2318 94.11 0.1776  0.315        0.1852         0.300
morph:col   0.0181  1.43 0.0042  2.535        0.0054         0.958
error       0.2137 80.00 0.1603  0.299        0.1678         0.283

> SD:
-----
          Estimate    DF CI LCL CI UCL One-Sided LCL One-Sided UCL
total        0.481 94.11 0.4214  0.562        0.4304         0.548
morph:col    0.135  1.43 0.0651  1.592        0.0732         0.979
error        0.462 80.00 0.4004  0.547        0.4096         0.532

> CV[%]:
--------
          Estimate    DF CI LCL CI UCL One-Sided LCL One-Sided UCL
total         52.4 94.11  45.85   61.1         46.83          59.6
morph:col     14.6  1.43   7.08  173.3          7.97         106.5
error         50.3 80.00  43.57   59.5         44.58          57.9


95% Confidence Level  
Satterthwaite methodology used for computing CIs 

Mixed model using lme4

This is an alternative approach where we specify fixed and random effects explicitly, and assess the model with a combination of REML and ML.

lesions.reml <- lmer(lsize~morph + (1|col), lesions)
summary(lesions.reml, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: lsize ~ morph + (1 | col)
   Data: lesions

REML criterion at convergence: 142

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.2303 -0.5402 -0.0075  0.4923  3.0487 

Random effects:
 Groups   Name        Variance Std.Dev.
 col      (Intercept) 0.0181   0.134   
 Residual             0.2137   0.462   
Number of obs: 100, groups:  col, 20

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)   0.9189     0.0551 18.0000   16.66  2.2e-12 ***
morph1       -0.3134     0.0551 18.0000   -5.68  2.2e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
       (Intr)
morph1 0.000 

anova(lesions.reml, ddf="Kenward-Roger")

Type III Analysis of Variance Table with Kenward-Roger's method
      Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)    
morph    6.9     6.9     1    18    32.3 2.2e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

lesions.ml <- lmer(lsize~morph + (1|col), lesions, REML=F)
summary(lesions.ml)

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: lsize ~ morph + (1 | col)
   Data: lesions

     AIC      BIC   logLik deviance df.resid 
   142.4    152.8    -67.2    134.4       96 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.1834 -0.5357 -0.0052  0.4999  3.0956 

Random effects:
 Groups   Name        Variance Std.Dev.
 col      (Intercept) 0.012    0.110   
 Residual             0.214    0.462   
Number of obs: 100, groups:  col, 20

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)   0.9189     0.0523 20.0000   17.57  1.3e-13 ***
morph1       -0.3134     0.0523 20.0000   -5.99  7.4e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
       (Intr)
morph1 0.000 

CI on variance components (remembering to square CIs from lmer with are in SD units)

lesions.ci <- confint.merMod(lesions.reml)

Computing profile confidence intervals ...

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

            2.5 % 97.5 %
.sig01      0.000 0.0675
.sigma      0.159 0.2955
(Intercept) 0.658 1.0539
morph1      0.177 0.0423

Test fixed effect using likelihood ratio tests based on ML

lesions.lme <- lmer(lsize~morph + (1|col), lesions, REML=F)
AICc(lesions.lme)

[1] 143

lesions.lme1 <- lmer(lsize~(1|col), lesions, REML=F)
AICc(lesions.lme1)

[1] 161

anova(lesions.lme,lesions.lme1)

Data: lesions
Models:
lesions.lme1: lsize ~ (1 | col)
lesions.lme: lsize ~ morph + (1 | col)
             npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)    
lesions.lme1    3 161 169  -77.5      155                        
lesions.lme     4 142 153  -67.2      134  20.6  1    5.8e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1