obuchowski.R at 641f842 ⋅ xrobin/pROC

# pROC: Tools Receiver operating characteristic (ROC curves) with
# (partial) area under the curve, confidence intervals and comparison. 
# Copyright (C) 2011-2014 Xavier Robin, Alexandre Hainard, Natacha Turck,
# Natalia Tiberti, Frédérique Lisacek, Jean-Charles Sanchez
# and Markus Müller
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

# Formula 3 from Obuchowski 2004, p. 1123
# Variance of an AUC given kappa
var.theta.obuchowski <- function(theta, kappa) {
    A <- qnorm(theta) * 1.414
    (0.0099 * exp(-A^2/2)) * ((5 * A^2 + 8) + (A^2 + 8)/kappa)
}

# Formulas from Obuchowski 1997, table 1 p. 1531
expr1 <- function(A, B) {
  exp(-A^2/(2 * (1 + B^2)))
}
expr2 <- function(B) {
  1 + B^2
}
expr3 <- function(cdagger1, cdagger2) {
  pnorm(cdagger1) - pnorm(cdagger2)
}
expr4 <- function(cddagger1, cddagger2) {
 exp(-cddagger1) - exp(- cddagger2)
}

cdagger <- function (A, B, FPRi) {
  (qnorm(FPRi) + A * B * (1 + B^2)^(-1) ) * (1 + B^2)^(1/2)
}

cddagger <- function(cdagger) {
  cdagger^2 / 2
}

f.full <- function(A, B) {
  expr1 <- expr1(A, B)
  expr2 <- expr2(B)
  expr1 * (2 * pi * expr2) ^ (-1/2)
}

f.partial <- function(A, B, FPR1, FPR2) {
  cdagger1 <- cdagger(A, B, FPR1)
  cdagger2 <- cdagger(A, B, FPR2)
  expr1 <- expr1(A, B)
  expr2 <- expr2(B)
  expr3 <- expr3(cdagger1, cdagger2)
  expr1 * (2 * pi * expr2) ^ (-1/2) * expr3
}

g.full <- function(A, B) {
  expr1 <- expr1(A, B)
  expr2 <- expr2(B)
  - expr1 * A * B * (2 * pi * expr2^3) ^ (-1/2)
}

g.partial <- function(A, B, FPR1, FPR2) {
  cdagger1 <- cdagger(A, B, FPR1)
  cdagger2 <- cdagger(A, B, FPR2)
  cddagger1 <- cddagger(cdagger1)
  cddagger2 <- cddagger(cdagger2)
  expr1 <- expr1(A, B)
  expr2 <- expr2(B)
  expr3 <- expr3(cdagger1, cdagger2)
  expr4 <- expr4(cddagger1, cddagger2)
  # WARNING: we have set (-expr4), in contradiction with Obuchowski paper
  expr1 * (2 * pi * expr2) ^ (-1) * (-expr4) - A * B * expr1 * (2 * pi * expr2^3) ^ (-1/2) * expr3
}

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

  if (!identical(attr(roc$auc, "partial.auc"), FALSE)) {
    FPR1 <- 1 - attr(roc$auc, "partial.auc")[2]
    FPR2 <- 1 - attr(roc$auc, "partial.auc")[1]
    va <- var.params.obuchowski(A, B, kappa, FPR1, FPR2)
  }
  else {
    va <- var.params.obuchowski(A, B, kappa)
  }
  return(va)
}

# Variance of a ROC curve given the parameters
# Obuchowski 1997, formula 4 p. 1530
# A and B: params of the binormal ROC curve
# kappa: proportion controls / cases
# FPR1, FPR2: the bottom (1) or top (2) bounds of the pAUC interval
var.params.obuchowski <- function(A, B, kappa, FPR1, FPR2) {
  if (!missing(FPR1) && !is.null(FPR1) && !missing(FPR1) && !is.null(FPR2)) {
    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)
  }
  else {
    f.full(A, B)^2 * (1 + B^2 / kappa + A^2/2) + g.full(A, B)^2 * B^2 * (1 + kappa) / (2*kappa)
  }
}

# Covariance of 2 given 'roc' objects (under the alternative hypothesis)
cov.roc.obuchowski <- function(roc1, roc2) {
  binormal1 <- smooth(roc1, method="binormal")$model
  A1 <- unname(coefficients(binormal1)[1])
  B1 <- unname(coefficients(binormal1)[2])
  binormal2 <- smooth(roc2, method="binormal")$model
  A2 <-unname(coefficients(binormal2)[1])
  B2 <- unname(coefficients(binormal2)[2])
  kappa <- length(roc1$controls) / length(roc1$cases)
  ra <- cor(roc1$cases, roc2$cases)
  rn <- cor(roc1$controls, roc2$controls)
  if (!identical(attr(roc1$auc, "partial.auc"), FALSE)) {
    FPR11 <- 1 - attr(roc1$auc, "partial.auc")[2]
    FPR12 <- 1 - attr(roc1$auc, "partial.auc")[1]
    FPR21 <- 1 - attr(roc2$auc, "partial.auc")[2]
    FPR22 <- 1 - attr(roc2$auc, "partial.auc")[1]
    co <- cov.params.obuchowski(A1, B1, A2, B2, rn, ra, kappa, FPR11, FPR12, FPR21, FPR22)
  }
  else {
    co <- cov.params.obuchowski(A1, B1, A2, B2, rn, ra, kappa)
  }
  return(co)
}

# Covariance under the null hypothesis
# roc1 is taken as null
cov0.roc.obuchowski <- function(roc1, roc2) {
  binormal <- smooth(roc, method="binormal")$model
  A <- unname(coefficients(binormal)[1])
  B <- unname(coefficients(binormal)[2])
  R <- length(roc1$controls) / length(roc1$cases)
  ra <- cor(roc1$cases, roc2$cases)
  rn <- cor(roc1$controls, roc2$controls)
  if (!identical(attr(roc1$auc, "partial.auc"), FALSE)) {
    FPR1 <- attr(roc1$auc, "partial.auc")[2]
    FPR2 <- attr(roc1$auc, "partial.auc")[1]
    co <- cov.params.obuchowski(A, B, A, B, rn, ra, kappa, FPR1, FPR2, FPR1, FPR2)
  }
  else {
    co <- cov.params.obuchowski(A, B, A, B, rn, ra, kappa)
  }
  return(co)
}


# Covariance of a ROC curve given the parameters
# Obuchowski 1997, formula 5 p. 1531
# (A|B)(1|2): A and B params of the binormal ROC curve
# rn, ra: correlation of the results in ROC curves 1 and 2 in controls (n) and cases (a) patients
# kappa: proportion controls / cases
# FPR(1|2)(1|2): the bounds of the pAUC interval:
#    ***** ROC curve 1 or 2
#         ***** bottom (1) or top (2) of the interval
cov.params.obuchowski <- function(A1, B1, A2, B2, rn, ra, kappa, FPR11, FPR12, FPR21, FPR22) {
  if (!missing(FPR11) && !is.null(FPR11) &&
      !missing(FPR12) && !is.null(FPR12) &&
      !missing(FPR21) && !is.null(FPR21) &&
      !missing(FPR22) && !is.null(FPR22)) {
    f1 <- f.partial(A1, B1, FPR11, FPR12)
    f2 <- f.partial(A2, B2, FPR21, FPR22)
    g1 <- g.partial(A1, B1, FPR11, FPR12)
    g2 <- g.partial(A2, B2, FPR21, FPR22)
  }
  else {
    f1 <- f.full(A1, B1)
    f2 <- f.full(A2, B2)
    g1 <- g.full(A1, B1)
    g2 <- g.full(A2, B2)
  }
  f1 * f2 * (ra + rn * B1 * B2 / kappa + ra^2 * A1 * A2  / 2) +
    g1 * g2 * (B1 * B2 * (rn^2 + kappa * ra^2) / (2 *  kappa)) + 
    f1 * g2 * (ra^2 * A1 * B2 / 2) + f2 * g1 * (ra^2 * A2 * B1 / 2)
}

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

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

Read our documentation on viewing source code .

Navigation	Overlay
`t` Navigate files	`h` Toggle hits
`y` Change url to tip of branch	`m` Toggle misses
`b / v` Jump to prev/next hit line	`p` Toggle partial
`z / x` Jump to prev/next missed or partial line	`1..9` Toggle flags
`shift + o` Open current page in GitHub	`a` Toggle all on
`/ or ?` Show keyboard shortcuts dialog	`c` Toggle context lines or commits

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	3	A <- qnorm(theta) * 1.414
24	3	(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	3	exp(-A^2/(2 * (1 + B^2)))
30		}
31		expr2 <- function(B) {
32	3	1 + B^2
33		}
34		expr3 <- function(cdagger1, cdagger2) {
35	3	pnorm(cdagger1) - pnorm(cdagger2)
36		}
37		expr4 <- function(cddagger1, cddagger2) {
38	3	exp(-cddagger1) - exp(- cddagger2)
39		}
40
41		cdagger <- function (A, B, FPRi) {
42	3	(qnorm(FPRi) + A * B * (1 + B^2)^(-1) ) * (1 + B^2)^(1/2)
43		}
44
45		cddagger <- function(cdagger) {
46	3	cdagger^2 / 2
47		}
48
49		f.full <- function(A, B) {
50	3	expr1 <- expr1(A, B)
51	3	expr2 <- expr2(B)
52	3	expr1 * (2 * pi * expr2) ^ (-1/2)
53		}
54
55		f.partial <- function(A, B, FPR1, FPR2) {
56	3	cdagger1 <- cdagger(A, B, FPR1)
57	3	cdagger2 <- cdagger(A, B, FPR2)
58	3	expr1 <- expr1(A, B)
59	3	expr2 <- expr2(B)
60	3	expr3 <- expr3(cdagger1, cdagger2)
61	3	expr1 * (2 * pi * expr2) ^ (-1/2) * expr3
62		}
63
64		g.full <- function(A, B) {
65	3	expr1 <- expr1(A, B)
66	3	expr2 <- expr2(B)
67	3	- expr1 * A * B * (2 * pi * expr2^3) ^ (-1/2)
68		}
69
70		g.partial <- function(A, B, FPR1, FPR2) {
71	3	cdagger1 <- cdagger(A, B, FPR1)
72	3	cdagger2 <- cdagger(A, B, FPR2)
73	3	cddagger1 <- cddagger(cdagger1)
74	3	cddagger2 <- cddagger(cdagger2)
75	3	expr1 <- expr1(A, B)
76	3	expr2 <- expr2(B)
77	3	expr3 <- expr3(cdagger1, cdagger2)
78	3	expr4 <- expr4(cddagger1, cddagger2)
79		# WARNING: we have set (-expr4), in contradiction with Obuchowski paper
80	3	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	3	binormal <- smooth(roc, method="binormal")$model
86	3	A <- unname(coefficients(binormal)[1])
87	3	B <- unname(coefficients(binormal)[2])
88	3	kappa <- length(roc$controls) / length(roc$cases)
89
90	3	if (!identical(attr(roc$auc, "partial.auc"), FALSE)) {
91	3	FPR1 <- 1 - attr(roc$auc, "partial.auc")[2]
92	3	FPR2 <- 1 - attr(roc$auc, "partial.auc")[1]
93	3	va <- var.params.obuchowski(A, B, kappa, FPR1, FPR2)
94		}
95		else {
96	3	va <- var.params.obuchowski(A, B, kappa)
97		}
98	3	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	3	if (!missing(FPR1) && !is.null(FPR1) && !missing(FPR1) && !is.null(FPR2)) {
108	3	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	3	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	3	binormal1 <- smooth(roc1, method="binormal")$model
118	3	A1 <- unname(coefficients(binormal1)[1])
119	3	B1 <- unname(coefficients(binormal1)[2])
120	3	binormal2 <- smooth(roc2, method="binormal")$model
121	3	A2 <-unname(coefficients(binormal2)[1])
122	3	B2 <- unname(coefficients(binormal2)[2])
123	3	kappa <- length(roc1$controls) / length(roc1$cases)
124	3	ra <- cor(roc1$cases, roc2$cases)
125	3	rn <- cor(roc1$controls, roc2$controls)
126	3	if (!identical(attr(roc1$auc, "partial.auc"), FALSE)) {
127	3	FPR11 <- 1 - attr(roc1$auc, "partial.auc")[2]
128	3	FPR12 <- 1 - attr(roc1$auc, "partial.auc")[1]
129	3	FPR21 <- 1 - attr(roc2$auc, "partial.auc")[2]
130	3	FPR22 <- 1 - attr(roc2$auc, "partial.auc")[1]
131	3	co <- cov.params.obuchowski(A1, B1, A2, B2, rn, ra, kappa, FPR11, FPR12, FPR21, FPR22)
132		}
133		else {
134	3	co <- cov.params.obuchowski(A1, B1, A2, B2, rn, ra, kappa)
135		}
136	3	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	3	if (!missing(FPR11) && !is.null(FPR11) &&
170	3	!missing(FPR12) && !is.null(FPR12) &&
171	3	!missing(FPR21) && !is.null(FPR21) &&
172	3	!missing(FPR22) && !is.null(FPR22)) {
173	3	f1 <- f.partial(A1, B1, FPR11, FPR12)
174	3	f2 <- f.partial(A2, B2, FPR21, FPR22)
175	3	g1 <- g.partial(A1, B1, FPR11, FPR12)
176	3	g2 <- g.partial(A2, B2, FPR21, FPR22)
177		}
178		else {
179	3	f1 <- f.full(A1, B1)
180	3	f2 <- f.full(A2, B2)
181	3	g1 <- g.full(A1, B1)
182	3	g2 <- g.full(A2, B2)
183		}
184	3	f1 * f2 * (ra + rn * B1 * B2 / kappa + ra^2 * A1 * A2 / 2) +
185	3	g1 * g2 * (B1 * B2 * (rn^2 + kappa * ra^2) / (2 * kappa)) +
186	3	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		}