add the wrappers for glm (#32)

* add the wrappers for glm

* [skip style] [skip vbump] Restyle files

* Empty

* update lint

* update wordlist

* [skip style] [skip vbump] Restyle files

* update docs and namespace

---------

Co-authored-by: github-actions <41898282+github-actions[bot]@users.noreply.github.com>

DESCRIPTION

@@ -14,6 +14,7 @@ Depends:
 Imports:
     checkmate,
     numDeriv,
+    prediction,
     sandwich,
     stats
 Suggests:

NAMESPACE

@@ -15,14 +15,21 @@ export(h_jac_ratio)
 export(h_odds_ratio)
 export(h_ratio)
 export(predict_counterfactual)
+export(robin_glm)
 export(treatment_effect)
 export(vcovANHECOVA)
 import(checkmate)
 importFrom(numDeriv,jacobian)
+importFrom(prediction,find_data)
 importFrom(sandwich,vcovHC)
 importFrom(stats,coefficients)
+importFrom(stats,family)
 importFrom(stats,fitted)
+importFrom(stats,gaussian)
+importFrom(stats,glm)
 importFrom(stats,model.matrix)
 importFrom(stats,model.response)
 importFrom(stats,predict)
 importFrom(stats,residuals)
+importFrom(stats,terms)
+importFrom(stats,var)

R/RobinCar2-package.R

@@ -7,6 +7,8 @@
 
 #' @import checkmate
 #' @importFrom numDeriv jacobian
-#' @importFrom stats predict residuals fitted model.response model.matrix coefficients
+#' @importFrom stats predict residuals fitted model.response model.matrix coefficients family
+#' gaussian terms glm var
 #' @importFrom sandwich vcovHC
+#' @importFrom prediction find_data
 NULL

R/predict_couterfactual.R

@@ -15,7 +15,7 @@ predict_counterfactual <- function(fit, treatment, data, unbiased) {
 }
 
 #' @export
-predict_counterfactual.lm <- function(fit, treatment, data, unbiased = TRUE) {
+predict_counterfactual.lm <- function(fit, treatment, data = find_data(fit), unbiased = TRUE) {
   treatment <- h_get_vars(treatment)
   assert_data_frame(data)
   assert_subset(unlist(treatment), colnames(data))
@@ -72,6 +72,6 @@ predict_counterfactual.lm <- function(fit, treatment, data, unbiased = TRUE) {
 }
 
 #' @export
-predict_counterfactual.glm <- function(fit, treatment, data = fit$data, unbiased = TRUE) {
+predict_counterfactual.glm <- function(fit, treatment, data = find_data(fit), unbiased = TRUE) {
   predict_counterfactual.lm(fit = fit, data = data, treatment = treatment, unbiased = unbiased)
 }

R/robin_glm.R

@@ -0,0 +1,55 @@
+#' Covariate adjusted glm model
+#' @param formula (`formula`) A formula of analysis.
+#' @param data (`data.frame`) Input data frame.
+#' @param treatment (`formula` or `character(1)`) A formula of treatment assignment or assignment by stratification,
+#' or a string name of treatment assignment.
+#' @param contrast (`function` or `character(1)`) A function to calculate the treatment effect, or character of
+#' "difference", "risk_ratio", "odds_ratio" for default contrasts.
+#' @param contrast_jac (`function`) A function to calculate the Jacobian of the contrast function. Ignored if using
+#' default contrasts.
+#' @param vcov (`function`) A function to calculate the variance-covariance matrix of the treatment effect,
+#' including `vcovHC` and `vcovANHECOVA`.
+#' @param family (`family`) A family object of the glm model.
+#' @param ... Additional arguments passed to `vcov`. For finer control of glm, refer to usage of `treatment_effect`,
+#' `difference`, `risk_ratio`, `odds_ratio`.
+#' @export
+#' @examples
+#' robin_glm(y ~ treatment * s1, data = dummy_data, treatment = treatment ~ s1, contrast = "difference")
+robin_glm <- function(
+    formula, data, treatment, contrast = "difference",
+    contrast_jac = NULL, vcov = vcovANHECOVA, family = gaussian, ...) {
+  attr(formula, ".Environment") <- environment()
+  fit <- glm(formula, family = family, data = data)
+  pc <- predict_counterfactual(fit, treatment, data, unbiased = TRUE)
+  has_interaction <- h_interaction(formula, treatment)
+  if (has_interaction && identical(vcov, vcovHC) && !identical(contrast, "difference")) {
+    stop(
+      "Huber-White standard error only works for difference contrasts in models without interaction term."
+    )
+  }
+  if (identical(contrast, "difference")) {
+    difference(pc)
+  } else if (identical(contrast, "risk_ratio")) {
+    risk_ratio(pc)
+  } else if (identical(contrast, "odds_ratio")) {
+    odds_ratio(pc)
+  } else {
+    assert_function(contrast)
+    assert_function(contrast_jac, null.ok = TRUE)
+    if (is.null(contrast_jac)) {
+      contrast_jac <- function(x) {
+        numDeriv::jacobian(contrast, x)
+      }
+    }
+    treatment_effect(pc, eff_measure = contrast, eff_jacobian = contrast_jac, variance = vcov, ...)
+  }
+}
+
+h_interaction <- function(formula, treatment) {
+  assert_formula(formula)
+  treatment <- h_get_vars(treatment)
+  assert_subset(treatment$treatment, all.vars(formula[[length(formula)]]))
+  tms <- terms(formula)
+  fct <- attr(tms, "factors")
+  any(fct[treatment$treatment, ] %in% c(1, 2) & colSums(fct != 0) > 1)
+}

R/treatment_effect.R

@@ -21,7 +21,6 @@ treatment_effect.prediction_cf <- function(
     test_null(variance)
   )
   assert_function(eff_measure)
-  assert_function(eff_jacobian)
   assert_vector(pair)
   assert(
     test_subset(pair, names(object)),
@@ -36,6 +35,7 @@ treatment_effect.prediction_cf <- function(
     if (missing(eff_jacobian)) {
       trt_jac <- numDeriv::jacobian(eff_measure, object[pair])
     } else {
+      assert_function(eff_jacobian)
       trt_jac <- eff_jacobian(object[pair])
     }
     trt_var <- trt_jac %*% inner_variance %*% t(trt_jac)
@@ -57,8 +57,8 @@ treatment_effect.prediction_cf <- function(
 #' @inheritParams predict_counterfactual
 treatment_effect.lm <- function(
     object, pair = names(object), variance = vcovANHECOVA, eff_measure, eff_jacobian,
-    treatment, data, unbiased = TRUE, ...) {
-  pc <- predict_counterfactual(object, treatment, data, unbiased)
+    treatment, data = find_data(object), unbiased = TRUE, ...) {
+  pc <- predict_counterfactual(object, data = data, treatment, unbiased)
   if (missing(pair)) {
     treatment_effect(pc, pair = , , variance = variance, eff_measure = eff_measure, eff_jacobian = eff_jacobian, ...)
   } else {
@@ -69,7 +69,7 @@ treatment_effect.lm <- function(
 #' @export
 treatment_effect.glm <- function(
     object, pair, variance = vcovANHECOVA, eff_measure, eff_jacobian,
-    treatment, data = object$data, unbiased = TRUE, ...) {
+    treatment, data = find_data(object), unbiased = TRUE, ...) {
   pc <- predict_counterfactual(object, treatment, data, unbiased)
   if (missing(pair)) {
     treatment_effect(pc, pair = , , variance = variance, eff_measure = eff_measure, eff_jacobian = eff_jacobian, ...)

inst/WORDLIST

@@ -1,5 +1,7 @@
 ANHECOVA
 RobinCar
+glm
 jacobian
+jacobians
 pbo
 trt

tests/testthat/test-robin_glm.R

@@ -0,0 +1,32 @@
+# h_interaction ----
+
+test_that("h_interaction works correctly", {
+  expect_false(h_interaction(y ~ trt + z, treatment = trt ~ x))
+  expect_true(h_interaction(y ~ trt:z, treatment = trt ~ x))
+  expect_true(h_interaction(trt * y ~ trt:z, treatment = trt ~ x))
+  expect_true(h_interaction(y ~ trt:z, treatment = "trt"))
+})
+
+# robin_glm ----
+
+test_that("robin_glm works correctly", {
+  expect_silent(
+    robin_glm(
+      y ~ treatment * s1,
+      data = dummy_data, treatment = treatment ~ s1,
+      contrast = "difference", vcov = vcovHC
+    )
+  )
+  expect_silent(robin_glm(y_b ~ treatment * s1, data = dummy_data, treatment = treatment ~ s1, contrast = "difference"))
+  expect_silent(robin_glm(y_b ~ treatment * s1, data = dummy_data, treatment = treatment ~ s1, contrast = "risk_ratio"))
+  expect_silent(robin_glm(y_b ~ treatment * s1, data = dummy_data, treatment = treatment ~ s1, contrast = "odds_ratio"))
+  expect_error(
+    robin_glm(
+      y_b ~ treatment * s1,
+      data = dummy_data, treatment = treatment ~ s1,
+      contrast = "odds_ratio", vcov = vcovHC
+    ),
+    "Huber-White standard error only works for difference contrasts in models without interaction term."
+  )
+  expect_silent(robin_glm(y_b ~ treatment * s1, data = dummy_data, treatment = treatment ~ s1, contrast = h_diff))
+})

tests/testthat/test-treatment_effect.R

@@ -135,3 +135,17 @@ test_that("treatment_effect works as expected", {
     tolerance = 1e-6
   )
 })
+
+test_that("treatment_effect works as expected for custom contrast", {
+  pc <- predict_counterfactual(fit_binom, treatment ~ s1)
+  h_contrast <- function(x) {
+    v <- outer(x, x, `-`)
+    v[lower.tri(v)]
+  }
+  expect_silent(treatment_effect(pc, eff_measure = h_contrast))
+})
+
+test_that("treatment_effect works for lm/glm object", {
+  expect_silent(treatment_effect(fit_binom, treatment = treatment ~ s1, eff_measure = h_diff))
+  expect_silent(treatment_effect(fit_lm, treatment = treatment ~ s1, eff_measure = h_diff))
+})