Оценка уровня избыточной смертности, ассоциированной с COVID-19 в Республике Беларусь

Главное в этой статье

• Выполнен статистический анализ многолетних месячных данных по смертности в Республике Беларусь с целью оценить уровень избыточной смертности в период пандемии COVID-19.
• Наиболее вероятная суммарная избыточная смертность в стране в апреле-июне 2020 г. оценена в 6737 человек, что значительно превышает ранее опубликованную оценку в 5500 человек.
• Наиболее вероятный удельный прирост смертности в апреле-июне 2020 г. оценен в 23%. В пересчете на 100 тыс. населения смертность за этот период времени возросла на 71 человека.
• Полученные оценки избыточной смертности в Беларуси сравнимы с таковыми в странах, существенно пострадавших от пандемии COVID-19. 

Введение

В отличие от предыдущих лет, Министерство Здравоохранения Республики Беларусь (РБ) пока не опубликовало официальные данные по количеству смертей в стране за 2020 г. Однако, как выяснили недавно журналисты новостного издания "Наша Нiва", такие месячные данные были переданы Минздравом РБ в ООН и их можно найти на сайте UNdata. Журналисты "Наша Нiва" проанализировали эти данные и сделали вывод о том, что в апреле-июне 2020 г. прирост смертности в стране составил примерно 5500 человек, по сравнению со средним уровнем смертности в эти же месяцы за последние 5 лет. Предположительно, эта цифра дает намного более реалистичную оценку избыточной смертности в апреле-июне 2020 г., связанной с пандемией новой короновирусной инфекции, чем официальная оценка Минздрава РБ в 398 человек.

Методологический недостаток опубликованного анализа

Хотя использованный журналистами "Наша Нiва" метод логичен, прост и часто применяется в подобных ситуациях, он имеет один существенный недостаток. Дело в том, что сравнение со средним показателем смертности за предыдущие годы имело бы смысл только при отсутствии тренда в данных. Однако, как можно видеть на рис. 1, в последнее время (примерно с 2005 г.) в Беларуси наблюдался довольно выраженный тренд на снижение смертности. В связи с этим более корректным было бы сравнение зарегистрированной смертности в апреле-июне 2020 г. с естественной смертностью, которая имела бы место без вспышки короновирусной инфекции (т.е. если бы 2020 г. был обычным годом). Ниже описана попытка оценить подобные естественные уровни смертности в апреле, мае и июне 2020 г. с использованием прогнозной модели, учитывающей тренд и другие свойства данных в период до наступления вспышки COVID-19.

Рис. 1. Динамика смертности в Республике Беларусь в период с января 1980 г. по июнь 2020 г. (включительно). Голубая линия (получена методом loess-регрессии) добавлена для иллюстрации тренда. Обратите внимание на пропуски в данных за 2000, 2001, 2006 и 2010 гг.

Прогнозирование естественного уровня смертности

Данные и код, использованные для получения описанных ниже результатов, можно найти в  Github-репозитории ranalytics/covid_by. С целью экономии места код для построения графиков ниже не приводится.

Разобьем имеющиеся данные (рис. 1) на следующие три части:

• обучающая выборка (с января 1980 г. по сентябрь 2019 г.) - для обучения альтернативных прогнозных моделей;
• проверочная выборка (с октября 2019 г. по март 2020 г.) - для проверки качества предсказаний моделей, подогнанных к данным из обучающей выборки;
• тестовая выборка (с апреля 2020 г. по июнь 2020 г.): данные за период пандемии COVID-19, которые будут использованы для расчета избыточной смертности.

require(dplyr)
require(readr)
require(tsibble)
require(feasts)
require(lubridate)
require(ggplot2)
require(tidyr)
require(RColorBrewer)

# Загрузка и очистка данных:
dat <- read_csv("./data/un_data.csv") %>% 
  setNames(., c("year", "month", "y")) %>% 
  dplyr::filter(!month %in% c("Total", "Unknown"))

# Дополнительная подготовка данных к анализу 
# (добавление дат в формате "год месяц", преобразование в
# объект класса tsibble и замена пропущенных значений на NA):
dat <- dat %>% 
  mutate(dm = paste(year, substr(month, 1, 3)) %>% yearmonth(.)) %>% 
  as_tsibble(., index = dm) %>% 
  fill_gaps() %>% 
  mutate(y = y /10000) %>% # смертность, выраженная в десятках тыс. человек
  dplyr::select(y, dm)

# Данные за период до апреля 2020 г.:
pre_covid <- dat %>% 
  filter(dat$dm < yearmonth("2020 Apr"))

# Создание обучающей, проверочной и тестовой выборок:
n_pre_covid <- nrow(pre_covid)
train <- pre_covid[1:(n_pre_covid - 6), ]
valid <- pre_covid[(n_pre_covid - 5):n_pre_covid, ]
test <- dat %>% dplyr::filter(dat$dm >= yearmonth("2020 Apr"))

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

Как видно на рис. 2, месячная смертность в Беларуси имеет четко выраженную сезонность с пиковыми уровнями, приходящимися обычно на декабрь, январь и март.

Рис. 2. Месячные уровни смертности в Республике Беларусь в разные годы наблюдений

Кроме того, в данных наблюдается умеренная  автокорреляция до сдвига, примерно равного трем (на последующих сдвигах она также имеет место, но становится заметно слабее; рис. 3).

Рис. 3. Автокорреляция в данных по месячной смертности в Республике Беларусь (выражена в десятках тыс. человек) на сдвигах от 1 до 6. Обозначения: y - смертность, lag(y, n) - смертность на сдвиге n

Таким образом, модель для предсказания уровня месячной смертности в Беларуси должна учитывать имеющиеся в данных тренд, сезонность и автокорреляцию. Для построения подобных моделей подходят многие методы. Однако учитывая наличие пропущенных значений в данных (см. рис. 1), мы воспользуемся методологией "байесовских структурных моделей временных рядов", для которой такие пропуски не являются критичными (подробнее об этом семействе моделей см. гл. 7 в моей новой книге  по анализу временных рядов).

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

require(bsts)

# Моделируемый отклик:
y <- train$y

# Спецификация компонент модели:
ss <- list()
ss <- AddLocalLinearTrend(ss, y)
ss <- AddSeasonal(ss, y, nseasons = 12)
ss <- AddAutoAr(ss, y, lag = 3)

# Подгонка модели:
m0 <- bsts(y, ss, niter = 2000, ping = 0, seed = 42)

Как следует из приведенных ниже количественных показателей качества полученной модели, она хорошо описывает обучающие данные (об этом, в частности, говорит высокий коэффициент детерминации rsquare):

summary(m0)

## $residual.sd
## [1] 0.01490257
## 
## $prediction.sd
## [1] 0.06201056
## 
## $rsquare
## [1] 0.9911943
## 
## $relative.gof
## [1] NA

О хорошей аппроксимации обучающих данных моделью m0 говорит также рис. 4.

 Рис. 4. Результат подгонки модели m0. Кружками показаны обучающие данные. Мелкими полупрозрачными черными точками показано апостериорное совместное распределение компонент модели на каждой временной отметке (на каждую отметку приходится около 2000 тыс. точек). Области плотного расположения этих точек выглядят темнее и соответствуют областям с большей апостериорной вероятностью

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

Рис. 5. Апостериорные распределения компонент модели m0 (слева направо: тренд, сезонная и авторегрессионная компоненты)

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

ss <- list()
ss <- AddLocalLinearTrend(ss, y)
ss <- AddSeasonal(ss, y, nseasons = 12)
m1 <- bsts(y, ss, niter = 2000, ping = 0, seed = 42)

summary(m1)

## $residual.sd
## [1] 0.05366628
## 
## $prediction.sd
## [1] 0.06683368
## 
## $rsquare
## [1] 0.8858062
## 
## $relative.gof
## [1] NA

Как видим, по всем количественным показателям модель m1 дает более низкое качество аппроксимации обучающих данных, чем модель m0. Это можно увидеть также из следующего графика (рис. 6):

   

Рис. 6. Сравнение качества аппроксимации обучающих данных моделями m0 и m1 с помощью т.н. "ошибок следующего шага" (англ. "one-step-ahead errors"; подробнее см. разд. 7.4 в книге Мастицкий, 2020). Вверху: накопленные средние абсолютные ошибки следующего шага. Внизу: обучающие данные. Более низкий уровень ошибок, накопленных моделью m0, говорит о ее более высоком качестве

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

Чтобы выполнить такой независимый тест воспользуемся моделями m0 и m1 и рассчитаем прогнозные значения смертности на следующие 6 месяцев, которые соответствуют временному периоду созданной нами ранее проверочной выборки (с октября 2019 г. по март 2020 г.; см. выше). Полученные прогнозные значения далее сравним с фактически наблюденными в этот период уровнями смертности. Для измерения качества предсказаний воспользуемся средней абсолютной ошибкой:

# Функция для расчета средней абсолютной ошибки:
mae <- function(observed, predicted){
  mean(abs(observed - predicted))
}

# Расчет средних абсолютных ошибок предсказаний для обеих моделей:
sapply(list("m0" = (predict(m0, horizon = 6) %>% .$median)*10000, 
            "m1" = (predict(m1, horizon = 6) %>% .$median)*10000), 
       mae, observed = valid$y * 10000) %>% round()
       
##  m0  m1 
## 415 508

Модель m0 ошиблась в своих предсказаниях для проверочной выборки в среднем на 415 смертей (см. также рис. 7), что намного меньше ошибки, допущенной моделью m1. За не имением дополнительных предикторов, которые потенциально могли бы еще больше повысить качество предсказаний m0, будем считать, что структура этой модели оптимальна для аппроксимации анализируемых данных.

Рис. 7. Иллюстрация качества предсказаний модели m0 для проверочной выборки (синяя линия). Светлоголубая полоса вокруг синей линии соответствует 95%-ной доверительной области полученных предсказаний. Черная кривая изображает исторические наблюдения (для удобства визуализации показаны только данные за последние 5 лет - с октября 2015 г. по сентябрь 2019 г.). 

Объединим теперь данные из обучающей и проверочной выборок и подгоним окончательную модель для прогнозирования смертности. Структура этой модели аналогична таковой у m0:

y_full <- c(y, valid$y)

ss <- list()
ss <- AddLocalLinearTrend(ss, y_full)
ss <- AddSeasonal(ss, y_full, nseasons = 12)
ss <- AddAutoAr(ss, y_full, lag = 3)
m_final <- bsts(y_full, ss, niter = 2000, ping = 0, seed = 42)

Вычислим прогнозные уровни смертности на период с апреля по июнь 2020 г. с помощью полученной модели m_final и изобразим результат графически (рис. 8):

# Вычисление прогнозных значений:
m_final_pred <- predict(m_final, horizon = 3)

Рис. 8. Прогноз естественного уровня месячной смертности в Республике Беларусь в апреле-июне 2020 г. (синяя линия), который имел бы место в отсутствие вспышки короновирусной инфекции. Светлоголубая полоса вокруг синей линии соответствует 95%-ной доверительной области полученных предсказаний. Черная кривая изображает исторические наблюдения (для удобства визуализации показаны только данные с апреля 2017 г. по март 2020 г.). Красными точками показаны фактически наблюденные уровни смертности во время пандемии

На рис. 8 хорошо видно, что зарегистрированные в действительности уровни смертности (показаны красными точками) намного превышают прогнозные значения (синяя линия). Рассчитаем теперь величину этого прироста.

Оценки уровней избыточной смертности

Важно подчеркнуть, что термин "избыточная смертность" ("excess mortality" в англоязычных источниках) используется здесь для описания всех дополнительных случаев смерти, имевших место в Беларуси в апреле-июне 2020 г. В отсутствие другой информации (в частности, разбивки по причинам смерти) невозможно утверждать, что все эти дополнительные случаи были непосредственно вызваны COVID-19: некоторые из них, например, могли быть обусловлены вторичными осложнениями от перенесенной короновирусной инфекции или какими-то совершенно другими тяжелыми заболеваниями (сердечно-сосудистые, рак, и т.п.), которым система здравоохранения не смогла уделить должного внимания в этот период из-за карантинных мероприятий и нехватки персонала и других ресурсов. Наконец, это могло быть и увеличение количества случаев, никак не связанных с медицинскими причинами (например, мгновенная смерть в результате участившихся тяжелых дорожно-транспортных происшествий, суициды и т.п.).

В простейшем случае избыточную смертность можно было бы рассчитать путем вычета спрогнозированных точечных оценок смертности (рис. 8, голубая линия) из фактически наблюденных уровней (рис. 8, красные точки). Однако благодаря байесовкой природе построенной модели у нас есть возможность рассчитать оценки избыточной смертности с учетом неопределенности в отношении прогнозных значений. Всего при подгонке модели было выполнено 2000 итераций по методу Монте-Карло с марковскими цепями (MCMC). После отброса итераций, потраченных на "разогрев" алгоритма (т.н. "burn-in period"), эффективное число итераций, доступное нам для выполнения расчетов, составило 1947:

# Для удобства вернемся с исходной шкале, на которой
# была зарегистрирована месячная смертность (т.е. умножим
# каждое прогнозное значение на 10000):
m_dist <- (m_final_pred$distribution*10000) %>% round()

dim(m_dist)
## [1] 1947    3

head(m_dist)
##       [,1]  [,2] [,3]
## [1,]  9754 10741 8785
## [2,] 10705 11378 8820
## [3,]  9663 11118 9161
## [4,]  9808  9861 9865
## [5,] 10049 10454 9739
## [6,] 10598 10333 9884

Столбцы в приведенной матрице соответствуют апрелю, маю и июню 2020 г., а каждая из 1947 строк - возможной реализации уровня смертности в эти месяцы согласно полученной прогнозной модели. Таким образом, для вычисления возможных реализаций избыточной смертности достаточно вычесть каждую строку матрицы из вектора с фактически наблюденными уровнями смертности:

covid_y <- test$y * 10000
excess_mortality <- apply(m_dist, 1, 
                          FUN = function(x){covid_y - x}) %>% t()
head(excess_mortality)
##      [,1] [,2] [,3]
## [1,] 1132 1225 4231
## [2,]  181  588 4196
## [3,] 1223  848 3855
## [4,] 1078 2105 3151
## [5,]  837 1512 3277
## [6,]  288 1633 3132

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

excess_mortality <- excess_mortality %>% 
  as_tibble(.name_repair = "minimal") %>% 
  setNames(., c("apr", "may", "jun")) %>% 
  rowwise() %>% 
  mutate(total = sum(c(apr, may, jun)))
  
excess_mortality
  
## # A tibble: 1,947 x 4
## # Rowwise: 
##      apr   may   jun total
##    <dbl> <dbl> <dbl> <dbl>
##  1  1131 1224.  4231 6586.
##  2   181  589.  4197 4967.
##  3  1222  849.  3855 5926.
##  4  1077 2105.  3151 6333.
##  5   836 1511.  3277 5624.
##  6   287 1633.  3131 5051.
##  7   542  838.  2766 4146.
##  8  1279 1631.  4255 7165.
##  9  1488 2440.  4497 8425.
## 10  1638 2245.  4313 8196.
## # ... with 1,937 more rows

Вычислим теперь наиболее вероятные (медианные) уровни избыточной смертности и их 95%-ные доверительные интервалы для каждого месяца в отдельности и по всем трем месяцам в сумме:

result_abs <- excess_mortality %>%
  dplyr::select(apr:total) %>% 
  pivot_longer(cols = apr:total) %>% 
  group_by(month = name) %>% 
  summarise(tibble(
    med = median(value),
    ll95 = quantile(value, p = 0.025),
    ul95 = quantile(value, p = 0.975)
  )) %>% 
  arrange(med) %>% 
  mutate(month = c("апрель", "май", "июнь", "всего"))

result_abs
## # A tibble: 4 x 4
##   month    med  ll95  ul95
##   <chr>  <dbl> <dbl> <dbl>
## 1 апрель  963  -113. 2018.
## 2 май    1874.  729. 3062.
## 3 июнь   3906  2722. 5090.
## 4 всего  6731. 4189. 9335.

Итак, наиболее вероятная избыточная смертность в Беларуси в апреле, мае и июне 2020 г. составила 963, 1874 и 3906 человек соответственно. Интересно, что нижний доверительный предел у апрельской оценки отрицательный (-113 человек). Это говорит о том, что зарегистрированная в апреле смертность статистически не отличается от прогнозного значения по этому месяцу. Другими словами, наблюденная в апреле смертность вполне укладывается в диапазон апрельских уровней, имевших место в предыдущие годы и учтенных в ходе построения прогнозной модели. Этот результат неудивителен, учитывая что на апрель 2020 г. пришлось начало пандемии, когда количество заболевших новым короновирусом было еще невелико. Однако со временем оно неизбежно возросло, что по-видимому привело также к значительному росту избыточной смертности в мае и июне. Полученные оценки избыточной смертности для апреля-июня и их 95%-ные доверительные интервалы представлены графически на рис. 9:

Рис. 9. Точечные оценки избыточной смертности в Беларуси в апреле-июне 2020 г. и их 95%-ные доверительные интервалы

В сумме за три рассматриваемых месяца наиболее вероятный уровень избыточной смертности составил 6731 человек, что более чем на 1000 смертей превышает величину, полученную журналистами издания "Наша Нiва" (5500 человек). Согласно приведенному выше доверительному интервалу этой оценки, истинный суммарный уровень избыточной смертности в апреле-мае с вероятностью 95% находился в диапазоне от 4189 (наиболее оптимистичная оценка) до 9335 человека (наиболее пессимистичная оценка). 

В пересчете на 100 тыс. населения получается, что истинный суммарный прирост смертности в апреле-июне 2020 г. с вероятностью 95% лежал в интервале от 44 до 99 человек, а наиболее вероятное значение составило 71 (этот вывод основан на допущении, что по состоянию на середину 2020 г. население Беларуси составляло 9449323 человека - оценка, полученная по данным ООН и взятая с сайта Worldometer). По этому показателю Беларусь оказывается сравнимой с несколькими странами, которые больше всего пострадали от пандемии COVID-19: Франция, Швеция, Португалия, США, Италия, Испания, Великобритания и др.

Помимо абсолютных величин прироста смертности интерес представляет также расчет удельного прироста (% от ожидаемой смертности в отсутствие COVID-19).  Такие оценки и их доверительные интервалы легко получить:

covid_sum <- sum(covid_y)

excess_mortality_rel <- apply(m_dist, 1, 
                              FUN = function(x){
                                monthly <- (covid_y - x) / x * 100
                                x_sum <- sum(x)
                                total <- (covid_sum - x_sum) / x_sum * 100
                                c(monthly, total)
                              }) %>% t() %>% 
  as_tibble(.name_repair = "minimal") %>% 
  setNames(., c("apr", "may", "jun", "total"))

result_rel <- excess_mortality_rel %>%
  dplyr::select(apr:total) %>% 
  pivot_longer(cols = apr:total) %>% 
  group_by(month = name) %>% 
  summarise(tibble(
    med = median(value),
    ll95 = quantile(value, p = 0.025),
    ul95 = quantile(value, p = 0.975)
  )) %>% 
  mutate(month = factor(month, ordered = TRUE,
                levels = c("apr", "may", "jun", "total"),
                labels = c("апрель", "май", "июнь", "всего"))) %>% 
  arrange(month)

result_rel
## # A tibble: 4 x 4
##   month    med  ll95  ul95
##   <ord>  <dbl> <dbl> <dbl>
## 1 апрель  9.70 -1.02  22.8
## 2 май    18.6   6.48  34.4
## 3 июнь   42.9  26.4   64.2
## 4 всего  23.1  13.2   35.2

Графическое изображение оценок удельного прироста смертности приведено на рис. 10.

Рис. 10. Точечные оценки удельного прироста смертности в Беларуси в апреле-июне 2020 г. и их 95%-ные доверительные интервалы

Таким образом, наиболее вероятный удельный прирост смертности в Беларуси в апреле, мае и июне 2020 г. составил 9.7%, 18.6% и 42.9% соответственно. В целом же за три месяца точечная оценка удельного прироста смертности оказалась равной 23.1%. При этом истинное значение показателя с вероятностью 95% находилось в пределах от 13.2% до 35.2%. Данные оценки доверительных границ ставят Беларусь по негативному воздействию COVID-19 в один ряд с такими европейскими странами, как Франция, Швеция и Италия, где удельный прирост смертности за 11 недель пандемии составил 20%, 28% и 35% соответственно.

Выполненный анализ можно улучшить

Использованные здесь данные по месячной смертности в Республике Беларусь в исходном варианте имеют несколько значений (всего их было 12), обозначенных как "Unknown", т.е. принадлежность этих значений к какому-то конкретному месяцу неизвестна. Такие значения были удалены из анализа. В большинстве случаев эти неизвестные значения не превышали 20 и их удаление вряд ли оказало сколь либо существенное влияния на полученные результаты. Однако одно значение (из 2005 г.) оказалось необычно высоким - 1446 человек, и его удаление потенциально могло сказаться на качестве построенной прогнозной модели. Если бы оказалось возможным понять природу этого "выброса" и скорректировать обучающие данные соответствующим образом, то приведенные выше оценки избыточной смертности могли бы несколько измениться.

Как было отмечено выше, оцененная в этом исследовании избыточная смертность включает все возможные причины смерти, в т.ч. не связанные с COVID-19. Более точную оценку избыточной смертности, непосредственно обусловленной короновирусной инфекцией, мог бы дать т.н. "причинно-следственный анализ" (англ. "causal impact analysis", реализован в пакете CausalImpact для R; Brodersen et al. 2015). В отличие от примененного здесь подхода, когда прогнозная модель была построена только с использованием информации о свойствах самого временного ряда, методология причинно-следственного анализа предполагает также включение в модель дополнительных предикторов, тесно коррелирующих с моделируемым откликом в период до наступления воздействия. Более того, важным условием применимости этой методологии является то, что наступление воздействия (в нашем случае это была бы вспышка короновирусной инфекции) никак не должно влиять на последующую динамику предикторов. Выполнение этого условия потенциально может дать намного более надежные прогнозные модели и, как результат, точнее оценить размер эффекта от изучаемого воздействия. К сожалению, мне пока не удалось обнаружить подходящие предикторы для выполнения подобного анализа смертности от COVID-19 в Беларуси.