Family.RmdThis vignette accompanies the Results section of the RobinCar Family paper.
library(RobinCar)
library(RobinCar2)
library(dplyr)
library(SuperLearner)
library(ranger)
library(xgboost)
library(glmnet)Data Manipulation
We use data from the Juraska et al. (2022) R package, which comes from a study of nucleoside treatment regimens for individuals with HIV-1 (Hammer et al. 1996) called the AIDS Clinical Trials Group Study 175. This trial had stratified permuted block randomization.
# Data are from the speff2trial package
data <- tibble(speff2trial::ACTG175)
data <- data %>% mutate(
# Create continuous and binary outcomes
y_cont = cd420 - cd40,
y_bin = as.numeric((y_cont / cd40) > 0.5)
) %>% filter(
# Focus on two treatment arms
arms %in% c("0", "1")
) %>% mutate(
# Factor variables for treatment
# and stratification variable
arms = factor(arms),
strat = factor(strat)
)Linear Adjustment
The following shows a linear model with treatment-by-covariate interactions, including strata as covariates, and using the four covariates of weight, hemophilia status, and prior use of non-zidovudine antiretroviral therapy. These specifications ensure that the covariate adjustment has guaranteed efficiency gain (asymptotically).
robin_lm(
y_cont ~ arms * (strat + wtkg + hemo + oprior),
treatment = arms ~ pb(strat),
data=data)## Model : y_cont ~ arms * (strat + wtkg + hemo + oprior)
## Randomization: arms ~ pb(strat) ( Permuted-Block )
## Variance Type: vcovG
## Marginal Mean:
## Estimate Std.Err 2.5 % 97.5 %
## 0 -16.9304 4.4915 -25.7336 -8.1273
## 1 54.0531 6.2848 41.7352 66.3710
##
## Contrast : h_diff
## Estimate Std.Err Z Value Pr(>|z|)
## 1 v.s. 0 70.984 7.673 9.2511 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Generalized Linear Models
For the binary outcome, we can use a logistic regression model instead. Rather than a linear contrast, we could instead estimate a log risk ratio. This does not change the marginal mean estimates, but it changes the contrast estimates (and p-values).
robin_glm(
y_bin ~ arms * (strat + wtkg + hemo + oprior),
treatment = arms ~ pb(strat),
family = binomial(link = "logit"),
contrast = "log_risk_ratio",
data=data
)## Model : y_bin ~ arms * (strat + wtkg + hemo + oprior)
## Randomization: arms ~ pb(strat) ( Permuted-Block )
## Variance Type: vcovG
## Marginal Mean:
## Estimate Std.Err 2.5 % 97.5 %
## 0 0.0493622 0.0093041 0.0311264 0.0676
## 1 0.1835664 0.0168944 0.1504539 0.2167
##
## Contrast : log_risk_ratio
## Estimate Std.Err Z Value Pr(>|z|)
## 1 v.s. 0 1.31339 0.20904 6.2831 3.318e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Mantel-Haenszel
RobinCar includes functions for robust Mantel-Haenszel estimation. Here we create a new variable that is the interaction of prior use of non-zidovudine antiretroviral therapy and hemophilia status, and use that as the stratification variable.
By default, the function uses the average treatment effect as the estimand of interest. Specifying the estimand argument as “MH” does not change the estimate, but it does change the standard error calculation.
The MH function in RobinCar is only valid for simple randomization.
data <- data %>% mutate(
X=interaction(oprior, hemo)
)
robincar_mh(
data, "arms", "y_bin",
strata_cols=c("oprior", "hemo")
)## Treatment group contrasts based on ATE
## Estimand: ATE
## Stratified by oprior, hemo
## SE calculated via modified Greenland's estimator
##
## Contrasts:
## # A tibble: 1 × 4
## contrast estimate se `pval (2-sided)`
## <chr[1d]> <dbl> <dbl> <dbl>
## 1 treat 1 - 0 0.134 0.0193 3.83e-12Survival Analysis
For a survival outcome with right-censoring, we can use the RobinCar2
function robin_surv to do a covariate-adjusted, stratified
log-rank test.
surv <- robin_surv(
Surv(days, cens) ~ wtkg + oprior + hemo + strat,
treatment = arms ~ pb(strat),
data = data
)
surv## Model : Surv(days, cens) ~ wtkg + oprior + hemo + strat
## Randomization: arms ~ pb(strat) (Permuted-Block)
## Covariates adjusted for: wtkg, oprior, hemo, strat (including interactions with arms)
##
## Contrast : Covariate-adjusted Log Hazard Ratio
##
## Estimate Std.Err Z Value Pr(>|z|)
## 1 v.s. 0 -0.68785 0.12237 -5.6209 1.9e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Test : Covariate-adjusted Log-Rank
##
## Test Stat. Pr(>|z|)
## 1 v.s. 0 -5.7444 9.223e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1You can use the table function to see the number of
events and number at-risk
table(surv)## Number of patients and events per treatment arm:
## arms Patients Events
## 1 0 532 181
## 2 1 522 103Data Manipulation
This case study uses a different dataset,
medicaldata::indo_rct, which was a randomized trial of
rectal indomethacin to prevent post-ERCP pancreatitis (Elmunzer,
Higgins, et al. 2012). This study used simple randomization stratified
by study site (this is not the same as stratified permuted block
randomization in the first example).
## Loading required package: medicaldata
data2 <- medicaldata::indo_rct %>%
mutate(
outcome2=as.numeric(outcome)-1
) %>%
select(
-c(outcome, status, type, bleed)
)First fit a really simple logistic regression working model only adjusting for site, and a risk difference as our contrast of interest.
robin_glm(
outcome2 ~ rx * site,
treatment = rx ~ sr(1),
data=data2,
family=binomial(link="logit"),
contrast = "difference"
)## Model : outcome2 ~ rx * site
## Randomization: rx ~ sr(1) ( Simple )
## Variance Type: vcovG
## Marginal Mean:
## Estimate Std.Err 2.5 % 97.5 %
## 0_placebo 0.167499 0.021227 0.125894 0.2091
## 1_indomethacin 0.092527 0.016766 0.059666 0.1254
##
## Contrast : difference
## Estimate Std.Err Z Value Pr(>|z|)
## 1_indomethacin v.s. 0_placebo -0.074971 0.026910 -2.786 0.005337 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There are 27 baseline covariates included in this dataset and 602 patients. If we include all of the covariates in a logistic regression working model, the model fails to converge and returns NA.
# Create formula with all baseline covariates included in dataset
covariates <- setdiff(colnames(data2), c("id", "outcome2", "rx"))
formula <- as.formula(
sprintf("outcome2 ~ rx * (%s)", paste0(covariates, collapse=" + "))
)
robin_glm(
as.formula(formula),
treatment = rx ~ sr(1),
data=data2,
family=binomial(link="logit"),
contrast = "difference"
)## Model : outcome2 ~ rx * (site + age + risk + gender + sod + pep + recpanc + psphinc + precut + difcan + pneudil + amp + paninj + acinar + brush + asa81 + asa325 + asa + prophystent + therastent + pdstent + sodsom + bsphinc + bstent + chole + pbmal + train)
## Randomization: rx ~ sr(1) ( Simple )
## Variance Type: vcovG
## Marginal Mean:
## Estimate Std.Err 2.5 % 97.5 %
## 0_placebo NA NA NA NA
## 1_indomethacin NA NA NA NA
##
## Contrast : difference
## Estimate Std.Err Z Value Pr(>|z|)
## 1_indomethacin v.s. 0_placebo NA NA NA NAIf you do want to include all of the baseline covariates because you aren’t sure which ones are the most prognostic, this is a use-case for the SuperLearner approach. We’ll test an ensemble of random forest (SL.ranger), elastic net regression (Sl.glmnet), and XGboost (SL.xgboost). SuperLearner will figure out the optimal combination of these approaches. See this link for documentation on the SuperLearner.
RobinCar includes SuperLearner working models with cross-fitting to ensure we aren’t over-fitting to the data. We will use 5 folds for cross-fitting. Since cross-fitting is a random process, we’ll also set a seed for reproducibility.
set.seed(1989)
sl <- robincar_SL(
df=data2,
response_col="outcome2",
treat_col="rx",
covariate_cols=covariates,
SL_libraries=c("SL.xgboost", "SL.glmnet", "SL.ranger"),
car_scheme="simple",
k_split=5
)## Done!## Warning in func(): Prediction unbiasedness does not hold.
# Make a linear contrast
contrast <- robincar_contrast(sl, contrast_h="diff")
contrast## Treatment group contrasts using linear contrast
##
## Contrasts:
## # A tibble: 1 × 4
## contrast estimate se `pval (2-sided)`
## <chr> <dbl> <dbl> <dbl>
## 1 treat 1_indomethacin - 0_placebo -0.0812 0.0269 0.00256
##
## Variance-Covariance Matrix for Contrasts:
## [,1]
## [1,] 0.0007245833