Pathogen Exposure Analysis

Questions

Question 1: Does method of fungal exposure affect TTD?
Question 2: Does pathogen exposure type influence space use?
Question 3: Does pathogen exposure type influence fragmentation?

Load libraries and data

#Stegodyphus-metarhizium analysis
####Load Libraries####
library(conflicted) # Conflicted functions
library(tidyverse) # Data carpentry
library(ggfortify) # plot survival curve
library(survival) # survival analysis
library(survminer) # plot survival
library(coxme) # proportional hazards model
library(emmeans) # Estimated marginal means
library(glmmTMB) # Another generalized linear mixed effects model package
library(DHARMa) # checking regression assumptions
library(ggstatsplot) # Chi Squared
library(rstatix) # Chi squared pairwise comparisons

# Fixed conflicted functions
conflicted::conflicts_prefer(dplyr::filter)

[conflicted] Removing existing preference.
[conflicted] Will prefer dplyr::filter over any other package.

conflicted::conflicts_prefer(dplyr::select)

[conflicted] Removing existing preference.
[conflicted] Will prefer dplyr::select over any other package.

####Load csv files from R project folder####
survival <- read.csv("data/ttd_2.csv")
location <- read.csv("data/exposure_location.csv") # 4=out and 5=dead
pilot <- read.csv("data/survival_pilot_2020.csv")

Body Condition

Calculate body condition

#Calculate SMI from Parthasarathy et al. 2018
SMI <- function(mass, prosoma){
  cw_bar <- mean(prosoma, na.rm=T)
  m <- lm(prosoma ~ mass)
  slope <- summary(m)[[4]][2]
  cor <- cor(mass, prosoma)
  
  smi <- mass*(cw_bar/prosoma)^(slope/cor)
  return(smi)
}

#Combine SMI vector with body size dataframe
smi <- SMI(survival$initial_mass, survival$initial_prosoma_width);survival <- cbind(survival, SMI=smi)

##Q0: Pilot study

Question 1: Does fungal exposure affect Time To Death (TTD)?

Model: Cox Proportional hazards model

Dependent Variable: TTD

Independent Variable:
(1) Treatment

Random Effect:
(1) Colony ID
(2) Spider age (juvenile or adult)

Run Cox model and polt survival curve

#### Q0: Analyze individual mortality with a proportional hazards mixed model####
# Cox model containing mixed (random and fixed) effects
# [ttdfit0]Average ttd when S. dumicola exposed to M. robertsii
ttdfit0 <- coxme(Surv(TTD) ~ treatment + (1|colony.ID) + (1|age), data = pilot)

print(ttdfit0)

Cox mixed-effects model fit by maximum likelihood
  Data: pilot
  events, n = 49, 49 (1 observation deleted due to missingness)
  Iterations= 6 34 
                    NULL Integrated  Fitted
Log-likelihood -144.5657  -129.7258 -128.01

                  Chisq   df          p   AIC   BIC
Integrated loglik 29.68 3.00 1.6115e-06 23.68 18.00
 Penalized loglik 33.11 1.92 5.5979e-08 29.28 25.65

Model:  Surv(TTD) ~ treatment + (1 | colony.ID) + (1 | age) 
Fixed coefficients
               coef exp(coef)  se(coef)    z     p
treatmentE 2.066782  7.899364 0.4224541 4.89 1e-06

Random effects
 Group     Variable  Std Dev      Variance    
 colony.ID Intercept 0.0197955955 0.0003918656
 age       Intercept 0.9273314196 0.8599435617

# Fit to model without random effects for plot
fit <- survfit(Surv(TTD) ~ treatment, data = pilot)

# Plot pilot survival curve
ggsurvplot(fit, data = pilot, size = 1, palette = c("#2E9FDF", "darkolivegreen3"), 
           risk.table = FALSE, 
           conf.int = TRUE,
           pval = TRUE,legend.labs = c("Control", "Exposed"),
           ggtheme = theme(panel.grid.major = element_blank(),
                           panel.grid.minor = element_blank(),
                           panel.background = element_blank(),
                           axis.line = element_line(colour ="black"),
           text = element_text(size=15)))

ggsave("figures/pilot_surv_curve.png")

Saving 7.29 x 4.51 in image

##Q1: Survival

Question 1: Does method of fungal exposure affect Time To Death (TTD)?

Model: Cox Proportional hazards model

Dependent Variable: TTD

Independent Variable:
(1) Treatment (2) SMI

Random Effect:
(1) Colony ID
(2) Block

Run Cox model

Raw data

#### Q1: Analyze individual mortality with a proportional hazards mixed model####
# Cox model containing mixed (random and fixed) effects
# [ttdfit1]Average ttd when S. dumicola exposed to M. robertsii
ttdfit1 <- coxme(Surv(ttd_2, censor) ~ treatment + SMI + (1|colony_id) + (1|block) + (1|source_retreat), data = survival)

print(ttdfit1)

Cox mixed-effects model fit by maximum likelihood
  Data: survival
  events, n = 74, 140
  Iterations= 9 58 
                    NULL Integrated    Fitted
Log-likelihood -341.6881  -291.3376 -275.7682

                   Chisq    df p   AIC   BIC
Integrated loglik 100.70  7.00 0  86.7 70.57
 Penalized loglik 131.84 14.92 0 102.0 67.62

Model:  Surv(ttd_2, censor) ~ treatment + SMI + (1 | colony_id) + (1 |      block) + (1 | source_retreat) 
Fixed coefficients
                         coef   exp(coef)    se(coef)     z       p
treatmentcricket  3.773037240  43.5120208 1.088314428  3.47 5.3e-04
treatmentpaper    2.982937042  19.7457256 1.096389045  2.72 6.5e-03
treatmentspider   5.063165038 158.0900864 1.089922271  4.65 3.4e-06
SMI              -0.002033213   0.9979689 0.005245269 -0.39 7.0e-01

Random effects
 Group          Variable  Std Dev      Variance    
 colony_id      Intercept 0.6835888057 0.4672936553
 block          Intercept 0.0199795322 0.0003991817
 source_retreat Intercept 0.0199861315 0.0003994455

# Pairwise comparisons for original data
pairwise_original <- pairwise_survdiff(Surv(ttd_2, censor) ~ treatment, data = survival)
print(pairwise_original)

    Pairwise comparisons using Log-Rank test 

data:  survival and treatment 

        control cricket paper  
cricket 7.0e-08 -       -      
paper   4.8e-05 0.031   -      
spider  1.4e-15 1.7e-05 8.3e-10

P value adjustment method: BH 

No primary cases

#### Re-analyze individual immunity with a proportional hazards mixed model with primary cases removed####
# Remove primary cases from curve
survival_no_primary <- filter(survival, primary_case=="n", died_before_meta == "n")

# Cox model containing mixed (random and fixed) effects
# [ttdfit3]Average ttd when Stegos exposed to metarhizium without primary cases
ttdfit2 <- coxme(Surv(ttd_2, censor) ~ treatment + SMI + (1|colony_id) + (1|block) + (1|source_retreat), data = survival_no_primary)
print(ttdfit2)

Cox mixed-effects model fit by maximum likelihood
  Data: survival_no_primary
  events, n = 66, 132
  Iterations= 9 58 
                    NULL Integrated    Fitted
Log-likelihood -302.3586  -257.5749 -241.6802

                   Chisq    df          p   AIC   BIC
Integrated loglik  89.57  7.00 1.1102e-16 75.57 60.24
 Penalized loglik 121.36 15.12 0.0000e+00 91.12 58.02

Model:  Surv(ttd_2, censor) ~ treatment + SMI + (1 | colony_id) + (1 |      block) + (1 | source_retreat) 
Fixed coefficients
                          coef  exp(coef)    se(coef)     z       p
treatmentcricket  3.8187701525  45.548157 1.097666555  3.48 5.0e-04
treatmentpaper    2.9400010268  18.915866 1.107455271  2.65 7.9e-03
treatmentspider   5.0157183003 150.764392 1.103419060  4.55 5.5e-06
SMI              -0.0009024413   0.999098 0.005537331 -0.16 8.7e-01

Random effects
 Group          Variable  Std Dev      Variance    
 colony_id      Intercept 0.7314800247 0.5350630266
 block          Intercept 0.0199960245 0.0003998410
 source_retreat Intercept 0.0199965029 0.0003998601

# Pairwise comparisons with no primaries
pairwise_no_primary <- pairwise_survdiff(Surv(ttd_2, censor) ~ treatment, data = survival_no_primary)
print(pairwise_no_primary)

    Pairwise comparisons using Log-Rank test 

data:  survival_no_primary and treatment 

        control cricket paper
cricket 7e-08   -       -    
paper   0.00010 0.01857 -    
spider  2e-14   0.00026 5e-09

P value adjustment method: BH 

Plot survival

Raw data

#### Plot Survival curve ####
# Filter out spiders that died before metarhizium exposure
surv <- survival %>% 
  filter(died_before_meta == "n")

# Fit to model without random effects for plot
fit <- survfit(Surv(ttd_2, censor) ~ treatment, data = surv)

# Plot with primary cases
plot_primary <- ggsurvplot(fit, data = surv, size = 1, palette = c("#2E9FDF", "orange4", "bisque3", "darkgreen"), 
                           risk.table = FALSE, 
                           conf.int = TRUE,
                           pval = FALSE,
                           legend.title = "Treatment",
                           legend.labs = c("Control", "Cricket", "Paper", "Spider"),
                           xlim = c(0, 30),
                           xlab = "Days",
                           ggtheme = theme(panel.grid.major = element_blank(),
                                           panel.grid.minor = element_blank(),
                                           panel.background = element_blank(),
                                           axis.line = element_line(colour ="black"),
                                           text = element_text(size=23)))

# Change x-axis scale to intervals of 5 days and add significance letters
plot_primary$plot <- plot_primary$plot + 
  scale_x_continuous(breaks = seq(0, 30, by = 5)) +
  annotate("text",
           x = 29, y = 0.98,
           label = "a", size = 7) +
  annotate("text",
           x = 29, y = 0.56,
           label = "b", size = 7) +
  annotate("text",
           x = 29, y = 0.3,
           label = "c", size = 7) +
  annotate("text",
           x = 29, y = 0.04,
           label = "d", size = 7)

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.

# View plot
plot_primary

# Save plot to figures folder
ggsave("figures/surv_curve.png", width = 9, height = 6, units = "in")

Polt with no primary cases

# Fit to model without primary cases and random effects for plot
fit2 <- survfit(Surv(ttd_2, censor) ~ treatment, data = survival_no_primary)

# Plot without primary cases
plot_no_primary <- ggsurvplot(fit2, data = survival_no_primary, size = 1, palette = c("#2E9FDF", "orange4", "bisque3", "darkgreen"), 
                              risk.table = FALSE, 
                              conf.int = TRUE,
                              pval = FALSE,
                              legend.title = "Treatment",
                              legend.labs = c("Control", "Cricket", "Paper", "Spider"),
                              xlim = c(0, 30),
                              xlab = "Days",
                              ggtheme = theme(panel.grid.major = element_blank(),
                                              panel.grid.minor = element_blank(),
                                              panel.background = element_blank(),
                                              axis.line = element_line(colour ="black"),
                                              text = element_text(size=18)))

# Change x-axis scale to intervals of 7 days
plot_no_primary$plot <- plot_no_primary$plot + 
  scale_x_continuous(breaks = seq(0, 30, by = 5)) +
  annotate("text",
           x = 29, y = 0.98,
           label = "a", size = 6) +
  annotate("text",
           x = 29, y = 0.56,
           label = "b", size = 6) +
  annotate("text",
           x = 29, y = 0.3,
           label = "c", size = 6) +
  annotate("text",
           x = 29, y = 0.04,
           label = "d", size = 6)

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.

# View plot
print(plot_no_primary)

# Save plot to figures folder
ggsave("figures/surv_curve_no_primary.png", width = 9, height = 6, units = "in")

##Q2: Space use

Question 2: Does pathogen exposure affect the proportion of time individuals spend inside vs outside their nest?

Model: GLMM

Dependent Variable: Proportion of time individuals spent inside nest

Independent Variable:
(1) Treatment

Random Effect:
(1) Colony ID
(2) Block

Wrangle TTD data

#### Q2: Pathogen exposure space use analysis ####
# Apply the pivot_longer() function to reshape the data
location_long <- location %>% 
  pivot_longer(cols = contains("location"), # Select columns containing "location"
               names_to = "spider_id", # Name for the new column with spider IDs 
               values_to = "location") # Name for the new column with location

# Filter rows by unique 'colony_id'
filtered_survival <- survival %>%
  distinct(colony_id, .keep_all = TRUE)

# Merge the datasets by 'colony_id'
merged_location <- location_long %>%
  left_join(dplyr::select(filtered_survival, colony_id, treatment), by = "colony_id") %>%
  filter(location != "5") %>%
  unite(subject_id, spider_id, colony_id, sep = "_", remove = FALSE) %>% 
  filter(subject_id != "yellow_location_13")

# Merge location and ttd data
merged_location_ttd <- merged_location %>%
  left_join(survival[,c("subject_id","ttd_2", "primary_case", "source_retreat", "censor")], by = c("subject_id"="subject_id"))

# Calculate the count of individuals at each location
location_in <- merged_location_ttd %>%
  filter(location < 4, date >= "2023-06-08") %>% 
  group_by(date, subject_id) %>%
  mutate(Count = n()) %>% 
  ungroup() %>%
  group_by(subject_id) %>%
  mutate(sumCount = sum(Count)/ttd_2) %>% 
  ungroup() %>% 
  filter(!duplicated(subject_id))

# Calculate summary statistics for space use by treatment
summary_space_use_stats <- location_in %>%
  group_by(treatment.y) %>%
  summarise(mean = mean(sumCount),
            median = median(sumCount),
            sd = sd(sumCount),
            min = min(sumCount),
            max = max(sumCount))

# Print the summary statistics
print(summary_space_use_stats)

Run GLMM

# Distribution of space use by treatment
ggplot(data = location_in) +
  geom_histogram(aes(x = sumCount, fill = treatment.y), bins = 10, color = "black") +
  scale_fill_manual(values = c("#2E9FDF", "orange4", "bisque3", "darkgreen")) +
  facet_wrap(~ treatment.y) 

# Fit a GLMM
space_glmm <- glmmTMB(sumCount ~ treatment.y + (1 | block) + (1 | colony_id) + (1|source_retreat), data = location_in)

# Display model summary
summary(space_glmm)

 Family: gaussian  ( identity )
Formula:          sumCount ~ treatment.y + (1 | block) + (1 | colony_id) + (1 |      source_retreat)
Data: location_in

     AIC      BIC   logLik deviance df.resid 
  -354.6   -331.1    185.3   -370.6      131 

Random effects:

Conditional model:
 Groups         Name        Variance  Std.Dev.
 block          (Intercept) 1.437e-05 0.003791
 colony_id      (Intercept) 1.132e-03 0.033648
 source_retreat (Intercept) 1.969e-04 0.014034
 Residual                   3.271e-03 0.057191
Number of obs: 139, groups:  block, 7; colony_id, 28; source_retreat, 4

Dispersion estimate for gaussian family (sigma^2): 0.00327 

Conditional model:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)         0.99095    0.01761   56.27  < 2e-16 ***
treatment.ycricket -0.03434    0.02281   -1.51    0.132    
treatment.ypaper   -0.03045    0.02286   -1.33    0.183    
treatment.yspider  -0.13275    0.02259   -5.88 4.21e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Inside: Check your model / regression assumptions
plot(simulateResiduals(space_glmm))

testDispersion(simulateResiduals(space_glmm))

    DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated

data:  simulationOutput
dispersion = 1.0213, p-value = 0.8
alternative hypothesis: two.sided

# Test significance
# Use type III sums of squares
Anova(space_glmm, type=3)

Analysis of Deviance Table (Type III Wald chisquare tests)

Response: sumCount
               Chisq Df Pr(>Chisq)    
(Intercept) 3165.958  1  < 2.2e-16 ***
treatment.y   38.904  3  1.818e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Posthoc Tukey
# Extract least squares means
pairs(emmeans(space_glmm, specs="treatment.y", type="response"))

 contrast          estimate     SE  df t.ratio p.value
 control - cricket  0.03434 0.0228 131   1.506  0.4370
 control - paper    0.03045 0.0229 131   1.332  0.5443
 control - spider   0.13275 0.0226 131   5.876  <.0001
 cricket - paper   -0.00389 0.0227 131  -0.172  0.9982
 cricket - spider   0.09841 0.0228 131   4.315  0.0002
 paper - spider     0.10230 0.0229 131   4.475  0.0001

P value adjustment: tukey method for comparing a family of 4 estimates 

Plot space use

All spiders

# Plot with raw data
raw_proportions <- ggplot(location_in, aes(x = factor(treatment.y, levels = c("control", "paper", "cricket", "spider")), y = sumCount, color = treatment.y, fill = treatment.y,
                                           shape = as.character(censor))) +
  geom_jitter(size = 2.75, width = 0.35, alpha = 0.5) +
  ylim(0.4,1.05) +
  stat_summary(fun = mean, size = 1.15, shape = 21, color = "black") + 
  stat_summary(aes(x = factor(treatment.y, levels = c("control", "paper", "cricket", "spider")), 
                    y = sumCount, color = treatment.y, color = treatment.y),
                geom="errorbar", width = 0.2, size = 1.25, inherit.aes = F) +
  scale_color_manual(values = c("#2E9FDF", "orange4", "bisque3", "darkgreen")) +
  scale_fill_manual(values = c("#2E9FDF", "orange4", "bisque3", "darkgreen")) +
  guides(color = "none") +
  labs(x = "Treatment", y = "Proportion of days inside nest", color = "", fill = "Treatments") +
  scale_x_discrete(labels=c("Control", "Paper", "Cricket", "Spider")) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour ="black"),
        text = element_text(size=23), legend.position = "none")

Warning: Duplicated aesthetics after name standardisation: colour

raw_proportions +
  annotate("text", x = unique(location_in$treatment.y), 
           y = rep(max(location_in$sumCount) + 0.04,length(unique(location_in$treatment.y))),
           label = c("a", "a", "a", "b"), color = "black", size = 7)

No summary function supplied, defaulting to `mean_se()`

Warning: Removed 4 rows containing missing values (`geom_segment()`).

ggsave("figures/space_use.png", width = 9, height = 6, units = "in")

No summary function supplied, defaulting to `mean_se()`

Warning: Removed 4 rows containing missing values (`geom_segment()`).

Only dead spiders

# Plot with raw data on only spiders that died
# Filter out living spiders
location_in_dead <- location_in %>% 
  filter(ttd_2 < 28)

# Plot
raw_proportions_dead <- ggplot(location_in_dead, aes(x = factor(treatment.y, 
                                                                levels = c("control", "paper", "cricket", "spider")), 
                                                     y = sumCount, color = treatment.y, fill = treatment.y)) +
  geom_jitter(width = 0.25, alpha = 0.5) +
  stat_summary(fun = mean, size = 1, shape = 21, color = "black") + 
  stat_summary(geom="errorbar", width = 0.2, size = 1.25) +
  scale_color_manual(values = c("#2E9FDF", "orange4", "bisque3", "darkgreen")) +
  scale_fill_manual(values = c("#2E9FDF", "orange4", "bisque3", "darkgreen")) +
  guides(color = "none") +
  labs(x = "Treatment", y = "Proportion of days inside nest", color = "", fill = "Treatments") +
  scale_x_discrete(labels=c("Control", "Paper", "Cricket", "Spider")) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour ="black"),
        text = element_text(size=18), legend.position = "none")

raw_proportions_dead +
  annotate("text", x = unique(location_in$treatment.y), 
           y = rep(max(location_in$sumCount) + 0.04,length(unique(location_in$treatment.y))),
           label = c("a", "a", "a", "b"), color = "black", size = 6)

No summary function supplied, defaulting to `mean_se()`

Warning: Removed 4 rows containing missing values (`geom_segment()`).

ggsave("figures/space_use_dead.png", width = 9, height = 6, units = "in")

No summary function supplied, defaulting to `mean_se()`

Warning: Removed 4 rows containing missing values (`geom_segment()`).

##Q3: Polydomy

Question 3: Does pathogen exposure mode influence the construction of new nests?

Model: Logistic regression/Chi squared

Dependent Variable: Polydomy y/n

Independent Variable:
(1) treatment

Wrangle data

#### Polydomy Analysis ####
# Filter out just last day of data
location_last <- location %>% 
  filter(date == "2023-07-06") %>% 
  dplyr::select(poly, treatment)

# Convert variables to factors
location_last$poly <- as.factor(location_last$poly)
location_last$treatment <- as.factor(location_last$treatment)

Run chi squared

# Chi squared comparing between treatments
contingency_table <- table(location_last$poly, location_last$treatment)

chi <- chisq.test(contingency_table)

Warning in chisq.test(contingency_table) :
  Chi-squared approximation may be incorrect

print(chi)

    Pearson's Chi-squared test

data:  contingency_table
X-squared = 3.1111, df = 3, p-value = 0.3748

ggplot(location_last, aes(treatment, fill = poly)) +
  geom_bar(colour="black", position = "fill") +
  scale_fill_manual(values = c("lightblue", "darkblue"),
                    labels = c("Monodomous", "Polydomous"),
                    name = "Fragmentation") +
  ylab("Proportion") +
  xlab("Treatment") +
  scale_x_discrete(labels=c("Control", "Cricket", "Paper", "Spider")) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour ="black"),
        text = element_text(size=15),
        legend.position = "top")


ggsave("figures/polydomy_treatment.png", width = 5, height = 5, units = "in")