In an interesting recent paper, Royle et al. (2012) presented a critical analysis of Maxent and proposed a method for obtaining probability of occurrence (rather than a vague suitability index) based on maximum likelihood from presence-background (rather than presence-absence) species distribution data. They provided an R package, maxlike, which can apply this method to a set of presence point coordinates and a raster stack of predictor variables. However, this implies not only that you need point coordinates of your species observations, but also that sampling units are square pixels, whereas we often need to model distribution data recorded on larger and more irregular units — such as provinces, river basins, or UTM cells, for example.
The maxlike.df function is a modification of maxlike that can be applied to a data frame instead of spatial objects. I added argument covariates (to use instead of rasters), which requires a matrix or data.frame of the variables (one in each column, with the same names used in argument formula), and argument presences (to use instead of points) for a vector/column with the species data (of the same size and in the same order as the covariates, with 1 for presence and 0 or NA for background or no records — mind that MaxLike, like MaxEnt, needs both presence and background data, so they’re not strictly “presence-only”). Everything else in maxlike.df works the same way as in maxlike (Chandler & Royle 2013), so please read the maxlike documentation for further information on how the method is implemented and interpreted. Note that they recommend that the variables be standardized (e.g. vars = scale(vars)) before modelling, and that you should select which variables to include in each model.
Note also that Royle et al.‘s (2012) approach has been supported by Fitzpatrick et al. (2013), but criticized by Hastie & Fithian (2013).
UPDATE: the latest version of the maxlike package (0.1-7) now allows models to be calibrated with data frames, via arguments x (data.frame with environmental conditions at presence point locations) and z (data.frame with environmental conditions at background locations). The function below is therefore deprecated, and you should use the maxlike package instead.
maxlike.df <- function (formula, rasters, points, presences, covariates, link = c("logit", "cloglog"), starts, hessian = TRUE, fixed, removeDuplicates = FALSE, savedata = FALSE, na.action = "na.omit", ...)
{
if (missing(rasters)) rasters <- covariates # new line
if (missing(points)) points <- presences # new line
if (identical(formula, ~1))
stop("At least one continuous covariate must be specified in the formula")
link <- match.arg(link)
varnames <- all.vars(formula)
call <- match.call()
npts <- sum(points == 1) # was nrow(points)
cd.class <- class(rasters)[1]
#if (cd.class != "rasterstack") stop("rasters must be a raster stack")
#pt.class <- class(points)[1]
#if (!pt.class %in% c("matrix", "data.frame")) stop("points must be a matrix or a data.frame")
#if (ncol(points) != 2) stop("points must have 2 columns containing the x- and y- coordinates")
pt.names <- colnames(points)
#if (identical(cd.class, "rasterstack")) {
cd.names <- names(rasters)
npix <- nrow(rasters) # was prod(dim(rasters)[1:2])
id <- 1:nrow(rasters) # new line
cellID <- id[points == 1] # was cellFromXY(rasters, points)
duplicates <- duplicated(cellID)
if (removeDuplicates) {
cellID <- unique(cellID)
npts <- length(cellID)
points.retained <- points[!duplicates, ]
}
x <- as.data.frame(rasters[points == 1, ]) # was as.data.frame(matrix(extract(rasters, cellID), npts))
z <- rasters # was as.data.frame(matrix(getValues(rasters), npix))
names(x) <- names(z) <- cd.names
#}
if (!all(varnames %in% cd.names)) stop("at least 1 covariate in the formula is not in names(rasters).")
X.mf <- model.frame(formula, x, na.action = na.action)
X.mf.a <- attributes(X.mf)
pts.removed <- integer(0)
points.retained <- points
if ("na.action" %in% names(X.mf.a)) {
pts.removed <- X.mf.a$na.action
npts.removed <- length(pts.removed)
if (npts.removed > 0) {
warning(paste(npts.removed, "points removed due to missing values"))
points.retained <- points.retained[-pts.removed] # I removed the comma before ']', as here 'points' is a vector rather than a two-column matrix/dataframe
}
}
X <- model.matrix(formula, X.mf)
Z.mf <- model.frame(formula, z, na.action = na.action)
Z.mf.a <- attributes(Z.mf)
pix.removed <- integer(0)
if ("na.action" %in% names(Z.mf.a)) {
pix.removed <- Z.mf.a$na.action
npix.removed <- length(pix.removed)
}
Z <- model.matrix(formula, Z.mf)
npars <- ncol(X)
parnames <- colnames(X)
if (!"(Intercept)" %in% parnames)
stop("The intercept must be estimated or fixed")
if (missing(starts)) {
starts <- rep(0, npars)
names(starts) <- parnames
}
else names(starts) <- parnames
if (identical(link, "logit")) {
nll <- function(pars) {
psix <- plogis(drop(X %*% pars))
psiz <- sum(plogis(drop(Z %*% pars)))
-1 * sum(log(psix/psiz))
}
}
else if (identical(link, "cloglog")) {
nll <- function(pars) {
psix <- 1 - exp(-exp(drop(X %*% pars)))
psiz <- sum(1 - exp(-exp(drop(Z %*% pars))))
-1 * sum(log(psix/psiz))
}
}
else stop("link function should be either 'logit' or 'cloglog'")
is.fixed <- rep(FALSE, npars)
if (!missing(fixed)) {
if (length(fixed) != length(starts))
stop("fixed should be a vector with the same length as the number of parameters to be estimated")
if (sum(is.double(fixed)) < 1)
stop("fixed must contain at least one real value")
is.fixed <- !is.na(fixed)
if (sum(!is.fixed) < 1)
stop("you cannot fix all parameters in the model")
npars <- sum(!is.fixed)
nll.fix <- function(p) {
p[is.fixed] <- fixed[is.fixed]
do.call("nll", list(pars = p))
}
fm <- optim(starts, nll.fix, hessian = hessian, ...)
fm$par[is.fixed] <- fixed[is.fixed]
}
else {
fm <- optim(starts, nll, hessian = hessian, ...)
}
not.fixed <- !is.fixed
par <- fm$par
if (hessian) {
vcTry <- try(solve(fm$hessian[not.fixed, not.fixed]))
if (identical(class(vcTry), "matrix")) {
vc <- matrix(0, length(par), length(par))
vc[not.fixed, not.fixed] <- vcTry
se <- sqrt(diag(vc))
}
else {
vc <- matrix(NA, npars, npars)
se <- rep(NA, npars)
}
}
else {
vc <- matrix(NA, npars, npars)
se <- rep(NA, npars)
}
dimnames(vc) <- list(parnames, parnames)
aic <- 2 * fm$value + 2 * npars
fitted <- plogis(Z %*% par)
fitted.values <- merge(data.frame(id = id), data.frame(id = rownames(fitted), fitted = fitted), all = TRUE) # new line
out <- list(Est = cbind(Est = par, SE = se), vcov = vc, AIC = aic, call = call, pts.removed = pts.removed, pix.removed = pix.removed, points.retained = points.retained, optim = fm, not.fixed = not.fixed, link = link, fitted.values = fitted.values$fitted) # 'fitted.values' added
if (class(rasters) %in% c("matrix", "data.frame")) savedata = TRUE # new line, otherwise 'predict' would fail
if (savedata)
out$rasters <- rasters
class(out) <- c("maxlikeFit", "list")
return(out)
}
# some changes to predict.maxlikeFit were also required to accomodate maxlike.df:
predict.maxlikeFit <- function (object, ...) {
e <- coef(object)
rasters <- object$rasters
if (is.null(rasters)) {
rasters <- try(get(as.character(object$call$rasters)))
if (identical(class(rasters)[1], "try-error"))
stop("could not find the raster data")
warning("raster data were not saved with object, using the data found in the workspace instead.")
}
link <- object$link
cd.names <- names(rasters)
cd.class <- class(rasters)[1] # new line
if (cd.class == "RasterStack") { # new line
npix <- prod(dim(rasters)[1:2])
z <- as.data.frame(matrix(getValues(rasters), npix))
} # new line
else if (cd.class %in% c("matrix", "data.frame")) { # new line
npix <- nrow(rasters) # new line
z <- rasters # new line
} else stop("'rasters' must be either a RasterStack, a matrix or a data.frame") # new line
names(z) <- cd.names
formula <- object$call$formula
varnames <- all.vars(formula)
if (!all(varnames %in% cd.names))
stop("at least 1 covariate in the formula is not in rasters.")
Z.mf <- model.frame(formula, z, na.action = "na.pass")
Z.terms <- attr(Z.mf, "terms")
Z <- model.matrix(Z.terms, Z.mf)
eta <- drop(Z %*% coef(object))
if (identical(link, "logit"))
psi.hat <- plogis(eta)
else if (identical(link, "cloglog"))
psi.hat <- 1 - exp(-exp(eta))
else stop("link function should be either 'logit' or 'cloglog'")
if (cd.class == "RasterStack") { # new line
psi.mat <- matrix(psi.hat, dim(rasters)[1], dim(rasters)[2], byrow = TRUE)
psi.raster <- raster(psi.mat)
extent(psi.raster) <- extent(rasters)
psi.raster
} # new line
else if (cd.class == "data.frame") return(psi.hat) # new line
}
[presented with Pretty R]
Usage example for a data frame with the species presence data in column 2 and the variables (named var1, var2, var3 and var4) in columns 3 to 6:
myspecies.maxlike <- maxlike.df(formula = ~ var1 + var2 + var3 + var4, covariates = mydata[ , 3:6], presences = mydata[ , 2], method = "BFGS")
If you then want to explore the model with functions such as summary or predict, you need to install and load the maxlike package (which defines the appropriate class for this type of model). If not, just with the maxlike.df function you can already get some useful information from the model:
# get the model predictions applied to the analysed data:
myspecies.maxlike$fitted.values
# get the coefficient estimates and their standard errors:
myspecies.maxlike$Est
# get the model's AIC:
myspecies.maxlike$AIC
You can then evaluate maxlike.df predictions with functions from the modEvA package, as in the following examples:
library(modEvA)
AUC(obs = mydata$myspecies, pred = myspecies.maxlike$fitted.values)
HLfit(obs = mydata$myspecies, pred = myspecies_maxlike$fitted.values)
plotGLM(obs = mydata$myspecies, pred = myspecies_maxlike$fitted.values)
References
Chandler, R.B. & Royle, J.A. (2013) Package ‘maxlike’: Model species distributions by estimating the probability of occurrence using presence-only data. Version 0.1-4, CRAN
Fitzpatrick, M.C., Gotelli, N.J. & Ellison, A.M. (2013) MaxEnt versus MaxLike: empirical comparisons with ant species distributions. Ecosphere 4, art55.
Hastie, T. & Fithian, W. (2013) Inference from presence-only data; the ongoing controversy. Ecography 36, 864-867.
Royle, J.A., Chandler, R.B., Yackulic, C. & Nichols, J.D. (2012) Likelihood analysis of species occurrence probability from presence-only data for modelling species distributions. Methods in Ecology and Evolution 3: 545–554
modTools 
