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 Jolly-Seber (Cjs) spatial capture-recapture (SCR) models can be used within an R6 object of the JsModel class. The JS model is formulated according to the Schwarz-Arnason parameterisation. I recommend you read the vignettes on the ScrData and ScrModel classes before reading this one. The folder inst/examples/ in the repository contains example code.

library(openpopscr)
library(secr)

data("jsmodel_example")

Simulating data

The openpopscr package uses the ScrData class to store data. See the vignette on the ScrData class to see how to create the data object from your own raw data. In this section, I simulate a JS SCR survey. The simulation function is useful when you want to try out a new method where you know what the correct parameter estimates are, you want to test how robust a method is to violations of its assumptions, or you want to see how data simulated from a model compare with real data.

To simulate a survey, I require true detection, survival parameters, and recruitment parameters, a detector layout, and a mesh. The specified value of beta is the proportion of the super-population (those ever to have lived) to be a live in the first occasion.

# set true parameters 
true_par <- list(lambda0 = 0.2, sigma = 20, phi = 0.8, beta = 0.2, D = 1000)

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

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

# set number of occasions to simulate
n_occasions <- 10

To simulate the data, I set a random seed so that the results will be the same every time the code is run. This is so that you, the reader, can re-run the code and produce the same output as me. In practice, you may not want to set the seed if you want the output from the simulation to change randomly.

scrdat <- simulate_js_openscr(true_par, n_occasions, detectors, mesh, seed = 52381)

You will notice the simulate_js_scr function produces useful output to let the user know what is going on. If you want to switch this off, for example if you are running lots of simulations and do not want this output being produced over and over, then you can set print = FALSE just as I have set the seed above.

In the simulated survey, r scrdat$n() unique individual were detected at least once. More summary information can be found by typing the name of the object into the console:

scrdat

Fitting a model

To create a JsModel, you need to specify a formula for each detection parameter the survival parameter, and the recruitment parameter; also, you need starting values for all parameters. It is best to choose starting values for the parameters that are biologically reasonable and match your data. You could simulate some data using proposed starting values as true parameters to see if they lead to similar data. Poor starting values can lead to spurious convergence in complicated models, models with many parameters. The function get_start_value can be used to produce a reasonable guess for each parameter.

# set each parameter to a be constant 
form <- list(lambda0 ~ 1, 
             sigma  ~ 1,
             phi ~ 1, 
             beta ~ 1)

# get some starting values based on data 
start <- get_start_values(scrdat, model = "JsModel")

# check starting values are reasonable, they seem ok 
start

Before fitting the model, I create the model object:

mod <- JsModel$new(form, scrdat, start)

Once the model object has been created, you can fit the model using the fit function:

mod$fit()

If you have the output printed, the "checking convergence" message reports whether convergence has been reached and if it has not it reports the code output by optim that tells you why it failed to converge.

Reading model output

Once a model has been fit, you can type its name into the console and some summary results are printed:

mod

The parameters are reported on the link scale; this is because when covariates are involved, it is easier to interpret their effects. The package openpopscr uses the log-link function for all of these parameters. So, you would take the exponential of these reported numbers to obtain the parameter values on the response scale.

Alternatively, you can use the get_par function:

mod$get_par("lambda0", k = 1, j = 1)
mod$get_par("sigma", k = 1, j = 1)
mod$get_par("phi", k = 1)
mod$get_par("beta", k = 1)
mod$get_par("D")

I specify $k = 1, j = 1$ in the functions above because by default get_par returns the value of these parameters for every occasion/detector. This is useful when they vary by occasion, but here they are constants, so I only want their value on a single occasion.

If you want to access the reported results, they are stored as a list which is returned by the function estimates:

mod$estimates()

The model AIC can be computed in the usual way:

AIC(mod)

Covariates

The parameters can change with observed covariate values. Covariates can change by detector, by occasion, or by mesh point. At present, individual-level covariates are not implemented in openpopscr. By default, the covariates $x, y, t$ are created to allow for parameters to depend on occasion or space.

As an example, one can fit a model where recruitment depends on time:

# simulate a wet/dry covariate 
set.seed(385638)
season <- as.factor(sample(c("wet", "dry"), size = scrdat$n_occasions(), replace = TRUE))

# add covariate 
scrdat$add_covariate("season", season, "k")
# copy formulae from previous model
form_seas <- form
# update formula for beta to include season 
form_seas[[4]]<- beta ~ season
# create new model object 
mod_seas <- JsModel$new(form_seas, scrdat, start)
# fit model
mod_seas$fit()

# look at results
mod_seas
# predict beta for each season
cbind(season, mod_seas$get_par("beta"))

Covariates that change over time, by detector, and by detector-time can be included in the same way. You can also include spatial covariates on the density parameter $D$. See ScrData and ScrModel vignettes for more details on loading covariates into the ScrData object and fitting models with covariates.

Time and Robust Design

There are two further arguments you can give when creating a ScrData object: time and primary.

The time option requires a vector that contains the start time of each occasion. The time unit is the user's choice. For example, if occasions occur in Weeks 1, 3, 4 of a particular month, this could be specified in a time vector as c(0, 3, 4) (in weeks) or c(0, 14, 21) (in days).

scrdat <- ScrData$new(capthist, mesh, time = time)

Once time is specified, the survival and estimates will be the estimated probability of survival in one time unit, so it will depend on what time unit you used.

The primary argument relates to robust design surveys. In a robust design, there are both primary and secondary occasions for sampling. Secondary occasions take place close together in time and it is assumed the population is closed between them. Population change occurs between primary occasions. The survival parameter between primary occasions is of interest.

For ScrData, you supply a single capthist object (combined across all primary occasions) and then specify in the primary argument which primary each secondary occasion belongs to. For example, suppose a survey has 3 primary occasions with 2, 4, 4 secondary occasions respectively. The capthist object would then have 10 occasions in it and the primary vector would be c(1,1,2,2,2,2,3,3,3,3)

scrdat <- ScrData$new(capthist, mesh, primary = primary)

An example of a robust design analysis is included inst/examples/js_robust_example.R

Transience

For JS models, transience can play an important role. When individuals are transient, their activity can pass in and out of range of detection; individuals that permanently pass out of range of detection are best explained, by the JSModel, as having died and so survival probability is underestimated. For individuals entering the survey area, they are best explained by increased recruitment. Furthermore, transience can cause overestimation of detection probability and substantial underestmation of density within occasion.

You can simulate Js SCR surveys with transience:

# set truth 
true_par <- list(D = 1000, lambda0 = 0.5, sigma = 20, phi = 0.8, beta = 0.2, sd = 20)
# make detectors array 
detectors <- make.grid(nx = 7, ny = 7, spacing = 20, detector = "count")
# make mesh 
mesh <- make.mask(detectors, buffer = 100, nx = 64, ny = 64, type = "trapbuffer")
# set number of occasions to simulate
n_occasions <- 10
# simulate ScrData 
scrtransdat <- simulate_js_openscr(true_par, n_occasions, detectors, mesh, move = TRUE, seed = 34583)

I will first fit a stationary JSModel to these data:

# set formulae and start values 
stat_form <- list(lambda0 ~ 1, sigma ~ 1, phi ~ 1, beta ~ 1)
start <- get_start_values(scrtransdat, model = "JsModel")
# create model object 
stat <- JsModel$new(stat_form, scrtransdat, start)
# fit model 
stat$fit()

# look at results 
stat
# look at parameters on response scale 
stat$get_par("lambda0", k = 1)
stat$get_par("sigma", k = 1)
stat$get_par("phi", k = 1)
stat$get_par("beta", k = 1)
stat$get_par("D")

I now fit the transient model. The transient model class is called JsTransientModel. We must also now specify a start value and formula for the movement parameter.

# specify formulas for parameters 
form <- list(lambda0 ~ 1, 
             sigma  ~ 1,
             phi ~ 1,
             beta ~ 1,
             sd ~ 1)

start <- get_start_values(scrtransdat, model = "JsTransientModel")

trans <- JsTransientModel$new(form, scrtransdat, start)

I now fit the model as usual:

# fit model 
trans$fit()

# look at results 
trans
# look at parameters on response scale 
trans$get_par("lambda0", k = 1)
trans$get_par("sigma", k = 1)
trans$get_par("sd", k = 1)
trans$get_par("phi", k = 1)
trans$get_par("beta", k = 1)
trans$get_par("D")