-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathMAIN-Revision.R
348 lines (289 loc) · 18.2 KB
/
MAIN-Revision.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
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
################################################################################################################
# DREAM CHALLENGES SUBMISSION (CHALLENGE 2.B)
#
# TEAM YODA | UNIVERSITY OF WASHINGTON TACOMA
# SUMMER - 2015
#
# Predicting Discontinuation of Docetaxel Treatment for Metastatic
# Castration-Resistant Prostate Cancer (mCRPC) with hill-climbing and Random Forest
#
# Authors Statement
# Kevin Anderson, Seyed Khankhajeh, Daniel Kristiyanto, Kaiyuan Shi, Seth West contributed equally as first authors.
#
# Kevin Anderson, Seyed Khankhajeh, Daniel Kristiyanto, Kaiyuan Shi, Seth West, Azu Lee, Qi Wei, Migao Wu, Yunhong Yin
# conducted empirical studies comparing the performance of various machine learning and feature selection algorithms
# on this dataset.
#
# Daniel Kristiyanto served as the team captain and managed the submission uploads.
# Ka Yee Yeung served as the PI, used this DREAM 9.5 challenge as a class project, provided guidance and co-ordinated
# this submission effort.
#
#
###########################################################################################
# ADJUST-ACCORDINGLY (IFNECESSARY)
path.script = "~/Dropbox/_SPACE/DREAM9.5/"
setwd(path.script)
library(neuralnet)
require(FSelector)
require(e1071)
require(timeROC)
require(caret)
require(unbalanced)
require(impute)
require(randomForest)
library(FSelector)
library(e1071)
library(impute)
library(caret)
library(randomForest)
source(paste0(path.script,"inc_datacleanup.r")) # Contains functions to clean up the data
source(paste0(path.script,"score.R")) # Contains Score.R from Synapse as scoring
source(paste0(path.script,"inc_functions.R")) # Contains various functions (evaluator, etc)
###### GLOBAL VARIALBLES ##############################################################
## GLOBAL VARIABLES AND PARAMETERS
CV = c(rep(1:10)) # Number of Crossfold Validation, each fold takes about 40 mins
#CV = c(rep(1:9), seq(10,100,10))
# Itterate through every study with following target
# (and use the other two as training data)
# 1. ACCENT2, 2. EFC6546, 3. CELGENE 4. ALL COMBINED
STUDY = c(1,2,3,4)
# STUDY = c(4) # Enable to perform prediction for the Core_Validation Dataset
## MODEL TUNING
FOLD.RATIO = 0.9 # How many goes as training data, for STUDY= 4 ONLY
baggingRatio = 0.9 # Used by Feature Selection Evaluation.
balace.ratio = 31
rf.ntree = 180
rf.mtry = 4
k = 6
# k = 43
x.axis = "Cross-fold Validation"
####### VARIABLE LIST ##############################################################
# IN GENERAL, HERE'S THE VARIABLES
core_table_training <- NULL # ORIGINAL CORE TRAINING (RAW)
core_table_validation <- NULL # ORIGINAL CORE TESTING (RAW)
table.for.model <- NULL # CLEANED, AUGMENTED CORE TRAINING
table.for.validation <- NULL # CLEANED, AUGMENTED CORE TESTING
curr.training.data <- NULL # A SUBSET OF CURRENT FOLD'S TRAINING DATA
curr.testing.data <- NULL # A SUBSET OF CURRENT FOLD'S TESTING DATA
SCORING.TABLE <- NULL # SCORE RESULTS
OUTPUT.TABLE <- NULL # SUBSET OF CURR.TESTING DATA WITH PROB RESULT
FINAL.TABLE <- NULL # CORE TESTING/VALIDATION FROM SYNAPSE WITH PREDICTION
#
####### PREPROCESS ##############################################################
## GET THE DATA
core_table_training <- read.csv("./DATA/CoreTable_training.csv" , stringsAsFactors = F)
core_table_validation1 <- read.csv("./DATA/CoreTable_validation.csv", stringsAsFactors = F)
core_table_validation2 <- read.csv("./DATA/CoreTable_leaderboard.csv", stringsAsFactors = F)
core_table_validation <- rbind(core_table_validation1, core_table_validation2)
## GET THE DATA CLEANED UP
#core_table_validation$DISCONT <- 0
combined_data <- rbind(core_table_training, core_table_validation)
clean_data <- dream9.cleanData(combined_data)
## SOME DATA CLEAN UP
clean_data$DEATH[clean_data$DEATH == "."] <- NA
clean_data$DEATH[clean_data$DEATH == "YES"] <- 1
clean_data$DEATH[clean_data$DEATH == ""] <- 0
clean_data$LKADT_P[clean_data$LKADT_P == "."] <- NA
clean_data$ENDTRS_C[clean_data$ENDTRS_C == "."] <- "UNKNOWN"
clean_data$ENTRT_PC[clean_data$ENTRT_PC == "."] <- NA
clean_data$PER_REF[clean_data$PER_REF == "."] <- NA
clean_data$LKADT_REF[clean_data$LKADT_REF == "."] <- NA
clean_data$ENTRT_PC[clean_data$ENTRT_PC == "."] <- NA
####### EXTRAPOLATE ##############################################################
dependant.cols <- c("DEATH", "LKADT_P", "ENDTRS_C", "ENTRT_PC", "PER_REF", "LKADT_REF","TSTAG_DX")
clean_data <- clean_data[,!(colnames(clean_data) %in% dependant.cols)]
## VARIABLE EXTRAPOLATE / AUGMENT
# to.explode <- c("AGEGRP2","RACE_C","REGION_C","SMOKE","SMOKFREQ","SMOKSTAT","ECOG_C","TRT2_ID","TRT3_ID")
to.explode <- c("TRT2_ID","SMOKE","SMOKFREQ","SMOKSTAT")
tmp <- clean_data[,to.explode]
f <- paste("~0+",paste(to.explode, collapse = "+"))
tmp <- model.matrix(as.formula(f), data=tmp)
tmp <- apply(tmp,2,factor)
clean_data <- cbind(clean_data,tmp)
factor.cols <- union(colnames(tmp),factor.cols)
binary.cols <- union(colnames(tmp),binary.cols)
to.drop <- c(to.explode,"TRT3_ID","GLEAS_DX","NON_TARGET", "ABDOMINAL","BONE","WEIGHTBL","HEIGHTBL")
clean_data <- clean_data[,!(colnames(clean_data) %in% to.drop)]
row.names(clean_data) <- clean_data$RPT
halabi <- c("RACE_C","AGEGRP2","BMI","PRIOR_RADIOTHERAPY","ANALGESICS","ECOG_C","GLEAS_DX","ALB","LDH", "WBC",
"AST","TBILI","PLT", "HB","ALT","TESTO","PSA","ALP")
halabi <- c("RACE_C","AGEGRP2","BMI","PRIOR_RADIOTHERAPY","ANALGESICS","ECOG_C","ALB","LDH", "WBC",
"AST","TBILI","PLT", "HB","ALT","PSA","ALP")
halabi <- c(halabi,
"RECTAL",
"LYMPH_NODES",
"KIDNEYS",
"LUNGS",
"LIVER",
"PLEURA",
"OTHER",
"PROSTATE",
"ADRENAL",
"BLADDER",
"PERITONEUM",
"COLON",
"SOFT_TISSUE"
)
# "HEAD_AND_NECK",
# "STOMACH",
# "PANCREAS",
# "THYROID"
## SPLIT THE DATA BACK TO TRAINING AND VALIDATION
removed.from.validation <- c("DISCONT","WEIGHTBL","HEIGHTBL")
table.for.validation <- clean_data[clean_data$STUDYID == "AZ",!(colnames(clean_data) %in% removed.from.validation)]
table.for.model <- clean_data[((clean_data$STUDYID != "AZ") & !is.na(clean_data$DISCONT)), ]
## DROP DEPENDANT VARIABLES AWAY
to.drop <- c("LKADT_P","DEATH","PER_REF","LKADT_REF","LKADT_PER", "TSTAG_DX")
table.for.model <- table.for.model[,!(colnames(table.for.model) %in% to.drop)]
# table.for.validation <- clean_validation
split.table <- split(table.for.model, table.for.model$STUDYID)
####### CROSSFOLD START ##############################################################
# CV <- rep(1:length(intersect(halabi, colnames(table.for.model))))
for(cv in CV) #Begin Cross-Fold for Validation or for Model Tuning
{
# Unless tuning is performed, these vaiables should be disabled
# rf.ntree = cv
# balace.ratio = cv
# rf.mtry = cv
k = cv
library(caret)
train.index <- createDataPartition(table.for.model$DISCONT, p = FOLD.RATIO,list = FALSE, times = 1)
detach(package:caret, unload=T, force=T)
for (curr.study in STUDY) ## LOOP THROUGH THE DATA FRAMES
{
## RETRIEVE THE CLEAN DATA SET
## THE FUNCTION ASSUMES THAT THE CORE DATA IS ALREADY LOADED
curr.training.data <- as.data.frame(get.training.data(fold=curr.study))
curr.testing.data <- as.data.frame(get.testing.data(fold=curr.study))
testing.name <- get.testing.name(curr.study)
testing.name <- get.testing.name(curr.study)
if (curr.study == 4 ){ # USE ALL FOR TESTING/TRAINING
set.seed(123) # Tuning? Normal: Disabled #####
curr.training.data <- table.for.model[train.index,]
curr.testing.data <- table.for.model[-train.index,]
testing.name <- paste("ALL-fold")
}
# ENABLE THIS TO SCORE AZ
#curr.training.data <- as.data.frame(rbind(curr.training.data,curr.testing.data))
## DROP VARIABLES THAT ARE NOT AVAILABLE IN TESTING FROM TRAINING AND SOME REDUDANT
# to.drop <- union(to.drop, c("SMOKFREQ", "SMOKSTAT", "HEIGHTBL", "WEIGHTBL", "WEIGHT", "NON_TARGET", "AGE", "X"))
# tmp <- curr.testing.data[curr.testing.data$DISCONT==1,]
# tmp <- tmp[,-c(1,2,3)]
# tmp <- sapply(tmp, var)
# tmp <- tmp[tmp==0]
# to.drop <- union(to.drop, attr(tmp,"names"))
# curr.training.data <- curr.training.data[,!names(curr.training.data) %in% to.drop]
#
## UPDATE THE LIST OF VARIABLE GROUPS
all.features <- intersect(colnames(curr.training.data),colnames(curr.testing.data))
binary.cols <- intersect(all.features, binary.cols)
numeric.cols <- intersect(all.features, numeric.cols)
cols.to.convert <- intersect(all.features, factor.cols)
####### BALANCE THE TRAINING DATA ##############################################################
library(mlr)
library(unbalanced)
# set.seed(10)
balancer.ubunder <- ubUnder(X=curr.training.data[,-c(3)], Y=curr.training.data[,'DISCONT'], perc = balace.ratio, method = "percUnder")
curr.training.data <- cbind(balancer.ubunder$Y, balancer.ubunder$X)
names(curr.training.data)[1] <- paste("DISCONT")
detach(package:unbalanced, unload=T, force=T)
detach(package:mlr, unload=T, force=T)
####### METADATA EDITOR ##############################################################
for(i in 1:length(cols.to.convert))
{
cname <- cols.to.convert[i]
# Convert to R Valid Names. Caret is very picky
curr.testing.data[,cname] <- as.factor(make.names(curr.testing.data[,cname]))
curr.training.data[,cname] <- as.factor(make.names(curr.training.data[,cname]))
table.for.validation[,cname] <- as.factor(make.names(table.for.validation[,cname]))
# Make sure that factor levels between training and testing match
level.list <- union(levels(curr.training.data[,cname]), levels(curr.testing.data[,cname]))
level.list <- union(levels(table.for.validation[,cname]),level.list)
levels(curr.training.data[,cname]) <- level.list
levels(curr.testing.data[,cname]) <- level.list
levels(table.for.validation[,cname]) <- level.list
}
##### FEATURE SELECTION ##############################################################
##### RF
# use.all.features <- colnames(curr.training.data[,setdiff(names(table.for.validation),c("RPT","STUDYID"))])
weights <- random.forest.importance(DISCONT ~., curr.training.data[,setdiff(names(curr.training.data),c(halabi,"RPT","STUDYID"))], importance.type = 1)
features <- c(halabi, cutoff.k(weights,k))
# ##### HC
# bag0 <- subset(curr.training.data, DISCONT=="0")
# bag1 <- subset(curr.training.data, DISCONT=="1")
# features <- setdiff(all.features,numeric.cols)
# features <- features[-c(1,2,3)]
# hill.climbing <- hill.climbing.search(features, evaluatorrf)
#
#### SAVED FEATURES FROM PREVIOUS RUN
HC.4 <- c("LYMPH_NODES","LIVER", "PLEURA", "PROSTATE", "GONADOTROPIN", "HMG_COA_REDUCT", "CEREBACC", "MHENDO", "AGEGRP2X..75", "ECOG_C2", "TRT3_IDUNKNOWN", "ALP", "NEU", "PLT", "TBILI", "CREACL", "MG", "GLU")
HC.Sanofi <- c("LYMPH_NODES", "LUNGS", "PLEURA", "PROSTATE", "BLADDER", "ORCHIDECTOMY", "BILATERAL_ORCHIDECTOMY", "GLUCOCORTICOID", "GONADOTROPIN", "CORTICOSTEROID", "IMIDAZOLE", "CEREBACC", "MHHEPATO", "MHINJURY", "MHRESP", "ALP", "AST", "CREAT", "HB", "NEU", "PSA", "CREACL", "MG", "ALB")
HC.Celgene <- c("NA.", "LYMPH_NODES", "OTHER", "PROSTATE", "TURP", "BILATERAL_ORCHIDECTOMY", "CORTICOSTEROID", "ACE_INHIBITORS", "BETA_BLOCKING", "GASTREFL", "PULMEMB", "COPD", "AGEGRP2X..75", "RACE_CWhite", "REGION_CWESTERN.EUROPE", "TRT2_IDUNKNOWN", "CREAT", "PSA", "TBILI", "CREACL", "ALB", "TPRO", "BUN")
HC.Ascent <- c("NA.", "LIVER", "OTHER", "ORCHIDECTOMY", "ANALGESICS", "ANTI_ANDROGENS", "HMG_COA_REDUCT", "ESTROGENS", "CEREBACC", "DVT", "PUD", "SPINCOMP", "AGEGRP2X65.74", "SMOKEYES", "SMOKSTATUNKNOWN", "TRT2_IDUNKNOWN","TRT3_IDUNKNOWN", "AST", "HB", "WBC", "MG", "ALB")
Halabi <- c("RACE_C", "AGEGRP2", "BMI", "PRIOR_RADIOTHERAPY", "ANALGESICS", "ECOG_C", "ALB", "LDH", "WBC", "AST", "TBILI", "PLT", "HB", "ALT", "PSA", "ALP", "RECTAL", "LYMPH_NODES", "KIDNEYS", "LUNGS", "LIVER", "PLEURA", "OTHER", "PROSTATE", "ADRENAL", "BLADDER", "PERITONEUM", "COLON", "SOFT_TISSUE")
RF.4 <- c("TRT2_IDUNKNOWN", "TRT2_IDDOCETAXEL", "REGION_C", "SMOKEUNKNOWN", "SPINCOMP", "GLU")
RF.Sanofi <- c("TRT2_IDUNKNOWN", "TRT2_IDDOCETAXEL", "REGION_C", "SMOKEUNKNOWN", "GLU", "SPINCOMP")
RF.Ascent <- c("TRT2_IDDOCETAXEL", "TRT2_IDUNKNOWN", "SMOKEUNKNOWN", "REGION_C", "CREACL", "MHBLOOD")
RF.Celgene <- c("TRT2_IDDOCETAXEL", "TRT2_IDUNKNOWN", "REGION_C", "SMOKEUNKNOWN", "COPD", "NA.")
HC <- union(union(HC.4,HC.Sanofi), union(HC.Celgene,HC.Ascent))
RF <- union(union(RF.4,RF.Sanofi), union(RF.Ascent, RF.Celgene))
features <- use.all.features
features <- setdiff(features, union(c("TRT2_IDDOCETAXEL","TRT2_IDUNKNOWN"),setdiff(features,colnames(curr.training.data))))
##### UNIVARIATE (DEFAULT=DISABLED)
# sub.fs <- curr.training.data
# sub.fs$RPT <- NULL
# weights <- chi.squared(DISCONT ~., sub.fs)
# features <- cutoff.k(weights,k)
####### CLASSIFICATION / MODEL ##############################################################
x <- as.data.frame(curr.training.data[,features])
y <- as.factor(make.names(curr.training.data$DISCONT,unique = F))
set.seed(123)
## RF
model.rf <- randomForest(x=x, y=y, mtry=round(length(features)/rf.mtry), na.action = na.omit, probability=T,
ntree = (length(features)*rf.ntree),type="classification",replace = T)
prob <- predict(model.rf, curr.testing.data, type = "prob")[,2] # nb: 'class', rf: 'prob'
val.prob <- predict(model.rf, table.for.validation,type="prob")[,2] # nb: 'class', rf: 'prob'
# NB
# model.nb <- naiveBayes(x=x, y=y, type='prob')
# prob <- predict(model.nb, curr.testing.data, type = "raw")[,2] # nb: 'class', rf: 'prob'
# val.prob <- predict(model.nb, table.for.validation,type="raw")[,2] # nb: 'class', rf: 'prob'
# library(class)
# model.knn <- knn(curr.training.data[, features], curr.testing.data[, features], cl=factor(y), k = length(features), l = 0, prob = TRUE, use.all = TRUE)
# prob <- attr(model.knn, 'prob')
# model.knn <- knn(curr.training.data[, features], table.for.validation[, features], cl=factor(y), k = length(features), l = 0, prob = TRUE, use.all = TRUE)
# val.prob <- attr(model.knn, 'prob')
#
# model.nn <- nnet(as.formula(paste("y ~ ", paste(features, collapse= "+"))), data=cbind(x,y), linear.output=FALSE, likelihood=TRUE, size=length(features))
# prob <- predict(model.nn, curr.testing.data)
# val.prob <- predict(model.nn, table.for.validation)
#
model.svm <- svm(as.formula(paste("y ~ ", paste(features, collapse= "+"))), data=cbind(x,y), probability=TRUE)
prob <- attr(predict(model.svm, curr.testing.data, probability=TRUE), 'probabilities')[,'X1']
val.prob <- attr(predict(model.svm, table.for.validation, probability=TRUE), 'probabilities')[, 'X1']
####### OUTPUT ##############################################################
OUTPUT.TABLE <- cbind(curr.testing.data[,features],prob, as.factor(round(prob)), curr.testing.data$DISCONT)
SCORING.TABLE <- rbind(SCORING.TABLE, c(testing.name, dream9.score(prob, curr.testing.data$DISCONT)))
ACC <- round(score_q2(prob, curr.testing.data$DISCONT)*100)
# setwd("~/Dropbox/DREAM-F1000/OUTPUT")
## WRITE OUTPUT AS CSV
# write.csv(OUTPUT.TABLE, file = paste("OUTPUT/SUBSETTEST-",testing.name,"-",ACC,".csv", sep=""))
# print(paste("---- FILE:", "SUBSETTEST-",testing.name,"-",ACC,".csv WRITTEN TO HARDDRIVE ---",sep=""))
# write.csv(SCORING.TABLE, file = paste("OUTPUT/SCORE-",testing.name,"-",ACC,".csv", sep=""))
# print(paste("---- FILE:", "SCORE-",testing.name,"-",ACC,".csv WRITTEN TO HARDDRIVE ---",sep=""))
## WRITE VALIDATION OUTPUT AS CSV
FINAL.TABLE <- as.data.frame(cbind(as.character(table.for.validation$RPT),val.prob, as.numeric(round(val.prob))))
colnames (FINAL.TABLE) <- c("RPT","RISK","DISCONT")
# write.csv(FINAL.TABLE, file = paste("OUTPUT/VALIDATION-",testing.name,ACC,".csv", sep=""),row.names = FALSE)
# print(paste("---- FILE:", "VALIDATION-",ACC,".csv WRITTEN TO HARDDRIVE ---", sep=""))
} # END OF FOLD PER STUDY
} # END OF CROSS-FOLD VALIDATION
scoring.rows <- SCORING.TABLE[,1]
SCORING.TABLE <- SCORING.TABLE[,-1]
SCORING.TABLE <- apply(SCORING.TABLE,2,as.numeric)
row.names(SCORING.TABLE) <- scoring.rows
print(SCORING.TABLE)
plot(cbind(CV,SCORING.TABLE[,1]),xlab="Fold", ylab="Average AUC", main="Random forest importance cutoff")
print(paste("MEAN AUC: ", mean(SCORING.TABLE[,1])))
CV[which.max(SCORING.TABLE[,1])]
write.csv(SCORING.TABLE, file = paste("Data/NB/SCORE-svm-RF4-",testing.name,"-",ACC,".csv", sep=""))