Add extract apollo

DESCRIPTION

@@ -19,7 +19,7 @@ Suggests:
   testthat (>= 2.0.0),
   lmtest (>= 0.9-34)
 Depends: R (>= 3.5)
-Imports: methods, stats, httr
+Imports: methods, stats, httr, apollo, janitor
 Enhances:
   AER,
   alpaca,

R/extract.R

@@ -8859,3 +8859,91 @@ extract.logitr <- function(model,
 #' @export
 setMethod("extract", signature = className("logitr", "logitr"),
           definition = extract.logitr)
+
+
+
+# -- extract.apollo_estimation (Apollo) ----------------------------------------
+
+#' @noRd
+extract.apollo_estimation <- function(model,
+                                      wtpest = NULL,
+                                      se = "rob",
+                                      ...) {
+  # validate 'se' argument
+  if (!se %in% c("rob", "normal", "bs")) {
+    stop("Invalid value for 'se'. Please use one of 'rob', 'normal', or 'bs'.")
+  }
+  # bootstrap SEs must exist in the model object
+  if (se == "bs" && !"bootse" %in% names(model)) {
+    stop("No bootstrapped SE found. The 'model' must contain 'bootse' for se = 'bs'.")
+  }
+
+  # pull out the standard apollo output table
+  settings <- list(printPVal = TRUE)
+  if (is.null(wtpest)) {
+    estimated <- janitor::clean_names(
+      as.data.frame(apollo::apollo_modelOutput(model, settings))
+    )
+    # pick which SE & p-value columns to use
+    switch(se,
+           rob = {
+             estimated$se <- estimated$rob_s_e
+             estimated$pv <- estimated$p_1_sided_2
+           },
+           bs = {
+             estimated$se <- estimated$bootstrap_s_e
+             estimated$pv <- estimated$p_1_sided_3
+           },
+           normal = {
+             estimated$se <- estimated$s_e
+             estimated$pv <- estimated$p_1_sided
+           }
+    )
+  } else {
+    # user-supplied WTP table
+    estimated <- wtpest
+    colnames(estimated) <- c("estimate", "se", "robt", "pv")
+  }
+
+  # clean up the coefficient names
+  coefnames <-  rownames(estimated)
+
+  # assemble into a texreg object
+  tr <- createTexreg(
+    coef.names   = coefnames,
+    coef         = estimated[["estimate"]],
+    se           = estimated[["se"]],
+    pvalues      = estimated[["pv"]],
+    gof.names    = c("Num. obs.", "Num. indiv.", "Log Likelihood (Null)", "Log Likelihood (Fit)"),
+    gof          = c(
+      model[["nObsTot"]],
+      model[["nIndivs"]],
+      model[["LL0"]][[1]],
+      model[["LLout"]][[1]]
+    ),
+    gof.decimal  = c(FALSE, FALSE, TRUE, TRUE)
+  )
+
+  return(tr)
+}
+
+#' \code{\link{extract}} method for \code{apollo} objects
+#'
+#' \code{\link{extract}} method for \code{apollo} objects created by
+#' \code{apollo::apollo_estimate()}.
+#'
+#' @param model   An object of class \code{apollo}.
+#' @param wtpest   Optional data.frame of willingness-to-pay estimates (and s.e.'s).
+#' @param se       Which standard errors to use: \code{"rob"}, \code{"normal"}, or \code{"bs"}.
+#' @param ...      Currently ignored.
+#'
+#' @method extract apollo
+#' @aliases extract.apollo
+#' @author Julian Sagebiel <julian.sagebiel@idiv.de>
+#' @export
+setMethod(
+  f          = "extract",
+  signature  = className("apollo", "apollo"),
+  definition = extract.apollo_estimation
+)
+

tests/testthat/test-extract.R

@@ -1231,3 +1231,131 @@ test_that("extract logitr objects from the logitr package", {
   expect_length(tr@gof.decimal, 4)
   expect_equivalent(dim(matrixreg(mnl_pref)), c(15, 2))
 })
+
+
+
+
+
+
+
+test_that("extract.apollo returns a proper texreg object for a tiny MNL model", {
+  testthat::skip_on_cran()
+  testthat::skip_if_not_installed("apollo", minimum_version = "0.3.0")
+
+  require(apollo)
+
+  # load and subset example data for speed
+  data("apollo_swissRouteChoiceData", package = "apollo")
+  database <- apollo_swissRouteChoiceData[1:200, ]
+
+  # initialise Apollo and control settings
+  apollo::apollo_initialise()
+  apollo_control <- list(
+    modelName  = "test_mnl",
+    modelDescr = "Tiny MNL for testing",
+    indivID    = "ID",
+    nCores     = 1
+  )
+  apollo_beta  <- c(ti = 0, cost = 0)
+  apollo_fixed <- c()
+
+  # validate inputs (uses database & apollo_control from global env)
+  apollo_inputs <- apollo_validateInputs(database = database,
+                                         apollo_control = apollo_control,
+                                         apollo_beta = apollo_beta,
+                                         apollo_fixed = apollo_fixed)
+
+  # define extremely simple MNL probabilities function
+  apollo_probabilities <- function(apollo_beta, apollo_inputs, functionality = "estimate") {
+
+
+
+
+
+       apollo_attach(apollo_beta, apollo_inputs)
+    on.exit(apollo_detach(apollo_beta, apollo_inputs))
+
+
+
+    ### Create list of probabilities P
+    P = list()
+
+    V <- list(
+      alt1 =  ti * tt1 + cost * tc1,
+      alt2 =  ti * tt2 + cost * tc2
+    )
+    mnl_settings <- list(
+      alternatives = c(alt1 = 1, alt2 = 2),
+      avail        = 1,
+      choiceVar    = choice,
+      V            = V
+    )
+    ### Compute probabilities using MNL model
+    P[["model"]] = apollo_mnl(mnl_settings, functionality)
+
+    ### Take product across observation for same individual
+    P = apollo_panelProd(P, apollo_inputs, functionality)
+
+    ### Prepare and return outputs of function
+    P = apollo_prepareProb(P, apollo_inputs, functionality)
+    return(P)
+
+  }
+
+  # estimate with very few iterations
+  model <- apollo_estimate(
+    apollo_beta,
+    apollo_fixed,
+    apollo_probabilities,
+    apollo_inputs,
+    estimate_settings = list(
+      estimationRoutine = "bfgs",
+      maxIterations    = 200,
+      silent           = TRUE
+    )
+  )
+
+unlink("test_mnl_iterations.csv")
+
+  # extract texreg object
+  tr <- extract(model)
+
+  # basic structure
+  expect_s4_class(tr, "texreg")
+  expect_equal(length(tr@coef), length(model$estimate))
+
+  # GOF names & values should match what our extract.apollo implementation uses
+  expect_equal(
+    tr@gof.names,
+    c("Num. obs.", "Num. indiv.", "Log Likelihood (Null)", "Log Likelihood (Fit)")
+  )
+  expect_equal(
+    tr@gof,
+    c(model$nObsTot, model$nIndivs, unname(model$LL0), unname(model$LLout))
+  )
+
+
+  wtpest <- data.frame(
+    estimate = c(0.5, 0.3),
+    se       = c(0.05, 0.03),
+    robt     = c(NA, NA),
+    pv       = c(0.10, 0.20)
+  )
+
+  # extract with wtpest
+
+  tr2 <- extract(model, wtpest = wtpest)
+
+  # coefficients, SEs and p-values come from our wtpest
+  expect_equal(tr2@coef, wtpest$estimate)
+  expect_equal(tr2@se,   wtpest$se)
+  expect_equal(tr2@pvalues, wtpest$pv)
+
+
+  expect_error(
+    extract(model = model, se = "invalid"),
+    "Invalid value for 'se'. Please use one of 'rob', 'normal', or 'bs'.")
+
+
+})
+

tests/testthat/test-plotreg.R

@@ -18,10 +18,10 @@ test_that("plotreg works", {
   p3 <- plotreg(list(m1, m2), custom.title = "My plot")
   p4 <- plotreg(list(m1, m2), custom.note = "My note")
 
-  expect_true(ggplot2::is.ggplot(p1))
-  expect_true(ggplot2::is.ggplot(p2))
-  expect_true(ggplot2::is.ggplot(p3))
-  expect_true(ggplot2::is.ggplot(p4))
+  expect_true(ggplot2::is_ggplot(p1))
+  expect_true(ggplot2::is_ggplot(p2))
+  expect_true(ggplot2::is_ggplot(p3))
+  expect_true(ggplot2::is_ggplot(p4))
 
   expect_s3_class(p1$data, "data.frame")
   expect_s3_class(p1, "gg")
@@ -45,7 +45,7 @@ test_that("plotreg works with confidence intervals using the biglm package", {
   gg <- log(Volume) ~ log(Girth) + log(Height) + offset(2 * log(Girth) + log(Height))
   b <- biglm::bigglm(gg, data = trees, chunksize = 10, sandwich = TRUE)
   p6 <- plotreg(list(a, b))
-  expect_true(ggplot2::is.ggplot(p6))
+  expect_true(ggplot2::is_ggplot(p6))
   expect_equal(dim(p6$data), dim(readRDS("../files/PlotDataFrameCI.RDS")))
   # saveRDS(p6$data, "../files/PlotDataFrameCI.RDS")
   expect_identical(p6$labels$y, "Bars denote SEs (95%). Circle points denote significance.")
@@ -63,7 +63,7 @@ test_that("plotreg -odds ratio", {
   m1 <- lm(weight ~ group)
   p7 <- plotreg(m1, override.coef = exp(m1$coefficients), ci.force = TRUE, ci.test = 1)
 
-  expect_true(ggplot2::is.ggplot(p7))
+  expect_true(ggplot2::is_ggplot(p7))
   expect_s3_class(p7$data, "data.frame")
   expect_s3_class(p7, "gg")
   expect_equal(dim(p7$data), dim(readRDS("../files/PlotDataFrameOR.RDS"))) # test if data frame is correctly constructed