evaluationFunctions.R

#preModelingFunctions.R
#Set of functions to evaluate species distribution models
#Author: Jorge Velásquez

#EvaluatePOModel
#Evaluation of maxen models. This uses k-fold partitioning to divide up occurrences 
#into a train and test sets, and fixed background/pseudoabsence train and test set. 
#Arguments:
# folds(numeric): folds to use in k-fold partitioning
# covs.pres(data frame): covariates associated with occurrences
# covs.bkg.train (data frame): covariates associated with background used for training
# covs.bkg.train (data frame): covariates associated with background used for testing
# mxnt.args(string): string to be passed to dismo's function maxent though the args argument
# path(string): string to be passed to dismo's function maxent though the path argument
#Returns:
# A data frame with values n.train(# training points), n.test(# testing points), 
# nparams (# parameters, based on the number of features with coefficients !=0),
# auc.train (AUC of training subset), auc.test (AUC of testing subset), for each
# fold.

EvaluatePOModel <- function(folds, covs.pres, covs.bkg.train, covs.bkg.test, mxnt.args, path){
  results<-data.frame(n.train=rep(0,folds), n.test=0, nparams=0, train.auc=0,
                      test.auc=0, stringsAsFactors=FALSE)
  kvector <- kfold(covs.pres, folds)
  
  for (k in 1:folds){
    n.train <- length(which(kvector!=k))
    n.test <- length(which(kvector==k))
    train.df <- rbind(covs.pres[kvector!=k, ], covs.bkg.train)
    test.df <- rbind(covs.pres[kvector==k, ], covs.bkg.test)
    y.train <- c(rep(1, n.train),rep(0,nrow(covs.bkg.train)))
    y.test <- c(rep(1, n.test),rep(0,nrow(covs.bkg.test)))
    mxnt.obj <- maxent(x=train.df, p=y.train, removeDuplicates=FALSE, args=mxnt.args, path=path)
    pocc.train <- predict(mxnt.obj, train.df)
    pocc.test <- predict(mxnt.obj, test.df)
    auc.train <- evaluate(pocc.train[y.train==1], pocc.train[y.train==0])@auc
    auc.test <- evaluate(pocc.test[y.test==1], pocc.test[y.test==0])@auc
    nparams <- sum(getLambdaTable(mxnt.obj@lambdas)[, 2] != 0)
    results[k, ]<-c(n.train, n.test, nparams, auc.train, auc.test)
  }
  return(results)
}

#getLambdaTable
#Convert the lambda object in maxent to a data frame
#Arguments:
# lambdas(character vector): character vector returned by maxent function. This is usually obtained
#  by accessing the lambda slot in a maxent object, e.g. mxnt.obj@lambdas
#Returns:
# A data frame of lambda values, with columns feature name, coefficient, min value, max value.

getLambdaTable<-function(lambdas){
  lambdas.list <- strsplit(lambdas,",")
  nparams = length(lambdas) - 4
  varnames=rep("NA",nparams)
  result<-data.frame(lambdas=rep(0,nparams))
  for (i in 1:nparams){
    varnames[i]<-lambdas.list[[i]][1]
    result[i,1]<-as.numeric(lambdas.list[[i]][2])
  }
  result<-data.frame(varnames,result,stringsAsFactors=F)
  return(result)
}

#OptimizeLambda
#Optimize regularization value in maxent
#Arguments:
# folds(numeric): folds to use in k-fold partitioning
# covs.pres(data frame): covariates associated with occurrences
# covs.bkg.train (data frame): covariates associated with background used for training
# covs.bkg.train (data frame): covariates associated with background used for testing
# mxnt.args(string): string to be passed to dismo's function maxent though the args argument
# wd(string): working directory where result will be saved
# sp.prefix(string): prefix for output csv file name.
# path(string): string to be passed to dismo's function maxent though the path argument
#Returns:
# Saves in working directory the evaluation statistics at several regularization values
# for each folds and in R returns a vector with best and optimum regularization values
# and associated statistics (see FindBestLambda).

OptimizeLambda <- function(folds, covs.pres, covs.bkg.train, covs.bkg.test, mxnt.args, 
                           wd=getwd(), sp.prefix="species", path){
  lambda.vector <- c(0.02,0.05,0.1,0.22,0.46,1,2.2,4.6)
  results <-data.frame()
  for(lambda in lambda.vector){
    mxnt.args <- c(mxnt.args,paste0("betamultiplier=",lambda))
    results <- rbind(results, 
                     EvaluatePOModel(folds, covs.pres, covs.bkg.train, covs.bkg.test, mxnt.args,path=path))
  }
  results <- cbind(lambda=rep(lambda.vector,each=folds), results)
  write.csv(results, paste0(wd, "/", sp.prefix, "_lambda.optimization.csv"), row.names=FALSE)
  lambda.params <- FindBestLambda(results)
  return(lambda.params)
}

#FindBestLambda
#Use Mann-Whitney test to identify the best minimum regularization value.
#Arguments:
# df(data frame): data frame returned within the OptimizeLambda function (results object). This could
#  also be obtained by reading the csv file generated by the OptimizeLambda function.
#Returns:
# Data frame object with columns: best.lambda (best or maximum performance regularization value), 
# best.nparams (number of parameters when using best regularization value), best.median.auc 
# (median auc value of models run with best regularization value), optimum lambda (regularization value
# that has the largest value and is not significantly different from the best regularization value),
# optimum.nparams (number of parameters associated with optimum regularization value), 
# optimum.median.auc (auc of models developed with optimum regularization value).

FindBestLambda<-function(df){
  summary.df <- ddply(df, "lambda", summarise, mean.auc=mean(test.auc),
                      median.auc=median(test.auc, na.rm=TRUE),mean.nparams=mean(nparams))
  best.lambda <- summary.df$lambda[which.max(summary.df$median.auc)]
  optimum.lambda <- NA
  opt.idx <- NA
  start.idx <- which.max(summary.df$median.auc) + 1
  if(start.idx > nrow(summary.df)){
    result=c(best.lambda = best.lambda, 
             best.nparams = summary.df$mean.nparams[(start.idx-1)],
             best.median.auc = max(summary.df$median.auc),
             optimum.lambda = NA,
             optimum.nparams = NA,
             optimum.median.auc = NA)
    return(result)
  }
  for(i in start.idx:nrow(summary.df)){
    pval <- with(df,
                 wilcox.test(test.auc[lambda == best.lambda], test.auc[lambda == summary.df$lambda[i]],
                             alternative="greater", paired=F)$p.value)
    if(pval<0.05){
      opt.idx <- (i-1)
      optimum.lambda <- summary.df$lambda[opt.idx]
      break
    }
  }
  if(is.na(opt.idx)){
    result=c(best.lambda = best.lambda, 
             best.nparams = summary.df$mean.nparams[(start.idx-1)],
             best.median.auc = max(summary.df$median.auc),
             optimum.lambda = NA,
             optimum.nparams = NA,
             optimum.median.auc = NA)
  } else {
    result=c(best.lambda = best.lambda,
             best.nparams = summary.df$mean.nparams[(start.idx-1)],
             best.median.auc = max(summary.df$median.auc),
             optimum.lambda = optimum.lambda,
             optimum.nparams = summary.df$mean.nparams[opt.idx],
             optimum.median.auc = summary.df$median.auc[opt.idx])
  }
  return(result)
}

#EvaluateBRTModel

EvaluateBRTModel<-function(folds, covs.pres, covs.bkg.train, covs.bkg.test,brt.params){
  results<-data.frame(n.train=rep(0,folds), n.test=0, ntrees=0, lr=brt.params[2], train.auc=0,
                      test.auc=0, stringsAsFactors=FALSE)
  kvector <- kfold(covs.pres, folds)
  lr.tmp=brt.params[2]
  for (k in 1:folds){
    print(paste0("Starting with learning rate ",lr.tmp))
    n.train <- length(which(kvector!=k))
    n.test <- length(which(kvector==k))
    train.df <- rbind(covs.pres[kvector!=k, ], covs.bkg.train)
    test.df <- rbind(covs.pres[kvector==k, ], covs.bkg.test)
    y.train <- c(rep(1, n.train),rep(0,nrow(covs.bkg.train)))
    y.test <- c(rep(1, n.test),rep(0,nrow(covs.bkg.test)))
    
    df <- data.frame(y.train, train.df)
    
    brt.obj <- gbm.step(data=df, gbm.x = 2:ncol(df), gbm.y = 1,
                        family = "bernoulli", tree.complexity = brt.params[1], 
                        learning.rate = lr.tmp, bag.fraction = brt.params[3],prev.stratify=FALSE,
                        site.weights=c(rep(1,n.train), rep(n.train/nrow(covs.bkg.train),nrow(covs.bkg.train))))
    
    while(is.null(brt.obj)){
      lr.tmp=lr.tmp*0.5
      print(paste0("\nTrying learning rate ",lr.tmp))
      brt.obj <- gbm.step(data=df, gbm.x = 2:ncol(df), gbm.y = 1,
                          family = "bernoulli", tree.complexity = brt.params[1], 
                          learning.rate = lr.tmp, bag.fraction = brt.params[3],prev.stratify=FALSE,
                          site.weights=c(rep(1,n.train), rep(n.train/nrow(covs.bkg.train),nrow(covs.bkg.train))))
    }
    
    pocc.train <- predict(brt.obj, train.df, n.trees=brt.obj$gbm.call$best.trees, type="response")
    pocc.test <- predict(brt.obj, test.df, n.trees=brt.obj$gbm.call$best.trees, type="response")
    
    auc.train <- evaluate(pocc.train[y.train==1], pocc.train[y.train==0])@auc
    auc.test <- evaluate(pocc.test[y.test==1], pocc.test[y.test==0])@auc
    ntrees <- brt.obj$n.trees
    results[k, ]<-c(n.train, n.test, ntrees, lr.tmp, auc.train, auc.test)
  }
  return(results)
}


#EvaluateBioclimModel
EvaluateBioclimModel<-function(folds, covs.pres, covs.bkg.train, covs.bkg.test){
  results<-data.frame(n.train=rep(0,folds), n.test=0, train.auc=0,
                      test.auc=0, stringsAsFactors=FALSE)
  kvector <- kfold(covs.pres, folds)
  
  for (k in 1:folds){
    n.train <- length(which(kvector!=k))
    n.test <- length(which(kvector==k))
    train.df <- rbind(covs.pres[kvector!=k, ], covs.bkg.train)
    test.df <- rbind(covs.pres[kvector==k, ], covs.bkg.test)
    y.train <- c(rep(1, n.train),rep(0,nrow(covs.bkg.train)))
    y.test <- c(rep(1, n.test),rep(0,nrow(covs.bkg.test)))
    
    bc.obj <- bioclim(covs.pres[kvector!=k, ])
    
    pocc.train <- predict(bc.obj, train.df)
    pocc.test <- predict(bc.obj, test.df)
    
    auc.train <- evaluate(pocc.train[y.train==1], pocc.train[y.train==0])@auc
    auc.test <- evaluate(pocc.test[y.test==1], pocc.test[y.test==0])@auc
    results[k, ]<-c(n.train, n.test, auc.train, auc.test)
  }
  return(results)
}