R/PBSperformance.R

In MichaelFolkes/PBSperformance: Performance analysis methods, with a selection of metrics.

#' @title Calculate Z-score (standardized) values.
#'
#' @param x A vector of numbers.
#' @param mean A vector of length one. User can supply a mean value. If NA (the default), then it is calculated.
#' @param sd A vector of length one. User can supply an sd value. If NA (the default), then it is calculated.
#'
#' @return A vector, the same length as x consisting of the z-score values.
#' @export
#'
#' @examples
calcZscore <- function(x, mean=NA, sd=NA){
	x.mean <- ifelse(is.na(mean), mean(x, na.rm=TRUE), mean )
	x.sd <- ifelse(is.na(sd), sd(x, na.rm=TRUE), sd )

	z <- (x-x.mean)/x.sd
	return(z)
}#END calcZscore




#' @title Jackknife data and run a function on the subset of data.
#'
#' @param dat A data frame.
#' @param FUN A function.
#'
#' @return A list.
#' @export
#'
#' @examples
#' \dontrun{
#'
#' fn.tmp <- function(list.obj){
#'  	df <- list.obj$data
#'  	newdata <- list.obj$newdata
#'  	lm.fit <- lm(b~a, data=df)
#'  	predict.val <- predict(object = lm.fit, newdata=newdata)
#'  	error <- predict.val-newdata$b
#'  	return(list(model.fit=lm.fit, predict.val=predict.val, error=error))
#'  	}#END fn.tmp
#'
#' df <- data.frame(a=1:10, b=rnorm(10))
#'
#' out.list <- jackknife(df, FUN = fn.tmp)
#' out.list[[1]]
#' }
jackknife <- function(dat, FUN=NULL){
	dat <- as.data.frame(dat)
	fun.output <- vector("list", nrow(dat))

	for(i in 1:nrow(dat)){
		fun.output[[i]]$data <- dat[-i,,drop=FALSE]
		fun.output[[i]]$newdata <- dat[i,,drop=FALSE]
		if(!is.null(FUN)){
			fun.output[[i]] <- c(fun.output[[i]], FUN(fun.output[[i]]))
			}
	}

	return(fun.output)
}#END jackknife


makePMdataframe <- function(modelNames){
  Performance.out<-data.frame(row.names=modelNames, mre=NA, mae=NA, mrse=NA, ustat2=NA, mase=NA)
  return   (Performance.out)
}#END MakePMdataframe




#' @title Performance measure wrapper
#'
#' @description Here is the description.
#'
#' @inheritParams PMs
#' @param metrics A vector of the performance metrics to be included.
#'
#' @return A matrix of named columns, one for each metric, and the number of
#' data pairs used in the calculations.
#' NA values are tolerated.
#' @export
performance <- function(expect,obs, metrics=c('all', 'mre', 'mae', 'mpe', 'mape', 'rmse', 'ustat2', 'ustat2.mean', 'mase'),...){
  if('all' %in% metrics) metrics <- c( 'mre', 'mae', 'mpe', 'mape', 'rmse', 'ustat2', 'ustat2mean', 'mase')

  pm.list <- lapply(metrics, do.call, list(expect, obs,...))
  names(pm.list) <- metrics

  pm.vec <- unlist(lapply(pm.list, "[[", 2))
  m <- data.frame(matrix(ncol=length(metrics)+1, data=c(pm.vec,  sum(complete.cases(expect,obs))  ),dimnames=list(NULL,c(metrics , "number.forecasts"))))

  return(list(errors=pm.list[[1]]$errors, m = m))
}#END Performance




#' @title Performance measures
#'
#' @param expect A number vector, representing the expected or forecasted
#'   values.
#' @param obs A number vector, representing the observed values.
#' @param layman A \code{logical}. Used in \code{\link{mre}} and \code{\link{mpe}}. If
#'   FALSE (default) then the numerator is defined as expected - observed. If
#'   TRUE the numerator is defined as observed - expected. The default method is
#'   the approach commonly used in forecasting literature (See Hyndman, 2006).
#'   The latter approach produces results that may be more intuitive to the
#'   layman...
#'
#' @return A \code{list}, comprising two elements is returned. The first
#'   element, named \code{errors}, is a vector of the index-specific errors. The
#'   second list element is the statistic. NA values are tolerated.
#' @seealso \url{https://en.wikipedia.org/wiki/Mean_absolute_error}
#' @name PMs
NULL




#' @title Calculate mean absolute error.
#' @inheritParams PMs
#'
#' @export
mae <- function(expect,obs,...) {
  #complete <- complete.cases(expect, obs)
  #expect <- expect[complete]
  #obs <- obs[complete]
  errors <- (abs(expect - obs))
  errors.mean <- mean(errors, na.rm = TRUE)
  return(list(errors=errors, mae=errors.mean))

}#END mae


#' @title Calculate mean absolute percent error.
#'
#' @inheritParams PMs
#' @export
mape <- function(expect,obs,...) {
  #complete <- complete.cases(expect, obs)
  #expect <- expect[complete]
  #obs <- obs[complete]
  errors <- (abs((expect - obs)/obs))
  errors.mean <- mean(errors[!is.infinite(errors)], na.rm = TRUE)
  list(errors=errors, mape=errors.mean)
}#END mape


#' @title Calculate mean absolute scaled error.
#'
#' @inheritParams PMs
#' @export
mase <- function(expect,obs,...){
  #Hyndman,RJ (2006)
  #complete <- complete.cases(expect, obs)
  #expect <- expect[complete]
  #obs <- obs[complete]
  res <- obs-expect #note it's opposite what is typical
  scale <- mean(abs(diff(obs)), na.rm = TRUE)
  mase <- mean(abs(res/scale), na.rm = TRUE)
  return(list(errors=NULL, mase=mase))
}#END mase



#' @title Calculate mean percent error.
#'
#' @inheritParams PMs
#'
#' @export
mpe <- function(expect,obs, layman=FALSE,...) {
  #complete <- complete.cases(expect, obs)
  #expect <- expect[complete]
  #obs <- obs[complete]
	if(layman){
		errors <- (expect - obs)/obs
	}else{
		errors <- (obs - expect)/obs
	}
  errors.mean <- mean(errors[!is.infinite(errors)], na.rm = TRUE)
  return(list(errors=errors, mpe=errors.mean))

}#END mpe



#' @title Calculate mean raw error.
#'
#' @inheritParams PMs
#'
#' @export
mre <- function(expect,obs, layman=FALSE,...) {
  #complete <- complete.cases(expect, obs)
  #expect <- expect[complete]
  #obs <- obs[complete]
	if(layman){
		errors <- (expect - obs)
	}else{
		errors <- (obs - expect)
	}
  errors.mean <- mean(errors[!is.infinite(errors)], na.rm = TRUE)
  return(list(errors=errors, mre=errors.mean))
}#END mre


#' @title Calculate mean squared error.
#'
#' @inheritParams PMs
#' @export
mse<-function(expect,obs,...) {
	#complete <- complete.cases(expect, obs)
	#expect <- expect[complete]
	#obs <- obs[complete]
	errors <- (expect - obs)^2
	errors.mean <- mean(errors, na.rm = TRUE)
	list(errors=errors, mse=errors.mean)
}#END mse





#' @title Calculate sqaure root of mean squared error.
#'
#' @inheritParams PMs
#' @export
rmse<-function(expect,obs,...) {
  #complete <- complete.cases(expect, obs)
  #expect <- expect[complete]
  #obs <- obs[complete]
  errors <- (expect - obs)^2
  errors.mean <- sqrt(mean(errors, na.rm = TRUE))
  list(errors=errors, rmse=errors.mean)
}#END rmse



#' @title Calculate Theil's U-statistic.
#' @inheritParams PMs
#'
#' @export
ustat2 <- function(expect,obs, benchmark=c("lly", "mean", "median"), window=NA,...){
  #Theil's U-statistic (Thiel, H 1966. Applied economic forecasting)
  #Hyndman,RJ (2006)
  # this fn uses one of three possible benchmark methods (like last year, mean, medain). The latter two can be based on a window size of chosen length (eg window=3 calcs mean or median on 3 obs points prior to current forecast). if window=NA, then use all data prior to current
  benchmark <- match.arg(benchmark)
  lly <- function(x,...){x[length(x)]}

  complete <- complete.cases(expect, obs)
  expect.bm <- expect
  for(ind in 1:length(expect.bm)){
    start <- ifelse(is.na(window), 1, max(1,ind-window ))
    X <- obs[start:(ind-1)]
    if(!is.na(expect.bm[ind])){
      expect.bm[ind] <- do.call(what = benchmark, args = list(x=X, na.rm=TRUE))
    }
  }

  expect.bm <- expect.bm[complete]
  expect <- expect[complete]
  obs <- obs[complete]
  numer <- rmse(expect, obs)
  denom <- rmse(expect =  expect.bm, obs =  obs )
  u2 <- numer/denom

  return(list(errors=NULL, u2=u2))
}#END ustat2



#' @title Calculate Theil's U-statistic standardized by benchmark method.
#'
#' @inheritParams PMs
#' @export
ustat2mean <- function(expect,obs,  window=NA,...){
  #Theil's U-statistic (Thiel, H 1966. Applied economic forecasting)
  #Hyndman,RJ (2006)
  # this fn uses one of three possible benchmark methods (like last year, mean, medain). The latter two can be based on a window size of chosen length (eg window=3 calcs mean or median on 3 obs points prior to current forecast). if window=NA, then use all data prior to current
  benchmark <- "mean"
  lly <- function(x,...){x[length(x)]}

  complete <- complete.cases(expect, obs)
  expect.bm <- expect
  for(ind in 1:length(expect.bm)){
    start <- ifelse(is.na(window), 1, max(1,ind-window ))
    X <- obs[start:(ind-1)]
    if(!is.na(expect.bm[ind])){
      expect.bm[ind] <- do.call(what = benchmark, args = list(x=X, na.rm=TRUE))
    }
  }

  expect.bm <- expect.bm[complete]
  expect <- expect[complete]
  obs <- obs[complete]
  numer <- rmse(expect, obs)
  denom <- rmse(expect =  expect.bm, obs =  obs )
  u2.mean <- numer/denom

  return(list(errors=NULL, u2.mean=u2.mean))
}#END ustat2.mean




#' @title Calculate Theil's U-statistic. Variant.
#'
#' @inheritParams PMs
#' @export
ustat2rmspe <- function(expect,obs,...){
  #Theil's U-statistic (Thiel, H 1966. Applied economic forecasting)
  #Hyndman,RJ (2006) describes that there are ambiguities to who defines ustat2. Makridakis et (1998) say U2 is the relative RMSPE as shown on this webpage:
  #http://docs.oracle.com/cd/E40248_01/epm.1112/cb_statistical/frameset.htm?ch07s02s03s04.html
  #relative RMSPE is what's shown here:
  complete <- complete.cases(expect, obs)
  expect <- expect[complete]
  obs <- obs[complete]
  numer <- sqrt(mean(((expect[-1] - obs[-1]) / obs[-length(obs)])^2))
  denom <- sqrt(mean(((obs[-1] - obs[-length(obs)]) / obs[-length(obs)])^2))
  u2RMSPE <- numer/denom

  return(list(errors=NULL, u2RMSPE=u2RMSPE))
}



#' @title Calculate Theil's U-statistic. Variant.
#'
#' @inheritParams PMs
#' @export
ustat2b <- function(expect,obs,...){
  #this version expects the obs vector to have one data point earlier than expect vector, so like last step forecast can work properly
  #first element of expect would be an NA
  #Theil's U-statistic (Thiel, H 1966. Applied economic forecasting)
  #Hyndman,RJ (2006)
  # complete <- complete.cases(expect, obs)

  numer <- rmse(expect, obs)
  denom <- rmse(obs[-length(obs)], obs[-1] ) #like last step model
  u2 <- numer/denom

  return(list(errors=NULL, ustat2b=ustat2b))
}#END ustat2b