## ============================================================================= ## Diagnosis-Associated Antimicrobial Prescribing Patterns in a Resource-Limited ## Veterinary District Setting: Evidence from Mongar, Bhutan (2023-2026) ## ## FULL ANALYSIS PIPELINE - Additional file 3 ## Author: Narayan Pokhrel, Cheda ## District Veterinary Hospital, Mongar, Bhutan ## ## HOW TO USE THIS SCRIPT: ## Run every step below, in order, from Step 1 through the end. Each step ## produces output described in the note above it. All numbers produced here ## should match the corresponding tables, figures, and in-text statistics in ## the manuscript. Update the file path in Step 1 if your data file is stored ## elsewhere. ## ============================================================================= ## ----------------------------------------------------------------------------- ## STEP 1: Load required packages and the raw dataset ## ----------------------------------------------------------------------------- ## Note: Installs and loads all packages used anywhere in this script. Update ## the file path below to point to your local copy of AMU_data.xlsx. required_packages <- c("readxl", "dplyr", "tidyr", "binom", "vcd", "ggplot2", "writexl") new_packages <- required_packages[!(required_packages %in% installed.packages()[,"Package"])] if (length(new_packages) > 0) install.packages(new_packages, repos = "https://cloud.r-project.org") library(readxl) library(dplyr) library(tidyr) library(binom) library(vcd) library(ggplot2) library(writexl) # --- UPDATE THIS PATH to your local file location --- df <- read_excel("D:/Paper3_Antimicrobialprescribing/AMU_data.xlsx", sheet = 1) cat("Step 1 complete: loaded", nrow(df), "medicine-dispensing line items.\n") ## ----------------------------------------------------------------------------- ## STEP 2: Build the visit-level dataset ## ----------------------------------------------------------------------------- ## Note: A single clinical visit (one animal, one treatment date) can generate ## more than one medicine line item. This step collapses the raw line-item ## data down to one row per visit, which is the correct unit of analysis for ## overall prescribing rate and diagnosis-condition-specific rate (Table 1, ## Table 4, Figure 3). Where a visit had more than one diagnosis condition ## recorded, the modal (most frequently repeated) value is used. ## ## NOTE ON TIE-BREAKING: for a small number of visits with a genuine tie between ## two equally-frequent diagnosis conditions, get_mode()'s tie-breaking behavior ## can vary very slightly by R version/locale. This affects at most a handful of ## visits and does not change any statistical conclusion, but if you see small ## (+/-1 to 5) differences in visit counts for a few diagnosis conditions compared ## to the published tables, this is the expected cause. get_mode <- function(x) { x <- x[!is.na(x)] if (length(x) == 0) return(NA) names(sort(table(x), decreasing = TRUE))[1] } df$visit_id <- paste(df$`Patient Id`, df$`Treatment Date`, sep = "_") visits <- df %>% group_by(visit_id) %>% summarise( species = get_mode(`Species Name`), diag_system = get_mode(`Final Diagnosis-System`), diag_condition= get_mode(`Final Diagnosis-Condition`), record_year = get_mode(`Record Year`), am_given = any(`Medicine Class` == "Antimicrobials", na.rm = TRUE), .groups = "drop" ) cat("Step 2 complete: built", nrow(visits), "unique clinical visits.\n") ## ----------------------------------------------------------------------------- ## STEP 3: Identify and characterise the "condition1" placeholder diagnosis ## ----------------------------------------------------------------------------- ## Note: A data-entry default value ("condition1") appears in place of a real ## diagnosis for a small number of visits. This step quantifies how many, ## confirms they are retained in the overall denominator, and characterises ## their species distribution and antimicrobial rate (reported in Methods ## and discussed as a limitation). visits_clean <- visits %>% filter(diag_condition != "condition1") condition1_visits <- visits %>% filter(diag_condition == "condition1") n_total <- nrow(visits) n_condition1 <- nrow(condition1_visits) n_clean <- nrow(visits_clean) am_total <- sum(visits$am_given) am_clean <- sum(visits_clean$am_given) am_condition1 <- sum(condition1_visits$am_given) cat("Step 3 complete:\n") cat(" Total visits:", n_total, "\n") cat(" condition1 (placeholder) visits:", n_condition1, sprintf("(%.1f%% of all visits)\n", 100 * n_condition1 / n_total)) cat(" Genuine-diagnosis visits:", n_clean, "\n") cat(" AM-treated, full sample:", am_total, sprintf("(%.2f%%)\n", 100 * am_total / n_total)) cat(" AM-treated, genuine-diagnosis only:", am_clean, sprintf("(%.2f%%)\n", 100 * am_clean / n_clean)) cat(" AM-treated among condition1 visits:", am_condition1, sprintf("(%.1f%%)\n", 100 * am_condition1 / n_condition1)) cat(" Species distribution of condition1 visits:\n") print(condition1_visits %>% count(species) %>% mutate(pct = round(100 * n / sum(n), 1))) ## ----------------------------------------------------------------------------- ## STEP 4: Overall antimicrobial prescribing rate, with 95% Wilson CI (Table 1) ## ----------------------------------------------------------------------------- ## Note: Produces the two headline prescribing-rate figures reported in the ## Abstract and Results: the rate across all visits, and the rate restricted ## to visits with a genuine (non-placeholder) diagnosis. ci_full <- binom.confint(am_total, n_total, method = "wilson") ci_clean <- binom.confint(am_clean, n_clean, method = "wilson") cat("Step 4 complete:\n") cat(sprintf(" Full sample: %.2f%% (95%% CI: %.1f-%.1f%%), n=%d/%d\n", 100*ci_full$mean, 100*ci_full$lower, 100*ci_full$upper, am_total, n_total)) cat(sprintf(" Clean sample: %.2f%% (95%% CI: %.1f-%.1f%%), n=%d/%d\n", 100*ci_clean$mean, 100*ci_clean$lower, 100*ci_clean$upper, am_clean, n_clean)) ## ----------------------------------------------------------------------------- ## STEP 5: Species and diagnosis-system distribution (Table 1) ## ----------------------------------------------------------------------------- ## Note: Descriptive breakdown of the full visit sample by species and by ## diagnosis system, as reported in Table 1 of the manuscript. species_table <- visits %>% count(species, sort = TRUE) %>% mutate(pct = round(100 * n / sum(n), 1)) diag_system_table <- visits %>% count(diag_system, sort = TRUE) %>% mutate(pct = round(100 * n / sum(n), 1)) cat("Step 5 complete. Species distribution:\n") print(species_table) cat("\nDiagnosis system distribution:\n") print(diag_system_table) write.csv(species_table, "Table1a_species_distribution.csv", row.names = FALSE) write.csv(diag_system_table, "Table1b_diagnosis_system_distribution.csv", row.names = FALSE) cat(" Saved: Table1a_species_distribution.csv and Table1b_diagnosis_system_distribution.csv\n") ## ----------------------------------------------------------------------------- ## STEP 6: Build the antimicrobial classification lookup table ## ----------------------------------------------------------------------------- ## Note: Assigns each of the 21 antimicrobial formulations to a pharmacological ## class and a WHO 2023 AWaRe tier (Access/Watch/Not classified). This ## classification is used throughout Tables 2, 3, and 4, Figure 4, and the ## logistic regression models. Corresponds to Supplementary Table S1. drug_lookup <- tibble::tribble( ~Medicine_Name, ~Pharm_Class, ~AWaRe, "Tetracycline Hcl W/S", "Tetracyclines", "Access", "Oxytetracycline LA Inj.", "Tetracyclines", "Watch", "Oxytetracycline SA", "Tetracyclines", "Watch", "Doxycycline hyclate tab", "Tetracyclines", "Access", "Amoxycillin Trihydrate bolus", "Penicillins", "Access", "Amoxycillin Trihydrate Susp. (250 mg)", "Penicillins", "Access", "Benzathine Penicillin LA Inj.", "Penicillins", "Access", "Benzathine Penicillin Inj.", "Penicillins", "Access", "Ampicillin + Cloxacillin Inj.", "Penicillins", "Access", "Sulphadiazine + Trimethoprim powder", "Sulfonamides+Trimethoprim", "Access", "Trimethoprim + Sulphadiazine bolus", "Sulfonamides+Trimethoprim", "Access", "Sulphadimidine Inj.", "Sulfonamides", "Access", "Sulfaquinoxaline, Diaveridine, Menadione", "Sulfonamides (combination)", "Not classified", "Enrofloxacin tab", "Fluoroquinolones", "Watch", "Enrofloxacin inj.", "Fluoroquinolones", "Watch", "Gentamicin inj.", "Aminoglycosides", "Access", "Cefotaxime Inj.", "Cephalosporins (3rd gen)", "Watch", "Cephalexin bolus", "Cephalosporins (1st gen)", "Access", "Strepto-penicillin Inj.", "Aminoglycoside+Penicillin", "Access", "Metronidazole inj.", "Nitroimidazoles", "Access", "Metronidazole Tab", "Nitroimidazoles", "Access" ) cat("Step 6 complete: classification lookup built for", nrow(drug_lookup), "formulations.\n") ## ----------------------------------------------------------------------------- ## STEP 7: Build the antimicrobial-treated line-item dataset ## ----------------------------------------------------------------------------- ## Note: Filters to antimicrobial line items only, excludes those linked to ## the condition1 placeholder diagnosis, and joins the classification lookup. ## This is the dataset (n=1008 line items) used for Tables 2, 3, and the ## diagnosis-drug-class and AWaRe-counseling analyses. Unlike the visit-level ## dataset, a visit contributing 2 antimicrobials contributes 2 rows here. amu <- df %>% filter(`Medicine Class` == "Antimicrobials", `Final Diagnosis-Condition` != "condition1") %>% left_join(drug_lookup, by = c("Medicine Name" = "Medicine_Name")) n_line_items <- nrow(amu) cat("Step 7 complete:", n_line_items, "antimicrobial line items (diagnosis-linked subset).\n") # Check for multi-antimicrobial visits and mixed-AWaRe-tier visits (reported in Methods) am_per_visit <- amu %>% count(visit_id, name = "n_am_lines") n_multi_visit <- sum(am_per_visit$n_am_lines > 1) tier_per_visit <- amu %>% group_by(visit_id) %>% summarise(n_tiers = n_distinct(AWaRe), .groups = "drop") n_mixed_tier <- sum(tier_per_visit$n_tiers > 1) cat(" Visits with >1 antimicrobial line item:", n_multi_visit, "\n") cat(" Visits with antimicrobials from >1 AWaRe tier:", n_mixed_tier, "\n") ## ----------------------------------------------------------------------------- ## STEP 8: Table 2 - Diagnosis system x pharmacological class cross-tabulation ## ----------------------------------------------------------------------------- ## Note: The core statistical test of the paper. Produces the chi-square ## statistic, a Monte Carlo simulated p-value (robust to the sparse cells in ## this table), and Cramer's V as an effect-size measure. t2_table <- table(amu$`Final Diagnosis-System`, amu$Pharm_Class) # Restrict the significance test to categories with >=10 line items, to avoid # an unstable test based on very sparse categories. Table 2 itself (the # descriptive counts reported in the manuscript) keeps all categories. sys_keep <- names(which(rowSums(t2_table) >= 10)) class_keep <- names(which(colSums(t2_table) >= 10)) t2_test <- t2_table[sys_keep, class_keep] chisq_result <- suppressWarnings(chisq.test(t2_test)) mc_result <- chisq.test(t2_test, simulate.p.value = TRUE, B = 10000) cramers_v <- assocstats(t2_test)$cramer n_sparse_cells <- sum(chisq_result$expected < 5) pct_sparse <- round(100 * n_sparse_cells / length(chisq_result$expected), 1) cat("Step 8 complete:\n") cat(sprintf(" Chi-square = %.2f, df = %d, p < 0.001\n", chisq_result$statistic, chisq_result$parameter)) cat(sprintf(" Cells with expected count <5: %d of %d (%.1f%%)\n", n_sparse_cells, length(chisq_result$expected), pct_sparse)) cat(sprintf(" Monte Carlo simulated p-value: %s\n", format.pval(mc_result$p.value))) cat(sprintf(" Cramer's V (effect size): %.3f\n", cramers_v)) write.csv(as.data.frame.matrix(t2_table), "Table2_diagnosis_by_drugclass.csv") cat(" Saved: Table2_diagnosis_by_drugclass.csv\n") ## ----------------------------------------------------------------------------- ## STEP 9: Sensitivity analysis - diagnosis-drug association within each species ## ----------------------------------------------------------------------------- ## Note: Tests whether the diagnosis-drug class association in Step 8 holds ## independently within cattle, dogs, and cats, to rule out species as a ## confound (since diagnosis system is itself associated with species). cat("Step 9: species-stratified diagnosis-drug association\n") for (sp in c("Bovine", "Canine", "Feline")) { sub <- amu %>% filter(`Species Name` == sp) tab <- table(sub$`Final Diagnosis-System`, sub$Pharm_Class) sk <- names(which(rowSums(tab) >= 5)) ck <- names(which(colSums(tab) >= 5)) if (length(sk) >= 2 && length(ck) >= 2) { tt <- tab[sk, ck, drop = FALSE] res <- suppressWarnings(chisq.test(tt, simulate.p.value = TRUE, B = 10000)) cat(sprintf(" %s: n=%d, chi-sq=%.2f, Monte Carlo p=%s\n", sp, sum(tab), res$statistic, format.pval(res$p.value))) } } ## ----------------------------------------------------------------------------- ## STEP 10: Sensitivity analysis - one antimicrobial per visit (independence check) ## ----------------------------------------------------------------------------- ## Note: Because some visits contribute more than one antimicrobial line item, ## this step repeats the Step 8 test restricted to one antimicrobial per visit ## (the first-listed line item, a prespecified rule), to confirm the result is ## not an artefact of non-independent observations within visits. visit_level_drug <- amu %>% group_by(visit_id) %>% slice(1) %>% ungroup() t2_visit <- table(visit_level_drug$`Final Diagnosis-System`, visit_level_drug$Pharm_Class) sk2 <- names(which(rowSums(t2_visit) >= 10)) ck2 <- names(which(colSums(t2_visit) >= 10)) t2_visit_test <- t2_visit[sk2, ck2] chisq_visit <- suppressWarnings(chisq.test(t2_visit_test)) mc_visit <- chisq.test(t2_visit_test, simulate.p.value = TRUE, B = 10000) cat("Step 10 complete:\n") cat(" Visit-level dataset n =", nrow(visit_level_drug), "(one antimicrobial per visit)\n") cat(sprintf(" Chi-square = %.2f, df = %d, Monte Carlo p = %s\n", chisq_visit$statistic, chisq_visit$parameter, format.pval(mc_visit$p.value))) ## ----------------------------------------------------------------------------- ## STEP 11: Table 3 - WHO AWaRe tier distribution by diagnosis system ## ----------------------------------------------------------------------------- ## Note: Cross-tabulates AWaRe tier (Access/Watch/Not classified) against ## diagnosis system, with percentages calculated among AWaRe-classified line ## items only. t3 <- amu %>% count(`Final Diagnosis-System`, AWaRe) %>% pivot_wider(names_from = AWaRe, values_from = n, values_fill = 0) %>% mutate( Total = Access + Watch + `Not classified`, # NOTE: percentages use Total (including any "Not classified" line items) as the # denominator, matching the published Table 3. Only 2 of 15 diagnosis systems # (Digestive System, General clinical examination) have a nonzero "Not classified" # count, so this differs from an Access/(Access+Watch)-only denominator by at most # ~0.3 percentage points on those two rows. Access_pct = round(100 * Access / Total, 1), Watch_pct = round(100 * Watch / Total, 1) ) %>% arrange(desc(Total)) cat("Step 11 complete. Table 3:\n") print(t3) write.csv(t3, "Table3_AWaRe_by_diagnosis_system.csv", row.names = FALSE) cat(" Saved: Table3_AWaRe_by_diagnosis_system.csv\n") ## ----------------------------------------------------------------------------- ## STEP 12: Table 4 - Antimicrobial prescribing rate by diagnosis condition ## ----------------------------------------------------------------------------- ## Note: Visit-level prescribing rate for each diagnosis condition with n>=20 ## visits, with 95% Wilson confidence intervals. This is the basis for Table 4 ## and Figure 3. Also computes the predominant species for each condition. diag_rates <- visits_clean %>% filter(!is.na(diag_condition)) %>% group_by(diag_condition) %>% summarise(n_visits = n(), n_am = sum(am_given), .groups = "drop") %>% filter(n_visits >= 20) ci_list <- binom.confint(diag_rates$n_am, diag_rates$n_visits, method = "wilson") diag_rates_ci <- diag_rates %>% mutate( am_rate = round(100 * n_am / n_visits, 1), ci_lower = round(100 * ci_list$lower, 1), ci_upper = round(100 * ci_list$upper, 1) ) %>% arrange(desc(am_rate)) # Predominant species per condition species_per_condition <- visits_clean %>% filter(diag_condition %in% diag_rates_ci$diag_condition) %>% count(diag_condition, species) %>% group_by(diag_condition) %>% mutate(pct = round(100 * n / sum(n), 1)) %>% slice_max(n, n = 1) %>% ungroup() %>% select(diag_condition, predominant_species = species, species_pct = pct) table4_full <- diag_rates_ci %>% left_join(species_per_condition, by = "diag_condition") cat("Step 12 complete. Table 4 (", nrow(table4_full), "conditions with n>=20 visits):\n") print(table4_full, n = 40) write.csv(table4_full, "Table4_AM_rate_by_diagnosis_with_species.csv", row.names = FALSE) cat(" Saved: Table4_AM_rate_by_diagnosis_with_species.csv\n") overall_rate <- round(mean(visits_clean$am_given) * 100, 1) cat(" Overall rate among visits with identifiable diagnosis (reference line):", overall_rate, "%\n") ## ----------------------------------------------------------------------------- ## STEP 13: Supplementary Table S2 - full diagnosis condition list (all n) ## ----------------------------------------------------------------------------- ## Note: Same as Step 12 but with no minimum sample-size restriction, covering ## all 129 diagnosis conditions recorded in the dataset. full_conditions <- visits_clean %>% filter(!is.na(diag_condition)) %>% group_by(diag_condition) %>% summarise(n_visits = n(), n_am = sum(am_given), am_rate = round(100 * n_am / n_visits, 1), .groups = "drop") %>% arrange(desc(n_visits)) cat("Step 13 complete:", nrow(full_conditions), "total diagnosis conditions.\n") write.csv(full_conditions, "SupplementaryTableS2_full_conditions.csv", row.names = FALSE) cat(" Saved: SupplementaryTableS2_full_conditions.csv\n") ## ----------------------------------------------------------------------------- ## STEP 14: Oxytetracycline classification sensitivity analysis (Supp. Table S3) ## ----------------------------------------------------------------------------- ## Note: Oxytetracycline accounts for a large share of antimicrobial use and ## is classified as Watch tier under WHO 2023 AWaRe. This step tests how much ## the food-producing-vs-companion-animal Watch-tier comparison depends on ## this specific classification choice (vs. classifying it as Access, as in ## some agricultural/WOAH-oriented frameworks). food_species <- c("Bovine", "Caprine", "Swine", "Poultry") companion_species <- c("Canine", "Feline") amu_scenario <- amu %>% mutate( species_group = case_when( `Species Name` %in% food_species ~ "Food-producing", `Species Name` %in% companion_species ~ "Companion", TRUE ~ "Other" ), AWaRe_scenario_B = ifelse(`Medicine Name` %in% c("Oxytetracycline LA Inj.", "Oxytetracycline SA"), "Access", AWaRe) ) oxy_share <- amu %>% filter(`Medicine Name` %in% c("Oxytetracycline LA Inj.", "Oxytetracycline SA")) %>% nrow() cat("Step 14 complete:\n") cat(sprintf(" Oxytetracycline: %d of %d line items (%.1f%%)\n", oxy_share, n_line_items, 100 * oxy_share / n_line_items)) for (grp in c("Food-producing", "Companion")) { sub <- amu_scenario %>% filter(species_group == grp) cat(sprintf("\n %s (n=%d):\n", grp, nrow(sub))) cat(" Scenario A (WHO 2023, oxytetracycline=Watch): ") print(round(100 * prop.table(table(sub$AWaRe)), 1)) cat(" Scenario B (oxytetracycline=Access): ") print(round(100 * prop.table(table(sub$AWaRe_scenario_B)), 1)) } ## ----------------------------------------------------------------------------- ## STEP 15: Course-completion counseling indicator ## ----------------------------------------------------------------------------- ## Note: Screens the free-text "Advice" field for course-completion-related ## keywords, producing a binary indicator for each antimicrobial line item. ## This indicator was separately validated against two independent human ## raters on a random sample of 150 entries (sensitivity 91.4%, specificity ## 97.0-98.5%, accuracy 93.9-94.6%; see Methods for full validation details - ## that exercise is not reproduced in this script as it required manual review). amu <- amu %>% mutate( advice_lower = tolower(as.character(Advice)), course_advice = grepl("complete|course|finish|full dose", advice_lower) ) overall_counsel_rate <- round(100 * mean(amu$course_advice), 1) cat("Step 15 complete: overall course-completion counseling rate =", overall_counsel_rate, "%\n") ## ----------------------------------------------------------------------------- ## STEP 16: Counseling by AWaRe tier (Figure 4) and by route of administration ## ----------------------------------------------------------------------------- ## Note: The two descriptive comparisons that motivate the adjusted logistic ## regression in Step 17. Produces the crude Access-vs-Watch counseling gap, ## the route-of-administration breakdown, and confirms route and AWaRe tier ## are themselves closely related (the basis for suspecting confounding). amu_classified <- amu %>% filter(AWaRe %in% c("Access", "Watch")) counsel_by_aware <- amu_classified %>% group_by(AWaRe) %>% summarise(n = n(), n_advice = sum(course_advice), pct = round(100 * n_advice / n, 1), .groups = "drop") chisq_aware <- chisq.test(table(amu_classified$AWaRe, amu_classified$course_advice), correct = TRUE) amu <- amu %>% mutate(Route2 = ifelse(Route == "Oral", "Oral", "Injectable/Other")) counsel_by_route <- amu %>% group_by(Route) %>% summarise(n = n(), pct_advice = round(100 * mean(course_advice), 1), .groups = "drop") route_by_aware <- amu_classified %>% mutate(Route2 = ifelse(Route == "Oral", "Oral", "Injectable/Other")) %>% count(Route2, AWaRe) %>% pivot_wider(names_from = AWaRe, values_from = n, values_fill = 0) cat("Step 16 complete.\n") cat("Counseling by AWaRe tier:\n"); print(counsel_by_aware) cat(sprintf(" Chi-square = %.1f, df = 1, p < 0.001\n", chisq_aware$statistic)) cat("\nCounseling by route of administration:\n"); print(counsel_by_route) cat("\nRoute by AWaRe tier (confound check):\n"); print(route_by_aware) ## ----------------------------------------------------------------------------- ## STEP 17: Adjusted logistic regression - counseling ~ AWaRe tier + route ## ----------------------------------------------------------------------------- ## Note: The key adjusted analysis. Tests whether AWaRe tier remains ## independently associated with counseling after accounting for route of ## administration. Reports unadjusted and adjusted odds ratios with 95% CIs. m_unadjusted <- glm(course_advice ~ AWaRe, data = amu_classified, family = binomial) m_adjusted <- glm(course_advice ~ AWaRe + Route2, data = amu_classified %>% mutate(Route2 = ifelse(Route == "Oral", "Oral", "Injectable/Other")), family = binomial) or_unadj <- exp(cbind(OR = coef(m_unadjusted), suppressMessages(confint(m_unadjusted)))) or_adj <- exp(cbind(OR = coef(m_adjusted), suppressMessages(confint(m_adjusted)))) cat("Step 17 complete.\n") cat("Unadjusted model (counseling ~ AWaRe tier):\n"); print(round(or_unadj, 3)) cat("\nAdjusted model (counseling ~ AWaRe tier + route):\n"); print(round(or_adj, 3)) cat("\nAIC comparison - full model vs. route-only model:\n") m_route_only <- glm(course_advice ~ Route2, data = amu_classified %>% mutate(Route2 = ifelse(Route == "Oral", "Oral", "Injectable/Other")), family = binomial) cat(sprintf(" AWaRe+Route AIC = %.1f | Route-only AIC = %.1f\n", AIC(m_adjusted), AIC(m_route_only))) ## ----------------------------------------------------------------------------- ## STEP 18: Visit-level sensitivity analysis for the counseling model ## ----------------------------------------------------------------------------- ## Note: Repeats Steps 16-17 restricted to one antimicrobial per visit, to ## confirm the counseling-gap finding is not an artefact of visits that ## contributed multiple antimicrobial line items. visit_level_counsel <- amu_classified %>% mutate(Route2 = ifelse(Route == "Oral", "Oral", "Injectable/Other")) %>% group_by(visit_id) %>% slice(1) %>% ungroup() chisq_visit_counsel <- chisq.test(table(visit_level_counsel$AWaRe, visit_level_counsel$course_advice), correct = TRUE) m_visit_adj <- glm(course_advice ~ AWaRe + Route2, data = visit_level_counsel, family = binomial) or_visit_adj <- exp(cbind(OR = coef(m_visit_adj), suppressMessages(confint(m_visit_adj)))) cat("Step 18 complete.\n") cat(" Visit-level dataset n =", nrow(visit_level_counsel), "\n") cat(sprintf(" Chi-square = %.2f, p < 0.001\n", chisq_visit_counsel$statistic)) cat(" Adjusted odds ratio (visit-level):\n"); print(round(or_visit_adj, 3)) ## ----------------------------------------------------------------------------- ## STEP 19: Facility-level comparison - DVH versus the RNR-Extension Centres ## ----------------------------------------------------------------------------- ## Note: Tests whether the three main outcomes (diagnosis-drug class ## association, AWaRe tier distribution, counseling rate) differ between DVH ## and the district's RNR-ECs, which serve markedly different caseloads. amu <- amu %>% mutate(facility = ifelse(`Centre Name` == "DVH", "DVH", "RNR-EC")) species_by_facility <- amu %>% count(facility, `Species Name`) %>% group_by(facility) %>% mutate(pct = round(100 * n / sum(n), 1)) cat("Step 19a: Species mix by facility\n") print(species_by_facility, n = 20) cat("\nStep 19b: Diagnosis-drug class association, tested separately per facility\n") for (fac in c("DVH", "RNR-EC")) { sub <- amu %>% filter(facility == fac) tab <- table(sub$`Final Diagnosis-System`, sub$Pharm_Class) sk <- names(which(rowSums(tab) >= 5)) ck <- names(which(colSums(tab) >= 5)) if (length(sk) >= 2 && length(ck) >= 2) { tt <- tab[sk, ck, drop = FALSE] res <- suppressWarnings(chisq.test(tt, simulate.p.value = TRUE, B = 10000)) cat(sprintf(" %s: n=%d, chi-sq=%.2f, Monte Carlo p=%s\n", fac, sum(tab), res$statistic, format.pval(res$p.value))) } } amu_classified <- amu_classified %>% mutate(facility = ifelse(`Centre Name` == "DVH", "DVH", "RNR-EC")) aware_by_facility <- amu_classified %>% count(facility, AWaRe) %>% group_by(facility) %>% mutate(pct = round(100 * n / sum(n), 1)) chisq_facility_aware <- chisq.test(table(amu_classified$facility, amu_classified$AWaRe == "Watch")) counsel_by_facility <- amu_classified %>% group_by(facility) %>% summarise(n = n(), pct_advice = round(100 * mean(course_advice), 1), .groups = "drop") chisq_facility_counsel <- chisq.test(table(amu_classified$facility, amu_classified$course_advice)) oxy_by_facility <- amu %>% mutate(is_oxy = `Medicine Name` %in% c("Oxytetracycline LA Inj.", "Oxytetracycline SA")) %>% group_by(facility) %>% summarise(n = n(), n_oxy = sum(is_oxy), pct_oxy = round(100 * n_oxy / n, 1), .groups = "drop") cat("\nStep 19c: AWaRe tier by facility\n"); print(aware_by_facility) cat(sprintf(" Chi-square = %.2f, p = %.4f\n", chisq_facility_aware$statistic, chisq_facility_aware$p.value)) cat("\nStep 19d: Counseling rate by facility\n"); print(counsel_by_facility) cat(sprintf(" Chi-square = %.2f, p = %.3f\n", chisq_facility_counsel$statistic, chisq_facility_counsel$p.value)) cat("\nStep 19e: Oxytetracycline share by facility (explains the AWaRe difference)\n"); print(oxy_by_facility) cat("Step 19 complete.\n") ## ----------------------------------------------------------------------------- ## STEP 20: Fuller multivariable model - AWaRe + route + species + facility ## ----------------------------------------------------------------------------- ## Note: Extends the Step 17 adjusted model to check whether species and ## facility type, both plausible confounders, change the AWaRe/route ## conclusion. Diagnosis system is deliberately excluded: several of its 15 ## categories have fewer than 15 line items, risking unstable estimates. amu_classified <- amu_classified %>% mutate(Route2 = ifelse(Route == "Oral", "Oral", "Injectable/Other"), species_grp = case_when( `Species Name` == "Bovine" ~ "Bovine", `Species Name` %in% c("Canine", "Feline") ~ "Companion", TRUE ~ "Other" )) m_full <- glm(course_advice ~ AWaRe + Route2 + species_grp + facility, data = amu_classified, family = binomial) or_full <- exp(cbind(OR = coef(m_full), suppressMessages(confint(m_full)))) lrt <- anova(m_adjusted, m_full, test = "Chisq") cat("Step 20 complete.\n") cat("Fuller model (counseling ~ AWaRe + route + species + facility):\n") print(round(or_full, 3)) cat("\nLikelihood ratio test vs. route-only-adjusted model:\n") print(lrt) ## ----------------------------------------------------------------------------- ## STEP 21: Within-animal clustering check ## ----------------------------------------------------------------------------- ## Note: A visit is defined by patient identifier + treatment date, so the ## same animal can contribute more than one visit. Quantifies repeat visits ## per animal and refits the Step 17 adjusted model with cluster-robust ## standard errors (clustered on animal identifier) to check whether this ## affects the precision or direction of the key adjusted estimate. if (!requireNamespace("sandwich", quietly = TRUE)) install.packages("sandwich", repos = "https://cloud.r-project.org") if (!requireNamespace("lmtest", quietly = TRUE)) install.packages("lmtest", repos = "https://cloud.r-project.org") library(sandwich); library(lmtest) all_visits_animals <- df %>% distinct(visit_id, `Patient Id`) %>% count(`Patient Id`, name = "n_visits") n_unique_animals <- nrow(all_visits_animals) pct_repeat_animals <- round(100 * mean(all_visits_animals$n_visits > 1), 1) pct_visits_from_repeats <- round(100 * sum(all_visits_animals$n_visits[all_visits_animals$n_visits > 1]) / sum(all_visits_animals$n_visits), 1) amu_classified <- amu_classified %>% mutate(animal_id = `Patient Id`) am_animals <- amu_classified %>% distinct(visit_id, animal_id) %>% count(animal_id, name = "n_am_visits") n_unique_am_animals <- nrow(am_animals) pct_repeat_am_animals <- round(100 * mean(am_animals$n_am_visits > 1), 1) cl_vcov <- vcovCL(m_adjusted, cluster = amu_classified$animal_id) cl_ci <- coefci(m_adjusted, vcov = cl_vcov) cl_or <- exp(cbind(OR = coef(m_adjusted), cl_ci)) cat("Step 21 complete.\n") cat(" Unique animals across all", nrow(visits), "visits:", n_unique_animals, "\n") cat(" % of animals with >1 visit:", pct_repeat_animals, "% (accounting for", pct_visits_from_repeats, "% of all visits)\n") cat(" Unique animals among", nrow(amu_classified), "AWaRe-classified AM line items:", n_unique_am_animals, "\n") cat(" % of those animals with >1 such visit:", pct_repeat_am_animals, "%\n") cat(" Cluster-robust OR (clustered on animal identifier):\n") print(round(cl_or, 3)) ## ----------------------------------------------------------------------------- ## STEP 22: Temporal trend - antimicrobial prescribing rate by year ## ----------------------------------------------------------------------------- ## Note: Tests whether the overall prescribing rate changed across the 3.5-year ## study period. yearly_rate <- visits %>% group_by(record_year) %>% summarise(n_visits = n(), n_am = sum(am_given), pct = round(100 * n_am / n_visits, 2), .groups = "drop") chisq_year <- chisq.test(table(visits$record_year, visits$am_given)) cat("Step 22 complete. Annual prescribing rate:\n") print(yearly_rate) cat(sprintf(" Chi-square = %.2f, df = %d, p = %.3f\n", chisq_year$statistic, chisq_year$parameter, chisq_year$p.value)) ## ----------------------------------------------------------------------------- ## STEP 23: Save all results to a single Excel workbook ## ----------------------------------------------------------------------------- ## Note: Consolidates every table produced above into one file for easy ## review and for archiving alongside the manuscript submission. results_list <- list( "Table1_species" = species_table, "Table1_diag_system" = diag_system_table, "Table2_diag_by_class" = as.data.frame(t2_table) %>% rename(DiagnosisSystem = Var1, PharmClass = Var2, n = Freq), "Table3_AWaRe_by_system"= t3, "Table4_diag_rates" = table4_full, "TableS2_full_conditions" = full_conditions, "Counseling_by_AWaRe" = counsel_by_aware, "Counseling_by_route" = counsel_by_route, "Yearly_AM_rate" = yearly_rate ) write_xlsx(results_list, "Paper3_All_Results.xlsx") cat("\nStep 24 complete: saved consolidated results to Paper3_All_Results.xlsx\n") ## ----------------------------------------------------------------------------- ## STEP 24: Generate Figure 2 - Top 10 diagnosis conditions driving AM use ## ----------------------------------------------------------------------------- top_conditions <- amu %>% count(`Final Diagnosis-Condition`, sort = TRUE) %>% slice_head(n = 10) %>% mutate(pct = round(100 * n / n_line_items, 1), label = paste0(n, " (", pct, "%)"), `Final Diagnosis-Condition` = factor(`Final Diagnosis-Condition`, levels = rev(`Final Diagnosis-Condition`))) fig1 <- ggplot(top_conditions, aes(x = `Final Diagnosis-Condition`, y = n)) + geom_col(fill = "#2C7FB8") + geom_text(aes(label = label), hjust = -0.05, size = 4) + coord_flip(clip = "off") + scale_y_continuous(limits = c(0, max(top_conditions$n) * 1.25), expand = c(0, 0)) + labs(title = "Top 10 diagnosis conditions driving antimicrobial prescribing", x = NULL, y = "Number of antimicrobial-treated line items") + theme_minimal(base_size = 13) + theme(plot.title = element_text(face = "bold", size = 15), plot.margin = margin(10, 60, 10, 10)) ggsave("Figure2_top_diagnosis_conditions.png", fig1, width = 10, height = 6, dpi = 300, bg = "white") cat("Step 24 complete: saved Figure2_top_diagnosis_conditions.png\n") ## ----------------------------------------------------------------------------- ## STEP 25: Generate Figure 3 - AM rate by diagnosis condition, with n labels ## ----------------------------------------------------------------------------- diag_rates_ci_fig <- diag_rates_ci %>% mutate(above_overall = am_rate > overall_rate, diag_condition = factor(diag_condition, levels = diag_condition[order(am_rate)])) fig2 <- ggplot(diag_rates_ci_fig, aes(x = diag_condition, y = am_rate, fill = above_overall)) + geom_col(width = 0.7) + geom_hline(yintercept = overall_rate, linetype = "dashed", color = "black", linewidth = 0.6) + geom_text(aes(label = paste0("n=", n_visits), y = am_rate + 2), hjust = 0, size = 3, color = "grey30") + annotate("text", x = 2, y = overall_rate + 6, label = paste0("Overall rate: ", overall_rate, "%"), size = 4, hjust = 0) + coord_flip(clip = "off") + scale_fill_manual(values = c("TRUE" = "#D73027", "FALSE" = "#A6C4E8"), labels = c("TRUE" = "Above overall rate", "FALSE" = "At or below overall rate"), name = NULL) + scale_y_continuous(limits = c(0, 92), expand = c(0, 0)) + labs(title = "Antimicrobial prescribing rate by diagnosis condition\n(n\u226520 visits per condition)", x = NULL, y = "Antimicrobial prescribing rate (%)") + theme_minimal(base_size = 13) + theme(plot.title = element_text(face = "bold", size = 14, hjust = 0), legend.position = "bottom", panel.grid.minor = element_blank(), plot.margin = margin(10, 60, 10, 10)) ggsave("Figure3_disproportionate_AM_use.png", fig2, width = 11, height = 11.5, dpi = 300, bg = "white") cat("Step 25 complete: saved Figure3_disproportionate_AM_use.png\n") ## ----------------------------------------------------------------------------- ## STEP 26: Generate Figure 4 - Counseling rate by AWaRe tier ## ----------------------------------------------------------------------------- fig3_data <- counsel_by_aware %>% mutate(label = paste0(pct, "%\n(n=", n, ")")) fig3 <- ggplot(fig3_data, aes(x = AWaRe, y = pct, fill = AWaRe)) + geom_col(width = 0.5) + geom_text(aes(label = label), vjust = -0.3, size = 5) + scale_fill_manual(values = c("Access" = "#1B9E4B", "Watch" = "#D73027")) + scale_y_continuous(limits = c(0, max(fig3_data$pct) * 1.25)) + labs(title = "Course-completion counseling by AWaRe tier", x = "WHO AWaRe tier", y = "Course-completion counseling rate (%)") + theme_minimal(base_size = 14) + theme(plot.title = element_text(face = "bold", size = 16), legend.position = "none") ggsave("Figure4_counseling_by_AWaRe.png", fig3, width = 7, height = 7, dpi = 300, bg = "white") cat("Step 26 complete: saved Figure4_counseling_by_AWaRe.png\n") ## ============================================================================= ## END OF SCRIPT ## All numbers printed above should match the corresponding tables, figures, ## and in-text statistics in the manuscript. If you obtain different visit or ## line-item counts than those reported in the paper, check that your ## AMU_data.xlsx file matches the version used for the published analysis. ## ============================================================================= cat("\n\n=== PIPELINE COMPLETE ===\n") cat("All tables saved as CSV files and consolidated in Paper3_All_Results.xlsx\n") cat("All figures saved as PNG files in your working directory:", getwd(), "\n")