Estimate the heterogeneous treatment effect

In PDRobust: Robust Longitudinal Effects Under Truncation by Death

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  dev = "ragg_png",
  dpi = 192,
  fig.width = 7,
  fig.height = 4.5,
  out.width = "90%",
  fig.align = "center",
  warning = TRUE,
  message = TRUE
)

Introduction

The complete workflow is illustrated as follows. This article focuses on the package's main analysis functions, including HTESepT() and HTEAllT(). Those are the most significant functions in this package.

data("BiSample") -> Mapping() -> DataCheck() -> DataStandard()
                 -> prediction / diagnostic / analysis functions

library(PDRobust)
data("BiSample", package = "PDRobust")
map <- Mapping(
  id = "id",
  time = "time",
  treatment = "A",
  survival = "S",
  outcome = "Y",
  baseline_time = 0,
  cutoff_time = 2,
  covariates = c("X1", "X2", "X3", "X4", "X5","X6"),
  interest_vars = c("X1", "X5"),
  y_type = "B"
)

pd_data <- DataStandard(BiSample, map)
head(pd_data)

ps_fo <- A ~ X1 + X3 + X4 + X5 + X6
prin_fo <- S ~ (X1 + X3 + X4 + X5 + X6 ) * A
out_fo <- Y ~ (X1 + X3 + X4 + X5 + X6) * A 

Time-varying heterogeneous treatment effects

HTESepT() estimate time-specific heterogeneous treatment effect at one or more specified time points conditional on variables of interest defined in interest_vars when mapping. It returns the point estimate and, when requested, subject-level bootstrap standard errors and confidence intervals.

The argument data specifies the standardized dataset used for analysis. The arguments ps_fo, prin_fo and out_fo specify model formulas for propensity score model, principal score model and conditional outcome model, respectively. These models are refitted internally for the original sample and for every bootstrap sample. The argument target_time specifies the time points for estimation, and it must be a numeric vector such as c(1, 2). The mapped baseline time and cutoff time can also be included. Although the results are reported only at the requested time points, the principal scores are accumulated over all times from baseline time to cutoff time points.

The argument B, conf_level, max_attempts and verbose control the bootstrap process. B specifies the number of successful subject-level bootstrap replications. If B = 0 , no bootstrap is performed; the function only returns point estimates, while bootstrap standard errors and confidence interval are reported as NA. conf_level is the confidence level for the Wald confidence interval and defaults to 0.95. The max_attempts argument specified the maximum number of bootstrap samples that are attempted to obtain B successful replications. When max_attempts = NULL, it defaults to B*10. This allows additional attempts when a resampled dataset can not product valid estimate, for example because of inadequate variation, model-fitting failure, or nonconvergence. The argument verbose is a logical value and determine whether to the bootstrap progress messages are displayed.

The argument for mapping is not required because the information is carried inside the attributes of standardized dataset and used automatically.

The five bootstrap replications below keep the example fast; they are insufficient for substantive standard errors or confidence intervals. Increase B and assess inference stability for an analysis. The seed fixes the random resampling for a given R and dependency environment.

set.seed(20260912)
separate <- HTESepT(
  data = pd_data,
  ps_fo = ps_fo,
  prin_fo = prin_fo,
  out_fo = out_fo,
  target_time = c(1, 2),
  B = 5,
  conf_level = 0.95,
  max_attempts = NULL,
  verbose = TRUE
)

names(separate)

If verbose = TRUE, these messages report the number of successful replications relative to the requested value of B, together with the total number of attempts made.

The returned object contains the following components:

names(separate)

The primary outputs are summary and forest_plot. The coefficients parameterize the working treatment-effect model at each requested time.

The intercept represents the reference component of the conditional treatment-effect model. The remaining coefficients describe how the treatment effect varies with the corresponding baseline effect modifiers.

For continuous outcomes, the working effect is the intercept plus the linear combination of baseline effect modifiers, on the outcome scale. For binary outcomes, if eta is this linear predictor, the working risk difference is 2 * plogis(eta) - 1. Binary-model coefficients are therefore on this link scale; they are not odds ratios or direct risk differences. The displayed intervals describe the coefficients.

separate$summary
separate$forest_plot

Supplementary components include bootstrap_info, which summarizes the requested and successful replications, total attempts, completion status, and failures.

names(separate$bootstrap_info)

separate$bootstrap_info

boot_mat, which contains the coefficient estimates from each successful bootstrap replication; and convergence, which provides information about the estimating-equation solver.

separate$boot_mat

Pooled heterogeneous treatment effects

HTEAllT() estimates pooled heterogeneous treatment effects using every observed standardized analysis time from the mapped baseline through the cutoff. Unlike HTESepT(), it does not accept a target_time argument.

The arguments data, ps_fo, prin_fo, out_fo, B, conf_level, max_attempts, and verbose have the same interpretations as in HTESepT().

pooled <- HTEAllT(
  data = pd_data,
  ps_fo = ps_fo,
  prin_fo = prin_fo,
  out_fo = out_fo,
  B = 0,
  conf_level = 0.95,
  max_attempts = NULL,
  verbose = FALSE
)
names(pooled)

The primary outputs are again summary and forest_plot. The summary component reports the pooled HTE-model coefficients, including the intercept, the mapped effect modifiers, and a Time Effect term when the dataset contains at least two analysis times. The Time Effect describes the linear change in the conditional treatment-effect function per one-unit increase in standardized time.

pooled$summary
pooled$forest_plot

Additional components include analysis_times, which identifies the time points included in the pooled analysis, and time_effect_estimable, which indicates whether the time effect could be estimated. When only one analysis time is available, the time-effect term is omitted, time_effect_estimable is FALSE, and an explanatory message is stored in note.

pooled$time_effect_estimable
pooled$analysis_times

Bootstrap results, convergence information, formulas, settings, and the mapping used for the analysis are also retained in the returned object.

PDRobust documentation built on Oct. 2, 2026, 5:09 p.m.