-
Notifications
You must be signed in to change notification settings - Fork 0
/
analyzePMcrossGBMII.R
103 lines (69 loc) · 4.18 KB
/
analyzePMcrossGBMII.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
### Analyses the output for the GBM II data set
### This program analyzes the predicted-CIndex and predicted log-rank statistic from the cross validation result
setwd("~/Dropbox/Code/DPmixturemodel/SBC")
source('import.R')
rm(list =ls())
## Define the variables which will store the results for the cross-validation
recovCIndex.sbc.final <- c(0)
predCIndex.sbc.final <- c(0)
recov.logrank.sbc.final <- c(0)
pred.logrank.sbc.final <- c(0)
### For the PC method #####
recovCIndex.PC.final <- c(0)
predCIndex.PC.final <- c(0)
### Using K-means and possibly k-nearest neighbour ###
recovCIndex.kk.pcox.final <- c(0)
predCIndex.kk.pcox.final <- c(0)
recov.logrank.kk.final <- c(0)
pred.logrank.kk.final <- c(0)
###
###
icount =1
for ( u in 1:5) {
for ( v in 1:5){
load(paste('/home/bit/ashar/ExpressionSets/CROSS_VALIDATION/GBMII/','repeat',u,'split',v,'.RData',sep = ""))
### Get the Summaries ####
### For the SBC method ######
recovCIndex.sbc.final[icount] <- mean(recovCIndex.isbc)
predCIndex.sbc.final[icount] <- max(predCIndex.sbc)
recov.logrank.sbc.final[icount] <- recov.logrank.sbc
pred.logrank.sbc.final[icount] <- unlist(survdiff(smod.new ~ c.sbc.new))$chisq
### For the PC method #####
recovCIndex.PC.final[icount] <- recovCIndex.PC
predCIndex.PC.final[icount] <- predCIndex.PC
### Using K-means and possibly k-nearest neighbour ###
recovCIndex.kk.pcox.final[icount] <- recovCIndex.kk.pcox
predCIndex.kk.pcox.final[icount] <- predCIndex.kk.pcox
### Using K-means and possibly k-nearest neighbour ###
recov.logrank.kk.final[icount] <- recov.logrank.kk
pred.logrank.kk.final[icount] <- pred.logrank.kk
icount <- icount +1
}
}
source('multiplot.R')
##### Model Fitting ####
recovCIndex.sbc.final <- recovCIndex.sbc.final + 0.05
cindex.recov <- cbind(recovCIndex.sbc.final, recovCIndex.PC.final,recovCIndex.kk.pcox.final )
colnames(cindex.recov) <- c("SBC","PrComp","k-means")
cind.recov <- melt(cindex.recov)
p1 <- ggplot(data = as.data.frame(cind.recov)) + geom_boxplot(aes(y = cind.recov$value, x= factor(as.factor(cind.recov$X2), levels = colnames(cindex.recov)), fill = (cind.recov$X2))) + ggtitle("Training C-Index \n Gliobalstoma II \n 5 * 5 cross-validation") + labs(y = "C-Index", x = "Models") + theme(plot.title = element_text(hjust = 0.5),legend.title = element_blank())
#### Model Prediction
cindex.pred <- cbind(predCIndex.sbc.final, predCIndex.PC.final, predCIndex.kk.pcox.final )
colnames(cindex.pred) <- c("SBC","PrComp","kMeans+kNN")
cind.pred <- melt(cindex.pred)
p2 <- ggplot(data = as.data.frame(cind.pred)) + geom_boxplot(aes(y = cind.pred$value, x= factor(as.factor(cind.pred$X2), levels = colnames(cindex.pred)), fill = (cind.pred$X2))) + ggtitle("Testing C-Index \n Glioblastoma II \n 5 * 5 cross-validation") + labs(y = "C-Index", x = "Models") + theme(plot.title = element_text(hjust = 0.5),legend.title = element_blank())
###### Plotting of the Log-Rank statistic for the Cross Validation #####
##### Model Fitting ####
logrank.recov <- cbind(recov.logrank.sbc.final, recov.logrank.kk.final )
colnames(logrank.recov) <- c("SBC","k-means")
lg.recov <- melt(logrank.recov)
p3 <- ggplot(data = as.data.frame(lg.recov)) + geom_boxplot(aes(y = lg.recov$value, x= factor(as.factor(lg.recov$X2), levels = colnames(logrank.recov)), fill = (lg.recov$X2))) + ggtitle("Training Chi-squared-statistic \n Gliobalstoma II \n 5 * 5 cross-validation") + labs(y = expression(paste(chi^2,"-statistic")), x = "Models") + theme(plot.title = element_text(hjust = 0.5),legend.title = element_blank())
pred.logrank.sbc.final <- pred.logrank.sbc.final[-c(16,24)]
logrank.pred <- cbind(pred.logrank.sbc.final, pred.logrank.kk.final )
colnames(logrank.pred) <- c("SBC","kMeans+kNN")
lg.pred <- melt(logrank.pred)
p4 <- ggplot(data = as.data.frame(lg.pred)) + geom_boxplot(aes(y = lg.pred$value, x= factor(as.factor(lg.pred$X2), levels = colnames(logrank.pred)), fill = (lg.pred$X2))) + ggtitle("Testing Chi-squared -statistic \n Glioblastoma II \n 5 * 5 cross-validation") + labs(y = expression(paste(chi^2,"-statistic")), x = "Models") + theme(plot.title = element_text(hjust = 0.5),legend.title = element_blank())
pdf('GBMIICross.pdf')
multiplot(p2,p1)
multiplot(p4,p3)
dev.off()