inst/doc/iterativeBMAsurv.R

In iterativeBMAsurv: The Iterative Bayesian Model Averaging (BMA) Algorithm For Survival Analysis

### R code from vignette source 'iterativeBMAsurv.Rnw'

###################################################
### code chunk number 1: Setup
###################################################
library(BMA)
library(iterativeBMAsurv)


###################################################
### code chunk number 2: getTrainData
###################################################
## Use the sample training data. The data matrix is called trainData.
data(trainData)
## The survival time vector for the training set is called trainSurv, where survival times are reported in years.
data(trainSurv)
## The censor vector for the training set is called trainCens, where 0 = censored and 1 = uncensored.
data(trainCens)


###################################################
### code chunk number 3: trainingStep
###################################################
## Training phase: select relevant genes
## In this example training set, the top 100 genes have already been sorted in decreasing order of their log likelihood
ret.list <- iterateBMAsurv.train.wrapper (x=trainData, surv.time=trainSurv, cens.vec=trainCens, nbest=5)

## Extract the {\tt bic.surv} object
ret.bic.surv <- ret.list$obj

## Extract the names of the genes from the last iteration of {\tt bic.surv}
gene.names <- ret.list$curr.names

## Get the selected genes with probne0 > 0
top.gene.names <- gene.names[ret.bic.surv$probne0 > 0]
top.gene.names

## Get the posterior probabilities for the selected models
ret.bic.surv$postprob


###################################################
### code chunk number 4: testStep
###################################################
## The test data matrix is called testData.
data(testData)
## The survival time vector for the test set is called testSurv, where survival times are reported in years
data(testSurv)
## The censor vector for the test set is called testCens, where 0 = censored and 1 = uncensored
data(testCens)

## Get the subset of test data with the genes from the last iteration of bic.surv
curr.test.dat <- testData[, top.gene.names]

## Compute the predicted risk scores for the test samples
y.pred.test <- apply (curr.test.dat, 1, predictBicSurv, postprob.vec=ret.bic.surv$postprob, mle.mat=ret.bic.surv$mle)

## Compute the risk scores for the training samples
y.pred.train <- apply (trainData[, top.gene.names], 1, predictBicSurv, postprob.vec=ret.bic.surv$postprob, mle.mat=ret.bic.surv$mle)

## Assign risk categories for test samples
## Argument {\tt cutPoint} is the percentage cutoff for separating the high-risk group from the low-risk group
ret.table <- predictiveAssessCategory (y.pred.test, y.pred.train, testCens, cutPoint=50) 

## Extract risk group vector and risk group table
risk.vector <- ret.table$groups
risk.table <- ret.table$assign.risk
risk.table

## Create a survival object from the test set
mySurv.obj <- Surv(testSurv, testCens)

## Extract statistics including p-value, chi-square, and variance matrix
stats <- survdiff(mySurv.obj ~ unlist(risk.vector))
stats


###################################################
### code chunk number 5: trainPredictTestStep
###################################################
## Use p=10 genes and nbest=5 for fast computation
ret.bma <- iterateBMAsurv.train.predict.assess (train.dat=trainData, test.dat=testData, surv.time.train=trainSurv, surv.time.test=testSurv, cens.vec.train=trainCens, cens.vec.test=testCens, p=10, nbest=5)

## Extract the statistics from this survival analysis run
number.genes <- ret.bma$nvar
number.models <- ret.bma$nmodel
evaluate.success <- ret.bma$statistics
evaluate.success


###################################################
### code chunk number 6: crossValidationStep
###################################################
## Perform 1 run of 2-fold cross validation on the training set, using p=10 genes and nbest=5 for fast computation
## Argument {\tt diseaseType} specifies the type of disease present the training samples, used for writing to file
cv <- crossVal (exset=trainData, survTime=trainSurv, censor=trainCens, diseaseType="DLBCL", noRuns=1, noFolds=2, p=10, nbest=5)


###################################################
### code chunk number 7: imageplot
###################################################
 imageplot.iterate.bma.surv (ret.bic.surv)