RobinCar Family

This vignette accompanies the Results section of the RobinCar Family paper.

library(RobinCar)
library(RobinCar2)

## 
## Attaching package: 'RobinCar2'

## The following object is masked from 'package:base':
## 
##     table

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

Data Manipulation

Creates continuous and binary outcomes, and does simple data manipulations.

# 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 ' ' 1

Generalized Linear Models

For the binary outcome, we can use a logistic regression model instead.

robin_glm(
  y_bin ~ arms * (wtkg + hemo + oprior),
  treatment = arms ~ pb(strat),
  family = binomial(link = "logit"),
  contrast = "risk_ratio",
  data = data
)

## Warning in robin_glm(y_bin ~ arms * (wtkg + hemo + oprior), treatment = arms ~
## : Consider using the log transformation `log_odds_ratio` and `log_risk_ratio`
## to replace `odds_ratio` and `risk_ratio` to improve the performance of normal
## approximation.

## Model        :  y_bin ~ arms * (wtkg + hemo + oprior) 
## Randomization:  arms ~ pb(strat)  ( Permuted-Block )
## Variance Type:  vcovG 
## Marginal Mean: 
##    Estimate   Std.Err     2.5 % 97.5 %
## 0 0.0496651 0.0092998 0.0314378 0.0679
## 1 0.1828118 0.0168988 0.1496908 0.2159
## 
## Contrast     :  risk_ratio
##          Estimate Std.Err Z Value  Pr(>|z|)    
## 1 v.s. 0  3.68089 0.76609  3.4994 0.0004662 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Rather than a linear contrast, we can instead estimate a 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 ' ' 1

Mantel-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-12

Survival 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 ' ' 1

You 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    103

Marlena Bannick

2026-01-19