library(dplyr)
library(sl3)
library(Rsolnp)
library(glmnet)
library(xgboost)
library(ranger)
library(earth)
library(unihtee)Background
Treatment effect modifiers (TEMs) are responsible for the disparate effects of a treatment on a population. They are confounders — meaning they influence both the likelihood of receiving treatment and the outcome — that modify the effect of treatment on the outcome. In precision medicine, these effect modifiers delineate patient subgroups which experience differing benefit from a given medical intervention. In public health, they are used to determine the effect of policy decisions on sub-populations. Identifying TEMs, if any exist, is required for a comprehensive understanding of an intervention’s effect on a population.
Traditional parametric modeling techniques, like generalized linear models (GLMs), define TEMs as the confounders with non-zero confounder-treatment interaction terms. In settings characterized by time-to-event outcomes, the Cox proportional hazards model can similarly be used. Inference about TEMs is therefore possible under stringent conditions about the data-generating process. When these (semi)parametric methods’ assumptions, like linearity of the confounder-outcome relationship in a linear model or proportional hazards in a Cox model, are violated, however, their results are invalid.
More flexible approaches instead focus on estimating the conditional average treatment effect (CATE) using interpretable modeling techniques, like the LASSO (Tibshirani 1996; Tian et al. 2014; S. Chen et al. 2017; Q. Zhao, Small, and Ertefaie 2018; Semenova and Chernozhukov 2020; Bahamyirou et al. 2022) or Random Forests (Breiman 2001; Wager and Athey 2018; Cui et al. 2022). Again, however, these methods require restrictive assumptions about the data-generating process, like sparsity of treatment effect modification (P. Zhao and Yu 2006), approximately uncorrelated confounders (P. Zhao and Yu 2006; Hastie, Tibshirani, and Friedman 2009) and sample sizes that are much larger than the number of confounders (P. Zhao and Yu 2006), to reliably recover TEMs.
Instead of identifying TEMs indirectly by modeling the confounder-treatment-outcome relationship, recent work has developed frameworks tailored to the task (Williamson et al. 2022; Boileau et al. 2022, 2023; Hines, Diaz-Ordaz, and Vansteelandt 2022). Relying on TEM variable importance parameters (TEM-VIP), it is possible to assess confounders’ influence on the treatment effect in nonparametric, algorithm-agnostic fashion that largely avoids the pitfalls of parametric and CATE-based methodologies. These TEM-VIPs are defined within nonparametric statistical models that may be augmented with causal interpretations, permitting formal inference about TEMs.
The unihtee R package implements the recent proposals of
Boileau et al. (2022) and Boileau et al. (2023). This TEM-VIP framework
relies on parameters assess the marginal effect of each confounder on
treatment effect heterogeneity, however treatment effect is defined.
Nonparametric estimators are provided in these works, and their
asymptotic properties established. In particular, these estimators are
shown to be asymptotically linear under minimal assumptions about the
data-generating process, and to recover treatment effect modifiers more
readily than competing methods. These TEM-VIPs and estimators are
introduced alongside worked examples in the following section.
unihtee in Action
The unihtee R package performs inference about TEM-VIPs
in data-generating processes with a binary exposure variable and either
continuous, binary, or right-censored time-to-event outcomes. No
restriction is placed on the confounding variables, though ordinal
confounders should be formatted as factors using ordered(),
and categorical confounders should be one-hot encoded.
Working with Continuous Outcomes
Let there be \(n\) independent and identically distributed (i.i.d.) random vectors \({X_i}_{i=1}^n\) such that \(X_i = (W_i, A_i, Y^{(0)}, Y_i^{(1)}) \sim P_{X,0} \in \mathcal{M}_X\). Dropping the index where possible for notational convenience for the remainder of the tutorial, \(W\) is defined as the set of \(p\) confounders, \(A\) as the binary treatment indicator, and \(Y^{(0)}, Y_i^{(1)}\) the continuous potential outcomes (Rubin 1974) produced under control and treatment conditions, respectively. Further, \(\mathcal{M}_X\) is the full-data nonparametric model of all possible data-generating processes, and contains the true but unknown full-data data-generating process \(P_{X,0}\).
Now, \(P_{X,0}\) is generally unknown because its random vectors are unobservable: only one potential outcome is ever measurable in practice. Still, \(P_{X,0}\) allows to us to readily define causal parameters on which inference may subsequently be performed. Examples are given in the subsections below.
We instead have access to \(n\) i.i.d. random observations \(O=(W,A,Y) \sim P_0 \in \mathcal{M}\), where \(W\) and \(A\) are defined as in the full-data, and where \(Y = AY^{(1)} + (1-A)Y{(0)}\). \(\mathcal{M}\) is the observed-data statistical model containing all possible observed-data data-generating processes, including the true data-generating process \(P_{0}\). Note too that the elements of \(\mathcal{M}\) are fully specified by their full-data counterparts in \(\mathcal{M}_{X}\).
Having established the requisite statistical objects, let’s define a data-generating process:
cont_outcome_dgp <- function(n_obs) {
# confounders
w_1 <- rnorm(n = n_obs)
w_2 <- rnorm(n = n_obs)
w_3 <- rnorm(n = n_obs)
w_4 <- rnorm(n = n_obs)
w_5 <- rnorm(n = n_obs)
# treatment
prop_score <- plogis(w_1 + w_2)
a <- rbinom(n_obs, 1, prob = prop_score)
# potential outcomes
y_1 <- rnorm(n = n_obs, mean = w_1 + 2 * w_2 + 2 * w_3 + 0.5 * a, sd = 0.1)
y_0 <- rnorm(n = n_obs, mean = w_1 + w_2, sd = 0.1)
# outcome
y <- a * y_1 + (1 - a) * y_0
# assemble the observations in a tibble
tibble(
w_1 = w_1,
w_2 = w_2,
w_3 = w_3,
w_4 = w_4,
w_5 = w_5,
a = a,
y = y
)
}It’s clear from the definition of the potential outcomes that
w_2 and w_3 are treatment effect modifiers:
they interact with the treatment indicator a. A TEM-VIP for
quantifying the strength of the treatment effect modification is
outlined next.
Absolute TEM-VIP
Indexing \(W\) by \(j=1,\ldots,p\), and assuming without loss of generality that \(\mathbb{E}_{P_{X,0}}[W_j] = 0, \mathbb{E}_{P_{X,0}}[W_j^2] > 0\), we can define our first TEM-VIP for the \(j^\text{th}\) confounder as follows:
\[ \Psi^F_j(P_{X,0}) = \frac{\mathbb{E}_{P_{X,0}}[(Y^{(1)}-Y^{(0)})W_j]}{\mathbb{E}_{P_{X,0}}[W_j^2]} =\frac{\mathbb{E}_{P_{X,0}}[(\mathbb{E}_{P_{X,0}}[Y^{(1)}|W]- \mathbb{E}_{P_{X,0}}[Y^{(0)}|W])W_j]} {\mathbb{E}_{P_{X,0}}[W_j^2]} \;. \]
This parameter is called the absolute TEM-VIP. Assuming the expectation of \(\mathbb{E}_{P_{X,0}}[Y^{(1)}|W]- \mathbb{E}_{P_{X,0}}[Y^{(0)}|W]\) conditional on \(W_j\) is linear in \(W_j\), \(\Psi^F_j(P_{X,0})\) is the simple linear regression coefficient produced by regressing the difference in expected potential outcomes against confounder \(W_j\). Even when this relationship is nonlinear, as is almost surely the case in most applications, \(\Psi^F_j(P_{X,0})\) corresponds to the correlation between the difference in potential outcomes and the \(j^\text{th}\) confounder, re-normalized to be on the same scale as the potential outcomes.
Of course, as stated earlier, we don’t generally have access to the full-data random vectors required to perform inference about the above parameter. Luckily, under the identification conditions outlined in Boileau et al. (2022) and Boileau et al. (2023) — namely, that there is no unmeasured confounding and no positivity violations — we can perform inference about the equivalent observed-data parameter:
\[ \Psi_j(P_0) =\frac{\mathbb{E}_{P_0}[(\mathbb{E}_{P_{0}}[Y|A=1,W]- \mathbb{E}_{P_0}[Y|A=0,W])W_j]}{\mathbb{E}_{P_0}[W_j^2]} = \Psi^F_j(P_{X,0})\;. \]
Boileau et al. (2022) and Boileau et al. (2023) derive two nonparametric estimators of this parameter: the one-step and targeted maximum likelihood (TML) estimators. Both require the estimation of two nuisance parameters: the propensity score and the expected outcome conditioned on the treatment and confounders. These estimators are double robust: only one nuisance parameter is required to be consistently estimated to ensure that the one-step and TML estimators are consistent. If both nuisance parameters converge to their true values at a fast enough rate, then the one-step and TML estimators are asymptotically normal. This permits hypothesis testing about \(\Psi_j(P_0)\) using Wald-type confidence intervals. Further details on these estimators and their asymptotic properties are provided in Boileau et al. (2022) and Boileau et al. (2023).
Having defined the parameter and briefly discussed the estimators, we
apply unicate() to recover the treatment effect modifiers,
as defined by the absolute TEM-VIP. Both the one-step and TML estimators
are showcased. The LASSO regression of Tibshirani
(1996) is used to estimate the nuisance parameters.
set.seed(514)
# simulate a random sample
sample_df <- cont_outcome_dgp(n_obs = 500)
# one-step estimator
unihtee(
data = sample_df,
confounders = c("w_1", "w_2", "w_3", "w_4", "w_5"),
modifiers = c("w_1", "w_2", "w_3", "w_4", "w_5"),
exposure = "a",
outcome = "y",
outcome_type = "continuous",
effect = "absolute",
estimator = "onestep",
cond_outcome_estimator = sl3::Lrnr_glmnet$new(
family = "gaussian",
formula = "~ a * w_1 + a * w_2 + a * w_3 + a * w_4 + a * w_5"
),
prop_score_estimator = sl3::Lrnr_glmnet$new(family = "binomial")
)
#> modifier estimate se z p_value ci_lower
#> 1: w_3 1.94918442 0.13370581 14.5781579 0.000000e+00 1.68712103
#> 2: w_2 0.84690336 0.10529721 8.0429801 8.881784e-16 0.64052083
#> 3: w_5 0.14521709 0.10425853 1.3928558 1.636634e-01 -0.05912962
#> 4: w_4 -0.13211440 0.10613296 -1.2448009 2.132049e-01 -0.34013499
#> 5: w_1 -0.05285389 0.08758479 -0.6034596 5.462030e-01 -0.22452009
#> ci_upper p_value_fdr
#> 1: 2.21124781 0.000000e+00
#> 2: 1.05328589 2.220446e-15
#> 3: 0.34956381 2.665062e-01
#> 4: 0.07590619 2.665062e-01
#> 5: 0.11881231 5.462030e-01
# TML estimator
unihtee(
data = sample_df,
confounders = c("w_1", "w_2", "w_3", "w_4", "w_5"),
modifiers = c("w_1", "w_2", "w_3", "w_4", "w_5"),
exposure = "a",
outcome = "y",
outcome_type = "continuous",
effect = "absolute",
estimator = "tmle",
cond_outcome_estimator = sl3::Lrnr_glmnet$new(
family = "gaussian",
formula = "~ a * w_1 + a * w_2 + a * w_3 + a * w_4 + a * w_5"
),
prop_score_estimator = sl3::Lrnr_glmnet$new(family = "binomial")
)
#> modifier estimate se z p_value ci_lower
#> 1: w_3 1.91206506 0.13386174 14.2838801 0.000000e+00 1.64969605
#> 2: w_2 0.81924540 0.10526619 7.7826068 7.105427e-15 0.61292365
#> 3: w_5 0.15243525 0.10422216 1.4625992 1.435771e-01 -0.05184018
#> 4: w_4 -0.11700272 0.10569020 -1.1070347 2.682789e-01 -0.32415551
#> 5: w_1 -0.06328782 0.08714635 -0.7262246 4.677011e-01 -0.23409467
#> ci_upper p_value_fdr
#> 1: 2.17443407 0.000000e+00
#> 2: 1.02556714 1.776357e-14
#> 3: 0.35671067 2.392951e-01
#> 4: 0.09015007 3.353486e-01
#> 5: 0.10751902 4.677011e-01unicate() outputs a table of results summarizing the
TEM-VIP inference procedure for each potential treatment effect
modifier. This table ordered by the nominal p-value of the inference
procedure. Both the one-step and TML estimators correctly identify the
TEM-VIP of w_2 and w_3 as being non-zero based
on the false discovery rate adjusted p-values. Each estimator’s point
estimates are also close to their true parameter values; for every unit
increase in w_2, the one-step and TML estimators suggest
that the average treatment effect increases by \(\approx 0.88\) and \(\approx 0.86\) units, respectively.
Similarly for w_3, the one-step and TML estimators indicate
that a unit increase in this confounder leads to an increase of \(\approx 2.04\) and \(\approx 1.98\) in the average treatment
effect, respectively.
Working with Binary Outcomes
In the binary outcome setting, the full-data and observed-data models and data-generating processes are identical to those of the continuous outcome setting, save that the outcome is binary. That means the absolute TEM-VIP can be used as is with binary-outcome data — the only practical difference is the choice of estimator used for the expected outcome conditional on the treatment and confounders. With that said, a relative TEM-VIP is generally more sensitive and informative in these scenarios. We propose such a parameter in the following section.
Before that, however, we define a data-generating process with two treatment effect modifiers:
bin_outcome_dgp <- function(n_obs) {
# confounders
w_1 <- rpois(n = n_obs, lambda = 3)
w_2 <- rnorm(n = n_obs)
w_3 <- rbinom(n = n_obs, size = 1, prob = 0.5)
# treatment
prop_score <- plogis(-0.5 + 0.25 * w_1 + w_2 + 0.5 * w_3)
a <- rbinom(n_obs, 1, prob = prop_score)
# potential outcomes
y_prob_1 <- plogis(-1 - 3 * w_3)
y_prob_0 <- plogis(1 + 2 * w_2)
y_1 <- rbinom(n_obs, 1, prob = y_prob_1)
y_0 <- rbinom(n_obs, 1, prob = y_prob_0)
# outcome
y <- a * y_1 + (1 - a) * y_0
# assemble the observations in a tibble
tibble(
w_1 = w_1,
w_2 = w_2,
w_3 = w_3,
a = a,
y = y
)
}Here, w_2 and w_3 are a treatment effect
modifier. They interact with the treatment variable, a.
Relative TEM-VIP
Again assuming that the confounders are centered at zero with non-zero variance, the relative TEM-VIP for the \(j^\text{th}\) confounder is \[ \Gamma^F_j(P_{X,0}) =\frac{\mathbb{E}_{P_{X,0}}[(\log\mathbb{E}_{P_{X,0}}[Y^{(1)}|W]- \log\mathbb{E}_{P_{X,0}}[Y^{(0)}|W])W_j]} {\mathbb{E}_{P_{X,0}}[W_j^2]} \;. \]
Assuming that the expectation of the log ratio of the expected conditional potential outcomes conditional on \(W_j\) is linear in \(W_j\), \(\Gamma^F_j(P_{X,0})\) is the simple linear regression coefficient obtained by regressing the log ratio of the expected conditional potential outcomes on \(W_j\). As with \(\Psi^F_j(P_{X,0})\), this parameter can be interpreted as standardized correlation coefficient of the log ratio of expected potential outcomes and the \(j^\text{th}\) confounder.
We highlight that this TEM-VIP isn’t restricted to data-generating processes with binary outcomes; it may also be used when the outcome is non-negative, continuous or discrete random variable.
Again, much like \(\Psi^F_j(P_{X,0})\), \(\Gamma^F_j(P_{X,0})\) isn’t generally ever observed since it is a parameter of the full-data model. Under the conditions of no unmeasured confounding and in the absence of positivity violations, however, an equivalent parameter can be estimated from the observed-data:
\[ \Gamma_j(P_0) =\frac{\mathbb{E}_{P_0}[(\log\mathbb{E}_{P_{0}}[Y|A=1,W]- \log\mathbb{E}_{P_0}[Y|A=0,W])W_j]}{\mathbb{E}_{P_0}[W_j^2]} = \Gamma^F_j(P_{X,0})\;. \]
Nonparametric one-step and TML estimators of this parameter are
derived by Boileau et al. (2023) and
implemented in the unicate package. Just like \(\Psi_j(P_0)\), inference about this
parameter depends on the estimation of the propensity score and the
expected outcome conditional on treatment indicator and confounders.
These estimators are consistent if both nuisance parameters are
consistently estimated, and are asymptotically normal if both converge
to their true values at rate of \(n^{-1/4}\). The latter permits hypothesis
testing based on the Gaussian. Additional details are provided by Boileau et al. (2023).
We now apply these estimators to a random sample generated by
bin_outcome_dgp(). The nuisance parameters are estimated
using Random Forests implemented in ranger package (Breiman 2001; Wright and Ziegler 2017). Note
too that these estimators are fit using cross-fitting (Zheng and van der Laan 2011; Chernozhukov et al.
2017).
# simulate a random sample
sample_df <- bin_outcome_dgp(n_obs = 500)
# one-step estimator
unihtee(
data = sample_df,
confounders = c("w_1", "w_2", "w_3"),
modifiers = c("w_1", "w_2", "w_3"),
exposure = "a",
outcome = "y",
outcome_type = "binary",
effect = "relative",
estimator = "onestep",
cond_outcome_estimator = sl3::Lrnr_ranger$new(),
prop_score_estimator = sl3::Lrnr_ranger$new(),
cross_fit = TRUE
)
#> modifier estimate se z p_value ci_lower ci_upper
#> 1: w_2 -1.1714214 0.2831918 -4.136495 3.526517e-05 -1.7264774 -0.6163655
#> 2: w_3 -2.7446161 1.0922870 -2.512724 1.198029e-02 -4.8854986 -0.6037336
#> 3: w_1 0.2347776 0.1871624 1.254406 2.096946e-01 -0.1320607 0.6016159
#> p_value_fdr
#> 1: 0.0001057955
#> 2: 0.0179704377
#> 3: 0.2096945985
# TML estimator
unihtee(
data = sample_df,
confounders = c("w_1", "w_2", "w_3"),
modifiers = c("w_1", "w_2", "w_3"),
exposure = "a",
outcome = "y",
outcome_type = "binary",
effect = "relative",
estimator = "tmle",
cond_outcome_estimator = sl3::Lrnr_ranger$new(),
prop_score_estimator = sl3::Lrnr_ranger$new(),
cross_fit = TRUE
)
#> modifier estimate se z p_value ci_lower
#> 1: w_3 -3.478340586 0.7746871 -4.48999399 7.122518e-06 -4.9967274
#> 2: w_2 -1.425906898 0.3203599 -4.45095295 8.549007e-06 -2.0538123
#> 3: w_1 0.009961001 0.1692693 0.05884707 9.530739e-01 -0.3218068
#> ci_upper p_value_fdr
#> 1: -1.9599538 1.282351e-05
#> 2: -0.7980015 1.282351e-05
#> 3: 0.3417288 9.530739e-01Both estimators correctly identify w_2 and
w_3 as treatment effect modifiers, as defined by the
relative TEM-VIP. The one-step estimator’s results suggest that for
every unit increase in w_2 and w_3, the log
ratio of the conditional expected potential outcomes is expected to
decrease by \(\approx 1.06\) and \(\approx 3.00\) units, respectively. For the
TML estimator, the log ration of the conditional expected outcomes is
expected to decrease by \(\approx
1.39\) and \(\approx 3.52\) for
every unit increase in w_2 and w_3,
respectively.
Working with Time-to-Event Outcomes
unihtee also provides functionality for estimating
TEM-VIPs for data with right-censored time-to-event outcomes. We detail
the assumed full-data and observed-data models here, then introduce an
absolute and a relative TEM-VIP in subsequent sections.
Consider \(n\) i.i.d. random vectors \(X = (W, A, C^{(0)}, C^{(1)}, T^{(0)}, T^{(1)}) \sim P_{X,0} \in \mathcal{M}_X\). \(W, A, P_{X,0},\) and \(\mathcal{M}_X\) are defined as in the continuous and binary outcome models. \(C^{(0)}\) and \(C^{(1)}\) correspond to the censoring times that would be observed under the control and treatment conditions, respectively. Similarly, \(T^{(0)}\) and \(T^{(1)}\) are the event times that would be observed under either treatment assignment.
Again, \(P_{X,0}\) is generally unobservable. We instead have access to random observations \(O = (W, A, \Delta, T) \sim P_0 \in \mathcal{M}\), where \(\Delta\) is a censoring indicator and \(T = A(\Delta C^{(1)} + (1-\Delta)T^{(1)}) + (1-A)(\Delta C^{(0)} + (1-\Delta) T^{(0)})\). \(W, A, P_0,\) and \(\mathcal{M}\) are defined as in the previous observed-data models.
Now, it’s worth introducing the full-data and observed-data conditional survival functions since both TEM-VIPs introduced in this section rely on them. The full-data survival function is defined as \(S_{P_{X,0}}(t|A,W) = \mathbb{P}_{P_{X,0}}[T^{(A)} > t | W]\). This parameter’s observed-data counterpart is defined analogously as \(S_0(t|A,W) = P_{P_0}[T > t |A, W]\).
Finally, we define a time-to-event data-generating process with nine
equidistant time points and three potential confounders,
w_1, w_2, and w_3. Of these
confounders, only w_1 is a treatment effect modifier. Note
that unihtee() requires a “wide” data format when working
time-to-event outcome data. That is, each observation is represented by
a single row in the data. An example is provided following this code
chunk.
# define hazard functions
cond_cens_hazard <- function(time, a, w_1, w_2, w_3) {
(time < 9) / (1 + exp(4 + 0.2 * (w_1 + w_2) - a))
}
cond_surv_hazard <- function(time, a, w_1, w_2, w_3) {
(time < 9) / (1 + exp(2 + 3 * a * w_1)) + (time == 9)
}
# failure time simulator
failure_time_sim <- function(n_obs, a, w_1, w_2, w_3) {
sapply(
seq_len(n_obs),
function(obs) {
failure_time <- NA
for (t in 1:9) {
prob <- cond_surv_hazard(t, a, w_1[obs], w_2[obs], w_3[obs])
status <- rbinom(1, 1, prob)
if (status == 1) {
failure_time <- t
break
}
}
return(failure_time)
}
)
}
# censoring time simulator
censor_time_sim <- function(n_obs, a, w_1, w_2, w_3) {
sapply(
seq_len(n_obs),
function(obs) {
censor_time <- NA
for (t in 1:9) {
prob <- cond_cens_hazard(t, a, w_1[obs], w_2[obs], w_3[obs])
status <- rbinom(1, 1, prob)
if (status == 1) {
censor_time <- t
break
}
}
if (is.na(censor_time)) censor_time <- 10
return(censor_time)
}
)
}
# time-to-event outcome data-generating process
tte_outcome_dgp <- function(n_obs) {
# confounders
w_1 <- rnorm(n = n_obs)
w_2 <- rnorm(n = n_obs)
w_3 <- rnorm(n = n_obs)
# treatment
prop_score <- plogis(-1 - w_1 + w_2 - w_3)
a <- rbinom(n_obs, 1, prob = prop_score)
# generate the failure events for t = 1 to 9
failure_time_1 <- failure_time_sim(n_obs, 1, w_1, w_2, w_3)
failure_time_0 <- failure_time_sim(n_obs, 0, w_1, w_2, w_3)
# generate the censoring events for t = 1 to 9
censor_time_1 <- censor_time_sim(n_obs, 1, w_1, w_2, w_3)
censor_time_0 <- censor_time_sim(n_obs, 0, w_1, w_2, w_3)
# compile the failure and censoring times
failure_time <- sapply(
seq_len(n_obs),
function(obs) {
if (a[obs] == 1) failure_time_1[obs] else failure_time_0[obs]
}
)
censor_time <- sapply(
seq_len(n_obs),
function(obs) {
if (a[obs] == 1) censor_time_1[obs] else censor_time_0[obs]
}
)
# determine the observed time-to-event and censoring indicator
time <- sapply(
seq_len(n_obs),
function(obs) {
if (censor_time[obs] < failure_time[obs]) {
censor_time[obs]
} else {
failure_time[obs]
}
}
)
censoring_indicator <- sapply(
seq_len(n_obs),
function(obs) if (time[obs] == censor_time[obs]) 1 else 0
)
# assemble the tibble
tibble(
w_1 = w_1,
w_2 = w_2,
w_3 = w_3,
a = a,
time = time,
censoring_indicator = censoring_indicator
)
}
tte_outcome_dgp(10)
#> # A tibble: 10 × 6
#> w_1 w_2 w_3 a time censoring_indicator
#> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
#> 1 0.185 0.369 0.753 0 9 0
#> 2 -1.19 1.13 -1.96 1 1 0
#> 3 -0.926 0.551 -1.23 1 1 0
#> 4 1.59 -0.0865 0.725 0 4 0
#> 5 -1.59 1.13 -1.88 1 1 0
#> 6 -1.07 1.31 0.492 0 7 0
#> 7 -0.861 -1.73 0.430 0 6 1
#> 8 -1.57 1.30 1.08 1 1 0
#> 9 -0.600 0.0367 -1.40 0 1 0
#> 10 -0.410 1.30 -0.238 0 6 0Absolute TEM-VIP
Boileau et al. (2023) Propose an absolute TEM-VIP based on the restricted mean survival time at some pre-specified time \(t\). Again assuming without loss of generality that the \(\mathbb{E}_{P_{X,0}}[W_j] = 0\) and that the confounders have variance bounded away from zero, this TEM-VIP is defined as \[ \Psi_j^F(P_{X,0};t) = \frac{\mathbb{E}_{P_{X,0}}\left[W_j \int_0^t \{ S_{P_{X,0}}(u|1, W) - S_{P_{X,0}}(u|0, W) \} du \right]} {\mathbb{E}_{P_{X,0}}[W_j^2]} \; . \]
“\(\Psi\)” is re-used to highlight that this is an absolute TEM-VIP. As in the continuous outcome scenario, this parameter captures the correlation between the conditional difference in the restricted mean survival time and the \(j^\text{th}\) confounder, re-normalized to be on the same scale as the outcome. Put another way, this parameter identifies confounders driving the largest difference in truncated survival times across treatment conditions.
As before, this full-data parameter is equal to the following observed-data parameter under the assumptions of (1) no unmeasured confounding, (2) censoring mechanism positivity, and (3) treatment assignment positivity: \[ \Psi_j(P_0;t) = \frac{\mathbb{E}_{P_0}\left[W_j \int_0^t \{ S_0(u|1, W) - S_0(u|0, W) \} du \right]} {\mathbb{E}_{P_0}[W_j^2]} \; . \]
As in the continuous and binary outcome scenarios, Boileau et al. (2023) derived two estimators of this absolute effect parameter: a one-step estimator and a TML estimator. Both rely on the accurate estimation of three nuisance parameters: the conditional event hazard function, the conditional censoring hazard function, and the propensity score. If either the conditional event hazard function or the conditional censoring hazard function and the propensity score are consistently estimated, then these nonparametric estimators are consistent. Further, if all of these nuisance parameters are converge to their true values at a rate of \(n^{-1/4}\), then these estimators are asymptotically normally distributed. This permits Gaussian-based hypothesis testing through the use of Wald-type confidence intervals. See Boileau et al. (2023) for more details.
We now apply these estimators to a random sample from the data-generating process defined above, choosing as target parameter \(\Psi(P_0; 8)\). LASSO regressions are used to estimate the propensity scores and censoring hazard functions, and, as demonstration of these methods’ flexibility, estimate the failure hazard function with XGBoost (T. Chen and Guestrin 2016).
# simulate a random sample
sample_df <- tte_outcome_dgp(n_obs = 250)
# one-step estimator
unihtee(
data = sample_df,
confounders = c("w_1", "w_2", "w_3"),
modifiers = c("w_1", "w_2", "w_3"),
exposure = "a",
outcome = "time",
censoring = "censoring_indicator",
time_cutoff = 8,
outcome_type = "time-to-event",
effect = "absolute",
estimator = "onestep",
prop_score_estimator = sl3::Lrnr_glmnet$new(),
failure_hazard_estimator = sl3::Lrnr_xgboost$new(),
censoring_hazard_estimator = sl3::Lrnr_glmnet$new()
)
#> modifier estimate se z p_value ci_lower ci_upper
#> 1: w_1 1.65578106 0.3100170 5.3409367 9.246754e-08 1.0481478 2.2634143
#> 2: w_2 -0.27580480 0.2378789 -1.1594338 2.462794e-01 -0.7420474 0.1904378
#> 3: w_3 0.04483202 0.2631764 0.1703497 8.647351e-01 -0.4709937 0.5606577
#> p_value_fdr
#> 1: 2.774026e-07
#> 2: 3.694191e-01
#> 3: 8.647351e-01
# TML estimator
unihtee(
data = sample_df,
confounders = c("w_1", "w_2", "w_3"),
modifiers = c("w_1", "w_2", "w_3"),
exposure = "a",
outcome = "time",
censoring = "censoring_indicator",
time_cutoff = 8,
outcome_type = "time-to-event",
effect = "absolute",
estimator = "tmle",
prop_score_estimator = sl3::Lrnr_glmnet$new(),
failure_hazard_estimator = sl3::Lrnr_xgboost$new(),
censoring_hazard_estimator = sl3::Lrnr_glmnet$new()
)
#> modifier estimate se z p_value ci_lower ci_upper
#> 1: w_1 1.24516252 0.3137506 3.9686383 7.228451e-05 0.6302114 1.86011363
#> 2: w_2 -0.37336556 0.2401380 -1.5547961 1.199946e-01 -0.8440360 0.09730484
#> 3: w_3 -0.07052526 0.2658042 -0.2653278 7.907569e-01 -0.5915015 0.45045101
#> p_value_fdr
#> 1: 0.0002168535
#> 2: 0.1799919603
#> 3: 0.7907569131These results suggest every unit increase in w_1 results
in an increase of \(\approx 1.66\) —
according to the one-step estimator — or \(\approx 1.25\) — according to the TML
estimator — time units for the difference in mean survival times,
truncated at time eight.
Relative TEM-VIP
We next consider a relative TEM-VIP for data with right-censored time-to-event outcomes. Again assuming that the potential treatment effect modifiers are centered and have variance bounded away from zero, we define this parameter as follows for the \(j^\text{th}\) confounder: \[ \Gamma_j^F(P_{X,0};t) = \frac{\mathbb{E}_{P_{X,0}}\left[W_j (\log S_{P_{X,0}}(t|1, W) - \log S_{P_{X,0}}(t|0, W))\right]} {\mathbb{E}_{P_{X,0}}[W_j^2]} \; . \]
For some pre-specified time point \(t\), this TEM-VIP represents the standardized correlation coefficient of \(W_j\) and the log ratio of conditional survival times under each treatment condition. As with “\(\Psi\)”, we re-use “\(\Gamma\)” to highlight that this is a relative effect parameter. Its observed-data counterpart is given by \[ \Gamma_j^F(P_0;t) = \frac{\mathbb{E}_{P_0}\left[W_j (\log S_0(t|1, W) - \log S_0(t|0, W))\right]} {\mathbb{E}_{P_0}[W_j^2]} \; , \] again assuming that there is no unmeasured confounding and that there are no positivity violations of the censoring or treatment assignment mechanisms.
Boileau et al. (2023) derived one-step and TML estimators of this parameter, again requiring the estimation of three nuisance parameters: the propensity score and the conditional failure and censoring event hazards. These one-step and TML estimators are consistent when all nuisance parameters are estimated consistently, and are asymptotically linear when these nuisance parameters converge to their true values at a fast-enough rate (Boileau et al. 2023).
We apply the nonparametric estimators of Boileau et al. (2023) to a random sample
generated by tte_outcome_dgp(), taking as target of
inference \(\Gamma_j(P_0; 5)\). The
nuisance parameters are estimated using the same strategy employed in
prior absolute TEM-VIP example.
# simulate a random sample
sample_df <- tte_outcome_dgp(n_obs = 1000)
# one-step estimator
unihtee(
data = sample_df,
confounders = c("w_1", "w_2", "w_3"),
modifiers = c("w_1", "w_2", "w_3"),
exposure = "a",
outcome = "time",
censoring = "censoring_indicator",
time_cutoff = 5,
outcome_type = "time-to-event",
effect = "relative",
estimator = "onestep",
prop_score_estimator = sl3::Lrnr_glmnet$new(),
failure_hazard_estimator = sl3::Lrnr_xgboost$new(),
censoring_hazard_estimator = sl3::Lrnr_glmnet$new(),
cross_fit = TRUE
)
#> modifier estimate se z p_value ci_lower ci_upper
#> 1: w_1 3.2868215 0.7571221 4.3412042 1.417039e-05 1.802862 4.7707807
#> 2: w_3 -1.6368035 1.2089918 -1.3538582 1.757816e-01 -4.006428 0.7328204
#> 3: w_2 -0.5309253 0.7692443 -0.6901908 4.900742e-01 -2.038644 0.9767935
#> p_value_fdr
#> 1: 4.251118e-05
#> 2: 2.636724e-01
#> 3: 4.900742e-01
# TML estimator
unihtee(
data = sample_df,
confounders = c("w_1", "w_2", "w_3"),
modifiers = c("w_1", "w_2", "w_3"),
exposure = "a",
outcome = "time",
censoring = "censoring_indicator",
time_cutoff = 5,
outcome_type = "time-to-event",
effect = "relative",
estimator = "tmle",
prop_score_estimator = sl3::Lrnr_glmnet$new(),
failure_hazard_estimator = sl3::Lrnr_xgboost$new(),
censoring_hazard_estimator = sl3::Lrnr_glmnet$new(),
cross_fit = TRUE
)
#> modifier estimate se z p_value ci_lower ci_upper
#> 1: w_1 3.2826140 1.139237 2.88141401 0.003958952 1.049709 5.515519
#> 2: w_2 0.2824440 1.057572 0.26706826 0.789416602 -1.790398 2.355286
#> 3: w_3 -0.2385342 6.600093 -0.03614104 0.971169901 -13.174716 12.697648
#> p_value_fdr
#> 1: 0.01187686
#> 2: 0.97116990
#> 3: 0.97116990Again, both the one-step and TML estimators successfully discern the
treatment effect modifier from the other potential confounders when
using a false discovery rate cutoff of \(5\%\). These results indicate that a unit
increase in w_1 results in an increase of the log ratio of
conditional survivals equal to either \(\approx 3.29\) or \(\approx 3.28\), depending on whether one
uses the one-step estimator or the TML estimator.
Additional Features and Notes
Working with Randomized Experiment Data
While unihtee’s estimators are applied to (simulated)
observational data in the previous section’s examples, they are also
suitable for data generated by randomized experiments. The propensity
scores used to randomize observations should be saved data object passed
to unihtee(), and the name of the variable containing these
scores should be indicated in the propensity_score_values
argument.
Providing the known propensity scores to the absolute TEM-VIP estimators produces some desirable results. In the continuous and binary outcome settings, providing known treatment assignment probabilities guarantees that the estimators are asymptotically linear, even when the expected conditional outcome estimator is misspecified. When the outcome is a right-censored time-to-event variable, asymptotic normality of the estimators is guaranteed so long as the conditional censoring hazard is estimated consistently at a rate of \(n^{-1/4}\). Additional details and discussions are provided in Boileau et al. (2022) and Boileau et al. (2023).
Cross-Fitting
Cross-fitting, popularized by Zheng and van
der Laan (2011) and Chernozhukov et al.
(2017), uses cross-validation to fit the nuisance parameters and
estimate the target parameter. Doing so relaxes the entropy constraints
required for asymptotic normality of the one-step and TML estimators
implemented in the unihtee package. This technique is also
known to improve these estimators finite-sample properties, such as
reducing their false discovery rate (see, for example, Hejazi et al.).
Users may employ the cross-fitted version of any estimator
implemented in unihtee by setting unihtee()’s
cross_fit argument to TRUE. Only K-fold
cross-validation is currently supported. The number of folds defaults to
five, though this can be modified using the cross_fit_folds
argument. Examples are provided in the binary outcome and time-to-event
outcome relative TEM-VIP sections of this vignette, as well as in the
subsection on Super Learners that follows.
sl3 Integration: Using Super Learners
unihtee uses the estimators implemented in the
sl3 (Coyle et al. 2021) to
estimate nuisance parameters. While individual estimators may be used,
as is done in the examples above, sl3’s main contribution
are facilities for constructing Super Learners (Laan, Polley, and Hubbard 2007). Super Learners
are ensemble algorithms, meaning they are the convex combination of
individual supervised learning algorithms. These individual algorithms
are combined on the basis of a user-defined risk, and this combination
is asymptotically optimal with respect to said risk under minimal
conditions about the data-generating process. A complete review on Super
Learners and sl3 is outside the scope of this vignette. We
invite the interested reader to the tutorials available on the
sl3 package’s website and in Chapter 6 of the tlverse
handbook (Phillips 2023).
Employing Super Learners to estimate nuisance parameters decreases the risk of model-misspecification — assuming, of course, that a diverse collection individual learning algorithms is considered. In turn, there is an increased likelihood that the resulting one-step and TML estimators are asymptotically normal.
We revisit our continuous outcome example, this time using Super Learners to estimate the propensity score and the conditional expected outcomes. Both Super Learners are composed of penalized GLMs, multivariate adaptive regression splines (Friedman 1991), and Random Forests. The propensity score’s Super Learner uses the log loss to form the convex combination of individual learning algorithms, while the conditional expected outcomes’ uses the squared error loss. Additionally, we employ 2-fold cross-fitting.
# simulate a random sample
sample_df <- cont_outcome_dgp(n_obs = 500)
# define the Super Learner for the propensity score
sl_bin <- Lrnr_sl$new(
learners = list(
Lrnr_glmnet$new(family = "binomial", alpha = 0),
Lrnr_glmnet$new(family = "binomial", alpha = 0),
Lrnr_glmnet$new(family = "binomial", alpha = 0.5),
Lrnr_earth$new(glm = list(family = "binomial")),
Lrnr_ranger$new()
),
metalearner = make_learner(
Lrnr_solnp, metalearner_logistic_binomial,
loss_squared_error
)
)
# define the Super Learner for the conditional expected outcomes
sl_cont <- Lrnr_sl$new(
learners = list(
Lrnr_glmnet$new(
family = "gaussian",
formula = "~ a * w_1 + a * w_2 + a * w_3 + a * w_4 + a * w_5"
),
Lrnr_glmnet$new(
family = "gaussian",
alpha = 0,
formula = "~ a * w_1 + a * w_2 + a * w_3 + a * w_4 + a * w_5"
),
Lrnr_glmnet$new(
family = "gaussian",
alpha = 0.5,
formula = "~ a * w_1 + a * w_2 + a * w_3 + a * w_4 + a * w_5"
),
Lrnr_earth$new(
glm = list(family = "gaussian"),
degree = 1
),
Lrnr_earth$new(
glm = list(family = "gaussian"),
degree = 2
),
Lrnr_earth$new(
glm = list(family = "gaussian"),
formula = "~ a * w_1 + a * w_2 + a * w_3 + a * w_4 + a * w_5"
),
Lrnr_ranger$new()
),
metalearner = make_learner(
Lrnr_solnp, metalearner_linear, sl3::loss_squared_error
)
)
# one-step estimator
unihtee(
data = sample_df,
confounders = c("w_1", "w_2", "w_3", "w_4", "w_5"),
modifiers = c("w_1", "w_2", "w_3", "w_4", "w_5"),
exposure = "a",
outcome = "y",
outcome_type = "continuous",
effect = "absolute",
estimator = "onestep",
cond_outcome_estimator = sl_cont,
prop_score_estimator = sl_bin,
cross_fit = TRUE,
cross_fit_folds = 2
)
#> modifier estimate se z p_value ci_lower ci_upper
#> 1: w_2 1.01322541 0.10199927 9.9336537 0.0000000 0.81330684 1.2131440
#> 2: w_3 1.98347638 0.13831712 14.3400640 0.0000000 1.71237482 2.2545779
#> 3: w_4 0.15980553 0.09755231 1.6381522 0.1013899 -0.03139699 0.3510081
#> 4: w_1 -0.01980657 0.10279772 -0.1926752 0.8472133 -0.22129010 0.1816770
#> 5: w_5 0.01434656 0.09377386 0.1529911 0.8784053 -0.16945020 0.1981433
#> p_value_fdr
#> 1: 0.0000000
#> 2: 0.0000000
#> 3: 0.1689832
#> 4: 0.8784053
#> 5: 0.8784053
# TML estimator
unihtee(
data = sample_df,
confounders = c("w_1", "w_2", "w_3", "w_4", "w_5"),
modifiers = c("w_1", "w_2", "w_3", "w_4", "w_5"),
exposure = "a",
outcome = "y",
outcome_type = "continuous",
effect = "absolute",
estimator = "tmle",
cond_outcome_estimator = sl_cont,
prop_score_estimator = sl_bin,
cross_fit = TRUE,
cross_fit_folds = 2
)
#> modifier estimate se z p_value ci_lower ci_upper
#> 1: w_2 0.98762264 0.10454871 9.4465311 0.0000000 0.7827072 1.1925381
#> 2: w_3 1.98247437 0.13797716 14.3681347 0.0000000 1.7120391 2.2529096
#> 3: w_4 0.15711411 0.09755985 1.6104382 0.1073022 -0.0341032 0.3483314
#> 4: w_1 -0.03473531 0.10365604 -0.3351016 0.7375484 -0.2379011 0.1684305
#> 5: w_5 0.01402747 0.09395342 0.1493024 0.8813150 -0.1701212 0.1981762
#> p_value_fdr
#> 1: 0.0000000
#> 2: 0.0000000
#> 3: 0.1788371
#> 4: 0.8813150
#> 5: 0.8813150Though this procedure is more robust, the resulting estimates and hypothesis tests are to those of the previous continuous outcome, absolute TEM-VIP example.
Session information
sessionInfo()
#> R version 4.2.2 (2022-10-31)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 22.04.2 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] unihtee_0.0.1 earth_5.3.2 plotmo_3.6.2 TeachingDemos_2.12
#> [5] plotrix_3.8-2 Formula_1.2-5 ranger_0.15.1 xgboost_1.7.5.1
#> [9] glmnet_4.1-7 Matrix_1.5-4 Rsolnp_1.16 sl3_1.4.4
#> [13] dplyr_1.1.1
#>
#> loaded via a namespace (and not attached):
#> [1] colorspace_2.1-0 class_7.3-21 rprojroot_2.0.3
#> [4] fs_1.6.1 listenv_0.9.0 prodlim_2023.03.31
#> [7] fansi_1.0.4 lubridate_1.9.2 codetools_0.2-19
#> [10] splines_4.2.2 R.methodsS3_1.8.2 cachem_1.0.7
#> [13] knitr_1.42 jsonlite_1.8.4 pROC_1.18.0
#> [16] caret_6.0-94 R.oo_1.25.0 compiler_4.2.2
#> [19] backports_1.4.1 assertthat_0.2.1 fastmap_1.1.1
#> [22] rstackdeque_1.1.1 cli_3.6.1 visNetwork_2.1.2
#> [25] htmltools_0.5.5 prettyunits_1.1.1 tools_4.2.2
#> [28] igraph_1.4.2 gtable_0.3.3 glue_1.6.2
#> [31] reshape2_1.4.4 Rcpp_1.0.10 jquerylib_0.1.4
#> [34] pkgdown_2.0.7 vctrs_0.6.1 nlme_3.1-162
#> [37] iterators_1.0.14 timeDate_4022.108 origami_1.0.7
#> [40] gower_1.0.1 xfun_0.38 stringr_1.5.0
#> [43] globals_0.16.2 rbibutils_2.2.13 delayed_0.4.0
#> [46] timechange_0.2.0 lifecycle_1.0.3 future_1.32.0
#> [49] MASS_7.3-58.3 scales_1.2.1 ipred_0.9-14
#> [52] ragg_1.2.5 hms_1.1.3 parallel_4.2.2
#> [55] BBmisc_1.13 yaml_2.3.7 memoise_2.0.1
#> [58] ggplot2_3.4.2 sass_0.4.5 rpart_4.1.19
#> [61] stringi_1.7.12 desc_1.4.2 foreach_1.5.2
#> [64] randomForest_4.7-1.1 checkmate_2.1.0 hardhat_1.3.0
#> [67] lava_1.7.2.1 truncnorm_1.0-9 shape_1.4.6
#> [70] Rdpack_2.4 rlang_1.1.0 pkgconfig_2.0.3
#> [73] systemfonts_1.0.4 evaluate_0.20 lattice_0.21-8
#> [76] ROCR_1.0-11 purrr_1.0.1 recipes_1.0.5
#> [79] htmlwidgets_1.6.2 tidyselect_1.2.0 parallelly_1.35.0
#> [82] plyr_1.8.8 magrittr_2.0.3 R6_2.5.1
#> [85] generics_0.1.3 pillar_1.9.0 withr_2.5.0
#> [88] survival_3.5-5 abind_1.4-5 nnet_7.3-18
#> [91] tibble_3.2.1 future.apply_1.10.0 crayon_1.5.2
#> [94] uuid_1.1-0 utf8_1.2.3 rmarkdown_2.21
#> [97] imputeMissings_0.0.3 progress_1.2.2 grid_4.2.2
#> [100] data.table_1.14.8 ModelMetrics_1.2.2.2 digest_0.6.31
#> [103] R.utils_2.12.2 textshaping_0.3.6 stats4_4.2.2
#> [106] munsell_0.5.0 bslib_0.4.2