From f9d5c12fd15db0f8dfffc2e37ebb9b9eb4d8c65d Mon Sep 17 00:00:00 2001 From: Liming <36079400+clarkliming@users.noreply.github.com> Date: Mon, 19 Aug 2024 15:26:45 +0800 Subject: [PATCH] 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 | 1 + NAMESPACE | 7 ++++ R/RobinCar2-package.R | 4 +- R/predict_couterfactual.R | 4 +- R/robin_glm.R | 55 ++++++++++++++++++++++++++ R/treatment_effect.R | 8 ++-- inst/WORDLIST | 2 + man/robin_glm.Rd | 45 +++++++++++++++++++++ tests/testthat/test-robin_glm.R | 32 +++++++++++++++ tests/testthat/test-treatment_effect.R | 14 +++++++ 10 files changed, 165 insertions(+), 7 deletions(-) create mode 100644 R/robin_glm.R create mode 100644 man/robin_glm.Rd create mode 100644 tests/testthat/test-robin_glm.R diff --git a/DESCRIPTION b/DESCRIPTION index e450e50..c2852fd 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -14,6 +14,7 @@ Depends: Imports: checkmate, numDeriv, + prediction, sandwich, stats Suggests: diff --git a/NAMESPACE b/NAMESPACE index c48614f..6377bd6 100644 --- a/NAMESPACE +++ b/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) diff --git a/R/RobinCar2-package.R b/R/RobinCar2-package.R index c4be6d6..dc1cbd9 100644 --- a/R/RobinCar2-package.R +++ b/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 diff --git a/R/predict_couterfactual.R b/R/predict_couterfactual.R index 445013d..441e9dd 100644 --- a/R/predict_couterfactual.R +++ b/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) } diff --git a/R/robin_glm.R b/R/robin_glm.R new file mode 100644 index 0000000..09302b6 --- /dev/null +++ b/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) +} diff --git a/R/treatment_effect.R b/R/treatment_effect.R index 06060f3..3cca457 100644 --- a/R/treatment_effect.R +++ b/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, ...) diff --git a/inst/WORDLIST b/inst/WORDLIST index 1085c9c..70f41ea 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -1,5 +1,7 @@ ANHECOVA RobinCar +glm jacobian +jacobians pbo trt diff --git a/man/robin_glm.Rd b/man/robin_glm.Rd new file mode 100644 index 0000000..342d299 --- /dev/null +++ b/man/robin_glm.Rd @@ -0,0 +1,45 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/robin_glm.R +\name{robin_glm} +\alias{robin_glm} +\title{Covariate adjusted glm model} +\usage{ +robin_glm( + formula, + data, + treatment, + contrast = "difference", + contrast_jac = NULL, + vcov = vcovANHECOVA, + family = gaussian, + ... +) +} +\arguments{ +\item{formula}{(\code{formula}) A formula of analysis.} + +\item{data}{(\code{data.frame}) Input data frame.} + +\item{treatment}{(\code{formula} or \code{character(1)}) A formula of treatment assignment or assignment by stratification, +or a string name of treatment assignment.} + +\item{contrast}{(\code{function} or \code{character(1)}) A function to calculate the treatment effect, or character of +"difference", "risk_ratio", "odds_ratio" for default contrasts.} + +\item{contrast_jac}{(\code{function}) A function to calculate the Jacobian of the contrast function. Ignored if using +default contrasts.} + +\item{vcov}{(\code{function}) A function to calculate the variance-covariance matrix of the treatment effect, +including \code{vcovHC} and \code{vcovANHECOVA}.} + +\item{family}{(\code{family}) A family object of the glm model.} + +\item{...}{Additional arguments passed to \code{vcov}. For finer control of glm, refer to usage of \code{treatment_effect}, +\code{difference}, \code{risk_ratio}, \code{odds_ratio}.} +} +\description{ +Covariate adjusted glm model +} +\examples{ +robin_glm(y ~ treatment * s1, data = dummy_data, treatment = treatment ~ s1, contrast = "difference") +} diff --git a/tests/testthat/test-robin_glm.R b/tests/testthat/test-robin_glm.R new file mode 100644 index 0000000..85b42d5 --- /dev/null +++ b/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)) +}) diff --git a/tests/testthat/test-treatment_effect.R b/tests/testthat/test-treatment_effect.R index 0c4bbab..96a7137 100644 --- a/tests/testthat/test-treatment_effect.R +++ b/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)) +})