Quatifying Spatial Pattern in Ecological Data

Accounting for Spatial Dependence in Ecological Data

Reference and Disclaimer: For this exercise, I have used explanation and major chunk of codes from the book Spatial Ecology and Conservation Modeling by Robert Fletcher & Marie-Josée Fortin.

Load the necessary packages

rm(list = ls())

## First specify the packages of interest
packages = c("stringr","rgdal", "spdep", "lme4","vegan", "mgcv", "MASS", "spaMM", "deldir", "dismo", "spatialreg", "huxtable","INLA")

## Now load or install&load all
package.check <- lapply(
  packages,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE)
      library(x, character.only = TRUE)
    }
  }
)

Load the raster data and the point data to be analysed

elev <- raster("./SpaStatAssign/CWT_regiona_10mDEM.tif")

point.data = readOGR("./SpaStatAssign/SpaStatData.gdb","BAWW2021_PointCounts800mgrid")

Checking the projection for the DEM layer

proj4string(elev)

## [1] "+proj=utm +zone=17 +datum=NAD83 +units=m +no_defs"

Checking the projection for the point data

proj4string(point.data)

## [1] "+proj=utm +zone=17 +datum=NAD83 +units=m +no_defs"

head(point.data)

ID	x1	x2	PointCountCapName	detectP1	detectP2	detectP3	detectP4	Detection
pcCoweeta001	2.76e+05	3.88e+06	PCCow001	0	0	1	0	1
pcCoweeta002	2.75e+05	3.88e+06	PCCow002	1	1	0	0	1
pcCoweeta002	2.75e+05	3.88e+06	PCCow002	1	1	1	0	1
pcCoweeta003	2.75e+05	3.88e+06	PCCow003	0	0	0	1	1
pcCoweeta004	2.76e+05	3.88e+06	PCCow004	1	0	1	0	1
pcCoweeta005	2.76e+05	3.88e+06	PCCow005	1	1	1	1	1

Using terrain() we derive other variables of interest: slope and aspect. We merge all the layers into a single raster stack using stack()

#create slope, aspect layers from elevation
elev.terr <- terrain(elev, opt=c("slope", "aspect"))
#make a multilayered file for extraction
layers <- stack(elev, elev.terr) #merge the slope and aspect layers into a single raster object that holds all raster layers
names(layers) <- c("elev", "slope", "aspect")

Plotting the raster stack

plot(layers)

Plotting the presence and absence of Black and White Warbler. We are concerned with the detection of the BAWW rather than the type of BAWW detected here

plot(elev)
points(point.data[point.data$Detection==1, c("x1","x2")], col="red")      #plots the presence points
points(point.data[point.data$Detection==0, c("x1","x2")], col="black")    #plots the absence points

We use extract() from raster package to extract covariate values from layers corresponding to the locations from the survey data of BAWW and merge them with the point data

#extract GIS data at sampling points
coords <- cbind(point.data$x1, point.data$x2)
colnames(coords) <- c("x", "y")
land.cov <- raster::extract(x=layers, y=coords) # If y represents points, extract returns the values of a Raster* object for the cells in which a set of points fall.
point.data <- cbind(point.data,land.cov)

Since centering and scaling can help improve model convergence and facilitates comparing coefficients for different parameter, we scale our predictors

point.data$elevs <- scale(point.data$elev, center = T, scale = T)
point.data$slopes <- scale(point.data$slope, center = T, scale = T)
point.data$aspects <- scale(point.data$aspect, center = T, scale = T)

Fitting Logistic Models of varying complexity

model_elev <- glm(Detection ~ elev, family="binomial", data=point.data)
model_slope <- glm(Detection ~ slope, family="binomial", data=point.data)
model_all <- glm(Detection ~ elev + slope + aspect, family="binomial", data=point.data)
model_slope2 <- glm(Detection ~ poly(slope,2),family="binomial", data=point.data)

model_tbl = list("Elevation only" = model_elev,
                 "Slope only" = model_slope,
                 "All" = model_all,
                 "Slope square" = model_slope2)

coefs <- names(coef(model_tbl[[1]]))[stringr::str_detect(names(coef(model_tbl[[1]])), "vdc")]

huxtable <- huxreg(model_tbl,number_format = 3, omit_coefs =coefs) %>%
 set_caption("Table: Logistic Regression Models")
huxtable

Table: Logistic Regression Models

	Elevation only	Slope only	All	Slope square
(Intercept)	0.676	-0.859	-1.536	1.014 ***
	(1.192)	(0.571)	(1.509)	(0.204)
elev	0.000		0.001
	(0.001)		(0.001)
slope		4.325 **	4.384 **
		(1.359)	(1.441)
aspect			-0.183
			(0.110)
poly(slope, 2)1				7.822 **
				(2.581)
poly(slope, 2)2				-0.032
				(2.494)
N	141	141	141	141
logLik	-84.070	-78.454	-76.608	-78.454
AIC	172.140	160.908	161.216	162.908
*** p < 0.001; ** p < 0.01; * p < 0.05.

From the above model comparison we have some evidence that the BAWW spotting increases with increasing slope. We see that the slope term is positive and significant at 1%. Considering the AIC along side the significance, I would prefer Slope only model as the number of variable used is smaller than other models.

Predicting the relationship using new dataset

We create a new dataset and use Slope only model to predict the relationship

Slope <- seq(min(point.data$slope), max(point.data$slope), length=15)
newdata <- data.frame(slope=Slope)
#Predict onto newdata
glm.pred <- predict(model_slope, newdata=newdata, type= "link", se=T) 
# Creating 95% confidence interval
ucl <- glm.pred$fit + 1.96*glm.pred$se.fit
lcl <- glm.pred$fit - 1.96*glm.pred$se.fit

Back-transforming the predictions from link to the probability scale

glm.newdata <- data.frame(newdata, pred=plogis(glm.pred$fit), lcl=plogis(lcl), ucl=plogis(ucl))

Plotting the predictions of the model onto new data in probability scale

plot(glm.newdata$slope, glm.newdata$pred, ylim=c(0,2), xlab="Slope", ylab="Prob. Occurrence")
lines(glm.newdata$slope, glm.newdata$lcl)
lines(glm.newdata$slope, glm.newdata$ucl)

Plotting the predictions of the model onto the raster stack layer

#Map the model
glm.raster <- predict(model=model_slope, object=layers,  type="response")
glm.raster <= exp(glm.raster)/(1+exp(glm.raster))

## class      : RasterLayer 
## dimensions : 952, 1285, 1223320  (nrow, ncol, ncell)
## resolution : 9.236852, 9.236852  (x, y)
## extent     : 269816.1, 281685.4, 3876588, 3885381  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=17 +datum=NAD83 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 1  (min, max)

plot(glm.raster, xlab = "Longitude", ylab = "Latitude")

We say that there exist spatial dependence: Spotting BMWW with increasing slope

Thanks!!!