The goal of this vignette is to show most of all possibilies with VT (for VirtualTwins) package.
+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 VT package can be used. Type ?sepsis to know more about this dataset. Anyway, the true subgroup is PRAPACHE <= 26 & AGE <= 49.80.
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.
NOTE: if you run VT with interactions, categorical covariables must be transformed into binary variables.
Type ?formatRCTDataset for details.
Related functions/classes in VirtualTwins package : VT.object(), vt.data(), formatRCTDataset.
Related functions/classes in aVirtualTwins package : VT.object(), vt.data(), formatRCTDataset.
From one of these methods you can estimate \(\hat{P_{1i}}\) and \(\hat{P_{0i}}\).
-Related functions/classes in VirtualTwins package : VT.difft(), VT.forest(), VT.forest.one(), VT.forest.double(), VT.forest.fold().
Related functions/classes in aVirtualTwins package : VT.difft(), VT.forest(), VT.forest.one(), VT.forest.double(), VT.forest.fold().
The idea is to identify which covariable from \(X\) described variation of \(Z\).
-Related functions/classes in VirtualTwins package : VT.tree(), VT.tree.class(), VT.tree.reg().
Related functions/classes in aVirtualTwins package : VT.tree(), VT.tree.class(), VT.tree.reg().
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 outcomeTHERAPY : 1 for active treatment, 0 for control treatmentTIMFIRST : Time from first sepsis-organ fail to start drugAGE : Patient age in yearsBLLPLAT : Baseline local plateletsblSOFA : Sum of baselin sofa (cardiovascular, hematology, hepaticrenal, and respiration scores)BLLCREAT : Base creatinineORGANNUM : Number of baseline organ failuresPRAPACHE : Pre-infusion apache-ii scoreBLGCS : Base GLASGOW coma scale scoreBLIL6 : Baseline serum IL-6 concentrationBLADL : Baseline activity of daily living scoreBLLBILI : Baseline local bilirubinSource: http://biopharmnet.com/wiki/Software_for_subgroup_identification_and_analysis
+In order to begin the two steps of VT method, VirtualTwins package need to be initialized with vt.data() function. type ?vt.data for more details.
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, ...){
@@ -145,18 +172,27 @@ vt.o.tmp <- vt.data(sepsis.tm
## Dummy variable cat.x_4 created
## Dummy variable cat.x_5 createdDummies variables are created for each category of cat.x variable. And cat.x is removed from dataset.
As described earlier, step 1 can be done via differents ways
-Following example used sepsis data created in previous part.
To perform simple random forest on VT.object, randomForest, caret and party package can be used.
Class VT.forest.one is used. It takes in arguments :
vt.object : return of vt.data() functionmodel : A random forest modelinteractions : logical, TRUE is default valuewith randomForest
# use randomForest::randomForest()
-library(randomForest, verbose = F)
-# Reproducibility
+library(randomForest, verbose = F)## randomForest 4.6-10
+## Type rfNews() to see new features/changes/bug fixes.# Reproducibility
set.seed(123)
# Fit rf model
# default params
@@ -164,7 +200,9 @@ vt.o.tmp <- vt.data(sepsis.tm
model.rf <- randomForest(x = vt.o$getX(interactions = T),
y = vt.o$getY())
# initialize VT.forest.one
-vt.f.rf <- VT.forest.one(vt.o, model.rf)with party
cforest() can be usefull however computing time is really long. I think there is an issue when giving cforest object in Reference Class parameter. Need to fix it.
# # use randomForest::randomForest()
@@ -176,7 +214,9 @@ vt.f.rf <- VT.forest.one(vt.o
# # 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.o, model.cf)with caret
Using caret can be usefull to deal with parallel computing for example.
NOTE: For caret levels of outcome can’t be 0, so i’ll change levels name into “n”/“y”
## "y" will be the favorable outcome# use caret::train()
-library(caret, verbose = F)
-# Reproducibility
+library(caret, verbose = F)## Loading required package: lattice
+## Loading required package: ggplot2# Reproducibility
set.seed(123)
# fit train model
fitControl <- trainControl(classProbs = T, method = "none")
@@ -201,12 +243,125 @@ model.tr <- train(tuneGrid = data.frame(mtry = 5),
trControl = fitControl)
# initialize VT.forest.one
-vt.f.tr <- VT.forest.one(vt.o.tr, model.tr)To perform double random forest on VT.object, same packages as simple random forest can be used.
Class VT.forest.double is used. It takes in arguments :
vt.object : return of vt.data() functionmodel_trt1 : a random forest model for \(T=1\)model_trt0 : a random forest model for \(T=0\)NOTE: use trt parameter in VT.object::getX() or VT.object::getY() methods to obtain part of data depending on treatment. See following example.
with randomForest
# 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.o, model.rf.trt1, model.rf.trt0)
+# Then, use run() to compute probabilities
+vt.doublef.rf$run()Follow the same structure for caret or cforest models.
This idea is taken from method 3 of Jared Foster paper :
+++A modification of [previous methods] is to obtain \(\hat{P_{1i}}\) and \(\hat{P_{0i}}\) via cross-validation. In this méthod the specific data for subject \(i\) is not used to obtain \(\hat{P_{1i}}\) and \(\hat{P_{0i}}\). Using k-fold cross-validation, we apply random forest regression approach to \(\frac{k-1}{k}\) of the data and use the resulting predictor to obtain estimates of \(P_{1i}\) and \(P_{0i}\) for the remaining \(\frac{1}{k}\) of the observations. This is repeated \(k\) times.
+
To use this approach, type VirtualTwins:::VT.forest.fold(). This class takes in argument :
vt.object : return of vt.data() functionfold : number of fold (e.g. \(5\))ratio : Control of sampsize balance. ratio of \(2\) means that there 2 times le highest level compared to the other. “Highest” means the level with larger observations. It’s in test.interactions : Logical. If TRUE, interactions between covariables and treatments are used. FALSE otherwise.NOTE: This function use only randomForest package.
# initialize k-fold RF
+model.fold <- VirtualTwins:::VT.forest.fold(vt.o, fold = 5, ratio = 1, interactions = T)
+# grow RF with randomForest package options
+# set do.trace option to see the 5 folds
+model.fold$run(ntree = 500, do.trace = 500)## ntree OOB 1 2
+## 500: 32.38% 16.46% 57.72%
+## ntree OOB 1 2
+## 500: 33.25% 18.50% 55.26%
+## ntree OOB 1 2
+## 500: 30.03% 16.52% 50.34%
+## ntree OOB 1 2
+## 500: 29.26% 13.56% 55.71%
+## ntree OOB 1 2
+## 500: 30.60% 13.79% 59.70%Random Forests are not the only models you can use to compute \(\hat{P_{1i}}\) and \(\hat{P_{0i}}\). Any prediction model can be used, as logitic regression, boosting …
+Anyway, aVirtualTwins package can be used. To do so, you can use VT.difft() class. It is important to note this the parent class of all “forests” classes. It takes in argument :
vt.object : return of vt.data() functiontwin1 : estimate of \(P(Y_{i} = 1 | T = T_{i})\) : meaning response probability under the correct treatment.twin1 : estimate of \(P(Y_{i} = 1 | T = 1-T_{i})\) : meaning response probability under the other treatment.method : absolute (default), relative or logit. See ?VT.difft() for details.# 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.2507360 -0.3905591 0.3102285 0.3986582 0.3747078 -0.2508240# Graph :
+# hist(model.difft$difft)NOTE: Also, you can clone repository, write your own child class of VT.difft() AND submit it !
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.
+We define a new variable \(Z^{*}\), \(Z^{*}_i=1\) if \(Z_i > c\) and \(Z^{*}_i=0\) otherwise. Classification tree’s goal is to explain the value \(Z^*=1\). \(c\) is a threshold given by the user. It’s the threshold for which the difference is “interesting”. One idea is to use quantiles of the difft distribution.
+To compute a classifiction tree, VT.tree.class() is used. Internally, rpart::rpart() is computed. It takes in argument:
vt.difft : VT.difft objectthreshold : \(0.05\) by default, corresponds to \(c\)sens : c(">", "<"). sens corresponds to the way \(Z^{*}\) is defined.
+">" (default) : \(Z^{*}\), \(Z^{*}_i=1\) if \(Z_i > c\) and \(Z^{*}_i=0\) otherwise."<" : \(Z^{*}\), \(Z^{*}_i=1\) if \(Z_i < c\) and \(Z^{*}_i=0\) otherwise.screening : NULL is default value. If TRUE only covariables in varimp vt.object ’s field is used.See ?VT.tree for details.
# initialize classification tree
+tr.class <- VT.tree.class(vt.f.rf, sens = ">", threshold = quantile(vt.f.rf$difft, 0.7))
+# compute tree with rpart option
+tr.class$run(cp = 0, maxdepth = 3, maxcompete = 1)Use regression tree to explain \(Z\) by covariables \(X\). Then some leafs have predicted \(Z_i\) greater than the threshold \(c\) (if \(sens\) is “>”), and it defines which covariables explain \(Z\).
+The class to use is VT.tree.reg(). It takes same parameters than classification mehod.
# initialize regression tree
+tr.reg <- VT.tree.reg(vt.f.rf, sens = ">", threshold = quantile(vt.f.rf$difft, 0.7))
+# compute tree with rpart option
+tr.reg$run(cp = 0, maxdepth = 3, maxcompete = 1)