SIR模型的贝叶斯估计

介绍

千叶大学/Nospare米仓是。这一次，我将解释传染病数学模型中经常使用的SIR模型，并介绍它的贝叶斯估计。

什么是 SIR 模型？

我将解释基于 Nishiura (2021) 的 SIR 模型。在这里，我们不考虑个体的异质性，并隐含地假设每个人都同样容易受到感染和康复。

SIR模型是一种称为区室模型的模型，是通过将传染病流行的人群（人口）分类为几个区室（区间）进行分析的模型。

具体来说，SIR 模型将总体分为

1. S。易感人群（易感人群）→未感染且未获得免疫力的人群。
2. 胃nfected → 感染传染病的人群。
3. R。ecovered（康复人群）→已获得免疫力的人群或已隔离的人群。

   分为三个状态。对此，我们做出以下假设。

   1. 感染后立即从 S 移动到 I。
   2. 只有 I 会在人群中引起继发感染。
   3. 感染持续时间呈指数分布。
   4. 感染期后我会转换为 R。
   5. 恢复后，获得免疫力并且不会发生再感染（不会从 R 转变为 S 或 I）。
      在这些假设下，我们将 S、I 和 R 的动态描述如下。

      S和I之间的接触引起了新的感染。在这种情况下，$\beta$是单位时间内S和I之间的接触率。换句话说，令$\beta$SI 是人口从S 移动到I 的频率，与每单位时间$\beta$SI 成正比。

      另一方面，假设总体从 I 移动到 R 与 $\gamma$ 成比例。换句话说，感染者每单位时间的恢复与$\gamma$成正比。

      综上所述，以下常微分方程
      $$\frac{dS(t)}{dt}=-\betaS(t)I(t)\
      $$
      $$\frac{dI(t)}{dt}=\beta S(t)I(t)-\gamma I(t)
      $$
      $$\frac{dR(t)}{dt}= \gammaI(t)
      $$
      ，称为 SIR 模型。如果我们取三个表达式的总和，
      $$
      \frac{d}{dt}(I(t)+S(t)+R(t))=0
      $$
      注意 。

      SIR 模型行为和繁殖数

      SIR 模型的行为因模型参数 $(\beta,\gamma)$ 的值而异。通过变换 $I(t)$ 的公式，

      \begin{align}
      \frac{dI(t)}{dt}&=\beta S(t)I(t)-\gamma I(t)\\
      &= I(t)(\beta S(t)-\gamma)
      \end{align}
      

      变成。根据定义，$I(t)$ 是非负的，所以 $(\beta S(t)-\gamma)$ 是正的，即 $$\frac{\beta S(t)}{\gamma}>1 $ 如果$，$I(t)$ 增加，这意味着感染人数增加。

      如果 N=1 是随时间变化的恒定标准化总体（$N=S(t)+I(t)+R(t)$），则 $S(0)=N= 当 $t=0$ 1$ 时，也就是说，在感染的早期，我们可以假设每个人都是易感的，所以 $\frac{\beta S(t)}{\gamma}$ 变成了 $\frac{\beta}{\gamma}$，这是我们称之为“基本再生数，$R_0$”。

      如果$R_{0}$为1，则1个感染者生育1个感染者，如果为2，则1个感染者生育2个感染者。t)$收敛到0，即感染收敛。反之，如果$R_0\geq 1$，则感染人数会增加。

      $R_{0}$ 是“没有其他感染者并且没有疫苗接种或封锁等政策干预”时的假设值。另一方面，对策时的再生数称为“执行再生数，$R_t$”，也就是前面提到的$\frac{\beta S(t)}{\gamma}$。

      通过这种方式通过SIR模型监测$R_0$和$R_t$，可以在一定程度上了解传染病的流行情况。

      SIR模型参数的贝叶斯估计

      一、R库“暴发”中1978年发生在英国一所寄宿学校的H1N1流感大流行的数据，流感_英格兰_1978_学校用来。

      该数据显示，从 1 月 22 日到 2 月 4 日，763 名学生中有 512 人被感染。内容为日期、当天卧床休息的学生人数、当天康复的学生人数。下图在横轴上显示日期。我经常在坐标轴上绘制当天在床上休息的学生人数。
      

      本节介绍如何从数据中统计估计 SIR 参数 $(\beta,\gamma)$。由于进行了统计估计，因此需要在某处添加概率项进行建模，但这次我们将解释如何基于 Grinsztajn et al. (2021) 进行贝叶斯估计。

      目的是从后验分布 $p(\theta\mid y)\propto p(y\mid \theta)p(\theta)$ 中采样，观察值为 $y$，参数为 $\theta$ .其中 $p(\theta)$ 是先验分布，$p(y\mid \theta)$ 是统计模型（似然函数）。

      如果适当地确定参数$S(t)、I(t)和R(t)$的初始值，则作为常微分方程的SIR可以在数值上求解为在这里，我们将由此获得的$I(t)$ 表示为$I(t)^M$。这里，这个$I(t)^M$和实际观察到的$I(t)$，$I(t)^O$=当天卧床休息的学生人数，在一个统计模型$p(我们将使用 y\mid \theta)$ 链接。

      特别是，这次我们采用以下负二项分布作为$ p(y\mid \theta)$。 $$I(t)^O\simNegBin(I(t)^M,\phi), $$
      其中$\phi$ 是负二项分布的色散参数。

      所以参数$\theta$变成$\theta=(\beta,\gamma,\phi)$。作为先验分布，$\beta$ 具有均值为 2，方差为 1 的正限制正态分布，$\gamma$ 具有均值为 0.4，方差为 0.5 的正限制正态分布，$\phi 假设指数分布为参数 5 为 $ 的倒数。

      使用估计结果，我们还可以估计 $R_0$ 的后验分布 $p(R_0\mid y)$。

      由斯坦实现

      首先，SIR模型，即常方程部分，在Stan中实现。一个普通方程通常可以使用函数 $f$ 表示为 $$\frac{dy(t)}{dt}=f(y(t))$$。对于 SIR 模型，$y(t)=(S(t),I(t),R(t))$。在 Stan 中，首先在功能块中

      real[] f(real time, real[] state, real[] theta, real[] x_r, int[] x_i)
      

      并且一般规定了常方程的内容。其中 f 是函数名，time 是时间，state 是 y 的每个分量，theta 是计算 f 时使用的参数，x_r 是计算 f 的实数值数据，x_i 是计算 f 数据的整数值。对于 SIR 模型，

      functions {
        real[] sir(real t, real[] y, real[] theta, 
                   real[] x_r, int[] x_i) {
      
            real S = y[1];
            real I = y[2];
            real R = y[3];
            real N = x_i[1];
            
            real beta = theta[1];
            real gamma = theta[2];
            
            real dS_dt = -beta * I * S / N;
            real dI_dt =  beta * I * S / N - gamma * I;
            real dR_dt =  gamma * I;
            
            return {dS_dt, dI_dt, dR_dt};
        }
      }
      

      可以描述 SIR 模型。模型写好后，SIR 模型就可以四阶龙格-库塔法使用 进行数值求解。 SIR 在 Stan

       y = integrate_ode_rk45(sir, y0, t0, ts, theta, x_r, x_i);
      

      可以使用四阶龙格-库塔法求解 SIR 模型。这里，sir 是求解常方程 sir 的规范，y0 是初始值，t0 是离散化时间的初始值，ts 是评估解的时间（普通方程离散化的精细程度或），其他变量同上。

      接下来，在数据部分描述相关的数据。

       data {
        int<lower=1> n_days;
        real y0[3];
        real t0;
        real ts[n_days];
        int N;
        int cases[n_days];
      }
      

      接下来，用转换后的数据转换数据以加快估计速度。请注意，此 SIR 估计中未使用 x_r，x_i 是总体大小 N，依此类推。

       real x_r[0];
        int x_i[1] = { N };
      }
      

      接下来，在参数块中声明参数 $\theta=(\beta,\gamma,\phi)$。

       parameters {
        real<lower=0> gamma;
        real<lower=0> beta;
        real<lower=0> phi_inv;
      }
      

      请注意，这里我们给出了 $\frac{1}{\phi}$ 而不是 $\phi$ 的先验，即过度分散的参数。
      在转换后的参数块中，按如下方式转换参数：

       transformed parameters{
        real y[n_days, 3];
        real phi = 1. / phi_inv;
        {
          real theta[2];
          theta[1] = beta;
          theta[2] = gamma;
      
          y = integrate_ode_rk45(sir, y0, t0, ts, theta, x_r, x_i);
        }
      }
      }
      

      像以前一样在模型块中指定先验和统计模型。请注意，我们为参数块中的每个参数指定了实数。

       model {
        #事前分布
        beta ~ normal(2, 1); 
        gamma ~ normal(0.4, 0.5); 
        phi_inv ~ exponential(5);
        
        #統計モデル
        #col(matrix x, int n) - 行列xの第nベクトル．ここでは2列目にIに関するベクトルがあるとしている．
        cases ~ neg_binomial_2(col(to_matrix(y), 2), phi);
      }
      

      最后，编写诸如 $R_0$、平均愈合时间 $\frac{1}{\gamma}$ 以及生成的数量块中感染人数的估计等计算。

       generated quantities {
        real R0 = beta / gamma;
        real recovery_time = 1 / gamma;
        real pred_cases[n_days];
        pred_cases = neg_binomial_2_rng(col(to_matrix(y), 2) + 1e-5, phi);
      }
      }
      

      将这些放在一起完成以下 Stan 代码（另存为 sir_negbin.stan）。

      functions {
        real[] sir(real t, real[] y, real[] theta, 
                   real[] x_r, int[] x_i) {
      
            real S = y[1];
            real I = y[2];
            real R = y[3];
            real N = x_i[1];
            
            real beta = theta[1];
            real gamma = theta[2];
            
            real dS_dt = -beta * I * S / N;
            real dI_dt =  beta * I * S / N - gamma * I;
            real dR_dt =  gamma * I;
            
            return {dS_dt, dI_dt, dR_dt};
        }
      }
      data {
        int<lower=1> n_days;
        real y0[3];
        real t0;
        real ts[n_days];
        int N;
        int cases[n_days];
      }
      transformed data {
        real x_r[0];
        int x_i[1] = { N };
      }
      parameters {
        real<lower=0> gamma;
        real<lower=0> beta;
        real<lower=0> phi_inv;
      }
      transformed parameters{
        real y[n_days, 3];
        real phi = 1. / phi_inv;
        {
          real theta[2];
          theta[1] = beta;
          theta[2] = gamma;
      
          y = integrate_ode_rk45(sir, y0, t0, ts, theta, x_r, x_i);
        }
      }
      model {
        beta ~ normal(2, 1);
        gamma ~ normal(0.4, 0.5);
        phi_inv ~ exponential(5);
        
        cases ~ neg_binomial_2(col(to_matrix(y), 2), phi);
      }
      generated quantities {
        real R0 = beta / gamma;
        real recovery_time = 1 / gamma;
        real pred_cases[n_days];
        pred_cases = neg_binomial_2_rng(col(to_matrix(y), 2), phi);}
      

      估计结果

      上面的 Stan 文件在 R 中执行如下。

      library(rstan)
      library(gridExtra)
      
      cases = influenza_england_1978_school$in_bed  # 静養中の学生の数(感染者数の観測値)
      
      # 学生の数（人口）
      N = 763;
      
      # 時間
      n_days = length(cases) #全期間
      t = seq(0, n_days, by = 1)　#タイムステップ
      t0 = 0 
      t = t[-1]
      
      #初期値の設定
      i0 = 1　#初期に一人感染者が人口内に出現
      s0 = N - i0 
      r0 = 0
      y0 = c(S = s0, I = i0, R = r0)
      
      #Stanにデータを渡す 
      data_sir = list(n_days = n_days, y0 = y0, t0 = t0, ts = t, N = N, cases = cases)
      
      # MCMCのイタレーションの数
      niter = 10000
      
      #コンパイル(hogehogeはスタンファイルがあるディレクトリ)
      model = stan_model("hogehoge/sir_negbin.stan")
      
      #4回の独立のチェインを求める．バーンインはそれぞれ5000
      fit_sir_negbin = sampling(model,
                                 data = data_sir,
                                 iter = niter,
                                 chains = 4, 
                                 seed = 0)
      
      pars=c('beta', 'gamma', "R0", "recovery_time")
      
      #結果のプリント
      print(fit_sir_negbin, pars = pars)
      
      #推定結果の事後密度
      stan_dens(fit_sir_negbin, pars = pars, separate_chains = TRUE)
      

      下面是后验分布等等。
      四个链的边缘后验分布具有相似的形状，可以看出 Rhat 收敛。

      

      

      其他高级主题

      本文介绍了最简单的 SIR 模型。还提出了在 R 之后返回 S 的 SIRS 模型，以及考虑处于潜伏期（E = Exposed）的患者的 SEIR 模型。了解更多信息这个网站请看。

      另一个类比Dureau 等人，（2012 年）在参数随时间变化的假设下估计随机微分方程类型的 SIR 模型。

      综上所述

      Nospare Co., Ltd. 拥有专门从事各种统计领域的研究人员。用于业务数据的统计咨询和分析诺斯派有限公司请联系我们。
      
      
      

      参考

      ・西浦“传染病流行病学数据分析导论”（2021年）
      ・Grinsztajn et al. (2021) “Stan 疾病传播建模的贝叶斯工作流程”

原创声明：本文系作者授权爱码网发表，未经许可，不得转载;

原文地址：https://www.likecs.com/show-308624550.html