-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy patharoc_compute.R
69 lines (69 loc) · 2.22 KB
/
aroc_compute.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
#' compute.threshold.YI.AROC.sp function modified to retrieve sensitivity and specificity
#'
#' @param data https://rdrr.io/cran/ROCnReg/src/R/compute.threshold.YI.AROC.sp.R
#'
#' @return compute.threshold.YI.AROC.sp with sensitivity and specificity
#' @export
#'
#' @examples
aroc_compute<-function (object, newdata, parallel = c("no", "multicore", "snow"),
ncpus = 1, cl = NULL)
{
if (class(object)[1] != "AROC.sp") {
stop(paste0("This function can not be used for this object class: ",
class(object)[1]))
}
names.cov <- all.vars(object$formula)[-1]
if (!missing(newdata) && !inherits(newdata, "data.frame"))
stop("Newdata must be a data frame")
if (!missing(newdata) && length(names.cov) != 0 && sum(is.na(match(names.cov,
names(newdata)))))
stop("Not all needed variables are supplied in newdata")
if (missing(newdata)) {
newdata <- cROCData(object$data, names.cov, object$group)
}
else {
newdata <- na.omit(newdata[, names.cov, drop = FALSE])
}
p <- seq(0, 1, length = 500)
np <- length(p)
npred <- nrow(newdata)
sigma0 <- summary(object$fit)$sigma
data.d <- (object$data[object$data[, object$group] != object$tag.h,
])[!object$missing.ind$d, ]
n1 <- nrow(data.d)
pre.placement.values <- (data.d[, object$marker] - predict(object$fit,
newdata = data.d))/sigma0
if (object$est.cdf.h == "normal") {
u1 <- 1 - pnorm(pre.placement.values)
}
else {
res0p <- object$fit$residuals/sigma0
F0res <- ecdf(res0p)
u1 <- 1 - F0res(pre.placement.values)
}
AROC <- numeric(np)
for (i in 1:np) {
AROC[i] <- sum(u1 <= p[i])/n1
}
difbb <- AROC - p
FPF <- mean(p[which(difbb == max(difbb))])
YI <- max(difbb)
pred0 <- predict(object$fit, newdata = newdata)
if (object$est.cdf == "normal") {
csf0_inv <- qnorm(1 - FPF)
}
else {
csf0_inv <- quantile(res0p, 1 - FPF, type = 1)
}
thresholds <- pred0 + sigma0 * csf0_inv
res <- list()
res$call <- match.call()
res$newdata <- newdata
res$thresholds <- thresholds
res$YI <- YI
res$FPF <- FPF
res$specificity <- (1-p)
res$sensitivity <- (AROC)
res
}