Loading and processing the data

Loading the package:

#setwd("c:/users/elie/box sync/seasonalsurvival/analysis")
require(cyclomort)

A nice “piping” workflow: load the data, make the Capture and End date POSIX dates, and create the “cycloSurv” object:

gazelle <- read.csv("../data/gazelle.csv") %>% 
  mutate(CaptureDate = ymd(CaptureDate),
         EndDate = ymd(EndDate),
         T = createCycloSurv(CaptureDate, EndDate, 
                             event = Fate == "Dead",
                             period = 1, timeunits = "years"))

here’s a portion of the relevant data:

sat_ID CaptureDate EndDate Region Fate T
26335a 2002-08-12 2002-10-31 Eastern Dead (0.6132920,0.8323249]
26337 2002-08-14 2003-08-15 Eastern Censored–Alive (0.6187678,1.6208432+]
41599 2003-09-20 2005-04-08 Eastern Censored–Alive (1.7194080,3.2690654+]
41602 2003-09-20 2005-01-08 Eastern Censored–Alive (1.7194080,3.0226535+]
41603 2003-09-21 2003-12-08 Eastern Dead (1.7221459,1.9357029]
41598 2003-09-22 2004-09-15 Eastern Censored–Alive (1.7248838,2.7077937+]

The last column is the cycloSurv object, which encodes the event timing (start, end and censoring) as the number of years from January 1, 2002.

Summary

There are 49 animals, of which:

## 
## Censored--Alive            Dead 
##              29              20

Plot the data:

require(ggplot2); require(ggtheme2)
ggplot(gazelle, aes(CaptureDate, sat_ID, col = Fate)) + 
  geom_errorbarh(aes(xmin = CaptureDate, xmax = EndDate))

There is some very weak evidence of a seasonal pattern (mainly, more deaths in March):

with(subset(gazelle, Fate == "Dead"), barplot(table(month(EndDate))))

Fit some survival models

null model (constant hazard):

f0 <- fit_cyclomort(gazelle$T, n.seasons = 0)
f0
## Multi-seasonal hazard function fit with 0 seasons with periodicity 1.
## 
## $meanhazard
##   estimate     CI.low    CI.high         se 
## 0.41556719 0.26810617 0.64413320 0.09292364 
## 
## Log-likelihood: -37.5622; AIC: 77.1244

The output of the model just gives an estimate of the mean hazard rate (0.41 / year) with confidence intervals

plot(f0, breaks = seq(0,1,length=12))

single season model

f1 <- fit_cyclomort(gazelle$T, n.seasons = 1)
f1
## Multi-seasonal hazard function fit with 1 seasons with periodicity 1.
## 
## $meanhazard
##    estimate    CI.low   CI.high        se
## 1 0.4139975 0.2256053 0.6023897 0.0941961
## 
## $`season 1`
##   parameter season  estimate       CI.low  CI.high
## 1      peak      1 0.1980790 -0.769357793 1.165516
## 2  duration      1 0.4687249  0.003530531 0.500000
## 
## Log-likelihood: -37.5215; AIC: 81.043
plot(f1, breaks = seq(0,1,length=12))

The single season model shows no good evidence for a significant sasonality. It estimates a peak at 0.2 years (March 13), but the confidence interval around that peak is huge. Note also, that the AIC model for the 1 season model is larger than the null model:

sapply(list(f0, f1), AIC)
## [1] 77.12444 81.04300