day2-2-behavior.Rmd

---
title: "Behavioral Complexities"
author: "Peter Solymos <solymos@ualberta.ca>"
---

## Introduction

We have reviewed so far how to fit _naive_ models to estimate the expected value of the observed counts, $\lambda$. So what is this $\lambda$? Here are some definitions for further discussion:

- **relative abundance**: $\lambda$ without any reference to nuisance variables, but possibly standardized by design, or nuisance variables used as fixed effects,
- **abundance**: $N=\lambda/C$, $C$ is a correction factor and $N$ refers to the number of individuals within the area surveyed -- the problem is that we cannot measure this directly (this is a latent variable), moreover the survey area is also often unknown (i.e. for unlimited distance counts),
- **occupancy**: the probability that the survey area is occupied, this is really equivalent to the indicator function $N>0$,
- **density** $D = N/A = \lambda/AC$, abundance per unit area -- same problems as above: $N$ and sometimes $A$ are unknowns.

Our objective in the following chapters is to work out the details of estimating abundance and density in some clever ways through learning about the nature of the mechanisms contributing to $C$.

## Prerequisites

```{r beh-libs,message=TRUE,warning=FALSE}
## mace a local copy of day 2 files
source("src/functions.R")
qpad_local(day=2)

## update bSims - some issues fixed
#remotes::install_github("psolymos/bSims")

library(bSims)              # simulations
library(detect)             # multinomial models
load("data/josm-data.rda")  # JOSM data
set.seed(1)
```

Just spend a bit of time admiring the package startup messages... 

## Simulations

The conditionally independent _layers_ of a **bSims** realization are the following, with the corresponding function:

1. landscape (`bsims_init`),
2. population (`bsims_populate`),
3. behavior with movement and vocalization events (`bsims_animate`),
4. the physical side of the observation process (`bsims_detect`), and
5. the human aspect of the observation process (`bsims_transcribe`).

See this example as a sneak peek that we'll explain in the following subsections:

```{r intro,fig.width=12,fig.height=8,out.width='100%'}
phi <- 0.5                 # singing rate
tau <- 1:3                 # detection distances by strata
tbr <- c(3, 5, 10)         # time intervals
rbr <- c(0.5, 1, 1.5)      # count radii

l <- bsims_init(extent=10, # landscape
  road=0.25, edge=0.5)
p <- bsims_populate(l,     # population
  density=c(1, 1, 0))
e <- bsims_animate(p,      # events
  vocal_rate=phi,
  move_rate=1, movement=0.2)
d <- bsims_detect(e,       # detections
  tau=tau)
x <- bsims_transcribe(d,   # transcription
  tint=tbr, rint=rbr)

get_table(x, "removal") # removal table
get_table(x, "visits")  # visits table

op <- par(mfrow=c(2,3), cex.main=2)
plot(l, main="Initialize")
plot(p, main="Populate")
plot(e, main="Animate")
plot(d, main="Detect")
plot(x, main="Transcribe")
par(op)
```

The layers allow us to fix some of the layers and simulate multiple realizations conditional on these fixed layers, e.g. fix the landscape and population layers and only change the behavior, etc.

### Landscape

The `bsims_ini` function sets up the geometry of a local landscape. The `extent` of the landscape determines the edge lengths of a square shaped area. With no argument values passed, the function assumes a homogeneous _habitat_ (H) in a 10 units x 10 units landscape, 1 unit is 100 meters. Having units this way allows
easier conversion to ha as area unit that is often used in the North American bird literature. As a result, our landscape has an area of 1 km$^2$.

The `road` argument defines the half-width of the road that is placed in a vertical position. The `edge` argument defines the width of the edge stratum on both sides of the road. Habitat (H), edge (E), and road (R) defines the 3 strata that we refer to by their initials (H for no stratification, HER for all 3 strata present).

The origin of the Cartesian coordinate system inside the landscape is centered at the middle of the square. The `offset` argument allows the road and edge strata to be shifted to the left (negative values) or to the right (positive values) of the horizontal axis. This makes it possible to create landscapes with only
two strata. The `bsims_init` function returns a landscape object (with class 'bsims_landscape')

```{r landscape,fig.width=10,fig.height=10}
(l1 <- bsims_init(extent = 10, road = 0, edge = 0, offset = 0))
(l2 <- bsims_init(extent = 10, road = 1, edge = 0, offset = 0))
(l3 <- bsims_init(extent = 10, road = 0.5, edge = 1, offset = 2))
(l4 <- bsims_init(extent = 10, road = 0, edge = 5, offset = 5))

op <- par(mfrow = c(2, 2))
plot(l1, main = "Habitat")
points(0, 0, pch=3)
plot(l2, main = "Habitat & road")
lines(c(0, 0), c(-5, 5), lty=2)
plot(l3, main = "Habitat, edge, road + offset")
arrows(0, 0, 2, 0, 0.1, 20)
lines(c(2, 2), c(-5, 5), lty=2)
points(0, 0, pch=3)
plot(l4, main = "2 habitats")
arrows(0, 0, 5, 0, 0.1, 20)
lines(c(5, 5), c(-5, 5), lty=2)
points(0, 0, pch=3)
par(op)
```


### Population

The `bsims_populate` function _populates_ the landscape we created by the `bsims_init` function, which is the first argument we have to pass to `bsims_populate`. The function returns a population object (with class 'bsims_population'). The most important argument that controls how many individuals will inhabit our landscape is `density` that defines the expected value of individuals per unit area (1 ha).

By default, `density = 1` ($D=1$) and we have 100 ha in the landscape ($A=100$) which translates into 100 individuals on average ($E[N]=\lambda=AD$).
The actual number of individuals in the landscape might deviate from this expectation, because $N$ is a random variable ($N \sim f(\lambda)$). The `abund_fun` argument controls this relationship between the expected ($\lambda$) and realized abundance ($N$). The default is a Poisson distribution:

```{r pop-pois}
bsims_populate(l1)
```

Changing `abund_fun` can be useful to make abundance constant or allow under- or over-dispersion, e.g.:

```{r pop-nb}
summary(rpois(100, 100)) # Poisson variation
summary(MASS::rnegbin(100, 100, 0.8)) # NegBin variation
negbin <- function(lambda, ...) MASS::rnegbin(1, lambda, ...)
bsims_populate(l1, abund_fun = negbin, theta = 0.8)
## constant abundance
bsims_populate(l1, abund_fun = function(lambda, ...) lambda)
```

Once we determine how many individuals will populate the landscape, we have control over the spatial arrangement of the nest location for each individual. The default is a homogeneous Poisson point process (complete spatial randomness). 
Deviations from this can be controlled by the `xy_fun`. This function takes distance as its only argument and returns a numeric value between 0 and 1. A function `function(d) reurn(1)` would be equivalent with the Poisson process, meaning that every new random location is accepted with probability 1 irrespective of the distance between the new location and the previously generated point locations in the landscape.

When this function varies with distance, it leads to a non-homogeneous point process via this accept-reject algorithm. The other arguments (`margin`, `maxit`, `fail`) are passed to the underlying `accepreject` function to remove edge effects and handle high rejection rates.

In the next example, we fix the abundance to be constant (i.e. not a random variable, $N=\lambda$) and with different spatial point processes:

```{r pop-xy,fig.height=9,fig.width=6}
D <- 0.5
f_abund <- function(lambda, ...) lambda

## systematic
f_syst <- function(d)
  (1-exp(-d^2/1^2) + dlnorm(d, 2)/dlnorm(exp(2-1),2)) / 2
## clustered
f_clust <- function(d)
  exp(-d^2/1^2) + 0.5*(1-exp(-d^2/4^2))

p1 <- bsims_populate(l1, density = D, abund_fun = f_abund)
p2 <- bsims_populate(l1, density = D, abund_fun = f_abund, xy_fun = f_syst)
p3 <- bsims_populate(l1, density = D, abund_fun = f_abund, xy_fun = f_clust)

distance <- seq(0,10,0.01)
op <- par(mfrow = c(3, 2))
plot(distance, rep(1, length(distance)), type="l", ylim = c(0, 1), 
  main = "random", ylab=expression(f(d)), col=2)
plot(p1)

plot(distance, f_syst(distance), type="l", ylim = c(0, 1), 
  main = "systematic", ylab=expression(f(d)), col=2)
plot(p2)

plot(distance, f_clust(distance), type="l", ylim = c(0, 1), 
  main = "clustered", ylab=expression(f(d)), col=2)
plot(p3)
par(op)
```

The `get_nests` function extracts the nest locations. `get_abundance` and `get_density` gives the total abundance ($N$) and density ($D=N/A$, where $A$ is `extent^2`) in the landscape, respectively.

If the landscape is stratified, that has no effect on density unless we specify different values through the `density` argument as a vector of length 3 referring to the HER strata:

```{r pop-dens,fig.width=10,fig.height=10,out.width='100%'}
D <- c(H = 2, E = 0.5, R = 0)

op <- par(mfrow = c(2, 2))
plot(bsims_populate(l1, density = D), main = "Habitat")
plot(bsims_populate(l2, density = D), main = "Habitat & road")
plot(bsims_populate(l3, density = D), main = "Habitat, edge, road + offset")
plot(bsims_populate(l4, density = D), main = "2 habitats")
par(op)
```

But birds don't just stay put in one place and do nothing. They move and vocalize. The `bsims_animate` function _animates_ the population created by the `bsims_populate` function. `bsims_animate` returns an events object (with class 'bsims_events'). The most important arguments are governing the `duration` of the simulation in minutes, the vocalization (`vocal_rate`), and the movement (`move_rate`) rates as average number of events per minute.

`bsims_animate` uses independent Exponential distributions with rates `vocal_rate` and `move_rate` to simulate vocalization and movement events, respectively.

```{r beh-events,fig.height=4,fig.width=6}
l <- bsims_init()
p <- bsims_populate(l, density = 0.5)
e1 <- bsims_animate(p, vocal_rate = 1)
```

There are no movement related events when `move_rate = 0`, the individuals are always located at the nest, i.e. there is no within territory movement. If we increase the movement rate, we also have to increase the value of `movement`, that is the standard deviation of bivariate Normal kernels centered around each nest location. This kernel is used to simulate new locations for the movement events. Increase the value of `movement` to see how that works.

Movement is illustrated by a line, crosses indicate nest locations, dots are the vocalization events

```{r beh-move,fig.width=10,fig.height=5,out.width='100%'}
e2 <- bsims_animate(p, move_rate = 1, movement = 0.25)

op <- par(mfrow = c(1, 2))
plot(e1, main = "Closure")
plot(e2, main = "Movement")
par(op)
```

## Exercise

Play time!

```{r eval=FALSE}
run_app("bsimsH")
```