Plotting covariate predictions by age class

[Andre23]

Hi, I ran a set of models in RMark using an individual covariate called "SFdist" (numerical value, range 0 - 5.5) as a potential impact on survival.

My hypothesis was that it could impact survival at all ages or only for juveniles.

Code: Select all

# 2. Process the data
cjs_proc <- process.data(rmark_data, model = "CJS", begin.time = start_year,groups="subpop")

# 3. Create design data
cjs_ddl <-  make.design.data(cjs_proc)

#Add age column to the design data (age classes)
cjs_ddl$Phi$ageclass3=cut(cjs_ddl$Phi$Age,c(0, 1, 4, Inf), right =FALSE)
levels(cjs_ddl$Phi$ageclass3)=c("Juveniles","Adults(1-3)","Adults(4+)")
# Create an indicator for juveniles (age0)
cjs_ddl$Phi$is_juvenile <- ifelse(cjs_ddl$Phi$ageclass3 == "Juveniles", 1, 0)
# Create an indicator for young adults (age1-3)
cjs_ddl$Phi$youngA <- ifelse(cjs_ddl$Phi$ageclass3 == "Adults(1-3)", 1, 0)
# Create an indicator for older adults (age4+)
cjs_ddl$Phi$olderA <- ifelse(cjs_ddl$Phi$ageclass3 == "Adults(4+)", 1, 0)

#modify another parameter (e.g. "p")
cjs_ddl$p$ageclass6=cut(cjs_ddl$p$Age,c(-Inf,2,3,4,5,6,Inf),right=FALSE)
levels(cjs_ddl$p$ageclass6) <- c("1","2","3","4","5","6+")

#Run set of models
cjs_time_age <- mark(cjs_proc, cjs_ddl, model.parameters = list(
  Phi = list(formula=~time:ageclass3),
  p = list(formula=~time:ageclass6+subpop,remove.intercept=TRUE,link="logit")))

cjs_time_age_SF <- mark(cjs_proc, cjs_ddl, model.parameters = list(
  Phi = list(formula=~time:ageclass3+SFdist),
  p = list(formula=~time:ageclass6+subpop,remove.intercept=TRUE,link="logit")))

cjs_time_age_SF_juv <- mark(cjs_proc, cjs_ddl, model.parameters = list(
  Phi = list(formula=~is_juvenile:time:SFdist + youngA:time + olderA:time), #seperate age classes, try with SFdist effect on juveniles only
  p = list(formula=~time:ageclass6+subpop,remove.intercept=TRUE,link="logit")))

How can I plot how SFdist is impacting:
1. Juveniles survival only?
2. Average survival estimates at any age class of the entire dataset?

I aim to achieve it the same way as in the RMark workshop (12. Covariate Predictions) with the weight effect on survival.

Thanks!

[jlaake]

You don't say whether you want to plot for an individual model or use model averaging. I assume it is the latter. Either use collect.models or change code to use mark.wrapper to get a set of model results in a marklist. Then use model average and specify the values of the SFdist that you want in a data frame passed as argument data to model.average. You can see specify computation for specific age parameters by specifying the indices from the design data for the survival parameter. See ?collect.models, ?model average.marklist.

[Andre23]

Hi @jlaake

Thank you, I managed to do it with covariate.predictions function

## Collect model results, adjust for c-hat value, and create table of models ranked by QAICc
model.list=collect.models(type="CJS")
adjusted.models<-adjust.chat(chat = 1.27, model.list)
my.c.models<-model.table(adjusted.models)

#Plot covariate effect on juveniles survival (plot for values of SFdist from 0 to 5)
Phi.SFdist.predictions=covariate.predictions(cjs_time_age_SF_juv,data=data.frame(index=rep(1,6),SFdist=0:5))

#plot
with(Phi.SFdist.predictions$estimates,
{
plot(SFdist, estimate,type="l",lwd=2,xlab="SFdist",
ylab="Survival",ylim=c(0,1))
lines(SFdist,lcl,lty=2)
lines(SFdist,ucl,lty=2)
})

[jlaake]

Yes it will work for a single model. If you choose to average the effect over models, use model average.