mcmcRocPrc.R is quite long
1 |
#'This function generates ROC and precision-recall curves
|
|
2 |
#'after fitting a Bayesian logit or probit model.
|
|
3 |
#'@title ROC and Precision-Recall Curves using Bayesian MCMC estimates generalized
|
|
4 |
#'@description This function generates ROC and Precision-Recall curves
|
|
5 |
#'after fitting a Bayesian logit or probit regression. For fast calculation for
|
|
6 |
#'from an "rjags" object use \code{\link{mcmcRocPrc}}
|
|
7 |
#'@param modelmatrix model matrix, including intercept (if the intercept is among the
|
|
8 |
#'parameters estimated in the model). Create with model.matrix(formula, data).
|
|
9 |
#'Note: the order of columns in the model matrix must correspond to the order of columns
|
|
10 |
#'in the matrix of posterior draws in the \code{mcmcout} argument. See the \code{mcmcout}
|
|
11 |
#'argument for more and Beger (2016) for background.
|
|
12 |
#'@param mcmcout posterior distributions of all logit coefficients,
|
|
13 |
#'in matrix form. This can be created from rstan, MCMCpack, R2jags, etc. and transformed
|
|
14 |
#'into a matrix using the function as.mcmc() from the coda package for \code{jags} class
|
|
15 |
#'objects, as.matrix() from base R for \code{mcmc}, \code{mcmc.list}, \code{stanreg}, and
|
|
16 |
#'\code{stanfit} class objects, and \code{object$sims.matrix} for \code{bugs} class objects.
|
|
17 |
#'Note: the order of columns in this matrix must correspond to the order of columns
|
|
18 |
#'in the model matrix. One can do this by examining the posterior distribution matrix and sorting the
|
|
19 |
#'variables in the order of this matrix when creating the model matrix. A useful function for sorting
|
|
20 |
#'column names containing both characters and numbers as
|
|
21 |
#'you create the matrix of posterior distributions is \code{mixedsort()} from the gtools package.
|
|
22 |
#'@param modelframe model frame in matrix form. Can be created using
|
|
23 |
#'as.matrix(model.frame(formula, data))
|
|
24 |
#'@param curves logical indicator of whether or not to return values to plot the ROC or Precision-Recall
|
|
25 |
#'curves. If set to \code{FALSE} (default), results are returned as a list without the extra
|
|
26 |
#'values.
|
|
27 |
#'@param link type of generalized linear model; a character vector set to \code{"logit"} (default) or \code{"probit"}.
|
|
28 |
#'@param fullsims logical indicator of whether full object (based on all MCMC draws
|
|
29 |
#'rather than their average) will be returned. Default is \code{FALSE}. Note: If \code{TRUE}
|
|
30 |
#'is chosen, the function takes notably longer to execute.
|
|
31 |
#'@references Beger, Andreas. 2016. “Precision-Recall Curves.” Available at SSRN:
|
|
32 |
#'https://ssrn.com/Abstract=2765419. http://dx.doi.org/10.2139/ssrn.2765419.
|
|
33 |
#'@return This function returns a list with 4 elements:
|
|
34 |
#'\itemize{
|
|
35 |
#'\item area_under_roc: area under ROC curve (scalar)
|
|
36 |
#'\item area_under_prc: area under precision-recall curve (scalar)
|
|
37 |
#'\item prc_dat: data to plot precision-recall curve (data frame)
|
|
38 |
#'\item roc_dat: data to plot ROC curve (data frame)
|
|
39 |
#'}
|
|
40 |
#'
|
|
41 |
#'@examples
|
|
42 |
#' \dontshow{.old_wd <- setwd(tempdir())}
|
|
43 |
#' \donttest{
|
|
44 |
#' # simulating data
|
|
45 |
#'
|
|
46 |
#' set.seed(123456)
|
|
47 |
#' b0 <- 0.2 # true value for the intercept
|
|
48 |
#' b1 <- 0.5 # true value for first beta
|
|
49 |
#' b2 <- 0.7 # true value for second beta
|
|
50 |
#' n <- 500 # sample size
|
|
51 |
#' X1 <- runif(n, -1, 1)
|
|
52 |
#' X2 <- runif(n, -1, 1)
|
|
53 |
#' Z <- b0 + b1 * X1 + b2 * X2
|
|
54 |
#' pr <- 1 / (1 + exp(-Z)) # inv logit function
|
|
55 |
#' Y <- rbinom(n, 1, pr)
|
|
56 |
#' df <- data.frame(cbind(X1, X2, Y))
|
|
57 |
#'
|
|
58 |
#' # formatting the data for jags
|
|
59 |
#' datjags <- as.list(df)
|
|
60 |
#' datjags$N <- length(datjags$Y)
|
|
61 |
#'
|
|
62 |
#' # creating jags model
|
|
63 |
#' model <- function() {
|
|
64 |
#'
|
|
65 |
#' for(i in 1:N){
|
|
66 |
#' Y[i] ~ dbern(p[i]) ## Bernoulli distribution of y_i
|
|
67 |
#' logit(p[i]) <- mu[i] ## Logit link function
|
|
68 |
#' mu[i] <- b[1] +
|
|
69 |
#' b[2] * X1[i] +
|
|
70 |
#' b[3] * X2[i]
|
|
71 |
#' }
|
|
72 |
#'
|
|
73 |
#' for(j in 1:3){
|
|
74 |
#' b[j] ~ dnorm(0, 0.001) ## Use a coefficient vector for simplicity
|
|
75 |
#' }
|
|
76 |
#'
|
|
77 |
#' }
|
|
78 |
#'
|
|
79 |
#' params <- c("b")
|
|
80 |
#' inits1 <- list("b" = rep(0, 3))
|
|
81 |
#' inits2 <- list("b" = rep(0, 3))
|
|
82 |
#' inits <- list(inits1, inits2)
|
|
83 |
#'
|
|
84 |
#' ## fitting the model with R2jags
|
|
85 |
#' set.seed(123)
|
|
86 |
#' fit <- R2jags::jags(data = datjags, inits = inits,
|
|
87 |
#' parameters.to.save = params, n.chains = 2, n.iter = 2000,
|
|
88 |
#' n.burnin = 1000, model.file = model)
|
|
89 |
#'
|
|
90 |
#' # processing the data
|
|
91 |
#' mm <- model.matrix(Y ~ X1 + X2, data = df)
|
|
92 |
#' xframe <- as.matrix(model.frame(Y ~ X1 + X2, data = df))
|
|
93 |
#' mcmc <- coda::as.mcmc(fit)
|
|
94 |
#' mcmc_mat <- as.matrix(mcmc)[, 1:ncol(xframe)]
|
|
95 |
#'
|
|
96 |
#' # using mcmcRocPrcGen
|
|
97 |
#' fit_sum <- mcmcRocPrcGen(modelmatrix = mm,
|
|
98 |
#' modelframe = xframe,
|
|
99 |
#' mcmcout = mcmc_mat,
|
|
100 |
#' curves = TRUE,
|
|
101 |
#' fullsims = FALSE)
|
|
102 |
#' }
|
|
103 |
#'
|
|
104 |
#' \dontshow{setwd(.old_wd)}
|
|
105 |
#'@export
|
|
106 |
|
|
107 |
mcmcRocPrcGen <- function(modelmatrix, |
|
108 |
mcmcout,
|
|
109 |
modelframe,
|
|
110 |
curves = FALSE, |
|
111 |
link = "logit", |
|
112 |
fullsims = FALSE){ |
|
113 |
|
|
114 |
if(link == "logit") { |
|
115 |
pred_prob <- plogis(t(modelmatrix %*% t(mcmcout))) |
|
116 |
} else if (link == "probit") { |
|
117 |
pred_prob <- pnorm(t(modelmatrix %*% t(mcmcout))) |
|
118 |
} else { |
|
119 |
stop("Please enter a valid link argument") |
|
120 |
}
|
|
121 |
|
|
122 |
if(missing(modelmatrix) | missing(mcmcout) | missing(modelframe)) { |
|
123 |
"Please enter the required arguments"
|
|
124 |
}
|
|
125 |
|
|
126 |
if(fullsims == FALSE){ |
|
127 |
y_pred <- apply(X = pred_prob, MARGIN = 2, FUN = function(x) median(x)) |
|
128 |
|
|
129 |
# Observed y and x
|
|
130 |
pred_obs <- data.frame(y_pred = y_pred, y_obs = modelframe[, 1]) |
|
131 |
|
|
132 |
auc_roc <- function(obs, pred) { |
|
133 |
pred <- prediction(pred, obs) |
|
134 |
auc <- performance(pred, "auc")@y.values[[1]] |
|
135 |
return(auc) |
|
136 |
}
|
|
137 |
|
|
138 |
auc_pr <- function(obs, pred) { |
|
139 |
xx.df <- prediction(pred, obs) |
|
140 |
perf <- performance(xx.df, "prec", "rec") |
|
141 |
xy <- data.frame(recall = perf@x.values[[1]], |
|
142 |
precision = perf@y.values[[1]]) |
|
143 |
|
|
144 |
# take out division by 0 for lowest threshold
|
|
145 |
xy <- subset(xy, !is.nan(xy$precision)) |
|
146 |
|
|
147 |
res <- caTools::trapz(xy$recall, xy$precision) |
|
148 |
res |
|
149 |
}
|
|
150 |
|
|
151 |
area_under_roc <- auc_roc(obs = pred_obs$y_obs, pred = pred_obs$y_pred) |
|
152 |
|
|
153 |
area_under_prc <- auc_pr(obs = pred_obs$y_obs, pred = pred_obs$y_pred) |
|
154 |
|
|
155 |
if(curves == FALSE){ |
|
156 |
# Results as a list
|
|
157 |
results <- list() |
|
158 |
results$area_under_roc <- area_under_roc |
|
159 |
results$area_under_prc <- area_under_prc |
|
160 |
|
|
161 |
return(results) |
|
162 |
}
|
|
163 |
|
|
164 |
if(curves == TRUE){ |
|
165 |
|
|
166 |
prediction_obj <- prediction(predictions = pred_obs$y_pred, |
|
167 |
labels = pred_obs$y_obs) |
|
168 |
|
|
169 |
prc_performance_obj <- performance(prediction.obj = prediction_obj, |
|
170 |
measure = "prec", |
|
171 |
x.measure = "rec") |
|
172 |
|
|
173 |
prc_dat <- data.frame(x = prc_performance_obj@x.values, |
|
174 |
y = prc_performance_obj@y.values) |
|
175 |
names(prc_dat) <- c("x", "y") |
|
176 |
|
|
177 |
roc_performance_obj <- performance(prediction.obj = prediction_obj, |
|
178 |
measure = "tpr", |
|
179 |
x.measure = "fpr") |
|
180 |
|
|
181 |
roc_dat <- data.frame(x = roc_performance_obj@x.values, |
|
182 |
y = roc_performance_obj@y.values) |
|
183 |
names(roc_dat) <- c("x", "y") |
|
184 |
|
|
185 |
# Results as a list
|
|
186 |
results <- list() |
|
187 |
results$area_under_roc <- area_under_roc |
|
188 |
results$area_under_prc <- area_under_prc |
|
189 |
results$prc_dat <- prc_dat |
|
190 |
results$roc_dat <- roc_dat |
|
191 |
|
|
192 |
return(results) |
|
193 |
}
|
|
194 |
}
|
|
195 |
|
|
196 |
if(fullsims == TRUE){ |
|
197 |
RocPrcOneDraw <- function(pred_prob_vector){ |
|
198 |
# run this function over each row (iteration) of the pred_prob matrix
|
|
199 |
|
|
200 |
# y_pred <- apply(X = pred_prob, MARGIN = 2, FUN = function(x) median(x))
|
|
201 |
|
|
202 |
# Observed y and x
|
|
203 |
pred_obs <- data.frame(y_pred = pred_prob_vector, y_obs = modelframe[, 1]) |
|
204 |
|
|
205 |
auc_roc <- function(obs, pred) { |
|
206 |
pred <- prediction(pred, obs) |
|
207 |
auc <- performance(pred, "auc")@y.values[[1]] |
|
208 |
return(auc) |
|
209 |
}
|
|
210 |
|
|
211 |
auc_pr <- function(obs, pred) { |
|
212 |
xx.df <- prediction(pred, obs) |
|
213 |
perf <- performance(xx.df, "prec", "rec") |
|
214 |
xy <- data.frame(recall = perf@x.values[[1]], |
|
215 |
precision = perf@y.values[[1]]) |
|
216 |
|
|
217 |
# take out division by 0 for lowest threshold
|
|
218 |
xy <- subset(xy, !is.nan(xy$precision)) |
|
219 |
|
|
220 |
res <- caTools::trapz(xy$recall, xy$precision) |
|
221 |
res |
|
222 |
}
|
|
223 |
|
|
224 |
area_under_roc <- auc_roc(obs = pred_obs$y_obs, pred = pred_obs$y_pred) |
|
225 |
|
|
226 |
area_under_prc <- auc_pr(obs = pred_obs$y_obs, pred = pred_obs$y_pred) |
|
227 |
|
|
228 |
if(curves == FALSE){ |
|
229 |
# Results as a list
|
|
230 |
one_result <- c(area_under_roc, area_under_prc) |
|
231 |
return(one_result) |
|
232 |
}
|
|
233 |
|
|
234 |
if(curves == TRUE){ |
|
235 |
prediction_obj <- prediction(predictions = pred_obs$y_pred, |
|
236 |
labels = pred_obs$y_obs) |
|
237 |
|
|
238 |
prc_performance_obj <- performance(prediction.obj = prediction_obj, |
|
239 |
measure = "prec", |
|
240 |
x.measure = "rec") |
|
241 |
|
|
242 |
prc_dat <- data.frame(x = prc_performance_obj@x.values, |
|
243 |
y = prc_performance_obj@y.values) |
|
244 |
names(prc_dat) <- c("x", "y") |
|
245 |
|
|
246 |
roc_performance_obj <- performance(prediction.obj = prediction_obj, |
|
247 |
measure = "tpr", |
|
248 |
x.measure = "fpr") |
|
249 |
|
|
250 |
roc_dat <- data.frame(x = roc_performance_obj@x.values, |
|
251 |
y = roc_performance_obj@y.values) |
|
252 |
names(roc_dat) <- c("x", "y") |
|
253 |
|
|
254 |
# Results as a list
|
|
255 |
one_result <- list() |
|
256 |
one_result$area_under_roc <- area_under_roc |
|
257 |
one_result$area_under_prc <- area_under_prc |
|
258 |
one_result$prc_dat <- prc_dat |
|
259 |
one_result$roc_dat <- roc_dat |
|
260 |
|
|
261 |
return(one_result) |
|
262 |
}
|
|
263 |
}
|
|
264 |
|
|
265 |
if(curves == FALSE){ |
|
266 |
all_results <- matrix(nrow = nrow(pred_prob), ncol = 2) |
|
267 |
for(i in 1:nrow(pred_prob)){ |
|
268 |
all_results[i, ] <- RocPrcOneDraw(pred_prob[i, ]) |
|
269 |
}
|
|
270 |
all_results <- as.data.frame(all_results) |
|
271 |
names(all_results) <- c("area_under_roc", "area_under_prc") |
|
272 |
}
|
|
273 |
|
|
274 |
if(curves == TRUE){ |
|
275 |
all_results <- list() |
|
276 |
for(i in 1:nrow(pred_prob)){ |
|
277 |
all_results[[i]] <- RocPrcOneDraw(pred_prob[i, ]) |
|
278 |
}
|
|
279 |
}
|
|
280 |
|
|
281 |
return(all_results) |
|
282 |
}
|
|
283 |
|
|
284 |
}
|
Read our documentation on viewing source code .