.Rhistory

#library(xlsx)
library(XML)
library(RCurl)
#library(multiMiR)
library(MASS)
library(matrixcalc)
library(stats)
library(Runuran)
library(MCMCpack)
library(VGAM)
library(statmod)
library(survcomp)
library(gtools)
library(ggplot2)
library(GenOrd)
library(plyr)
library(iCluster)
library(CCA)
#library(PReMiuM)
library(caret)
library(mcclust)
library(Biobase)
library(mixtools)
rm(list = ls())
#################################### SIMULATED DATA PROPERTIES ####################################################
## Number of points
N.test = 500
N.train = 500
## Number of Clusters
F = 2
k =F
N = N.train
## Distribution of the points within three clusters
p.dist = c(0.5,0.5)
## Total Number of features D
D = 20
## Total Percentage of irrelevant feature
prob.noise.feature = 0.50
## Overlap between Cluster of molecular Data of the relevant features
prob.overlap = 0.05
###### Get the Data #####################################
## Initialize the Training Data
source('simulate.R')
simulate()
####### Assign training and testing data ###############
Y <- Y.dat
Y.new <- Y.new.dat
smod <-  Surv(exp(time), censoring)
smod.new <- Surv(exp(time.new), censoring.new)
source('Comparisonx.R')
Comparisonx()
smod <-  Surv(exp(time), censoring)
smod.new <-  Surv(exp(time.new), censoring.new)
### Fitting A Penalized Cox Proportional Hazard's Model
reg.pcox <- cv.glmnet(x = Y, y = smod, family = "cox")
lp <- predict(object =reg.pcox, newx = Y, s= "lambda.min")
recovCIndex.na.cox <<- as.numeric(survConcordance(smod ~lp)[1])
recovCIndex.na.cox
linear.pred.cox <- c(0)
### see if we can use glmnet
reg.pcox <- cv.glmnet(x = Y, y = Surv(exp(time), censoring), family = "cox")
linear.pred.cox <- predict(object =reg.pcox, newx = Y.new, s= "lambda.min")
predCIndex.na.cox <<- as.numeric(survConcordance(smod.new ~ linear.pred.cox)[1])
predCIndex.na.cox
reg.pcox
plot(reg.pcox)
plot(reg.pcox$cvm)
reg.pcox$nzero
dim(Y)
smod <-  Surv(exp(time), censoring)
smod.new <-  Surv(exp(time.new), censoring.new)
### Fitting A Penalized Cox Proportional Hazard's Model
reg.pcox <- cv.glmnet(x = Y, y = smod, family = "cox")
lp <- predict(object =reg.pcox, newx = Y, s= "lambda.min")
recovCIndex.na.cox <<- as.numeric(survConcordance(smod ~lp)[1])
######## Prediction with penalized Cox Proportional Hazards Model ###########################################
linear.pred.cox <- c(0)
### see if we can use glmnet
reg.pcox <- cv.glmnet(x = Y, y = Surv(exp(time), censoring), family = "cox")
linear.pred.cox <- predict(object =reg.pcox, newx = Y.new, s= "lambda.min")
predCIndex.na.cox <<- as.numeric(survConcordance(smod.new ~ linear.pred.cox)[1])
reg <- cv.glmnet(x = Y, y = time, family = "gaussian")
linear.aft <- predict(object = reg, newx = Y, s = "lambda.min")
cindex.paft <- as.numeric(survConcordance(smod ~ exp(-linear.aft))[1])
recovCIndex.na.aft <<- as.numeric(cindex.paft)
recovCIndex.na.aft
linear.pred.paft <- c(0)
### see if we can use glmnet
reg.paft <- cv.glmnet(x = Y, y = time, family = "gaussian")
linear.pred.paft <- predict(object = reg.paft, newx = Y.new, s= "lambda.min")
predCIndex.na.aft <<- as.numeric(survConcordance(smod.new ~ exp(-linear.pred.paft))[1])
predCIndex.na.aft
gr.km <- kmeans(Y, F, nstart =10)
recovRandIndex.km <<- adjustedRandIndex(c.true,as.factor(gr.km$cluster))
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)
### Tune the parameter k
knnFit <- train(x = Y, y = label.train, method =  "knn", trControl = fitControl, preProcess = c("center","scale"), tuneLength = 5)
knnPredict <- predict(knnFit,newdata = Y.new )
label.test <- knnPredict
predRandIndex.knear <<- adjustedRandIndex(c.true.new, label.test)
colnames(Y)
dim(Y)
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)
### Tune the parameter k
knnFit <- train(x = Y, y = label.train, method =  "knn", trControl = fitControl, preProcess = c("center","scale"), tuneLength = 5)
true.knn <- predict(knnFit,newdata = Y.new )
predRandIndex.true.knear <<- adjustedRandIndex(c.true.new, true.knn)
gr.km <- kmeans(Y, F, nstart =10)
adjustedRandIndex(c.true,as.factor(gr.km$cluster))
gr.km <- kmeans(Y, F, nstart =10)
recovRandIndex.km <<- adjustedRandIndex(c.true,as.factor(gr.km$cluster))
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)
gr.km <- kmeans(Y, F, nstart =10)
recovRandIndex.km <<- adjustedRandIndex(c.true,as.factor(gr.km$cluster))
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)
### Tune the parameter k
knnFit <- caret::train(x = Y, y = label.train, method =  "knn", trControl = fitControl, preProcess = c("center","scale"), tuneLength = 5)
caret::train()
caret::predict.train()
### Tune the parameter k
knnFit <- caret::predict.train(x = Y, y = label.train, method =  "knn", trControl = fitControl, preProcess = c("center","scale"), tuneLength = 5)
caret::train()
### Tune the parameter k
knnFit <- caret::train(x = Y, y = label.train, method =  "knn", trControl = fitControl, preProcess = c("center","scale"), tuneLength = 5)
### Tune the parameter k
knnFit <- caret::train(x = as.data.frame(Y), y = label.train, method =  "knn", trControl = fitControl, preProcess = c("center","scale"), tuneLength = 5)
knnPredict <- predict(knnFit,newdata = as.data.frame(Y.new ))
knnPredict
label.test <- knnPredict
adjustedRandIndex(c.true.new, label.test)
recovRandIndex.km
linear.cox <- c(0)
for ( q in 1:F){
ind <- which((gr.km$cluster) == q)
time.tmp <- time[ind]
censoring.tmp <- censoring[ind]
Y.tmp <- Y[ind,]
coxreg <- list(0)
coxreg$x <- Y.tmp
coxreg$time <- exp(time.tmp)
coxreg$status <- censoring.tmp
reg.pcox <- cv.glmnet(x = Y.tmp, y = Surv(coxreg$time, coxreg$status), family = "cox")
linear.cox[ind] <- predict(object =reg.pcox, newx = Y.tmp, s= "lambda.min")
}
recovCIndex.km.pcox <<- as.numeric(survConcordance(smod ~ linear.cox)[1])
### Prediction with k-means + k-nearest neghbour
linear.kkpcox.prediction <- c(0)
for ( q in 1:F){
ind <- which(label.train == q)
ind.new <- which(label.test == q)
reg.aft <- cv.glmnet(x = Y[ind,], y = Surv(exp(time[ind]),censoring[ind]), family = "cox")
linear.kkpcox.prediction[ind.new] <- predict(object =reg.aft, newx = Y.new[ind.new,], s= "lambda.min")
}
predCIndex.kk.pcox <<- as.numeric(survConcordance(smod.new ~ linear.kkpcox.prediction)[1])
linear.aft <- c(0)
for ( q in 1:F){
ind <- which((gr.km$cluster) == 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")
rel.coeff <- coeff.pred[2:(D+1)]
ind.rel <- which(rel.coeff !=0)
linear.aft[ind] <- predict(object = reg, newx = Y.tmp, s = "lambda.min")
}
recovCIndex.km.paft <<- as.numeric(survConcordance(smod ~ exp(-linear.aft))[1])
#### prediction with penAFT ###################################################################
linear.kkpaft.prediction <- c(0)
for ( q in 1:F){
ind <- which(label.train == q)
ind.new <- which(label.test == q)
reg.aft <- cv.glmnet(x = Y[ind,], y = time[ind], family = "gaussian")
linear.kkpaft.prediction[ind.new] <- predict(object =reg.aft, newx = Y.new[ind.new,], s= "lambda.min")
}
predCIndex.kn.paft <<- as.numeric(survConcordance(smod.new ~ exp(-linear.kkpaft.prediction))[1])
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)
### Tune the parameter k
knnFit <- train(x = Y, y = label.train, method =  "knn", trControl = fitControl, preProcess = c("center","scale"), tuneLength = 5)
true.knn <- predict(knnFit,newdata = Y.new )
predRandIndex.true.knear <<- adjustedRandIndex(c.true.new, true.knn)
linear.pred <- c(0)
cox.pred <- c(0)
for ( q in 1:F){
ind <- which((c.true) == q)
ind.new <- which(true.knn == q)
time.tmp <- time[ind]
censoring.tmp <- censoring[ind]
Y.tmp <- Y[ind,1:D]
Y.tmp.new <- Y.new[ind.new,1:D]
coxreg <- list(0)
coxreg$x <- Y.tmp
coxreg$time <- exp(time.tmp)
coxreg$status <- censoring.tmp
reg.pcox <- cv.glmnet(x = Y.tmp, y = Surv(coxreg$time, coxreg$status), family = "cox")
linear.pred[ind] <- predict(object =reg.pcox, newx = Y.tmp, s= "lambda.min")
cox.pred[ind.new] <- predict(object =reg.pcox, newx = Y.tmp.new, s= "lambda.min")
}
recovCIndex.true.pcox <<- as.numeric(survConcordance(smod ~ linear.pred)[1])
predCIndex.true.knn.pcox <<- as.numeric(survConcordance(smod.new ~ cox.pred)[1])
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)
### Tune the parameter k
knnFit <- caret::train(x = as.data.frame(Y), y = label.train, method =  "knn", trControl = fitControl, preProcess = c("center","scale"), tuneLength = 5)
true.knn <- caret::predict(knnFit,newdata = as.data.frame(Y.new ))
predRandIndex.true.knear <<- adjustedRandIndex(c.true.new, true.knn)
true.knn <- predict(knnFit,newdata = as.data.frame(Y.new ))
adjustedRandIndex(c.true.new, true.knn)
predRandIndex.true.knear <<- adjustedRandIndex(c.true.new, true.knn)
linear.pred <- c(0)
cox.pred <- c(0)
for ( q in 1:F){
ind <- which((c.true) == q)
ind.new <- which(true.knn == q)
time.tmp <- time[ind]
censoring.tmp <- censoring[ind]
Y.tmp <- Y[ind,1:D]
Y.tmp.new <- Y.new[ind.new,1:D]
coxreg <- list(0)
coxreg$x <- Y.tmp
coxreg$time <- exp(time.tmp)
coxreg$status <- censoring.tmp
reg.pcox <- cv.glmnet(x = Y.tmp, y = Surv(coxreg$time, coxreg$status), family = "cox")
linear.pred[ind] <- predict(object =reg.pcox, newx = Y.tmp, s= "lambda.min")
cox.pred[ind.new] <- predict(object =reg.pcox, newx = Y.tmp.new, s= "lambda.min")
}
recovCIndex.true.pcox <<- as.numeric(survConcordance(smod ~ linear.pred)[1])
predCIndex.true.knn.pcox <<- as.numeric(survConcordance(smod.new ~ cox.pred)[1])
km.perm <- KMeansSparseCluster.permute(x = Y, K= k ,wbounds=c(1.5,2:6),nperms=5)
km.out <- KMeansSparseCluster(x = Y, K=k,wbounds=km.perm$bestw)
recovRandIndexSKM <<- adjustedRandIndex(km.out[[1]]$Cs, c.true)
###################################################
########### sparse hierarchical clustering #########################
#############################################
#############################################
perm.out <- HierarchicalSparseCluster.permute(x = Y, wbounds=c(1.5,2:6), nperms=5)
# Perform sparse hierarchical clustering
sparsehc <- HierarchicalSparseCluster(dists=perm.out$dists, wbound=perm.out$bestw, method="complete")
recovRandIndexSHC <<- adjustedRandIndex(cutree(sparsehc$hc, k = k), c.true)
rm(list = ls())
#################################### SIMULATED DATA PROPERTIES ####################################################
## Number of points
N.test = 500
N.train = 500
## Number of Clusters
F = 2
k =F
N = N.train
## Distribution of the points within three clusters
p.dist = c(0.5,0.5)
## Total Number of features D
D = 20
## Total Percentage of irrelevant feature
prob.noise.feature = 0.50
## Overlap between Cluster of molecular Data of the relevant features
prob.overlap = 0.05
###### Get the Data #####################################
## Initialize the Training Data
source('simulate.R')
simulate()
####### Assign training and testing data ###############
Y <- Y.dat
Y.new <- Y.new.dat
smod <-  Surv(exp(time), censoring)
smod.new <- Surv(exp(time.new), censoring.new)
##################### STATE OF THE ART TECHNIQUES #################################
##################### BASIC METHODS + SOME ADVANCED METHODS ########################
source('Comparisonx.R')
Comparisonx()
######################### Initialize the Parameters ##############################
source('initialize.R')
initialize()
###################### Start with a good configuration ###########################
source('startSBC.R')
startSBC()
install.packages('flexmix')
library(flexmix)
source('ComparisionFLX.R')
ComparisionFLX()
source('startSBC.R')
startSBC()
### This is the Main Function and contains a simulation case
### Also CHECK THE TIME REQUIRED FOR THE MODEL
setwd("~/Dropbox/Code/DPmixturemodel/SBC")
rm(list = ls())
#################################### SIMULATED DATA PROPERTIES ####################################################
## Number of points
N.test = 500
N.train = 500
## Number of Clusters
F = 2
k =F
N = N.train
## Distribution of the points within three clusters
p.dist = c(0.5,0.5)
## Total Number of features D
D = 20
## Total Percentage of irrelevant feature
prob.noise.feature = 0.50
## Overlap between Cluster of molecular Data of the relevant features
prob.overlap = 0.05
###### Get the Data #####################################
## Initialize the Training Data
source('simulate.R')
simulate()
####### Assign training and testing data ###############
Y <- Y.dat
Y.new <- Y.new.dat
smod <-  Surv(exp(time), censoring)
smod.new <- Surv(exp(time.new), censoring.new)
##################### STATE OF THE ART TECHNIQUES #################################
##################### BASIC METHODS + SOME ADVANCED METHODS ########################
source('Comparisonx.R')
Comparisonx()
source('ComparisionFLX.R')
ComparisionFLX()
#
# source('ComparisionPReMiuM.R')
# ComparisionPReMiuM()
# setwd("~/Dropbox/Code/DPmixturemodel/SBC")
######################### Initialize the Parameters ##############################
source('initialize.R')
initialize()
###################### Start with a good configuration ###########################
source('startSBC.R')
startSBC()
table(c)
c
adjustedRandIndex(c,c.true)
### This is the Main Function and contains a simulation case
### Also CHECK THE TIME REQUIRED FOR THE MODEL
setwd("~/Dropbox/Code/DPmixturemodel/SBC")
rm(list = ls())
#################################### SIMULATED DATA PROPERTIES ####################################################
## Number of points
N.test = 500
N.train = 500
## Number of Clusters
F = 2
k =F
N = N.train
## Distribution of the points within three clusters
p.dist = c(0.5,0.5)
## Total Number of features D
D = 20
## Total Percentage of irrelevant feature
prob.noise.feature = 0.50
## Overlap between Cluster of molecular Data of the relevant features
prob.overlap = 0.05
###### Get the Data #####################################
## Initialize the Training Data
source('simulate.R')
simulate()
####### Assign training and testing data ###############
Y <- Y.dat
Y.new <- Y.new.dat
smod <-  Surv(exp(time), censoring)
smod.new <- Surv(exp(time.new), censoring.new)
##################### STATE OF THE ART TECHNIQUES #################################
##################### BASIC METHODS + SOME ADVANCED METHODS ########################
source('Comparisonx.R')
Comparisonx()
source('ComparisionFLX.R')
ComparisionFLX()
#
# source('ComparisionPReMiuM.R')
# ComparisionPReMiuM()
# setwd("~/Dropbox/Code/DPmixturemodel/SBC")
######################### Initialize the Parameters ##############################
source('initialize.R')
initialize()
c
adjustedRandIndex(c,c.true)
recovCIndex.km.paft
recovCIndex.flx.paft
data <- data.frame(time, Y)
fo <- sample(rep(seq(10), length = nrow(data)))
gr.flx <- flexmix(time ~ ., data = data, k = F, cluster = gr.km$cluster, model = FLXMRglmnet(foldid = fo, adaptive= FALSE, family = c("gaussian")), control = list(iter.max = 500))
gr.km <- kmeans(Y, F, nstart =10)
## FlexMix Clustering
data <- data.frame(time, Y)
fo <- sample(rep(seq(10), length = nrow(data)))
gr.flx <- flexmix(time ~ ., data = data, k = F, cluster = gr.km$cluster, model = FLXMRglmnet(foldid = fo, adaptive= FALSE, family = c("gaussian")), control = list(iter.max = 500))
linear.flx <- c(0)
beta.flx <- matrix(0, nrow = D, ncol = F)
for ( q in 1:F){
ind <- which(clusters(gr.flx) == 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")
rel.coeff <- coeff.pred[2:(D+1)]
beta.flx[1:D,q] <- rel.coeff
linear.flx[ind] <- predict(object = reg, newx = Y.tmp, s = "lambda.min")
}
recovCIndex.flx.paft  <- as.numeric(survConcordance(smod ~ exp(-linear.flx))[1])
recovCIndex.flx.paft
recovCIndex.km.paft
### This is the Main Function and contains a simulation case
### Also CHECK THE TIME REQUIRED FOR THE MODEL
setwd("~/Dropbox/Code/DPmixturemodel/SBC")
rm(list = ls())
#################################### SIMULATED DATA PROPERTIES ####################################################
## Number of points
N.test = 500
N.train = 500
## Number of Clusters
F = 2
k =F
N = N.train
## Distribution of the points within three clusters
p.dist = c(0.5,0.5)
## Total Number of features D
D = 20
## Total Percentage of irrelevant feature
prob.noise.feature = 0.50
## Overlap between Cluster of molecular Data of the relevant features
prob.overlap = 0.05
###### Get the Data #####################################
## Initialize the Training Data
source('simulate.R')
simulate()
####### Assign training and testing data ###############
Y <- Y.dat
Y.new <- Y.new.dat
smod <-  Surv(exp(time), censoring)
smod.new <- Surv(exp(time.new), censoring.new)
##################### STATE OF THE ART TECHNIQUES #################################
##################### BASIC METHODS + SOME ADVANCED METHODS ########################
source('Comparisonx.R')
Comparisonx()
source('ComparisionFLX.R')
ComparisionFLX()
#
# source('ComparisionPReMiuM.R')
# ComparisionPReMiuM()
# setwd("~/Dropbox/Code/DPmixturemodel/SBC")
######################### Initialize the Parameters ##############################
source('initialize.R')
initialize()
###################### Start with a good configuration ###########################
source('startSBC.R')
startSBC()
### This is the Main Function and contains a simulation case
### Also CHECK THE TIME REQUIRED FOR THE MODEL
setwd("~/Dropbox/Code/DPmixturemodel/SBC")
rm(list = ls())
#################################### SIMULATED DATA PROPERTIES ####################################################
## Number of points
N.test = 500
N.train = 500
## Number of Clusters
F = 2
k =F
N = N.train
## Distribution of the points within three clusters
p.dist = c(0.5,0.5)
## Total Number of features D
D = 40
## Total Percentage of irrelevant feature
prob.noise.feature = 0.50
## Overlap between Cluster of molecular Data of the relevant features
prob.overlap = 0.15
###### Get the Data #####################################
## Initialize the Training Data
source('simulate.R')
simulate()
####### Assign training and testing data ###############
Y <- Y.dat
Y.new <- Y.new.dat
smod <-  Surv(exp(time), censoring)
smod.new <- Surv(exp(time.new), censoring.new)
##################### STATE OF THE ART TECHNIQUES #################################
##################### BASIC METHODS + SOME ADVANCED METHODS ########################
source('Comparisonx.R')
Comparisonx()
source('ComparisionFLX.R')
ComparisionFLX()
#
# source('ComparisionPReMiuM.R')
# ComparisionPReMiuM()
# setwd("~/Dropbox/Code/DPmixturemodel/SBC")
######################### Initialize the Parameters ##############################
source('initialize.R')
initialize()
###################### Start with a good configuration ###########################
source('startSBC.R')
startSBC()
############################# PARAMETERS for GIBB's SAMPLING ####
iter = 20
iter.burnin = 50
iter.thin  = 2
Nps = as.integer(iter/iter.thin)
########### Train the Model #########################################
source('burninDPMM.R')
burninDPMM()
source('gibbsDPMM.R')
gibbsDPMM()