xrobin / pROC
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 8
power.roc.test <- function(...)
21 8
  UseMethod("power.roc.test")
22

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

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

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

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

102
      # Restore sig.level if two.sided
103 8
      if (alternative == "two.sided") {
104 8
        sig.level <- sig.level * 2
105
      }
106 8
      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"))
107
    }
108
    else {
109 0
      stop("'roc2' must be an object of class 'roc'.")
110
    }
111
  }
112
  else {
113 8
  	ncontrols <- length(roc1$controls)
114 8
  	ncases <- length(roc1$cases)
115 8
    if (! is.null(sig.level) && ! is.null(power)) {
116 8
      if (is.null(kappa)) {
117 8
        kappa <- ncontrols / ncases
118
      }
119 8
      ncontrols <- ncases <- NULL
120
    }
121 8
    auc <- auc(roc1)
122
    # TODO: implement this with var() and cov() for the given ROC curve
123 8
    return(power.roc.test.numeric(ncontrols = ncontrols, ncases = ncases, auc = auc, sig.level = sig.level, power = power, alternative = alternative, kappa = kappa, ...))
124
  }
125
}
126

127
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"), ...) {
128
  # basic sanity checks
129 8
  if (!is.null(ncases) && ncases < 0)
130 0
    stop("'ncases' must be positive")
131 8
  if (!is.null(ncontrols) && ncontrols < 0)
132 0
    stop("'ncontrols' must be positive")
133 8
  if (!is.null(kappa) && kappa < 0)
134 0
    stop("'kappa' must be positive")
135 8
  if (!is.null(power) && (power < 0 || power > 1))
136 0
    stop("'power' must range from 0 to 1")
137 8
  if (!is.null(sig.level) && (sig.level < 0 || sig.level > 1))
138 0
    stop("'sig.level' must range from 0 to 1")
139
  
140
  # Complete ncontrols and ncases with kappa
141 8
  if (is.null(ncontrols) && ! is.null(ncases) && !is.null(kappa))
142 0
    ncontrols <- kappa * ncases
143 8
  else if (is.null(ncases) && ! is.null(ncontrols) && !is.null(kappa))
144 0
    ncases <- ncontrols / kappa
145
  
146 8
  alternative <- match.arg(alternative)
147 8
  if (alternative == "two.sided" && !is.null(sig.level)) {
148 8
    sig.level <- sig.level / 2
149
  }
150

151
  # determine AUC
152 8
  if (is.null(auc)) {
153 8
    if (is.null(ncontrols) || is.null(ncases))
154 0
      stop("'ncontrols' and 'ncases' (or one of these with 'kappa') or 'auc' must be provided.")
155 8
    else if (is.null(power))
156 0
      stop("'power' or 'auc' must be provided.")
157 8
    else if (is.null(sig.level))
158 0
      stop("'sig.level' or 'auc' must be provided.")
159 8
    kappa <- ncontrols / ncases
160 8
    zalpha <- qnorm(1 - sig.level)
161 8
    zbeta <- qnorm(power)
162

163 8
    tryCatch(
164 8
             root <- uniroot(power.roc.test.optimize.auc.function, interval=c(0.5, 1-1e-16), ncontrols=ncontrols, ncases=ncases, zalpha=zalpha, zbeta=zbeta),
165 0
             error=function(e) {stop(sprintf("AUC could not be solved:\n%s", e))}
166
             )
167 8
    auc <- root$root
168
  }
169

170
  # Determine number of patients (sample size)
171 8
  else if (is.null(ncases) && is.null(ncontrols)) {
172 8
    if (is.null(power))
173 0
      stop("'ncontrols' and 'ncases' (or one of these with 'kappa') or 'auc' must be provided.")
174 8
    else if (is.null(kappa))
175 0
      stop("'kappa' must be provided.")
176 8
    else if (is.null(auc))
177 0
      stop("'auc' or 'ncases' and 'ncontrols' must be provided.")
178 8
    else if (is.null(sig.level))
179 0
      stop("'sig.level' or 'ncases' and 'ncontrols' must be provided.")
180

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

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

243
power.roc.test.list <- function(parslist, ncontrols = NULL, ncases = NULL, sig.level = 0.05, power = NULL,  kappa = 1, alternative = c("two.sided", "one.sided"), ...) {
244
  # basic sanity checks
245 8
  if (!is.null(ncases) && ncases < 0)
246 0
    stop("'ncases' must be positive")
247 8
  if (!is.null(ncontrols) && ncontrols < 0)
248 0
    stop("'ncontrols' must be positive")
249 8
  if (!is.null(kappa) && kappa < 0)
250 0
    stop("'kappa' must be positive")
251 8
  if (!is.null(power) && (power < 0 || power > 1))
252 0
    stop("'power' must range from 0 to 1")
253 8
  if (!is.null(sig.level) && (sig.level < 0 || sig.level > 1))
254 0
    stop("'sig.level' must range from 0 to 1")
255

256
  
257
  # Complete ncontrols and ncases with kappa
258 8
  if (is.null(ncontrols) && ! is.null(ncases) && !is.null(kappa))
259 8
    ncontrols <- kappa * ncases
260 8
  else if (is.null(ncases) && ! is.null(ncontrols) && !is.null(kappa))
261 0
    ncases <- ncontrols / kappa
262

263
  # Warn if anything is passed with ...
264 8
  if (length(list(...)) > 0) {
265 0
    warning(paste("The following arguments were ignored:", paste(names(list(...)), collapse=", ")))
266
  }
267
  
268 8
  alternative <- match.arg(alternative)
269 8
  if (alternative == "two.sided" && !is.null(sig.level)) {
270 8
    sig.level <- sig.level / 2
271
  }
272

273
  # Check required elements of parslist
274 8
  required <- c("A1", "B1", "A2", "B2", "rn", "ra", "delta")
275 8
  if (any(! required %in% names(parslist))) {
276 0
    stop(paste("Missing parameter(s):", paste(required[! required %in% names(parslist) ], collapse=", ")))
277
  }
278

279
  # Determine number of patients (sample size)
280 8
  if (is.null(ncases) && is.null(ncontrols)) {
281 8
    if (is.null(power))
282 0
      stop("'ncontrols' and 'ncases' (or one of these with 'kappa') or 'auc' must be provided.")
283 8
    else if (is.null(kappa))
284 0
      stop("'kappa' must be provided.")
285 8
    else if (is.null(sig.level))
286 0
      stop("'sig.level' or 'ncases' and 'ncontrols' must be provided.")
287

288 8
    zalpha <- qnorm(1 - sig.level)
289 8
    zbeta <- qnorm(power)
290 8
    ncases <- ncases.obuchowski.params(parslist, zalpha, zbeta, kappa)
291 8
    ncontrols <- kappa * ncases
292
  }
293
  
294
  # Determine power
295 8
  else if (is.null(power)) {
296 8
    if (is.null(ncontrols) || is.null(ncases))
297 0
      stop("'ncontrols' and 'ncases' (or one of these with 'kappa') or 'auc' must be provided.")
298 8
    else if (is.null(sig.level))
299 0
      stop("'sig.level' or 'power' must be provided.")
300 8
    kappa <- ncontrols / ncases
301

302 8
    zalpha <- qnorm(1 - sig.level)
303 8
    zbeta <- zbeta.obuchowski.params(parslist, zalpha, ncases, kappa)
304 8
    power <- pnorm(zbeta)
305
  }
306

307
  # Determine sig.level
308 8
  else  if (is.null(sig.level)) { 
309 8
    if (is.null(ncontrols) || is.null(ncases))
310 0
      stop("'ncontrols' and 'ncases' (or one of these with 'kappa') or 'auc' must be provided.")
311 8
    else if (is.null(power))
312 0
      stop("'power' or 'sig.level' must be provided.")
313 8
    kappa <- ncontrols / ncases
314

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

328
}
329

330

331
#### HIDDEN FUNCTIONS ####
332

333
# A function to 'optimize' auc
334
power.roc.test.optimize.auc.function <- function(x, ncontrols, ncases, zalpha, zbeta) {
335 8
  kappa <- ncontrols / ncases
336 8
  Vtheta <- var.theta.obuchowski(x, kappa)
337 8
  (zalpha * sqrt(0.0792 * (1 + 1/kappa)) + zbeta * sqrt(Vtheta))^2 / (x - 0.5)^2 - ncases
338
}
339

340
# Compute variance of a delta from a 'covvar' list (see 'covvar' below)
341
var.delta.covvar <- function(covvar) {
342 8
  covvar$var1 + covvar$var2 - 2 * covvar$cov12
343
}
344

345
# Compute variance of a delta from a 'covvar' list (see 'covvar' below)
346
# under the null hypothesis
347
# roc1 taken as reference.
348
var0.delta.covvar <- function(covvar) {
349 8
  2 * covvar$var1 - 2 * covvar$cov12
350
}
351

352
# Compute the number of cases with Obuchowski formula and var(... method=method)
353
ncases.obuchowski <- function(roc1, roc2, zalpha, zbeta, method, ...) {
354 8
  delta <- roc1$auc - roc2$auc
355 8
  covvar <- covvar(roc1, roc2, method, ...)
356 8
  v0 <- var0.delta.covvar(covvar)
357 8
  va <- var.delta.covvar(covvar)
358 8
  nd <- solve.nd(zalpha = zalpha,
359 8
  			   zbeta = zbeta,
360 8
  			   v0 = v0, va = va,
361 8
  			   delta = delta)
362 8
  return(nd)
363
}
364

365
# Compute the number of cases with Obuchowski formula from params
366
ncases.obuchowski.params <- function(parslist, zalpha, zbeta, kappa) {
367 8
  covvar <- list(
368 8
                 var1 = var.params.obuchowski(parslist$A1, parslist$B1, kappa, parslist$FPR11, parslist$FPR12),
369 8
                 var2 = var.params.obuchowski(parslist$A2, parslist$B2, kappa, parslist$FPR21, parslist$FPR22),
370 8
                 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)
371
                 )
372 8
  v0 <- var0.delta.covvar(covvar)
373 8
  va <- var.delta.covvar(covvar)
374 8
  nd <- solve.nd(zalpha = zalpha,
375 8
  			   zbeta = zbeta,
376 8
  			   v0 = v0, va = va,
377 8
  			   delta = parslist$delta)
378 8
  return(nd)
379
}
380

381
# Compute the z alpha with Obuchowski formula and var(... method=method)
382
zalpha.obuchowski <- function(roc1, roc2, zbeta, method, ...) {
383 8
  delta <- roc1$auc - roc2$auc
384 8
  ncases <- length(roc1$cases)
385 8
  covvar <- covvar(roc1, roc2, method, ...)
386 8
  v0 <- var0.delta.covvar(covvar)
387 8
  va <- var.delta.covvar(covvar)
388 8
  zalpha <- solve.zalpha(nd=ncases,
389 8
  					   zbeta = zbeta,
390 8
  					   v0 = v0, va = va,
391 8
  					   delta = delta)
392 8
  return(zalpha)
393
}
394

395
# Compute the z alpha with Obuchowski formula from params
396
zalpha.obuchowski.params <- function(parslist, zbeta, ncases, kappa) {
397 8
  covvar <- list(
398 8
                 var1 = var.params.obuchowski(parslist$A1, parslist$B1, kappa, parslist$FPR11, parslist$FPR12),
399 8
                 var2 = var.params.obuchowski(parslist$A2, parslist$B2, kappa, parslist$FPR21, parslist$FPR22),
400 8
                 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)
401
                 )
402 8
  v0 <- var0.delta.covvar(covvar)
403 8
  va <- var.delta.covvar(covvar)
404 8
  zalpha <- solve.zalpha(nd=ncases,
405 8
  					   zbeta = zbeta,
406 8
  					   v0 = v0, va = va,
407 8
  					   delta = parslist$delta)
408 8
  return(zalpha)
409
}
410

411
# Compute the z beta with Obuchowski formula and var(... method=method)
412
zbeta.obuchowski <- function(roc1, roc2, zalpha, method, ...) {
413 8
  delta <- roc1$auc - roc2$auc
414 8
  ncases <- length(roc1$cases)
415 8
  covvar <- covvar(roc1, roc2, method, ...)
416 8
  v0 <- var0.delta.covvar(covvar)
417 8
  va <- var.delta.covvar(covvar)
418 8
  zbeta <- solve.zbeta(nd=ncases,
419 8
  					 zalpha = zalpha,
420 8
  					 v0 = v0, va = va,
421 8
  					 delta = delta)
422 8
  return(zbeta)
423
}
424

425
# Compute the z beta with Obuchowski formula from params
426
zbeta.obuchowski.params <- function(parslist, zalpha, ncases, kappa) {
427 8
	covvar <- list(
428 8
		var1 = var.params.obuchowski(parslist$A1, parslist$B1, kappa, parslist$FPR11, parslist$FPR12),
429 8
		var2 = var.params.obuchowski(parslist$A2, parslist$B2, kappa, parslist$FPR21, parslist$FPR22),
430 8
		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)
431
	)
432 8
	v0 <- var0.delta.covvar(covvar)
433 8
	va <- var.delta.covvar(covvar)
434 8
	a <- va
435 8
	zbeta <- solve.zbeta(nd=ncases,
436 8
						 zalpha = zalpha,
437 8
						 v0 = v0, va = va,
438 8
						 delta = parslist$delta)
439 8
	return(zbeta)
440
}
441

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

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

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

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

Read our documentation on viewing source code .

Loading