library(metafor)
library(jsonlite)
a <- commandArgs(trailingOnly=TRUE)
e <- read.csv(file.path(a[1], 'results/effects.csv'))
early <- subset(e, outcome=='depression' & window=='early')
m <- subset(early, group=='masked')
f <- rma(yi=m$yi, vi=m$vi, method='REML', test='knha')
s <- subset(early, setting=='surgery' & group %in% c('masked','awake'))
g <- rma(yi=s$yi, vi=s$vi, mods=~I(s$group=='awake')+I(s$selected), method='REML', test='knha')
reg <- fromJSON(file.path(a[2], 'data/registry.json'))$registrations
r <- subset(reg, eligible & due)
x <- fromJSON(file.path(a[2], 'data/extraction.json'), simplifyVector=FALSE)
t <- Filter(function(v) v$trial=='T097',x)[[1]]
y <- Filter(function(v) identical(v$measure,'PHQ-9') && v$day==1 && v$form=='mean_sd',t$depression)[[1]]$values
pre <- unlist(y[['Pre-ESK group']]); post <- unlist(y[['Post-ESK group']])
n1 <- as.numeric(pre['n']); n2 <- as.numeric(post['n']); df <- n1+n2-2
sp <- sqrt(((n1-1)*as.numeric(pre['sd'])^2+(n2-1)*as.numeric(post['sd'])^2)/df)
J <- exp(lgamma(df/2)-log(sqrt(df/2))-lgamma((df-1)/2))
h <- J*(as.numeric(pre['mean'])-as.numeric(post['mean']))/sp
v <- 1/n1+1/n2+h^2/(2*(n1+n2))
summary <- list(primary=list(k=f$k,n=sum(m$n),g=unname(coef(f)[1]),lower=f$ci.lb,upper=f$ci.ub),
 timing_regression=list(k=g$k,g=unname(coef(g)[2]),lower=g$ci.lb[2],upper=g$ci.ub[2]),
 registry=list(due=nrow(r),unreported=sum(!r$published & !r$results_posted)),
 within_trial=list(n=n1+n2,g=h,lower=h-1.96*sqrt(v),upper=h+1.96*sqrt(v)),
 early_trials=length(unique(early$parent)),early_failed=length(unique(early$parent[!early$trustworthy])))
stopifnot(f$k==23,sum(m$n)==3202,abs(coef(f)[1]-(-0.61))<0.01,
 g$k==37,abs(coef(g)[2]-0.12)<0.01,nrow(r)==39,sum(!r$published & !r$results_posted)==18,
 abs(h-(-0.72))<0.01,length(unique(early$parent))==51,length(unique(early$parent[!early$trustworthy]))==21)
cat(toJSON(summary,auto_unbox=TRUE,pretty=TRUE,digits=12),'\n')