Heatwave Foraging Analysis
Steph Tran
2026-09-06
Loading in the data
data <- read.csv("~/Desktop/forage.csv") %>%
mutate(foraging_time = lubridate::ms(time),
foraging_time_sec = as.numeric(time),
latency_surv = Surv(latency, peck>0),
population = relevel(factor(population), ref = "LU"),
stream = factor(stream),
sex = factor(sex),
treatment = factor(treatment))## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `foraging_time = lubridate::ms(time)`.
## Caused by warning in `.parse_hms()`:
## ! Some strings failed to parse
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.#check to see if data loaded in correctly
head(data)## stream population sex treatment mass sl time latency peck event
## 1 Maracas HU F C 0.313 20.497 83 142 1
## 2 Maracas HU F C 0.144 17.051 14:56 204 127 1
## 3 Maracas HU F C 0.182 18.035 11:39 59 106 1
## 4 Maracas HU F C 0.206 18.325 12:10 58 98 1
## 5 Maracas HU F C 0.306 19.653 13:02 44 96 1
## 6 Maracas HU F C 0.224 17.606 13:25 44 92 1
## date initial comments foraging_time foraging_time_sec latency_surv
## 1 12/14/25 NV <NA> NA 83
## 2 7/27/25 SMT 14M 56S NA 204
## 3 4/12/25 SMT 11M 39S NA 59
## 4 4/12/25 SMT 12M 10S NA 58
## 5 11/23/25 NV 13M 2S NA 44
## 6 11/23/25 NV 13M 25S NA 44table(data$sex,data$treatment, data$population, data$stream)## , , = LU, = Maracas
##
##
## C H
## F 17 17
## M 15 15
##
## , , = HU, = Maracas
##
##
## C H
## F 19 19
## M 20 20
##
## , , = LU, = Tacarigua
##
##
## C H
## F 19 18
## M 19 19
##
## , , = HU, = Tacarigua
##
##
## C H
## F 16 15
## M 20 20Checking which model is better to run
model_pecks_nb <- glmmTMB(peck ~ mass + (population + treatment + sex)^3 +
(1|stream),
family = nbinom2,
data = data)
model_pecks_p <- glmmTMB(peck ~ mass + (population + treatment + sex)^3 +
(1|stream),
family = genpois,
data = data)
AIC(model_pecks_p, model_pecks_nb)## df AIC
## model_pecks_p 11 3003.402
## model_pecks_nb 11 2954.505The negative binomial model is better to run.
Peck and Latency Models
model_peck <- glmmTMB(peck ~ log(mass) + (population + treatment + sex)^3 +
(1|stream),
family = nbinom2,
data = data)
# coxme coef are opposite from glmmtmb models - looks at mortality rather than survival
#i.e., (-) is actually an increase in latency ;;; slower
model_latency <- coxph(latency_surv ~ log(mass) + (population + treatment + sex)^3 +
cluster(stream),
data = data)
summary(model_peck)## Family: nbinom2 ( log )
## Formula:
## peck ~ log(mass) + (population + treatment + sex)^3 + (1 | stream)
## Data: data
##
## AIC BIC logLik -2*log(L) df.resid
## 2956.6 2996.9 -1467.3 2934.6 277
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## stream (Intercept) 0.02829 0.1682
## Number of obs: 288, groups: stream, 2
##
## Dispersion parameter for nbinom2 family (): 1.31
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.17977 0.38923 8.169 3.1e-16 ***
## log(mass) -0.47144 0.17612 -2.677 0.00743 **
## populationHU 0.31064 0.21221 1.464 0.14323
## treatmentH -0.01280 0.20991 -0.061 0.95136
## sexM -0.11927 0.25508 -0.468 0.64010
## populationHU:treatmentH -0.44208 0.30033 -1.472 0.14103
## populationHU:sexM -0.63025 0.29755 -2.118 0.03416 *
## treatmentH:sexM -0.07112 0.29941 -0.238 0.81224
## populationHU:treatmentH:sexM 0.33319 0.42061 0.792 0.42827
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(model_latency)## Call:
## coxph(formula = latency_surv ~ log(mass) + population + treatment +
## sex + population:treatment + population:sex + treatment:sex +
## population:treatment:sex, data = data, cluster = stream)
##
## n= 288, number of events= 271
##
## coef exp(coef) se(coef) robust se z
## log(mass) -0.11846 0.88829 0.19897 0.11251 -1.053
## populationHU -0.01330 0.98679 0.24662 0.47471 -0.028
## treatmentH -0.21103 0.80975 0.24596 0.12258 -1.722
## sexM 0.49708 1.64392 0.31841 0.23420 2.122
## populationHU:treatmentH -0.27224 0.76167 0.35292 0.20748 -1.312
## populationHU:sexM -0.32604 0.72177 0.33988 0.19056 -1.711
## treatmentH:sexM -0.48444 0.61604 0.34766 0.02468 -19.625
## populationHU:treatmentH:sexM 0.36981 1.44746 0.48830 0.24520 1.508
## Pr(>|z|)
## log(mass) 0.2924
## populationHU 0.9776
## treatmentH 0.0852 .
## sexM 0.0338 *
## populationHU:treatmentH 0.1895
## populationHU:sexM 0.0871 .
## treatmentH:sexM <2e-16 ***
## populationHU:treatmentH:sexM 0.1315
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## log(mass) 0.8883 1.1258 0.7125 1.1074
## populationHU 0.9868 1.0134 0.3892 2.5021
## treatmentH 0.8098 1.2349 0.6368 1.0297
## sexM 1.6439 0.6083 1.0388 2.6015
## populationHU:treatmentH 0.7617 1.3129 0.5072 1.1439
## populationHU:sexM 0.7218 1.3855 0.4968 1.0486
## treatmentH:sexM 0.6160 1.6233 0.5869 0.6466
## populationHU:treatmentH:sexM 1.4475 0.6909 0.8951 2.3406
##
## Concordance= 0.603 (se = 0.031 )
## Likelihood ratio test= 26.78 on 8 df, p=8e-04
## Wald test = 1.11 on 8 df, p=1
## Score (logrank) test = 29.8 on 8 df, p=2e-04, Robust = 2 p=1
##
## (Note: the likelihood ratio and score tests assume independence of
## observations within a cluster, the Wald and robust score tests do not).Predicted values for plotting
Peck
pred_peck <- ggpredict(model_peck, terms= c("treatment", "population","sex"), back_transform = TRUE) %>%
mutate(streampair = paste(x, group, sep = "_"),
sexpair = paste(group, facet, sep = "_")) %>%
drop_na()## You are calculating adjusted predictions on the population-level (i.e.
## `type = "fixed"`) for a *generalized* linear mixed model.
## This may produce biased estimates due to Jensen's inequality. Consider
## setting `bias_correction = TRUE` to correct for this bias.
## See also the documentation of the `bias_correction` argument.streampairlevels <-c("C_LU","H_LU", "C_HU", "H_HU")
pred_peck$streampair = factor(pred_peck$streampair, levels=streampairlevels)Latency
pred_latency <- ggpredict(model_latency, terms= c("treatment", "population","sex"), back_transform = TRUE) %>%
mutate(streampair = paste(x, group, sep = "_"),
sexpair = paste(group, facet, sep = "_")) %>%
drop_na()
pred_latency$streampair = factor(pred_latency$streampair, levels=streampairlevels)PLOTS
Figure 1. Peck
Figure 2. Latency
Test Plots
Test plots, by sex - note: the latency plots do not have the same values as the original plots, they may be back transformed. The plots are also not aligned by Y axis, so you cannot compare them visually side by side.
