Analysis of lifetime mortality trajectories in wildlife disease research: Basta and beyond

R Code

library(tidyverse)
library(BaSTA)
library(data.table)
library(snowfall)

Load Data

The following code shows the steps to analyse the entire badger data set. In the journal article we first analyse ‘Cub-positive’ and ‘Never positive’ individuals separately. The steps are the same as below.

CH.data<- read.table("Hudson_2019_MTA_Data.txt", header = TRUE)
CH.data<- as.data.frame(CH.data)

Prepare matrices for BaSTA

Create capture history

Using BaSTA function CensusToCaptureHist and specifing ID and date column.

Y.full<- CensusToCaptHist(ID=CH.data[,1], d=CH.data[,2]) 

Create birth and death year matrix

Keep only single entry per badger and create birth and death year matrix.

CH.data<- data.table(CH.data)
CH.data<- distinct(CH.data, CH.data$ID, .keep_all = TRUE)
CH.data <- droplevels(CH.data)
BirthDeath.full<- CH.data[,c(1,8:9)]

Create co-variate matrix

Specify model co-variate as status. (This is a 4 category co-variate: 1. Cub positive female; 2. Cub positive male; 3. Never positive female; 4. Never positive male)

CovMat.full<- MakeCovMat(x=~status, CH.data)

Create input matrix

Combine the birth/death matirx, co-variate matrix and full capture history.

inputMat.aim <- as.data.frame(cbind(BirthDeath.full, Y.full[,-1], CovMat.full[,-1]))

Data check

newData.aim <- DataCheck(inputMat.aim, studyStart = 1,studyEnd = 163, autofix = rep(1, 7),silent = FALSE)

## No problems were detected with the data.
## 
## *DataSummary*
## - Number of individuals         =    2,222 
## - Number with known birth year  =    2,196 
## - Number with known death year  =     336 
## - Number with known birth
##  AND death years                =     326 
## 
## - Total number of detections
##  in recapture matrix            =   10,048 
## 
## - Earliest detection time       =       1 
## - Latest detection time         =     163 
## - Earliest recorded birth year  =       3 
## - Latest recorded birth year    =     161 
## - Earliest recorded death year  =       1 
## - Latest recorded death year    =     164

Begin BaSTA analysis

All Badgers:

Analyse all badgers together specifying all models to compare. Recapture probabilities allowed to vary on each occasion. The following code takes approximately 10 days to run.

All.badgers.out <- multibasta(object = inputMat.aim, studyStart = 2, studyEnd = 173, nsim = 4, niter = 1000000, burnin = 10001, thinning = 100, parallel = TRUE, ncpus = 4, models = c("EX", "GO", "LO"), shapes = c("simple","Makeham", "bathtub"), recaptTrans = 2:173)

Most supported model

The most supported model for both the ‘Cub-positive’ and ‘Never positive’ datum was the Gompertz bathtub or Siler function. The complete data was then analysed together, an example of the summary output from this analysis is shown below (Some lines have been omitted).

summary(All.badgers.out$runs$GO.Bt)

Plot output

View the survival and mortality trajectories together with posterior density distributions of survival parameters.

plot(All.badgers.out$runs$GO.Bt, fancy = TRUE)