QK Box 12.5

Parkinson et al. (2002) examined the differences in the abundance of birds between different habitats on a floodplain. Habitat was a between-subjects factor, and there were three groups: permanent billabongs, temporary billabongs, and non-wetland habitat. Four sites (subjects) were chosen within each habitat group. The authors were also interested in whether the pattern between habitats was consistent over months (austral summer: November, December, January, February) and time of the day (morning and evening chorus), so each site was sampled morning and evening on one day each month; month and time of day were the crossed within-subjects factors. The relative abundance of birds was the response variable.

Australian Pelican, Pelecanus conspicillatus. Mick Keough, CC BY-SA 4.0.

Australasian Darter, Anhinga novaehollandiae. Mick Keough, CC BY-SA 4.0.

Link to paper and data

Parkinson, A., Mac Nally, R. & Quinn, G. P. (2002). Differential macrohabitat use by birds on the unregulated Ovens River floodplain of southeastern Australia. River Research and Applications, 18, 495-506.

Preliminaries

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

Import parkinson data file

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

Set contrasts using afex

set_sum_contrasts()

Check residuals - uneven spread related to mean

park1.aov <- aov(abund~habitat*ampm*month, parkinson)
plot(park1.aov)

Unequal spread

Original authors used 4th root transform but we will use log10 here

By now, you should be able to run your own version of this analysis with 4th root transform.

parkinson$labund <- log10(parkinson$abund)
park2.aov <- aov(labund~habitat*ampm*month, parkinson)
plot(park2.aov)

Log10 improves spread in residual plot

Run as split-plot

park3.aov <- aov(labund~habitat*month*ampm+Error(bong/(month*ampm)), parkinson)
summary(park3.aov)

Error: bong
          Df Sum Sq Mean Sq F value   Pr(>F)    
habitat    2  3.941  1.9705    36.6 4.76e-05 ***
Residuals  9  0.485  0.0538                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: bong:month
              Df Sum Sq Mean Sq F value Pr(>F)   
month          3 0.0597 0.01990   0.717 0.5502   
habitat:month  6 0.7201 0.12002   4.326 0.0035 **
Residuals     27 0.7490 0.02774                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: bong:ampm
             Df  Sum Sq Mean Sq F value Pr(>F)  
ampm          1 0.24433 0.24433   9.273 0.0139 *
habitat:ampm  2 0.01888 0.00944   0.358 0.7084  
Residuals     9 0.23714 0.02635                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: bong:month:ampm
                   Df Sum Sq Mean Sq F value Pr(>F)  
month:ampm          3 0.1755 0.05849   2.534 0.0779 .
habitat:month:ampm  6 0.1432 0.02386   1.034 0.4253  
Residuals          27 0.6231 0.02308                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Run with ezANOVA to get GG and HF adjustments

ezpark1 <- ezANOVA(data=parkinson, dv=labund, wid=bong, within=month*ampm, between=habitat, type=3, detailed=TRUE)

Warning: Converting "bong" to factor for ANOVA.Warning: Converting "month" to factor for ANOVA.Warning: Converting "ampm" to factor for ANOVA.Warning: Converting "habitat" to factor for ANOVA.

print(ezpark1)

$ANOVA
              Effect DFn DFd          SSn       SSd            F            p p<.05         ges
1        (Intercept)   1   9 185.82587994 0.4845966 3451.1861758 6.044855e-13     * 0.988857661
2            habitat   2   9   3.94109825 0.4845966   36.5973345 4.756581e-05     * 0.653044220
3              month   3  27   0.05970429 0.7489888    0.7174187 5.502357e-01       0.027723404
5               ampm   1   9   0.24432786 0.2371379    9.2728783 1.390178e-02     * 0.104494288
4      habitat:month   6  27   0.72010127 0.7489888    4.3264409 3.495319e-03     * 0.255902542
6       habitat:ampm   2   9   0.01888359 0.2371379    0.3583408 7.083704e-01       0.008937925
7         month:ampm   3  27   0.17547497 0.6231424    2.5343723 7.791341e-02       0.077324211
8 habitat:month:ampm   6  27   0.14317702 0.6231424    1.0339477 4.252937e-01       0.064002814

$`Mauchly's Test for Sphericity`
              Effect         W          p p<.05
3              month 0.7890302 0.87321554      
4      habitat:month 0.7890302 0.87321554      
7         month:ampm 0.2455222 0.05655546      
8 habitat:month:ampm 0.2455222 0.05655546      

$`Sphericity Corrections`
              Effect       GGe       p[GG] p[GG]<.05       HFe       p[HF] p[HF]<.05
3              month 0.8526384 0.530946777           1.2200581 0.550235681          
4      habitat:month 0.8526384 0.006128013         * 1.2200581 0.003495319         *
7         month:ampm 0.5380760 0.121213461           0.6382668 0.110170806          
8 habitat:month:ampm 0.5380760 0.411016520           0.6382668 0.415745122          

Draw bar graph

Pool times as time didn’t interact with anything

park_stats <- summarySE(parkinson,measurevar='labund', groupvars= c('habitat','month'))
ggplot(park_stats, aes(x=month, y=labund, fill=habitat))+
  geom_bar(stat="identity", position=position_dodge())+
  geom_errorbar(aes(ymin=labund-se, ymax=labund+se), position=position_dodge(0.9), width=0.3, color="darkblue")+
  scale_fill_uchicago()

Run as mixed effects with lme4

park1.lmer <- lmer(labund~habitat+month+ampm+habitat*month+habitat*ampm+month*ampm+habitat*month*ampm+(1|bong)+(1|bong:month)+(1|bong:ampm), parkinson)
summary(park1.lmer, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: labund ~ habitat + month + ampm + habitat * month + habitat *  
    ampm + month * ampm + habitat * month * ampm + (1 | bong) +      (1 | bong:month) + (1 | bong:ampm)
   Data: parkinson

REML criterion at convergence: 30.9

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.32462 -0.45772 -0.02583  0.43892  1.99277 

Random effects:
 Groups     Name        Variance  Std.Dev.
 bong:month (Intercept) 0.0023301 0.04827 
 bong:ampm  (Intercept) 0.0008171 0.02858 
 bong       (Intercept) 0.0028547 0.05343 
 Residual               0.0230796 0.15192 
Number of obs: 96, groups:  bong:month, 48; bong:ampm, 24; bong, 12

Fixed effects:
                       Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)            1.391289   0.023683  9.000000  58.746 6.05e-13 ***
habitat1               0.078617   0.033493  9.000000   2.347 0.043500 *  
habitat2               0.199321   0.033493  9.000000   5.951 0.000215 ***
month1                -0.032860   0.029443 27.000000  -1.116 0.274232    
month2                -0.012602   0.029443 27.000000  -0.428 0.672037    
month3                 0.033117   0.029443 27.000000   1.125 0.270581    
ampm1                  0.050449   0.016567  9.000000   3.045 0.013901 *  
habitat1:month1       -0.095749   0.041638 27.000000  -2.300 0.029437 *  
habitat2:month1        0.048929   0.041638 27.000000   1.175 0.250218    
habitat1:month2        0.181415   0.041638 27.000000   4.357 0.000171 ***
habitat2:month2       -0.141554   0.041638 27.000000  -3.400 0.002112 ** 
habitat1:month3        0.009678   0.041638 27.000000   0.232 0.817961    
habitat2:month3        0.057023   0.041638 27.000000   1.369 0.182134    
habitat1:ampm1        -0.015329   0.023429  9.000000  -0.654 0.529298    
habitat2:ampm1        -0.003236   0.023429  9.000000  -0.138 0.893183    
month1:ampm1          -0.008753   0.026856 27.000000  -0.326 0.746987    
month2:ampm1           0.067314   0.026856 27.000000   2.506 0.018519 *  
month3:ampm1          -0.051518   0.026856 27.000000  -1.918 0.065702 .  
habitat1:month1:ampm1  0.024475   0.037980 27.000000   0.644 0.524745    
habitat2:month1:ampm1  0.031960   0.037980 27.000000   0.841 0.407460    
habitat1:month2:ampm1 -0.042139   0.037980 27.000000  -1.110 0.276997    
habitat2:month2:ampm1 -0.024581   0.037980 27.000000  -0.647 0.522957    
habitat1:month3:ampm1 -0.035500   0.037980 27.000000  -0.935 0.358222    
habitat2:month3:ampm1  0.012819   0.037980 27.000000   0.338 0.738333    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation matrix not shown by default, as p = 24 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

anova(park1.lmer, ddf="Kenward-Roger")

Type III Analysis of Variance Table with Kenward-Roger's method
                    Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
habitat            1.68926 0.84463     2     9 36.5963 4.757e-05 ***
month              0.04967 0.01656     3    27  0.7174  0.550230    
ampm               0.21402 0.21402     1     9  9.2731  0.013901 *  
habitat:month      0.59912 0.09985     6    27  4.3265  0.003495 ** 
habitat:ampm       0.01654 0.00827     2     9  0.3583  0.708365    
month:ampm         0.17547 0.05849     3    27  2.5343  0.077916 .  
habitat:month:ampm 0.14318 0.02386     6    27  1.0339  0.425301    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get variance components

park1.ci <- confint.merMod(park1.lmer, oldNames=FALSE)

Computing profile confidence intervals ...

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

                                 2.5 %       97.5 %
sd_(Intercept)|bong:month 0.0000000000 9.205968e-03
sd_(Intercept)|bong:ampm  0.0000000000 7.721817e-03
sd_(Intercept)|bong       0.0000000000 1.000865e-02
sigma                     0.0112715453 2.597718e-02
(Intercept)               1.8161607435 2.059021e+00
habitat1                  0.0002856626 1.969313e-02
habitat2                  0.0189354222 6.814016e-02
month1                    0.0070893803 3.414726e-04
month2                    0.0040883970 1.500544e-03
month3                    0.0003320442 7.132707e-03
ampm1                     0.0003968621 6.557028e-03
habitat1:month1           0.0283428708 5.357182e-04
habitat2:month1           0.0005605165 1.477017e-02
habitat1:month2           0.0118397560 6.452538e-02
habitat2:month2           0.0458635686 4.754101e-03
habitat1:month3           0.0039597034 6.770270e-03
habitat2:month3           0.0002427766 1.680304e-02
habitat1:ampm1            0.0034223958 7.751915e-04
habitat2:ampm1            0.0021537385 1.594814e-03
month1:ampm1              0.0030904695 1.450516e-03
month2:ampm1              0.0004192180 1.303078e-02
month3:ampm1              0.0096741040 2.189676e-05
habitat1:month1:ampm1     0.0017443317 8.229206e-03
habitat2:month1:ampm1     0.0011751096 9.643298e-03
habitat1:month2:ampm1     0.0117459953 5.808691e-04
habitat2:month2:ampm1     0.0082484957 1.735467e-03
habitat1:month3:ampm1     0.0103510775 9.449429e-04
habitat2:month3:ampm1     0.0028537899 6.250372e-03

Evaluate whether to keep site within habitat by ampm or by month interactions (use ML)

park1a.lmer <- lmer(labund~habitat+month+ampm+habitat*month+habitat*ampm+month*ampm+habitat*month*ampm+(1|bong)+(1|bong:month)+(1|bong:ampm), REML=FALSE,parkinson)
park2.lmer <- lmer(labund~habitat+month+ampm+habitat*month+habitat*ampm+month*ampm+habitat*month*ampm+(1|bong)+(1|bong:month), REML=FALSE, parkinson)
park3.lmer <- lmer(labund~habitat+month+ampm+habitat*month+habitat*ampm+month*ampm+habitat*month*ampm+(1|bong)+(1|bong:ampm), REML=FALSE, parkinson)
anova(park1a.lmer, park2.lmer)

Data: parkinson
Models:
park2.lmer: labund ~ habitat + month + ampm + habitat * month + habitat * ampm + month * ampm + habitat * month * ampm + (1 | bong) + (1 | bong:month)
park1a.lmer: labund ~ habitat + month + ampm + habitat * month + habitat * ampm + month * ampm + habitat * month * ampm + (1 | bong) + (1 | bong:month) + (1 | bong:ampm)
            npar     AIC    BIC logLik deviance  Chisq Df Pr(>Chisq)
park2.lmer    27 -44.529 24.708 49.265  -98.529                     
park1a.lmer   28 -42.610 29.192 49.305  -98.610 0.0807  1     0.7763

anova(park1a.lmer, park3.lmer)

Data: parkinson
Models:
park3.lmer: labund ~ habitat + month + ampm + habitat * month + habitat * ampm + month * ampm + habitat * month * ampm + (1 | bong) + (1 | bong:ampm)
park1a.lmer: labund ~ habitat + month + ampm + habitat * month + habitat * ampm + month * ampm + habitat * month * ampm + (1 | bong) + (1 | bong:month) + (1 | bong:ampm)
            npar     AIC    BIC logLik deviance  Chisq Df Pr(>Chisq)
park3.lmer    27 -44.306 24.931 49.153  -98.306                     
park1a.lmer   28 -42.610 29.192 49.305  -98.610 0.3041  1     0.5813

AICc(park1a.lmer, park2.lmer, park3.lmer)

Neither contribute to model fit, refit model without them

park4.lmer <- lmer(labund~habitat+month+ampm+habitat*month+habitat*ampm+month*ampm+habitat*month*ampm+(1|bong), REML=FALSE, parkinson)
anova(park1a.lmer, park4.lmer)

Data: parkinson
Models:
park4.lmer: labund ~ habitat + month + ampm + habitat * month + habitat * ampm + month * ampm + habitat * month * ampm + (1 | bong)
park1a.lmer: labund ~ habitat + month + ampm + habitat * month + habitat * ampm + month * ampm + habitat * month * ampm + (1 | bong) + (1 | bong:month) + (1 | bong:ampm)
            npar     AIC    BIC logLik deviance  Chisq Df Pr(>Chisq)
park4.lmer    26 -46.299 20.374 49.150  -98.299                     
park1a.lmer   28 -42.610 29.192 49.305  -98.610 0.3109  2      0.856

AICc(park1a.lmer, park2.lmer, park3.lmer, park4.lmer)

Refit simpler model with REML

park4a.lmer <- lmer(labund~habitat+month+ampm+habitat*month+habitat*ampm+month*ampm+habitat*month*ampm+(1|bong), parkinson)
summary(park4a.lmer, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: labund ~ habitat + month + ampm + habitat * month + habitat *  
    ampm + month * ampm + habitat * month * ampm + (1 | bong)
   Data: parkinson

REML criterion at convergence: 31.1

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.45047 -0.51655 -0.04362  0.53175  1.97897 

Random effects:
 Groups   Name        Variance Std.Dev.
 bong     (Intercept) 0.003538 0.05948 
 Residual             0.025544 0.15982 
Number of obs: 96, groups:  bong, 12

Fixed effects:
                       Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)            1.391289   0.023683  9.000000  58.747 6.04e-13 ***
habitat1               0.078617   0.033493  9.000000   2.347 0.043497 *  
habitat2               0.199321   0.033493  9.000000   5.951 0.000215 ***
month1                -0.032860   0.028253 63.000000  -1.163 0.249199    
month2                -0.012602   0.028253 63.000000  -0.446 0.657104    
month3                 0.033117   0.028253 63.000000   1.172 0.245557    
ampm1                  0.050449   0.016312 63.000000   3.093 0.002954 ** 
habitat1:month1       -0.095749   0.039956 63.000000  -2.396 0.019539 *  
habitat2:month1        0.048929   0.039956 63.000000   1.225 0.225301    
habitat1:month2        0.181415   0.039956 63.000000   4.540 2.59e-05 ***
habitat2:month2       -0.141554   0.039956 63.000000  -3.543 0.000752 ***
habitat1:month3        0.009678   0.039956 63.000000   0.242 0.809405    
habitat2:month3        0.057023   0.039956 63.000000   1.427 0.158481    
habitat1:ampm1        -0.015329   0.023069 63.000000  -0.664 0.508801    
habitat2:ampm1        -0.003236   0.023069 63.000000  -0.140 0.888885    
month1:ampm1          -0.008753   0.028253 63.000000  -0.310 0.757725    
month2:ampm1           0.067314   0.028253 63.000000   2.383 0.020225 *  
month3:ampm1          -0.051518   0.028253 63.000000  -1.823 0.072981 .  
habitat1:month1:ampm1  0.024475   0.039956 63.000000   0.613 0.542385    
habitat2:month1:ampm1  0.031960   0.039956 63.000000   0.800 0.426787    
habitat1:month2:ampm1 -0.042139   0.039956 63.000000  -1.055 0.295623    
habitat2:month2:ampm1 -0.024581   0.039956 63.000000  -0.615 0.540635    
habitat1:month3:ampm1 -0.035500   0.039956 63.000000  -0.888 0.377662    
habitat2:month3:ampm1  0.012819   0.039956 63.000000   0.321 0.749402    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation matrix not shown by default, as p = 24 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

anova(park4a.lmer, ddf="Kenward-Roger")

Type III Analysis of Variance Table with Kenward-Roger's method
                    Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
habitat            1.86968 0.93484     2     9 36.5973 4.757e-05 ***
month              0.05970 0.01990     3    63  0.7791 0.5100272    
ampm               0.24433 0.24433     1    63  9.5650 0.0029545 ** 
habitat:month      0.72010 0.12002     6    63  4.6984 0.0005205 ***
habitat:ampm       0.01888 0.00944     2    63  0.3696 0.6924788    
month:ampm         0.17547 0.05849     3    63  2.2898 0.0868987 .  
habitat:month:ampm 0.14318 0.02386     6    63  0.9342 0.4769777    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1