Poisson (PNE) mark-resight models in RMark

[jslewis]

Greetings,

I am running the Poisson (PNE) mark-resight models in RMark to estimate abundance. Jeff Laake and Brett McClintock were very helpful in understanding the models and figuring out the R code for RMark -- their assistance was greatly appreciated! They recommended that I post the code for other people that were interested. The following code produces the same results as when the models are run in program MARK. Where you see a "**" in the code you can fill in text/values that are specific to your analyses.

Note that there are some options included below that are specific to a particular model or my personal preferences. For example, the option of simulated annealing is included in the model statement because it was the only way to successfully run some models (and sometimes a particular model had to be run a couple of times to work with this option), but for other models it was unnecessary and the option had to be removed from the model statement for some models to run. Also, for collecting models and viewing the AIC table, Jeff helped with code so that my user-defined model names (e.g., alpha.tsog_sigma.zero) were used in the AIC table (especially useful when parameters were set to 0 and wanted to name models accordingly), but this might not be of interest to everybody.


For the code below, we evaluated 2 covariates "tsog" and "weight", in addition to an "id" label for each animal.

An example input file is:

ch weight tsog id
06 59 0.15 M6
+0 44 0.14 F16
26 70 0.55 M55
04 38 0.09 F72
04 36 0.73 F94
04 54 0.08 M32
02 41 0.25 F3
02 43 0.11 F96
02 65 0.11 M71


# Mark-resight Poisson (PNE) models using RMark
# Jesse Lewis 2013

rm(list=ls())

library(RMark)

# Set working directory -- MARK files will be saved here
# file.choose()
setwd( ** )

# import data file and define field.types as numeric "n" or factor "f"
Input_file = import.chdata( ** ,
header=TRUE, field.names = NULL, field.types = c("n", "n", "f"))
Input_file

n = length(Input_file[,1]) # calculates number of animals in input file
n

pois.proc = process.data(Input_file, model = "PoissonMR",
counts = list("Unmarked Seen" = **,
"Marked Unidentified" = **,
"Known Marks" = n))
pois.ddl = make.design.data(pois.proc)


# To Run Multiple Models

# FIRST: estimate starting values from a(.) sig(.) model

alpha.dot_sig.dot = mark(pois.proc, pois.ddl,
model.parameters = list(alpha = list(formula = ~1),
U = list(formula = ~1),
sigma = list(formula = ~1)),
threads = 2, options = "SIMANNEAL")
alpha.dot_sig.dot$results$derived


# then use these values as initial start values for alpha, U, sig for multiple models below

# define values of alpha, sigma, and U
alpha = **
sigma = **
U = **

# Here are an additional 9 models for examples

alpha.dot_sigma.zero = mark(pois.proc, pois.ddl,
model.parameters = list(alpha = list(formula = ~1),
U = list(formula = ~1),
sigma = list(formula = ~1, fixed = 0)),
initial = c(alpha = alpha, U = U, sigma = 0),
threads = 2, options = "SIMANNEAL")
alpha.dot_sigma.zero$results$derived

alpha.dot_sigma.tsog = mark(pois.proc, pois.ddl,
model.parameters = list(alpha = list(formula = ~1),
U = list(formula = ~1),
sigma = list(formula = ~tsog)),
initial = c(alpha = alpha, U = U, sigma = sigma),
threads = 2, options = "SIMANNEAL")
alpha.dot_sigma.tsog$results$derived

alpha.dot_sigma.weight = mark(pois.proc, pois.ddl,
model.parameters = list(alpha = list(formula = ~1),
U = list(formula = ~1),
sigma = list(formula = ~weight)),
initial = c(alpha = alpha, U = U, sigma = sigma),
threads = 2, options = "SIMANNEAL")
alpha.dot_sigma.weight$results$derived

alpha.dot_sigma.tsog.weight = mark(pois.proc, pois.ddl,
model.parameters = list(alpha = list(formula = ~1),
U = list(formula = ~1),
sigma = list(formula = ~tsog + weight)),
initial = c(alpha = alpha, U = U, sigma = sigma),
threads = 2, options = "SIMANNEAL")
alpha.dot_sigma.tsog.weight$results$derived

alpha.tsog_sigma.dot = mark(pois.proc, pois.ddl,
model.parameters = list(alpha = list(formula = ~tsog),
U = list(formula = ~1),
sigma = list(formula = ~1)),
initial = c(alpha = alpha, U = U, sigma = sigma),
threads = 2, options = "SIMANNEAL")
alpha.tsog_sigma.dot$results$derived

alpha.tsog_sigma.zero = mark(pois.proc, pois.ddl,
model.parameters = list(alpha = list(formula = ~tsog),
U = list(formula = ~1),
sigma = list(formula = ~1, fixed = 0)),
initial = c(alpha = alpha, U = U, sigma = 0),
threads = 2, options = "SIMANNEAL")
alpha.tsog_sigma.zero$results$derived

alpha.tsog_sigma.tsog = mark(pois.proc, pois.ddl,
model.parameters = list(alpha = list(formula = ~tsog),
U = list(formula = ~1),
sigma = list(formula = ~tsog)),
initial = c(alpha = alpha, U = U, sigma = sigma),
threads = 2, options = "SIMANNEAL")
alpha.tsog_sigma.tsog$results$derived

alpha.tsog_sigma.weight = mark(pois.proc, pois.ddl,
model.parameters = list(alpha = list(formula = ~tsog),
U = list(formula = ~1),
sigma = list(formula = ~weight)),
initial = c(alpha = alpha, U = U, sigma = sigma),
threads = 2, options = "SIMANNEAL")
alpha.tsog_sigma.weight$results$derived

alpha.tsog_sigma.tsog.weight = mark(pois.proc, pois.ddl,
model.parameters = list(alpha = list(formula = ~tsog),
U = list(formula = ~1),
sigma = list(formula = ~tsog + weight)),
initial = c(alpha = alpha, U = U, sigma = sigma),
threads = 2, options = "SIMANNEAL")
alpha.tsog_sigma.tsog.weight$results$derived


multiple.model.results = collect.models() # collect all models from above
multiple.model.results$model.table=model.table(multiple.model.results,model.name=F)
# call results to see summary table
options(width = 200) # creates wider window to display results so that they do not wrap in R console
multiple.model.results

# model average N -- derived parameter

model_names = as.vector(as.character(multiple.model.results$model$model)) # get model names
# note that the "model_names" statement can be removed above and below in "N_table" to
# produce a table without model names and quotation marks around values
weights = multiple.model.results$model.table$weight
model_nums = as.numeric(row.names(multiple.model.results$model.table))
N_estimates = sapply(model_nums, function(x) multiple.model.results[[x]]$results$derived$estimate)
N_stderrors = sapply(model_nums, function (x) multiple.model.results[[x]]$results$derived$se)
N_table = cbind(model_nums, model_names, N_estimates, N_stderrors)
N_table

# for model averaged estimate of N and se
N = model.average(x=list(estimate = N_estimates, weight = weights, se = N_stderrors))
N # derived estimate of N with se


Thanks again to Jeff and Brett for their help! Hopefully other people find this code useful for their work. Good luck.

Jesse Lewis