跳轉到主要內容
工程的博客

使用貝葉斯分層模型來推斷疾病COVID-19參數

Srijith Rajamohan博士。
通過Srijith Rajamohan博士。

2021年6月29日 在工程的博客

分享這篇文章

在以前的文章中,我們看看如何使用PyMC3 COVID-19的疾病動態建模。這篇文章是建立在這個用例並探索如何使用貝葉斯分層模型來推斷COVID-19疾病參數和帶來的好處相比,彙集或未共享模型。我們裝一個爵士模型合成數據,生成的常微分方程(ODE),以評估疾病等參數R0。然後我們展示這個框架可以應用到一個真實的數據集(例如每天感染不同國家)的數量。我們結束這個模型的局限性和輪廓的步驟改善推理過程。

我也推出了一係列的課程Coursera覆蓋這一主題的貝葉斯建模和推理,2和3的課程尤其與這篇文章有關。檢查出來的Coursera磚計算統計課程頁麵。

爵士模型

爵士模型,如我們之前所示COVID-19發布模型,包括一組三個的常微分方程(常微分方程)。在這個模型中有三個隔間:年代,我和R。

層次模型COVID博客圖片1

這裏‘S’,‘我’和‘R’指的是易感,感染和恢複部分的人口規模的“N”等

S +我+ R = N

這裏的假設是,一旦你已經從疾病中恢複過來,終生免疫是授予一個獨立的個體。這不是很多的疾病,因此,可能不是一個有效的模型。

λ的感染率,μ的速度從疾病中康複。分數的人從感染中恢複是由“f”但是對於這項工作的目的,“f”設置為1。最後一個初值問題(IVP)為我們設置的詩賦,我(0)被認為是已知的案例數量開始大流行和S(0)可以被估計為N -我(0)。這裏我們假設整個人口是敏感。我們的目標是做到以下幾點:

  • 利用貝葉斯推理估計λ和μ
  • 使用上述參數估計我任何時間(t) ' t '
  • 計算R0

池,未共享和分層模型

假設你有關於感染的數量的信息從不同的州在美國。一種方法使用這些數據來推斷COVID-19的疾病參數(如。R0)是和所有單個參數估計。這被稱為池模型。然而,這種方法的問題是,細粒度的信息可能包含在這些單個國家或團體。另一個極端是估計單個參數R0每一個國家。這種方法的結果未共享模型。然而,考慮到我們要估計的參數對應於相同的病毒,必須有一種方式來執行這個集體,這就是我們所說的層次模型。這是特別有用,當沒有足夠的信息在某些州創建精確的估計。分層模型允許我們分享來自其他州的信息使用共享“hyperprior”。讓我們更詳細地看看這個配方使用λ的例子:

池模型,我們可以從一個分布固定參數λλμ,λσ。

層次模型COVID-symbol-4

對於一個未共享模型,我們可以得出每個使用固定參數λλμ,λι。

層次模型,我們之前有一個參數化的不恒定參數來自其他發行版。這裏,我們畫λs對於每一個國家,然而,他們通過共享連接hyperprior分布如下所示。

查看課程3介紹PyMC3貝葉斯建模和推理在最近的Coursera專門化層次模型。

分層模型合成數據

實現和演示使用分層模型,我們使用的一組常微分方程生成數據定義先生模型。這些值在預設生成步伐;這裏的時間間隔是0.25。我們也選擇兩組為了便於說明,然而,一個可以有任意多個組。λ和μ的值設置為[4.0,3.0]和[1.0,2.0]分別對兩組。這個代碼來生成以及由此產生的時間序列曲線如下所示。

生成合成數據

def先生(y, t, p):ds=- - - - - -p (0]*y (0]*y (1]迪=p (0]*y (0]*y (1]- - - - - -p (1]*y (1]返回(ds, di)次=np.arange (0,5,0.25)cases_obs=(0]*2lamμ=4.0,1.0y=odeint(先生,t=次,y0=(0.99,0.01),參數=((lamμ)),rtol=1 e-8)社會毒瘤=np.random.lognormal(平均=np。日誌(y [1::]),σ=(0.1,0.1])cases_obs [0]=:社會毒瘤,1]plt.plot(倍[1::],社會毒瘤,標記=“o”,線型=“沒有”)plt。情節(倍,y (:,0),顏色=“C0”α=0.5、標簽=f' S (t)美元)plt。情節(倍,y (:,1),顏色=“C1”α=0.5、標簽=f“我(t)美元”)plt.legend ()plt.show ()lamμ=3.0,2.0y=odeint(先生,t=次,y0=(0.99,0.01),參數=((lamμ)),rtol=1 e-8)社會毒瘤=np.random.lognormal(平均=np。日誌(y [1::]),σ=(0.1,0.1])cases_obs [1]=:社會毒瘤,1]plt.plot(倍[1::],社會毒瘤,標記=“o”,線型=“沒有”)plt。情節(倍,y (:,0),顏色=“C0”α=0.5、標簽=f' S (t)美元)plt。情節(倍,y (:,1),顏色=“C1”α=0.5、標簽=f“我(t)美元”)plt.legend ()plt.show ()

實現和演示使用分層模型,我們使用的一組常微分方程生成數據定義先生模型。

使用分層模型進行推理

defSIR_sunode(t y p):返回{“年代”:- p。林* y。* y.I,“我”:p。林* y。* y。I - p。μ* y.I,}
              sample_period = covid_data.sample_periodcases_obs = covid_data.cases_obstime_range = np.arange (0,len(covid_data.cases_obs [0)))* covid_data。sample_period I0 = covid_data.data [0]#數據是按比例縮小的S0 =1——來自S_init = S0I_init =錢數cases_obs_scaled = covid_data.data與pm.Model ()作為model4:# - - - - - - - - - - - - - - - - - - -設置先驗和hyperpriors - - - - - - - - - - - - - - - #prior_lam = pm.Lognormal (“prior_lam”,0.75,2)prior_mu = pm.Lognormal (“prior_mu”,0.75,2)prior_lam_std = pm.HalfNormal (“prior_lam_std”,1.0)prior_mu_std = pm.HalfNormal (“prior_mu_std”,1.0)林= pm.Lognormal (“λ”=,prior_lam_std prior_lam形狀2)# 1.5,1.5μ= pm.Lognormal (“畝”=,prior_mu_std prior_mu形狀2)# 1.5,1.5# - - - - - - - - - - - - - - - - - - - - - - ODE模型- - - - - - - - - - - - - - - #
              res, _,問題解決者,_,_ = sunode.wrappers.as_theano.solve_ivp (y0 = {“年代”:(S_init, (2)),“我”:(I_init, (2)},params = {“林”(lam (2)),“畝”(μ,(2)),“_dummy”:(np.array (1。),())},rhs = SIR_sunode,#的時間點,我們要訪問的解決方案tvals = time_range [1:,t0 = time_range [0),)# - - - - - - - - - - - - - - - - - - -設置可能的觀測數據- - - - - - - - - - - - - - - #我= pm.Normal (“我”μ= res (“我”),σ=0.01觀察= cases_obs_scaled [1:])R0 = pm.Deterministic (“R0”林/畝)# - - - - - - - - - - - - - - - - - - -樣本分布- - - - - - - - - - - - - - - ##如果你增加方差和分布看上去波濤洶湧,增加更有效地優化樣本容量的樣本空間#,增加樣品的總數跟蹤= pm.sample (8000年調=2000年鏈=4核=4)data = az.from_pymc3(跟蹤=跟蹤)az.plot_posterior(數據,point_estimate =“模式”round_to =2)az.summary(跟蹤)traceplot(數據)

變量的後驗分布密度最高的間隔(HDI)

Traceplots和密度圖的變量

層次模型COVID博客圖片13

現實生活中的COVID-19數據

這裏使用的數據獲得的約翰霍普金斯CSSE Github頁麵案件數量的定期更新。這裏我們情節和使用infections-per-day病例數的兩個國家,美國和巴西。然而,並沒有限製的選擇或一些國家可以使用層次模型。的情況下,下麵是3月1日,2020年1月1日2021年。圖似乎遵循類似的軌跡,雖然上的鱗片軸這些國家是不同的。考慮到這些病例都來自COVID-19病毒,這是合理的。然而,有差異需要考慮,如不同的變體,不同的地理結構和社會距離規則,醫療基礎設施等等。

情節COVID-19病例數的兩個國家

推斷的參數

采樣後驗分布如下所示,連同他們的94%最高密度區間(HDI)。

采樣後驗分布以及它們的94%最高密度區間(HDI)。

我們也可以檢查traceplots收斂,顯示了良好的混合所有的變量——一個好跡象,取樣器探索的空間。有很好的協議所有的痕跡。這種行為可以被證實相當狹窄的人類發展指數在上麵的情節。

traceplots和密度R0和其它變量

下表總結了各種推斷變量和參數的分布,隨著采樣器統計數據。而估計的變量是至關重要的,這張桌子是特別有用的通知我們取樣器的質量和效率。例如,小紅帽都等於1,說明良好之間的協議所有的連鎖店。有效樣本大小是另一個關鍵指標。如果這是小樣本的總數相比,這就是一個明確信號與采樣器的麻煩。即使小紅帽值看起來不錯,一定要檢查的有效樣本大小!

表推斷變量的分布和抽樣統計

雖然這取得了令人滿意的對我們的參數估計,我們經常碰到的問題取樣器沒有有效執行。在本係列的下一篇文章中,我們將看一些診斷問題,提高建模過程。難度遞增的順序,這些列出如下:

  1. 增加調整大小和樣本的數量。
  2. 減少target_accept參數采樣器,以減少之間的自相關樣本。利用自相關圖證實了這一點。
  3. 添加更多的樣本觀測數據,即增加采樣頻率。
  4. 使用更好的先驗和hyperpriors參數。
  5. 使用一個替代模型的參數化。
  6. 改變社會距離等措施納入模型。

您可以了解更多關於這些主題在我Coursera專業化,包含以下課程:

  1. 介紹貝葉斯統計
  2. 貝葉斯推理與密度
  3. 介紹PyMC3貝葉斯建模和推理

看到課程清單

免費試著磚
開始
看到所有工程的博客的帖子
Baidu
map