mallard | R Documentation |
A nest survival data set on mallards. The data set and analysis is described by Rotella et al.(2004).
A data frame with 565 observations on the following 13 variables.
the day the nest was first found
the last day that chicks were present
the last day the nest was checked
the fate of the nest; 0=hatch and 1 depredated
the frequency of nests with this data; always 1 in this example
Robel reading of vegetation thickness
proportion grass in vicinity of nest
dummy 0/1 variable; 1 if native vegetation
dummy 0/1 variable; 1 if planted vegetation
dummy 0/1 variable; 1 if wetland vegetation
dummy 0/1 variable; 1 if roadside vegetation
age of nest in days the day the nest was found
age of nest at beginning of study
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.
Jay Rotella
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.
# 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) # 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.