Read data
data <- read.csv("~/Desktop/forage.csv") %>%
mutate(foraging_time_of_day = lubridate::hm(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_of_day = lubridate::hm(time)`.
## Caused by warning in `.parse_hms()`:
## ! Some strings failed to parse
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
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_of_day foraging_time_sec latency_surv
## 1 12/14/25 NV <NA> NA 83
## 2 7/27/25 SMT 14H 56M 0S NA 204
## 3 4/12/25 SMT 11H 39M 0S NA 59
## 4 4/12/25 SMT 12H 10M 0S NA 58
## 5 11/23/25 NV 13H 2M 0S NA 44
## 6 11/23/25 NV 13H 25M 0S NA 44
table(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 20
f_data <- data %>%
filter(sex == 'F')
m_data <- data %>%
filter(sex == 'M')
#checking which model is better to run by looking for lower AIC
model_pecks_nb <- glmmTMB(peck ~ mass + (population + treatment)^2 +
(1|stream),
family = nbinom2,
data = data)
model_pecks_p <- glmmTMB(peck ~ mass + (population + treatment)^2 +
(1|stream),
family = genpois,
data = data)
AIC(model_pecks_p, model_pecks_nb)
## df AIC
## model_pecks_p 7 3017.472
## model_pecks_nb 7 2955.959
#negative binomial model is better to run
ANALYSES
Male Peck Count
m_model_pecks <- glmmTMB(peck ~ log(mass) + population + treatment +
(1|stream),
family = nbinom2,
data = m_data)
summary(m_model_pecks)
## Family: nbinom2 ( log )
## Formula: peck ~ log(mass) + population + treatment + (1 | stream)
## Data: m_data
##
## AIC BIC logLik -2*log(L) df.resid
## 1521.3 1539.3 -754.6 1509.3 142
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## stream (Intercept) 0.04096 0.2024
## Number of obs: 148, groups: stream, 2
##
## Dispersion parameter for nbinom2 family (): 1.42
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.68896 0.99180 4.728 2.27e-06 ***
## log(mass) 0.09605 0.34627 0.277 0.78148
## populationHU -0.40334 0.14276 -2.825 0.00472 **
## treatmentH -0.17602 0.14170 -1.242 0.21417
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m_pred_peck <- predict_response(m_model_pecks, terms= c("treatment", "population"), margin = "empirical", type = "link")
Female Peck Count
f_model_pecks <- glmmTMB(peck ~ log(mass) + population + treatment +
(1|stream),
family = nbinom2,
data = f_data)
summary(f_model_pecks)
## Family: nbinom2 ( log )
## Formula: peck ~ log(mass) + population + treatment + (1 | stream)
## Data: f_data
##
## AIC BIC logLik -2*log(L) df.resid
## 1435.6 1453.3 -711.8 1423.6 134
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## stream (Intercept) 1.63e-08 0.0001277
## Number of obs: 140, groups: stream, 2
##
## Dispersion parameter for nbinom2 family (): 1.22
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.1874 0.4001 7.966 1.64e-15 ***
## log(mass) -0.5291 0.2004 -2.640 0.00829 **
## populationHU 0.1031 0.1576 0.654 0.51294
## treatmentH -0.2272 0.1555 -1.461 0.14409
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
f_pred_peck <- predict_response(f_model_pecks, terms= c("treatment", "population"), margin = "empirical", type = "link")
Male Foraging Latency
# coxme coef are opposite from glmmtmb models - looks at mortality rather than survival. doctors study MORTALITY not survival -- how fast you die
#i.e., (-) is actually an increase in latency ;;; slower
#latency assessment for males and females
m_coxmodel_latency <- coxme(latency_surv ~ log(mass) + population + treatment +
(1|stream),
data = m_data)
m_coxmodel_latency
## Cox mixed-effects model fit by maximum likelihood
## Data: m_data
## events, n = 141, 148
## Iterations= 2 10
## NULL Integrated Fitted
## Log-likelihood -586.4804 -576.5247 -576.4992
##
## Chisq df p AIC BIC
## Integrated loglik 19.91 4.00 0.00051996 11.91 0.12
## Penalized loglik 19.96 3.02 0.00017690 13.92 5.02
##
## Model: latency_surv ~ log(mass) + population + treatment + (1 | stream)
## Fixed coefficients
## coef exp(coef) se(coef) z p
## log(mass) 0.8383360 2.3125157 0.3507707 2.39 0.01700
## populationHU -0.3410949 0.7109914 0.1720869 -1.98 0.04700
## treatmentH -0.6071461 0.5449038 0.1724338 -3.52 0.00043
##
## Random effects
## Group Variable Std Dev Variance
## stream Intercept 0.019990474 0.000399619
Female Foraging Latency
f_coxmodel_latency <- coxme(latency_surv ~ log(mass) + population + treatment +
(1|stream),
data = f_data)
f_coxmodel_latency
## Cox mixed-effects model fit by maximum likelihood
## Data: f_data
## events, n = 130, 140
## Iterations= 5 22
## NULL Integrated Fitted
## Log-likelihood -540.1159 -534.2502 -534.2473
##
## Chisq df p AIC BIC
## Integrated loglik 11.73 4 0.0194650 3.73 -7.74
## Penalized loglik 11.74 3 0.0083754 5.73 -2.89
##
## Model: latency_surv ~ log(mass) + population + treatment + (1 | stream)
## Fixed coefficients
## coef exp(coef) se(coef) z p
## log(mass) -0.5945128 0.5518314 0.2318170 -2.56 0.010
## populationHU -0.1185763 0.8881840 0.1780441 -0.67 0.510
## treatmentH -0.3637813 0.6950432 0.1776084 -2.05 0.041
##
## Random effects
## Group Variable Std Dev Variance
## stream Intercept 8.902317e-03 7.925125e-05
FIGURES
Fig 1a. Male Peck
fig1a <- ggplot(m_pred_peck, aes(x, exp(predicted), group = group, color=group)) +
geom_point(position = position_dodge(width=0.5), size = 2) +
geom_line(linewidth=0.75, position = position_dodge(width = 0.5)) +
geom_errorbar(aes(ymin=(exp(predicted - std.error)), ymax=(exp(predicted + std.error))),
linewidth = 0.75, width = 0.25, position = position_dodge(width = 0.5)) +
scale_color_manual(values=c("chartreuse4", "gray25")) +
labs(x = NULL,
y = "Peck Count",
title = "Peck count, males",
col = "Population") +
ylim(35,100) +
theme_bw() +
theme(legend.position="none")
fig1a
Fig 1b. Female Peck
fig1b <- ggplot(f_pred_peck, aes(x, exp(predicted), group = group, color=group)) +
geom_point(position = position_dodge(width=0.5), size = 2) +
geom_line(linewidth=0.75, position = position_dodge(width = 0.5)) +
geom_errorbar(aes(ymin=exp(predicted - std.error), ymax=exp(predicted + std.error)), linewidth=0.75, width = 0.1, position = position_dodge(width = 0.5)) +
scale_color_manual(values=c("chartreuse4", "gray25")) +
labs(x = NULL,
y = "Peck Count",
title = "Peck count, females",
col = "Population") +
ylim(35,100) +
theme_bw() +
theme(legend.position="none",
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank())
fig1b
Fig 2a. Male Latency
fig2a <- m_data %>%
ggplot(aes(y = latency, x = treatment, group = population, color = population)) +
scale_color_manual(values=c("chartreuse4", "gray25")) +
stat_summary(fun = mean, geom = "point", size = 1.5, position = position_dodge(width = 0.5)) +
stat_summary(fun = mean, geom = "line",linewidth= 0.75, position = position_dodge(width = 0.5)) +
stat_summary(fun.data = "mean_se", geom = "errorbar", width = 0.1, position = position_dodge(width = 0.5)) +
theme_bw() +
labs(x = NULL,
y = "Latency (s)",
title = "Latency (s), males",
col = "Population") +
theme_bw() +
theme(legend.position="none")
fig2a
Fig 2b. Female Latency
fig2b <- f_data %>%
ggplot(aes(y = latency, x = treatment, group = population, color = population)) +
scale_color_manual(values=c("chartreuse4", "gray25")) +
stat_summary(fun = mean, geom = "point", size = 1.5, position = position_dodge(width = 0.5)) +
stat_summary(fun = mean, geom = "line",linewidth= 0.75,position = position_dodge(width = 0.5)) +
stat_summary(fun.data = "mean_se", geom = "errorbar", width = 0.1,position = position_dodge(width = 0.5)) +
theme_bw() +
labs(x = NULL,
y = "Latency (s)",
title = "Latency (s), females",
col = "Population") +
theme_bw() +
theme(legend.position="none",
axis.title.y = element_blank())
fig2b
Supplementary
# Female Survival Curve by Population
f_surv_fit <- surv_fit(Surv(latency, event = event) ~ population,
data = f_data,
na.action = "na.fail")
# Female Survival Curve Plot
ggsurvplot(f_surv_fit, conf.int = TRUE, fun = 'event',
palette = c('black', 'grey'),
cumcensor = F,
xlab = 'Time (s)',
legend.title = "Population",
legend.labs = c("HU","LU"))
# Male Survival Curve by Population
m_surv_fit <- surv_fit(Surv(latency, event = event) ~ population,
data = m_data,
na.action = "na.fail")
# Male Survival Curve Plot
ggsurvplot(m_surv_fit, conf.int = TRUE, fun = 'event',
palette = c('black', 'grey'),
cumcensor = F,
xlab = 'Time (s)',
legend.title = "Population",
legend.labs = c("HU","LU"))
