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

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

## -----------------------------------------------------------------------------
ps_diag <- PSDiag(data = pd_data, 
                  ps_fo = ps_fo)
names(ps_diag)

## ----fig.alt = "Absolute standardized mean differences before and after propensity-score weighting."----
print(ps_diag)

## -----------------------------------------------------------------------------
ps_diag$weight_type

## -----------------------------------------------------------------------------
prin_diag <- PrinSDiag(data = pd_data, 
                       ps_fo = ps_fo, 
                       prin_fo = prin_fo)
names(prin_diag)

## ----fig.alt = "Standardized principal-score balance statistics for the selected covariates."----
print(prin_diag)

## -----------------------------------------------------------------------------
head(prin_diag$p0)
head(prin_diag$p1)

## -----------------------------------------------------------------------------
set.seed(20160878)
sa <- SA(
  data = pd_data,
  ps_fo = ps_fo,
  prin_fo = prin_fo,
  out_fo = out_fo,
  ratiovec = c(0, 0.05, 0.1)
)
names(sa)

## ----fig.alt = "Estimated effect-modification coefficients over time at different outcome-noise variance ratios."----
print(sa)

## -----------------------------------------------------------------------------
sa$variance_by_time
sa$warnings

## -----------------------------------------------------------------------------
profile <- QR(
  data = pd_data,
  prin_fo = prin_fo,
  quantile_level = c(0.5, 0.95)
)

names(profile)

## -----------------------------------------------------------------------------
print(profile)

## -----------------------------------------------------------------------------
head(profile$weights)

## -----------------------------------------------------------------------------
or_fo <- S ~ X1 + X2 + X4
or0 <- ORCI(data = pd_data, 
            formula = or_fo,
            a = 0,
            conf_level = 0.95)

names(or0)

## ----fig.alt = "Cutoff survival odds ratios and confidence intervals within treatment group zero."----
print(or0)

## -----------------------------------------------------------------------------
or0$model_diagnostics