Survival analysis models time-to-event as an outcome (e.g. time to death, how long do you wait for an animal to appear in a camera trap, when will an earthquake strike, when will a mango fall, etc). There are, in short, lots of ecological applications. But the biggest by far is, in fact, modeling survival.
The goals in this lecture are to get an intuition for why time-to-event data is fundamentally different than other kinds of observations, and - on the basis of a simple experiment - get intuition for the hazard function and censored data and how to estimate those using maximum likelihood.
Discrete Hazard process
We rolled dice - 52 to be precise - each until we got a 6.
Each roll represented a year. Once a 6 was rolled, the die was dead and removed. Like with these
lethal dice
Here are the data:
X <- read.csv("survival.csv")[,1]
mean(X)
## [1] 7.211538
plot(table(X), type = "h")
If you like, you can also simulate this process:
X <- rep(NA, 50)
for(i in 1:length(X)){
rolls <- 0
x <- 0
while(x != 6){
rolls <- rolls+1
x <- sample(1:6,1)
}
X[i] <- rolls
}
Probability distribution
The probability mass function of mortality time is:
\[Pr(X = x) = f(x|p) = p(1-p)^{x-1}\]
This can be reasoned out, but also has a name, the negative binomial distribution: \(X \sim \text{NegBinom}(p)\)1Important note: the standard definition of the negative binomial is number of failures before first success, so it includes 0, whereas in our definition if the event dies in the first attempt, that’s a 1, not a 0, so technically: \(X+1 \sim \text{NegBinom}(p)\). This becomes important later..
But now we can also get a confidence interval, using the all important Hessian. Recall that for a maximum likelihood estimate \(\widehat{\theta}\) , an asymptotic estimate of its variance is given by \[\sigma^2_{\hat\theta} = {-1\over d^2l(\hat\theta)/d\theta^2}\]
i.e., the greater the (negative) curvature, the smaller the standard error. Here’s one way to compute this:
require(Rdistance)
p.dd <- secondDeriv(p.hat$maximum, p.ll, X = X)
se = sqrt(-1/p.dd)
(CI = p.hat$maximum + c(-2,2) * se[1,1])
## [1] 0.1029695 0.1743538
optim, rather than optimize, does this all at once if you ask for the “Hessian”. Some notes about the somewhat fussy (but very important) optim function:
- by default it minimizes, so to make it maximize you have to add the control = list(fnscale = -1) OR you can just make your function return the negative log-likelihood,
- method = "L-BFGS-B" gives “box constraints” - limiting the search to a known range, but because the log-likelihood is infinity at 0 and 1, we have to
- the par argument is an “initial guess” - the place where the optimizer will start looking for a maximum.
This is a rare likelihood that you can actually solve analytically - which is a satisfying problem to solve. To get the estimate you have to take the derivative of the log-likelihood against \(p\) and set it to 0.
The standard error of this estimate can also be calculated analytically! That’s an exercise for a lazy Sunday afternoon, or a train ride, or some other time.
The solution is \(\sigma = {n(C-n) \over C^2(C-n) + C^2n}\) where \(C = \sum X_i\). Check it out:
n <- length(X)
C <- sum(X)
list(sigma_math = sqrt(n*(C-n)/(C^2 * (C-n) + C^2*n)),
sigma_optim = p.se)
Censoring is a very important feature of survival data. Basically, you can’t follow your population forever, so at some point objects “fall out” of the study. You don’t know when the event happened, just that is had to happen after we were done making observations.
Here, we’ll do something very simple: we’ll simply say we can only follow our individuals for 6 years. Afterwards, we don’t track them anymore.
Let’s generate more data (this time I just use the negative binomial shortcut, and manually censor):
X <- read.csv("survival.csv")[,1]
X.full <- X
Y <- rep(6, sum(X > 6))
X <- X[X <= 6]
The likelihood function now asks what is \(p\) given \(\bf X\), which contains \(n\) known fate observations, and \(\bf Y\), which contains \(m\) censored observations is:
Ok, that’s kind of a cop out. It’s not so easy or necessary to write out, but we can code it up using negative binomial short cuts. Namely that \(f(x,p)\) is (very close to) the negative binomial probability mass function, and \(g(y,p)\) is one minus the negative binomial cumulative mass function.2 THIS is where that differnece is important!
p.censored.ll <- function(p, Z, Y){
# true observations are the probability function dnbinom
log.Z <- dnbinom(Z-1, prob = p, size = 1, log = TRUE)
# censored observations are the 1 - the cumulative probability pnbinom
log.Y <- log(1-pnbinom(Y-1, prob = p, size = 1))
# here, the log has to bo on the outside, can you see why?
sum(log.Z) + sum(log.Y)
}
Note that in R, all distribution functions have a log option … now you maybe understand why:
Let’s compare the likelihood function for both the full \(X\) dataset, and the censored version:
