In r-glennie/openpopscr: Open Population Spatial Capture-Recapture

knitr::opts_chunk$set(
  eval = FALSE, 
  collapse = TRUE,
  fig.width = 7, 
  fig.height = 7, 
  comment = "#>"
)

In this vignette, I describe how other types of model classes (ScrModel, CjsModel, JsTransientModel, etc.) can be fit to multiple strata using the StrataModel class. This is one way to include individual-level effects that are factor variables when these variables are known for all detected individuals (e.g. sex). Each group of individuals is called a stratum. Strata could also be created for different sites or species. The difference between a stratified analysis and analysing these groups separetely is that the stratified model can share some parameters across strata.

library(openpopscr)
library(secr)

data("stratamodel_example")

Simulate data

I will simulate a survey that is conducted over three strata. The idea is strata may share some parameters (i.e., parameters are equal among strata) and for other parameters, each stratum will have its own particular value. Strata are treated independently conditional on the parameters.

Here, all strata have the same value of lambda0, but different sigma and D parameters. I create a ScrData object for each stratum seperately and then store these in a list. I use the same detector array setup and mesh in each stratum, but these could differ.

true_par <- list(s1 = list(D = 1000, lambda0 = 2, sigma = 20), 
                 s2= list(D = 500, lambda0 = 2, sigma = 20), 
                 s3 = list(D = 1000, lambda0 = 2, sigma = 10))

# make detectors array 
detectors <- make.grid(nx = 7, ny = 7, spacing = 20, detector = "multi")

# make mesh 
mesh <- make.mask(detectors, buffer = 100, nx = 64, ny = 64, type = "trapbuffer")

# set number of occasions to simulate
n_occasions <- 5 

# simulate for each stratum
scrdat1 <- simulate_scr(true_par[[1]], n_occasions, detectors, mesh, seed = 15483)
scrdat2 <- simulate_scr(true_par[[2]], n_occasions, detectors, mesh, seed = 79696)
scrdat3 <- simulate_scr(true_par[[3]], n_occasions, detectors, mesh, seed = 43523)

scrdat <- list(s1 = scrdat1, s2 = scrdat2, s3 = scrdat3)

scrdat1
scrdat2
scrdat3

Fitting the model

The StrataModel class works by specifying a formula for each parameter that is to be common among all strata, then separately specifying a formula for each stratum for each parameter that is to be particular to strata. This allows stratum-specific parameters to depend on different covariate combinations.

When creating the StrataModel you must also specify what type of model you wish to fit to each stratum. Here, I fit the ScrModel. You can also specify any other model class available within the openpopscr package.

shared_form <- list(lambda0 ~ 1)

private_form <- list(s1 = list(sigma ~ 1),
                     s2 = list(sigma ~ 1), 
                     s3 = list(sigma ~ 1))


start <- list(lambda0 = 2, sigma = 20, D = 1000)

obj <- StrataModel$new(scrdat, "ScrModel", shared_form, private_form, start)

obj$fit()

Once fitted, the object's name can be typed into the terminal to return the results for each stratum.

obj

Furthermore, the get_object function can be used to extract the model object that corresponds to each stratum. You can then use the functions that model object has to extract estimates or other quantities.

obj$get_object(1)$estimates()

You can extract parameters for each stratum using the get_par function:

# should be the same for all strata
obj$get_par("lambda0", k = 1)
# differ between strata (in simulation strata 1 & 2 have same, 3 differs)
obj$get_par("sigma", k = 1)
# density by stratum
obj$get_par("D")

To compare, let's fit a model where all parameters are shared. This is equivalent to just a single ScrModel. Notice we must create empty lists to specify that no parameters are shared.

shared_form2 <- list(lambda0 ~ 1, sigma ~ 1)

private_form2 <- list(s1 = list(), s2 = list(), s3 = list())

obj2 <- StrataModel$new(scrdat, "ScrModel", shared_form2, private_form2, start)

obj2$fit()