Purpose: demonstrate methods for simulating space use from fitted RSF and SSFs.

Preamble

Load libraries

library(knitr)
library(lubridate)
library(raster)
library(amt) 
library(broom)
library(tidyverse)
library(tictoc)
library(here)
library(glmmTMB)

options(width=165)
opts_chunk$set(fig.width=12,fig.height=4.5, error=TRUE,cache = FALSE)

Record time for running all code

tic("total")

RSF models

Read in annotated available data for RSF modeling

rsfdat <- read_rds("data/rsf_dat.rds")

Forest versus non-forest

rsfdat$forest <- ifelse(rsfdat$landC == "forest", 1, 0)

Weight available data

rsfdat$w <- ifelse(rsfdat$case_ == 1, 1, 5e3)

We will use only one animal to illustrate simulations. Simulations could be extendend to multiple animals, but not all the details are worked out/implemented yet.

Fit an RSF for F1

dat <- rsfdat %>% filter(id == "F1")
rsffit <- dat %>% mutate(elev = elev / 1e3) %>% 
  glm(case_ ~ forest + elev, data = ., weight = w, family = binomial())

summary(rsffit)
## 
## Call:
## glm(formula = case_ ~ forest + elev, family = binomial(), data = ., 
##     weights = w)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.4408  -0.3886  -0.3622  -0.1858   5.0715  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -11.79414    0.25173 -46.852  < 2e-16 ***
## forest        1.46158    0.07646  19.116  < 2e-16 ***
## elev        937.24059  291.71346   3.213  0.00131 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 33760  on 28325  degrees of freedom
## Residual deviance: 33249  on 28323  degrees of freedom
## AIC: 33255
## 
## Number of Fisher Scoring iterations: 7

Next we prepare the environmental covariates: forest and elevation.

landuse <- raster("data/landuse/landuse_study_area.tif") 
landuse <- projectRaster(landuse, crs = sp::CRS("+init=epsg:5070"), method = "ngb")
forest <- landuse %in% 41:43

elev <- raster("data/elevation/ASTER ASTGTM2 Elevation-20100101120000000-0-0.tif")

We have to resample the elevation raster, so that it matches the forests.

elev <- projectRaster(elev, crs = sp::CRS("+init=epsg:5070"))
elev <- resample(elev, landuse)
elev <- elev / 1e3

Next, we crop the raster to the 1 km buffer of a bounding box of the observed data.

forest_small <- crop(forest, extent(c(range(dat$x_) + c(-1e3, 1e3), 
                                       range(dat$y_) + c(-1e3, 1e3))))

elev_small <- crop(elev, extent(c(range(dat$x_) + c(-1e3, 1e3), 
                                    range(dat$y_) + c(-1e3, 1e3))))

plot(forest_small)

plot(elev_small)

hull_polygon <- make_track(dat, x_, y_) %>% hr_mcp(level = 1) %>% hr_isopleths()
## .t missing, creating `track_xy`.
# mask the elev and forest rasters
elev_mcp <- mask(elev_small, hull_polygon)
forest_mcp<-mask(forest_small, hull_polygon)

Now we can calcualte UD from the RSF. Note, we are ignoring the intercept since it cancels out.

ud_lp <- coef(rsffit)[2] * forest_small + 
  coef(rsffit)[3] * elev_small 
ud <- exp(ud_lp)
ud <- ud / sum(ud[], na.rm = TRUE)
plot(ud)
with(filter(dat, case_ == 1), points(x_, y_, pch="."))

Better to plot the log of the UD.

plot(log(ud))
with(filter(dat, case_ == 1), points(x_, y_, pch="."))

Now, if limit just to mcp

ud_lp <- coef(rsffit)[2] * forest_mcp + 
  coef(rsffit)[3] * elev_mcp 
ud <- exp(ud_lp)
ud <- ud / sum(ud[], na.rm = TRUE)
plot(ud)
with(filter(dat, case_ == 1), points(x_, y_, pch="."))

Better to plot the log of the UD.

plot(log(ud))
with(filter(dat, case_ == 1), points(x_, y_, pch="."))

Simulate UDs from fitted SSF models

Read in annotated available data for SSF modeling

ssfdat <- read_rds("data/ssf_dat.rds")

Forest versus non-forest

ssfdat$forest <- ifelse(ssfdat$landC == "forest", 1, 0)

We will use only one animal to illustrate simulations. Simulations could be extendend to multiple animals, but not all the details are worked out/implemented yet.

dat <- ssfdat %>% filter(id == "F1")
class(dat) <- c("stepx_xyt", "steps_xy", class(dat)) # class is lost due to the nesting. 
ssffit <- dat %>% 
  fit_issf(case_ ~ forest + strata(step_id_))

summary(ssffit)
## Call:
## coxph(formula = Surv(rep(1, 12848L), case_) ~ forest + strata(step_id_), 
##     data = data, method = "exact")
## 
##   n= 12848, number of events= 1168 
## 
##          coef exp(coef) se(coef)     z Pr(>|z|)    
## forest 0.8496    2.3388   0.1019 8.335   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##        exp(coef) exp(-coef) lower .95 upper .95
## forest     2.339     0.4276     1.915     2.856
## 
## Concordance= 0.545  (se = 0.005 )
## Rsquare= 0.006   (max possible= 0.353 )
## Likelihood ratio test= 77.71  on 1 df,   p=<2e-16
## Wald test            = 69.47  on 1 df,   p=<2e-16
## Score (logrank) test = 71.29  on 1 df,   p=<2e-16
landuse <- raster("data/landuse/landuse_study_area.tif") 
landuse <- projectRaster(landuse, crs = sp::CRS("+init=epsg:5070"), method = "ngb")
forest <- landuse %in% 41:43
names(forest) <- "forest"

# get movement characteristics
sl_char <- fit_sl_dist_base(filter(dat, case_ == 1) %>% pull(sl_))
sl_char$fit$estimate["shape"]
##     shape 
## 0.6495437
mk <- movement_kernel(
  shape = sl_char$fit$estimate["shape"], 
  scale = 1 / sl_char$fit$estimate["rate"], 
  template = forest
)
plot(mk)

hk <- habitat_kernel(
  coef = coef(ssffit), 
  resources = forest
)
  
plot(hk)

We want to restrict the area again to the area where the animal actually occured.

forest_small <- crop(forest, extent(c(range(dat$x1_) + c(-5e3, 5e3), 
                                      range(dat$y1_) + c(-5e3, 5e3))))
plot(forest_small)

hk_small <- habitat_kernel(
  coef = coef(ssffit), 
  resources = forest_small
)

Now lets simulate an UD

ud <- simulate_ud(mk, hk_small, start = as.numeric(dat[1, c("x1_", "y1_")]), 
                   n = 5e5)
plot(ud)

It is probably more realistic to add home ranging behaviour. We can do this by including x and y and x^2 and y^2 in the model.

y <- x <- hk_small
x[] <- xFromCell(hk_small, 1:ncell(hk_small))
y[] <- yFromCell(hk_small, 1:ncell(hk_small))

x <- (x - mean(x[])) / 1e3
y <- (y - mean(y[])) / 1e3

plot(x)

plot(y)

xy <- stack(x, y, x^2, y^2)
names(xy) <- c("x", "y", "x2", "y2")
plot(xy)

dat <- dat %>% extract_covariates(xy)

ssffit1 <- dat %>% 
  fit_issf(case_ ~ forest + x + y + x2 + y2 + strata(step_id_))

summary(ssffit1)
## Call:
## coxph(formula = Surv(rep(1, 12848L), case_) ~ forest + x + y + 
##     x2 + y2 + strata(step_id_), data = data, method = "exact")
## 
##   n= 12848, number of events= 1168 
## 
##            coef exp(coef) se(coef)      z Pr(>|z|)    
## forest  0.84014   2.31670  0.10321  8.140 3.96e-16 ***
## x      -0.24257   0.78461  0.40323 -0.602    0.547    
## y      -0.45747   0.63289  0.40390 -1.133    0.257    
## x2     -0.17047   0.84327  0.28068 -0.607    0.544    
## y2     -0.06656   0.93561  0.28724 -0.232    0.817    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##        exp(coef) exp(-coef) lower .95 upper .95
## forest    2.3167     0.4316    1.8924     2.836
## x         0.7846     1.2745    0.3560     1.729
## y         0.6329     1.5801    0.2868     1.397
## x2        0.8433     1.1859    0.4865     1.462
## y2        0.9356     1.0688    0.5328     1.643
## 
## Concordance= 0.545  (se = 0.01 )
## Rsquare= 0.006   (max possible= 0.353 )
## Likelihood ratio test= 79.57  on 5 df,   p=1e-15
## Wald test            = 71.01  on 5 df,   p=6e-14
## Score (logrank) test = 72.76  on 5 df,   p=3e-14
hk_small <- habitat_kernel(
  coef = coef(ssffit1), 
  resources = stack(forest_small, xy)
)

ud1 <- simulate_ud(mk, hk_small, start = as.numeric(dat[1, c("x1_", "y1_")]), 
                   n = 1e6)

plot(ud1)

Probably it is a good idea to increase n to 1e7 (this will take a bit longer for the simulations to run)

ud2 <- simulate_ud(mk, hk_small, start = as.numeric(dat[1, c("x1_", "y1_")]), 
                   n = 1e7)

plot(ud2)
with(filter(dat, case_ == 1), points(x1_, y1_, pch="."))