library(BatchExperiments)
library(EpiGLM)
library(plyr)

data(estData)
load("estBoot.RData")
resL <- reduceResultsList(reg, ids=1)
#fit1ll <- sapply(resL, function(x) {b <- x$fit$admb$loglik
#                                    if(!is.null(b)) b else NA})
#F <- ecdf(fit1ll)
#F(estData$fit$admb$loglik)

resdf <- reduceResultsExperiments(reg, fun=est2df)
test <- !is.na(resdf$hat_ctr.1)
resdf <- resdf[test,]

theta <- c(estData$fit$admb$coefficients, coef(estData$ip_fit$m))              
names(theta) <- paste("theta", names(theta), sep="_")
names(theta) <- make.names(names(theta))

resdf <- cbind(t(theta), resdf)

is_covered <- function(x, y, z) x < y && y < z
parnames <- names(theta)
for(parname in parnames){
  y <- resdf[, parname]
  hatname <- gsub("theta", "hat", parname)
  lname <- paste(hatname, "lower", sep="_")
  x <- resdf[, lname]
  uname <- paste(hatname, "upper", sep="_")
  z <- resdf[, uname]
  nm <- paste(parname, "is_covered", sep="_")
  resdf[, nm] <- mapply(is_covered, x, y, z)
}

tmpf <- function(piece){
  res <- list()
  for(parname in parnames){
    y <- piece[[parname]]
    hatname <- gsub("theta", "hat", parname)
    yh <- piece[[hatname]]    
    nm <- paste(hatname, "bias", sep="_")
    res[[nm]] <- mean(yh - y, na.omit=TRUE)
    nm <- paste(nm, "se", sep="_")
    res[[nm]] <- sd(yh, na.rm=TRUE)/sqrt(length(na.omit(yh)))

    is_covered <- paste(parname, "is_covered", sep="_")
    yic <- piece[[is_covered]]
    nm <- paste(hatname, "coverage", sep="_")
    res[[nm]] <- p <- mean(yic, na.rm=TRUE)
    nm <- paste(nm, "se", sep="_")
    res[[nm]] <- p * (1 - p) / sqrt(length(na.omit(yic)))
    
    nm <- paste(hatname, "bootLower", sep="_")
    res[[nm]] <- quantile(yh, .025, na.rm=TRUE)
    nm <- paste(nm, "se", sep="_")
    res[[nm]] <- tmpff(yh, piece$index, .025)
    
    nm <- paste(hatname, "bootUpper", sep="_")
    res[[nm]] <- quantile(yh, .975, na.rm=TRUE)
    nm <- paste(nm, "se", sep="_")
    res[[nm]] <- tmpff(yh, piece$index, .975)
  }
  as.data.frame(res)
}

k <- 10
resdf$index <- rep(1:k, length.out=nrow(resdf))
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(resdf, c("prob", "floor_times", "maxinftime",
                                  "limit_staff_to_one_day"), tmpf)

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