diff --git a/DESCRIPTION b/DESCRIPTION index 7eab876..a2d15d9 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -36,7 +36,8 @@ Imports: future (>= 1.17.0), progressr, data.table (>= 1.13.0), - checkmate (>= 2.1.0) + checkmate (>= 2.1.0), + SuperRiesz URL: https://github.com/nt-williams/lmtp BugReports: https://github.com/nt-williams/lmtp/issues Suggests: diff --git a/R/density_ratios.R b/R/density_ratios.R index 4e7e0e1..325861b 100644 --- a/R/density_ratios.R +++ b/R/density_ratios.R @@ -12,7 +12,7 @@ cf_r <- function(Task, learners, mtp, control, pb) { seed = TRUE) } - trim_ratios(recombine_ratios(future::value(out), Task$folds), control$.trim) + trim_ratios(recombine_ratios(future::value(out), task$folds), control$.trim) } estimate_r <- function(natural, shifted, trt, cens, risk, tau, node_list, learners, pb, mtp, control) { @@ -56,7 +56,7 @@ estimate_r <- function(natural, shifted, trt, cens, risk, tau, node_list, learne ratios <- density_ratios(pred, irv, drv, frv, mtp) densratios[, t] <- ratios - pb() + progress_bar() } list(ratios = densratios, fits = fits) diff --git a/R/estimators.R b/R/estimators.R index 9267228..ee4049c 100644 --- a/R/estimators.R +++ b/R/estimators.R @@ -54,6 +54,8 @@ #' @param learners_trt \[\code{character}\]\cr A vector of \code{mlr3superlearner} algorithms for estimation #' of the outcome regression. Default is \code{c("mean", "glm")}. #' \bold{Only include candidate learners capable of binary classification}. +#' @param trt_method \[\code{character}\]\cr +#' Method for estimating treatment assignment mechanism (default or riesz) #' @param folds \[\code{integer(1)}\]\cr #' The number of folds to be used for cross-fitting. #' @param weights \[\code{numeric(nrow(data))}\]\cr @@ -96,6 +98,7 @@ lmtp_tmle <- function(data, trt, outcome, baseline = NULL, time_vary = NULL, id = NULL, bounds = NULL, learners_outcome = c("mean", "glm"), learners_trt = c("mean", "glm"), + trt_method = "default", folds = 10, weights = NULL, control = lmtp_control()) { @@ -146,7 +149,13 @@ lmtp_tmle <- function(data, trt, outcome, baseline = NULL, time_vary = NULL, pb <- progressr::progressor(Task$tau*folds*2) - ratios <- cf_r(Task, learners_trt, mtp, control, pb) + if (trt_method == "default") { + ratios <- cf_r(Task, learners_trt, mtp, control, pb) + } + else { + ratios <- cf_rr(Task, learners_trt, mtp, control, pb) + } + estims <- cf_tmle(Task, "tmp_lmtp_scaled_outcome", ratios$ratios, learners_outcome, control, pb) @@ -155,9 +164,10 @@ lmtp_tmle <- function(data, trt, outcome, baseline = NULL, time_vary = NULL, estimator = "TMLE", m = list(natural = estims$natural, shifted = estims$shifted), r = ratios$ratios, + cumulated = trt_method == "riesz", tau = Task$tau, folds = Task$folds, - id = Task$id, + id = Task$natural$lmtp_id, outcome_type = Task$outcome_type, bounds = Task$bounds, weights = Task$weights, @@ -265,7 +275,6 @@ lmtp_tmle <- function(data, trt, outcome, baseline = NULL, time_vary = NULL, lmtp_sdr <- function(data, trt, outcome, baseline = NULL, time_vary = NULL, cens = NULL, shift = NULL, shifted = NULL, k = Inf, mtp = FALSE, - # intervention_type = c("static", "dynamic", "mtp"), outcome_type = c("binomial", "continuous", "survival"), id = NULL, bounds = NULL, learners_outcome = c("mean", "glm"), @@ -329,9 +338,10 @@ lmtp_sdr <- function(data, trt, outcome, baseline = NULL, time_vary = NULL, estimator = "SDR", m = list(natural = estims$natural, shifted = estims$shifted), r = ratios$ratios, + cumulated = FALSE, tau = Task$tau, folds = Task$folds, - id = Task$id, + id = Task$natural$lmtp_id, outcome_type = Task$outcome_type, bounds = Task$bounds, weights = Task$weights, @@ -536,6 +546,8 @@ lmtp_sub <- function(data, trt, outcome, baseline = NULL, time_vary = NULL, cens #' @param learners \[\code{character}\]\cr A vector of \code{mlr3superlearner} algorithms for estimation #' of the outcome regression. Default is \code{c("mean", "glm")}. #' \bold{Only include candidate learners capable of binary classification}. +#' @param trt_method \[\code{character}\]\cr +#' Method for estimating treatment assignment mechanism (default or riesz) #' @param folds \[\code{integer(1)}\]\cr #' The number of folds to be used for cross-fitting. #' @param weights \[\code{numeric(nrow(data))}\]\cr @@ -568,10 +580,10 @@ lmtp_sub <- function(data, trt, outcome, baseline = NULL, time_vary = NULL, cens #' @example inst/examples/ipw-ex.R lmtp_ipw <- function(data, trt, outcome, baseline = NULL, time_vary = NULL, cens = NULL, shift = NULL, shifted = NULL, mtp = FALSE, - # intervention_type = c("static", "dynamic", "mtp"), k = Inf, id = NULL, outcome_type = c("binomial", "continuous", "survival"), learners = c("mean", "glm"), + trt_method = "default", folds = 10, weights = NULL, control = lmtp_control()) { @@ -619,16 +631,21 @@ lmtp_ipw <- function(data, trt, outcome, baseline = NULL, time_vary = NULL, cens ) pb <- progressr::progressor(Task$tau*folds) - - ratios <- cf_r(Task, learners, mtp, control, pb) - - theta_ipw( - eta = list( - r = matrix( + + if (trt_method == "default") { + ratios <- cf_r(Task, learners, mtp, control, pb) + ratios$ratios <- matrix( t(apply(ratios$ratios, 1, cumprod)), nrow = nrow(ratios$ratios), ncol = ncol(ratios$ratios) - ), + ) + } else { + ratios <- cf_rr(Task, learners, mtp, control, pb) + } + + theta_ipw( + eta = list( + r = ratios$ratios, y = if (Task$survival) { convert_to_surv(data[[final_outcome(outcome)]]) } else { diff --git a/R/gcomp.R b/R/gcomp.R index 28ef367..088b826 100644 --- a/R/gcomp.R +++ b/R/gcomp.R @@ -17,7 +17,7 @@ cf_sub <- function(Task, outcome, learners, control, pb) { out <- future::value(out) list( - m = recombine_outcome(out, "m", Task$folds), + m = recombine_outcome(out, "m", task$folds), fits = lapply(out, function(x) x[["fits"]]) ) } @@ -74,7 +74,7 @@ estimate_sub <- function(natural, shifted, trt, outcome, node_list, cens, risk, natural$train[!rt, pseudo] <- 0 m[!rv, t] <- 0 - pb() + progress_bar() } list(m = m, fits = fits) diff --git a/R/lmtp_options.R b/R/lmtp_options.R new file mode 100644 index 0000000..d1240ba --- /dev/null +++ b/R/lmtp_options.R @@ -0,0 +1,29 @@ +#' Set LMTP Estimation Parameters +#' +#' @param .trim \[\code{numeric(1)}\]\cr +#' Determines the amount the density ratios should be trimmed. +#' The default is 0.999, trimming the density ratios greater than the 0.999 percentile +#' to the 0.999 percentile. A value of 1 indicates no trimming. +#' @param .learners_outcome_folds \[\code{integer(1)}\]\cr +#' The number of cross-validation folds for \code{learners_outcome}. +#' @param .learners_trt_folds \[\code{integer(1)}\]\cr +#' The number of cross-validation folds for \code{learners_trt}. +#' @param .return_full_fits \[\code{logical(1)}\]\cr +#' Return full \code{mlr3superlearner} fits? Default is \code{FALSE}. +#' +#' @return A list of parameters controlling the estimation procedure. +#' @export +#' +#' @examples +#' lmtp_control(.trim = 0.975) +lmtp_control <- function(...) { + change <- list(...) + control <- list(.trim = 0.999, + .learners_outcome_folds = NULL, + .learners_trt_folds = NULL, + .return_full_fits = FALSE) + if (length(change) == 0) return(control) + change <- change[names(change) %in% names(control)] + control[names(change)] <- change + control +} diff --git a/R/riesz_representer.R b/R/riesz_representer.R new file mode 100644 index 0000000..d2fad95 --- /dev/null +++ b/R/riesz_representer.R @@ -0,0 +1,69 @@ +cf_rr <- function(task, learners, mtp, control, progress_bar) { + out <- list() + for (fold in seq_along(task$folds)) { + out[[fold]] <- future::future({ + estimate_rr(get_folded_data(task$natural, task$folds, fold), + get_folded_data(task$shifted, task$folds, fold), + task$trt, + task$cens, + task$risk, + task$tau, + task$node_list$trt, + learners, + mtp, + control, + progress_bar) + }, + seed = TRUE) + } + + recombine_ratios(future::value(out), task$folds) +} + + +estimate_rr <- function(natural, shifted, trt, cens, risk, tau, node_list, learners, mtp, control, progress_bar) { + representers <- matrix(nrow = nrow(natural$valid), ncol = tau) + fits <- list() + + for (t in 1:tau) { + jrt <- censored(natural$train, cens, t)$j + drt <- at_risk(natural$train, risk, t) + irv <- censored(natural$valid, cens, t)$i + jrv <- censored(natural$valid, cens, t)$j + drv <- at_risk(natural$valid, risk, t) + + trt_t <- ifelse(length(trt) > 1, trt[t], trt) + + frv <- followed_rule(natural$valid[[trt_t]], shifted$valid[[trt_t]], mtp) + + vars <- c(node_list[[t]], cens[[t]]) + + conditional_indicator_train <- matrix(1, ncol = 1, nrow = nrow(natural$train)) + conditional_indicator_valid <- matrix(1, ncol = 1, nrow = nrow(natural$valid)) + fit <- run_riesz_ensemble( + learners, + natural$train[jrt & drt, vars, drop = FALSE], + shifted$train[jrt & drt, vars, drop = FALSE], + conditional_indicator_train[jrt & drt,,drop = FALSE], + natural$valid[jrv & drv, vars, drop = FALSE], + shifted$valid[jrv & drv, vars, drop = FALSE], + conditional_indicator_valid[jrv & drv, ,drop = FALSE], + folds = control$.learners_trt_folds + ) + + if (control$.return_full_fits) { + fits[[t]] <- fit + } else { + fits[[t]] <- extract_sl_weights(fit) + } + + pred <- matrix(-999L, nrow = nrow(natural$valid), ncol = 1) + pred[jrv & drv, ] <- fit$predictions + + representers[, t] <- pred + + progress_bar() + } + + list(ratios = representers, fits = fits) +} diff --git a/R/sdr.R b/R/sdr.R index e221532..fddd3d3 100644 --- a/R/sdr.R +++ b/R/sdr.R @@ -16,18 +16,21 @@ cf_sdr <- function(Task, outcome, ratios, learners, control, pb) { out <- future::value(out) - list(natural = recombine_outcome(out, "natural", Task$folds), - shifted = recombine_outcome(out, "shifted", Task$folds), + list(natural = recombine_outcome(out, "natural", task$folds), + shifted = recombine_outcome(out, "shifted", task$folds), fits = lapply(out, function(x) x[["fits"]])) } estimate_sdr <- function(natural, shifted, trt, outcome, node_list, cens, risk, tau, outcome_type, ratios, learners, control, pb) { - m_natural_train <- m_shifted_train <- - cbind(matrix(nrow = nrow(natural$train), ncol = tau), natural$train[[outcome]]) - m_natural_valid <- m_shifted_valid <- - cbind(matrix(nrow = nrow(natural$valid), ncol = tau), natural$valid[[outcome]]) + m_natural_train <- m_shifted_train <- cbind(matrix(nrow = nrow(natural$train), + ncol = tau), + natural$train[[outcome]]) + + m_natural_valid <- m_shifted_valid <- cbind(matrix(nrow = nrow(natural$valid), + ncol = tau), + natural$valid[[outcome]]) fits <- vector("list", length = tau) @@ -57,12 +60,13 @@ estimate_sdr <- function(natural, shifted, trt, outcome, node_list, cens, risk, } if (t < tau) { - tmp <- transform_sdr(compute_weights(ratios, t + 1, tau), - t, tau, - m_shifted_train, - m_natural_train) + densratio <- transform_sdr(compute_weights(ratios, t + 1, tau), + t, + tau, + m_shifted_train, + m_natural_train) - natural$train[, pseudo] <- shifted$train[, pseudo] <- tmp + natural$train[, pseudo] <- shifted$train[, pseudo] <- densratio fit <- run_ensemble(natural$train[i & rt, c("lmtp_id", vars, pseudo)], pseudo, @@ -100,7 +104,7 @@ estimate_sdr <- function(natural, shifted, trt, outcome, node_list, cens, risk, m_natural_valid[!rv, t] <- 0 m_shifted_valid[!rv, t] <- 0 - pb() + progress_bar() } list(natural = m_natural_valid, diff --git a/R/sl_riesz.R b/R/sl_riesz.R new file mode 100644 index 0000000..16e07a3 --- /dev/null +++ b/R/sl_riesz.R @@ -0,0 +1,31 @@ +riesz_superlearner_weights <- function(learners, task_valid) { + risks <- lapply(learners, function(x) { + x$loss(task_valid) + }) + + weights <- numeric(length(learners)) + weights[which.min(risks)] <- 1 + list(weights = weights, risk = risks) +} + +#' @importFrom SuperRiesz super_riesz +run_riesz_ensemble <- function(learners, natural_train, shifted_train, conditional_train, + natural_valid, shifted_valid, conditional_valid, folds) { + + if(is.null(folds)) folds <- 5 + sl <- SuperRiesz::super_riesz( + natural_train, + list(shifted = shifted_train, weight = data.frame(weight = conditional_train / mean(conditional_train))), + library = learners, + folds = folds, + m = \(alpha, data) alpha(data("shifted")) * data("weight")[,1] + ) + predictions = predict(sl, shifted_valid) * mean(conditional_valid[, 1]) + + list( + predictions = predictions, + fits = sl, + coef = sl$weights, + risk = sl$risk + ) +} diff --git a/R/theta.R b/R/theta.R index 09bdcfd..e0a56f0 100644 --- a/R/theta.R +++ b/R/theta.R @@ -58,15 +58,22 @@ theta_ipw <- function(eta) { out } -eif <- function(r, tau, shifted, natural) { +eif <- function(r, cumulated, tau, shifted, natural) { natural[is.na(natural)] <- -999 shifted[is.na(shifted)] <- -999 m <- shifted[, 2:(tau + 1), drop = FALSE] - natural[, 1:tau, drop = FALSE] - rowSums(compute_weights(r, 1, tau) * m, na.rm = TRUE) + shifted[, 1] + if(cumulated == TRUE) { + weights <- r + } + else { + weights <- compute_weights(r, 1, tau) + } + rowSums(weights * m, na.rm = TRUE) + shifted[, 1] } theta_dr <- function(eta, augmented = FALSE) { inflnce <- eif(r = eta$r, + cumulated = eta$cumulated, tau = eta$tau, shifted = eta$m$shifted, natural = eta$m$natural) diff --git a/R/tmle.R b/R/tmle.R index ca84a0f..168f2bf 100644 --- a/R/tmle.R +++ b/R/tmle.R @@ -1,10 +1,13 @@ -cf_tmle <- function(Task, outcome, ratios, learners, control, pb) { +cf_tmle <- function(Task, outcome, cumulated, ratios, learners, control, pb) { out <- vector("list", length = length(Task$folds)) - ratios <- matrix(t(apply(ratios, 1, cumprod)), - nrow = nrow(ratios), - ncol = ncol(ratios)) + + if(cumulated == FALSE) { + ratios <- matrix(t(apply(ratios, 1, cumprod)), + nrow = nrow(ratios), + ncol = ncol(ratios)) + } - for (fold in seq_along(Task$folds)) { + for (fold in seq_along(task$folds)) { out[[fold]] <- future::future({ estimate_tmle( get_folded_data(Task$natural, Task$folds, fold), @@ -20,11 +23,10 @@ cf_tmle <- function(Task, outcome, ratios, learners, control, pb) { out <- future::value(out) - list( - natural = recombine_outcome(out, "natural", Task$folds), - shifted = cbind(recombine_outcome(out, "shifted", Task$folds), Task$natural[["tmp_lmtp_scaled_outcome"]]), - fits = lapply(out, function(x) x[["fits"]]) - ) + list(natural = recombine_outcome(out, "natural", task$folds), + shifted = cbind(recombine_outcome(out, "shifted", task$folds), + task$natural[["tmp_lmtp_scaled_outcome"]]), + fits = lapply(out, function(x) x[["fits"]])) } estimate_tmle <- function(natural, shifted, trt, outcome, node_list, cens, @@ -94,12 +96,10 @@ estimate_tmle <- function(natural, shifted, trt, outcome, node_list, cens, m_natural_valid[!rv, t] <- 0 m_shifted_valid[!rv, t] <- 0 - pb() + progress_bar() } - list( - natural = m_natural_valid, - shifted = m_shifted_valid, - fits = fits - ) + list(natural = m_natural_valid, + shifted = m_shifted_valid, + fits = fits) } diff --git a/README.Rmd b/README.Rmd index 2ce2e86..6b50559 100644 --- a/README.Rmd +++ b/README.Rmd @@ -35,22 +35,8 @@ For an in-depth look at the package's functionality, please consult the accompan ## Installation -**lmtp** can be installed from CRAN with: - -``` r -install.packages("lmtp") -``` - -The stable, development version can be installed from GitHub with: - -``` r -devtools::install_github("nt-williams/lmtp@devel") -``` - -A version allowing for different covariates sets for the treatment, censoring, and outcome regressions: - -``` r -devtools::install_github("nt-williams/lmtp@separate-variable-sets") +```r +remotes::install_github("nt-williams/lmtp@riesz") ``` ## What even is a modified treatment policy? diff --git a/README.md b/README.md index 8374a0f..1c15b74 100644 --- a/README.md +++ b/README.md @@ -35,7 +35,7 @@ treatment weighting estimator are provided for the sake of being thorough but their use is recommended against in favor of the TML and SDR estimators). Both binary and continuous outcomes (both with censoring) are allowed. **lmtp** is built atop the -[`SuperLearner`](https://CRAN.R-project.org/package=SuperLearner) +[`mlr3superlearner`](https://github.com/nt-williams/mlr3superlearner) package to utilize ensemble machine learning for estimation. The treatment mechanism is estimated via a density ratio classification procedure irrespective of treatment variable type providing decreased @@ -51,23 +51,8 @@ Studies](https://muse.jhu.edu/article/883479). ## Installation -**lmtp** can be installed from CRAN with: - -``` r -install.packages("lmtp") -``` - -The stable, development version can be installed from GitHub with: - -``` r -devtools::install_github("nt-williams/lmtp@devel") -``` - -A version allowing for different covariates sets for the treatment, -censoring, and outcome regressions: - ``` r -devtools::install_github("nt-williams/lmtp@separate-variable-sets") +remotes::install_github("nt-williams/lmtp@mlr3superlearner") ``` ## What even is a modified treatment policy? @@ -117,9 +102,7 @@ regimes for binary, continuous, and survival outcomes. library(lmtp) #> Loading required package: mlr3superlearner #> Loading required package: mlr3learners -#> Warning: package 'mlr3learners' was built under R version 4.2.3 #> Loading required package: mlr3 -#> Warning: package 'mlr3' was built under R version 4.2.3 # the data: 4 treatment nodes with time varying covariates and a binary outcome head(sim_t4) diff --git a/cran-comments.md b/cran-comments.md old mode 100644 new mode 100755 diff --git a/inst/examples/ipw-ex.R b/inst/examples/ipw-ex.R index 83e06aa..cdc832d 100644 --- a/inst/examples/ipw-ex.R +++ b/inst/examples/ipw-ex.R @@ -74,7 +74,7 @@ # Example 2.4 # Using the same data as examples 2.1, 2.2, and 2.3, but now treating the exposure - # as an ordered categorical variable. To account for the exposure being a + # as a categorical variable. To account for the exposure being a # factor we just need to modify the shift function (and the original data) # so as to respect this. tmp <- sim_t4 diff --git a/inst/examples/sdr-ex.R b/inst/examples/sdr-ex.R index 282120a..904324b 100644 --- a/inst/examples/sdr-ex.R +++ b/inst/examples/sdr-ex.R @@ -74,7 +74,7 @@ # Example 2.4 # Using the same data as examples 2.1, 2.2, and 2.3, but now treating the exposure - # as an ordered categorical variable. To account for the exposure being a + # as a categorical variable. To account for the exposure being a # factor we just need to modify the shift function (and the original data) # so as to respect this. tmp <- sim_t4 diff --git a/inst/examples/sub-ex.R b/inst/examples/sub-ex.R index 5537df8..4413d8d 100644 --- a/inst/examples/sub-ex.R +++ b/inst/examples/sub-ex.R @@ -72,7 +72,7 @@ # Example 2.4 # Using the same data as examples 2.1, 2.2, and 2.3, but now treating the exposure - # as an ordered categorical variable. To account for the exposure being a + # as a categorical variable. To account for the exposure being a # factor we just need to modify the shift function (and the original data) # so as to respect this. tmp <- sim_t4 diff --git a/inst/examples/tmle-ex.R b/inst/examples/tmle-ex.R index 3314d30..f5907b3 100644 --- a/inst/examples/tmle-ex.R +++ b/inst/examples/tmle-ex.R @@ -73,7 +73,7 @@ # Example 2.4 # Using the same data as examples 2.1, 2.2, and 2.3, but now treating the exposure - # as an ordered categorical variable. To account for the exposure being a + # as a categorical variable. To account for the exposure being a # factor we just need to modify the shift function (and the original data) # so as to respect this. tmp <- sim_t4 diff --git a/man/lmtp_ipw.Rd b/man/lmtp_ipw.Rd index 951dce5..c89fd42 100644 --- a/man/lmtp_ipw.Rd +++ b/man/lmtp_ipw.Rd @@ -196,7 +196,7 @@ by some amount, use \code{mtp = TRUE}}. # Example 2.4 # Using the same data as examples 2.1, 2.2, and 2.3, but now treating the exposure - # as an ordered categorical variable. To account for the exposure being a + # as a categorical variable. To account for the exposure being a # factor we just need to modify the shift function (and the original data) # so as to respect this. tmp <- sim_t4 diff --git a/man/lmtp_sdr.Rd b/man/lmtp_sdr.Rd index 0f882b4..01b3450 100644 --- a/man/lmtp_sdr.Rd +++ b/man/lmtp_sdr.Rd @@ -212,7 +212,7 @@ by some amount, use \code{mtp = TRUE}}. # Example 2.4 # Using the same data as examples 2.1, 2.2, and 2.3, but now treating the exposure - # as an ordered categorical variable. To account for the exposure being a + # as a categorical variable. To account for the exposure being a # factor we just need to modify the shift function (and the original data) # so as to respect this. tmp <- sim_t4 diff --git a/man/lmtp_sub.Rd b/man/lmtp_sub.Rd index 9809582..1a2e85c 100644 --- a/man/lmtp_sub.Rd +++ b/man/lmtp_sub.Rd @@ -185,7 +185,7 @@ continuous, or time-to-event outcomes. Supports binary, categorical, and continu # Example 2.4 # Using the same data as examples 2.1, 2.2, and 2.3, but now treating the exposure - # as an ordered categorical variable. To account for the exposure being a + # as a categorical variable. To account for the exposure being a # factor we just need to modify the shift function (and the original data) # so as to respect this. tmp <- sim_t4 diff --git a/man/lmtp_tmle.Rd b/man/lmtp_tmle.Rd index 37bbf17..7abb0b5 100644 --- a/man/lmtp_tmle.Rd +++ b/man/lmtp_tmle.Rd @@ -211,7 +211,7 @@ by some amount, use \code{mtp = TRUE}}. # Example 2.4 # Using the same data as examples 2.1, 2.2, and 2.3, but now treating the exposure - # as an ordered categorical variable. To account for the exposure being a + # as a categorical variable. To account for the exposure being a # factor we just need to modify the shift function (and the original data) # so as to respect this. tmp <- sim_t4 diff --git a/tests/testthat/test-censoring.R b/tests/testthat/test-censoring.R index c7b2ae1..9ad4184 100644 --- a/tests/testthat/test-censoring.R +++ b/tests/testthat/test-censoring.R @@ -8,6 +8,7 @@ truth <- 0.8 sub <- sw(lmtp_sub(sim_cens, A, "Y", time_vary = L, cens = C, k = 0, shift = NULL, folds = 1)) ipw <- sw(lmtp_ipw(sim_cens, A, "Y", time_vary = L, cens = C, k = 0, shift = NULL, folds = 1)) tmle <- sw(lmtp_tmle(sim_cens, A, "Y", time_vary = L, cens = C, k = 0, shift = NULL, folds = 1)) +tmle_riesz <- sw(lmtp_tmle(sim_cens, A, "Y", time_vary = L, cens = C, k = 0, shift = NULL, folds = 1, trt_method = "riesz", learners_trt = c("glm"))) sdr <- sw(lmtp_sdr(sim_cens, A, "Y", time_vary = L, cens = C, k = 0, shift = NULL, folds = 1)) # tests @@ -15,5 +16,6 @@ test_that("estimator fidelity with censoring", { expect_equal(truth, sub$theta, tolerance = 0.15) expect_equal(truth, ipw$theta, tolerance = 0.15) expect_equal(truth, tmle$theta, tolerance = 0.15) + expect_equal(truth, tmle_riesz$theta, tolerance = 0.15) expect_equal(truth, sdr$theta, tolerance = 0.15) }) diff --git a/tests/testthat/test-point_treatment.R b/tests/testthat/test-point_treatment.R index 76bbab5..806ba37 100644 --- a/tests/testthat/test-point_treatment.R +++ b/tests/testthat/test-point_treatment.R @@ -6,7 +6,7 @@ n <- 1e4 W1 <- rbinom(n, size = 1, prob = 0.5) W2 <- rbinom(n, size = 1, prob = 0.65) A <- rbinom(n, size = 1, prob = plogis(-0.4 + 0.2 * W2 + 0.15 * W1)) -Y.1 <-rbinom(n, size = 1, prob = plogis(-1 + 1 - 0.1 * W1 + 0.3 * W2)) +Y.1 <- rbinom(n, size = 1, prob = plogis(-1 + 1 - 0.1 * W1 + 0.3 * W2)) Y.0 <- rbinom(n, size = 1, prob = plogis(-1 + 0 - 0.1 * W1 + 0.3 * W2)) Y <- Y.1 * A + Y.0 * (1 - A) @@ -16,6 +16,7 @@ truth <- mean(tmp$Y.1) sub <- sw(lmtp_sub(tmp, "A", "Y", baseline = c("W1", "W2"), shift = static_binary_on, folds = 1)) ipw <- sw(lmtp_ipw(tmp, "A", "Y", baseline = c("W1", "W2"), shift = static_binary_on, folds = 1)) tmle <- sw(lmtp_tmle(tmp, "A", "Y", baseline = c("W1", "W2"), shift = static_binary_on, folds = 1)) +tmle_riesz <- sw(lmtp_tmle(tmp, "A", "Y", baseline = c("W1", "W2"), shift = static_binary_on, folds = 1, trt_method = "riesz", learners_trt = c("glm"))) sdr <- sw(lmtp_sdr(tmp, "A", "Y", baseline = c("W1", "W2"), shift = static_binary_on, folds = 1)) # tests @@ -23,5 +24,6 @@ test_that("point treatment fidelity", { expect_equal(truth, sub$theta, tolerance = 0.025) expect_equal(truth, ipw$theta, tolerance = 0.025) expect_equal(truth, tmle$theta, tolerance = 0.025) + expect_equal(truth, tmle_riesz$theta, tolerance = 0.025) expect_equal(truth, sdr$theta, tolerance = 0.025) }) diff --git a/tests/testthat/test-shifted.R b/tests/testthat/test-shifted.R index b85f90b..7c6b368 100644 --- a/tests/testthat/test-shifted.R +++ b/tests/testthat/test-shifted.R @@ -22,6 +22,12 @@ tmle <- cens, k = 0, shifted = sc, outcome_type = "binomial", folds = 2, mtp = TRUE)) +tmle_riesz <- + sw(lmtp_tmle(sim_cens, a, "Y", nodes, baseline = NULL, + cens, k = 0, shifted = sc, trt_method = "riesz", + learners_trt = c("constant"), + outcome_type = "binomial", folds = 2, intervention_type = "mtp")) + sdr <- sw(lmtp_sdr(sim_cens, a, "Y", nodes, baseline = NULL, cens, k = 0, shifted = sc, @@ -32,5 +38,6 @@ test_that("estimator fidelity with shifted data supplied", { expect_equal(truth, sub$theta, tolerance = 0.05) expect_equal(truth, ipw$theta, tolerance = 0.025) expect_equal(truth, tmle$theta, tolerance = 0.025) + expect_equal(truth, tmle_riesz$theta, tolerance = 0.025) expect_equal(truth, sdr$theta, tolerance = 0.025) }) diff --git a/tests/testthat/test-survey.R b/tests/testthat/test-survey.R index 1cf934b..90b7c78 100644 --- a/tests/testthat/test-survey.R +++ b/tests/testthat/test-survey.R @@ -27,6 +27,10 @@ ipw <- sw(lmtp_ipw(tmp, "A", "Y", baseline = c("W1", "W2"), shift = static_binar tmle <- sw(lmtp_tmle(tmp, "A", "Y", baseline = c("W1", "W2"), shift = static_binary_on, weights = wts, folds = 2)) +tmle_riesz <- sw(lmtp_tmle(tmp, "A", "Y", baseline = c("W1", "W2"), shift = static_binary_on, + trt_method = "riesz", learners_trt = c("glm"), + weights = wts, folds = 2)) + sdr <- sw(lmtp_sdr(tmp, "A", "Y", baseline = c("W1", "W2"), shift = static_binary_on, weights = wts, folds = 2)) @@ -36,5 +40,6 @@ test_that("survey weight fidelity", { expect_equal(truth, sub$theta, tolerance = 0.025) expect_equal(truth, ipw$theta, tolerance = 0.025) expect_equal(truth, tmle$theta, tolerance = 0.025) + expect_equal(truth, tmle_riesz$theta, tolerance = 0.025) expect_equal(truth, sdr$theta, tolerance = 0.025) }) diff --git a/tests/testthat/test-time_varying_treatment.R b/tests/testthat/test-time_varying_treatment.R index 5903b52..72e772f 100644 --- a/tests/testthat/test-time_varying_treatment.R +++ b/tests/testthat/test-time_varying_treatment.R @@ -23,14 +23,15 @@ d <- function(data, trt) { truth <- 0.305 -sub <- sw(lmtp_sub(tmp, a, "Y", time_vary = time_varying, shift = d, folds = 1)) ipw <- sw(lmtp_ipw(tmp, a, "Y", time_vary = time_varying, shift = d, mtp = T, folds = 1)) tmle <- sw(lmtp_tmle(tmp, a, "Y", time_vary = time_varying, shift = d, mtp = T, folds = 1)) +tmle_riesz <- sw(lmtp_tmle(tmp, a, "Y", time_vary = time_varying, shift = d, mtp = T, folds = 1, trt_method = "riesz", learners_trt = c("constant"))) sdr <- sw(lmtp_sdr(tmp, a, "Y", time_vary = time_varying, shift = d, mtp = T, folds = 1)) test_that("time varying treatment fidelity, t = 4", { expect_equal(truth, sub$theta, tolerance = 0.025) expect_equal(truth, ipw$theta, tolerance = 0.05) expect_equal(truth, tmle$theta, tolerance = 0.025) + expect_equal(truth, tmle_riesz$theta, tolerance = 0.025) expect_equal(truth, sdr$theta, tolerance = 0.025) })