在前面我們已經介紹了基本再生數（R0）的原理和計算，詳情請點擊以下連結：
《COVID-19專題：關於R0，你想知道的都在這裡》
由於R0通常只用於估計傳染病早期的傳播力，為了遏制傳染病的流行，制定更為合理的防控措施，需要有一個指標能夠實時反映傳染病在人群中的傳播力，這個指標通常使用實時再生數（Time-dependent reproduction number, Rt）。
Part 1
什麼是Rt？
Rt是指t時刻平均一個發病的病例能夠導致的新發病例數。t表示一個時間單位，可以是小時、天、周等，最常使用的是天，也就是計算每一天的再生數。若Rt小於1，表示維持當前防控措施可以逐漸控制傳染病流行；若Rt大於1，表示傳染病將持續擴散，提示防控措施需要優化和加強。
Part 2
怎麼計算Rt？
Rt的概念發展了很多年，計算的方法在不斷改進，也在不同的統計軟體開發了多種算法。下面，我們利用R語言軟體中的「EpiEstim」包，基於貝葉斯框架，對Rt進行計算。
01
數據準備
數據的結構根據不同的情況有所差異。若是第一個時間點之後的病例均來自本地傳播（如COVID-19前期在武漢的流行），則僅需要2列數據，一列為發病（症狀出現）時間，一列為每天的病例數。但是，多數地方可能存在持續的輸入病例、由輸入病例導致本地病例的情況，這種情景至少需要3列數據，一列為發病（症狀出現）時間，一列為每天的輸入病例數，一列為每天的本地病例數。如表1所示，date表示發病日期，imported表示輸入病例數，local表示本地病例數。
表1 數據結構示例
02
計算Rt
計算Rt還需要一個特別重要的數據，就是代際時間分布（詳見「COVID-19專題：關於R0，你想知道的都在這裡」）。我們根據有無個案代際時間，將Rt的計算分為2種情況。前者可利用部分或全部病例的代際時間獲得分布；後者根據文獻已有的參數形成分布。
1
有個案代際時間的Rt計算
通過流行病學調查等方法，可以獲得病例的傳播關係，獲得案例之間的代際時間（如果只有部分關係，可理解為從整體中抽樣）。在「EpiEstim」包中，可以利用馬爾科夫鏈蒙特卡洛方法(Markov Chain Monte Carlo，MCMC)取得對這些代際時間進行重採樣。計算代際時間前需要將數據資料整理成類似表2的數據框。表2中每一行數據代表一個由感染者和被感染者組成對子。EL表示感染者發病時刻（整數型）的下邊界，ER表示感染者發病時刻的上邊界（ER≥EL，感染者在ER至EL期間發病）；SL表示被感染者發病時刻（整數型）的下邊界，SR表示被感染者發病時刻的上邊界（SR≥SL，被感染者在SR至SL期間發病）；type包括0,1,2，分別表示雙側不確定（感染者和被感染者的確切發病時間均未知）、單側不確定（感染者和被感染者只有一方有確切發病時間）和確切觀測（感染者和被感染者均有確切發病時間）。
表2 個案代際時間數據
由於公開數據中可獲得的COVID-19的代際時間對子較少，我們模擬了一些數據進行分析。
首先，基於MCMC對個案代際時間進行重抽樣。該過程依賴於利用dic.fit.mcmc函數，dat參數指定數據(si_data表示代際時間數據，結構如表2)，dist參數指定初始分布為gamma分布，init.pars為初始設置，burnin為去除初始抽樣的次數（初始抽樣不穩定），n.samples為重抽樣次數，seed為種子數。
code
library(EpiEstim)
library(coarseDataTools)
MCMC_seed
#使用MCMC進行重抽樣#
SI.fit
init.pars = init_mcmc_params(si_data, "G"), burnin = 1000,
n.samples = 5000, seed = MCMC_seed)
然後，將重抽樣結果通過coarse2estim函數轉化為所需要的格式，再放入estimate_R函數中進行計算Rt。estimate_R函數中，incid表示病例數據；method為代際時間的確定方法；si_sample為MCMC重抽樣之後的數據；config為貝葉斯框架的參數設置，n2表示後驗分布的採樣次數，seed表示種子數。
code
si_sample
overall_seed
R_si_from_sample
method = "si_from_sample",
si_sample = si_sample,
config = make_config(list(n2 = 50,
seed = overall_seed)))
plot(R_si_from_sample, legend=FALSE)
結果可以通過三個圖進行展示，分別為流行曲線圖、Rt圖，代際時間分布圖。
2
無個案代際時間的Rt計算
若無個案代際分布時間，可以參考一些已有的研究，如中國疾控中心團隊在《新英格蘭醫學雜誌》發表的文章發現，COVID-19的代際時間符合gamma分布，均值、標準差分別為7.5和3.4。
以某市公開的COVID-19個案數據為例，計算Rt。由於某市的病例持續存在輸入病例，我們需要區分本地和輸入病例。Incid表示病例數據，需要指定本地病例和輸入病例；代際時間分布的方法選擇參數法parametric_si；config確實參數，mean_si和std_si分別表示代際時間分布的均值和標準差，t_start和t_end用於設置平滑時間。由於每天的Rt波動較大，通常需要對其進行滑動計算，下例表示滑動時間為10天，若不設置t_start和t_end則默認平滑7天。
code
covid19.R
method = "parametric_si",
config = make_config(
mean_si = 7.5,
std_si = 3.4,
t_start = 2:(nrow(data)-10),
t_end = (2+10):nrow(data)))
plot(covid19.R,legend = FALSE, add_imported_cases = TRUE,
options_I = list(col = c("imported" = "red", "local" = "black"),
interval = 1,
ylab = "Incidence"), options_R = list(ylab = "COVID-19 R"))
從結果上看，某市Rt有波動，曾一度接近1，隨後快速下降，全程Rt均保持在1以下。這表明, 某市儘管面臨巨大疫情輸入壓力，但及時採取了一系列強有力防控的措施，有效地控制COVID-19本地疫情傳播，使之在低位水平發展。
出於分析或展示的需要，我們可以從計算結果中將Rt的數據提取出來，保存為數據框。提取數據的代碼如下：
code
eR
dataRgd
RL=covid19.R$R$`Quantile.0.025(R)`,
RH=covid19.R$R$`Quantile.0.975(R)`,
date=seq.Date(from = as.Date(data$date[nrow(data)])-length(eR)+1,
to = as.Date(data$date[nrow(data)]),by="day"))
我們將數據提取出來後，用ggplot2進行重新繪圖。後面我們將針對COVID-19推出可視化的專題，歡迎大家繼續關注。
參考文獻
[1] Thompson R N, Stockwin J E, van Gaalen R D, et al. Improved inference of time-varying reproduction numbers during infectious disease outbreaks[J]. Epidemics, 2019, 29: 100356.
[2] LiQ, Guan X, Wu P, et al. Early transmissiondynamics in Wuhan, China, of novel coronavirus–infected pneumonia[J]. NewEngland Journal of Medicine,2020.
[3] 何冠豪,容祖華,胡建雄,劉濤,肖建鵬,郭凌川,曾韋霖,朱志華,龔德鑫,殷李華,萬東華,吳君樂,康敏,宋鐵,何劍峰,馬文軍.新型冠狀病毒肺炎兩種不同流行模式及其防控效果比較：基於廣州和溫州市的分析[J].中華流行病學雜誌,2020.
製作：胡建雄、何冠豪
審核：劉濤、肖建鵬
指導：馬文軍
關於我們
《數據分析和應用》公眾號隸屬於廣東省公共衛生研究院，旨在分享數據分析方法、案例應用，並總結實踐經驗。團隊學科背景包括衛生統計學、流行病學、環境衛生學、衛生管理學等，具有紮實的統計基礎，歡迎共同開展各類研究。
掃碼關注我們
郵箱：statistic@gdiph.org.cn