library(glmmTMB) library(MuMIn) library(dplyr) library(DHARMa) whole <- read.csv("Data2.csv") str(whole) whole <- whole %>% mutate( across(c(grpid, sex, individ, op), factor), across(c(mass, score, trialfirst, grpday, avtemp, scorefirst, trialfirst, latency), as.numeric), maxgrp = as.numeric(grp) ) dat_pass <- whole %>% mutate( y_log = log(score), y_sqrt = sqrt(score), opF = factor(op), ztemp = as.numeric(scale(avtemp)) ) %>% add_count(grpid, name = "n_g") %>% mutate( grpid_disp = factor(ifelse(n_g >= 6, as.character(grpid), "other")) ) ###########BASE MODEL#################### base <- glmmTMB(y_log ~ (1|individ) + (1|grpid), family = gaussian(),dispformula = ~ 0 + grpid_disp, data = dat_pass) summary(base) confint(base, method = "Wald") AICc(base) #108.3945 res <- simulateResiduals(base) plot(res) testDispersion(res) #dispersion = 1.0085, p-value = 0.896 ###########PHASE#################### op <- glmmTMB(y_log ~ op + (1|individ)+ (1|grpid), family = gaussian(),dispformula = ~ 0 + grpid_disp, data = dat_pass) summary(op) confint(op, method = "Wald") AICc(op) #AICc 47.5093 #op2 0.12573 0.04413 2.85 0.00439 ** #op3 -0.42539 0.04029 -10.56 < 2e-16 *** r1 <- DHARMa::simulateResiduals(op, n=2000); plot(r1) plot(r1) testDispersion(r1) #dispersion = 1.0155, p-value = 0.883 ###########MASS#################### mass <- glmmTMB(y_log ~ scale(mass) + (1|individ)+ (1|grpid), family = gaussian(), dispformula = ~ 0 + grpid_disp, data = dat_pass) summary(mass) confint(mass, method = "Wald") AICc(mass) #AICc 108.1614 cale(mass) -0.06321 p 0.077 . res <- simulateResiduals(mass) plot(res) testDispersion(res) #dispersion = 1.01, p-value = 0.88 ###########GROUPSIZE############## grp <- glmmTMB(y_log ~ scale(grp) + (1|individ)+ (1|grpid), family = gaussian(), dispformula = ~ 0 + grpid_disp, data = dat_pass) summary(grp) confint(grp, method = "Wald") AICc(grp) #AICc 110.844 scale(grp) -0.02453 p 0.525 res <- simulateResiduals(grp) plot(res) testDispersion(res) #dispersion = 0.97493, p-value = 0.96 ###########GROUPDAY############## grpday <- glmmTMB(y_log ~ scale(grpday) + (1|individ)+ (1|grpid), family = gaussian(), dispformula = ~ 0 + grpid_disp, data = dat_pass) summary(grpday) confint(grpday, method = "Wald") AICc(grpday) #AICc 111.101 scale(grpday) 0.01573 p 0.706 res <- simulateResiduals(grpday) plot(res) testDispersion(res) #dispersion = 1.0086, p-value = 0.88 ###########AGE############## age <- glmmTMB(y_log ~ scale(age) + (1|individ) + (1|grpid), family = gaussian(), dispformula = ~ 0 + grpid_disp, data = dat_pass) summary(age) confint(age, method = "Wald") AICc(age) #AICc 111.2434 scale(age) -0.0006439 p 0.987 res <- simulateResiduals(age) plot(res) testDispersion(res) #dispersion = 1.0085, p-value = 0.896 ###########SEX############## sex <- glmmTMB(y_log ~ sex + (1|individ)+ (1|grpid), family = gaussian(), dispformula = ~ 0 + grpid_disp, data = dat_pass) summary(sex) confint(sex, method = "Wald") AICc(sex) #AICc 113.8874 #sexM -0.01593 0.07717 -0.21 0.836 #sexU 0.07145 0.17473 0.41 0.683 res <- simulateResiduals(sex) plot(res) testDispersion(res) #dispersion = 1.009, p-value = 0.912 ###########SEXRATIO############## sexratio <- glmmTMB(y_log ~ scale(sexratio) + (1|individ)+ (1|grpid), family = gaussian(), dispformula = ~ 0 + grpid_disp, data = dat_pass) summary(sexratio) confint(sexratio, method = "Wald") AICc(sexratio) #AICc 110.454 scale(sexratio) -0.03365 p 0.367 res <- simulateResiduals(sexratio) plot(res) testDispersion(res) #dispersion = 1.0071, p-value = 0.896 ###########TEMP############## temp <- glmmTMB(y_log ~ scale(avtemp) + (1|individ)+ (1|grpid), family = gaussian(), dispformula = ~ 0 + grpid_disp, data = dat_pass) summary(temp) confint(temp, method = "Wald") AICc(temp) #AICc 106.673 scale(avtemp) 0.08492 p 0.0203 * --> significance should not be considered (phase 3 #had significantly lower temperatures and was different in pass rates so when phase is added into the model, #the significance disappears) res <- simulateResiduals(temp) plot(res) testDispersion(res) #dispersion = 1.0159, p-value = 0.88 ###########LATENCY############## dat_pass$zlat <- as.numeric(scale(dat_pass$latency)) latency <- glmmTMB(y_log ~ zlat + I(zlat^2) + (1|individ)+ (1|grpid), family = gaussian(), dispformula = ~ 0 + grpid_disp, data = dat_pass) summary(latency) confint(latency, method = "Wald") AICc(latency) #AICc 101.7683 #zlat 0.19040 0.04815 3.95 7.68e-05 *** #I(zlat^2) -0.04699 0.01630 -2.88 0.00394 ** r1 <- DHARMa::simulateResiduals(latency, n=2000); plot(r1) plot(r1) testDispersion(r1) #dispersion = 1.0059, p-value = 0.92 ###########NEO############## neo <- glmmTMB(y_log ~ scale(avneo) + (1|individ)+ (1|grpid), family = gaussian(), dispformula = ~ 0 + grpid_disp, data = dat_pass) summary(neo) confint(neo, method = "Wald") AICc(neo) #AICc 109.3524 scale(avneo) 0.05709 p 0.164 res <- simulateResiduals(neo) plot(res) testDispersion(res) #dispersion = 1.0104, p-value = 0.896 ###########FULL############## full <- glmmTMB(y_log ~ zlat + I(zlat^2) + op + (1 |individ)+ (1|grpid), family = gaussian(), dispformula = ~ 0 + grpid_disp, data = dat_pass) summary(full) confint(full, method = "Wald") AICc(full) #AICc 42.62806 #zlat 0.066458 0.041668 1.59 0.11073 #I(zlat^2) -0.028579 0.009835 -2.91 0.00366 ** #op2 0.213204 0.032526 6.55 5.57e-11 *** #op3 -0.349694 0.036728 -9.52 < 2e-16 *** r1 <- DHARMa::simulateResiduals(full, n=2000); plot(r1) plot(r1) testDispersion(r1) #dispersion = 1.0239, p-value = 0.852 #########AICc############### AICc(base, mass, op, age, sex, grp, grpday, temp, sexratio, latency, neo, full) # df AICc #base 14 108.39452 #mass 15 108.16145 #op 16 47.50930 #age 15 111.24342 #sex 16 113.88739 #grp 15 110.84400 #grpday 15 111.10103 #temp 15 106.67297 #sexratio 15 110.45400 #latency 15 106.03012 #neo 15 109.35241 #full 18 42.62806