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.
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.
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")
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
\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
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
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)
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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.