Rotella wrote:Here is an updated version of the script that includes code for working with the output.
- Code: Select all
# 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. #
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
library(RMark)
# 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)
}
# 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
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)
[unparseable or potentially dangerous latex formula]
Hi, Rotella
I try to obtain estimates of DSR for various values of 'NestAge' and 'Year' in my analysis following your code showed above. While the result I got was strange, please see the code and results as below.
Thanks.
Chenyang Liu
AgeYear <- RWB.results$AgeYear
fc <- find.covariates(AgeYear,RWB)
fc$value[1:24] <- 1:24 # assign 1:24 to 1st 24 nest ages
fc$value[fc$var == "Year"] <- 2017 # assign new value to AGE
design <- fill.covariates(AgeYear,fc) # fill design matrix with values
RWB.Age.survival <- compute.real(AgeYear,design=design)[1:24,]
RWB.Age.survival <- cbind(design[1:24,c(1:2)],RWB.Age.survival)
colnames(RWB.Age.survival) <- c("Age","Year","DSR","seDSR","lclDSR","uclDSR")
RWB.Age.survival
Age Year DSR seDSR lclDSR uclDSR NA
1 1 0 0.9804642 0.004109337 0.9705531 0.9870839
2 1 0 0.9801073 0.004139078 0.9701485 0.9867890
3 1 0 0.9797440 0.004171609 0.9697295 0.9864915
4 1 0 0.9793742 0.004207171 0.9692954 0.9861917
5 1 0 0.9789978 0.004246011 0.9688455 0.9858899
6 1 0 0.9786146 0.004288385 0.9683791 0.9855863
7 1 0 0.9782247 0.004334554 0.9678954 0.9852812
8 1 0 0.9778278 0.004384784 0.9673935 0.9849749
9 1 0 0.9774238 0.004439346 0.9668727 0.9846676
10 1 0 0.9770126 0.004498512 0.9663322 0.9843597
11 1 0 0.9765941 0.004562555 0.9657712 0.9840514
12 1 0 0.9761682 0.004631750 0.9651888 0.9837430
13 1 0 0.9757347 0.004706368 0.9645842 0.9834348
14 1 0 0.9752936 0.004786679 0.9639565 0.9831271
15 1 0 0.9748446 0.004872950 0.9633050 0.9828200
16 1 0 0.9743877 0.004965440 0.9626287 0.9825139
17 1 0 0.9739227 0.005064406 0.9619268 0.9822089
18 1 0 0.9734495 0.005170096 0.9611985 0.9819052
19 1 0 0.9729679 0.005282753 0.9604430 0.9816031
20 1 0 0.9724779 0.005402609 0.9596594 0.9813026
21 1 0 0.9719792 0.005529893 0.9588469 0.9810039
22 1 0 0.9714718 0.005664822 0.9580047 0.9807072
23 1 0 0.9709554 0.005807607 0.9571319 0.9804126
24 1 0 0.9704300 0.005958451 0.9562279 0.9801200
> fc$value[1:9] <- 12 # assign 1:24 to 1st 24 nest ages
> fc$value[fc$var == "Year"] <- seq(2012,2020, length=9) # assign new value to AGE
> design <- fill.covariates(AgeYear,fc) # fill design matrix with values
> RWB.Age.survival <- compute.real(AgeYear,design=design)[1:9,]
> RWB.Age.survival <- cbind(design[1:9,c(2:3)],RWB.Age.survival)
> colnames(RWB.Age.survival) <- c("Age","Year","DSR","seDSR","lclDSR","uclDSR")
> RWB.Age.survival
Age Year DSR seDSR lclDSR uclDSR NA
1 0 0 0.9761682 0.00463175 0.9651888 0.983743
2 0 0 0.9761682 0.00463175 0.9651888 0.983743
3 0 0 0.9761682 0.00463175 0.9651888 0.983743
4 0 0 0.9761682 0.00463175 0.9651888 0.983743
5 0 0 0.9761682 0.00463175 0.9651888 0.983743
6 0 0 0.9761682 0.00463175 0.9651888 0.983743
7 0 0 0.9761682 0.00463175 0.9651888 0.983743
8 0 0 0.9761682 0.00463175 0.9651888 0.983743
9 0 0 0.9761682 0.00463175 0.9651888 0.983743