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

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

## -----------------------------------------------------------------------------
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)

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

## -----------------------------------------------------------------------------
check$valid
check$ready_for_analysis
check$manual_resolution_required
check$can_standardize

## -----------------------------------------------------------------------------


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

## -----------------------------------------------------------------------------
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)

## -----------------------------------------------------------------------------
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

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

## ----ps_dgn, fig.alt = "Absolute standardized mean differences before and after propensity-score weighting."----
ps_diagnostic <- PSDiag(data = pd_data,
                        ps_fo = ps_fo)

print(ps_diagnostic)

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

head(p0)

## ----pps_dgn, fig.alt = "Standardized principal-score balance statistics for the selected covariates."----
principal_diagnostic <- PrinSDiag(
  data = pd_data, 
  ps_fo = ps_fo, 
  prin_fo = prin_fo)

print(principal_diagnostic)

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

head(mu1)

## ----sa, fig.alt = "Estimated effect-modification coefficients over time at different outcome-noise variance ratios."----
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)

## ----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

## ----or_ci, fig.alt = "Cutoff survival odds ratios and confidence intervals within treatment group zero."----
or_control <- ORCI(
  data = pd_data,
  formula = S ~ X1 + X3 + X4,
  a = 0,
  conf_level = 0.95
)

print(or_control)

## ----htesept, fig.alt = "Time-specific treatment-effect model coefficients and demonstration bootstrap confidence intervals."----
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)

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