Отслеживание случаев Covid-19 в Швейцарии в режиме реального времени

Рюта Йошимацу, Саймон Хефти

Введение

Когда мы говорим о Covid-19, мы иногда ссылаемся на действующий номер воспроизведения, Rt. Rt — один из ключевых показателей, позволяющих оценить наше положение с точки зрения роста эпидемии. Проще говоря, он показывает, сколько людей в среднем заражает зараженный человек в момент времени t. Когда это значение выше 1, число заражений будет расти экспоненциально, а когда оно ниже 1, ожидается обратное. Получить оценку Rt в режиме реального времени сложно из-за неэффективности отчетности. Например, последняя оценка Rt, доступная сегодня, предоставленная Swiss Science Task Force, была сделана 2 недели назад.

В этой записи блога мы рассмотрим один из возможных способов восполнить этот пробел и дать оценку сегодняшнему Rt. Сначала мы рассмотрим общее количество зарегистрированных случаев Covid-19 в Швейцарии [данные] с начала эпидемии и изучим динамику волн заражения. Затем мы смоделируем временной ряд и сделаем прогнозы на ближайшее прошлое (описано ниже). Поскольку мы знаем, что Rt является функцией ежедневного количества случаев вместе с предположениями модели, мы будем использовать наш прогноз для оценки Rt пропущенных дней, включая сегодняшний день.

Мы признаем, что количество зарегистрированных случаев, которые мы использовали в этой работе [данные], служит только приблизительным показателем фактического количества инфекций. Эти цифры иногда в значительной степени занижались из-за перегрузки возможностей тестирования, а иногда завышались из-за ошибок и предвзятости тестирования. Более того, оценки Rt, опубликованные Швейцарской научной целевой группой, по которым мы проверили наш прогноз, являются оценкой, основанной на еще одной оценке. Мы ни в коем случае не пытаемся заявить, что разработали мощную модель, которая фиксирует полную динамику инфекций Covid-19, или идеальную модель прогнозирования для оценки неизвестного Rt. Мы всего лишь два любопытных специалиста по данным. пытаясь осмыслить имеющиеся у нас данные и поделиться тем, что мы узнали. Наслаждаться!

Понимание зарегистрированных случаев Covid-19 в Швейцарии

На верхнем рисунке показано изменение во времени общего (кумулятивного) числа зарегистрированных случаев Covid-19, а внизу — общее количество зарегистрированных смертельных случаев из-за Covid-19 в Швейцарии [данные]. Глядя на графики, мы замечаем периоды роста и плато, которые следуют друг за другом.

Этот волнообразный рост населения часто моделируется с помощью логистического уравнения. Логистические кривые могут удобно отображать динамику экспоненциального роста, за которым следует линейный рост и, в конечном итоге, насыщение (плато). Для Covid-19 уже было опубликовано множество исследований с использованием того или иного варианта этого уравнения, другими словами, обобщенных логистических уравнений: напр. бумага1, бумага2.

В данной работе мы использовали растянутое логистическое уравнение (уравнение ниже), которое также относится к семейству обобщенных логистических уравнений. По сути, растянутое логистическое уравнение кодирует временную зависимость в экспоненциальной скорости роста. Здесь предполагается, что темпы роста со временем снижаются из-за внешних воздействий, таких как вмешательство государства (изоляция) или повышение осведомленности людей. Исследователи показали, что растянутая логистическая кривая довольно хорошо подходит для случаев Covid-19, а также для распространения других эпидемий.

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

Подгонка выполняется методом наименьших квадратов, а количество растянутых логистических условий, которые должны быть включены (i в приведенном выше уравнении) в модель в момент времени t, определяется путем сравнения квадратичная ошибка, создаваемая отдельными моделями, включающими от 1 до m терминов. Наилучшая подходящая модель по состоянию на 15 марта 2021 г. содержит три растянутые логистические кривые. Анимация ниже показывает процесс выбора динамической модели. Серая и синяя кривые представляют наблюдаемое общее количество случаев, а оранжевая кривая лучше всего соответствует приведенному выше уравнению. Оранжевые вертикальные пунктирные линии обозначают t1, t2 и t3. Обратите внимание, что подгонка больше отклоняется от фактического наблюдения в начале каждого роста. Общеизвестно, что хороший прогноз можно сделать только во второй половине сигмоиды. К счастью, сегодня (15 марта 2021 г.) совпадение кажется стабильным, но в начале нового периода роста нам, возможно, придется найти более надежную модель.

Ниже показаны снимки подгонки по состоянию на 15 марта 2021 г. с использованием линейной комбинации трех растянутых логистических кривых. Эта модель неплохо справляется с общей тенденцией.

Оценка последнего Rt

Rt на сегодняшний день может быть выведен из модели передачи болезни, начального распределения вероятностей Rt и истории ежедневного количества случаев, включая сегодняшнее, с использованием Байесовский каркас. Черная кривая на графике ниже — это оценочное Rt, полученное с использованием этого подхода, а красная кривая — логарифм ежедневного числа случаев. На графике видно, что как только Rt превышает пороговое значение, количество случаев начинает увеличиваться, и наоборот, когда оно становится ниже.

Трудность в оценке последнего Rt заключается в получении ежедневного количества случаев на сегодняшний день (красный прямоугольник на рисунке выше). Существует неизбежный временной лаг между фактическим случаем заражения и моментом сообщения об этом событии. По этой причине Швейцарская научная рабочая группа обычно ждет около 10–14 дней, пока сообщаемое число не стабилизируется, а затем оценивает Rt для данного дня. Вот почему Rt, представленный на сайте сегодня, отражает динамику заражения только 10–14 дней назад.

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

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

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

Мы используем красную кривую как тренд и вычитаем ее из черной кривой. Оставшаяся часть представляет собой временной ряд без тренда (рисунок выше). Здесь вы также можете четко увидеть сильную недельную сезонность (колебания) и волатильность (зависимая от времени дисперсия). Чтобы дать хороший прогноз, нам также необходимо зафиксировать сезонные компоненты в нашей модели.

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

Мы выполнили обобщенную линейную регрессию на временных рядах без тренда с ошибками ARIMA и членами Фурье базовой периодичности 7 (неделя) и 30,4 (месяц) в качестве дополнительных внешних регрессоров. Выбор модели (порядок АР, порядок МА, количество слагаемых Фурье на базовую частоту, амплитуды гармонических слагаемых и т. д.) осуществлялся с использованием информационных критериев Акаике. Затем мы использовали модель наилучшего соответствия, чтобы сделать прогноз для временных рядов без тренда за последние две недели (вверху справа), наложили эти значения на общий тренд, полученный из растянутых логистических кривых, и, наконец, пришли к нашему прогнозу ближайшего прошлого. (красная волнистая кривая на нижнем графике). Посетите наш репозиторий на GitHub, чтобы узнать подробности реализации.

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

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

Выше приведен снимок набора данных, который мы создали с использованием данных временных рядов: то есть время, Rt и количество ежедневных случаев. Мы разработали функции задержки, которые представляют собой Rt и ежедневные наблюдения на предыдущих временных шагах, а также функции скользящего среднего, которые представляют собой совокупность Rt и ежедневных наблюдений за фиксированное окно предыдущих временные шаги. Эти функции предоставляют информацию о последовательных корреляциях между значениями в наборе данных. Rt — это наша переменная ответа, а все остальное, кроме времени, — это наш предиктор. Мы выполнили регрессию Пуассона, используя алгоритм повышения градиента. Посетите наш репозиторий на GitHub, чтобы узнать подробности реализации.

Мы проверили модель повышения градиента, используя метод продвижения вперед. Верхний график на приведенном выше рисунке показывает среднеквадратичную ошибку (RMSE) предсказанного Rt для следующих 14 дней со дня x, где x находится в диапазоне от 1 апреля, с 2020 г. по 1 марта 2021 г. В день x модель обучается с использованием данных, доступных только до этого дня. Мы видим, что в начале эпидемии модель изо всех сил пытается сделать хороший прогноз со RMSE > 0,2, но по мере того, как проходит время и становится доступным больше данных, она начинает изучать функцию, и метрика стабилизируется. Он снова пострадал в начале второго роста, но в течение третьего периода роста и после этого прогнозы сработали со среднеквадратичной ошибкой ‹0,005.

И, наконец, наш долгий путь восполнения недостающего Rt подошёл к концу. На верхнем рисунке красная волнистая кривая на вершине черной кривой, которая представляет собой эффективную Rt, рассчитанную Швейцарской научной группой, является нашим прогнозом на последние 14 дней, включая сегодняшний день (15 марта). , 2021). Чтобы проверить этот прогноз на данных в реальном времени, мы терпеливо ждали около 2 недель и проверили эффективную Rt, опубликованную Швейцарской научной группой за эти дни. Результаты были положительными. Оценочный Rt вырос (хотя и не так сильно, как мы предсказывали), и наш прогноз находится в пределах доверительного интервала: наш прогноз 1,25, фактическая оценка 1,13 (1–1,27) на 15 марта 2021 г.

Вывод

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

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

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

Наконец, мы хотели бы повторить, что цель этой работы состоит не в том, чтобы объяснить науку об эпидемиях или предложить более совершенные модели для оценки Rt, а скорее в том, чтобы продемонстрировать, как мы объединили устоявшиеся и новые методы моделирования. на предоставленных нам данных.

Подтверждение

Мы благодарим Lucas Böttcher за его экспертное мнение по теме и Jacqueline, Vito и Marcel из d-one за отзывы о нашей работе.