machineLearning/cindex.R at master · EGates1/machineLearning · GitHub

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
library(compareC, quietly=TRUE)

dataset <- read.csv("file:///home/gpauloski/git-repos/TACE/gmmdatamatrix_July_28.csv")
models <- read.csv("file:///home/gpauloski/git-repos/TACE/modelPredictions.csv")

#dataset[,"liver_BCLC"] <- as.numeric(factor(dataset[,"liver_BCLC"], levels=levels(factor(dataset[,"liver_BCLC"]))))
#dataset[,"liver_CLIP"] <- as.numeric(factor(dataset[,"liver_CLIP"], levels=levels(factor(dataset[,"liver_CLIP"]))))
#dataset[,"liver_TNM"] <- as.numeric(factor(dataset[,"liver_TNM"], levels=levels(factor(dataset[,"liver_TNM"]))))
#dataset[,"liver_Okuda"] <- as.numeric(factor(dataset[,"liver_Okuda"], levels=levels(factor(dataset[,"liver_Okuda"]))))

#dataset[,"liver_BCLC"] <- -1*dataset[,"liver_BCLC"]
#dataset[,"liver_CLIP"] <- -1*dataset[,"liver_CLIP"]
#dataset[,"liver_TNM"] <- -1*dataset[,"liver_TNM"]
#dataset[,"liver_Okuda"] <- -1*dataset[,"liver_Okuda"]

#dataset[,"liver_BCLC"] <- ifelse(dataset[,"liver_BCLC"] < 4, 2, 1)
#dataset[,"liver_CLIP"] <- ifelse(dataset[,"liver_CLIP"] < 3, 2, 1)
#dataset[,"liver_TNM"] <- ifelse(dataset[,"liver_TNM"] < 3, 2, 1)
#dataset[,"liver_Okuda"] <- ifelse(dataset[,"liver_Okuda"] < 2, 2, 1)

cindex <- function(time, pred) {
	nobs <- length(time)
	pairs <- combn(nobs,2)
	conc <- 0
	tie <- 0

	for(i in 1:ncol(pairs)) {
		if(time[pairs[1,i]] > time[pairs[2,i]]) {
			if(pred[pairs[1,i]] == pred[pairs[2,i]])
				tie = tie + 1
			if(pred[pairs[1,i]] > pred[pairs[2,i]])
				conc = conc + 1
		} else if(time[pairs[1,i]] < time[pairs[2,i]]) {
			if(pred[pairs[1,i]] == pred[pairs[2,i]])
				tie = tie + 1
			if(pred[pairs[1,i]] < pred[pairs[2,i]])
				conc = conc + 1
		} else tie = tie + 1
	}

	return((conc+(0.5*tie))/ncol(pairs))
}

#cat("\nBCLC:", cindex(dataset[,"liver_TTP"],dataset[,"liver_BCLC"]))
#cat("\nCLIP:", cindex(dataset[,"liver_TTP"],dataset[,"liver_CLIP"]))
#cat("\nTNM:", cindex(dataset[,"liver_TTP"],dataset[,"liver_TNM"]))
#cat("\nOkuda:", cindex(dataset[,"liver_TTP"],dataset[,"liver_Okuda"]),"\n")

#cat("\nBCLC:", estC(dataset[,"liver_TTP"],dataset[,"liver_CensorStatus"],dataset[,"liver_BCLC"]))
#cat("\nCLIP:", estC(dataset[,"liver_TTP"],dataset[,"liver_CensorStatus"],dataset[,"liver_CLIP"]))
#cat("\nTNM:", estC(dataset[,"liver_TTP"],dataset[,"liver_CensorStatus"],dataset[,"liver_TNM"]))
#cat("\nOkuda:", estC(dataset[,"liver_TTP"],dataset[,"liver_CensorStatus"],dataset[,"liver_Okuda"]),"\n")

output <- data.frame(model=as.character(seq(1,ncol(models[,-1]),1)))
output <- cbind(output, c_index=seq(1,ncol(models[,-1]),1))
#output <- cbind(output, c_index2 = output$c_index)
for(i in 1:ncol(models[,-1])) {
	obs <- sort(c(which(models[,i+1]==1),which(models[,i+1]==2)))
	#output$c_index[i] <- estC(dataset[obs,"liver_TTP"],
	#	dataset[obs,"liver_CensorStatus"],models[obs,i+1])
	output$model[i] <- colnames(models)[i+1]
	output$c_index[i] <- cindex(dataset[obs,"liver_TTP"],models[obs,i+1])
}

output$model <- colnames(models)[2:ncol(models)]
write.csv(output, file="model_c_index.csv")
cat("\nC-Index values saved to model_c_index.csv\n\n")