library(BatchExperiments)
library(EpiGLM)
library(plyr)
load("simStudy.RData")

res_df <- reduceResultsExperiments(reg, fun=est2df)
is_covered <- function(x, y, z) x < y && y < z
res_df$is_covered <- with(res_df, mapply(is_covered, x=lower, y=ip, z=upper))

tmpf <- function(x) {x <- mean(x, na.rm=TRUE); x*(1-x)}
k <- 10
res_df$index <- rep(1:k, length.out=nrow(res_df))
tmpff <- function(hats, index, prob){
  splits <- split(hats, index)
  x <- sapply(splits, quantile, prob=prob, na.rm=TRUE)
  sd(x, na.rm=TRUE) / sqrt(length(na.omit(x)))
}
  
res_tab <- ddply(res_df, c("nreps", "ip", "ip_dispersion", "floor_times", "propFreeImp"),
  summarize, bias=mean(hat_scale - ip),
  bias_se=sd(hat_scale)/sqrt(length(hat_scale)),
  avCIwidth=mean(width, na.rm=TRUE),
  avCIwidth_se=sd(width, na.rm=TRUE)/sqrt(length(na.omit(width))),
  Pcoverage=mean(is_covered, na.rm=TRUE),
  Pcoverage_se=tmpf(is_covered) / sqrt(length(is_covered)),
  bootLower=quantile(hat_scale, .025), bootUpper=quantile(hat_scale, .975),
  bootLower_se=tmpff(hat_scale, index, .025),
  bootUpper_se=tmpff(hat_scale, index, .975),
  bias_disp=mean(ip_dispersion_hat - ip_dispersion, na.rm=TRUE),
  bias_disp_se=sd(ip_dispersion_hat, na.rm=TRUE)/sqrt(length(na.omit(ip_dispersion_hat))),
  bootLower_disp=quantile(ip_dispersion_hat, .025, na.rm=TRUE),
  bootUpper_disp=quantile(ip_dispersion_hat, .975, na.rm=TRUE),
  bootLower_se_disp=tmpff(ip_dispersion_hat, index, .025),
  bootUpper_se_disp=tmpff(ip_dispersion_hat, index, .975))

cleanf <- function(x){
  ## source: http://stackoverflow.com/q/5465314/548147
  oldx <- c(FALSE, x[-1]==x[-length(x)])
  ## is the value equal to the previous?
  res <- x
  res[oldx] <- NA
  return(res)
}

T <- res_tab[, c("ip", "ip_dispersion", "nreps", "propFreeImp", "floor_times", 
                 "bias", "avCIwidth", "Pcoverage", "bias_disp",
                 "bootLower_disp", "bootUpper_disp")]
ord <- with(T, order(ip, ip_dispersion, nreps, propFreeImp, floor_times))
T <- T[ord, ]
T$propFreeImp <- cleanf(T$propFreeImp)
T$floor_times <- as.integer(T$floor_times)
T$nreps <- cleanf(T$nreps)
T$ip <- cleanf(T$ip)
T$ip_dispersion <- cleanf(T$ip_dispersion)
colnames(T) <- c("$\\mu$", "$\\rho$", "$n$", "imputed", "rounded", 
                 "$\\textrm{bias}(\\hat{\\mu})$",
                 "$\\textrm{av. s.e.}(\\hat{\\mu})$", "cover. (\\%)",
                 "$\\textrm{bias}(\\hat{\\rho})$", "lower $\\hat{\\rho}$",
                 "upper $\\hat{\\rho}$")
xT <- xtable::xtable(T, digits=c(0, 2, 2, 0, 2, 0, 2, 2,
                          2, 2, 2, 2))
print(xT,sanitize.text.function = function(x){x}, include.rownames=FALSE,
      file="tab1.tex")
save.image(file="simStudy-post.RData")
