Fisher SSF

Preamble

Load libraries

library(ezknitr)
library(knitr)
library(lubridate)
library(raster)
library(move)
library(amt) 
library(tidyverse)
library(here)
options(width=150)
opts_chunk$set(fig.width=12,fig.height=4.5)

Read in the fisher tracks

ssfdat <- read_rds("data/trk.Rdata")

Prepare environmental data

Add environmental covariates, this is very similar to the first part of the RSF script. Again, we will use three covariates: land use class, elevation and population density

landuse <- raster(here("data/landuse/landuse_study_area.tif"))
elevation <- raster(here("data/elevation/ASTER ASTGTM2 Elevation-20100101120000000-0-0.tif"))
popden <- raster(here("data/pop_den/pop_den.tif"))

First, we need to bring all thre layers to the same CRS.

get_crs(ssfdat)

## CRS arguments:
##  +init=epsg:5070 +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80
## +towgs84=0,0,0

landuse <- projectRaster(landuse, crs = get_crs(ssfdat), method = "ngb")
elevation <- projectRaster(elevation, crs = get_crs(ssfdat))
popden <- projectRaster(popden, crs = get_crs(ssfdat))

Plot the raster

plot(landuse)

plot(elevation)

plot(popden)

Check the resolution

res(landuse)

## [1] 30 30

res(elevation)

## [1] 44.0 60.4

res(popden)

## [1] 659 905

Assign each raster a meaningful name

names(landuse) <- "landclass"
names(elevation) <- "ele"
names(popden) <- "popden"

Resolutions (and also extents) are different. This means, that we can not stack the rasters. There two options:

Resample rasters to the coarsest raster (popden) using the function raster::resample.
Extract covariate values from each raster idendependetly. We will continue with second option. ## Create steps and annotate with environmental data sampling rate

summarize_sampling_rate_many(ssfdat, id)

## # A tibble: 8 x 10
##   id        min    q1 median  mean    q3   max    sd     n unit 
##   <fct>   <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 M1    13.3    14.7   15.0  32.7  15.4   977. 104.    918 min  
## 2 M4     0.1000  1.93   2.03  8.04  2.57 1209.  44.0  8957 min  
## 3 F2     0.4     1.97   2.03  8.93  3.98 1080.  48.1  3003 min  
## 4 M3     0.417   1.97   2.17 13.3   5.42 2889.  90.2  2435 min  
## 5 F3     0.417   1.97   2.07 15.7   4.08 2167. 104.   1500 min  
## 6 M2     1.35    9.78  10.0  21.7  10.3  2111.  87.2  1637 min  
## 7 F1     3.92    9.80  10.0  20.7  10.3  1650.  95.2  1348 min  
## 8 M5     0.250   1.93   2.02  7.41  2.37  976.  25.0 13098 min

The approach is as follows:

Resample track and filter bursts
Convert track to steps
Create random steps
Extract covariate values

For M4

ssfdat %>% filter(id == "M4") %>% 
  track_resample(rate = minutes(10), tolerance = minutes(1)) %>% 
  filter_min_n_burst() %>% 
  steps_by_burst() %>% 
  random_steps() %>% 
  extract_covariates(landuse) %>%  
  extract_covariates(elevation) %>% 
  extract_covariates(popden)

## Warning in random_steps.steps_xy(.): Step-lengths or turning angles contained NA, which were removed.

## # A tibble: 13,024 x 15
##    burst_ step_id_ case_      x1_      y1_      x2_      y2_ t1_                 t2_                 dt_             sl_    ta_ landclass   ele popden
##  *  <dbl>    <int> <lgl>    <dbl>    <dbl>    <dbl>    <dbl> <dttm>              <dttm>              <time>        <dbl>  <dbl>     <dbl> <dbl>  <dbl>
##  1     16        1 TRUE  1780469. 2412224. 1780468. 2412230. 2010-02-10 03:10:39 2010-02-10 03:20:10 9.500017 mi~   6.01  2.96         90  112.   309.
##  2     16        1 FALSE 1780469. 2412224. 1780496. 2412005. 2010-02-10 03:10:39 2010-02-10 03:20:10 9.500017 mi~ 221.   -0.183        21  111.   309.
##  3     16        1 FALSE 1780469. 2412224. 1780454. 2412189. 2010-02-10 03:10:39 2010-02-10 03:20:10 9.500017 mi~  37.9  -0.718        90  111.   309.
##  4     16        1 FALSE 1780469. 2412224. 1780505. 2412369. 2010-02-10 03:10:39 2010-02-10 03:20:10 9.500017 mi~ 149.    2.60         90  114.   309.
##  5     16        1 FALSE 1780469. 2412224. 1780484. 2412245. 2010-02-10 03:10:39 2010-02-10 03:20:10 9.500017 mi~  25.7   2.22         90  112.   309.
##  6     16        1 FALSE 1780469. 2412224. 1780474. 2412232. 2010-02-10 03:10:39 2010-02-10 03:20:10 9.500017 mi~   9.13  2.29         90  112.   309.
##  7     16        1 FALSE 1780469. 2412224. 1780524. 2412257. 2010-02-10 03:10:39 2010-02-10 03:20:10 9.500017 mi~  64.3   1.81         90  114.   309.
##  8     16        1 FALSE 1780469. 2412224. 1780441. 2412018. 2010-02-10 03:10:39 2010-02-10 03:20:10 9.500017 mi~ 208.   -0.439        21  113.   309.
##  9     16        1 FALSE 1780469. 2412224. 1780340. 2412185. 2010-02-10 03:10:39 2010-02-10 03:20:10 9.500017 mi~ 135.   -1.58         90  118.   309.
## 10     16        1 FALSE 1780469. 2412224. 1780500. 2412066. 2010-02-10 03:10:39 2010-02-10 03:20:10 9.500017 mi~ 161.   -0.108        21  111.   309.
## # ... with 13,014 more rows

Because SSF coefficients depend on the time interval between observations, it is advantageous for all individuals to have the same sampling rate. #’ Here, all individuasl except M1 had 10 minute or 2 minute sampling rates. Lets resample all individuals except M1 to a 10 minute sampling rate.

ssfdat <- ssfdat %>% filter(id != "M1") %>% group_by(id) %>% nest() %>% 
  mutate(data = map(
    data, ~ .x %>% 
      track_resample(rate = minutes(10), tolerance = minutes(1)) %>%
      filter_min_n_burst() %>% 
      steps_by_burst() %>% 
      random_steps() %>% 
      extract_covariates(landuse) %>%  
      extract_covariates(elevation) %>% 
      extract_covariates(popden))) %>% unnest()

## Warning in random_steps.steps_xy(.): Step-lengths or turning angles contained NA, which were removed.

## Warning in random_steps.steps_xy(.): Step-lengths or turning angles contained NA, which were removed.

## Warning in random_steps.steps_xy(.): Step-lengths or turning angles contained NA, which were removed.

## Warning in random_steps.steps_xy(.): Step-lengths or turning angles contained NA, which were removed.

## Warning in random_steps.steps_xy(.): Step-lengths or turning angles contained NA, which were removed.

## Warning in random_steps.steps_xy(.): Step-lengths or turning angles contained NA, which were removed.

## Warning in random_steps.steps_xy(.): Step-lengths or turning angles contained NA, which were removed.

Create landcover classes (as suggested by Scott Lapoint :)

ssfdat <- ssfdat %>%  mutate(
  landclass = as.character(landclass), 
  landC = fct_collapse(landclass,
      agri = c("81", "82"),
      forest =c("41", "42", "43"),
      shrub = c("52"),
      grass = c("31", "71"),
      wet = c("90", "95"),
      other = c("11", "21", "22", "23", "24")))

Center and scale variables, make response numeric

ssfdat <- ssfdat %>% mutate(elev = scale(ele), 
                            popD = scale(popden), 
                            case_ = as.numeric(case_))

Save SSF data for later (multiple animal)

write_rds(ssfdat, "data/ssf_dat.rds")

Explore the data

Look Distribution of habitat class for used and available data

ggplot(ssfdat, aes(x = landC, y = ..prop.., group = case_, colour = case_))+
  geom_bar(position = "dodge", aes(fill = case_)) +
  facet_wrap( ~ id, scales = "free")

SSF model fitting

Now, fit an SSF model to data from each animal. Since not all animals experience all habitat types, lets just explore forest versus non-forest.

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

Fit an SSF to a single animal

summary(fit_issf(case_ ~ elev + popD + forest + sl_ + log(sl_) + strata(step_id_), 
                 data = filter(ssfdat, id == "M4")))

## Call:
## coxph(formula = Surv(rep(1, 13024L), case_) ~ elev + popD + forest + 
##     sl_ + log(sl_) + strata(step_id_), data = data, method = "exact")
## 
##   n= 13024, number of events= 1184 
## 
##                coef  exp(coef)   se(coef)      z Pr(>|z|)    
## elev      0.7444981  2.1053846  0.6208104  1.199   0.2304    
## popD     -0.1166513  0.8898955  0.0674116 -1.730   0.0836 .  
## forest    0.4554875  1.5769419  0.0735100  6.196 5.78e-10 ***
## sl_       0.0002094  1.0002094  0.0004180  0.501   0.6165    
## log(sl_) -0.0063753  0.9936450  0.0374310 -0.170   0.8648    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## elev        2.1054     0.4750    0.6236     7.108
## popD        0.8899     1.1237    0.7798     1.016
## forest      1.5769     0.6341    1.3653     1.821
## sl_         1.0002     0.9998    0.9994     1.001
## log(sl_)    0.9936     1.0064    0.9234     1.069
## 
## Concordance= 0.551  (se = 0.01 )
## Rsquare= 0.004   (max possible= 0.353 )
## Likelihood ratio test= 46.12  on 5 df,   p=9e-09
## Wald test            = 45.44  on 5 df,   p=1e-08
## Score (logrank) test = 45.82  on 5 df,   p=1e-08

Fit an SSF model to data from each animal

fitted_ssf <- function(data){
    fit_issf(case_ ~ elev + popD + forest + sl_ + log(sl_) + strata(step_id_), data=data)
}
ssffits <-ssfdat %>%  group_by(id) %>% nest() %>% 
  mutate(mod = map(data, fitted_ssf))

Look at first model

ssffits$mod[[1]]

## $model
## Call:
## survival::clogit(formula, data = data, ...)
## 
##                coef  exp(coef)   se(coef)      z        p
## elev      0.7444981  2.1053846  0.6208104  1.199   0.2304
## popD     -0.1166513  0.8898955  0.0674116 -1.730   0.0836
## forest    0.4554875  1.5769419  0.0735100  6.196 5.78e-10
## sl_       0.0002094  1.0002094  0.0004180  0.501   0.6165
## log(sl_) -0.0063753  0.9936450  0.0374310 -0.170   0.8648
## 
## Likelihood ratio test=46.12  on 5 df, p=8.581e-09
## n= 13024, number of events= 1184 
## 
## $sl_
## NULL
## 
## $ta_
## NULL
## 
## $more
## NULL
## 
## attr(,"class")
## [1] "fit_clogit" "list"

Now, use tidy to extract information about the model fits

ssffits <- ssffits %>%
  mutate(tidy = map(mod, ~ broom::tidy(.x$model)),
         n = map_int(data, nrow))

ssffits$tidy

## [[1]]
## # A tibble: 5 x 7
##   term      estimate std.error statistic  p.value  conf.low conf.high
##   <chr>        <dbl>     <dbl>     <dbl>    <dbl>     <dbl>     <dbl>
## 1 elev      0.744     0.621        1.20  2.30e- 1 -0.472      1.96   
## 2 popD     -0.117     0.0674      -1.73  8.36e- 2 -0.249      0.0155 
## 3 forest    0.455     0.0735       6.20  5.78e-10  0.311      0.600  
## 4 sl_       0.000209  0.000418     0.501 6.16e- 1 -0.000610   0.00103
## 5 log(sl_) -0.00638   0.0374      -0.170 8.65e- 1 -0.0797     0.0670 
## 
## [[2]]
## # A tibble: 5 x 7
##   term       estimate std.error statistic  p.value conf.low conf.high
##   <chr>         <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 elev      5.87       1.64        3.58   0.000340  2.66      9.08   
## 2 popD     -0.0136     0.129      -0.106  0.916    -0.266     0.239  
## 3 forest    0.353      0.145       2.43   0.0151    0.0683    0.638  
## 4 sl_      -0.0000190  0.000582   -0.0327 0.974    -0.00116   0.00112
## 5 log(sl_) -0.00176    0.0732     -0.0240 0.981    -0.145     0.142  
## 
## [[3]]
## # A tibble: 5 x 7
##   term      estimate std.error statistic p.value conf.low conf.high
##   <chr>        <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 elev      3.65      3.12         1.17  0.242   -2.46      9.76   
## 2 popD     -0.120     0.182       -0.663 0.507   -0.477     0.236  
## 3 forest    0.501     0.192        2.60  0.00921  0.124     0.878  
## 4 sl_       0.000153  0.000780     0.196 0.845   -0.00138   0.00168
## 5 log(sl_)  0.0337    0.0948       0.356 0.722   -0.152     0.219  
## 
## [[4]]
## # A tibble: 5 x 7
##   term      estimate std.error statistic p.value conf.low conf.high
##   <chr>        <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 elev      8.60       4.57       1.88    0.0601 -0.366    17.6    
## 2 popD     -0.571      0.367     -1.55    0.120  -1.29      0.149  
## 3 forest    0.282      0.252      1.12    0.264  -0.212     0.777  
## 4 sl_       0.000138   0.00138    0.0997  0.921  -0.00257   0.00284
## 5 log(sl_)  0.0123     0.103      0.120   0.905  -0.189     0.213  
## 
## [[5]]
## # A tibble: 5 x 7
##   term      estimate std.error statistic   p.value  conf.low conf.high
##   <chr>        <dbl>     <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
## 1 elev      2.20      0.837       2.63   0.00843    0.564     3.84    
## 2 popD     -0.0265    0.105      -0.253  0.800     -0.232     0.179   
## 3 forest    0.389     0.0903      4.30   0.0000168  0.212     0.565   
## 4 sl_       0.000175  0.000346    0.506  0.613     -0.000503  0.000853
## 5 log(sl_) -0.000959  0.0221     -0.0433 0.965     -0.0443    0.0424  
## 
## [[6]]
## # A tibble: 5 x 7
##   term      estimate std.error statistic  p.value  conf.low conf.high
##   <chr>        <dbl>     <dbl>     <dbl>    <dbl>     <dbl>     <dbl>
## 1 elev      0.169     0.716        0.236 8.14e- 1 -1.23       1.57   
## 2 popD      0.274     0.144        1.89  5.82e- 2 -0.00953    0.557  
## 3 forest    0.881     0.106        8.28  1.26e-16  0.673      1.09   
## 4 sl_       0.000800  0.000472     1.70  9.00e- 2 -0.000125   0.00172
## 5 log(sl_) -0.00959   0.0251      -0.382 7.02e- 1 -0.0588     0.0396 
## 
## [[7]]
## # A tibble: 5 x 7
##   term       estimate std.error statistic  p.value  conf.low conf.high
##   <chr>         <dbl>     <dbl>     <dbl>    <dbl>     <dbl>     <dbl>
## 1 elev      0.0525     0.189        0.277 0.782    -0.318     0.423   
## 2 popD      2.92       0.792        3.69  0.000227  1.37      4.47    
## 3 forest   -0.334      0.0911      -3.67  0.000246 -0.513    -0.155   
## 4 sl_      -0.0000754  0.000311    -0.243 0.808    -0.000684  0.000534
## 5 log(sl_)  0.00862    0.0299       0.289 0.773    -0.0499    0.0671

Now, create data frame w/ the coefficients, etc

ssf_coefs <- ssffits %>%
  unnest(tidy) %>%
  select(-(std.error:conf.high)) 
  
ssf_coefs %>% spread(term, estimate)

## # A tibble: 7 x 7
##   id        n   elev forest `log(sl_)`    popD        sl_
##   <fct> <int>  <dbl>  <dbl>      <dbl>   <dbl>      <dbl>
## 1 F1    12848 0.169   0.881  -0.00959   0.274   0.000800 
## 2 F2     3817 5.87    0.353  -0.00176  -0.0136 -0.0000190
## 3 F3     1815 8.60    0.282   0.0123   -0.571   0.000138 
## 4 M2    14113 2.20    0.389  -0.000959 -0.0265  0.000175 
## 5 M3     2233 3.65    0.501   0.0337   -0.120   0.000153 
## 6 M4    13024 0.744   0.455  -0.00638  -0.117   0.000209 
## 7 M5    21681 0.0525 -0.334   0.00862   2.92   -0.0000754

Plot coefficients

ssf_coefs %>% 
  ggplot(., aes(x=1, y=estimate)) + 
  geom_dotplot(binaxis="y", stackdir="center")+geom_hline(yintercept=0)+
  facet_wrap(~term, scales="free")

## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

Write out coefficients for MultipleAnimals.R

save(ssf_coefs, file="data/ssfcoefs.Rdata")

Document Footer

Document spun with: ezspin(“FisherSSF.R”, fig_dir = “figures”, keep_md=FALSE)

Session Information:

sessionInfo()

## R version 3.5.3 (2019-03-11)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 16299)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252    LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                           LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] here_0.1        forcats_0.4.0   stringr_1.4.0   dplyr_0.8.0.1   purrr_0.3.2     readr_1.3.1     tidyr_0.8.3     tibble_2.0.1    ggplot2_3.1.0  
## [10] tidyverse_1.2.1 amt_0.0.6       move_3.1.0      rgdal_1.4-2     geosphere_1.5-7 raster_2.8-19   sp_1.3-1        lubridate_1.7.4 knitr_1.22     
## [19] ezknitr_0.6    
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.1          mvtnorm_1.0-10      lattice_0.20-38     utf8_1.1.4          assertthat_0.2.1    rprojroot_1.3-2     digest_0.6.18      
##  [8] R6_2.4.0            cellranger_1.1.0    plyr_1.8.4          backports_1.1.4     evaluate_0.13       httr_1.4.0          highr_0.7          
## [15] pillar_1.3.1        rlang_0.3.1         lazyeval_0.2.2      readxl_1.3.0        rstudioapi_0.9.0    Matrix_1.2-17       rmarkdown_1.11     
## [22] labeling_0.3        splines_3.5.3       munsell_0.5.0       broom_0.5.1         compiler_3.5.3      modelr_0.1.4        xfun_0.5           
## [29] pkgconfig_2.0.2     htmltools_0.3.6     tidyselect_0.2.5    fitdistrplus_1.0-14 codetools_0.2-16    fansi_0.4.0         crayon_1.3.4       
## [36] withr_2.1.2         MASS_7.3-51.1       grid_3.5.3          nlme_3.1-137        jsonlite_1.6        gtable_0.2.0        npsurv_0.4-0       
## [43] magrittr_1.5        scales_1.0.0        lsei_1.2-0          cli_1.0.1           stringi_1.1.7       xml2_1.2.0          generics_0.0.2     
## [50] boot_1.3-20         tools_3.5.3         glue_1.3.1          hms_0.4.2           parallel_3.5.3      survival_2.43-3     yaml_2.2.0         
## [57] colorspace_1.4-0    circular_0.4-93     rvest_0.3.2         memoise_1.1.0       haven_2.1.0