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.94   1.971    36.6 4.8e-05 ***
Residuals  9   0.48   0.054                    
---
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.060  0.0199    0.72 0.5502   
habitat:month  6  0.720  0.1200    4.33 0.0035 **
Residuals     27  0.749  0.0277                  
---
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.2443  0.2443    9.27  0.014 *
habitat:ampm  2 0.0189  0.0094    0.36  0.708  
Residuals     9 0.2371  0.0263                 
---
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.175  0.0585    2.53  0.078 .
habitat:month:ampm  6  0.143  0.0239    1.03  0.425  
Residuals          27  0.623  0.0231                 
---
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.825880 0.48460 3451.18618 6.0449e-13     * 0.9888577
2            habitat   2   9   3.941098 0.48460   36.59733 4.7566e-05     * 0.6530442
3              month   3  27   0.059704 0.74899    0.71742 5.5024e-01       0.0277234
5               ampm   1   9   0.244328 0.23714    9.27288 1.3902e-02     * 0.1044943
4      habitat:month   6  27   0.720101 0.74899    4.32644 3.4953e-03     * 0.2559025
6       habitat:ampm   2   9   0.018884 0.23714    0.35834 7.0837e-01       0.0089379
7         month:ampm   3  27   0.175475 0.62314    2.53437 7.7913e-02       0.0773242
8 habitat:month:ampm   6  27   0.143177 0.62314    1.03395 4.2529e-01       0.0640028

$`Mauchly's Test for Sphericity`
              Effect       W        p p<.05
3              month 0.78903 0.873216      
4      habitat:month 0.78903 0.873216      
7         month:ampm 0.24552 0.056555      
8 habitat:month:ampm 0.24552 0.056555      

$`Sphericity Corrections`
              Effect     GGe    p[GG] p[GG]<.05     HFe     p[HF] p[HF]<.05
3              month 0.85264 0.530947           1.22006 0.5502357          
4      habitat:month 0.85264 0.006128         * 1.22006 0.0034953         *
7         month:ampm 0.53808 0.121213           0.63827 0.1101708          
8 habitat:month:ampm 0.53808 0.411017           0.63827 0.4157451          

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.3246 -0.4577 -0.0258  0.4389  1.9928 

Random effects:
 Groups     Name        Variance Std.Dev.
 bong:month (Intercept) 0.002330 0.0483  
 bong:ampm  (Intercept) 0.000817 0.0286  
 bong       (Intercept) 0.002855 0.0534  
 Residual               0.023080 0.1519  
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.39129    0.02368  9.00000   58.75    6e-13 ***
habitat1               0.07862    0.03349  9.00000    2.35  0.04350 *  
habitat2               0.19932    0.03349  9.00000    5.95  0.00022 ***
month1                -0.03286    0.02944 27.00000   -1.12  0.27423    
month2                -0.01260    0.02944 27.00000   -0.43  0.67204    
month3                 0.03312    0.02944 27.00000    1.12  0.27058    
ampm1                  0.05045    0.01657  9.00000    3.05  0.01390 *  
habitat1:month1       -0.09575    0.04164 27.00000   -2.30  0.02944 *  
habitat2:month1        0.04893    0.04164 27.00000    1.18  0.25022    
habitat1:month2        0.18141    0.04164 27.00000    4.36  0.00017 ***
habitat2:month2       -0.14155    0.04164 27.00000   -3.40  0.00211 ** 
habitat1:month3        0.00968    0.04164 27.00000    0.23  0.81796    
habitat2:month3        0.05702    0.04164 27.00000    1.37  0.18213    
habitat1:ampm1        -0.01533    0.02343  9.00000   -0.65  0.52930    
habitat2:ampm1        -0.00324    0.02343  9.00000   -0.14  0.89318    
month1:ampm1          -0.00875    0.02686 27.00000   -0.33  0.74699    
month2:ampm1           0.06731    0.02686 27.00000    2.51  0.01852 *  
month3:ampm1          -0.05152    0.02686 27.00000   -1.92  0.06570 .  
habitat1:month1:ampm1  0.02447    0.03798 27.00000    0.64  0.52475    
habitat2:month1:ampm1  0.03196    0.03798 27.00000    0.84  0.40746    
habitat1:month2:ampm1 -0.04214    0.03798 27.00000   -1.11  0.27700    
habitat2:month2:ampm1 -0.02458    0.03798 27.00000   -0.65  0.52296    
habitat1:month3:ampm1 -0.03550    0.03798 27.00000   -0.93  0.35822    
habitat2:month3:ampm1  0.01282    0.03798 27.00000    0.34  0.73833    
---
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.689   0.845     2     9   36.60 4.8e-05 ***
month               0.050   0.017     3    27    0.72  0.5502    
ampm                0.214   0.214     1     9    9.27  0.0139 *  
habitat:month       0.599   0.100     6    27    4.33  0.0035 ** 
habitat:ampm        0.017   0.008     2     9    0.36  0.7084    
month:ampm          0.175   0.058     3    27    2.53  0.0779 .  
habitat:month:ampm  0.143   0.024     6    27    1.03  0.4253    
---
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.00000000 9.2066e-03
sd_(Intercept)|bong:ampm  0.00000000 7.7209e-03
sd_(Intercept)|bong       0.00000000 1.0009e-02
sigma                     0.01127155 2.5977e-02
(Intercept)               1.81616074 2.0590e+00
habitat1                  0.00028566 1.9693e-02
habitat2                  0.01893542 6.8140e-02
month1                    0.00708938 3.4147e-04
month2                    0.00408840 1.5005e-03
month3                    0.00033204 7.1327e-03
ampm1                     0.00039686 6.5570e-03
habitat1:month1           0.02834287 5.3572e-04
habitat2:month1           0.00056052 1.4770e-02
habitat1:month2           0.01183976 6.4525e-02
habitat2:month2           0.04586357 4.7541e-03
habitat1:month3           0.00395970 6.7703e-03
habitat2:month3           0.00024278 1.6803e-02
habitat1:ampm1            0.00342240 7.7519e-04
habitat2:ampm1            0.00215374 1.5948e-03
month1:ampm1              0.00309047 1.4505e-03
month2:ampm1              0.00041922 1.3031e-02
month3:ampm1              0.00967410 2.1897e-05
habitat1:month1:ampm1     0.00174433 8.2292e-03
habitat2:month1:ampm1     0.00117511 9.6433e-03
habitat1:month2:ampm1     0.01174600 5.8087e-04
habitat2:month2:ampm1     0.00824850 1.7355e-03
habitat1:month3:ampm1     0.01035108 9.4494e-04
habitat2:month3:ampm1     0.00285379 6.2504e-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.5 24.7   49.3    -98.5                    
park1a.lmer   28 -42.6 29.2   49.3    -98.6  0.08  1       0.78

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.3 24.9   49.2    -98.3                    
park1a.lmer   28 -42.6 29.2   49.3    -98.6   0.3  1       0.58

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.3 20.4   49.1    -98.3                    
park1a.lmer   28 -42.6 29.2   49.3    -98.6  0.31  2       0.86

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.4505 -0.5166 -0.0436  0.5317  1.9790 

Random effects:
 Groups   Name        Variance Std.Dev.
 bong     (Intercept) 0.00354  0.0595  
 Residual             0.02554  0.1598  
Number of obs: 96, groups:  bong, 12

Fixed effects:
                      Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)            1.39129    0.02368  9.00000   58.75  6.0e-13 ***
habitat1               0.07862    0.03349  9.00000    2.35  0.04350 *  
habitat2               0.19932    0.03349  9.00000    5.95  0.00022 ***
month1                -0.03286    0.02825 63.00000   -1.16  0.24920    
month2                -0.01260    0.02825 63.00000   -0.45  0.65710    
month3                 0.03312    0.02825 63.00000    1.17  0.24556    
ampm1                  0.05045    0.01631 63.00000    3.09  0.00295 ** 
habitat1:month1       -0.09575    0.03996 63.00000   -2.40  0.01954 *  
habitat2:month1        0.04893    0.03996 63.00000    1.22  0.22530    
habitat1:month2        0.18141    0.03996 63.00000    4.54  2.6e-05 ***
habitat2:month2       -0.14155    0.03996 63.00000   -3.54  0.00075 ***
habitat1:month3        0.00968    0.03996 63.00000    0.24  0.80941    
habitat2:month3        0.05702    0.03996 63.00000    1.43  0.15848    
habitat1:ampm1        -0.01533    0.02307 63.00000   -0.66  0.50880    
habitat2:ampm1        -0.00324    0.02307 63.00000   -0.14  0.88889    
month1:ampm1          -0.00875    0.02825 63.00000   -0.31  0.75772    
month2:ampm1           0.06731    0.02825 63.00000    2.38  0.02023 *  
month3:ampm1          -0.05152    0.02825 63.00000   -1.82  0.07298 .  
habitat1:month1:ampm1  0.02447    0.03996 63.00000    0.61  0.54239    
habitat2:month1:ampm1  0.03196    0.03996 63.00000    0.80  0.42679    
habitat1:month2:ampm1 -0.04214    0.03996 63.00000   -1.05  0.29562    
habitat2:month2:ampm1 -0.02458    0.03996 63.00000   -0.62  0.54064    
habitat1:month3:ampm1 -0.03550    0.03996 63.00000   -0.89  0.37766    
habitat2:month3:ampm1  0.01282    0.03996 63.00000    0.32  0.74940    
---
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.870   0.935     2     9   36.60 4.8e-05 ***
month               0.060   0.020     3    63    0.78 0.51004    
ampm                0.244   0.244     1    63    9.56 0.00295 ** 
habitat:month       0.720   0.120     6    63    4.70 0.00052 ***
habitat:ampm        0.019   0.009     2    63    0.37 0.69248    
month:ampm          0.175   0.058     3    63    2.29 0.08690 .  
habitat:month:ampm  0.143   0.024     6    63    0.93 0.47698    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1