1
# pROC: Tools Receiver operating characteristic (ROC curves) with
2
# (partial) area under the curve, confidence intervals and comparison. 
3
# Copyright (C) 2011-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
# Formula 3 from Obuchowski 2004, p. 1123
21
# Variance of an AUC given kappa
22
var.theta.obuchowski <- function(theta, kappa) {
23 29
    A <- qnorm(theta) * 1.414
24 29
    (0.0099 * exp(-A^2/2)) * ((5 * A^2 + 8) + (A^2 + 8)/kappa)
25
}
26

27
# Formulas from Obuchowski 1997, table 1 p. 1531
28
expr1 <- function(A, B) {
29 29
  exp(-A^2/(2 * (1 + B^2)))
30
}
31
expr2 <- function(B) {
32 29
  1 + B^2
33
}
34
expr3 <- function(cdagger1, cdagger2) {
35 29
  pnorm(cdagger1) - pnorm(cdagger2)
36
}
37
expr4 <- function(cddagger1, cddagger2) {
38 29
 exp(-cddagger1) - exp(- cddagger2)
39
}
40

41
cdagger <- function (A, B, FPRi) {
42 29
  (qnorm(FPRi) + A * B * (1 + B^2)^(-1) ) * (1 + B^2)^(1/2)
43
}
44

45
cddagger <- function(cdagger) {
46 29
  cdagger^2 / 2
47
}
48

49
f.full <- function(A, B) {
50 29
  expr1 <- expr1(A, B)
51 29
  expr2 <- expr2(B)
52 29
  expr1 * (2 * pi * expr2) ^ (-1/2)
53
}
54

55
f.partial <- function(A, B, FPR1, FPR2) {
56 29
  cdagger1 <- cdagger(A, B, FPR1)
57 29
  cdagger2 <- cdagger(A, B, FPR2)
58 29
  expr1 <- expr1(A, B)
59 29
  expr2 <- expr2(B)
60 29
  expr3 <- expr3(cdagger1, cdagger2)
61 29
  expr1 * (2 * pi * expr2) ^ (-1/2) * expr3
62
}
63

64
g.full <- function(A, B) {
65 29
  expr1 <- expr1(A, B)
66 29
  expr2 <- expr2(B)
67 29
  - expr1 * A * B * (2 * pi * expr2^3) ^ (-1/2)
68
}
69

70
g.partial <- function(A, B, FPR1, FPR2) {
71 29
  cdagger1 <- cdagger(A, B, FPR1)
72 29
  cdagger2 <- cdagger(A, B, FPR2)
73 29
  cddagger1 <- cddagger(cdagger1)
74 29
  cddagger2 <- cddagger(cdagger2)
75 29
  expr1 <- expr1(A, B)
76 29
  expr2 <- expr2(B)
77 29
  expr3 <- expr3(cdagger1, cdagger2)
78 29
  expr4 <- expr4(cddagger1, cddagger2)
79
  # WARNING: we have set (-expr4), in contradiction with Obuchowski paper
80 29
  expr1 * (2 * pi * expr2) ^ (-1) * (-expr4) - A * B * expr1 * (2 * pi * expr2^3) ^ (-1/2) * expr3
81
}
82

83
# Variance of a ROC curve given a 'roc' object
84
var.roc.obuchowski <- function(roc) {
85 29
  binormal <- smooth(roc, method="binormal")$model
86 29
  A <- unname(coefficients(binormal)[1])
87 29
  B <- unname(coefficients(binormal)[2])
88 29
  kappa <- length(roc$controls) / length(roc$cases)
89

90 29
  if (!identical(attr(roc$auc, "partial.auc"), FALSE)) {
91 29
    FPR1 <- 1 - attr(roc$auc, "partial.auc")[2]
92 29
    FPR2 <- 1 - attr(roc$auc, "partial.auc")[1]
93 29
    va <- var.params.obuchowski(A, B, kappa, FPR1, FPR2)
94
  }
95
  else {
96 29
    va <- var.params.obuchowski(A, B, kappa)
97
  }
98 29
  return(va)
99
}
100

101
# Variance of a ROC curve given the parameters
102
# Obuchowski 1997, formula 4 p. 1530
103
# A and B: params of the binormal ROC curve
104
# kappa: proportion controls / cases
105
# FPR1, FPR2: the bottom (1) or top (2) bounds of the pAUC interval
106
var.params.obuchowski <- function(A, B, kappa, FPR1, FPR2) {
107 29
  if (!missing(FPR1) && !is.null(FPR1) && !missing(FPR1) && !is.null(FPR2)) {
108 29
    f.partial(A, B, FPR1, FPR2)^2 * (1 + B^2 / kappa + A^2/2) + g.partial(A, B, FPR1, FPR2)^2 * B^2 * (1 + kappa) / (2*kappa)
109
  }
110
  else {
111 29
    f.full(A, B)^2 * (1 + B^2 / kappa + A^2/2) + g.full(A, B)^2 * B^2 * (1 + kappa) / (2*kappa)
112
  }
113
}
114

115
# Covariance of 2 given 'roc' objects (under the alternative hypothesis)
116
cov.roc.obuchowski <- function(roc1, roc2) {
117 29
  binormal1 <- smooth(roc1, method="binormal")$model
118 29
  A1 <- unname(coefficients(binormal1)[1])
119 29
  B1 <- unname(coefficients(binormal1)[2])
120 29
  binormal2 <- smooth(roc2, method="binormal")$model
121 29
  A2 <-unname(coefficients(binormal2)[1])
122 29
  B2 <- unname(coefficients(binormal2)[2])
123 29
  kappa <- length(roc1$controls) / length(roc1$cases)
124 29
  ra <- cor(roc1$cases, roc2$cases)
125 29
  rn <- cor(roc1$controls, roc2$controls)
126 29
  if (!identical(attr(roc1$auc, "partial.auc"), FALSE)) {
127 29
    FPR11 <- 1 - attr(roc1$auc, "partial.auc")[2]
128 29
    FPR12 <- 1 - attr(roc1$auc, "partial.auc")[1]
129 29
    FPR21 <- 1 - attr(roc2$auc, "partial.auc")[2]
130 29
    FPR22 <- 1 - attr(roc2$auc, "partial.auc")[1]
131 29
    co <- cov.params.obuchowski(A1, B1, A2, B2, rn, ra, kappa, FPR11, FPR12, FPR21, FPR22)
132
  }
133
  else {
134 29
    co <- cov.params.obuchowski(A1, B1, A2, B2, rn, ra, kappa)
135
  }
136 29
  return(co)
137
}
138

139
# Covariance under the null hypothesis
140
# roc1 is taken as null
141
cov0.roc.obuchowski <- function(roc1, roc2) {
142 0
  binormal <- smooth(roc, method="binormal")$model
143 0
  A <- unname(coefficients(binormal)[1])
144 0
  B <- unname(coefficients(binormal)[2])
145 0
  R <- length(roc1$controls) / length(roc1$cases)
146 0
  ra <- cor(roc1$cases, roc2$cases)
147 0
  rn <- cor(roc1$controls, roc2$controls)
148 0
  if (!identical(attr(roc1$auc, "partial.auc"), FALSE)) {
149 0
    FPR1 <- attr(roc1$auc, "partial.auc")[2]
150 0
    FPR2 <- attr(roc1$auc, "partial.auc")[1]
151 0
    co <- cov.params.obuchowski(A, B, A, B, rn, ra, kappa, FPR1, FPR2, FPR1, FPR2)
152
  }
153
  else {
154 0
    co <- cov.params.obuchowski(A, B, A, B, rn, ra, kappa)
155
  }
156 0
  return(co)
157
}
158

159

160
# Covariance of a ROC curve given the parameters
161
# Obuchowski 1997, formula 5 p. 1531
162
# (A|B)(1|2): A and B params of the binormal ROC curve
163
# rn, ra: correlation of the results in ROC curves 1 and 2 in controls (n) and cases (a) patients
164
# kappa: proportion controls / cases
165
# FPR(1|2)(1|2): the bounds of the pAUC interval:
166
#    ***** ROC curve 1 or 2
167
#         ***** bottom (1) or top (2) of the interval
168
cov.params.obuchowski <- function(A1, B1, A2, B2, rn, ra, kappa, FPR11, FPR12, FPR21, FPR22) {
169 29
  if (!missing(FPR11) && !is.null(FPR11) &&
170 29
      !missing(FPR12) && !is.null(FPR12) &&
171 29
      !missing(FPR21) && !is.null(FPR21) &&
172 29
      !missing(FPR22) && !is.null(FPR22)) {
173 29
    f1 <- f.partial(A1, B1, FPR11, FPR12)
174 29
    f2 <- f.partial(A2, B2, FPR21, FPR22)
175 29
    g1 <- g.partial(A1, B1, FPR11, FPR12)
176 29
    g2 <- g.partial(A2, B2, FPR21, FPR22)
177
  }
178
  else {
179 29
    f1 <- f.full(A1, B1)
180 29
    f2 <- f.full(A2, B2)
181 29
    g1 <- g.full(A1, B1)
182 29
    g2 <- g.full(A2, B2)
183
  }
184 29
  f1 * f2 * (ra + rn * B1 * B2 / kappa + ra^2 * A1 * A2  / 2) +
185 29
    g1 * g2 * (B1 * B2 * (rn^2 + kappa * ra^2) / (2 *  kappa)) + 
186 29
    f1 * g2 * (ra^2 * A1 * B2 / 2) + f2 * g1 * (ra^2 * A2 * B1 / 2)
187
}
188

189
# Variance of a difference between two ROC curves given the parameters
190
# Obuchowski 1997, formula 4 and 5 p. 1530--1531
191
# (A|B)(1|2): A and B params of the binormal ROC curve
192
# rn, ra: correlation of the results in ROC curves 1 and 2 in controls (n) and cases (a) patients
193
# kappa: proportion controls / cases
194
# FPR(1|2)(1|2): the bounds of the pAUC interval:
195
#    ***** ROC curve 1 or 2
196
#         ***** bottom (1) or top (2) of the interval
197
vardiff.params.obuchowski <- function(A1, B1, A2, B2, rn, ra, kappa, FPR11, FPR12, FPR21, FPR22) {
198 0
  var(A1, B1, kappa, FPR11, FPR12) + var(A2, B2, kappa, FPR21, FPR22) +
199 0
    2 * cov(A1, B1, A2, B2, rn, ra, kappa, FPR11, FPR12, FPR21, FPR22)
200
}
201

202
# Variance of a difference between two ROC curves given the parameters
203
# under the null hypothesis. ROC curve 1 is taken as null
204
# Obuchowski 1997, formula 4 and 5 p. 1530--1531
205
# (A|B)(1|2): A and B params of the binormal ROC curve
206
# rn, ra: correlation of the results in ROC curves 1 and 2 in controls (n) and cases (a) patients
207
# kappa: proportion controls / cases
208
# FPR(1|2)(1|2): the bounds of the pAUC interval:
209
#    ***** ROC curve 1 or 2
210
#         ***** bottom (1) or top (2) of the interval
211
vardiff0.params.obuchowski <- function(A1, B1, A2, B2, rn, ra, kappa, FPR11, FPR12, FPR21, FPR22) {
212 0
  2 * var(A1, B1, kappa, FPR11, FPR12) + 
213 0
    2 * cov(A1, B1, A2, B2, rn, ra, kappa, FPR11, FPR12, FPR21, FPR22)
214
}

Read our documentation on viewing source code .

Loading