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

reg$file.dir <- file.path(getwd(), "simStudy2-files")


#resL <- reduceResultsList(reg)

rre <- function(...){
  # Character vectors are returned separately to prevent coercion to
  # factors
  foo <- reduceResultsExperiments(...)
  bar <- reduceResultsExperiments(..., returnch=TRUE)
  cbind(foo, bar)
}
#res_df_ad <- rre(reg, fun=trest2df, method="admb")
#res_df_uni <- rre(reg, fun=trest2df, method="uniroot")
#res_df_tr <- rbind(res_df_ad, res_df_uni)
system.time(res_df_tr <- rre(reg, fun=trest2df, method="uniroot"))

system.time(res_df_ip <- reduceResultsExperiments(reg, fun=ipest2df))

is_covered <- function(x, y, z) x < y && y < z
res_df_ip$is_covered <- with(res_df_ip,
                             mapply(is_covered, x=lower, y=ip, z=upper))

res_df_tr$is_covered_gr <- with(res_df_tr,
                                mapply(is_covered, x=hat_cgr_lower,
                                       y=gr, z=hat_cgr_upper))
res_df_tr$is_covered_tr <- with(res_df_tr,
                                mapply(is_covered, x=hat_ctr_lower,
                                       y=tr, z=hat_ctr_upper))

tmpf <- function(x) {x <- mean(x, na.rm=TRUE); x*(1-x)}
k <- 10
res_df_ip$index <- rep(1:k, length.out=nrow(res_df_ip))
res_df_tr$index <- rep(1:k, length.out=nrow(res_df_tr))
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)))
}

test <- res_df_ip$propOnsetNA == 0
res_tab_ip <- ddply(res_df_ip[test, ],
                    c("nreps", "ip", "ip_dispersion",
                      "propFreeImp", "floor_times"),
  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(na.omit(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))

res_tab_tr <- ddply(res_df_tr, c("nreps", "tr", "gr", "floor_times", "propFreeImp", "propOnsetNA"),
  summarize, bias_tr=mean(hat_ctr - tr, na.rm=TRUE),
  bias_gr=mean(hat_cgr - gr, na.rm=TRUE),
  bias_tr_se=sd(hat_ctr, na.rm=TRUE)/sqrt(length(na.omit(hat_ctr))),
  bias_gr_se=sd(hat_cgr, na.rm=TRUE)/sqrt(length(na.omit(hat_cgr))),
  avCIwidth_tr=mean(hat_ctr_width, na.rm=TRUE),
  avCIwidth_gr=mean(hat_cgr_width, na.rm=TRUE),
  avCIwidth_se_tr=sd(hat_ctr_width, na.rm=TRUE)/sqrt(length(na.omit(hat_ctr_width))),
  avCIwidth_se_gr=sd(hat_cgr_width, na.rm=TRUE)/sqrt(length(na.omit(hat_cgr_width))),
  Pcoverage_tr=mean(is_covered_tr, na.rm=TRUE),                    
  Pcoverage_se_tr=tmpf(is_covered_tr) / sqrt(length(na.omit(is_covered_tr))),
  Pcoverage_gr=mean(is_covered_gr, na.rm=TRUE),                    
  Pcoverage_se_gr=tmpf(is_covered_gr) / sqrt(length(na.omit(is_covered_gr))),
  bootLower_tr=quantile(hat_ctr, .025, na.rm=TRUE),
  bootUpper_tr=quantile(hat_ctr, .975, na.rm=TRUE),
  bootLower_gr=quantile(hat_cgr, .025, na.rm=TRUE),
  bootUpper_gr=quantile(hat_cgr, .975, na.rm=TRUE),
  bootLower_se_tr=tmpff(hat_ctr, index, .025),
  bootUpper_se_tr=tmpff(hat_ctr, index, .975),
  bootLower_se_gr=tmpff(hat_cgr, index, .025),
  bootUpper_se_gr=tmpff(hat_cgr, index, .975),
  meancor=mean(rho_ctr_gr, na.rm=TRUE),
  meancor_se=sd(rho_ctr_gr, na.rm=TRUE)/sqrt(length(na.omit(rho_ctr_gr))),
  Paderr=mean(note=="admb_error"),
  Pnodat=mean(note=="no_transmission_data"),
  Pdiverg=mean(note=="divergent"),
  Puniden=mean(note=="unidentifiable"),
  Pexpzer=mean(note=="exposure_rounded_to_zero"),
  Pminx=mean(note=="minimum_feasible_X0"))

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)
}

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

Ttr <- res_tab_tr[, c("tr", "gr", "nreps", "propFreeImp", "propOnsetNA",
                      "floor_times", "bias_tr", "avCIwidth_tr", "Pcoverage_tr",
                      "bias_gr", "avCIwidth_gr", "Pcoverage_gr", "meancor",
                      "Paderr", "Pnodat", "Pdiverg", "Puniden", "Pexpzer",
                      "Pminx")]
ord <- with(Ttr, order(tr, gr, nreps, propFreeImp, propOnsetNA, floor_times))
Ttr <- Ttr[ord, ]

Ttr$floor_times <- as.integer(Ttr$floor_times)
Ttr$nreps <- cleanf(Ttr$nreps)
Ttr$tr <- cleanf(Ttr$tr)
Ttr$gr <- cleanf(Ttr$gr)
Ttr$propFreeImp <- cleanf(Ttr$propFreeImp)
Ttr$propOnsetNA <- cleanf(Ttr$propOnsetNA)

colnames(Ttr)[1:13] <- c("$\\beta$", "$r$", "$n$", "imputed", "missing",
                         "rounded", "$\\textrm{bias}(\\hat{\\beta})$",
                         "$\\textrm{av. s.e.}(\\hat{\\beta})$",
                         "$\\beta$ cover. (\\%)", "$\\textrm{bias}(\\hat{r})$",
                         "$\\textrm{av. s.e.}(\\hat{r})$", "$r$ cover. (\\%)",
                         "av. corr.")
xTtr <- xtable::xtable(Ttr, digits=c(0, 4, 2, 0, 2, 2, 0, 4, 4, 2, 2,
                              2, 2, 2, 2, 2, 2, 2, 2, 2))
print(xTtr,sanitize.text.function = function(x){x}, include.rownames=FALSE,
      file="tab_tr.tex")

save.image(file="simStudy-post.RData")
