
lmtp: An R Package for Estimating the Causal Effects of Modified Treatment Policies
We present the lmtp R package for causal inference from longitudinal observational or randomized studies. This package implements the estimators of Díaz et al. (2021) for estimating general nonparametric causal effects based on modified treatment policies. Modified treatment policies generalize static and dynamic interventions, making lmtp and allpurpose package for nonparametric causal inference in observational studies. The methods provided can be applied to both pointtreatment and longitudinal settings, and can account for timevarying exposure, covariates, and right censoring thereby providing a very general tool for causal inference. Additionally, two of the provided estimators are based on flexible machine learning regression algorithms, and avoid bias due to parametric model misspecification while maintaining valid statistical inference.
Causal inference, longitudinal data, nonparametric, modified treatment policies
1. Introduction
Modified treatment policies (MTPs) are a recently developed class of interventions that generalize many interesting and widely used parameters in causal inference. They generalize static and dynamic interventions for categorical exposure variables, but also provide a framework to define interesting parameters for continuous or multivariate exposures. The fundamental ideas date back to Robins et al. (2004). In a series of recent articles, these methods have been further developed to allow for more flexible data types and estimation algorithms (Díaz and van der Laan, 2012; Haneuse and Rotnitzky, 2013; Young et al., 2014). Specifically, two recent papers (Díaz et al., 2021, 2022) have generalized the MTP framework to account for (1) timevarying treatment and covariates, (2) informative rightcensoring, (3) competing risks, (4) machine learning estimation of nuisance parameters.
In their general formulation, MTPs are hypothetical interventions where the postintervention value of exposure can depend on the actual observed exposure level and the [End Page 103] unit’s history. As such, MTPs are useful to assess the effect of continuous treatments. For example, Haneuse and Rotnitzky (2013) assess the effect of reducing surgery time by a predetermined amount (e.g., 5 minutes) for lung cancer patients, where the reduction is carried out only for those patients for whom the intervention is feasible. Furthermore, as explained above, MTPs generalize many important causal inference estimands, such as the effect of a dynamic treatment rule in which the treatment level is assigned as a function of a unit’s history. For example, dynamic treatment rules, a particular case of MTPs, may be used to estimate the effect of policies such as switching HIV treatment once the CD4 Tcell count goes below a predetermined threshold (Petersen et al., 2014). MTPs also generalize many interesting causal effects such as the average treatment effects, the causal risk ratio, and causal odds ratio. In this article we describe how lmtp can be used for estimating the causal effects of MTPs, and present examples on the use of the software for several of the above cases.
Flexible and general software for estimation of causal effects in longitudinal studies is currently limited. The ipw and gfoRmula (van der Wal and Geskus, 2011; McGrath et al., 2020) packages provide routines for estimating causal effects using inverse probability weighting (IPW) and the parametric gformula respectively. The validity of these methods, however, requires making unrealistic parametric assumptions of the underlying data generating distributions and are not doublyrobust. The ltmle and survtmle (Lendle et al., 2017; Benkeser and Hejazi, 2019) packages implement a doublyrobust method for estimating causal effects from longitudinal data, but these packages do not support continuous valued exposures. The txshift (Hejazi and Benkeser, 2020) package provides a method for estimating effects of continuous, pointtreatment exposures using MTPs defined as additive shifts. In contrast, lmtp can be used for estimating causal effects for binary, categorical, and continuous valued exposures from longitudinal data using stateoftheart methodologies that may incorporate machine learning and avoid parametric modeling assumptions.
The goal of this paper is to introduce the lmtp package and to illustrate its capabilities. The package implements four methods for estimating the effects of MTPs. The first two estimators are simple estimation strategies based on gcomputation and inverse probability weighting. The sampling distributions of these estimators under machine learning of the nuisance parameters are unknown, and therefore we do not recommend their use in practice. The other two estimators, a targeted minimumloss based estimator (TMLE)(Laan and Rose, 2011; Laan and Rubin, 2006) and a sequentially doublyrobust estimator (SDR) (Buckley and James, 1979; Fan and Gijbels, 1994; van der Laan and Dudoit, 2003; Rotnitzky et al., 2006; Rubin and Laan, 2006; Kennedy et al., 2017), are multiplyrobust and asymptotically normal under conditions, even when the nuisance parameters are fitted with machine learning. TMLE and SDR are implemented using crossfitting to allow for the use of flexible machine learning regression methodology (Díaz et al., 2021). The package may be download from CRAN at cran.rproject.org/package=lmtp.
2. Notation and modified treatment policies
2.1 Data structure
We start the paper with a brief introduction to MTPs. The interested reader is encouraged to consult Díaz et al. (2021) and Díaz et al. (2022) for more technical details and theoretical [End Page 104] results. We will use the notation of Díaz et al. (2021) with only slight modifications. Let i be the index of an observation from a data set with n total units and t be the index of time for a total number of time points τ . The observed data for observation i may be denoted as
where W denotes baseline covariates, L_{t} denotes timevarying covariates, A_{t} denotes a vector of exposure and/or censoring variables and Y denotes an outcome measured at the end of study followup. We observe n i.i.d. copies of Z with distribution P. Upper case reference (A_{t}) to a variable indicates a random variable while lower case reference (a_{t}) indicates a realization of that random variable. If rightcensoring exists, A_{t} can be adapted so that A_{t} = (A_{1,}_{t}, A_{2,}_{t}) where A_{1,}_{t} equals one if an observation is still in the study at time t and zero otherwise, and A_{2,}_{t} denotes the exposure at time t. We use an overbar to indicate the history of a variable up until time t. We then use H_{t} = (L̄_{t},Ā_{t−}_{1}) to denote the history of all variables up until just before A_{t}.
2.2 Modified treatment policies
We use the potential outcomes framework to define the causal effect of interest using our established data structure. We consider a hypothetical policy where Ā is set to a regime d defined as A^{d}_{t} = d_{t}(A_{t}, H^{d}_{t}), where H^{d}_{t} = (L̄_{t},Ā^{d}_{t−}_{1}), for a set of usergiven regimes d_{t} : t ∈ {1, ..., τ }. The defining characteristic that makes regime d_{t} a modified treatment policy is that it depends on the natural value of treatment A_{t}, that is, the value that the treatment would have taken under no intervention. However, when the function d_{t} only depends on H_{t}, the LMTP reduces to the dynamic treatment regimes studied in the literature. Furthermore, when d_{t} is a constant and does not depend on either A_{t} or H_{t}, then LMTPs reduce to the conventional static rules studied in the causal inference literature (e.g., Bang and Robins, 2005; van der Laan and Gruber, 2011). Below we present examples of all these interventions.
First, consider a study of the effect of physical activity on mortality in the elderly. Assume that each patient is monitored at several time points, and that a measure of physical activity such as the metabolic equivalent of task (MET) (Mendes et al., 2018) is measured together with a number of lifestyle, health status, and demographic variables. In this setup, a natural question to ask would be “what is the effect on mortality of an intervention that increases physical activity by δ units for patients whose socioeconomic and health status allows it?” Formally, consider a longitudinal study with losstofollowup. Let A_{t} = (A_{1,}_{t}, A_{2,}_{t}) where A_{1,}_{t} equals one if an observation is still in the study at time t and zero otherwise, and A_{2,}_{t} denote a continuous exposure at time t that can be changed through some intervention. A modified treatment policy that increases A_{2,}_{t}, whenever it is feasible to do so, can be defined as
where u_{t}(h_{t}) defines the maximum level of physical activity allowed for a patient with characteristics h_{t}. Note that we also consider an intervention on A_{1,}_{t} because we are interested [End Page 105] in a hypothetical world where there is no losstofollowup. In this case the hypothetical exposure after intervention, A^{d}_{t} depends on the actually observed exposure, A_{t}. This is in contrast to a deterministic intervention where A^{d}_{t} would be set to some prespecified value with probability one.
For dynamic treatment rules, consider a hypothetical longitudinal study where two different antiviral treatments are administered to HIV positive patients. Sometimes an antiviral drug works at first, until the virus develops resistance, at which point it is necessary to change the treatment regime. Assume we are interested in assessing a policy with two treatments encoded as A_{t} ∈ {0, 1}, and we want to assess the effect of a regime that would switch the antiviral treatment as soon as the CD4 T cell count drops bellow 300. Let A_{t} = (A_{1,}_{t}, A_{2,}_{t}) where A_{1,}_{t} equals one if an observation is still in the study at time t and zero otherwise, and A_{2,}_{t} denotes the treatment arm at time t. Let L_{t} denote the CD4 T cell count at time t. In this case, one may decide to assess the effect of the rule
In contrast to the previous rule (2), the dynamic treatment rule (3) does not depend on the natural value of treatment at time t, it only depends on the history. The software and methods presented here handle both cases seamlessly.
In the case of a single time point setting where the data structure is Z = (W, A, Y ), it follows trivially from the above definitions that the average treatment effect from a crosssectional study, defined as E[Y (1) − Y (0)], can be estimated using MTPs by simply letting τ = 1 and contrasting two MTPs d(A) = 1 and d(A) = 0. The lmtp package presented in this article allows the contrast of different MTPs using differences, ratios, and odds ratios. We provide examples of these contrasts in §4 below.
In what follows we focus on estimating the the causal effect of MTP d on outcome Y, using lmtp, through the causal parameter
where Y (Ā^{d}) is the counterfactual outcome in a world, where possibly contrary to fact, each entry ofĀ was modified according to the MTP d. When Y is continuous, θ is the mean population value of Y under MTP d; when Y is dichotomous, θ is the population proportion of event Y under MTP d. Similarly, when Y is the indicator of an event by end of the study, θ is defined as the cumulative incidence of Y under MTP d.
2.3 Identification
The ability to estimate θ depends on the ability to identify an expression for θ as a function of the data generating distribution P using only the observed data Z and not counterfactual variables Y^{d}. A full review of these identification assumptions is outside the scope of this article but is presented in our technical paper (Díaz et al., 2021). Briefly, the following standard assumptions must hold
Assumption 1 (Consistency)Ā =ā =⇒ Y = Y (ā) for allā ∈ supp Ā [End Page 106]
Assumption 2 (Exchangeability) A_{t} ╨ Y (ā)H_{t} for allā ∈ suppĀ and t ∈ {1, ..., τ }
Assumption 3 (Positivity) If (a_{t}, h_{t}) ∈ supp{A_{t}, H_{t}} then (d(a_{t}, h_{t}), h_{t}) ∈ supp{A_{t}, H_{t}} for t ∈ {1, ..., τ }
The consistency assumption states that the potential outcome for an observation under their observed exposure is the value of the outcome that we did actually observe. The exchangeability assumption is often also referred to as the nounmeasured confounding assumption; it is satisfied if all common causes of the exposure and outcome are measured and adjusted for.
Of particular importance to this article is the positivity assumption which states that the distribution of the exposure under the MTP is supported in the data. Concretely, in a study with a continuous exposure and losstofollowup, the positivity assumption states that if an observation with covariate history h_{t} and exposure a_{t} who was not losttofollowup at time t exists then there is also an observation with covariate history h_{t} who was not losttofollowup at time t but whose exposure was observed as d(a_{t}, h_{t}) that also exists. This assumption is often an issue when working with continuous exposures and/or multiple time points in the context of static or dynamic interventions. Note that, if the treatment is binary and d_{t}(a_{t}, h_{t}) = 1, the positivity assumption reduces to “if h_{t} ∈ supp{H_{t}} then (1, h_{t}) ∈ supp{A_{t}, H_{t}}”, which is the standard positivity assumption requiring that g_{t}(1  h_{t}) > 0 whenever p(h_{t}) > 0.
One strength of MTPs is that if the support of the exposure variable is known, then they may be formulated to avoid violations of the positivity assumption. For example,Haneuse and Rotnitzky (2013) consider the effect of an MTP reducing surgery operating times on postsurgery complications. There are real world constraints–ethical considerations, surgery characteristics, limits of the surgeon–that imply a lower bound on surgery operation time and therefore would imply positivity violations for certain regimes. Taking these constraints into consideration, the authors were able to formulate interventions that were in the support of the data. Another example and further discussion is given in Rudolph et al. (2022). Furthermore, the MTP methods presented here can be used to solve positivity problems in other interesting ways, such as defining incremental propensity score interventions that shift the propensity score in a riskratio scale Wen et al. (2021); Díaz et al. (2022).
Under these identification assumptions, the causal parameter θ is identified as follows. Set m_{τ}_{+1} = Y . For t = τ, . . ., 1, recursively define
Then θ is identified as E[m_{1}(A^{d}_{1}, L_{1})]. To illustrate this recursive definition, consider an example when τ = 2. Then we have
[End Page 107]
which easily generalizes to many time points. In addition to the function m_{t}, the description of the estimation methods below will benefit from the definition of the ratio between the density of the postintervention exposure, A^{d}_{t}, and the density of the observed exposure, A_{t}. That is, define g^{d}_{t}(a_{t}  h_{t}) as the density of A^{d}_{t} conditional on H_{t} = h_{t} evaluated at h_{t}. Likewise define g_{t}(a_{t}  h_{t}) as the conditional density of A_{t}. Our estimation procedure will make use of the density ratio defined as:
This density ratio reduces to the wellknown inverse probability weights used in estimation of the ATE in the singletime point setting where Z = (W, A, Y ) and A is binary.
3. Estimating modified treatment policy effects
The literature on estimation of the effects of dynamic treatment regimes d_{t}(h_{t}) that do not depend on a_{t} is rich. Compared to this literature, estimation of the effect of modified treatment policies d(a_{t}, h_{t}) that also depend on a_{t} introduces certain technical complexities that make the generalization of the methods nontrivial. Specifically, the parameter functional identifying the average effect of a treatment regime d_{t}(h_{t}) does not depend on the treatment mechanism g_{t}(a_{t}  h_{t}). This means that for estimation purposes, the treatment mechanism can often be treated as a nuisance parameter, which leads to important and useful properties such as doublyrobust estimation, etc. When the treatment policies are d(a_{t}, h_{t}) and do depend on a_{t}, the functional identifying the counterfactual average does depend on g_{t}(a_{t}  h_{t}). This means that the treatment mechanism cannot be treated as a nuisance parameter, and it means that it is unclear whether doublyrobust estimators can be constructed. In this paper and in our companion technical paper Díaz et al. (2021), we have made an assumption called piecewise smooth invertibility on the treatment policy d_{t}(a_{t}, h_{t}), and have shown that this assumption is sufficient to allow the construction of doublyrobust estimators. The methods implemented in the package require that the analyst specifies treatment policies that satisfy this assumption. Examples of treatment policies d(a_{t}, h_{t}) where doublyrobust estimation is unachievable are given, e.g., in Díaz and van der Laan (2013); Kennedy (2019); Díaz et al. (2022).
3.1 Estimation methods
The lmtp package implements four estimation methods: a targeted minimumloss based estimator (TMLE), a sequential doublyrobust estimator (SDR), an estimator based on the parametric Gformula, and an inverse probability weighted (IPW) estimator. We will only describe the use of TMLE, lmtp_tmle, and SDR, lmtp_sdr, as their use is strongly encouraged over the others based on their advantageous theoretical properties which allow for machine learning while maintaining the ability to compute valid confidence intervals and pvalues.
Targeted minimumloss based estimation is a general framework for constructing asymptotically linear estimators leveraging machine learning, with an optimal biasvariance tradeoff for the target causal parameter (van der Laan and Rose, 2011, 2018). In general, TMLE is constructed from a factorization of observed data likelihood into an outcome regression [End Page 108] and an intervention mechanism. Using the outcome regression, an initial estimate of the target parameter is constructed and then debiased by a fluctuation that depends on a function of the intervention mechanism. The sequential doublyrobust estimator is based on a unbiased transformation of the efficient influence function of the target estimand. For a thorough discussion of TMLE and SDR for static, dynamic, and modified treatment policies, we refer the reader to van der Laan and Gruber (2011); Luedtke et al. (2017); Rotnitzky et al. (2017); Díaz et al. (2021).
TMLE and SDR require estimation of two parameters at each time point: an outcome mechanism and an intervention mechanism. Both TMLE and SDR are multiplyrobust in that they allow certain configurations of nuisance parameters to be inconsistently estimated. Specifically, TMLE is considered τ + 1multiply robust in that it allows for inconsistent estimation of all the intervention mechanisms prior to any time point t, as long as all outcome mechanisms after time t are consistently estimated. SDR is 2^{τ} robust in that at each time point, estimation of at most either the intervention mechanism or outcome mechanism is allowed to be inconsistent. Both TMLE and SDR are efficient when all the treatment mechanism and outcome regression are consistently estimated at a given consistency rate, but the SDR has better protection against model misspecification (see Luedtke et al., 2017; Rotnitzky et al., 2017; Díaz et al., 2021, for more details).
It is important to note that the SDR estimator can produce an estimate θˆ outside of the bounds of the parameter space (e.g., probability estimates outside [0, 1]), while the TMLE guarantees that the estimate is within bounds of the parameter space. With this in mind and because for a single timepoint TMLE and SDR are equally robust, we recommend use of TMLE for the case of a single timepoint, while we recommend use of SDR for the longitudinal setting.
3.2 Required data structure
Data is passed to lmtp estimators through the data argument. Data should be in wide format with one column per variable per time point under study (i.e., there should be one column for every variable in Z). These columns do not have to be in any specific order and the data set may contain variables that are not used in estimation. We refer the reader to Hoffman et al. (2022) for a complete example of preprocessing data for lmtp. The names of treatment variables, censoring variables, baseline covariates, and timevarying covariates are specified using the trt, cens, baseline, and time_vary arguments respectively. The trt, cens, and baseline arguments accept character vectors and the trt and cens arguments should be ordered according to the timeordering of the data generating mechanism. The time_vary argument accepts an unnamed list sorted according to the timeordering of the model with each index containing the name of the timevarying covariates for the given time. The outcome variable is specified through the outcome argument.
Estimators are compatible with continuous, dichotomous and survival outcomes. In the case of a dichotomous or continuous outcome, only a single variable name should be passed to the outcome argument. For survival outcomes, a vector containing the names of the intermediate outcome and final outcome variables, ordered according to time, should be specified with the outcome argument. Dichotomous and survival outcomes should be coded using zero’s and one’s where one indicates the occurrence of an event and zero otherwise. [End Page 109]
If working with a survival outcome, once an observation experiences an outcome, all future outcome variables should also be coded with a one. The outcome_type argument should be set to “continuous” for continuous outcomes, “binomial” for dichotomous outcomes, and “survival” for timetoevent outcomes. If the study is subject to losstofollowup, the cens argument must be provided. Censoring indicators should be coded using zero’s and one’s where one indicates an observation is observed at the next time and zero indicates losstofollowup. Once an observation’s censoring status is switched to zero it cannot change back to one. Missing data before an observation is losttofollowup is not allowed. While multiple solutions exist for the problem of missing data, we recommend the use of multiple imputation pending future work on incorporating missingness into the causal model. Assuming the missingness only depends on past variables, multiple imputation is an appropriate approach.
The k argument controls a Markov assumption on the data generating mechanism. When k = Inf, the history H_{t} will be constructed using all previous timepoint variables while setting k to any other value will restrict H_{t} to timevarying covariates from time t − k − 1 until t−1. Baseline confounders are always included in H_{t}. Making a Markov assumption is particularly helpful when the sample size is smaller than the number of variables in the the full covariate history. The create_node_list function may be used to inspect how variables will be used for estimation. It is specified with the same trt, baseline, time_vary, and k arguments as lmtp estimators and is used internally to create a “node list” that encodes [End Page 110] which variables should be used at each time point of estimation. For example, consider a study with the observed data structure
We can translate this data structure to R with
A list of lists is returned with the names of the variables in H_{t} to be used for estimation of the outcome regression and the treatment mechanism at every time t. Notice that variables A_{1} and A_{2} are included in the list of variables used for estimation of the treatment mechamism (trt). This is due to the fact that the nuisance parameter for the treatment mechanism is the density ratio r_{t}, which is a function of A_{1} and A_{2}.
The density ratio is estimated based on a classification trick using an auxiliary variable Λ as a pseudo outcome and the treatment as a predictor. We now briefly describe how this density ratio estimation is done; the process is fully automated and hidden from the software user. Specifically, the TMLE and SDR estimation methods require estimation of the ratio of the densities of A^{d}_{t} and A_{t}, conditional on the history H_{t}, defined as r_{t} above. This is achieved through computing the odds in a classification problem in an augmented dataset with 2n observations where the outcome is the auxiliary variable Λ (defined below) and the predictors are the variables A_{t} and H_{t}. In the 2n augmented data set, the data structure at time t is redefined as
where Λ_{λ}_{,}_{i} = λ_{i} indexes duplicate values. For all duplicated observations λ ∈ {0, 1} with the same i, H_{λ}_{,}_{i}_{,}_{t} is the same. For λ = 0, A_{λ}_{,}_{i}_{,}_{t} equals the observed exposure values A_{i}_{,}_{t}, whereas for λ = 1, A_{λ}_{,}_{i}_{,}_{t} equals the exposure values under the MTP d, namely A^{d}_{t}. The [End Page 111] classification approach to density ratio estimation proceeds by estimating the conditional probability that ∆ = 1 in this dataset, and dividing it by the corresponding estimate of the conditional probability that ∆ = 0. Specifically, denoting p^{λ} the distribution of the data in the augmented dataset, we have:
Further details on this algorithm may be found in our technical paper (Díaz et al., 2021).
3.3 Specifying Interventions
The treatment rule d_{t} can be specified using one of two arguments: shift or shifted. If using shift, the analyst provides a twoargument function where the first argument should correspond to a dataset containing the variables in Z and the second argument should accept a string corresponding to the variable A_{t} in Z. This function should return a size n vector of A_{t} modified according to d_{t}. We provide examples of these functions in §4 below. Alternatively, the analyst can supply a modified version of data to shifted where A_{t} are replaced with A^{d}_{t}.
3.4 Machine Learnning
An attractive property of multiplyrobust estimators is that they can incorporate flexible machinelearning algorithms for the estimation of nuisance parameters m_{t} and r_{t} while remaining consistent. The super learner algorithm is an ensemble learner than incorporates a set of candidate models through a weighted convexcombination based on crossvalidation (van der Laan et al., 2007). Asymptotically, this weighted combination of models, called the metalearner, will outperform any single one of its components.
Our package uses the implementation of the super learner provided by the SuperLearner package (Polley et al., 2021). The algorithms to be used in the super learner with lmtp_tmle and lmtp_sdr calls are specified with the lrnrs_trt and lrnrs_outcome arguments. The outcome variable type should guide users on selecting the appropriate candidate learners for use with the lrnrs_outcome argument. Regardless of whether an exposure is continuous, dichotomous, or categorical, the exposure mechanism is estimated using classification as discussed above, users should thus only include candidate learners capable of binary classification with the lrnrs_trt argument. If the data is highdimensional, the analyst should restrict the SuperLearner library to candidate algorithms suitable for highdimensional data.
Candidate learners that rely on crossvalidation for the tuning of hyperparameters should support grouped data if used with lrnrs_trt. Because estimation of the treatment mechanism relies on the augmented 2n duplicated data set, duplicated observations must be put into the same fold during samplesplitting. This is done automatically by the package.
3.5 Additional arguments
Samplesplitting and crossfitting is used with all methods to avoid certain technical conditions that may not hold for machine learning estimators (Zheng and van der Laan, 2011; [End Page 112] Chernozhukov et al., 2018). Specifically, the data are split in V volds and each nuisance parameter is estimated V times, excluding data in one fold each time. The number of folds V can be set with the folds argument. If data has a hierarchical structure, the id argument is used to indicate the name of a variable in the data set indicating unique groups. These identifiers will be used for generation of crossvalidation folds and will be accounted for in standard error calculations. If a continuous outcome has known bounds, these bounds may be specified using the bounds argument with a length two numeric vector where the first index is the lower bound and the second index is the upper bound. Parallel processing is supported using the future package Bengtsson (2021). The analyst can set a seed for reproducibility using the typical set.seed() approach (the L’Ecuyer CMRG algorithm, an algorithm that’s suitable for random number generation with parallel processing is used by future).
4. Examples
Example 1 Longitudinal MTP with no losstofollowup
We have simulated data on n = 5000 observations over a 5month period. Each observation has a continuous exposure (A_1, A_2, A_3, A_4) and covariate (L_1, L_2, L_3, L_4) recorded at months one through four and a dichotomous outcome (Y) at month five. We assume no losstofollowup and no Markov assumption. This data set is installed with the package and is stored in the object sim_t4.
For this example, we are interested in the effect of a longitudinal MTP where at each month an observation’s exposure decreases by one only if their observed exposure wouldn’t be less than one if modified. Our data structure has no baseline confounders and we will use only GLMs for estimation so the only objects we must specify are the treatment variables, the timevarying covariates, the outcome variable, and the MTP shift function. Because we are only using parametric models for estimation, crossfitting is not necessary and we may set folds = 1.
Example 2 Longitudinal MTP with right censoring
For this example, we have a simulated dataset of n = 1000 observations. Data was simulated for three time points with a continuous timevarying exposure at times t ∈ {1, 2} (A1, A2), a dichotomous timevarying covariate at times t ∈ {1, 2} (L1, L2), and a dichotomous outcome (Y) at time τ + 1 = 3. Losstofollowup is present after time t = 1 so the data set contains censoring indicators (C1, C2). This data is installed with the package and is stored in the object sim_cens.
Suppose we are interested in the additive effect of an MTP where exposure is increased by 0.5 at every time point for all observations. Instead of using a linear model, we will estimate the outcome regression and treatment mechanism using a super learner composed of a generalized linear model, a random forest (Wright and Ziegler, 2017), and multivariate adaptive regression splines (Milborrow, 2019).
If losstofollowup exists, we can estimate the population mean outcome under the observed exposures adjusting for informative losstofollowup by specifying shift = NULL.
Example 3 Survival analysis and deterministic effects
The lmtp package may also be used to estimate deterministic causal effects, such as the causal relative risk. Suppose we have timetoevent data on n = 2000 observations with a timeinvariant dichotomous exposure followed for a period of seven days. We wish to estimate the causal relative risk of experiencing the event by day seven. This data is installed with the package and is stored in the object sim_point_surv.
Example 4 Dynamic treatment regimes
Dynamic treatment regimes are treatment rules where treatment is applied based on a fixed rule that depends on covariate history. The lmtp package is capable of estimating the effects of deterministic dynamic treatment rules and modified treatment policies that depend on covariate history.
For this example we will use the same data from example 1. However, we will extend the longitudinal MTP used in that example so that a shift in the exposure also depends on covariate history and time. Specifically, if t = 1 exposure decreases by one if an observation’s realized exposure wouldn’t be less than one if modified. If t > 1 and the value of the timevarying covariate at that timepoint equals one then exposure should decrease by one if an observation’s realized exposure wouldn’t be less than one if modified; otherwise exposure shouldn’t be modified. We also make a Markov assumption and set k = 1, restricting the timevarying covariate history to only variables at t − 1.
Example 5 The shifted argument
In some scenarios, shift functions may be very complex and it may be simpler for the analyst to create a copy of the data with the treatment variables modified by d_{t} themselves. The modified copy of the data can then be supplied as an argument to the shifted function inlieu of a shift function to the shift argument. Only the variables in trt can be different between the data in data and shifted.
As an example, consider the same analysis from example 4. The example 4 shift function required: identifying the current timepoint using regular expressions, checking the value of a timevarying covariate, and defining a function inside of itself. Modifying the data according to the d_{t} prior to calling lmtp_sdr greatly simplifies the code.
5. Summary
The lmtp package provides a unified implementation of four estimators for general parameters in a causal inference setting with longitudinal data. The package is very flexible: it accommodates continuous, multivariate, and categorical exposures, continuous, binary, timetoevent outcomes, and rightcensoring. Two estimators in the package allow the use of machine learning model ensembles that mitigate the risk of model misspecification bias, providing confidence intervals and standard errors. We demonstrated the ease of use and flexibility of the package in four examples with data structures of varying complexity. The package provides further flexibility and allows the user to contrast the average of the counterfactual outcomes in any desirable scale: difference, risk ratio, or odds ratio. [End Page 118]
Mailman School of Public Health
Columbia University
New York, NY 10032, USA
Department of Population Health
New York University Grossman School of Medicine
New York, NY 10016, USA