QK Box 12.3

Garcia et al. (2015) studied the effect of the appetite-regulating hormone leptin on appetite and mating preferences in the spadefoot toad Spea bombifrons. Eighteen female toads collected from the wild were allocated to a treatment group (n = 9) which received a subcutaneous injection of leptin once per day for six days, and a control group (n = 9), which received saline injections with the same frequency. One hour after the day 6 injections, each toad was presented with approximately 50 crickets. The response variable was the cumulative number of attacks by each toad over three-minute intervals for 15 minutes. Treatment (leptin versus control was the fixed between-subject factor and toads were the subjects. The within-subjects fixed factor was time with five groups representing 3, 6, 9, 12 and 15 minutes after the introduction of crickets.

Plains spadefoot toad. USFWS Mountain-Prairie, Public domain, via Wikimedia Commons

Garcia, N. W., Pfennig, K. S. & Burmeister, S. S. (2015). Leptin manipulation reduces appetite and causes a switch in mating preference in the Plains Spadefoot Toad (Spea bombifrons). PLoS One, 10, e0125981.

Link to the paper

Preliminaries

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

Import garcia data file

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

set contrasts using afex

set_sum_contrasts()

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

set toad as factor and create timefac as categorical version of time

garcia$toad <- factor(garcia$toad)
garcia$timefac <- factor(garcia$time, ordered=TRUE)

Diagnostics

Check residuals - uneven spread related to mean

garcia1.aov <- aov(cumattack~treatment*time, garcia)
plot(garcia1.aov)

Check boxplots - unequal spread related to mean

boxplot(cumattack~treatment*time,data=garcia)

Check homogeneity of within-group variances - very different variances

garcia_stats <- summarySE(data=garcia,measurevar="cumattack", groupvars=c("treatment","time"))
garcia_stats

Original authors did not use log10 transform but worth tryings

garcia$lcumattack <- log10(garcia$cumattack)
garcia1log.aov <- aov(lcumattack~treatment*time, garcia)
plot(garcia1log.aov)

boxplot(lcumattack~treatment*time,data=garcia)

garcialog_stats <- summarySE(data=garcia,measurevar="lcumattack", groupvars=c("treatment","time"))
garcialog_stats

Log10 overcorrects - analyse original data and try different variances with lme First analyse with time as categorical factor

Run as simple split-plot with polynomial contrasts

fully balanced so all SS are OK

garcia2.aov <- aov(cumattack~treatment*timefac+Error(toad), garcia)
summary(garcia2.aov)

Error: toad
          Df Sum Sq Mean Sq F value  Pr(>F)   
treatment  1   1647  1646.9   8.593 0.00978 **
Residuals 16   3066   191.7                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: Within
                  Df Sum Sq Mean Sq F value   Pr(>F)    
timefac            4   3525   881.2   91.81  < 2e-16 ***
treatment:timefac  4    579   144.7   15.08 9.82e-09 ***
Residuals         64    614     9.6                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

summary(garcia2.aov, split=list(timefac=list(linear=1, quadratic=2, cubic=3)))

Error: toad
          Df Sum Sq Mean Sq F value  Pr(>F)   
treatment  1   1647  1646.9   8.593 0.00978 **
Residuals 16   3066   191.7                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: Within
                               Df Sum Sq Mean Sq F value   Pr(>F)    
timefac                         4   3525     881  91.813  < 2e-16 ***
  timefac: linear               1   3494    3494 364.023  < 2e-16 ***
  timefac: quadratic            1     27      27   2.848   0.0963 .  
  timefac: cubic                1      2       2   0.255   0.6151    
treatment:timefac               4    579     145  15.077 9.82e-09 ***
  treatment:timefac: linear     1    572     572  59.647 9.92e-11 ***
  treatment:timefac: quadratic  1      1       1   0.070   0.7924    
  treatment:timefac: cubic      1      4       4   0.422   0.5183    
Residuals                      64    614      10                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

emmeans(garcia2.aov, ~timefac|treatment)

treatment = Leptin:
 timefac emmean   SE   df lower.CL upper.CL
 3         6.89 2.26 22.8     2.21     11.6
 6        10.78 2.26 22.8     6.10     15.5
 9        13.67 2.26 22.8     8.99     18.3
 12       15.89 2.26 22.8    11.21     20.6
 15       17.44 2.26 22.8    12.77     22.1

treatment = Saline:
 timefac emmean   SE   df lower.CL upper.CL
 3         8.89 2.26 22.8     4.21     13.6
 6        14.78 2.26 22.8    10.10     19.5
 9        22.44 2.26 22.8    17.77     27.1
 12       28.22 2.26 22.8    23.54     32.9
 15       33.11 2.26 22.8    28.43     37.8

Warning: EMMs are biased unless design is perfectly balanced 
Confidence level used: 0.95 

Run with ez to get GG and HF adjustments

ezgarcia1 <- ezANOVA(data=garcia, dv=cumattack, wid=toad, within=timefac, between=treatment, type=3, detailed=TRUE)

Warning: Converting "treatment" to factor for ANOVA.

print(ezgarcia1)

$ANOVA
             Effect DFn DFd        SSn       SSd          F            p p<.05       ges
1       (Intercept)   1  16 26660.0111 3066.4444 139.105790 2.642036e-09     * 0.8786887
2         treatment   1  16  1646.9444 3066.4444   8.593376 9.780788e-03     * 0.3091338
3           timefac   4  64  3524.6000  614.2222  91.813025 8.657540e-26     * 0.4891700
4 treatment:timefac   4  64   578.7778  614.2222  15.076700 9.818711e-09     * 0.1358810

$`Mauchly's Test for Sphericity`
             Effect          W            p p<.05
3           timefac 0.03247946 1.720133e-07     *
4 treatment:timefac 0.03247946 1.720133e-07     *

$`Sphericity Corrections`
             Effect      GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
3           timefac 0.371348 6.293100e-11         * 0.4004872 1.275427e-11         *
4 treatment:timefac 0.371348 1.881095e-04         * 0.4004872 1.182037e-04         *

Run as REML mixed effects

Use lme4 using time as factor - use anova and Anova to look at fixed effects (Type 3 is OK to decide whether to keep intyeraction)

First random intercepts only to match aov model

garcia1.lmer <- lmer(cumattack~treatment*timefac+(1|toad), garcia)
summary(garcia1.lmer, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: cumattack ~ treatment * timefac + (1 | toad)
   Data: garcia

REML criterion at convergence: 488

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.49115 -0.47848  0.02588  0.47055  2.87421 

Random effects:
 Groups   Name        Variance Std.Dev.
 toad     (Intercept) 36.411   6.034   
 Residual              9.597   3.098   
Number of obs: 90, groups:  toad, 18

Fixed effects:
                     Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)           17.2111     1.4593 16.0000  11.794 2.64e-09 ***
treatment1            -4.2778     1.4593 16.0000  -2.931  0.00978 ** 
timefac.L             13.9316     0.7302 64.0000  19.079  < 2e-16 ***
timefac.Q             -1.2324     0.7302 64.0000  -1.688  0.09633 .  
timefac.C             -0.3689     0.7302 64.0000  -0.505  0.61512    
timefac^4              0.2590     0.7302 64.0000   0.355  0.72401    
treatment1:timefac.L  -5.6394     0.7302 64.0000  -7.723 9.92e-11 ***
treatment1:timefac.Q  -0.1930     0.7302 64.0000  -0.264  0.79236    
treatment1:timefac.C   0.4743     0.7302 64.0000   0.650  0.51827    
treatment1:timefac^4  -0.2988     0.7302 64.0000  -0.409  0.68375    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) trtmn1 tmfc.L tmfc.Q tmfc.C tmfc^4 tr1:.L tr1:.Q tr1:.C
treatment1  0.000                                                         
timefac.L   0.000  0.000                                                  
timefac.Q   0.000  0.000  0.000                                           
timefac.C   0.000  0.000  0.000  0.000                                    
timefac^4   0.000  0.000  0.000  0.000  0.000                             
trtmnt1:t.L 0.000  0.000  0.000  0.000  0.000  0.000                      
trtmnt1:t.Q 0.000  0.000  0.000  0.000  0.000  0.000  0.000               
trtmnt1:t.C 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000        
trtmnt1:t^4 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 

anova(garcia1.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)    
treatment           82.5   82.47     1    16  8.5934  0.009781 ** 
timefac           3524.6  881.15     4    64 91.8130 < 2.2e-16 ***
treatment:timefac  578.8  144.69     4    64 15.0767 9.819e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get var comps

garcia1.ci <- confint.merMod(garcia1.lmer)

Computing profile confidence intervals ...

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

                          2.5 %       97.5 %
.sig01                17.178221  69.20280818
.sigma                 6.258163  12.05099182
(Intercept)          206.327227 402.32844979
treatment1            50.762687   2.04708022
timefac.L            157.856396 234.06210257
timefac.Q              6.759322   0.01825933
timefac.C              3.015192   0.99713403
timefac^4              1.228844   2.64538702
treatment1:timefac.L  49.096557  18.24910212
treatment1:timefac.Q   2.435224   1.37939474
treatment1:timefac.C   0.797729   3.39237469
treatment1:timefac^4   2.776574   1.14210094

Try random intercepts and slopes - doesn’t work as too many parameters (too many df used for toad*time interaction)

garcia1a.lmer <- lmer(cumattack~treatment*timefac+(timefac|toad), garcia)

Error: number of observations (=90) <= number of random effects (=90) for term (timefac | toad); the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable

Run as mixed effects with nlme using time as factor - random intercepts only to match above - matches above

garcia1.lme <- lme(cumattack~treatment*timefac, random=~1|toad, method="REML",garcia)
summary(garcia1.lme)

Linear mixed-effects model fit by REML
  Data: garcia 

Random effects:
 Formula: ~1 | toad
        (Intercept) Residual
StdDev:    6.034162 3.097938

Fixed effects:  cumattack ~ treatment * timefac 

 Correlation: 
                     (Intr) trtmn1 tmfc.L tmfc.Q tmfc.C tmfc^4 tr1:.L tr1:.Q tr1:.C
treatment1           0                                                             
timefac.L            0      0                                                      
timefac.Q            0      0      0                                               
timefac.C            0      0      0      0                                        
timefac^4            0      0      0      0      0                                 
treatment1:timefac.L 0      0      0      0      0      0                          
treatment1:timefac.Q 0      0      0      0      0      0      0                   
treatment1:timefac.C 0      0      0      0      0      0      0      0            
treatment1:timefac^4 0      0      0      0      0      0      0      0      0     

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.49114675 -0.47847867  0.02588476  0.47055427  2.87421248 

Number of Observations: 90
Number of Groups: 18 

anova(garcia1.lme, type="marginal")

Compare above model to one that allows group-specific variances to differ - refit with ML

garcia2.lme <- lme(cumattack~treatment*timefac, random=~1|toad, method="ML",garcia)
garcia3.lme <- lme(cumattack~treatment*timefac, random=~1|toad, weights = varIdent(form= ~ 1 | treatment*timefac),method="ML",garcia)
anova(garcia2.lme, garcia3.lme)

AICc(garcia2.lme, garcia3.lme)

Refit unequal variances model with REML

garcia4.lme <- lme(cumattack~treatment*timefac, random=~1|toad, weights = varIdent(form= ~ 1 | treatment*timefac),method="REML",garcia)
anova(garcia4.lme, type="marginal")

Second, analyse with time as continuous covariate Plot cumattack-time relationships for each toad - all relationships positive but some differences in slopes

xyplot(cumattack~time|toad, groups=treatment, type=c("p","r"), auto.key=T, garcia)

Now examine individual slopes (and intercepts) - again, variation in slopes

list_garcia <- lmList(cumattack~time|toad, garcia)
summary(list_garcia)

Call:
  Model: cumattack ~ time | toad 
   Data: garcia 

Coefficients:
   (Intercept) 
   Estimate Std. Error     t value     Pr(>|t|)
1       3.3   1.767138  1.86742599 6.727360e-02
2      -3.0   1.767138 -1.69765999 9.532860e-02
3       2.8   1.767138  1.58448266 1.189234e-01
4       9.9   1.767138  5.60227798 7.342449e-07
5      -2.5   1.767138 -1.41471666 1.628923e-01
6       9.1   1.767138  5.14956865 3.765460e-06
7       0.1   1.767138  0.05658867 9.550817e-01
8       6.7   1.767138  3.79144065 3.795788e-04
9       7.2   1.767138  4.07438399 1.522204e-04
10      1.7   1.767138  0.96200733 3.403341e-01
11      3.0   1.767138  1.69765999 9.532860e-02
12      2.3   1.767138  1.30153933 1.985994e-01
13     11.5   1.767138  6.50769665 2.590922e-08
14      5.1   1.767138  2.88602199 5.597604e-03
15      3.8   1.767138  2.15036933 3.602079e-02
16      7.9   1.767138  4.47050465 4.042883e-05
17      4.0   1.767138  2.26354666 2.764355e-02
18     -1.0   1.767138 -0.56588666 5.738150e-01
   time 
    Estimate Std. Error   t value     Pr(>|t|)
1  1.2333333  0.1776041  6.944285 5.074414e-09
2  1.1333333  0.1776041  6.381235 4.149348e-08
3  3.0666667  0.1776041 17.266870 1.206819e-23
4  1.7666667  0.1776041  9.947219 8.244742e-14
5  2.3000000  0.1776041 12.950153 3.392342e-18
6  0.7666667  0.1776041  4.316718 6.804628e-05
7  2.0333333  0.1776041 11.448686 4.591067e-16
8  1.1666667  0.1776041  6.568918 2.062162e-08
9  0.6000000  0.1776041  3.378301 1.359576e-03
10 0.3666667  0.1776041  2.064517 4.378554e-02
11 2.0666667  0.1776041 11.636369 2.446958e-16
12 3.3000000  0.1776041 18.580654 4.030655e-25
13 0.9666667  0.1776041  5.442818 1.310882e-06
14 0.4333333  0.1776041  2.439884 1.800572e-02
15 0.8666667  0.1776041  4.879768 9.796650e-06
16 1.7000000  0.1776041  9.571852 3.144289e-13
17 1.2000000  0.1776041  6.756601 1.023403e-08
18 1.4666667  0.1776041  8.258068 3.789301e-11

Residual standard error: 1.6849 on 54 degrees of freedom

Run as mixed effects with lme4 using time as continuous covariate

First ML model with random slopes and intercepts

garcia2.lmer <- lmer(cumattack~treatment+time+treatment*time+(time|toad), REML=FALSE, garcia)

Second just random intercepts

garcia3.lmer <- lmer(cumattack~treatment*time+(1|toad), REML=FALSE, garcia)

Test random slopes term - note the comparison is based on ML fits

anova(garcia3.lmer,garcia2.lmer)

Data: garcia
Models:
garcia3.lmer: cumattack ~ treatment * time + (1 | toad)
garcia2.lmer: cumattack ~ treatment + time + treatment * time + (time | toad)
             npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
garcia3.lmer    6 518.48 533.48 -253.24   506.48                         
garcia2.lmer    8 465.00 485.00 -224.50   449.00 57.483  2  3.294e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

AICc(garcia3.lmer,garcia2.lmer)

Use random slopes model - refit using REML

garcia4.lmer <- lmer(cumattack~treatment+time+treatment*time+(time|toad), REML=TRUE, garcia)
summary(garcia4.lmer, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: cumattack ~ treatment + time + treatment * time + (time | toad)
   Data: garcia

REML criterion at convergence: 449.9

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.89438 -0.62632  0.07927  0.48621  1.75457 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 toad     (Intercept) 13.9359  3.7331        
          time         0.3144  0.5607   -0.04
 Residual              2.8390  1.6849        
Number of obs: 90, groups:  toad, 18

Fixed effects:
                Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)       3.9944     0.9735 16.0000   4.103 0.000831 ***
treatment1        1.0722     0.9735 16.0000   1.101 0.287014    
time              1.4685     0.1386 16.0000  10.592 1.23e-08 ***
treatment1:time  -0.5944     0.1386 16.0000  -4.288 0.000565 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) trtmn1 time  
treatment1   0.000              
time        -0.154  0.000       
tretmnt1:tm  0.000 -0.154  0.000

anova(garcia4.lmer, 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)    
treatment        3.44    3.44     1    16   1.2131 0.2870140    
time           318.53  318.53     1    16 112.2005 1.226e-08 ***
treatment:time  52.19   52.19     1    16  18.3848 0.0005649 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

get var comps

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

Computing profile confidence intervals ...

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

                              2.5 %     97.5 %
sd_(Intercept)|toad       5.1156549 28.4594487
cor_time.(Intercept)|toad 0.2677897  0.2562486
sd_time|toad              0.1382368  0.6086524
sigma                     1.9907550  4.2456171
(Intercept)               4.3897541 34.7358522
treatment1                0.6840088  8.8297571
time                      1.4352649  3.0241614
treatment1:time           0.7481169  0.1049445

Compare to random intercepts for interest - refit using REML

garcia5.lmer <- lmer(cumattack~treatment*time+(1|toad), REML=TRUE, garcia)
summary(garcia5.lmer, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: cumattack ~ treatment * time + (1 | toad)
   Data: garcia

REML criterion at convergence: 508.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.7723 -0.5665  0.1216  0.5436  2.6657 

Random effects:
 Groups   Name        Variance Std.Dev.
 toad     (Intercept) 36.469   6.039   
 Residual              9.308   3.051   
Number of obs: 90, groups:  toad, 18

Fixed effects:
                Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)       3.9944     1.6109 23.5011   2.480   0.0207 *  
treatment1        1.0722     1.6109 23.5011   0.666   0.5121    
time              1.4685     0.0758 70.0000  19.374  < 2e-16 ***
treatment1:time  -0.5944     0.0758 70.0000  -7.842 3.56e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) trtmn1 time  
treatment1   0.000              
time        -0.423  0.000       
tretmnt1:tm  0.000 -0.423  0.000

anova(garcia5.lmer, 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)    
treatment         4.1     4.1     1 23.501   0.4431    0.5121    
time           3493.6  3493.6     1 70.000 375.3426 < 2.2e-16 ***
treatment:time  572.4   572.4     1 70.000  61.5023 3.565e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Get CIs as well
garcia5.ci <- confint.merMod(garcia5.lmer, oldNames=FALSE)

Computing profile confidence intervals ...

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

                         2.5 %     97.5 %
sd_(Intercept)|toad 17.0727070 69.0998289
sigma                6.6384482 12.7832490
(Intercept)          0.7734961 50.5436071
treatment1           4.1727702 17.5324806
time                 1.7425498  2.6146253
treatment1:time      0.5519097  0.1989005

Now allow covariances to differ using nlme - compare default to AR(1) - use random slopes model

lmeControl(maxIter = 100)  #Controls lme optimization process

$maxIter
[1] 100

$msMaxIter
[1] 50

$tolerance
[1] 1e-06

$niterEM
[1] 25

$msMaxEval
[1] 200

$msTol
[1] 1e-07

$msVerbose
[1] FALSE

$returnObject
[1] FALSE

$gradHess
[1] TRUE

$apVar
[1] TRUE

$.relStep
[1] 6.055454e-06

$opt
[1] "nlminb"

$optimMethod
[1] "BFGS"

$minAbsParApVar
[1] 0.05

$natural
[1] TRUE

$sigma
[1] 0

$allow.n.lt.q
[1] FALSE

garcia5.lme <- lme(cumattack~treatment*time, random=~time|toad, method="ML",garcia)
garcia6.lme <- lme(cumattack~treatment*time, random=~time|toad, correlation=corAR1(, form=~1|toad),method="ML",garcia)

Error in lme.formula(cumattack ~ treatment * time, random = ~time | toad,  : 
  nlminb problem, convergence error code = 1
  message = iteration limit reached without convergence (10)

couldn’t get convergence with AR(1) with random slopes even when increasing max iterations to 100 so use random intercept model. For the rand. int. model, we still need to keep the number of iterations above the default

garcia7.lme <- lme(cumattack~treatment*time, random=~1|toad, method="ML",garcia)
garcia8.lme <- lme(cumattack~treatment*time, random=~1|toad, correlation=corAR1(, form=~1|toad),method="ML",garcia)

Compare models

anova(garcia7.lme,garcia8.lme)

AICc(garcia7.lme,garcia8.lme)