## code to reproduce Table 4
## THIS FILE: combining n=1000 and n=5000 simulations into table
options(width=150)
library(xtable)
sessionInfo()

## set parameters
lambda <- seq(from=0, to=2, length=20)
nlambda <- length(lambda)
nparam <- 1
.beta <- -1.5

#########################################
## set input and output files
#########################################

## input: SIMEX results
infile1a <- "~/simex-vfunc/code-files-for-psychometrika-paper-aug2016/code-for-table4-n1000.Robj"
infile5a <- "~/simex-vfunc/code-files-for-psychometrika-paper-aug2016/code-for-table4-n5000.Robj"
## output: combined table
tabfilec <- "~/simex-vfunc/code-files-for-psychometrika-paper-aug2016/code-for-table4-combine.tex"

#########################################
## function to make table from output
#########################################
doit <- function(infilea){

  load(infilea)
  res <- do.call("rbind", lapply(bigres, function(x){ x$res }))
  
  tapply(res$b1,  res$method, summary)
  
  ## get mean, SD, RMSE by estimator by averaging across sims
  res$bmean   <- ave(res$b1, res$method)
  res$bbias   <- res$bmean - .beta
  res$bsd     <- ave(res$b1, res$method, FUN=sd)
  res$rmse    <- ave( (res$b1 - .beta)^2, res$method, FUN=function(x){ sqrt(mean(x))})
  res         <- unique(res[,c("method","bmean","bbias","bsd","rmse")])
  
  ## put the pieces together for table, wrapping quartic to wide
  rownames(res) <- 1:nrow(res)
  res2 <- res[c(1:3,4:10),]
  res4 <- res[c(1:3,11:17),]
  tab <- data.frame(res2, bmean4 = res4$bmean, bbias4 = res4$bbias, bsd4 = res4$bsd, rmse4 = res4$rmse)
  
  nv <- setdiff(names(tab),"method")
  for(v in nv){
    tab[,v] <- round(tab[,v],digits=2)
  }
  
  ## order according to paper and output
  print(tab <- tab[c(1,2,3,4,9,5,10,7,6,8),])
  return(tab)
}

###############################################
## create everything for the output files
###############################################
tab1k <- doit(infile1a)
tab5k <- doit(infile5a)

## put the RMSE for 5000 onto the 1000 table
stopifnot(all(tab1k$method == tab5k$method))
tab1k$rmse_5000  <- tab5k$rmse
tab1k$rmse4_5000 <- tab5k$rmse4
o <- c("method","bmean","bbias","bsd","rmse","rmse_5000","bmean4","bbias4","bsd4","rmse4","rmse4_5000")
stopifnot(all(sort(o) == sort(names(tab1k))))
tab1k <- tab1k[,o]

sink(tabfilec)
print(xtable(tab1k))
sink(NULL)