# Optional independent validation, run from the bundle root with metafor/jsonlite.
# This script reselects source records and recombines arms in base R.
# It creates data/metafor-reference.json for the Python runner's cross-checks.
library(metafor)
library(jsonlite)
Sys.setlocale('LC_COLLATE', 'C')
d <- read.csv('data/source.csv', check.names=FALSE, stringsAsFactors=FALSE)
walk <- 'Walking / Jogging'
controls <- c('Educational','Social','Social or educational control',
              'Usual care','Placebo pill','Stretching')
combine <- function(x) {
  total <- sum(x$n)
  m <- weighted.mean(x$mean, x$n)
  s <- sqrt((sum((x$n-1)*x$sd^2) + sum(x$n*(x$mean-m)^2))/(total-1))
  c(n=total,mean=m,sd=s)
}
select <- function(control_labels=controls, clinician=TRUE) {
  all <- list()
  for (study in sort(unique(d$studyID))) {
    x <- d[d$studyID==study & d$weeks_from_end_of_treatment_to_measurement==0 &
           (d$trt==walk | d$trt %in% control_labels),]
    if (!any(x$trt==walk) || !any(x$trt %in% control_labels)) next
    stopifnot(!anyDuplicated(x[,c('arm_number','outcome')]))
    common <- Reduce(intersect, split(x$outcome,x$arm_number))
    if (!length(common)) next
    common <- sort(common)
    first <- common[1]
    if (clinician) {
      rated <- common[vapply(common,function(m) all(x$reported[x$outcome==m]=='Clinician'),logical(1))]
      if (length(rated)) first <- rated[1]
    }
    x <- x[x$outcome==first,]
    a <- combine(x[x$trt==walk,]); b <- combine(x[x$trt %in% control_labels,])
    es <- escalc(measure='SMD', m1i=a['mean'],sd1i=a['sd'],n1i=a['n'],
                 m2i=b['mean'],sd2i=b['sd'],n2i=b['n'],vtype='LS')
    all[[study]] <- list(studyID=study,outcome=first,yi=as.numeric(es$yi),vi=as.numeric(es$vi),
                        risk_random=all(x$random_sequence_generation_selection_bias=='Low risk'),
                        risk_allocation=all(x$allocation_concealment_selection_bias=='Low risk'),
                        risk_assessor=all(x$blinded_outcome_assessor=='Low risk'))
  }
  all
}
fit <- function(items,method='REML',test='adhoc') {
  if (length(items)<2) return(list(k=length(items),status='insufficient'))
  model <- rma.uni(yi=vapply(items,function(x) x$yi,numeric(1)),
                   vi=vapply(items,function(x) x$vi,numeric(1)),method=method,test=test,
                   control=list(threshold=1e-10,maxiter=1000,stepadj=.5))
  list(k=length(items),g=as.numeric(model$b),se_g=model$se,
       ci95_low_g=model$ci.lb,ci95_high_g=model$ci.ub,tau2=model$tau2,
       Q=model$QE,Q_p_value=model$QEp,p_two_sided=model$pval)
}
selected <- select()
reference <- list(
  software=list(R=as.character(getRversion()),metafor=as.character(packageVersion('metafor')),
                jsonlite=as.character(packageVersion('jsonlite'))),
  studies=lapply(selected,function(x) x[c('studyID','outcome','yi','vi')]),
  primary=fit(selected),
  sensitivity=list(REML_normal=fit(selected,test='z'),REML_unmodified_HK=fit(selected,test='knha'),
                   fixed_normal=fit(selected,method='FE',test='z'),DL_normal=fit(selected,method='DL',test='z'),
                   low_randomization_and_allocation=fit(Filter(function(x) x$risk_random && x$risk_allocation,selected)),
                   low_blinded_assessor=fit(Filter(function(x) x$risk_assessor,selected)),
                   lexicographic_measure=fit(select(clinician=FALSE)),usual_care_only=fit(select('Usual care'))),
  leave_one_out=lapply(names(selected),function(name) c(list(omitted=name),fit(selected[names(selected)!=name])))
)
write_json(reference,'data/metafor-reference.json',pretty=TRUE,auto_unbox=TRUE,digits=NA)
cat('Independent reference generated using R',reference$software$R,'and metafor',reference$software$metafor,'\n')