# Reported viral non-suppression, missed ART doses and interest in
# long-acting injectable treatment among people living with HIV in Uganda
# Short report analysis | Stephen Charles Olet | Version 2.0 | 7 October 2026
#
# Standalone base R. No contributed packages required.
# Run from this folder:  source("ART_Adherence_Analysis.R")
# Or:                    Rscript ART_Adherence_Analysis.R /absolute/path/to/folder
#
# Main question: how common is interest in long-acting injectable ART within
# subgroups defined by reported viral status and missed ART doses?
# The analysis is descriptive. It does not test whether any subgroup is MORE
# interested than another, and age is described, not modelled, because age and
# interview session cannot be separated among PLHIV in this survey.

args <- commandArgs(trailingOnly = TRUE)
root <- if (length(args)) args[1] else getwd()
data_file   <- file.path(root, "data", "lai_art_acceptance_v1.csv")
roster_file <- file.path(root, "data", "parent_232_roster.csv")
out <- file.path(root, "outputs")
dir.create(out, recursive = TRUE, showWarnings = FALSE)
options(stringsAsFactors = FALSE, width = 120)

save_csv <- function(x, name) write.csv(x, file.path(out, paste0(name, ".csv")), row.names = FALSE, na = "")
log_lines <- character()
note <- function(x) { log_lines <<- c(log_lines, x); message(x) }
check <- function(ok, msg) if (!isTRUE(ok)) stop(msg, call. = FALSE)

# ---------------------------------------------------------------- 1. Inputs
raw <- read.csv(data_file, check.names = FALSE, na.strings = c("", "NA"))
roster <- read.csv(roster_file, check.names = FALSE)

needed <- c("instanceID", "Q01", "Q02", "Q03", "Q04", "Q15", "Q16", "Q17",
            "Q18", "Q19", "Q20", "Q29", "Q36")
check(all(needed %in% names(raw)), "Required columns are missing. Use the original coded CSV export.")
check(!anyNA(raw$instanceID) && !anyDuplicated(raw$instanceID), "Missing or duplicate identifiers in the export.")
check(nrow(roster) == 232 && !anyDuplicated(roster$instanceID), "Expected the locked 232-person roster without duplicates.")
check(all(roster$instanceID %in% raw$instanceID), "Some roster identifiers are absent from the export.")

as_code <- function(v, allowed) {
  z <- suppressWarnings(as.numeric(as.character(raw[[v]])))
  check(!any(!is.na(raw[[v]]) & is.na(z)), paste("Non-numeric coding in", v))
  check(all(z[!is.na(z)] %in% allowed), paste("Unexpected code in", v))
  z
}
raw$Q01 <- as_code("Q01", 1:5); raw$Q03 <- as_code("Q03", 0:1); raw$Q04 <- as_code("Q04", 1:2)
raw$Q16 <- as_code("Q16", 1:3); raw$Q18 <- as_code("Q18", 1:3); raw$Q19 <- as_code("Q19", 0:3)
raw$Q29 <- as_code("Q29", 1:5); raw$Q36 <- as_code("Q36", 1:5)
for (v in c("Q02", "Q15", "Q17")) raw[[v]] <- suppressWarnings(as.numeric(as.character(raw[[v]])))
check(all(raw$Q02 >= 15, na.rm = TRUE), "Participant younger than 15 years found.")

# ---------------------------------------------------------------- 2. Variables
raw$in_roster  <- raw$instanceID %in% roster$instanceID
raw$eligible   <- raw$Q01 %in% c(1, 3, 4, 5) & !is.na(raw$Q02) & raw$Q02 >= 15
raw$on_art     <- raw$Q01 == 1
raw$interest   <- as.integer(raw$Q29 >= 4)            # willing or very willing
raw$benefit    <- as.integer(raw$Q36 >= 4)            # agrees injectables would help adherence
raw$viral      <- factor(raw$Q18, 1:3, c("Reported suppressed", "Reported not suppressed", "Not tested"))
raw$missed     <- factor(ifelse(raw$Q19 == 0, "No missed doses", ifelse(raw$Q19 == 1, "1-2 missed doses", "3 or more missed doses")),
                         levels = c("No missed doses", "1-2 missed doses", "3 or more missed doses"))
raw$any_missed <- as.integer(raw$Q19 > 0)
raw$age_group  <- factor(ifelse(raw$Q02 < 25, "15-24 years", "25 years or older"), c("15-24 years", "25 years or older"))
raw$sex        <- factor(raw$Q04, 1:2, c("Male", "Female"))
raw$regimen    <- factor(raw$Q16, 1:3, c("Single-tablet regimen", "Multiple pills", "Not sure"))

challenge_labels <- c("Pill fatigue", "Side effects", "Stigma", "Travel", "Food insecurity", "Forgetfulness", "Other")
answered_q20 <- !is.na(raw$Q20) & nzchar(trimws(raw$Q20))
q20_tokens <- strsplit(ifelse(answered_q20, trimws(raw$Q20), ""), "\\s+")
check(all(unlist(q20_tokens) %in% c("", as.character(1:7))), "Unexpected Q20 code.")
for (k in 1:7) raw[[paste0("challenge_", k)]] <- ifelse(answered_q20, vapply(q20_tokens, function(t) as.integer(as.character(k) %in% t), integer(1)), NA)

main     <- raw[raw$in_roster & raw$on_art, ]
expanded <- raw[raw$eligible & raw$on_art, ]
check(nrow(main) == 116, "Expected 116 PLHIV on ART in the roster.")
cohorts <- list("Main cohort (roster, n=116)" = main, "All eligible PLHIV" = expanded)

# ---------------------------------------------------------------- 3. Participant flow
outside <- raw[!raw$in_roster, ]
save_csv(data.frame(
  step = c("Records in export", "Eligible records", "Records outside the parent roster",
           "  of which eligible PLHIV on ART", "  of which eligible PrEP-service clients",
           "PLHIV on ART in the export", "PLHIV on ART in the parent roster (main cohort)"),
  n = c(nrow(raw), sum(raw$eligible), nrow(outside),
        sum(outside$eligible & outside$on_art), sum(outside$eligible & outside$Q01 %in% 3:5),
        sum(raw$eligible & raw$on_art), nrow(main))), "01_participant_flow")

# ---------------------------------------------------------------- 4. Helpers
wilson <- function(e, n) {
  if (n == 0) return(c(NA, NA, NA))
  z <- qnorm(0.975); p <- e / n; d <- 1 + z^2 / n
  half <- z * sqrt(p * (1 - p) / n + z^2 / (4 * n^2))
  c(p, (p + z^2 / (2 * n) - half) / d, (p + z^2 / (2 * n) + half) / d)
}
prop_row <- function(x, label) {
  n <- sum(!is.na(x)); e <- sum(x == 1, na.rm = TRUE); ci <- wilson(e, n)
  data.frame(subgroup = label, events = e, n = n, percent = round(100 * ci[1], 1),
             lower95 = round(100 * ci[2], 1), upper95 = round(100 * ci[3], 1),
             note = ifelse(n < 10, "Fewer than 10 participants: report the count only", ""))
}
pct <- function(k, n) sprintf("%d (%.1f)", k, 100 * k / n)
fmt <- function(v) format(round(v, 1), nsmall = 0, trim = TRUE)
mid <- function(x) sprintf("%s (%s-%s)", fmt(median(x, na.rm = TRUE)), fmt(quantile(x, .25, na.rm = TRUE)), fmt(quantile(x, .75, na.rm = TRUE)))

subgroups <- function(d) list(
  "All PLHIV on ART"                       = rep(TRUE, nrow(d)),
  "Reported suppressed"                    = d$viral == "Reported suppressed",
  "Reported not suppressed"                = d$viral == "Reported not suppressed",
  "Not tested"                             = d$viral == "Not tested",
  "No missed doses"                        = d$missed == "No missed doses",
  "1-2 missed doses"                       = d$missed == "1-2 missed doses",
  "3 or more missed doses"                 = d$missed == "3 or more missed doses",
  "Suppressed, no missed doses"            = d$viral == "Reported suppressed" & d$any_missed == 0,
  "Suppressed, any missed doses"           = d$viral == "Reported suppressed" & d$any_missed == 1,
  "Not suppressed, no missed doses"        = d$viral == "Reported not suppressed" & d$any_missed == 0,
  "Not suppressed, any missed doses"       = d$viral == "Reported not suppressed" & d$any_missed == 1)

# ---------------------------------------------------------------- 5. Table 1: characteristics
table1 <- function(d) {
  n <- nrow(d)
  rows <- list(
    c("Age in years, median (IQR)", mid(d$Q02)),
    c("Aged 15-24 years, n (%)", pct(sum(d$age_group == "15-24 years"), n)),
    c("Female, n (%)", pct(sum(d$sex == "Female"), n)),
    c("Years since HIV diagnosis, median (IQR)", mid(d$Q15)),
    c("Years on ART, median (IQR)", mid(d$Q17)))
  for (lv in levels(d$regimen)) rows[[length(rows) + 1]] <- c(paste0("Regimen: ", lv, ", n (%)"), pct(sum(d$regimen == lv), n))
  for (lv in levels(d$viral))   rows[[length(rows) + 1]] <- c(paste0("Viral status (past 12 months): ", lv, ", n (%)"), pct(sum(d$viral == lv), n))
  for (lv in levels(d$missed))  rows[[length(rows) + 1]] <- c(paste0("Missed doses (past 30 days): ", lv, ", n (%)"), pct(sum(d$missed == lv), n))
  na <- sum(!is.na(d$challenge_1))
  for (k in 1:7) rows[[length(rows) + 1]] <- c(paste0("Adherence challenge: ", challenge_labels[k], ", n (%)"), pct(sum(d[[paste0("challenge_", k)]] == 1, na.rm = TRUE), na))
  t <- do.call(rbind, rows); data.frame(characteristic = t[, 1], value = t[, 2])
}
t1 <- table1(main); names(t1)[2] <- names(cohorts)[1]
t1[[names(cohorts)[2]]] <- table1(expanded)$value
save_csv(t1, "02_Table1_characteristics")

# ---------------------------------------------------------------- 6. Table 2: interest by subgroup (MAIN RESULT)
interest_table <- function(var, file) {
  res <- do.call(rbind, lapply(names(cohorts), function(cn) {
    d <- cohorts[[cn]]; g <- subgroups(d)
    r <- do.call(rbind, lapply(names(g), function(s) prop_row(d[[var]][g[[s]]], s)))
    r$cohort <- cn; r
  }))
  save_csv(res, file); res
}
t2 <- interest_table("interest", "03_Table2_interest_by_subgroup")  # assignment, not printed
invisible(interest_table("benefit", "04_TableS1_perceived_adherence_benefit_by_subgroup"))

# Full five-level willingness distribution, so that "interest" is not only shown as a cut-point
lv29 <- c("Very unwilling", "Unwilling", "Neutral", "Willing", "Very willing")
dist <- do.call(rbind, lapply(c("viral", "missed"), function(v) {
  tab <- table(main[[v]], factor(main$Q29, 1:5, lv29))
  data.frame(grouping = v, subgroup = rownames(tab), as.data.frame.matrix(tab), n = rowSums(tab), check.names = FALSE)
}))
save_csv(dist, "05_TableS2_full_willingness_distribution")

# ---------------------------------------------------------------- 7. Table 3: overlap of the two clinical measures
overlap <- as.data.frame.matrix(table(main$viral, main$missed))
overlap <- cbind(viral_status = rownames(overlap), overlap, Total = rowSums(overlap))
save_csv(overlap, "06_Table3_viral_status_by_missed_doses")

# ---------------------------------------------------------------- 8. Supplementary descriptive contrasts
# Differences in proportions with Newcombe hybrid-score intervals and Fisher exact p-values.
# These show the size of any difference and its uncertainty. They are not tests of a
# hypothesis that one subgroup is more interested.
newcombe <- function(e1, n1, e2, n2) {
  a <- wilson(e1, n1); b <- wilson(e2, n2); diff <- a[1] - b[1]
  c(diff, diff - sqrt((a[1] - a[2])^2 + (b[3] - b[1])^2), diff + sqrt((a[3] - a[1])^2 + (b[1] - b[2])^2))
}
contrast <- function(d, cn, label, g1, g0) {
  e1 <- sum(d$interest[g1]); n1 <- sum(g1); e0 <- sum(d$interest[g0]); n0 <- sum(g0)
  nc <- newcombe(e1, n1, e0, n0)
  data.frame(cohort = cn, contrast = label, group = paste0(e1, "/", n1), comparison = paste0(e0, "/", n0),
             difference_pp = round(100 * nc[1], 1), lower95 = round(100 * nc[2], 1), upper95 = round(100 * nc[3], 1),
             fisher_p = signif(fisher.test(matrix(c(e1, n1 - e1, e0, n0 - e0), 2))$p.value, 3))
}
save_csv(do.call(rbind, lapply(names(cohorts), function(cn) { d <- cohorts[[cn]]; rbind(
  contrast(d, cn, "Reported not suppressed minus reported suppressed", d$viral == "Reported not suppressed", d$viral == "Reported suppressed"),
  contrast(d, cn, "Any missed doses minus no missed doses", d$any_missed == 1, d$any_missed == 0)) })),
  "07_TableS3_descriptive_differences")

# ---------------------------------------------------------------- 9. Retained model (HC0): overlap check only
# Modified Poisson for reported non-suppression by any missed doses, among those tested.
# Describes how closely the two self-reports agree. It is not a causal model.
hc0_pr <- function(d, cn) {
  d <- d[d$viral != "Not tested", ]; y <- as.integer(d$viral == "Reported not suppressed")
  X <- cbind(1, d$any_missed); fit <- glm.fit(X, y, family = poisson(link = "log"))
  check(fit$converged, "Overlap model did not converge.")
  mu <- fit$fitted.values; bread <- solve(crossprod(X, X * mu)); U <- X * (y - mu)
  V <- bread %*% crossprod(U) %*% bread; b <- fit$coefficients[2]; se <- sqrt(V[2, 2])
  data.frame(cohort = cn, model = "Reported non-suppression by any missed doses (tested participants)",
             n = length(y), events = sum(y), PR = round(exp(b), 2), lower95 = round(exp(b - 1.96 * se), 2),
             upper95 = round(exp(b + 1.96 * se), 2), variance = "HC0 robust", max_fitted = round(max(mu), 3))
}
save_csv(do.call(rbind, lapply(names(cohorts), function(cn) hc0_pr(cohorts[[cn]], cn))), "08_TableS4_overlap_model_HC0")

# ---------------------------------------------------------------- 10. Age, descriptive only
age_rows <- do.call(rbind, lapply(names(subgroups(main))[1:7], function(s) {
  g <- subgroups(main)[[s]]; x <- main[g, ]
  data.frame(subgroup = s, n = nrow(x), median_age = if (nrow(x)) median(x$Q02) else NA,
             aged_15_24 = if (nrow(x)) pct(sum(x$age_group == "15-24 years"), nrow(x)) else "")
}))
save_csv(age_rows, "09_TableS5_age_by_subgroup_descriptive")

# ---------------------------------------------------------------- 11. Figure 1
fig_data <- t2[t2$cohort == names(cohorts)[1] & t2$n >= 10 & t2$subgroup %in% names(subgroups(main))[1:7], ]
draw_fig <- function() {
  par(mar = c(5, 15, 1, 6)); y <- rev(seq_len(nrow(fig_data)))
  plot(fig_data$percent, y, xlim = c(0, 100), ylim = c(0.5, nrow(fig_data) + 0.5), yaxt = "n", ylab = "",
       xlab = "Participants willing or very willing to use injectable ART (%)", pch = 19, bty = "l")
  segments(fig_data$lower95, y, fig_data$upper95, y)
  axis(2, at = y, labels = fig_data$subgroup, las = 1, cex.axis = 0.85)
  axis(4, at = y, labels = paste0(fig_data$events, "/", fig_data$n), las = 1, tick = FALSE, cex.axis = 0.8)
}
pdf(file.path(out, "Figure1_interest_by_subgroup.pdf"), width = 8, height = 4.5, family = "sans"); draw_fig(); invisible(dev.off())
png(file.path(out, "Figure1_interest_by_subgroup.png"), width = 2400, height = 1350, res = 300); draw_fig(); invisible(dev.off())

# ---------------------------------------------------------------- 12. Notes and provenance
writeLines(c(
  "Main question: how common is interest in long-acting injectable ART within subgroups defined by reported viral status and missed ART doses?",
  "Interest means willing or very willing (Q29 = 4 or 5). The full five-level distribution is in Table S2.",
  "Main cohort: the 116 PLHIV on ART in the locked 232-person parent roster. Sensitivity cohort: all eligible PLHIV on ART in the export.",
  "Participants outside the roster were eligible and were removed only by the parent study's equal-group selection; confirm the PLHIV selection method with the study team.",
  "Viral status is self-reported for the past 12 months, with no laboratory value or defined threshold. Missed doses refer to the past 30 days.",
  "Subgroups with fewer than 10 participants are reported as counts only and are left out of Figure 1.",
  "Results are descriptive. Do not describe any subgroup as more or less interested; Table S3 gives differences with intervals for transparency.",
  "Age is described only. Among PLHIV, younger and older participants were interviewed in separate sessions.",
  "The only model (Table S4) checks agreement between the two self-reported clinical measures and uses HC0 robust variance, matching the other papers."),
  file.path(out, "00_methods_notes.txt"))
writeLines(capture.output(sessionInfo()), file.path(out, "sessionInfo.txt"))
save_csv(data.frame(file = basename(c(data_file, roster_file)), md5 = unname(tools::md5sum(c(data_file, roster_file)))), "input_fingerprints")
writeLines(log_lines, file.path(out, "analysis_log.txt"))
note(paste("Analysis complete. Outputs in", out))