-
Notifications
You must be signed in to change notification settings - Fork 3
/
plot_cv.R
110 lines (104 loc) · 5.15 KB
/
plot_cv.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
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
#' Visualizations for the prediction diagnostics of an estimated ecological niche
#'
#' Create multiple plots of output from the \code{\link{lrren}} function, specifically for the internal k-fold cross-validation diagnostics.
#'
#' @param input An object of class 'list' from the \code{\link{lrren}} function.
#' @param alpha Numeric. The two-tailed alpha level for the significance threshold (default is 0.05).
#'
#' @return This function produces two plots: 1) area under the receiver operating characteristic curve and 2) precision-recall curve. Each plot shows predictions for the log relative risk surface. The red-colored lines are the average curves.
#'
#' @importFrom cvAUC ci.cvAUC cvAUC
#' @importFrom fields image.plot
#' @importFrom graphics abline layout legend lines mtext par plot plot.new title
#' @importFrom ROCR performance prediction
#' @export
#'
#' @examples
#' if (interactive()) {
#' set.seed(1234) # for reproducibility
#'
#' # Using the 'bei' and 'bei.extra' data within {spatstat.data}
#'
#' # Covariate data (centered and scaled)
#' elev <- spatstat.data::bei.extra[[1]]
#' grad <- spatstat.data::bei.extra[[2]]
#' elev$v <- scale(elev)
#' grad$v <- scale(grad)
#' elev_raster <- terra::rast(elev)
#' grad_raster <- terra::rast(grad)
#'
#' # Presence data
#' presence <- spatstat.data::bei
#' spatstat.geom::marks(presence) <- data.frame("presence" = rep(1, presence$n),
#' "lon" = presence$x,
#' "lat" = presence$y)
#' spatstat.geom::marks(presence)$elev <- elev[presence]
#' spatstat.geom::marks(presence)$grad <- grad[presence]
#'
#' # (Pseudo-)Absence data
#' absence <- spatstat.random::rpoispp(0.008, win = elev)
#' spatstat.geom::marks(absence) <- data.frame("presence" = rep(0, absence$n),
#' "lon" = absence$x,
#' "lat" = absence$y)
#' spatstat.geom::marks(absence)$elev <- elev[absence]
#' spatstat.geom::marks(absence)$grad <- grad[absence]
#'
#' # Combine into readable format
#' obs_locs <- spatstat.geom::superimpose(presence, absence, check = FALSE)
#' obs_locs <- spatstat.geom::marks(obs_locs)
#' obs_locs$id <- seq(1, nrow(obs_locs), 1)
#' obs_locs <- obs_locs[ , c(6, 2, 3, 1, 4, 5)]
#'
#' # Run lrren
#' test_lrren <- lrren(obs_locs = obs_locs,
#' cv = TRUE)
#'
#' # Run plot_cv
#' plot_cv(input = test_lrren)
#' }
#'
plot_cv <- function(input, alpha = 0.05) {
if (alpha >= 1 | alpha <= 0) {
stop("The argument 'alpha' must be a numeric value between 0 and 1")
}
op <- graphics::par(no.readonly = TRUE)
on.exit(graphics::par(op))
kfold <- length(input$cv$cv_predictions_rr)
nsamp <- input$out$presence$n
out_cv_rr <- cvAUC::cvAUC(input$cv$cv_predictions_rr, input$cv$cv_labels)
out_ci_rr <- cvAUC::ci.cvAUC(input$cv$cv_predictions_rr, input$cv$cv_labels,
confidence = 1 - alpha)
pred_rr <- ROCR::prediction(input$cv$cv_predictions_rr, input$cv$cv_labels)
perf_rr <- ROCR::performance(pred_rr, "prec", "rec") # PRREC same as "ppv", "tpr"
graphics::layout(matrix(c(1, 2, 3, 3), ncol = 2, byrow = TRUE), heights = c(4, 1))
graphics::par(oma = c(0, 1, 0, 0), mar = c(0.1, 4.1, 4.1, 2.1), pty = "s")
graphics::plot(out_cv_rr$perf, col = "black", lty = 3,
xlab = "False Positive Rate (FPR)\n",
ylab = "\nTrue Positive Rate (TPR)") #Plot fold AUCs
graphics::abline(0, 1, col = "black", lty = 2)
graphics::plot(out_cv_rr$perf, col = "red", avg = "vertical", add = TRUE, lwd = 2) #Plot CV AUC
graphics::title(paste("Area Under the ROC Curve\nAUC = ",
round(out_cv_rr$cvAUC, digits = 3), " (", floor((1-alpha)*100), "% CI: ",
round(out_ci_rr$ci[1], digits = 3), " - ",
round(out_ci_rr$ci[2], digits = 3), ")", sep = ""),
cex.main = 1.1)
graphics::plot(perf_rr, ylim = c(0, 1), xlim = c(0, 1), lty = 3,
xlab = "True Positive Rate (Sensitivity or Recall)\n",
ylab = "\nPositive Predictive Value (Precision)")
graphics::abline((nsamp / kfold) / length(input$cv$cv_labels[[1]]), 0, lty = 2, col = "black")
suppressWarnings(graphics::lines(colMeans(do.call(rbind, perf_rr@x.values)),
colMeans(do.call(rbind, perf_rr@y.values)),
col = "red", lty = 1, lwd = 2)) # mean PRREC
graphics::title("Precision-Recall Curve", cex.main = 1.1)
graphics::par(mai = c(0, 0, 0, 0), mar = c(5.1, 4.1, 0.1, 2.1) / 5, pty = "m")
graphics::plot.new()
graphics::legend(x = "top", inset = 0, title = "Legend",
legend = c("individual k-fold",
"average",
"luck (reference)"),
lty = c(3, 1, 2), bty = "n",
col = c("black", "red", "black"))
graphics::mtext(paste("Internal ", kfold,
"-fold cross-validation, alpha = ", alpha, sep = ""),
side = 3, line = -4, outer = TRUE, cex = 1.25)
}