Simulating movement from a recharge dynamics model

[meryltheng]

Hi Brett,

I've been using your R package momentuHMM for simple movement simulations for a project I'm working on and it's been really great! I had some questions, I hope you don't mind.

For consistency, I was looking to simulate more complex habitat-driven movement using the same package and wondering if what I'm envisioning is possible. Specifically, would it be possible to simulate from a recharge dynamics model to replicate the empirical movement patterns seen in Fig. 26 of this document (African buffalo example; https://cran.r-project.org/web/packages/momentuHMM/vignettes/momentuHMM.pdf)? What I'm looking for is for short, tortuous movement in resource areas (charged state) and for long, directed movements towards the nearest resource area in the recharge state.

I've only managed this so far,

where the yellow dots represent habitat patches and coloured lines represent movement trajectories from two animals. I'm not quite sure if the animals are switching states, it appears the blue one might be. If so, is there a way to get the interpatch movement (during the recharge state) to be more directed?

For reference, here is the code I used:

"-------- Create simulation landscape --------"
r <- raster(ncol=1000, nrow=1000)
extent(r)<- c(0,1000,0,1000)
# create randomly located circular patches
patch_xy <- cbind(x=runif(50,100,900),y=runif(50,100,900))
patch_xy <- SpatialPoints(patch_xy)
patch_pol <- gBuffer(patch_xy, width = 10 ,quadsegs=round(1/10,0))
r <- rasterize(patch_pol,r)
proj4string(r) <- sp::CRS("+proj=utm +zone=53H +south +datum=WGS84") # assign dummy CRS
plot(r)

"-------- Prepare spatial covariate for model input --------"
# create distance to resource spatial covariate
dist2r <- distance(r)
plot(dist2r)
dist2r_scaled <- dist2r / mean(raster::values(terrain(dist2r,opt = "slope")),na.rm = T)

# calculate gradient
D_scaled <- ctmcmove::rast.grad(dist2r_scaled)
intercept <- raster(dist2r)
intercept[] <- 1
W <- stack(list("intercept" = intercept,
                "near_r" = dist2r < 5 ))
W_names <- names(W)

## orthogonalize W based on locations (adapted from vignette, not quite sure if I'm doing it right)
W_ortho <- W
W_path <- extract(x = W, y = cbind(x=runif(200,0,1000),y=runif(200,0,1000)))
W_tilde <- apply(W_path * c(0, rep(1,199)), 2, cumsum)
W_tilde_svd <- svd(W_tilde)
W_tilde_proj_mat <- W_tilde_svd$v %*% diag(W_tilde_svd$d^(-1))
W_mat <- as.matrix(W)
W_mat_proj <- W_mat %*% W_tilde_proj_mat
for(layer in 1:ncol(W_mat)){
  values(W_ortho[[layer]]) <- W_mat_proj[, layer]
  names(W_ortho[[layer]]) <- paste0("svd", layer)
}

spatialCovs <- list(W_intercept=W_ortho$svd1,
                    W_near_r=W_ortho$svd2,
                    dist2r=dist2r,
                    D.x=D_scaled$rast.grad.x,
                    D.y=D_scaled$rast.grad.y)


"-------- Excecute movement simulations --------"
# states
nbStates <- 2
stateNames <- c("charged", "discharged")

# draw random initial locations for each individual
N = 2
initialPositions <- array(500, c(N,2)) # start from centre of landscape
posits.df <- split(initialPositions, seq(nrow(initialPositions)))

cores=1
#### RUN MOMENTUHMM ####
simDat <- simData(nbAnimals=N, nbStates=2,
                  obsPerAnimal=1000, 
                  dist=list(mu="rw_mvnorm2"),  # bivariate normal random walk
                  DM=list(mu=list(mean.x=~mu.x_tm1+state2(D.x),
                                  mean.y=~mu.y_tm1+state2(D.y),
                                  sigma.x=~1,sigma.xy=~1,sigma.y=~1)),
                  # regression coeffs for mu params
                  Par=list(mu=c( 1, # mu.x_tm1 s1
                                 1, # mu.x_tm1 s2
                                 -2, # x strength to D (s2)
                                 1, # mu.y_tm1 s1
                                 1, # mu.y_tm1 s2
                                 -2, # y strength to D (s2)
                                 1, 5, # sigma.x (s1, s2)
                                 0, 0, # sigma.xy
                                 1, 5)), # sigma.y (s1, s2) 
                  ## specify recharge formula
                  formula = ~recharge(g0=~1,theta=~W_intercept+W_near_r),
                  beta = list(beta=matrix(c(0,-1,0,-1),2,2),g0=-0.3,theta=c(0,1,5)),
                  ## remove Markov property
                 # betaRef <- c(1,1), # make state 1 the reference state
                  #delta=c(0.5,0.5),
                  spatialCovs = spatialCovs,
                  #covs = covs,
                  mvnCoords = "mu",
                  initialPosition = posits.df,
                  ncores = cores,retrySims = 10,states=TRUE)

Let me know if I need to provide more information. Many thanks in advance!

Cheers,
Meryl

[bmcclintock]

The recharge model of Hooten et al is non-Markov (i.e. the probability of switching to the "charged" or "discharged" state is not dependent on the current state and is solely a function of the current value of the recharge function). The beta=matrix(c(0,-1,0,-1),2,2) constraints can be used to accomplish this, but they only make sense if the reference state is the same for the "charged" and "discharged" states. This can be done by setting betaRef = c(1,1).

Orthogonalization of the recharge covariates is not necessary for simulation; this is more just a "trick" for better behavior during model fitting. I'd suggest skipping this step, in which case theta=~near_r could be used in the recharge formula. Then theta[1] should be negative (i.e. discharge) and theta[2]>abs(theta[1]) should be positive (i.e. recharge).

There is a lot of interplay between the parameters, covariates values (and the magnitude of their gradients), and spatial scale in this model. It therefore requires a fair bit tweaking to achieve the desired behavior. For example, if the "charged" state is just a simple random walk with short steps as described above, then an individual is likely to never leave the current patch (and its recharge function will "blow up" to a large positive value).

Using your code with betaRef=c(1,1) we get:

plotSpatialCov(simDat, spatialCovs$dist2r, states=simDat$states, col=stateCols, stateNames = stateNames, size=1.25, alpha=0.35)

# plot true recharge function
for(i in unique(simDat$ID)){
  ind <- which(simDat$ID==i)
  g_t <- cumsum(c(beta0$g0,cbind(rep(1,length(ind)),simDat$W_near_r[ind])%*%beta0$theta[2:3]))
  plot(1:(length(ind)+1),g_t,ylab="g(t)",xlab="t",main=paste("ID",i,"recharge function"),pch=20,col=stateCols[simDat$states[ind]])
  segments(y0=g_t[-length(g_t)],x0=1:(length(ind)-1),y1=g_t[-1],x1=2:length(ind),col=stateCols[simDat$states[ind]])
  legend("topright",legend=stateNames,pch=20,col=stateCols,bty="n")
  abline(h=0,lty=2)
}

So both individuals start in the "discharged" state, promptly move towards the nearest patch, and the recharge function "blows up". Thus the scale and relative values of the recharge function parameters need to be "balanced" to allow for both discharge and recharge as the individual moves through the landscape, and movement in the "charged" state needs to allow for movement to other patches (otherwise they are likely to get "stuck" in the nearest patch). Although not necessarily an issue, also note that in this model the gradient is zero within patches, so movement for both the "charged" and "discharged" states is a simple random walk when individuals are within a patch.

Here is an example where the recharge function is more balanced, but the "charged" state is still just a simple random walk:

library(raster)
library(momentuHMM)
library(rgeos)
"-------- Create simulation landscape --------"
r <- raster(ncol=1000, nrow=1000)
extent(r)<- c(0,1000,0,1000)
# create randomly located circular patches
set.seed(6, kind="L'Ecuyer-CMRG", normal.kind="Inversion")
patch_xy <- cbind(x=runif(50,100,900),y=runif(50,100,900))
patch_xy <- SpatialPoints(patch_xy)
patch_pol <- gBuffer(patch_xy, width = 10 ,quadsegs=round(1/10,0))
r <- rasterize(patch_pol,r)
proj4string(r) <- sp::CRS("+proj=utm +zone=53H +south +datum=WGS84") # assign dummy CRS
plot(r)

"-------- Prepare spatial covariate for model input --------"
# create distance to resource spatial covariate
dist2r <- distance(r)
plot(dist2r)
dist2r_scaled <- dist2r / mean(raster::values(terrain(dist2r,opt = "slope")),na.rm = T)

# calculate gradient
D_scaled <- ctmcmove::rast.grad(dist2r_scaled)
intercept <- raster(dist2r)
intercept[] <- 1
W <- stack(list("intercept" = intercept,
                "near_r" = dist2r < 5 ))
W_names <- names(W)

spatialCovs <- list(W_intercept=W$intercept,
                    W_near_r=W$near_r,
                    dist2r=dist2r,
                    D.x=D_scaled$rast.grad.x,
                    D.y=D_scaled$rast.grad.y)


"-------- Excecute movement simulations --------"
# states
nbStates <- 2
stateNames <- c("charged", "discharged")

beta0 <- list(beta=matrix(c(0,-1,0,-1),2,2),g0=-0.3,theta=c(-0.25,0.35))

cores=1

# draw random initial locations for each individual
N = 2
initialPositions <- array(500, c(N,2)) # start from centre of landscape
posits.df <- split(initialPositions, seq(nrow(initialPositions)))

#### RUN MOMENTUHMM ####
set.seed(1, kind="L'Ecuyer-CMRG", normal.kind="Inversion")
simDat <- simData(nbAnimals=N, nbStates=2,
                  obsPerAnimal=5000, 
                  dist=list(mu="rw_mvnorm2"),  # bivariate normal random walk
                  DM=list(mu=list(mean.x=~mu.x_tm1+state2(D.x),
                                  mean.y=~mu.y_tm1+state2(D.y),
                                  sigma.x=~1,sigma.xy=~1,sigma.y=~1)),
                  # regression coeffs for mu params
                  Par=list(mu=c( 1, # mu.x_tm1 s1
                                 1, # mu.x_tm1 s2
                                 -2, # x strength to D (s2)
                                 1, # mu.y_tm1 s1
                                 1, # mu.y_tm1 s2
                                 -2, # y strength to D (s2)
                                 0.25, 1.25, # sigma.x (s1, s2)
                                 0, 0, # sigma.xy
                                 0.25, 1.25)), # sigma.y (s1, s2) 
                  ## specify recharge formula
                  formula = ~recharge(g0=~1,theta=~W_near_r),
                  beta = beta0,
                  ## remove Markov property
                  betaRef = c(1,1), # make state 1 the reference state
                  delta=c(0.5,0.5),
                  spatialCovs = spatialCovs,
                  #covs = covs,
                  mvnCoords = "mu",
                  initialPosition = posits.df,
                  ncores = cores,#retrySims = 10,
                  states=TRUE,
                  stateNames = stateNames)

# plot simulated data
stateCols <- c("#009E73", "#F0E442")
plot(simDat,dataNames=c("mu.x","mu.y"),compact=TRUE,ask=FALSE)
plotSpatialCov(simDat, spatialCovs$dist2r, states=simDat$states, col=stateCols, stateNames = stateNames, size=1.25, alpha=0.35)

# plot true recharge function
for(i in unique(simDat$ID)){
  ind <- which(simDat$ID==i)
  g_t <- cumsum(c(beta0$g0,cbind(rep(1,length(ind)),simDat$W_near_r[ind])%*%beta0$theta))
  plot(1:(length(ind)+1),g_t,ylab="g(t)",xlab="t",main=paste("ID",i,"recharge function"),pch=20,col=stateCols[simDat$states[ind]])
  segments(y0=g_t[-length(g_t)],x0=1:(length(ind)-1),y1=g_t[-1],x1=2:length(ind),col=stateCols[simDat$states[ind]])
  legend("topright",legend=stateNames,pch=20,col=stateCols,bty="n")
  abline(h=0,lty=2)
}

So here the recharge function is better behaved, but we don't get much movement between patches. But now if we were to add a potential function for the "charged" state that tends to move "charged" individuals away from patches:

# add potential function away from patches for state 1
set.seed(1, kind="L'Ecuyer-CMRG", normal.kind="Inversion")
simDat2 <- simData(nbAnimals=N, nbStates=2,
                  obsPerAnimal=5000, 
                  dist=list(mu="rw_mvnorm2"),  # bivariate normal random walk
                  DM=list(mu=list(mean.x=~mu.x_tm1+D.x,
                                  mean.y=~mu.y_tm1+D.y,
                                  sigma.x=~1,sigma.xy=~1,sigma.y=~1)),
                  # regression coeffs for mu params
                  Par=list(mu=c( 1, # mu.x_tm1 s1
                                 1, # x strength away from D (s1)
                                 1, # mu.x_tm1 s2
                                 -2, # x strength to D (s2)
                                 1, # mu.y_tm1 s1
                                 1, # x strength away from D (s1)
                                 1, # mu.y_tm1 s2
                                 -2, # y strength to D (s2)
                                 0.25, 1.25, # sigma.x (s1, s2)
                                 0, 0, # sigma.xy
                                 0.25, 1.25)), # sigma.y (s1, s2) 
                  ## specify recharge formula
                  formula = ~recharge(g0=~1,theta=~W_near_r),
                  beta = beta0,
                  ## remove Markov property
                  betaRef = c(1,1), # make state 1 the reference state
                  delta=c(0.5,0.5),
                  spatialCovs = spatialCovs,
                  #covs = covs,
                  mvnCoords = "mu",
                  initialPosition = posits.df,
                  ncores = cores,#retrySims = 10,
                  states=TRUE,
                  stateNames = stateNames)

# plot simulated data
plot(simDat2,dataNames=c("mu.x","mu.y"),compact=TRUE,ask=FALSE)
plotSpatialCov(simDat2, spatialCovs$dist2r, states=simDat2$states, col=stateCols, stateNames = stateNames, size=1.25, alpha=0.35)

# plot true recharge function
for(i in unique(simDat2$ID)){
  ind <- which(simDat2$ID==i)
  g_t <- cumsum(c(beta0$g0,cbind(rep(1,length(ind)),simDat2$W_near_r[ind])%*%beta0$theta))
  plot(1:(length(ind)+1),g_t,ylab="g(t)",xlab="t",main=paste("ID",i,"recharge function"),pch=20,col=stateCols[simDat2$states[ind]])
  segments(y0=g_t[-length(g_t)],x0=1:(length(ind)-1),y1=g_t[-1],x1=2:length(ind),col=stateCols[simDat2$states[ind]])
  legend("topright",legend=stateNames,pch=20,col=stateCols,bty="n")
  abline(h=0,lty=2)
}

So we now see some more movement among the patches and frequent switches between the "charged" and "discharged" states. Hopefully you can play around with this code and get a feel for the interplay between the movement models for the "charged" and "discharged" states, the covariates, and the model parameters. There are lots of other options for movement of the "charged" state, including adding correlated movement (via the crw special function) or alternative potential functions. One could also potentially have patches turn "on" and "off" by making spatialCovs time dependent.

[meryltheng]

Hi Brett,

This is fantastic and exactly what I was looking for. I'll keep playing around with what you've provided and try to add in more features (such as the crw).

Thank you so much for the solutions and helpful explanations!

Cheers,
Meryl