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.

Method

VT is a two steps method but with many possibilities

let \(\hat{P_{1i}} = P(Y_i = 1|T_i = 1, X_i)\)
let \(\hat{P_{0i}} = P(Y_i = 1|T_i = 0, X_i)\)
let \(X = \left\{X_1, \ldots, X_{p-2}\right\}\)

Simple Random Forest

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() function
• model : A random forest model
• interactions : logical, TRUE is default value

with randomForest

# use randomForest::randomForest()
library(randomForest, verbose = F)
# 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())
# initialize VT.forest.one
vt.f.rf <- VT.forest.one(vt.o, model.rf)
# Then, use run() to compute probabilities
vt.f.rf$run()

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()
# 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.o, model.cf)
# # Then, use run() to compute probabilities
# vt.f.cf$run()

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”

# 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)
# 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.tr)
# Then, use run() to compute probabilities
vt.f.tr$run()

Double Random Forest

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() function
• model_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.

K Fold Random Forest

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 aVirtualTwins:::VT.forest.fold(). This class takes in argument :

• vt.object : return of vt.data() function
• fold : 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 <- aVirtualTwins:::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%

Build Your Own Model

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() function
• twin1 : 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 !

---

Classification

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 object
• threshold : \(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)

Regression

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)

Subgroups and results