關於疫情控制的經濟模型:SIR

雖然說是去年的論文了,我也以為不會寫這東西,但是最近第二波疫情似乎來勢洶洶,可以寫一下這篇的code以及變數方法(所有的code都是用Python寫成的,我的環境是Anaconda 3),提供給有興趣研究的人參考,同時也科普一下大家常在說的控制Rt到底是怎麼一回事。

這篇論文出自於Andrew Atkeson(2020)的這一篇”What Will Be the Economic Impact of COVID-19 in the US? Rough Estimates of Disease Scenarios“,後續則有許多人引用這篇發展出更深入的探討,但我們可以把重心放在這邊,而我自己則是把疫苗施打率放進去裡面當作一個變數,試圖分析第二波如何在疫苗施打率之下有效控制,目前沒看過有人寫台灣的就自己寫一個。

基本上他的概念是我們要如何透過社會隔離得到最少的染疫人口比例,而這其實真的是可以計算的,事實上去年甚至到今年許多國家都還是採用這個模型作為開端去做符合當地情況的變化,而Robert Shiller於2017在”Narrative Economics“這篇論文裡面也有用這樣的模型來說明報紙以及口接交談的概念或者謠言如何改變我們的思維,所以他不僅僅是一篇一病相關的,未來還有很多可以探討的空間,而這就讓我們開始吧。

SIR模型是什麼?

SIR,其實是疫病傳染的三種途徑,susceptible(可能受感染),infection(感染),remove(移除),但其實詳細來講中間從S->I的過程中還有一個expose(暴露感染源)的階段。

在假設目前沒有任何疫苗的情況下,我們假設總人口最終都會到達remove的階段,亦即死亡或者康復,而從S->E->I->R這過程中會取決於:

1.疫情的傳遞速率(簡稱為\beta(t))

2.病毒淺伏時間(簡稱為\sigma(t))

3.復原速度,可能是康復或者死亡(簡稱為\gamma(t))

而我們好奇的是從每個階段到下一個階段隨著時間的變化,因此可以寫成簡稱為\frac{dS}{dt},也就是隨著時間的變化S是如何改變的,而我們可以把這樣的符號寫成\dot{s},所以你如果看\dot{s},\dot{e}之類的我在講的是隨著時間的變化,S,E等變數的人口數是怎麼改變的。

那我們可以簡單得出幾條公式:

(1)   \begin{equation*} \dot{s}=-\beta(t)*s(t)*i(t)\end{equation*}

(2)   \begin{equation*}\dot{e}=\beta(t)*s(t)*i(t)-\sigma*e(t)\end{equation*}

(3)   \begin{equation*}\dot{i}=\sigma*e(t)-\gamma*i(t)\end{equation*}

看起來很複雜?其實不會,讓我們一條接著一條講述:

第一條公式講的是每日未受感染的人口數,但這些人他有可能被感染,於是健康人口數會隨著疫情的傳遞速率(\beta(t))而減少。

第二條公式講的則是是已經暴露於感染者當中的人口數每日的變化量,他當然是跟著前面的感染人口數增加增加,但我們要扣掉那些淺伏期間過去的人口,那些是確定已經受到感染的。

第三條則是每日確定受感染人口數的每日變化量,他取決於暴露於病毒的人過去淺伏期而發作的人口數卻要扣掉完全康復的。

而我們可以在額外加上一個變數r,也就是完全免疫恢復或者因為疫情死亡的人口數,因此我們可以得到r=1-s-e-i。而我們還要多定義一個變數為R(t),稱呼其為有效感染率(Effective Reproduction Rate),而\frac{\beta(t)}{\gamma}=R(t),因為真正的感染速度其實還要考慮到一個人從真正染病到完全康復的時間,而這段時間居處的人才會是真正的傳染率,而這也是你會常常在網路上面看到人家說要降低R裡面的那個R,以及很多人在說R=3或者R=1之類的那個R。而每種不同的病毒都有不同的R,但是透過社交隔離等政策我們可以有效降低這個數字。

好了,那我們有充分的知識了,現在我們來看幾種不同的情況下模擬的疫情狀況,而我要聲明的有幾點:

1.這只是一個模型,沒有模型是完美的,不要拿模型的數字去想精準預估現實。

2.這是一個簡單模型,我們沒有算進去跟家人再一起染疫的機率,也沒有算進去情突變,也沒算上每個人染病以後可能根據體質而有不同淺伏期等等

3.假設已經染病康復的不會在染病,但現實上我們看過二次感染的例子。

4.以下所有的案例只是說明不同情境會有什麼變化,詳細的數字都還是要已各國政府公告為主,是學術研究論文的產物而已,而後續也有許多人改良這模型,有興趣的可以去閱讀例如Acemoglu et al(2020)”Optimal Targeted Lockdowns in a Multi-Group SIR Model“或者是Farboodi et al(2020)”Internal and External Effects of Social Distancing in a Pandemic“,而我相信如果能結合更多醫學以及生物學的專業會是更可靠的模型,不要只聽經濟學的一面之詞。

所以R到底有啥作用?

先前那一段我們最後提到了R,那我們先寫個code模擬一下這個疫情的發展情況,首先先架設好這一次需要用到的外插軟體,以下只是簡單的環境設定而已,跟模型無關:

%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (11, 5)  #set default figure size
import numpy as np
from numpy import exp
from scipy.integrate import odeint

接下來,設定一下我們有興趣的幾個變數:

pop_size = 24000000
γ = 1 / 18
σ = 1 / 5.2
def F(x, t, R0=1.6):
    """
    Time derivative of the state vector.

        * x is the state vector (array_like)
        * t is time (scalar)
        * R0 is the effective transmission rate, defaulting to a constant

    """
    s, e, i = x

    # New exposure of susceptibles
    β = R0(t) * γ if callable(R0) else R0 * γ
    ne = β * s * i

    # Time derivatives
    ds = - ne
    de = ne - σ * e
    di = σ * e - γ * i

    return ds, de, di

首先第一行是人口數,我設定為2千四百萬人,而我們從文獻上知道COVID-19從染病到康復總共是18天,平均淺伏期是5.2天,因此把\sigma\gamma設定下來。

而下面則是定義我們剛剛提到的SIR模型,這麼方程式(F)裡面總共有主要的需要考慮的:主要變數(x,也就是我們提到的s,e,i)、時間(t)以及有效傳染率R,而我們先把R定為常數為1.6。

接下來則是根據有效傳染率回推到傳染率\beta,而接下來則是跟上述講的3條公式一模一樣。接下來設定一些一開始的變數:

# initial conditions
i_0 = 800 / pop_size
e_0 = 2400 / pop_size
s_0 = 1 - i_0 - e_0
x_0 = s_0, e_0, i_0

在這邊我們假設染病者有800人,已經受到感染但是還沒有病狀的有2400人,也就是前面兩條講的,後面則是可以得知s=1-i-e以及把s,i,e設為基本的變數向量x。

而記得我們考慮的都是隨著時間的變化而改變的不同階段人口數嗎?而我們需要一個叫做odeint的小工具幫我們計算時間內的變化量,而對於數學有點概念的這則是幫我們積分不同系統的ODE(ordinary differential equation),也就是上述的那幾條,而code是這樣寫的:

def solve_path(R0, t_vec, x_init=x_0):
    """
    Solve for i(t) and c(t) via numerical integration,
    given the time path for R0.

    """
    G = lambda x, t: F(x, t, R0)
    s_path, e_path, i_path = odeint(G, x_init, t_vec).transpose()

    c_path = 1 - s_path - e_path       
    return i_path, c_path

透過這一條可以解出不同路徑,進而達到我們想知道的比例,而在這邊我想知道的是染病人數以及總感染人數(染病人數+康復人數),也就是最後一行return出來的東西。

而上面有一個我解釋道的是t_vec,也就是我們預期模擬的時間區塊,我設定為300,也就是300天,而我把模擬的次數模擬為1000次,而下面這條code可以幫我們做到這一點:

t_length = 300
grid_size = 1000
t_vec = np.linspace(0, t_length, grid_size)

接下來我們可以進入一開始的問題:所以R有什麼用?下面的code就是把R從1.6到3區分為6個區間,讓你看看到底不同R的模擬:

R0_vals = np.linspace(1.6, 3.0, 6)
labels = [f'R0 = {r:.2f}' for r in R0_vals]
i_paths, c_paths = [], []

for r in R0_vals:
    i_path, c_path = solve_path(r, t_vec)
    i_paths.append(i_path)
    c_paths.append(c_path)

前面兩行只是不同的R(t)切割而已,後面則是要針對不同的R去作解,沒什麼特別的,最後我們設定一個幫助我們畫圖的函數就好:

def plot_paths(paths, labels, times=t_vec):

    fig, ax = plt.subplots()

    for path, label in zip(paths, labels):
        ax.plot(times, path, label=label)

    ax.legend(loc='upper left')

    plt.show()

好了,萬事俱備,接著請你輕鬆地使用下面的code把你畫圖的函數叫出來:

plot_paths(i_paths, labels)

你應該會得到圖一的結果,y軸是感染的人數佔總人口比例,x軸則是天數,你可以看到如果R=1.6的話感染人數是可以有效壓制的,而當R越來越大的時候,感染的人數高鋒會越大,而爆發也會略早導致醫院很可能負荷不了。

圖一:不同R的影響

而R怎麼控制?事實上每種病毒的R並不一樣,但是如果人類自己有效跟他人隔離的話,R是可以有效降低的。那些下來,我們來模擬一下如果有校舍交距離控制的情況下能夠如何改變R。

社交距離控制時間

我們假設一個變數\eta他代表的是不同強度的社交距離控制時間,你可想像成為數字要押多久,或者要經過多長的觀測時間來決定是否實行隔離,而對於研究經濟與封城有興趣的學者會對不同的\eta給予不同的經濟成本進而去分析最大效益,但這邊我們就單純看不同社交距離控制的影響,而我們假設\eta介於0到1之間,越接近1越接近全境封鎖的程度。

而R就會變成一個\eta以及時間(t)的函數:

(4)   \begin{equation*}R(t)=R0*e^{-\eta*t}+(1-e^{-\eta*t})*\bar{R}\end{equation*}

R0代表的是什麼也不做的情況下的疫情傳染率,而\bar{R}則是隨著時間最後可以達到的疫情傳染率,而我們用指數型的e來衡量不同時間的變化可以改變的程度。

根據上面的公式我們可以定義以下的不同情況的隔離政策的函數,我假設最初的R是3,而最後的是1.6:

def R0_mitigating(t, r0=3, η=1, r_bar=1.6):
    R0 = r0 * exp(- η * t) + (1 - exp(- η * t)) * r_bar
    return R0

接著,我們來測試不同的\eta會產生的變化,如同下面的code:

η_vals = 1/5, 1/10, 1/20, 1/50, 1/100
labels = [fr'\eta = {η:.2f}' for η in η_vals]

接著,我們可以用下面的code把它畫出來,並得到圖2的結果,可以看到不同的\eta產生不同的R隨著時間的變化量

fig, ax = plt.subplots()

for η, label in zip(η_vals, labels):
    ax.plot(t_vec, R0_mitigating(t_vec, η=η), label=label)

ax.legend()
plt.show()
圖二:不同\eta產生不同的R

而從上一章節我們應該已經知道這會對於染疫者數量有很大的差別,因此下面的code可以讓我們看清楚染疫者的數量,得到圖三的結果:

i_paths, c_paths = [], []

for η in η_vals:
    R0 = lambda t: R0_mitigating(t, η=η)
    i_path, c_path = solve_path(R0, t_vec)
    i_paths.append(i_path)
    c_paths.append(c_path)
plot_paths(i_paths, labels)
圖三:不同\eta的染疫者數量

我們從圖三可以看到,光是同0.01升級到0.02就會差距很大,但是更值得注意的是從0.05到0.1或者是0.2的差別不到,如果說考慮到每增加0.01單位的\eta經濟成本都是一樣的,那其實並非越快宣布越好,而是可以抓到一個適當的值只得經濟衝擊最小,因為即便控管變的嚴厲但是得到的效果其實沒有差別到哪裡去。

不同封城時間的效果

接下來,我們來模擬不同封城時間的效果,我們可以模擬到底封城30天跟60天的效果比較,而我們假設封城可以把R拉到剩下0.5好了,其實我們唯一要改的只有solve_path的函數而已,而具體如下:

t_length = 500
grid_size = 1000
t_vec = np.linspace(0, t_length, grid_size)
R0_paths = (lambda t: 0.5 if t < 30 else 2,
            lambda t: 0.5 if t < 60 else 2)

labels = [f'scenario {i}' for i in (1, 2)]

i_paths, c_paths = [], []

for R0 in R0_paths:
    i_path, c_path = solve_path(R0, t_vec, x_init=x_0)
    i_paths.append(i_path)
    c_paths.append(c_path)

最上面的三行code是把時間天數拉長而已,而其餘這只是在說如果封城30天內傳染力為0.5以後完全開放又變成2的情況,以及封城60天之後完全開放的情況下如此而已,沒什麼特別的,唯一值得注意的我只比較兩個極端事件例如0.5到2,當然你可以加更多條件例如循序開放之類的可能會更為現實。

而接下來我們同樣用下面的code繪圖,會呈現圖四:

plot_paths(i_paths, labels)

圖四:不同封城效果

從圖四可以看到,其實都是會到頂峰的而封城其實最大的效果是在於買時間,買注射疫苗以及使得醫院可以維持校去的時間,但這是因為我們假設只有封城跟解封兩種情況,而我們也沒有把疫苗的案例放進去,而下面,我們可以嘗試把疫苗注射人口放進去,試圖模擬如果疫苗有效施打的話,疫情會如何。

疫苗問世

原本Atkeson(2020)的論文沒有把疫苗這個變數放進去,而下面的假設是我自己做的,所以可能有很大的問題,但我是架構在原本的模型上多了一個疫苗而已,而我的假設是如果一個人施打疫苗以後他就會跳過S->E->I->R的階段,直接從S->R,而我同樣是上面的公式但是多了一條疫苗施打速度:

(5)   \begin{equation*} \dot{v}=\alpha*v(t)\end{equation*}

(6)   \begin{equation*} \dot{s}=-\beta(t)*s(t)*i(t)*(1-v(t))-v(t)\end{equation*}

(7)   \begin{equation*}\dot{e}=\beta(t)*s(t)*i(t)-\sigma*e(t)\end{equation*}

(8)   \begin{equation*}\dot{i}=\sigma*e(t)-\gamma*i(t)\end{equation*}

所以我透過\alpha這個值來模擬疫苗施打人數的變化率,而同樣的施打疫苗的人應該從可能被感染者,S裡面移除,但是沒施打的人還是會依照原本的流程走,而我可以重寫一段code如下,基本的道理都一樣,我就稍微帶過:

pop_size = 24000000
γ = 1 / 18
σ = 1 / 5.2
α=1/1000
def F(x, t, R0=1.6):
   
    s, e, i, v = x

    
    β = R0(t) * γ if callable(R0) else R0 * γ
    ne = β * s * i
    

    
    dv=α*v
    ds = - ne*(1-v)-v
    de = ne - σ * e
    di = σ * e - γ * i

    return ds, de, di, dv

上面的code只是多了\alpha作為疫苗施打速度而已,而其他則如同上述的公式一樣,沒有太大的改變。

# initial conditions
i_0 = 800 / pop_size
e_0 = 2400 / pop_size
v_0=7000 / pop_size
s_0 = 1 - i_0 - e_0-v_0
x_0 = s_0, e_0, i_0,v_0

在這邊,我設定打疫苗人口數在一開始為7000人,其他不便。

def solve_path(R0, t_vec, x_init=x_0):
    """
    Solve for i(t) and c(t) via numerical integration,
    given the time path for R0.

    """
    G = lambda x, t: F(x, t, R0)
    s_path, e_path, i_path,v_path = odeint(G, x_init, t_vec).transpose()

    c_path = 1 - s_path - e_path-v_path       # cumulative cases
    return i_path, c_path

同樣解ODE,只是增加疫苗變數而已。

R0_paths = (lambda t: 0.5 if t < 30 else 2,
            lambda t: 0.5 if t < 60 else 2)

labels = [f'scenario {i}' for i in (1, 2)]

i_paths, c_paths = [], []

for R0 in R0_paths:
    i_path, c_path = solve_path(R0, t_vec, x_init=x_0)
    i_paths.append(i_path)
    c_paths.append(c_path)
plot_paths(i_paths, labels)

上面的code就同樣是解理種不同的封城情境,而你可以看到結果如同圖五,就會看到波峰不再是一樣的,而是明顯有減少,而比較圖四你也會發現即便是30天的封城他的波峰也從原本的12%的染病人口數變成了8%。

圖五:施打疫苗以後的兩種封城情況

而如果我們不是7000人有打疫苗,而是有20000人呢?下面的code更改了基本條件,重新跑一次,可以得到圖六:

# initial conditions
i_0 = 800 / pop_size
e_0 = 2400 / pop_size
v_0=14000 / pop_size
s_0 = 1 - i_0 - e_0-v_0
x_0 = s_0, e_0, i_0,v_0
R0_paths = (lambda t: 0.5 if t < 30 else 2,
            lambda t: 0.5 if t < 60 else 2)

labels = [f'scenario {i}' for i in (1, 2)]

i_paths, c_paths = [], []

for R0 in R0_paths:
    i_path, c_path = solve_path(R0, t_vec, x_init=x_0)
    i_paths.append(i_path)
    c_paths.append(c_path)
plot_paths(i_paths, labels)

圖六:疫苗初始施打人數為20000人

比較圖六以及圖五你可以看到更明顯的差距,疫苗的施打人口基數會造成疫情二次爆發有著重大的影響,甚至可以說如果疫苗的人數增加了,那麼封城的經濟考量其實可以放寬一點。

結論

這麼模型只是一個經濟學上的模擬模型而已,他其實很需要生物或者疫病學上面的知識補充,而且他並不能代表真實的社會情況,事實上原作者Andrew Atkeson的那篇論文無論是染病人數或者是死亡人數的模擬結果都遠遠都高出真實數據,因此這裡面的圖片y軸代表的人數絕對不代表真實人數。

但這篇文章可以告訴我們什麼:

1.不同的傳染速率R的確有差,而我們可以靠著人為控制降低R的值,進而使得感染數量下降。

2.快速反應使得R降低的確能有效壓制疫情,但是並非永遠都是同樣的效果,事實上,我們應該可以抓到最大的效果,使得經濟損失控制在一定程度裡面,而這是經濟學模型可以幫忙的地方。

3.單純的封城的確可以買到時間,但真正有效的還是疫苗的問世。

4.疫苗問世以後,施打的人口基數會大幅決定二次感染的疫情影響程度,但這是假定策略一樣的情況下,我並沒有模擬疫苗問世但是又同時放鬆款管制的結果。

而我們其實可以把這個模型做延伸,如果這不是病毒而是假消息呢?有多少人口數會被假消息影響,進而從不影響->懷疑自己->深信假消息->被打臉後不相信,這樣的過程裡面做出選擇以及動態的均衡,而假消息的疫苗假設是科普教育的話,那麼科普教育的普及勢必會對假消息的散播有著重大的影響,而這也是我寫這部落格的初衷,雖然更新得比較慢,但至少我們希望散播的能夠有效對抗假消息的疫苗,進而降低感染率。

發佈留言

發佈留言必須填寫的電子郵件地址不會公開。 必填欄位標示為 *