############################################################
# Statistical analysis of Bd prevalence in frogs
# 
# This script analyses Bd infection prevalence using
# binomial generalized linear models (GLMs) with a logit
# link function. Effects of host traits (adult ecology
# and reproductive mode) on Bd presence are assessed.
#
# Software:
# R version 4.2.2 (R Core Team, 2025)
#
# Packages:
# dplyr, emmeans, multcomp
############################################################

# Load required packages
library(dplyr)
library(emmeans)
library(multcomp)

# Load data
Ecol <- read.csv("Amphibian_Ecology.csv")

# Convert categorical predictors to factors
Ecol$Ecology <- as.factor(Ecol$Ecology)
Ecol$Reproductive_Mode <- as.factor(Ecol$Reproductive_Mode)

############################################################
# Binomial GLMs: Bd presence ~ host ecology
############################################################

# Model 1: Effect of adult ecology on Bd presence
m_ecology <- glm(
  Bd_presence ~ Ecology,
  family = binomial(link = "logit"),
  data = Ecol
)

# Null model for likelihood ratio test
m_ecology_null <- glm(
  Bd_presence ~ 1,
  family = binomial(link = "logit"),
  data = Ecol
)

# Likelihood ratio test
anova(m_ecology_null, m_ecology, test = "Chisq")

# Post-hoc Tukey pairwise comparisons
ecology_emm <- emmeans(m_ecology, ~ Ecology)
pairs(ecology_emm, adjust = "tukey")

# Predicted probabilities
summary(ecology_emm, type = "response")

############################################################

# Model 2: Effect of reproductive mode on Bd presence
m_repro <- glm(
  Bd_presence ~ Reproductive_Mode,
  family = binomial(link = "logit"),
  data = Ecol
)

summary(m_repro)

# Null model for likelihood ratio test
m_repro_null <- glm(
  Bd_presence ~ 1,
  family = binomial(link = "logit"),
  data = Ecol
)

# Likelihood ratio test
anova(m_repro_null, m_repro, test = "Chisq")

# Post-hoc Tukey pairwise comparisons
repro_emm <- emmeans(m_repro, ~ Reproductive_Mode)
pairs(repro_emm, adjust = "tukey")

# Predicted probabilities
summary(repro_emm, type = "response")