Detailed Function Presentation

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 = FALSE,
  message = FALSE
)

PDRobust

This vignette demonstrates validation, prediction, diagnostics, and treatment effect estimation using the bundled data.

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

1. Load the built-in package data

library(PDRobust)
data("BiSample", package = "PDRobust")
head(BiSample)

2. Define roles and analysis settings with Mapping()

mapping <- Mapping(
  id = "id",
  time = "time",
  treatment = "A",
  survival = "S",
  outcome = "Y",
  baseline_time = 0,
  cutoff_time = 2,
  covariates = c("X1", "X3", "X4", "X5", "X6"),
  interest_vars = c("X1", "X4"),
  y_type = "B"
)

print(mapping)

3. Validate the raw data with DataCheck()

check <- DataCheck(BiSample, mapping, strict = FALSE)
names(check)           

check$valid
check$ready_for_analysis
check$manual_resolution_required
check$can_standardize

The itemized report is in check$checks; supporting details are in check$diagnostics.

4. Standardize the panel with DataStandard()

pd_data <- DataStandard(BiSample, mapping, drop =TRUE)
head(pd_data)

Imperfect dataset

For imperfect dataset, we have:

data("ImperfectConSample", package = "PDRobust")
head(ImperfectConSample)

con_mapping <- Mapping(
  id = "patient_id",
  time = "visit_month",
  treatment = "treatment",
  survival = "alive_status",
  outcome = "clinical_outcome",
  baseline_time = 0,
  cutoff_time = 12,
  covariates = c("X1", "X2", "X3", "X4", "X5", "X6"),
  interest_vars = c("X1", "X2"),
  y_type = "C"
)

con_check <- DataCheck(ImperfectConSample, con_mapping, strict = FALSE)

con_check$valid
con_check$ready_for_analysis
con_check$manual_resolution_required
con_check$can_standardize

con_data <- DataStandard(ImperfectConSample, con_mapping, drop = TRUE)
head(con_data)

print(dim(ImperfectConSample))
print(dim(con_data))

names(attributes(con_data))

attr_standard <- attributes(con_data)
attr_standard$pd_standardization$time_map
head(attr_standard$pd_standardization$id_map)

5.1 Prediction functions and Diagnostics

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 + S

Propensity score model

ps <- PSPred(
  ps_fo = ps_fo,
  fit_dat = pd_data,
  pred_dat = pd_data,
  mapping = mapping
)

head(ps)

ps_diagnostic <- PSDiag(data = pd_data,
                        ps_fo = ps_fo)

print(ps_diagnostic)

Principal score model

p0 <- PrinPred(
  prin_fo = prin_fo,
  fit_dat = pd_data,
  pred_dat = pd_data,
  a = 0,
  mapping = mapping
)

head(p0)

principal_diagnostic <- PrinSDiag(
  data = pd_data, 
  ps_fo = ps_fo, 
  prin_fo = prin_fo)

print(principal_diagnostic)

Outcome model

mu1 <- OutPred(
  out_fo = out_fo,
  fit_dat = pd_data,
  pred_dat = pd_data,
  a = 1,
  mapping = mapping
)

head(mu1)

set.seed(12345)
sensitivity <- SA(
  data  = pd_data,
  ps_fo = ps_fo,
  prin_fo = prin_fo,
  out_fo = out_fo,
  ratiovec = c(0.05,0.1, 0.2)
)
print(sensitivity)

Principal-stratum profiling with QR()

principal_profile <- QR(
  data = pd_data,
  prin_fo = prin_fo,
  quantile_level = c(0.25, 0.50, 0.75)
)

print(principal_profile)
principal_profile$data

Treatment-group odds ratios

or_control <- ORCI(
  data = pd_data,
  formula = S ~ X1 + X3 + X4,
  a = 0,
  conf_level = 0.95
)

print(or_control)

5.2 Heterogeneous treatment effect

The five bootstrap replications below are only for a fast demonstration. Substantive standard errors and confidence intervals require more replications and an assessment of their stability. Use B = 0 for point estimates alone.

set.seed(12345)
separate_hte <- 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
)

separate_hte$summary
separate_hte$forest_plot

head(separate_hte$boot_mat)

pooled_hte <- HTEAllT(
  data = pd_data,
  ps_fo = ps_fo,
  prin_fo = prin_fo,
  out_fo = out_fo,
  B = 0,
  verbose = FALSE
)
pooled_hte$summary
pooled_hte$forest_plot

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