Semiparametric Estimation of the Impacts of Longitudinal Interventions on Adolescent Obesity using Targeted Maximum-Likelihood: Accessible Estimation with the ltmle Package

Data structure

For this analysis we used a subset of the NGHS data consisting of the participants recruited by the University of California at Berkeley center, one of three recruitment sites for NGHS (n = 530). The participants were 9–10 years of age at study entry and followed for 10 years; 887 girls enrolled at baseline. The study collected anthropometric measurements annually and an extensive set of variables potentially relevant to weight gain, including physical, behavioral, socioeconomic and mental health factors such as pubertal maturation stage, diet, physical activity, parental education, and perceived self-worth [10, 15–20].

The time-dependent outcome variable of interest was BMI-for-age z-score, which indicates BMI relative to other girls of the same age on a standard deviation scale [21]. We focused on physical activity and total calories as the exposures/potential interventions of interest, with the hypotheses that increasing physical activity would result in a lower 10-year BMI and that increasing total calories would have the opposite effect. Physical activity was measured using a Habitual Activity Questionnaire that was adapted from a questionnaire developed by Ku et al. (1981) [22], and compared against two other assessment methods [18]. Total calories were estimated from 3-day food diaries [7]. Potential confounding variables were selected based primarily upon previously reported associations with BMI [10, 16–24]. These included baseline race (white or African American), and the time-dependent variables pubertal maturation stage (four levels: prepuberty, early maturity, midpuberty, maturity), number of hours of television watched per week, perceived stress scale [23], global self-worth score (an indicator of self-esteem measured using a Harter’s Self Perception Profile for Children) [18, 24], as well as the outcome, BMI, at previous measurements.

We focused on three time points after enrollment: years 0 (age 9–10 years), 5 and 10 (19–20 years); selection of these time points was based on preliminary work indicating that these were the most relevant for capturing BMI trajectories. We restricted the sample to the subset of participants who had total calories, physical activity, and BMI z-score measured at years 0, 5, and 10 (n=530) since this is the level of resolution such that there would be relatively little missing data. However, even at this resolution, there were still missing data for some of the time-varying covariates (confounders), and we imputed missing values by using a local average of the years around the time point of interest within the subject, rather than omitting the girls with missing values completely. This did not appreciably change the overall summary statistics for the covariates when compared to a complete case analysis.

Parameter of interest

The research question we endeavored to answer concerned the potential causal effect of intervening to set longitudinal profiles of physical activity level or caloric intake amount on BMI z-score at year 10 of the study, when the girls were aged 19–20 years. An ideal experiment that would answer this question would be to randomize a cohort of girls at baseline to each potential longitudinal pattern of physical activity and diet, then follow up with the girls, ensuring perfect adherence and no attrition. The causal inference estimation framework makes transparent the identifiability assumptions necessary to estimate such parameters from observational data. The structural causal model (SCM) reflecting our belief about the time-ordering and relationships between the exposure, covariates, and outcome of interest was

where the subscripts denote the three time points of interest (years 0, 5 and 10), A_j, denotes the exposure of interest (physical activity or total calories) within time interval j, L_j are the vector of time-varying confounders (variables that affect both future exposure and the outcome of interest, including intervening measurements of BMI) within interval j, Ydenotes the outcome process of interest (BMI at year 10), and the Us denote unmeasured, independent (exogenous) variables, so that the variables that make up our data set are deterministic (but unknown) functions of the measured history, and some unmeasured error term. Now, we can represent the data as an independent, and identically distributed sample of O=(Aˉ,Lˉ,Y), where the over-bar notation denotes the history of the variable, e.g. Aˉ=(A0,A5,A10), and Aˉ5=(A0,A5).

Intervening on the exposure of interest at a certain time point corresponds to deterministically setting A0,A5,A10=a0,a5,a10, resulting in a modified set of structural equations. The counterfactual outcome is denoted Yaˉ and is interpreted as the value thatY would have taken under universal application in the population of the hypothetical intervention Aˉ=aˉ. The modified SCM, based on intervening by setting nodes to fixed values, is

where the (1) notation represents counterfactuals indexed by specific treatment regimes, a. Once the SCM is specified, the next step is to specify the counterfactuals indexed by interventions on A. The interventions simulated in this analysis were ones that deterministically set values of the vector aˉ=a0,a5,a10. The variables intervened on were physical activity level (1 denotes a high level defined as > 20 METS-times/wk, 0 denotes a low level defined as ≤ 20 METS-times/wk), caloric intake (1 denotes a total caloric intake of 2,000 kcal/day or less and 0 denotes a total caloric intake more than 2,000 kcal). The cutoff for the total calories intervention was determined based on average energy requirement recommendations for women [25]. The physical activity cutoff was determined based on the distribution of physical activity during the first five study years, before physical activity began to decline steeply [26]. The counterfactuals of interest were therefore denoted as Y(a₁,a₅,a₁₀), which can be interpreted as the BMI z-score that participants would have had under interventions on physical activity and total calories. A value a_j = 1 means that a person can be considered “treated” at time point j, while a value a_j = 0 means that they were “untreated” at that timepoint j. The intervention aˉ = {1, 1, 1} was chosen because it represents good physical activity behavior and/or eating habits being instilled at a young age and continued through adolescence. In contrast, the early intervention of aˉ = {1, 0, 0} was chosen to represent a decline in physical activity and/or eating habits with age. The late intervention aˉ = {0, 0, 1} represents the opposite.

Identifiability assumptions – G-computation

In order to estimate the marginal distribution of different counterfactuals from observed data, identifiability assumptions must be asserted. In longitudinal data, one such assumption, the sequential randomization assumption [14, 27] corresponds to

, where Parents indicates all preceding measured variables, and in this case Y(aˉ) are the set of counterfactuals defined by the combinations of possible interventions. Under this assumption (as well as others outlined below), we can write the counterfactual as an estimand that is a function only of the observed data-generating distribution, P₀:

This result is the longitudinal G-computation formula, which allows specification of parameters of the observed data distribution. The parameters we specified depended on the different intervention patterns contained in a¯. The parameters of interest are functions of the marginal mean under longitudinal treatment regimes, E[Y(a₀, a₅, a₁₀)], or:

where ψ1 can be interpreted as the mean difference of BMI z-scores within a population of universal, early sustained intervention versus the mean BMI z-score in the population without any intervention; ψ2 is a comparison of early sustained intervention to late intervention; ψ3, compares early intervention, that is not sustained, versus late intervention; ψ4 is equivalent to ψ1, but now comparing early, but not sustained intervention to a population without any intervention. Finally, ψ5 compares sustained intervention to a population where everyone follows relatively high caloric intake or relatively low physical activity.

In order for one to be able to estimate these causal parameters, there needs to be sufficient natural experimentation of the variable of interest (the A_j) within groups defined by the history of covariates, Lˉj and previous levels of the A_j–. This is the so-called positivity assumption. There cannot be combinations of covariate and exposure histories, for which all participants are only “treated” or “untreated”, since this would make the right-hand side of eq. (2) undefined. Positivity is especially problematic when the covariates and treatment are continuous and in these data, there were many covariates. Thus we discretized the covariates by creating groupings defined by reasonable cutoffs. Specifically, we dichotomized the number of hours of television according to the recommendations from the American Academy of Pediatrics as “high” if the individual watched more than 2 hours/wk and “low” if they watched two or fewer hours per week. Stress was dichotomized as “high” if the stress score was above 25 and “low” if the score was equal to or less than 25. Race, pubertal stage, and self-worth score did not require discretization to address the positivity requirement. Energy intake was used as a confounder in the analysis of an intervention on physical activity and vice versa, using the cutoffs defined above.

The substitution estimator is then given by

This estimator requires the estimates of the relevant regressions and joint distribution of intermediates, given the past. Typically, one derives the mean via time-sequential simulation of the L-process, setting the treatment history to the desired rule. Thus, the standard formulation of this substitution estimator requires often challenging estimation, since it requires estimation of joint densities.

Van der Laan and Gruber [28] present a different representation of eq. (3) that avoids having to estimate the joint conditional density of intermediate variables, e.g. the PLt=lt|Aˉt−1=aˉt−1,Lˉt−1=lˉt−1. Specifically, by invoking the tower rule, one can represent E[Yaˉ] as an iterative conditional expectation (ICE). The approach is based on the G-computation representation of the density of the covariate process, L, which can be written as:

, where Qjalˉj=Plj|L¯j−1=lˉj−1,A¯j−1=aˉj−1 is shorthand notation for defining the conditional distribution of a covariate process given a treatment process at a particular time j; equivalently define Q¯ja to be the associated conditional expectation implied by Qja. In addition, introduce notation for a set of random variables, La=L0,La1,…,LaJ,Y=LaJ+1 which has distribution, Pa, that is counterfactual random vectors of the covariate process generated under a specific treatment regime, a. The ICE representation is based on an equality which is formed from ICEs relative to these counterfactual distributions, or:

The substitution estimator, following eq. (6), starts with estimating the innermost conditional expectation, E[QˉJ+1a(Y)|LaJ)], and then moves outward, estimating the relevant regressions to generate the estimated conditional expectations, until at the end, one simply gets an average across all observations. The obvious computational virtue is that unlike the estimator (5), one does not have to estimate conditional densities of the covariate process, L, but only a set of conditional expectations (regressions) that can be done much more straightforwardly, using data-adaptive (semiparametric) techniques. This is an enormous advantage over the estimator based on eq. (5). The one small disadvantage, is these set of regressions must be run for each desired treatment rule, a, whereas using the original formulation, one only has to estimate the conditional densities once. However, this is a very small price to pay when the covariate process is high dimensional, i.e. consists of many variables, so any practical joint conditional density estimation requires either very strong assumptions (conditional independence, or joint normality, etc.) or a very large degree of smoothing (e.g. histogram density with very large bins). Table 1 contains the steps for defining the parameter of interest, as well as the associated estimation steps, for our specific obesity study.

Because the Qˉn,L,taAˉt−,Lˉ10 are defined via particular treatment regime of interest, a, this procedure is repeated for each particular regime of interest. This formulation offers a very compelling alternative to a substitution based on the standard representation (4). However, using the estimator based on the ICE approach still requires estimates of very high-dimensional regression, outcomes versus often many variables, potentially measured at several times in the past), in a very big, semi-parametric model. We discuss below how to target the estimation towards our particular parameter of interest, but before we do so, we discuss a general data-adaptive approach for deriving the regression models required for the estimate listed in Table 1.

Targeted maximum likelihood estimation

The data-adaptive fitting of the models that make up Step 1 above are targeted to optimizing the ability of models to predict the outcome, and not towards estimating the parameter. One can improve the estimator by “targeting” these fits towards the parameter of interest. Specifically, we used Targeted Maximum-Likelihood Estimation (TMLE), a two-stage estimator that augments the estimated models that comprise the ICE-inspired algorithm discussed in Table 1. This results in a bias-reduction step, which will remove residual bias relative to the substitution estimator if the treatment mechanism can be estimated consistently [14]. The TMLE augmentation [28] involves adding a so-called clever covariates to each of these regressions, which requires estimation of the treatment mechanism. Specifically, for each time point, it requires estimates of the probability of being in the treatment group (e.g. high physical activity) given the past (that is, all past covariates, treatments, and outcomes) [36].

Specifically, we use the same formulation as discussed above, but now based on Qn,L,J∗a, which is obtained, as a regression, treating the initial estimator, Qn,L,Ja as an offset, of Y against covariate:

where gk≡P(Ak=ak|Parents(Ak)). To derive the clever covariate one needs an estimate of g_k (g_n,k) and in this case, we used simple main terms logistic regression, both for the outcome of current treatment given the past. This is done as a compromise to ameliorate the residual bias that could result from the original SL fits, which are designed to minimize the risk of the prediction but avoid adding problems with creating large outliers via these clever covariates that can result from models that estimate probabilities of censoring or treatment either close to 1 or 0 (see chapter on CTMLE in Van der Laan et al. for a less ad hoc approach to model selection for the treatment/censoring models in the context of TMLE).

As discussed above, one of the virtues of TMLE over the initial substitution estimator is that it reduces bias due to the fact that the original statistical models are chosen based on minimizing expected loss with regards to prediction of the outcomes (intermediate and final), and not with regards to minimizing the mean-squared-error of the estimate of the parameter of interest. Relatedly, the TMLE estimate is consistent if either the original models or the treatment/censoring models are consistently estimated (doubly-robust). If they both are consistently estimated, the estimator is semi-parametrically (locally) efficient. Finally, the TMLE can be thought of a smoothing of the original substitution estimator, and thus can have more predictable sampling distributional properties; it is an asymptotically linear estimator with a known influence function, and this can be used to derive robust asymptotic inference [14]. Thus, one can derive an approximate SE by simply taking the sample standard deviation of the estimated IC of a subject divided by the sample size, or SEΨˆ=varˆ(IC)n (see van der Laan and Gruber [28] for the form of the IC).

The required calculations would seem a daunting task to code, and thus the effort required could inhibit taking such an approach. However, an R package, ltmle, provides a one-stop, user-friendly implementation of the ICE, TMLE estimator. One can either use the defaults available for estimating the prediction and treatment models, or one can specify particular algorithms, including SL with an associated library of learners. Though not discussed here, the package incorporates the possibility of missing data and not just fixed treatments of interest, but treatment rules. The package also returns, in addition to the TMLE estimate, the ICE based on non-augmented model, the so-called inverse probability of treatment weighted (IPTW) estimators, as well as a “naïve” estimate based on no dependent confounding. Finally, the package will also return estimates of marginal structural models if one wants to model the treatment-specific mean as a potentially smooth function of the history of intervention (e.g. the total number of time intervals with treatment). Thus, the powerful ltmle packages makes estimating a relatively complex suite of estimators amazingly straightforward, which opens the door to estimation of targeted casual parameters based on potentially complex interventions to a much wider audience.