RobinCar2 provides robust covariate adjustment methods for estimating and inferring treatment effects through generalized linear models (glm) under different randomization schema.
Common Usage
A minimal call of robin_lm()
and
robin_glm()
, consisting of only formula, data arguments and
the randomization scheme, will produce an object of class
treatment_effect
.
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 Outcomes Using the General Variance
For the continuous outcome y
, the linear model includes
covar
as a covariate, s1
as a stratification
factor. The randomization scheme is a permuted-block randomization
stratified by s1
. The model formula also includes the
treatment by stratification interaction as
y ~ treatment * s1 + covar
.
robin_lm(y ~ treatment * s1 + covar,
data = dummy_data,
treatment = treatment ~ pb(s1)
)
#> Model : y ~ treatment * s1 + covar
#> Randomization: treatment ~ pb(s1) ( Permuted-Block )
#> Variance Type: vcovG
#> Marginal Mean:
#> Estimate Std.Err 2.5 % 97.5 %
#> pbo 0.200321 0.067690 0.067651 0.3330
#> trt1 0.763971 0.075929 0.615152 0.9128
#> trt2 0.971250 0.076543 0.821228 1.1213
#>
#> Contrast : h_diff
#> Estimate Std.Err Z Value Pr(>|z|)
#> trt1 v.s. pbo 0.56365 0.10074 5.5952 2.203e-08 ***
#> trt2 v.s. pbo 0.77093 0.10133 7.6082 2.779e-14 ***
#> trt2 v.s. trt1 0.20728 0.10683 1.9402 0.05235 .
#> ---
#> 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 the treatment by stratification (covariate)
interaction.
robin_lm(y ~ treatment + s1 + covar,
data = dummy_data,
treatment = treatment ~ pb(s1),
vcov = "vcovHC"
)
#> Model : y ~ treatment + s1 + covar
#> Randomization: treatment ~ pb(s1) ( Permuted-Block )
#> Variance Type: vcovG
#> Marginal Mean:
#> Estimate Std.Err 2.5 % 97.5 %
#> pbo 0.200449 0.067690 0.067779 0.3331
#> trt1 0.763978 0.075930 0.615158 0.9128
#> trt2 0.971285 0.076539 0.821271 1.1213
#>
#> Contrast : h_diff
#> Estimate Std.Err Z Value Pr(>|z|)
#> trt1 v.s. pbo 0.56353 0.10074 5.5941 2.218e-08 ***
#> trt2 v.s. pbo 0.77084 0.10133 7.6074 2.796e-14 ***
#> trt2 v.s. trt1 0.20731 0.10683 1.9405 0.05232 .
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note that robin_glm
can also handle continuous outcomes
using the default family and the default link function
family = gaussian()
.
Obtain Treatment Effect for Binary Outcomes
For binary outcomes, the logistic model includes covar
as a covariate, s1
as a stratification factor. The
randomization scheme is a permuted-block randomization stratified by
s1
. The model formula also includes the treatment by
stratification interaction as 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")
)
#> Model : y_b ~ treatment * s1 + covar
#> Randomization: treatment ~ pb(s1) ( Permuted-Block )
#> Variance Type: vcovG
#> Marginal Mean:
#> Estimate Std.Err 2.5 % 97.5 %
#> pbo 0.356097 0.033599 0.290243 0.4219
#> trt1 0.580696 0.034418 0.513238 0.6482
#> trt2 0.621386 0.034019 0.554711 0.6881
#>
#> Contrast : h_diff
#> Estimate Std.Err Z Value Pr(>|z|)
#> trt1 v.s. pbo 0.224599 0.047711 4.7075 2.508e-06 ***
#> trt2 v.s. pbo 0.265290 0.047534 5.5810 2.391e-08 ***
#> trt2 v.s. trt1 0.040691 0.047941 0.8488 0.396
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Obtain Treatment Effect for Counts
For counts, the log link model includes covar
as a
covariate, s1
as a stratification factor. The randomization
scheme is a permuted-block randomization stratified by s1
.
The model formula also includes the treatment by stratification
interaction as y_count ~ treatment * s1 + covar
. Note here
we need to specify family
to be
poisson(link = "log")
to use the Poisson model or to be
MASS::negative.binomial(theta = NA)
to use the negative
binomial model. A fixed theta
could be provided if it is
known.
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)
)
#> Model : y_count ~ treatment * s1 + covar
#> Randomization: treatment ~ pb(s1) ( Permuted-Block )
#> Variance Type: vcovG
#> Marginal Mean:
#> Estimate Std.Err 2.5 % 97.5 %
#> pbo 19.15272 0.30047 18.56381 19.742
#> trt1 20.06258 0.33673 19.40260 20.723
#> trt2 20.59901 0.32085 19.97016 21.228
#>
#> Contrast : h_diff
#> Estimate Std.Err Z Value Pr(>|z|)
#> trt1 v.s. pbo 0.90986 0.45255 2.0105 0.0443766 *
#> trt2 v.s. pbo 1.44629 0.43915 3.2934 0.0009898 ***
#> trt2 v.s. trt1 0.53643 0.46552 1.1523 0.2491880
#> ---
#> 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 RobinCar2 supports
sp
for the simple randomization, pb
for the
permuted-block randomization, and ps
for the Pocock-Simon
randomization.
robin_glm(y_b ~ treatment * s1 + covar,
data = dummy_data,
treatment = treatment ~ ps(s1),
family = binomial(link = "logit")
)
#> Model : y_b ~ treatment * s1 + covar
#> Randomization: treatment ~ ps(s1) ( Pocock-Simon )
#> Variance Type: vcovG
#> Marginal Mean:
#> Estimate Std.Err 2.5 % 97.5 %
#> pbo 0.356097 0.033599 0.290243 0.4219
#> trt1 0.580696 0.034418 0.513238 0.6482
#> trt2 0.621386 0.034019 0.554711 0.6881
#>
#> Contrast : h_diff
#> Estimate Std.Err Z Value Pr(>|z|)
#> trt1 v.s. pbo 0.224599 0.047711 4.7075 2.508e-06 ***
#> trt2 v.s. pbo 0.265290 0.047534 5.5810 2.391e-08 ***
#> trt2 v.s. trt1 0.040691 0.047941 0.8488 0.396
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1