############################################################################################################ #################################### ABOVE CHANCE ANALYSIS ################################################# ############################################################################################################ library(glmmTMB); library(emmeans) chance <- read.csv("Data.csv") chance$op <- factor(chance$op) m <- glmmTMB(cbind(scorefirst, trialfirst - scorefirst) ~ op + (1|grpid) + (1|individ), family = binomial, data = chance) probcorrect <- emmeans(m, ~ op, type = "response") test(probcorrect, null = 0, side = "greater") #op prob SE df null z.ratio p.value #1 0.956 0.0160 Inf 0.5 8.088 <0.0001 #2 0.742 0.0628 Inf 0.5 3.221 0.0013 #3 0.811 0.0560 Inf 0.5 3.987 <0.0001 table(chance$op) # 1 2 3 #38 33 32 ############################################################################################################ ####################################### FIRST ATTEMPTS ANALYSIS ############################################ ############################################################################################################ install.packages("lme4") install.packages("glmmTMB") install.packages("DHARMa") install.packages("MuMIn") library(lme4) library(MuMIn) library(lme4) library(glmmTMB) library(DHARMa) first <- read.csv("Data.csv") str(first) hist(first$proportion) first$grpid <- factor(first$grpid) first$individ <- factor(first$individ) first$op <- factor(first$op) first$mass <- as.numeric(first$mass) first$score <- as.numeric(first$score) first$scorefirst <- as.numeric(first$scorefirst) first$trialfirst <- as.numeric(first$trialfirst) first$latency <- as.numeric(first$latency) first$maxgrp <- as.numeric(first$grp) ###########BASE MODEL#################### basefirst <- glmmTMB(cbind(scorefirst, trialfirst - scorefirst) ~ (1 | grpid) + (1 |individ), family = binomial, data = first) summary(basefirst) confint(basefirst, method = "Wald") AICc(basefirst) #AICc 411.5321 res <- simulateResiduals(basefirst) plot(res) testDispersion(res) #dispersion = 1.2742, p-value < 2.2e-16 ##############PHASE#################### opfirst <- glmmTMB(cbind(scorefirst, trialfirst - scorefirst) ~ op + (1 | grpid) + (1 |individ), family = binomial, data = first) summary(opfirst) confint(opfirst, method = "Wald") AICc(opfirst) #AICc 347.5744 #op2 -2.0233 0.2858 -7.080 1.44e-12 *** #op3 -1.6240 0.3240 -5.013 5.37e-07 *** res <- simulateResiduals(opfirst) plot(res) testDispersion(res) #dispersion = 1.0655, p-value = 0.288 ###########MASS#################### massfirst <- glmmTMB(cbind(scorefirst, trialfirst - scorefirst) ~ scale(mass) + (1 | grpid) + (1 | individ), family = binomial, data = first) summary(massfirst) confint(massfirst, method = "Wald") AICc(massfirst) # AICc 409.0076 #scale(mass) 0.5149 0.2348 2.193 0.0283 * res <- simulateResiduals(massfirst) plot(res) testDispersion(res) #dispersion = 1.2686, p-value < 2.2e-16 ###########GROUPSIZE############## groupfirst <- glmmTMB(cbind(scorefirst, trialfirst - scorefirst) ~ scale(grp) + (1 | grpid) + (1 | individ), family = binomial, data = first) summary(groupfirst) confint(groupfirst, method = "Wald") AICc(groupfirst) # AICc 413.5801 scale(grp) -0.08673 0.25296 -0.343 0.732 res <- simulateResiduals(groupfirst) plot(res) testDispersion(res) #dispersion = 1.2736, p-value < 2.2e-16 ###########GROUPDAY############## groupdayfirst <- glmmTMB(cbind(scorefirst, trialfirst - scorefirst) ~ scale(grpday) + (1 | grpid) + (1 | individ), family = binomial, data = first) summary(groupdayfirst) confint(groupdayfirst, method = "Wald") AICc(groupdayfirst) # AICc 412.8071 scale(grpday) 0.1449 0.1541 0.94 0.347 res <- simulateResiduals(groupdayfirst) plot(res) testDispersion(res) #dispersion = 1.2631, p-value < 2.2e-16 ###########AGE############## agefirst <- glmmTMB(cbind(scorefirst, trialfirst - scorefirst) ~ scale(age) + (1 | grpid) + (1 | individ), family = binomial, data = first) summary(agefirst) confint(agefirst, method = "Wald") AICc(agefirst) # AICc 412.9044 pv 0.373 -0.2725997 0.7273652 res <- simulateResiduals(agefirst) plot(res) testDispersion(res) #dispersion = 1.1292, p-value = 0.28 ###########SEX############## sexfirst <- glmmTMB(cbind(scorefirst, trialfirst - scorefirst) ~ sex + (1 | grpid) + (1 | individ), family = binomial, data = first) summary(sexfirst) confint(sexfirst, method = "Wald") AICc(sexfirst) # AICc 415.5144 #sexM 0.533 #sexU 0.934 res <- simulateResiduals(sexfirst) plot(res) testDispersion(res) #dispersion = 1.2746, p-value < 2.2e-16 ###########SEXRATIO############## sexratiofirst <- glmmTMB(cbind(scorefirst, trialfirst - scorefirst) ~ scale(sexratio) + (1 | grpid) + (1 | individ), family = binomial, data = first) summary(sexratiofirst) confint(sexratiofirst, method = "Wald") AICc(sexratiofirst) # AICc 413.4941 scale(sexratio) -0.1102 0.2435 -0.453 0.651 res <- simulateResiduals(sexratiofirst) plot(res) testDispersion(res) #dispersion = 1.2741, p-value < 2.2e-16 ###########TEMP############## tempfirst <- glmmTMB(cbind(scorefirst, trialfirst - scorefirst) ~ scale(avtemp) + (1 | grpid) + (1 | individ), family = binomial, data = first) summary(tempfirst) confint(tempfirst, method = "Wald") AICc(tempfirst) # AICc 411.1933 scale(avtemp) 0.2052 0.1301 1.578 0.115 res <- simulateResiduals(tempfirst) plot(res) testDispersion(res) #dispersion = 1.2328, p-value < 2.2e-16 ###########LATENCY############## first$zlat <- as.numeric(scale(first$latency)) latencyfirst <- glmmTMB(cbind(scorefirst, trialfirst - scorefirst) ~ scale(latency) + (1 | grpid) + (1 | individ), family = binomial, data = first) summary(latencyfirst) confint(latencyfirst, method = "Wald") AICc(latencyfirst) # AIC 404.8047 scale(latency) 0.4760 0.1813 2.626 0.00863 ** res <- simulateResiduals(latencyfirst) plot(res) testDispersion(res) #dispersion = 1.1886, p-value < 2.2e-16 ###########NEO############## neofirst <- glmmTMB(cbind(scorefirst, trialfirst - scorefirst) ~ scale(avneo) + (1 | grpid) + (1 | individ), family = binomial, data = first) summary(neofirst) confint(neofirst, method = "Wald") AICc(neofirst) # AIC 409.9793 scale(avneo) -0.4056 0.2017 -2.011 0.0443 * res <- simulateResiduals(neofirst) plot(res) testDispersion(res) #dispersion = 1.2787, p-value < 2.2e-16 ###########FULL############## fullfirst <- glmmTMB(cbind(scorefirst, trialfirst - scorefirst) ~ op + zlat + I(zlat^2) + (1 | grpid) + (1 |individ), family = binomial, data = first) summary(fullfirst) confint(fullfirst, method = "Wald") AICc(fullfirst) #AICc 350.0102 #op2 -2.12180 0.31035 -6.837 8.10e-12 *** #op3 -1.77167 0.37590 -4.713 2.44e-06 *** res <- simulateResiduals(fullfirst) plot(res) testDispersion(res) #dispersion = 1.0661, p-value = 0.264 ############################################################################################################ ################################# OVERALL ATTEMPTS ANALYSIS / TRIALS TO CRITERION ########################## ############################################################################################################ library(glmmTMB) library(MuMIn) library(dplyr) library(DHARMa) whole <- read.csv("Data.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) #178.1018 res <- simulateResiduals(base) plot(res) testDispersion(res) #dispersion = 0.99916, p-value = 0.936 ###########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 132.2989 #op2 0.26128 0.07121 3.67 p-value= 0.000243 *** #op3 -0.32779 0.06160 -5.32 p-value= 1.03e-07 *** r1 <- DHARMa::simulateResiduals(op, n=2000); plot(r1) plot(r1) testDispersion(r1) #dispersion = 0.95679, p-value = 0.909 ###########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 179.7761 scale(mass) -0.04573 p-value= 0.284 res <- simulateResiduals(mass) plot(res) testDispersion(res) #dispersion = 0.99575, p-value = 0.944 ###########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 180.91 scale(grp) 0.0005271 p-value= 0.992 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 180.8616 scale(grpday) -0.01009 p-value= 0.826 res <- simulateResiduals(grpday) plot(res) testDispersion(res) #dispersion = 0.97539, p-value = 0.96 ###########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 170.795 scale(age) -0.009607 p-value= 0.814 res <- simulateResiduals(age) plot(res) testDispersion(res) #dispersion = 0.97501, p-value = 0.96 ###########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 183.7791 #sexM 0.002531 p-value= 0.977 #sexU -0.018073 p-value= 0.950 res <- simulateResiduals(sex) plot(res) testDispersion(res) #dispersion = 0.97465, p-value = 0.96 ###########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 180.8903 scale(sexratio) 0.007158 p-value= 0.888 res <- simulateResiduals(sexratio) plot(res) testDispersion(res) #dispersion = 0.97618, p-value = 0.96 ###########TEMP############## temp <- glmmTMB(y_log ~ scale(avtemp) + op + (1|individ)+ (1|grpid), family = gaussian(), dispformula = ~ 0 + grpid_disp, data = dat_pass) summary(temp) confint(temp, method = "Wald") AICc(temp) #AICc 171.554 scale(avtemp) 0.11790 p-value= 0.00049 *** --> 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 = 0.9837, p-value = 0.952 ###########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 174.2795 #zlat 0.19157 p-value= 0.000198 *** #I(zlat^2) -0.05240 p-value= 0.001867 ** r1 <- DHARMa::simulateResiduals(latency, n=2000); plot(r1) plot(r1) testDispersion(r1) #dispersion = 0.98509, p-value = 0.996 ###########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 179.5961 scale(avneo) 0.05264 p-value= 0.249 res <- simulateResiduals(neo) plot(res) testDispersion(res) #dispersion = 0.97447, p-value = 0.952 ###########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 120.1263 #zlat 0.09912 0.04402 2.25 0.024330 * #I(zlat^2) -0.03447 0.01002 -3.44 0.000582 *** #op2 0.24706 0.03089 8.00 1.25e-15 *** #op3 -0.32430 0.03820 -8.49 < 2e-16 *** r1 <- DHARMa::simulateResiduals(full, n=2000); plot(r1) plot(r1) testDispersion(r1) #dispersion = 0.95653, p-value = 0.904 #########AICc############### AICc(base, mass, op, age, sex, grp, grpday, temp, sexratio, latency, neo, full) # df AICc # df AICc #base 14 168.0925 #mass 15 169.8165 #op 16 121.6136 #age 15 170.7950 #sex 16 173.6575 #grp 15 170.8390 #grpday 15 170.7972 #temp 15 162.0459 #sexratio 15 168.2651 #latency 16 164.4859 #neo 15 169.6782 #full 18 109.5283 ############################################################################################################ ################################# PERFORMANCE CONSISTANCY/ASSOCIATIONS ##################################### ############################################################################################################ library(glmmTMB) library(dplyr) library(tidyr) dat <- first %>% mutate(op = factor(op, levels = c("1","2","3")), prop = scorefirst / trialfirst) %>% select(individ, grpid, op, score = scorefirst, trial = trialfirst, prop) wide <- dat %>% select(individ, grpid, op, score, trial, prop) %>% pivot_wider(names_from = op, values_from = c(score, trial, prop), names_vary = "slowest") #OP1 VS OP2 op12 <- wide %>% filter(!is.na(score_1), !is.na(score_2)) model12 <- glmmTMB( cbind(score_2, trial_2 - score_2) ~ scale(prop_1) + (1 | grpid) + (1 | individ), family = binomial, data = op12 ) summary(model12) # Estimate Std. Error z value Pr(>|z|) #scale(prop_1) 0.6638 0.1308 5.074 3.89e-07 *** #OP1 VS OP3 op13 <- wide %>% filter(!is.na(score_1), !is.na(score_3)) model13 <- glmmTMB( cbind(score_3, trial_3 - score_3) ~ scale(prop_1) + (1 | grpid) + (1 | individ), family = binomial, data = op13 ) summary(model13) # Estimate Std. Error z value Pr(>|z|) #scale(prop_1) 0.2145 0.1985 1.081 0.28 #OP2 VS OP3 op23 <- wide %>% filter(!is.na(score_2), !is.na(score_3)) model23 <- glmmTMB( cbind(score_3, trial_3 - score_3) ~ scale(prop_2) + (1 | grpid) + (1 | individ), family = binomial, data = op23 ) summary(model23) # Estimate Std. Error z value Pr(>|z|) #scale(prop_2) 0.1516 0.2109 0.719 0.472 #HOLM adjustments pval <- c( summary(model12)$coefficients$cond["scale(prop_1)", "Pr(>|z|)"], summary(model13)$coefficients$cond["scale(prop_1)", "Pr(>|z|)"], summary(model23)$coefficients$cond["scale(prop_2)", "Pr(>|z|)"] ) p.adjust(pval, method = "holm") #[1] 1.168021e-06 5.596806e-01 5.596806e-01 #0.000001168021 #0.5596806 #0.5596806 ############################################################################################################ ########################################## TASK DIFFICULTY ################################################# ############################################################################################################ levels(first$op) m <- glmmTMB(cbind(scorefirst, trialfirst - scorefirst) ~ op + (1|grpid) + (1|individ), family = binomial, data = first) emm <- emmeans(m, ~ op) con <- pairs(emm, adjust = "holm") sc <- summary(con, infer = c(TRUE, TRUE)) contrast <- data.frame( Contrast = sc$contrast, OR = exp(sc$estimate), `95% CI (lower)` = exp(sc$asymp.LCL), `95% CI (upper)` = exp(sc$asymp.UCL), z = sc$z.ratio, `p (Holm)` = sc$p.value ) contrast # Contrast OR X95..CI..lower. X95..CI..upper. z p..Holm. #1 op1 - op2 7.5631610 3.8158673 14.990407 7.080198 4.318474e-12 #2 op1 - op3 5.0733674 2.3358885 11.018958 5.012600 1.073991e-06 #3 op2 - op3 0.6707999 0.3749785 1.199995 -1.643527 1.002739e-01 ############################################################################################################ ########################################## POWER TESTING ################################################### ############################################################################################################ install.packages("pwr") library(pwr) #small pwr.f2.test( u = 3, v = 103 - 3 - 1, f2 = 0.02, sig.level = 0.05 ) # power = 0.1906988 #medium pwr.f2.test( u = 3, v = 103 - 3 - 1, f2 = 0.15, sig.level = 0.05 ) # power = 0.9141335 #large pwr.f2.test( u = 3, v = 103 - 3 - 1, f2 = 0.35, sig.level = 0.05 ) #power = 0.9995714