---
title: "Estimate the heterogeneous treatment effect"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Estimate the heterogeneous treatment effect}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  markdown: 
    wrap: 80
---

```{r, include = FALSE}
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.

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

```{r}
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"
)
```

```{r}
pd_data <- DataStandard(BiSample, map)
head(pd_data)
```

```{r}
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.

```{r}
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:

```{r}
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.

```{r, fig.alt = "Time-specific treatment-effect model coefficients and demonstration bootstrap confidence intervals."}
separate$summary
separate$forest_plot
```

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

```{r}
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.

```{r}
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()`.

```{r}
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.

```{r, fig.alt = "Pooled treatment-effect model point estimates; bootstrap intervals are not calculated in this example."}
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`.

```{r}
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.