Package Introduction • RobinCar2

RobinCar2 provides robust inference on treatment effect estimates obtained through (generalized) linear models(glm) with variance link functions.

Here with the dummy data provided in the package, we demonstrate how to use the package.

Common Usage

library(RobinCar2)
head(dummy_data)
#>   id treatment s1 s2      covar           y y_b
#> 1  1       pbo  b  c  0.5119022 -0.02761963   0
#> 2  2       pbo  a  c -0.7941720  0.49919508   0
#> 3  3      trt1  b  d  0.8988804  0.48037375   0
#> 4  4      trt2  b  c -0.4821317  0.67490126   1
#> 5  5      trt1  a  d -0.2285514  0.55499267   0
#> 6  6      trt2  a  c  0.2742069  2.39830584   1

Data Introduction

In the dummy_data, we have the following columns:

id is the patient identifier.
treatment is the treatment assignment.
s1 is the first stratification factor.
s2 is the second stratification factor.
covar is the continuous covariate.
y is the continuous outcome.
y_b is the binary outcome.

Obtain Treatment Effect for Continuous Outcome

In simplest case of a continuous outcome, and we want to adjust for covar, s1, where covar is a covariate and s1 is a stratification factor. Randomization is permuted-block randomization stratified by s1. A proper model to use is y ~ treatment * s1 + covar.

robin_glm(y ~ treatment * s1 + covar, data = dummy_data, treatment = treatment ~ pb(s1))
#> Treatment Effect
#> -------------
#> Model        :  y ~ treatment * s1 + covar 
#> Randomization:  treatment ~ pb(s1) 
#> Marginal Mean: 
#> counter-factual prediction
#> 
#>       pbo      trt1      trt2 
#> 0.2003208 0.7639709 0.9712499 
#> 
#> Variance Type:  gvcov 
#>             Estimate Std.Err Z Value Pr(>|z|)    
#> trt1 - pbo     0.564   0.101    5.60  2.2e-08 ***
#> trt2 - pbo     0.771   0.101    7.61  2.8e-14 ***
#> trt2 - trt1    0.207   0.107    1.94    0.052 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We can also use the Huber-White variance estimator by setting vcov = "vcovHC". Please note that in this case, the model formula should not contain treatment interaction with stratification factors.

robin_glm(y ~ treatment + s1 + covar, data = dummy_data, treatment = treatment ~ pb(s1), vcov = "vcovHC")
#> Treatment Effect
#> -------------
#> Model        :  y ~ treatment + s1 + covar 
#> Randomization:  treatment ~ pb(s1) 
#> Marginal Mean: 
#> counter-factual prediction
#> 
#>       pbo      trt1      trt2 
#> 0.2004488 0.7639780 0.9712848 
#> 
#> Variance Type:  vcovHC 
#>             Estimate Std.Err Z Value Pr(>|z|)    
#> trt1 - pbo     0.564   0.101    5.60  2.1e-08 ***
#> trt2 - pbo     0.771   0.102    7.57  3.7e-14 ***
#> trt2 - trt1    0.207   0.108    1.92    0.055 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Obtain Treatment Effect for Binary Outcome

If the outcome is binary, and we want to adjust for covar, s1, where covar is a covariate and s1 is a stratification factor. Randomization is permuted-block randomization stratified by s1. A proper model to use is y_b ~ treatment * s1 + covar. Note here we need to specify family to be binomial(link = "logit").

robin_glm(y_b ~ treatment * s1 + covar, data = dummy_data, treatment = treatment ~ pb(s1), family = binomial(link = "logit"))
#> Treatment Effect
#> -------------
#> Model        :  y_b ~ treatment * s1 + covar 
#> Randomization:  treatment ~ pb(s1) 
#> Marginal Mean: 
#> counter-factual prediction
#> 
#>       pbo      trt1      trt2 
#> 0.3560965 0.5806957 0.6213865 
#> 
#> Variance Type:  gvcov 
#>             Estimate Std.Err Z Value Pr(>|z|)    
#> trt1 - pbo    0.2246  0.0477    4.71  2.5e-06 ***
#> trt2 - pbo    0.2653  0.0475    5.58  2.4e-08 ***
#> trt2 - trt1   0.0407  0.0479    0.85      0.4    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Obtain Treatment Effect for Counts

If the outcome is count, and we want to adjust for the same covariates, we can use family = MASS::negative.binomial(theta = NA), to correctly optimize over theta. Otherwise we can also have a fixed theta if we do know the truth.

Other than negative binomial link function, it is also possible that we use Poisson link function.

dummy_data$y_count <- rpois(nrow(dummy_data), lambda = 20)
robin_glm(
  y_count ~ treatment * s1 + covar,
  data = dummy_data,
  treatment = treatment ~ pb(s1), family = MASS::negative.binomial(theta = 1)
)
#> Treatment Effect
#> -------------
#> Model        :  y_count ~ treatment * s1 + covar 
#> Randomization:  treatment ~ pb(s1) 
#> Marginal Mean: 
#> counter-factual prediction
#> 
#>      pbo     trt1     trt2 
#> 19.15274 20.06251 20.59907 
#> 
#> Variance Type:  gvcov 
#>             Estimate Std.Err Z Value Pr(>|z|)    
#> trt1 - pbo     0.910   0.453    2.01  0.04440 *  
#> trt2 - pbo     1.446   0.439    3.29  0.00099 ***
#> trt2 - trt1    0.537   0.466    1.15  0.24907    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Using Different Covariate-Adaptive Randomization Schema

If the randomization schema is not permuted-block randomization, we can use other randomization schema. Currently we have also sp(simple), pb(permuted-block), ps(Pocock-Simon).

robin_glm(y_b ~ treatment * s1 + covar, data = dummy_data, treatment = treatment ~ ps(s1), family = binomial(link = "logit"))
#> Treatment Effect
#> -------------
#> Model        :  y_b ~ treatment * s1 + covar 
#> Randomization:  treatment ~ ps(s1) 
#> Marginal Mean: 
#> counter-factual prediction
#> 
#>       pbo      trt1      trt2 
#> 0.3560965 0.5806957 0.6213865 
#> 
#> Variance Type:  gvcov 
#>             Estimate Std.Err Z Value Pr(>|z|)    
#> trt1 - pbo    0.2246  0.0477    4.71  2.5e-06 ***
#> trt2 - pbo    0.2653  0.0475    5.58  2.4e-08 ***
#> trt2 - trt1   0.0407  0.0479    0.85      0.4    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1