- 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