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 1Data Introduction
In the dummy_data, we have the following columns:
-
idis the patient identifier. -
treatmentis the treatment assignment. -
s1is the first stratification factor. -
s2is the second stratification factor. -
covaris the continuous covariate. -
yis the continuous outcome. -
y_bis 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
#>
#> Marginal Mean Variance:
#> pbo trt1 trt2
#> 0.06768998 0.07592944 0.07654319
#>
#>
#> Variance Type: vcovG
#> 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 ' ' 1We 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
#>
#> Marginal Mean Variance:
#> pbo trt1 trt2
#> 0.06693369 0.07540157 0.07688832
#>
#>
#> 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 ' ' 1Obtain 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
#>
#> Marginal Mean Variance:
#> pbo trt1 trt2
#> 0.03359913 0.03441801 0.03401864
#>
#>
#> Variance Type: vcovG
#> 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 ' ' 1Obtain 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
#>
#> Marginal Mean Variance:
#> pbo trt1 trt2
#> 0.3004728 0.3367318 0.3208476
#>
#>
#> Variance Type: vcovG
#> 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 ' ' 1Using 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
#>
#> Marginal Mean Variance:
#> pbo trt1 trt2
#> 0.03359913 0.03441801 0.03401864
#>
#>
#> Variance Type: vcovG
#> 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