add a helper function for bootstrapping

DESCRIPTION

@@ -20,6 +20,7 @@ Imports:
     parallel,
     stats
 Suggests:
+    boot (>= 1.3-10),
     coala (>= 0.2.1),
     knitr,
     rmarkdown,

NAMESPACE

@@ -8,6 +8,7 @@ S3method(create_jaatha_model,"function")
 S3method(create_jaatha_model,coalmodel)
 S3method(fit_glm,jaatha_model)
 S3method(fit_glm,jaatha_stat_basic)
+export(boot_jaatha)
 export(create_jaatha_data)
 export(create_jaatha_model)
 export(create_jaatha_stat)

R/boot_jaatha.R

@@ -0,0 +1,60 @@
+#' Parametric Bootstrapping of Jaatha Estimates
+#' 
+#' This function is a helper function for using the \code{\link[boot]{boot}}
+#' function to bootstrap Jaatha estimates. Each bootstap replication requires
+#' a complete jaatha estimation on data simulated with the original parameter
+#' estimates. Therefore, bootstrapping is normally computationally demanding and
+#' should be executed on a computing cluster. The jaatha analyses are
+#' resticted to a single CPU cores, so that as many replicas as possible can 
+#' be executed in parallel using the corresponding options of 
+#' \code{\link[boot]{boot}}.
+#' 
+#' @param model The jaatha model
+#' @param data The jaatha data
+#' @param results The results of an \code{\link{jaatha}} analysis performed with 
+#'   the same \code{model} and \code{data} as passed to this function. 
+#' @param R The number of bootstrapping replicates that are performed.
+#' @param ... Additional arguments that are passed on \code{\link[boot]{boot}}.
+#'   It is highly recommended to use its \code{parallel} and \code{ncpus} 
+#'   options to parallelize the bootstrap replicates.
+#' @return The result of \code{\link[boot]{boot}}. This object can be used to
+#'   estimate standard errors or confidence intervals of the estimates using
+#'   the functions available in package \pkg{boot}.
+#' 
+#' @export
+boot_jaatha <- function(model, data, results, R, ...) {
+  require_package("boot")
+
+  args <- results$args
+  sim_func <- model$get_sim_func()
+
+  log_folder <- tempfile("logs_")
+  message("Logfiles in folder: ", tempdir())
+  message("This might take a while...")
+  dir.create(log_folder)
+
+  jaatha_stat <- function(data) {
+    null <- capture.output({
+      results <- jaatha(model, data,
+                        repetitions = args$repetition,
+                        sim = args$sim,
+                        max_steps = args$max_steps,
+                        init_method = args$init_method,
+                        cores = 1)
+    }, 
+    file = tempfile(paste0("boot_log_", Sys.getpid(), "_")), 
+    type = "message")
+
+    results$estimate
+  }
+
+  simulate_data <- function(data, param) {
+    sim_data <- sim_func(param)
+    create_jaatha_data.default(sim_data, model)
+  }
+
+  boot::boot(data, jaatha_stat, R = R, 
+             sim = "parametric", 
+             ran.gen = simulate_data, 
+             mle = results$estimate, ...)
+}

R/jaatha_log.R

@@ -88,7 +88,7 @@ jaatha_log_class <- R6Class("jaatha_log",
       "creates the results list the main function returns"
       best_estimate <- self$get_best_estimates(1, TRUE)
       param <- as.numeric(best_estimate[1, -(1:3)])
-      res <- list(param = private$par_ranges$denormalize(param),
+      res <- list(estimate = private$par_ranges$denormalize(param),
                   loglikelihood = as.numeric(best_estimate[1, 3]),
                   converged = all(private$converged),
                   args = list(repetitions = private$reps,

R/jaatha_model.R

@@ -90,7 +90,8 @@ jaatha_model_class <- R6Class("jaatha_model",
       invisible(sim_data)
     },
     get_scaling_factor = function() private$scaling_factor,
-    get_par_number = function() private$par_ranges$get_par_number()
+    get_par_number = function() private$par_ranges$get_par_number(),
+    get_sim_func = function() private$sim_func
   )
 )
 

tests/testthat/test-zzz-bootstrapping.R

@@ -0,0 +1,19 @@
+context("Bootstrapping")
+
+test_that("the front-end for boot works", {
+  skip_if_not_installed("boot")
+
+  model <- create_test_model()
+  data <- create_test_data(model)
+
+  results <- list(estimate = c(5.05, 5.05),
+                  loglikelihood = -5,
+                  converged = TRUE,
+                  args = list(repetitions = 1,
+                              sim = 10,
+                              max_steps = 10,
+                              init_method = "middle"))
+
+  capture.output(boot_values <- boot_jaatha(model, data, results, 10))
+  expect_is(boot_values, "boot")
+})