############################################################
# Risk Factors for Mild and Severe Subclinical Mastitis
# Supplementary Materials R code
# Author: Muhammad Aktaruzzaman, M. Ariful Islam, Solama Akter Shanta
#Mst. Tahomina Akter, Muhammad Tofazzal Hossain, Md. Nazmul Islam, 
#Md. Shaffiul Alam, A K M Anisur Rahman

# Date: 19/11/2025
############################################################

# -----------------------------
# 1. Load libraries
# -----------------------------
library(readr)
library(dplyr)
library(tidyr)
library(stringr)
library(ggplot2)
library(purrr)
library(broom)
library(scales)
library(brms)
library(DescTools)
library(car)

# -----------------------------
# 2. Load dataset
# -----------------------------
scc_dat <- read_csv("SCC_data.csv")
# Check dimensions and column names
dim(scc_dat)
names(scc_dat)

# -----------------------------
# 3. Herd size summary
# -----------------------------
herd_size_by_system <- scc_dat %>%
  group_by(Management_system, Herd) %>%
  summarise(Herd_size = n(), .groups = "drop")

herd_size_summary <- herd_size_by_system %>%
  group_by(Management_system) %>%
  summarise(
    n_herds = n(),
    Mean_herd_size = mean(Herd_size),
    Median_herd_size = median(Herd_size),
    Min_herd_size = min(Herd_size),
    Max_herd_size = max(Herd_size),
    Range = paste0(Min_herd_size, "–", Max_herd_size)
  )

print(herd_size_summary)

# -----------------------------
# 4. SCC classification
# -----------------------------
scc_dat <- scc_dat %>%
  mutate(
    SCC_class = case_when(
      SCC <= 147292 ~ "Healthy",
      SCC <= 799051 ~ "Mild",
      TRUE ~ "Severe"
    ),
    SCC_class = factor(SCC_class, levels = c("Healthy", "Mild", "Severe"))
  )

# -----------------------------
# 5. Prevalence by Management System (for Figure 2)
# -----------------------------
# Compute prevalence
prev_core_ms <- scc_dat %>%
  group_by(Management_system) %>%
  mutate(n_total_ms = n()) %>%
  group_by(Management_system, SCC_class, n_total_ms) %>%
  summarise(n = n(), .groups = "drop")

# Overall SCM prevalence (Mild + Severe)
overall_ms <- prev_core_ms %>%
  filter(SCC_class %in% c("Mild", "Severe")) %>%
  group_by(Management_system, n_total_ms) %>%
  summarise(n = sum(n), Prevalence = sum(n)/n_total_ms,
            SCC_class = "Overall", .groups = "drop")

# Combine
prev_core_ms2 <- bind_rows(prev_core_ms, overall_ms) %>%
  mutate(SCC_class = factor(SCC_class, levels = c("Healthy", "Mild", "Severe", "Overall")))

# Compute 95% CI
prev_data_ms <- prev_core_ms2 %>%
  mutate(
    CI = map2(n, n_total_ms, ~ tidy(binom.test(.x, .y))),
    CI_lower = map_dbl(CI, ~ .x$conf.low),
    CI_upper = map_dbl(CI, ~ .x$conf.high),
    Prev_pct = n / n_total_ms * 100,
    CI_lower_pct = CI_lower * 100,
    CI_upper_pct = CI_upper * 100
  ) %>%
  select(-CI)

# -----------------------------
# Figure 2: Prevalence by Management System
# -----------------------------
my_cols <- c(
  "Healthy" = "#1b9e77",
  "Mild"    = "#d95f02",
  "Severe"  = "#7570b3",
  "Overall" = "#e7298a"
)

prev_faceted <- ggplot(prev_data_ms, aes(x = SCC_class, y = Prev_pct, fill = SCC_class)) +
  geom_col(width = 0.65, color = "black") +
  geom_segment(aes(x = SCC_class, xend = SCC_class, y = CI_lower_pct, yend = CI_upper_pct),
               size = 0.9, color = "black") +
  geom_segment(aes(x = as.numeric(SCC_class) - 0.1, xend = as.numeric(SCC_class) + 0.1,
                   y = CI_lower_pct, yend = CI_lower_pct), size = 0.9, color = "black") +
  geom_segment(aes(x = as.numeric(SCC_class) - 0.1, xend = as.numeric(SCC_class) + 0.1,
                   y = CI_upper_pct, yend = CI_upper_pct), size = 0.9, color = "black") +
  geom_text(aes(label = sprintf("%.1f%%", Prev_pct), y = CI_upper_pct + 6), size = 4, fontface = "bold") +
  geom_text(aes(label = sprintf("95%% CI: %.1f–%.1f", CI_lower_pct, CI_upper_pct),
                y = CI_upper_pct + 12), size = 3.8, fontface = "bold", color = "gray20") +
  facet_grid(~ Management_system, switch = "x") +
  scale_fill_manual(values = my_cols) +
  labs(x = "Subclinical Mastitis Status", y = "Prevalence (%)", fill = "SCM status") +
  theme_bw() +
  theme(strip.placement = "outside", strip.background = element_blank(),
        strip.text.x = element_text(size = 14, face = "bold"),
        axis.title = element_text(size = 15, face = "bold"),
        axis.text = element_text(size = 13, face = "bold"),
        legend.position = "top", legend.title = element_text(size = 13, face = "bold"),
        legend.text = element_text(size = 12))

ggsave("Figure_2.tiff", prev_faceted, device = "tiff", dpi = 700, compression = "lzw",
       width = 11, height = 6, units = "in")

print(prev_faceted)

# -----------------------------
# 6. Descriptive statistics
# -----------------------------
num_vars <- c("Age", "Parity", "Milk_yield", "Lact_stage","SCC")

summarize_mean_sd <- function(x) {
  sprintf("%.1f ± %.1f", mean(x, na.rm = TRUE), sd(x, na.rm = TRUE))
}

summary_table <- scc_dat %>%
  group_by(SCC_class) %>%
  summarise(across(all_of(num_vars), summarize_mean_sd)) %>%
  ungroup()

overall_summary <- scc_dat %>%
  summarise(across(all_of(num_vars), summarize_mean_sd)) %>%
  mutate(SCC_class = "Overall")

final_summary <- bind_rows(overall_summary, summary_table) %>%
  relocate(SCC_class)

print(final_summary)

# -----------------------------
# 7. Univariable and multivariable models
# -----------------------------
# Prepare dataset with lowest-prevalence reference for categorical variables
set_lowest_ref <- function(var) {
  var <- factor(var)
  tbl <- table(var)
  ref_level <- names(tbl)[which.min(tbl)]
  factor(var, levels = c(ref_level, setdiff(levels(var), ref_level)))
}

scc_dat_prepped <- scc_dat %>%
  mutate(
    Parity_qtile = set_lowest_ref(factor(ntile(Parity, 5))),
    Lact_stage_qtile = set_lowest_ref(factor(ntile(Lact_stage, 5))),
    Hist_cli_mast = set_lowest_ref(factor(Hist_cli_mast)),
    Pregcat = set_lowest_ref(factor(Pregcat)),
    SCC_class = factor(SCC_class, levels = c("Healthy", "Mild", "Severe")),
    Milk_yield = as.numeric(scale(Milk_yield))
  )

# -----------------------------
# Multivariable Bayesian multinomial regression
# -----------------------------
multi_formula <- as.formula(
  "SCC_class ~ Milk_yield + Hist_cli_mast + Parity_qtile + Lact_stage_qtile + (1 | Herd)"
)

multi_fit <- brm(
  formula = multi_formula,
  data = scc_dat_prepped,
  family = categorical(),
  chains = 3,
  iter = 33000,
  warmup = 5000,
  cores = 2,
  seed = 123,
  control = list(adapt_delta = 0.95, max_treedepth = 15)
)

# -----------------------------
# Extract fixed effects posterior
# -----------------------------
fe_post <- fixef(multi_fit, summary = FALSE)

posterior_summary <- map_dfr(colnames(fe_post), function(par) {
  post_par <- fe_post[, par]
  outcome <- strsplit(par, "_")[[1]][2]  
  tibble(
    Predictor = strsplit(par, "_")[[1]][1],
    Outcome = outcome,
    Parameter = par,
    Posterior_Mean = mean(post_par),
    CI_2.5 = quantile(post_par, 0.025),
    CI_97.5 = quantile(post_par, 0.975),
    Prob_Positive = mean(post_par > 0)
  )
})

# -----------------------------
# Convert log-odds to RRR
# -----------------------------
posterior_summary_rrr <- posterior_summary %>%
  mutate(
    RRR = exp(Posterior_Mean),
    RRR_lower = exp(CI_2.5),
    RRR_upper = exp(CI_97.5),
    RRR_summary = paste0(round(RRR, 3), " (", round(RRR_lower, 3), "–", round(RRR_upper, 3), ")")
  )

write.csv(posterior_summary_rrr, "posterior_table_rrr.csv", row.names = FALSE)
print(posterior_summary_rrr)

# -----------------------------
# Figure 3: Significant RRR plot
# -----------------------------
# Parse predictors and remove intercepts
df_sig <- posterior_summary_rrr %>%
  filter(!str_detect(Parameter, "Intercept")) %>%
  mutate(
    Predictor_group = case_when(
      str_detect(Predictor, "Milk_yield") ~ "Milk yield",
      str_detect(Predictor, "Parity_qtile") ~ "Parity",
      str_detect(Predictor, "Lact_stage_qtile") ~ "Lactation stage",
      str_detect(Predictor, "Hist_cli_mast") ~ "History of clinical mastitis",
      TRUE ~ Predictor
    ),
    AxisLabel = Predictor_group
  ) %>%
  filter(!(RRR_lower < 1 & RRR_upper > 1))

p_sig <- ggplot(df_sig, aes(x = AxisLabel, y = RRR, color = Outcome, shape = Outcome)) +
  geom_point(size = 3, position = position_dodge(0.7)) +
  geom_errorbar(aes(ymin = RRR_lower, ymax = RRR_upper),
                width = 0.25, position = position_dodge(0.7)) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "gray50") +
  coord_flip() +
  theme_minimal(base_size = 14) +
  labs(y = "Relative Risk Ratio (RRR, 95% CI)", x = NULL) +
  scale_color_manual(values = c("Mild" = "#1f77b4", "Severe" = "#d62728")) +
  scale_shape_manual(values = c("Mild" = 16, "Severe" = 17)) +
  scale_y_log10(breaks = c(0.25, 0.5, 1, 2, 4, 8),
                labels = c("0.25", "0.5", "1", "2", "4", "8"))

ggsave("Figure_3.tiff", p_sig, width = 9, height = 7, dpi = 900, compression = "lzw")
print(p_sig)