|
1
|
|
# pROC: Tools Receiver operating characteristic (ROC curves) with
|
|
2
|
|
# (partial) area under the curve, confidence intervals and comparison.
|
|
3
|
|
# Copyright (C) 2010-2014 Xavier Robin, Alexandre Hainard, Natacha Turck,
|
|
4
|
|
# Natalia Tiberti, Frédérique Lisacek, Jean-Charles Sanchez
|
|
5
|
|
# and Markus Müller
|
|
6
|
|
#
|
|
7
|
|
# This program is free software: you can redistribute it and/or modify
|
|
8
|
|
# it under the terms of the GNU General Public License as published by
|
|
9
|
|
# the Free Software Foundation, either version 3 of the License, or
|
|
10
|
|
# (at your option) any later version.
|
|
11
|
|
#
|
|
12
|
|
# This program is distributed in the hope that it will be useful,
|
|
13
|
|
# but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
14
|
|
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
15
|
|
# GNU General Public License for more details.
|
|
16
|
|
#
|
|
17
|
|
# You should have received a copy of the GNU General Public License
|
|
18
|
|
# along with this program. If not, see <http://www.gnu.org/licenses/>.
|
|
19
|
|
|
|
20
|
|
venkatraman.paired.test <- function(roc1, roc2, boot.n, ties.method="first", progress, parallel) {
|
|
21
|
3
|
X <- roc1$predictor
|
|
22
|
3
|
Y <- roc2$predictor
|
|
23
|
3
|
R <- rank(X, ties.method = ties.method)
|
|
24
|
3
|
S <- rank(Y, ties.method = ties.method)
|
|
25
|
3
|
D <- roc1$response # because roc1&roc2 are paired
|
|
26
|
|
|
|
27
|
3
|
E <- venkatraman.paired.stat(R, S, D, roc1$levels)
|
|
28
|
3
|
EP <- laply(seq_len(boot.n), venkatraman.paired.permutation, R=R, S=S, D=D, levels=roc1$levels, ties.method=ties.method, .progress=progress, .parallel=parallel)
|
|
29
|
3
|
return(list(E, EP))
|
|
30
|
|
}
|
|
31
|
|
|
|
32
|
|
venkatraman.unpaired.test <- function(roc1, roc2, boot.n, ties.method="first", progress, parallel) {
|
|
33
|
3
|
X <- roc1$predictor
|
|
34
|
3
|
Y <- roc2$predictor
|
|
35
|
3
|
R <- rank(X, ties.method = ties.method)
|
|
36
|
3
|
S <- rank(Y, ties.method = ties.method)
|
|
37
|
3
|
D1<- roc1$response
|
|
38
|
3
|
D2 <- roc2$response
|
|
39
|
3
|
mp <- (sum(D1 == roc1$levels[2]) + sum(D2 == roc2$levels[2]))/(length(D1) + length(D1)) # mixing proportion, kappa
|
|
40
|
|
|
|
41
|
3
|
E <- venkatraman.unpaired.stat(R, S, D1, D2, roc1$levels, roc2$levels, mp)
|
|
42
|
3
|
EP <- laply(seq_len(boot.n), venkatraman.unpaired.permutation, R=R, S=S, D1=D1, D2=D2, levels1=roc1$levels, levels2=roc2$levels, mp=mp, ties.method=ties.method, .progress=progress, .parallel=parallel)
|
|
43
|
3
|
return(list(E, EP))
|
|
44
|
|
}
|
|
45
|
|
|
|
46
|
|
venkatraman.paired.permutation <- function(n, R, S, D, levels, ties.method) {
|
|
47
|
|
# Break ties
|
|
48
|
3
|
R2 <- R + runif(length(D)) - 0.5 # Add small amount of random but keep same mean
|
|
49
|
3
|
S2 <- S + runif(length(D)) - 0.5
|
|
50
|
|
|
|
51
|
|
# Permutation
|
|
52
|
3
|
q <- 1 - round(runif(length(D)))
|
|
53
|
3
|
R3 <- R2 * q + (1 - q) * S
|
|
54
|
3
|
S3 <- S2 * q + (1 - q) * R
|
|
55
|
|
|
|
56
|
3
|
return(venkatraman.paired.stat(rank(R3, ties.method=ties.method), rank(S3, ties.method=ties.method), D, levels))
|
|
57
|
|
}
|
|
58
|
|
|
|
59
|
|
venkatraman.unpaired.permutation <- function(n, R, S, D1, D2, levels1, levels2, mp, ties.method) {
|
|
60
|
|
# Break ties
|
|
61
|
3
|
R <- R + runif(length(D1)) - 0.5 # Add small amount of random but keep same mean
|
|
62
|
3
|
S <- S + runif(length(D2)) - 0.5
|
|
63
|
|
|
|
64
|
3
|
R.controls <- R[D1==levels1[1]]
|
|
65
|
3
|
R.cases <- R[D1==levels1[2]]
|
|
66
|
3
|
S.controls <- S[D2==levels2[1]]
|
|
67
|
3
|
S.cases <- S[D2==levels2[2]]
|
|
68
|
|
|
|
69
|
|
# Permutation
|
|
70
|
3
|
controls <- sample(c(R.controls, S.controls))
|
|
71
|
3
|
cases <- sample(c(R.cases, S.cases))
|
|
72
|
3
|
R[D1==levels1[1]] <- controls[1:length(R.controls)]
|
|
73
|
3
|
S[D2==levels2[1]] <- controls[(length(R.controls)+1):length(controls)]
|
|
74
|
3
|
R[D1==levels1[2]] <- cases[1:length(R.cases)]
|
|
75
|
3
|
S[D2==levels2[2]] <- cases[(length(R.cases)+1):length(cases)]
|
|
76
|
|
|
|
77
|
3
|
return(venkatraman.unpaired.stat(rank(R, ties.method=ties.method), rank(S, ties.method=ties.method), D1, D2, levels1, levels2, mp))
|
|
78
|
|
}
|
|
79
|
|
|
|
80
|
|
venkatraman.paired.stat <- function(R, S, D, levels) {
|
|
81
|
3
|
R.controls <- R[D==levels[1]]
|
|
82
|
3
|
R.cases <- R[D==levels[2]]
|
|
83
|
3
|
S.controls <- S[D==levels[1]]
|
|
84
|
3
|
S.cases <- S[D==levels[2]]
|
|
85
|
3
|
n <- length(D)
|
|
86
|
|
|
|
87
|
3
|
R.fn <- sapply(1:n, function(x) sum(R.cases <= x))
|
|
88
|
3
|
R.fp <- sapply(1:n, function(x) sum(R.controls > x))
|
|
89
|
3
|
S.fn <- sapply(1:n, function(x) sum(S.cases <= x))
|
|
90
|
3
|
S.fp <- sapply(1:n, function(x) sum(S.controls > x))
|
|
91
|
|
|
|
92
|
3
|
return(sum(abs((S.fn + S.fp) - (R.fn + R.fp))))
|
|
93
|
|
}
|
|
94
|
|
|
|
95
|
|
venkatraman.unpaired.stat <- function(R, S, D1, D2, levels1, levels2, mp) {
|
|
96
|
3
|
R.controls <- R[D1==levels1[1]]
|
|
97
|
3
|
R.cases <- R[D1==levels1[2]]
|
|
98
|
3
|
S.controls <- S[D2==levels2[1]]
|
|
99
|
3
|
S.cases <- S[D2==levels2[2]]
|
|
100
|
3
|
n <- length(D1)
|
|
101
|
3
|
m <- length(D2)
|
|
102
|
|
|
|
103
|
3
|
R.fx <- sapply(1:n, function(x) sum(R.cases <= x)) / length(R.cases)
|
|
104
|
3
|
R.gx <- sapply(1:n, function(x) sum(R.controls <= x)) / length(R.controls)
|
|
105
|
3
|
S.fx <- sapply(1:m, function(x) sum(S.cases <= x)) / length(S.cases)
|
|
106
|
3
|
S.gx <- sapply(1:m, function(x) sum(S.controls <= x)) / length(S.controls)
|
|
107
|
3
|
R.p <- mp*R.fx + (1 - mp)*R.gx
|
|
108
|
3
|
S.p <- mp*S.fx + (1 - mp)*S.gx
|
|
109
|
3
|
R.exp <- mp*R.fx + (1 - mp)*(1-R.gx)
|
|
110
|
3
|
S.exp <- mp*S.fx + (1 - mp)*(1-S.gx)
|
|
111
|
|
|
|
112
|
|
# Do the integration
|
|
113
|
3
|
x <- sort(c(R.p, S.p))
|
|
114
|
3
|
R.f <- approxfun(R.p, R.exp)
|
|
115
|
3
|
S.f <- approxfun(S.p, S.exp)
|
|
116
|
3
|
f <- function(x) abs(R.f(x)-S.f(x))
|
|
117
|
3
|
y <- f(x)
|
|
118
|
|
#trapezoid integration:
|
|
119
|
3
|
idx <- 2:length(x)
|
|
120
|
3
|
integral <- sum(((y[idx] + y[idx-1]) * (x[idx] - x[idx-1])) / 2, na.rm=TRUE) # remove NA that can appear in the borders
|
|
121
|
3
|
return(integral)
|
|
122
|
|
}
|
xrobin