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 0
  if(link == "logit") {
115 0
    pred_prob <- plogis(t(modelmatrix %*% t(mcmcout)))    
116 0
  } else if (link == "probit") {
117 0
    pred_prob <- pnorm(t(modelmatrix %*% t(mcmcout)))
118
  } else {
119 0
    stop("Please enter a valid link argument")
120
  }
121
  
122 0
  if(missing(modelmatrix) | missing(mcmcout) | missing(modelframe)) {
123 0
    "Please enter the required arguments"
124
  }
125
  
126 0
  if(fullsims == FALSE){
127 0
    y_pred <- apply(X = pred_prob, MARGIN = 2, FUN = function(x) median(x))
128
    
129
    # Observed y and x
130 0
    pred_obs <- data.frame(y_pred = y_pred, y_obs = modelframe[, 1])             
131
    
132 0
    auc_roc <- function(obs, pred) {
133 0
      pred <- prediction(pred, obs)
134 0
      auc  <- performance(pred, "auc")@y.values[[1]]
135 0
      return(auc)
136
    }
137
    
138 0
    auc_pr <- function(obs, pred) {
139 0
      xx.df <- prediction(pred, obs)
140 0
      perf  <- performance(xx.df, "prec", "rec")
141 0
      xy    <- data.frame(recall = perf@x.values[[1]], 
142 0
                          precision = perf@y.values[[1]])
143
      
144
      # take out division by 0 for lowest threshold
145 0
      xy <- subset(xy, !is.nan(xy$precision))
146
      
147 0
      res   <- caTools::trapz(xy$recall, xy$precision)
148 0
      res
149
    }
150
    
151 0
    area_under_roc <- auc_roc(obs = pred_obs$y_obs, pred = pred_obs$y_pred)
152
    
153 0
    area_under_prc <- auc_pr(obs = pred_obs$y_obs, pred = pred_obs$y_pred)
154
    
155 0
    if(curves == FALSE){
156
      # Results as a list
157 0
      results <- list()
158 0
      results$area_under_roc <- area_under_roc
159 0
      results$area_under_prc <- area_under_prc
160
      
161 0
      return(results)
162
    }
163
    
164 0
    if(curves == TRUE){
165
      
166 0
      prediction_obj <- prediction(predictions = pred_obs$y_pred,
167 0
                                         labels = pred_obs$y_obs)
168
      
169 0
      prc_performance_obj <- performance(prediction.obj = prediction_obj,
170 0
                                               measure = "prec",
171 0
                                               x.measure = "rec")
172
      
173 0
      prc_dat <- data.frame(x = prc_performance_obj@x.values,
174 0
                            y = prc_performance_obj@y.values)
175 0
      names(prc_dat) <- c("x", "y")
176
      
177 0
      roc_performance_obj <- performance(prediction.obj = prediction_obj,
178 0
                                               measure = "tpr",
179 0
                                               x.measure = "fpr")
180
      
181 0
      roc_dat <- data.frame(x = roc_performance_obj@x.values,
182 0
                            y = roc_performance_obj@y.values)
183 0
      names(roc_dat) <- c("x", "y")
184
      
185
      # Results as a list
186 0
      results <- list()
187 0
      results$area_under_roc <- area_under_roc
188 0
      results$area_under_prc <- area_under_prc
189 0
      results$prc_dat <- prc_dat
190 0
      results$roc_dat <- roc_dat
191
      
192 0
      return(results)
193
    }
194
  }
195
  
196 0
  if(fullsims == TRUE){
197 0
    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 0
      pred_obs <- data.frame(y_pred = pred_prob_vector, y_obs = modelframe[, 1])						 
204
      
205 0
      auc_roc <- function(obs, pred) {
206 0
        pred <- prediction(pred, obs)
207 0
        auc  <- performance(pred, "auc")@y.values[[1]]
208 0
        return(auc)
209
      }
210
      
211 0
      auc_pr <- function(obs, pred) {
212 0
        xx.df <- prediction(pred, obs)
213 0
        perf  <- performance(xx.df, "prec", "rec")
214 0
        xy    <- data.frame(recall = perf@x.values[[1]], 
215 0
                            precision = perf@y.values[[1]])
216
        
217
        # take out division by 0 for lowest threshold
218 0
        xy <- subset(xy, !is.nan(xy$precision))
219
        
220 0
        res   <- caTools::trapz(xy$recall, xy$precision)
221 0
        res
222
      }
223
      
224 0
      area_under_roc <- auc_roc(obs = pred_obs$y_obs, pred = pred_obs$y_pred)
225
      
226 0
      area_under_prc <- auc_pr(obs = pred_obs$y_obs, pred = pred_obs$y_pred)
227
      
228 0
      if(curves == FALSE){
229
        # Results as a list
230 0
        one_result <- c(area_under_roc, area_under_prc)
231 0
        return(one_result)
232
      }
233
      
234 0
      if(curves == TRUE){
235 0
        prediction_obj <- prediction(predictions = pred_obs$y_pred,
236 0
                                           labels = pred_obs$y_obs)
237
        
238 0
        prc_performance_obj <- performance(prediction.obj = prediction_obj,
239 0
                                                 measure = "prec",
240 0
                                                 x.measure = "rec")
241
        
242 0
        prc_dat <- data.frame(x = prc_performance_obj@x.values,
243 0
                              y = prc_performance_obj@y.values)
244 0
        names(prc_dat) <- c("x", "y")
245
        
246 0
        roc_performance_obj <- performance(prediction.obj = prediction_obj,
247 0
                                                 measure = "tpr",
248 0
                                                 x.measure = "fpr")
249
        
250 0
        roc_dat <- data.frame(x = roc_performance_obj@x.values,
251 0
                              y = roc_performance_obj@y.values)
252 0
        names(roc_dat) <- c("x", "y")
253
        
254
        # Results as a list
255 0
        one_result <- list()
256 0
        one_result$area_under_roc <- area_under_roc
257 0
        one_result$area_under_prc <- area_under_prc
258 0
        one_result$prc_dat <- prc_dat
259 0
        one_result$roc_dat <- roc_dat
260
        
261 0
        return(one_result)
262
      }
263
    }
264
    
265 0
    if(curves == FALSE){
266 0
      all_results <- matrix(nrow = nrow(pred_prob), ncol = 2)
267 0
      for(i in 1:nrow(pred_prob)){
268 0
        all_results[i, ] <- RocPrcOneDraw(pred_prob[i, ])
269
      }
270 0
      all_results <- as.data.frame(all_results)
271 0
      names(all_results) <- c("area_under_roc", "area_under_prc")
272
    }
273
    
274 0
    if(curves == TRUE){
275 0
      all_results <- list()
276 0
      for(i in 1:nrow(pred_prob)){
277 0
        all_results[[i]] <- RocPrcOneDraw(pred_prob[i, ])
278
      }
279
    }
280
    
281 0
    return(all_results)
282
  }
283
  
284
}

Read our documentation on viewing source code .

Loading