# ================================================================
# District-Level Trends in Veterinary Antimicrobial Prescribing
# in a Resource-Limited Setting: Evidence from Mongar, Bhutan (2023-2026)
#
# Analysis script - Paper 1
# Data: DVH Mongar + 16 RNR-Extension Centres, VIS clinical records
#       Jan 2023 - Jun 2026 (merged and cleaned in Excel prior to
#       import: 4 annual exports combined, columns standardized,
#       Medicine Class capitalization/spelling corrected)
# Author: Narayan Pokhrel, Cheda Cheda
# ================================================================


# ---- Step 1: Load required packages ----
library(tidyverse)
library(readxl)

# ---- Step 2: Load the cleaned AMU dataset ----
setwd("D:/AMR Manuscript")
df <- read_excel("AMU data.xlsx")
dim(df)
colnames(df)

# ---- Step 3: Filter to Antimicrobials only ----
amu <- df %>% filter(`Medicine Class` == "Antimicrobials")
nrow(amu)
table(amu$`Record Year`)

# ---- Step 4: Year-wise summary with annualized 2026 ----
yearly_summary <- amu %>%
  count(`Record Year`, name = "n_prescriptions") %>%
  mutate(months_covered = ifelse(`Record Year` == 2026, 6, 12),
         annualized_n = round(n_prescriptions / months_covered * 12, 1))
yearly_summary

# ---- Step 5: Plot yearly trend [Figure 1] ----
fig1 <- ggplot(yearly_summary, aes(x = `Record Year`)) +
  geom_col(aes(y = n_prescriptions), fill = "steelblue") +
  geom_point(aes(y = annualized_n), color = "red", size = 3) +
  geom_line(aes(y = annualized_n, group = 1), color = "red", linetype = "dashed") +
  labs(title = "Antimicrobial Prescriptions by Year, Mongar (2023-2026)",
       subtitle = "Bars = raw count | Red dashed line = annualized estimate",
       x = "Year", y = "Number of Prescriptions") +
  theme_minimal()
fig1
ggsave("Figure1_Yearly_Trend.png", fig1, width = 8, height = 5, dpi = 300)

# ---- Step 6: Yearly-level trend test ----
trend_model <- lm(n_prescriptions ~ `Record Year`, data = yearly_summary)
summary(trend_model)

# ---- Step 7: Monthly-level trend test ----
monthly_summary <- amu %>%
  count(`Record Year`, `Record Month`, name = "n_prescriptions") %>%
  mutate(month_index = row_number())
monthly_summary
monthly_trend_model <- lm(n_prescriptions ~ month_index, data = monthly_summary)
summary(monthly_trend_model)

# ---- Step 8: Prescriptions by year and species ----
species_yearly <- amu %>%
  count(`Record Year`, `Species Name`, name = "n_prescriptions")
species_yearly

# ---- Step 9: Restrict to main 3 species, annualize 2026 ----
species_yearly_adj <- species_yearly %>%
  filter(`Species Name` %in% c("Bovine", "Canine", "Feline")) %>%
  mutate(months_covered = ifelse(`Record Year` == 2026, 6, 12),
         annualized_n = round(n_prescriptions / months_covered * 12, 1))
species_yearly_adj

# ---- Step 9b: Plot species trend [Figure 2] ----
fig2 <- ggplot(species_yearly_adj, aes(x = `Record Year`, y = annualized_n, color = `Species Name`)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  labs(title = "Antimicrobial Prescriptions by Species and Year, Mongar (2023-2026)",
       subtitle = "Bovine, Canine, and Feline only (2026 annualized)",
       x = "Year", y = "Annualized Number of Prescriptions", color = "Species") +
  theme_minimal()
fig2
ggsave("Figure2_Species_Trend.png", fig2, width = 8, height = 5, dpi = 300)

# ---- Step 10: Bovine monthly trend test ----
bovine_monthly <- amu %>%
  filter(`Species Name` == "Bovine") %>%
  count(`Record Year`, `Record Month`, name = "n_prescriptions") %>%
  mutate(month_index = row_number())
bovine_monthly
bovine_trend_model <- lm(n_prescriptions ~ month_index, data = bovine_monthly)
summary(bovine_trend_model)

# ---- Step 11: Non-linear (quadratic) test for Bovine ----
bovine_quad_model <- lm(n_prescriptions ~ month_index + I(month_index^2), data = bovine_monthly)
summary(bovine_quad_model)

# ---- Step 12: Canine and Feline monthly trend tests ----
canine_monthly <- amu %>%
  filter(`Species Name` == "Canine") %>%
  count(`Record Year`, `Record Month`, name = "n_prescriptions") %>%
  mutate(month_index = row_number())
canine_trend <- lm(n_prescriptions ~ month_index, data = canine_monthly)
summary(canine_trend)

feline_monthly <- amu %>%
  filter(`Species Name` == "Feline") %>%
  count(`Record Year`, `Record Month`, name = "n_prescriptions") %>%
  mutate(month_index = row_number())
feline_trend <- lm(n_prescriptions ~ month_index, data = feline_monthly)
summary(feline_trend)

# ---- Step 13: Top 10 antimicrobials ----
top_drugs <- amu %>%
  count(`Medicine Name`, name = "n_prescriptions") %>%
  arrange(desc(n_prescriptions)) %>%
  slice_head(n = 10)
top_drugs

# ---- Step 13b: Plot top drugs [Figure 3] ----
fig3 <- ggplot(top_drugs, aes(x = reorder(`Medicine Name`, n_prescriptions), y = n_prescriptions)) +
  geom_col(fill = "steelblue") +
  coord_flip() +
  labs(title = "Top 10 Antimicrobials Prescribed, Mongar (2023-2026)",
       x = NULL, y = "Number of Prescriptions") +
  theme_minimal()
fig3
ggsave("Figure3_Top_Drugs.png", fig3, width = 8, height = 5, dpi = 300)

# ---- Step 14: Group into broader drug classes ----
amu_classed <- amu %>%
  mutate(drug_class = case_when(
    str_detect(`Medicine Name`, regex("tetracycline", ignore_case = TRUE)) ~ "Tetracyclines",
    str_detect(`Medicine Name`, regex("penicillin|amoxycillin|ampicillin|cloxacillin", ignore_case = TRUE)) ~ "Penicillins",
    str_detect(`Medicine Name`, regex("sulpha|trimethoprim", ignore_case = TRUE)) ~ "Sulphonamides",
    str_detect(`Medicine Name`, regex("enrofloxacin|ciprofloxacin|fluoroquinolone", ignore_case = TRUE)) ~ "Fluoroquinolones",
    str_detect(`Medicine Name`, regex("gentamicin|streptomycin", ignore_case = TRUE)) ~ "Aminoglycosides",
    TRUE ~ "Other"
  ))
class_summary <- amu_classed %>%
  count(drug_class, name = "n_prescriptions") %>%
  arrange(desc(n_prescriptions))
class_summary

# ---- Step 15: Species x drug class association ----
amu_main3 <- amu_classed %>%
  filter(`Species Name` %in% c("Bovine", "Canine", "Feline"))
species_class_table_main3 <- table(amu_main3$`Species Name`, amu_main3$drug_class)
species_class_table_main3
chisq.test(species_class_table_main3, simulate.p.value = TRUE, B = 10000)

# ---- Step 16: Facility type - DVH Mongar vs RNR-EC ----
amu_main3 <- amu_main3 %>%
  mutate(facility_type = ifelse(Jurisdiction == "Mongar", "DVH Mongar", "RNR-EC"))
table(amu_main3$facility_type)
facility_class_table <- table(amu_main3$facility_type, amu_main3$drug_class)
facility_class_table
chisq.test(facility_class_table, simulate.p.value = TRUE, B = 10000)

# ---- Step 17: Proportional breakdown by facility type ----
facility_class_prop <- amu_main3 %>%
  count(facility_type, drug_class) %>%
  group_by(facility_type) %>%
  mutate(percent = round(n / sum(n) * 100, 1)) %>%
  arrange(facility_type, desc(percent))
print(facility_class_prop, n = 20)

# ---- Step 18: Species distribution by facility type ----
species_facility_table <- table(amu_main3$facility_type, amu_main3$`Species Name`)
species_facility_table
prop.table(species_facility_table, margin = 1) * 100

# ---- Step 19: Facility effect within Bovine only (controls for species) ----
bovine_only <- amu_main3 %>% filter(`Species Name` == "Bovine")
bovine_facility_table <- table(bovine_only$facility_type, bovine_only$drug_class)
bovine_facility_table
chisq.test(bovine_facility_table, simulate.p.value = TRUE, B = 10000)

# ---- Step 20: Seasonal (month-of-year) pattern ----
monthly_pattern <- amu %>%
  count(`Record Month`, name = "n_prescriptions")
monthly_pattern

monthly_pattern_avg <- amu %>%
  count(`Record Year`, `Record Month`, name = "n") %>%
  group_by(`Record Month`) %>%
  summarise(avg_prescriptions = round(mean(n), 1), n_years = n())
monthly_pattern_avg

chisq.test(monthly_pattern_avg$avg_prescriptions)

# ---- Step 20b: Plot seasonal pattern [Figure 4] ----
fig4 <- ggplot(monthly_pattern_avg, aes(x = factor(`Record Month`), y = avg_prescriptions)) +
  geom_col(fill = "darkorange") +
  labs(title = "Seasonal Pattern in Antimicrobial Prescribing, Mongar (2023-2026)",
       subtitle = "Year-adjusted monthly averages",
       x = "Month", y = "Average Number of Prescriptions") +
  theme_minimal()
fig4
ggsave("Figure4_Seasonality.png", fig4, width = 8, height = 5, dpi = 300)

# ---- Step 21: Export a de-identified version of the dataset ----
# (for data-sharing purposes only; not used in any analysis above)

deidentified <- df %>%
  select(-`Owner Name`, -CID, -`Mobile No.`, -Nationality)

# Confirm the identifiable columns are gone
colnames(deidentified)

# Save as a separate file - keep this distinct from your working data
write.csv(deidentified, "AMU_data_deidentified.csv", row.names = FALSE)

# ---- Step 22: Build visit-level dataset and test the % antimicrobial metric ----

# A "visit" = unique combination of Patient Id + Treatment Date
# First, check how many visits have MULTIPLE rows (i.e., multiple drugs
# prescribed in one visit) - this tells us if visit-level counting is clean

visit_summary <- df %>%
  group_by(`Patient Id`, `Treatment Date`) %>%
  summarise(
    n_line_items = n(),
    any_antimicrobial = any(`Medicine Class` == "Antimicrobials", na.rm = TRUE),
    .groups = "drop"
  )

# Total unique visits - should be around 8,881
nrow(visit_summary)

# How many visits had more than one drug/line item recorded?
table(visit_summary$n_line_items > 1)

# The key number: % of visits where at least one antimicrobial was prescribed
mean(visit_summary$any_antimicrobial) * 100

# ---- Step 23: Investigate the facility-level prescribing difference ----
# Testing three plausible explanations within Bovine cases only:
# (1) Route - injectable use may signal more severe/systemic disease
# (2) Diagnostic testing - more test-based diagnosis may signal referral of complex cases
# (3) Duration - longer courses may signal more complex/chronic disease

# Rebuild bovine_only if needed (from amu_main3, filtered to Bovine)
bovine_only <- amu_main3 %>% filter(`Species Name` == "Bovine")

# ---- (1) Route: injectable vs non-injectable, by facility ----
bovine_only <- bovine_only %>%
  mutate(route_type = ifelse(str_detect(Route, regex("inj|intra|parenteral", ignore_case = TRUE)),
                             "Injectable", "Non-injectable"))

route_table <- table(bovine_only$facility_type, bovine_only$route_type)
route_table
chisq.test(route_table, simulate.p.value = TRUE, B = 10000)

# ---- (2) Diagnostic testing: was a test requested, by facility ----
bovine_only <- bovine_only %>%
  mutate(test_requested = ifelse(is.na(`Diagnostic Test Request`) | `Diagnostic Test Request` == "",
                                 "No test", "Test requested"))

test_table <- table(bovine_only$facility_type, bovine_only$test_requested)
test_table
chisq.test(test_table, simulate.p.value = TRUE, B = 10000)

# ---- (3) Duration: average treatment course length, by facility ----
bovine_only %>%
  group_by(facility_type) %>%
  summarise(
    mean_duration = mean(Duration, na.rm = TRUE),
    median_duration = median(Duration, na.rm = TRUE),
    n = n()
  )

# t-test comparing durations
t.test(Duration ~ facility_type, data = bovine_only)

# ---- Step 24: Cramér's V with bootstrapped 95% CI for chi-square results ----
# install.packages("rcompanion")  # only needed once
install.packages("rcompanion")
library(rcompanion)


# Species x drug class (Bovine/Canine/Feline)
cramerV(species_class_table_main3, ci = TRUE, R = 1000)

# Facility type x drug class (all species)
cramerV(facility_class_table, ci = TRUE, R = 1000)

# Facility type x drug class (Bovine only, controlling for species)
cramerV(bovine_facility_table, ci = TRUE, R = 1000)

# Route of administration x facility type (Bovine only)
cramerV(route_table, ci = TRUE, R = 1000)
# ================================================================
# End of analysis - Figures 1-4 saved to D:/AMR Manuscript
# ================================================================