Bootstrap inference for a prespecified personalization / recommendation model

Runs B bootstrap samples using a prespecified model then computes the two I estimates based on cross validation. p values of the two I estimates are computed for a given \(H_0: \mu_{I_0} = \mu_0\) and confidence intervals are provided using the basic, percentile methods by default and the BCa method as well if desired.

Usage

PTE_bootstrap_inference(
  X,
  y,
  regression_type = "continuous",
  incidence_metric = "odds_ratio",
  personalized_model_build_function = NULL,
  censored = NULL,
  predict_function = NULL,
  difference_function = NULL,
  cleanup_mod_function = NULL,
  y_higher_is_better = TRUE,
  verbose = FALSE,
  full_verbose = FALSE,
  H_0_mu_equals = NULL,
  pct_leave_out = 0.1,
  m_pow_of_n = 0.75,
  B = 3000,
  alpha = 0.05,
  run_bca_bootstrap = FALSE,
  display_adversarial_score = FALSE,
  num_cores = NULL
)

Arguments

X

A \(n \times p\) dataframe of covariates where one column is labeled "treatment" and it is a binary vector of treatment allocations in the study.

y

An \(n\)-length numeric vector which is the response

regression_type

A string indicating the regression problem. Legal values are "continous" (the response y is a real number with no missing data, the default), "incidence" (the reponse y is either 0 or 1) and "survival". If the type is "survival", the user must also supply additional data via the parameter censored.

incidence_metric

Ignored unless the regression_type is "incidence" and difference_function is set to NULL (in the latter case, you have specified a more custom metric). Then, this parameter allows the user to select which of the three standard metrics to use for comparison: "probability_difference", "risk_ratio", "odds_ratio" where the default is "odds_ratio".

personalized_model_build_function

An R function that will be evaluated to construct the personalized medicine / recommendation model. In the formula for the model, the response is "y", the treatment vector is "treatment" and the data is "Xytrain". This function must return some type of object that can be used for prediction later via predict_function. Here are the defaults for each regression_type. They are linear models with first order interactions:


                                    personalized_model_build_function = switch(regression_type,
                                        continuous = function(Xytrain) { #defalt is OLS regression
                                            lm(y ~ . * treatment,
                                                data = Xytrain)
                                        },
                                        incidence = function(Xytrain) { #default is logistic regression
                                            glm(y ~ . * treatment,
                                                data = Xytrain,
                                                family = "binomial")
                                        },
                                        survival = function(Xytrain) { #default is Weibull regression
                                            survreg(Surv(Xytrain$y, Xytrain$censored) ~ (. - censored) * treatment,
                                                data = Xytrain,
                                                dist = "weibull")
                                        }
                                    )

censored

Only required if the regression_type is "survival". In this case, this vector is of length \(n\) and is binary where 0 indicates censored and 1 indicates uncensored. In a clinical trial, someone who is still alive at the end of the study or was lost to follow up will receive a censor value of 0, while someone who died during the study will receive a censor value of 1. \(n\) and is binary where 0 indicates censorship (e.g. the patient died).

predict_function

An R function that will be evaluated on left out data after the model is built with the training data. This function uses the object "mod" that is the result of the personalized_model_build_function and it must make use of "Xyleftout", a subset of observations from X. This function must return a scalar numeric quantity for comparison. The default function is predict(mod, obs_left_out) e.g. the default looks like:


                                    function(mod, Xyleftout) {
                                        predict(mod, Xyleftout)
                                    }

difference_function

A function which takes the result of one out of sample experiment (boostrap or not) of all n samples and converts it into a difference that will be used as a score in a score distribution to determine if the personalization model is statistically significantly able to distinguish subjects. The function looks as follows:


                                    function(results, indices_1_1, indices_0_0, indices_0_1, indices_1_0) {
                                        ...
                                        c(rec_vs_non_rec_diff_score, rec_vs_all_score, rec_vs_best_score)
                                    }

where results is a matrix consisting of columns of the estimated response of the treatment administered, the estimated response of the counterfactual treatment, the administered treatment, the recommended treatment based on the personalization model, the real response, and if this subject was censored (0 if so). Here are a couple of example entries:

est_true est_counterfactual given_tx rec_tx real_y censored 166.8 152.2 1 1 324 1 1679.1 2072.0 1 0 160 0

The arguments indices_1_1, indices_0_0, indices_0_1, indices_1_0 give the indices of the subjects whose treatment was administered as 1 and whose optimal was 1, whose treatment was administered 0 and whose optimal was 0, etc.

This function should return three numeric scores: the recommend vs. the non-recommended (adversarial), the recommended vs. all (all) and the recommended vs. the best average treatment (best) as a 3-dimensional vector as illustrated above.

By default, this parameter is NULL which means for continuous and incidence the average difference is used and for survival, the median Kaplan-Meier survival is used.

cleanup_mod_function

A function that is called at the end of a cross validation iteration to cleanup the model in some way. This is used for instance if you would like to release the memory your model is using but generally does not apply. The default is NA for "no function."

y_higher_is_better

True if a response value being higher is clinically "better" than one that is lower (e.g. cognitive ability in a drug trial for the mentally ill). False if the response value being lower is clinically "better" than one that is higher (e.g. amount of weight lost in a weight-loss trial). Default is TRUE.

verbose

Prints out a dot for each bootstrap sample. This only works on some platforms.

full_verbose

Prints out full information for each cross validation model for each bootstrap sample. This only works on some platforms.

H_0_mu_equals

The \(\mu_{I_0}\) value in \(H_0\). Default is NULL which specifies 0 for regression types continuous, survival and incidence (with incidence metric "probability_difference") or 1 if the regression type is incidence and the incidence metric is "risk_ratio" or "odds_ratio". These defaults essentially answer the question: does my allocation procedure do better than the business-as-usual / naive allocation procedure?

pct_leave_out

In the cross-validation, the proportion of the original dataset left out to estimate out-of-sample metrics. The default is 0.1 which corresponds to 10-fold cross validation.

m_pow_of_n

Within each bootstrap sample, the exponent \(\gamma\) used to set the resample size (sampled with replacement) to \(m = \lfloor n^{\gamma} \rfloor\), implementing the "m-out-of-n bootstrap." Politis, Romano & Wolf (Subsampling, 1999) show that bootstrap consistency for non-smooth (non-regular) functionals – as arises here because the treatment recommendation is a hard threshold on the estimated treatment effect – requires \(m \to \infty\) and \(m / n \to 0\), i.e. \(\gamma < 1\); \(\gamma = 1\) recovers the classical (and here potentially inconsistent) \(n\)-out-of-\(n\) bootstrap. The default \(\gamma = 0.75\) follows the Politis-Romano-Wolf recommended rate as a practical compromise between bias (\(\gamma\) too close to 1, reintroducing the non-regularity) and excess bootstrap Monte Carlo variance (\(\gamma\) too small); see their book for the data-driven "minimum volatility" method for calibrating \(\gamma\) to a specific dataset. When \(\gamma < 1\), the reported confidence intervals are rescaled by \(\sqrt{m / n}\) as required by m-out-of-n bootstrap theory; this rescaling is a no-op when \(\gamma = 1\) (i.e. \(m = n\)).

B

The number of bootstrap samples to take. We recommend making this as high as you can tolerate given speed considerations. The default is 3000.

alpha

Defines the confidence interval size (1 - alpha). Defaults to 0.05.

run_bca_bootstrap

Do the BCA bootstrap as well. This takes double the time. It defaults to FALSE.

display_adversarial_score

The adversarial score records the personalization metric versus the deliberate opposite of the personalization. This does not correspond to any practical situation but it is useful for debugging. Default is FALSE.

num_cores

The number of cores to use in parallel to run the bootstrap samples more rapidly. Defaults to NULL which uses max(cores - 1, 1) where cores is detected once when the PTE package is loaded and cached in options(mc.cores = ...) (see ?options). Pass an explicit value here to override the cached default for a single call.

Value

A results object of type "PTE_bootstrap_results" that contains much information about the observed results and the bootstrap runs, including hypothesis testing and confidence intervals.

Author

Adam Kapelner

Examples

library(PTE)

##B and num_cores are kept small here so the examples run quickly; in practice make B as
##large as you can tolerate and omit num_cores to use all-but-one core.

##response: continuous
data(continuous_example)
X = continuous_example$X
y = continuous_example$y
pte_results = PTE_bootstrap_inference(X, y, regression_type = "continuous", B = 5, num_cores = 1)
pte_results
#>     I_random observed_est = 0.077,  pctile p-val = 0.167, 
#>       95% CI's: basic = [-0.02, 0.186], pctile = [-0.033, 0.173], 
#>     I_best observed_est = 0.065,  p val = 0.333, 
#>       95% CI's: basic = [0.014, 0.253], pctile = [-0.124, 0.116], 

##response: incidence. incidence_metric also accepts "risk_ratio" and "probability_difference"
##(see the incidence_metric parameter above); odds_ratio is illustrated here for speed.
##force incidence and pretend y came to you this way
y_incidence = ifelse(y > quantile(y, 0.75), 1, 0)
pte_results = PTE_bootstrap_inference(X, y_incidence,
    regression_type = "incidence",
    B = 5, num_cores = 1)
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
pte_results
#>     I_random observed_est = 1.155,  pctile p-val = 0.167, 
#>       95% CI's: basic = [0.506, 1.337], pctile = [0.974, 1.805], 
#>     I_best observed_est = 1.04,  p val = 0.333, 
#>       95% CI's: basic = [0.454, 1.264], pctile = [0.816, 1.626], 

##response: survival
data(survival_example)
pte_results = PTE_bootstrap_inference(survival_example$X, survival_example$y,
    censored = survival_example$censored,
    regression_type = "survival",
    B = 5, num_cores = 1)
#> Warning: 'varlist' has changed (from nvar=6) to new 7 after EncodeVars() -- should no longer happen!
#> Warning: 'varlist' has changed (from nvar=6) to new 7 after EncodeVars() -- should no longer happen!
pte_results
#>     I_random observed_est = 0.148,  pctile p-val = 0.167, 
#>       95% CI's: basic = [0.082, 0.222], pctile = [0.074, 0.214], 
#>     I_best observed_est = -0.041,  p val = 0.833, 
#>       95% CI's: basic = [-0.027, 0.012], pctile = [-0.094, -0.055],