2 Mapping Occurrence Probability
In this example, we use the occu function to fit the single-season occupancy model of MacKenzie et al. (2002) to data on the European crossbill (Loxia curvirostra) collected in 267 1-km\(^2\) sample quadrats in Switzerland, 1999 (Schmid et al. 2004). We then use the model to compute the expected probability of occurrence at each pixel in a raster defining the Swiss landscape.
First we load the crossbill data, which contains the detection/non-detection data and covariates. The dataset actually contains data from multiple years, but we are only going to analyze data from the first year, 1999. A multi-year analysis of occupancy dynamics could be accomplished using the colext function, and in addition to mapping occurrence probability, it would be possible to also map colonization and extinction probabilities. The following commands format the data.
library(unmarked)
data(crossbill)
umf <- unmarkedFrameOccu(
y=as.matrix(crossbill[,c("det991", "det992", "det993")]),
siteCovs=crossbill[,c("ele", "forest")],
obsCovs=list(date=crossbill[,c("date991", "date992", "date993")]))
sc <- scale(siteCovs(umf))
siteCovs(umf) <- sc
head(umf)## Data frame representation of unmarkedFrame object.
## y.1 y.2 y.3 ele forest date.1 date.2 date.3
## 1 0 0 0 -1.15390803 -1.14710440 34 59 65
## 2 0 0 0 -1.15390803 -0.49668159 17 33 65
## 3 NA NA NA -0.21745379 -0.09920098 NA NA NA
## 4 0 0 0 -0.37352950 -0.93029679 29 59 65
## 5 0 0 0 -0.06137809 0.00920282 24 45 65
## 6 NA NA NA -0.99783232 -1.18323900 NA NA NA
## 7 0 0 0 -0.68568091 -1.03870060 26 54 74
## 8 0 0 0 -0.84175662 0.91256783 23 43 71
## 9 0 0 0 -0.99783232 -1.07483520 21 36 56
## 10 0 0 0 -0.99783232 -0.78575839 37 62 75Notice that the site covariates, elevation and forest, were standardized using the scale function. Standardization isn’t always necessary, but it can make it easier to find the maximum likelihood estimates. When standardizing covariates and then making predictions, it is important to retain the original sample mean and standard deviation. The reason for this is explained below.
Fitting a model is now straight-forward. In many cases, we would fit several models corresponding to competing hypotheses, but for simplicity, we stick with this single model.
(fm.occu <- occu(~date ~ele + I(ele^2) + forest, umf))## Warning: 22 sites have been discarded because of missing data.##
## Call:
## occu(formula = ~date ~ ele + I(ele^2) + forest, data = umf)
##
## Occupancy:
## Estimate SE z P(>|z|)
## (Intercept) 0.464 0.541 0.857 0.39137
## ele 1.228 0.319 3.846 0.00012
## I(ele^2) -1.333 0.417 -3.194 0.00140
## forest 0.654 0.274 2.385 0.01706
##
## Detection:
## Estimate SE z P(>|z|)
## (Intercept) -2.2356 0.52283 -4.28 0.000019
## date 0.0254 0.00752 3.39 0.000711
##
## AIC: 455.3892Now that we have our fitted model, we can use it to predict occurrence probability at each pixel in the Swiss landscape. The Switzerland dataset contains country-wide data. There are many ways to display it—here is an example of mapping elevation using the levelplot function in the lattice package (Sarkar 2008).
library(lattice)
data(Switzerland)
print(levelplot(elevation ~ x + y, Switzerland, aspect="iso",
xlab="Easting (m)", ylab="Northing (m)",
col.regions=terrain.colors(100)))
Figure 1. Elevation in Switzerland
The raster package (Jacob van Etten 2012) provides another alternative. Here we create two raster objects and specify the coordinate system.
library(raster)elevation <- rasterFromXYZ(Switzerland[,c("x","y","elevation")],
crs="+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs")
forest <- rasterFromXYZ(Switzerland[,c("x","y","forest")],
crs="+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs")Since we standardized the covariates during the model fitting process, we need to transform the country-wide data using the same values. Note, we don’t want to use the mean and SD of the rasters themselves, we want to use the mean and SD of the original covariates used to fit the models, which are stored as attributes of the sc object. The following commands display the original means and SDs and then transform the rasters and join them in a raster “stack”.
attr(sc, "scaled:center")## ele forest
## 1189.32584 34.74532attr(sc, "scaled:scale")## ele forest
## 640.71471 27.67431ele.s <- (elevation-1189)/640
forest.s <- (forest-34.7)/27.7
ef <- stack(ele.s, forest.s)
names(ef) <- c("ele", "forest")
plot(ef, col=terrain.colors(100))
Figure 2. Elevation and forest cover, standardized
It is important to assign names that exactly match the covariate names used to fit the model. This is required by the predict function as demonstrated later. The predict function is useful for computing spatially-referenced model predictions, standard errors, and confidence intervals, but it is computationally demanding when there are many pixels in the raster. Thus, if measures of uncertainty are not required, the following code can be used to quickly produce the species distribution map shown in Fig. 3.
(beta <- coef(fm.occu, type="state"))## psi(Int) psi(ele) psi(I(ele^2)) psi(forest)
## 0.4638233 1.2276426 -1.3327186 0.6544976logit.psi <- beta[1] + beta[2]*ele.s + beta[3]*ele.s^2 + beta[4]*forest.s
psi <- exp(logit.psi) / (1 + exp(logit.psi))
print(spplot(psi, col.regions=terrain.colors(100)))
Figure 3. A species distribution map for the European crossbill in Switzerland. The colors represent occurrence probability.
As of version 0.9-6, the predict method in unmarked can make predictions using an object of class RasterStack from the package. As mentioned previously, the rasters must be named, perhaps by using the names(someraster) <- somename method. The object returned by predict is another raster stack with rasters for the expected values of the parameter of interest, the standard errors, and the upper and lower confidence intervals. The following example is very slow because there are many of pixels in the raster. The resulting map is shown in Fig. 4.
E.psi <- predict(fm.occu, type="state", newdata=ef)
plot(E.psi, axes=FALSE, col=terrain.colors(100))
Figure 4. Expected occurrence probability along with standard errors and the limits of the asymptotic 95% confidence interval.
Users should be cautious when predicting from models that have categorical predictor variables, i.e. factors. The package does not have advanced methods for handling factors, and thus it is not easy to automatically create dummy variables from them as can typically be done using model.matrix. The safest option is to create the dummy variables manually before fitting the models, and to use the same variables as rasters for prediction.
A more important consideration when creating species distribution maps based upon occurrence probability is that of spatial scale. Occurrence probability will typically depend upon the area of the “site” in question. Thus, in our crossbill example, it would not be appropriate to use our model to predict occurrence probability for 10-km\(^2\) pixels since the surveys were done in 1-km\(^2\) quadrats. In some cases it might be possible to directly model the effect of site area on occurrence probability, in which case the effect could be accounted for in the predictions.
