# Last updated September 28, 2019
# The original mallard example shown below uses idividual calls to mark function.
# This is not as efficient as using mark.wrapper and can cause difficulties if different
# groups arguments are used and model averaging is attempted. Below this original example
# the more efficient approach is shown.
# Read in data, which are in a simple text file that
# looks like a MARK input file but (1) with no comments or semicolons and
# (2) with a 1st row that contains column labels
# mallard=read.table("mallard.txt",header=TRUE)
# The mallard data set is also incuded with RMark and can be retrieved with
# data(mallard)
# \donttest{
# This example is excluded from testing to reduce package check time
# ggplot commands have been commented out so RMark doesn't require ggplot
# scripted analysis of mallard nest-survival data in RMark
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
# Example of use of RMark for modeling nest survival data - #
# Mallard nests example #
# The example runs the 9 models that are used in the Nest #
# Survival chapter of the Gentle Introduction to MARK and that#
# appear in Table 3 (page 198) of #
# Rotella, J.J., S. J. Dinsmore, T.L. Shaffer. 2004. #
# Modeling nest-survival data: a comparison of recently #
# developed methods that can be implemented in MARK and SAS. #
# Animal Biodiversity and Conservation 27:187-204. #
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
# The mallard data set is also incuded with RMark and can be retrieved with
data(mallard)
# use the indicator variables for the 4 habitat types to yield
# 1 variable with habitat as a factor with 4 levels that
# can be used for a group variable in RMark
mallard$habitat <- ifelse(mallard$Native == 1, "Native",
ifelse(mallard$Planted == 1, "Planted",
ifelse(mallard$Roadside == 1, "Roadside",
"Wetland")))
# make the new variable a factor
mallard$habitat <- as.factor(mallard$habitat)
mallard.pr <- process.data(mallard,
nocc=90,
model="Nest",
groups=("habitat"))
# Write a function for evaluating a set of competing models
run.mallard <- function()
{
# 1. A model of constant daily survival rate (DSR)
S.Dot = list(formula = ~1)
# 2. DSR varies by habitat type - treats habitats as factors
# and the output provides S-hats for each habitat type
S.Hab = list(formula = ~habitat)
# 3. DSR varies with vegetation thickness (Robel reading)
S.Robel = list(formula = ~Robel)
# 4. DSR varies with the amount of native vegetation in the surrounding area
S.PpnGr = list(formula = ~PpnGrass)
# 5. DSR follows a trend through time
S.TimeTrend = list(formula = ~Time)
# 6. DSR varies with nest age
S.Age = list(formula = ~NestAge)
# 7. DSR varies with nest age & habitat type
S.AgeHab = list(formula = ~NestAge + habitat)
# 8. DSR varies with nest age & vegetation thickness
S.AgeRobel = list(formula = ~NestAge + Robel)
# 9. DSR varies with nest age & amount of native vegetation in
# surrounding area
S.AgePpnGrass = list(formula = ~NestAge + PpnGrass)
# Return model table and list of models
mallard.model.list = create.model.list("Nest")
mallard.results = mark.wrapper(mallard.model.list,
data = mallard.pr,
adjust = FALSE,delete=TRUE)
}
# The next line runs the 9 models above and takes a minute or 2
mallard.results <- run.mallard()
mallard.results
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
# Examine table of model-selection results #
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
# next line exports files; emove # and delete=TRUE in mark.wrapper above
#export.MARK(mallard.results$S.Age$data,
# "MallDSR",
# mallard.results,
# replace = TRUE,
# ind.covariates = "all")
mallard.results # print model-selection table to screen
options(width = 100) # set page width to 100 characters
#sink("results.table.txt") # capture screen output to file
#print(mallard.results) # send output to file
#sink() # return output to screen
# remove "#" on next line to see output in notepad in Windows
# system("notepad results.table.txt",invisible=FALSE,wait=FALSE)
# remove "#" on next line to see output in texteditor editor on Mac
# system("open -t results.table.txt", wait = FALSE)
names(mallard.results)
mallard.results$S.Dot$results$beta
mallard.results$S.Dot$results$real
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
# Examine output for 'DSR by habitat' model #
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
# Remove "#" on next line to see output
# mallard.results$S.Hab # print MARK output to designated text editor
mallard.results$S.Hab$design.matrix # view the design matrix that was used
mallard.results$S.Hab$results$beta # view estimated beta for model in R
mallard.results$S.Hab$results$beta.vcv # view variance-covariance matrix for beta's
mallard.results$S.Hab$results$real # view the estimates of Daily Survival Rate
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
# Examine output for best model #
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
# Remove "#" on next line to see output
# mallard.results$AgePpnGrass # print MARK output to designated text editor
mallard.results$S.AgePpnGrass$results$beta # view estimated beta's in R
mallard.results$S.AgePpnGrass$results$beta.vcv # view estimated var-cov matrix in R
# To obtain estimates of DSR for various values of 'NestAge' and 'PpnGrass'
# some work additional work is needed.
# Store model results in object with simpler name
AgePpnGrass <- mallard.results$S.AgePpnGrass
# Build design matrix with ages and ppn grass values of interest
# Relevant ages are age 1 to 35 for mallards
# For ppngrass, use a value of 0.5
fc <- find.covariates(AgePpnGrass,mallard)
fc$value[1:35] <- 1:35 # assign 1:35 to 1st 35 nest ages
fc$value[fc$var == "PpnGrass"] <- 0.1 # assign new value to PpnGrass
design <- fill.covariates(AgePpnGrass, fc) # fill design matrix with values
# extract 1st 35 rows of output
AgePpn.survival <- compute.real(AgePpnGrass, design = design)[1:35, ]
# insert covariate columns
AgePpn.survival <- cbind(design[1:35, c(2:3)], AgePpn.survival)
colnames(AgePpn.survival) <- c("Age", "PpnGrass","DSR", "seDSR", "lclDSR",
"uclDSR")
# view estimates of DSR for each age and PpnGrass combo
AgePpn.survival
#library(ggplot2)
#ggplot(AgePpn.survival, aes(x = Age, y = DSR)) +
# geom_line() +
# geom_ribbon(aes(ymin = lclDSR, ymax = uclDSR), alpha = 0.3) +
# xlab("Nest Age (days)") +
# ylab("Estimated DSR") +
# theme_bw()
# assign 17 to 1st 50 nest ages
fc$value[1:89] <- 17
# assign range of values to PpnGrass
fc$value[fc$var == "PpnGrass"] <- seq(0.01, 0.99, length = 89)
# fill design matrix with values
design <- fill.covariates(AgePpnGrass,fc)
AgePpn.survival <- compute.real(AgePpnGrass, design = design)
# insert covariate columns
AgePpn.survival <- cbind(design[ , c(2:3)], AgePpn.survival)
colnames(AgePpn.survival) <-
c("Age", "PpnGrass", "DSR", "seDSR", "lclDSR", "uclDSR")
# view estimates of DSR for each age and PpnGrass combo
AgePpn.survival
# Plot results
#ggplot(AgePpn.survival, aes(x = PpnGrass, y = DSR)) +
# geom_line() +
# geom_ribbon(aes(ymin = lclDSR, ymax = uclDSR), alpha = 0.3) +
# xlab("Proportion Grass on Site") +
# ylab("Estimated DSR") +
# theme_bw()
# If you want to clean up the mark*.inp, .vcv, .res and .out
# and .tmp files created by RMark in the working directory,
# execute 'rm(list = ls(all = TRUE))' - see 2 lines below.
# NOTE: this will delete all objects in the R session.
# rm(list = ls(all=TRUE))
# Then, execute 'cleanup(ask = FALSE)' to delete orphaned output
# files from MARK. Execute '?cleanup' to learn more
# cleanup(ask = FALSE)
# }
Run the code above in your browser using DataLab