marked R package by Laake et al 2013

The package Rmark by Jeff Lakke for interfacing R and Program MARK has been extended to be a stand-alone mark-recapture package called marked.

The package is detailed in the paper Laake et al. 2013. marked: an R package for maximumlikelihood and Markov Chain Monte Carlo analysis of capture–recapture data. Methods in Ecology & Evolution 4, 885–890

Package marked appears to function similarly to Rmarked but uses Markov-chain Monte Carlo methods and AD Model Builder to do its calculations instead of calling Program MARK. It appears to also further integrate R functions into mark-recapture analysis.

The Laake et al. (2013) article provides some tutorials using the classic European Dipper dataset. However, it for the sake of brevity it skips some code that is avaiable through a vignette for the package.

I have not located the vignette but I think similar code is available in the documentation for the original RMark package in the Program MARK: AGentle Introduction text by Cooch and White. Here I get the examples in Laake et al. (2013) – hereafter LJC (2013) after the authors initials – up and running.

Some preliminaries

# Load the package
library(marked)
## This is marked 1.1.3

# Load the classic dipper data
data(dipper)

The code in LJC (2013) uses individual, time-varying covariates to model trap dependence. However, the original dipper data only contains covariates for sex.

# The stock dipper dataframe
summary(dipper)
##       ch                sex     
##  Length:294         Female:153  
##  Class :character   Male  :141  
##  Mode  :character

Laake & Rexstad (2014) provide code for generating individual time-varying covariates for the dipper data in “Appendix C: RMark - an alternative approach tobuilding linear models in MARK” in the MARK book. This code can be found on page 73 of the 2014 version of the appendix, in the section “C.16. Individual covariates.”

Below is their code, modified to feed directly into package marked instead of Rmark, with some additional annotations by me.

First LR (2014) create a single weight covariate. Then they create a time-varying covariate for trap-dependence. The LJC (2013) article also uses a factor varible for flood, which as of the first version of this script I have not created. The information regrading this covariate and the code for generating it I believe are also in LR (2014).

Original Code from LJC (2013)

Using the stock dipper dataframe this code results in an error because there is no “Flood” or trap-dependence (“td”) covariates in the original data.


### Process data ###

# Process dipper data with process.data() function from marked package
dipper.proc <- process.data(dipper, model = "cjs", begin.time = 1981)
## 255 capture histories collapsed into  53

# Create design data with static and time-varying covariates
design.Phi = list(static = c("weight"), time.varying = c("Flood"))
design.p = list(static = c("sex"), time.varying = c("td"))

design.parameters <- list(Phi = design.Phi, p = design.p)

# This results in an error because there are not any covariates!
ddl <- make.design.data(dipper.proc, parameters = design.parameters)
## Error: Missing time varying variable
## Flood1981,Flood1982,Flood1983,Flood1984,Flood1985,Flood1986

Simulate Individual “Weight” Covariate

This code from LR (2014) simulates a single individual covariate, weight, with mean = 10 and SD = 2.

# First, create a single 'weight' column create some imaginary weight data
# using a random normal.
dipper$weight <- rnorm(294, 10, 2)

Create time-varying Individual Covariates

This code creates a “trap dependence” (“td”) individual covariates. If a bird was captured the previous year it gets coded as “1”. If it was not capture every before, or its been more than 1 years since capture, it gets a “0”.


# Code from page 76 of Appendix C

### Create a function for processing the data ###
create.td <- function(ch, varname = "td", begin.time = 1) 
# Arguments: ch = capture history vector (0/1 values only) varname = prefix
# for variable name begin.time = time for first occasion Value: returns a
# dataframe with trap-dependent variables named
# varnamet+1,...,varnamet+(nocc-1) where t is begin.time and nocc is the
# number of occasions

# Function starts here:
{
    # turn vector of capture history strings into a vector of characters
    char.vec <- unlist(strsplit(ch, ""))

    # test to make sure they only contain 0 or 1
    if (!all(char.vec %in% c(0, 1))) 
        stop("Function only valid for CJS model without missing values") else {
        # get number of occasions (nocc) and change it into a matrix of numbers
        nocc <- nchar(ch[1])
        tdmat <- matrix(as.numeric(char.vec), ncol = nocc, byrow = TRUE)

        ### remove the >> last column << of individual covariates, which is not used
        ### ### NB: this is important. There are always (nocc - 1) individual
        ### covariates.  This is because the #covaraiates relate specifically to the
        ### transitions BETWEEN capture occassions.  So if a dipper study was
        ### conducted from 2003 to 2006, nocc = 4 and ther would be 3 'td' individual
        ### covariates, td2004, td2005,td2006.  The covariate for 2006 is not
        ### releavant because that data would only be useful for the non-existent
        ### transition from 2006 to 2007.
        tdmat <- tdmat[, 1:(nocc - 1)]

        # turn it into a dataframe and assign the field (column) names
        tdmat <- as.data.frame(tdmat)
        names(tdmat) <- paste(varname, (begin.time + 1):(begin.time + nocc - 
            1), sep = "")
        return(tdmat)
    }
}

Test the function create.td() function scipted above. Note that you can't just give it the dipper objectbut have to specifically give it dipper$ch

dipper.td <- create.td(dipper$ch, begin.time = 1981)  #This must be properly set for everything else to run. 

The dipper.td object contains just the trap-dependence dummy variables

head(dipper.td)
##   td1982 td1983 td1984 td1985 td1986 td1987
## 1      0      0      0      0      0      0
## 2      0      0      0      0      0      0
## 3      0      0      0      0      0      0
## 4      0      0      0      0      0      0
## 5      0      0      0      0      0      0
## 6      0      0      0      0      0      0

I'll add the dipper.td object to the original dipper object

dipper2 <- cbind(dipper, dipper.td)

head(dipper2)
##        ch    sex weight td1982 td1983 td1984 td1985 td1986 td1987
## 1 0000001 Female  8.105      0      0      0      0      0      0
## 2 0000001 Female  9.605      0      0      0      0      0      0
## 3 0000001 Female  9.559      0      0      0      0      0      0
## 4 0000001 Female 12.426      0      0      0      0      0      0
## 5 0000001 Female 10.732      0      0      0      0      0      0
## 6 0000001 Female  8.648      0      0      0      0      0      0

Some notes on the structure of the individual covariate data. The “begin.time” for these data is set as 1981 and there the number of occassion (“nocc”) is 7. So, the study ran from 1981 until 1987. The “td” columns in the dipper2 dataframe I've created are named starting with td1982 and run only through td1987.

The covariates are structured and named this way because the td1982 covariate is used in the calculation of the parameters for survival from 1981 to 1982 and the probablity of being observed in 1982. I get a little confused about this naming convention when, for example, the individual covariate is something like length. If these covariates were the mass of the bird, the covariate would have been measured in 1981 when the bird was caught but the covariate would be named “wt1982” because it relates to the survival and observation probablities in 1982.

###Running LJC (2013) examples with individual covariates

Now we have a simluated weight and trap dependece covariates. I'll re-run the data from LJC (2013) again. In this case, I've removed the code for the “flood” covariate by commenting it out


### Process data ###

# Process dipper data with process.data() function from marked package
dipper.proc <- process.data(dipper2, #my new dipper2 object
                            model='cjs', 
                            begin.time=1981)
## Error: replacement has 254 rows, data has 255

#Create design data with static and time-varying covariates
design.Phi=list(static=c('weight')#,
                #time.varying=c('Flood') #flood removed
                )

design.p=list(static=c('sex'),
              time.varying=c('td')) #td

design.parameters <- list(Phi=design.Phi, p=design.p)

#The make.design.data() function from the marked package now runs properly
ddl <- make.design.data(dipper.proc, 
                        parameters=design.parameters)
## Error: undefined columns selected

Run models found in Methods in Ecology & Evolution Paper

This runs the models found in the paper.


### Define parameters ###
#Define the survival (phi) model
Phi.sfw <- list(formula= ~ weight) #NB: the original code include "flood"

#Define the p (observation) model
p.ast <- list(formula= ~ sex+td)


### Run Models ###

#Call the crm() function to run the model with AD Model builder (admb)
#This requires the R2admb package, which is loaded automaticall, ADMB itself, and a C++ compiler.  See 
#the helpf file for the crm (eg "?crm") for information on setting ADMB up.

# model_admb <- crm(dipper.proc, 
#                   ddl,
#                   hessian=TRUE,
#                   model.parameters = list(Phi=Phi.sfw,p=p.ast), 
#                   use.admb=TRUE)

#Call crm() to run the model with markov-chain monte-carlo (MCMC)
model_mcmc <- crm(dipper.proc,             #There might be typo in original code; the object called is    
                                           # dipper, not dipper.proc
                  model="probitCJS",
                  begin.time=1981,
                  design.parameters=design.parameters,
                  model.parameters=list(Phi=Phi.sfw,p=p.ast),
                  burnin=1000,
                  iter=10000,
                  debug = FALSE)
## Creating design data.
## Error: undefined columns selected

The MCMC model fits almost instantly!

Fitting splines

The authors provide an example fitting polynomial splines to a numeric (continous) time covariate using R splines package.

library(splines)

# define a function for the spline magic
do.example <- function() {

    # Phi models
    Phi.1 <- list(formula = ~1)
    Phi.2 <- list(formula = ~bs(Time))  #'bs' must relate to splines.  'Time' is continous time variable

    # P models
    p.1 <- list(formula = ~1)
    p.2 <- list(formula = ~time)  #'time' with lower-case 't' is normal time, as a factor

    # Use the built-in mark wrapper function create.model.list() (produces
    # copious output)
    crmodel.list <- create.model.list(c("Phi", "p"))
    return(crm.wrapper(crmodel.list, data = dipper2, begin.time = 1981, debug = FALSE, 
        hessian = TRUE, use.admb = FALSE, external = FALSE))  #debug does not seem to supress output
}

example <- do.example()
## Phi.1.p.1
## Model: CJS
## 
## Processing data
## Error: replacement has 254 rows, data has 255

All of the models and their AIC valeus can be called up using model.table

model.table(example)
## Error: object of type 'closure' is not subsettable

Out put of one the models can be accessed using print.crm

print.crm(example[[2]])
## Error: object of type 'closure' is not subsettable

Extract one of the models

ex1 <- example[[2]]
## Error: object of type 'closure' is not subsettable
ex1$results$real$p
## Error: object 'ex1' not found

Package marked Utility functions

# Phi.mean(x,age=0,time=NULL,age.bins=NULL,age.levels=NULL)
# p.mean(x,age=0,time=NULL,age.bins=NULL,age.levels=NULL)
# Phi.boxplot(x,age=0,time=NULL,sex=NULL)
# p.boxplot(x,age=0,time=NULL,sex=NULL)

Additional Rmark code from appendix

Additional code from the appendix for running models with just the weight covariate. This code creates a “mark wrapper” function. In the packaed marked there is the function crm.wrapper() to take care of this.

# #### Create a 'mark wrapper' function to do the analyses. #### #This
# 'wrapper' packages up multiple models and processes them together for
# comparison and model- #selection purposes
# do.dipper.covariate.example=function() { # process the data for CJS model
# and make default design data
# dipper.processed=process.data(dipper,begin.time=1980,groups='sex')
# dipper.ddl=make.design.data(dipper.processed) # create additional Phi
# fields adult and young dipper.ddl$Phi$adult=0
# dipper.ddl$Phi$adult[dipper.ddl$Phi$Age>=1]=1 dipper.ddl$Phi$young=1
# dipper.ddl$Phi$adult ### Define several models for Phi ### #Intercept
# model Phi.dot=list(formula=~1) # weight only for all survivals
# Phi.weight=list(formula=~weight) # sex and sex dependent slope for weight
# Phi.weight.x.sex=list(formula=~weight*sex) # same intercept for
# male/female with a sex dependent slope for weight
# Phi.weight.sex=list(formula=~weight:sex) # effect of weight only for first
# time post capture; if you exclude the # adult term, then an adult would
# have the intercept survival which would # be the value for weight=0
# Phi.weight.plus.sex=list(formula=~adult + young:weight) ### define models
# for p ### p.dot=list(formula=~1) p.time=list(formula=~time) # create model
# list cml=create.model.list('CJS') # run and return models
# return(mark.wrapper(cml,data=dipper.processed,ddl=dipper.ddl)) }
# #Commented out: running the function:
# covariate.results=do.dipper.covariate.example()

Additional Rmark code from appendix for creating a mark wrapper function to run model with individual covariates. In the packaed marked there is the function crm.wrapper() to take care of this.


# ### Create a mark-wrapper function ### # This function 1)generates the
# trap dependence (td) covaraites and 2)runs the models
# do.dipper.td=function() { # get data and add the td time varying
# covariate, process the data # and create the design data dipper <-
# cbind(dipper,create.td(dipper$ch,begin.time=1980)) dipper.processed <-
# process.data(dipper,begin.time=1980) dipper.ddl <-
# make.design.data(dipper.processed) # create additional p field adult
# dipper.ddl$p$adult <- 0 dipper.ddl$p$adult[dipper.ddl$p$Age > 1] <- 1 #
# define models for Phi Phi.dot=list(formula=~1) # define models for p p.td
# <- list(formula=~td) p.td.adult <- list(formula=~td:adult) p.td.time <-
# list(formula=~td:time) p.time.plus.td <- list(formula=~time+td) # create
# model list cml <- create.model.list('CJS') # run and return models
# return(mark.wrapper(cml,data=dipper.processed,ddl=dipper.ddl)) }

td.results <- do.dipper.td()
## Error: could not find function "do.dipper.td"