Learn R Programming

RMark (version 3.0.0)

mallard: Mallard nest survival example

Description

A nest survival data set on mallards. The data set and analysis is described by Rotella et al.(2004).

Arguments

Format

A data frame with 565 observations on the following 13 variables.

FirstFound

the day the nest was first found

LastPresent

the last day that chicks were present

LastChecked

the last day the nest was checked

Fate

the fate of the nest; 0=hatch and 1 depredated

Freq

the frequency of nests with this data; always 1 in this example

Robel

Robel reading of vegetation thickness

PpnGrass

proportion grass in vicinity of nest

Native

dummy 0/1 variable; 1 if native vegetation

Planted

dummy 0/1 variable; 1 if planted vegetation

Wetland

dummy 0/1 variable; 1 if wetland vegetation

Roadside

dummy 0/1 variable; 1 if roadside vegetation

AgeFound

age of nest in days the day the nest was found

AgeDay1

age of nest at beginning of study

Author

Jay Rotella

Details

The first 5 fields must be named as they are shown for nest survival models. Freq and the remaining fields are optional. See killdeer for more description of the nest survival data structure and the use of the special field AgeDay1. The habitat variables Native,Planted,Wetland,Roadside were originally coded as 0/1 dummy variables to enable easier modelling with MARK. A better alternative in RMark is to create a single variable habitat with values of "Native","Planted", "Wetland", and "Roadside". Then the Hab model in the example below would become:

 Hab=mark(mallard,nocc=90,model="Nest",
model.parameters=list(S=list(formula=~habitat)), groups="habitat") 

For this example, that doesn't make a big difference but if you have more than one factor and you want to construct an interaction or you have a factor with many levels, then it is more efficient to use factor variables rather than dummy variables.

Examples

Run this code

# 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