QK Box 10.7

We will analyse a modified dataset from Schlegel et al. (2012) who examined swimming speed of sea urchin sperm under three different seawater pH levels [pH 8.1(control) and two more acidic treatments of pH 7.8 and 7.5] and from 19 different randomly chosen individual animals. There were between nine and ten runs per individual and pH combination and the response variable was the average swimming speed of sperm.

Sea urchin, Heliocidaris erythrogramma. Mick Keough, CC BY 4.0

The paper is here

Schlegel, P., Havenhand, J. N., Gillings, M. R. & Williamson, J. E. (2012). Individual variability in reproductive success determines winners and losers under ocean acidification: a case study with sea urchins. PLoS One, 7, e53118.

Preliminaries

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

Import schlegel data file (schlegel.csv)

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

Set contrasts to sum from afex

set_sum_contrasts()

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

Convert pH and indiv to factors

schlegel$ph <- factor(schlegel$ph)
schlegel$indiv <- factor(schlegel$indiv)

Check assumptions with boxplots

boxplot(avspeedtot~ph, data=schlegel)

Fit model with default aov type I SS

schlegel.aov <- aov(avspeedtot~ph*indiv, data=schlegel)

Check residuals - look OK

plot(schlegel.aov)

Examine results

options(digits = 4)
summary(schlegel.aov)

             Df Sum Sq Mean Sq F value  Pr(>F)    
ph            2    469     235    15.9 1.9e-07 ***
indiv        18  11166     620    42.2 < 2e-16 ***
ph:indiv     36    688      19     1.3    0.12    
Residuals   509   7488      15                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get correct F-ratio for ph

234.7/19.1

[1] 12.29

21-pf(234.7/19.1, 2, 36, lower.tail = TRUE, log.p = FALSE)

[1] 20

Fit models with lm to get type III SS using car

options(digits=10)
schlegel.lm1 <- lm(avspeedtot~ph*indiv, data=schlegel)
Anova(schlegel.lm1, type='III')

Anova Table (Type III tests)

Response: avspeedtot
                Sum Sq  Df     F value     Pr(>F)    
(Intercept) 265761.549   1 18064.53988 < 2.22e-16 ***
ph             470.394   2    15.98697 1.8463e-07 ***
indiv        11170.259  18    42.18184 < 2.22e-16 ***
ph:indiv       687.936  36     1.29891    0.11855    
Residuals     7488.296 509                           
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get correct F-ratio and P value for ph

235.197/19.109

[1] 12.30817939

1-pf(235.197/19.109, 2, 36, lower.tail = TRUE, log.p = FALSE)

[1] 8.449437684e-05

Using VCA package to get anova var comps (with CIs that can be -ve)

schlegel.vca <- anovaMM(avspeedtot~ph+(indiv)+(ph*indiv), schlegel)
schlegel.vca

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

    [Fixed Effects]

      int     ph7.6     ph7.8     ph8.1 
22.932291 -2.153381 -1.604031  0.000000 


    [Variance Components]

  Name     DF        SS          MS         VC        %Total    SD       CV[%]    
1 total    51.006764                        35.337903 100       5.944569 27.349571
2 indiv    18        11166.04084 620.335602 20.183217 57.114925 4.492574 20.669279
3 ph:indiv 36        687.936231  19.10934   0.442906  1.253344  0.665512 3.061863 
4 error    509       7488.296372 14.711781  14.711781 41.631731 3.835594 17.646669

Mean: 21.735512 (N = 566) 

Experimental Design: unbalanced  |  Method: ANOVA

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

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

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

    [Fixed Effects]

    int   ph7.6   ph7.8   ph8.1 
22.9323 -2.1534 -1.6040  0.0000 


    [Variance Components]

  Name     DF      SS         MS       VC      %Total  SD     CV[%]   Var(VC)
1 total    51.0068                     35.3379 100     5.9446 27.3496        
2 indiv    18      11166.0408 620.3356 20.1832 57.1149 4.4926 20.6693 48.2322
3 ph:indiv 36      687.9362   19.1093  0.4429  1.2533  0.6655 3.0619  0.2144 
4 error    509     7488.2964  14.7118  14.7118 41.6317 3.8356 17.6467 0.8504 

Mean: 21.7355 (N = 566) 

Experimental Design: unbalanced  |  Method: ANOVA


> VC:
-----
         Estimate  CI LCL  CI UCL One-Sided LCL One-Sided UCL
total     35.3379 24.8192 54.3450       26.2456       50.6234
indiv     20.1832  6.5714 33.7951        8.7598       31.6066
ph:indiv   0.4429 -0.4647  1.3505       -0.3188        1.2046
error     14.7118 13.0593 16.7008       13.3103       16.3614

> SD:
-----
         Estimate  CI LCL CI UCL One-Sided LCL One-Sided UCL
total      5.9446  4.9819 7.3719        5.1230        7.1150
indiv      4.4926  2.5635 5.8134        2.9597        5.6220
ph:indiv   0.6655 -0.6817 1.1621       -0.5646        1.0975
error      3.8356  3.6138 4.0867        3.6483        4.0449

> CV[%]:
--------
         Estimate  CI LCL  CI UCL One-Sided LCL One-Sided UCL
total     27.3496 22.9205 33.9164       23.5699       32.7345
indiv     20.6693 11.7939 26.7459       13.6169       25.8654
ph:indiv   3.0619 -3.1362  5.3466       -2.5975        5.0495
error     17.6467 16.6261 18.8018       16.7851       18.6097


95% Confidence Level  
SAS PROC MIXED method used for computing CIs 

Fit random “slope” and intercept (i.e. full) model using lme4 and REML to compare models with LR test

schlegel.lmer1 <- lmer(avspeedtot~ph + (1|indiv) + (1|ph:indiv), REML=TRUE, schlegel)
summary(schlegel.lmer1)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: avspeedtot ~ ph + (1 | indiv) + (1 | ph:indiv)
   Data: schlegel

REML criterion at convergence: 3206.4

Scaled residuals: 
       Min         1Q     Median         3Q        Max 
-3.2039636 -0.6556270 -0.0297379  0.6068584  3.8945501 

Random effects:
 Groups   Name        Variance   Std.Dev. 
 ph:indiv (Intercept)  0.4451441 0.6671912
 indiv    (Intercept) 20.5554324 4.5338099
 Residual             14.7103546 3.8354080
Number of obs: 566, groups:  ph:indiv, 57; indiv, 19

Fixed effects:
              Estimate Std. Error         df  t value   Pr(>|t|)    
(Intercept) 21.6797931  1.0562577 17.9831797 20.52510 6.2745e-14 ***
ph1         -0.9009115  0.2602038 36.2020559 -3.46233   0.001392 ** 
ph2         -0.3515580  0.2599264 36.0546682 -1.35253   0.184629    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
    (Intr) ph1   
ph1  0.000       
ph2  0.000 -0.501

AICc(schlegel.lmer1)

[1] 3218.569462

Fit random intercept-only mixed model using lme4 and REML to compare models with LR test

schlegel.lmer2 <- lmer(avspeedtot~ph + (1|indiv), REML=TRUE, schlegel)
summary(schlegel.lmer2)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: avspeedtot ~ ph + (1 | indiv)
   Data: schlegel

REML criterion at convergence: 3207.7

Scaled residuals: 
       Min         1Q     Median         3Q        Max 
-3.2111046 -0.6685941 -0.0356255  0.6208264  3.9690798 

Random effects:
 Groups   Name        Variance Std.Dev.
 indiv    (Intercept) 20.69177 4.548821
 Residual             15.00250 3.873306
Number of obs: 566, groups:  indiv, 19

Fixed effects:
               Estimate  Std. Error          df  t value   Pr(>|t|)    
(Intercept)  21.6802670   1.0562017  17.9828315 20.52664 6.2688e-14 ***
ph1          -0.9013241   0.2304615 544.9859217 -3.91095 0.00010354 ***
ph2          -0.3517865   0.2301462 544.9836170 -1.52853 0.12696002    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
    (Intr) ph1   
ph1  0.000       
ph2  0.000 -0.501

AICc(schlegel.lmer2)

[1] 3217.784777

Compare model fit to test random interaction fit

anova(schlegel.lmer1, schlegel.lmer2)

refitting model(s) with ML (instead of REML)

Data: schlegel
Models:
schlegel.lmer2: avspeedtot ~ ph + (1 | indiv)
schlegel.lmer1: avspeedtot ~ ph + (1 | indiv) + (1 | ph:indiv)
               npar       AIC       BIC     logLik  deviance   Chisq Df Pr(>Chisq)
schlegel.lmer2    5 3217.1076 3238.8006 -1603.5538 3207.1076                      
schlegel.lmer1    6 3218.2849 3244.3164 -1603.1424 3206.2849 0.82275  1    0.36438

ranova(schlegel.lmer2, reduce.terms=TRUE)

ANOVA-like table for random-effects: Single term deletions

Model:
avspeedtot ~ ph + (1 | indiv)
            npar     logLik       AIC       LRT Df Pr(>Chisq)    
<none>         5 -1603.8388 3217.6776                            
(1 | indiv)    4 -1803.4214 3614.8427 399.16508  1 < 2.22e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get CIs

schlegel.ci1 <- confint.merMod(schlegel.lmer1)

Computing profile confidence intervals ...

schlegel.vc1 <- (schlegel.ci1)^2
print(schlegel.vc1)

                     2.5 %          97.5 %
.sig01        0.0000000000   1.49493825337
.sig02       10.6793375185  40.29003392955
.sigma       13.0418238075  16.67571388336
(Intercept) 382.5301594921 566.48711669828
ph1           1.9881649319   0.15338949529
ph2           0.7399354733   0.02469286792

Refit simpler model and change reference group to 8.1 for meaningful contrasts

schlegel$ph <- relevel(schlegel$ph, ref="8.1")
schlegel.lmer3 <- lmer(avspeedtot~ph + (1|indiv), REML=TRUE, schlegel)
summary(schlegel.lmer3, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: avspeedtot ~ ph + (1 | indiv)
   Data: schlegel

REML criterion at convergence: 3207.7

Scaled residuals: 
       Min         1Q     Median         3Q        Max 
-3.2111046 -0.6685941 -0.0356255  0.6208264  3.9690798 

Random effects:
 Groups   Name        Variance Std.Dev.
 indiv    (Intercept) 20.69177 4.548821
 Residual             15.00250 3.873306
Number of obs: 566, groups:  indiv, 19

Fixed effects:
               Estimate  Std. Error          df  t value   Pr(>|t|)    
(Intercept)  21.6802670   1.0562018  17.9999924 20.52663  6.147e-14 ***
ph1           1.2531106   0.2301462 545.0007765  5.44485  7.854e-08 ***
ph2          -0.9013241   0.2304616 545.0030792 -3.91095 0.00010354 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
    (Intr) ph1   
ph1  0.000       
ph2  0.000 -0.501

ranova(schlegel.lmer3, reduce.terms=TRUE)

ANOVA-like table for random-effects: Single term deletions

Model:
avspeedtot ~ ph + (1 | indiv)
            npar     logLik       AIC       LRT Df Pr(>Chisq)    
<none>         5 -1603.8388 3217.6776                            
(1 | indiv)    4 -1803.4214 3614.8427 399.16508  1 < 2.22e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

F tests of fixed effects (Type III SS) - matches t-test from summary command

anova(schlegel.lmer3, ddf="Kenward-Roger")

Type III Analysis of Variance Table with Kenward-Roger's method
      Sum Sq   Mean Sq NumDF     DenDF  F value     Pr(>F)    
ph 472.87444 236.43722     2 545.00154 15.75986 2.2192e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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

schlegel.ci3 <- confint.merMod(schlegel.lmer3)

Computing profile confidence intervals ...

schlegel.vc3 <- (schlegel.ci3)^2
print(schlegel.vc3)

                     2.5 %         97.5 %
.sig01       10.7911510977  40.3923816785
.sigma       13.3075963683  16.8686962959
(Intercept) 382.5485943306 566.5011081919
ph1           0.6433013719   2.9041293937
ph2           1.8305473550   0.2021863353