Lecture 12

class: center, middle, inverse, title-slide

.title[
# Lecture 12
]
.subtitle[
## Simple Movement Models
]
.author[
### Spring 2025
]

---

# Readings:

McClintock, B. T., B. Abrahms, R. B. Chandler, P. B. Conn, S. J. Converse, R. L. Emmet, B. Gardner, N. J. Hostetter, and D. S. Johnson. 2022. An integrated path for spatial capture–recapture and animal movement modeling. Ecology 103(10): e3473. 10.1002/ecy.3473

---
## Individual-Level Models

So far in this class we've discussed models that group individuals together to model population-level processes.

But not all individuals are created equal.

Individual-level models are often more informative about population *mechanisms* because population dynamics ultimately occur at the individual level.

---
## Individual-Level Models

Central to the study of population dynamics is an argument over the importance of individual-level complexity in model structure.

Does individual variation:

a) create noise that blocks the identification of universal ecological patterns

+ John H. Lawton. Are There General Laws in Ecology?  *Oikos*, Vol. 84, No. 2. (Feb., 1999), pp. 177-192.

+ To discover general patterns, ecology may need to pay less attention to the 'middle ground' of community ecology and rely less on reductionism and experimental manipulation

---
## Individual-Level Models

Central to the study of population dynamics is an argument over the importance of individual-level complexity in model structure.

Does individual variation:

a) create noise that blocks the identification of universal ecological patterns

+  Hubbell, S.P. (2005), Neutral theory in community ecology and the hypothesis of functional equivalence. Functional Ecology, 19: 166-172.

+ All individual level patterns are caused by chance. Population patterns are caused by niche requirements and interactions with other species at the ecosystem level

---
## Individual-Level Models

Central to the study of population dynamics is an argument over the importance of individual-level complexity in model structure.

Does individual variation:

a) create noise that blocks the identification of universal ecological patterns

b) individual-level behaviors and demographic differences fundamentally important in their own right

+ Michael Huston, Donald DeAngelis, Wilfred Post, New Computer Models Unify Ecological Theory: Computer simulations show that many ecological patterns can be explained by interactions among individual organisms, BioScience, Volume 38, Issue 10, November 1988, Pages 682–691

+ Individuals are the reason population patterns exist. We can't truly understand the reason for the patterns if we ignore the individual

---
## Individual-Level Models

Central to the study of population dynamics is an argument over the importance of individual-level complexity in model structure.

Does individual variation:

a) create noise that blocks the identification of universal ecological patterns

b) individual-level behaviors and demographic differences fundamentally important in their own right

+ DeAngelis, D. L. 1996. Individual-based models and approaches in ecology: Populations, communities and ecosystems,Chapman & Hall, New York.

+ Most experimental data is initially collected at the individual level. It becomes harder to test models for validity if we have to transform our data to fit the model

---
## Individual-Level Models

Central to the study of population dynamics is an argument over the importance of individual-level complexity in model structure.

Does individual variation:

a) create noise that blocks the identification of universal ecological patterns

b) individual-level behaviors and demographic differences fundamentally important in their own right

This argument has been flourishing since the mid-1950s with both sides bemoaning the failure of mathematical models to realistically reproduce complex dynamics in ecology. This issue remains hotly contested to this day.

---
## Movement Models

As noted by DeAngelis (1996), we often collect data at the individual-level, so it makes sense to analyse it that way.

A great example of this is movement data.

One of the simplest movement models we can use is a random walk model, also called Brownian Bridge movement.

`$$s_{i,x,t=1} \sim Uniform(x_1, x_2)$$`

`$$s_{i,y,t=1} \sim Uniform(y_1, y_2)$$`

`$$s_{i, x, t+1} \sim Normal(s_{i, x, t}, \sigma ^2 )$$`

`$$s_{i, y, t+1} \sim Normal(s_{i, y, t}, \sigma^2 )$$`

---
## Movement Models

Under this model, animals move entirely randomly. The only parameter of interest is `$\sigma$`, which controls the average distance between points.

`$$\Large s_{i, x, t+1} \sim Normal(s_{i, x, t}, \sigma ^2 )$$`

`$$\Large s_{i, y, t+1} \sim Normal(s_{i, y, t}, \sigma^2 )$$`

You will often see these models written in the literature as:

`$$\Large \mu_t = Normal(\mu_{t-1}, \sigma^2)$$`
or:

`$$\Large \mu_t = \mu_{t-1} + \epsilon_{t} \sim Normal(0, \sigma^2)$$`
These are all equivalent

---
## Random Walk Movement Models

Here's a simple simulation of a random walk we used long, long ago in Lab 2.

``` r
sim_move <- function(sigma, n, nocc){
s <- array(NA, c(n,2,nocc))
for(i in 1:n){ #each individual
 s[i,, 1] <- runif(2, 0, 1) #coords time 1 
 for(t in 2:nocc){
 s[i,, t] <- rnorm(2, s[i,, t-1], sd = sigma)
 }
}
return(s)
}
```

---
## Random Walk Movement Models

If `$\sigma$` is very small, the animals don't move that much:

``` r
nocc <- 50
s <- sim_move(sigma = .02, n = 1, nocc = nocc)
plot(s[1,1,],s[1,2,], asp=1, , xlim = c(-2, 2), ylim = c(-2, 2), type = 'b')
arrows(s[1,1,1:(nocc-1)], s[1,2, 1:(nocc-1)],
 s[1,1, 2:nocc], s[1,2, 2:nocc], length=0.1, col="blue")
```

---
## Random Walk Movement Models

If `$\sigma$` is large the animals make larger movements:

``` r
nocc <- 50
s <- sim_move(sigma = .5, n = 1, nocc = nocc)
plot(s[1,1,],s[1,2,], asp=1, xlim = c(-2, 2), ylim = c(-2, 2), type = 'b')
arrows(s[1,1,1:(nocc-1)], s[1,2, 1:(nocc-1)],
 s[1,1, 2:nocc], s[1,2, 2:nocc], length=0.1, col="blue")
```

---
## Correlated Random Walk

As with all simple models, they are useful if you don't have any idea about movement, but it is unlikely to fully capture the true movement patterns of animals in nature.

For one thing, animals don't just wander around the world aimlessly (usually).

What would our model look like if we accounted for correlation ('short-term directional persistence')?

---
## Correlated Random Walk

Consider this model:

`$$\Large \mu_t = \mu_{t-1} + \epsilon_{t} \sim Normal(0, \sigma^2)$$`

--

Under this formulation, each location is only dependent on the location one time step before it. To change that, we can calculate something like this:

`$$\Large\mu_t = \mu_{t-1} + \gamma(\mu_{t-1} - \mu_{t-2}) + \epsilon_{t} \sim Normal(0, \sigma^2)$$`
 
where `$\gamma \in [0, 1]$` is a 'damping parameter' that impacts how strongly the past locations affect the current location. 
---
## Correlated Random Walk

Let's simulate how that looks:

``` r
sim_move2 <- function(sigma, n, nocc, gamma){
s <- array(NA, c(n,2,nocc))
for(i in 1:n){
 s[i,, 1] <- runif(2, 0, 1) #coords time 1 
 s[i,, 2] <- rnorm(2, s[i,, 1], sd = sigma) #coords time 2 
 for(t in 3:nocc){
 s[i,, t] <- rnorm(2, s[i,, t-1] + gamma*(s[i,, t-1] - s[i,, t-2]), sd = sigma)
 }
}
return(s)
}
```

---
## Correlated Random Walk

When `$\gamma$` is high, movement is highly correlated. The animal walks in a path:

``` r
s <- sim_move2(sigma = .5, n = 1, nocc = 50, gamma = 1)
plot(s[1,1,],s[1,2,], asp=1, type = 'b')
arrows(s[1,1,1:(nocc-1)], s[1,2, 1:(nocc-1)],
 s[1,1, 2:nocc], s[1,2, 2:nocc], length=0.1, col="blue")
```

<img src="12_MovementModels_files/figure-html/unnamed-chunk-6-1.png" width="288" style="display: block; margin: auto;" />
---
## Correlated Random Walk

When `$\gamma$` is close to 0, movement reverts to a random walk.

``` r
s <- sim_move2(sigma = .5, n = 1, nocc = 50, gamma = .01)
plot(s[1,1,],s[1,2,], asp=1, type = 'b')
arrows(s[1,1,1:(nocc-1)], s[1,2, 1:(nocc-1)],
 s[1,1, 2:nocc], s[1,2, 2:nocc], length=0.1, col="blue")
```

---
## Correlated Random Walk with Rotation

What if want to include things like turn angles?

--
We can easily accommodate this by adding in a 'rotational component' parameter, `$R$`.

`$$\mathbf{R} = \gamma
\begin{bmatrix}
cos(\beta) & -sin(\beta)\\
sin(\beta) & cos(\beta)
\end{bmatrix}$$`

- `$\gamma$` represents the strength of the correlation between past movements and the current position.

- `$\beta \in [-\pi, \pi)$` is the turn angle of the individual

If `$\mathbf{R} = 0$`, we have a random walk model

If `$\beta = 0$` we have a tendancy to walk in a straight line

---
## Correlated Random Walk  with Rotation

This is all well and good, but how do we code this? This requires a little bit of matrix math.

`$$\Large \gamma
\begin{bmatrix}
cos(\beta) & -sin(\beta)\\
sin(\beta) & cos(\beta)
\end{bmatrix} * (\mu_{t-1} - \mu_{t-2}) =$$`

`$$\Large\begin{bmatrix}
\gamma*cos(\beta) & \gamma*-sin(\beta)\\
\gamma*sin(\beta) & \gamma*cos(\beta)
\end{bmatrix}*\begin{bmatrix} s_{i, x, t-1} - s_{i, x, t-2} \\ s_{i, y, t-1} - s_{i, y, t-2}\end{bmatrix} =$$`

`$$\Large\begin{bmatrix}
\gamma*cos(\beta)*(s_{i, x, t-1} - s_{i, x, t-2}) + \gamma*-sin(\beta)*(s_{i, y, t-1} - s_{i, y, t-2})\\
\gamma*sin(\beta)*(s_{i, x, t-1} - s_{i, x, t-2}) + \gamma*cos(\beta)*(s_{i, y, t-1} - s_{i, y, t-2})\end{bmatrix}$$`
---
## Correlated Random Walk with Rotation

Luckily for us, we don't have to worry about this math! We can just use the `%*%` function in R and it will figure out the matrix math.

Here's what that looks like in code:

``` r
sim_move3 <- function(sigma, n, nocc, gamma, beta){
s <- array(NA, c(n,2,nocc))
R <- gamma*matrix(c(cos(beta), sin(beta),
 -sin(beta), cos(beta)),
 2, 2) #rotational matrix
for(i in 1:n){
 s[i,, 1] <- runif(2, 0, 1) #coords time 1 
 s[i,, 2] <- rnorm(2, s[i,, 1], sd = sigma) #coords time 2 
 for(t in 3:nocc){
 s[i,, t] <- rnorm(2, s[i,, t-1] + R %*% (s[i,, t-1] - s[i,, t-2]), sd = sigma)
 }
}
return(s)
}
```

---
## Correlated Random Walk with Rotation

Beta is the turning angle in radians, so we can make some fun patterns if we chose specific betas. Here's some 90° turns to the left:

``` r
s <- sim_move3(sigma = .5, n = 1, nocc = 50, gamma = 1, beta = 3.14/2) #90° turns
plot(s[1,1,],s[1,2,], asp=1, type = 'b')
arrows(s[1,1,1:(nocc-1)], s[1,2, 1:(nocc-1)],
 s[1,1, 2:nocc], s[1,2, 2:nocc], length=0.1, col="blue")
```

---
## Correlated Random Walk with Rotation

Beta is the turning angle in radians, so we can make some fun patterns if we chose specific betas. Here's some 30° left turns when the correlation between locations is weaker:

``` r
s <- sim_move3(sigma = .5, n = 1, nocc = 50, gamma = .5, beta = 3.14/6) #30° turns
plot(s[1,1,],s[1,2,], asp=1, type = 'b')
arrows(s[1,1,1:(nocc-1)], s[1,2, 1:(nocc-1)],
 s[1,1, 2:nocc], s[1,2, 2:nocc], length=0.1, col="blue")
```

---
## In Class Exercise

Simulate an individual across 100 time steps with the following constraints:

- `set.seed(10)`

- `$\gamma$` = .85

- `$\sigma$` = .1

- The individual only makes right turns at 45 degree angles (it has a deformed right wing, it's very sad)

---
## In Class Exercise

``` r
set.seed(10)
nocc <- 100
s <- sim_move3(sigma = .1, n = 1, nocc = nocc, gamma = .85, beta = 3.14*7/4)
plot(s[1,1,],s[1,2,], asp=1, type = 'b')
arrows(s[1,1,1:(nocc-1)], s[1,2, 1:(nocc-1)],
 s[1,1, 2:nocc], s[1,2, 2:nocc], length=0.1, col="blue")
```

Simulate an individual across 100 time steps with the following constraints:

- `set.seed(10)`

- `$\gamma$` = .8

- The individual only makes left turns at 45 degree angles (it has a deformed right wing, it's very sad)

---
## Biased Random Walks

Movement ecologists, are often interested in home ranges - where does the animal spend most of its time? How much space does that take up? Where is the home range center?

Right now we have some sense of parameters that might control movement direction and correlation, but we haven't written anything into our model that encourages the individual to stay in one general area. 
--

Similar to what we have done before, we can add an attraction parameter `$\zeta \in [0, 1]$` that controls attraction to some specific location,  `$\mathbf{a}$`

`$$\Large\mu_t = \mu_{t-1} + \zeta(\mathbf{a} - \mu_{t-1}) + \epsilon_{t} \sim Normal(0, \sigma^2)$$`

---
## Biased Random Walks

``` r
sim_move4 <- function(sigma, n, nocc, zeta){
a <- runif(2) #attraction point
s <- array(NA, c(n,2,nocc))
for(i in 1:n){
 s[i,, 1] <- a #coords time 1 
 for(t in 2:nocc){
 s[i,, t] <- rnorm(2, s[i,, t-1] + zeta*(a - s[i,, t-1]), sd = sigma)
 }
}
return(s)
}
```

---
## Biased Random Walks

When `$\zeta$` is large, we see more obvious central tendency.

``` r
set.seed(10)
nocc <- 250
s <- sim_move4(sigma = .5, n = 1, nocc = nocc, zeta = .97)
plot(s[1,1,],s[1,2,], asp=1, type = 'b', ylim = c(-4, 4))
arrows(s[1,1,1:(nocc-1)], s[1,2, 1:(nocc-1)],
 s[1,1, 2:nocc], s[1,2, 2:nocc], length=0.1, col="blue")
```

---
## Biased Random Walks

When `$\zeta$` is small, we see get behavior closer to the random walk.

``` r
set.seed(10)
nocc <- 250
s <- sim_move4(sigma = .5, n = 1, nocc = nocc, zeta = .005)
plot(s[1,1,],s[1,2,], asp=1, type = 'b', ylim = c(-4, 4))
arrows(s[1,1,1:(nocc-1)], s[1,2, 1:(nocc-1)],
 s[1,1, 2:nocc], s[1,2, 2:nocc], length=0.1, col="blue")
```

---
## Random Walks in Nimble

Let's see how we might use these models to analyze actual data.

In the WILD8370 package, you'll find some location data from 6 collared white tailed deer in Arkansas. Note that in practice you would want a larger sample size.

``` r
data(collars)
colnames(collars)
```

```
##  [1] "Collar.ID"        "SCTS..UTC."       "Acq..Time..UTC."  "Origin"          
##  [5] "Latitude.deg."    "Longitude.deg."   "Altitude.m."      "DOP"             
##  [9] "Fix.Type"         "Mortality.Status" "Temp..C."         "Main.V."         
## [13] "Beacon..V."       "x"                "y"                "date1"           
## [17] "year"
```

---
## Random Walks in Nimble

Let's start by plotting out the deer, just to see what our data looks like.

``` r
ggplot(collars, aes(x = x/1e6, y= y/1e6))+
    geom_line()+
    facet_wrap(~Collar.ID, scales = 'free')+
    theme_bw()
```

---
## Random Walks in Nimble

We'll first fit the most basic movement model we can to this data - the random walk.

`$$s_{i,x,t=1} \sim Uniform(x_1, x_2)$$`

`$$s_{i,y,t=1} \sim Uniform(y_1, y_2)$$`

`$$s_{i, x, t+1} \sim Normal(s_{i, x, t}, \sigma ^2 )$$`

`$$s_{i, y, t+1} \sim Normal(s_{i, y, t}, \sigma^2 )$$`
---
## Random Walks in Nimble

We'll first fit the most basic movement model we can to this data - the random walk. The only complicated part of this model is that not all animals have the same number of locations, so we have to make sure our `t` loops accounts for this:

``` r
r_walk <- nimbleCode({
 for(i in 1:nind){
 s[i,1,1] ~ dunif(xmin, xmax)
 s[i,2,1] ~ dunif(ymin, ymax)
 for(t in 2:nTelemetry[i]){
 s[i,1,t] ~ dnorm(s[i,1,t-1], sd = sigma)
 s[i,2,t] ~ dnorm(s[i,2,t-1], sd= sigma)
 }
 }
 sigma ~ dgamma(1,1)
})
```

---
## Random Walks in Nimble

We need to clean up the data. We need the coordinates to be in meters (UTMs) and I find it easier to use small numbers, so I divide the UTMS by 1e6.

``` r
collars$Collar.ID_f <- as.numeric(as.factor(collars$Collar.ID))
nind <- length(unique(collars$Collar.ID_f))
nTelemetry <- as.vector(table(collars$Collar.ID_f))
s <- array(NA, c(nind, 2, max(nTelemetry)))
for(i in 1:nind){
 me <- collars[collars$Collar.ID_f == i,]
 s[i,1,1:nTelemetry[i]] <- me$x/1e6
 s[i,2,1:nTelemetry[i]] <- me$y/1e6
}
xmin <- min(s[,1,], na.rm = T)-1 #add some space around locations
xmax <- max(s[,1,], na.rm = T)+1 #add some space around locations
ymin <- min(s[,2,], na.rm = T)-1 #add some space around locations
ymax <- max(s[,2,], na.rm = T)+1 #add some space around locations
```

---
## Random Walks in Nimble

``` r
consts <- list(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax, nind = nind, nTelemetry = nTelemetry)
dat <- list(s = s)
inits <- function(){list(sigma = rgamma(1,1,1))}
params <- 'sigma'
```

``` r
rw_out <- nimbleMCMC(code = r_walk,
 data = dat,
 constants = consts,
 inits = inits(),
 monitors = params,
 thin = 1,
 niter = 5000,
 nburnin = 2500,
 nchains = 3,
 samplesAsCodaMCMC = TRUE
 )
MCMCvis::MCMCtrace(rw_out, pdf = F, Rhat = T, n.eff = T)
```

---
## Random Walks in Nimble

---
## Random Walks in Nimble

What have we learned?

``` r
2*summary(rw_out)$quantiles*1e6
```

```
##  2.5%   25%   50%   75% 97.5% 
## 489.9 493.8 495.5 497.5 500.6
```

Our model tells us that on average, the deer moved around 495 meters from its previous location (remember these are locations taken every 4 hours).

--
 
This model is fine, but not especially useful.

---
## Home Range Estimation

Let's try to fit a different type of movement model so we can learn about the approximate home range of the individual. In this case, we are interested in determining 2 latent parameters - the 'activity center' of the individual, and the average distance an individual might be at any given time from that central point.  This model ignores any temporal correlation between locations.

Note that this model looks very similar to the random walk model, but we really aren't explicitly modeling 'movement.'

``` r
central_walk <- nimbleCode({
 for(i in 1:nind){
 mu[i,1] ~ dunif(xmin, xmax)
 mu[i,2] ~ dunif(ymin, ymax)
 for(t in 1:nTelemetry[i]){
 s[i,1,t] ~ dnorm(mu[i,1], sd = sigma[i])
 s[i,2,t] ~ dnorm(mu[i,2], sd= sigma[i])
 }
 sigma[i] ~ dgamma(1,1)
 }
})
```

---
## Home Range Estimation

``` r
consts <- list(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax, nind = nind, nTelemetry = nTelemetry)
dat <- list(s = s)
```

We will need to give initial values for mu and sigma

``` r
mu.init <- apply(s, c(1,2), function(x) mean(x, na.rm= T))
inits <- function(){list(sigma = rgamma(6,1,1), mu = mu.init)}
params <- c('mu', 'sigma')
```

---
## Home Range Estimation

Here's a quick visual of our initial value for `$mu$` for the first individual relative to the recorded locations.

``` r
plot(s[1,1,], s[1,2,], type = 'b', pch = 19, cex= .5)
points(mu.init[1,1], mu.init[1,2], col = 'red', pch= 19)
```

---
## Home Range Estimation

``` r
central_out <- nimbleMCMC(code = central_walk,
 data = dat,
 constants = consts,
 inits = inits(),
 monitors = params,
 thin = 1,
 niter = 5000,
 nburnin = 2500,
 nchains = 3,
 samplesAsCodaMCMC = TRUE
 )
MCMCvis::MCMCtrace(central_out, params = c('sigma'), pdf = F, Rhat = T, n.eff = T)
```

---
## Home Range Estimation

---
## Home Range Estimation

A fun fact, that you do not need to memorize, is that home range size can be estimated from this model as a circle with radius `$2.45*\sigma$`.

Let's see what that says about our deer and then compare the estimate home range with the raw data.

``` r
sigma <- MCMCvis::MCMCsummary(central_out, params = 'sigma', Rhat = F, n.eff = F)
2.45*sigma*1e3 #radius in km
```

```
##            mean       sd   2.5%    50%  97.5%
## sigma[1] 0.6767 0.009373 0.6579 0.6764 0.6952
## sigma[2] 0.6971 0.010118 0.6769 0.6972 0.7180
## sigma[3] 0.9731 0.012304 0.9493 0.9729 0.9987
## sigma[4] 0.7703 0.010488 0.7499 0.7702 0.7912
## sigma[5] 1.1073 0.013913 1.0802 1.1067 1.1352
## sigma[6] 1.8658 0.023238 1.8223 1.8660 1.9120
```

``` r
pi*(2.45*sigma$mean*1e3)^2 #~home range in km^2
```

```
## [1]  1.439  1.527  2.975  1.864  3.852 10.937
```

Remember that our data is a combination of deer of multiple ages/sexes (and CWD statuses), so the variation here is not entirely unexpected.

---
## Home Range Estimation - Graph

``` r
library(ggforce) #lets us draw circles
mux <- MCMCsummary(central_out, params = 'mu')$`50%`[1:6] 
muy <- MCMCsummary(central_out, params = 'mu')$`50%`[7:12] 
mus <- cbind(mux, muy)

collars$Collar.ID <- as.character(collars$Collar.ID)
ggplot(collars)+
 geom_line(dat = collars, aes(x = x/1e6, y= y/1e6, group = Collar.ID, col = Collar.ID))+
 geom_point(dat =mus, aes(x = mux, y =muy), cex = .5)+
 geom_circle(dat = mus, aes(x0 = mux, y0 = muy, r = 2.45*sigma$mean))+
 scale_color_manual(values = rainbow(6))+
 theme_bw()+
 coord_equal()+
 xlab('UTM W (km)')+
 ylab('UTM N (km)')+
 theme(axis.title = element_text(size = 15))
```

---
## Home Range Estimation - Graph

Clearly our model is not suuuuuuper accurate, but it's a good starting point.