Model selection function: Perform model selection on a specified full model (LM, LMM, GLM, GLMM) by calculating the AICs and producing a table of the best models within a given deltaAIC. It also gives the final/reduced model as the model, within those, that has the most parameters and lowest AIC. Model selection is done sequentially using function MuMIn::dredge.

model_selection <- function (model, delta = 2){

  # Model selection
  AICtable_full <- MuMIn::dredge(model)

  top_models_index <- which(AICtable_full[,'delta'] <= delta)
  AICtable <- AICtable_full[top_models_index,]

  model_list <-MuMIn::get.models(AICtable, subset = delta <= delta)

  final_model_index <- which.max(unlist(lapply(lapply(model_list, coefficients), length)))

  final_model <- model_list[[final_model_index]]

  return(list(AIC_table = AICtable,
              final_model = final_model))
}

Import data:

#make_forage()
load(file = paste0(here::here('data'),'/data_forage.Rda'))

kable(head(data_forage))
## Warning in `[<-.data.frame`(`*tmp*`, , isn, value = structure(list(Litter =
## structure(c("3", : provided 9 variables to replace 8 variables
Population ID Sex Mom_ID Litter Treatment Pred_Tutor Pred Tutor_pop Birth.Day Day Date Time Water Weight Run Latency Attempts Notes Tutor event time2peck Predator
AH AHF2_1G4G.3.3 F AHF2_1G4G 3 pred-/same pred-/AH pred- AH 2-Feb 2 22-Mar 15:47 pred- 0.133 1 16 17 AH 1 16 pred-
AH AHF2_1G4G.3.3 F AHF2_1G4G 3 pred-/same pred-/AH pred- AH 2-Feb 1 21-Mar 15:06 pred+ 0.133 1 73 99 AH 1 73 pred-
AH AHF2_1G4G.6.2 F AHF2_1G4G 6 pred-/same pred-/AH pred- AH 14-Aug 1 29-Sep 15:17 pred- 0.144 1 35 166 first few in middle AH 1 35 pred-
AH AHF2_1G4G.6.2 F AHF2_1G4G 6 pred-/same pred-/AH pred- AH 14-Aug 2 30-Sep 15:06 pred+ 0.144 1 58 7 AH 1 58 pred-
AH AHF2_1G4P.2.1 F AHF2_1G4P 2 pred-/same pred-/AH pred- AH 15-Jan 1 2-Mar 15:30 pred+ 0.128 1 39 55 AH 1 39 pred-
AH AHF2_1G4P.2.1 F AHF2_1G4P 2 pred-/same pred-/AH pred- AH 15-Jan 2 3-Mar 16:04 pred- 0.128 1 310 72 AH 1 310 pred-

Remove AHF2_1O3P.1.3 and AHF2_1O3P.1.3 because they were assessed at 76 and 77 days old

data_forage <- data_forage[!data_forage$ID =="AHF2_1O3P.1.3",]

data_forage <- data_forage[!data_forage$ID =="AHF2_1O3P.1.4",]

#data_forage <- data_forage[data_forage$ID =="AHF2_1O3P.1.3",]

#data_forage <- data_forage[data_forage$ID =="AHF2_1O3P.1.4",]

#data_forage %>%  filter(ID=='AHF2_1O3P.1.3' | ID=='AHF2_1O3P.1.4')

# THis isn't working, I'll just manually delete these guys!

Ben Bolker (creater of LME4) recommends using glmmTMB instead of lme4::glmer.nb for negative binomial GLMMs. Nevertheless, he said that is smart to check model parameters with different implementations/algorithms to make sure the answers are consistent - he calls it the gold standard for addressing convergence warnings, so that is what we are going to do!

# data_f <- data_forage %>%
#   filter(Sex == 'F') %>%
#   na.omit
# 
# library(lme4)
# m1 <- glmer.nb(Attempts ~ (Population + Pred + Tutor_pop + Water)^2 + (1|ID) + (1|Mom_ID), data = data_f,  na.action = "na.fail", family = poisson(), control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)))
# # Fails to converge, but outpute suggests that it could be false +
# 
# library(glmmTMB)
# m2 <- glmmTMB(Attempts ~ (Population + Pred + Tutor_pop + Water)^2 + (1|ID) + (1|Mom_ID), na.action = "na.fail", family = 'nbinom2', data = data)
# 
# library(ggplot2)
# library(dotwhisker)
# library(broom.mixed) 
# dwplot(list(glmer.nb=model1,glmmTMB=model),effect="fixed") + theme_bw() +
#     geom_vline(xintercept=0,lty=2)

Coefficients are similar, so either either output could be used. But let’s continue with the glmmTMB since Bolker recommends that…

Foraging Attempts Analyses

Female Foraging Attempts
  1. Subset female data
data <- data_forage %>% filter(Sex == 'F') %>% na.omit
  1. Create generalized linear mixed model with all interactions
library(glmmTMB)
model_pecks <- glmmTMB(Attempts ~ (Population + Pred + Tutor_pop + Water)^2 + (1|ID) + (1|Mom_ID), na.action = "na.fail", family = 'nbinom2', data = data)
  1. Run model selection (Andres’s model selection function does not work with the glmmTMB, so need to do model selection manually by look at the AIC_table)
mp <- glmmTMB(Attempts ~ Population + Pred + Tutor_pop + Water + Population:Pred + Population:Water + Pred:Water + Tutor_pop:Water + (1|ID) + (1|Mom_ID), na.action = "na.fail", family = 'nbinom2', data = data)
  1. Summary of final generalized linear mixed model for female foraging attempts
a1 <- car::Anova(mp, type = 'III')

kable(a1)
Chisq Df Pr(>Chisq)
(Intercept) 222.3902196 1 0.0000000
Population 17.1216186 1 0.0000351
Pred 0.5657917 1 0.4519360
Tutor_pop 3.9085832 1 0.0480401
Water 4.8073447 1 0.0283387
Population:Pred 1.6992990 1 0.1923797
Population:Water 1.4278858 1 0.2321098
Pred:Water 6.2026415 1 0.0127560
Tutor_pop:Water 2.5092898 1 0.1131769 
kable(summary(mp)$coefficients)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.8994624 0.2614851 14.9127536 0.0000000
PopulationAL -1.3074853 0.3159835 -4.1378278 0.0000351
Predpred+ -0.2274158 0.3023378 -0.7521913 0.4519360
Tutor_popAL -0.4878996 0.2467861 -1.9770137 0.0480401
Waterpred+ -0.6624155 0.3021189 -2.1925658 0.0283387
PopulationAL:Predpred+ 0.5041141 0.3867176 1.3035716 0.1923797
PopulationAL:Waterpred+ -0.3813682 0.3191521 -1.1949418 0.2321098
Predpred+:Waterpred+ 0.7881241 0.3164508 2.4905103 0.0127560
Tutor_popAL:Waterpred+ 0.4976624 0.3141662 1.5840738 0.1131769

|| || || ||

|| || || ||
Male Foraging Attempts
  1. Subset male data
data2 <- data_forage %>% filter(Sex == 'M') %>% na.omit
  1. Create generalized linear mixed model with all interactions
model_pecks2 <- glmmTMB(Attempts ~ (Population + Pred + Tutor_pop + Water)^2 + (1|ID) + (1|Mom_ID), na.action = "na.fail", family = 'nbinom2', data = data2)
  1. Run model selection (Andres’s model selection function does not work with the glmmTMB, so need to do model selection manually by look at the AIC_table)
mp2 <- glmmTMB(Attempts ~ Population + Tutor_pop + Water + Population:Water + Tutor_pop:Water + (1|ID) + (1|Mom_ID), na.action = "na.fail", family = 'nbinom2', data = data2)
  1. Summary of final generalized linear mixed model for male foraging attempts
a2 <- car::Anova(mp2, type = 'III')

kable(a2)
Chisq Df Pr(>Chisq)
(Intercept) 125.9727200 1 0.0000000
Population 2.4194959 1 0.1198335
Tutor_pop 0.8703434 1 0.3508602
Water 2.7989793 1 0.0943243
Population:Water 0.4179748 1 0.5179493
Tutor_pop:Water 7.2645051 1 0.0070331 
kable(summary(mp2)$coefficients)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.9535167 0.2631487 11.2237569 0.0000000
PopulationAL -0.4703049 0.3023549 -1.5554729 0.1198335
Tutor_popAL -0.2455058 0.2631579 -0.9329220 0.3508602
Waterpred+ -0.4888291 0.2921845 -1.6730150 0.0943243
PopulationAL:Waterpred+ -0.2135711 0.3303448 -0.6465097 0.5179493
Tutor_popAL:Waterpred+ 0.8975093 0.3329936 2.6952746 0.0070331

|| || || ||

|| || || ||

Note that it is important to remove NAs (with na.omit) and set na.action = "na.fail" so that the AIC works with models that have the exact same number of datapoints.

Now we can use our home-made function model_selection.

Latency Analyses

Female Latency
  1. Subset female data
data <- data_forage %>% filter(Sex == 'F') %>% na.omit
  1. Create Cox proportional hazards model with all interactions
library(coxme)

ml1 <- coxme(Surv(Latency, event = event) ~ (Population + Pred + Tutor_pop + Water)^2 + (1|ID) + (1|Mom_ID), data = data, na.action = 'na.fail')
  1. Run model selection
models_ml1 <- model_selection(ml1)
## Registered S3 methods overwritten by 'MuMIn':
##   method        from 
##   formula.coxme coxme
##   logLik.coxme  coxme
##   logLik.lmekin coxme
kable(models_ml1$AIC_table)
Population Pred Tutor_pop Water Population:Pred Population:Tutor_pop Population:Water Pred:Tutor_pop Pred:Water Tutor_pop:Water df logLik AICc delta weight
6 + NA + NA NA NA NA NA NA NA 59 -1242.297 2633.595 0.000000 0.3105946
38 + NA + NA NA + NA NA NA NA 59 -1242.447 2634.997 1.401849 0.1540942
152 + + + NA + NA NA + NA NA 59 -1243.166 2635.017 1.422390 0.1525197
136 + + + NA NA NA NA + NA NA 59 -1242.315 2635.291 1.696209 0.1330046
24 + + + NA + NA NA NA NA NA 59 -1243.192 2635.329 1.734152 0.1305051
8 + + + NA NA NA NA NA NA NA 60 -1242.157 2635.509 1.914002 0.1192817
  1. Summary of final cox proportional hazards mixed model summary for female foraging latency
anova <- car::Anova(models_ml1$final_model, type = 'III')
summary(models_ml1$final_model)
## Cox mixed-effects model fit by maximum likelihood
##   Data: data
##   events, n = 261, 299
##   Iterations= 18 95 
##                     NULL Integrated    Fitted
## Log-likelihood -1306.234  -1279.864 -1214.361
## 
##                    Chisq    df          p   AIC     BIC
## Integrated loglik  52.74  7.00 4.1745e-09 38.74   13.79
##  Penalized loglik 183.75 59.34 1.2101e-14 65.07 -146.45
## 
## Model:  Surv(Latency, event = event) ~ Population + Pred + Tutor_pop +      (1 | ID) + (1 | Mom_ID) + Population:Pred + Pred:Tutor_pop 
## Fixed coefficients
##                              coef exp(coef)  se(coef)     z       p
## PopulationAL           -1.0762220 0.3408809 0.2495925 -4.31 1.6e-05
## Predpred+               0.1038197 1.1094004 0.2677995  0.39 7.0e-01
## Tutor_popAL            -0.1779397 0.8369929 0.2361161 -0.75 4.5e-01
## PopulationAL:Predpred+  0.2095678 1.2331449 0.3324783  0.63 5.3e-01
## Predpred+:Tutor_popAL  -0.3081240 0.7348242 0.3300455 -0.93 3.5e-01
## 
## Random effects
##  Group  Variable  Std Dev    Variance  
##  ID     Intercept 0.60680167 0.36820827
##  Mom_ID Intercept 0.13055081 0.01704351
kable(anova)
Df Chisq Pr(>Chisq)
Population 1 18.5926208 0.0000162
Pred 1 0.1502933 0.6982552
Tutor_pop 1 0.5679292 0.4510830
Population:Pred 1 0.3973035 0.5284852
Pred:Tutor_pop 1 0.8715727 0.3505203
Male Latency
  1. Subset male data
data2 <- data_forage %>% filter(Sex == 'M') %>% na.omit
  1. Create Cox proportional hazards model with all interactions
ml2 <- coxme(Surv(Latency, event = event) ~ (Population + Pred + Tutor_pop + Water)^2 + (1|ID) + (1|Mom_ID), data = data2, na.action = 'na.fail')
  1. Run model selection
models_ml2 <- model_selection(ml2)

kable(models_ml2$AIC_table)
Population Pred Tutor_pop Water Population:Pred Population:Tutor_pop Population:Water Pred:Tutor_pop Pred:Water Tutor_pop:Water df logLik AICc delta weight
28 + + NA + + NA NA NA NA NA 4 -1041.389 2091.063 0.0000000 0.2556563
92 + + NA + + NA + NA NA NA 5 -1040.654 2091.676 0.6127019 0.1881958
544 + + + + + NA NA NA NA + 6 -1039.616 2091.730 0.6668614 0.1831679
20 + + NA NA + NA NA NA NA NA 3 -1042.875 2091.965 0.9014543 0.1628952
608 + + + + + NA + NA NA + 7 -1039.069 2092.761 1.6977908 0.1093921
284 + + NA + + NA NA NA + NA 5 -1041.280 2092.927 1.8635246 0.1006926
  1. Summary of final cox proportional hazards mixed model summary for male foraging latency
anova2 <- car::Anova(models_ml2$final_model, type = 'III')
summary(models_ml2$final_model)
## Cox mixed-effects model fit by maximum likelihood
##   Data: data2
##   events, n = 217, 252
##   Iterations= 7 37 
##                     NULL Integrated    Fitted
## Log-likelihood -1052.964   -1039.11 -1039.029
## 
##                   Chisq   df          p   AIC    BIC
## Integrated loglik 27.71 9.00 0.00106710  9.71 -20.71
##  Penalized loglik 27.87 7.08 0.00024696 13.72 -10.20
## 
## Model:  Surv(Latency, event = event) ~ Population + Pred + Tutor_pop +      Water + (1 | ID) + (1 | Mom_ID) + Population:Pred + Population:Water +      Tutor_pop:Water 
## Fixed coefficients
##                               coef exp(coef)  se(coef)     z      p
## PopulationAL            -0.1256642 0.8819109 0.2263486 -0.56 0.5800
## Predpred+                0.5551293 1.7421662 0.1918494  2.89 0.0038
## Tutor_popAL             -0.2330738 0.7920951 0.1903984 -1.22 0.2200
## Waterpred+              -0.3343386 0.7158113 0.2408074 -1.39 0.1700
## PopulationAL:Predpred+  -0.5462639 0.5791094 0.2754232 -1.98 0.0470
## PopulationAL:Waterpred+ -0.2857480 0.7514520 0.2739933 -1.04 0.3000
## Tutor_popAL:Waterpred+   0.4868620 1.6272021 0.2742041  1.78 0.0760
## 
## Random effects
##  Group  Variable  Std Dev      Variance    
##  ID     Intercept 1.899319e-02 3.607414e-04
##  Mom_ID Intercept 4.456235e-03 1.985803e-05
kable(anova2)
Df Chisq Pr(>Chisq)
Population 1 0.3082247 0.5787716
Pred 1 8.3727392 0.0038089
Tutor_pop 1 1.4985114 0.2209005
Water 1 1.9276731 0.1650132
Population:Pred 1 3.9337263 0.0473270
Population:Water 1 1.0876433 0.2969940
Tutor_pop:Water 1 3.1525630 0.0758078

Post-hoc analyses

Male Latency

library(emmeans)

ML <- coxph(models_ml2$final_model$formulaList$fixed, data = data2, na.action = 'na.fail')

emmML1 <- emmeans(ML, specs = pairwise ~ Population:Pred)
emmML2 <- emmeans(ML, specs = pairwise ~ Population|Pred)
emmML3 <- emmeans(ML, specs = pairwise ~ Pred|Population)

emmML1$emmeans
##  Population Pred  emmean    SE  df asymp.LCL asymp.UCL
##  AH         pred- -0.162 0.138 Inf   -0.4319    0.1080
##  AL         pred- -0.430 0.234 Inf   -0.8889    0.0283
##  AH         pred+  0.393 0.233 Inf   -0.0637    0.8499
##  AL         pred+ -0.422 0.249 Inf   -0.9092    0.0660
## 
## Results are averaged over the levels of: Tutor_pop, Water 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95
emmML2$emmeans
## Pred = pred-:
##  Population emmean    SE  df asymp.LCL asymp.UCL
##  AH         -0.162 0.138 Inf   -0.4319    0.1080
##  AL         -0.430 0.234 Inf   -0.8889    0.0283
## 
## Pred = pred+:
##  Population emmean    SE  df asymp.LCL asymp.UCL
##  AH          0.393 0.233 Inf   -0.0637    0.8499
##  AL         -0.422 0.249 Inf   -0.9092    0.0660
## 
## Results are averaged over the levels of: Tutor_pop, Water 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95
emmML3$emmeans
## Population = AH:
##  Pred  emmean    SE  df asymp.LCL asymp.UCL
##  pred- -0.162 0.138 Inf   -0.4319    0.1080
##  pred+  0.393 0.233 Inf   -0.0637    0.8499
## 
## Population = AL:
##  Pred  emmean    SE  df asymp.LCL asymp.UCL
##  pred- -0.430 0.234 Inf   -0.8889    0.0283
##  pred+ -0.422 0.249 Inf   -0.9092    0.0660
## 
## Results are averaged over the levels of: Tutor_pop, Water 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95

Male Foraging Attempts:

emmMF1 <- emmeans(model_pecks2, specs = pairwise~Tutor_pop:Water)
emmMF2 <- emmeans(model_pecks2, specs = pairwise~Tutor_pop|Water)
emmMF3 <- emmeans(model_pecks2, specs = pairwise~Water|Tutor_pop)

emmMF1$emmeans
##  Tutor_pop Water emmean    SE  df lower.CL upper.CL
##  AH        pred-   2.73 0.202 238     2.33     3.13
##  AL        pred-   2.49 0.208 238     2.08     2.90
##  AH        pred+   2.11 0.212 238     1.69     2.53
##  AL        pred+   2.77 0.209 238     2.35     3.18
## 
## Results are averaged over the levels of: Population, Pred 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95
emmMF2$emmeans
## Water = pred-:
##  Tutor_pop emmean    SE  df lower.CL upper.CL
##  AH          2.73 0.202 238     2.33     3.13
##  AL          2.49 0.208 238     2.08     2.90
## 
## Water = pred+:
##  Tutor_pop emmean    SE  df lower.CL upper.CL
##  AH          2.11 0.212 238     1.69     2.53
##  AL          2.77 0.209 238     2.35     3.18
## 
## Results are averaged over the levels of: Population, Pred 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95
emmMF3$emmeans
## Tutor_pop = AH:
##  Water emmean    SE  df lower.CL upper.CL
##  pred-   2.73 0.202 238     2.33     3.13
##  pred+   2.11 0.212 238     1.69     2.53
## 
## Tutor_pop = AL:
##  Water emmean    SE  df lower.CL upper.CL
##  pred-   2.49 0.208 238     2.08     2.90
##  pred+   2.77 0.209 238     2.35     3.18
## 
## Results are averaged over the levels of: Population, Pred 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95

Female Foraging Attempts:

emmFF1 <- emmeans(model_pecks, specs = pairwise~Water:Pred)
emmFF2 <- emmeans(model_pecks, specs = pairwise~Water|Pred)
emmFF3 <- emmeans(model_pecks, specs = pairwise~Pred|Water)

emmFF1$emmeans
##  Water Pred  emmean    SE  df lower.CL upper.CL
##  pred- pred-   3.00 0.177 285     2.65     3.35
##  pred+ pred-   2.40 0.193 285     2.02     2.78
##  pred- pred+   3.02 0.190 285     2.65     3.39
##  pred+ pred+   3.21 0.181 285     2.85     3.56
## 
## Results are averaged over the levels of: Population, Tutor_pop 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95
emmFF2$emmeans
## Pred = pred-:
##  Water emmean    SE  df lower.CL upper.CL
##  pred-   3.00 0.177 285     2.65     3.35
##  pred+   2.40 0.193 285     2.02     2.78
## 
## Pred = pred+:
##  Water emmean    SE  df lower.CL upper.CL
##  pred-   3.02 0.190 285     2.65     3.39
##  pred+   3.21 0.181 285     2.85     3.56
## 
## Results are averaged over the levels of: Population, Tutor_pop 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95
emmFF3$emmeans
## Water = pred-:
##  Pred  emmean    SE  df lower.CL upper.CL
##  pred-   3.00 0.177 285     2.65     3.35
##  pred+   3.02 0.190 285     2.65     3.39
## 
## Water = pred+:
##  Pred  emmean    SE  df lower.CL upper.CL
##  pred-   2.40 0.193 285     2.02     2.78
##  pred+   3.21 0.181 285     2.85     3.56
## 
## Results are averaged over the levels of: Population, Tutor_pop 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95

Figures

Female Latency

Only population origin had a significant effect on female latency to forage. However, tutor population and the interaction between population and tutor population were included in the ‘best’ model.

Option 1

Option 2 Currently in manuscript version

Option 3

Option 4

Male Latency

Predation and there interaction between Predation and Population origin had a significant effect on male latency to forage.

Option 1 Currently in manuscript version

Option 2

Option 3

Option 4

Option 5

Option 6
This includes all the variables for latency that are in the model, but it is just so hard to comprehend, that I think we go with the simplified version above that shows the significant relationships

Female Attempts

Population, Tutor Population, and Test Water, and the interaction between Rearing Water and Test Water had significant effects on female foraging attempts.

Option 1 Currently in manuscript version

Option 2

Option 3

Male Attempts

Only the interaction between Tutor Population and Test Water had significant effects on male foraging attempts.

Option 1 Currently in manuscript version