R/ecospat.eval.R

In ecospat: Spatial Ecology Miscellaneous Methods

################ MEVA.TABLE: Function originally from A. Guisan (Unil-ECOSPAT, Switzerland)
################ Modified from L. Maiorano, O. Broennimann and F. Collart (Unil-ECOSPAT, Switzerland)
ecospat.meva.table <- function(Pred, Sp.occ, th) # Pred: vector of predicted probabilities Sp.occ: vector of binary observations
# th: threshold used to cut the probability to binary values
{
  # eva <- list()
  tab <- table(Pred >= th, Sp.occ, dnn = list("Predicted values", "Observed values"))
  a <- tab[4]
  b <- tab[2]
  c <- tab[3]
  d <- tab[1]
  N <- a + b + c + d
  prev <- (a + c)/N
  ccr <- (a + d)/N
  mcr <- (b + c)/N
  se <- a/(a + c)
  sp <- d/(b + d)
  fpr <- b/(b + d)
  fnr <- c/(a + c)
  ppp <- a/(a + b)
  npp <- d/(c + d)
  or <- (a * d)/(c * b)
  kappa <- ecospat.cohen.kappa(tab)$kap
  nmi <- (-(a * log(a)) - (b * log(b)) - (c * log(c)) - (d * log(d)) + ((a + b) *
    log(a + b)) + ((c + d) * log(c + d)))/((N * log(N)) - ((a + c) * log(a +
    c)) - ((b + d) * log(b + d)))
  tss <- se + sp - 1
  orss <- (a*d - b*c)/(a*d + b*c)
  sed <- ( log(fpr) - log(se) -  log(1-fpr) + log(1-se) ) /( log(fpr) + log(se) +  log(1-fpr) + log(1-se))
  
  mtrcNAME <- c("Prevalence", "Correct classification rate", "Mis-classification rate",
    "Sensitivity", "Specificity", "False positive rate", "False negative rate",
    "Positive predictive power", "Negative predictive power", "Odds Ratio", "Kappa", "Normalized mutual information",
    "True skill statistic","Odds Ratio Skill Score","Symmetric Extremal Dependence Index")
  mtrcVALUE <- c(round(prev, digits = 4), round(ccr, digits = 4), round(mcr, digits = 4),
    round(se, digits = 4), round(sp, digits = 4), round(fpr, digits = 4), round(fnr,
      digits = 4), round(ppp, digits = 4), round(npp, digits = 4), round(or,
      digits = 4), round(kappa, digits = 4), round(nmi, digits = 4), round(tss,
      digits = 4),round(orss, digits = 4),round(sed, digits = 4))
  mtrcMATR <- cbind.data.frame(Metric=mtrcNAME, Value =mtrcVALUE)
  return(list(CONTINGENCY_TABLE = tab, EVALUATION_METRICS = mtrcMATR))
}


############### MAX-KAPPA ## Function originally from A. Guisan (Unil-ECOSPAT, Switzerland)
############### Modified from L. Maiorano (Unil-ECOSPAT, Switzerland) and Olivier Broennimann
ecospat.max.kappa <- function(Pred, Sp.occ) # Pred: vector of predicted probabilities Sp.occ: vector of binary observations
{
  FUN<-function(thresh){
    xtab<-table(Pred >= thresh, Sp.occ)
    return(ecospat.cohen.kappa(xtab)$kap )
  }
  
  threshold<-(1:100)/100
  k<-sapply(threshold,FUN)
  table<-data.frame(cbind(threshold,k))
  
  max.Kappa <- max(table$k,na.rm = TRUE)
  max.threshold <- table$threshold[which.max(table$k)][1]
  return(list(table=table,max.Kappa=max.Kappa,max.threshold=max.threshold))
}


############### MAX-TSS ## Function originally from L. Maiorano (Unil-ECOSPAT, Switzerland)
############### modyfying MAX-KAPPA from A. Guisan (Unil-ECOSPAT, Switzerland), modified by Olivier Broennimann and Flavien Collart
ecospat.max.tss <- function(Pred, Sp.occ) # Pred: vector of predicted probabilities Sp.occ: vector of binary observations
{
  FUN<-function(thresh){
    xtab<-table(Pred >= thresh, Sp.occ)
    if(nrow(xtab)!=2){ ## if the model is bad the Pred values are not changing across the whole dataset
      return(0)
    }else{
      a <- xtab[4]
      b <- xtab[2]
      c <- xtab[3]
      d <- xtab[1]
      se <- a/(a + c)
      sp <- d/(b + d)
      return(se + sp - 1)
    }
  }
  threshold<-(1:100)/100
  tss<-sapply(threshold,FUN)
  table<-data.frame(cbind(threshold,tss))
  
  max.TSS <- max(table$tss,na.rm = TRUE)
  max.threshold <- table$threshold[which.max(table$tss)][1]
  return(list(table=table,max.TSS=max.TSS,max.threshold=max.threshold))
}

################## PLOT-K ## Function originally from L. Maiorano (Unil-ECOSPAT, Switzerland)
ecospat.plot.kappa <- function(Pred, Sp.occ) # Pred: vector of predicted probabilities Sp.occ: vector of binary observations
{
  k <- 0.01
  i <- 1
  evak <- data.frame(k)
  while (k < 1) {
    a <- table(Pred >= k, Sp.occ)[4]
    b <- table(Pred >= k, Sp.occ)[2]
    c <- table(Pred >= k, Sp.occ)[3]
    d <- table(Pred >= k, Sp.occ)[1]
    N <- a + b + c + d
    tab <- table(Pred >= k, Sp.occ)
    evak[i, "THRESHOLD"] <- k
    evak[i, "k"] <- ecospat.cohen.kappa(tab)$kap
    # evak[i,'k'] <- ecospat.cohen.kappa(table(Pred >= k, Sp.occ))$kap
    k <- k + 0.01
    i <- i + 1
  }
  plot(evak$THRESHOLD, evak$k, type = "l", xlab = "THRESHOLD", ylab = "KAPPA")
}


################## PLOT-TSS ## Function originally from L. Maiorano (Unil-ECOSPAT, Switzerland)
ecospat.plot.tss <- function(Pred, Sp.occ) # Pred: vector of predicted probabilities Sp.occ: vector of binary observations
{
  tss <- 0.01
  i <- 1
  evatss <- data.frame(tss)
  while (tss < 1) {
    a <- table(Pred >= tss, Sp.occ)[4]
    b <- table(Pred >= tss, Sp.occ)[2]
    c <- table(Pred >= tss, Sp.occ)[3]
    d <- table(Pred >= tss, Sp.occ)[1]
    N <- a + b + c + d
    se <- a/(a + c)
    sp <- d/(b + d)
    evatss[i, "THRESHOLD"] <- tss
    evatss[i, "tss"] <- se + sp - 1
    tss <- tss + 0.01
    i <- i + 1
  }
  plot(evatss$THRESHOLD, evatss$tss, type = "l", xlab = "THRESHOLD", ylab = "TRUE SKILL STATISTIC")
}

############################################################################################################################


ecospat.cohen.kappa <- function(xtab) {
  # computes Cohen's kappa and variance estimates, 95% CI from a symmetric
  # agreement table xtab return value is a list with elements kap and vark See
  # Bishop et al. Disc. Mult. Analysis pp. 395-397 1975

  # Originally coded by C. Randin, UNIL Adjusted by N.E.Zimmermann, WSL

  ci <- vector(length = 2)
  kap <- 0
  vark <- 0
  totn <- 0
  if (nrow(xtab) != ncol(xtab)) {
    # cat('\n freq. table not symmetric.\n')
    result = list(kap = kap, vark = vark, totn = totn)
    return(result)
  }
  # compute obs props:
  totn <- sum(xtab)
  # calc marginals, expected diag props.
  obs <- xtab/totn
  pi <- apply(xtab, 1, sum)/totn
  pj <- apply(xtab, 2, sum)/totn
  exp <- outer(pi, pj)
  theta.2 <- sum(diag(exp))
  theta.1 <- sum(diag(obs))
  kap <- (theta.1 - theta.2)/(1 - theta.2)
  theta.3 <- diag(obs) %*% (pi + pj)
  theta.4 <- sum(obs * outer(pi, pj, FUN = "+")^2)
  # cat('\n thetas 1,2,3,4 \n', theta.1, theta.2, theta.3, theta.4, '\n') calc
  # the var est
  vark <- (theta.1 * (1 - theta.1))/(1 - theta.2)^2
  vark <- vark + (2 * (1 - theta.1) * (2 * theta.1 * theta.2 - theta.3))/(1 - theta.2)^3
  vark <- vark + ((1 - theta.1)^2 * (theta.4 - 4 * theta.2^2))/(1 - theta.2)^4
  vark <- as.vector(vark/totn)
  ci[1] <- kap - 1.96 * sqrt(vark)
  ci[2] <- kap + 1.96 * sqrt(vark)
  result = list(kap = kap, vark = vark, totn = totn, ci = ci)
  return(result)
}