# Better Than New#

Click here to run this notebook on Colab.

## Show code cell content

```
# Install empiricaldist if we don't already have it
try:
import empiricaldist
except ImportError:
!pip install empiricaldist
```

## Show code cell content

```
# download utils.py
from os.path import basename, exists
def download(url):
filename = basename(url)
if not exists(filename):
from urllib.request import urlretrieve
local, _ = urlretrieve(url, filename)
print("Downloaded " + local)
download(
"https://github.com/AllenDowney/ProbablyOverthinkingIt/raw/book/notebooks/utils.py"
)
```

## Show code cell content

```
DATA_PATH = "https://github.com/AllenDowney/ProbablyOverthinkingIt/raw/book/data/"
```

## Show code cell content

```
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(font_scale=1.1)
sns.set(rc={'figure.figsize':(6, 3.5)})
sns.set_style('white')
# Set the random seed so we get the same results every time
np.random.seed(17)
```

```
import utils
def decorate(**options):
utils.decorate(**options)
sns.despine(left=True, bottom=True)
```

```
from utils import fit_normal
from utils import make_cdf
from utils import normal_error_bounds
def make_plot(series, model_label=None, plot_bounds=True, qs=None, **options):
"""Plot the data and a normal model.
series: Series
model_label: string
plot_bounds: boolean, whether to plot the expected upper bound
qs: sequence of quantities where the model should be evaluated
options: passed to plt.plot
"""
cdf = make_cdf(series)
dist = fit_normal(series)
p = 1 - (1 / len(series))
upper = dist.ppf(p)
if plot_bounds:
plt.plot(upper, p * 100, "+", color="black", ms=12)
n = len(series)
if qs is None:
q_max = max(cdf.qs.max(), upper)
qs = np.linspace(cdf.qs.min(), q_max)
low, high = normal_error_bounds(dist, n, qs)
plt.fill_between(qs, low * 100, high * 100, lw=0, color="gray", alpha=0.2)
cdf.plot(**options)
```

```
def savefig(filename, **options):
if 'dpi' not in options:
options['dpi'] = 300
plt.savefig(filename, **options)
```

Suppose you work in a hospital, and one day you have lunch with three of your colleagues. One is a facilities engineer working on a new lighting system, one is an obstetrician who works in the maternity ward, and one is an oncologist who works with cancer patients. While you all enjoy the hospital food, each of them poses a statistical puzzle.

The engineer says they are replacing old incandescent light bulbs with LED bulbs, and they’ve decided to replace the oldest bulbs first. According to previous tests, the bulbs last 1400 hours on average. So, they ask, which do you think will last longer: a new bulb or one that has already been lit for 1000 hours?

Sensing a trick question, you ask if the new bulb might be defective. The engineer says, no, let’s assume we’ve confirmed that it works. In that case, you say, I think the new bulb will last longer.

“That’s right,” says the engineer. “Light bulbs behave as you expect; they wear out over time, so the longer they’ve been in use, the sooner they burn out, on average.”

“However,” says the obstetrician, “not everything works that way. For example, most often, pregnancy lasts 39 or 40 weeks. Today I saw three patients who are all pregnant; the first is at the beginning of week 39, the second is at the beginning of week 40, and the third is at the beginning of week 41. Which one do you think will deliver her baby first?”

Now you are sure it’s a trick question, but just to play along, you say the third patient is likely to deliver first.

The obstetrician says no, the remaining duration of the three pregnancies is nearly the same, about four days. Even taking medical intervention into account, all three have the same chance of delivering first.

“That’s surprising,” says the oncologist. “But in my field things are even stranger. For example, today I saw two patients with glioblastoma, which is a kind of brain cancer. They are about the same age, and the stage of their cancers is about the same, but one of them was diagnosed a week ago and one was diagnosed a year ago. Unfortunately, the average survival time after diagnosis is only about a year. So you probably expect the first patient to live longer.”

By now you know better than to guess, so you wait for the answer.

The oncologist explains that many patients with glioblastoma live only a few months after diagnosis.
So, it turns out, a patient who survives one year after diagnosis is then *more* likely to survive a second year.

Based on this conversation, we can see that there are three ways survival times can go:

Many things wear out over time, like light bulbs, so we expect something new to last longer than something old.

But there are some situations, like patients after a cancer diagnosis, that are the other way around: the longer someone has survived, the longer we expect them to survive.

And there are some situations, like women expecting babies, where the average remaining time doesn’t change, at least for a while.

In this chapter I’ll demonstrate and explain each of these effects, starting with light bulbs.

## Light Bulbs#

The data are available in a gist:

```
download(
"https://gist.github.com/epogrebnyak/7933e16c0ad215742c4c104be4fbdeb1/raw/c932bc5b6aa6317770c4cbf43eb591511fec08f9/lamps.csv"
)
```

```
df = pd.read_csv("lamps.csv", index_col=0)
df.tail()
```

h | f | K | |
---|---|---|---|

i | |||

28 | 1812 | 1 | 4 |

29 | 1836 | 1 | 3 |

30 | 1860 | 1 | 2 |

31 | 1980 | 1 | 1 |

32 | 2568 | 1 | 0 |

The `h`

column contains the quantities; the `k`

column contains the frequencies (counts).

```
from empiricaldist import Pmf
pmf_lifetimes = Pmf(df["f"].values, index=df["h"])
pmf_lifetimes.index.name = "t"
```

We can use a `Counter`

to reconstitute the 50 individual lifetimes.

```
from collections import Counter
lifetimes = pd.Series(Counter(pmf_lifetimes.to_dict()).elements())
lifetimes.describe()
```

```
count 50.000000
mean 1413.840000
std 346.520628
min 840.000000
25% 1176.000000
50% 1446.000000
75% 1653.000000
max 2568.000000
dtype: float64
```

Here’s a normal distribution that has the same mean and standard deviation as the data.

```
from scipy.stats import norm
mu, sigma = lifetimes.mean(), lifetimes.std()
dist = norm(mu, sigma)
mu, sigma
```

```
(1413.84, 346.52062753337276)
```

The following figure shows the distribution of survival times for these light bulbs, plotted as a cumulative distribution function (CDF), along with a Gaussian model.

```
from empiricaldist import Cdf
qs = np.linspace(500, 2500)
ps = dist.cdf(qs)
plt.plot(qs, ps, color="gray", alpha=0.3, label="Gaussian model")
n = len(lifetimes)
low, high = normal_error_bounds(dist, n, qs)
plt.fill_between(qs, low, high, lw=0, color="gray", alpha=0.1)
surv = Cdf.from_seq(lifetimes)
surv.plot(ls="--", label="Data")
decorate(
xlabel="Duration (hours)",
ylabel="CDF",
title="Distribution of light bulb lifetimes",
)
```

```
from scipy.stats import gaussian_kde
xs = np.linspace(0, 3000, 201)
ys = gaussian_kde(lifetimes).evaluate(xs)
density = Pmf(ys, xs, name='')
density.normalize()
```

```
0.06666298081917303
```

```
mean = lifetimes.mean()
mean
```

```
1413.84
```

```
def decorate_kde(title=''):
decorate(
xlabel="Duration (hours)",
ylabel="Density",
xlim=[0, 3000],
yticks=[],
title=title
)
```

```
def plot_density(density):
mean = density.mean()
density.plot()
plt.axvline(mean, ls=':')
# plt.text(mean-20, 0.002, f'mean {mean:0.0f}', ha='center')
decorate_kde()
```

```
plot_density(density)
savefig('nbue1a.png')
```

```
def plot_conditional(density, age):
left = density[density.qs < age]
right = density[density.qs >= age]
left.plot()
right.plot(color='C0')
plt.fill_between(right.index, 0, right.values, alpha=0.1)
decorate_kde()
return right
```

```
age = 1000
conditional = plot_conditional(density, age)
savefig('nbue1b.png')
```

```
shifted = Pmf(conditional.values, conditional.index - age)
shifted.normalize()
plot_density(shifted)
savefig('nbue1c.png')
```

```
age = 500
conditional = plot_conditional(shifted, age)
savefig('nbue1d.png')
```

```
shifted = Pmf(conditional.values, conditional.index - age)
shifted.normalize()
plot_density(shifted)
savefig('nbue1e.png')
```

The shaded area shows how much variation we expect from the Gaussian model. Except for one unusually long-lasting bulb, the lifetimes fall within the bounds, which shows that the data are consistent with the model.

In this dataset, the average lifetime for a new light bulb is about 1414 hours. For a bulb that has been used for 1000 hours, the average lifetime is higher, about 1495 hours; however, since it has already burned 1000 hours, its average remaining lifetime is only 495 hours. So we would rather have the new bulb.

```
np.mean(lifetimes)
```

```
1413.84
```

```
lifetimes[lifetimes > 1000].mean() - 1000
```

```
495.2558139534883
```

```
lifetimes[lifetimes > 1500].mean() - 1500
```

```
213.391304347826
```

We can do the same calculation for a range of elapsed times from 0 to 2568 hours (the lifespan of the oldest bulb). At each point in time, \(t\), we can compute the average lifetime for bulbs that survive past \(t\) and the average remaining lifetime we expect.

```
def remaining_lifetimes_seq(seq, qs=None):
"""Compute average remaining lifetimes.
seq: sequence of quantities
qs: ages to evaluate
returns: Series
"""
if qs is None:
qs = np.linspace(0, seq.max(), 200)
series = pd.Series(index=qs, dtype=float)
for q in qs:
conditional = Pmf.from_seq(seq[seq >= q])
if conditional.sum() > 0:
conditional.normalize()
series[q] = conditional.mean() - q
return series.dropna()
```

To compute error bounds, we’ll use resampling.

```
def resample(xs):
"""Resample a sequence of quantities.
xs: quantities
returns: sample
"""
return np.random.choice(xs, size=len(xs), replace=True)
```

```
mu, sigma = lifetimes.mean(), lifetimes.std()
dist = norm(mu, sigma)
n = len(lifetimes)
qs = pmf_lifetimes.index
```

The following loop generates multiple resamplings and compute the remaining lifetimes for each.

```
np.random.seed(17)
from scipy.stats import gaussian_kde
kde = gaussian_kde(lifetimes)
res = []
qs = np.arange(0, 2550, 50)
for i in range(201):
sample = kde.resample(len(lifetimes))
rem = remaining_lifetimes_seq(sample, qs)
res.append(rem)
```

```
x = 1000
y = lifetimes[lifetimes > 1000].mean() - x
x, y
```

```
(1000, 495.2558139534883)
```

The following two functions compile the results and compute error bounds.

```
def percentile_rows(series_seq, ps):
"""Computes percentiles from aligned series.
series_seq: list of sequences
ps: cumulative probabilities
returns: Series of x-values, NumPy array with selected rows
"""
# NOTE: This is a specialized version of this function that
# fills NaNs with 0
df = pd.concat(series_seq, axis=1).fillna(0)
xs = df.index
array = df.values.transpose()
array = np.sort(array, axis=0)
nrows, ncols = array.shape
ps = np.asarray(ps)
indices = (ps * nrows).astype(int)
rows = array[indices]
return xs, rows
```

```
def plot_percentiles(series_seq, label=None, **options):
"""Plot the low, median, and high percentiles.
series_seq: sequence of Series
ps: percentiles to use for low, medium and high
label: string label for the median line
options: options passed plt.plot and plt.fill_between
"""
ps = [0.05, 0.5, 0.95]
xs, rows = percentile_rows(series_seq, ps)
low, med, high = rows
plt.plot(xs, med, alpha=0.5, label=label, **options)
plt.fill_between(xs, low, high, linewidth=0, alpha=0.2, **options)
```

```
ps = [0.5]
xs, rows = percentile_rows(res, ps)
(med,) = rows
x = 1000
y = pd.Series(med, xs)[x]
```

The following figure shows the result.

```
plt.plot([0, x, x], [y, y, 0], ":", color="gray", alpha=0.5)
plot_percentiles(res, color="C0")
decorate(
xlabel="Elapsed time in hours",
ylabel="Average remaining lifetime",
# title="Average remaining lifetime of light bulbs",
)
savefig('nbue2.png')
```

The x-axis shows elapsed times since the installation of a hypothetical light bulb. The y-axis shows the average remaining lifetime. In this example, a bulb that has been burning for 1000 hours is expected to last about 495 hours more, as indicated by the dotted lines.

So a new bulb generally lasts longer than a used bulb.

```
(lifetimes > 1450).sum(), (lifetimes > 1710).sum(), (lifetimes > 1800).sum()
```

```
(25, 11, 5)
```

## Any day now#

Now let’s consider the question posed by the imaginary obstetrician at lunch. Suppose you visit a maternity ward and meet women who are starting their 39th, 40th, and 41st week of pregnancy. Which one do you think will deliver first?

To answer this question, we need to know the distribution of gestation times, which we can get from the National Survey of Family Growth (NSFG), which was the source of the birth weights in Chapter~\ref{extremes-outliers-and-goats}. I gathered data collected between 2002 and 2017, which includes information about 43,939 live births.

I’ve made a subset of the data available in an HDF file.

```
download(DATA_PATH + "nsfg.hdf5")
```

```
# data cleaning is in FirstLateNSFG repo
# includes data from 2002, 2006-2010, 2011-2013, 2013-2015, 2015-2017
nsfg = pd.read_hdf("nsfg.hdf5", "nsfg")
nsfg.shape
```

```
(62539, 7)
```

```
nsfg.head()
```

pregend1 | nbrnaliv | prglngth | outcome | birthord | finalwgt | cycle | |
---|---|---|---|---|---|---|---|

0 | 6.0 | 1.0 | 40 | 1 | 1.0 | 19877.457610 | 10 |

1 | 1.0 | NaN | 14 | 4 | NaN | 19877.457610 | 10 |

2 | 6.0 | 1.0 | 39 | 1 | 2.0 | 19877.457610 | 10 |

3 | 6.0 | 1.0 | 39 | 1 | 1.0 | 4221.017695 | 10 |

4 | 6.0 | 1.0 | 39 | 1 | 2.0 | 4221.017695 | 10 |

I selected live births with pregnancy lengths less than 45 weeks (the outliers are probably errors).

```
live = nsfg["outcome"] == 1
reasonable = nsfg["prglngth"] < 45
lengths = nsfg.loc[live & reasonable, "prglngth"]
pmf_length = Pmf.from_seq(lengths)
len(lengths)
```

```
43246
```

```
(lengths < 28).mean() * 100
```

```
0.9272533875965407
```

```
gt = pmf_length.qs >= 28
pmf_length_gt = Pmf(pmf_length[gt] * 100)
```

The following figure shows the distribution of their durations (except for the 1% of babies born prior to 28 weeks).

```
pmf_length_gt.plot(color="C1", marker="o", alpha=0.4, label='')
decorate(
xlabel="Pregnancy length (weeks)",
ylabel="Percentage",
# title="Distribution of pregnancy length",
xticks=np.arange(28, 46, 2),
)
savefig('nbue3.png')
```

About 41% of these births were during the 39th week of pregnancy, and another 18% during the 40th week.

```
lengths.mean(), pmf_length_gt[39], pmf_length_gt[40]
```

```
(38.505341534477175, 42.357674698237986, 17.400453221107153)
```

I’ll use resampling to compute error bounds. Since the sampling weights are specific to each cycle, I’ll do the resampling within the cycles.

```
for name, group in nsfg.groupby("cycle"):
print(name, len(group))
```

```
6 13593
7 20492
8 9543
9 9358
10 9553
```

```
def sample_by_cycle(df):
sample = df.groupby("cycle").sample(frac=1, replace=True, weights="finalwgt")
return sample
```

The following loop computes the remaining lifetimes for each resampling.

```
qs = np.arange(36, 44)
mrt_seq = []
for i in range(101):
sample = sample_by_cycle(nsfg)
live = sample["outcome"] == 1
reasonable = sample["prglngth"] < 50
lengths = sample.loc[live & reasonable, "prglngth"]
mrt = remaining_lifetimes_seq(lengths, qs)
mrt_seq.append(mrt)
```

```
def decorate_pregnancy(title=''):
decorate(
xlabel="Week of pregnancy",
ylabel="Average remaining time",
title=title,
ylim=[0, 3.5],
)
```

```
# the previous version of plot_percentiles was specialized;
# for this analysis we want the general one (not that it matters)
from utils import plot_percentiles
plot_percentiles(mrt_seq, color="C1")
#decorate_pregnancy(title="Average remaining duration of pregancy")
decorate_pregnancy(title="")
savefig('nbue4.png')
```

Between weeks 36 and 39, the curve behaves as we expect: as time goes on, we get closer to the finish line. For example, at the beginning of week 36, the average remaining time is 3.2 weeks. At the beginning of week 37, it is down to 2.3 weeks. So far, so good: a week goes by and the finish line is just about one week closer.

But then, cruelly, the curve levels off. At the beginning of week 39, the average remaining time is 0.68 weeks, so the end seems near. But if we get to the beginning of week 40, and the baby has not been born, the average remaining time is 0.63 weeks. A week has passed and the finish line is only 8 hours closer.

And if we get to the beginning of week 41, and the baby has not been born, the average remaining time is 0.59 weeks. Another week has passed and the finish line is only 7 hours closer!

Here’s the computation for the average remaining time at the end of each week.

```
from utils import percentile_rows
xs, ys = percentile_rows(mrt_seq, [0.025, 0.5, 0.975])
df = pd.DataFrame(np.transpose(ys), xs, columns=[2.5, 50, 97.5])
df
```

2.5 | 50.0 | 97.5 | |
---|---|---|---|

36 | 3.142822 | 3.153736 | 3.165931 |

37 | 2.292390 | 2.303369 | 2.316928 |

38 | 1.465161 | 1.474933 | 1.488573 |

39 | 0.669659 | 0.679651 | 0.692118 |

40 | 0.612146 | 0.629809 | 0.647226 |

41 | 0.560671 | 0.586682 | 0.611351 |

42 | 0.320918 | 0.352373 | 0.382273 |

43 | 0.220000 | 0.271242 | 0.340611 |

```
df[50].diff().loc[[40, 41]] * 7 * 24
```

```
40 -8.373551
41 -7.245228
Name: 50.0, dtype: float64
```

```
from scipy.interpolate import interp1d
from empiricaldist import Hazard
series = df[50]
interp = interp1d(series.index, series.values, fill_value="extrapolate")
x = 39
for letter in 'abcdef':
y = interp(x)
print(x, y*7)
plt.figure()
plot_percentiles(mrt_seq, color="C1")
plt.plot([36, x, x], [y, y, 0], ":", color="gray", alpha=0.5)
plt.text(37, y+0.1, f'{y*7:0.2f} days')
decorate_pregnancy(title='')
savefig(f'nbue4{letter}.png')
x += y
```

```
39 4.757559429477013
39.679651347068145 4.520430465974382
40.32542712792163 4.310420075240543
40.941201424384566 4.124527367596406
41.530419619755484 3.2367992226202267
41.99281950870123 2.4783853029123057
```

## Cancer survival times#

Finally, let’s consider the surprising result reported by the imaginary oncologist at lunch: for many cancers, a patient who has survived a year after diagnosis is expected to live longer than someone who has just been diagnosed.

To demonstrate and explain this result, I’ll use data from the the Surveillance, Epidemiology, and End Results (SEER) program, which is run by the U.S. National Institutes of Health (NIH). Starting in 1973, SEER has collected data on cancer cases from registries in several locations in the United States. In the most recent datasets, these registries cover about one third of the U.S. population.

NOTE: I am definitely not allowed to redistribute SEER data, so the data in this example in synthetic; that is, I have generated random data that is statistically similar to the actual data I reported in the book. The results below will differ from what’s in the book, but the conclusions are qualitatively similar.

```
filename = "brain.hdf"
download(DATA_PATH + filename)
```

```
brain = pd.read_hdf("brain.hdf", "brain")
brain.shape
```

```
(16202, 2)
```

```
brain.head()
```

duration | observed | |
---|---|---|

0 | 4.50 | True |

1 | 15.75 | True |

2 | 10.75 | True |

3 | 6.25 | True |

4 | 29.50 | True |

From the SEER data I selected the 16,202 cases of glioblastoma, diagnosed between 2000 and 2016, where the survival times are available. We can use this data to infer the distribution of survival times, but first we have to deal with a statistical obstacle: some of the patients are still alive, or were alive the last time they appeared in a registry.

For these cases, the time until their death is unknown, which is a good thing. However, in order to work with data like this, we need a special method, called Kaplan-Meier estimation, to compute the distribution of lifetimes.

`duration`

is the observed duration in months, for both complete and ongoing cases.

```
duration = brain["duration"]
duration.describe()
```

```
count 16202.000000
mean 13.457783
std 19.088951
min 0.250000
25% 3.250000
50% 7.750000
75% 16.500000
max 201.750000
Name: duration, dtype: float64
```

`observed`

is a boolean that indicates whether each case is complete or ongoing.

```
observed = brain["observed"]
observed.sum()
```

```
14548
```

```
(duration[observed] == 0).mean(), (duration[~observed] == 0).mean()
```

```
(0.0, 0.0)
```

As we should expect, the distributions are different for the complete and ongoing cases.

```
Cdf.from_seq(duration[observed]).plot(label="complete")
Cdf.from_seq(duration[~observed]).plot(label="ongoing")
decorate(xscale="linear")
```

With Kaplan-Meier estimation, we can use all cases to estimate the survival function.

```
import lifelines
from empiricaldist import Surv
def km_fit(duration, observed):
"""Compute a survival function by Kaplan-Meier estimation.
duration: sequence of durations for complete and incomplete cases
observed: sequence booleans indicating whether each case is complete
returns: tuple of Surv (estimate, low, high)
"""
fit = lifelines.KaplanMeierFitter().fit(duration, observed)
surv_km = Surv(fit.survival_function_["KM_estimate"])
ci_fit = fit.confidence_interval_
surv_low = Surv(ci_fit["KM_estimate_lower_0.95"])
surv_high = Surv(ci_fit["KM_estimate_upper_0.95"])
return surv_km, surv_low, surv_high
```

```
surv_km, surv_low, surv_high = km_fit(duration, observed)
surv_km.head()
```

probs | |
---|---|

timeline | |

0.00 | 1.00000 |

0.25 | 0.97667 |

0.50 | 0.95130 |

```
def plot_fit(surv_km, surv_low, surv_high, **options):
"""Plot an estimated survival curve with error bounds.
surv_km: estimated Surv
surv_low, surv_high: lower and upper bounds
options: passed to Series.plot
"""
xs = surv_low.index
plt.fill_between(xs, surv_low, surv_high, color="gray", alpha=0.3)
surv_km.plot(**options)
```

The following figure shows the result on a log scale, plotted in the familiar form of a CDF and in a new form called a survival curve.

```
ax1 = plt.gca()
surv_km.plot(ls="--", color='C2', label="Survival curve")
h1, l1 = plt.gca().get_legend_handles_labels()
plt.ylabel("Probability of survival")
plt.xlabel("Survival time (months)")
plt.xscale("log")
labels = [1, 2, 5, 10, 20, 50, 100]
ticks = labels
plt.xticks(ticks, labels)
ax2 = plt.twinx()
ax1.tick_params(left=False, right=False)
cdf = surv_km.make_cdf() * 100
cdf.plot(color="gray", ls=":", label="CDF")
h2, l2 = plt.gca().get_legend_handles_labels()
plt.legend(h1 + h2, l1 + l2, loc="center right")
plt.ylabel("Percentile rank")
plt.title("Survival curve and CDF for patients with glioblastoma")
ax2.tick_params(left=False, right=False)
savefig('nbue5.png')
```

```
surv_km.plot(ls="--", color='C2', label="Survival curve")
plt.xscale("log")
labels = [1, 2, 5, 10, 20, 50, 100]
ticks = labels
plt.xticks(ticks, labels);
decorate(ylabel="Probability of survival",
xlabel="Survival time (months)")
savefig('nbue5a.png')
```

The survival curve shows the probability of survival past a given time on a scale from 0 to 1. It is the complement of the CDF, so as the CDF increases from left to right, the survival curve decreases. The two curves contain the same information; the only reason to use one or the other is convention. Survival curves are used more often in medicine and reliability engineering, CDFs in many other fields.

One thing that is apparent – from either curve – is that glioblastoma is a serious diagnosis. The median survival time after diagnosis is less than 9 months, and only about 16% of patients survive more than two years.

Please keep in mind that this curve lumps together people of different ages with different health conditions, diagnosed at different stages of disease over a period of about 16 years. Survival times depend on all of these factors, so this curve does not provide a prognosis for any specific patient. In particular, as treatment has gradually improved, the prognosis is better for someone with a more recent diagnosis. If you or someone you know is diagnosed with glioblastoma, you should get a prognosis from a doctor, based on specifics of the case, not from aggregated data in a book demonstrating basic statistical methods.

```
surv_km(12), surv_km(24), surv_km(24)/surv_km(12)
```

```
(array(0.38549537), array(0.16368232), 0.42460256589027323)
```

As we did with light bulbs and pregnancy lengths, we can use this distribution to compute the average remaining survival time for patients at each time after diagnosis. The following figure shows the result.

```
from utils import remaining_lifetimes_pmf
np.random.seed(17)
qs = np.linspace(0, 165, 101)
rem_seq = []
for i in range(101):
sample = brain.sample(frac=1, replace=True)
observed = sample["observed"]
duration = sample["duration"]
surv_km, surv_high, surv_low = km_fit(duration, observed)
rem = remaining_lifetimes_pmf(surv_km.make_pmf(), qs)
# rem.plot(color="gray", alpha=0.01)
rem_seq.append(rem)
```

```
def gray_box(x, x_max=170):
"""Make a gray box that spans the y-axis.
x_max: location of the right bound
"""
plt.axvspan(x, x_max, facecolor="gray", edgecolor='none', alpha=0.1)
plt.xlim([-5, x_max])
```

```
plot_percentiles(rem_seq, color='C2')
gray_box(125)
for x in [14, 32, 60]:
plt.axvline(x, color="gray", ls=":", alpha=0.5)
decorate(
xlabel="Time since diagnosis (months)",
ylabel="Average remaining survival",
title="Average remaining survival for patients with glioblastoma",
)
savefig('nbue6.png')
```

```
qs = np.arange(0, 165)
series = remaining_lifetimes_pmf(surv_km.make_pmf(), qs)
milestones = series[[0, 14, 32, 60]]
milestones
```

```
0 13.773990
14 17.069424
32 24.735119
60 27.839032
dtype: float64
```

```
ps = [0.05, 0.5, 0.95]
xs, rows = percentile_rows(rem_seq, ps)
low, med, high = rows
```

```
for i, letter in enumerate('abcd'):
plt.figure()
plt.plot(xs, med, alpha=0.7, label='', color='C2')
gray_box(125)
for x in milestones.cumsum()[:i+1]:
print(x)
plt.axvline(x, color="gray", ls=":", alpha=0.5)
decorate(
xlabel="Time since diagnosis (months)",
ylabel="Average remaining survival",
)
savefig(f'nbue6{letter}.png')
```

```
13.773989705091937
13.773989705091937
30.843413403818417
13.773989705091937
30.843413403818417
55.578532081470954
13.773989705091937
30.843413403818417
55.578532081470954
83.41756386272576
```

At the time of diagnosis, the average survival time is about 14 months. That is certainly a bleak prognosis, but there is some good news to follow. If a patient survives the first 14 months, we expect them to survive another 18 months, on average. If they survive those 18 months, for a total of 32, we expect them to survive another 28 months. And if they survive those 28 months, for a total of 60 months (five years), we expect them to survive another 35 months (almost three years). The vertical lines indicate these milestones. It’s like running a race where the finish line keeps moving, and the farther you go, the faster it retreats.

If you hear that the average survival time after diagnosis is 14 months, you might imagine a Gaussian distribution where 14 months is the most common value and an equal number of patients live more or less than 14 months. But that would be a very misleading picture.

To show how bad it would be, I chose a Gaussian distribution that matches the distribution of survival times as well as possible – which is not very well – and used it to compute average remaining survival times.

```
from utils import make_normal_model
pmf_km = surv_km.make_pmf()
pmf_normal = make_normal_model(pmf_km)
pmf_km.plot(label='data')
pmf_normal.plot(label='model')
decorate(
xlabel="Time since diagnosis (months)",
ylabel="PMF",
title="The normal model does not fit the data",
)
```

```
series = remaining_lifetimes_pmf(pmf_normal)
series.tail()
```

```
197.694724 0.972382
198.708543 0.829894
199.722362 0.630459
200.736181 0.358620
201.750000 0.000000
dtype: float64
```

As a result, the remaining lifetimes from the model don’t resemble the remaining lifetimes from the data.

```
plot_percentiles(rem_seq, color='C2', label="Actual data")
series.plot(color="C4", ls=":", label="Gaussian model")
gray_box(125)
decorate(
xlabel="Time since diagnosis (months)",
ylabel="Average remaining time",
title="Average remaining time compared to Gaussian model",
loc="upper right",
)
```

With the Gaussian model, the average remaining survival time starts around 20 months, drops quickly at first, and levels off around 5 months. So it behaves nothing like the actual averages.

On the other hand, if your mental model of the distribution is lognormal, you would get it about right. To demonstrate, I chose a lognormal distribution that fits the actual distribution of survival times and used it to compute average remaining lifetimes.

Here’s the function that fits the distribution.

```
from scipy.optimize import least_squares
def fit_lognormal(surv, xs=None):
"""Fit a lognormal distribution to a survival function.
surv: Surv object
xs: places to evaluate the survival function
returns: Surv
"""
def error_func(params):
mu, sigma = params
# just fit over the range from 0 to 120 months
xs = np.logspace(0, np.log10(120))
ps = norm.sf(np.log(xs), mu, sigma)
error = ps - surv_km(xs)
return error
pmf = surv.make_pmf()
pmf.normalize()
params = pmf.mean(), pmf.std()
res = least_squares(error_func, x0=params, xtol=1e-3)
assert res.success
mu, sigma = res.x
xs = surv.index if xs is None else xs
ps = norm.sf(np.log(xs), mu, sigma)
return Surv(ps, xs)
```

```
surv_lognormal = fit_lognormal(surv_km)
pmf_lognormal = surv_lognormal.make_pmf()
```

```
series = remaining_lifetimes_pmf(pmf_lognormal)
```

The following figure shows the result.

```
plot_percentiles(rem_seq, color='C2', label="Actual data")
series.plot(color="C4", ls=":", label="Lognormal model")
gray_box(125)
decorate(
xlabel="Time since diagnosis (months)",
ylabel="Average remaining time",
title="Average remaining time compared to lognormal model",
)
plt.legend(loc="lower left")
None
```

During the first 24 months, the model is a little too optimistic, and after 120 months it is much too optimistic. But the lognormal model gets the shape of the curve right: if your mental model of the distribution is lognormal, you would have a reasonably accurate understanding of the situation.

## Life Expectancy At Birth#

In 2012, a team of demographers at the University of Southern California estimated life expectancy for people born in Sweden in the early 1800s and 1900s. They chose Sweden because it “has the deepest historical record of high-quality [demographic] data.”

For ages from 0 to 91 years, they estimated the mortality rate, which is the fraction of people at each age who die.

I used an online graph digitizer to get the data from the figure in their paper and store it in a CSV file.

```
filename = "mortality_rates_beltran2012.csv"
download(DATA_PATH + filename)
```

```
mortality = pd.read_csv(filename, header=[0, 1])
mortality.head()
```

1905 | 1800 | |||
---|---|---|---|---|

X | Y | X | Y | |

0 | 0.056633 | -2.685744 | 0.066515 | -1.541222 |

1 | 0.175223 | -2.866213 | 0.224635 | -1.692583 |

2 | 0.293812 | -3.043189 | 0.362990 | -1.890516 |

3 | 0.412402 | -3.216672 | 0.481580 | -2.052938 |

4 | 0.530992 | -3.390155 | 0.623887 | -2.239811 |

```
def make_series(x, y):
"""Make a pandas Series.
x: values for the index
y: values of the series
returns: Series
"""
return pd.Series(y, x).dropna()
```

```
mort1800 = make_series(mortality["1800", "X"], mortality["1800", "Y"].values)
mort1800.head()
```

```
(1800, X)
0.066515 -1.541222
0.224635 -1.692583
0.362990 -1.890516
0.481580 -2.052938
0.623887 -2.239811
dtype: float64
```

```
mort1905 = make_series(mortality["1905", "X"], mortality["1905", "Y"].values)
mort1905.head()
```

```
(1905, X)
0.056633 -2.685744
0.175223 -2.866213
0.293812 -3.043189
0.412402 -3.216672
0.530992 -3.390155
dtype: float64
```

```
np.exp(mort1800.iloc[0]), np.exp(mort1905.iloc[0])
```

```
(0.21411938757305218, 0.06817045483264002)
```

The following function interpolates the data from the figure and puts it in a Hazard function.

```
from scipy.interpolate import interp1d
from empiricaldist import Hazard
def make_hazard(series):
"""Make a Hazard function by interpolating a Series.
series: Series
returns: Hazard
"""
interp = interp1d(series.index, series.values, fill_value="extrapolate")
xs = np.arange(0, 108)
ys = np.exp(interp(xs))
return Hazard(ys, xs)
```

```
haz1800 = make_hazard(mort1800)
haz1905 = make_hazard(mort1905)
```

The following figure shows the results for two cohorts: people born between 1800 and 1810, and people born between 1905 and 1915.

```
haz1800[:91].plot(label="Born 1800–1810", ls="--")
haz1905[:91].plot(label="Born 1905–1915")
decorate(
xlabel="Age", ylabel="Mortality rate", title="Historic mortality rates in Sweden"
)
savefig('nbue7.png')
```

The notable feature is that mortality rates were lower for the later cohort at every age.

```
haz1800[0], haz1905[0]
```

```
(0.2281961996457684, 0.07430619990025686)
```

```
haz1800.min(), haz1800.idxmin(), haz1905[14]
```

```
(0.004624153294009557, 14, 0.0019040788545011136)
```

```
haz1800[80], haz1905[80]
```

```
(0.13156208289057392, 0.06441158462898641)
```

```
pmf1800 = haz1800.make_pmf()
pmf1800.normalize()
```

```
0.9999999998644409
```

```
pmf1905 = haz1905.make_pmf()
pmf1905.normalize()
```

```
0.9999711879585407
```

```
pmf1800.mean(), pmf1905.mean()
```

```
(36.473535530889244, 65.67277539967878)
```

The following figure shows life expectancy as a function of age for people born in Sweden around 1800 and 1905.

```
qs = np.arange(0, 91)
series1905 = remaining_lifetimes_pmf(pmf1905, qs)
series1905.plot(ls="--", color="C0", label="Sweden 1905")
series1800 = remaining_lifetimes_pmf(pmf1800, qs)
series1800.plot(ls=":", color="C0", label="Sweden 1800")
decorate(
xlabel="Age", ylabel="Average remaining lifetime", title="Remaining lifetime vs age"
)
savefig('nbue8a.png')
```

```
def draw_lines(pmf, series):
x = pmf.mean()
y = series[int(round(x))]
plt.plot([0, x, x], [y, y, 0], ":", color="gray", alpha=0.5)
```

```
qs = np.arange(0, 91)
series1905 = remaining_lifetimes_pmf(pmf1905, qs)
draw_lines(pmf1905, series1905)
series1905.plot(ls="--", color="C0", label="Sweden 1905")
series1800 = remaining_lifetimes_pmf(pmf1800, qs)
draw_lines(pmf1800, series1800)
series1800.plot(ls=":", color="C0", label="Sweden 1800")
decorate(
xlabel="Age", ylabel="Average remaining lifetime", title="Remaining lifetime vs age"
)
savefig('nbue8b.png')
```

In both cohorts, used was better than new, at least for the first few years of life.

```
series1800[0], series1800.idxmax(), series1800.max()
```

```
(36.473535530889244, 5, 49.9171982587444)
```

```
series1905[0], series1905.idxmax(), series1905.max()
```

```
(65.67277539967881, 2, 70.18860681489375)
```

```
series1800[0], series1800[36], series1800[65]
```

```
(36.473535530889244, 29.241600197695945, 11.018435234417836)
```

```
series1905[5]
```

```
68.80456362586426
```

```
series1905[0], series1905[66], series1905[82]
```

```
(65.67277539967881, 16.058634133594552, 6.471530454841982)
```

## Child mortality#

The following figure shows the percentage of children who die before age 5 for four geographical regions, from 1900 to 2019. These data were combined from several sources by Gapminder, a foundation based in Sweden that “promotes sustainable global development […] by increased use and understanding of statistics.”

Documentation of the data is at https://www.gapminder.org/data/documentation/gd005

The data is in this online spreadsheet, which I downloaded in August 2022.

```
filename = "GM-ChildMortality-Dataset-v11.xlsx"
download(DATA_PATH + filename)
```

```
gm = pd.read_excel(filename, sheet_name="data-for-regions-by-year")
gm.tail()
```

geo | name | time | Child mortality | |
---|---|---|---|---|

1199 | americas | The Americas | 2096 | 2.70 |

1200 | americas | The Americas | 2097 | 2.67 |

1201 | americas | The Americas | 2098 | 2.63 |

1202 | americas | The Americas | 2099 | 2.60 |

1203 | americas | The Americas | 2100 | 2.59 |

```
grouped = gm.groupby("geo")
line_styles = {"africa": ":", "asia": "-", "americas": "--", "europe": "-."}
for x in [1918, 1942, 1960]:
plt.axvline(x, color="gray", alpha=0.2)
for name in ["africa", "asia", "americas", "europe"]:
group = grouped.get_group(name)
series = pd.Series(group["Child mortality"].values / 10, group["time"])
series.loc[1900:2019].plot(label=name.title(), ls=line_styles[name])
decorate(
xlabel="Year",
ylabel="Percent who die before age 5",
title="Child mortality by region",
)
```

In every region, child mortality has decreased consistently and substantially.

I collected recent mortality data from the Global Health Observatory of the World Health Organization (WHO). For people born in 2019, we don’t know what their future lifetimes will be, but we can estimate it if we assume that the mortality rate in each age group will not change over their lifetimes.

```
def read_hazard(filename):
"""Read age-specific death rates and make a Hazard object.
filename: string
returns: Hazard
"""
df = pd.read_csv(filename, header=[1])
query = 'Indicator == "nMx - age-specific death rate between ages x and x+n"'
nmx = df.query(query)
index = nmx["Age Group"].values
ages = np.arange(-5, 90, 5)
age_series = pd.Series(ages, index)
age_series.iloc[[0, 1]] = [0, 1]
rate = nmx["Both sexes"]
rate.index = age_series.values
hazard = rate.reindex(np.arange(111)).ffill()
return Hazard(hazard)
```

```
download(DATA_PATH + "life_table_sweden_2019.csv")
download(DATA_PATH + "life_table_nigeria_2019.csv")
```

```
filename = "life_table_sweden_2019.csv"
hazard = read_hazard(filename)
pmf_swe = hazard.make_pmf()
pmf_swe.normalize()
pmf_swe.mean()
```

```
81.6599453119953
```

```
filename = "../data/life_table_nigeria_2019.csv"
hazard = read_hazard(filename)
hazard.head()
```

probs | |
---|---|

0 | 0.077958 |

1 | 0.011613 |

2 | 0.011613 |

```
pmf_nga = hazard.make_pmf()
pmf_nga.normalize()
pmf_nga.mean()
```

```
61.59620754538476
```

Based on that simplification, the following figure shows average remaining lifetime as a function of age for Sweden and Nigeria in 2019, compared to Sweden in 1905.

```
rem_swe = remaining_lifetimes_pmf(pmf_swe)
rem_nga = remaining_lifetimes_pmf(pmf_nga)
rem_swe[0]
```

```
81.6599453119953
```

```
rem_swe.loc[:91].plot(label="Sweden 2019", ls="-")
series1905.plot(ls="--", color="C0", label="Sweden 1905")
rem_nga.loc[:91].plot(ls="-.", label="Nigeria 2019", color="C2")
decorate(
xlabel="Age",
ylabel="Average remaining lifetime",
title="Average remaining lifetime vs age",
)
```

```
rem_swe.loc[:91].plot(label="Sweden 2019", ls="-")
series1905.plot(ls="--", color="C0", label="Sweden 1905")
series1800.plot(ls=":", color="C0", label="Sweden 1800")
decorate(
xlabel="Age", ylabel="Average remaining lifetime", title="Remaining lifetime vs age"
)
savefig('nbue9.png')
```

```
rem_nga.head()
```

```
0.000000 61.596208
0.552764 66.256747
1.105528 66.478013
1.658291 65.925249
2.211055 66.143866
dtype: float64
```

## The Immortal Swede#

Going back to the data from Sweden, the following figure shows the mortality rate for each age group, updated every ten years from 2000 to 2019.

```
download(DATA_PATH + "life_table_sweden.csv")
```

```
columns = ["Indicator", "Age Group", 2019, 2015, 2010, 2005, 2000]
df = pd.read_csv("life_table_sweden.csv", header=[1])
df.columns = columns
df.head()
```

Indicator | Age Group | 2019 | 2015 | 2010 | 2005 | 2000 | |
---|---|---|---|---|---|---|---|

0 | nMx - age-specific death rate between ages x a... | <1 year | 0.002042 | 0.002452 | 0.002579 | 0.002448 | 0.003451 |

1 | nMx - age-specific death rate between ages x a... | 1-4 years | 0.000121 | 0.000127 | 0.000150 | 0.000211 | 0.000123 |

2 | nMx - age-specific death rate between ages x a... | 5-9 years | 0.000061 | 0.000068 | 0.000063 | 0.000102 | 0.000097 |

3 | nMx - age-specific death rate between ages x a... | 10-14 years | 0.000081 | 0.000110 | 0.000090 | 0.000101 | 0.000134 |

4 | nMx - age-specific death rate between ages x a... | 15-19 years | 0.000217 | 0.000230 | 0.000271 | 0.000262 | 0.000355 |

```
index = df["Age Group"].values
ages = np.arange(-5, 90, 5)
age_series = pd.Series(ages, index)
age_series.iloc[[0, 1]] = [0, 1]
age_series
```

```
<1 year 0
1-4 years 1
5-9 years 5
10-14 years 10
15-19 years 15
20-24 years 20
25-29 years 25
30-34 years 30
35-39 years 35
40-44 years 40
45-49 years 45
50-54 years 50
55-59 years 55
60-64 years 60
65-69 years 65
70-74 years 70
75-79 years 75
80-84 years 80
85+ years 85
dtype: int64
```

```
df.drop(columns=["Indicator", "Age Group"], inplace=True)
df.index = age_series.values
df
```

2019 | 2015 | 2010 | 2005 | 2000 | |
---|---|---|---|---|---|

0 | 0.002042 | 0.002452 | 0.002579 | 0.002448 | 0.003451 |

1 | 0.000121 | 0.000127 | 0.000150 | 0.000211 | 0.000123 |

5 | 0.000061 | 0.000068 | 0.000063 | 0.000102 | 0.000097 |

10 | 0.000081 | 0.000110 | 0.000090 | 0.000101 | 0.000134 |

15 | 0.000217 | 0.000230 | 0.000271 | 0.000262 | 0.000355 |

20 | 0.000407 | 0.000457 | 0.000458 | 0.000472 | 0.000511 |

25 | 0.000549 | 0.000575 | 0.000498 | 0.000500 | 0.000499 |

30 | 0.000598 | 0.000563 | 0.000505 | 0.000498 | 0.000524 |

35 | 0.000631 | 0.000604 | 0.000570 | 0.000710 | 0.000827 |

40 | 0.000836 | 0.000860 | 0.000916 | 0.001102 | 0.001240 |

45 | 0.001285 | 0.001348 | 0.001581 | 0.001838 | 0.002161 |

50 | 0.002087 | 0.002457 | 0.002690 | 0.003153 | 0.003394 |

55 | 0.003508 | 0.004089 | 0.004410 | 0.005092 | 0.005415 |

60 | 0.006164 | 0.006585 | 0.007485 | 0.007953 | 0.008854 |

65 | 0.010054 | 0.010798 | 0.011625 | 0.013320 | 0.014281 |

70 | 0.016627 | 0.018914 | 0.019199 | 0.021721 | 0.024041 |

75 | 0.028767 | 0.031977 | 0.034546 | 0.037426 | 0.041425 |

80 | 0.053991 | 0.058500 | 0.062250 | 0.068268 | 0.074020 |

85 | 0.156954 | 0.156715 | 0.156301 | 0.161742 | 0.167294 |

```
from utils import Re70, Or70, Gr60, Bl50, Pu40
def plot_mortality(df):
"""Plot columns from a mortality DataFrame.
df: DataFrame
"""
df[2000].plot(lw=1, style="x-", color=Re70)
# df[2005].plot(lw=1, ms=4, style="^--", color=Or70)
df[2010].plot(lw=1, ms=4, style="s--", color=Gr60)
# df[2015].plot(lw=1, ms=4, style="*--", color=Bl50)
df[2019].plot(lw=1, ms=4, style="o-", color=Pu40)
decorate(xlabel="Age", ylabel="Mortality rate", yscale="log")
```

```
plot_mortality(df)
labels = ["1/10", "1/100", "1/1000", "1/10000"]
ticks = [1e-1, 1e-2, 1e-3, 1e-4]
plt.yticks(ticks, labels)
#decorate(title="Mortality rates in Sweden")
plt.tight_layout()
savefig('nbue10.png')
```

The straight-line increase after age 35 was described by Benjamin Gompertz in 1825, so this phenomenon is called the Gompertz Law. It is an empirical law, which is to say that it names a pattern we see in nature, but at this point we don’t have an explanation of why it’s true, or whether it is certain to be true in the future. Nevertheless, the data in this example fall in a remarkably straight line.

```
df.loc[[0, 50], 2019] * 100
```

```
0 0.204202
50 0.208688
Name: 2019, dtype: float64
```

```
df.loc[[10, 40, 65, 85], 2019] * 100
```

```
10 0.008100
40 0.083590
65 1.005406
85 15.695439
Name: 2019, dtype: float64
```

The previous figure also shows that mortality rates decreased between 2000 and 2019 in almost every age group. If we zoom in on the age range from 40 to 80, we can see the changes in adult mortality more clearly.

```
subset = df.loc[40:80]
plot_mortality(subset)
def add_yticks():
labels = ["1/10", "3/100", "1/100", "3/1000", "1/1000"]
ticks = [1e-1, 3e-2, 1e-2, 3e-3, 1e-3]
plt.yticks(ticks, labels)
add_yticks()
#decorate(title="Mortality rates for adults in Sweden")
plt.tight_layout()
savefig('nbue11.png')
```

In these age groups, the decreases in mortality have been remarkably consistent. By fitting a model to this data, we can estimate the rate of change as a function of both age and time. According to the model, as you move from one age group to the next, your mortality rate increases by about 11% per year. At the same time, for the reasons I just mentioned, the mortality rate in every age group decreases by about 2% per year.

```
data = subset.stack().reset_index()
data.columns = ["age", "year", "rate"]
data["log_rate"] = np.log10(data["rate"])
data.head()
```

age | year | rate | log_rate | |
---|---|---|---|---|

0 | 40 | 2019 | 0.000836 | -3.077846 |

1 | 40 | 2015 | 0.000860 | -3.065354 |

2 | 40 | 2010 | 0.000916 | -3.038094 |

3 | 40 | 2005 | 0.001102 | -2.957944 |

4 | 40 | 2000 | 0.001240 | -2.906429 |

```
import statsmodels.formula.api as smf
formula = "log_rate ~ age + year"
results = smf.ols(formula, data=data).fit()
results.summary()
```

Dep. Variable: | log_rate | R-squared: | 0.999 |
---|---|---|---|

Model: | OLS | Adj. R-squared: | 0.999 |

Method: | Least Squares | F-statistic: | 1.581e+04 |

Date: | Sun, 18 Feb 2024 | Prob (F-statistic): | 3.78e-61 |

Time: | 17:29:03 | Log-Likelihood: | 109.74 |

No. Observations: | 45 | AIC: | -213.5 |

Df Residuals: | 42 | BIC: | -208.1 |

Df Model: | 2 | ||

Covariance Type: | nonrobust |

coef | std err | t | P>|t| | [0.025 | 0.975] | |
---|---|---|---|---|---|---|

Intercept | 13.5134 | 0.964 | 14.017 | 0.000 | 11.568 | 15.459 |

age | 0.0446 | 0.000 | 176.799 | 0.000 | 0.044 | 0.045 |

year | -0.0091 | 0.000 | -19.003 | 0.000 | -0.010 | -0.008 |

Omnibus: | 2.326 | Durbin-Watson: | 1.177 |
---|---|---|---|

Prob(Omnibus): | 0.313 | Jarque-Bera (JB): | 2.195 |

Skew: | 0.507 | Prob(JB): | 0.334 |

Kurtosis: | 2.623 | Cond. No. | 5.95e+05 |

Notes:

[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

[2] The condition number is large, 5.95e+05. This might indicate that there are

strong multicollinearity or other numerical problems.

These results imply that the life expectancies we computed in the previous section are too pessimistic because they take into account only the first effect – the increase with age – and not the second – the decrease over time. So let’s see what happens if we include the second effect as well, that is, if we assume that mortality rates will continue to decrease.

The following figure shows the actual mortality rates for 2000 and 2019 again, along with predictions for 2040 and 2060.

```
ratios = 10**results.params
ratios["age"] - 1, ratios["year"] - 1
```

```
(0.10821618126745203, -0.02076728510750292)
```

The following function uses the results from the model to generate predictions – assuming future trends at current rates.

```
def pred_hazard(results, pred_frame):
"""Predict hazard rates
"""
hazard = Hazard(10 ** results.predict(pred_frame))
hazard.index = pred_frame["age"]
return hazard
```

```
pred_frame = pd.DataFrame(dtype=float)
pred_frame["age"] = np.arange(40, 105)
pred_frame["year"] = 2019
hazard_no_progress = pred_hazard(results, pred_frame)
```

```
pred_frame = pd.DataFrame(dtype=float)
pred_frame["age"] = np.arange(40, 120)
pred_frame["year"] = 2019 + pred_frame["age"] - 40
hazard_progress = pred_hazard(results, pred_frame)
```

```
def make_pred(subset, results, year):
"""Use a regression model to generate predictions.
subset: DataFrame
results: RegressionResults
year: int
returns: Hazard
"""
pred_frame = pd.DataFrame(dtype=float)
pred_frame["age"] = subset.index.values
pred_frame["year"] = year
pred = 10 ** results.predict(pred_frame)
pred.index = subset.index
return Hazard(pred)
```

```
subset[2000].plot(lw=1, style="x-", color=Re70)
subset[2019].plot(lw=1, ms=4, style="o-", color=Pu40)
pred = make_pred(subset, results, 2040)
pred.plot(style=":", lw=1, color="gray")
pred = make_pred(subset, results, 2060)
pred.plot(style=":", lw=1, color="gray")
for x, marker in zip([60, 80], ["^", "s"]):
y = hazard_progress[x]
plt.plot(x, y, marker, color="C0")
hazard_progress.loc[:80].plot(lw=1, label="With current progress")
decorate(xlabel="Age", ylabel="Mortality rate", yscale="log")
add_yticks()
#decorate(title="Mortality rate in Sweden with continued progress", loc="lower right")
plt.tight_layout()
savefig('nbue12.png')
```

The line labeled “With current progress” indicates the mortality rates we expect for a hypothetical Swede who was 40 in 2020. When they are 60, it will be 2040, so we expect them to have the mortality rate of a 60-year-old in 2040, indicated by a triangle. And when they are 80, it will be 2060, so we expect them to have the mortality rate of an 80-year-old in 2060, indicated by a square.

We can use these mortality rates to compute survival curves, as shown in the following figure.

```
hazard_no_progress.make_surv().plot(label="Without progress")
decorate(
xlabel="Age",
ylabel="Probability of survival",
#title="Survival curve, Sweden 2019",
)
savefig('nbue12a.png')
```

```
hazard_no_progress.make_surv().plot(label="Without progress")
hazard_progress.make_surv().plot(ls="--", label="With current progress")
decorate(
xlabel="Age",
ylabel="Probability of survival",
#title="Survival curves with and without progress",
)
savefig('nbue12b.png')
```

```
def hypothetical_results(b2):
"""Make RegressionResults with a counterfactual slope.
b2: slope
returns: RegressionResults
"""
results = smf.ols(formula, data=data).fit()
i1, a1, b1 = results.params
i2 = i1 + (b1 - b2) * 2019
results.params["Intercept"] = i2
results.params["year"] = b2
return results
```

```
def hazard_hypo_progress(b2, upper):
"""Compute a Hazard function with hypothetical progress.
b2: slope
upper: upper end of the qs
returns: Hazard
"""
hypo = hypothetical_results(b2)
pred_frame = pd.DataFrame(dtype=float)
pred_frame["age"] = np.arange(40, upper)
pred_frame["year"] = 2019 + pred_frame["age"] - 40
return pred_hazard(hypo, pred_frame)
```

```
b2 = results.params["year"] * 2
hazard_progress2 = hazard_hypo_progress(b2, 145)
```

```
b3 = results.params["year"] * 3
hazard_progress3 = hazard_hypo_progress(b3, 180)
```

```
b4 = results.params["year"] * 4
hazard_progress4 = hazard_hypo_progress(b4, 300)
```

```
results.params["age"] / results.params["year"]
```

```
-4.89621102405056
```

```
b49 = -results.params["age"]
hazard_progress49 = hazard_hypo_progress(b49, 8000)
```

```
hazard_progress.plot()
hazard_progress2.plot()
hazard_progress3.plot()
hazard_progress4.plot()
hazard_progress49.loc[:140].plot()
decorate(yscale="log")
```

```
hazard_no_progress.make_surv().plot(ls="-", label="Without progress")
hazard_progress.make_surv().plot(ls="--", label="With current progress")
hazard_progress2.make_surv().plot(ls=":", label="With progress x2")
decorate(
xlabel="Age",
ylabel="Probability of survival",
#title="Survival curves with different rates of progress",
)
savefig('nbue13.png')
```

The dashed line on the left shows the survival curve we expect if there is no further decrease in mortality rates; in that scenario, life expectancy at age 40 is 82 years and the probability of living to 100 is only 1.4%.

The solid line shows the survival curve if mortality continues to decrease at the current rate; in that case, life expectancy for the same 40-year-old is 90 years and the chance of living to 100 is 25%.

```
def mean(hazard):
pmf = hazard.make_pmf()
pmf.normalize()
return pmf.mean()
```

```
mean(hazard_no_progress), mean(hazard_progress), mean(hazard_progress2)
```

```
(81.88119575997098, 90.01328388938188, 102.72010650635067)
```

```
mean(hazard_progress3), mean(hazard_progress4)
```

```
(125.65619219578866, 185.38884205632206)
```

```
hazard_no_progress.make_surv()(100)
```

```
array(0.01353305)
```

```
hazard_progress.make_surv()(100)
```

```
array(0.24964025)
```

```
hazard_progress2.make_surv()([100, 120, 130, 140])
```

```
array([0.60379629, 0.17113289, 0.03569291, 0.00159352])
```

Finally, the dotted line on the right shows the survival curve if mortality decreases at twice the current rate. Life expectancy at age 40 would be 102 and the probability of living to 100 would be 60%.

The following figure shows the survival curve under the assumption that mortality rates, starting in 2019, decrease at four times the current rate.

```
hazard_progress4.make_surv().plot(color="C3", ls="-.", label="With progress ×4")
decorate(
xlabel="Age",
ylabel="Probability of survival",
# title="Survival curve with progress ×4",
)
```

```
hazard_no_progress.make_surv().plot(ls="-", label="Without progress")
hazard_progress.make_surv().plot(ls="--", label="With current progress")
hazard_progress2.make_surv().plot(ls=":", label="With progress ×2")
hazard_progress4.make_surv().plot(color="C3", ls="-.", label="With progress ×4")
decorate(
xlabel="Age",
ylabel="Probability of survival",
# title="Survival curve with 4x rate of progress",
)
savefig('nbue14.png')
```

In this scenario, some people live to be 300 years old! However, even with these optimistic assumptions, the shape of the curve is similar to what we get with slower rates of progress. And, as shown in the following figure, average remaining lifetimes decrease with almost the same slope.

```
qs = np.arange(40, hazard_progress4.qs.max())
rem_progress4 = remaining_lifetimes_pmf(hazard_progress4.make_pmf(), qs)
qs = np.arange(40, hazard_progress3.qs.max())
rem_progress3 = remaining_lifetimes_pmf(hazard_progress3.make_pmf(), qs)
qs = np.arange(40, 140)
rem_progress2 = remaining_lifetimes_pmf(hazard_progress2.make_pmf(), qs)
qs = np.arange(40, 120)
rem_progress = remaining_lifetimes_pmf(hazard_progress.make_pmf(), qs)
qs = np.linspace(40, 100)
rem_no_progress = remaining_lifetimes_pmf(hazard_no_progress.make_pmf(), qs)
```

```
rem_progress4.plot(color="C3", ls="-.", label="With progress ×4")
# rem_progress3.plot(ls=':', label='With progress ×3')
rem_progress2.plot(ls=":", label="With progress ×2")
rem_progress.plot(ls="-", label="With current progress")
rem_no_progress.plot(ls="--", label="Without progress")
decorate(
xlabel="Age",
ylabel="Average remaining lifetime",
# title="Average remaining lifetime with different rates of progress",
)
savefig('nbue15.png')
```

With faster progress, people live longer, but they still have the new-better-than-used property: as each year passes, they get closer to the grave, on average.

However, something remarkable happens if progress accelerates by a factor of 4.9. At that speed, the increase in mortality due to aging by one year is exactly offset by the decrease due to progress. That means that the probability of dying is the same from one year to the next, forever. The result is a survival curve that looks like this.

```
dashdotdot = (0, (3, 1, 1, 1, 1, 1))
hazard_no_progress.make_surv().plot(ls="-", label="Without progress")
hazard_progress.make_surv().plot(ls="--", label="With current progress")
hazard_progress2.make_surv().plot(ls=":", label="With progress ×2")
hazard_progress4.make_surv().plot(color="C3", ls="-.", label="With progress ×4")
hazard_progress49.make_surv().plot(
color="C4", ls=dashdotdot, label="With progress ×4.9"
)
decorate(
xlabel="Age",
ylabel="Probability of survival",
# title="Survival curves with progress ×4 and ×4.9 progress",
)
savefig('nbue15.png')
```

The difference between 4 and 4.9 is qualitative. At a factor of 4, the survival curve plummets toward an inevitable end; at a factor of 4.9, it extends out to ages that are beyond biblical.

```
lam = hazard_no_progress[40]
lam * 1e4
```

```
7.888830027209969
```

```
1 / lam, 2 / lam
```

```
(1267.6150919094762, 2535.2301838189524)
```

```
from scipy.stats import expon
dist = expon(scale=1 / lam)
dist.mean(), dist.median()
```

```
(1267.6150919094762, 878.6438269922893)
```

```
dist.sf([2000, 4000])
```

```
array([0.20643576, 0.04261572])
```

```
dist.ppf(1 - 1 / 7e9)
```

```
28735.7893553912
```

In the near future, if child mortality continues to decrease, everyone will be better new than used. But eventually, if adult mortality decreases fast enough, everyone will have the same remaining time, on average: new, used, or in between.

Probably Overthinking It

Copyright 2022 Allen Downey

The code in this notebook and `utils.py`

is under the MIT license.