Q & K Box 10.3

Bleeker et al. (2017) examined variability in reproductive characteristics of a goby (Neogonius melanostomus) along the Rhine River in Europe. They had five sampling sites along the Rhine and imagine that the five sites represented a random sample from a population of possible sites. At each site they had between seven and 22 fish; the response variable was gonad mass and the covariate was body length (centered as a body length of zero is nonsensical). Our interest is how the relationship between gonad mass and body length varied across the random effect of site.

Eric Engbretson, U.S. Fish and Wildlife Service, Public domain, via Wikimedia Commons

The paper is here

Bleeker, K., de Jong, K., van Kessel, N., Hinde, C. A. & Nagelkerke, L. A. J. (2017). Evidence for ontogenetically and morphologically distinct alternative reproductive tactics in the invasive Round Goby Neogobius melanostomus. PLoS One, 12.

Preliminaries

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

Import bleeker data file (bleeker1.csv)

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

Set contrasts from afex

set_sum_contrasts()

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

Focus on models with covariate

Check linearity and then variances

panel.smoother <- function(x, y) {
  panel.xyplot(x, y) 
  panel.lmline(x, y)
}
xyplot(gm~tl|river,main="Scatterplots by River", ylab="gonad weight", xlab="total length", panel=panel.smoother, bleeker1)

# check boxplots - variance difference for gm but not related to mean
boxplot(gm~river, bleeker1)

boxplot(tl~river, bleeker1)

Linearity seems OK and variance difference for gonad mass but not related to mean

Fit random coefficients model

bleeker.lmer1 <- lmer(gm~tl +(tl|river), REML=TRUE, bleeker1)

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

summary(bleeker.lmer1)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: gm ~ tl + (tl | river)
   Data: bleeker1

REML criterion at convergence: 8.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6723 -0.6896 -0.2498  0.5697  2.6659 

Random effects:
 Groups   Name        Variance  Std.Dev.  Corr 
 river    (Intercept) 1.267e-05 0.0035598      
          tl          2.036e-07 0.0004513 -1.00
 Residual             5.687e-02 0.2384758      
Number of obs: 69, groups:  river, 5

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)    
(Intercept) -0.002204   0.108021 64.778750  -0.020 0.983788    
tl           0.042661   0.010310 60.038508   4.138 0.000111 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
   (Intr)
tl -0.964
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

plot(bleeker.lmer1)

isSingular(bleeker.lmer1)

[1] TRUE

Remove correlation between random effects

bleeker.lmer2 <- lmer(gm~tl +(tl||river), REML=TRUE, bleeker1)

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

summary(bleeker.lmer2)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: gm ~ tl + (tl || river)
   Data: bleeker1

REML criterion at convergence: 8.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6723 -0.6894 -0.2498  0.5698  2.6659 

Random effects:
 Groups   Name        Variance Std.Dev.
 river    (Intercept) 0.00000  0.0000  
 river.1  tl          0.00000  0.0000  
 Residual             0.05687  0.2385  
Number of obs: 69, groups:  river, 5

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)    
(Intercept) -0.002245   0.108005 67.000000  -0.021    0.983    
tl           0.042666   0.010307 67.000000   4.139 9.95e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
   (Intr)
tl -0.964
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

plot(bleeker.lmer2)

isSingular(bleeker.lmer2)

[1] TRUE

Removing correlation makes little difference to model’s fit

Singular warning so try with centered covariate

bleeker1$ctl <- scale(bleeker1$tl, center=TRUE, scale=FALSE)
bleeker.lmer3 <- lmer(gm~ctl +(ctl||river), REML=TRUE, bleeker1)

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

summary(bleeker.lmer3)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: gm ~ ctl + (ctl || river)
   Data: bleeker1

REML criterion at convergence: 8.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6723 -0.6894 -0.2498  0.5698  2.6659 

Random effects:
 Groups   Name        Variance Std.Dev.
 river    (Intercept) 0.00000  0.0000  
 river.1  ctl         0.00000  0.0000  
 Residual             0.05687  0.2385  
Number of obs: 69, groups:  river, 5

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)  0.42874    0.02871 67.00000  14.934  < 2e-16 ***
ctl          0.04267    0.01031 67.00000   4.139 9.95e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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

isSingular(bleeker.lmer3)

[1] TRUE

Correlation between random effects reduced but still warning about singularity

Proceeded to compare random coefficients and random intercepts model with nlme. Also generate AICc for each model

bleeker.lme1 <- lme(fixed=gm~ctl, random=~1|river/ctl, method="REML", bleeker1)
bleeker.lme2 <- lme(fixed=gm~ctl, random=~1|river, method="REML", bleeker1)
anova(bleeker.lme1, bleeker.lme2)

#AICc
AICc(bleeker.lme1)

[1] 18.73185

AICc(bleeker.lme2)

[1] 17.19429

Difference in fit of the two models is not great, so we will ignore the random coefficients (river x length)

Fit random intercept model with lmer

bleeker.lmer4 <- lmer(gm~ctl +(1|river), REML=TRUE, bleeker1)

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

plot(bleeker.lmer4)

Residuals look fine

summary(bleeker.lmer4, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: gm ~ ctl + (1 | river)
   Data: bleeker1

REML criterion at convergence: 8.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6723 -0.6894 -0.2498  0.5698  2.6659 

Random effects:
 Groups   Name        Variance Std.Dev.
 river    (Intercept) 0.00000  0.0000  
 Residual             0.05687  0.2385  
Number of obs: 69, groups:  river, 5

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)  0.42874    0.03074  2.77327  13.946 0.001190 ** 
ctl          0.04267    0.01170 36.24997   3.648 0.000825 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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

anova(bleeker.lmer4, ddf="Kenward-Roger")

Type III Analysis of Variance Table with Kenward-Roger's method
     Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
ctl 0.75688 0.75688     1 36.25  13.308 0.0008248 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get profile confidence intervals from lmer fit

bleeker.ci <- confint.merMod(bleeker.lmer4)

Computing profile confidence intervals ...

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

                   2.5 %      97.5 %
.sig01      0.0000000000 0.008510938
.sigma      0.0402567780 0.078627247
(Intercept) 0.1382150218 0.236989230
ctl         0.0005035901 0.003950525