fit outcome model method (#156)

* fitter and test

* another test

* add pkgdown entries

Author: gravesti

.pre-commit-config.yaml

@@ -18,6 +18,7 @@ repos:
           - parglm
           - sandwich
           - duckdb
+          - lmtest
     # codemeta must be above use-tidy-description when both are used
     # -   id: codemeta-description-updated
     -   id: use-tidy-description

DESCRIPTION

@@ -29,6 +29,7 @@ Imports:
     DBI,
     duckdb,
     formula.tools,
+    lmtest,
     methods,
     mvtnorm,
     parglm,

R/generics.R

@@ -174,12 +174,14 @@ setGeneric("read_expanded_data", function(object, period = NULL) standardGeneric
 #' @param object A `te_outcome_fitter` object
 #' @param data `data.frame` containing outcomes and covariates as defined in `formula`.
 #' @param formula `formula` describing the model.
+#' @param weights `numeric` vector of weights.
 #'
 #' @return An object of class `te_outcome_fitted`
 #' @export
 #' @keywords internal
 #' @examples
-setGeneric("fit_outcome_model", function(object, data, formula) standardGeneric("fit_outcome_model"))
+#' fit_outcome_model
+setGeneric("fit_outcome_model", function(object, data, formula, weights = NULL) standardGeneric("fit_outcome_model"))
 
 #' Method for fitting weight models
 #'

R/te_model_fitter.R

@@ -1,3 +1,6 @@
+#' @include generics.R
+NULL
+
 #' Outcome Model Fitter Class
 #'
 #' This is a virtual class which other outcome model fitter classes should inherit from. Objects of these class exist to
@@ -65,3 +68,51 @@ setMethod(
     )
   }
 )
+
+#' @rdname fit_outcome_model
+setMethod(
+  f = "fit_outcome_model",
+  signature = "te_stats_glm_logit",
+  function(object, data, formula, weights) {
+    data$weights <- if (is.null(weights)) rep(1, nrow(data)) else weights
+    model <- glm(
+      formula = formula,
+      data = data,
+      family = binomial("logit"),
+      x = FALSE,
+      y = FALSE,
+      weights = weights
+    )
+    if (!is.na(object@save_path)) {
+      if (!dir.exists(object@save_path)) dir.create(object@save_path, recursive = TRUE)
+      file <- tempfile(pattern = "model_", tmpdir = object@save_path, fileext = ".rds")
+      saveRDS(model, file = file)
+      save_path <- data.frame(save = file)
+    }
+
+    vcov <- sandwich::vcovCL(
+      model,
+      cluster = data[["id"]],
+      type = NULL,
+      sandwich = TRUE,
+      fix = FALSE
+    )
+
+    model_list <- list(
+      model = model,
+      vcov = vcov
+    )
+
+    coef_obj <- lmtest::coeftest(model, vcov. = vcov, save = TRUE)
+    summary_list <- list()
+    summary_list[["tidy"]] <- broom::tidy(coef_obj, conf.int = TRUE)
+    summary_list[["glance"]] <- broom::glance(coef_obj)
+    if (!is.na(object@save_path)) summary_list[["save_path"]] <- save_path
+
+    new(
+      "te_stats_glm_logit_outcome_fitted",
+      model = model_list,
+      summary = summary_list
+    )
+  }
+)

R/te_outcome_model.R

@@ -1,8 +1,35 @@
+#' Fitted Outcome Model Object
+#'
+#' @slot model list containing fitted model objects.
+#' @slot summary list of data.frames. Tidy model summaries a la `broom()` and `glance()`
+setClass(
+  "te_outcome_fitted",
+  slots = c(
+    model = "list",
+    summary = "list"
+  )
+)
+
+#' Fitted Outcome Model Object
+#'
+#' @slot formula `formula` object for the model fitting
+#' @slot treatment_var character. The treatment variable
+#' @slot adjustment_vars character. Adjustment variables
+#' @slot model_fitter Model fitter object
+#' @slot fitted list. Saves the model objects
 setClass("te_outcome_model",
   slots = c(
     formula = "formula",
     treatment_var = "character",
     adjustment_vars = "character",
-    model_fitter = "te_model_fitter"
+    model_fitter = "te_model_fitter",
+    fitted = "te_outcome_fitted"
   )
 )
+
+# te_stats_glm_logit_outcome_fitted -----
+
+setClass(
+  "te_stats_glm_logit_outcome_fitted",
+  contains = "te_outcome_fitted"
+)

R/te_weights.R

@@ -2,7 +2,7 @@
 
 #' Fitted Weights Object
 #'
-#' @slot specification list. The parameters specifying how the model should be fit
+#' @slot label string. A short description of the model
 #' @slot summary list of data.frames. Tidy model summaries a la `broom()` and `glance()`
 #' @slot fitted list. Saves the model objects or at least summaries if large
 setClass("te_weights_fitted",

_pkgdown.yml

@@ -44,6 +44,8 @@ reference:
   - te_data-class
   - te_datastore-class
   - te_model_fitter-class
+  - te_outcome_fitted-class
+  - te_outcome_model-class
 - title: Internal Methods
 - contents:
   - save_expanded_data

man/fit_outcome_model.Rd

@@ -27,3 +27,62 @@ test_that("fit_weights_model works for stats_glm_logit", {
   saved_model <- readRDS(result@summary$save_path$path)
   expect_class(saved_model, "glm")
 })
+
+
+test_that("fit_outcome_model works for stats_glm_logit", {
+  object <- stats_glm_logit(NA)
+  result <- fit_outcome_model(
+    object,
+    data = vignette_switch_data,
+    formula = outcome ~ assigned_treatment + followup_time + nvarC,
+    weights = vignette_switch_data$weight
+  )
+
+  expect_class(result, "te_outcome_fitted")
+  expect_class(result, "te_stats_glm_logit_outcome_fitted")
+
+  expect_equal(
+    result@summary[["tidy"]]$estimate,
+    c(-3.23428770614918, 0.0429082229813285, -0.000849656189863374, -0.0370118224213905)
+  )
+  expect_equal(
+    as.data.frame(result@summary[["tidy"]][2, c("conf.low", "conf.high")]),
+    data.frame(conf.low = -0.501620572971525, conf.high = 0.587437018934182)
+  )
+  expect_equal(result@summary[["glance"]]$df.null, 1939052)
+
+  expect_matrix(result@model$vcov, nrows = 4, ncols = 4)
+  expect_equal(
+    diag(result@model$vcov),
+    c(
+      `(Intercept)` = 0.145504452458149, assigned_treatment = 0.0771872414783003,
+      followup_time = 1.59771524178973e-06, nvarC = 4.40825112262011e-05
+    )
+  )
+})
+
+test_that("fit_outcome_model works for stats_glm_logit with save_dir", {
+  dir <- withr::local_tempdir(pattern = "glm_test")
+  object <- stats_glm_logit(save_path = dir)
+  data <- data.frame(
+    y = rep(c(1, 0, 1, 0), times = c(15, 5, 5, 15)),
+    x = rep(c(1, 0), times = c(20, 20))
+  )
+  result <- fit_outcome_model(
+    object,
+    data = data,
+    formula = y ~ x
+  )
+
+  expect_class(result, "te_outcome_fitted")
+  expect_class(result, "te_stats_glm_logit_outcome_fitted")
+
+  expect_equal(result@summary[["tidy"]]$estimate, c(-1.0986123, 2.1972246))
+  fitted_model <- readRDS(result@summary$save_path$save)
+  expect_equal(
+    fitted_model,
+    result@model$model,
+    ignore_function_env = TRUE,
+    ignore_formula_env = TRUE
+  )
+})