Preliminary analyses
2025-08-03
DATA ENTRY
data <- read.csv("../data/forage.csv") %>%
mutate(foraging_time = lubridate::ms(foraging_time),
foraging_time_sec = as.numeric(foraging_time),
latency_surv = Surv(latency),
stream = relevel(factor(stream), ref = "LU"),
sex = factor(sex),
treatment = factor(treatment))
head(data)## id stream population sex treatment mass sl foraging_time latency
## 1 MHU.F.1 HU Maracas F H 0.202 18.651 12M 20S 110
## 2 MHU.F.10 HU Maracas F C 0.198 18.635 11M 27S 109
## 3 MHU.F.11 HU Maracas F C 0.182 18.035 11M 39S 59
## 4 MHU.F.12 HU Maracas F C 0.158 16.192 11M 49S 24
## 5 MHU.F.13 HU Maracas F C 0.244 19.827 11M 59S 67
## 6 MHU.F.14 HU Maracas F C 0.206 18.325 12M 10S 58
## pecks rate forage_date foraging_time_sec latency_surv
## 1 53 17.66667 12-Apr 740 110
## 2 67 22.33333 12-Apr 687 109
## 3 106 35.33333 12-Apr 699 59
## 4 47 15.66667 12-Apr 709 24
## 5 80 26.66667 12-Apr 719 67
## 6 98 32.66667 12-Apr 730 58PECKS
model_pecks_nb <- glmmTMB(pecks ~ (stream + treatment + sex)^3 +
(1|population),
family = nbinom2,
data = data)
model_pecks_p <- glmmTMB(pecks ~ (stream + treatment + sex)^3 +
(1|population),
family = genpois,
data = data)
AIC(model_pecks_p, model_pecks_nb)## df AIC
## model_pecks_p 10 1151.713
## model_pecks_nb 10 1133.617The negative binomial model is better, so we will work with it by removing non significant interactions
model_pecks <- glmmTMB(pecks ~ (sex + stream + treatment) +
(1|population),
family = nbinom2,
data = data)
summary(model_pecks)## Family: nbinom2 ( log )
## Formula: pecks ~ (sex + stream + treatment) + (1 | population)
## Data: data
##
## AIC BIC logLik deviance df.resid
## 1126.1 1142.3 -557.1 1114.1 104
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## population (Intercept) 0.005965 0.07724
## Number of obs: 110, groups: population, 2
##
## Dispersion parameter for nbinom2 family (): 2.21
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.248792 0.148132 28.682 <2e-16 ***
## sexM -0.112320 0.130567 -0.860 0.390
## streamHU 0.008874 0.130831 0.068 0.946
## treatmentH 0.003045 0.131441 0.023 0.982
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1LATENCY
model_latency <- coxme(latency_surv ~ (stream + treatment + sex)^2 +
(1|population),
data = data)
model_latency## Cox mixed-effects model fit by maximum likelihood
## Data: data
## events, n = 110, 110
## Iterations= 5 31
## NULL Integrated Fitted
## Log-likelihood -410.3228 -400.8392 -400.8372
##
## Chisq df p AIC BIC
## Integrated loglik 18.97 7 0.0082907 4.97 -13.94
## Penalized loglik 18.97 6 0.0042235 6.96 -9.25
##
## Model: latency_surv ~ (stream + treatment + sex)^2 + (1 | population)
## Fixed coefficients
## coef exp(coef) se(coef) z p
## streamHU -0.2223944 0.8005995 0.3411493 -0.65 0.51000
## treatmentH -0.2917965 0.7469205 0.3427208 -0.85 0.39000
## sexM 1.1769982 3.2446198 0.3529072 3.34 0.00085
## streamHU:treatmentH 0.5647093 1.7589364 0.3920254 1.44 0.15000
## streamHU:sexM -1.0983232 0.3334297 0.3981083 -2.76 0.00580
## treatmentH:sexM -0.6470157 0.5236060 0.3923138 -1.65 0.09900
##
## Random effects
## Group Variable Std Dev Variance
## population Intercept 8.923498e-03 7.962881e-05model_latency2 <- glmmTMB(log(latency) ~ (stream + treatment + sex)^2 +
(1|population),
data = data)
summary(model_latency2)## Family: gaussian ( identity )
## Formula:
## log(latency) ~ (stream + treatment + sex)^2 + (1 | population)
## Data: data
##
## AIC BIC logLik deviance df.resid
## 351.6 375.9 -166.8 333.6 101
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## population (Intercept) 7.629e-10 2.762e-05
## Residual 1.215e+00 1.102e+00
## Number of obs: 110, groups: population, 2
##
## Dispersion estimate for gaussian family (sigma^2): 1.21
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.82442 0.29455 12.984 < 2e-16 ***
## streamHU 0.09362 0.37163 0.252 0.80111
## treatmentH 0.18186 0.37163 0.489 0.62460
## sexM -1.08216 0.37163 -2.912 0.00359 **
## streamHU:treatmentH -0.31415 0.42079 -0.747 0.45532
## streamHU:sexM 0.90719 0.42079 2.156 0.03109 *
## treatmentH:sexM 0.76694 0.42079 1.823 0.06836 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1pred_lat <- ggpredict(model_latency2, terms= c("stream", "treatment", "sex"),
back_transform = F)## Model has log transformed response. Predictions are on transformed
## scale.ggplot(pred_lat, aes(y = exp(predicted), x = group, col = facet, group = facet)) +
geom_point(position = position_dodge(width = 0.5)) +
geom_errorbar(
aes(ymin = exp(predicted - std.error), ymax = exp(predicted + std.error)),
width = 0,
position = position_dodge(width = 0.5)
) +
geom_line(position = position_dodge(width = 0.5)) +
labs(y = "Peck latency (sec)", x = "Treatment", col = "Sex") +
facet_grid(~x) +
theme_bw()