Глава 13. Анализ выживаемости
Анализ выживаемости (survival analysis) — это способ описания того, как долго что-либо длится. Он часто используется для изучения продолжительности жизни человека, но также подходит для моделирования «выживаемости» механических и электронных компонентов, или, в более общем плане, интервала времени до наступления некоторого события, или даже расстояния в пространстве.
Начнем с простого примера — срока службы электрических лампочек, а затем рассмотрим другой, более существенный пример — возраст вступления в первый брак и то, как он изменился в США за последние 50 лет.
Функция выживания
Фундаментальным понятием в анализе выживаемости является функция выживания (survival function). Она представляет собой долю генеральной совокупности, которая сохраняется во времени дольше заданной продолжительности. В качестве первого примера мы рассчитаем функцию выживания для времени работы лампы накаливания.
Используем данные эксперимента, проведенного в 2007 году. Исследователи установили 50 новых ламп накаливания и оставили их постоянно включенными. Они проверяли лампочки каждые 12 часов и записывали срок службы каждой перегоревшей. Эксперимент продолжался, пока из строя не вышли все 50 лампочек. Инструкции по загрузке данных приведены в Jupyter-блокноте этой главы.
Считаем данные:
df = pd.read_csv("lamps.csv", index_col=0)
df.tail()
| i | h | f | K |
| 28 | 1812 | 1 | 4 |
| 29 | 1836 | 1 | 3 |
| 30 | 1860 | 1 | 2 |
| 31 | 1980 | 1 | 1 |
| 32 | 2568 | 1 | 0 |
В столбце h приведена продолжительность работы в часах (hours). В столбце f указано количество перегоревших лампочек для каждого значения h. Чтобы представить распределение срока службы, поместим эти значения в объект Pmf и нормализуем его:
from empiricaldist import Pmf
pmf_bulblife = Pmf(df["f"].values, index=df["h"])
pmf_bulblife.normalize()
50
С помощью метода make_cdf построим интегральную функцию распределения (ИФР), показывающую долю лампочек, вышедших из строя до или в момент времени h. Например, к моменту времени 1656 часов перегорает 78 % ламп:
cdf_bulblife = pmf_bulblife.make_cdf()
cdf_bulblife[1656]
0.7800000000000002
Функция выживания — это доля лампочек, перегорающих после какого-либо момента времени h, являющаяся дополнением (complement) к ИФР. Ее можно вычислить следующим образом:
complementary_cdf = 1 - cdf_bulblife
complementary_cdf[1656]
0.21999999999999975
После 1656 часов выходит из строя 22 % лампочек.
В библиотеке empiricaldist имеется класс Surv, реализующий функцию выживания. Объекты этого класса можно создавать при помощи метода make_surv:
surv_bulblife = cdf_bulblife.make_surv()
surv_bulblife[1656]
0.21999999999999997
Если построить графики ИФР и функции выживания, то станет видно, что они дополняют друг друга, то есть их сумма равна 1 при любых значениях h:
cdf_bulblife.plot(ls="--", label="ИФР")
surv_bulblife.plot(label="Функция выживания")
decorate(xlabel="Продолжительность работы лампы (часы)", ylabel="Вероятность")
В этом смысле ИФР и функция выживания эквивалентны: если дана любая из них, можно вычислить другую. Но в рамках анализа выживаемости чаще всего приходится работать с кривыми выживания. А вычисление кривой выживания — это шаг к следующему важному понятию, функции риска (hazard function).
Функция риска
В датасете о лампах накаливания каждое значение h представляет собой 12-часовой интервал, заканчивающийся в час h. Его я буду называть «интервал h». Предположим, мы знаем, что лампочка проработала вплоть до интервала h, и хотели бы узнать вероятность того, что она перегорит в течение интервала h. Для ответа на этот вопрос можно использовать функцию выживания, показывающую долю лампочек, «переживших» интервал h, и функцию вероятности (ФВ), которая дает долю лампочек, перегоревших в течение интервала h. Сумма этих показателей равна доле ламп, которые могут перегореть в течение интервала h и поэтому называются находящимися «в группе риска» (at risk). В данном примере 26 % лампочек находились в группе риска в течение интервала 1656:
at_risk = pmf_bulblife + surv_bulblife
at_risk[1656]
0.25999999999999995
И 4 % всех лампочек перегорело в течение интервала 1656:
pmf_bulblife[1656]
0.04
Риск (hazard), или интенсивность отказов, в данном случае — это отношение pmf_bulblife и at_risk:
hazard = pmf_bulblife / at_risk
hazard[1656]
0.15384615384615388
Из всех ламп, которые остались в строю до интервала 1656, около 15 % перегорели в течение этого интервала.
Функцию риска можно вычислить самостоятельно или же воспользоваться библиотекой empiricaldist. В ней представлен объект Hazard, реализующий функцию риска, и метод make_hazard, который ее вычисляет:
hazard_bulblife = surv_bulblife.make_hazard()
hazard_bulblife[1656]
0.15384615384615397
Вот как выглядит функция риска для ламп накаливания:
hazard_bulblife.plot()
decorate(xlabel="Продолжительность работы лампы (часы)", ylabel="Риск")
Видны всплески уровня риска в некоторых местах графика, но такой способ представления функции риска может вводить в заблуждение, особенно в тех частях диапазона, где мало данных. Лучше всего построить график функции накопленного риска (cumulative hazard function), которая представляет собой накопленную сумму (cumulative sum) риска:
cumulative_hazard = hazard_bulblife.cumsum()
cumulative_hazard.plot()
decorate(xlabel="Продолжительность работы лампы (часы)", ylabel="Накопленный
риск")
Там, где вероятность выхода из строя высока, функция накопленного риска показывает крутой подъем. Там же, где эта вероятность мала, кумулятивная функция риска растет полого. В нашем примере можно заметить, что наивысший уровень риска наблюдается между 1500 и 2000 часами. После этого риск снижается — хотя это только благодаря одной необычно долгоживущей лампочке, поэтому в другом датасете он может выглядеть иначе.
Теперь, когда у вас появилось общее представление о функциях выживания и риска, применим их к более объемному датасету.
Данные о браке
Во многих странах люди сегодня вступают в брак позже, чем это было в прежние годы, и все больше людей остаются одинокими. Чтобы изучить эти тенденции на примере США, воспользуемся инструментами анализа выживаемости и данными Национального исследования роста семьи (NSFG) США.
Использованный в предыдущих главах датасет NSFG был посвящен беременностям и содержал по одной строке для каждой упомянутой респонденткой беременности. В этой главе мы будем работать с файлом, который содержит информацию о самих респондентках.
Я собрал ответы, полученные в ходе девяти этапов опроса, проводившихся в период с 1982 по 2019 год, и отобрал данные, касающиеся брака. Инструкции по загрузке этой выборки приведены в блокноте этой главы.
Считаем данные:
resp = pd.read_csv("marriage_nsfg_female.csv.gz")
resp.shape
(70183, 34)
В выборке каждая строка соответствует одной из более чем 70 000 респонденток и содержит следующие переменные, связанные с возрастом и браком:
cmbirth
Известная для всех респонденток дата рождения.
cmintvw
Известная для всех респонденток дата проведения опроса.
cmmarrhx
Дата вступления респондентки в первый брак, если брак имел место и дата известна.
evrmarry
1 — если респондентка состояла в браке до даты проведения опроса, 0 — в противном случае.
Первые три переменные записаны в кодировке «век — месяц» (Century-Month Code, CMC), то есть как целое число месяцев, прошедших с декабря 1899 года. Так что 1 в этом формате записи — это январь 1900 года.
Чтобы изучить эффект смены поколений, сгруппируем респонденток по десятилетию их даты рождения. В этом нам поможет следующая функция, которая принимает значение переменной cmbirth и вычисляет десятилетие даты рождения. В ней использован оператор целочисленного деления // для деления на 10 и округления в меньшую сторону:
month0 = pd.to_datetime("1899-12-31")
def decade_of_birth(cmbirth):
date = month0 + pd.DateOffset(months=cmbirth)
return date.year // 10 * 10
При помощи метода apply применим эту функцию и вычислим десятилетие даты рождения всех респонденток, а затем поместим полученные значения в новый столбец cohort. В данном контексте когорта (cohort) — это группа людей, объединенных чем-то общим (например, десятилетием своего рождения), для целей анализа рассматриваемая как единое целое.
Результатом функции value_counts является количество женщин в каждой когорте:
from thinkstats import value_counts
resp["cohort"] = resp["cmbirth"].apply(decade_of_birth)
value_counts(resp["cohort"])
cohort
1930 325
1940 3608
1950 10631
1960 14953
1970 16438
1980 14271
1990 8552
2000 1405
Name: count, dtype: int64
В этом датасете содержится информация о более 10 000 женщин, родившихся с 1950-х по 1980-е годы, и меньшее количество женщин из десятилетий до и после этого промежутка времени.
Далее рассчитаем возраст каждой респондентки на момент вступления в брак (если он имел место) и на дату проведения опроса:
resp["agemarr"] = (resp["cmmarrhx"] - resp["cmbirth"]) / 12
resp["age"] = (resp["cmintvw"] - resp["cmbirth"]) / 12
Прежде чем приступить к работе с этими данными, определим следующую функцию, которая принимает в качестве аргументов объект DataFrame и список когорт, а возвращает словарь, сопоставляющий каждую когорту с объектом Surv. Для каждой когорты функция выбирает возраст респондентки на момент вступления в первый брак и использует Surv.from_seq для вычисления функции выживания. Аргумент dropna=False позволяет учесть значения NaN в функции выживания, поэтому результат включает женщин, которые не выходили замуж:
from empiricaldist import Surv
def make_survival_map(resp, cohorts):
surv_map = {}
grouped = resp.groupby("cohort")
for cohort in cohorts:
group = grouped.get_group(cohort)
surv_map[cohort] = Surv.from_seq(group["agemarr"], dropna=False)
return surv_map
Вот пример ее использования:
cohorts = [1980, 1960, 1940]
surv_map = make_survival_map(resp, cohorts)
А вот результаты для женщин, родившихся в 1940-х, 1960-х и 1980-х годах:
for cohort, surv in surv_map.items():
surv.plot(label=f"{cohort}-е")
ylim = [-0.05, 1.05]
decorate(xlabel="Возраст", ylabel="P(не состояла в браке)", ylim=ylim)
Если принять эти результаты за истину, то они показывают, что женщины предыдущих поколений вступали в брак в более молодом возрасте и доля хотя бы раз выходивших замуж была больше. Однако мы не будем заниматься их интерпретацией, поскольку они некорректны по двум причинам:
• Как указано в отчете «Национального исследования роста семьи» на с. 3, NSFG использует стратифицированную выборку, это означает, что в нем намеренно избыточно представлены некоторые социальные группы.
• Кроме того, наш способ вычисления функции выживания не учитывает должным образом людей, которые пока не состоят в браке.
Для решения первой проблемы воспользуемся методом ресемплинга, называемым взвешенным бутстрепингом (weighted bootstrap resampling). Для решения второй применим метод, называемый оценкой Каплана — Мейера. Начнем с ресемплинга.
Взвешенный бутстрепинг
В переменной finalwgt датасета NSFG содержится выборочный вес (sampling weight) каждого респондента, то есть количество людей из генеральной совокупности, которое он представляет. Эти веса можно использовать в процессе ресемплинга, чтобы сделать поправку на стратифицированность выборки. Функция ниже принимает на вход объект DataFrame и имя столбца, содержащего выборочные веса. Она выполняет ресемплинг строк датафрейма с учетом выборочных весов и возвращает новый объект DataFrame:
def resample_rows_weighted(df, column="finalwgt"):
n = len(df)
weights = df[column]
return df.sample(n, weights=weights, replace=True)
Текущий датасет включает респондентов из нескольких итераций опроса, называемых циклами. Поэтому для ресемплинга нужно сначала сгруппировать респондентов по циклам, провести ресемплинг каждой группы, а затем снова объединить все группы. Эти шаги выполняет следующая функция:
def resample_cycles(resp):
grouped = resp.groupby("cycle")
samples = [resample_rows_weighted(group) for _, group in grouped]
return pd.concat(samples)
Сначала проведем ресемплинг данных один раз:
sample = resample_cycles(resp)
Позже мы проведем множественный ресемплинг, чтобы посмотреть, насколько велика изменчивость, обусловленная случайной выборкой.
На рисунке ниже показаны результаты после проведения ресемплинга в сравнении с результатами из предыдущего раздела, где он не проводился, обозначенными пунктиром:
for label, surv in surv_map.items():
surv.plot(ls=":", color="gray", alpha=0.6)
survs_resampled = make_survival_map(sample, cohorts)
for label, surv in survs_resampled.items():
surv.plot(label=label)
decorate(xlabel="Возраст", ylabel="P(не состояла в браке)", ylim=ylim)
Разница результатов с проведением ресемплинга и без него существенна, и это доказывает, что для получения верных результатов нужно вводить поправку на стратифицированность выборки.
Теперь перейдем ко второй проблеме, связанной с неполнотой данных.
Оценка функции риска
В примере с лампами накаливания мы знаем продолжительность работы всех 50 ламп, поэтому можем напрямую вычислить функцию выживания, которую затем использовать для вычисления функции риска.
В примере со статистикой браков мы знаем возраст вступления в первый брак для тех респонденток, которые в него вступали к моменту опроса. Но для респонденток, которые никогда не выходили замуж, мы не знаем, в каком возрасте они сделают это в будущем и вступят ли в брак вообще.
Такого рода недостающие данные называются цензурированными (censored). Этот термин может показаться странным, поскольку информация, подвергшаяся цензуре, обычно скрывается преднамеренно, но в данном случае она отсутствует только потому, что мы не знаем будущего.
Однако имеется неполная информация, с которой можно работать: если на момент опроса женщина не состоит в браке, то мы знаем, что возраст ее вступления в брак (если это вообще произойдет) должен превышать ее текущий возраст.
Эту неполную информацию можно использовать для оценки функции риска, а затем с помощью функции риска вычислить функцию выживания. Этот процесс называется оценкой Каплана — Мейера (Kaplan-Meier estimation).
В демонстративных целях я выберу только одну когорту из повторной выборки данных:
resp60 = sample.query("cohort == 1960")
Для респонденток, которые состояли в браке на момент проведения опроса, выберем их возраст на момент вступления в первый брак. Всего имеется 9921 такой случай, назовем их «завершенными» (complete):
complete = resp60.query("evrmarry == 1")["agemarr"]
complete.count()
9921
Для респонденток, которые не состояли в браке, выберем их возраст на момент опроса. Всего есть 5468 таких случаев; их мы назовем «незавершенными» (ongoing):
ongoing = resp60.query("evrmarry == 0")["age"]
ongoing.count()
5468
Теперь, чтобы оценить функцию риска, вычислим общее количество случаев, которые были «в группе риска» для всех значений возраста, включая всех тех, кто не состоял в браке до этого возраста. Будет удобно создать объект FreqTab, который подсчитывает количество завершенных и незавершенных случаев для всех возрастов:
from empiricaldist import FreqTab
ft_complete = FreqTab.from_seq(complete)
ft_ongoing = FreqTab.from_seq(ongoing)
Например, 58 респонденток сообщили, что они впервые вышли замуж в возрасте 25 лет:
ft_complete[25]
58
А 5 респонденток, которые были опрошены в возрасте 25 лет, сообщили, что никогда не состояли в браке:
ft_ongoing[25]
5
Из этих объектов класса FreqTab можно вычислить ненормированные объекты класса Surv, которые содержат количество завершенных и незавершенных случаев для возраста, превышающего заданный:
surv_complete = ft_complete.make_surv()
surv_ongoing = ft_ongoing.make_surv()
Например, 2848 женщин сообщили, что вступили в брак после 25 лет:
surv_complete[25]
2848
И 2273 опрошенных в возрасте после 25 лет сообщили, что никогда не состояли в браке:
surv_ongoing[25]
2273
Сумма только что вычисленных четырех чисел представляет собой число респондентов, которые находились в группе риска, то есть могли бы вступить в брак в возрасте 25 лет. Термин «группа риска» является наследием анализа выживаемости в медицине, где он зачастую связан с угрозой заболевания или смерти. Он может показаться неуместным применительно к браку, обычно считающемуся позитивным этапом жизни. Как бы то ни было, вот наши выкладки:
at_risk = ft_complete[25] + ft_ongoing[25] + surv_complete[25] +
surv_ongoing[25]
at_risk
5184
Из этого числа количество тех, кто действительно вступил в брак в возрасте 25 лет, равно ft_complete[25]. И теперь мы можем вычислить функцию риска в возрасте 25 лет:
hazard = ft_complete[25] / at_risk
hazard
0.011188271604938271
Мы показали, как вычислить функцию риска для определенного возраста. Теперь вычислим ее целиком, для всех возрастов. Воспользуемся методом union класса Index библиотеки Pandas, чтобы получить индекс, содержащий все значения возраста из объектов ft_complete и ft_ongoing по порядку:
ts = pd.Index.union(ft_complete.index, ft_ongoing.index)
Теперь посчитаем количество женщин группы риска для всех возрастов, сделав выборку значений в каждом из объектов FreqTab и Surv по возрастам из переменной ts:
at_risk = ft_complete(ts) + ft_ongoing(ts) + surv_complete(ts) +
surv_ongoing(ts)
Наконец, вычислим функцию риска для всех значений возраста и поместим результаты в объект Hazard:
from empiricaldist import Hazard
hs = ft_complete(ts) / at_risk
hazard = Hazard(hs, ts)
Вот как выглядит функция накопленного риска:
hazard.cumsum().plot()
decorate(xlabel="Возраст", ylabel="Накопленный риск")
В диапазоне от 20 до 30 лет наблюдается наибольшая крутизна графика, это означает, что в этом возрасте одинокая женщина подвергается наибольшему «риску» вступления в брак. После этого накопленный риск выходит на плато, что означает постепенное снижение риска.
Оценка функции выживания
Зная функцию выживания, можно вычислить функцию риска. А если наоборот?
Вот один из способов. Функция риска показывает вероятность вступить в брак в зависимости от возраста, если брака раньше не было. Следовательно, комплементарной к функции риска является вероятность остаться не замужем как функция возраста.
Чтобы «пережить» определенный возраст t, нужно оставаться не замужем в любом возрасте до t включительно. И вероятность этого события является накопленным произведением (cumulative product) функции, комплементарной к функции риска, которую можно вычислить следующим образом:
ps = (1 - hazard).cumprod()
У объекта Hazard есть метод make_surv, который выполняет этот подсчет:
surv = hazard.make_surv()
Вот как выглядит полученный результат по сравнению с предыдущим (пунктирная линия), в котором делалась поправка на стратифицированность выборки с помощью ресемплинга, но не учитывались цензурированные данные:
survs_resampled[1960].plot(ls=":", color="gray", label="повторная выборка")
surv.plot(label="оценка Каплана — Мейера")
decorate(xlabel="Возраст", ylabel="P(не состояла в браке)", ylim=ylim)
Это показывает, насколько важно правильно обрабатывать цензурированные данные.
Подобная функция выживания легла в основу известной статьи 1986 года. Журнал Newsweek сообщал, что у 40-летней незамужней женщины «больше шансов быть убитой террористом», чем выйти замуж. Это утверждение получило широкую огласку и стало частью массовой культуры, но оно изначально было неверным (потому что основывалось на ошибочном анализе) и оказалось еще более неверным (из-за культурных изменений, которые уже начались тогда). В 2006 году Newsweek опубликовал другую статью, в которой признал свою ошибку.
Я рекомендую узнать больше об этой статье, статистических данных, на которых она была основана, и реакции на нее. Она напоминает о наших моральных обязательствах тщательно проводить статистический анализ, интерпретировать результаты с должным уровнем скептицизма, а также аккуратно и честно представлять их на суд общественности.
Следующая функция реализует все эти этапы оценки Каплана — Мейера. Она принимает на вход последовательности значений времени выживания для завершенных и незавершенных случаев, а затем возвращает объект Hazard:
def estimate_hazard(complete, ongoing):
"""Оценка Каплана — Мейера."""
ft_complete = FreqTab.from_seq(complete)
ft_ongoing = FreqTab.from_seq(ongoing)
surv_complete = ft_complete.make_surv()
surv_ongoing = ft_ongoing.make_surv()
ts = pd.Index.union(ft_complete.index, ft_ongoing.index)
at_risk = (
ft_complete(ts) + ft_ongoing(ts) +
surv_complete(ts) + surv_ongoing(ts)
)
hs = ft_complete(ts) / at_risk
return Hazard(hs, ts)
А функция ниже принимает на вход данные группы респондентов, извлекает значения времени выживания, получает функцию риска в результате вызова estimate_hazard, а затем вычисляет соответствующую функцию выживания:
def estimate_survival(group):
"""Оценка функции выживания."""
complete = group.query("evrmarry == 1")["agemarr"]
ongoing = group.query("evrmarry == 0")["age"]
hf = estimate_hazard(complete, ongoing)
sf = hf.make_surv()
return sf
Нам вскоре потребуются эти функции для вычисления доверительного интервала для функции выживания. Но сначала рассмотрим еще один способ вычисления оценок Каплана — Мейера.
Библиотека Lifelines
Python-пакет lifelines предоставляет инструменты для анализа выживаемости, в том числе функции для вычисления оценок Каплана — Мейера.
Воспроизведем с его помощью результат, полученный в предыдущем разделе, тем самым доказав его правильность. Сначала вычислим функцию выживания, вызвав estimate_survival:
surv = estimate_survival(resp60)
Далее получим ее с помощью инструментов lifelines. Сначала преобразуем данные в формат, совместимый с lifelines:
complete = complete.dropna()
durations = np.concatenate([complete, ongoing])
event_observed = np.concatenate([np.ones(len(complete)),
np.zeros(len(ongoing))])
Теперь можно создать объект KaplanMeierFitter и подогнать модель под данные:
from lifelines import KaplanMeierFitter
kmf = KaplanMeierFitter()
kmf.fit(durations=durations, event_observed=event_observed)
<lifelines.KaplanMeierFitter:"KM_estimate", fitted with 15389 total
observations,
5468 right-censored observations>
После подгонки к данным можно отобразить результаты с помощью функции plot. Результаты содержат оценку функции выживания и доверительный интервал, хотя в данном случае последний неверен, поскольку не учитывает стратифицированность выборки:
kmf.plot()
decorate(xlabel="Возраст", ylabel="P(не состояла в браке)", ylim=ylim)
В отличие от функции выживания, которую мы вычислили, функция из библиотеки lifelines определена, начиная с 0. Но остальная часть значений функции совпадает с нашими — в пределах погрешности арифметики с плавающей точкой:
ps = kmf.survival_function_["KM_estimate"].drop(0)
np.allclose(ps, surv)
True
В следующем разделе воспользуемся взвешенным ресемплингом для вычисления доверительного интервала, который учитывает стратифицированность выборки.
Доверительный интервал
Вычисленная нами оценка Каплана — Мейера основана на однократном ресемплинге датасета. Чтобы получить представление о том, насколько велика вариативность из-за случайной выборки, проведем расчеты с несколькими повторными выборками и построим график результатов.
Применим следующую функцию, которая принимает объект DataFrame и список когорт, оценивает функцию выживания для каждой когорты и возвращает словарь, сопоставляющий каждую целочисленную когорту с объектом класса Surv.
Эта функция идентична функции make_survival_map, за исключением того что она вызывает функцию estimate _survival, использующую оценку Каплана — Мейера, вместо метода Surv.from_seq, который подходит только при отсутствии цензурированных данных:
def estimate_survival_map(resp, cohorts):
"""Создать словарь, сопоставляющий когортам объекты Surv."""
surv_map = {}
grouped = resp.groupby("cohort")
for cohort in cohorts:
group = grouped.get_group(cohort)
surv_map[cohort] = estimate_survival(group)
return surv_map
В цикле ниже случайно генерируется 101 повторная выборка нашего датасета и создается список из 101 словаря, содержащего оценки функции выживания:
cohorts = [1940, 1950, 1960, 1970, 1980, 1990]
surv_maps = [estimate_survival_map(resample_cycles(resp), cohorts)
for i in range(101)]
Чтобы отобразить результаты, воспользуемся функцией ниже, которая принимает полученный список словарей, целочисленное значение когорты и строку с названием цвета графика (color). Функция перебирает словари, выбирает функцию выживаемости для указанной когорты и отображает ее полупрозрачной линией — это один из способов визуализировать изменчивость в повторных выборках:
def plot_cohort(surv_maps, cohort, color):
"""Построить диаграмму результатов для одной когорты."""
survs = [surv_map[cohort] for surv_map in surv_maps]
for surv in survs:
surv.plot(color=color, alpha=0.05)
x, y = surv.index[-1], surv.iloc[-1]
plt.text(x + 1, y, f"{cohort}-е", ha="left", va="center")
Вот результаты по когортам рождения с 1940-х по 1990-е годы:
colors = [f"C{i}" for i in range(len(cohorts))]
for cohort, color in zip(cohorts, colors):
plot_cohort(surv_maps, cohort, color)
xlim = [8, 55]
decorate(xlabel="Возраст", ylabel="P(не состояла в браке)",
xlim=xlim, ylim=ylim)
Полученная диаграмма вполне подходит для изучения, но линии выглядят размытыми, а некоторые надписи перекрывают друг друга. Для получения рисунка, подходящего для публикации, потребуется больше усилий, но мы пока не будем ничего усложнять.
На диаграмме видны несколько закономерностей:
• Женщины, родившиеся в 1940-х годах, выходили замуж раньше всех, а женщины, родившиеся в 1950-х и 1960-х годах, хотя и выходили замуж позже, но среди них доля оставшихся незамужними примерно одинакова.
• Женщины, родившиеся в 1970-х годах, позже выходили замуж и чаще оставались незамужними по сравнению с женщинами в предыдущих когортах.
• Родившиеся в 1980-х и 1990-х годах вступают в брак еще позже и еще чаще остаются незамужними, хотя в будущем эти тенденции могут измениться.
Придется дождаться следующего выпуска данных NSFG, чтобы продолжить наши исследования.
Ожидаемое остаточное время жизни
Имея некоторое распределение, можно рассчитать ожидаемое остаточное время жизни (expected remaining lifetime) как функцию от прошедшего времени. Например, если дано распределение продолжительности беременности, можно рассчитать ожидаемое время до наступления родов. В демонстрационных целях будем использовать данные о беременностях из NSFG.
С помощью функции get_nsfg_groups считаем данные и разделим их на соответствующие первым (firsts) и последующим (others) детям:
from nsfg import get_nsfg_groups
live, firsts, others = get_nsfg_groups()
Начнем с однократного ресемплинга данных:
sample = resample_rows_weighted(live, "finalwgt")
Вот функция вероятности продолжительности беременности:
pmf_durations = Pmf.from_seq(sample["prglngth"])
Теперь предположим, что сейчас начало 36-й недели беременности. Зная, что чаще всего срок беременности составляет 39 недель, естественно ожидать, что остаточное время составит 3–4 недели. Чтобы уточнить эту оценку, отберем значения в распределении, которые больше или равны 36 неделям:
t = 36
is_remaining = pmf_durations.qs >= t
Далее создадим новый объект Pmf, содержащий только эти значения со сдвигом влево — чтобы текущее время было равно 0:
ps = pmf_durations.ps[is_remaining]
qs = pmf_durations.qs[is_remaining] — t
pmf_remaining = Pmf(ps, qs)
Поскольку было выбрано подмножество значений исходного объекта Pmf, вероятности больше не суммируются в 1, но это можно исправить при помощи нормировки Pmf:
pmf_remaining.normalize()
0.9155006558810669
Вот диаграмма полученного распределения остаточного времени на начало 36-й недели:
pmf_remaining.bar(label="36-я неделя")
decorate(xlabel="Оставшийся срок (недели)", ylabel="Вероятность")
Среднее значение этого распределения — ожидаемое остаточное время:
pmf_remaining.mean()
3.2145671641791043
Следующая функция объединяет эти шаги и вычисляет распределение остаточного времени для заданного объекта Pmf в данный момент времени t:
def compute_pmf_remaining(pmf, t):
"""Распределение остаточного времени."""
is_remaining = pmf.qs >= t
ps = pmf.ps[is_remaining]
qs = pmf.qs[is_remaining] — t
pmf_remaining = Pmf(ps, qs)
pmf_remaining.normalize()
return pmf_remaining
Следующая функция принимает на вход объект Pmf с продолжительностями беременности и вычисляет ожидаемое остаточное время в начале каждой недели с 36-й по 43-ю:
def expected_remaining(pmf):
index = range(36, 44)
expected = pd.Series(index=index)
for t in index:
pmf_remaining = compute_pmf_remaining(pmf, t)
expected[t] = pmf_remaining.mean()
return expected
Вот результаты однократного ресемплинга данных:
expected = expected_remaining(pmf_durations)
expected
36 3.214567
37 2.337714
38 1.479095
39 0.610133
40 0.912517
41 0.784211
42 0.582301
43 0.589372
dtype: float64
Чтобы оценить, насколько велика вариативность, обусловленная случайной выборкой, проведем эксперимент с несколькими ресемплингами и построим график результатов:
for i in range(21):
sample = resample_rows_weighted(live, "finalwgt")
pmf_durations = Pmf.from_seq(sample["prglngth"])
expected = expected_remaining(pmf_durations)
expected.plot(color="C0", alpha=0.1)
decorate(
xlabel="Неделя беременности",
ylabel="Ожидаемое остаточное время (недели)", ylim=[0, 3.4]
)
С 36-й по 39-ю неделю ожидаемое остаточное время уменьшается и в начале 39-й недели достигает уровня около 0.6 недели. Но далее кривая стабилизируется. В начале 40-й недели ожидаемое остаточное время все еще близко к 0.6 недели, как и в начале 41-й, 42-й и 43-й недель. Для родителей, с нетерпением ожидающих появления на свет ребенка, такая закономерность может показаться довольно мучительной.
Глоссарий
Анализ выживаемости (survival analysis)
Набор методов для описания и предсказания времени до наступления определенного события, часто касающийся анализа продолжительности жизни или срока службы.
Функция выживания (survival function)
Функция, которая сопоставляет времени t вероятность «пережить» момент t.
Функция риска (hazard function)
Функция, которая сопоставляет времени t долю случаев, в которых событие наступило в момент t, среди всех случаев, «доживших» до момента t.
Функция накопленного риска (cumulative hazard function)
Накопленная сумма функции риска, часто используемая для визуализации.
Взвешенный бутстрепинг (weighted bootstrap resampling)
Форма ресемплинга, которая использует выборочные веса для корректировки стратифицированной выборки с целью моделирования репрезентативной выборки.
Цензурированные данные (censored data)
Данные, которые известны лишь частично, поскольку определенное событие еще не произошло или осталось ненаблюдаемым.
Оценка Каплана — Мейера (Kaplan-Meier estimation)
Метод оценки функции выживания и функции риска для датасетов с цензурированными наблюдениями.
Когорта (cohort)
Группа испытуемых с общими характеристиками, например диагноз или десятилетие даты рождения.
Упражнения
Упражнение 13.1
Методы, описанные в этой главе, можно использовать для оценки функции риска и функции выживания применительно к продолжительности брака. Для простоты мы рассмотрим только первые браки и ограничимся разводом в качестве конечной точки, не рассматривая такие сценарии, как раздельное проживание или смерть.
В данных NSFG переменная cmdivorcx содержит дату развода для первого брака каждого респондента, если таковой имел место, записанную в кодировке «век — месяц». Подсчитайте продолжительность браков, завершившихся разводом, и продолжительность текущих браков.
• Для завершенных случаев вычислите время, прошедшее между cmdivorcx и cmmarrhx. Если оба значения корректны, то есть не равны NaN, это означает, что первый брак респондента закончился разводом.
• Чтобы выявить незавершенные случаи, выберите женщин, которые были замужем только один раз и все еще состоят в браке. Используйте переменную fmarno, в которой записано количество браков каждой респондентки, и переменную fmarital, кодирующую ее семейное положение, — значение 1 указывает на то, что респондентка замужем.
Иногда значения этих переменных приблизительные, поэтому вы можете встретить небольшое количество случаев с отрицательной разницей, но она не должна превышать одного года.
Проведите оценку функции риска и функции выживания для продолжительности брака. Постройте график функции накопленного риска — когда риск развода наиболее высок? Постройте график функции выживания — какая доля браков заканчивается разводом?
Упражнение 13.2
В 2012 году группа демографов из Университета Южной Калифорнии подсчитала ожидаемую продолжительность жизни людей, родившихся в Швеции в начале XIX и XX веков. Для людей в возрасте от 0 до 91 года они рассчитали повозрастной коэффициент смертности (age-specific mortality rate), представляющий собой долю людей, умирающих в определенном возрасте, из всех, кто доживает до этого возраста, — зависимость, в которой вы могли узнать функцию риска.
Я оцифровал онлайн-данные из статьи и сохранил их в CSV-файле. Инструкции по загрузке данных приведены в Jupyter-блокноте этой главы.
Данные можно загрузить так:
mortality = pd.read_csv("mortality_rates_beltran2012.csv",
header=[0, 1]).dropna()
Следующая функция с помощью интерполяции данных вычисляет функцию риска с приближенными значениями коэффициента смертности для всех возрастов от 0 до 99 лет:
from scipy.interpolate import interp1d
from empiricaldist import Hazard
def make_hazard(ages, rates):
interp = interp1d(ages, rates, fill_value="extrapolate")
xs = np.arange(0, 100)
ys = np.exp(interp(xs))
return Hazard(ys, xs)
Теперь можно создать объект Hazard:
ages = mortality["1800", "X"].values
rates = mortality["1800", "Y"].values
hazard = make_hazard(ages, rates)
Вот график коэффициента смертности:
hazard.plot()
decorate(xlabel="Возраст (годы), ylabel="Риск")
На основе этих данных с помощью метода make_surv создайте функцию выживания и вызовите make_cdf, чтобы вычислить соответствующую ИФР. Постройте график результатов.
Затем с помощью метода make_pmf создайте объект Pmf, представляющий распределение продолжительности жизни, и изобразите его графически. Наконец, используйте функцию compute_pmf_remaining для вычисления среднего остаточного времени жизни для всех возрастов от 0 до 99 лет. Постройте диаграмму полученных результатов.
На кривой остаточного времени жизни вы должны заметить парадоксальную закономерность: в течение первых нескольких лет жизни остаточное время жизни увеличивается. Так как в начале XIX века детская смертность была очень высокой, можно было ожидать, что дети более старшего возраста проживут дольше, чем дети младшего возраста. Примерно после 5 лет ожидаемая продолжительность жизни возвращается к предсказуемому поведению — ожидается, что молодые люди проживут дольше, чем пожилые.
Если вас заинтересовала эта тема, рекомендую ознакомиться с главой 5 моей книги Probably Overthinking It (University of Chicago Press, 2023). В ней вы найдете другие столь же непредсказуемые результаты из разных областей статистики.
fitted with 15389 total observations — подогнано к 15 389 наблюдениям. — Примеч. пер.
5468 right-censored observations — 5468 цензурированных справа наблюдений. — Примеч. пер.
