www.phidot.org

Hi, I'm trying to run some CJS models in marked to look at apparent survival in a local American kestrel population. There's a network of ~150 nest boxes where we band adult males/females during incubation, and nestlings that we're able to sex at ~25 days old. I would like to look at the effect of winter severity on survival and how survival differs between the two age classes (birds marked as adults and birds marked as juveniles). My question is whether I've properly created an age-specific data structure for the analysis. Here's my code and data below; age is indicated as 'A' for adults and 'Y' for birds banded as nestlings.

Code: Select all

library(marked)
library(tidyverse)

# read in the data --------------------------------------------------------

# kestrel capture histories
kestrel <- read_csv('https://raw.githubusercontent.com/jaymwin/kestrelSurvival/master/kestrel_ch.csv',
                    col_types = cols(
                      ch = col_character(),
                      sex = col_factor(),
                      age = col_factor()
                      )
                    ) %>%
  as.data.frame() # leading 0's disappear when I read this data with read.csv
kestrel

# winter severity data
anoms <- read_csv('https://raw.githubusercontent.com/jaymwin/kestrelSurvival/master/ID_daymet_anoms.csv') %>%
  filter(winter_year > 2007 & winter_year < 2018)

# manipulate and bind with capture histories
anoms <- anoms %>%
  mutate(
    winter_year = str_c('tmin', winter_year)
  ) %>%
  spread(winter_year, tmin_anom) %>%
  as.data.frame()
anoms

kestrel <- cbind(kestrel, anoms)
kestrel


# process data ------------------------------------------------------------

kestrel.proc <- process.data(kestrel,
                             model = 'cjs',
                             begin.time = 2008)

# create design data with static and time-varying covariates
design.Phi = list(static = c('sex'),
                  time.varying = c('tmin') # winter severity
)

design.p=list(static = c('sex'))

design.parameters <- list(Phi = design.Phi,
                          p = design.p)

kestrel.ddl <- make.design.data(kestrel.proc,
                                parameters = design.parameters)
kestrel.ddl

# now deal with adults and hatch-year birds; is this right?
kestrel.ddl$Phi$adult=0
kestrel.ddl$Phi$adult[kestrel.ddl$Phi$Age>=1]=1
kestrel.ddl$Phi$young=1


# run models --------------------------------------------------------------

# fit the models; not an exhaustive list
fit.models = function()
{
  Phi.sex = list(formula = ~ sex)
  Phi.sex.time = list(formula = ~ sex + time)
  Phi.time = list(formula = ~ time)
  Phi.tmin = list(formula = ~ tmin)
  Phi.sex.tmin = list(formula = ~ sex + tmin)
  Phi.sex.tmin.time = list(formula = ~ sex + tmin + time)
  Phi.age = list(formula = ~ adult) # age model
  p.sex = list(formula = ~ sex)
  p.dot = list(formula = ~ 1)
  p.time = list(formula = ~ time)
  cml = create.model.list(c("Phi","p"))
  results = crm.wrapper(cml, data = kestrel.proc, ddl = kestrel.ddl, hessian = TRUE,
                        external = FALSE, accumulate = FALSE)
  return(results)
}
kestrel.models = fit.models()
kestrel.models

Does this look correct?

Thanks,

Jay

I think I've dealt correctly now with HY birds transitioning to AHY by re-formatting the data and using RMark instead (see code below), but I'm still uncertain how to do the same thing in marked.

Code: Select all

library(tidyverse)
library(RMark)

CH_final <- read_csv('https://raw.githubusercontent.com/jaymwin/kestrelSurvival/master/kestrel_ch.csv',
                     col_types = cols(
                       ID = col_factor(),
                       ch = col_character(),
                       sex = col_factor(),
                       year = col_factor(),
                       age = col_factor()
                     )
) %>%
  as.data.frame()
CH_final

# Set up the basic model structure -------
kestrel_process <- process.data(data = CH_Final,
                                model = "CJS",
                                begin.time = 2008,
                                initial.ages = c(0,1),
                                groups = c("sex","age"))

kestrel.ddl = make.design.data(kestrel_process,
                               parameters = list(Phi = list(age.bins = c(0, 1, 11)),
                                                 p = list(age.bins = c(1, 2, 11))), right = FALSE)

Thanks,

Jay

Sorry but I have been away from my computer for a few days. One of the advantages of being retired.

What you did with marked was fine if all you did was band nestlings. But that was not the case based on what you said and what you did with RMark. I think you'll find that the setup approach in marked is fairly identical to what you did with RMark. When I wrote marked, I wanted them to be as similar as possible. Here is an example using the dipper data. Note I had to use the age.var argument because the groups variable I chose isn't named age.

Code: Select all

library(marked)
data(dipper)
# create phony age data
dipper$ahy="Adult"
dipper$ahy[1:100]="HY"
dipper$ahy=factor(dipper$ahy)
#process data assigning initial ages
dp=process.data(dipper,groups="ahy",initial.age=c(1,0),age.var=1,accumulate=FALSE)
# make design data
ddl=make.design.data(dp)
# create 0, 1+ age variable
ddl$Phi$ya=cut(ddl$Phi$Age,c(0,1,10),right=FALSE,labels=c("HY","AHY"))
# fit model
model=crm(dp,ddl,model.parameters=list(Phi=list(formula=~ya)))
# append design data to design matrix to show model
xx=cbind(ddl$Phi,as.matrix(model$results$model_data$Phi.dm))
head(xx)
tail(xx)


Meant to mention that initial.ages values are in order of the levels of the factor variable used for age (ahy in this example). They are alphabetical by default so AHY is before HY so initial.ages are 1,0 rather than 0,1. Always best to look at your design data and see if the created values are what you want.