I've been thinking about ways to tackle an epidemic modelling problem I've been working on, and I've come up against a conceptual difficulty over the way survival analysis works. Here's a really simplified version with all the extraneous details stripped out.
There is a collection of $n$ individuals. At the beginning (time $t = 0$), 1 individual is infectious with a disease, and the other $n-1$ individuals are susceptible. As time progresses, susceptible people can get the disease through contact with infectious individuals, and pass the disease on to other people. To keep focus on the core of my question we'll use these (unrealistic) simplifying assumptions:
- There are no births, deaths, immigration, or emigration
- Individuals mix uniformly (e.g. no preferential mixing by age or sex or location, etc)
- No delay between contracting the disease and ability to infect others (i.e. no incubation period)
- Once infected, an individual remains infectious forever (no recovery)
Here is how transmission works: we use a continuous-time additive hazard model. This means that we let each uninfected individual $i$ have their own hazard function $h_i(t)$. This function tells us "given that the individual $i$ has survived until the time $t$, what is the infinitesimal rate of failure (i.e. infection) at time $t$?"
We define it as follows:
$$h_i(t ;\lambda) = \lambda \sum_{j \neq i}^n I_j(t)$$
Where $I_j(t)$ is simply an indicator function that is 1 when individual $j$ is infectious at time $t$, and 0 otherwise. The $\lambda$ is a strictly-positive parameter governing how transmissible the disease is: the higher the value, the faster the disease spreads. Given our assumptions, $h_i(t)$ is thus piecewise constant and non-decreasing: it looks like a staircase as the epidemic spreads to more people and the risk of infection increases.
The goal is this: given a dataset consisting of the times $t_1, t_2, \dots, t_n$ that each individual became infected, infer the value of $\lambda$. At first this seems simple: since we have complete data, we can easily calculate what the hazard functions must have been for each individual over time, and we can work out the pdf for the infection times from the hazard functions:
$$f_i(t | \lambda) = h_i(t)\exp{\left(- \int_0^t h_i(u)du \right)}$$
(The integral in the exponential is just the cumulative hazard function, see the Wikipedia link for more details.)
With this, it's a "simple" matter of writing down the likelihood of the entire data and maximizing it with respect to $\lambda$:
$$\hat{\lambda} = \mathrm{argmax}_\lambda \prod_i^n f_i(t_i | \lambda)$$
But here's the problem:
Writing the product of probabilities in that way is assuming that the events are independent, and I'm really not so sure that's the case. If we were to change one of the infection times (e.g. $t_5$ becomes $t'_5$), that would propagate through and potentially change many of the hazard functions, because they are dependent on when the other individuals become infectious (via $I_j(t)$). This would go on to change the other probability densities. It would be as if, in a linear regression model, some units' $X$-covariates causally depended on other units' $Y$-observations.
On the other hand, I can't think of how else to write down the likelihood. And well, it's damned tempting to just ignore the issue entirely because it gives me a headache to think about otherwise.
My questions are:
- Is this a big problem, or can I just ignore it? Will I get anything meaningful by just maximizing that "likelihood"?
- If it is a big problem, is this approach salvageable or do I just need to scrap the whole idea?
- If the latter, what are some other ways of going about this (this being analysing epidemics in continuous time rather than discrete)?
If anyone is wondering about the XY-problem, let me briefly explain why I'm formulating it this way. I want to make a stochastic continuous-time model for my epidemic data, and since it can be understood as "time-to-failure", survival analysis seems like a logical way of tackling it. Having it be an additive (rather than proportional) risk model works for this case because there is no such thing as a "baseline hazard" in this kind of epidemic (i.e. if there are no infected people around, the chance of transmission is exactly zero, which cannot be represented in a proportional hazard model). It also lets me incorporate more complexity in a very natural way: I can have changing population size and incubation periods just by making straightforward changes to the infection indicator functions, and I can have preferential mixing by introducing more $\lambda$-like parameters that act on covariate information of the individuals.
