From 8d4ccd8d17e7e950508e0208576f0d6d57943e44 Mon Sep 17 00:00:00 2001 From: nt-williams Date: Fri, 26 Apr 2024 10:20:46 -0700 Subject: [PATCH 01/26] Modified TMLE for computationally effiecient bootstrap --- DESCRIPTION | 2 +- R/estimators.R | 38 +++++++++++++++++++++++++++--- R/gcomp.R | 14 ++++++----- R/lmtp_control.R | 10 ++++++-- R/theta.R | 35 ++++++++++++++++++++++++++++ R/tmle2.R | 57 +++++++++++++++++++++++++++++++++++++++++++++ man/lmtp_control.Rd | 8 ++++++- man/lmtp_tmle.Rd | 8 ++++++- 8 files changed, 158 insertions(+), 14 deletions(-) create mode 100644 R/tmle2.R diff --git a/DESCRIPTION b/DESCRIPTION index 14c2e99..977052d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -24,7 +24,7 @@ License: AGPL-3 Encoding: UTF-8 LazyData: true Roxygen: list(markdown = TRUE) -RoxygenNote: 7.2.3 +RoxygenNote: 7.3.1 Imports: stats, nnls, diff --git a/R/estimators.R b/R/estimators.R index 8f4b7c5..e3e80fd 100644 --- a/R/estimators.R +++ b/R/estimators.R @@ -41,6 +41,9 @@ #' @param mtp \[\code{logical(1)}\]\cr #' Is the intervention of interest a modified treatment policy? #' Default is \code{FALSE}. If treatment variables are continuous this should be \code{TRUE}. +#' @param boot \[\code{logical(1)}\]\cr +#' Compute standard errors using the bootstrap? Default is \code{TRUE}. If \code{FALSE}, standard +#' errors will be calculated using the empirical variance of the efficient influence function. #' @param outcome_type \[\code{character(1)}\]\cr #' Outcome variable type (i.e., continuous, binomial, survival). #' @param id \[\code{character(1)}\]\cr @@ -77,7 +80,8 @@ #' \item{standard_error}{The estimated standard error of the LMTP effect.} #' \item{low}{Lower bound of the 95% confidence interval of the LMTP effect.} #' \item{high}{Upper bound of the 95% confidence interval of the LMTP effect.} -#' \item{eif}{The estimated, un-centered, influence function of the estimate.} +#' \item{eif}{The estimated, un-centered, influence function of the estimate, +#' \code{NULL} if \code{boot = TRUE}.} #' \item{shift}{The shift function specifying the treatment policy of interest.} #' \item{outcome_reg}{An n x Tau + 1 matrix of outcome regression predictions. #' The mean of the first column is used for calculating theta.} @@ -92,7 +96,8 @@ #' @export lmtp_tmle <- function(data, trt, outcome, baseline = NULL, time_vary = NULL, cens = NULL, shift = NULL, shifted = NULL, k = Inf, - mtp = FALSE, outcome_type = c("binomial", "continuous", "survival"), + mtp = FALSE, boot = TRUE, + outcome_type = c("binomial", "continuous", "survival"), id = NULL, bounds = NULL, learners_outcome = "SL.glm", learners_trt = "SL.glm", @@ -125,6 +130,7 @@ lmtp_tmle <- function(data, trt, outcome, baseline = NULL, time_vary = NULL, checkmate::assertNumber(control$.bound) checkmate::assertNumber(control$.trim, upper = 1) checkmate::assertLogical(control$.return_full_fits, len = 1) + checkmate::assertLogical(boot, len = 1) task <- lmtp_task$new( data = data, @@ -146,6 +152,32 @@ lmtp_tmle <- function(data, trt, outcome, baseline = NULL, time_vary = NULL, pb <- progressr::progressor(task$tau*folds*2) + if (isTRUE(boot)) { + ratios <- cf_r(task, learners_trt, mtp, control, pb) + Qn <- cf_sub(task, "tmp_lmtp_scaled_outcome", learners_outcome, control, pb) + Qnb_eps <- cf_tmle2(task, ratios$ratios, Qn, control) + + ans <- theta_boot( + list( + estimator = "TMLE", + m = Qnb_eps$psi, + r = ratios$ratios, + boots = Qnb_eps$booted, + tau = task$tau, + folds = task$folds, + id = task$id, + outcome_type = task$outcome_type, + bounds = task$bounds, + weights = task$weights, + shift = if (is.null(shifted)) deparse(substitute((shift))) else NULL, + fits_m = Qn$fits, + fits_r = ratios$fits, + outcome_type = task$outcome_type + ) + ) + return(ans) + } + ratios <- cf_r(task, learners_trt, mtp, control, pb) estims <- cf_tmle(task, "tmp_lmtp_scaled_outcome", @@ -484,7 +516,7 @@ lmtp_sub <- function(data, trt, outcome, baseline = NULL, time_vary = NULL, cens theta_sub( eta = list( - m = estims$m, + m = estims$ms, outcome_type = task$outcome_type, bounds = task$bounds, folds = task$folds, diff --git a/R/gcomp.R b/R/gcomp.R index ad0b718..a5feeb9 100644 --- a/R/gcomp.R +++ b/R/gcomp.R @@ -23,15 +23,15 @@ cf_sub <- function(task, outcome, learners, control, pb) { out <- future::value(out) list( - m = recombine_outcome(out, "m", task$folds), + ms = recombine_outcome(out, "ms", task$folds), + mn = recombine_outcome(out, "mn", task$folds), fits = lapply(out, function(x) x[["fits"]]) ) } estimate_sub <- function(natural, shifted, trt, outcome, node_list, cens, risk, tau, outcome_type, learners, control, pb) { - - m <- matrix(nrow = nrow(natural$valid), ncol = tau) + ms <- mn <- matrix(nrow = nrow(natural$valid), ncol = tau) fits <- vector("list", length = tau) for (t in tau:1) { @@ -79,13 +79,15 @@ estimate_sub <- function(natural, shifted, trt, outcome, node_list, cens, risk, under_shift_valid[, trt_t] <- shifted$valid[jv & rv, trt_t] natural$train[jt & rt, pseudo] <- bound(SL_predict(fit, under_shift_train), 1e-05) - m[jv & rv, t] <- bound(SL_predict(fit, under_shift_valid), 1e-05) + ms[jv & rv, t] <- bound(SL_predict(fit, under_shift_valid), 1e-05) + mn[jv & rv, t] <- bound(SL_predict(fit, natural$valid[jv & rv, vars]), 1e-05) natural$train[!rt, pseudo] <- 0 - m[!rv, t] <- 0 + ms[!rv, t] <- 0 + mn[!rv, t] <- 0 pb() } - list(m = m, fits = fits) + list(ms = ms, mn = mn, fits = fits) } diff --git a/R/lmtp_control.R b/R/lmtp_control.R index 4c8ad6b..86a3237 100644 --- a/R/lmtp_control.R +++ b/R/lmtp_control.R @@ -13,6 +13,8 @@ #' The number of cross-validation folds for \code{learners_trt}. #' @param .return_full_fits \[\code{logical(1)}\]\cr #' Return full SuperLearner fits? Default is \code{FALSE}, return only SuperLearner weights. +#' @param .B description +#' @param .boot_seed description #' #' @return A list of parameters controlling the estimation procedure. #' @export @@ -23,10 +25,14 @@ lmtp_control <- function(.bound = 1e5, .trim = 0.999, .learners_outcome_folds = 10, .learners_trt_folds = 10, - .return_full_fits = FALSE) { + .return_full_fits = FALSE, + .B = 1000, + .boot_seed = NULL) { list(.bound = .bound, .trim = .trim, .learners_outcome_folds = .learners_outcome_folds, .learners_trt_folds = .learners_trt_folds, - .return_full_fits = .return_full_fits) + .return_full_fits = .return_full_fits, + .B = .B, + .boot_seed = .boot_seed) } diff --git a/R/theta.R b/R/theta.R index 66582c5..56a156f 100644 --- a/R/theta.R +++ b/R/theta.R @@ -111,3 +111,38 @@ theta_dr <- function(eta, augmented = FALSE) { class(out) <- "lmtp" out } + +theta_boot <- function(eta) { + theta <- weighted.mean(eta$m[, 1], eta$weights) + + if (eta$outcome_type == "continuous") { + theta <- rescale_y_continuous(theta, eta$bounds) + } + + # NEED TO FIGURE OUT HOW THIS WOULD WORK WITH CLUSTERING + se <- sqrt(var(eta$boots)) + ci_low <- theta - (qnorm(0.975) * se) + ci_high <- theta + (qnorm(0.975) * se) + + out <- list( + estimator = eta$estimator, + theta = theta, + standard_error = se, + low = ci_low, + high = ci_high, + id = eta$id, + shift = eta$shift, + outcome_reg = switch( + eta$outcome_type, + continuous = rescale_y_continuous(eta$m, eta$bounds), + binomial = eta$m + ), + density_ratios = eta$r, + fits_m = eta$fits_m, + fits_r = eta$fits_r, + outcome_type = eta$outcome_type + ) + + class(out) <- "lmtp" + out +} diff --git a/R/tmle2.R b/R/tmle2.R new file mode 100644 index 0000000..87517ff --- /dev/null +++ b/R/tmle2.R @@ -0,0 +1,57 @@ +cf_tmle2 <- function(task, ratios, m_init, control) { + psi <- estimate_tmle2(task$natural, + task$cens, + task$risk, + task$tau, + m_init$mn, + m_init$ms, + ratios, + task$weights) + + if (!is.null(control$.boot_seed)) set.seed(control$.boot_seed) + + boots <- replicate(control$.B, + sample(1:nrow(task$natural), nrow(task$natural), replace = TRUE), + simplify = FALSE) + + Qnb <- sapply(boots, function(i) { + psis <- estimate_tmle2(task$natural[i, , drop = FALSE], + task$cens, + task$risk, + task$tau, + m_init$mn[i, , drop = FALSE], + m_init$ms[i, , drop = FALSE], + ratios[i, , drop = FALSE], + task$weights[i]) + weighted.mean(psis[,1], task$weights[i]) + }) + + list(psi = psi, booted = Qnb) +} + +estimate_tmle2 <- function(data, cens, risk, tau, mn, ms, ratios, weights) { + m_eps <- matrix(nrow = nrow(data), ncol = tau + 1) + m_eps[, tau + 1] <- data$tmp_lmtp_scaled_outcome + + fits <- vector("list", length = tau) + for (t in tau:1) { + i <- censored(data, cens, t)$i + j <- censored(data, cens, t)$j + r <- at_risk(data, risk, t) + + wts <- ratios[i & r, t] * weights[i & r] + + fit <- sw( + glm( + m_eps[i & r, t + 1] ~ offset(qlogis(mn[i & r, t])), + weights = wts, + family = "binomial" + ) + ) + + m_eps[j & r, t] <- bound(plogis(qlogis(ms[j & r, t]) + coef(fit))) + m_eps[!r, t] <- 0 + } + + m_eps +} diff --git a/man/lmtp_control.Rd b/man/lmtp_control.Rd index 2e1d10a..8172476 100644 --- a/man/lmtp_control.Rd +++ b/man/lmtp_control.Rd @@ -9,7 +9,9 @@ lmtp_control( .trim = 0.999, .learners_outcome_folds = 10, .learners_trt_folds = 10, - .return_full_fits = FALSE + .return_full_fits = FALSE, + .B = 1000, + .boot_seed = NULL ) } \arguments{ @@ -30,6 +32,10 @@ The number of cross-validation folds for \code{learners_trt}.} \item{.return_full_fits}{[\code{logical(1)}]\cr Return full SuperLearner fits? Default is \code{FALSE}, return only SuperLearner weights.} + +\item{.B}{description} + +\item{.boot_seed}{description} } \value{ A list of parameters controlling the estimation procedure. diff --git a/man/lmtp_tmle.Rd b/man/lmtp_tmle.Rd index f790bab..1d19454 100644 --- a/man/lmtp_tmle.Rd +++ b/man/lmtp_tmle.Rd @@ -15,6 +15,7 @@ lmtp_tmle( shifted = NULL, k = Inf, mtp = FALSE, + boot = TRUE, outcome_type = c("binomial", "continuous", "survival"), id = NULL, bounds = NULL, @@ -73,6 +74,10 @@ all time points.} Is the intervention of interest a modified treatment policy? Default is \code{FALSE}. If treatment variables are continuous this should be \code{TRUE}.} +\item{boot}{[\code{logical(1)}]\cr +Compute standard errors using the bootstrap? Default is \code{TRUE}. If \code{FALSE}, standard +errors will be calculated using the empirical variance of the efficient influence function.} + \item{outcome_type}{[\code{character(1)}]\cr Outcome variable type (i.e., continuous, binomial, survival).} @@ -108,7 +113,8 @@ A list of class \code{lmtp} containing the following components: \item{standard_error}{The estimated standard error of the LMTP effect.} \item{low}{Lower bound of the 95\% confidence interval of the LMTP effect.} \item{high}{Upper bound of the 95\% confidence interval of the LMTP effect.} -\item{eif}{The estimated, un-centered, influence function of the estimate.} +\item{eif}{The estimated, un-centered, influence function of the estimate, +\code{NULL} if \code{boot = TRUE}.} \item{shift}{The shift function specifying the treatment policy of interest.} \item{outcome_reg}{An n x Tau + 1 matrix of outcome regression predictions. The mean of the first column is used for calculating theta.} From a73f4866aa8409b710d34bcddae73d8778fa9255 Mon Sep 17 00:00:00 2001 From: nt-williams Date: Fri, 26 Apr 2024 10:28:03 -0700 Subject: [PATCH 02/26] Pre-allocated for loop with future for bootstrap --- R/tmle2.R | 31 ++++++++++++++++++------------- 1 file changed, 18 insertions(+), 13 deletions(-) diff --git a/R/tmle2.R b/R/tmle2.R index 87517ff..75bd45b 100644 --- a/R/tmle2.R +++ b/R/tmle2.R @@ -14,19 +14,24 @@ cf_tmle2 <- function(task, ratios, m_init, control) { sample(1:nrow(task$natural), nrow(task$natural), replace = TRUE), simplify = FALSE) - Qnb <- sapply(boots, function(i) { - psis <- estimate_tmle2(task$natural[i, , drop = FALSE], - task$cens, - task$risk, - task$tau, - m_init$mn[i, , drop = FALSE], - m_init$ms[i, , drop = FALSE], - ratios[i, , drop = FALSE], - task$weights[i]) - weighted.mean(psis[,1], task$weights[i]) - }) - - list(psi = psi, booted = Qnb) + Qnb <- vector("list", control$.B) + for (i in seq_along(boots)) { + ib <- boots[[i]] + Qnb[[i]] <- future::future({ + psis <- estimate_tmle2(task$natural[ib, , drop = FALSE], + task$cens, + task$risk, + task$tau, + m_init$mn[ib, , drop = FALSE], + m_init$ms[ib, , drop = FALSE], + ratios[ib, , drop = FALSE], + task$weights[ib]) + weighted.mean(psis[, 1], task$weights[ib]) + }, + seed = TRUE) + } + + list(psi = psi, booted = unlist(future::value(Qnb))) } estimate_tmle2 <- function(data, cens, risk, tau, mn, ms, ratios, weights) { From 327c0bf1ceff17d7ea4c9e96d282b8ae70797ff8 Mon Sep 17 00:00:00 2001 From: nt-williams Date: Fri, 26 Apr 2024 10:32:35 -0700 Subject: [PATCH 03/26] Updating tests --- tests/testthat/test-censoring.R | 2 +- tests/testthat/test-dynamic.R | 6 +++--- tests/testthat/test-point_treatment.R | 2 +- tests/testthat/test-shifted.R | 2 +- tests/testthat/test-survey.R | 2 +- tests/testthat/test-time_varying_treatment.R | 2 +- 6 files changed, 8 insertions(+), 8 deletions(-) diff --git a/tests/testthat/test-censoring.R b/tests/testthat/test-censoring.R index c7b2ae1..3b1547d 100644 --- a/tests/testthat/test-censoring.R +++ b/tests/testthat/test-censoring.R @@ -7,7 +7,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 <- sw(lmtp_tmle(sim_cens, A, "Y", time_vary = L, cens = C, k = 0, shift = NULL, boot = F, folds = 1)) sdr <- sw(lmtp_sdr(sim_cens, A, "Y", time_vary = L, cens = C, k = 0, shift = NULL, folds = 1)) # tests diff --git a/tests/testthat/test-dynamic.R b/tests/testthat/test-dynamic.R index 94bf0c2..27e7e51 100644 --- a/tests/testthat/test-dynamic.R +++ b/tests/testthat/test-dynamic.R @@ -32,18 +32,18 @@ cens <- c("C1", "C2") nodes <- list(c(NULL), c("L")) # truth = 0.308 -tml.stc <- sw(lmtp_tmle(sim, a, "Y", baseline, nodes, cens, shift = static_binary_on, folds = 1)) +tml.stc <- sw(lmtp_tmle(sim, a, "Y", baseline, nodes, cens, shift = static_binary_on, folds = 1, boot = F)) # truth = 0.528 sdr.stc <- sw(lmtp_sdr(sim, a, "Y", baseline, nodes, cens, shift = static_binary_off, folds = 1)) # truth = 0.433 -tml.tv <- sw(lmtp_tmle(sim, a, "Y", baseline, nodes, cens, shift = time_vary_on, folds = 1)) +tml.tv <- sw(lmtp_tmle(sim, a, "Y", baseline, nodes, cens, shift = time_vary_on, folds = 1, boot = F)) sdr.tv <- sw(lmtp_sdr(sim, a, "Y", baseline, nodes, cens, shift = time_vary_on, folds = 1)) # time varying and covariate dynamic # truth = 0.345 -tml.dyn <- sw(lmtp_tmle(sim, a, "Y", baseline, nodes, cens, shift = dynamic_vec, folds = 1)) +tml.dyn <- sw(lmtp_tmle(sim, a, "Y", baseline, nodes, cens, shift = dynamic_vec, folds = 1, boot = F)) sdr.dyn <- sw(lmtp_sdr(sim, a, "Y", baseline, nodes, cens, shift = dynamic_vec, folds = 1)) test_that("Dynamic intervention fidelity", { diff --git a/tests/testthat/test-point_treatment.R b/tests/testthat/test-point_treatment.R index 638e103..5f2b8d3 100644 --- a/tests/testthat/test-point_treatment.R +++ b/tests/testthat/test-point_treatment.R @@ -15,7 +15,7 @@ truth <- mean(tmp$Y.1) sub <- lmtp_sub(tmp, "A", "Y", baseline = c("W1", "W2"), shift = static_binary_on, folds = 1) ipw <- lmtp_ipw(tmp, "A", "Y", baseline = c("W1", "W2"), shift = static_binary_on, folds = 1) -tmle <- lmtp_tmle(tmp, "A", "Y", baseline = c("W1", "W2"), shift = static_binary_on, folds = 1) +tmle <- lmtp_tmle(tmp, "A", "Y", baseline = c("W1", "W2"), shift = static_binary_on, folds = 1, boot = F) sdr <- lmtp_sdr(tmp, "A", "Y", baseline = c("W1", "W2"), shift = static_binary_on, folds = 1) # tests diff --git a/tests/testthat/test-shifted.R b/tests/testthat/test-shifted.R index b85f90b..161af27 100644 --- a/tests/testthat/test-shifted.R +++ b/tests/testthat/test-shifted.R @@ -20,7 +20,7 @@ ipw <- tmle <- sw(lmtp_tmle(sim_cens, a, "Y", nodes, baseline = NULL, cens, k = 0, shifted = sc, - outcome_type = "binomial", folds = 2, mtp = TRUE)) + outcome_type = "binomial", folds = 2, mtp = TRUE, boot = F)) sdr <- sw(lmtp_sdr(sim_cens, a, "Y", nodes, baseline = NULL, diff --git a/tests/testthat/test-survey.R b/tests/testthat/test-survey.R index 2cdd6d5..b8021bf 100644 --- a/tests/testthat/test-survey.R +++ b/tests/testthat/test-survey.R @@ -25,7 +25,7 @@ ipw <- lmtp_ipw(tmp, "A", "Y", baseline = c("W1", "W2"), shift = static_binary_o weights = wts, folds = 2) tmle <- lmtp_tmle(tmp, "A", "Y", baseline = c("W1", "W2"), shift = static_binary_on, - weights = wts, folds = 2) + weights = wts, folds = 2, boot = F) sdr <- lmtp_sdr(tmp, "A", "Y", baseline = c("W1", "W2"), shift = static_binary_on, weights = wts, folds = 2) diff --git a/tests/testthat/test-time_varying_treatment.R b/tests/testthat/test-time_varying_treatment.R index 7e7cd1e..4b289e0 100644 --- a/tests/testthat/test-time_varying_treatment.R +++ b/tests/testthat/test-time_varying_treatment.R @@ -25,7 +25,7 @@ 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 <- sw(lmtp_tmle(tmp, a, "Y", time_vary = time_varying, shift = d, mtp = T, folds = 1, boot = F)) 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", { From 6e63f60b70db9a27570cebcd3e21d75a5dfdb0c6 Mon Sep 17 00:00:00 2001 From: nt-williams Date: Fri, 26 Apr 2024 10:33:47 -0700 Subject: [PATCH 04/26] Updated NEWS --- NEWS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/NEWS.md b/NEWS.md index f512846..99c17a0 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,6 +5,7 @@ - Can now estimate the effects of simultaneous interventions on multiple variables. - New pre-packaged shift function, `ipsi()` for estimating IPSI effects using the risk ratio. - `lmtp_control()` now replaces extra estimator arguments. +- Bootstrap for TMLE with the `boot` argument using a modified TMLE algorithm (https://arxiv.org/abs/1810.03030). ### Buf Fixes From 2a72e6640fbd0dc8300cf2b05618a223791682ee Mon Sep 17 00:00:00 2001 From: nt-williams Date: Fri, 26 Apr 2024 10:37:31 -0700 Subject: [PATCH 05/26] Todo notes --- R/theta.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/theta.R b/R/theta.R index 56a156f..a0065c8 100644 --- a/R/theta.R +++ b/R/theta.R @@ -112,6 +112,7 @@ theta_dr <- function(eta, augmented = FALSE) { out } +# TODO: NEED TO SAVE THE SEED FOR THE REPLICATES AND THE BOOTED ESTIMATES FOR ESTIMATNG CONTRASTS theta_boot <- function(eta) { theta <- weighted.mean(eta$m[, 1], eta$weights) @@ -119,7 +120,7 @@ theta_boot <- function(eta) { theta <- rescale_y_continuous(theta, eta$bounds) } - # NEED TO FIGURE OUT HOW THIS WOULD WORK WITH CLUSTERING + # TODO: NEED TO FIGURE OUT HOW THIS WOULD WORK WITH CLUSTERING se <- sqrt(var(eta$boots)) ci_low <- theta - (qnorm(0.975) * se) ci_high <- theta + (qnorm(0.975) * se) From 6b1a70745fe1a02c443f9656df10b6e8630bd477 Mon Sep 17 00:00:00 2001 From: nt-williams Date: Thu, 27 Jun 2024 15:12:09 -0700 Subject: [PATCH 06/26] Basic implementation for lmtp_survival --- R/lmtp_survival.R | 55 +++++++++++++++++++++++++++++++++++++++++++++++ R/print.R | 7 ++++++ 2 files changed, 62 insertions(+) create mode 100644 R/lmtp_survival.R diff --git a/R/lmtp_survival.R b/R/lmtp_survival.R new file mode 100644 index 0000000..60730be --- /dev/null +++ b/R/lmtp_survival.R @@ -0,0 +1,55 @@ +lmtp_survival <- function(data, trt, outcomes, baseline = NULL, time_vary = NULL, + cens = NULL, shift = NULL, shifted = NULL, estimator = c("lmtp_tmle", "lmtp_sdr"), k = Inf, + mtp = FALSE, + id = NULL, + learners_outcome = "SL.glm", + learners_trt = "SL.glm", + folds = 10, weights = NULL, + control = lmtp_control()) { + + checkmate::assertCharacter(outcomes, min.len = 2, null.ok = FALSE, unique = TRUE, any.missing = FALSE) + + estimator <- match.arg(estimator) + tau <- length(outcomes) + estimates <- vector("list", tau) + + for (t in 1:tau) { + args <- list( + data = data, + trt = trt, + outcome = outcomes[1:t], + baseline = baseline, + time_vary = time_vary, + cens = cens[1:t], + shift = shift, + shifted = shifted, + k = k, + mtp = mtp, + outcome_type = ifelse(t == 1, "binomial", "survival"), + id = id, + learners_outcome = learners_outcome, + learners_trt = learners_trt, + folds = folds, + weights = weights, + control = control + ) + + estimates[[t]] <- future::future({ + if (estimator == "lmtp_tmle") do.call(lmtp_tmle, args) + else do.call(lmtp_sdr, args) + }, + seed = TRUE) + } + + estimates <- future::value(estimates) + + estimates[[1]]$theta <- 1 - estimates[[1]]$theta + high <- estimates[[1]]$high + low <- estimates[[1]]$low + estimates[[1]]$high <- 1 - low + estimates[[1]]$low <- 1 - high + estimates[[1]]$eif <- 1 - estimates[[1]]$eif + + class(estimates) <- "lmtp_survival" + estimates +} \ No newline at end of file diff --git a/R/print.R b/R/print.R index 0fac95f..eb7bf40 100644 --- a/R/print.R +++ b/R/print.R @@ -21,3 +21,10 @@ print.lmtp_contrast <- function(x, ...) { x$vals$p.value <- format.pval(x$vals$p.value, digits = 3, eps = 0.001) print(format(x$vals, digits = 3)) } + +#' @export +print.lmtp_survival <- function(x, ...) { + out <- do.call("rbind", lapply(x, tidy)) + out$t <- 1:length(x) + out[, c(ncol(out), 1:ncol(out) - 1)] +} From c989ed030d3cbbbc317fe0971eca9b232e7b0e9c Mon Sep 17 00:00:00 2001 From: nt-williams Date: Thu, 27 Jun 2024 15:53:47 -0700 Subject: [PATCH 07/26] Isotonic project to keep estimates monotone --- DESCRIPTION | 3 ++- R/lmtp_survival.R | 21 ++++++++++++++------- R/print.R | 2 +- R/utils.R | 10 ++++++++++ 4 files changed, 27 insertions(+), 9 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index d0403f7..14eebb3 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -37,7 +37,8 @@ Imports: data.table (>= 1.13.0), checkmate (>= 2.1.0), SuperLearner, - schoolmath + schoolmath, + isotone URL: https://beyondtheate.com/, https://github.com/nt-williams/lmtp BugReports: https://github.com/nt-williams/lmtp/issues Suggests: diff --git a/R/lmtp_survival.R b/R/lmtp_survival.R index 60730be..7e69c57 100644 --- a/R/lmtp_survival.R +++ b/R/lmtp_survival.R @@ -42,14 +42,21 @@ lmtp_survival <- function(data, trt, outcomes, baseline = NULL, time_vary = NULL } estimates <- future::value(estimates) - - estimates[[1]]$theta <- 1 - estimates[[1]]$theta - high <- estimates[[1]]$high - low <- estimates[[1]]$low - estimates[[1]]$high <- 1 - low - estimates[[1]]$low <- 1 - high - estimates[[1]]$eif <- 1 - estimates[[1]]$eif + estimates <- fix_surv_time1(estimates) + estimates <- isotonic_projection(estimates) class(estimates) <- "lmtp_survival" estimates +} + +isotonic_projection <- function(x, alpha = 0.05) { + cv <- abs(qnorm(p = alpha / 2)) + estim <- do.call("rbind", lapply(x, tidy)) + iso_fit <- isotone::gpava(1:length(x), estim$estimate) + for (i in seq_along(x)) { + x[[i]]$theta <- iso_fit$y[i] + x[[i]]$low <- x[[i]]$theta - (qnorm(0.975) * x[[i]]$standard_error) + x[[i]]$high <- x[[i]]$theta + (qnorm(0.975) * x[[i]]$standard_error) + } + x } \ No newline at end of file diff --git a/R/print.R b/R/print.R index eb7bf40..af7061f 100644 --- a/R/print.R +++ b/R/print.R @@ -26,5 +26,5 @@ print.lmtp_contrast <- function(x, ...) { print.lmtp_survival <- function(x, ...) { out <- do.call("rbind", lapply(x, tidy)) out$t <- 1:length(x) - out[, c(ncol(out), 1:ncol(out) - 1)] + as.data.frame(out[, c(ncol(out), 1:ncol(out) - 1)]) } diff --git a/R/utils.R b/R/utils.R index 92becb1..91b60fd 100644 --- a/R/utils.R +++ b/R/utils.R @@ -227,3 +227,13 @@ is_normalized <- function(x, tolerance = .Machine$double.eps^0.5) { # Check if the mean is approximately 1 within the given tolerance abs(mean(x) - 1) < tolerance } + +fix_surv_time1 <- function(x) { + x[[1]]$theta <- 1 - x[[1]]$theta + high <- x[[1]]$high + low <- x[[1]]$low + x[[1]]$high <- 1 - low + x[[1]]$low <- 1 - high + x[[1]]$eif <- 1 - x[[1]]$eif + x +} \ No newline at end of file From d30f1f5ddec5aa11e6d7358af23f5cf009ccc1aa Mon Sep 17 00:00:00 2001 From: nt-williams Date: Thu, 27 Jun 2024 16:04:21 -0700 Subject: [PATCH 08/26] lmtp_survival documentation --- NAMESPACE | 2 + R/lmtp_survival.R | 112 +++++++++++++++++++++++------- inst/examples/lmtp_survival-ex.R | 10 +++ man/lmtp_survival.Rd | 113 +++++++++++++++++++++++++++++++ 4 files changed, 212 insertions(+), 25 deletions(-) create mode 100644 inst/examples/lmtp_survival-ex.R create mode 100644 man/lmtp_survival.Rd diff --git a/NAMESPACE b/NAMESPACE index 102394e..5706430 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,6 +2,7 @@ S3method(print,lmtp) S3method(print,lmtp_contrast) +S3method(print,lmtp_survival) S3method(tidy,lmtp) export(create_node_list) export(event_locf) @@ -11,6 +12,7 @@ export(lmtp_control) export(lmtp_ipw) export(lmtp_sdr) export(lmtp_sub) +export(lmtp_survival) export(lmtp_tmle) export(static_binary_off) export(static_binary_on) diff --git a/R/lmtp_survival.R b/R/lmtp_survival.R index 7e69c57..6733f47 100644 --- a/R/lmtp_survival.R +++ b/R/lmtp_survival.R @@ -1,12 +1,74 @@ +#' LMTP Survival Curve Estimator +#' +#' Wrapper around \code{lmtp_tmle} and \code{lmtp_sdr} for survival outcomes to estimate the entire survival curve. +#' Estimates are reconstructed using isotonic regression to enforce monotonicity of the survival curve. +#' +#' @param data \[\code{data.frame}\]\cr +#' A \code{data.frame} in wide format containing all necessary variables +#' for the estimation problem. Must not be a \code{data.table}. +#' @param trt \[\code{character}\] or \[\code{list}\]\cr +#' A vector containing the column names of treatment variables ordered by time. +#' Or, a list of vectors, the same length as the number of time points of observation. +#' Vectors should contain column names for the treatment variables at each time point. The list +#' should be ordered following the time ordering of the model. +#' @param outcomes \[\code{character}\]\cr +#' A vector containing the columns names of intermediate outcome variables and the final +#' outcome variable ordered by time. Only numeric values are allowed. Variables should be coded as 0 and 1. +#' @param baseline \[\code{character}\]\cr +#' An optional vector containing the column names of baseline covariates to be +#' included for adjustment at every time point. +#' @param time_vary \[\code{list}\]\cr +#' A list the same length as the number of time points of observation with +#' the column names for new time-varying covariates introduced at each time point. The list +#' should be ordered following the time ordering of the model. +#' @param cens \[\code{character}\]\cr +#' An optional vector of column names of censoring indicators the same +#' length as the number of time points of observation. If missingness in the outcome is +#' present or if time-to-event outcome, must be provided. +#' @param shift \[\code{closure}\]\cr +#' A two argument function that specifies how treatment variables should be shifted. +#' See examples for how to specify shift functions for continuous, binary, and categorical exposures. +#' @param shifted \[\code{data.frame}\]\cr +#' An optional data frame, the same as in \code{data}, but modified according +#' to the treatment policy of interest. If specified, \code{shift} is ignored. +#' @param estimator \[\code{character(1)}\]\cr +#' The estimator to use. Either \code{"lmtp_tmle"} or \code{"lmtp_sdr"}. +#' @param k \[\code{integer(1)}\]\cr +#' An integer specifying how previous time points should be +#' used for estimation at the given time point. Default is \code{Inf}, +#' all time points. +#' @param mtp \[\code{logical(1)}\]\cr +#' Is the intervention of interest a modified treatment policy? +#' Default is \code{FALSE}. If treatment variables are continuous this should be \code{TRUE}. +#' @param id \[\code{character(1)}\]\cr +#' An optional column name containing cluster level identifiers. +#' @param learners_outcome \[\code{character}\]\cr A vector of \code{SuperLearner} algorithms for estimation +#' of the outcome regression. Default is \code{"SL.glm"}, a main effects GLM. +#' @param learners_trt \[\code{character}\]\cr A vector of \code{SuperLearner} algorithms for estimation +#' of the exposure mechanism. Default is \code{"SL.glm"}, a main effects GLM. +#' \bold{Only include candidate learners capable of binary classification}. +#' @param folds \[\code{integer(1)}\]\cr +#' The number of folds to be used for cross-fitting. +#' @param weights \[\code{numeric(nrow(data))}\]\cr +#' An optional vector containing sampling weights. +#' @param control \[\code{list()}\]\cr +#' Output of \code{lmtp_control()}. +#' +#' @return A list of class \code{lmtp_survival} containing \code{lmtp} objects for each time point. +#' +#' @example inst/examples/lmtp_survival-ex.R +#' @export lmtp_survival <- function(data, trt, outcomes, baseline = NULL, time_vary = NULL, - cens = NULL, shift = NULL, shifted = NULL, estimator = c("lmtp_tmle", "lmtp_sdr"), k = Inf, - mtp = FALSE, - id = NULL, - learners_outcome = "SL.glm", - learners_trt = "SL.glm", - folds = 10, weights = NULL, - control = lmtp_control()) { - + cens = NULL, shift = NULL, shifted = NULL, + estimator = c("lmtp_tmle", "lmtp_sdr"), + k = Inf, + mtp = FALSE, + id = NULL, + learners_outcome = "SL.glm", + learners_trt = "SL.glm", + folds = 10, weights = NULL, + control = lmtp_control()) { + checkmate::assertCharacter(outcomes, min.len = 2, null.ok = FALSE, unique = TRUE, any.missing = FALSE) estimator <- match.arg(estimator) @@ -15,22 +77,22 @@ lmtp_survival <- function(data, trt, outcomes, baseline = NULL, time_vary = NULL for (t in 1:tau) { args <- list( - data = data, - trt = trt, - outcome = outcomes[1:t], - baseline = baseline, - time_vary = time_vary, - cens = cens[1:t], - shift = shift, - shifted = shifted, - k = k, - mtp = mtp, - outcome_type = ifelse(t == 1, "binomial", "survival"), - id = id, - learners_outcome = learners_outcome, - learners_trt = learners_trt, - folds = folds, - weights = weights, + data = data, + trt = trt, + outcome = outcomes[1:t], + baseline = baseline, + time_vary = time_vary, + cens = cens[1:t], + shift = shift, + shifted = shifted, + k = k, + mtp = mtp, + outcome_type = ifelse(t == 1, "binomial", "survival"), + id = id, + learners_outcome = learners_outcome, + learners_trt = learners_trt, + folds = folds, + weights = weights, control = control ) @@ -59,4 +121,4 @@ isotonic_projection <- function(x, alpha = 0.05) { x[[i]]$high <- x[[i]]$theta + (qnorm(0.975) * x[[i]]$standard_error) } x -} \ No newline at end of file +} diff --git a/inst/examples/lmtp_survival-ex.R b/inst/examples/lmtp_survival-ex.R new file mode 100644 index 0000000..0722fed --- /dev/null +++ b/inst/examples/lmtp_survival-ex.R @@ -0,0 +1,10 @@ +# Time-to-event analysis with a binary time-invariant exposure. Interested in +# the effect of treatment being given to all observations on the cumulative +# incidence of the outcome. +A <- "trt" +Y <- paste0("Y.", 1:6) +C <- paste0("C.", 0:5) +W <- c("W1", "W2") + +lmtp_survival(sim_point_surv, A, Y, W, cens = C, folds = 1, + shift = static_binary_on, estimator = "lmtp_tmle") diff --git a/man/lmtp_survival.Rd b/man/lmtp_survival.Rd new file mode 100644 index 0000000..8cae122 --- /dev/null +++ b/man/lmtp_survival.Rd @@ -0,0 +1,113 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/lmtp_survival.R +\name{lmtp_survival} +\alias{lmtp_survival} +\title{LMTP Survival Curve Estimator} +\usage{ +lmtp_survival( + data, + trt, + outcomes, + baseline = NULL, + time_vary = NULL, + cens = NULL, + shift = NULL, + shifted = NULL, + estimator = c("lmtp_tmle", "lmtp_sdr"), + k = Inf, + mtp = FALSE, + id = NULL, + learners_outcome = "SL.glm", + learners_trt = "SL.glm", + folds = 10, + weights = NULL, + control = lmtp_control() +) +} +\arguments{ +\item{data}{[\code{data.frame}]\cr +A \code{data.frame} in wide format containing all necessary variables +for the estimation problem. Must not be a \code{data.table}.} + +\item{trt}{[\code{character}] or [\code{list}]\cr +A vector containing the column names of treatment variables ordered by time. +Or, a list of vectors, the same length as the number of time points of observation. +Vectors should contain column names for the treatment variables at each time point. The list +should be ordered following the time ordering of the model.} + +\item{outcomes}{[\code{character}]\cr +A vector containing the columns names of intermediate outcome variables and the final +outcome variable ordered by time. Only numeric values are allowed. Variables should be coded as 0 and 1.} + +\item{baseline}{[\code{character}]\cr +An optional vector containing the column names of baseline covariates to be +included for adjustment at every time point.} + +\item{time_vary}{[\code{list}]\cr +A list the same length as the number of time points of observation with +the column names for new time-varying covariates introduced at each time point. The list +should be ordered following the time ordering of the model.} + +\item{cens}{[\code{character}]\cr +An optional vector of column names of censoring indicators the same +length as the number of time points of observation. If missingness in the outcome is +present or if time-to-event outcome, must be provided.} + +\item{shift}{[\code{closure}]\cr +A two argument function that specifies how treatment variables should be shifted. +See examples for how to specify shift functions for continuous, binary, and categorical exposures.} + +\item{shifted}{[\code{data.frame}]\cr +An optional data frame, the same as in \code{data}, but modified according +to the treatment policy of interest. If specified, \code{shift} is ignored.} + +\item{estimator}{[\code{character(1)}]\cr +The estimator to use. Either \code{"lmtp_tmle"} or \code{"lmtp_sdr"}.} + +\item{k}{[\code{integer(1)}]\cr +An integer specifying how previous time points should be +used for estimation at the given time point. Default is \code{Inf}, +all time points.} + +\item{mtp}{[\code{logical(1)}]\cr +Is the intervention of interest a modified treatment policy? +Default is \code{FALSE}. If treatment variables are continuous this should be \code{TRUE}.} + +\item{id}{[\code{character(1)}]\cr +An optional column name containing cluster level identifiers.} + +\item{learners_outcome}{[\code{character}]\cr A vector of \code{SuperLearner} algorithms for estimation +of the outcome regression. Default is \code{"SL.glm"}, a main effects GLM.} + +\item{learners_trt}{[\code{character}]\cr A vector of \code{SuperLearner} algorithms for estimation +of the exposure mechanism. Default is \code{"SL.glm"}, a main effects GLM. +\bold{Only include candidate learners capable of binary classification}.} + +\item{folds}{[\code{integer(1)}]\cr +The number of folds to be used for cross-fitting.} + +\item{weights}{[\code{numeric(nrow(data))}]\cr +An optional vector containing sampling weights.} + +\item{control}{[\code{list()}]\cr +Output of \code{lmtp_control()}.} +} +\value{ +A list of class \code{lmtp_survival} containing \code{lmtp} objects for each time point. +} +\description{ +Wrapper around \code{lmtp_tmle} and \code{lmtp_sdr} for survival outcomes to estimate the entire survival curve. +Estimates are reconstructed using isotonic regression to enforce monotonicity of the survival curve. +} +\examples{ +# Time-to-event analysis with a binary time-invariant exposure. Interested in +# the effect of treatment being given to all observations on the cumulative +# incidence of the outcome. +A <- "trt" +Y <- paste0("Y.", 1:6) +C <- paste0("C.", 0:5) +W <- c("W1", "W2") + +lmtp_survival(sim_point_surv, A, Y, W, cens = C, folds = 1, + shift = static_binary_on, estimator = "lmtp_tmle") +} From c7284332dff2fada8948c66e5ca9eb4314c685b6 Mon Sep 17 00:00:00 2001 From: nt-williams Date: Thu, 27 Jun 2024 16:04:46 -0700 Subject: [PATCH 09/26] Version bump --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 14eebb3..1cc8e67 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: lmtp Title: Non-Parametric Causal Effects of Feasible Interventions Based on Modified Treatment Policies -Version: 1.4.0 +Version: 1.4.1 Authors@R: c(person(given = "Nicholas", family = "Williams", From 6434a860ae64cdbba2f699f1480fdc24bab62c7f Mon Sep 17 00:00:00 2001 From: nt-williams Date: Thu, 27 Jun 2024 16:12:01 -0700 Subject: [PATCH 10/26] Updateds NEWS.md --- NEWS.md | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/NEWS.md b/NEWS.md index 10cdb3a..3b6a0d1 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,13 @@ +# lmtp 1.4.1 + +### New Features + +- Added `lmtp_survival()` function for estimating the entire survival curve. Enforces monotonicity using isotonic regression (see issue \#140. + +### Bug Fixes + +### General + # lmtp 1.4.0 ### New Features From 609cb610f254635988ed2c859b75cb806e7a18b9 Mon Sep 17 00:00:00 2001 From: nt-williams Date: Fri, 28 Jun 2024 09:51:04 -0700 Subject: [PATCH 11/26] Tidy method for lmtp_survival --- NAMESPACE | 1 + R/lmtp_survival.R | 2 +- R/print.R | 4 +--- R/tidy.R | 14 ++++++++++++++ inst/examples/lmtp_survival-ex.R | 8 ++++++-- man/lmtp_survival.Rd | 8 ++++++-- man/tidy.lmtp_survival.Rd | 32 ++++++++++++++++++++++++++++++++ 7 files changed, 61 insertions(+), 8 deletions(-) create mode 100644 man/tidy.lmtp_survival.Rd diff --git a/NAMESPACE b/NAMESPACE index 5706430..bbeee75 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,6 +4,7 @@ S3method(print,lmtp) S3method(print,lmtp_contrast) S3method(print,lmtp_survival) S3method(tidy,lmtp) +S3method(tidy,lmtp_survival) export(create_node_list) export(event_locf) export(ipsi) diff --git a/R/lmtp_survival.R b/R/lmtp_survival.R index 6733f47..b0c4685 100644 --- a/R/lmtp_survival.R +++ b/R/lmtp_survival.R @@ -113,7 +113,7 @@ lmtp_survival <- function(data, trt, outcomes, baseline = NULL, time_vary = NULL isotonic_projection <- function(x, alpha = 0.05) { cv <- abs(qnorm(p = alpha / 2)) - estim <- do.call("rbind", lapply(x, tidy)) + estim <- tidy.lmtp_survival(x) iso_fit <- isotone::gpava(1:length(x), estim$estimate) for (i in seq_along(x)) { x[[i]]$theta <- iso_fit$y[i] diff --git a/R/print.R b/R/print.R index af7061f..432815b 100644 --- a/R/print.R +++ b/R/print.R @@ -24,7 +24,5 @@ print.lmtp_contrast <- function(x, ...) { #' @export print.lmtp_survival <- function(x, ...) { - out <- do.call("rbind", lapply(x, tidy)) - out$t <- 1:length(x) - as.data.frame(out[, c(ncol(out), 1:ncol(out) - 1)]) + as.data.frame(tidy.lmtp_survival(x)) } diff --git a/R/tidy.R b/R/tidy.R index 8ff520b..12f0e96 100644 --- a/R/tidy.R +++ b/R/tidy.R @@ -28,3 +28,17 @@ tidy.lmtp <- function(x, ...) { class(out) <- c("tbl_df", "tbl", "data.frame") out } + +#' Tidy a(n) lmtp_survival object +#' +#' @param x A `lmtp_survival` object produced by a call to [lmtp::lmtp_survival()]. +#' @param ... Unused, included for generic consistency only. +#' +#' @example inst/examples/lmtp_survival-ex.R +#' +#' @export +tidy.lmtp_survival <- function(x, ...) { + out <- do.call("rbind", lapply(x, tidy)) + out$t <- 1:length(x) + out[, c(ncol(out), 1:ncol(out) - 1)] +} diff --git a/inst/examples/lmtp_survival-ex.R b/inst/examples/lmtp_survival-ex.R index 0722fed..8a01cd7 100644 --- a/inst/examples/lmtp_survival-ex.R +++ b/inst/examples/lmtp_survival-ex.R @@ -1,3 +1,4 @@ +\donttest{ # Time-to-event analysis with a binary time-invariant exposure. Interested in # the effect of treatment being given to all observations on the cumulative # incidence of the outcome. @@ -6,5 +7,8 @@ Y <- paste0("Y.", 1:6) C <- paste0("C.", 0:5) W <- c("W1", "W2") -lmtp_survival(sim_point_surv, A, Y, W, cens = C, folds = 1, - shift = static_binary_on, estimator = "lmtp_tmle") +curve <- lmtp_survival(sim_point_surv, A, Y, W, cens = C, folds = 1, + shift = static_binary_on, estimator = "lmtp_tmle") + +tidy(curve) +} diff --git a/man/lmtp_survival.Rd b/man/lmtp_survival.Rd index 8cae122..6e634d0 100644 --- a/man/lmtp_survival.Rd +++ b/man/lmtp_survival.Rd @@ -100,6 +100,7 @@ Wrapper around \code{lmtp_tmle} and \code{lmtp_sdr} for survival outcomes to est Estimates are reconstructed using isotonic regression to enforce monotonicity of the survival curve. } \examples{ +\donttest{ # Time-to-event analysis with a binary time-invariant exposure. Interested in # the effect of treatment being given to all observations on the cumulative # incidence of the outcome. @@ -108,6 +109,9 @@ Y <- paste0("Y.", 1:6) C <- paste0("C.", 0:5) W <- c("W1", "W2") -lmtp_survival(sim_point_surv, A, Y, W, cens = C, folds = 1, - shift = static_binary_on, estimator = "lmtp_tmle") +curve <- lmtp_survival(sim_point_surv, A, Y, W, cens = C, folds = 1, + shift = static_binary_on, estimator = "lmtp_tmle") + +tidy(curve) +} } diff --git a/man/tidy.lmtp_survival.Rd b/man/tidy.lmtp_survival.Rd new file mode 100644 index 0000000..a356f0c --- /dev/null +++ b/man/tidy.lmtp_survival.Rd @@ -0,0 +1,32 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/tidy.R +\name{tidy.lmtp_survival} +\alias{tidy.lmtp_survival} +\title{Tidy a(n) lmtp_survival object} +\usage{ +\method{tidy}{lmtp_survival}(x, ...) +} +\arguments{ +\item{x}{A \code{lmtp_survival} object produced by a call to \code{\link[=lmtp_survival]{lmtp_survival()}}.} + +\item{...}{Unused, included for generic consistency only.} +} +\description{ +Tidy a(n) lmtp_survival object +} +\examples{ +\donttest{ +# Time-to-event analysis with a binary time-invariant exposure. Interested in +# the effect of treatment being given to all observations on the cumulative +# incidence of the outcome. +A <- "trt" +Y <- paste0("Y.", 1:6) +C <- paste0("C.", 0:5) +W <- c("W1", "W2") + +curve <- lmtp_survival(sim_point_surv, A, Y, W, cens = C, folds = 1, + shift = static_binary_on, estimator = "lmtp_tmle") + +tidy(curve) +} +} From 0c31ee385603651b4dd318246407955933c5e568 Mon Sep 17 00:00:00 2001 From: nt-williams Date: Fri, 28 Jun 2024 10:01:52 -0700 Subject: [PATCH 12/26] Progress ticker for lmtp_survival, updated documentation --- R/lmtp_survival.R | 5 +++++ R/print.R | 2 +- man/lmtp_survival.Rd | 1 + 3 files changed, 7 insertions(+), 1 deletion(-) diff --git a/R/lmtp_survival.R b/R/lmtp_survival.R index b0c4685..313d474 100644 --- a/R/lmtp_survival.R +++ b/R/lmtp_survival.R @@ -2,6 +2,7 @@ #' #' Wrapper around \code{lmtp_tmle} and \code{lmtp_sdr} for survival outcomes to estimate the entire survival curve. #' Estimates are reconstructed using isotonic regression to enforce monotonicity of the survival curve. +#' \bold{Confidence intervals correspond to marginal confidence intervals for the survival curve, not simultaneous intervals.} #' #' @param data \[\code{data.frame}\]\cr #' A \code{data.frame} in wide format containing all necessary variables @@ -75,6 +76,8 @@ lmtp_survival <- function(data, trt, outcomes, baseline = NULL, time_vary = NULL tau <- length(outcomes) estimates <- vector("list", tau) + t <- 1 + cli::cli_progress_step("Working on time {t}/{tau}...") for (t in 1:tau) { args <- list( data = data, @@ -101,8 +104,10 @@ lmtp_survival <- function(data, trt, outcomes, baseline = NULL, time_vary = NULL else do.call(lmtp_sdr, args) }, seed = TRUE) + cli::cli_progress_update() } + cli::cli_progress_done() estimates <- future::value(estimates) estimates <- fix_surv_time1(estimates) estimates <- isotonic_projection(estimates) diff --git a/R/print.R b/R/print.R index 432815b..7f931f1 100644 --- a/R/print.R +++ b/R/print.R @@ -24,5 +24,5 @@ print.lmtp_contrast <- function(x, ...) { #' @export print.lmtp_survival <- function(x, ...) { - as.data.frame(tidy.lmtp_survival(x)) + print(as.data.frame(tidy.lmtp_survival(x))) } diff --git a/man/lmtp_survival.Rd b/man/lmtp_survival.Rd index 6e634d0..9f693bd 100644 --- a/man/lmtp_survival.Rd +++ b/man/lmtp_survival.Rd @@ -98,6 +98,7 @@ A list of class \code{lmtp_survival} containing \code{lmtp} objects for each tim \description{ Wrapper around \code{lmtp_tmle} and \code{lmtp_sdr} for survival outcomes to estimate the entire survival curve. Estimates are reconstructed using isotonic regression to enforce monotonicity of the survival curve. +\bold{Confidence intervals correspond to marginal confidence intervals for the survival curve, not simultaneous intervals.} } \examples{ \donttest{ From 102002bfa1beb9290ba9e9a34c2a24ffadd16014 Mon Sep 17 00:00:00 2001 From: nt-williams Date: Sat, 29 Jun 2024 16:31:59 -0700 Subject: [PATCH 13/26] Refactoring for loop in lmtp_survival --- R/lmtp_survival.R | 53 ++++++++++++++++++++++++----------------------- 1 file changed, 27 insertions(+), 26 deletions(-) diff --git a/R/lmtp_survival.R b/R/lmtp_survival.R index 313d474..42c985d 100644 --- a/R/lmtp_survival.R +++ b/R/lmtp_survival.R @@ -67,7 +67,8 @@ lmtp_survival <- function(data, trt, outcomes, baseline = NULL, time_vary = NULL id = NULL, learners_outcome = "SL.glm", learners_trt = "SL.glm", - folds = 10, weights = NULL, + folds = 10, + weights = NULL, control = lmtp_control()) { checkmate::assertCharacter(outcomes, min.len = 2, null.ok = FALSE, unique = TRUE, any.missing = FALSE) @@ -76,34 +77,34 @@ lmtp_survival <- function(data, trt, outcomes, baseline = NULL, time_vary = NULL tau <- length(outcomes) estimates <- vector("list", tau) + args <- list( + data = data, + trt = trt, + baseline = baseline, + time_vary = time_vary, + shift = shift, + shifted = shifted, + k = k, + mtp = mtp, + id = id, + learners_outcome = learners_outcome, + learners_trt = learners_trt, + folds = folds, + weights = weights, + control = control + ) + + if (estimator == "lmtp_tmle") expr <- expression(do.call(lmtp_tmle, args)) + else expr <- expression(do.call(lmtp_sdr, args)) + t <- 1 cli::cli_progress_step("Working on time {t}/{tau}...") for (t in 1:tau) { - args <- list( - data = data, - trt = trt, - outcome = outcomes[1:t], - baseline = baseline, - time_vary = time_vary, - cens = cens[1:t], - shift = shift, - shifted = shifted, - k = k, - mtp = mtp, - outcome_type = ifelse(t == 1, "binomial", "survival"), - id = id, - learners_outcome = learners_outcome, - learners_trt = learners_trt, - folds = folds, - weights = weights, - control = control - ) + args$outcome <- outcomes[1:t] + args$cens <- cens[1:t] + args$outcome_type <- ifelse(t == 1, "binomial", "survival") - estimates[[t]] <- future::future({ - if (estimator == "lmtp_tmle") do.call(lmtp_tmle, args) - else do.call(lmtp_sdr, args) - }, - seed = TRUE) + estimates[[t]] <- future::future(eval(expr), seed = TRUE) cli::cli_progress_update() } @@ -122,7 +123,7 @@ isotonic_projection <- function(x, alpha = 0.05) { iso_fit <- isotone::gpava(1:length(x), estim$estimate) for (i in seq_along(x)) { x[[i]]$theta <- iso_fit$y[i] - x[[i]]$low <- x[[i]]$theta - (qnorm(0.975) * x[[i]]$standard_error) + x[[i]]$low <- x[[i]]$theta - (qnorm(0.975) * x[[i]]$standard_error) x[[i]]$high <- x[[i]]$theta + (qnorm(0.975) * x[[i]]$standard_error) } x From b3c826e07b71fda971679613c42e94819d17416b Mon Sep 17 00:00:00 2001 From: nt-williams Date: Sat, 29 Jun 2024 16:44:58 -0700 Subject: [PATCH 14/26] Capturing boot call in lmtp_survival --- R/lmtp_survival.R | 11 +++++++++-- man/lmtp_survival.Rd | 6 ++++++ 2 files changed, 15 insertions(+), 2 deletions(-) diff --git a/R/lmtp_survival.R b/R/lmtp_survival.R index 42c985d..ffaac84 100644 --- a/R/lmtp_survival.R +++ b/R/lmtp_survival.R @@ -41,6 +41,10 @@ #' @param mtp \[\code{logical(1)}\]\cr #' Is the intervention of interest a modified treatment policy? #' Default is \code{FALSE}. If treatment variables are continuous this should be \code{TRUE}. +#' @param boot \[\code{logical(1)}\]\cr +#' Compute standard errors using the bootstrap? Default is \code{TRUE}. If \code{FALSE}, standard +#' errors will be calculated using the empirical variance of the efficient influence function. +#' Ignored if \code{estimator = "lmtp_sdr"}. #' @param id \[\code{character(1)}\]\cr #' An optional column name containing cluster level identifiers. #' @param learners_outcome \[\code{character}\]\cr A vector of \code{SuperLearner} algorithms for estimation @@ -64,6 +68,7 @@ lmtp_survival <- function(data, trt, outcomes, baseline = NULL, time_vary = NULL estimator = c("lmtp_tmle", "lmtp_sdr"), k = Inf, mtp = FALSE, + boot = TRUE, id = NULL, learners_outcome = "SL.glm", learners_trt = "SL.glm", @@ -94,8 +99,10 @@ lmtp_survival <- function(data, trt, outcomes, baseline = NULL, time_vary = NULL control = control ) - if (estimator == "lmtp_tmle") expr <- expression(do.call(lmtp_tmle, args)) - else expr <- expression(do.call(lmtp_sdr, args)) + if (estimator == "lmtp_tmle") { + args$boot <- boot + expr <- expression(do.call(lmtp_tmle, args)) + } else expr <- expression(do.call(lmtp_sdr, args)) t <- 1 cli::cli_progress_step("Working on time {t}/{tau}...") diff --git a/man/lmtp_survival.Rd b/man/lmtp_survival.Rd index 9f693bd..906c8c1 100644 --- a/man/lmtp_survival.Rd +++ b/man/lmtp_survival.Rd @@ -16,6 +16,7 @@ lmtp_survival( estimator = c("lmtp_tmle", "lmtp_sdr"), k = Inf, mtp = FALSE, + boot = TRUE, id = NULL, learners_outcome = "SL.glm", learners_trt = "SL.glm", @@ -73,6 +74,11 @@ all time points.} Is the intervention of interest a modified treatment policy? Default is \code{FALSE}. If treatment variables are continuous this should be \code{TRUE}.} +\item{boot}{[\code{logical(1)}]\cr +Compute standard errors using the bootstrap? Default is \code{TRUE}. If \code{FALSE}, standard +errors will be calculated using the empirical variance of the efficient influence function. +Ignored if \code{estimator = "lmtp_sdr"}.} + \item{id}{[\code{character(1)}]\cr An optional column name containing cluster level identifiers.} From dbe5f7f75c9f7b51efb1d3a18da8bd48b8a1bc7f Mon Sep 17 00:00:00 2001 From: nt-williams Date: Sun, 30 Jun 2024 09:41:02 -0600 Subject: [PATCH 15/26] bug fixes for lmtp_survival with time varying treatment --- R/lmtp_survival.R | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/R/lmtp_survival.R b/R/lmtp_survival.R index 42c985d..6ef1808 100644 --- a/R/lmtp_survival.R +++ b/R/lmtp_survival.R @@ -79,9 +79,7 @@ lmtp_survival <- function(data, trt, outcomes, baseline = NULL, time_vary = NULL args <- list( data = data, - trt = trt, baseline = baseline, - time_vary = time_vary, shift = shift, shifted = shifted, k = k, @@ -94,12 +92,17 @@ lmtp_survival <- function(data, trt, outcomes, baseline = NULL, time_vary = NULL control = control ) + if (length(trt) == 1) args$trt <- trt + if (length(time_vary) == 1) args$time_vary <- time_vary + if (estimator == "lmtp_tmle") expr <- expression(do.call(lmtp_tmle, args)) else expr <- expression(do.call(lmtp_sdr, args)) t <- 1 cli::cli_progress_step("Working on time {t}/{tau}...") for (t in 1:tau) { + if (length(trt) > 1) args$trt <- trt[1:t] + if (length(args$time_vary) > 1) args$time_vary <- time_vary[1:t] args$outcome <- outcomes[1:t] args$cens <- cens[1:t] args$outcome_type <- ifelse(t == 1, "binomial", "survival") From 3fb77331da5f7877f02f2fedcdacbe3976a1a63c Mon Sep 17 00:00:00 2001 From: Nicholas Williams Date: Wed, 11 Sep 2024 14:44:44 -0700 Subject: [PATCH 16/26] Removing schoolmath dependency --- DESCRIPTION | 3 +-- R/utils.R | 7 ++++++- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 1cc8e67..09a318b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -24,7 +24,7 @@ License: AGPL-3 Encoding: UTF-8 LazyData: true Roxygen: list(markdown = TRUE) -RoxygenNote: 7.3.1 +RoxygenNote: 7.3.2 Imports: stats, nnls, @@ -37,7 +37,6 @@ Imports: data.table (>= 1.13.0), checkmate (>= 2.1.0), SuperLearner, - schoolmath, isotone URL: https://beyondtheate.com/, https://github.com/nt-williams/lmtp BugReports: https://github.com/nt-williams/lmtp/issues diff --git a/R/utils.R b/R/utils.R index 91b60fd..a2889b8 100644 --- a/R/utils.R +++ b/R/utils.R @@ -236,4 +236,9 @@ fix_surv_time1 <- function(x) { x[[1]]$low <- 1 - high x[[1]]$eif <- 1 - x[[1]]$eif x -} \ No newline at end of file +} + +is_decimal <- function(x) { + test <- floor(x) + !(x == test) +} From 257bfd0d89d9dbec13a22d203f1788ea79a0336f Mon Sep 17 00:00:00 2001 From: Nicholas Williams Date: Wed, 11 Sep 2024 14:47:21 -0700 Subject: [PATCH 17/26] Removing schoolmath dependency --- R/checks.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/checks.R b/R/checks.R index 10e7649..f6ae1ce 100644 --- a/R/checks.R +++ b/R/checks.R @@ -177,7 +177,7 @@ check_trt_type <- function(data, trt, mtp) { for (i in seq_along(trt)) { a <- data[[trt[i]]] if (is.character(a) | is.factor(a)) next - is_decimal[i] <- any(schoolmath::is.decimal(a[!is.na(a)])) + is_decimal[i] <- any(is_decimal(a[!is.na(a)])) } if (any(is_decimal) & isFALSE(mtp)) { cli::cli_warn("Detected decimalish `trt` values and {.code mtp = FALSE}. Consider setting {.code mtp = TRUE} if getting errors.") From 818efea8e3554b680fda1e01cdd67c206164d05c Mon Sep 17 00:00:00 2001 From: Nicholas Williams Date: Wed, 11 Sep 2024 14:47:49 -0700 Subject: [PATCH 18/26] version bump --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 09a318b..1318698 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: lmtp Title: Non-Parametric Causal Effects of Feasible Interventions Based on Modified Treatment Policies -Version: 1.4.1 +Version: 1.4.2 Authors@R: c(person(given = "Nicholas", family = "Williams", From db81fcceff8d7d2191f0f51f6e2f4cca0ce9752c Mon Sep 17 00:00:00 2001 From: Nicholas Williams Date: Wed, 11 Sep 2024 14:48:54 -0700 Subject: [PATCH 19/26] Updated NEWS --- NEWS.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/NEWS.md b/NEWS.md index 3b6a0d1..6acbb75 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,9 @@ +# lmtp 1.4.2 + +### General + +- Removed dependency on `schoolmath` which used a very slow function for testing if a vector was "decimalish". + # lmtp 1.4.1 ### New Features From 98fd3833c40a19323ce363d1490d65e7cabbad78 Mon Sep 17 00:00:00 2001 From: Nicholas Williams Date: Tue, 17 Sep 2024 10:35:47 -0700 Subject: [PATCH 20/26] Changing default for boot to false --- R/estimators.R | 4 ++-- R/lmtp_survival.R | 4 ++-- man/lmtp_survival.Rd | 4 ++-- man/lmtp_tmle.Rd | 4 ++-- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/R/estimators.R b/R/estimators.R index 2d13277..5d8c8d3 100644 --- a/R/estimators.R +++ b/R/estimators.R @@ -42,7 +42,7 @@ #' Is the intervention of interest a modified treatment policy? #' Default is \code{FALSE}. If treatment variables are continuous this should be \code{TRUE}. #' @param boot \[\code{logical(1)}\]\cr -#' Compute standard errors using the bootstrap? Default is \code{TRUE}. If \code{FALSE}, standard +#' Compute standard errors using the bootstrap? Default is \code{FALSE}. If \code{FALSE}, standard #' errors will be calculated using the empirical variance of the efficient influence function. #' @param outcome_type \[\code{character(1)}\]\cr #' Outcome variable type (i.e., continuous, binomial, survival). @@ -96,7 +96,7 @@ #' @export lmtp_tmle <- function(data, trt, outcome, baseline = NULL, time_vary = NULL, cens = NULL, shift = NULL, shifted = NULL, k = Inf, - mtp = FALSE, boot = TRUE, + mtp = FALSE, boot = FALSE, outcome_type = c("binomial", "continuous", "survival"), id = NULL, bounds = NULL, learners_outcome = "SL.glm", diff --git a/R/lmtp_survival.R b/R/lmtp_survival.R index ac20e83..5f12759 100644 --- a/R/lmtp_survival.R +++ b/R/lmtp_survival.R @@ -42,7 +42,7 @@ #' Is the intervention of interest a modified treatment policy? #' Default is \code{FALSE}. If treatment variables are continuous this should be \code{TRUE}. #' @param boot \[\code{logical(1)}\]\cr -#' Compute standard errors using the bootstrap? Default is \code{TRUE}. If \code{FALSE}, standard +#' Compute standard errors using the bootstrap? Default is \code{FALSE}. If \code{FALSE}, standard #' errors will be calculated using the empirical variance of the efficient influence function. #' Ignored if \code{estimator = "lmtp_sdr"}. #' @param id \[\code{character(1)}\]\cr @@ -68,7 +68,7 @@ lmtp_survival <- function(data, trt, outcomes, baseline = NULL, time_vary = NULL estimator = c("lmtp_tmle", "lmtp_sdr"), k = Inf, mtp = FALSE, - boot = TRUE, + boot = FALSE, id = NULL, learners_outcome = "SL.glm", learners_trt = "SL.glm", diff --git a/man/lmtp_survival.Rd b/man/lmtp_survival.Rd index 906c8c1..5baf4e0 100644 --- a/man/lmtp_survival.Rd +++ b/man/lmtp_survival.Rd @@ -16,7 +16,7 @@ lmtp_survival( estimator = c("lmtp_tmle", "lmtp_sdr"), k = Inf, mtp = FALSE, - boot = TRUE, + boot = FALSE, id = NULL, learners_outcome = "SL.glm", learners_trt = "SL.glm", @@ -75,7 +75,7 @@ Is the intervention of interest a modified treatment policy? Default is \code{FALSE}. If treatment variables are continuous this should be \code{TRUE}.} \item{boot}{[\code{logical(1)}]\cr -Compute standard errors using the bootstrap? Default is \code{TRUE}. If \code{FALSE}, standard +Compute standard errors using the bootstrap? Default is \code{FALSE}. If \code{FALSE}, standard errors will be calculated using the empirical variance of the efficient influence function. Ignored if \code{estimator = "lmtp_sdr"}.} diff --git a/man/lmtp_tmle.Rd b/man/lmtp_tmle.Rd index 1d19454..e8dcc50 100644 --- a/man/lmtp_tmle.Rd +++ b/man/lmtp_tmle.Rd @@ -15,7 +15,7 @@ lmtp_tmle( shifted = NULL, k = Inf, mtp = FALSE, - boot = TRUE, + boot = FALSE, outcome_type = c("binomial", "continuous", "survival"), id = NULL, bounds = NULL, @@ -75,7 +75,7 @@ Is the intervention of interest a modified treatment policy? Default is \code{FALSE}. If treatment variables are continuous this should be \code{TRUE}.} \item{boot}{[\code{logical(1)}]\cr -Compute standard errors using the bootstrap? Default is \code{TRUE}. If \code{FALSE}, standard +Compute standard errors using the bootstrap? Default is \code{FALSE}. If \code{FALSE}, standard errors will be calculated using the empirical variance of the efficient influence function.} \item{outcome_type}{[\code{character(1)}]\cr From 387070a59126c499bc0efae4284c599de6805947 Mon Sep 17 00:00:00 2001 From: Nicholas Williams Date: Thu, 19 Sep 2024 09:18:47 -0700 Subject: [PATCH 21/26] Accounting for ignored boot call in lmtp_survival --- R/lmtp_survival.R | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/R/lmtp_survival.R b/R/lmtp_survival.R index 5f12759..bca37b2 100644 --- a/R/lmtp_survival.R +++ b/R/lmtp_survival.R @@ -100,8 +100,12 @@ lmtp_survival <- function(data, trt, outcomes, baseline = NULL, time_vary = NULL if (length(trt) == 1) args$trt <- trt if (length(time_vary) == 1) args$time_vary <- time_vary - if (estimator == "lmtp_tmle") expr <- expression(do.call(lmtp_tmle, args)) - else expr <- expression(do.call(lmtp_sdr, args)) + if (estimator == "lmtp_tmle") { + args$boot <- boot + expr <- expression(do.call(lmtp_tmle, args)) + } else { + expr <- expression(do.call(lmtp_sdr, args)) + } t <- 1 cli::cli_progress_step("Working on time {t}/{tau}...") From 15241d85811d220983b40b35f9a34c1cbf4c11d2 Mon Sep 17 00:00:00 2001 From: Nicholas Williams Date: Thu, 19 Sep 2024 09:19:41 -0700 Subject: [PATCH 22/26] Updated NEWS --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 41795f0..bb90160 100644 --- a/NEWS.md +++ b/NEWS.md @@ -3,13 +3,13 @@ ### General - Removed dependency on `schoolmath` which used a very slow function for testing if a vector was "decimalish". +- Bootstrap for TMLE with the `boot` argument using a modified TMLE algorithm (https://arxiv.org/abs/1810.03030). # lmtp 1.4.1 ### New Features - Added `lmtp_survival()` function for estimating the entire survival curve. Enforces monotonicity using isotonic regression (see issue \#140). -- Bootstrap for TMLE with the `boot` argument using a modified TMLE algorithm (https://arxiv.org/abs/1810.03030). ### Bug Fixes From 7be1f431789a1d4ae05c676849da4a31cc7d953a Mon Sep 17 00:00:00 2001 From: Nicholas Williams Date: Thu, 26 Sep 2024 11:13:27 -0700 Subject: [PATCH 23/26] Bug fix in isotonic regression --- DESCRIPTION | 2 +- NEWS.md | 6 ++++++ R/lmtp_survival.R | 4 ++-- 3 files changed, 9 insertions(+), 3 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 1318698..7fb8908 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: lmtp Title: Non-Parametric Causal Effects of Feasible Interventions Based on Modified Treatment Policies -Version: 1.4.2 +Version: 1.4.3 Authors@R: c(person(given = "Nicholas", family = "Williams", diff --git a/NEWS.md b/NEWS.md index bb90160..56f4315 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,9 @@ +# lmtp 1.4.3 + +### Bug Fixes + +- Using fitted values from isotonic regression in `lmtp_survival()` instead of the original values. + # lmtp 1.4.2 ### General diff --git a/R/lmtp_survival.R b/R/lmtp_survival.R index bca37b2..afb484c 100644 --- a/R/lmtp_survival.R +++ b/R/lmtp_survival.R @@ -132,9 +132,9 @@ lmtp_survival <- function(data, trt, outcomes, baseline = NULL, time_vary = NULL isotonic_projection <- function(x, alpha = 0.05) { cv <- abs(qnorm(p = alpha / 2)) estim <- tidy.lmtp_survival(x) - iso_fit <- isotone::gpava(1:length(x), estim$estimate) + iso_fit <- isotone::gpava(1:length(x), 1 - estim$estimate) for (i in seq_along(x)) { - x[[i]]$theta <- iso_fit$y[i] + x[[i]]$theta <- (1 - iso_fit$x[i]) x[[i]]$low <- x[[i]]$theta - (qnorm(0.975) * x[[i]]$standard_error) x[[i]]$high <- x[[i]]$theta + (qnorm(0.975) * x[[i]]$standard_error) } From 74db0f08eafbbfbb4a63a02a11046927879edaa0 Mon Sep 17 00:00:00 2001 From: Nicholas Williams Date: Thu, 26 Sep 2024 11:14:39 -0700 Subject: [PATCH 24/26] updates NEWS --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 56f4315..27e4d9d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,7 +2,7 @@ ### Bug Fixes -- Using fitted values from isotonic regression in `lmtp_survival()` instead of the original values. +- Using fitted values from isotonic regression in `lmtp_survival()` instead of the original values (see issue \#149). # lmtp 1.4.2 From 65e16c004a328269c42f44de4c9cefaa1059ac97 Mon Sep 17 00:00:00 2001 From: Nicholas Williams Date: Thu, 26 Sep 2024 11:16:52 -0700 Subject: [PATCH 25/26] Fixing version discrepancies --- DESCRIPTION | 2 +- NEWS.md | 18 +++++------------- 2 files changed, 6 insertions(+), 14 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 7fb8908..09a318b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: lmtp Title: Non-Parametric Causal Effects of Feasible Interventions Based on Modified Treatment Policies -Version: 1.4.3 +Version: 1.4.1 Authors@R: c(person(given = "Nicholas", family = "Williams", diff --git a/NEWS.md b/NEWS.md index 27e4d9d..9ea0cd3 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,26 +1,18 @@ -# lmtp 1.4.3 - -### Bug Fixes - -- Using fitted values from isotonic regression in `lmtp_survival()` instead of the original values (see issue \#149). - -# lmtp 1.4.2 - -### General - -- Removed dependency on `schoolmath` which used a very slow function for testing if a vector was "decimalish". -- Bootstrap for TMLE with the `boot` argument using a modified TMLE algorithm (https://arxiv.org/abs/1810.03030). - # lmtp 1.4.1 ### New Features - Added `lmtp_survival()` function for estimating the entire survival curve. Enforces monotonicity using isotonic regression (see issue \#140). +- Bootstrap for TMLE with the `boot` argument using a modified TMLE algorithm (https://arxiv.org/abs/1810.03030). ### Bug Fixes +- Using fitted values from isotonic regression in `lmtp_survival()` instead of the original values (see issue \#149). + ### General +- Removed dependency on `schoolmath` which used a very slow function for testing if a vector was "decimalish". + # lmtp 1.4.0 ### New Features From 9bec66196739e619d27adfb7bd3d106b454a7dfb Mon Sep 17 00:00:00 2001 From: Nicholas Williams Date: Mon, 4 Nov 2024 09:38:00 -0800 Subject: [PATCH 26/26] Fix for issue 151 --- NEWS.md | 1 + R/tmle2.R | 4 ++++ 2 files changed, 5 insertions(+) diff --git a/NEWS.md b/NEWS.md index 9ea0cd3..92942cd 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,6 +8,7 @@ ### Bug Fixes - Using fitted values from isotonic regression in `lmtp_survival()` instead of the original values (see issue \#149). +- Bootstrap TMLE uses cumulative density ratios (see issue \#151). ### General diff --git a/R/tmle2.R b/R/tmle2.R index d429489..c4c81e5 100644 --- a/R/tmle2.R +++ b/R/tmle2.R @@ -1,4 +1,8 @@ cf_tmle2 <- function(task, ratios, m_init, control) { + ratios <- matrix(t(apply(ratios, 1, cumprod)), + nrow = nrow(ratios), + ncol = ncol(ratios)) + psi <- estimate_tmle2(task$natural, task$cens, task$risk,