Plotting daily survival rates | nest survival.

[lbrow43]

Hello.
My top model for American oystercatcher nest survival was predicted by age alone. I am trying to plot the daily survival rate estimates with their confidence intervals. I am having trouble calculating the confidence intervals.

I followed the above example with the formula DSR = exp(log-odds of DSR)/(1+exp(log-odds of DSR)) to calculate the daily survival rate of American oystercatchers across ages 0–35. What would the formula be if I wanted to hand calculate the confidence intervals associated with the daily survival rates? I recognize that I can pull the value of the standard error for the age term from my top model, after naming the top model, and plug this value into a typical standard error formula of lcl = dsr – 1.96*se, but I'm not sure if it makes sense to have a confidence interval that is a consistent amount across all ages so that the graph appears to have a single unchanging width representing the confidence interval.

My main question is if I am pulling the correct value for the standard error to use?

I have previously tried using the function covariate.predictions to get my daily survival rate estimates and confidence intervals in order to plot them, but was not able to successfully do this. However, the above formula for calculating the daily survival rate by hand worked to get my daily survival rate values. I am hoping to use a similar formula to manually calculate the confidence intervals surrounding my daily survival rates, in order to plot them.

Code is below.
Thanks in advance for the help.

Code: Select all

m0 <- mark(NestSurvivalData, nocc = 115, model = "Nest", model.parameters = list(S=list(formula=~NestAge)))
m0

age <- c(rep(1:35)) #get ages
#make predictions
##betas
#intercept
#4.2755890
#age
#-0.0344769
#se of age
#0.0060471

#estimated log odds dsr
logodds = 4.2755890 +(-0.0344769*age)

#DSR = exp(log-odds of DSR)/(1+exp(log-odds of DSR))
dsr = exp(logodds)/(1+exp(logodds))

#calculate ci
#m0 says se of age= 0.0060471. probs need to do calculation, not use raw
lcl = dsr - (1.96*0.0060471)
ucl = dsr + (1.96*0.0060471)

ageplot <- cbind.data.frame(age, dsr, lcl, ucl)

library(ggplot2)
ggplot(ageplot, aes(x = age, y = dsr)) +
  geom_line(size = 2) +
  geom_ribbon(aes(ymin = lcl, ymax = ucl), alpha = 0.2)


[jlaake]

If you want to do manually you will want to read the appendix in Cooch and White that discusses the delta method.

If you want to use covariate.predictions please reply with the problem you encountered.

Jeff

[jlaake]

Here is a simple example with the killdeer data. I had to create a dummy AgeDay1 field because it didn't contain one. This creates a sequence of covariates names NestAge1,NestAge2,...,NestAge40 for the occasions. All of the values are the same in a column because I set all AgeDay1 values to 1. The following shows how to use covariate.predictions from a model. My guess is that you assigned NestAge to the data but that will not work because the covariates are named NestAgex where x is the time for the occasion.

Code: Select all

library(RMark)
data(killdeer)
killdeer$AgeDay1=1
model=mark(killdeer,model="Nest",model.parameters=list(S=list(formula=~NestAge)),nocc=40)
est=covariate.predictions(model,data=data.frame(NestAge1=1:40,index=1))
estimates=est$estimates
plot(estimates$NestAge1,estimates$estimate,xlab="Nest Age in days",ylab="Daily survival rate")


If you only have NestAge in the model then it is that simple. It is more complicated if you have other covariates or time effects.

[lbrow43]

Thank you so much Jeff! This worked.

In a separate model of nest survival to look at the number of American oystercatcher nests that survived to hatching across a nocc = 118, my top model of daily survival rate was predicted by year. I had three years of data covering 2022, 2023, and 2024. When I tried to develop my plot of daily survival rate predicted by year, where I am imagining 3 separate lines of daily survival rate, one line for each year, I ran into issues.

My top model outputs a single estimate of daily survival rate for each year. To get the estimates across the entire season of 118 days for each year, I tried calling up the design matrix and PIMS structure.
This allowed me to get 118 rows of data for my covariates, but lack 118 rows of data for my daily survival rate. I keep returning to the problem that I only have a single estimate for daily survival rate for each year rather than estimates across 118 season days.

Is it possible to form this plot? Or am I only able to obtain a single estimate for each year?

Code: Select all

##Prep for DSR + period survival estimates
m4 =  mark(NestSurvivalData, nocc = 118, model = "Nest", groups = "Year", model.parameters = list(S=list(formula=~Year)))
#first time that I called the input file below in my RMARK code chunk within R, I saw each year have nest survival predicted out for 117 days. How do I access this data?
m4
#But when I called the object m4 as a separate external window, I am only given a single estimate for each year. What is happening?
m4$results


m4df = data.frame(m4$results$real)
  #usually works in past to get real vals across nocc

#SO close (has 351 datapoints 117 x 3 yrs). just need survival estimates added
m4$design.data$S
m4$design.data$pimtypes$S # i am using all pims, so what's problem?

#closer
PIMS(m4, "S", simplified = F)
PIMS(m4, "S", simplified = T)

PIMS(amoy.results1[[3]], "S")
PIMS(amoy.results1[[1]], "S")

as.numeric(row.names(m4$Phi[m4$Phi$group=="Year",]))

#make sequence of yr vals
yrvalues=c(2022, 2023, 2024)
yrvalues = data.frame(Year = yrvalues)

#predict DSR over yr
yr.predict <- covariate.predictions(m4, data=yrvalues, indices = c(1))
#this should make the prediction matrix

#next test that the values in the prediction matrix make sense
yr.predict #print it console so you can see
#look -- estimates appear in R below this chunk!
yr.predict$estimates



[jlaake]

The model ~Year means that you have a constant DSR within each year but they differ across years. So there are only 3 different estimates. Did you mean to fit ~Year*NestAge?

[lbrow43]

OK thanks for confirming. I only have year as a predictor. That is what I thought but I wanted to double check that I should not be receiving estimates across the 118 days or plotting the data.
Thanks Jeff