xrobin / pROC
1
# pROC: Tools Receiver operating characteristic (ROC curves) with
2
# (partial) area under the curve, confidence intervals and comparison. 
3
# Copyright (C) 2010-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
smooth <- function(...)
21 8
  UseMethod("smooth")
22

23
smooth.default <- function(...) {
24 8
  stats::smooth(...)
25
}
26

27
smooth.smooth.roc <- function(smooth.roc, ...) {
28 0
  roc <- attr(smooth.roc, "roc")
29 0
  if (is.null(roc))
30 0
    stop("Cannot smooth a ROC curve generated directly with numeric 'density.controls' and 'density.cases'.")
31 0
  smooth.roc(roc, ...)
32
}
33

34
smooth.roc <- function(roc, method = c("binormal", "density", "fitdistr", "logcondens", "logcondens.smooth"), n = 512, bw = "nrd0",
35
                       density = NULL, density.controls = density, density.cases = density,
36
                       start = NULL, start.controls = start, start.cases = start,
37
                       reuse.auc = TRUE, reuse.ci = FALSE,
38
                       ...) {
39 8
  method <- match.arg(method)
40
  
41 8
  if (is.ordered(roc$original.predictor) && (method == "density" || method =="fitidstr"))
42 0
    stop("ROC curves of ordered predictors can be smoothed only with binormal smoothing.")
43

44 8
  if (method == "binormal") {
45 8
    sesp <- smooth.roc.binormal(roc, n)
46
  }
47 8
  else if (method == "fitdistr") {
48 8
    sesp <- smooth.roc.fitdistr(roc, n, density.controls, density.cases, start.controls, start.cases, ...)
49
  }
50 8
  else if (method == "density") {
51 8
    sesp <- smooth.roc.density(roc, n, density.controls, density.cases, bw, ...)
52
  }
53 8
  else if (method == "logcondens") {
54 8
    sesp <- smooth.roc.logcondens(roc, n)
55
  }
56 8
  else if (method == "logcondens.smooth") {
57 8
    sesp <- smooth.roc.logcondens.smooth(roc, n)
58
  }
59
  else {
60 0
    stop(sprintf("Impossible smooth method value '%s'. Please report this bug to <%s>.",
61 0
                 method, utils::packageDescription("pROC")$BugReports))
62
  }
63

64 8
  class(sesp) <- "smooth.roc"
65 8
  sesp <- sort(sesp) # sort se and sp
66
  # anchor SE/SP at 0/100
67 8
  sesp$specificities <- c(0, as.vector(sesp$specificities), ifelse(roc$percent, 100, 1))
68 8
  sesp$sensitivities <- c(ifelse(roc$percent, 100, 1), as.vector(sesp$sensitivities), 0)
69 8
  attr(sesp, "roc") <- roc # keep the original roc. May be useful in CI.
70 8
  sesp$percent <- roc$percent # keep some basic roc specifications
71 8
  sesp$direction <- roc$direction
72 8
  sesp$call <- match.call()
73
  # keep smoothing arguments (for print and bootstrap)
74 8
  sesp$smoothing.args <- list(...)
75 8
  sesp$smoothing.args$method <- method
76 8
  sesp$smoothing.args$n <- n
77 8
  sesp$smoothing.args$start.controls <- start.controls
78 8
  sesp$smoothing.args$start.cases <- start.cases
79 8
  sesp$smoothing.args$density.controls <- density.controls
80 8
  sesp$smoothing.args$density.cases <- density.cases
81 8
  sesp$smoothing.args$bw <- bw
82
  # complete fit.controls/cases if a function was passed as densfun
83 8
  if (method == "fitdistr") {
84 8
    if (is.null(sesp$fit.controls$densfun)) {
85 8
      if (missing(density.controls))
86 0
        sesp$fit.controls$densfun <- deparse(substitute(density))
87
      else
88 8
        sesp$fit.controls$densfun <- deparse(substitute(density.controls))
89
    }
90 8
    if (is.null(sesp$fit.cases$densfun)) {
91 8
      if (missing(density.cases))
92 0
        sesp$fit.cases$densfun <- deparse(substitute(density))
93
      else
94 8
        sesp$fit.cases$densfun <- deparse(substitute(density.cases))
95
    }
96
  }
97

98
  # if there was an auc and a ci, re-do them
99 8
  if (!is.null(roc$auc) && reuse.auc) {
100 8
    args <- attributes(roc$auc)
101 8
    args$roc <- NULL
102 8
    args$smooth.roc <- sesp
103 8
    sesp$auc <- do.call("auc.smooth.roc", args)
104
  }
105 8
  if (!is.null(roc$ci) && reuse.ci){
106 0
    args <- attributes(roc$ci)
107 0
    args$roc <- NULL
108 0
    args$smooth.roc <- sesp
109 0
    sesp$ci <- do.call(paste(class(roc$ci), "smooth.roc", sep="."), args)
110
  }
111

112 8
  return(sesp)
113
}
114

115
smooth.roc.density <- function(roc, n, density.controls, density.cases, bw,
116
                               # catch args for density
117
                               cut = 3, adjust = 1, kernel = window, window = "gaussian",
118
                               percent = roc$percent, direction = roc$direction,
119
                               ...) {
120 8
  if (!is.numeric(density.controls) || !is.numeric(density.cases)) {
121 8
  	predictor <- c(roc$controls, roc$cases)
122 8
    if (is.character(bw)) {
123 8
    	bw <- match.fun(paste("bw", bw, sep="."))(predictor)
124
    }
125 8
    bw <- bw * adjust
126 8
    from <- min(predictor) - (cut * bw)
127 8
    to <- max(predictor) + (cut * bw)
128
  }
129 8
  if (mode(density.controls) == "function") {
130 8
    density.controls <- density.controls(roc$controls, n=n, from=from, to=to, bw=bw, kernel=kernel, ...)
131 8
    if (! is.numeric(density.controls)) {
132 0
      if (is.list(density.controls) && !is.null(density.controls$y) && is.numeric(density.controls$y))
133 0
        density.controls <- density.controls$y
134
      else
135 0
        stop("The 'density' function must return a numeric vector or a list with a 'y' item.")
136
    }
137
  }
138 8
  else if (is.null(density.controls))
139 8
    density.controls <- suppressWarnings(density(roc$controls, n=n, from=from, to=to, bw=bw, kernel=kernel, ...))$y
140 8
  else if (! is.numeric(density.controls))
141 0
    stop("'density.controls' must be either NULL, a function or numeric values of density (over the y axis).")
142 8
  if (mode(density.cases) == "function") {
143 8
    density.cases <- density.cases(roc$cases, n=n, from=from, to=to, bw=bw, kernel=kernel, ...)
144 8
    if (! is.numeric(density.cases)) {
145 0
      if (is.list(density.cases) && !is.null(density.cases$y) && is.numeric(density.cases$y))
146 0
        density.cases <- density.cases$y
147
      else
148 0
        stop("The 'density' function must return a numeric vector or a list with a 'y' item.")
149
    }
150
  }
151 8
  else if (is.null(density.cases))
152 8
    density.cases <- suppressWarnings(density(roc$cases, n=n, from=from, to=to, bw=bw, kernel=kernel, ...))$y
153 8
  else if (! is.numeric(density.cases))
154 0
    stop("'density.cases' must be either NULL, a function or numeric values of density (over the y axis).")
155 8
  if (length(density.controls) != length(density.cases))
156 0
    stop("Length of 'density.controls' and 'density.cases' differ.")
157

158 8
  perfs <- sapply((1:length(density.controls))+.5, roc.utils.perfs.dens, x=(1:length(density.controls))+.5, dens.controls=density.controls, dens.cases=density.cases, direction=direction)
159

160 8
  return(list(sensitivities = perfs[2,] * ifelse(percent, 100, 1),
161 8
              specificities = perfs[1,] * ifelse(percent, 100, 1)))
162
}
163

164
smooth.roc.binormal <- function(roc, n) {
165 8
  df <- data.frame(sp=qnorm(roc$sp * ifelse(roc$percent, 1/100, 1)), se=qnorm(roc$se * ifelse(roc$percent, 1/100, 1)))
166 8
  df <- df[apply(df, 1, function(x) all(is.finite(x))),]
167 8
  if (dim(df)[1] <= 1) # ROC curve or with only 1 point
168 0
    stop("ROC curve not smoothable (not enough points).")
169 8
  model <- lm(sp~se, df)
170 8
  if(any(is.na(model$coefficients[2])))
171 0
    stop("ROC curve not smoothable (not enough points).")
172 8
  se <- qnorm(seq(0, 1, 1/(n-1)))
173 8
  sp <- predict(model, data.frame(se))
174

175 8
  return(list(sensitivities = pnorm(se) * ifelse(roc$percent, 100, 1),
176 8
              specificities = pnorm(sp) * ifelse(roc$percent, 100, 1),
177 8
              model = model))
178
}
179

180
smooth.roc.fitdistr <- function(roc, n, densfun.controls, densfun.cases, start.controls, start.cases, ...) {
181 8
  load.suggested.package("MASS")
182

183 8
  densfuns.list <- list(beta = "dbeta", cauchy = "dcauchy", "chi-squared" = "dchisq", exponential = "dexp", f = "df",
184 8
                        gamma = "dgamma", geometric = "dgeom", "log-normal" = "dlnorm", lognormal = "dlnorm",
185 8
                        logistic = "dlogis", "negative binomial" = "dnbinom", normal = "dnorm", poisson = "dpois",
186 8
                        t = "dt", weibull = "dweibull")
187

188 8
  if (is.null(densfun.controls))
189 8
    densfun.controls <- "normal"
190 8
  else if (is.character(densfun.controls))
191 8
    densfun.controls <- match.arg(densfun.controls, names(densfuns.list))
192

193 8
  if (is.null(densfun.cases))
194 8
    densfun.cases <- "normal"              
195 8
  else if (is.character(densfun.cases))
196 8
    densfun.cases <- match.arg(densfun.cases, names(densfuns.list))
197

198 8
  fit.controls <- MASS::fitdistr(roc$controls, densfun.controls, start.controls, ...)
199 8
  fit.cases <- MASS::fitdistr(roc$cases, densfun.cases, start.cases, ...)
200

201
  # store function name in fitting results
202 8
  if (mode(densfun.controls) != "function")
203 8
    fit.controls$densfun <- densfun.controls
204 8
  if (mode(densfun.cases) != "function")
205 8
    fit.cases$densfun <- densfun.cases
206

207 8
  x <- seq(min(c(roc$controls, roc$cases)), max(c(roc$controls, roc$cases)), length.out=n)
208

209
  # get the actual function name for do.call
210 8
  if (is.character(densfun.controls))
211 8
    densfun.controls <- match.fun(densfuns.list[[densfun.controls]])
212 8
  if (is.character(densfun.cases))
213 8
    densfun.cases <- match.fun(densfuns.list[[densfun.cases]])
214

215
  # ... that should be passed to densfun
216 8
  if (length(list(...)) > 0 && any(names(list(...)) %in% names(formals(densfun.controls))))
217 0
    dots.controls <- list(...)[names(formals(densfun.controls))[match(names(list(...)), names(formals(densfun.controls)))]]
218
  else
219 8
    dots.controls <- list()
220 8
  if (length(list(...)) > 0 && any(names(list(...)) %in% names(formals(densfun.cases))))
221 0
    dots.cases <- list(...)[names(formals(densfun.cases))[match(names(list(...)), names(formals(densfun.cases)))]]
222
  else
223 8
    dots.cases <- list()
224
   
225 8
  density.controls <- do.call(densfun.controls, c(list(x=x), fit.controls$estimate, dots.controls))
226 8
  density.cases <- do.call(densfun.cases, c(list(x=x), fit.cases$estimate, dots.cases))
227

228 8
  perfs <- sapply(x, roc.utils.perfs.dens, x=x, dens.controls=density.controls, dens.cases=density.cases, direction=roc$direction)
229

230 8
  return(list(sensitivities = perfs[2,] * ifelse(roc$percent, 100, 1),
231 8
              specificities = perfs[1,] * ifelse(roc$percent, 100, 1),
232 8
              fit.controls=fit.controls, fit.cases=fit.cases))
233
}
234

235
smooth.roc.logcondens <- function(roc, n) {
236 8
  load.suggested.package("logcondens")
237

238 8
  sp <- seq(0, 1, 1/(n-1))
239 8
  logcondens <- logcondens::logConROC(roc$cases, roc$controls, sp)
240 8
  se <- logcondens$fROC
241

242 8
  return(list(sensitivities = se * ifelse(roc$percent, 100, 1),
243 8
              specificities = (1 - sp) * ifelse(roc$percent, 100, 1),
244 8
              logcondens = logcondens))
245
}
246

247
smooth.roc.logcondens.smooth <- function(roc, n) {
248 8
  load.suggested.package("logcondens")
249

250 8
  sp <- seq(0, 1, 1/(n-1))
251 8
  logcondens <- logcondens::logConROC(roc$cases, roc$controls, sp)
252 8
  se <- logcondens$fROC.smooth
253

254 8
  return(list(sensitivities = se * ifelse(roc$percent, 100, 1),
255 8
              specificities = (1 - sp) * ifelse(roc$percent, 100, 1),
256 8
              logcondens = logcondens))
257
}

Read our documentation on viewing source code .

Loading