Title
seqtte -- Sequential target trial emulation
Syntax
seqtte outcomevar [if] [in], id(varname) time(varname) treatment(varname) [options]
options Description ---------------------------------------------------------------------------------------------- Required id(varname) individual identifier variable time(varname) integer calendar time variable treatment(varname) binary treatment indicator (0/1) ---------------------------------------------------------------------------------------------- Optional covariates(varlist) adjustment covariates for the outcome model estimator(string) itt (default) or pp wdenominator(varlist) denominator weight model covariates; if omitted with pp, an unweighted per-protocol analysis is performed wnumerator(varlist) numerator weight model covariates; if supplied, stabilized weights are used truncation(#) truncation threshold for cumulative weights; default 25 selectionrandom randomly subsample the control-arm (id, trial) pairs selectionsample(#) proportion of control-arm pairs to retain under selectionrandom; default 0.5 bootstrap(#) number of bootstrap replicates for the standard error and 95% percentile CI; default 0 seed(#) random-number seed; default -1 (seed not set) plot plot cumulative incidence curves for each treatment arm survivalmax(#) cap the follow-up shown in the cumulative incidence curves; default is adaptive expandonly return the expanded sequential-trial dataset and skip the analysis ----------------------------------------------------------------------------------------------
Description
seqtte estimates the causal effect of a sustained treatment strategy using sequential target trial emulation (Hernán and Robins, 2016). Two estimators are available.
Intent-to-treat (ITT) (estimator(itt), the default). Estimates the effect of being assigned to treatment at trial entry, regardless of subsequent treatment changes. No weighting is used.
Per-protocol (PP) (estimator(pp)). Estimates the effect of sustained adherence to the assigned treatment strategy. Individuals are censored at the period in which they deviate from their assigned treatment. If wdenominator() is supplied, inverse probability of censoring weights (IPCW) are applied to adjust for informative censoring (weighted PP). If wdenominator() is omitted, censoring is applied but no weights are used (unweighted PP).
Input data. The data should be in long (person-period) format with one row per individual per time period. outcomevar should equal 1 only in the period the event first occurs and 0 otherwise. The time variable should take consecutive integer values.
Algorithm. For each eligible person-period (i.e. periods in which the individual has not yet received treatment), a trial is initiated. The individual is then followed from that trial entry time to the end of their observed follow-up. Data are expanded so that each person can contribute to multiple trials. A pooled logistic regression model is then fitted with quadratic polynomial terms for follow-up time within trial and trial number, with standard errors clustered by individual. The pooled logistic regression approximates a discrete-time hazard model, so the exponentiated treatment coefficient is reported as a hazard ratio.
Weight models (weighted PP only). When wdenominator() is supplied with estimator(pp), four logistic regression models are fitted on the pre-expansion data, stratified by prior treatment status (A_lag = 0 or 1): a denominator model including wdenominator() covariates and a cubic polynomial in calendar time, and (if wnumerator() is supplied) a numerator model including wnumerator() covariates and the same time polynomial. Unstabilized weights are used when wnumerator() is omitted; stabilized weights (numerator/denominator) are used when it is supplied. Cumulative products of per-period weights are formed within each trial, then truncated at truncation(). When wdenominator() is omitted, no weight models are fitted and the outcome model is fitted on the censored data without weighting.
Options
+----------+ ----+ Required +------------------------------------------------------------------------------
id(varname) specifies the variable identifying individuals.
time(varname) specifies the integer calendar time variable. Consecutive integer values represent consecutive periods.
treatment(varname) specifies the binary treatment indicator (0 = untreated, 1 = treated).
+----------+ ----+ Optional +------------------------------------------------------------------------------
covariates(varlist) specifies adjustment covariates for the outcome model. In the expanded dataset each person-trial record takes covariate values from the trial entry period, so both time-fixed and time-varying covariates are included at their trial-entry values. Factor-variable notation (e.g. i.group) is allowed.
estimator(string) specifies the estimator: itt (default) for the intent-to-treat effect, or pp for the per-protocol effect.
wdenominator(varlist) specifies the covariates for the denominator weight models. These are fitted on the pre-expansion data and should include all time-varying confounders of the treatment–outcome relationship. When supplied with estimator(pp), IPCW weights are applied (weighted PP); when omitted, censoring is applied without weighting (unweighted PP). Factor-variable notation (e.g. i.group) is allowed.
wnumerator(varlist) specifies the covariates for the numerator (stabilization) weight models. These should be a subset of wdenominator(), typically restricted to baseline (study-entry) values of covariates. When omitted, unstabilized weights are used. Factor-variable notation (e.g. i.group) is allowed.
truncation(#) specifies the upper truncation threshold applied to the cumulative IPW weights. Default is 25.
selectionrandom randomly subsamples the control-arm (id, trial) pairs (Bernoulli sampling), retaining all treated-arm pairs. This reduces the size of the expanded dataset and is useful when the expanded data are very large. By default no subsampling is applied (all pairs are used).
selectionsample(#) specifies the proportion of control-arm pairs to retain when selectionrandom is specified. Must be in (0, 1]; default is 0.5. It has no effect unless selectionrandom is also specified.
bootstrap(#) requests # bootstrap replicates (resampling individuals with replacement) to compute a bootstrap standard error and 95% percentile confidence interval for the log-hazard ratio, reported below the regression table. Default is 0 (no bootstrap).
seed(#) sets the random-number seed, for reproducibility of the bootstrap resampling and of selectionrandom. Default is -1, meaning the seed is not set.
plot produces cumulative incidence (1 - survival) curves for each treatment arm by g-computation from the fitted outcome model. Cannot be combined with expandonly.
survivalmax(#) caps the follow-up time shown in the cumulative incidence curves (with plot). The g-computation projects every trial over the follow-up grid, so at long follow-up times, where few trials contributed observed data, the curves are driven by model extrapolation. By default the follow-up is capped at the largest time where at least 10% of the baseline trials still contribute observed (uncensored) data (minimum 5 trials); specify survivalmax(#) to set the cap explicitly.
expandonly performs the data expansion only and leaves the expanded sequential-trial dataset in memory, skipping the weight models, outcome model, bootstrap, and cumulative-incidence steps. The returned data contain the original variables together with trial (calendar time of trial entry), followup (time since trial entry), period (calendar time, equal to trial + followup), and event (the period-specific outcome indicator); for estimator(pp) a censored indicator (and, for weighted PP, the cumulative weight) are also included. This option cannot be combined with bootstrap() or plot.
Examples
Setup: generate a synthetic person-period dataset with a baseline covariate (age) and a time-varying confounder (bmi)
. clear . set seed 42 . set obs 500 . gen id = _n . gen age = rnormal(50, 10) . expand 12 . bysort id (age): gen time = _n - 1 . gen bmi = . . bysort id (time): replace bmi = rnormal(25, 4) if time == 0 . bysort id (time): replace bmi = bmi[_n-1] + rnormal(0, 1) if time > 0 . gen treatment = . . bysort id (time): replace treatment = (runiform() < invlogit(-2 + 0.1 * (bmi - 25))) if time == 0 . bysort id (time): replace treatment = cond(treatment[_n-1] == 0, runiform() < invlogit(-2 + 0.1 * (bmi - 25)), runiform() < 0.7) if time > 0 . gen outcome = (runiform() < invlogit(-3 + 0.05 * (bmi - 25) - 0.4 * treatment)) . bysort id (time): gen cumev = sum(outcome) . drop if cumev > 1 . drop cumev
ITT estimator
. seqtte outcome, id(id) time(time) treatment(treatment) covariates(age)
Unweighted PP estimator (censoring applied, no IPCW weights)
. seqtte outcome, id(id) time(time) treatment(treatment) covariates(age) estimator(pp)
PP estimator with unstabilized weights
. seqtte outcome, id(id) time(time) treatment(treatment) covariates(age) estimator(pp) wdenominator(age bmi)
PP estimator with stabilized weights. Note that the numerator model should be a submodel of the denominator model; if the same covariates are given to both the models coincide, every weight is 1, and the fit is identical to the unweighted per-protocol analysis.
. seqtte outcome, id(id) time(time) treatment(treatment) covariates(age) estimator(pp) wdenominator(age bmi) wnumerator(age)
Cumulative incidence curves by treatment arm
. seqtte outcome, id(id) time(time) treatment(treatment) covariates(age) plot
Stored results
seqtte stores the following in e():
Scalars e(N) number of observations in the pooled logistic regression e(N_indiv) number of individuals in the original data e(N_orig) number of observations in the original data e(N_exp) number of observations in the expanded dataset e(r2_p) pseudo R-squared e(ll) log likelihood e(N_uniq_arm0) number of individuals contributing follow-up to arm 0 e(N_uniq_arm1) number of individuals contributing follow-up to arm 1 e(N_nonuniq_arm0) number of follow-up intervals in arm 0 e(N_nonuniq_arm1) number of follow-up intervals in arm 1 e(N_sel) observations after random selection (if selectionrandom) e(selection_sample) control-arm sampling proportion (if selectionrandom) e(N_boot) number of successful bootstrap replicates (if bootstrap()) e(bs_se) bootstrap standard error of the log-hazard ratio (if bootstrap()) e(bs_ll) lower bound of the bootstrap 95% percentile CI for the hazard ratio (if bootstrap()) e(bs_ul) upper bound of the bootstrap 95% percentile CI for the hazard ratio (if bootstrap())
Macros e(cmd) seqtte e(estimator) itt or pp e(depvar) name of the outcome variable e(clustvar) name of the cluster variable e(vcetype) Robust
Matrices e(b) coefficient vector (log scale; exponentiated values reported as hazard ratios) e(V) variance-covariance matrix e(bs_b) bootstrap log-hazard-ratio replicates (if bootstrap()) e(cif) cumulative incidence by arm and follow-up time (if plot)
References
Hernán MA, Robins JM. 2016. Using Big Data to Emulate a Target Trial When a Randomized Trial Is Not Available. American Journal of Epidemiology 183(8): 758–764.
Danaei G, Rodríguez LAG, Cantero OF, Logan R, Hernán MA. 2013. Observational data for comparative effectiveness research: an emulation of randomised trials of statins and primary prevention of coronary heart disease. Statistical Methods in Medical Research 22(1): 70–96.
Maringe C, Benitez Majano S, Exarchakou A, et al. 2020. Reflection on modern methods: trial emulation in the presence of immortal-time bias. Assessing the benefit of major surgery for elderly lung cancer patients using observational data. International Journal of Epidemiology 49(5): 1719–1729.
Author
Tom Palmer, University of Bristol, Bristol, UK. remlapmot@hotmail.com
Michalis Katsoulis, UCL, London, UK.
Please report any bugs or feature requests at https://github.com/remlapmot/seqtte/issues.