A short tutorial to analyse stress in wild rodents.
rodent, Maxomys surifer, cortisol, cmr (capture, mark, recapture)
Summary (according to google AI): In wild rodent populations, cortisol (or corticosterone in rodents) levels vary based on several factors, including age, sex, social status (Schradin 2008), and environmental conditions (Beery 2015). Levels can fluctuate seasonally, with breeding males often exhibiting lower corticosterone levels during the breeding season, while other classes show declines in corticosterone during periods of low food abundance. Furthermore, urban populations may experience different stress hormone profiles compared to rural populations (Łopucki 2019).
Factors Influencing Cortisol Levels:
Age: Cortisol levels can vary with age, with juveniles potentially showing different levels than adults.
Sex: Studies have shown differences in cortisol levels between males and females.
Social Status: In some species, subordinate animals may experience higher cortisol levels due to factors like social instability or limited social contact.
Seasonal Changes: Corticosterone levels in some rodent species can fluctuate seasonally, with breeding females and non-breeding males and females showing a decline in levels during the non-breeding season.
Environmental Conditions: Studies have shown that urban and rural populations may exhibit different stress hormone profiles, with urban populations sometimes having narrower ranges and lower median values.
Stressors: Various stressors, such as capture, handling, and exposure to predators or competitors, can trigger acute increases in cortisol levels.
Reproductive Status: Pregnancy and other reproductive events can also influence cortisol levels.
Animals were trapped using living trap cages. Their feces were collected and prepared for estimating the level of cortisol.
load("data_for_cortisol.Rdata")
data <- data |>
dplyr::mutate(meanOD = (duplicate1 + duplicate2) / 2) |>
dplyr::mutate(netOD = meanOD - dplyr::first(meanOD)) |>
dplyr::mutate(BB0 = (netOD/dplyr::last(netOD))*100) |>
dplyr::mutate(meanOD = round(meanOD, 3),
netOD = round(netOD, 3),
BB0 = round(BB0, 2))Concentrations of cortisol are in picogram / mL
DT::datatable(data, class = 'cell-border stripe',
caption = 'Table 1: "Data calibration for cortisol"',
options = list(columnDefs =
list(list(className = 'dt-center',
targets = "_all"))))library(drc)
model <- drc::drm(netOD ~ concentration, data = data, fct = LL.4())# Tidy results
tidy_model <- broom::tidy(model)
knitr::kable(tidy_model, caption = "Concentration – Optical Measure Model Summary")| term | curve | estimate | std.error | statistic | p.value |
|---|---|---|---|---|---|
| b | (Intercept) | 1.3333070 | 0.2475773 | 5.385417 | 0.0057478 |
| c | (Intercept) | 0.1017402 | 0.0844127 | 1.205272 | 0.2945235 |
| d | (Intercept) | 0.9328456 | 0.0259820 | 35.903500 | 0.0000036 |
| e | (Intercept) | 702.2786678 | 146.5903989 | 4.790755 | 0.0087068 |
#Model residuals (exact ones used by summary(model))
res <- residuals(model)
#Fitted values
fit <- fitted(model)
#Observed values
obs <- model$y # guaranteed aligned
#Counts
n <- length(res)
p <- length(coef(model))
df <- n - p
# Residual standard error (same formula as summary(model))
RSE <- sqrt(sum(res^2) / df)
#RMSE (divide by n)
RMSE <- sqrt(mean(res^2))
#MSE
MSE <- mean(res^2)
#MAE
MAE <- mean(abs(res))
#Pseudo-R² (variance explained on response scale)
R2 <- 1 - (sum(res^2) / sum((obs - mean(obs))^2))
#Model statistics from glance()
gl <- broom::glance(model)
diagnostics <- dplyr::tibble(
n = n,
p = p,
df_residual = df,
AIC = gl$AIC,
BIC = gl$BIC,
logLik = gl$logLik,
RSE = RSE,
RMSE = RMSE,
MSE = MSE,
MAE = MAE
)
knitr::kable(diagnostics, caption = "Model Concentration – Optical Measure Diagnostic Metrics")| n | p | df_residual | AIC | BIC | logLik | RSE | RMSE | MSE | MAE |
|---|---|---|---|---|---|---|---|---|---|
| 8 | 4 | 4 | -26.89533 | -26.49812 | 18.44766 | 0.0341046 | 0.0241156 | 0.0005816 | 0.0171099 |
plot(model)
response <- response |>
dplyr::mutate(meanOD = (duplicate1 + duplicate2) / 2) |>
dplyr::mutate(netOD = meanOD - dplyr::first(data$meanOD)) |>
dplyr::mutate(BB0 = (netOD/dplyr::last(data$netOD))*100) |>
dplyr::mutate(meanOD = round(meanOD, 3),
netOD = round(netOD, 3),
BB0 = round(BB0, 2))
DT::datatable(response, class = 'cell-border stripe',
caption = 'Table 2: "Data lecture of samples"',
options = list(columnDefs =
list(list(className = 'dt-center',
targets = "_all"))))DOSEx<-ED(model,response$netOD,type="absolute",display=F)
plot(model)
points(y=response$netOD,x=DOSEx[,1],col="blue",pch=19,cex=2)
#arrows(DOSEx[,1],response$netOD,DOSEx[,1]+DOSEx[,2]*1.96,response$netOD,length=0.1,angle=90,lwd=3,col="blue")
#arrows(DOSEx[,1],response$netOD,DOSEx[,1]-DOSEx[,2]*1.96,response$netOD,length=0.1,angle=90,lwd=3,col="blue")data_merge <- cbind(response,DOSEx) |>
dplyr::left_join(data_rodent,
by = c("id_indiv","trappingDate")) |>
dplyr::mutate(body_condition = bodyWeight / length_mm) |>
dplyr::mutate(body_condition = round(body_condition, 2),
length_mm = round(length_mm, 2),
Estimate = round(Estimate,0),
`Std. Error` = round(`Std. Error`,0))
DT::datatable(data_merge, class = 'cell-border stripe',
caption = 'Table 3: "Data on individual rodents with estimated concentration of cortisol"',
options = list(columnDefs =
list(list(className = 'dt-center',
targets = "_all"))))model_bd <- lm(bodyWeight ~ length_mm,data_merge)
plot(effects::allEffects(model_bd),
col = 2,
ylab = "Body weight",
#ylim = c(-1, 1),
type = "response")
bmi_res <- as.data.frame(resid(model_bd ))
names(bmi_res) = c("bmi_res")
data_merge <- cbind(data_merge,bmi_res)data_merge <- data_merge |>
dplyr::mutate(trappingDate = as.Date(trappingDate)) |>
dplyr::group_by(id_indiv) |>
dplyr::mutate(day = trappingDate - dplyr::first(trappingDate) +1) |>
dplyr::mutate(day = as.integer(day))
DT::datatable(data_merge, class = 'cell-border stripe',
caption = 'Table 3: "Data on individual rodents with estimated concentration of cortisol"',
options = list(columnDefs =
list(list(className = 'dt-center',
targets = "_all"))))random factor = individuals (id_indiv) as several measures on the same individuals explanatory variables = number of days after the first recapture (day 1) sex BMI (residuals)
library(lme4)
model <- glmer(Estimate ~ sex + day + (1|id_indiv), data = data_merge)library(sjPlot)
library(sjstats)
library(sjlabelled)
library(sjmisc)
#plot_model(model, vline.color = "red",show.values = TRUE, value.offset = .3) +
# theme_sjplot2() +
# scale_color_sjplot("simply")
tab_model(model, show.se = TRUE, show.std = TRUE, show.stat = TRUE,
show.aic = TRUE)| Estimate | ||||||||
|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | std. Beta | standardized std. Error | CI | standardized CI | Statistic | p |
| (Intercept) | 472.43 | 63.58 | 0.38 | 0.19 | 342.20 – 602.67 | -0.01 – 0.76 | 7.43 | <0.001 |
| sex [M] | -205.26 | 63.63 | -0.96 | 0.30 | -335.61 – -74.91 | -1.57 – -0.35 | -3.23 | 0.003 |
| day | 29.63 | 11.29 | 0.36 | 0.14 | 6.51 – 52.74 | 0.08 – 0.65 | 2.62 | 0.014 |
| Random Effects | ||||||||
| σ2 | 26142.76 | |||||||
| τ00 id_indiv | 1028.42 | |||||||
| ICC | 0.04 | |||||||
| N id_indiv | 10 | |||||||
| Observations | 33 | |||||||
| Marginal R2 / Conditional R2 | 0.421 / 0.443 | |||||||
| AIC | 412.019 | |||||||
library(effects)
effects = allEffects(model)
plot(effects,
#col = 3,
ylab = "Cortisol level",
type = "response")
library(ggplot2)
ggplot(data_merge, aes(x = day, y = Estimate)) +
geom_point(aes(colour = sex)) +
geom_line(color = "blue") +
facet_wrap( ~ id_indiv) +
theme_classic()
library(leaflet)
library(leafpop)
map <- leaflet() %>%
addTiles() %>%
addCircleMarkers(data = data_merge,
~long, ~lat,
label = ~ as.character(id_indiv),
popup = popupTable(data_merge |>
dplyr::select(sex,Estimate,day)),
fillColor = "blue", weight = 0.2,
group = "Maxomys surifer",
stroke = FALSE, fillOpacity = 0.4) %>%
addProviderTiles(providers$Esri.WorldStreetMap,group ="WSM") %>%
addTiles(group = "OSM") %>%
addTiles(urlTemplate = "https://mts1.google.com/vt/lyrs=s&hl=en&src=app&x={x}&y={y}&z={z}&s=G",
attribution = 'Google',group ="Google") %>%
addLayersControl(baseGroups = c("WSM","OSM","Google"),
overlayGroups = c("Maxomys surifer"),
options = layersControlOptions(collapsed = FALSE)) %>%
addScaleBar(position = "bottomright",
scaleBarOptions(maxWidth = 100, metric = TRUE, imperial = FALSE,
updateWhenIdle = TRUE))
maplongevity_df <- readxl::read_excel("/Users/serge/Documents/DATA/data_muka/muka_rodents/Rodent_MUKA_database_20250731.xlsx",
sheet="trapRoudup") |>
dplyr::filter(id_indiv %in% unique(data_merge$id_indiv)) |>
dplyr::group_by(id_indiv) |>
dplyr::summarise(
first_seen = min(trappingDate),
last_seen = max(trappingDate),
longevity_days = as.numeric(difftime(last_seen, first_seen, units = "days"))) |>
dplyr::arrange(desc(longevity_days)) |>
dplyr::left_join(data_merge |>
dplyr::filter(day == "1") |>
dplyr::select(id_indiv,id_grid_line,sex,bodyWeight, Estimate),
by="id_indiv") |>
dplyr::rename(longevity = longevity_days,
grid = id_grid_line,
weight = bodyWeight) |>
dplyr::mutate(first_seen= as.Date(first_seen),
last_seen= as.Date(last_seen),
Estimate= round(Estimate, 1)) |>
dplyr::left_join(read.csv("data_area_range.csv"), by = "id_indiv") |>
dplyr::mutate(area_sq_meters = round(area_sq_meters, 0)) |>
dplyr::rename(area = area_sq_meters,
cortisol = Estimate) |>
dplyr::relocate(cortisol, .after = last_col())
DT::datatable(longevity_df, class = 'cell-border stripe',
caption = 'Table 4: "Individual longevity"',
options = list(columnDefs =
list(list(className = 'dt-center',
targets = "_all"))))# Kruskal Wallis Test One Way Anova by Ranks
KW_test <- kruskal.test(Estimate ~ id_grid_line, data = data_merge |>
dplyr::filter(day == "1"))
# Tidy results
tidy_model2 <- broom::tidy(KW_test)
knitr::kable(tidy_model2, caption = "Effect of location: Kruskal Wallis Test One Way Anova by Ranks")| statistic | p.value | parameter | method |
|---|---|---|---|
| 3.345455 | 0.1877344 | 2 | Kruskal-Wallis rank sum test |
data <- data_merge |>
dplyr::filter(day == "1")
SR_test <- cor.test(data$bmi_res, data$Estimate, method = "spearman")
# Tidy results
tidy_model3 <- broom::tidy(SR_test)
knitr::kable(tidy_model3, caption = "Effect of of body mass index (BMI): Spearman rank correlation")| estimate | statistic | p.value | method | alternative |
|---|---|---|---|---|
| -0.0060606 | 166 | 1 | Spearman’s rank correlation rho | two.sided |
no statistical effect of home range on the level of cortisol
SR_test2 <- cor.test(longevity_df$area, longevity_df$cortisol, method = "spearman")
# Tidy results
tidy_model4 <- broom::tidy(SR_test2)
knitr::kable(tidy_model4, caption = "Effect of of home range: Spearman rank correlation")| estimate | statistic | p.value | method | alternative |
|---|---|---|---|---|
| -0.1428571 | 64 | 0.7825397 | Spearman’s rank correlation rho | two.sided |