-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathSimulation_MAIN.R
152 lines (96 loc) · 4.13 KB
/
Simulation_MAIN.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
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
### 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()
############################# 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()
########## Analyze the fit ##########################################
### Good feature selection from heatmap plus cindex plus randindex
source('MCMCanalyze.R')
MCMCanalyze()
recovRandIndex.sbc <<- as.numeric(adjustedRandIndex(c.true, c.sbc))
pc <- prcomp(Y)
pc.pred <- predict(pc,newdata = Y)
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2], colour= as.factor(c.sbc))) + ggtitle(" SBC Clustering \n Test Set") + geom_point(shape=19) + labs(y = "PC1", x = "PC2", colour = "Classes")
######## Predict on New Data Set BASED ON JUST THE MOLECULAR DATA #####################################
source('predictCLASS.R')
predictCLASS(Y.new)
## Check the predicted Rand Index
predRandIndex.sbc <<- as.numeric(adjustedRandIndex(c.true.new, c.sbc.new))
c.sbc.new.knn <<- knn(train = Y, test = Y.new, cl = c.sbc, k = F)
pc <- prcomp(Y.new)
pc.pred <- predict(pc,newdata = Y.new)
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2], colour= as.factor(c.sbc.new))) + ggtitle(" SBC Clustering \n Test Set") + geom_point(shape=19) + labs(y = "PC1", x = "PC2", colour = "Classes")
surv.fit <- survfit(smod.new ~ c.sbc.new)
logrank <- survdiff(smod.new ~ c.sbc.new)
p5 <- ggsurv(surv.fit, main = " DPMM \n Kaplan Meier Estimators") + ggplot2::guides(linetype = FALSE) + ggplot2::scale_colour_discrete(name = 'Classes',breaks = c(1,2),labels = c('1', '2'))
#### Choose that configuration which has the highest difference in survival curves
predRandIndex.sbc <- c(0)
for (j in 1:Nps){
predRandIndex.sbc[j] <- adjustedRandIndex(c.true.new,c.matrix.new[,j])
}
### Use SBC + knn to make predictions on future probabilities ####
cl.old <- c.sbc
cl.new <- c.sbc.new.knn
### Use PAFT now to build the corresponding cluster specific PAFT models
pre.sbc <- c(0)
for ( q in 1:F){
ind <- which(cl.old == q)
ind.new <- which(cl.new == q)
time.tmp <- time[ind]
Y.tmp <- Y[ind,]
Y.tmp.new <- Y.new[ind.new,]
reg <- cv.glmnet(x = Y.tmp, y = time.tmp, family = "gaussian")
pre.sbc[ind.new] <- predict(object = reg, newx = Y.tmp.new, s = "lambda.min")
}
predCIndex.sbc.knn.aft <<- as.numeric(survConcordance(smod.new ~ exp(-pre.sbc))[1])
source('predictTIME.R')
predictchineseAFTtime(Y.new)