crossBreastCancer.R

##################################################################################################
#### THE SBC model in Cross Validation setting ####################################################
##### This file Runs a cross validation 5 times cross validation on the Breast Cancer Data Set ####
setwd("~/Dropbox/Code/DPmixturemodel/SBC")
rm(list =ls())
source('import.R')

### Load the data
data(vanDeVijver, package="seventyGeneData")
## Load the Pheno Data for VIJVER
pheno.vijver <- read.csv('/home/bit/ashar/ExpressionSets/ONE_VIEW/breastcancer/VIJVERpheno.csv')
## Load the 70 gene siganture
source('loadVIJVERsignature.R')
loadVIJVERsignature()

########## Preparing the Data Matrix ######################################
Y.full <- t(impute(exprs(vanDeVijver)))

######## Defining the Pheno type ###########
time.pre <- pheno.vijver$Follow_up_time_or_metastasis 
censoring.pre <- pheno.vijver$event_metastasis
### Define the Poor vs Good discrimination from VIJVER
c.vijver <- pheno.vijver$X70_genes
levels(c.vijver)[1:2] <- c(1:2)



#### Prepare for the parallel computation ###
library(parallel)
number.of.cores <- detectCores()
cluster.parallel <- makeCluster(2, type= "FORK")
library(doParallel)
registerDoParallel(cluster.parallel)
library(foreach)
getDoParWorkers()

############################# PARAMETERS for GIBB's SAMPLING ####
iter = 20
iter.burnin = 20
iter.thin  = 5
Nps = as.integer(iter/ iter.thin)

### Some parameters for the DP mixture model ###################
k = 2
F =k


#### Define key variables which store the results ####
# recov.logrank.vijver <- c(0)
# pred.logrank.vijver <- c(0)
# 
# recovCIndex.vv.pcox <- c(0)
# predCIndex.vv.pcox <- c(0)
# 
# predCIndex.vv.kk.pcox <- c(0)
# pred.logrank.vv.kk <- c(0)
# 
# 
# recov.logrank.kk <- c(0)
# pred.logrank.kk <- c(0)
# 
# 
# recovCIndex.kk.pcox <- c(0)
# predCIndex.kk.pcox <- c(0)
# 
# recov.logrank.sbc <- c(0)
# pred.logrank.sbc <- c(0)
# 
# recovCIndex.PC <- c(0)
# predCIndex.PC <- c(0)
# 
# 
# recovCIndex.sbc.matrix <- matrix(0, nrow =  Nps, ncol = 5)
# recovCIndex.sbc.paft.matrix <- matrix(0, nrow =  Nps, ncol = 5)
# 
# predCIndex.sbc.matrix <- matrix(0, nrow =  Nps, ncol = 5)
# predCIndex.sbc.pcox <- c(0)

# recovCIndex.PC <- c(0)
# predCIndex.PC <- c(0)

#signature.sbc.list <- list(0)



##Define the key functions to be used
source('rchinese.R')
source('priorPARAMETERS.R')
source('posteriorCensoredTime.R')
source('loadSBCvijversignature.R')
source('posteriorCLASS.R')
source('posteriorGMM.R')
source('posteriorAFT.R')
source('posteriorCensoredTime.R')
source('priorPARAMETERS.R')
source('calculateLIKELIHOOD.R')
source('posteriorhyperGMM.R')
source('posterioralpha.R')



#############################################################################################
########################## BEGIN THE CROSS-VALIDATION FOLD ###################################
##############################################################################################

## Set the initial conditions
j =1
set.seed(42*j)


#### Define the folds
folds <- createFolds(c.vijver, k = 2, list = TRUE, returnTrain = FALSE)
 
#### The actual cross-validation loop contains around 200 lines of code
### Begins here #####
results <- foreach(i = 1:1, .export = ls(),.packages = list.p) %dopar%{
test.index <- folds[[i]]
####### If we define the splits ourselves ####
############### Defining Our data without prefiltering #########################
Y.pre.train <- Y.full[-test.index, ]
Y.pre.test <-  Y.full[test.index,]
c.true <- c.vijver[-test.index]
c.true.new <- c.vijver[test.index]
time <- time.pre[-test.index]
censoring <- censoring.pre[-test.index]
time.new <- time.pre[test.index]
censoring.new <- censoring.pre[test.index]


### setting up the survival objects
smod <-  Surv(time, censoring)
smod.new <- Surv(time.new, censoring.new)


#### Use VIJVER's classification to come up with separability and c-Index##
#### Using VIJVER training and Testing Data Sets ##########################
Y.vijver.train <- Y.pre.train[,colnames(Y.pre.train) %in% vijver.signature]
Y.vijver.test <- Y.pre.test[,colnames(Y.pre.test) %in% vijver.signature]

###
recov.logrank.vijver<- unlist(survdiff(smod ~ c.true))$chisq
pred.logrank.vijver <-  unlist(survdiff(smod.new ~ c.true.new))$chisq

####### Recovering C-Indexes ##############
### Using L1 penalized Cox-PH ############
linear.vv.recovery <- c(0)
linear.vv.prediction <- c(0)
for ( q in 1:2){
  ind <- which(c.true == q)
  ind.new <- which(c.true.new == q)
  reg.cox <- cv.glmnet(x = Y.vijver.train[ind,], y = smod[ind], family = "cox", maxit = 10000000)
  linear.vv.recovery[ind] <- predict(object =reg.cox, newx = Y.vijver.train[ind,], s= "lambda.min")
  linear.vv.prediction[ind.new] <- predict(object =reg.cox, newx = Y.vijver.test[ind.new,], s= "lambda.min")
  }
recovCIndex.vv.pcox <- as.numeric(survConcordance(smod ~ linear.vv.recovery)[1])
predCIndex.vv.pcox <- as.numeric(survConcordance(smod.new ~ linear.vv.prediction)[1])


### Use k-NN to predict using the VIJVER classification
label.train <- as.factor(c.true)
### One has to to tune the k-NN classifier for k ###
fitControl <- trainControl(method = "repeatedcv", number = 5,repeats = 5, allowParallel = FALSE)
### Tune the parameter k
knnFit <- train(x = Y.vijver.train, y = label.train, method =  "knn", trControl = fitControl, tuneLength = 5)
knnPredict <- predict(knnFit,newdata = Y.vijver.test)
label.test <- knnPredict
linear.vv.kk.prediction <- c(0)
for ( q in 1:2){
  ind <- which(label.train == q)
  ind.new <- which(label.test == q)
  reg.cox <- cv.glmnet(x = Y.vijver.train[ind,], y = smod[ind], family = "cox", maxit = 10000000)
  linear.vv.kk.prediction[ind.new] <- predict(object =reg.cox, newx = Y.vijver.test[ind.new,], s= "lambda.min")
}
predCIndex.vv.kk.pcox<- as.numeric(survConcordance(smod.new ~ linear.vv.kk.prediction)[1])
pred.logrank.vv.kk <- unlist(survdiff(smod.new ~ label.test))$chisq




#### Use SBC signature ##
rel <- loadSBCvijversignature(Y.pre.train, time, censoring)
signature.sbc <- rel$signature.sbc


###### We get the signature and then we can define the folds #####
########## Creating Training and Test Data ##########################
Y <- Y.pre.train[,colnames(Y.pre.train) %in% signature.sbc ]
Y.new <- Y.pre.test[,colnames(Y.pre.test) %in% signature.sbc]

### Setting the important parameters
D <- ncol(Y)
N <- nrow(Y)
N.new <- nrow(Y.new)
K <-  as.integer(N/5)

####### Use HC clustering + k-NN prediction  #########################################
###### Use SBC signature        ##################################################
d <- dist(as.matrix(Y))
hc <- hclust(d)
label.train <- as.factor(cutree(hc, k = F))
#### Use k-means clustering
# gr.km <- kmeans(Y, F, nstart =10)
# label.train <- as.factor(gr.km$cluster)
### One has to to tune the k-NN classifier for k ###
fitControl <- trainControl(method = "repeatedcv", number = 5,repeats = 5, allowParallel = FALSE)
### Tune the parameter k
knnFit <- train(x = Y, y = label.train, method =  "knn", trControl = fitControl, tuneLength = 5)
knnPredict <- predict(knnFit,newdata = Y.new )
label.test <- knnPredict

recov.logrank.kk <- unlist(survdiff(smod ~ label.train))$chisq
pred.logrank.kk <- unlist(survdiff(smod.new ~ label.test))$chisq



###### penCox on top of the clustering ###################################################################
######## Penalized Cox PH with k-means clustering###########################################
linear.kk.recovery <- c(0)
linear.kk.prediction <- c(0)
for ( q in 1:2){
  ind <- which(label.train == q)
  ind.new <- which(label.test == q)
  reg.cox <- cv.glmnet(x = Y[ind,], y = smod[ind], family = "cox", maxit = 10000000)
  linear.kk.recovery[ind] <- predict(object =reg.cox, newx = Y[ind,], s= "lambda.min")
  linear.kk.prediction[ind.new] <- predict(object =reg.cox, newx = Y.new[ind.new,], s= "lambda.min")
}
recovCIndex.kk.pcox <- as.numeric(survConcordance(smod ~ linear.kk.recovery)[1])
predCIndex.kk.pcox  <- as.numeric(survConcordance(smod.new ~ linear.kk.prediction)[1])


# ######################### Initialize the Parameters ##############################
Time <- cbind(time, censoring)
## HYPER PRIORS
## Hyper parameters of the DP
shape.alpha <- 2
rate.alpha <- 1
## Hyperparameters for the GMM
beta  = D+1
ro = 0.5



alpha  = rgamma(1, shape = shape.alpha, rate = rate.alpha )
c <-  rchinese(N,alpha)
f <- table(factor(c, levels = 1:max(c)))

## Empirical Bayes Estimate of the Hyperparameters
epsilon = as.vector(apply(Y,2,mean))
W = diag(diag(cov(Y)))


## Initialization of the parameters for Gaussian Mixture
mu = matrix(data = NA, nrow = K, ncol = D)
S = array(data = NA, dim =c(K,D,D))


#Sparsity controlling hyperparameter of the BAYESIAN LASSO MODEL
r =1
si = 1.78


##Actual parameters
lambda2 <- numeric(K)
tau2 = matrix(data = NA, nrow = K, ncol = D)
betahat = matrix(data = NA, nrow = K, ncol = D)
sigma2 <- rep(NA, K)
beta0 <- rep(NA, K)
That <-  c(0)

## Fitting a linear model to the whole model
Ysc <- scale(Y[1:N,1:D], center = TRUE, scale =TRUE)
lm.data <- lm(time ~ Ysc)
sig2.dat <-  var(lm.data$residuals)


## Set Some Initial Values for the Cluster Parameters


disclass <- table(factor(c, levels = 1:K))
activeclass <- which(disclass!=0)
for ( j in 1:length(activeclass)){

  priorone <- priordraw(beta, W, epsilon, ro, r, si, N, D, sig2.dat)
  mu[activeclass[j],] <- (priorone$mu)
  S[activeclass[j],1:D,1:D]  <- priorone$Sigma
  beta0[activeclass[j]] <- priorone$beta0
  sigma2[activeclass[j]] <- priorone$sigma2
  betahat[activeclass[j],1:D] <- priorone$betahat
  lambda2[activeclass[j]] <- priorone$lambda2
  tau2[activeclass[j], 1:D] <- priorone$tau2
}

# The Time has to be initialized
ti <- updatetime(c, Y, Time,That, beta0, betahat, sigma2)
That <- ti$time


################# K-Means BLASSO INITIALIZATION ############################################
G <- F

k.data <- kmeans(Y,F,nstart =10)

c <- k.data$cluster

c.kmeans <- c
#### Under special cases
###### c <- c.true


prior.numclust <- table(factor(c, levels = 1:K))
prior.activeclass <- which(prior.numclust!=0)

### The means are set using the k-means
for ( i in 1:length(prior.activeclass)){
  mu[prior.activeclass[i],1:D] <-  k.data$centers[i,1:D]

  S[prior.activeclass[i],1:D,1:D] <-  priordraw(beta, W, epsilon, ro, r, si,N,D, sig2.dat)$Sigma

  lclust <- which(c == prior.activeclass[i])

  reg.blas <- 0

  sum <- c(0)

  coeff <- 0

  Ytemp <-  matrix(NA, nrow = length(lclust), ncol = D)

  Ytemp <- scale(Y[lclust,1:D], center = TRUE, scale = TRUE)


  ### Part where I use the MONOMVN PACKAGE

  Ttemp <- as.vector(That[lclust])

  ntemp <- length(lclust)

  reg.blas <- blasso(Ytemp, Ttemp, T = 300,thin =  50, RJ = TRUE, mprior = 0.0 ,normalize = TRUE, verb = 0)

  sum <- summary(reg.blas, burnin= 100)

  ## Selecting those features which are relevant

  coeff <- unlist(lapply(strsplit(sum$coef[3,], split = ":"), function(x) as.numeric(unlist(x)[2])))

  beta0[prior.activeclass[1]] <- coeff[1]

  indexplusone <- D+1

  ind <- 2:indexplusone

  betahat[prior.activeclass[i], ] <- coeff[ind]

  ta <- unlist(lapply(strsplit(sum$tau2i[3,], split = ":"), function(x) as.numeric(unlist(x)[2])))

  tau2[prior.activeclass[i],] <- ta

  sigma2[prior.activeclass[i]] <- sum$s2[3]

  lambda2[prior.activeclass[i]] <- sum$lambda2[3]

}

## Deleting those values which are no longer relevant
g <- table(factor(c, levels = 1:K))
inactive <- which(g==0)

for ( i in 1:length(inactive)){
  mu[inactive[i],1:D]  <- NA
  S[inactive[i],1:D,1:D]  <- NA
  beta0[inactive[i]] <- NA
  sigma2[inactive[i]] <- NA
  betahat[inactive[i],1:D] <- NA
  lambda2[inactive[i]] <- NA
  tau2[inactive[i], 1:D] <- NA
}

############# Train the Model #########################################

iter.burnin = iter.burnin
cognate <- NA
param <- NA
paramtime <- NA
loglike<- rep(0, iter)
timeparam <- NA
time.predicted <- c(0)
cindex <- c(0)


likli <- c(0)
gmm.likli <- c(0)
aft.likli <- c(0)


cog <- loglikelihood(c,Y,mu,S,alpha,That, beta0, betahat, sigma2, lambda2, tau2, K, epsilon, W, beta, ro,D, r, si, Time,N, sig2.dat)
likli[1] <- cog$loglikelihood
gmm.likli[1] <- cog$GMMlikelihood
aft.likli[1] <- cog$AFTlikelihood

rmse <- c(0)
randy <- c(0)



o =1
#################### BURNIN PHASE ###################################################
o.iter = o

for (o in o.iter:iter.burnin) {


  ################## PARAMETERS OF THE DP Mixture Model ######################################################
  ## Updating the parameters based on the observations
  param <- posteriorGMMparametrs(c,Y,mu,S, alpha,K, epsilon, W, beta, ro,N,D )
  mu <- param$mean
  S <- param$precision
  paramtime <- posteriortimeparameters(c, That, lambda2,tau2,sigma2,beta0, betahat, Y, K, epsilon, W, beta, ro,D, r, si, Time,N, sig2.data)
  beta0 <- paramtime$beta0
  betahat <- paramtime$betahat
  sigma2 <- paramtime$sigma2
  lambda2 <- paramtime$lambda2
  tau2 <- paramtime$tau2

  ########################## THE HYPERPARAMETERS OF THE GMM #################################
  #  Updating the hyper paramters
  hypercognate <- posteriorhyperPLUS (c, Y, mu, S, epsilon, W, beta, ro )
  epsilon <- hypercognate$epsilon
  tmpW <- hypercognate$W
  W <- matrix(as.matrix(tmpW),nrow = D, ncol =D)
  ro <- hypercognate$ro

  #   if( o%%10 == 0){
  #     res <- try(posteriorbeta(c, beta, D, S, W))
  #     if (class(res) == "try-error"){
  #       beta = beta
  #     } else{
  #       beta <- posteriorbeta(c, beta, D, S, W)
  #
  #     }
  #   }
  #
  ################# INDICATOR VARIABLE ##################################################################
  ## Updating the indicator variables and the parameters
  cognate <- posteriorchineseAFT(c,Y,mu,S,alpha,That, beta0, betahat, sigma2, lambda2, tau2, K, epsilon, W, beta, ro,D, r, si, Time,N, sig2.dat)
  c <- cognate$indicator
  mu <- cognate$mean
  S <- cognate$precision
  beta0 <- cognate$beta0
  betahat <- cognate$betahat
  sigma2 <- cognate$sigma2
  lambda2 <- cognate$lambda2
  tau2 <- cognate$tau2

  ########################### The Concentration Parameter #################################################################
  # Updating the concentration parameter
  alpha <- posterioralpha(c, N, alpha, shape.alpha, rate.alpha)

  ####################### The Censored Times ###########################################################
  # Updating the Time Variable
  ti <- NA
  ti <- updatetime(c, Y, Time,That, beta0, betahat, sigma2)
  That <- ti$time


  ##################### Print SOME Statistics #####################################################
  randy[o] <- adjustedRandIndex(c.kmeans,as.factor(c))
  print(randy[o])
  cog <- loglikelihood(c,Y,mu,S,alpha,That, beta0, betahat, sigma2, lambda2, tau2, K, epsilon, W, beta, ro,D, r, si, Time,N, sig2.dat)
  likli[o+1] <- cog$loglikelihood
  gmm.likli[o+1] <- cog$GMMlikelihood
  aft.likli[o+1] <- cog$AFTlikelihood

  print(likli[o+1])
  print(gmm.likli[o+1])
  print(aft.likli[o+1])

  print(o/iter.burnin)


  ##### Print the status bar
  Sys.sleep(0.1)
  #setTxtProgressBar(pb, o)


  ####### If At all the W get's NA
  if (sum(is.na(diag(W))+ 0) > 0){
    W <- diag(diag(cov(Y)))
  }
}


# ### Gibb's Sampling
# ############## GIBBS SAMPLING WITH THINNING ######################################################
# nmrse <- c(0)
# mu.list <- list(0)
# S.list <- list(0)
# beta0.list <- list(0)
# betahat.list <- list(0)
# sigma2.list <- list(0)
# lambda2.list <- list(0)
# tau2.list <- list(0)
# c.list <- list(0)
# That.list <- list(0)
# likli.gibbs <- c(0)
# W.list <- list(0)
# 
# 
# 
# 
# 
# count = 1
# for (o in 1:iter) {
#   ################## PARAMETERS OF THE DP Mixture Model ######################################################
#   ## Updating the parameters based on the observations
#   param <- posteriorGMMparametrs(c,Y,mu,S, alpha,K, epsilon, W, beta, ro,N,D )
#   mu <- param$mean
#   S <- param$precision
#   paramtime <- posteriortimeparameters(c, That, lambda2,tau2,sigma2,beta0, betahat, Y, K, epsilon, W, beta, ro,D, r, si, Time,N, sig2.data)
#   beta0 <- paramtime$beta0
#   betahat <- paramtime$betahat
#   sigma2 <- paramtime$sigma2
#   lambda2 <- paramtime$lambda2
#   tau2 <- paramtime$tau2
# 
#   ########################## THE HYPERPARAMETERS OF THE GMM #################################
#   #  Updating the hyper paramters
#   hypercognate <- posteriorhyperPLUS (c, Y, mu, S, epsilon, W, beta, ro )
#   epsilon <- hypercognate$epsilon
#   tmpW <- hypercognate$W
#   W <- matrix(as.matrix(tmpW),nrow = D, ncol =D)
#   ro <- hypercognate$ro
# 
#   ################# INDICATOR VARIABLE ##################################################################
#   ## Updating the indicator variables and the parameters
#   cognate <- posteriorchineseAFT(c,Y,mu,S,alpha,That, beta0, betahat, sigma2, lambda2, tau2, K, epsilon, W, beta, ro,D, r, si, Time,N, sig2.dat)
#   c <- cognate$indicator
#   mu <- cognate$mean
#   S <- cognate$precision
#   beta0 <- cognate$beta0
#   betahat <- cognate$betahat
#   sigma2 <- cognate$sigma2
#   lambda2 <- cognate$lambda2
#   tau2 <- cognate$tau2
# 
#   ########################### The Concentration Parameter #################################################################
#   # Updating the concentration parameter
#   alpha <- posterioralpha(c, N, alpha, shape.alpha, rate.alpha)
# 
# 
#   ######################## The Censored Times ###########################################################
#   # Updating the Time Variable
#   ti <- NA
#   ti <- updatetime(c, Y, Time,That, beta0, betahat, sigma2)
#   That <- ti$time
# 
# 
#   ################## PARAMETERS OF THE DP Mixture Model ######################################################
#   ## Updating the parameters based on the observations
# 
#   if(o%% iter.thin == 0 ){
#     mu.list[[count]] <- mu
#     S.list[[count]] <- S
#     W.list[[count]] <- W
#     beta0.list[[count]] <- beta0
#     betahat.list[[count]] <- betahat
#     sigma2.list[[count]] <- sigma2
#     lambda2.list[[count]] <- lambda2
#     tau2.list[[count]] <- tau2
#     c.list[[count]] <- c
#     That.list[[count]] <- That
# 
#     #     nmrse[count] <- calcrmse(time.real,time.predicted)$rmse
#     count <- count +1
#   }
# 
#   # likli.gibbs[o] <- loglikelihood(c,Y,mu,S,alpha,That, beta0, betahat, sigma2, lambda2, tau2, K, epsilon, W, beta, ro,D, r, si, Time,N, sig2.dat)
# 
# 
# 
#   #   print(o/iter)
#   #   print(loglike[o])
#   #   print(cindex)
#   #   Sys.sleep(0.1)
#   }




# 
# # ########## Analyze the fit ##########################################
# # ### Good feature selection from heatmap plus cindex plus randindex
# ########## ANLAYSING THE MCMC samples AND CALCULATING METRICES #######################################################
# source('linearprediction.R')
# count <- Nps
# final.rand <- c(0)
# 
# ############ The Matrices that will store the results #################################################
# ##### Class Assignments ########################
# c.matrix <- matrix(NA, nrow = N, ncol = count)
# recovCIndex.sbc <- c(0)
# recovCIndex.sbc.paft <- c(0)
# ################ Begin Analysig the MCMC samples #######################################################
# 
# require(Hmisc) 
# 
# for (h in 1:Nps){
#   ### See C-Index (concordance index) returned by the model
#   surv.aft <- Surv(exp(time),censoring)
#   ### Predict Time from the model
#   tem.tim <- predicttime(c.list[[h]], Y, That.list[[h]], Time, beta0.list[[h]], betahat.list[[h]], sigma2.list[[h]])$predicttime
#   library(Hmisc)
#   recovCIndex.sbc[h] <-  survConcordance(surv.aft ~ exp(-tem.tim))[[1]]
#   
#   ## try to fit a Cluster-specific AFT model on top of it ##
#   
#   ###### penAFT ###################################################################
#   ######## Penalized AFT with k-means clustering ######################################################
#   sbc.aft <- c(0)
#   for ( q in 1:F){
#     ind <- which((c.list[[h]]) == q)
#     L= length(ind)
#     
#     time.tmp <- time[ind]
#     censoring.tmp <- censoring[ind]
#     Y.tmp <- Y[ind,]
#     
#     reg <- cv.glmnet(x = Y.tmp, y = time.tmp, family = "gaussian")
#     coeff.pred <- coef(object =reg, newx = Y.tmp, s= "lambda.min")
#     sbc.aft[ind] <- predict(object = reg, newx = Y.tmp, s = "lambda.min") 
#   }
#   recovCIndex.sbc.paft[[h]] <- as.numeric(survConcordance(smod ~ exp(-sbc.aft))[1])
#   
#   c.matrix[,h] <- c.list[[h]]
#   
# }
# 
# 
# 
# 
# ###############################################
# ###### Calculating POINT ESTIMATES ############
# ###############################################
# psm <- comp.psm(t(c.matrix))
# mpear <- maxpear(psm)
# 
# ### If we build a cluster specific sbc approach
# c.final <- mpear$cl
# c.sbc <- mpear$cl
# 
# active <- as.numeric(rownames(table(c.final)))
# 
# 
# ############ Time Covariate Slopes FOR Relevant Clusters ############
# list.betahat <- list(0)
# 
# count = Nps
# 
# for ( i in 1:count){
#   list.betahat[[i]] <- (betahat.list[[i]][active,] != 0) +0
# }
# 
# Q <- length(active)
# matrix.betahat <- array(data = NA, dim =c(Q,count,D))
# 
# for ( z in 1:Q){
#   for ( x  in 1:count){
#     matrix.betahat[z,x,] <- list.betahat[[x]][z,]
#   }
# }
# 
# final.betahat <- apply(matrix.betahat,c(1,3),mean)
# ### Probability of betahat of genes FOR ONE SIMULATION
# ##colnames(final.betahat) =  c(rep("relevant",rel.D),rep("irrelevant",irrel.D))
# heatmapdata <- as.data.frame(final.betahat)
# heatmap.2(t(as.matrix(heatmapdata)),dendrogram="none", margins=c(6,10), main = "Posterior prob. \n for Selection \n ", cexCol = 0.85, cexRow = 0.7, Rowv = FALSE)
# 
# recov.logrank.sbc <- unlist(survdiff(smod ~ c.sbc))$chisq
# 
#  
# ######## Predict CLASS MEMBERSHIP on  New Data Set  BASED ON JUST THE MOLECULAR DATA #####################################
# source('priorPARAMETERS.R')
# Y.new  <- Y.input
# N.new <- nrow(Y.new)
# c.new.list <- list(0)
# ## The number of posterior samples
# 
# print("GOING THROUGH MCMC Samples")
# pb <- txtProgressBar(min = 1, max = Nps , style = 3)
# 
# 
# 
# ctemp.new <- c(0)
# modelweights <- c(0)
# 
# 
# for (count in 1:Nps){
#   
#   ## Assign the parameters to the posterior sample
#   ctemp <- c.list[[count]]
#   mu <- mu.list[[count]]
#   S <- S.list[[count]]
#   
#   g <- table(factor(ctemp, levels = 1:K))
#   activeclass <- which(g!=0)
#   ## The table function helps converting the data point specific indicator variables to class specific indicator variables
#   kminus <- length(activeclass)
#   # active <- activeclass
#   #Two Auxilary Variables
#   #The name of the auxilary variables are taken to be one and two more than the maximum value in the already active cluster set
#   activeclass <- append(activeclass, max(activeclass)+1)
#   activeclass <- append(activeclass, max(activeclass)+1)
#   active <- activeclass 
#   
#   
#   
#   ### Assigning values to parameters 
#   
#   priortwo <- NA
#   priortwo <- priordraw(beta, W, epsilon, ro, r, si,N,D, sig2.dat)
#   mu[active[kminus+1],1:D]  <- priortwo$mu  
#   S[active[kminus+1],1:D,1:D]  <- priortwo$Sigma[1:D,1:D]  
#   
#   
#   
#   
#   
#   priorthree <- NA
#   priorthree <- priordraw(beta, W, epsilon, ro, r, si,N,D, sig2.dat)
#   mu[active[kminus+2],1:D]  <- priorthree$mu  
#   S[active[kminus+2],1:D,1:D]  <- priorthree$Sigma[1:D,1:D]  
#   
#   
#   
#   ###### Some quantities used to store probabilities  
#   posteriorweight <- matrix(0, nrow = length(active), ncol = N.new)
#   weights <- matrix(0, nrow = length(active), ncol = N.new)
#   weights.final <- c(0)
#   ctemp.new <- c(0)
#   ## This can't be parallelized !!!!!
#   for(l in 1:N.new)  {
#     
#     ## Calculating the Expectations and also the normalization constant for the Expectation
#     for (j in 1:kminus) {
#       posteriorweight[j,l] <- log(g[active[j]]/ (N-1+alpha))  +  dMVN(as.vector(t(Y.new[l,1:D])), mean = mu[active[j],1:D], Q = S[active[j],1:D,1:D], log =TRUE) 
#     }
#     
#     res <- try(dMVN(as.vector(t(Y.new[l,1:D])), mean = mu[active[kminus+1],1:D], Q= S[active[kminus+1],1:D,1:D]), silent=TRUE)
#     if (class(res) == "try-error"){
#       posteriorweight[kminus+1,l] <- -Inf
#     } else{
#       posteriorweight[kminus+1,l] <-  log(alpha/ (N-1+alpha)) +  dMVN(as.vector(t(Y.new[l,1:D])), mean = mu[active[kminus+1],1:D], Q= S[active[kminus+1],1:D,1:D], log = TRUE) 
#     }
#     
#     res2 <- try(dMVN(as.vector(t(Y.new[l,1:D])), mean = mu[active[kminus+2],1:D], Q= S[active[kminus+2],1:D,1:D]), silent=TRUE)
#     if (class(res) == "try-error"){
#       posteriorweight[kminus+2,l] <- -Inf
#       
#     } else{
#       posteriorweight[kminus+2,l] <-  log(alpha/ (N-1+alpha))  +  dMVN(as.vector(t(Y.new[l,1:D])), mean = mu[active[kminus+2],1:D], Q= S[active[kminus+2],1:D,1:D], log = TRUE)      
#     }
#     
#     weights[,l] <- exp(posteriorweight[,l])/sum(exp(posteriorweight[,l]))
#     
#     
#     if (sum(exp(posteriorweight[,l])) < 1e-200){
#       ctemp.new[l] <- sample(active, 1, prob = rep(1,length(active)), replace = TRUE)
#     } else {  
#       ctemp.new[l] <- sample(active, 1, prob= weights[,l], replace = TRUE)
#       
#     }
#     
#     weights.final[l] <- dMVN(as.vector(t(Y.new[l,1:D])), mean = mu[ctemp.new[l],1:D], Q = S[ctemp.new[l],1:D,1:D], log =TRUE) 
#     
#   }
#   
#   
#   modelweights[count] <- sum(weights.final)
#   
#   c.new.list[[count]] <- ctemp.new
#   Sys.sleep(0.1)
#   setTxtProgressBar(pb, count)
#   
# }
# 
# 
# 
# 
# ## Converting the list to a matrix
# 
# c.matrix.new <- matrix(NA, nrow = N.new, ncol = Nps)
# for( h in 1:Nps){
#   c.matrix.new[,h] <- c.new.list[[h]]
# }
# 
# # ####### Use Ad-hoc methods to calculate the actual separability and C-Index ##################
# ### This function takes the posterior parameters AND predicts CLUSTER MEMEBERSHIP for the new points
# ### The posterior parameters being the RECOVERED CLUSTER MEMBERSHIP
# #### THEN IT FITS A LOGISTIC/MULTINOMIAL regression FOR LINEAR SEPARABILITY
# #### WE CAN ALSO EXPLORE NON-LINEAR Classifiers as k-SVM (RBF kernel)
# ### Check to see if this gives differences to state of the art
# 
# ### POSSIBLE APPROACHES TO GET POINT ESTIMATES #####
# 
# 
# 
# ### POSSIBILITY 1
# ### Use MPEAR APPROACH to get POINT ESTIMATES #############################################################################
# psm2 <- comp.psm(t(c.matrix.new))
# mpear2 <- maxpear(psm2)
# ### Generally the MPEAR output needs post-processing
# c.sbc.new.mpear <- mpear2$cl
# 
# 
# 
# 
# ### POSSIBILITY 2
# ### Use Logistic regression to get the labels
# reg <- cv.glmnet(x = Y, y = as.factor(c.sbc), family = "binomial")
# reg.new <- predict(object = reg, newx = Y.new, s = "lambda.min", type="class") 
# c.sbc.new.log <- as.numeric(reg.new)
# 
# 
# ### POSSIBILITY 3
# ### Use k-Nearest neighbour
# label.train <- as.factor(c.sbc)
# ### One has to to tune the k-NN classifier for k ###
# fitControl <- trainControl(method = "repeatedcv", number = 5,repeats = 5)
# ### Tune the parameter k 
# knnFit <- train(x = Y, y = label.train, method =  "knn", trControl = fitControl, tuneLength = 5)
# c.sbc.new.knn <- predict(knnFit,newdata = Y.new )
# #### Choose that configuration which has the highest difference in survival curves ####################
# 
# ### POSSIBILITY 4
# ### Use k-svm 
# 
# library(e1071)
# library(kernlab)
# df <- data.frame(cbind(c.sbc,Y))
# sig <- sigest(c.sbc ~.,data = df)
# range.sig <-  seq(from = sig[1],to = sig[3], length.out = 25)
# obj.base <- tune(svm, train.x = Y, train.y = factor(c.sbc), ranges = list(gamma = 2^(-10:14), cost = 2^(-10:14)), tunecontrol = tune.control(sampling = "cross"))
# 
# sigma.final <- obj.base$best.parameters$gamma
# cost.final <- obj.base$best.parameters$cost
# 
# # Training the Classifier
# ksvm.verk = ksvm(Y, factor(c.sbc), kernel = "rbfdot",kpar = list(sigma =sigma.final),C= cost.final, prob.model =TRUE)
# ## Predicting
# labels.ksvm <- predict(ksvm.verk ,Y.new, type = "response")
# pred.logrank.rbf <- survdiff(smod.new ~ labels.ksvm)
# c.sbc.new.rbf <- labels.ksvm
# 
# 
# ### POSSIBILITY 5
# ### Use one of the actual cluster assignments
# 
# c.matrix.new.sbc <- matrix(0, nrow = N.new, ncol = Nps + 4)
# c.matrix.new.sbc[,1:Nps] <- c.matrix.new[,1:Nps]
# c.matrix.new.sbc[,(Nps+1)] <- c.sbc.new.mpear
# c.matrix.new.sbc[,(Nps+2)] <- c.sbc.new.log
# c.matrix.new.sbc[,(Nps+3)] <- c.sbc.new.knn 
# c.matrix.new.sbc[,(Nps+4)] <- c.sbc.new.rbf
# 
# 
# Nps.mod <- Nps +4
# lr <- c(0)
# for (j in 1:Nps.mod){
#   lr[j] <-  1 - pchisq(unlist(survdiff(smod.new ~ c.matrix.new.sbc[,j]))$chisq,df = nlevels(as.factor(c.matrix.new.sbc[,i])) -1 )
# }
# 
# c.sbc.new <- c.matrix.new.sbc[,which.min(lr)]
# 
# 
# 
# 
# # ######## Predict C_INDEX on  New Data Set  BASED ON JUST THE MOLECULAR DATA #####################################
# source('priorPARAMETERS.R')
# 
# Y.new  <- Y.input
# N.new <- nrow(Y.new)
# c.new.list <- list(0)
# ## The number of posterior samples
# #That.new <- time.new 
# post.time  = matrix(NA,nrow = nrow(Y.new), ncol = Nps)
# print("GOING THROUGH MCMC Samples")
# pb <- txtProgressBar(min = 1, max = Nps , style = 3)
# 
# cind <- c(0)
# 
# # modelweights <- c(0)
# 
# predCIndex.sbc <- c(0)
# 
# 
# for (count in 1:Nps){
#   
#   ## Assign the parameters to the posterior sample
#   ctemp <- c.list[[count]]
#   mu <- mu.list[[count]]
#   S <- S.list[[count]]
#   beta0 <- beta0.list[[count]]
#   betahat  <- betahat.list[[count]] 
#   sigma2  <- sigma2.list[[count]] 
#   g <- table(factor(ctemp, levels = 1:K))
#   activeclass <- which(g!=0)
#   ## The table function helps converting the data point specific indicator variables to class specific indicator variables
#   kminus <- length(activeclass)
#   # active <- activeclass
#   #Two Auxilary Variables
#   #The name of the auxilary variables are taken to be one and two more than the maximum value in the already active cluster set
#   activeclass <- append(activeclass, max(activeclass)+1)
#   activeclass <- append(activeclass, max(activeclass)+1)
#   active <- activeclass 
#   
#   
#   
#   ### Assigning values to parameters 
#   
#   priortwo <- NA
#   priortwo <- priordraw(beta, W, epsilon, ro, r, si,N,D, sig2.dat)
#   mu[active[kminus+1],1:D]  <- priortwo$mu  
#   S[active[kminus+1],1:D,1:D]  <- priortwo$Sigma[1:D,1:D]  
#   beta0[active[kminus+1]] <- priortwo$beta0 
#   sigma2[active[kminus+1]] <- priortwo$sigma2
#   betahat[active[kminus+1],1:D] <- priortwo$betahat 
#   
#   
#   
#   
#   priorthree <- NA
#   priorthree <- priordraw(beta, W, epsilon, ro, r, si,N,D, sig2.dat)
#   mu[active[kminus+2],1:D]  <- priorthree$mu  
#   S[active[kminus+2],1:D,1:D]  <- priorthree$Sigma[1:D,1:D]  
#   beta0[active[kminus+2]] <- priorthree$beta0 
#   sigma2[active[kminus+2]] <- priorthree$sigma2
#   betahat[active[kminus+2],1:D] <- priorthree$betahat 
#   
#   
#   #######################################################
#   ########### Scaling the data ########################
#   Y.new.scaled.list <- list(0)
#   for (j in 1:kminus) {
#     clust <- which(ctemp == active[j])
#     
#     if(length(clust) > 1){
#       obj.t <- scale(Y[clust,1:D], center = TRUE, scale = TRUE)
#       Y.new.scaled.list[[j]] <- scale(Y.new, center = attr(obj.t,"scaled:center"), scale = (attr(obj.t,"scaled:scale")))
#     } else {
#       obj.t <- scale(Y[,1:D], center = TRUE, scale = TRUE)
#       Y.new.scaled.list[[j]] <- scale(Y.new, center = attr(obj.t,"scaled:center"), scale = (attr(obj.t,"scaled:scale")))
#       
#     }
#     
#     
#   }
#   
#   for (j in (kminus+1):(kminus+2)) {
#     obj.t <- scale(Y[,1:D], center = TRUE, scale = TRUE)
#     
#     Y.new.scaled.list[[j]] <- scale(Y.new, center = attr(obj.t,"scaled:center"), scale = (attr(obj.t,"scaled:scale")))
#   }
#   
#   
#   ###### Some quantities used to store probabilities  
#   posteriortime <- matrix(0, nrow = length(active), ncol = N.new)
#   posteriortimeweight <- matrix(0, nrow = length(active), ncol = N.new)
#   weights <- matrix(0, nrow = length(active), ncol = N.new)
#   
#   
#   
#   ## This can't be parallelized !!!!!
#   for(l in 1:N.new)  {
#     
#     ## Calculating the Expectations and also the normalization constant for the Expectation
#     for (j in 1:kminus) {
#       
#       posteriortime[j,l] <-  beta0[active[j]] + betahat[active[j],1:D] %*% as.vector(t(Y.new.scaled.list[[j]][l,])) 
#       
#       posteriortimeweight[j,l] <- log(g[active[j]]/ (N-1+alpha))  +  dMVN(as.vector(t(Y.new[l,1:D])), mean = mu[active[j],1:D], Q = S[active[j],1:D,1:D], log =TRUE) 
#     }
#     
#     res <- try(dMVN(as.vector(t(Y.new[l,1:D])), mean = mu[active[kminus+1],1:D], Q= S[active[kminus+1],1:D,1:D]), silent=TRUE)
#     if (class(res) == "try-error"){
#       posteriortime[kminus+1,l] <- 0
#       posteriortimeweight[kminus+1,l] <- -Inf
#     } else{
#       posteriortime[kminus+1,l] <-  beta0[active[kminus+1]] + betahat[active[kminus+1],1:D] %*% as.vector(t(Y.new.scaled.list[[j]][l,]))
#       posteriortimeweight[kminus+1,l] <-  log(alpha/ (N-1+alpha)) +  dMVN(as.vector(t(Y.new[l,1:D])), mean = mu[active[kminus+1],1:D], Q= S[active[kminus+1],1:D,1:D], log = TRUE) 
#     }
#     
#     res2 <- try(dMVN(as.vector(t(Y.new[l,1:D])), mean = mu[active[kminus+2],1:D], Q= S[active[kminus+2],1:D,1:D]), silent=TRUE)
#     if (class(res) == "try-error"){
#       posteriortime[kminus+2,l] <- 0
#       posteriortimeweight[kminus+2,l] <- -Inf
#       
#     } else{
#       posteriortime[kminus+2,l] <- beta0[active[kminus+2]] + betahat[active[kminus+2],1:D] %*% as.vector(t(Y.new.scaled.list[[j]][l,])) 
#       posteriortimeweight[kminus+2,l] <-  log(alpha/ (N-1+alpha))  +  dMVN(as.vector(t(Y.new[l,1:D])), mean = mu[active[kminus+2],1:D], Q= S[active[kminus+2],1:D,1:D], log = TRUE)      
#     }
#     
#     
#     
#     weights[,l] <- exp(posteriortimeweight[,l])/sum(exp(posteriortimeweight[,l]))
#   }
#   
#   
#   
#   for ( l in 1:N.new){
#     post.time[l,count] <- as.numeric(t(posteriortime[,l]) %*% weights[,l])
#     
#   } 
#   
#   # modelweights[count] <- sum(exp((1/N.new) *apply(posteriortimeweight,1,sum)))
#   
#   predCIndex.sbc[count] <- as.numeric(survConcordance(smod.new ~ exp(-post.time[,count]))[1]) 
#   # 
#   #    print(cind[count])
#   
#   ## Calculting the Model Weight
#   
#   Sys.sleep(0.1)
#   setTxtProgressBar(pb, count)
#   
#   
# }
# 

#### Use PC method of the reviewer #####
pX <- prcomp(Y.pre.train)
pc.recov <- predict(pX,newdata = Y.pre.train)
pc.pred <- predict(pX,newdata = Y.pre.test)
cox.fit <- coxph(smod ~ ., data=as.data.frame(pc.recov[,1:20]))
linear.PC.recovery <- predict(cox.fit, newdata=as.data.frame(pc.recov[,1:20]))
linear.PC.prediction <- predict(cox.fit,  newdata= as.data.frame(pc.pred[,1:20]))

#### Getting the C-Indices ####
recovCIndex.PC <- as.numeric(survConcordance(smod ~ linear.PC.recovery)[1])
predCIndex.PC <- as.numeric(survConcordance(smod.new ~ linear.PC.prediction)[1])


list('recov.logrank.vijver'= recov.logrank.vijver,
       'pred.logrank.vijver' = pred.logrank.vijver,
       'recovCIndex.vv.pcox' = recovCIndex.vv.pcox, 
       'predCIndex.vv.pcox'  = predCIndex.vv.pcox,
       'predCIndex.vv.kk.pcox' = predCIndex.vv.kk.pcox,
        'pred.logrank.vv.kk'  = pred.logrank.vv.kk,
        'recov.logrank.kk' = recov.logrank.kk,
        'pred.logrank.kk' = pred.logrank.kk,
        'recovCIndex.kk.pcox' = recovCIndex.kk.pcox,
       'predCIndex.kk.pcox' = predCIndex.kk.pcox,
        # 'recov.logrank.sbc' = recov.logrank.sbc,
        #  'pred.logrank.sbc' = pred.logrank.sbc,
         'recovCIndex.PC' = recovCIndex.PC,
        'predCIndex.PC' = predCIndex.PC,
       # 'recovCIndex.sbc' = recovCIndex.sbc,
       # 'recovCIndex.sbc.paft' = recovCIndex.sbc.paft,
       # 'predCIndex.sbc' = predCIndex.sbc,
      'signature.sbc' = signature.sbc)



}


### Stop the cluster
stopCluster(cluster.parallel)
stopImplicitCluster()
### End of program