Introduction

The goal of this vignette is to show most of all possibilies with aVT (for aVirtualTwins meaning adaptation of Virtual Twins method) package.

VT method (Jared Foster and al, 2011) has been created to find subgroup of patients with enhanced treatment effect, if it exists. Theorically, this method can be used for binary and continous outcome. This package only deals with binary outcome in a two arms clinical trial.

VT method is based on random forests and regression/classification trees.

I decided to use a simulated dataset called sepsis in order to show how aVT package can be used. Type ?sepsis to know more about this dataset. Anyway, the true subgroup is PRAPACHE <= 26 & AGE <= 49.80.

NOTE: This true subgroup is defined with the lower event rate (survival = 1) in treatement arm. Therefore in following examples we’ll search the subgroup with the highest event rate, and we know it is PRAPACHE > 26 & AGE > 49.80.

Quick preview

## Dataset

Data used in VT are modelized by \(\left\{Y, T, X_1, \ldots, X_{p-2}\right\}\). \(p\) is the number of variables.

• \(Y\) is a binary outcome. In R, \(Y\) is a factor. Second level of this factor will be the desirable event. (\(Y=1\))
• \(T\) is treatment variable, \(T=1\) means active treatement, \(T=0\) means control treatment. In R, \(T\) is numeric.
• \(X_i\) is covariables, \(X_i\) can be categorical, continous, binary.

NOTE: if you run VT with interactions, categorical covariables must be transformed into binary variables.

Type ?formatRCTDataset for details.

Related functions/classes in aVirtualTwins package : VT.object(), vt.data(), formatRCTDataset.

Sepsis dataset

See Introduction.

Sepsis dataset is a simulated clinical trial with two groups treatment about sepsis desease. See details. This dataset is taken from SIDES method

Sepsis contains simulated data on 470 subjects with a binary outcome survival, that stores survival status for patient after 28 days of treatment, value of 1 for subjects who died after 28 days and 0 otherwise. There are 11 covariates, listed below, all of which are numerical variables.

Note that contrary to the original dataset used in SIDES, missing values have been imputed by random forest randomForest::rfImpute(). See file data-raw/sepsis.R for more details.

True subgroup is PRAPACHE <= 26 & AGE <= 49.80. NOTE: This subgroup is defined with the lower event rate (survival = 1) in treatement arm.

470 patients and 13 variables:

• survival : binary outcome
• THERAPY : 1 for active treatment, 0 for control treatment
• TIMFIRST : Time from first sepsis-organ fail to start drug
• AGE : Patient age in years
• BLLPLAT : Baseline local platelets
• blSOFA : Sum of baselin sofa (cardiovascular, hematology, hepaticrenal, and respiration scores)
• BLLCREAT : Base creatinine
• ORGANNUM : Number of baseline organ failures
• PRAPACHE : Pre-infusion apache-ii score
• BLGCS : Base GLASGOW coma scale score
• BLIL6 : Baseline serum IL-6 concentration
• BLADL : Baseline activity of daily living score
• BLLBILI : Baseline local bilirubin

Source: http://biopharmnet.com/subgroup-analysis-software/

Create object VirtualTwins

In order to begin the two steps of VT method, aVirtualTwins package needs to be initialized with vt.data() function. type ?vt.data for more details.

NOTE: if running VT with interactions between \(T\) and \(X\), set interactions = TRUE.

Code of vt.data() :

vt.data <- function(dataset, outcome.field, treatment.field, interactions = TRUE, ...){
  data <- formatRCTDataset(dataset, outcome.field, treatment.field, interactions = TRUE)
  VT.object(data = data, ...)
}

Example with Sepsis

# load library VT
library(aVirtualTwins)
# load data sepsis
data(sepsis)
# initialize VT.object
vt.o <- vt.data(sepsis, "survival", "THERAPY", TRUE)
## "1" will be the favorable outcome

1 will be the favorable outcome because 1 is the second level of "survival" column. It means that \(P(Y=1)\) is the probability of interest. Anyway, it’s still possible to compute \(P(Y=0)\).

Quick example

Sepsis does not have any categorical variable, following example show how vt.data deals with categorical values depending on interactions parameter

# Creation of categorical variable
cat.x <- rep(1:5, (nrow(sepsis))/5)
cat.x <- as.factor(cat.x)
sepsis.tmp <- cbind(sepsis, cat.x)
vt.o.tmp <- vt.data(sepsis.tmp, "survival", "THERAPY", TRUE)
## "1" will be the favorable outcome 
## Creation of dummy variables for cat.x 
## Dummy variable cat.x_1 created 
## Dummy variable cat.x_2 created 
## Dummy variable cat.x_3 created 
## Dummy variable cat.x_4 created 
## Dummy variable cat.x_5 created

Dummies variables are created for each category of cat.x variable. And cat.x is removed from dataset.

Step 1 : compute \(\hat{P_{1i}}\) and \(\hat{P_{0i}}\)

As described earlier, step 1 can be done via differents ways

# use randomForest::randomForest()
library(randomForest, verbose = F)
## randomForest 4.7-1.2
## Type rfNews() to see new features/changes/bug fixes.
# Reproducibility
set.seed(123)
# Fit rf model 
# default params
# set interactions to TRUE if using interaction between T and X
model.rf <- randomForest(x = vt.o$getX(interactions = T),
                         y = vt.o$getY(),
                         ntree = 500)
# initialize VT.forest.one
vt.f.rf <- vt.forest("one", vt.data = vt.o, model = model.rf, interactions = T)
### or you can use randomForest inside vt.forest()
vt.f.rf <- vt.forest("one", vt.data = vt.o, interactions = T, ntree = 500)

# # use randomForest::randomForest()
# library(party, verbose = F)
# # Reproducibility
# set.seed(123)
# # Fit cforest model 
# # default params
# # set interactions to TRUE if using interaction between T and X
# model.cf <- cforest(formula = vt.o$getFormula(), data = vt.o$getData(interactions = T))
# # initialize VT.forest.one
# vt.f.cf <- vt.forest("one", vt.data = vt.o, model = model.cf)

# Copy new object
vt.o.tr <- vt.o$copy()
# Change levels
tmp <- ifelse(vt.o.tr$data$survival == 1, "y", "n")
vt.o.tr$data$survival <- as.factor(tmp)
rm(tmp)
# Check new data to be sure
formatRCTDataset(vt.o.tr$data, "survival", "THERAPY")
## "y" will be the favorable outcome
# use caret::train()
library(caret, verbose = F)
## Loading required package: ggplot2
## 
## Attaching package: 'ggplot2'
## The following object is masked from 'package:randomForest':
## 
##     margin
## Loading required package: lattice
# Reproducibility
set.seed(123)
# fit train model
fitControl <- trainControl(classProbs = T, method = "none")
model.tr <- train(x = vt.o.tr$getX(interactions = T),
                  y = vt.o.tr$getY(),
                  method = "rf",
                  tuneGrid = data.frame(mtry = 5),
                  trControl = fitControl)
# initialize VT.forest.one
vt.f.tr <- vt.forest("one", vt.o.tr, model = model.tr)

# grow RF for T = 1
model.rf.trt1 <- randomForest(x = vt.o$getX(trt = 1),
                              y = vt.o$getY(trt = 1))
# grow RF for T = 0
model.rf.trt0 <- randomForest(x = vt.o$getX(trt = 0),
                              y = vt.o$getY(trt = 0))
# initialize VT.forest.double()
vt.doublef.rf <- vt.forest("double",
                           vt.data = vt.o, 
                           model_trt1 = model.rf.trt1, 
                           model_trt0 = model.rf.trt0)
### Or you can use randomForest() inside
vt.doublef.rf <- vt.forest("double",
                           vt.data = vt.o,
                           ntree = 200)

# initialize k-fold RF
# you can use randomForest options
model.fold <- vt.forest("fold", vt.data = vt.o, fold = 5, ratio = 1, interactions = T, ntree = 200)

# you get twin1 and twin2 by your own method
# here, i'll use random number between 0 and 1 :
twin1_random <- runif(470)
twin2_random <- runif(470)

# then you can initialize VT.difft class : 
model.difft <- VT.difft(vt.o, twin1 = twin1_random, twin2 = twin2_random, "absolute")
# compute difference of twins : 
model.difft$computeDifft()
# See results
head(model.difft$difft)
## [1]  0.73871597  0.26876747  0.45043814 -0.04273889 -0.03562053  0.05085339
# Graph :
# hist(model.difft$difft)

Step 2 : Estimate a Regression or Classification Tree

As described in the method, we define \(Z_i = \hat{P_{1i}} - \hat{P_{0i}}\). It’s the difference in term of response of the active treatments compared to the control treatment. The idea is to try to explain this difference by few covariables.

# initialize classification tree
tr.class <- vt.tree("class",
                    vt.difft = vt.f.rf,
                    sens = ">",
                    threshold = quantile(vt.f.rf$difft, seq(.5, .8, .1)),
                    maxdepth = 3,
                    cp = 0,
                    maxcompete = 2) 
# tr.class is a list if threshold is a vectoor
class(tr.class)
## [1] "list"
# acce trees with treeXX
class(tr.class$tree1)
## [1] "VT.tree.class"
## attr(,"package")
## [1] "aVirtualTwins"

# initialize regression tree
tr.reg <- vt.tree("reg",
                  vt.difft = vt.f.rf,
                  sens = ">",
                  threshold = quantile(vt.f.rf$difft, seq(.5, .8, .1)))
# tr.class is a list if threshold is a vectoor
class(tr.reg)
## [1] "list"
# access trees with treeXX
class(tr.reg$tree1)
## [1] "VT.tree.reg"
## attr(,"package")
## [1] "aVirtualTwins"

Subgroups and results

Once trees have been computed, you surely want to see what are the subgroups. This package provides a wrapper function of intern methods of VT.tree class : vt.subgroups().

This function takes in argument :

• vt.tree : object or list of class VT.tree. Return of the vt.tree() function.
• only.leaf : logical. Set TRUE (default) to visualize only terminal nodes.
• only.fav : logical. Set TRUE (default) to visualize only class 1 nodes. (\(\hat{A}\))
• tables : logical. Set FALSE (default) to prevent tables of incidences from being printed.
• verbose : logical. Set FALSE (default) to prevent detailed stuffs from being printed.
• compete : logical. Set TRUE to print competitors rules thanks to competitors split. FALSE is default value.

If vt.tree is a list, unique subgroups are printed.

# use tr.class computed previously
vt.sbgrps <- vt.subgroups(tr.class)
# print tables with knitr package
library(knitr)
knitr::kable(vt.sbgrps)

You can plot one tree with package rpart.plot

## Loading required package: rpart
## Warning: Cannot retrieve the data used to build the model (so cannot determine roundint and is.binary for the variables).
## To silence this warning:
##     Call rpart.plot with roundint=FALSE,
##     or rebuild the rpart model with model=TRUE.

If you want to see competitors split :

tr.class$tree2$createCompetitors()
head(tr.class$tree2$competitors)
##    count ncat     improve   index      var path          string
## 1    470   -1 123.0361124 26.5000 PRAPACHE    2 PRAPACHE < 26.5
## 2    470    1  30.3361464 11.5000    BLGCS    2   BLGCS >= 11.5
## 3    470   -1  29.0631063 61.6710      AGE    2    AGE < 61.671
## 9    157   -1   6.6549865 51.7395    AGE.1    6   AGE.1 < 51.74
## 10   157   -1   1.7771209  1.5000 ORGANNUM    6  ORGANNUM < 1.5
## 11   157    1   0.8252186  0.6500  BLLBILI    6 BLLBILI >= 0.65

If you want to print incidence of a subgroup :

vt.o$getIncidences("PRAPACHE >= 26 & AGE >= 52")
## $table.selected
## $table.selected$table
##            trt
## resp              0       1     sum
##   0          30.000  18.000  48.000
##   1          15.000  72.000  87.000
##   sum        45.000  90.000 135.000
##   Incidence   0.333   0.800   0.644
## 
## $table.selected$rr
## [1] 2.402402
## 
## 
## $table.not.selected
## $table.not.selected$table
##            trt
## resp              0       1     sum
##   0          71.000 170.000 241.000
##   1          37.000  57.000  94.000
##   sum       108.000 227.000 335.000
##   Incidence   0.343   0.251   0.281
## 
## $table.not.selected$rr
## [1] 0.7317784
# or
# tr.class$tree2$getIncidences("PRAPACHE >= 26 & AGE >= 52")

If you want to get infos about the tree

tr.class$tree2$getInfos()
## 
## Threshold = 0.0747
## Delta = 0.0671
## Sens : >
## Size of Ahat : 144
# access Ahat
# tr.class$tree2$Ahat

You can re-run rpart computation:

tr.class$tree2$run(maxdepth = 2)

Type ?VT.tree for details.

Questions

This vignette is a bit messy right now, therefore feel free to ask anything to the repository issue reports :

1. https://github.com/prise6/aVirtualTwins
2. read documentation with ?aVirtualTwins