Calculate Cohen's Kappa from a vector of binary observations and a vector of continuous or binary predictions

kappa.R

confusion <- function(obs, pred, thr=NULL) {
  # obs: a vector of observed binary values (0, 1)
  # pred: a vector of continuous predictions (range: 0-1)
  # thr: a numeric scalar specifying the value at which to threshold 
  #      `pred`, or NULL if pred is already binary.
  if(is.null(thr) && !all(pred %in% 0:1)) {
    stop('If thr=NULL, all pred must be either 0 or 1.')
  }
  if(any(pred > 1) || any(pred < 0)) {
    stop('pred must all be between 0 and 1 (inclusive).')
  }
  if(!is.null(thr)) {
    pred <- as.numeric(pred >= thr)
  }
  table(obs=factor(obs, levels=0:1), pred=factor(pred, levels=0:1))
}

kappa <- function(obs, pred, thr=NULL) {
  # obs: a vector of observed binary values (0, 1)
  # pred: a vector of continuous predictions (range: 0-1)
  # thr: a numeric scalar specifying the value at which to threshold 
  #      `pred`, or NULL if pred is already binary.
  m <- confusion(obs, pred, thr)
  n <- sum(m)
  po <- sum(diag(m))/n
  p0 <- sum(m[2, ])/n * sum(m[, 2])/n
  p1 <- sum(m[1, ])/n * sum(m[, 1])/n
  pe <- p0 + p1
  (po - pe)/(1 - pe)
}

max_kappa <- function(obs, pred, precision=0.001, plot=FALSE) {
  # obs: a vector of observed binary values (0, 1)
  # pred: a vector of continuous predictions (range: 0-1)
  # precision: controls the number thresholds to evaluate. Precision of 0.001
  #            corresponds to 1/0.001 = 1000 thresholds evaluated b/w 0 and 1.
  # plot: logical. Plot a trace of Kappa wrt threshold?
  if(all(pred %in% 0:1)) stop('pred are all 0/1.')
  thrs <- seq(0, 1, length.out=1/precision)
  kk <- sapply(thrs, function(x) kappa(obs, pred, x))
  k_max <- max(kk)
  thr_max <- thrs[which.max(kk)]
  if(isTRUE(plot)) {
    plot(kk~thrs, xlab='Threshold', ylab='Kappa', type='l', las=1)
    abline(v=thr_max, h=k_max, lwd=2, col='#00000050')
    axis(3, at=thr_max, cex.axis=0.8, line=-1, xpd=NA, tick=FALSE, 
         labels=bquote(
           max~kappa*':'~.(round(k_max, 3))~'('*
             threshold==.(round(thr_max, 3))*')'))
  }
  list(kappa=k_max, threshold=thr_max)
}

Example:

set.seed(1)
pred <- runif(1000)
obs <- rbinom(1000, 1, pred)
max_kappa(obs, pred, precision=0.001, plot=TRUE)

$kappa
[1] 0.4921514

$threshold
[1] 0.5445445