Lecture 19

class: center, middle, inverse, title-slide

.title[
# Lecture 19
]
.subtitle[
## Open Population Mark Recapture Part 1
]
.author[
### Spring 2025
]

---

## Readings

> Analysis of Capture-Recapture Data By Rachel S. McCrea, Byron J. T. Morgan

> Kéry & Schaub Chapter 10

---
## Closed Capture-mark-recapture

#### Key to closed CMR is that because individuals do not enter or leave the population, we assume any 0 in the capture history is non-detection `$\large (1-p)$`

- in simplest model, easy to translate capture history into probabilistic statements based only on `$p$`:

.pull-left[
##### Individual 1: `$\large 101101$`  
##### Individual 2: `$\large 011001$`  
]

.pull-right[
#### `$\large p(1-p)pp(1-p)p$`  
#### `$\large (1-p)pp(1-p)(1-p)p$`
]

---
## Closed CMR models

#### In all closed CMR models, `$\large N$` is unknown

`$$\Large z_i \sim Bernoulli(\psi)$$`

`$$\Large y_{ik} \sim Bernoulli(z_ip)$$`

`$$\Large \psi \sim beta(1,1)$$`

`$$\Large p \sim beta(1,1)$$`

`$$\Large N = \sum_{i=1}^M z_i$$`

---
## Open Capture-mark-recapture

#### In open populations, we **cannot** assume any 0 in the capture history is non-detection.

Consider the following capture history

`$$\large 101$$`

+ In a closed model, `$\large p(1-p)p$` 
 
 
 
 + In an open model: `$\large (p_1*\psi_1)*((1-p_2)*\psi_2 + (1-\psi_2))*(p_3*\psi_3)$`

If we want to model `$\large \psi$` and `$\large p$` as changing over time, we will need a *robust design* framework.

---
## Robust Design

You have seen the robust design before, we just didn't call it anything fancy.

There is no disadvantage to using a robust design other than the cost of sampling, and most modern studies use this framework.

---
## Robust Design

CJS models were first used in the late 60's/early 70's to look at survival across time. The robust design didn't appear until the 80s.

If observations only happen once per primary period, CJS models cannot not be modified for heterogeneity or permanent trap response
among individuals. Remember the assumptions of the original CJS models:

+ Every animal has the same chance of capture, p
  
  + Every animal has same probability of surviving `$\large \phi$`

Closed CJS models are robust to this un-modeled variation; Open models are *NOT*

In his original papers, the likelihood was estimated ah-hoc, by first estimating each primary period as a closed population model and then separately estimating the likelihood of moving from one time period to the next.

We no longer estimate the likelihood in this exact way, but we have 40+ years of technology and research on our side :)

---
## Robust Design

Advantages of the Robust Design:

+ Estimates of `$p$` (and thus, `$N$` and `$\gamma$`) are less biased by heterogeneity in capture probability

+ Temporary emigration can be estimated (in most cases)

+ If there is no temporary emigration, can estimate `$N$`, `$\phi$` and `$\gamma$` in all time periods (not possible with classic JS models)

+ Better precision on estimates

---
## The Open Population
<img src="figs/flow.png" width="75%" style="display: block; margin: auto;" />

---
## Assumptions of the Open C-R Models

+ Tagged animals are a random sample of the population of interest

+ Numbers of releases are known

+ Tagging is accurate, no tag loss, no misread tags, no data entry errors

+ The fate of individual animals and the fates of animals in differing cohorts are independent

+ All animals, both marked and unmarked have the *same* capture process/probability**

** unless you use the robust design

---
## Open C-R Models - Restricted Dynamic Occupancy
#### Very similar to open occupancy models, except you can't come back from the dead
`$$\Large z_{i1} \sim Bernoulli(\psi)$$`
`$$\Large z_{i,t} \sim Bernoulli(z_{i,t-1} * \phi+(1-\prod_{t= 1}^{t-1} (1- z_{it}))\gamma)$$`

`$$\Large y_{it} \sim Binomial(z_{it}p, J)$$`

`$$\Large \psi \sim beta(1,1)$$`

`$$\Large N_t = \sum_{i=1}^M z_{it}$$`
---
## Open C-R Models - Restricted Dynamic Occupancy

What does this mean in the equation?

`$$\Large (1-\prod_{t= 1}^{t-1} (1- z_{it}))\gamma$$`

In these models, `$z_{it}$` = 0 means one of two things:

+ This individual was previously alive and has died
    
  + This individual is availble to enter the population but hasn't entered it yet

We only want individuals that have not yet been born (or immigrated in) to enter the population. No zombies.

---
## Open C-R Models - Restricted Dynamic Occupancy

One problem with this formulation is that we only have a fixed number of individuals (`$M$`) in our augmented population  (that's why we augmented in the first place).

In time period 1, the expected number of newly arriving individuals is `$M\gamma_1$` (everyone is new). In time period 2, it is `$M(1-\gamma_1)\gamma_2$`.  More generally, it is `$M \prod_{i=1}^{t-1}(1-\gamma_i)\gamma_t$`.

Thus even if per capita recruitment is stable over time, the expected number of recruits will tend to **decrease** and `$\gamma$` becomes a biologically meaningless number.

This is not necessarily a problem and can be partially avoided by allowing `$\gamma$` to change over time.

---
## Open C-R Models - Restricted Dynamic Occupancy

Imagine a population that start with 20 individuals and `$\gamma$` = .2. Our augmented population brings us up to 200 individuals. Survival does not impact recruitment.

```
##       Avail Not Avail Recruits Recruitment Probability
##  [1,]   180        20       20                  0.2000
##  [2,]   144        56       36                  0.2500
##  [3,]   115        85       29                  0.2522
##  [4,]    92       108       23                  0.2500
##  [5,]    74       126       18                  0.2432
##  [6,]    59       141       15                  0.2542
##  [7,]    47       153       12                  0.2553
##  [8,]    38       162        9                  0.2368
##  [9,]    30       170        8                  0.2667
## [10,]    24       176        6                  0.2500
```

---
## Open C-R Models - Restricted Dynamic Occupancy

If the number of recruits is stable over time (rather than `$\gamma$`), how does entry probability have to change to accommodate that?

```
##       Avail Not Avail Recruits Entry Probability
##  [1,]   180        20       20            0.2000
##  [2,]   160        40       20            0.1250
##  [3,]   140        60       20            0.1429
##  [4,]   120        80       20            0.1667
##  [5,]   100       100       20            0.2000
##  [6,]    80       120       20            0.2500
##  [7,]    60       140       20            0.3333
##  [8,]    40       160       20            0.5000
##  [9,]    20       180       20            1.0000
## [10,]     0       200       20               Inf
```
---
## Parameters of Interest

Lots of parameters of potential interest in these models

- `$\phi$` = probability of survival

- `$\gamma$` = biologically meaningless parameter indicating if you are 'removed' from the available population

- `$N_t$` = realized abundance in time `$t$`

- `$N_s$` = total number of individuals ever alive (superpopulation)

- `$B_t$` = realized recruits in time `$t$`

- `$b_t$` = probability of being recruited in time `$t$` given that you ARE recruited at some point from time `$t = 1$` to time `$t = T$` (semi-meaningless biologically)

- `$\psi$` = proportion of the total augmented population (real + imaginary) that are real

- `$f = \frac{B_t}{N_t}$` = per capita realized recruitment

---
## Restricted Dynamic Occupancy Pros and Cons

Pros:

+ Relatively fast 
  
  + Fairly simple code (all things considered)
  
  + Can be marginalized to remove `$z$`

Cons:

+ Gamma parameter is biologically meaningless 
  
  + Can't put a prior on the superpopulation (total ever alive)
  
  + Can't put a prior on fecundity 
  
  + Recruitment probability remains a black box

---
## Restricted Dynamic Occupancy in Nimble

State Process

``` r
openCR1 <- nimbleCode({
 for(t in 1:nocc){
 gamma[t] ~ dbeta(1,1) #some probability of entry
 }
 phi ~ dbeta(1,1)
 for(i in 1:M){
 z[i,1] ~ dbern(gamma[1]) #alive/real in time 1
 a[i,1] <- 1-z[i,1] #available for recruitment

for(t in 2:nocc){
 phi.p[i,t] <- phi*z[i,t-1] + a[i,t-1]*gamma[t] #effective survival or recruitment 
 a[i,t] <- (1-z[i,t])*a[i,t-1] #available for recruitment
 z[i,t] ~ dbern(phi.p[i,t]) #survival/recruitment determine status
 } #end t
 }# end i
```

---
## Restricted Dynamic Occupancy in Nimble

Detection Process

``` r
  ##----- DETECTION PROCESS -----##
  for(t in 1:nocc){ #each year
    p[t] ~ dbeta(1,1) 
    for(i in 1:M){
      y[i,t] ~ dbinom(size = ntrials[t], p = p[t]*z[i,t])
  }#i
  } #t
```

---
## Restricted Dynamic Occupancy in Nimble

Derived Parameters

``` r
 ### Derived parameters
 for(t in 1:nocc){
 N[t] <- sum(z[1:M, t]) #realized abundance in time t
 A[t] <- sum(a[1:M, t]) #available for recruitment in time t
 qgamma[t] <- 1-gamma[t] #not entering this time period
 } #end t
 cprob[1] <- gamma[1] #expected recruitment probability 
 for(t in 2:nocc){
 cprob[t] <- gamma[t]*prod(qgamma[1:(t-1)]) #expected recruitment prob
 B[t] <- A[t-1]-A[t] #how many actually got recruited
 f[t] <- B[t]/N[t] #realized per-capita recruitment
 }
 psi <- sum(cprob[1:nocc]) #overall inclusion probability of M
 for(t in 1:nocc){
 b[t] <- cprob[t]/psi #entry probability, given I entered at some point
 }
 Nsuper <- M-A[nocc] #Superpopulation (ever alive)
})
```

---
## Open C-R Models - Per Capita

Let `$\Large A_{t}$` represent the sum of all the individuals available for recruitment in time `$t$`. Then instead of introducing recruitment via:

`$$(1-\prod_{t= 1}^{t-1} (1- z_{it}))\gamma$$`

We use:

`$$(1-\prod_{t= 1}^{t-1} (1- z_{it}))\gamma_t'$$`
where 
`$$\gamma'_t = \frac{N_{t-1}}{A_{t-1}}\gamma$$`
Thus we are weighting the value of `$\gamma$` by the number of individuals that are alive in the last time period and the total number still available to be recruited.

---
## Open C-R Models - Per Capita

Imagine that same population that started with 20 individuals and `$\gamma$` = .2. Our augmented population brings us up to 200 individuals. Let `$\phi$` = 0.95.

```
##       Avail Not Avail Recruits Now Alive Entry Probability
##  [1,]   180        20       20        20           0.20000
##  [2,]   178        22        2        18           0.02222
##  [3,]   174        26        4        17           0.02022
##  [4,]   169        31        5        16           0.01954
##  [5,]   167        33        2        15           0.01893
##  [6,]   164        36        3        15           0.01796
##  [7,]   163        37        1        15           0.01829
##  [8,]   162        38        1        15           0.01840
##  [9,]   160        40        2        14           0.01852
## [10,]   156        44        4        13           0.01750
```

---
## Open C-R Models - Per Capita in Nimble

In Nimble, we only need to change a few lines:

``` r
openCR2 <- nimbleCode({
 for(t in 1:nocc){
 gamma[t] ~ dbeta(1,1) #some probability of entry
 }
 phi ~ dbeta(1,1)
 for(i in 1:M){
 z[i,1] ~ dbern(gamma[1]) #alive/real in time 1
 a[i,1] <- 1-z[i,1] #available for recruitment
 } #end i
 for(t in 2:nocc){
 gammaPrime[t] <- N[t-1]*gamma[t-1]/A[t-1] #weight gamma by number alive in previous t
 for(i in 1:M){
 phi.p[i,t] <- phi*z[i,t-1] + a[i,t-1]*gammaPrime[t] #effective survival or recruitment 
 a[i,t] <- (1-z[i,t])*a[i,t-1] #available for recruitment
 z[i,t] ~ dbern(phi.p[i,t]) #survival/recruitment determine status
 } #end t
 }# end i
```

---
## Open C-R Models - Per Capita in Nimble

Detection Process (doesn't change)

``` r
  ##----- DETECTION PROCESS -----##
  for(t in 1:nocc){ #each year
    p[t] ~ dbeta(1,1) 
    for(i in 1:M){
      y[i,t] ~ dbinom(size = ntrials[t], p = p[t]*z[i,t])
  }#i
  } #t
```

---
## Open C-R Models - Per Capita in Nimble
Derived Parameters (small change)

``` r
 ### Derived parameters
 for(t in 1:nocc){
 N[t] <- sum(z[1:M, t]) #realized abundance in time t
 A[t] <- sum(a[1:M, t]) #available for recruitment in time t
 } #end t
 cprob[1] <- gamma[1] #expected recruitment probability 
 qgamma[1] <- 1 - gamma[1]
 for(t in 2:nocc){
 qgamma[t] <- 1-gammaPrime[t] #not entering this time period
 cprob[t] <- gammaPrime[t]*prod(qgamma[1:(t-1)]) #expected recruitment prob
 B[t] <- A[t-1]-A[t] #how many actually got recruited
 f[t] <- B[t]/N[t] #realized per-capita recruitment
 }
 psi <- sum(cprob[1:nocc]) #overall inclusion probability of M
 for(t in 1:nocc){
 b[t] <- cprob[t]/psi #entry probability, given I entered at some point
 }
 Nsuper <- M-A[nocc] #Superpopulation (ever alive)
})
```

---
## Parameters of Interest

To review:

- `$\phi$` = probability of survival

- `$\gamma$` = biologically meaningless parameter indicating if you are 'removed' from the available population

- `$N_t$` = realized abundance in time `$t$`

- `$N_s$` = total number of individuals ever alive (superpopulation)

- `$B_t$` = realized recruits in time `$t$`

- `$b_t$` = probability of being recruited in time `$t$` given that you ARE recruited at some point from time `$t = 1$` to time `$t = T$`

- `$\psi$` = proportion of the total augmented population (real + imaginary) that are real

- `$f = \frac{B_t}{N_t}$` = per capita realized recruitment

---
## All parameters are connected