Q & K Box 13.8

Ryeland et al (2017) studied the roosting behaviour of several species of shorebirds. They recorded the proportion of time (number of minutes as a proportion of total minutes in a video bout) individuals of various species spent in the backrest position while roosting. They used a binomial model with a logit link for proportions with four fixed predictors recorded for each video bout: ambient temperature, wind speed, size of group focal bird was in, and distance focal bird was from the observer. We will analyze the data for a single species, the sharp-tailed sandpiper (Calidris acuminata). This would be a standard binomial GLM except that more than one bird was sometimes recorded in each bout so bout was included as a random effect since birds closer together may be correlated in their behaviour. The resulting model is a binomial GLMM.

Sharp-tailed sandpiper. patrickkavanagh, , via Wikimedia Commons

The paper is here

Ryeland, J., Weston, M. A., Symonds, M. R. E. & Overgaard, J. (2017). Bill size mediates behavioural thermoregulation in birds. Functional Ecology, 31, 885-93.

Preliminaries

First, load the required packages (car, performance, MuMIn, lme4, glmmTMB, lmtest, ggplot2)

Import ryeland data file (ryeland.csv)

ryeland <- read.csv("../data/ryeland.csv")
head(ryeland)

Create proportion of time facing back

ryeland$prophb<-ryeland$timehb/ryeland$filmp

Create success and fail columns

ryeland$success<-as.integer(ryeland$timehb)
ryeland$fail<-as.integer(ryeland$filmp-ryeland$timehb)

Create response variable

ryeland.prop<-cbind(ryeland$success,ryeland$fail)

Preliminary checks

Scatterplots

plot(prophb~wind, data=ryeland)

plot(prophb~dist, data=ryeland)

plot(prophb~groupsize, data=ryeland)

plot(prophb~temp, data=ryeland)

Boxplots

boxplot(ryeland$dist)

boxplot(log10(ryeland$dist))

boxplot(ryeland$groupsize)

boxplot(ryeland$temp)

boxplot(ryeland$wind)

boxplot(log10(ryeland$wind))

Check collinearity

vif(lm(ryeland$prophb ~ dist+groupsize+temp+wind,data=ryeland))

     dist groupsize      temp      wind 
     1.23      2.03      1.85      1.16 

cor(ryeland[,c('dist','groupsize','temp','wind')])

             dist groupsize    temp   wind
dist       1.0000    -0.203  0.0162 -0.274
groupsize -0.2030     1.000  0.6525 -0.157
temp       0.0162     0.652  1.0000 -0.047
wind      -0.2738    -0.157 -0.0470  1.000

Fit GLMM with binomial distribution and test predictor effects

ryeland1.glmer<-glmer(ryeland.prop~dist+groupsize+temp+wind+(1|vbout),data=ryeland,family=binomial)
summary(ryeland1.glmer)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: ryeland.prop ~ dist + groupsize + temp + wind + (1 | vbout)
   Data: ryeland

     AIC      BIC   logLik deviance df.resid 
     224      234     -106      212       36 

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.529 -0.697 -0.221  0.408  3.338 

Random effects:
 Groups Name        Variance Std.Dev.
 vbout  (Intercept) 1.43     1.19    
Number of obs: 42, groups:  vbout, 21

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  4.46783    1.23192    3.63  0.00029 ***
dist         0.00364    0.00770    0.47  0.63666    
groupsize    0.00316    0.01311    0.24  0.80938    
temp        -0.24929    0.05531   -4.51  6.6e-06 ***
wind        -0.03106    0.09391   -0.33  0.74084    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
          (Intr) dist   gropsz temp  
dist      -0.375                     
groupsize -0.155  0.453              
temp      -0.551 -0.362 -0.597       
wind      -0.395  0.387  0.345 -0.200

AIC(ryeland1.glmer)

[1] 224

AICc(ryeland1.glmer)

[1] 226

anova(ryeland1.glmer, type="lrtest")

Warning: additional arguments ignored: ‘type’

Analysis of Variance Table
          npar Sum Sq Mean Sq F value
dist         1   0.08    0.08    0.08
groupsize    1   8.33    8.33    8.33
temp         1  22.00   22.00   22.00
wind         1   0.11    0.11    0.11

Test wind

ryeland2.glmer<-glmer(ryeland.prop~dist+groupsize+temp+(1|vbout),data=ryeland,family=binomial)
lrtest(ryeland1.glmer, ryeland2.glmer)

Likelihood ratio test

Model 1: ryeland.prop ~ dist + groupsize + temp + wind + (1 | vbout)
Model 2: ryeland.prop ~ dist + groupsize + temp + (1 | vbout)
  #Df LogLik Df Chisq Pr(>Chisq)
1   6   -106                    
2   5   -106 -1  0.11       0.74

Test temp

ryeland3.glmer<-glmer(ryeland.prop~dist+groupsize+wind+(1|vbout),data=ryeland,family=binomial)
lrtest(ryeland1.glmer, ryeland3.glmer)

Likelihood ratio test

Model 1: ryeland.prop ~ dist + groupsize + temp + wind + (1 | vbout)
Model 2: ryeland.prop ~ dist + groupsize + wind + (1 | vbout)
  #Df LogLik Df Chisq Pr(>Chisq)    
1   6   -106                        
2   5   -115 -1  17.4      3e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Test groupsize

ryeland4.glmer<-glmer(ryeland.prop~dist+temp+wind+(1|vbout),data=ryeland,family=binomial)
lrtest(ryeland1.glmer, ryeland4.glmer)

Likelihood ratio test

Model 1: ryeland.prop ~ dist + groupsize + temp + wind + (1 | vbout)
Model 2: ryeland.prop ~ dist + temp + wind + (1 | vbout)
  #Df LogLik Df Chisq Pr(>Chisq)
1   6   -106                    
2   5   -106 -1  0.06       0.81

Test dist

ryeland5.glmer<-glmer(ryeland.prop~groupsize+wind+temp+wind+(1|vbout),data=ryeland,family=binomial)
lrtest(ryeland1.glmer, ryeland5.glmer)

Likelihood ratio test

Model 1: ryeland.prop ~ dist + groupsize + temp + wind + (1 | vbout)
Model 2: ryeland.prop ~ groupsize + wind + temp + wind + (1 | vbout)
  #Df LogLik Df Chisq Pr(>Chisq)
1   6   -106                    
2   5   -106 -1  0.22       0.64

Check residuals

plot(ryeland1.glmer, resid(., type = "deviance") ~ fitted(.))

residuals(ryeland1.glmer, type="deviance")

      1       2       3       4       5       6       7       8       9      10      11      12      13      14      15      16      17      18      19      20      21      22 
-3.7634  2.8380  0.8204 -0.6598 -0.5124 -0.5124  3.8535 -3.4811 -1.7306  2.0263  2.7415 -1.5641  0.2951  0.1486 -0.9874 -0.9874  0.0823 -0.4955  0.4147 -1.1464  2.3104 -1.8960 
     23      24      25      26      27      28      29      30      31      32      33      34      35      36      37      38      39      40      41      42 
-0.5897 -0.5897 -1.1729  1.3381  0.0043  0.0043 -0.2598 -0.2598 -0.3127 -0.3127 -0.2564 -0.2564 -1.8043  2.4378 -0.3901  0.5396  0.4258 -1.1398 -0.6372 -0.6372 

Check overdispersion

(code from Zuur et al 2013)

presid1 <- resid(ryeland1.glmer, type="pearson")
ssize1 <- nrow(ryeland)
params1 <- length(fixef(ryeland1.glmer)+1)
disp1 <- sum(presid1^2)/(ssize1-params1)
disp1

[1] 2.23

Add Observation Level Random Effect into model and retest the predictors

ryeland6.glmer<-glmer(ryeland.prop~dist+groupsize+temp+wind+(1|vbout)+(1|olre),data=ryeland,family=binomial)

Warning: Model failed to converge with max|grad| = 0.00479232 (tol = 0.002, component 1)

summary(ryeland6.glmer)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: ryeland.prop ~ dist + groupsize + temp + wind + (1 | vbout) +      (1 | olre)
   Data: ryeland

     AIC      BIC   logLik deviance df.resid 
     176      188      -81      162       35 

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.129 -0.365 -0.110  0.146  0.806 

Random effects:
 Groups Name        Variance Std.Dev.
 olre   (Intercept) 2.282    1.511   
 vbout  (Intercept) 0.653    0.808   
Number of obs: 42, groups:  olre, 42; vbout, 21

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  5.03538    1.39144    3.62   0.0003 ***
dist         0.00488    0.00859    0.57   0.5695    
groupsize    0.00664    0.01469    0.45   0.6511    
temp        -0.28931    0.06361   -4.55  5.4e-06 ***
wind        -0.02193    0.10482   -0.21   0.8343    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
          (Intr) dist   gropsz temp  
dist      -0.355                     
groupsize -0.125  0.448              
temp      -0.567 -0.368 -0.607       
wind      -0.387  0.391  0.344 -0.203
optimizer (Nelder_Mead) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.00479232 (tol = 0.002, component 1)

AIC(ryeland6.glmer)

[1] 176

AICc(ryeland6.glmer)

[1] 179

Test wind

ryeland7.glmer<-glmer(ryeland.prop~dist+groupsize+temp+(1|vbout)+(1|olre),data=ryeland,family=binomial)
lrtest(ryeland6.glmer, ryeland7.glmer)

Likelihood ratio test

Model 1: ryeland.prop ~ dist + groupsize + temp + wind + (1 | vbout) + 
    (1 | olre)
Model 2: ryeland.prop ~ dist + groupsize + temp + (1 | vbout) + (1 | olre)
  #Df LogLik Df Chisq Pr(>Chisq)
1   7    -81                    
2   6    -81 -1  0.04       0.83

Test temp

ryeland8.glmer<-glmer(ryeland.prop~dist+groupsize+wind+(1|vbout)+(1|olre),data=ryeland,family=binomial)
lrtest(ryeland6.glmer, ryeland8.glmer)

Likelihood ratio test

Model 1: ryeland.prop ~ dist + groupsize + temp + wind + (1 | vbout) + 
    (1 | olre)
Model 2: ryeland.prop ~ dist + groupsize + wind + (1 | vbout) + (1 | olre)
  #Df LogLik Df Chisq Pr(>Chisq)    
1   7  -81.0                        
2   6  -90.1 -1  18.2      2e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Test groupsize

ryeland9.glmer<-glmer(ryeland.prop~dist+temp+wind+(1|vbout)+(1|olre),data=ryeland,family=binomial)
lrtest(ryeland6.glmer, ryeland9.glmer)

Likelihood ratio test

Model 1: ryeland.prop ~ dist + groupsize + temp + wind + (1 | vbout) + 
    (1 | olre)
Model 2: ryeland.prop ~ dist + temp + wind + (1 | vbout) + (1 | olre)
  #Df LogLik Df Chisq Pr(>Chisq)
1   7  -81.0                    
2   6  -81.1 -1   0.2       0.65

Test dist

ryeland10.glmer<-glmer(ryeland.prop~groupsize+wind+temp+wind+(1|vbout)+(1|olre),data=ryeland,family=binomial)
lrtest(ryeland6.glmer, ryeland10.glmer)

Likelihood ratio test

Model 1: ryeland.prop ~ dist + groupsize + temp + wind + (1 | vbout) + 
    (1 | olre)
Model 2: ryeland.prop ~ groupsize + wind + temp + wind + (1 | vbout) + 
    (1 | olre)
  #Df LogLik Df Chisq Pr(>Chisq)
1   7  -81.0                    
2   6  -81.1 -1  0.32       0.57

Check residuals

residuals(ryeland6.glmer, type="deviance")

        1         2         3         4         5         6         7         8         9        10        11        12        13        14        15        16        17        18 
-1.544739  0.517286  0.119681 -0.032263 -0.468487 -0.468487  0.532565 -0.413694 -0.252082  0.424475  0.888607  0.099123  0.205676  0.194426 -0.863413 -0.863413 -0.092399 -0.205755 
       19        20        21        22        23        24        25        26        27        28        29        30        31        32        33        34        35        36 
 0.129297 -0.810856  0.465602  0.081393 -0.519161 -0.519161 -0.556220  0.670993  0.008321  0.008321 -0.193177 -0.193177 -0.229764 -0.229764 -0.180583 -0.180583  0.021693  0.501009 
       37        38        39        40        41        42 
 0.000936  0.219131  0.147782 -0.798193 -0.502843 -0.502843 

plot(ryeland6.glmer, resid(., type = "deviance") ~ fitted(.))

Fit model with only temp to match paper

ryeland11.glmer<-glmer(ryeland.prop~temp+(1|vbout)+(1|olre),data=ryeland,family=binomial)
summary(ryeland11.glmer)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: ryeland.prop ~ temp + (1 | vbout) + (1 | olre)
   Data: ryeland

     AIC      BIC   logLik deviance df.resid 
   170.6    177.5    -81.3    162.6       38 

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.128 -0.350 -0.108  0.163  0.812 

Random effects:
 Groups Name        Variance Std.Dev.
 olre   (Intercept) 2.264    1.505   
 vbout  (Intercept) 0.781    0.884   
Number of obs: 42, groups:  olre, 42; vbout, 21

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   5.1380     1.2625    4.07  4.7e-05 ***
temp         -0.2701     0.0508   -5.31  1.1e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
     (Intr)
temp -0.953

AIC(ryeland11.glmer)

[1] 171

AICc(ryeland11.glmer)

[1] 172

Plot temp to match paper Need to add fit logistic curve?

ggplot(ryeland, aes(temp, prophb)) + geom_point() +
  geom_smooth(method = "loess", span = 0.5, se = FALSE) 

Rerun original GLMM with glmmTMB

ryeland.glmmb <- glmmTMB(ryeland.prop~groupsize+wind+temp+wind+(1|vbout), family=binomial, data=ryeland)
summary(ryeland.glmmb)

 Family: binomial  ( logit )
Formula:          ryeland.prop ~ groupsize + wind + temp + wind + (1 | vbout)
Data: ryeland

     AIC      BIC   logLik deviance df.resid 
     222      231     -106      212       37 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 vbout  (Intercept) 1.45     1.2     
Number of obs: 42, groups:  vbout, 21

Conditional model:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  4.69865    1.14625    4.10  4.1e-05 ***
groupsize    0.00034    0.01174    0.03     0.98    
wind        -0.04793    0.08722   -0.55     0.58    
temp        -0.24056    0.05164   -4.66  3.2e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Now fit betabinomial to compare to OLRE

ryeland.glmmbb <- glmmTMB(ryeland.prop~dist+groupsize+temp+wind+(1|vbout),data=ryeland,family=betabinomial(link="logit"))
summary(ryeland.glmmbb)

 Family: betabinomial  ( logit )
Formula:          ryeland.prop ~ dist + groupsize + temp + wind + (1 | vbout)
Data: ryeland

     AIC      BIC   logLik deviance df.resid 
   178.7    190.8    -82.3    164.7       35 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 vbout  (Intercept) 0.124    0.353   
Number of obs: 42, groups:  vbout, 21

Dispersion parameter for betabinomial family ():  2.6 

Conditional model:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  3.39744    1.02782    3.31  0.00095 ***
dist         0.00215    0.00601    0.36  0.72045    
groupsize    0.00685    0.00984    0.70  0.48647    
temp        -0.19311    0.04846   -3.98  6.8e-05 ***
wind        -0.01632    0.06696   -0.24  0.80738    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

AICc(ryeland.glmmbb)

[1] 182

Test temperature

ryeland1.glmmbb <- glmmTMB(ryeland.prop~dist+groupsize+wind+(1|vbout),data=ryeland,family=betabinomial(link="logit"))
lrtest(ryeland.glmmbb, ryeland1.glmmbb)

Likelihood ratio test

Model 1: ryeland.prop ~ dist + groupsize + temp + wind + (1 | vbout)
Model 2: ryeland.prop ~ dist + groupsize + wind + (1 | vbout)
  #Df LogLik Df Chisq Pr(>Chisq)    
1   7  -82.3                        
2   6  -91.9 -1  19.1    1.2e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1