Доверительный интервал для числа заболевших коронавирусом (расчёт по летальности)

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

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

И вот перед нами число 3. Три умерших, не вектор какой-нибудь, не таблица и не выборка. Это факт. Три умерших попали к нам совершенно не случайно. Они умерли.

Итак, рассмотрим один из самых простых методов определения количества заболевших — по коэффициенту летальности и количеству смертей. Допустим, что нам известна летальность и она равна 1%. В этой ситуации логичным и правильным будет считать, что количество выздоровевших будет равно 297. Но какова же надёжность этого суждения? Можем ли мы просто отмахнуться от того, что у нас есть трое умерших, заявив, что три — не статистика?



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

qnbinom(p=c(.025,.975),size=3, prob=0.01)

Здесь p — это 2.5% снизу и 2.5% сверху, между которыми и находится искомый диапазон.

Результат — доверительный интервал от 60 до 717. Не так уж и плохо! Вполне вероятно, что трое умерших означают совсем даже не 297 выздоровевших, а всего шестьдесят! Но может быть и за семьсот. :-(

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

Напишем функцию random_infected, которое моделирует ситуацию болезней и смертей.

random_infected <- function(deaths, fatality_rate)
{
  dead = 0
  all = 1
  while (dead < deaths) {
    if (runif(1) < fatality_rate) {
      all = all + 1
      dead = dead + 1
    } else
      all = all + 1
  }
  return(all)
}

Эта функция делает следующее — бросает «n-гранный» кубик (используя равномерное распределение). Если выпала единичка, то увеличивает число dead и число all на единицу. А если не выпала — то только число all. Каждый бросок этого кубика — это заболевший человек, который может либо умереть, либо выздороветь. Как только у нас набирается заданное параметром deaths количетсво умерших, мы останавливаемся и сообщаем, сколько раз бросали кубик (число all). Вероятность выпадения единички на нашем воображаемом кубике — это и есть летальность, в нашем случае параметр fatality_rate.

infected_sizes<-replicate(100000,random_infected(deaths=3,fatality_rate=0.01))

А теперь посчитаем это число 100 тысяч раз. У меня ноутбук старенький, так что неохота ждать, пока посчитается миллион.

После этого можно посчитать среднее арифметическое полученных чисел. У меня получилось 301.2 — очень похоже на ожидаемое число 300. Вот так выглядит распределение количеств бросаний нашего кубика смерти:

library(ggplot2)
theme_set(theme_classic())

g <- ggplot(data.frame(infected_sizes=infected_sizes), aes(infected_sizes))
g + geom_density(alpha=0.8,fill="plum")

Вот оно — отрицательное биномиальное распределение, прошу любить и жаловать. Можно относительно легко по таким данным дать примерные ответы на вопросы — какова вероятность, что общее количество заболевших меньше пятидесяти (1.2%) или больше 1000 (0.3%).

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

Напомню, что это кубик мы бросаем мгновенно. К модели расчёта по смертности, которая использовалась в нашумевшей статье Томаса Пуэйо у меня есть небольшая претензия. Там мы предполагаем, исходя из 3 умерших в день X, летальности в 1% и знания о том, что среднее время между заражением и смертью составляет 17 дней, что в день X-17 заразилось 300 человек. Однако такой расчёт справедлив только в случае, если количество заболевающих каждый день одинаковое. Так как 17 дней — это не строгое число, у него тоже есть доверительные интервалы и ошибки. Если у нас бурный рост количества заболевших, то среди умерших в день X мы имеем некоторое количество заразившихся не 17 дней назад, а 16 или 15 дней, а может и 10 дней назад. Возможно, их даже больше, чем тех, кто заразился 17 дней назад. Таким образом, в ситуации быстрого роста количества заболевших такой обратный расчёт может приводить к завышенным оценкам распространённости заболевания. В общем, всё сложно.

P.S. Спасибо Григорию Демину за подсказку про тип распределения.