# Hazard function for a growing population

I want to model the expected failure rate in a growing population.

Each individual has a high infant mortality (e.g. Weibull(shape=0.5, scale=100)).

So the failure rate (or hazard function) at time t, H(t), is the density function divided by the survival function: H(t) = pdf(t)/sf(t)

The complication is that at every timestep, I have many individuals in my population. They have different age, and therefore different expected failure rates, which is given by their hazard function h(t). The expected failure rate of the population as a whole is a mix of these hazard functions. How do I compute that?

My first instinct was to sum up densities and survival functions individually, and then divide one by the other.

pdfs = numpy.zeros(max_time)
sfs = numpy.zeros(max_time)
for t in individual_start_times:
ts = numpy.arange(max_time-t)
pdfs[t:] += scipy.stats.weibull_min.pdf(ts, shape=0.5, scale=100)
sfs[t:] += scipy.stats.weibull_min.sf(ts, shape=0.5, scale=100)
hazard = pdfs/sfs


Does that make sense, or is that bogus? I'm struggling to verify this.

• As an epidemiologist reading "high infant mortality" is quite confusing: "high" implies that "infant mortality" which "high" intensifies is an infant mortality rate. However, individuals do not have mortality rates: they are either dead or alive (or sometimes "kinda dead" or "mostly dead", c.f. S. Morgenstern) Oct 2, 2018 at 15:43
• @Alexis Well, an individual can be not even born yet, and still we can have expectations about their lifespan. In other words, consider an individual as a random variable, not as a particular outcome. Suppose you know the dates at which kids will be born in an area, and you know the distribution of their lifespan. Then you have a probability of death at day t for a particular kid. But how do you get the expected death rate at day t of the entire population? Oct 2, 2018 at 17:10
• I would use "expected mortality risk" (or, if you care about life span "life expectancy"). Oct 2, 2018 at 17:30
• Where you say "expected failure rate of the population as a whole" are you trying to calculate the distribution of the time to extinction (i.e. the expected time until the entire population is dead)? Oct 4, 2018 at 10:19
• No, I care about the expected failure rate. The expected number of deaths over expected number alive. Oct 4, 2018 at 10:28

If I understood you correctly, you could fit a Weibull model to your data with Age as the covariate/predictor. And then from this model extract the estimated hazard function. The following code illustrates how this can be done in R using the survival package:

library("survival")

# Accelerated Failure Time (AFT) Weibull model
fm <- survreg(Surv(time, status) ~ age, data = lung, dist = "weibull")
summary(fm)

# transform coefficients to Proportional Hazards (PH) scale
# (Weibull only model capable of having both AFT and PH formulations)
shape <- 1 / fm\$scale
X <- model.matrix(fm)
scale <-  exp(- shape * c(X %*% coef(fm)))

# function to calculate hazard function for
# every subject in the data for a set of time points
ht <- function (time, shape, scale) {
outer(shape * time^(shape - 1), scale)
}

# A figure of the estimated hazard for the
# first 5 subjects
times <- seq(0.001, 1000, length = 500)
matplot(times, ht(times, shape, scale[1:5]), type = "l", lty = 1,
xlab = "Follow-up Time (days)", ylab = "Estimated Hazard Function")
legend("bottomright", paste("Patients", 1:5), col = 1:5, lty = 1, bty = "n")


• Thanks for your answer. This isn't what I'm looking for though. I have the fitted Weibull distribution already, and I know the hazard function of this distribution. But for any time t, I have multiple things alive, and each thing has a different age and thus a different chance of failure at that time. So I have the expected chance of failure of each instance, and what I'm looking for is the expected failure rate of the entire population at time t. Oct 4, 2018 at 19:08
• Thanks for the clarification. Then would it work if you fitted the Weibull model without age as a predictor? This would give you one hazard function for the whole sample that would be an estimate of the hazard function of the entire population. Oct 4, 2018 at 19:13
• I don't think so. I don't have a sample I'm trying to fit to. I'm doing predictions about the future, so I have a rate at which new lifetimes start, each life follows the same Weibull distribution, and I want to draw conclusions about the global failure rate at any future time. Oct 4, 2018 at 19:19