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 33
power.roc.test <- function(...)
21 33
  UseMethod("power.roc.test")
22

23
power.roc.test.roc <- function(roc1, roc2, sig.level = 0.05, power = NULL, alternative = c("two.sided", "one.sided"), reuse.auc=TRUE, method=c("delong", "bootstrap", "obuchowski"), ...) {
24
  # Basic sanity checks
25 33
  if (!is.null(power) && (power < 0 || power > 1))
26 0
    stop("'power' must range from 0 to 1")
27 33
  if (!is.null(sig.level) && (sig.level < 0 || sig.level > 1))
28 0
    stop("'sig.level' must range from 0 to 1")
29
  
30
  # check that the AUC of roc1 was computed, or do it now
31 33
  if (is.null(roc1$auc) | !reuse.auc) {
32 0
    roc1$auc <- auc(roc1, ...)
33
  }
34 33
  if (!is.null(attr(roc1$auc, "partial.auc.correct")) && attr(roc1$auc, "partial.auc.correct")) {
35 0
    stop("Cannot compute power with corrected partial AUCs")
36
  }
37 33
  roc1 <- roc.utils.unpercent(roc1)
38

39 33
  if (!missing(roc2) && !is.null(roc2)) {
40 33
    alternative <- match.arg(alternative)
41 33
    if (!is.null(sig.level) && alternative == "two.sided") {
42 33
      sig.level <- sig.level / 2
43
    }
44
    
45 33
    if ("roc" %in% class(roc2)) {
46
      # check that the AUC of roc2 was computed, or do it now
47 33
      if (is.null(roc2$auc) | !reuse.auc) {
48 0
        roc2$auc <- auc(roc2, ...)
49
      }
50 33
      if (!is.null(attr(roc2$auc, "partial.auc.correct")) && attr(roc2$auc, "partial.auc.correct")) {
51 0
        stop("Cannot compute power with corrected partial AUCs")
52
      }
53 33
      roc2 <- roc.utils.unpercent(roc2)
54

55
      # Make sure the ROC curves are paired
56 33
      rocs.are.paired <- are.paired(roc1, roc2)
57 33
      if (!rocs.are.paired) {
58 0
        stop("The sample size for a difference in AUC cannot be applied to unpaired ROC curves yet.")
59
      }
60
      # Make sure the AUC specifications are identical
61 33
      attr1 <- attributes(roc1$auc); attr1$roc <- NULL
62 33
      attr2 <- attributes(roc2$auc); attr2$roc <- NULL
63 33
      if (!identical(attr1, attr2)) {
64 0
        stop("Different AUC specifications in the ROC curves.")
65
      }
66
    
67
      # check that the same region was requested in auc. Otherwise, issue a warning
68 33
      if (!identical(attributes(roc1$auc)[names(attributes(roc1$auc))!="roc"], attributes(roc2$auc)[names(attributes(roc2$auc))!="roc"]))
69 0
        warning("Different AUC specifications in the ROC curves. Enforcing the inconsistency, but unexpected results may be produced.")
70
 
71 33
      ncontrols <- length(roc1$controls)
72 33
      ncases <- length(roc1$cases)
73 33
      kappa <- ncontrols / ncases
74

75
      # Power test
76 33
      if (is.null(power)) {
77 33
        if (is.null(sig.level))
78 0
          stop("'sig.level' or 'power' must be provided.")
79 33
        zalpha <- qnorm(1 - sig.level)
80 33
        zbeta <- zbeta.obuchowski(roc1, roc2, zalpha, method=method, ...)
81 33
        power <- pnorm(zbeta)
82
      }
83
      # sig.level
84 33
      else if (is.null(sig.level)) {
85 33
        zbeta <- qnorm(power)
86 33
        zalpha <- zalpha.obuchowski(roc1, roc2, zbeta, method=method, ...)
87 33
        sig.level <- 1 - pnorm(zalpha)
88
      }
89
      # Sample size
90
      else {
91 33
        zalpha <- qnorm(1 - sig.level)
92 33
        zbeta <- qnorm(power)
93 33
        ncases <- ncases.obuchowski(roc1, roc2, zalpha, zbeta, method=method, ...)
94 33
        ncontrols <- kappa * ncases
95
      }
96

97
      # Restore sig.level if two.sided
98 33
      if (alternative == "two.sided") {
99 33
        sig.level <- sig.level * 2
100
      }
101 33
      return(structure(list(ncases=ncases, ncontrols=ncontrols, auc1=roc1$auc, auc2=roc2$auc, sig.level=sig.level, power=power, alternative=alternative, method="Two ROC curves power calculation"), class="power.htest"))
102
    }
103
    else {
104 0
      stop("'roc2' must be an object of class 'roc'.")
105
    }
106
  }
107
  else {
108 33
    if (is.null(sig.level) || is.null(power)) {
109 33
      ncontrols <- length(roc1$controls)
110 33
      ncases <- length(roc1$cases)
111
    }
112
    else {
113 0
      ncontrols <- ncases <- NULL
114
    }
115 33
    auc <- auc(roc1)
116
    # TODO: implement this with var() and cov() for the given ROC curve
117 33
    return(power.roc.test.numeric(ncontrols = ncontrols, ncases = ncases, auc = auc, sig.level = sig.level, power = power, alternative = alternative, ...))
118
  }
119
}
120

121
power.roc.test.numeric <- function(auc = NULL, ncontrols = NULL, ncases = NULL, sig.level = 0.05, power = NULL,  kappa = 1, alternative = c("two.sided", "one.sided"), ...) {
122
  # basic sanity checks
123 33
  if (!is.null(ncases) && ncases < 0)
124 0
    stop("'ncases' must be positive")
125 33
  if (!is.null(ncontrols) && ncontrols < 0)
126 0
    stop("'ncontrols' must be positive")
127 33
  if (!is.null(kappa) && kappa < 0)
128 0
    stop("'kappa' must be positive")
129 33
  if (!is.null(power) && (power < 0 || power > 1))
130 0
    stop("'power' must range from 0 to 1")
131 33
  if (!is.null(sig.level) && (sig.level < 0 || sig.level > 1))
132 0
    stop("'sig.level' must range from 0 to 1")
133
  
134
  # Complete ncontrols and ncases with kappa
135 33
  if (is.null(ncontrols) && ! is.null(ncases) && !is.null(kappa))
136 0
    ncontrols <- kappa * ncases
137 33
  else if (is.null(ncases) && ! is.null(ncontrols) && !is.null(kappa))
138 0
    ncases <- ncontrols / kappa
139
  
140 33
  alternative <- match.arg(alternative)
141 33
  if (alternative == "two.sided" && !is.null(sig.level)) {
142 33
    sig.level <- sig.level / 2
143
  }
144

145
  # determine AUC
146 33
  if (is.null(auc)) {
147 33
    if (is.null(ncontrols) || is.null(ncases))
148 0
      stop("'ncontrols' and 'ncases' (or one of these with 'kappa') or 'auc' must be provided.")
149 33
    else if (is.null(power))
150 0
      stop("'power' or 'auc' must be provided.")
151 33
    else if (is.null(sig.level))
152 0
      stop("'sig.level' or 'auc' must be provided.")
153 33
    kappa <- ncontrols / ncases
154 33
    zalpha <- qnorm(1 - sig.level)
155 33
    zbeta <- qnorm(power)
156

157 33
    tryCatch(
158 33
             root <- uniroot(power.roc.test.optimize.auc.function, interval=c(0.5, 1-1e-16), ncontrols=ncontrols, ncases=ncases, zalpha=zalpha, zbeta=zbeta),
159 0
             error=function(e) {stop(sprintf("AUC could not be solved:\n%s", e))}
160
             )
161 33
    auc <- root$root
162
  }
163

164
  # Determine number of patients (sample size)
165 33
  else if (is.null(ncases) && is.null(ncontrols)) {
166 33
    if (is.null(power))
167 0
      stop("'ncontrols' and 'ncases' (or one of these with 'kappa') or 'auc' must be provided.")
168 33
    else if (is.null(kappa))
169 0
      stop("'kappa' must be provided.")
170 33
    else if (is.null(auc))
171 0
      stop("'auc' or 'ncases' and 'ncontrols' must be provided.")
172 33
    else if (is.null(sig.level))
173 0
      stop("'sig.level' or 'ncases' and 'ncontrols' must be provided.")
174

175 33
    theta <- as.numeric(auc)
176 33
    Vtheta <- var.theta.obuchowski(theta, kappa)
177
    
178 33
    ncases <- solve.nd(zalpha = qnorm(1 - sig.level),
179 33
    				   zbeta = qnorm(power),
180 33
    				   v0 = 0.0792 * (1 + 1 / kappa),
181 33
    				   va = Vtheta,
182 33
    				   delta = theta - 0.5)
183 33
    ncontrols <- kappa * ncases
184
  }
185
  
186
  # Determine power
187 33
  else if (is.null(power)) { 
188 33
    if (is.null(ncontrols) || is.null(ncases))
189 0
      stop("'ncontrols' and 'ncases' (or one of these with 'kappa') or 'auc' must be provided.")
190 33
    else if (is.null(auc))
191 0
      stop("'auc' or 'power' must be provided.")
192 33
    else if (is.null(sig.level))
193 0
      stop("'sig.level' or 'power' must be provided.")
194
  	
195 33
    kappa <- ncontrols / ncases
196 33
    theta <- as.numeric(auc)
197 33
    Vtheta <- var.theta.obuchowski(theta, kappa)
198
    
199 33
    zbeta <- solve.zbeta(nd = ncases,
200 33
    					 zalpha = qnorm(1 - sig.level),
201 33
    					 v0 = 0.0792 * (1 + 1 / kappa),
202 33
    					 va = Vtheta,
203 33
    					 delta = theta - 0.5)
204 33
    power <- pnorm(zbeta)
205
  }
206

207
  # Determine sig.level
208 33
  else  if (is.null(sig.level)) { 
209 33
    if (is.null(ncontrols) || is.null(ncases))
210 0
      stop("'ncontrols' and 'ncases' (or one of these with 'kappa') or 'auc' must be provided.")
211 33
    else if (is.null(auc))
212 0
      stop("'auc' or 'sig.level' must be provided.")
213 33
    else if (is.null(power))
214 0
      stop("'power' or 'sig.level' must be provided.")
215
  	
216 33
    kappa <- ncontrols / ncases
217 33
    theta <- as.numeric(auc)
218 33
    Vtheta <- var.theta.obuchowski(theta, kappa)
219
    
220 33
    zalpha <- solve.zalpha(nd = ncases,
221 33
    					  zbeta = qnorm(power),
222 33
    					  v0 = 0.0792 * (1 + 1 / kappa),
223 33
    					  va = Vtheta,
224 33
    					  delta = theta - 0.5)
225 33
    sig.level <- 1 - pnorm(zalpha)
226
  }
227
  else {
228 0
    stop("One of 'power', 'sig.level', 'auc', or both 'ncases' and 'ncontrols' must be NULL.")
229
  }
230
  # Restore sig.level if two.sided
231 33
  if (alternative == "two.sided") {
232 33
    sig.level <- sig.level * 2
233
  }
234 33
  return(structure(list(ncases=ncases, ncontrols=ncontrols, auc=auc, sig.level=sig.level, power=power, method="One ROC curve power calculation"), class="power.htest"))
235
}
236

237
power.roc.test.list <- function(parslist, ncontrols = NULL, ncases = NULL, sig.level = 0.05, power = NULL,  kappa = 1, alternative = c("two.sided", "one.sided"), ...) {
238
  # basic sanity checks
239 33
  if (!is.null(ncases) && ncases < 0)
240 0
    stop("'ncases' must be positive")
241 33
  if (!is.null(ncontrols) && ncontrols < 0)
242 0
    stop("'ncontrols' must be positive")
243 33
  if (!is.null(kappa) && kappa < 0)
244 0
    stop("'kappa' must be positive")
245 33
  if (!is.null(power) && (power < 0 || power > 1))
246 0
    stop("'power' must range from 0 to 1")
247 33
  if (!is.null(sig.level) && (sig.level < 0 || sig.level > 1))
248 0
    stop("'sig.level' must range from 0 to 1")
249

250
  
251
  # Complete ncontrols and ncases with kappa
252 33
  if (is.null(ncontrols) && ! is.null(ncases) && !is.null(kappa))
253 33
    ncontrols <- kappa * ncases
254 33
  else if (is.null(ncases) && ! is.null(ncontrols) && !is.null(kappa))
255 0
    ncases <- ncontrols / kappa
256

257
  # Warn if anything is passed with ...
258 33
  if (length(list(...)) > 0) {
259 0
    warning(paste("The following arguments were ignored:", paste(names(list(...)), collapse=", ")))
260
  }
261
  
262 33
  alternative <- match.arg(alternative)
263 33
  if (alternative == "two.sided" && !is.null(sig.level)) {
264 33
    sig.level <- sig.level / 2
265
  }
266

267
  # Check required elements of parslist
268 33
  required <- c("A1", "B1", "A2", "B2", "rn", "ra", "delta")
269 33
  if (any(! required %in% names(parslist))) {
270 0
    stop(paste("Missing parameter(s):", paste(required[! required %in% names(parslist) ], collapse=", ")))
271
  }
272

273
  # Determine number of patients (sample size)
274 33
  if (is.null(ncases) && is.null(ncontrols)) {
275 33
    if (is.null(power))
276 0
      stop("'ncontrols' and 'ncases' (or one of these with 'kappa') or 'auc' must be provided.")
277 33
    else if (is.null(kappa))
278 0
      stop("'kappa' must be provided.")
279 33
    else if (is.null(sig.level))
280 0
      stop("'sig.level' or 'ncases' and 'ncontrols' must be provided.")
281

282 33
    zalpha <- qnorm(1 - sig.level)
283 33
    zbeta <- qnorm(power)
284 33
    ncases <- ncases.obuchowski.params(parslist, zalpha, zbeta, kappa)
285 33
    ncontrols <- kappa * ncases
286
  }
287
  
288
  # Determine power
289 33
  else if (is.null(power)) {
290 33
    if (is.null(ncontrols) || is.null(ncases))
291 0
      stop("'ncontrols' and 'ncases' (or one of these with 'kappa') or 'auc' must be provided.")
292 33
    else if (is.null(sig.level))
293 0
      stop("'sig.level' or 'power' must be provided.")
294 33
    kappa <- ncontrols / ncases
295

296 33
    zalpha <- qnorm(1 - sig.level)
297 33
    zbeta <- zbeta.obuchowski.params(parslist, zalpha, ncases, kappa)
298 33
    power <- pnorm(zbeta)
299
  }
300

301
  # Determine sig.level
302 33
  else  if (is.null(sig.level)) { 
303 33
    if (is.null(ncontrols) || is.null(ncases))
304 0
      stop("'ncontrols' and 'ncases' (or one of these with 'kappa') or 'auc' must be provided.")
305 33
    else if (is.null(power))
306 0
      stop("'power' or 'sig.level' must be provided.")
307 33
    kappa <- ncontrols / ncases
308

309 33
    zbeta <- qnorm(power)
310 33
    zalpha <- zalpha.obuchowski.params(parslist, zbeta, ncases, kappa)
311 33
    sig.level <- 1 - pnorm(zalpha)
312
  }
313
  else {
314 0
    stop("One of 'power', 'sig.level', 'auc', or both 'ncases' and 'ncontrols' must be NULL.")
315
  }
316
  # Restore sig.level if two.sided
317 33
  if (alternative == "two.sided") {
318 33
    sig.level <- sig.level * 2
319
  }
320 33
  return(structure(list(ncases=ncases, ncontrols=ncontrols, sig.level=sig.level, power=power, method="Two ROC curves power calculation"), class="power.htest"))
321

322
}
323

324

325
#### HIDDEN FUNCTIONS ####
326

327
# A function to 'optimize' auc
328
power.roc.test.optimize.auc.function <- function(x, ncontrols, ncases, zalpha, zbeta) {
329 33
  kappa <- ncontrols / ncases
330 33
  Vtheta <- var.theta.obuchowski(x, kappa)
331 33
  (zalpha * sqrt(0.0792 * (1 + 1/kappa)) + zbeta * sqrt(Vtheta))^2 / (x - 0.5)^2 - ncases
332
}
333

334
# Compute variance of a delta from a 'covvar' list (see 'covvar' below)
335
var.delta.covvar <- function(covvar) {
336 33
  covvar$var1 + covvar$var2 - 2 * covvar$cov12
337
}
338

339
# Compute variance of a delta from a 'covvar' list (see 'covvar' below)
340
# under the null hypothesis
341
# roc1 taken as reference.
342
var0.delta.covvar <- function(covvar) {
343 33
  2 * covvar$var1 - 2 * covvar$cov12
344
}
345

346
# Compute the number of cases with Obuchowski formula and var(... method=method)
347
ncases.obuchowski <- function(roc1, roc2, zalpha, zbeta, method, ...) {
348 33
  delta <- roc1$auc - roc2$auc
349 33
  covvar <- covvar(roc1, roc2, method, ...)
350 33
  v0 <- var0.delta.covvar(covvar)
351 33
  va <- var.delta.covvar(covvar)
352 33
  nd <- solve.nd(zalpha = zalpha,
353 33
  			   zbeta = zbeta,
354 33
  			   v0 = v0, va = va,
355 33
  			   delta = delta)
356 33
  return(nd)
357
}
358

359
# Compute the number of cases with Obuchowski formula from params
360
ncases.obuchowski.params <- function(parslist, zalpha, zbeta, kappa) {
361 33
  covvar <- list(
362 33
                 var1 = var.params.obuchowski(parslist$A1, parslist$B1, kappa, parslist$FPR11, parslist$FPR12),
363 33
                 var2 = var.params.obuchowski(parslist$A2, parslist$B2, kappa, parslist$FPR21, parslist$FPR22),
364 33
                 cov12 = cov.params.obuchowski(parslist$A1, parslist$B1, parslist$A2, parslist$B2, parslist$rn, parslist$ra, kappa, parslist$FPR11, parslist$FPR12, parslist$FPR21, parslist$FPR22)
365
                 )
366 33
  v0 <- var0.delta.covvar(covvar)
367 33
  va <- var.delta.covvar(covvar)
368 33
  nd <- solve.nd(zalpha = zalpha,
369 33
  			   zbeta = zbeta,
370 33
  			   v0 = v0, va = va,
371 33
  			   delta = parslist$delta)
372 33
  return(nd)
373
}
374

375
# Compute the z alpha with Obuchowski formula and var(... method=method)
376
zalpha.obuchowski <- function(roc1, roc2, zbeta, method, ...) {
377 33
  delta <- roc1$auc - roc2$auc
378 33
  ncases <- length(roc1$cases)
379 33
  covvar <- covvar(roc1, roc2, method, ...)
380 33
  v0 <- var0.delta.covvar(covvar)
381 33
  va <- var.delta.covvar(covvar)
382 33
  zalpha <- solve.zalpha(nd=ncases,
383 33
  					   zbeta = zbeta,
384 33
  					   v0 = v0, va = va,
385 33
  					   delta = delta)
386 33
  return(zalpha)
387
}
388

389
# Compute the z alpha with Obuchowski formula from params
390
zalpha.obuchowski.params <- function(parslist, zbeta, ncases, kappa) {
391 33
  covvar <- list(
392 33
                 var1 = var.params.obuchowski(parslist$A1, parslist$B1, kappa, parslist$FPR11, parslist$FPR12),
393 33
                 var2 = var.params.obuchowski(parslist$A2, parslist$B2, kappa, parslist$FPR21, parslist$FPR22),
394 33
                 cov12 = cov.params.obuchowski(parslist$A1, parslist$B1, parslist$A2, parslist$B2, parslist$rn, parslist$ra, kappa, parslist$FPR11, parslist$FPR12, parslist$FPR21, parslist$FPR22)
395
                 )
396 33
  v0 <- var0.delta.covvar(covvar)
397 33
  va <- var.delta.covvar(covvar)
398 33
  zalpha <- solve.zalpha(nd=ncases,
399 33
  					   zbeta = zbeta,
400 33
  					   v0 = v0, va = va,
401 33
  					   delta = parslist$delta)
402 33
  return(zalpha)
403
}
404

405
# Compute the z beta with Obuchowski formula and var(... method=method)
406
zbeta.obuchowski <- function(roc1, roc2, zalpha, method, ...) {
407 33
  delta <- roc1$auc - roc2$auc
408 33
  ncases <- length(roc1$cases)
409 33
  covvar <- covvar(roc1, roc2, method, ...)
410 33
  v0 <- var0.delta.covvar(covvar)
411 33
  va <- var.delta.covvar(covvar)
412 33
  zbeta <- solve.zbeta(nd=ncases,
413 33
  					 zalpha = zalpha,
414 33
  					 v0 = v0, va = va,
415 33
  					 delta = delta)
416 33
  return(zbeta)
417
}
418

419
# Compute the z beta with Obuchowski formula from params
420
zbeta.obuchowski.params <- function(parslist, zalpha, ncases, kappa) {
421 33
	covvar <- list(
422 33
		var1 = var.params.obuchowski(parslist$A1, parslist$B1, kappa, parslist$FPR11, parslist$FPR12),
423 33
		var2 = var.params.obuchowski(parslist$A2, parslist$B2, kappa, parslist$FPR21, parslist$FPR22),
424 33
		cov12 = cov.params.obuchowski(parslist$A1, parslist$B1, parslist$A2, parslist$B2, parslist$rn, parslist$ra, kappa, parslist$FPR11, parslist$FPR12, parslist$FPR21, parslist$FPR22)
425
	)
426 33
	v0 <- var0.delta.covvar(covvar)
427 33
	va <- var.delta.covvar(covvar)
428 33
	a <- va
429 33
	zbeta <- solve.zbeta(nd=ncases,
430 33
						 zalpha = zalpha,
431 33
						 v0 = v0, va = va,
432 33
						 delta = parslist$delta)
433 33
	return(zbeta)
434
}
435

436
solve.zbeta <- function(nd, zalpha, v0, va, delta) {
437
	# Solve for z_\beta in Obuchowski formula:
438
	# See formula 2 in Obuchowsk & McClish 1997  (2 ROC curves)
439
	# or formula 2 in Obuchowski et al 2004 (1 ROC curve)
440
	# The formula is of the form:
441
	# nd = (z_alpha * sqrt(v0) - z_beta * sqrt(va)) / delta ^ 2
442
	# Re-organized:
443
	# z_beta = (sqrt(nd * delta ^ 2) - z_alpha * sqrt(v0)) / sqrt(va)
444
	# @param nd: number of diseased patients (or abornmal, N_A in Obuchowsk & McClish 1997)
445
	# @param zalpha: upper \alpha (sig.level) percentile of the standard normal distribution
446
	# @param v0 the null variance associated with z_alpha
447
	# @param va: the alternative variance associated with z_beta
448
	# @param delta: the difference in AUC
449 33
	return((sqrt(nd * delta ^ 2) - zalpha * sqrt(v0)) / sqrt(va))
450
}
451

452
solve.nd <- function(zalpha, zbeta, v0, va, delta) {
453
	# Solve for number of diseased (abnormal) patients in Obuchowski formula:
454
	# See formula 2 in Obuchowsk & McClish 1997  (2 ROC curves)
455
	# or formula 2 in Obuchowski et al 2004 (1 ROC curve)
456
	# nd = (z_alpha * sqrt(v0) - z_beta * sqrt(va)) / delta ^ 2
457
	# @param zalpha: upper \alpha (sig.level) percentile of the standard normal distribution
458
	# @param zbeta: upper \beta (power) percentile of the standard normal distribution
459
	# @param v0 the null variance associated with z_alpha
460
	# @param va: the alternative variance associated with z_beta
461
	# @param delta: the difference in AUC
462 33
	return((zalpha * sqrt(v0) + zbeta * sqrt(va)) ^ 2 / delta ^ 2)
463
}
464

465
solve.zalpha <- function(nd, zbeta, v0, va, delta) {
466
	# Solve for z_\alpha in Obuchowski formula:
467
	# See formula 2 in Obuchowsk & McClish 1997  (2 ROC curves)
468
	# or formula 2 in Obuchowski et al 2004 (1 ROC curve)
469
	# The formula is of the form:
470
	# nd = (z_alpha * sqrt(v0) - z_beta * sqrt(va)) / delta ^ 2
471
	# Re-organized:
472
	# z_alpha = (sqrt(nd * delta ^ 2) - z_beta * sqrt(va)) / sqrt(v0)
473
	# @param nd: number of diseased patients (or abornmal, N_A in Obuchowsk & McClish 1997)
474
	# @param zbeta: upper \beta (power) percentile of the standard normal distribution
475
	# @param v0 the null variance associated with z_alpha
476
	# @param va: the alternative variance associated with z_beta
477
	# @param delta: the difference in AUC
478 33
	return((sqrt(nd * delta ^ 2) - zbeta * sqrt(va)) / sqrt(v0))
479
}
480

481
# Compute var and cov of two ROC curves by bootstrap in a single bootstrap run
482
covvar <- function(roc1, roc2, method, ...) {
483 33
  cov12 <- cov(roc1, roc2, boot.return=TRUE, method=method, ...)
484 33
  if (!is.null(attr(cov12, "resampled.values"))) {
485 0
    var1 <- var(attr(cov12, "resampled.values")[,1])
486 0
    var2 <- var(attr(cov12, "resampled.values")[,2])
487 0
    attr(cov12, "resampled.values") <- NULL
488
  }
489
  else {
490 33
    var1 <- var(roc1, method=method, ...)
491 33
    var2 <- var(roc2, method=method, ...)
492
  }
493 33
  ncases <- length(roc1$cases)
494 33
  return(list(var1 = var1 * ncases, var2 = var2 * ncases, cov12 = cov12 * ncases))
495
}

Read our documentation on viewing source code .

Loading