# NOT RUN {
## load Western Arctic Herd data and convert to cycloSurv
data(wah_morts)
wah <- with(wah_morts, create_cycloSurv(start = start, end = end,
event = fate == "dead", period = 365))
# censor and trim
cutoff = "2016-01-01"
wah_pre = censor_cycloSurv(wah, censor.time = cutoff)
wah_post = trim_cycloSurv(wah, trim.time = cutoff)
# combine into dataframe
par.init <- par(no.readonly = TRUE)
par(mfrow = c(1,2))
plot(wah_pre[,1], 1:length(wah_pre), xlim = range(wah_pre[,1:2]), type= "n", main = "pre")
segments(wah_pre[,1], 1:length(wah_pre), wah_pre[,2], 1:length(wah_pre), col = wah_pre[,3]+1)
plot(wah_post[,1], 1:length(wah_post), xlim = range(wah_post[,1:2]), type= "n", main = "post")
segments(wah_post[,1], 1:length(wah_post), wah_post[,2], 1:length(wah_post), col = wah_pre[,3]+1)
# fit seasonal model before and after
wah_fit_pre <- fit_cyclomort(wah_pre, n.seasons = 1)
wah_fit_post <- fit_cyclomort(wah_post, n.seasons = 1)
# some evidence of a shift, though confidence intervals are wide
summary(wah_fit_pre)
summary(wah_fit_post)
# }
# NOT RUN {
par(mfrow = c(1,2))
plot(wah_fit_pre, plotCI = TRUE, breaks = 10); title("pre cut-off")
plot(wah_fit_post, plotCI = TRUE, breaks = 10); title("post cut-off")
# }
# NOT RUN {
par(par.init)
# }
Run the code above in your browser using DataLab