QK Box 12.4

Skrip et al (2016) studied the effects of dietary antioxidants and exercise on the allocation of nutrients to eggs in zebra finches (Taeniopygia guttata). We will focus on one aspect of their study, how the order of egg laying affected the mass of the eggs. Mating pairs of finches were allocated randomly to two diet groups (supplemented with high antioxidant food, no supplement) and two exercise groups (additional flight exercise, no additional exercise). The finches laid between one and nine eggs but few laid more than five and pairs that laid only one egg don’t provide enough information for fitting a regression between order and egg mass. So we will analyze egg mass as the response variable from pairs that laid between two and five eggs. The between-subjects component was a two factor (diet and exercise) crossed design with pairs as subjects, and egg order (within each pair) was the within-subjects factor.

A pair of Australian zebra finches. PotMart186, , via Wikimedia Commons.

Link to paper

Preliminaries

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

Import skrip data file

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

set contrasts using afex

set_sum_contrasts()

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

Diagnostics

check individual regressions

xyplot(wmass~order|pair, groups=group, type=c("p","r"), auto.key=T, skrip)

list_skrip <- lmList(wmass~order|pair, skrip, na.action=na.omit)
summary(list_skrip)

Call:
  Model: wmass ~ order | pair 
   Data: skrip 

Coefficients:
   (Intercept) 
   Estimate Std. Error  t value     Pr(>|t|)
a 1.1207500 0.08654913 12.94929 1.308318e-16
b 0.9963171 0.08530375 11.67964 4.502662e-15
c 0.9457000 0.08654913 10.92674 4.047996e-14
d 1.2784300 0.07411624 17.24899 3.364372e-21
e 1.0951667 0.10794573 10.14553 4.270621e-13
f 1.1296500 0.08654913 13.05212 9.910952e-17
g 1.1619600 0.07411624 15.67753 1.265078e-19
h 1.2429500 0.08654913 14.36121 3.233641e-18
i 1.0043300 0.07411624 13.55074 2.626501e-17
j 1.2924600 0.07411624 17.43828 2.209927e-21
k 1.1215000 0.07411624 15.13164 4.738873e-19
l 1.0617000 0.07411624 14.32480 3.546833e-18
m 0.9854500 0.08654913 11.38602 1.051081e-14
n 0.9360900 0.07411624 12.63002 3.124596e-16
o 1.0614000 0.08654913 12.26356 8.622791e-16
p 0.9439400 0.07411624 12.73594 2.337511e-16
q 1.0535300 0.07411624 14.21456 4.696690e-18
r 1.0016300 0.07411624 13.51431 2.891131e-17
   order 
      Estimate Std. Error    t value     Pr(>|t|)
a -0.013050000 0.03271249 -0.3989301 6.918743e-01
b  0.008848571 0.02388983  0.3703907 7.128685e-01
c  0.010560000 0.03160327  0.3341426 7.398594e-01
d -0.013730000 0.02234689 -0.6144032 5.421131e-01
e  0.055050000 0.04996916  1.1016795 2.765911e-01
f  0.004920000 0.03160327  0.1556801 8.769972e-01
g  0.014480000 0.02234689  0.6479649 5.203747e-01
h -0.013040000 0.03160327 -0.4126155 6.818920e-01
i  0.030930000 0.02234689  1.3840853 1.733130e-01
j -0.017040000 0.02234689 -0.7625223 4.498170e-01
k  0.040600000 0.02234689  1.8168078 7.606196e-02
l  0.009220000 0.02234689  0.4125854 6.819139e-01
m -0.016890000 0.03160327 -0.5344383 5.957287e-01
n  0.089510000 0.02234689  4.0054794 2.349747e-04
o  0.005820000 0.03160327  0.1841581 8.547355e-01
p  0.123080000 0.02234689  5.5077020 1.775475e-06
q  0.094630000 0.02234689  4.2345941 1.148112e-04
r  0.049090000 0.02234689  2.1967264 3.335257e-02

Residual standard error: 0.07066707 on 44 degrees of freedom

some evidence of different slopes

First model with random slopes and intercepts

skrip1.lmer <- lmer(wmass~diet*exercise*order+(order|pair), REML=FALSE, skrip, na.action=na.omit)

Second just random intercepts

skrip2.lmer <- lmer(wmass~diet*exercise*order+(1|pair), REML=FALSE, skrip, na.action=na.omit)

Test random slopes term

rand(skrip1.lmer)

ANOVA-like table for random-effects: Single term deletions

Model:
wmass ~ diet + exercise + order + (order | pair) + diet:exercise + diet:order + exercise:order + diet:exercise:order
                        npar logLik     AIC    LRT Df Pr(>Chisq)
<none>                    12 81.191 -138.38                     
order in (order | pair)   10 79.831 -139.66 2.7195  2     0.2567

anova(skrip2.lmer,skrip1.lmer)

Data: skrip
Models:
skrip2.lmer: wmass ~ diet * exercise * order + (1 | pair)
skrip1.lmer: wmass ~ diet * exercise * order + (order | pair)
            npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)
skrip2.lmer   10 -139.66 -115.84 79.831  -159.66                     
skrip1.lmer   12 -138.38 -109.80 81.191  -162.38 2.7195  2     0.2567

AICc(skrip2.lmer,skrip1.lmer)

No strong evidence for random slopes so use random intercepts.

Refit model with REML

skrip3.lmer <- lmer(wmass~diet*exercise*order+(1|pair), REML=TRUE, skrip, na.action=na.omit)

Check residual plot - OK

plot(skrip3.lmer)

summary(skrip3.lmer, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: wmass ~ diet * exercise * order + (1 | pair)
   Data: skrip

REML criterion at convergence: -103.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.1386 -0.4493  0.0028  0.4177  2.5362 

Random effects:
 Groups   Name        Variance Std.Dev.
 pair     (Intercept) 0.008465 0.09200 
 Residual             0.005614 0.07493 
Number of obs: 80, groups:  pair, 18

Fixed effects:
                       Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)            1.070290   0.029803 31.437992  35.912  < 2e-16 ***
diet1                  0.062936   0.029803 31.437992   2.112   0.0428 *  
exercise1              0.010650   0.029803 31.437992   0.357   0.7232    
order                  0.032613   0.006371 58.803754   5.119 3.55e-06 ***
diet1:exercise1        0.018304   0.029803 31.437992   0.614   0.5435    
diet1:order           -0.029235   0.006371 58.803754  -4.589 2.39e-05 ***
exercise1:order        0.014902   0.006371 58.803754   2.339   0.0228 *  
diet1:exercise1:order -0.012184   0.006371 58.803754  -1.913   0.0607 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) diet1  exrcs1 order  dt1:x1 dt1:rd exrc1:
diet1       -0.115                                          
exercise1    0.115 -0.145                                   
order       -0.600  0.028 -0.028                            
diet1:xrcs1 -0.145  0.115 -0.115  0.084                     
diet1:order  0.028 -0.600  0.084 -0.009 -0.028              
exercs1:rdr -0.028  0.084 -0.600  0.009  0.028 -0.137       
dt1:xrcs1:r  0.084 -0.028  0.028 -0.137 -0.600  0.009 -0.009

anova(skrip3.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)    
diet                0.025037 0.025037     1 31.438  4.4594   0.04275 *  
exercise            0.000717 0.000717     1 31.438  0.1277   0.72323    
order               0.147133 0.147133     1 58.804 26.2061 3.548e-06 ***
diet:exercise       0.002118 0.002118     1 31.438  0.3772   0.54353    
diet:order          0.118230 0.118230     1 58.804 21.0582 2.389e-05 ***
exercise:order      0.030718 0.030718     1 58.804  5.4712   0.02275 *  
diet:exercise:order 0.020536 0.020536     1 58.804  3.6578   0.06068 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get var comps with CIs

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

Computing profile confidence intervals ...

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

                             2.5 %       97.5 %
sd_(Intercept)|pair   2.952456e-03 1.455114e-02
sigma                 3.768687e-03 7.636959e-03
(Intercept)           1.031731e+00 1.265179e+00
diet1                 7.192387e-05 1.381911e-02
exercise1             1.924024e-03 4.250565e-03
order                 4.169228e-04 2.016701e-03
diet1:exercise1       1.318915e-03 5.293660e-03
diet1:order           1.732100e-03 2.929702e-04
exercise1:order       6.808788e-06 7.343093e-04
diet1:exercise1:order 5.905052e-04 4.097410e-08

Draw line plot

skrip_stats <- summarySE(skrip,measurevar='wmass', groupvars= c('group','order'), na.rm=TRUE)
skrip_stats

ggplot(skrip_stats, aes(x=order, y=wmass, fill=group, colour=group))+
  geom_line()+
  theme_qk()+
  scale_colour_uchicago()