Resource and Step-Selection Functions
Elie Gurarie
August 26, 2016
- 1 Resource Selection
- 1.1 Important Packages
- 1.2 Resource selection: Mule Deer Habitat
- 1.3 Determine available space (simplistically: MCP)
- 1.4 Create data frame of covariates
- 1.5 Obtain some random points within mcp
- 1.6 Calculate the covariates for the null set
- 1.7 Combine…
- 1.8 … and Model!
- 1.9 Model Selection
- 1.10 Predict from model
- 2 Step Selection
1 Resource Selection
- Some habitats are preferred over other habitats
- Resource Selection Functions quantify that preference and avoidance.
- The strategy is:
- Find used locations
- Find available locations
- Perform logistic regression
- The RSF: \[w(X) = \exp(\beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + ...)\]
- Positive \(\beta\)’s = preference
- Negative \(\beta\)’s = avoidance
Nicely … generate a map.
1.1 Important Packages
require(sp) # - for spatial points
require(adehabitatHR) # - for MCP's
require(raster) # - for loading and working with rasters
require(plyr)
require(magrittr)1.2 Resource selection: Mule Deer Habitat
deer.locs <- read.csv("rsf/deer.csv")
elevation <- raster("rsf/elevation.asc")
shrub <- raster("rsf/shrub.asc")
barren <- raster("rsf/barren.asc")
roads <- raster("rsf/rds.asc")
wells <- raster("rsf/wells.asc")
save(deer.locs, elevation, shrub, barren, roads, wells, file = "deerrasters.robj")p <- function(alpha=0.3) points(deer.locs, cex=0.5, col=rgb(0,0,1,alpha), pch=19)
par(mfrow=c(2,3), cex.lab = 0.5, mar = c(1,1,1,1), oma = c(5,3,1,1), xpd=NA, bty="n")
plot(elevation, main = "elevation", xaxt="n", xlab=""); p()
plot(wells, main = "wells", xaxt="n", xlab="", yaxt="n", ylab=""); p()
plot(shrub, main = "shrub", xaxt="n", xlab="", yaxt="n", ylab=""); p()
plot(barren, main = "barren"); p()
plot(roads, main = "roads", yaxt="n", ylab=""); p()
1.3 Determine available space (simplistically: MCP)
deer.mcp <- mcp(SpatialPoints(deer.locs), 100)# OLD VERSION
deer.mcp.xy <- deer.mcp@polygons[[1]]@Polygons[[1]]@coords
# FORTIYFY
require(ggplot2)
deer.mcp.xy <- fortify(deer.mcp)## Regions defined for each Polygonsplot(elevation)
points(deer.locs, asp=1, pch=19, col=rgb(0,0,.4,.5), cex=0.5)
lines(deer.mcp, col="darkred", lwd=2)
1.4 Create data frame of covariates
deer.sp<-SpatialPoints(deer.locs)
used.df <- data.frame(
used = rep(1, length(deer.sp)),
elevation = extract(elevation, deer.sp),
shrub = ceiling(extract(shrub, deer.sp)),
barren = ceiling(extract(barren, deer.sp)),
wells = ceiling(extract(wells, deer.sp)),
roads = ceiling(extract(roads, deer.sp)))
head(used.df)## used elevation shrub barren wells roads
## 1 1 2052.134 1 0 1 180
## 2 1 2050.978 1 0 1 130
## 3 1 2039.769 0 0 0 434
## 4 1 2044.926 1 0 0 308
## 5 1 2010.022 1 0 0 114
## 6 1 2023.509 1 0 0 71.5 Obtain some random points within mcp
random.points <- spsample(deer.mcp, 1000, "random")
1.6 Calculate the covariates for the null set
avail.df <- data.frame(
used = rep(0, length(random.points)),
elevation = extract(elevation, random.points),
shrub = ceiling(extract(shrub, random.points)),
barren = ceiling(extract(barren, random.points)),
wells = ceiling(extract(wells, random.points)),
roads = ceiling(extract(roads, random.points)))
head(avail.df)## used elevation shrub barren wells roads
## 1 0 2009.945 0 0 1 6
## 2 0 2016.945 1 0 1 46
## 3 0 2026.820 1 0 1 81
## 4 0 2055.548 1 0 1 228
## 5 0 2013.272 1 0 0 446
## 6 0 2025.959 1 0 0 1181.7 Combine…
deer.data <- rbind(avail.df, used.df)
with(deer.data, table(used, barren))## barren
## used 0 1
## 0 913 87
## 1 494 311.8 … and Model!
Simplest imaginable logistic regression model:
s <- function(x) x/sd(x)
deer.rsf <- glm(used ~ barren + shrub + s(elevation) + s(roads) + s(wells), family="binomial",
data = deer.data)## Estimate Std. Error z value P.value
## (Intercept) -8.18798873 3.75163161 -2.1825141 0.029
## barren -0.44857934 0.22538468 -1.9902832 0.047
## shrub -0.04090364 0.12841080 -0.3185374 0.750
## s(elevation) 0.12147586 0.06047401 2.0087282 0.045
## s(roads) 0.19206420 0.05579811 3.4421276 0.001
## s(wells) -0.27720721 0.06185451 -4.4816006 0.0001.9 Model Selection
For intense model selection …. Dredge!
#install.packages("MuMIn")
require(MuMIn)
deer.dredge <- dredge(deer.rsf)## Fixed term is "(Intercept)"options(na.action = "na.fail")
deer.rsf <- glm(used ~ shrub * barren * elevation, family="binomial", data = deer.data)
deer.dredge <- dredge(deer.rsf)## Fixed term is "(Intercept)"head(deer.dredge)## Global model call: glm(formula = used ~ shrub * barren * elevation, family = "binomial",
## data = deer.data)
## ---
## Model selection table
## (Int) brr elv shr brr:elv brr:shr elv:shr brr:elv:shr
## 39 -21.78 0.010410 27.60 -0.01360
## 40 -19.84 -0.2933 0.009469 25.98 -0.01282
## 56 -18.61 -0.5022 0.008874 24.51 0.4287 -0.01211
## 48 -21.81 14.6300 0.010440 27.56 -0.007452 -0.01360
## 64 -20.32 12.1700 0.009719 25.91 -0.006312 0.3651 -0.01280
## 128 -22.58 28.3300 0.010830 28.69 -0.014370 -36.5000 -0.01417 0.01843
## df logLik AICc delta weight
## 39 4 -975.307 1958.6 0.00 0.305
## 40 5 -974.436 1958.9 0.27 0.266
## 56 6 -973.984 1960.0 1.38 0.153
## 48 6 -973.999 1960.1 1.41 0.150
## 64 7 -973.685 1961.4 2.80 0.075
## 128 8 -973.049 1962.2 3.55 0.052
## Models ranked by AICc(x)Top model:
deer.rsf <- glm(used ~ barren + elevation + shrub +
barren*(elevation + shrub) + elevation:shrub,
family="binomial", data = deer.data)1.10 Predict from model
deer.rsf <- glm(used ~ barren + elevation + shrub + roads + wells,
family="binomial", data = deer.data)
coefs <- deer.rsf$coef
rsf<-exp(elevation*coefs["elevation"] + roads*coefs["roads"] + wells*coefs["wells"] +
shrub*coefs["shrub"] + barren*coefs["barren"])plot(rsf, main = "RSF"); p(alpha = 0.1)
2 Step Selection
A useful review:
The work flow:
The schematic:
Henrik Thurfjell et al. 2014, Movement Ecology
2.1 Generate null steps
Recalling the fuzzy catterpillar polar bear code (applied to the mule deer)…
require(scales)
load("deerrasters.robj")
Z <- (deer.locs[,1] + 1i*deer.locs[,2])[1:100]
# get pieces
n <- length(Z)
S <- Mod(diff(Z))
Phi <- Arg(diff(Z))
Theta <- diff(Phi)
RelSteps <- complex(mod = S[-1], arg=Theta)
# calculate null set
Z0 <- Z[-((n-1):n)]
Z1 <- Z[-c(1,n)]
Rotate <- complex(mod = 1, arg=Arg(Z1-Z0))
n.pseudo <- 3
palette(rainbow(length(Z)))
Z.null <- matrix(0, ncol=n.pseudo, nrow=n-2)
for(i in 1:length(Z1))
Z.null[i,] <- Z1[i] + sample(RelSteps,n.pseudo) * Rotate[i]
# plot
plot(Z, type="o", pch=19, asp=1, col=alpha(1:length(Z), .3))
for(i in 1:nrow(Z.null))
segments(rep(Re(Z1[i]), n-2), rep(Im(Z1[i]), n-2),
Re(Z.null[i,]), Im(Z.null[i,]), col=i)
2.2 Habitat “used”
Vital: to include “stratum” - just an indicator of the point against which you compare all null points!
deer.sp<-SpatialPoints(deer.locs)
used.df <- data.frame(
used = rep(1, length(deer.sp)),
elevation = extract(elevation, deer.sp),
shrub = ceiling(extract(shrub, deer.sp)),
barren = ceiling(extract(barren, deer.sp)),
wells = ceiling(extract(wells, deer.sp)),
roads = ceiling(extract(roads, deer.sp)),
stratum = 1:length(deer.sp))
head(used.df)## used elevation shrub barren wells roads stratum
## 1 1 2052.134 1 0 1 180 1
## 2 1 2050.978 1 0 1 130 2
## 3 1 2039.769 0 0 0 434 3
## 4 1 2044.926 1 0 0 308 4
## 5 1 2010.022 1 0 0 114 5
## 6 1 2023.509 1 0 0 7 62.3 Habitat available
Reshape the “Z.null” and add stratum (code hidden here)
stratumStack <- function(x){
df <- data.frame(stratum = rep(1:nrow(x), ncol(x)), as.vector(x))
names(df)[2] <- "Z"
return(df)
}
avail.df <- stratumStack(Z.null) %>% mutate(X = Re(Z), Y = Im(Z))
avail.points <- SpatialPoints(avail.df[,c("X","Y")])
avail.df <- data.frame(
used = rep(0, length(avail.points)),
elevation = extract(elevation, avail.points),
shrub = ceiling(extract(shrub, avail.points)),
barren = ceiling(extract(barren, avail.points)),
wells = ceiling(extract(wells, avail.points)),
roads = ceiling(extract(roads, avail.points)),
stratum = avail.df$stratum)## Warning in .doExtract(x, i, ..., drop = drop): some indices are invalid (NA
## returned)
## Warning in .doExtract(x, i, ..., drop = drop): some indices are invalid (NA
## returned)
## Warning in .doExtract(x, i, ..., drop = drop): some indices are invalid (NA
## returned)
## Warning in .doExtract(x, i, ..., drop = drop): some indices are invalid (NA
## returned)
## Warning in .doExtract(x, i, ..., drop = drop): some indices are invalid (NA
## returned)head(avail.df)## used elevation shrub barren wells roads stratum
## 1 0 2043.044 1 0 1 26 1
## 2 0 2030.998 1 0 0 518 2
## 3 0 2036.394 1 0 0 298 3
## 4 0 2012.110 1 0 1 320 4
## 5 0 1996.990 1 0 0 261 5
## 6 0 2048.271 0 0 0 137 62.4 Conditional Regression Model
\[{\rm Pr}\Bigl[(y_{i1},\ldots,y_{i{n_i}})|\sum_{j=1}^{n_i} y_{ij} = k_{1i}\Bigr] = \frac{\exp(\sum_{j=1}^{n_i} y_{ij} x_{ij}'\beta)}{\sum_{{\bf d}_i\in S_i}\exp(\sum_{j=1}^{n_i} y_{ij} x_{ij}'\beta)}\]
deer.ssf.model <- clogit(used~elevation+shrub+barren+wells+roads +
+ strata(stratum), data = deer.ssf.data,
na.action = "na.fail")
summary(deer.ssf.model)## Call:
## coxph(formula = Surv(rep(1, 817L), used) ~ elevation + shrub +
## barren + wells + roads + +strata(stratum), data = deer.ssf.data,
## na.action = "na.fail", method = "exact")
##
## n= 817, number of events= 525
##
## coef exp(coef) se(coef) z Pr(>|z|)
## elevation 0.0133872 1.0134772 0.0074636 1.794 0.0729 .
## shrub 0.4205956 1.5228683 0.3230715 1.302 0.1930
## barren -0.9855931 0.3732178 0.8110667 -1.215 0.2243
## wells -0.5545708 0.5743187 0.3179635 -1.744 0.0811 .
## roads 0.0002976 1.0002976 0.0008469 0.351 0.7253
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## elevation 1.0135 0.9867 0.99876 1.028
## shrub 1.5229 0.6567 0.80847 2.869
## barren 0.3732 2.6794 0.07613 1.830
## wells 0.5743 1.7412 0.30797 1.071
## roads 1.0003 0.9997 0.99864 1.002
##
## Rsquare= 0.013 (max possible= 0.282 )
## Likelihood ratio test= 10.79 on 5 df, p=0.05571
## Wald test = 9.41 on 5 df, p=0.09388
## Score (logrank) test = 10.13 on 5 df, p=0.07162.5 SSF to think about: Scale!
2.6 SSF to think about: Edges!
2.7 SSF to think about: Path vs. Points!
