R version 3.2.4 (2016-03-10) -- "Very Secure Dishes"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> ## code to reproduce Table 4
> ## THIS FILE: combining n=1000 and n=5000 simulations into table
> options(width=150)
> library(xtable)
> sessionInfo()
R version 3.2.4 (2016-03-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux Server release 6.8 (Santiago)

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8   
 [6] LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] xtable_1.8-2
> 
> ## 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)
         method bmean bbias  bsd rmse bmean4 bbias4 bsd4 rmse4
1         ideal -1.50  0.00 0.11 0.10  -1.50   0.00 0.11  0.10
2         naive -1.31  0.19 0.10 0.22  -1.31   0.19 0.10  0.22
3      regcalib -1.43  0.07 0.10 0.12  -1.43   0.07 0.10  0.12
4     simex2_gx -1.51 -0.01 0.12 0.12  -1.51  -0.01 0.13  0.13
9     simex2_aw -1.47  0.03 0.11 0.12  -1.45   0.05 0.12  0.13
5  simex2_avggx -1.59 -0.09 0.13 0.16  -1.65  -0.15 0.16  0.22
10  simex2_rcal -1.51 -0.01 0.12 0.12  -1.50   0.00 0.13  0.13
7     simex2_gw -1.48  0.02 0.11 0.11  -1.45   0.05 0.13  0.14
6  simex2_avggw -1.63 -0.13 0.13 0.18  -1.70  -0.20 0.16  0.26
8   simex2_gexw -1.49  0.01 0.12 0.12  -1.49   0.01 0.13  0.13
> tab5k <- doit(infile5a)
         method bmean bbias  bsd rmse bmean4 bbias4 bsd4 rmse4
1         ideal -1.50  0.00 0.05 0.05  -1.50   0.00 0.05  0.05
2         naive -1.31  0.19 0.04 0.20  -1.31   0.19 0.04  0.20
3      regcalib -1.44  0.06 0.05 0.08  -1.44   0.06 0.05  0.08
4     simex2_gx -1.51 -0.01 0.05 0.05  -1.51  -0.01 0.06  0.06
9     simex2_aw -1.47  0.03 0.05 0.06  -1.46   0.04 0.06  0.07
5  simex2_avggx -1.59 -0.09 0.06 0.11  -1.65  -0.15 0.07  0.16
10  simex2_rcal -1.51 -0.01 0.05 0.05  -1.51  -0.01 0.06  0.06
7     simex2_gw -1.48  0.02 0.05 0.05  -1.45   0.05 0.06  0.08
6  simex2_avggw -1.63 -0.13 0.06 0.14  -1.70  -0.20 0.07  0.21
8   simex2_gexw -1.49  0.01 0.05 0.05  -1.49   0.01 0.06  0.06
> 
> ## 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)
> 
> proc.time()
   user  system elapsed 
  0.740   0.054   1.092