PropensitySub (version 0.2.0)

bootstrap_propen: Calculate bootstrap CI for treatment effect estimate

Description

Calculate bootstrap CI for treatment effect estimate

Usage

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
)

Arguments

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

Value

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

Examples

# NOT RUN {
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 

# }