QK Box 12.1

Steiger et al. (2008) studied the Coolidge effect, the decline in males’ interest in mating with the same female compared to novel females, using the burying beetle Nicrophorus vespilloides. Eighteen male beetles were presented with the same female beetle four times, and then a novel female on the fifth occasion. This was a repeated measures design as the same individual males were repeatedly presented with females. There was no evidence that physical exhaustion affected time to mating as a separate control group of males were presented with novel, unmated females five times in succession, and there was no change in time to mating. The within-subjects factor was the order of presented females, and while this could have been treated as a continuous covariate, we treated it as a fixed factor with five groups. The response variable recorded on each occasion was time to mating.

Adult burying beetle. Francisco Welter-Schultes. Wikimedia Commons Creative Commons CC0 1.0 Universal Public Domain Dedication

Francisco Welter-Schultes, CC0, via Wikimedia Commons

Steiger, S., Franz, R., Eggert, A. K. & Muller, J. K. (2008). The Coolidge effect, individual recognition and selection for distinctive cuticular signatures in a burying beetle. Proceedings of the Royal Society B, 275, 1831-8.

Link to the paper: doi: 10.1098/rspb.2008.0375 and data

Preliminaries

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

Import steiger data file (steiger.csv)

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

Set contrasts from afex. Make individual a factor Rearrange matingnumber order to override default alphabetical order

set_sum_contrasts()

steiger$individual <- factor(steiger$individual)
steiger$matingnumber <- factor(steiger$matingnumber, levels=c("first","second","third","fourth","fifth"))

Quick boxplot to look at variances for different mating numbers

boxplot(time~matingnumber, data=steiger)

Fit OLS model with default aov SS

steiger.aov <- aov(time~matingnumber+individual, data=steiger)

Check residuals - some evidence for interaction

plot(steiger.aov)

Do interaction plot

interaction.plot(steiger$matingnumber, steiger$individual, steiger$time)

The rank order of treatments consistent but the sizes of the differences vary greatly between individuals. Quickly try boxplot with log-transformed data to see if any improvement

boxplot(log10(time)~matingnumber, data=steiger)

Transform to logs due to variance heterogeneity and to minimise interaction

steiger$ltime <- log10(steiger$time)

Recheck diagnostics - much better with less interaction

boxplot(ltime~matingnumber, data=steiger)

interaction.plot(steiger$matingnumber, steiger$individual, steiger$ltime)

Fit OLS model with default aov SS

steiger.aov1 <- aov(ltime~matingnumber+individual, data=steiger)

check residuals - look OK

plot(steiger.aov1)

Examine results

Once we’re happy with the model that we’re fitting, we can look at the results.

tidy(steiger.aov1)

emmeans(steiger.aov1, "matingnumber")

 matingnumber emmean    SE df lower.CL upper.CL
 first          1.47 0.136 68     1.20     1.74
 second         1.83 0.136 68     1.56     2.10
 third          2.01 0.136 68     1.74     2.28
 fourth         2.14 0.136 68     1.87     2.41
 fifth          1.35 0.136 68     1.08     1.62

Results are averaged over the levels of: individual 
Confidence level used: 0.95 

Get Greenhouse Geiser adjusted results

ezsteiger1 <- ezANOVA(data=steiger, dv=ltime, wid=individual, within=matingnumber, type=3, detailed=TRUE)
print(ezsteiger1)

$ANOVA
        Effect DFn DFd    SSn   SSd      F        p p<.05   ges
1  (Intercept)   1  17 278.55  6.72 704.49 2.81e-15     * 0.905
2 matingnumber   4  68   8.36 22.50   6.32 2.20e-04     * 0.222

$`Mauchly's Test for Sphericity`
        Effect     W     p p<.05
2 matingnumber 0.431 0.166      

$`Sphericity Corrections`
        Effect   GGe   p[GG] p[GG]<.05   HFe    p[HF] p[HF]<.05
2 matingnumber 0.678 0.00155         * 0.819 0.000657         *

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

steiger.vca <- anovaMM(ltime~matingnumber+(individual), steiger)
steiger.vca

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

    [Fixed Effects]

               int  matingnumberfifth  matingnumberfirst matingnumberfourth matingnumbersecond  matingnumberthird 
             2.009             -0.658             -0.544              0.131             -0.180              0.000 


    [Variance Components]

  Name       DF        SS        MS       VC       %Total    SD       CV[%]    
1 total      84.524647                    0.34383  100       0.58637  33.330258
2 individual 17        6.721791  0.395399 0.012892 3.749614  0.113544 6.454044 
3 error      68        22.503779 0.330938 0.330938 96.250386 0.575272 32.69941 

Mean: 1.76 (N = 90) 

Experimental Design: balanced  |  Method: ANOVA

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

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

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

    [Fixed Effects]

               int  matingnumberfifth  matingnumberfirst matingnumberfourth matingnumbersecond  matingnumberthird 
             2.009             -0.657             -0.543              0.131             -0.180              0.000 


    [Variance Components]

  Name       DF      SS      MS     VC     %Total  SD     CV[%]   Var(VC)
1 total      84.5246                0.3438 100     0.5864 33.3303 0.0028 
2 individual 17      6.7218  0.3954 0.0129 3.7496  0.1135 6.454   9e-04  
3 error      68      22.5038 0.3309 0.3309 96.2504 0.5753 32.6994 0.0032 

Mean: 1.76 (N = 90) 

Experimental Design: balanced  |  Method: ANOVA


> VC:
-----
           Estimate  CI LCL CI UCL One-Sided LCL One-Sided UCL
total        0.3438  0.2598 0.4765        0.2716        0.4517
individual   0.0129 -0.0447 0.0705       -0.0355        0.0613
error        0.3309  0.2428 0.4779        0.2550        0.4499

> SD:
-----
           Estimate CI LCL CI UCL One-Sided LCL One-Sided UCL
total         0.586  0.510  0.690         0.521         0.672
individual    0.114 -0.211  0.266        -0.188         0.247
error         0.575  0.493  0.691         0.505         0.671

> CV[%]:
--------
           Estimate CI LCL CI UCL One-Sided LCL One-Sided UCL
total         33.33     29   39.2          29.6          38.2
individual     6.45    -12   15.1         -10.7          14.1
error         32.70     28   39.3          28.7          38.1


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

Fit random intercept model using lme4 and REML

steiger.lmer1 <- lmer(ltime~matingnumber + (1|individual), REML=TRUE, steiger)
summary(steiger.lmer1, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: ltime ~ matingnumber + (1 | individual)
   Data: steiger

REML criterion at convergence: 168

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.921 -0.815  0.205  0.698  1.927 

Random effects:
 Groups     Name        Variance Std.Dev.
 individual (Intercept) 0.0129   0.114   
 Residual               0.3309   0.575   
Number of obs: 90, groups:  individual, 18

Fixed effects:
              Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)     1.7593     0.0663 17.0000   26.54  2.8e-15 ***
matingnumber1  -0.2934     0.1213 68.0000   -2.42   0.0182 *  
matingnumber2   0.0700     0.1213 68.0000    0.58   0.5657    
matingnumber3   0.2501     0.1213 68.0000    2.06   0.0430 *  
matingnumber4   0.3809     0.1213 68.0000    3.14   0.0025 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) mtngn1 mtngn2 mtngn3
matingnmbr1  0.000                     
matingnmbr2  0.000 -0.250              
matingnmbr3  0.000 -0.250 -0.250       
matingnmbr4  0.000 -0.250 -0.250 -0.250

anova(steiger.lmer1, type=3, ddf="Kenward-Roger")

Type III Analysis of Variance Table with Kenward-Roger's method
             Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)    
matingnumber   8.36    2.09     4    68    6.32 0.00022 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

emmeans(steiger.lmer1, ~matingnumber)

 matingnumber emmean    SE   df lower.CL upper.CL
 first          1.47 0.138 84.5     1.19     1.74
 second         1.83 0.138 84.5     1.55     2.10
 third          2.01 0.138 84.5     1.74     2.28
 fourth         2.14 0.138 84.5     1.86     2.42
 fifth          1.35 0.138 84.5     1.08     1.63

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

Note different CIs compared to OLS model fitting due to K-R adjustment

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

steiger.ci1 <- confint.merMod(steiger.lmer1, oldNames=FALSE)

Computing profile confidence intervals ...

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

                             2.5 %  97.5 %
sd_(Intercept)|individual 0.000000 0.09445
sigma                     0.229286 0.43705
(Intercept)               2.643805 3.58182
matingnumber1             0.278325 0.00352
matingnumber2             0.026936 0.09249
matingnumber3             0.000254 0.23443
matingnumber4             0.021535 0.37820

Bar graph for time main effect

steiger_sum <- summarySE(steiger, measurevar= 'ltime', groupvars= 'matingnumber')
ggplot(steiger_sum, aes(x=matingnumber, y=ltime))+
  geom_bar(stat="identity", position="dodge", fill="lightblue")+
  geom_errorbar(aes(ymin=ltime-se, ymax=ltime+se), width=0.3, color="darkblue")