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 . ``` 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 .