評估新型冠狀病毒地區防控效果的一種近似方法

李冀鵬 洪峰 白薇 廖敬儀 張彥如 周濤?

1) (電子科技大學大數據研究中心, 成都 611731)2) (清華大學深圳國際研究生院, 深圳 518055)(2020年3月14日收到; 2020年5月15日收到修改稿)

我們觀察到地區累計確診的病例數目和武漢封城前流入的人口總數高度相關, 且本地第三代感染者占比很小. 基于此, 提出了一種考慮輸入病例和地區人口效應的定量化評估新型冠狀病毒地區防控效果的近似方法, 并將其用于評估武漢流出人口前50 的城市防控的成效. 防控效果最顯著的10個城市依次是石家莊、洛陽、恩施、周口、廈門、貴陽、咸寧、安慶、信陽、南寧.

1 引 言

新型冠狀病毒（COVID-19）在世界范圍內迅速流行，對人類健康和全球經濟發展造成了不可估量的影響[1?5]. 估計一種傳染病傳染能力最重要的參 數 是 基 本 再 生 數R0(basic reproduction number), 指在一個全是易感染態個體構成的群體中, 一個感染態的個體在恢復之前平均能感染的人數[6]:R0> 1 表示有流行病爆發的可能,R0< 1 則表示疾病難以傳播. 早期的分析顯示, COVID-19 的基本再生數所屬區間約為[2.0, 4.0], 例如2.47—2.86[1], 1.4—3.9[7], 2.39—4.13[8], 2.0—3.3[9],1.5—3.5[10], 2.88—3.67[11], 1.4—3.8[12]. Zhou等[13]

假設COVID-19 早期傳播可以用一個SEIR 動力學模型近似刻畫, 且自由傳播時呈指數趨勢, 估計基本再生數為2.8—3.9, 后根據最新的流行病學特征參數修正為2.2—3.0. 宋倩倩等[14]用指數增長法、最大似然法和SEIR 模型分別對R0進行了估計, 其估計區間分別為3.63—3.87, 2.90—3.43和3.71—4.11. 總體來說, 新型冠狀病毒肺炎屬于具有中-高等級傳播能力的傳染性疾病, 最近全球范圍內的快速爆發也驗證了這一結論.

2020年1月20日, 中華人民共和國國家衛生健康委員會發布公告, 將新型冠狀病毒肺炎納入《中華人民共和國傳染病防治法》規定的乙類傳染病, 并采取甲類傳染病的預防、控制措施. 隨后,各省市自治區出臺了有力的防控措施, COVID-19 的傳播經歷了短暫的高峰之后迅速在中國得到遏制. 定量化評估一個地區防控的效果具有非常重要的意義—既有助于我國流行病學專業人士和相關管理人員比較性地回顧總結得失, 又可以為依然飽受病毒肆虐的海外地區提供參考. 為了評價防控措施的效果, 通常采用有效再生數對傳染病在采取了防控措施(非自由傳播狀態)下的傳播能力進行實時評估. 有效再生數Rt定義為在某時刻t, 一個感染者平均而言可以感染其他個體的數量[6], 如果有效再生數下降到1 以下, 就說明該傳染病已經處于可控的狀態. Zhang 等[15]分析了湖北之外若干地區的有效再生數, 發現均在一月底之前降到了1 以下. Chen 和Zhou[16]以1月21日為全國強有力防控措施的起點, 通過計算有效再生數, 發現絕大部分省份在1月21日之后一周以內, 有效再生數就下降到1, 湖北的有效再生數也在2月2日下降到1 以下. 陳端兵等[17]進一步的分析顯示, 經歷了短暫的爆發增長后, 韓國的有效再生數也于2月25日下降到1 以下. 盡管有效再生數可以在疾病傳播過程中反映當前控制的效果, 但依然不適于作為一個地區防控效果的總體指標. 這是因為有效再生數是一個時間序列, 而我們很難直接比較兩個時間序列, 除非某序列的取值一直低于或者高于另外一個序列. 文獻[15?17]暗示可以用有效再生數降到1 的時間來衡量防控效果. 這有明顯的實用價值, 但是作為事后評估是不公平的, 因為不同地區人口密度、人口結構和生活工作習慣可能存在重大差異, 因此在傳染病自由傳播時期的基本再生數就可能存在顯著差異. 另外一種直觀而又簡單的方式, 就是直接用地區的累計確診人數或者累計確診者占地區總人口的比例(確診比例)作為防控效果的指標. 前者顯然是不合適的, 因為不同地區人口總數差異很大; 后者表面合理, 實際上也不公平,因為不同地區與武漢“親疏有別”, 輸入病例數目相差很大, 如果只看確診比例, 那么與武漢人員來往越少, 防控效果就越好.

本文擬提出一種可以量化評估一個地區防控效果的方法. 我們的基本思路是根據COVID-19 在中國傳播的特征和流行病動力學的機制, 推斷影響一個地區流行病傳播的普適規律并把它們看作不可控的因素. 利用數據分析的方法剝離這些不可控因素后, 剩余的差異部分就被認為是不同控制效果導致的, 從而可以用這些差異來評估一個地區防控的效果.

2 方 法

在介紹具體的方法和假設之前, 我們先指出兩個重要的事實. 第一, 各地區確診人數很大程度上受到了武漢封城前流入人口的影響[18?20]. 表1 列出了武漢封城前2 周流入人口前50 名的城市的流入人口數(按照流入人口數排序)、常住人口數和截止到2020年2月25日的累計確診數(后面零星確診案例可以忽略). 圖1 顯示了這50個城市武漢流入人口數和累計確診數的關聯, 可以看到, 兩者幾乎是線性關系, 相應的Pearson 關聯系數高達0.9705. 第二, 不同地區確診病例中的相當一部分是輸入病例, 例如Zhang 等[15]采集的有詳細信息的湖北外8579 例病例中, 截止到2月17日, 有42.80%在武漢或除武漢外湖北其他地區有過暴露史. 能夠獲取較完整病例數據的地區, 例如深圳、四川等地, 均有1/3—1/2 病例系輸入病例—絕大多數來源于武漢和除武漢外湖北其他地區, 少數來自于其他省、市、自治區. 如果將輸入病例視作第一代, 輸入病例在本地的感染者視作第二代, 第二代感染的病例視作第三代, 則第三代及以上的感染者是少見的.

圖1 武漢流入人口最多的50個城市的累計確診數(人)與武漢流入人口數(萬人)之間的關系. 其中數據點代表城市, 紅色直線代表線性擬合的結果(R2 = 0.942)Fig. 1. The relationship between the cumulative number of confirmed cases in the top-50 cities with maximum inflow population from Wuhan and the number of corresponding inflow people. Cities are represented by data points and the linear fitting is represented by the red line (R2 = 0.942).

基于以上事實, 我們忽略一個地區三代以及三代以上的感染者, 則一個地區最終的累計確診人數主要由兩部分組成: 一是輸入病例(包括武漢和其他地區, 但以武漢為主), 二是由輸入病例直接引起的本地傳播. 根據倉室模型, 后者應該與輸入病例和本地常住人口(易感人群數量)的乘積成正比[6,21].記地區i的累計確診人數為Ni, 可得:

其中u是待定參數, 大約是武漢封城前人口感染比例(武漢封城前出現了人口流出的峰值, 具體請參見文獻[18]),Ai是從武漢流入地區i的人口數,Pi是地區i的常住人口數. 需要計算的關鍵量是控制參數ci, 它標志了地區i防控措施的效果.ci值越小, 說明地區i的防控成效越顯著. 顯然, 這種全接觸的假設太過簡單, 例如對于千萬人口的大城市, 任何人都不能接觸到所有人口, 而只能接觸到很少一部分人[22]—如“頓巴數”所暗示的[23,24], 一個人具有緊密關系的人是非常有限的. 事實上, 在很多關于流行病傳播動力學的研究中, 每一個染病者單位時間能夠接觸的易感者人數被假設與該個體活動區域的易感人口無關, 而是一個常數[25?27].在這種情況下, 一個地區i的確診人數可以近似為

表1 武漢封城前2 周流入人口前50 的城市名稱、流入人口數、常住人口數、截止到2020年2月25日的累計確診人數、控制因子c 的取值以及利用本文方法所得到的防控效果的排名. 武漢流入人口統計的是2020年1月10日至1月22日離開武漢的人數, 數據來源于“百度遷徙”. 常住人口數取自第六次全國人口普查數據(http://www.stats.gov.cn/tjsj/pcsj/rkpc/6rp/indexch.htm)Table 1. The six columns respectively show the names of the top-50 cities with maximum inflow population from Wuhan during the two weeks before the closure, the inflow population, the permanent resident population, the cumulative number of confirmed cases as of February 25, 2020, the values of the controlling parameter c, and the rankings of the control efficacy by the present method. The inflow population accounts for the number of people who left Wuhan from January 10 to January 22, 2020, obtained from “Baidu Migration”. The number of permanent residents is taken from the sixth national census (http://www.stats.gov.cn/tjsj/pcsj/rkpc/6rp/indexch.htm).

表1 （續） 武漢封城前2 周流入人口前50 的城市名稱、流入人口數、常住人口數、截止到2020年2月25日的累計確診人數、控制因子c 的取值以及利用本文方法所得到的防控效果的排名. 武漢流入人口統計的是2020年1月10日至1月22日離開武漢的人數, 數據來源于“百度遷徙”. 常住人口數取自第六次全國人口普查數據(http://www.stats.gov.cn/tjsj/pcsj/rkpc/6rp/indexch.htm)Table 1 (continued). The six columns respectively show the names of the top-50 cities with maximum inflow population from Wuhan during the two weeks before the closure, the inflow population, the permanent resident population, the cumulative number of confirmed cases as of February 25, 2020, the values of the controlling parameter c, and the rankings of the control efficacy by the present method. The inflow population accounts for the number of people who left Wuhan from January 10 to January 22,2020, obtained from “ Baidu Migration” . The number of permanent residents is taken from the sixth national census(http://www.stats.gov.cn/tjsj/pcsj/rkpc/6rp/indexch.htm).

(1)式和(2)式中的ci都是待計算的控制參數, 但是取值顯然不同. (1)式和(2)式都是極端的簡化,盡管在人口密集的大城市中一個感染者不可能和所有易感人口接觸, 但相比小城市或者鄉鎮, 他有更多機會通過聚集居住的小區、公共交通和商場超市等場所和更多人接觸. 事實上, 在流行病仿真模型中, 人口密集的大城市中平均每個人產生的接觸也更多[28,29], 因此我們認為真實情況應該介于兩者之間, 即:

其中a是一個[0, 1]區間中的自由參數, 當a=1時方程(3)退化到方程(1), 當a= 0時方程(3)退化到方程(2). 顯然, 如果能確定u和a的取值, 就可以通過公式

計算出每個地區i的控制參數ci, 進而對不同地區防控的效果進行評價.

3 結 果

首先估計u的取值. 2020年1月26日, 湖北省人民政府新聞辦公室就新型冠狀病毒感染的肺炎疫情防控工作召開了新聞發布會, 武漢市長周先旺表示: 受春節和疫情影響, 目前有 500 多萬人離開武漢, 還有 900 多萬人留在城里. 武漢2019年統計的常住人口是908 萬人, 考慮到本文所有計算都選用常住人口數(春節大量流動人口返鄉), 且武漢常住人口數和周先旺的表述基本一致, 本文選擇武漢封城后人口為908 萬. Li 等[30]估計武漢封城前兩周, 即從1月10日到1月23日, 武漢累計新增感染人數為13118 人(95%置信區間為[2974,23435]), 則封城前后感染率大約為0.1445%, 所以可以認為u≈0.001445(所有變量以“人”為單位).根據Zhang 等[15]的報告, 截止到2月17日, 湖北省外在武漢或者除武漢外的湖北其他地區暴露過的輸入病例占比為42.80%. 因為武漢1月23日后封城, 故2月17日之后出現癥狀的武漢輸入病例可以忽略不計. 根據國家衛健委的公開數據, 2月17日湖北省外累計確診病例數為12447 人, 2月25日湖北省外累計確診病例數為12877 人. 因此在本文數據分析的時間點2月25日, 在武漢或者除武漢外湖北其他地區暴露過的輸入病例占比約為0.428 × 12447/12877 = 41.37%. 考慮到本文分析的是武漢的輸入病例數, 要從上述數字中去掉在除武漢外湖北其他地區暴露的病例, 因此我們估計這個比例應該在30%—40%之間. 如果根據從武漢流入人口最多的50個城市的數據為基礎(表1),以公式

來估計, 則u的取值應該在區間[0.001436, 0.001915].兩種方法(數據來源于兩篇不同文獻[15,30])獲得的u的取值具有相當的一致性. 考慮到估計比較粗糙,我們取兩位有效數字, 即u= 0.0014.

圖2 雙對數坐標下(Ni – uAi)/uAi 與Pi 的關聯關系. 其中數據點代表城市, 紅色直線代表線性擬合的結果(R2 =0.141), 擬合的斜率是0.234Fig. 2. The correlation between (Ni – uAi)/uAi and Pi in the log-log coordinates, where data points represent cities and the red line represents the linear fitting (R2 = 0.141).The fitting slope is 0.234.

至此, 我們認為, 對于各地區最終的確診人數而言, 輸入病例數是最主要的因素, 有很強的解釋作用; 當地的人口數(大約是人口數的1/4次方)是次要的因素, 也具有一定的解釋作用. 控制掉這兩個因素后, 各地區存在的差異, 也就是ci取值的大小, 就可以反映各地區防控的效果. 基于估計的u值和a值, 利用(4)式計算得到各地區的ci值, 并按從小到大排序, 得到的控制效果排名即為表1 的最后一列. 根據我們的計算結果, 防控效果最好的10個城市(從最好往下排序)依次是石家莊、洛陽、恩施、周口、廈門、貴陽、咸寧、安慶、信陽、南寧; 防控效果最差的10個城市(從最差的往上排序)依次是溫州、深圳、廣州、仙桃、天津、鄂州、南昌、阜陽、杭州、十堰.

以上通過假設各ci相同來估計a值的方法是很粗糙的(實際上我們又要用這個估計的a值來算出不同取值的ci), 因此以a= 0.234 為中心, 就a的不同取值對最終排序結果的影響進行了敏感性分析. 給定兩個不同的a值a1,a2, 以及對應的兩個排序X1,X2, 我們用Kendall’s Tau 值[32]來刻畫排序X1和X2之間的相近程度. 本文考慮50個城市, 則會形成1225個城市對. 依次考察這1225個城市對, 如果某對城市i,j在X1和X2中是同序的(i在X1和X2中排名都高于j, 或i在X1和X2中排名都低于j), 則加1分; 如果它們在X1和X2中是逆序的(i在X1中排名高于j但在X2中排名低于j, 或i在X1中排名低于j但在X2中排名高于j), 則減1分. Kendall‘s Tau 值就是總分除以1225. 如果Tau 值為1, 則說明兩個排序完全一致; 如果Tau 值為–1, 則說明兩個排序完全相反.從圖3 可以看出, 在以a= 0.234 為中心的一個相當大的范圍內, 這個排序是相當穩定的——所有其他排序與a= 0.234 對應排序之間的Tau 值都大于0.8. 基于此, 我們認為a= 0.234 對應的排序結果具有一定的參考價值. 如果考慮每一個輸入病例接觸的易感者數目是一個常數[25?27], 即a= 0, 則得到的城市排名和a= 0.234 對應的排名的Kendall’s Tau 值為0.827, 具有高度的一致性. 舉例而言,a= 0 對應排名前10 的城市為石家莊、恩施、洛陽、廈門、貴陽、咸寧、周口、潛江、安慶和信陽, 其中有9個城市都和a= 0.234 對應排名重合.

圖3 不 同a 值所得到的50個城市防控效果排序的Kendall’s Tau值Fig. 3. Kendall's Tau between rankings of control efficacy for different a.

4 討 論

本文的分析顯示石家莊的防控效果最佳. 我們認為取得如此出色的成績有四方面主要的原因. 一是石家莊具有發達的制藥業和紡織業, 提供了疫情攻堅戰亟需的基礎防疫物資. 石家莊有一百多家藥企, 抗生素和維生素產量全球領先, 石家莊也是全國紡織基地之一, 其中僅匯康日用品有限公司的防護口罩日產量在疫情初期就達到了80 萬只. 二是石家莊是武漢外最早開展全民摸排的城市, 在1月21日就對全體居民進行了摸排, 特別關注湖北遷入人員和出現新冠肺炎類似感染癥狀的居民. 石家莊抓住了“防止外來輸入病例, 控制傳染源”這一疫情初期防控措施的關鍵, 有效阻斷了這些潛在的傳染鏈條. 三是石家莊老百姓對疫情非常重視, 配合程度高. 在“信用石家莊”網站公布的疫情處罰公示和疫情失信名單上未出現一例個人, 均為超市和藥房. 2月25日前石家莊新聞報導中僅出現一位女士因為未佩戴口罩和小區保安爭吵. 四是石家莊特別重視重點地區(例如交通樞紐和人口密度大的居住區域)的消毒和醫廢處理, 動用了無人機等新技術手段提高消毒的覆蓋率和頻率. 與石家莊相比, 溫州之所以防疫效果不佳, 除了基礎防疫物資條件不足外, 還歸咎于政府前期重視程度和宣傳力度不足, 造成市民對防疫工作的配合度也相對較差—溫州本地新聞中有多人多起暴力反抗檢疫和拒不配合隔離的惡性事件, 以及基層干部自身不配合檢疫措施的事件. 根據“信用溫州網”公布的信息, 截止2月20日, 溫州共有254 人因不配合疫情相關政策而失信. 在2003年非典疫情階段, 石家莊是重災區, 確診患者達33 人, 而溫州僅有1 人, 因此政府和老百姓對于類似公共衛生事件的應急處置經驗以及重視程度或有不同. 從流行病學的角度看, 政府早期的重視, 特別是全面摸排(盡早發現可能的感染者), 以及通過宣傳和防控措施加強對老百姓的引導(減少人與人之間的接觸), 是快速高效控制疫情的關鍵. 希望石家莊的經驗能夠成為未來我國及其他國家和地區處理相似公共衛生事件的參考.

本文的方法也具有一些局限性. 首先, 沒有考慮三代以上的感染者, 因此這種方法并不適用于評價武漢與國外其他疫情嚴重的地區. 事實上湖北部分受疫情影響嚴重的城市, 三代及以上的感染者數目不在少數. 其次, 不管是(1)式的倉室模型還是修正后的(3)式, 都暗含了平均場的假設, 而實際的傳播是具有非常明顯的聚集性[33]. 我們注意到圖2 的關系暗示一個人能夠在城市中產生的密切接觸數與城市人口規模正相關, 大體上和城市人口之間是1/4次方的關系. 這是合理的, 因為大城市提供了更多大量人群聚集和接觸的設施, 但是這些設施的數目不是和城市人口規模成正比, 而是呈亞線性的冪律關系[34]. 盡管這個數值本身具有直觀的合理性, 但是我們估計a值的方法很粗糙, 實際上, 應該通過外在的方法確定a的取值. 使用健康碼形成的帶有時空位置信息的“亮碼/掃碼”數據來量化分析不同規模和社會經濟發展水平下城市人口的流動和接觸模式, 從而確定a值. 我們多次呼吁政府在進行數據脫敏(人員和地點匿名化)后開放這部分數據, 并已經獲得了多個城市的支持—這部分數據不僅可以用于確定a值, 還可以用于建立流行病動力學的仿真模型, 從而更好預測傳播趨勢和預估政策效果. 總體來說, 本工作是以傳播動力學機制為出發點, 嘗試定量化評估地區傳染病防控效果的有益嘗試, 盡管還存在一些局限, 但我們相信本文報導的結果對于評價城市防控效果是有一定借鑒意義的. 與此同時, 需要提醒讀者, 本文結果尚屬于學術探索的階段性成果, 不宜作為政府或者相關管理結構對地區防控情況評估的官方參考.