Байесовская модель точки изменения для оценки даты, когда количество новых случаев COVID-19 начнет сглаживаться в разных странах.

Проблема

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

Прежде чем мы перейдем к этому, важно напомнить, что все модели неверны, но некоторые полезны. Эта модель полезна для оценки даты изменения, а не для прогнозирования того, что произойдет с COVID-19. Его не следует принимать за удивительную эпидемиологическую модель, которая сообщит нам, когда закончится карантин, а за способ описания того, что мы уже наблюдали, с помощью распределений вероятностей.

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

Модель

Мы хотим описать y, журнал количества новых случаев COVID-19 в конкретной стране каждый день как функцию от t, количества дней с момента заражения вирусом. началось в этой стране. Мы сделаем это с помощью модели сегментированной регрессии. Точка, в которой мы сегментируем, будет определяться параметром обучения τ, как описано ниже:

Другими словами, y будет моделироваться как w ₁ + b ₁ для дней до дня τ. После этого он будет смоделирован как w ₂ + b ₂.

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

Данные

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

Вирус также запускается в разное время в разных странах. Поскольку у нас есть регрессионная модель, нецелесообразно включать данные до того, как вирус оказался в конкретной стране. Эта дата выбирается вручную для каждой страны в зависимости от развития новых случаев и никогда не является датой регистрации первого пациента. Дату «начала» лучше интерпретировать как дату начала последовательного роста вируса, в отличие от даты регистрации пациента 0.

Предварительная спецификация

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

Начиная с w ₁ и w ₂, эти параметры можно условно интерпретировать как скорость роста вируса до и после изменения даты. Мы знаем, что вначале рост будет положительным и вряд ли будет больше 1 (в конце концов, w ₁ - это градиент). С этими предположениями, w ₁ ~ N (0,5, 0,25) является подходящим предварительным решением.
Мы будем использовать аналогичную логику для p (w ₂ ), но нужно помнить о гибкости. Без достаточно гибкого предварительного решения модель не будет работать в тех случаях, когда в данных нет реальной точки изменения. В этих случаях w ₂ ≈ w ₁, и мы увидим и пример этого в разделе результатов. На данный момент мы хотим, чтобы p (w ₂ ) был симметричным относительно 0, с большинством значений, лежащих между (-0,5, 0,5). Мы будем использовать w ₂ ~ N (0, 0,25).

Далее следуют условия смещения, b ₁ и b ₂. Приоры по этим параметрам особенно чувствительны к загородным характеристикам. Страны, которые более подвержены COVID-19 (по какой-либо причине), будут иметь больше подтвержденных случаев на пике, чем страны, которые менее подвержены. Это напрямую повлияет на апостериорное распределение для b ₂ (что является смещением для второй регрессии). Чтобы автоматически адаптировать этот параметр к разным странам, мы используем среднее значение первого и четвертого квартилей y в качестве предшествующих средних значений b ₁ и b ₂ соответственно. Стандартное отклонение для b ₁ принимается равным 1, что делает p (b ₁ ) относительно ровным априорным значением. Стандартное отклонение p (b ₂ ) принимается как четверть его предыдущего среднего значения, так что предыдущие масштабируются с более высокими средними значениями.

Что касается τ, поскольку в настоящее время у нас нет доступа ко всем данным (вирус продолжается), мы не можем иметь полностью фиксированное априорное значение и рассчитывать его с помощью модели. Вместо этого предполагается, что изменение, скорее всего, произойдет во второй половине текущего диапазона дат, поэтому мы используем τ ~ Beta (4, 3).

Ниже представлена ​​модель, написанная на Pyro.

class COVID_change(PyroModule):
    def __init__(self, in_features, out_features, b1_mu, b2_mu):
        super().__init__()
        self.linear1 = PyroModule[nn.Linear](in_features, out_features, bias = False)
        self.linear1.weight = PyroSample(dist.Normal(0.5, 0.25).expand([1, 1]).to_event(1))
        self.linear1.bias = PyroSample(dist.Normal(b1_mu, 1.))
        self.linear2 = PyroModule[nn.Linear](in_features, out_features, bias = False)
        self.linear2.weight = PyroSample(dist.Normal(0., 0.25).expand([1, 1])) #.to_event(1))
        self.linear2.bias = PyroSample(dist.Normal(b2_mu, b2_mu/4))
   def forward(self, x, y=None):
        tau = pyro.sample(“tau”, dist.Beta(4, 3))
        sigma = pyro.sample(“sigma”, dist.Uniform(0., 3.))
        # fit lm’s to data based on tau
        sep = int(np.ceil(tau.detach().numpy() * len(x)))
        mean1 = self.linear1(x[:sep]).squeeze(-1)
        mean2 = self.linear2(x[sep:]).squeeze(-1)
        mean = torch.cat((mean1, mean2))
        obs = pyro.sample(“obs”, dist.Normal(mean, sigma), obs=y)
        return mean

Гамильтониан Монте-Карло используется для апостериорной выборки. Код для этого показан ниже.

model = COVID_change(1, 1,
                     b1_mu = bias_1_mean,
                     b2_mu = bias_2_mean)
num_samples = 800
# mcmc
nuts_kernel = NUTS(model)
mcmc = MCMC(nuts_kernel,
            num_samples=num_samples,
            warmup_steps = 100,
            num_chains = 4)
mcmc.run(x_data, y_data)
samples = mcmc.get_samples()

Полученные результаты

Поскольку я живу в Канаде и знаю даты, когда были приняты меры предосторожности, моделирование начнется здесь. Мы будем использовать 27 февраля в качестве даты «запуска» вируса.

Приоры:

Последующее распространение:

Начиная с апостериорных значений для w ₁ и w ₂, если бы не было изменений в данных, мы ожидали бы увидеть эти два распределения близко друг к другу, поскольку они управляют ростом скорость вируса. Это хороший знак, что эти распределения вместе с апостериорными для b ₁ и b ₂ не пересекаются. Это свидетельство того, что точка изменения, оцененная нашей моделью, верна.

Предполагаемая точка изменения: 2020–03–28.

В качестве примечания (предупреждение: без науки) моя компания 16 марта выпустила обязательную политику работы из дома. Примерно в этот день у большинства компаний в Торонто возникнут проблемы с обязательной работой из дома, где это применимо. Если предположить, что инкубационный период вируса составляет до 14 дней, это предполагаемое изменение даты имеет смысл, поскольку это 12 дней после начала широкомасштабных мер социального дистанцирования!

Соответствие модели и 90% вероятных интервалов можно увидеть на графике ниже. Слева находится журнал количества ежедневных обращений, которое мы использовали для соответствия модели, а справа - истинное количество ежедневных обращений. Очень сложно визуально определить точку изменения, просто посмотрев на количество ежедневных случаев, и еще сложнее, посмотрев на общее количество подтвержденных случаев.

Оценка конвергенции

При проведении этих экспериментов наиболее важным шагом является диагностика MCMC на предмет сходимости. Я использую 3 способа оценки сходимости для этой модели, наблюдая за перемешиванием и стационарностью цепочек и R_hat. R_hat - коэффициент, на который каждое апостериорное распределение будет уменьшаться по мере того, как количество выборок стремится к бесконечности. Идеальное значение R_hat равно 1, но значения меньше 1,1 являются сильным признаком сходимости. Мы наблюдаем смешивание и стационарность цепей Маркова, чтобы знать, производит ли HMC соответствующие апостериорные образцы.

Ниже приведены графики трассировки для каждого параметра. Каждая цепочка неподвижна и хорошо перемешивается. Кроме того, все значения R_hat меньше 1,1.

После конвергенции последнее, что нужно проверить, прежде чем переходить к другим примерам, - насколько модель подходит для данных. Соответствует ли это предположениям, сделанным ранее? Чтобы проверить это, мы будем использовать график остатков и график QQ, как показано ниже.
Я обрисовал предполагаемую точку изменения, чтобы сравнить остатки до и после изменения для проверки гомоскедастичности. Остатки соответствуют нормальному распределению с нулевым средним и не зависят от времени до и после даты изменения.

А как насчет без изменений?

Ага! Что, если количество новых случаев еще не начало сокращаться ?! Чтобы проверить надежность модели в этом случае, мы рассмотрим данные по Канаде до 28 марта. Это день, когда в Канаде началось сглаживание кривой, рассчитанное по модели.

Тот факт, что нет истинной даты изменения, не означает, что модель выдаст «Без изменений» (смеется). Нам придется использовать апостериорные распределения, чтобы обосновать, что дата изменения, указанная в модели, неуместна, и, следовательно, в данных нет изменений.

Приоры:

Последующее распространение:

Задние элементы для w ₁ и w ₂ имеют значительное перекрытие, что указывает на то, что скорость роста вируса существенно не изменилась. Постеры для b ₁ и b ₂ также перекрываются. Они показывают, что модель изо всех сил пытается оценить разумное значение τ, что является для нас хорошим подтверждением того, что априорные значения не слишком сильны.

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

Как и в предыдущем примере, MCMC сошлась. Графики следа ниже показывают достаточное перемешивание и стационарность цепочек, и большинство значений R_hat меньше 1,1.

Следующие шаги и открытые вопросы

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

Следите за сообщениями по теме в моем личном блоге!

Спасибо, что прочитали, и обязательно свяжитесь со мной по электронной почте или другим способом, если у вас есть предложения или рекомендации, или даже просто поболтать!