Чем грозит Москве британский штамм COVID-19? Отвечаем с помощью Python и дифуров

Всем привет! Меня зовут Борис, я выпускник программы Науки о данных ФКН ВШЭ, работаю ML Инженером и преподаю в OTUS на курсах ML Professional, DL Basic, Computer Vision.

В первых числах января 2021 я узнал про британский штамм коронавируса,прогнозы о новой волне в США. Я подумал: аналитик данных я или кто? Мне захотелось забить гвоздик своим микроскопом и узнать, вызовет ли британский штамм волну заражений в Москвеи стоит ли покупать авиабилеты на лето.

Выглядело как приключение на две недели, но превратилось в исследование на три месяца. В процессе я выяснил, что хороших материалов по созданию эпидемиологических моделей практически нет. Банально авторы статей по моделированию COVID-19 в топовых журналах даже не делают train-test split.

Я предлагаю туториал на основе своего исследования. В нём я постарался передать все важные детали, которые сэкономили бы мне много недель, если бы о них кто-то писал.

Ноутбук к туториалу.

Краткая презентация с результатами исследования.

Препринт статьи: TBA

Какой британский штамм?

COVID-19 в представлении не нуждается, к нему мы уже привыкли. Однако новые мутации (штаммы) вируса способны изменить ход эпидемии. Сейчас ВОЗвыделяеттри таких штамма (variants of concern). Среди них штамм B.1.1.7, так называемый британский штамм.

Почему меня заинтересовал именно он? Во-первых,исследования показывают, что этот штамм распространяется на 40% - 90% быстрее, чем обычный коронавирус (для продвинутых: R0 на 40% - 90% больше). Во-вторых, случай заражения этим штаммом уже былобнаруженв России. В-третьих, непохоже, чтобы кто-то в России готовился к его появлению.

По результатам исследования пришел к выводу, что B.1.1.7 действительно способен привести к новой волне заражений COVID-19 и может унести 30 тысяч жизней.

Исходные данные

Нам понадобится статистика по заражениям, смертям и выздоровлениям в Москве. Её можно скачать со страницы.
На момент исследования был доступен промежуток дат: 2020.03.10 - 2021.03.23.

Изначально данные даны для всей России. Модель, которую мы будем использовать, работает лучше для замкнутых популяций, так что ограничимся моделированием для Москвы. Кроме того поправим индексы и переименуем колонки.

df.columns = ['date', 'region',             'total_infected', 'total_recovered', 'total_dead',            'deaths_per_day', 'infected_per_day', 'recovered_per_day']df = df[df.region == 'Москва'].reset_index()df['date'] = pd.to_datetime(df['date'], format='%d.%m.%Y')df['infected'] = df['total_infected'] - df['total_recovered'] - df['total_dead']df = df.drop(columns=['index', 'region'])df = df.sort_values(by='date')df.index = df.datedf['date'] = df.date.asfreq('D')

Добавим сглаженные скользящим окном в 7 дней версии всех колонок. Они пригодятся нам позже.

df_smoothed = df.rolling(7).mean().round(5)df_smoothed.columns = [col + '_ma7' for col in df_smoothed.columns]full_df = pd.concat([df, df_smoothed], axis=1)for column in full_df.columns:    if column.endswith('_ma7'):        original_column = column.strip('_ma7')        full_df[column] = full_df[column].fillna(full_df[original_column])

В итоге работаем со следующими колонками, а так же их сглаженными версиями:

• infected_per_day новых заражений в данный день.

• recovered_per_day новых выздоровлений в данный день.

• deaths_per_day погибших от инфекции в данный день.

• total_infected всего заражений к этому моменту, кумулятивная суммаinfected_per_day.

• total_dead всего погибших к этому моменту, кумулятивная сумма отdeaths_per_day.

• total_recovered всего выздоровлений к этому моменту, кумулятивная сумма отrecovered_per_day.

Основа модели: SEIRD

SIRэто очень простая модель, которая симулирует развитие эпидемии во времени. Популяция разделяется на три группы: Susceptible, Infected, Recovered.
Идея на пальцах:

• Есть замкнутая популяция и изначальное количество зараженных (Infected).

• В течение дня каждый зараженный с некоторой вероятностью инфицирует кого-то из группы Susceptible. Новые зараженные переходят в группу Infected.

• Находящиеся в группе Infected выздоравливают, когда проходит срок длительности болезни и переходят в группу Recovered.

• Запускаем этот процесс на много дней вперед и смотрим, как меняется численность групп.

SIR не учитывает наличие инкубационного периода и погибших, так что сразу перейдем к расширению SEIRD и рассмотрим его подробнее. SEIRD работает точно так же, но разделяет популяцию на следующие группы:

• Susceptible могут быть заражены.

• Exposed заражены вирусом, но не распространяют его, то есть проходят инкубационный период.

• Infectious распространяют вирус.

• Recovered переболели и не могут быть заражены.

• Dead погибли.

Изменение численности групп за день определяют дифференциальные уравнения:

Например, третье уравнение описывает изменение группы Infected за день. С вероятностьюу человека из группыEзаканчивается инкубационный период и он переходит в группуI. Так же с вероятностьючеловек из группыIвыздоравливает или погибает, переходя в группуRилиD. Куда он попадет определяется параметром, который соответствует доле смертности (Infection Fatality Rate). Таким образом, каждый день группу пополняют новые зараженные, а переболевшие уходят из группы.

Параметры модели:
case fatality rate.
количество человек, которых один переносчик заражает за день.
1 делить на среднюю длину инкубационного периода.
y 1 делить на среднюю длительность болезни.
R0 = /y базовое репродуктивное число, то есть количество человек, которых один переносчик заражает за время своей болезни.

Есть три причины использовать именно SEIRD. Во-первых, параметры SEIRD это просто характеристики болезни, а значит мы сможем прикинуть их значения на основании медицинских исследований COVID-19 как болезни. Во-вторых, эту модель очень легко модифицировать, что мы и будем делать дальше, добавляя меры карантина, неполноту статистики и второй штамм. В-третьих, эта модель позволяет легко покрутить разные сценарии. Например, можно смоделировать ситуацию: Что будет, если в этой же популяции появится второй штамм и параметрR0у него на 90% больше? Сделать подобное с помощью, например, градиентного бустинга, гораздо сложнее.

Главный минус SEIRD: это очень простая и полностью детерминированная модель. Она чувствительна к изначальным значениям и параметрам. Предсказания такой модели могут на порядки отличаться от реальных значений. Для нас это не страшно: мы ведь не пытаемся сделать самый точный прогноз по количеству зараженных, мы просто пытаемся понятьможет ли штамм B.1.1.7 вызвать новую волну. Интересно, что, несмотря на это, лучшая на текущий момент модель по прогнозированию COVID-19 этоочень навороченная SEIR модель.

Реализуем классический SEIRD

Привожу минимальную реализацию SEIRD.

class BarebonesSEIR:    def init(self, params=None):        self.params = params<span class="hljs-function" style="box-sizing: border-box;"><span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">def</span> <span class="hljs-title" style="box-sizing: border-box; color: rgb(121, 93, 163);">get_fit_params</span>(<span class="hljs-params" style="box-sizing: border-box;">self</span>):</span>    params = lmfit.Parameters()    params.add(<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">"population"</span>, value=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">12_000_000</span>, vary=<span class="hljs-literal" style="box-sizing: border-box; color: rgb(0, 134, 179);">False</span>)    params.add(<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">"epidemic_started_days_ago"</span>, value=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">10</span>, vary=<span class="hljs-literal" style="box-sizing: border-box; color: rgb(0, 134, 179);">False</span>)    params.add(<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">"r0"</span>, value=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">4</span>, <span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">min</span>=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">3</span>, <span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">max</span>=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">5</span>, vary=<span class="hljs-literal" style="box-sizing: border-box; color: rgb(0, 134, 179);">True</span>)    params.add(<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">"alpha"</span>, value=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">0.0064</span>, <span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">min</span>=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">0.005</span>, <span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">max</span>=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">0.0078</span>, vary=<span class="hljs-literal" style="box-sizing: border-box; color: rgb(0, 134, 179);">True</span>)  <span class="hljs-comment" style="box-sizing: border-box; color: rgb(150, 152, 150);"># CFR</span>    params.add(<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">"delta"</span>, value=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>/<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">3</span>, <span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">min</span>=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>/<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">14</span>, <span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">max</span>=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>/<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">2</span>, vary=<span class="hljs-literal" style="box-sizing: border-box; color: rgb(0, 134, 179);">True</span>)  <span class="hljs-comment" style="box-sizing: border-box; color: rgb(150, 152, 150);"># E -&gt; I rate</span>    params.add(<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">"gamma"</span>, value=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>/<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">9</span>, <span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">min</span>=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>/<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">14</span>, <span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">max</span>=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>/<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">7</span>, vary=<span class="hljs-literal" style="box-sizing: border-box; color: rgb(0, 134, 179);">True</span>)  <span class="hljs-comment" style="box-sizing: border-box; color: rgb(150, 152, 150);"># I -&gt; R rate</span>    params.add(<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">"rho"</span>, expr=<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'gamma'</span>, vary=<span class="hljs-literal" style="box-sizing: border-box; color: rgb(0, 134, 179);">False</span>)  <span class="hljs-comment" style="box-sizing: border-box; color: rgb(150, 152, 150);"># I -&gt; D rate</span>    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">return</span> params<span class="hljs-function" style="box-sizing: border-box;"><span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">def</span> <span class="hljs-title" style="box-sizing: border-box; color: rgb(121, 93, 163);">get_initial_conditions</span>(<span class="hljs-params" style="box-sizing: border-box;">self, data</span>):</span>    <span class="hljs-comment" style="box-sizing: border-box; color: rgb(150, 152, 150);"># Simulate such initial params as to obtain as many deaths as in data</span>    population = self.params[<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'population'</span>]    epidemic_started_days_ago = self.params[<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'epidemic_started_days_ago'</span>]    t = np.arange(epidemic_started_days_ago)    (S, E, I, R, D) = self.predict(t, (population - <span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>, <span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">0</span>, <span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>, <span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">0</span>, <span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">0</span>))    I0 = I[-<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>]    E0 = E[-<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>]    Rec0 = R[-<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>]    D0 = D[-<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>]    S0 = S[-<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>]    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">return</span> (S0, E0, I0, Rec0, D0)<span class="hljs-function" style="box-sizing: border-box;"><span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">def</span> <span class="hljs-title" style="box-sizing: border-box; color: rgb(121, 93, 163);">step</span>(<span class="hljs-params" style="box-sizing: border-box;">self, initial_conditions, t</span>):</span>    population = self.params[<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'population'</span>]    delta = self.params[<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'delta'</span>]    gamma = self.params[<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'gamma'</span>]    alpha = self.params[<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'alpha'</span>]    rho = self.params[<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'rho'</span>]        rt = self.params[<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'r0'</span>].value    beta = rt * gamma    S, E, I, R, D = initial_conditions    new_exposed = beta * I * (S / population)    new_infected = delta * E    new_dead = alpha * rho * I    new_recovered = gamma * (<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span> - alpha) * I    dSdt = -new_exposed    dEdt = new_exposed - new_infected    dIdt = new_infected - new_recovered - new_dead    dRdt = new_recovered    dDdt = new_dead    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">assert</span> S + E + I + R + D - population &lt;= <span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1e10</span>    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">assert</span> dSdt + dIdt + dEdt + dRdt + dDdt &lt;= <span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1e10</span>    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">return</span> dSdt, dEdt, dIdt, dRdt, dDdt<span class="hljs-function" style="box-sizing: border-box;"><span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">def</span> <span class="hljs-title" style="box-sizing: border-box; color: rgb(121, 93, 163);">predict</span>(<span class="hljs-params" style="box-sizing: border-box;">self, t_range, initial_conditions</span>):</span>    ret = odeint(self.step, initial_conditions, t_range)    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">return</span> ret.T

Давайте разберемся.

Методget_fit_paramsпросто создает хранилище парамтеров. Мы используем библиотекуlmfit, но пока что не оптимизируем параметры, а просто используем структуруParametersкак продвинутый словарь. Для задания диапазонов параметров я использовалмета-анализы характеристик COVID-19.

Параметрepidemic_started_days_agoпозволяет задать дату появления первого зараженного. В моём предположении это2 марта 2020.

Методstepсодержит дифференциальные уравнения модели. Он получает на вход изначальные численности всех групп в виде кортежаinitial_conditionsи номер дняt, возвращает изменения всех групп за день.

Методpredictполучает на вход начальные условия и массивt_range. Он просто запускаетscipy.integrate.odeintна методеstep, чтобы проинтегрировать дифференциальные уравнения и получить численность групп в каждый день из массиваt_range.

Методget_initial_conditionsзапускает модель с момента появления первого зараженнного наepidemic_started_days_agoдней и возвращает полученные численности групп. Далее эти численности используются как начальные значения для моделирования последующих дней.

Запустим модель и сравним её предсказания с реальностью:

model = BarebonesSEIR()model.params = model.get_fit_params()train_initial_conditions = model.get_initial_conditions(train_subset)train_t = np.arange(len(train_subset))(S, E, I, R, D) = model.predict(train_t, train_initial_conditions)plt.figure(figsize=(10, 7))plt.plot(train_subset.date, train_subset['total_dead'], label='ground truth')plt.plot(train_subset.date, D, label='predicted', color='black', linestyle='dashed' )plt.legend()plt.title('Total deaths')plt.show()

Модель предсказывает, что все уже умерли или переболели, а эпидемия давно закончилась. Результат ужасный, но ожидаемый.

Взглянем на динамику всех групп SEIRD:

Модель описывает одну волну заражений, после которой почти все выздоравливают, достигается групповой иммунитет и эпидемия сходит на нет. Даже если мы подберем самые лучшие параметры, эта модель просто не способна описать несколько волн эпидемии. Нам необходимо изменить систему уравнений, чтобы получить более полезный результат.

Добавляем сдерживающие меры

Одна из причин появления нескольких волн это сдерживающие меры. Предположим, что карантинные меры в деньt, снижаютR0болезни на процентq(t). Тогда для каждого дня можно расчитать репродуктивное число с учетом карантина:Rt = R0 - R0 * q(t). Отсюда мы можем вычислить(t) = Rt * yи использовать её в уравнениях.

Осталось только задать функцииq(t). В своём исследовании я использовал простую ступенчатую функцию. Делим весь промежуток на отрезки по 60 дней, для каждого отрезка добавляем новый параметр: уровень карантина в данный период. Переходы между отрезками стоит сделать плавными с помощью, иначе функция не будет гладкой и подбирать параметры будет сложно.

Реализуем небольшой пример для наглядности.

def sigmoid(x, xmin, xmax, a, b, c, r):    x_scaled = (x - xmin) / (xmax - xmin)    out = (a * np.exp(c * r) + b * np.exp(r * x_scaled)) / (np.exp(c * r) + np.exp(x_scaled * r))    return outdef stepwise_soft(t, coefficients, r=20, c=0.5):    t_arr = np.array(list(coefficients.keys()))min_index = np.<span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">min</span>(t_arr)max_index = np.<span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">max</span>(t_arr)<span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">if</span> t &lt;= min_index:    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">return</span> coefficients[min_index]<span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">elif</span> t &gt;= max_index:    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">return</span> coefficients[max_index]<span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">else</span>:    index = np.<span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">min</span>(t_arr[t_arr &gt;= t])<span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">if</span> <span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">len</span>(t_arr[t_arr &lt; index]) == <span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">0</span>:    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">return</span> coefficients[index]prev_index = np.<span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">max</span>(t_arr[t_arr &lt; index])<span class="hljs-comment" style="box-sizing: border-box; color: rgb(150, 152, 150);"># sigmoid smoothing</span>q0, q1 = coefficients[prev_index], coefficients[index]out = sigmoid(t, prev_index, index, q0, q1, c, r)<span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">return</span> outt_range = np.arange(100)coefficients = {    0: 0,    30: 0.5,    60: 1,    100: 0.4,}plt.title('Пример функции уровня карантина')plt.scatter(coefficients.keys(), coefficients.values(), label='Моменты изменения уровня карантина')plt.plot(t_range, [stepwise_soft(t, coefficients, r=20, c=0.5) for t in t_range], label='Ступенчатая функция с плавными переходами')plt.xlabel('t')plt.ylabel('Уровень канартина')plt.legend(loc='lower center', bbox_to_anchor=(0.5, -0.6),)plt.show()

Чтобы добавить это в нашу SEIRD модель нужно сделать несколько шагов:

1. Добавляем мультипликаторы карантина вget_fit_params. Изначально карантин находится на уровне нуля, но для остальных отрезков его уровень будет определен оптимизацией. Так же появляется новый гиперпараметр, задающий длину отрезков:stepwise_size, по умолчанию 60 дней.

def get_fit_params(self, data):    ...    params.add(f"t0_q", value=0, min=0, max=0.99, brute_step=0.1, vary=False)    piece_size = self.stepwise_size    for t in range(piece_size, len(data), piece_size):        params.add(f"t{t}_q", value=0.5, min=0, max=0.99, brute_step=0.1, vary=True)    return params

2. Добавляемновый методget_step_rt_beta, который вычисляет(t)иRtиспользуя ступенчатую функцию.

3. В методеstepиспользуемget_step_rt_beta,(t)с учетом карантина в деньt.

Учитываем неполноту статистики

Известно, что не все случаи заражения попадают в официальную статистику. Многие люди просто не обращаются в больницы, а часть болеет безсимптомно. Мы можем учесть это в нашей модели.

Изменим группы Infected, Recovered и Dead следующим образом:

• Iv(t) распространяют вирус, видны в статистике.

• I(t) распространяют вирус, не видны в статистике.

• Rv(t) переболели, видны в статистике.

• R(t) переболели, не видны в статистике.

• Dv(t) погибли, видны в статистике.

• D(t) погибли, не видны в статистике.

Добавим два новых параметра:

• pi вероятность попадания зараженного в статистику, то есть в группу Iv.

• pd вероятность регистрации смерти зараженного, невидимого в статистике, в группу Dv. Зараженные из группы Iv всегда попадают в Dv, но этот параметр позволяет учесть, что часть смертей незамеченных зараженных так же попадает в статистику по смертности от COVID-19.

В конце концов мы приходим к такой схеме модели, где вершины графа это группы, а стрелки определяют как люди перемещаются между группами:

Для удобства будем называть её SEIRD со скрытыми состояниями. Она отличается от обычного SEIRD только дифференциальными уравнениями. Подробности реализации можно посмотретьздесь.

Оптимизируем параметры

Модель готова, осталось только подобрать оптимальные значения параметров так, чтобы выдаваемые моделью Iv(t), Rv(t) и Dv(t) были как можно ближе к историческим данным. Это можно сделать с помощью метода наименьших квадратов.

Конкретнее, используем для этогоlmfit.minimize. Эта функция получает на вход callable, возвращающий список отклонений предсказаний от истины (остатки, residuals), и подбирает такие параметры, чтобы сумма отклонений была минимальна. Используем алгоритм Levenberg-Marquardt (method=leastsq), который двигает каждый параметр на очень маленькое значение и проверяет, уменьшается ли ошибка. Важно знать, что по этой причине алгоритм не может оптимизировать дискретные параметры, поэтому, например, не сможет подобратьstepwise_size.

Я написалнебольшую оберткувокруг этой функции под названиемBaseFitter, которая сохраняет историю параметров и делает другие служебные вещи. Но в конце концов всё сводится к вызовуminimize.

Рассмотрим получение остатков подробнее:

def smape_resid_transform(true, pred, eps=1e-5):    return (true - pred) / (np.abs(true) + np.abs(pred) + eps)class HiddenCurveFitter(BaseFitter):...    def residual(self, params, t_vals, data, model):        model.params = params    initial_conditions = model.get_initial_conditions(data)    (S, E, I, Iv, R, Rv, D, Dv), history = model.predict(t_vals, initial_conditions, history=<span class="hljs-literal" style="box-sizing: border-box; color: rgb(0, 134, 179);">False</span>)    (new_exposed,     new_infected_invisible, new_infected_visible,     new_recovered_invisible,     new_recovered_visible,     new_dead_invisible, new_dead_visible) = model.compute_daily_values(S, E, I, Iv, R, Rv, D, Dv)    new_infected_visible = new_infected_visible    new_dead_visible = new_dead_visible    new_recovered_visible = new_recovered_visible    true_daily_cases = data[self.new_cases_col].values[<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>:]    true_daily_deaths = data[self.new_deaths_col].values[<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>:]    true_daily_recoveries = data[self.new_recoveries_col].values[<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>:]    resid_I_new = smape_resid_transform(true_daily_cases, new_infected_visible)    resid_D_new = smape_resid_transform(true_daily_deaths, new_dead_visible)    resid_R_new = smape_resid_transform(true_daily_recoveries, new_recovered_visible)    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">if</span> self.weights:        residuals = np.concatenate([            self.weights[<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'I'</span>] * resid_I_new,            self.weights[<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'D'</span>] * resid_D_new,            self.weights[<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'R'</span>] * resid_R_new,        ]).flatten()    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">else</span>:        residuals = np.concatenate([            resid_I_new,            resid_D_new,            resid_R_new,        ]).flatten()    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">return</span> residuals

Основная идея простая. Для промежуткаt_valsполучаем предсказания заражений в день, смертей в день и выздоровлений в день. Вычисляем отклонения предсказанных значений от реальных и кладем всё это в один массив.

Также применяется несколько трюков, до которых пришлось дойти методом эмпирического тыка.

Во-первых, количество новых зараженных и количество новых смертей имеют абсолютно разный масштаб: нормально видеть 1000 новых зараженных за день и лишь 1 смерть в тот же день. Если ориентироваться на абсолютные остатки, то оптимизация будет обращать внимание только на остатки по зараженным. Чтобы привести все остатки к одному масштабу я применил трансформацию котносительной ошибкес помощью функцииsmape_resid_transform.

Во-вторых, можно предположить, что статистика по смертям достовернее статистики по заражениям и выздоровлениям. Чтобы внести это предположение в оптимизацию вводятся веса (self.weights) для остатков. Хорошо сработали такие веса:0.5для остатков по смертям и0.25для остальных.

Тренируем модель

Наконец, натренируем SEIRD со скрытыми состояниями. Статистика вносится с некоторым лагом по времени. Чтобы это меньше влияло на нашу модель будем использовать сглаженные версии колонок для расчета остатков.

from sir_models.fitters import HiddenCurveFitterfrom sir_models.models import SEIRHiddenstepwize_size = 60weights = {    'I': 0.25,    'R': 0.25,    'D': 0.5,}model = SEIRHidden(stepwise_size=stepwize_size)fitter = HiddenCurveFitter(     new_deaths_col='deaths_per_day_ma7',     new_cases_col='infected_per_day_ma7',     new_recoveries_col='recovered_per_day_ma7',     weights=weights,     max_iters=1000,)fitter.fit(model, train_subset)result = fitter.resultresult

Полученные параметры:

Проверим фит модели на тренировочных данных:

Вот это другое дело! Модель показывает, что почти половина смертей не попадают в статистику. По заражениям другая картина: в статистику попадает примерно одна пятая.

Полезно взглянуть на все предсказываемые моделью компоненты: ежедневные смерти, заражения и выздоровления. Модель может выдавать абсолютно верные предсказания по смертям, но при этом прогнозировать в 10 раз больше выздоровлений, или наоборот.

Отлично, здесь тоже всё выглядит нормально. В некоторых местах модель достаточно сильно ошибается, но нас это устраивает: нам не нужен идеальный фит на тренировочных данных. Модель верно захватывает общий тренд. Значит, она нормально описывает механизм эпидемии, даже если ошибается в конкретные дни.

Верифицируем модель c помощью кросс-валидации

Хорошие графики на тренировочных данных это одно, а качество на реальных данных это совсем другое. Хоть мы и не строим прогнозную модель, но единственный способ понять, что она способна описать эпидемию это проверить качество прогноза.

В статьях по моделированию COVID-19 стандартом является проверять модель на предсказаниях кумулятивного числа погибших, поэтому сделаем так же.

Процесс верификации:

• Выбираем несколько дат из тренировочных данных, например каждый 20-ый день.

• Для каждой даты:

  • Обучаем модель на всех данных до этой даты.

  • Делаем прогноз суммарного количества погибших на 30 дней вперед.

  • Считаем ошибку предсказания.

• Вычисляем среднюю ошибку модели.

• Сравниваем ошибку с ошибкой бейзлайна: предсказанием предыдущего дня.

from sir_models.utils import eval_on_select_dates_and_k_days_aheadfrom sir_models.utils import smapefrom sklearn.metrics import mean_absolute_errorK = 30last_day = train_subset.date.iloc[-1] - pd.to_timedelta(K, unit='D')eval_dates = pd.date_range(start='2020-06-01', end=last_day)[::20]def eval_hidden_moscow(train_df, t, train_t, eval_t):    weights = {        'I': 0.25,        'R': 0.25,        'D': 0.5,    }    model = SEIRHidden()    fitter = HiddenCurveFitter(        new_deaths_col='deaths_per_day_ma7',        new_cases_col='infected_per_day_ma7',        new_recoveries_col='recovered_per_day_ma7',        weights=weights,        max_iters=1000,        save_params_every=500)    fitter.fit(model, train_df)train_initial_conditions = model.get_initial_conditions(train_df)train_states, history = model.predict(train_t, train_initial_conditions, history=<span class="hljs-literal" style="box-sizing: border-box; color: rgb(0, 134, 179);">False</span>)test_initial_conds = [compartment[-<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>] <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">for</span> compartment <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">in</span> train_states]test_states, history = model.predict(eval_t, test_initial_conds, history=<span class="hljs-literal" style="box-sizing: border-box; color: rgb(0, 134, 179);">False</span>)    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">return</span> model, fitter, test_states(models, fitters, model_predictions, train_dfs, test_dfs) = eval_on_select_dates_and_k_days_ahead(train_subset,                                        eval_func=eval_hidden_moscow,                                        eval_dates=eval_dates,                                        k=K)model_pred_D = [pred[7] for pred in model_predictions]true_D = [tdf.total_dead.values for tdf in test_dfs]baseline_pred_D = [[tdf.iloc[-1].total_dead]*K for tdf in train_dfs]overall_errors_model = [mean_absolute_error(true, pred) for true, pred in zip(true_D, model_pred_D)]overall_errors_baseline = [mean_absolute_error(true, pred) for true, pred in zip(true_D, baseline_pred_D)]print('Mean overall error baseline', np.mean(overall_errors_baseline).round(3))print('Mean overall error model', np.mean(overall_errors_model).round(3))overall_smape_model = [smape(true, pred) for true, pred in zip(true_D, model_pred_D)]np.median(overall_smape_model)

Результат:

• Mean Absolute Arror бейзлайна: 714.

• Mean Absolute Arror модели: 550.

• Symmetric Mean Absolute Percentage Error: 4.6%

Модель побила бейзлайн и ошибка в пределах 5% Можно заключить, что модель подходит для оценки ситуации с двумя штаммами.

Расширение модели для двух штаммов

Ситуацию с двумя штаммами можно смоделировать как две параллельные эпидемии, которые заражают людей из одной популяции. При этом зараженные одним штаммом не могут заразиться другим. В контексте SEIR моделей это означает, например, что вместо группыI(t)для всех зараженных у нас теперь будет две группы для разных штаммов:I1(t)иI2(t).

Модель эпидемии обычного SARS-CoV-2 мы только что обучили. Осталось получить модель для британского штамма. СогласноисследованиямB.1.1.7 отличается только параметромR0: он в 1.4 - 1.9 раз больше. Значит мы можем взять обученную модель для обычного SARS-CoV-2, увеличитьR0и получить модель для B.1.1.7.

Полная схема модели с двумя штаммами выглядит так:

Главное отличиереализации этой модели в кодеот предыдущих в том, что эта модель не обучается, а создается на основе обученной модели для одного штамма. Добавляется параметрbeta2_mult, который задает, на сколько надо умножить1(t)первого штамма, чтобы получить2(t)второго.

class SEIRHiddenTwoStrains(SEIRHidden):    ...    @classmethod    def from_strain_one_model(cls, model):        strain1_params = model.params        strain1_params.add("beta2_mult", value=1.5, min=1, max=2, vary=False)        return cls(params=deepcopy(strain1_params))

Мы предполагаем, что меры карантина и всё остальное действуют на новый штамм так же, как на старый, а так же что уровень карантина в будущем сохраняется такой же, как в последний день тренировочных данных.

Сценарии развития эпидемии при появлении B.1.1.7

Наконец, мы получили модель для двух штаммов. Неизвестно какое именноR0у B.1.1.7, а так же сколько сейчас в Москве носителей этого штамма. Рассмотрим несколько сценариев.

ЗначениеR0влияет на дату и высоту пика заражений новым штаммом. В лучшем случае B.1.1.7 просто продлевает эпидемию в текущем состоянии на полгода-год. В худшем случае случается большая волна заражений.

Изначальное количество носителей британского штамма влияет только на дату пика. Во всех случаях новая волна случается в течение года после появления нового штамма.

Ослабление карантина приводит к значительному количеству заражений.

Автор статьи: Борис Цейтлин

Узнайте подробнее о курсе "Machine Learning. Professional".

Также приглашаем всех желающих на бесплатный интенсив по теме: Деплой ML модели: от грязного кода в ноутбуке к рабочему сервису.