Supply Chain Analytics
for Health Programmes

library(tidyverse)
library(readxl)

# Import the dataset
data <- read_excel("mpumalanga_stock_data.xlsx")

# Inspect immediately after import
glimpse(data)          # column names, types, first values
nrow(data)             # expected: 24 facilities × 6 months × 8 products
colSums(is.na(data))   # check for missing values by column

# Add MOS and risk flag to every row
data <- data |>
  mutate(
    MOS = ifelse(AMC > 0, SOH / AMC, NA),
    risk_flag = case_when(
      MOS < 2  ~ "HIGH RISK",
      MOS < 3  ~ "WATCH",
      MOS > 6  ~ "OVERSTOCK",
      TRUE     ~ "OK"
    )
  )

# Stockout rate by product — most recent month only
stockout_rate <- data |>
  filter(month == max(month)) |>
  group_by(product) |>
  summarise(
    n_facilities   = n(),
    n_high_risk    = sum(risk_flag == "HIGH RISK", na.rm = TRUE),
    pct_high_risk  = round(n_high_risk / n_facilities * 100, 1),
    avg_MOS        = round(mean(MOS, na.rm = TRUE), 2)
  ) |>
  arrange(desc(pct_high_risk))

library(ggplot2)

# Monthly average MOS for TDF/3TC/DTG across all facilities
arv_trend <- data |>
  filter(product == "TDF/3TC/DTG") |>
  group_by(month) |>
  summarise(avg_MOS = mean(MOS, na.rm = TRUE))

ggplot(arv_trend, aes(x = month, y = avg_MOS)) +
  geom_line(colour = "#0f766e", linewidth = 1.2) +
  geom_point(colour = "#0f766e", size = 3) +
  geom_hline(yintercept = 2,
    colour = "#dc2626", linetype = "dashed") +
  annotate("text", x = min(arv_trend$month), y = 2.15,
    label = "2-month minimum", colour = "#dc2626",
    size = 3.2, hjust = 0) +
  labs(
    title = "TDF/3TC/DTG — Average Months of Stock",
    subtitle = "All 24 facilities, Mpumalanga district",
    x = NULL, y = "Months of stock"
  ) +
  theme_minimal(base_size = 12)

library(writexl)

# Export multiple outputs to one workbook
write_xlsx(
  list(
    "KPI Summary"      = stockout_rate,
    "MOS by Facility"  = data,
    "Risk Flags"       = filter(data, risk_flag == "HIGH RISK")
  ),
  path = "mpumalanga_analysis_output.xlsx"
)

# R Markdown: render report for one facility/district
# In your .Rmd YAML header:
#   params:
#     facility: "Nhlazatshe CHC"
# Then filter data with: filter(data, facility == params$facility)
# Knit once per facility in a loop:

purrr::walk(unique(data$facility), ~
  rmarkdown::render("facility_report.Rmd",
    params = list(facility = .x),
    output_file = paste0(.x, "_report.docx")
  )
)

# 1. Completeness — % of facility-months reporting
completeness <- data |>
  group_by(facility, month) |>
  summarise(n_products = n(), .groups = "drop") |>
  mutate(complete = n_products >= 8)
cat("Reporting completeness:", round(mean(completeness$complete)*100,1), "%\n")

# 2. Range check — negative values
negatives <- data |> filter(SOH < 0 | consumption < 0)
cat("Records with negative values:", nrow(negatives), "\n")

# 3. Consistency — stock accounting equation
data <- data |>
  mutate(
    calc_cb = opening_balance + receipts - consumption - losses,
    discrepancy = abs(calc_cb - closing_balance),
    dq_flag = ifelse(discrepancy > 0, "Flag", "OK")
  )
cat("Records with accounting discrepancy:",
  sum(data$dq_flag == "Flag", na.rm=TRUE), "\n")