Skip to content

What stresses out California ground squirrels?​

Code Block 1

#LINEAR MIXED EFFECTS MODEL: CORT DISTURBANCE DATA
#install packages for running statistical models, assessing repeatability, and making graphs
#install.packages('glmmTMB')

#load libraries running statistical models, assessing repeatability, and making graphs
library(glmmTMB)
library(ggplot2)

# Compare "stress levels" at the two study sites
CORT_full=read.csv("/data/groups/UB_Workshop/Day_1/Session_3/CORT data for full analysis minus disturbance_5-28-24.csv")
m_full_area <- glmmTMB(lnCort ~ day + hour + site*mass_lhf + sex*mass_lhf + min_age_master*mass_lhf + (1|uid) + (1|area) + (1|year), family = gaussian, data = CORT_full)
summary(m_full_area)

# Create a violin plot to compare "stress levels" at the two study sites
dp <- ggplot(CORT_full, aes(x=disturbance, y=lnCort, fill=disturbance)) + 
  geom_violin(trim=FALSE)+
  geom_boxplot(width=0.1, fill="white") +
  labs(x="Disturbance at study site", y = "Ln fecal glucocorticoid metabolites (ng/g feces)")
# Change color by groups, choose a palette: "Dark2" or "Reds" or "Greens" or "Oranges"
dp + scale_fill_brewer(palette="Dark2") + theme_minimal()
dp + scale_fill_brewer(palette="Reds") + theme_classic()

# Create scatterplot to compare effects of sex, mass and body size on stress
p <- ggplot(CORT_full, aes(x = mass_lhf, y = lnCort, color=sex, shape = sex)) + 
  geom_point(colour = "grey32", alpha=I(0.3)) + theme_bw() + 
  stat_smooth(method = "lm") +
theme(axis.text=element_text(size=20)) + 
theme(axis.title=element_text(size=15)) + 
theme(legend.text=element_text(size=20)) + 
theme(legend.title=element_text(size=20)) +
  labs(x = "Mass/left hind foot (g/mm)", y = "Ln fecal glucocorticoid metabolites (ng/g feces)", legend = "Sex")
p+scale_color_brewer(palette="Dark2")

Code Block 2

###ASSESS EFFECTS OF INDIVIDUAL IDENTITY AND YEAR ON STRESS
library(glmmTMB)
library(rptR)
CORT_full=read.csv("/data/groups/UB_Workshop/Day_1/Session_3/CORT data for full analysis minus disturbance_5-28-24.csv")

## Random significance tests
### COMPARING MODELS WITH AND WITHOUT RANDOM EFFECTS
m_all_uid <- glmmTMB(lnCort ~ day + hour + site*mass_lhf + sex*mass_lhf + min_age_master*mass_lhf + (1|uid), family = gaussian, data = CORT_full)
summary(m_all_uid)

m_all_area <- glmmTMB(lnCort ~ day + hour + site*mass_lhf + sex*mass_lhf + min_age_master*mass_lhf + (1|area), family = gaussian, data = CORT_full)
summary(m_all_area)

m_all_year <- glmmTMB(lnCort ~ day + hour + site*mass_lhf + sex*mass_lhf + min_age_master*mass_lhf + (1|year), family = gaussian, data = CORT_full)
summary(m_all_year)

m_all_no_random <- glmmTMB(lnCort ~ day + hour + site*mass_lhf + sex*mass_lhf + min_age_master*mass_lhf, family = gaussian, data = CORT_full)
summary(m_all_no_random)

## USE THIS TO LOOK AT SIGNIFICANCE OF RANDOM EFFECTS
anova(m_all_no_random, m_all_uid)
anova(m_all_no_random, m_all_area)
anova(m_all_no_random, m_all_year)

## USE THIS TO LOOK AT SIGNIFICANCE OF RANDOM EFFECTS
anova(m_all_no_random, m_all_uid)
anova(m_all_no_random, m_all_area)
anova(m_all_no_random, m_all_year)
#All random variables are statistically significant

##### Repeatbility of CORT 

#Repeatability of squirrel identity
rep1 <- rpt(lnCort ~ (1|uid), grname = "uid", data = CORT_full, datatype = "Gaussian", nboot = 10, npermut = 10)
rep2 <- rpt(lnCort ~ day + hour + site*mass_lhf + sex*mass_lhf + min_age_master*mass_lhf + (1|uid), grname = "uid", data = CORT_full, datatype = "Gaussian", nboot = 1000, npermut = 1000)

#Repeatability of trapping area
rep3 <- rpt(lnCort ~ (1|area), grname = "area", data = CORT_full, datatype = "Gaussian", nboot = 10, npermut = 10)
rep4 <- rpt(lnCort ~ day + hour + site*mass_lhf + sex*mass_lhf + min_age_master*mass_lhf + (1|area), grname = "area", data = CORT_full, datatype = "Gaussian", nboot = 1000, npermut = 1000)

#GRAPH repeatability estimates from small vs large number of simulations
#repeatability estimates of stress for squirrel identity from only 10 simulations 
print(rep1)
plot(rep1, cex.main = 1)
#repeatability estimates of stress for squirrel identity from 1000 simulations 
print(rep2)
plot(rep2, cex.main = 1)

#repeatability estimates of stress for area of capture from only 10 simulations 
print(rep3)
plot(rep3, cex.main = 1)
#repeatability estimates of stress for area of capture from only 1000 simulations 
print(rep4)
plot(rep4, cex.main = 1)

Code Block 3

####EFFECTS OF HUMAN & DOG ACTIVITY ON CORT 2018-2023
#load libraries running statistical models, assessing repeatability, and making graphs
library(glmmTMB)
library(ggplot2)

#LET'S LOAD THE REDUCED DATA SET FOR HUMAN & DOG ACTIVITY from cluster
#CORT_crow18_23=read.csv("CORT data for full analysis_AND_disturb_day_scaled_AB_corrected_18_23_6-19-24.csv", h=T)

CORT_crow18_23=read.csv("/data/groups/UB_Workshop/Day_1/Session_3/CORT data for full analysis_AND_disturb_day_scaled_AB_corrected_18_23_6-19-24.csv")

#LET'S RUN THE STATISTICAL MODEL:
m_crow18_23 <- glmmTMB(lnCort ~ hour + day.s + mass_lhf + sex*mass_lhf + min_age_master*mass_lhf + disttotalscale + (1|uid) + (1|year), family = gaussian, data = CORT_crow18_23)
summary(m_crow18_23)

#LET'S PLOT THE DATA TO SEE HOW THESE EFFECTS VARY BY YEAR:
ggplot(CORT_crow18_23, aes(x = disttotalscale, y = lnCort)) +  geom_point(colour = "grey32", alpha=I(0.3)) + stat_smooth(method = "lm") + 
theme_bw() + labs(x = "Human and dog activity (scaled)", y = "Ln fecal glucocorticoid metabolites(ng/g feces)") + facet_wrap(~year, strip.position = "bottom")