Estimate Survival Difference Using LBC-Net and IPW Nelson–Aalen

`lbc_net_surv()` extends lbc_net to time-to-event outcomes. It first estimates propensity scores via LBC-Net and then combines them with inverse probability weights (IPW) to construct IPW-weighted Nelson–Aalen estimators for the marginal survival functions under treatment and control. The primary estimand is the survival difference \(S_1(t) - S_0(t)\) at one or more time points, where the survival estimators are based on the IPW-weighted Nelson–Aalen approach of Deng and Wang (2025).

Usage

lbc_net_surv(
  data = NULL,
  formula = NULL,
  Z = NULL,
  Tr = NULL,
  time = NULL,
  delta = NULL,
  K = 99,
  rho = 0.15,
  na.action = na.fail,
  gpu = 0,
  show_progress = TRUE,
  t_star = NULL,
  t_grid = NULL,
  ...,
  setup_lbcnet_args = list()
)

Arguments

data

an optional data frame, list, or environment (or an object coercible by `as.data.frame` to a data frame) containing the variables specified in `formula`. If a variable is not found in `data`, it is taken from `environment(formula)`, typically the environment from which `lbc_net()` is called. If `formula` is not provided, `Z` (a matrix of covariates) and `Tr` (a numeric treatment assignment vector) must be explicitly supplied.

formula

An object of class `"formula"` (or one that can be coerced to that class), specifying a symbolic description of the model in the form `Tr ~ X1 + X2 + ...`. If `formula` is provided, `Z` and `Tr` are extracted from `data`. If omitted, `Z` and `Tr` must be supplied explicitly.

Z

A numeric matrix or data frame of covariates. Required if `formula` is not provided. Each row represents an observation, and each column represents a covariate. If `formula` is used, `Z` is extracted from `data` automatically.

Tr

A numeric vector representing treatment assignment (typically 0/1). Required if `formula` is not provided. Must have the same number of rows as `Z`. If `formula` is used, `Tr` is extracted from `data` automatically.

time

Event or censoring time. This can be:

• a numeric vector of length equal to the sample size, or

• a character string naming the column in `data` that contains the event/censoring times.

delta

Event indicator (1 = event, 0 = censored). This can be:

• a numeric or integer vector of length equal to the sample size, or

• a character string naming the column in `data` that contains the event indicator.

K

an integer specifying the number of grid center points used to compute kernel weights in the local neighborhood. These weights are used to assess balance and calibration conditions. If specified, `ck` is automatically computed as `k / (K + 1)`, where `k = 1, ..., K`. The default is `99`, generating grid points from `0.01` to `0.99`. See **Details** for more information on kernel weights.

rho

a numeric value specifying the span used to determine the adaptive bandwidth `h` when `h` is not provided. The span controls the proportion of data included in the local neighborhood, ensuring a sufficient sample size for accurate training. The choice of `rho` influences local balance assessment and should be selected based on the data structure. While cross-validation can be used to approximate an optimal span, user discretion is advised. The default is `0.15`.

na.action

A function to specify the action to be taken if NAs are found. The default action is for the procedure to fail. An alternative is `na.omit`, which leads to rejection of cases with missing values on any required variable. (NOTE: If given, this argument must be named.)

gpu

An integer specifying the GPU device ID for computation if CUDA is available. If set to `0`, the function will attempt to use the default GPU. If CUDA is unavailable, computations will automatically fall back to the CPU.

show_progress

A logical value indicating whether to display a progress bar during training. If `TRUE` (default), displays elapsed time, remaining time, loss values, and training speed (iterations per second). This helps monitor training progress efficiently. Set to `FALSE` to disable the display.

t_star

Optional numeric value giving a single time point \(t^*\) at which to highlight the survival difference \(S_1(t^*) - S_0(t^*)\) in the output. If `t_star = NULL` (the default), the function sets \(t^*\) to the median survival time of the control group based on the Kaplan–Meier estimator; if that median is undefined (e.g., no events in control), the pooled-sample median is used instead.

t_grid

Optional numeric vector of time points at which to evaluate the survival curves and survival difference. If `t_grid` is `NULL` and `t_star` is specified or imputed (default), the function uses a single-point grid equal to `t_star`. If `t_grid` is non-`NULL`, it overrides `t_star` for evaluation purposes while `t_star` is still recorded in the output object.

...

Additional parameters for model tuning, including:

`ck`

a numeric vector of kernel center points. Values should be strictly between `0` and `1`. If `NULL`, `ck` is automatically computed using the default `K`. If provided, user-defined grid points should adhere to the constraint `0 < ck < 1`.

`h`

A numeric vector of bandwidth values for kernel weighting. By default, an adaptive bandwidth is automatically computed via preliminary probabilities estimated using logistic regression (`glm`), given `ck` and the span parameter `rho`. If `NULL`, `h` is computed using span_bw.

`kernel`

A character string specifying the kernel function used for local weighting. The default is `"gaussian"`. Available options include:

• `"gaussian"` (default): Ensures smooth weighting and continuity.

• `"uniform"`: Assigns equal weight to all observations within a fixed bandwidth.

• `"epanechnikov"`: Gives higher weight to observations closer to the kernel center.

See **Details** for more information on kernel weighting and its role in local balance estimation.

`seed`

An integer specifying the random seed for reproducibility in training the neural network. This seed is applied to PyTorch using `torch.manual_seed(seed)`, ensuring consistent results across runs when using stochastic optimization methods.

`hidden_dim`

An integer specifying the number of hidden units in the LBC-Net model. A rule of thumb is to set this as two to three times the number of covariates, but it should be significantly smaller than the sample size to prevent overfitting. Default is `100`.

`num_hidden_layers`

An integer specifying the number of hidden layers in the LBC-Net model. Default is `1`, which results in a three-layer network overall (input layer, one hidden layer, and output layer).

`vae_epochs`

An integer specifying the number of training epochs for the Variational Autoencoder (VAE). This determines how long the VAE component of the model is trained before being used in LBC-Net. Default is `250`.

`vae_lr`

A numeric value specifying the learning rate for training the VAE using the Adam optimizer. Controls how quickly the model updates its parameters during training. Default is `0.01`.

`max_epochs`

An integer specifying the maximum number of training epochs for LBC-Net. Early stopping is applied based on `lsd_threshold` to prevent unnecessary training. Default is `5000`.

`lr`

A numeric value specifying the initial learning rate for LBC-Net training using the Adam optimizer. The learning rate controls how much the model updates during each training step. Default is `0.05`.

`weight_decay`

A numeric value specifying the regularization parameter in the Adam optimizer for LBC-Net. Helps prevent overfitting by penalizing large weights in the model. Default is `1e-5`.

`balance_lambda`

A numeric value controlling the relative contributions of the local balance loss (\(Q_1\)) and calibration loss (\(Q_2\)) in the objective function, where the total loss is defined as \(Q = Q_1 + \lambda Q_2\). Default is `1.0`.

`epsilon`

A small numeric value controlling the lower and upper bounds of the estimated propensity scores. The default is `0.001`, ensuring scores remain within \([\epsilon, 1 - \epsilon]\) for numerical stability, particularly in cases of poor overlap. Setting `epsilon = 0` reverts to the standard logit link function. See **Details** for more on its role in model stabilization.

`lsd_threshold`

A numeric value defining the stopping criterion based on the Local Standardized mean Difference (LSD). Training stops when the rolling average of the maximum local balance falls below this threshold. The default `lsd_threshold = 2` balances efficiency and precision. In cases of poor overlap or small sample sizes, a more relaxed threshold (e.g., `5%` or `10%`) may be used to allow more flexibility in training.

`rolling_window`

An integer specifying the number of recent epochs used to compute the rolling average of the maximum local balance. Default is `5`. The early stopping mechanism is triggered when the rolling average of the maximum LSD over the most recent `rolling_window` epochs falls below `lsd_threshold`. Specifically, at every 200-epoch step, the maximum local balance is calculated, and a rolling average over the last `rolling_window` steps is updated. Training halts when this rolling average drops below `lsd_threshold`, or when the predefined maximum epochs is reached, ensuring sufficient learning capacity.

setup_lbcnet_args

List. Optional arguments passed to setup_lbcnet for configuring the Python environment. If Python is not set up, `setup_lbcnet()` is automatically called using these parameters. Default is `list(envname = "r-lbcnet", create_if_missing = TRUE)`, meaning it will attempt to use the virtual environment `"r-lbcnet"` and create it if missing.

Value

An object of class "lbc_net_surv" (and "lbc_net"), containing:

• fitted.values: Estimated propensity scores \(\hat e(X)\).

• weights: IPW weights

• survival: A list with components:

  • times: Time grid at which survival functions are evaluated.

  • S1: Estimated survival function under treatment.

  • S0: Estimated survival function under control.

  • diff: Estimated survival difference S1 - S0 at each time.

  • se_diff: Standard error of the survival difference.

  • ci_lower: Lower bound of the 95% Wald CI for the survival difference.

  • ci_upper: Upper bound of the 95% Wald CI for the survival difference.

  • t_star: The highlighted time point (median control survival by default, or user-supplied).

Details

Let \(T\) be the event or censoring time, \(\Delta\) the event indicator (1 = event, 0 = censored), and \(A \in \{0,1\}\) the treatment. LBC-Net first estimates the propensity score \(e(X) = P(A = 1 \mid X)\) using the same deep learning framework as lbc_net, with local balance and calibration constraints.

The resulting propensity scores \(\hat e(X)\) are used to form Inverse Probability Weights $$ W_i = \frac{1}{A_i \hat e(X_i) + (1 - A_i)\{1 - \hat e(X_i)\}}, $$ corresponding to the frequency weight function \(\omega^{*}(p) = 1\) (no additional tilting). Using these weights, Deng and Wang's (2025) IPW–Nelson–Aalen estimator is applied separately for the treated and control groups to obtain cumulative hazards \(\hat\Lambda_1(t)\), \(\hat\Lambda_0(t)\), and survival functions $$ \hat S_a(t) = \exp\{-\hat\Lambda_a(t)\}, \quad a \in \{0,1\}. $$

The function returns the estimated survival difference $$ \hat\Delta(t) = \hat S_1(t) - \hat S_0(t) $$ evaluated at a user-specified grid of time points `t_grid`. If the user does not provide a grid or a specific `t_star`, the default evaluation time is set to the median survival time of the control arm, estimated from a standard Kaplan–Meier curve based on observed data.

At this stage, `lbc_net_surv()` provides point estimates of \(\hat S_1(t)\), \(\hat S_0(t)\), and their difference \(\hat\Delta(t)\) along with variance estimation and influence-function-based standard errors for the survival difference at each evaluation time.

Note

Compared with lbc_net, this function:

• uses time and delta instead of an outcome Y,

• ignores outcome-related arguments such as Y, outcome_type, and estimand flags (ATE/ATT/Y).

See also

lbc_net for propensity score estimation in non-survival settings.

Examples

if (FALSE) { # \dontrun{
  set.seed(123)
  n  <- 1000
  X1 <- rnorm(n)
  X2 <- rnorm(n)
  Tr <- rbinom(n, 1, plogis(0.5 * X1 - 0.3 * X2))

  lambda0 <- 0.02
  hr      <- 0.7
  rate    <- lambda0 * ifelse(Tr == 1, hr, 1)
  T_true  <- rexp(n, rate = rate)
  C       <- rexp(n, rate = 0.01)

  time  <- pmin(T_true, C)
  delta <- as.integer(T_true <= C)

  dat <- data.frame(Tr = Tr, X1 = X1, X2 = X2,
                    time = time, delta = delta)

  fit_surv <- lbc_net_surv(
    data    = dat,
    formula = Tr ~ X1 + X2,
    time    = "time",
    delta   = "delta"
  )

  # Estimated survival difference and SE at evaluation times
  head(cbind(
    time = fit_surv$survival$times,
    diff = fit_surv$survival$diff,
    se   = fit_surv$survival$se_diff
  ))
} # }