# Manuscript 1: Oral PrEP history and stated interest in injectable prevention.
# Base R only. Run from directory containing lai_art_acceptance_v1.csv:
# source('01_PrEP_history_revised.R')
# Or Rscript 01_PrEP_history_revised.R path/to/raw.csv path/to/output

args <- commandArgs(trailingOnly=TRUE)
input <- if(length(args)>=1) args[1] else 'lai_art_acceptance_v1.csv'
out <- if(length(args)>=2) args[2] else 'Manuscript_1_outputs'
dir.create(out,recursive=TRUE,showWarnings=FALSE)
if(!file.exists(input)) stop('Input CSV not found: ',input)
d <- read.csv(input,check.names=FALSE,stringsAsFactors=FALSE,na.strings=c('','NA'))
required <- c('instanceID','Q01','Q02','Q04','Q05','Q06','Q08','Q09','Q10','Q21','Q25','Q29','Q30','Q33')
if(length(setdiff(required,names(d)))) stop('Missing columns: ',paste(setdiff(required,names(d)),collapse=', '))
if(anyDuplicated(d$instanceID)) stop('Duplicate instanceID in raw file')
num <- function(x) suppressWarnings(as.numeric(x))
save <- function(x,n) write.csv(x,file.path(out,paste0(n,'.csv')),row.names=FALSE,na='')
s <- d[d$Q01 %in% c('3','4','5') & !is.na(num(d$Q02)) & num(d$Q02)>=15 & num(d$Q29) %in% 1:5,,drop=FALSE]
s$history <- factor(s$Q01,levels=c('3','4','5'),labels=c('Current oral PrEP','Previous oral PrEP','Never used oral PrEP'))
s$age <- num(s$Q02)
s$age_decade <- s$age/10
s$sex <- factor(s$Q04,levels=c('1','2'),labels=c('Male','Female'))
s$willing_score <- num(s$Q29)
s$willing <- as.integer(s$willing_score>=4)
stopifnot(nrow(s)==120L, all(table(s$history)==c(56,34,30))) # Study-specific audit; revise only after documented raw-data change.
save(data.frame(stage=c('Raw','PrEP service eligible','Current','Previous','Never'),n=c(nrow(d),nrow(s),as.integer(table(s$history)))),'cohort_flow')
save(data.frame(variable=names(s),missing=vapply(s,function(x) sum(is.na(x)),integer(1))),'missingness_prep_population')

# Wilson binomial intervals, stated willingness is scores 4-5 of 5.
wilson <- function(k,n,z=qnorm(.975)) {
  p<-k/n; den<-1+z*z/n
  c(lower=(p+z*z/(2*n)-z*sqrt((p*(1-p)+z*z/(4*n))/n))/den,
    upper=(p+z*z/(2*n)+z*sqrt((p*(1-p)+z*z/(4*n))/n))/den)
}
groups<-levels(s$history)
summ<-do.call(rbind,lapply(groups,function(g) {
  x<-s$willing[s$history==g]; ci<-wilson(sum(x),length(x))
  data.frame(group=g,n=length(x),willing=sum(x),percent=100*mean(x),lower=100*ci[1],upper=100*ci[2])
}))
save(summ,'Table_1_primary_willingness_with_95CI')
score<-as.data.frame.matrix(table(s$history,factor(s$willing_score,levels=1:5)))
save(data.frame(group=rownames(score),score,check.names=FALSE),'Table_2_five_level_willingness_counts')
save(as.data.frame(table(history=s$history,sex=s$sex,useNA='ifany')),'Table_S1_sex_by_history')
save(data.frame(history=groups,mean_age=sapply(groups,function(g)mean(s$age[s$history==g])),sd_age=sapply(groups,function(g)sd(s$age[s$history==g]))),'Table_S1_age_by_history')
save(as.data.frame(table(history=s$history,Q21_code=s$Q21,useNA='ifany')),'Table_S2_Q21_status_codes')
save(as.data.frame(table(history=s$history,timing_code=s$Q30,useNA='ifany')),'Table_S3_timing_codes')
save(as.data.frame(table(history=s$history,care_mode_code=s$Q33,useNA='ifany')),'Table_S4_care_mode_codes')
baseline <- function(v,label) {
  x<-ifelse(is.na(v),'Missing',as.character(v))
  as.data.frame(table(variable=rep(label,length(x)),history=s$history,category=x),stringsAsFactors=FALSE)
}
baseline_table<-do.call(rbind,list(
  baseline(s$sex,'Sex'),baseline(s$Q05,'Marital status code'),
  baseline(s$Q06,'Education code'),baseline(s$Q08,'Employment code'),
  baseline(s$Q10,'Residence code'),baseline(s$Q25,'Prior awareness code')))
save(baseline_table,'Table_0_characteristics_by_history_codes')
save(data.frame(history=groups,
  median_age=sapply(groups,function(g)median(s$age[s$history==g])),
  age_IQR=sapply(groups,function(g)IQR(s$age[s$history==g])),
  median_income_UGX=sapply(groups,function(g)median(num(s$Q09[s$history==g]),na.rm=TRUE))),
  'Table_0_age_income_summaries')

# Crude and adjusted modified Poisson models with HC0 robust covariance.
model_data<-s[complete.cases(s[,c('willing','history','age_decade','sex')]),,drop=FALSE]
save(data.frame(n=nrow(model_data),events=sum(model_data$willing),excluded=nrow(s)-nrow(model_data)),'model_sample')
fit <- function(formula,dat) {
  m<-glm(formula,data=dat,family=poisson(link='log'))
  X<-model.matrix(m); u<-as.vector(residuals(m,type='response'))
  B<-solve(crossprod(X,m$fitted.values*X)); S<-crossprod(X*u)
  V<-B%*%S%*%B
  se<-sqrt(diag(V)); if(any(!is.finite(se)|se<=0)) stop('Invalid robust standard errors')
  b<-coef(m)
  data.frame(term=names(b),PR=exp(b),lower=exp(b-1.96*se),upper=exp(b+1.96*se),p=2*pnorm(-abs(b/se)),row.names=NULL)
}
crude<-fit(willing~history,model_data)
adjusted<-fit(willing~history+age_decade+sex,model_data)
save(crude,'Table_3_crude_PR_HC0')
save(adjusted,'Table_3_adjusted_PR_HC0')
if(any(adjusted$lower[grepl('history',adjusted$term)]<=0)) stop('Invalid confidence interval')
diagnostic_model<-glm(willing~history+age_decade+sex,data=model_data,family=poisson(link='log'))
design<-model.matrix(diagnostic_model)
save(data.frame(n=nrow(model_data),events=sum(model_data$willing),
  design_rank=qr(design)$rank,design_columns=ncol(design),
  max_fitted_mean=max(fitted(diagnostic_model)),
  max_Cooks_distance=max(cooks.distance(diagnostic_model)),
  influential_Cooks_over_4_per_n=sum(cooks.distance(diagnostic_model)>4/nrow(model_data))),
  'model_diagnostics')
if(qr(design)$rank<ncol(design)) stop('Design matrix is rank deficient')
if(max(fitted(diagnostic_model))>1) warning('Poisson model predicts a mean above one for some observations; inspect model diagnostics')

# Absolute contrasts and exact tests complement relative estimates.
cont<-do.call(rbind,lapply(c('Previous oral PrEP','Never used oral PrEP'),function(g){
  x<-s$willing[s$history=='Current oral PrEP'];y<-s$willing[s$history==g]
  mat<-matrix(c(sum(x),length(x)-sum(x),sum(y),length(y)-sum(y)),nrow=2,byrow=TRUE)
  data.frame(comparison=paste(g,'versus current'),difference_percentage_points=100*(mean(y)-mean(x)),fisher_p=fisher.test(mat)$p.value)
}))
save(cont,'Table_S5_absolute_differences_exploratory')

# Stratified nonparametric bootstrap for percentage-point differences.
set.seed(20260924)
boot<-do.call(rbind,lapply(c('Previous oral PrEP','Never used oral PrEP'),function(g){
  x<-s$willing[s$history=='Current oral PrEP'];y<-s$willing[s$history==g]
  draws<-replicate(3000,mean(sample(y,length(y),replace=TRUE))-mean(sample(x,length(x),replace=TRUE)))
  data.frame(comparison=paste(g,'minus current'),
    difference_percentage_points=100*(mean(y)-mean(x)),
    bootstrap_lower=100*unname(quantile(draws,.025)),
    bootstrap_upper=100*unname(quantile(draws,.975)))
}))
save(boot,'Table_S5_bootstrap_absolute_differences')
tests<-data.frame(test=c('Pearson chi-square 3 by 2 willingness','Kruskal-Wallis five-level willingness'),
  p=c(chisq.test(table(s$history,s$willing),correct=FALSE)$p.value,
      kruskal.test(willing_score~history,data=s)$p.value))
save(tests,'Table_S7_exploratory_global_tests')

# Original five-level willingness: descriptive sensitivity, avoiding a new arbitrary cutpoint.
score_prop<-prop.table(table(s$history,factor(s$willing_score,levels=1:5)),1)
save(data.frame(group=rownames(score_prop),100*score_prop,check.names=FALSE),'Table_S6_five_level_percent')

navy<-'#152D45'; teal<-'#127A86'; orange<-'#D08345'; muted<-'#738398'
cols<-c(teal,orange,navy)
render <- function(stem,plotfun,width=9,height=6) {
  png(file.path(out,paste0(stem,'.png')),width=width,height=height,units='in',res=300,bg='white');plotfun();dev.off()
  pdf(file.path(out,paste0(stem,'.pdf')),width=width,height=height,useDingbats=FALSE);plotfun();dev.off()
}
render('Figure_1_Willingness_by_PrEP_history',function(){
  par(mar=c(5.3,13,3.4,2),family='sans',fg=navy)
  y<-3:1;plot(NA,xlim=c(0,125),ylim=c(.5,3.5),axes=FALSE,xlab='',ylab='',main='Interest in injectable PrEP by oral PrEP history')
  abline(v=seq(0,100,25),col='#E5EAF0',lwd=1)
  for(i in 1:3){rect(0,y[i]-.22,summ$percent[i],y[i]+.22,col=cols[i],border=NA)
    segments(summ$lower[i],y[i],summ$upper[i],y[i],lwd=2,col=navy)
    segments(summ$lower[i],y[i]-.10,summ$lower[i],y[i]+.10,col=navy,lwd=2)
    segments(summ$upper[i],y[i]-.10,summ$upper[i],y[i]+.10,col=navy,lwd=2)
    text(summ$upper[i]+2,y[i],sprintf('%d/%d  %.1f%%',summ$willing[i],summ$n[i],summ$percent[i]),adj=0,cex=.83,col=navy)
  }
  axis(1,at=seq(0,100,25),labels=paste0(seq(0,100,25),'%'),col.axis=navy)
  axis(2,at=y,labels=groups,las=1,tick=FALSE,col.axis=navy)
  mtext('Bars show percentages; whiskers show Wilson 95% confidence intervals',side=1,line=3.3,cex=.78,col=muted)
  box(bty='l',col='#BEC8D1')
},width=11.5,height=6.0)

render('Figure_2_Five_level_willingness',function(){
  layout(matrix(c(1,2),ncol=1),heights=c(5,1))
  par(mar=c(4.8,13,3.5,2),family='sans',fg=navy)
  scorecols<-c('#B35D59','#DBA194','#E4E9EE','#72AEB1','#167D89')
  vals<-as.matrix(score_prop[rev(groups),,drop=FALSE])*100
  barplot(t(vals),horiz=TRUE,beside=FALSE,col=scorecols,border=NA,space=.45,xlim=c(0,100),axes=FALSE,names.arg=rev(groups),las=1,main='Full distribution of willingness')
  axis(1,at=seq(0,100,25),labels=paste0(seq(0,100,25),'%'),col.axis=navy)
  mtext('Responses to the five-level willingness question',side=1,line=3.2,cex=.8,col=muted)
  box(bty='l',col='#BEC8D1')
  par(mar=c(0,0,0,0));plot.new()
  legend('center',legend=c('Very unwilling','Unwilling','Neutral','Willing','Very willing'),fill=scorecols,bty='n',cex=.76,horiz=TRUE)
  layout(1)
},width=11.5,height=6.8)

render('Figure_3_Adjusted_prevalence_ratios',function(){
  p<-adjusted[grepl('^history',adjusted$term),,drop=FALSE]
  yy<-2:1;par(mar=c(5.5,16,3.3,2.5),family='sans',fg=navy)
  plot(NA,xlim=c(.2,1.7),ylim=c(.45,2.6),axes=FALSE,xlab='',ylab='',main='Adjusted prevalence of willingness')
  abline(v=1,lty=2,col='#A3AEB8');axis(1,at=seq(.2,1.2,.2),col.axis=navy)
  axis(2,at=yy,labels=c('Previous vs current oral PrEP','Never vs current oral PrEP'),las=1,tick=FALSE,col.axis=navy)
  for(i in seq_len(nrow(p))){segments(p$lower[i],yy[i],p$upper[i],yy[i],lwd=2,col=navy)
    points(p$PR[i],yy[i],pch=21,bg=cols[i+1],col=navy,cex=1.7)
    text(1.26,yy[i],sprintf('%.2f (%.2f–%.2f)',p$PR[i],p$lower[i],p$upper[i]),adj=0,cex=.8,col=navy)
  }
  mtext('Prevalence ratio, adjusted for age and sex; robust 95% confidence interval',side=1,line=3.5,cex=.78,col=muted)
  box(bty='l',col='#BEC8D1')
},width=12.5,height=6.0)

writeLines(c('Manuscript 1: oral PrEP history among 120 prevention-service respondents.',
  'Primary outcome: Q29 scores 4-5; comparison is hypothetical willingness, not actual uptake.',
  'Three treatment-history categories reflect Q01=3 current, Q01=4 previous, Q01=5 never oral PrEP.',
  'The adjusted model contains age per decade and sex. HC0 sandwich standard errors.',
  'Q21, Q30 and Q33 code labels require questionnaire verification before publication.',
  'Table S5 tests are exploratory; avoid causal interpretation or claims based on isolated p values.'),file.path(out,'METHODS_AND_INTERPRETATION.txt'))
cat('Finished Manuscript 1 analysis:',nrow(s),'participants. Outputs in',normalizePath(out),'\n')