jmsigner / amt

37 stars 13 forks source link

Parameters for the gamma and vonmises in make_issf()? #111

Closed MarieAugerMethe closed 3 months ago

MarieAugerMethe commented 3 months ago

Hi Johannes and team,

Thank you for creating such a wonderfully useful package!

I'm trying to create steady-state UDs from an issf modelled with glmmTMB. I have been following the instructions of your various papers on the topic (Signer et al. 2017, 2019, 2024), and I think I have it mostly sorted out. However, I'm not sure what to use for the parameters of the movement kernel. My understanding from the Avgar et al. 2016 issf paper in MEE was that the movement kernel was using the observed steps to calculate the base shape and scale parameters of the step length and then using the coefficients associated with the sl and ln(sl) to modify it, and similarly for the turning angle. But I think I must be misunderstanding something as I can't replicate the issf object with make_issf(). Here is an example. This is based on the data from your recent paper (Signer et al. 2024), which is available on zenodo: https://zenodo.org/records/10160168

# Based on data from Signer et al. 2024 MEE
library(terra)
library(tidyverse)
library(amt) # dev version needed: remotes::install_github("jmsigner/amt")
library(here) # here::here() should point to the folder "supplement2" containing 
library(raster)
# data/cilla.rds and data/env_covar.rds
#https://zenodo.org/records/10160168

set.seed(123333)

cilla <- read_rds(here::here("cilla.rds"))
env <- rast(read_rds(here::here("env_covar.rds")))
env

# Create random steps
ssf_cilla <- cilla %>% steps_by_burst() %>% random_steps(n_control = 100) |> 
  extract_covariates(env, where = "both") %>% 
  mutate(hour = hour(t1_) + minute(t1_) / 60) %>% 
  filter(complete.cases(.))

m_0 <- fit_clogit(ssf_cilla, case_ ~ cos(ta_) + sl_ + log(sl_) + 
                    water_dist_end + strata(step_id_))

# Shape parameter from the fitted model
m_0$sl_$params$shape

# Shape parameter from the observed step length (from movement data) - I expected this to be the same as above
fit_distr(subset(ssf_cilla, case_)$sl_, "gamma")$params$shape

# Shape parameter using the used and available - just looked in case
fit_distr(ssf_cilla$sl_, "gamma")$params$shape

m_0$ta_$params$kappa
fit_distr(subset(ssf_cilla, case_)$ta_, "vonmises")$params$kappa
fit_distr(ssf_cilla$ta_, "vonmises")$params$kappa

As expected, the kappa parameter estimated from the observed data appears to be the same as that for the fitted model. However, the shape parameters are not the same.

I'm sure I'm missing something obvious, and my apologies in advance.

Many thanks!

Marie

jmsigner commented 3 months ago

Hi @MarieAugerMethe, I think I figured out what is going on:

library(amt) # dev version needed: remotes::install_github("jmsigner/amt")
library(tidyverse)

set.seed(123333)

cilla <- read_rds("cilla.rds")

# Create random steps
ssf_cilla <- cilla %>% steps_by_burst() %>% random_steps(n_control = 100) %>% 
    filter(complete.cases(.))

m_0 <- fit_clogit(ssf_cilla, case_ ~ cos(ta_) + sl_ + log(sl_) + strata(step_id_))

# Shape parameter from the fitted model: This are the same
attr(ssf_cilla, "sl_")$params
m_0$sl_$params

# Also the same
cilla_bursts <- cilla %>% steps_by_burst()
fit_distr(cilla_bursts$sl_, "gamma")$params

# this is different
fit_distr(filter(ssf_cilla, case_)$sl_, "gamma")$params

# Wher does this difference come from: 
cilla_bursts |> nrow()

# We have two bursts
unique(ssf_cilla$burst_)

# When creating random steps, the first step of each burst are lost (because we do not know the initial direction). 
filter(ssf_cilla, case_) |> nrow()
MarieAugerMethe commented 3 months ago

Thank you!! Much appreciated!

In summary: the distribution used to sample the available/random steps is based on all steps, but the first step of each burst gets removed when creating the random steps. Thanks!