In papayoun/Rhabit: Estimation of animal habitat selection using the Langevin movement model

rm(list = ls())
knitr::opts_chunk$set(
  message = F, error = F, warning = F,
  comment = NA
)

library(Rhabit)
library(viridis)
library(ggplot2)
library(gridExtra)
library(abind)
library(dplyr)

This vignette illustrates the use of the Rhabit package to simulate data from a Langevin model and to estimate the parameters.

Simulation of the Langevin model

In the Langevin model, we assume that the movement is directed by a field of covariates. Those fields of covariates are mandatory to simulate a path according to the model.

Simulating covariate fields

In the presented example, the path is observed on a square, centered on the origine and whose sides equal 30.

In order to g a enerate a field of covariates, we propose to simulate two Mattern processes with the same parameters parameters $(\rho, \nu, \sigma^2)$ The covriates are obtained by applying a log transformation to the square of the Mattern process.

set.seed(1)# repeatability
lim <- c(-15, 15, -15, 15) # limits of map
resol <- 0.5 # grid resolution
rho <- 4; nu <- 1.5; sigma2 <- 10# Matern covariance parameters
mean_function <- function(z){# mean function
  -log(3 + sum(z^2))
}
J <- 2 # number of covariates

Then, simulating

covariates <- replicate(J, 
                        simSpatialCov(lim, nu, rho, sigma2, resol = resol,
                                      mean_function = mean_function,
                                      raster_like = T),
                        simplify = F)

The resulting covariates are available in the package~:

data("covariates", package = "Rhabit")

Representation of the fileds of covariates

ggopts <- theme_light()+theme(text = element_text(size = 20),
                axis.text = element_blank(),
                axis.title = element_blank(),
                legend.key.height= unit(3, "line"),
                strip.background = element_rect(fill="white",colour = "grey") ,
                strip.text = element_text(size = 25, color = "black"))
levels <- factor(paste("Covariate", 1:J), levels = paste("Covariate", 1:J))
cov_df <- do.call(rbind.data.frame,
                 lapply(1:J,
                        function(j){
                          Rhabit::rasterToDataFrame(covariates[[j]], levels[j])
                        }))
ggplot(cov_df, aes(x,y)) + geom_raster(aes(fill = val)) +
  coord_equal() + scale_fill_viridis(name = "Value") + facet_wrap(level~.) +
  ggopts

The resulting Utilization Distribution

\begin{equation} \rmd\X_t = \dfrac{\gamma^2}{2} \nabla \log \ud(\X_t) \rmd t + \gamma \rmd \bm{W}_t,\quad \X_0 =\bm{x}_0. \end{equation}

beta_true <- c(10, 5)
ud_rast <- covariates[[1]] # initialization
dx <- diff(ud_rast$x)[1]
dy <- diff(ud_rast$x)[1]
ud_rast$z <- exp(Reduce("+",
                    lapply(1:J, function(j)
                      dx * dy * beta_true[j] * covariates[[j]]$z)))
ud_rast$z <- ud_rast$z / sum(ud_rast$z)
ud_df <- rasterToDataFrame(ud_rast)
ggplot(ud_df, aes(x,y)) + geom_raster(aes(fill = val)) +
  coord_equal() + scale_fill_viridis(name = "Value") +
  ggopts

Simulation

Setting simulation parameters:

t_final <- 10
n_obs <- 10001
times <- seq(0, t_final, length.out = n_obs)

Then, one can simulate using cov_list

set.seed(1) #repeatability
tracks <- simLangevinMM(beta = beta_true, gamma2 = 10,
                        times = times, loc0 = c(0, 0),
                        cov_list = covariates, keep_grad = T)

These simulated tracks are a data set that can be loaded

data("tracks", package = "Rhabit")

And can be plotted

ggplot(ud_df, aes(x,y)) + geom_raster(aes(fill = val)) +
  coord_equal() + scale_fill_viridis(name = "Value") +
  geom_path(data = tracks[, 1:2]) +
  ggopts

Estimation

The dataset is expected to contain columns - x and y which specify the relocations - t , for the corresponding times - for each variable of interest, let's say V1, two columns V1_x and V1_y which contain the computed gradient of covariates V1 at the recorded location.

## test with the formula
fitted_langevin <- fit_langevin_ud(
  cbind(x,y) ~ grad_c1 + grad_c2, data = tracks)

The effects of the covariates can be obtained via

coef(fitted_langevin)

while the speed parameter is obtained through

speed_coef(fitted_langevin)

Testing thhe effects

The function \code{summary} tests for every covariate, if the corresponding parameter equals 0 or not.

The Null hypothesis for the speed parameter is $\gamma^2 =0$.

summary(fitted_langevin)

The confidence intervals for each parameter are available with.

fitted_langevin$CI

The expected trend in the movement at each observed location is given by

head(fitted_langevin$predicted, n = 10)

papayoun/Rhabit documentation built on July 19, 2023, 8:04 p.m.