xrobin / pROC
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 4
  X <- roc1$predictor
22 4
  Y <- roc2$predictor
23 4
  R <- rank(X, ties.method = ties.method)
24 4
  S <- rank(Y, ties.method = ties.method)
25 4
  D <- roc1$response # because roc1&roc2 are paired
26

27 4
  E <- venkatraman.paired.stat(R, S, D, roc1$levels)
28 4
  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 4
  return(list(E, EP))
30
}
31

32
venkatraman.unpaired.test <- function(roc1, roc2, boot.n, ties.method="first", progress, parallel) {
33 4
  X <- roc1$predictor
34 4
  Y <- roc2$predictor
35 4
  R <- rank(X, ties.method = ties.method)
36 4
  S <- rank(Y, ties.method = ties.method)
37 4
  D1<- roc1$response
38 4
  D2 <- roc2$response
39 4
  mp <- (sum(D1 == roc1$levels[2]) + sum(D2 == roc2$levels[2]))/(length(D1) + length(D1)) # mixing proportion, kappa
40

41 4
  E <- venkatraman.unpaired.stat(R, S, D1, D2, roc1$levels, roc2$levels, mp)
42 4
  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 4
  return(list(E, EP))
44
}
45

46
venkatraman.paired.permutation <- function(n, R, S, D, levels, ties.method) {
47
  # Break ties
48 4
  R2 <- R + runif(length(D)) - 0.5 # Add small amount of random but keep same mean
49 4
  S2 <- S + runif(length(D)) - 0.5
50

51
  # Permutation
52 4
  q <- 1 - round(runif(length(D)))
53 4
  R3 <- R2 * q + (1 - q) * S
54 4
  S3 <- S2 * q + (1 - q) * R
55

56 4
  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 4
  R <- R + runif(length(D1)) - 0.5 # Add small amount of random but keep same mean
62 4
  S <- S + runif(length(D2)) - 0.5
63

64 4
  R.controls <- R[D1==levels1[1]]
65 4
  R.cases <- R[D1==levels1[2]]
66 4
  S.controls <- S[D2==levels2[1]]
67 4
  S.cases <- S[D2==levels2[2]]
68

69
  # Permutation
70 4
  controls <- sample(c(R.controls, S.controls))
71 4
  cases <- sample(c(R.cases, S.cases))
72 4
  R[D1==levels1[1]] <- controls[1:length(R.controls)]
73 4
  S[D2==levels2[1]] <- controls[(length(R.controls)+1):length(controls)]
74 4
  R[D1==levels1[2]] <- cases[1:length(R.cases)]
75 4
  S[D2==levels2[2]] <- cases[(length(R.cases)+1):length(cases)]
76

77 4
  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 4
  R.controls <- R[D==levels[1]]
82 4
  R.cases <- R[D==levels[2]]
83 4
  S.controls <- S[D==levels[1]]
84 4
  S.cases <- S[D==levels[2]]
85 4
  n <- length(D)
86

87 4
  R.fn <- sapply(1:n, function(x) sum(R.cases <= x))
88 4
  R.fp <- sapply(1:n, function(x) sum(R.controls > x))
89 4
  S.fn <- sapply(1:n, function(x) sum(S.cases <= x))
90 4
  S.fp <- sapply(1:n, function(x) sum(S.controls > x))
91

92 4
  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 4
  R.controls <- R[D1==levels1[1]]
97 4
  R.cases <- R[D1==levels1[2]]
98 4
  S.controls <- S[D2==levels2[1]]
99 4
  S.cases <- S[D2==levels2[2]]
100 4
  n <- length(D1)
101 4
  m <- length(D2)
102

103 4
  R.fx <- sapply(1:n, function(x) sum(R.cases <= x)) / length(R.cases)
104 4
  R.gx <- sapply(1:n, function(x) sum(R.controls <= x)) / length(R.controls)
105 4
  S.fx <- sapply(1:m, function(x) sum(S.cases <= x)) / length(S.cases)
106 4
  S.gx <- sapply(1:m, function(x) sum(S.controls <= x)) / length(S.controls)
107 4
  R.p <- mp*R.fx + (1 - mp)*R.gx
108 4
  S.p <- mp*S.fx + (1 - mp)*S.gx
109 4
  R.exp <- mp*R.fx + (1 - mp)*(1-R.gx)
110 4
  S.exp <- mp*S.fx + (1 - mp)*(1-S.gx)
111

112
  # Do the integration
113 4
  x <- sort(c(R.p, S.p))
114 4
  R.f <- approxfun(R.p, R.exp)
115 4
  S.f <- approxfun(S.p, S.exp)
116 4
  f  <- function(x) abs(R.f(x)-S.f(x))
117 4
  y <- f(x)
118
  #trapezoid integration:
119 4
  idx <- 2:length(x)
120 4
  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 4
  return(integral)
122
}

Read our documentation on viewing source code .

Loading