The Fisher movement data was obtained from LaPoint, S., Gallery, P., Wikelski, M. and Kays, R. (2013). Animal behavior, cost-based corridor models, and real corridors. Landscape Ecology, 28, 1615-1630.
library(tidyverse)
library(amt)
library(copula)
library(cylcop) #Some functions are only implemented in the development version.
#Available at https://github.com/floo66/cylcop
dat <- read_csv("Fisher.csv") %>%
filter(id %in% c(1465, 1466, 1072, 1078, 1016, 1469))
dat_all <- dat %>% nest(data=c(x,y,t))
dat_all$sex <- c("f", "f", "f", "m", "m", "m")We then create a track from the data and convert it to an accurate projection for that geographic region with unit meters (epsg:5070).
dat_all <- dat_all %>%
mutate(trk = map(data, function(d) {
make_track(d, x, y, t, crs = sp::CRS("+init=epsg:4326")) %>%
transform_coords(sp::CRS("+init=epsg:5070"))
}))Next, we show information on the sampling rate of the different individuals,
dat_all %>% mutate(sr = lapply(trk, summarize_sampling_rate)) %>%
select(id, sr) %>% unnest(cols=c(sr))## # A tibble: 6 x 10
## id min q1 median mean q3 max sd n unit
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 1072 3.92 9.80 10.0 20.7 10.3 1650. 95.2 1348 min
## 2 1465 0.4 1.97 2.03 8.93 3.98 1080. 48.1 3003 min
## 3 1466 0.417 1.97 2.07 15.7 4.08 2167. 104. 1500 min
## 4 1078 1.35 9.78 10.0 21.7 10.3 2111. 87.2 1637 min
## 5 1469 0.417 1.97 2.17 13.3 5.42 2889. 90.2 2435 min
## 6 1016 0.100 1.93 2.03 8.04 2.57 1209. 44.0 8957 minresample tracks at a rate based on the sampling rates above,
dat_all <- dat_all %>% mutate(dat_clean = map(trk, ~ {
.x %>% track_resample(rate = minutes(10), tolerance = seconds(60))
}))and convert to steps.
dat_all <-
dat_all %>%
mutate(stps = map(dat_clean, ~ .x %>% steps_by_burst()))For fitting, we will pool together the step lengths and turn angles of all animals. The step length distributions are fitted to all available step lengths, not only to those where the corresponding turning angle is not NA.
dat_pool <- dat_all %>% select(stps) %>% unnest(cols=c(stps))
x <- dat_pool$sl_
theta <- dat_pool$ta_[!is.na(dat_pool$ta_)]We can now fit distributions to the step lengths.
cylcop::fit_steplength(x, parametric = "gamma")
cylcop::fit_steplength(x, parametric = "weibull")
cylcop::fit_steplength(x, parametric = "lognorm")
cylcop::fit_steplength(x, parametric = "weibullmix",ncomp = 2)
cylcop::fit_steplength(x, parametric = "lnormmix",ncomp = 2)
marg_dist <- cylcop::fit_steplength(x, parametric = "gammamix",ncomp = 2)The lowest AIC is achieved with a log-normal distribution.
marg_dist## $coef
## $coef$prop
## [1] 0.4116987 0.5883013
##
## $coef$shape
## [1] 1.698242 1.473340
##
## $coef$scale
## [1] 10.60563 110.62237
##
##
## $se
## list()
##
## $logL
## [1] -30583.77
##
## $AIC
## [1] 61177.54
##
## $name
## [1] "gammamix"Next, we fit circular distributions to the turn angles.
cylcop::fit_angle(theta, parametric = "vonmisesmix", ncomp = 2)
cylcop::fit_angle(theta, parametric = "vonmises")
cylcop::fit_angle(theta, parametric = "wrappedcauchy")
marg_angle <- cylcop::fit_angle(theta,
parametric = "vonmisesmix",
mu = c(0, pi),
ncomp = 2)The lowest AIC is achieved with a mixed von Mises distribution with fixed location parameters.
marg_angle## $coef
## $coef$mu
## [1] 0.000000 3.141593
##
## $coef$kappa
## [1] 0.4839508 6.3070649
##
## $coef$prop
## [1] 0.7757562 0.2242438
##
##
## $se
## list()
##
## $logL
## [1] -7791.917
##
## $AIC
## [1] 15589.83
##
## $name
## [1] "vonmisesmix"First, we remove step lengths that correspond to NA in turn angles from the data.
x <- dat_pool$sl_[!is.na(dat_pool$ta_)]Next, we investigate circular-linear correlation measures of the bivariate data.
cylcop::cor_cyl(theta = theta, x = x)## [1] 0.1024313cylcop::mi_cyl(
theta = theta,
x = x,
normalize = T,
symmetrize = T
)## [1] 0.03732642After this, we roughly estimate copula parameters using a grid search.
start_quadsec <-
cylcop::optCor(
copula = cyl_quadsec(),
theta = theta,
x = x,
method = "cor_cyl"
)
start_quadsec <-
cylcop::optCor(
copula = cyl_quadsec(),
theta = theta,
x = x,
method = "mi_cyl"
)
start_cubsec <-
cylcop::optCor(
copula = cyl_cubsec(),
theta = theta,
x = x,
method = "cor_cyl",
acc = 0.05,
n = 10000
)For cyl_rect_combine copulae with symmetric rectangles spanning the entire unit square, we can calculate parameters from Kendall’s tau.
start_rect_frank <-
cylcop::optCor(
copula = cyl_rect_combine(
copula = frankCopula(),
low_rect = c(0, 0.5),
up_rect = "symmetric"
),
theta = theta,
x = x,
method = "tau"
)
start_rect_clayton <-
cylcop::optCor(
copula = cyl_rect_combine(
copula = claytonCopula(),
low_rect = c(0, 0.5),
up_rect = "symmetric"
),
theta = theta,
x = x,
method = "tau"
)
start_rect_gumbel <-
cylcop::optCor(
copula = cyl_rect_combine(
copula = gumbelCopula(),
low_rect = c(0, 0.5),
up_rect = "symmetric"
),
theta = theta,
x = x,
method = "tau"
)Using these parameter values as starting points, we optimize copula parameters with a maximum pseudo likelihood approach.
copula_quadsec <-
cylcop::optML(
copula = cyl_quadsec(),
theta = theta,
x = x,
parameters = "a",
start = start_quadsec
)
copula_cubsec <-
cylcop::optML(
copula = cyl_cubsec(),
theta = theta,
x = x,
parameters = c("a", "b"),
start = start_cubsec,
traceOpt = T
)
copula_rect_frank <-
cylcop::optML(
copula = cyl_rect_combine(
copula = frankCopula(param = 2),
low_rect = c(0, 0.5),
up_rect = "symmetric"
),
theta = theta,
x = x,
parameters = "alpha",
start = start_rect_frank,
traceOpt = T
)With a Clayton copula with negative parameter the rectangular patchwork copula fits the data very badly, and we therefore set the lower bound close to 0.
copula_rect_clayton <-
cylcop::optML(
copula = cyl_rect_combine(
copula = claytonCopula(),
low_rect = c(0, 0.5),
up_rect = "symmetric"
),
theta = theta,
x = x,
parameters = "alpha",
start = start_rect_clayton,
traceOpt = T,
lower = 0.0001
)
copula_rect_gumbel <-
cylcop::optML(
copula = cyl_rect_combine(
copula = gumbelCopula(),
low_rect = c(0, 0.5),
up_rect = "symmetric"
),
theta = theta,
x = x,
parameters = "alpha",
start = start_rect_gumbel,
traceOpt = T
)With the following copula below, optimizing the upper bound of the lower rectangle (and with it the lower bound of the upper rectangle, because up_rect=“symmetric”) results either in rectangles of area 0 and a cyl_quadsec copula or in rectangles spanning the entire unit square, i.e. copula_rect_frank.
cylcop::optML(
copula = cyl_rect_combine(
copula = frankCopula(param = 2),
low_rect = c(0, 0.5),
up_rect = "symmetric",
background = cyl_quadsec()
),
theta = theta,
x = x,
parameters = c("alpha", "low_rect2", "bg_a"),
start = c(1.415, 0.357, 0.113),
lower = c(0.2, 0, 0),
upper = c(10, 0.45, 1 / (2 * pi)),
traceOpt = T
)Instead, we arbitrarily set the rectangles to \(R_1=[0.0,0.3]\times[0,1]\) and \(R_2=[0.7,1.0]\times[0,1]\)
copula_rect_frank_quadsec <- cylcop::optML(
copula = cyl_rect_combine(
copula = frankCopula(param = 2),
low_rect = c(0, 0.3),
up_rect = "symmetric",
background = cyl_quadsec()
),
theta = theta,
x = x,
parameters = c("alpha", "bg_a"),
start = c(1.415, 0.113),
lower = c(0.0001, 0),
traceOpt = T
)The copula giving the lowest AIC score was the cubic sections copula.
copula_cubsec## $copula
## Cub. sect. copula
## a = 0.06604285
## b = -0.1591549
##
## $logL
## [1] 274.8908
##
## $AIC
## [1] -545.7815 cylcop::wasserstein(copula_cubsec$copula, data=cbind(theta,x))## [1] 0.01894369sessionInfo()## R version 4.0.5 (2021-03-31)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19044)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=German_Switzerland.1252 LC_CTYPE=German_Switzerland.1252
## [3] LC_MONETARY=German_Switzerland.1252 LC_NUMERIC=C
## [5] LC_TIME=German_Switzerland.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] cylcop_0.1.0 copula_1.0-1 amt_0.1.4 forcats_0.5.1 stringr_1.4.0
## [6] dplyr_1.0.7 purrr_0.3.4 readr_2.0.2 tidyr_1.1.4 tibble_3.1.5
## [11] ggplot2_3.3.5 tidyverse_1.3.1
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.0-2 ellipsis_0.3.2 rgdal_1.5-27 fs_1.5.0
## [5] rstudioapi_0.13 circular_0.4-93 bit64_4.0.5 gsl_2.1-6
## [9] fansi_0.5.0 mvtnorm_1.1-1 lubridate_1.8.0 xml2_1.3.2
## [13] codetools_0.2-18 splines_4.0.5 doParallel_1.0.16 knitr_1.36
## [17] jsonlite_1.7.2 broom_0.7.9 dbplyr_2.1.1 GoFKernel_2.1-1
## [21] stabledist_0.7-1 shiny_1.7.1 compiler_4.0.5 httr_1.4.2
## [25] backports_1.2.1 assertthat_0.2.1 Matrix_1.3-4 fastmap_1.1.0
## [29] lazyeval_0.2.2 cli_3.0.1 later_1.3.0 htmltools_0.5.2
## [33] tools_4.0.5 gtable_0.3.0 glue_1.4.2 movMF_0.2-5
## [37] doRNG_1.8.2 tinytex_0.34 Rcpp_1.0.7 slam_0.1-48
## [41] jquerylib_0.1.4 cellranger_1.1.0 raster_3.4-13 vctrs_0.3.8
## [45] crawl_2.2.1 iterators_1.0.13 xfun_0.26 rbibutils_2.2.3
## [49] rvest_1.0.1 mime_0.12 lifecycle_1.0.1 rngtools_1.5.2
## [53] MASS_7.3-54 scales_1.1.1 vroom_1.5.5 hms_1.1.1
## [57] promises_1.2.0.1 Brobdingnag_1.2-6 parallel_4.0.5 yaml_2.2.1
## [61] momentuHMM_1.5.4 gridExtra_2.3 sass_0.4.0 stringi_1.7.5
## [65] highr_0.9 pcaPP_1.9-74 foreach_1.5.1 boot_1.3-28
## [69] Rdpack_2.1.2 rlang_0.4.11 pkgconfig_2.0.3 rgl_0.107.14
## [73] evaluate_0.14 lattice_0.20-45 htmlwidgets_1.5.4 bit_4.0.4
## [77] cowplot_1.1.1 transport_0.12-2 tidyselect_1.1.1 magrittr_2.0.1
## [81] R6_2.5.1 generics_0.1.0 ADGofTest_0.3 DBI_1.1.1
## [85] pillar_1.6.3 haven_2.4.3 mixR_0.2.0 withr_2.4.2
## [89] survival_3.2-13 sp_1.4-5 pspline_1.0-18 modelr_0.1.8
## [93] crayon_1.4.1 KernSmooth_2.23-20 utf8_1.2.2 plotly_4.10.0
## [97] tzdb_0.1.2 rmarkdown_2.11 viridis_0.6.1 grid_4.0.5
## [101] readxl_1.3.1 data.table_1.14.2 infotheo_1.2.0 CircStats_0.2-6
## [105] reprex_2.0.1 digest_0.6.28 xtable_1.8-4 httpuv_1.6.3
## [109] extraDistr_1.9.1 numDeriv_2016.8-1.1 stats4_4.0.5 munsell_0.5.0
## [113] viridisLite_0.4.0 bslib_0.3.1