Parametric g-formula implementations for causal survival analyses

Racial and ethnic differences in the risk of dementia diagnosis under hypothetical blood pressure-lowering interventions: The Multi-Ethnic Study of Atherosclerosis

INTRODUCTION

Substantial racial and ethnic disparities in hypertension and dementia exist in the United States. We evaluated the effect of maintaining systolic blood pressure (SBP) below clinical thresholds on dementia incidence.

METHODS

We included 6806 Multi-Ethnic Study of Atherosclerosis participants (44 to 84 years old). We implemented the parametric g-formula to simulate the hypothetical interventions to reduce SBP below 120 and 140 mmHg over time, accounting for time-varying confounding. We estimated risk ratios (RRs) and risk differences for dementia incidence at 19 years.

RESULTS

The RRs (95% confidence intervals [CIs]) comparing an intervention reducing SBP below 120 mmHg to no intervention were 0.93 (0.87 to 0.99) for total sample, 0.95 (0.88 to 1.02) for White, 0.90 (0.79 to 1.02) for Black, 0.90 (0.78 to 1.05) for Latino, and 1.16 (0.83 to 1.55) for Chinese American participants. Results for lowering SBP below 140 mmHg and with death as competing event were attenuated.

DISCUSSION

The reduction of SBP below 120 mmHg over time has modest effects on reducing dementia incidence. More work is needed to understand the heterogeneity across racial and ethnic groups.

HIGHLIGHTS

There is a potential beneficial effect in lowering SBP to reduce the risk of dementia, which may vary by race and ethnicity. The percentage of participants who would need intervention on blood pressure to meet clinical thresholds is greater for Black and Latino communities. Results are sensitive to the way that death is specified in the research question and analysis.

Semiparametric g-computation for survival outcomes with time-fixed exposures: An illustration

PURPOSE

Generalized (g-) computation is a useful tool for causal inference in epidemiology. However, in settings when the outcome is a survival time subject to right censoring, the standard pooled logistic regression approach to g-computation requires arbitrary discretization of time, parametric modeling of the baseline hazard function, and the need to expand one's dataset. We illustrate a semiparametric Breslow estimator for g-computation with time-fixed treatments and survival outcomes that is not subject to these limitations.

METHODS

We compare performance of the Breslow g-computation estimator to the pooled logistic g-computation estimator in simulations and illustrate both approaches to estimate the effect of a 3-drug vs 2-drug antiretroviral therapy regimen among people with HIV.

RESULTS

In simulations, both approaches performed well at the end of follow-up. The pooled logistic approach was biased at times between the endpoints of the discrete time intervals used, while the Breslow approach was not. In the example, both approaches estimated a 1-year risk difference of about 6 % in favor of the 3-drug regimen, but the shape of the survival curves differed.

CONCLUSIONS

The Breslow g-computation estimator of counterfactual risk functions does not rely on strong parametric assumptions about the time-to-event distribution or onerous dataset expansions.

Evaluating Effects of Multilevel Interventions on Disparity in Health and Healthcare Decisions

In this paper, we introduce an analytic approach for assessing effects of multilevel interventions on disparity in health outcomes and health-related decision outcomes (i.e., a treatment decision made by a healthcare provider). We outline common challenges that are encountered in interventional health disparity research, including issues of effect scale and interpretation, choice of covariates for adjustment and its impact on effect magnitude, and the methodological challenges involved with studying decision-based outcomes. To address these challenges, we introduce total effects of interventions on disparity for the entire sample and the treated sample, and corresponding direct effects that are relevant for decision-based outcomes. We provide weighting and g-computation estimators in the presence of study attrition and sketch a simulation-based procedure for sample size determinations based on precision (e.g., confidence interval width). We validate our proposed methods through a brief simulation study and apply our approach to evaluate the RICH LIFE intervention, a multilevel healthcare intervention designed to reduce racial and ethnic disparities in hypertension control.

Bias due to coarsening of time intervals in the inference for the effectiveness of colorectal cancer screening

BACKGROUND

Observational studies are frequently used to estimate the comparative effectiveness of different colorectal cancer (CRC) screening methods due to the practical limitations and time needed to conduct large clinical trials. However, time-varying confounders, e.g. polyp detection in the last screening, can bias statistical results. Recently, generalized methods, or G-methods, have been used for the analysis of observational studies of CRC screening, given their ability to account for such time-varying confounders. Discretization, or the process of converting continuous functions into discrete counterparts, is required for G-methods when the treatment and outcomes are assessed at a continuous scale.

DEVELOPMENT

This paper evaluates the interplay between time-varying confounding and discretization, which can induce bias in assessing screening effectiveness. We investigate this bias in evaluating the effect of different CRC screening methods that differ from each other in typical screening frequency.

APPLICATION

First, using theory, we establish the direction of the bias. Then, we use simulations of hypothetical settings to study the bias magnitude for varying levels of discretization, frequency of screening and length of the study period. We develop a method to assess possible bias due to coarsening in simulated situations.

CONCLUSIONS

The proposed method can inform future studies of screening effectiveness, especially for CRC, by determining the choice of interval lengths where data are discretized to minimize bias due to coarsening while balancing computational costs.

Tree-based subgroup discovery using electronic health record data: heterogeneity of treatment effects for DTG-containing therapies

The rich longitudinal individual level data available from electronic health records (EHRs) can be used to examine treatment effect heterogeneity. However, estimating treatment effects using EHR data poses several challenges, including time-varying confounding, repeated and temporally non-aligned measurements of covariates, treatment assignments and outcomes, and loss-to-follow-up due to dropout. Here, we develop the subgroup discovery for longitudinal data algorithm, a tree-based algorithm for discovering subgroups with heterogeneous treatment effects using longitudinal data by combining the generalized interaction tree algorithm, a general data-driven method for subgroup discovery, with longitudinal targeted maximum likelihood estimation. We apply the algorithm to EHR data to discover subgroups of people living with human immunodeficiency virus who are at higher risk of weight gain when receiving dolutegravir (DTG)-containing antiretroviral therapies (ARTs) versus when receiving non-DTG-containing ARTs.

RKHS-based covariate balancing for survival causal effect estimation

Survival causal effect estimation based on right-censored data is of key interest in both survival analysis and causal inference. Propensity score weighting is one of the most popular methods in the literature. However, since it involves the inverse of propensity score estimates, its practical performance may be very unstable, especially when the covariate overlap is limited between treatment and control groups. To address this problem, a covariate balancing method is developed in this paper to estimate the counterfactual survival function. The proposed method is nonparametric and balances covariates in a reproducing kernel Hilbert space (RKHS) via weights that are counterparts of inverse propensity scores. The uniform rate of convergence for the proposed estimator is shown to be the same as that for the classical Kaplan-Meier estimator. The appealing practical performance of the proposed method is demonstrated by a simulation study as well as two real data applications to study the causal effect of smoking on survival time of stroke patients and that of endotoxin on survival time for female patients with lung cancer respectively.

Evaluating Model Specification When Using the Parametric G-Formula in the Presence of Censoring

The noniterative conditional expectation (NICE) parametric g-formula can be used to estimate the causal effect of sustained treatment strategies. In addition to identifiability conditions, the validity of the NICE parametric g-formula generally requires the correct specification of models for time-varying outcomes, treatments, and confounders at each follow-up time point. An informal approach for evaluating model specification is to compare the observed distributions of the outcome, treatments, and confounders with their parametric g-formula estimates under the "natural course." In the presence of loss to follow-up, however, the observed and natural-course risks can differ even if the identifiability conditions of the parametric g-formula hold and there is no model misspecification. Here, we describe 2 approaches for evaluating model specification when using the parametric g-formula in the presence of censoring: 1) comparing factual risks estimated by the g-formula with nonparametric Kaplan-Meier estimates and 2) comparing natural-course risks estimated by inverse probability weighting with those estimated by the g-formula. We also describe how to correctly compute natural-course estimates of time-varying covariate means when using a computationally efficient g-formula algorithm. We evaluate the proposed methods via simulation and implement them to estimate the effects of dietary interventions in 2 cohort studies.

Longitudinal plasmode algorithms to evaluate statistical methods in realistic scenarios: an illustration applied to occupational epidemiology

INTRODUCTION

Plasmode simulations are a type of simulations that use real data to determine the synthetic data-generating equations. Such simulations thus allow evaluating statistical methods under realistic conditions. As far as we know, no plasmode algorithm has been proposed for simulating longitudinal data. In this paper, we propose a longitudinal plasmode framework to generate realistic data with both a time-varying exposure and time-varying covariates. This work was motivated by the objective of comparing different methods for estimating the causal effect of a cumulative exposure to psychosocial stressors at work over time.

METHODS

We developed two longitudinal plasmode algorithms: a parametric and a nonparametric algorithms. Data from the PROspective Québec (PROQ) Study on Work and Health were used as an input to generate data with the proposed plasmode algorithms. We evaluated the performance of multiple estimators of the parameters of marginal structural models (MSMs): inverse probability of treatment weighting, g-computation and targeted maximum likelihood estimation. These estimators were also compared to standard regression approaches with either adjustment for baseline covariates only or with adjustment for both baseline and time-varying covariates.

RESULTS

Standard regression methods were susceptible to yield biased estimates with confidence intervals having coverage probability lower than their nominal level. The bias was much lower and coverage of confidence intervals was much closer to the nominal level when considering MSMs. Among MSM estimators, g-computation overall produced the best results relative to bias, root mean squared error and coverage of confidence intervals. No method produced unbiased estimates with adequate coverage for all parameters in the more realistic nonparametric plasmode simulation.

CONCLUSION

The proposed longitudinal plasmode algorithms can be important methodological tools for evaluating and comparing analytical methods in realistic simulation scenarios. To facilitate the use of these algorithms, we provide R functions on GitHub. We also recommend using MSMs when estimating the effect of cumulative exposure to psychosocial stressors at work.

Propensity score weighting methods for causal subgroup analysis with time-to-event outcomes

Evaluating causal effects of an intervention in pre-specified subgroups is a standard goal in comparative effectiveness research. Despite recent advancements in causal subgroup analysis, research on time-to-event outcomes has been lacking. This article investigates the propensity score weighting method for causal subgroup survival analysis. We introduce two causal estimands, the subgroup marginal hazard ratio and subgroup restricted average causal effect, and provide corresponding propensity score weighting estimators. We analytically established that the bias of subgroup-restricted average causal effect is determined by subgroup covariate balance. Using extensive simulations, we compare the performance of various combinations of propensity score models (logistic regression, random forests, least absolute shrinkage and selection operator, and generalized boosted models) and weighting schemes (inverse probability weighting, and overlap weighting) for estimating the causal estimands. We find that the logistic model with subgroup-covariate interactions selected by least absolute shrinkage and selection operator consistently outperforms other propensity score models. Also, overlap weighting generally outperforms inverse probability weighting in terms of balance, bias and variance, and the advantage is particularly pronounced in small subgroups and/or in the presence of poor overlap. We applied the methods to the observational Comparing Options for Management: PAtient-centered REsults for Uterine Fibroids study to evaluate the causal effects of myomectomy versus hysterectomy on the time to disease recurrence in a number of pre-specified subgroups of patients with uterine fibroids.

Intervention treatment distributions that depend on the observed treatment process and model double robustness in causal survival analysis

The generalized g-formula can be used to estimate the probability of survival under a sustained treatment strategy. When treatment strategies are deterministic, estimators derived from the so-called efficient influence function (EIF) for the g-formula will be doubly robust to model misspecification. In recent years, several practical applications have motivated estimation of the g-formula under non-deterministic treatment strategies where treatment assignment at each time point depends on the observed treatment process. In this case, EIF-based estimators may or may not be doubly robust. In this paper, we provide sufficient conditions to ensure the existence of doubly robust estimators for intervention treatment distributions that depend on the observed treatment process for point treatment interventions and give a class of intervention treatment distributions dependent on the observed treatment process that guarantee model doubly and multiply robust estimators in longitudinal settings. Motivated by an application to pre-exposure prophylaxis (PrEP) initiation studies, we propose a new treatment intervention dependent on the observed treatment process. We show there exist (1) estimators that are doubly and multiply robust to model misspecification and (2) estimators that when used with machine learning algorithms can attain fast convergence rates for our proposed intervention. Finally, we explore the finite sample performance of our estimators via simulation studies.

Comparative Effectiveness and Safety of Seizure Prophylaxis Among Adults After Acute Ischemic Stroke

BACKGROUND

Older adults occasionally receive seizure prophylaxis in an acute ischemic stroke (AIS) setting, despite safety concerns. There are no trial data available about the net impact of early seizure prophylaxis on post-AIS survival.

METHODS

Using a stroke registry (American Heart Association's Get With The Guidelines) individually linked to electronic health records, we examined the effect of initiating seizure prophylaxis (ie, epilepsy-specific antiseizure drugs) within 7 days of an AIS admission versus not initiating in patients ≥65 years admitted for a new, nonsevere AIS (National Institutes of Health Stroke Severity score ≤20) between 2014 and 2021 with no recorded use of epilepsy-specific antiseizure drugs in the previous 3 months. We addressed confounding by using inverse-probability weights. We performed standardization accounting for pertinent clinical and health care factors (eg, National Institutes of Health Stroke Severity scale, prescription counts, seizure-like events).

RESULTS

The study sample included 151 patients who received antiseizure drugs and 3020 who did not. The crude 30-day mortality risks were 219 deaths per 1000 patients among epilepsy-specific antiseizure drugs initiators and 120 deaths per 1000 among noninitiators. After standardization, the estimated mortality was 251 (95% CI, 190-307) deaths per 1000 among initiators and 120 (95% CI, 86-144) deaths per 1000 among noninitiators, corresponding to a risk difference of 131 (95% CI, 65-200) excess deaths per 1000 patients. In the prespecified subgroup analyses, the risk difference was 52 (95% CI, 11-72) among patients with minor AIS and 138 (95% CI, 52-222) among moderate-to-severe AIS patients. Similarly, the risk differences were 86 (95% CI, 18-118) and 157 (95% CI, 57-219) among patients aged 65 to 74 years and ≥75 years, respectively.

CONCLUSIONS

There was a higher risk of 30-day mortality associated with initiating versus not initiating seizure prophylaxis within 7 days post-AIS. This study does not support the role of seizure prophylaxis in reducing 30-day poststroke mortality.

A Simulation Study Comparing the Performance of Time-Varying Inverse Probability Weighting and G-Computation in Survival Analysis

Inverse probability weighting (IPW) and g-computation are commonly used in time-varying analyses. To inform decisions on which to use, we compared these methods using a plasmode simulation based on data from the Effects of Aspirin in Gestation and Reproduction (EAGeR) Trial (June 15, 2007-July 15, 2011). In our main analysis, we simulated a cohort study of 1,226 individuals followed for up to 10 weeks. The exposure was weekly exercise, and the outcome was time to pregnancy. We controlled for 6 confounding factors: 4 baseline confounders (race, ever smoking, age, and body mass index) and 2 time-varying confounders (compliance with assigned treatment and nausea). We sought to estimate the average causal risk difference by 10 weeks, using IPW and g-computation implemented using a Monte Carlo estimator and iterated conditional expectations (ICE). Across 500 simulations, we compared the bias, empirical standard error (ESE), average standard error, standard error ratio, and 95% confidence interval coverage of each approach. IPW (bias = 0.02; ESE = 0.04; coverage = 92.6%) and Monte Carlo g-computation (bias = -0.01; ESE = 0.03; coverage = 94.2%) performed similarly. ICE g-computation was the least biased but least precise estimator (bias = 0.01; ESE = 0.06; coverage = 93.4%). When choosing an estimator, one should consider factors like the research question, the prevalences of the exposure and outcome, and the number of time points being analyzed.

MULTIPLY ROBUST ESTIMATORS OF CAUSAL EFFECTS FOR SURVIVAL OUTCOMES

Multiply robust estimators of the longitudinal g-formula have recently been proposed to protect against model misspecification better than the standard augmented inverse probability weighted estimator (Rotnitzky et al., 2017; Luedtke et al., 2018). These multiply robust estimators ensure consistency if one of the models for the treatment process or outcome process is correctly specified at each time point. We study the multiply robust estimators of Rotnitzky et al. (2017) in the context of a survival outcome. Specifically, we compare various estimators of the g-formula for survival outcomes in order to 1) understand how the estimators may be related to one another, 2) understand each estimator's robustness to model misspecification, and 3) construct estimators that can be more efficient than others in certain model misspecification scenarios. We propose a modification of the multiply robust estimators to gain efficiency under misspecification of the outcome model by using calibrated propensity scores over non-calibrated propensity scores at each time point. Theoretical results are confirmed via simulation studies, and a practical comparison of these estimators is conducted through an application to the US Veterans Aging Cohort Study.

Revisiting the g-null Paradox

The (noniterative conditional expectation) parametric g-formula is an approach to estimating causal effects of sustained treatment strategies from observational data. An often-cited limitation of the parametric g-formula is the g-null paradox: a phenomenon in which model misspecification in the parametric g-formula is guaranteed in some settings consistent with the conditions that motivate its use (i.e., when identifiability conditions hold and measured time-varying confounders are affected by past treatment). Many users of the parametric g-formula acknowledge the g-null paradox as a limitation when reporting results but still require clarity on its meaning and implications. Here, we revisit the g-null paradox to clarify its role in causal inference studies. In doing so, we present analytic examples and a simulation-based illustration of the bias of parametric g-formula estimates under the conditions associated with this paradox. Our results highlight the importance of avoiding overly parsimonious models for the components of the g-formula when using this method.

G-computation and doubly robust standardisation for continuous-time data: A comparison with inverse probability weighting

In time-to-event settings, g-computation and doubly robust estimators are based on discrete-time data. However, many biological processes are evolving continuously over time. In this paper, we extend the g-computation and the doubly robust standardisation procedures to a continuous-time context. We compare their performance to the well-known inverse-probability-weighting estimator for the estimation of the hazard ratio and restricted mean survival times difference, using a simulation study. Under a correct model specification, all methods are unbiased, but g-computation and the doubly robust standardisation are more efficient than inverse-probability-weighting. We also analyse two real-world datasets to illustrate the practical implementation of these approaches. We have updated the R package RISCA to facilitate the use of these methods and their dissemination.