bootstrap_propen: Calculate bootstrap CI for treatment effect estimate
In PropensitySub: Treatment Effect Estimate in Strata with Missing Data

Description Usage Arguments Value Note Examples

Calculate bootstrap CI for treatment effect estimate

bootstrap_propen(
  data.in,
  indicator.var = "indicator",
  formula,
  indicator.next = NULL,
  seed = 100,
  class.of.int,
  estimate.res,
  n.boot = 1000,
  method = "ipw",
  wild.boot = FALSE,
  tte = "AVAL",
  event = "event",
  trt = "trt",
  response = NULL,
  caliper = NULL,
  pairs = NULL,
  hr.ratio.ref = NULL,
  ref.denom = TRUE,
  model = "plain",
  max.num.run = 5000,
  non.converge.check = FALSE,
  multinom.maxit = 100,
  non.converge.check.thresh = 1
)

`data.in`	(`data.frame`) input data
`indicator.var`	(`string`) column name of the strata indicator variable which must be numeric. Assume arm1 has strata labeling and arm2 does not have strata labeling. pts without strata labeling should be indicated as -1 (e.g. pts in the arm1, or pts in arm2 but with missing label). within arm1 (the arm with strata labeling), subclasss should be indicated as 0,1,2...
`formula`	(`formula`) to input to the logistic or multinomial logistic model (in the form of strata~features)
`indicator.next`	(`string`) column name of the column which indicates status at a different measurement. It should be coded in the same way as in indicator.var (e.g. -1, 0, 1). Patients who have both missing current status and missing next status should be excluded in the modeling.
`seed`	seed
`class.of.int`	(`list`) classes (stratum) of interest. Request to be in list format. It could be subset of classes in arm1; it could also define combined classes. For example: class.of.int = list("class1"=0, "class2"=1, "class3"=2, "class2or3"="c(1,2)"). for "class2or3", Prob(class 2 or 3) will be calculated as Prob(class2) + Prob(class3)
`estimate.res`	result object from ipw_strata() or ps_match_strata()
`n.boot`	number of bootstraps to run; note only runs without warning/error msg will be considered when calculating bootstrap CI; so it is possible more that n.boot runs are performed (but capped by max.num.run)
`method`	"ipw" or "ps". If "ipw", ipw_strata() will be used. If "ps", ps_match_strata() will be used.
`wild.boot`	whether wild bootstrap should be used. If so, weights will be generated using rexp(1)
`tte`	(`string`) column name of the time to event variable
`event`	(`string`) column name of the event variable (1: event, 0: censor)
`trt`	(`string`) column name of the treatment variable. The variable is expected to be in `factor` format and the first level will be considered as reference (control) arm when calculating summary statistics.
`response`	(`string`) column name of the response variable. 1 indicates responder and 0 indicates non responder. if response is not NULL, tte and event will be ignored and the function will assume binary outcome.
`caliper`	A scalar or vector denoting the caliper(s) which should be used when matching. A caliper is the distance which is acceptable for any match. Observations which are outside of the caliper are dropped. If a scalar caliper is provided, this caliper is used for all covariates in `X`. If a vector of calipers is provided, a caliper value should be provided for each covariate in `X`. The caliper is interpreted to be in standardized units. For example, `caliper=.25` means that all matches not equal to or within .25 standard deviations of each covariate in `X` are dropped. Note that dropping observations generally changes the quantity being estimated.
`pairs`	pairs of interest when calculating ratio of HR (delta of delta for OR). this should be a matrix whose rows are names of strata, 1st column indicates the stratum to be used as numerator (HR or ORR diff); 2nd column indicates denominator. If pairs is NULL, ratio of HR (difference of OR difference) will not be calculated.
`hr.ratio.ref`	no longer to be used, please use pairs instead
`ref.denom`	no longer to be used, please use pairs instead
`model`	(`string`) one of (`plain`, `dwc`, `wri`). "plain" when 2 levels are specified in indicator variable, a binomial glm will be fitted; when more than 2 levels are specified, a multinomial glm will be fitted; "dwc" Doubly weighted control: Two separated models will be fitted: one is binomial glm of 2 vs. (1, 0), the other one is bionomial glm of 1 vs 0. The probability of being each class is then calculated by aggregating these two models. Note this is similar to the plain method but with different (less rigid) covariance assumptions. "wri" Weight regression imputation: the current status is going to be learned from the next status. Indicator of the next status should be specified using indicator.next. Currently "wri" only support the case where there are only two non-missing strata. In indicator variable, the two nonmissing strata should be coded as 0 and 1, the missing group should be coded as 2.
`max.num.run`	max number of bootstraps to run (include both valid and not valid runs)
`non.converge.check`	whether to output number of time each level of each categorical variable for each stratum specified in indicator having N<non.converge.check.thresh when non-convergence occurs
`multinom.maxit`	see parameter `maxit` in nnet::multinom, default is 100
`non.converge.check.thresh`	see above

return a list containing the following components:

boot.out.est a matrix with rows as estimates such as Coefficient and Variance in strata and columns as summary statistics such as Mean and Median of the estimates.
est.ci.mat a matrix with rows as strata and columns as Estimate and Bootstrap CIs.
eff.comp.ci.mat a matrix with rows as strata comparisons and columns as Estimate and Bootstrp CIs.
conv.est a logical vector to indicate whether model in each bootstrap converges.
error.est numeric to indicate the total number of models in bootstrap which gives errors.
boot.results a matrix with rows as each bootstrap and columns as model results such as Coefficient in strata.
glm.warn.est a logical vector to indicate whether glm model gives warning in each bootstrap.
num.valid.boots numeric to indicate the total number of valid bootstraps.
num.total.boots numeric for the total number of bootstrap runs.
warning.msg a list to capture warning messages from models.
error.msg a list to capture error messages from models.
non.converge.dat a matrix with rows as each level of each categorical variable for each stratum specified in indicator having N less than non.converge.check.thresh and columns as treatment groups

only estimates from runs without error or warning will be considered when calculating bootstrap CI. If none of the bootstrap runs is error/warning free, CI of est.ci.mat will be NA

library(dplyr)
clinical_1 <- clinical %>% mutate( 
  indicator = case_when(
    STRATUM == "strata_1" ~ 0, 
    STRATUM == "strata_2" ~ 1,
    is.na(STRATUM) & ARM == "experimental" ~ 1,
    TRUE ~ -1 
  ),
  ARM  = factor(ARM, levels = c("control","experimental")),
  BNLR = case_when(
    is.na(BNLR) ~ median(BNLR, na.rm = TRUE),
    TRUE ~ BNLR
  )
)
# ipw: default model
ipw_res <- ipw_strata(
  data.in = clinical_1, formula = indicator ~ BECOG + SEX + BNLR,
  indicator.var = "indicator", tte = "OS_MONTH", event = "OS_EVENT", trt = "ARM",
  class.of.int = list("strata_1" = 1, "strata_2" = 0)
 )
boot_ipw <- bootstrap_propen(
  data.in = clinical_1, formula = indicator ~ BECOG + SEX + BNLR,
  indicator.var = "indicator", tte = "OS_MONTH", event = "OS_EVENT", trt = "ARM",
  class.of.int = list("strata_1" = 1, "strata_2" = 0),
  estimate.res = ipw_res, method = "ipw", n.boot = 5
)
boot_ipw$est.ci.mat
boot_ipw$boot.out.est 
# ps: DWC model
clinical_2 <- clinical %>% mutate( 
  indicator = case_when(
    STRATUM == "strata_1" ~ 0, 
    STRATUM == "strata_2" ~ 1,
    is.na(STRATUM) & ARM == "experimental" ~ 2,
    TRUE ~ -1
  ),
  ARM  = factor(ARM, levels = c("control","experimental")),
  BNLR = case_when(
    is.na(BNLR) ~ median(BNLR, na.rm = TRUE),
    TRUE ~ BNLR
  )
)
ps_res <- ps_match_strata(
  data.in = clinical_2, formula = indicator ~ BECOG + SEX + BNLR, model = "dwc",
  indicator.var = "indicator", tte = "OS_MONTH", event = "OS_EVENT", trt = "ARM",
  class.of.int = list("strata_1" = 0, "strata_2" = 1, "missing" = 2)
 ) 
boot_ps <- bootstrap_propen(
  data.in = clinical_2, formula = indicator ~ BECOG + SEX + BNLR, model = "dwc",
  indicator.var = "indicator", tte = "OS_MONTH", event = "OS_EVENT", trt = "ARM",
  class.of.int = list("strata_1" = 0, "strata_2" = 1, "missing" = 2),
  estimate.res = ps_res, method = "ps", n.boot = 5
)
boot_ps$est.ci.mat
boot_ps$boot.out.est