基于拟合优度检验的无症状感染者统计模型分析——以上海疫情为例

Statistical Model Analysis of Asymptomatic Patients Based on Goodness of Fit Test—Taking the Shanghai Epidemic as an Example

摘要: 新冠病毒的不断变异和大规模人群接种疫苗等原因，使得疫情爆发地区中的无症状感染者增多。分析和掌握无症状感染者的增长规律，对制定防疫政策和精准防疫措施有重要的参考意义。本文对2022年上海疫情中无症状感染者数据进行统计分布拟合，分析拟合分布的偏态、峰值和分位点，发现伽马分布和对数正态分布与真实分布较为近似，并通过EDF统计量对拟合分布作拟合优度检验，结果显示伽马分布检验p值较大，检验通过率较高，其拟合效果较好。

Abstract: Due to the continuous mutation of novel coronavirus and mass vaccination of the population, the number of asymptomatic patients in regional outbreaks has increased. Analyzing and mastering the growth pattern of asymptomatic infections has important reference significance for formulating epidemic prevention policies and precise epidemic prevention measures. This paper performs sta-tistical distribution fitting on the data of asymptomatic patients in the Shanghai epidemic in 2022, analyzes skewness, peaks, and quantiles of fitted distributions, it is found that the gamma distribu-tion and lognormal distribution are more approximate to the true distribution. In this paper, EDF statistic is used to test the goodness of fit of the fitting distributions. The results show that the p-value of the gamma distribution test is large, and the test pass rate is high, indicating that the fit-ting effect is good.

1. 引言

2022年新型冠状病毒(2019-nCoV)在多地区爆发，但随着病毒的不断变异、疫苗接种率的提升以及疫情管控措施的加强，相对于2020年疫情发展，无症状感染者成为主要病毒感染人群。现阶段新冠疫情中，无症状感染者也具有较强的传染性，同时具有隐匿性，使得无症状感染者的人数增长速度要快于确诊病例，对于疫情防控而言是一个新的挑战。因此分析无症状感染者增长的统计分布，研究其传播规律，对当前疫情走势判断、防控政策的制定有着重要的意义。

统计分布，是用来描述随机现象的基本工具，是数理统计学的基础，任何统计方法都离不开统计分布的概念和各种具体分布的性质 [1]。根据确定的统计分布，可以研究事件发展变化情况，对未知事件的发生进行预测。统计分布分包括离散型统计分布和连续型统计分布，本文将利用连续型分布研究无症状感染者的增长情况。常见的连续型统计分布有均匀分布、正态分布、对数正态分布、伽马分布和威布尔分布等，本文通过参数估计的方法拟合出无症状感染者数据的多种连续型统计分布，并对每一种分布进行拟合优度检验，找到最优的统计模型。拟合优度检验是用来检验观测数与依照某种假设或分布模型计算得到的理论数之间一致性的一种统计假设检验，以便判断该假设或模型是否与实际观测数相吻合，是评价用已知分布拟合现实数据优劣的一种方法，有着广泛的实际应用 [2]。拟合优度检验包括作图法与回归方法、 ${\chi }^2$ 型检验和EDF型检验等，借助EDF统计量检验拟合出的分布是本文的主要研究方法。

2. 统计分布拟合

2.1. 分布模型

2.1.1. 正态分布模型

若随机变量X的概率密度函数为

$f\left(x\right)=\frac{1}{\sqrt{2\pi }\sigma }{\text{e}}^{-\frac{{\left(x-\mu \right)}^2}{2{\sigma }^2}}$

则称随机变量X服从均值为 $\mu$，方差为 ${\sigma }^2$ 的正态分布，记作 $X\sim N\left(\mu ,{\sigma }^2\right)$，其中

$\hat{\mu }=\bar{X}$， $\hat{\sigma }=\frac{1}{n-1}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(X_i-\bar{X}\right)}^2}$ (1)

分别为 $\mu$ 和 ${\sigma }^2$ 的一个无偏估计，那么可以用分布 $N\left(\hat{\mu },{\hat{\sigma }}^2\right)$ 近似随机变量X的真实分布。

2.1.2. 对数正态分布模型

若随机变量X的取值为正数，且 $\ln X\sim N\left(\mu ,{\sigma }^2\right)$，则称随机变量X服从对数正态分布，其概率密度函数为

$f\left(x\right)=\{\begin{array}{l}\frac{1}{x\sqrt{2\pi }\sigma }{\text{e}}^{-\frac{{\left(\ln x-\mu \right)}^2}{2{\sigma }^2}},x>0\\ 0,x\le 0\end{array}$

在对数正态分布的实际应用中，一般取 $\mu$ 和 ${\sigma }^2$ 的极大似然估计作为 $\mu$ 和 ${\sigma }^2$ 估计量，分别为

${\hat{\mu }}_{MLE}=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\ln X_i}$， ${\hat{\sigma }}_{MLE}=\frac{1}{n}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(\ln X_i-\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\ln X_i}\right)}}^2$ (2)

2.1.3. 伽马分布模型

若随机变量X的概率密度函数为

$f\left(x\right)=\{\begin{array}{l}\frac{{\beta }^{\alpha }}{\Gamma \left(\alpha \right)}{x}^{\alpha -1}{\text{e}}^{-\beta x},x>0\\ 0,x\le 0\end{array}$

其中 $\Gamma \left(\alpha \right)={\displaystyle {\int }_0^{+\infty }{x}^{\alpha -1}{\text{e}}^{-x}\text{d}x}$， $\alpha$ 为形状参数， $\beta$ 为尺度参数，则称随机变量X服从伽马分布，记作 $X\sim G\left(\alpha ,\beta \right)$，伽马分布的期望和方差分别为

$E\left(X\right)=\frac{\alpha }{\beta }$， $Var\left(X\right)=\frac{\alpha }{{\beta }^2}$ (3)

用样本均值 $\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{X}_i}$ 代替真实均值，样本方差 $\frac{1}{n-1}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(X_i-\bar{X}\right)}^2}$ 代替真实方差，根据公式(3)可以计算出参数 $\alpha$ 和 $\beta$ 的估计量

$\hat{\alpha }=\frac{{\left(\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{X}_i}\right)}^2}{\frac{1}{n-1}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(X_i-\bar{X}\right)}^2}}$， $\hat{\beta }=\frac{\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{X}_i}}{\frac{1}{n-1}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(X_i-\bar{X}\right)}^2}}$ (4)

2.1.4. 威布尔分布模型

若随机变量X的概率密度函数为

$f\left(x\right)=\{\begin{array}{l}\frac{k}{\lambda }{\left(\frac{x}{\lambda }\right)}^{k-1}{\text{e}}^{-{\left(\frac{x}{\lambda }\right)}^k},x>0\\ 0,x\le 0\end{array}$

其中 $\lambda >0$ 为尺度参数， $k>0$ 为形状参数，则称随机变量X服从威布尔分布，记作 $X\sim W\left(\lambda ,k\right)$。用实际数据进行分布拟合时，一般使用的参数估计为

$\hat{k}=\frac{1.2}{\sqrt{\frac{1}{n-1}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(\ln X_i-\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\ln X_i}\right)}}}$， $\hat{\lambda }=\exp \left(\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\ln X_i}+\frac{0.572}{\hat{k}}\right)$ (5)

2.2. 数据拟合

本文收集了上海2022年3月1日至6月30日共122天的无症状感染者日增长数据，数据来源于上海市卫生健康委官方网站 [3]。为了研究上海新冠疫情无症状感染者的增长规律，本节对日增长数据进行统计分布拟合。在疫情蔓延过程中，日增长数据会出现偶然的波动，为了消除这种偶然因素，提高模型的拟合优度，本节还采用移动平均的方法拟合数据，其中以7天作为平均长度。在数据拟合过程中，将无症状感染者出现的日期作为随机变量，为了便于描述，把3月1日记作上海疫情第1天，6月30日记作上海疫情第122天，则随机变量X的取值为 $1,2,3,\cdots ,122$。由公式(1)、(2)、(4)、(5)，根据无症状感染者日增长数据和7天移动平均数据，计算出四个统计分布的参数估计，具体如表1。

Table 1. Probability distribution parameter estimates for two classes of data

表1. 两类数据的统计分布参数估计值

为了便于观察无症状感染者增长数，本文的分布图采用的是频数分布图。图1和图2的直方图分别展示了无症状感染者日增长数和7天移动平均增长数，图中曲线表示四种拟合分布的频数变化。由图可以发现，无症状感染者前期增长速度较快，呈正偏态分布。在拟合的分布曲线中，对数正态分布和伽马分布是正偏态的，而威布尔分布是负偏态的；同时，对数正态分布和伽马分布对尾部数据拟合较好，能较好地表示疫情后期无症状感染者的增长情况；从峰值上看，四个分布的峰值都低于真实峰值，但对数正态分布和伽马分布的峰值与真实数据的峰值最为接近；从峰值出现的时间上看，真实数据峰值出现的时间与对数正态分布和伽马分布峰值出现的时间近似，而正态分布和威布尔分布的峰值出现的时间要晚。通过以上分析可知，对数正态分布和伽马分布对两类数据拟合的效果较好，而正态分布和威布尔分布拟合效果较差。

(a) 正态分布拟合 (b) 对数正态分布拟合 (c) 伽马分布拟合 (d) 威布尔分布拟合

Figure 1. Frequency plots of fitted distributions for daily growth data

图1. 日增长数据拟合分布频数图

(a) 正态分布拟合 (b) 对数正态分布拟合 (c) 伽马分布拟合 (d) 威布尔分布拟合

Figure 2. Frequency plots of fitted distributions for 7-day moving average growth data

图2. 7天移动平均增长数据拟合分布频数图

中国精算师协会传染病数学模型研究小组在研究2003年非典疫情北京数据数学模型时，利用分位点评价分布拟合效果 [4]，本文运用了同样的方法对分布进行拟合评价。表2和表3中的数据是两类数据真实分布和拟合分布的分位点，通过表格数据可以分析疫情发展情况与时间的关系，例如，真实数据的50%分位点为39，表示第1天至39天无症状感染者增长总数占122天总感染人数的50%，此时疫情发展处于中期。从表中数据可以发现，四个拟合分布的25%、50%、75%分位点与真实分布的分位点都很接近，伽马分布的这三个分位点与真实分布分位点完全相同。若99.9%分位点作为此次疫情即将结束的时间点，则真实分布中第85天进入疫情结束阶段，正态分布和威布尔分布将此时间点提前7至8天，对数正态分布推迟8至9天，而伽马分布仅推迟1至2天。在95%、99%分位点上，同样是伽马分布表现得最好。综上所述，对于疫情发展情况，四个分布在疫情前中期都能够较好地判断，而在疫情发展后期，伽马分布判断较为准确。

Table 2. Quantile of the fitted distributions for daily growth data

表2. 日增长数据拟合分布的分位点

Table 3. Quantile of the fitted distributions for 7-day moving average growth data

表3. 7天移动平均增长数据拟合分布的分位点

3. EDF拟合优度检验

EDF型拟合优度检验是建立在经验分布函数基础上的检验方法，常见的检验方法有KS检验、CvM检验和AD检验。设 $X_1,X_2,\cdots ,X_n$ 是一组来源于连续分布F的独立样本，EDF型检验统计量检验的问题为

$H_0:F=F_0;H_1:F\ne F_0$

其中， $F_0$ 为确定的理论分布函数。

3.1. 检验方法

3.1.1. KS检验

Kolmogorov和Smirnov提出了KS检验 [5] [6]，KS检验统计量为

$D_n=\sqrt{n}\underset{-\infty <x<\infty }{\sup }\left|F_n\left(x\right)-F_0\left(x\right)\right|$，

$D_n^{+}=\sqrt{n}\underset{-\infty <x<\infty }{\sup }\left(F_n\left(x\right)-F_0\left(x\right)\right),D_n^{-}=\sqrt{n}\underset{-\infty <x<\infty }{\sup }(F_0\left(x\right)-F_n(\; x\; ))$

其中 $F_n$ 为经验分布函数。当样本给定时，将样本按照从小到大排序得到 $X_{\left(1\right)},X_{\left(2\right)},\cdots ,X_{\left(n\right)}$，其统计量可以表示为

$D_n=\sqrt{n}\underset{1\le i\le n}{\max }\left[\left(\frac{i}{n}-F_0\left(X_{\left(i\right)}\right)\right),\left(F_0\left(X_{\left(i\right)}\right)-\frac{i-1}{n}\right)\right]$ (6)

3.1.2. CvM检验

CvM型检验统计量 [7] 的一般形式为

${\omega }_n^2=n{\displaystyle {\int }_{-\infty }^{\infty }{\left[F_n\left(x\right)-F_0\left(x\right)\right]}^{2}g\left(x\right)\text{d}{F}_0\left(x\right)}$ (7)

当取 $g\left(x\right)=1$ 时，便得到了CvM检验统计量

$W_n^2=n{\displaystyle {\int }_{-\infty }^{\infty }{\left[F_n\left(x\right)-F_0\left(x\right)\right]}^2\text{d}{F}_0\left(x\right)}$ (8)

当样本给定时，CvM检验统计量可以表示为

$W_n^2={{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left[F_0\left(X_{\left(i\right)}\right)-\frac{2i-1}{2n}\right]}}^2+\frac{1}{12n}$ (8)

3.1.3. AD检验

Anderson和Darling提出了AD检验统计量 [8]，当(7)式中 $g\left(x\right)=\frac{1}{x\left(x-1\right)}$ 时，可以得到AD检验统计量

$A_n^2=n\frac{{\displaystyle {\int }_{-\infty }^{\infty }{\left[F_n\left(x\right)-F_0\left(x\right)\right]}^2}}{F_0\left(x\right)\left[1-F_0\left(x\right)\right]}\text{d}{F}_0(\; x\; )$

当样本给定时，AD检验统计量可以表示为

$A_n^2=-n-\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left[\left(2i-1\right)\ln \left[F_0\left(X_{\left(i\right)}\right)\right]+\left(2n+1-2i\right)\ln \left[1-F_0\left(X_{\left(i\right)}\right)\right]\right]}$ (9)

3.2. 检验p值的计算方法

在进行拟合优度检验时，本文利用p值来确定是否拒绝原假设。对于样本 $X_1,X_2,\cdots ,X_n$ 和拟合分布 $F_0$，EDF型检验统计量检验p值的Monte-Carlo模拟计算步骤如下：

1) 将样本按升序排列得到 $X_{\left(1\right)},X_{\left(2\right)},\cdots ,X_{\left(n\right)}$，再带入公式(6)、(8)、(9)中，分别计算出KS统计量 $D^0$，CvM统计量 $W^0$，AD统计量 $A^0$。

2) 运用Monte-Carl模拟方法，生成n个服从 $U\left[0,1\right]$ 分布的随机变量 $U_1,U_2,\cdots ,U_n$。

3) 随机变量 $U_1,U_2,\cdots ,U_n$ 升序排列得到 $U_{\left(1\right)},U_{\left(2\right)},\cdots ,U_{\left(n\right)}$，令 $F_0\left(X_{\left(i\right)}\right)=U_{\left(i\right)}$，带入公式(6)、(8)、(9)中，计算三种EDF统计量 $D,C,W$。

4) 重复步骤2)和3) M次(一般M大于1000)，得到 $D_1,D_2,\cdots ,D_n;W_1,W_2,\cdots ,W_n;A_1,A_2,\cdots ,A_n$。则KS、CvM和AD统计量的检验p值分别为

$p_{KS}=\frac{\#\left(D>D^0\right)}{M},p_{CvM}=\frac{\#\left(W>W^0\right)}{M},p_{AD}=\frac{\#\left(A>A^0\right)}{M}$

3.3. 检验分析

本文收集的样本量超过59万，这些样本不能严格地服从拟合出的四个分布，若将所有的样本进行拟合优度检验，那么四个分布都不能通过检验，而本文所研究的问题是选择一个近似的分布代替真实分布，所以需要减少检验的样本量。本文采用的方法是将所有的样本分组平均，每组的样本量为N，那么总共分为 $H=\left[N/T\right]$ 组，其中T为所有样本的数量，每组样本取均值生成新的样本 $Y_1,Y_2,\cdots ,Y_H$，则 $Y_k$ $\left(1\le k\le H\right)$ 表示122天中前 $N\left(k-1\right)+1$ 到 $Nk$ 个被感染者出现时间的期望。

表4和表5中显示的是两类数据拟合分布的EDF拟合优度检验的p值，在p值计算过程中取 $M=10000$。两表中数据表明，当N越小时，检验p值越小，这是因为用作检验的样本点越多，拟合优度检验越苛刻，要求分布拟合的优度越高，则拟合分布越不容易通过检验。而无论N取何值，伽马分布的三种检验p值是最大的，而威布尔分布是最小的。取检验水平 $\alpha =0.05$，当 $N=500$ 时，从表4中数据可以发现，对数正态分布通过CvM检验，伽马分布通过了CvM和AD检验，而另外两个分布没有通过任何检验，表5中只有伽马分布通过了所有检验。通过以上分析可知，伽马分布的拟合效果最好。

Table 4. EDF tests p-values of fitted distributions for daily growth data

表4. 日增长数据拟合分布的EDF检验p值

Table 5. EDF tests p-values of fitted distributions for 7-day moving average growth data

表5. 7天移动平均增长数据拟合分布的EDF检验p值

4. 总结

对拟合分布图进行分析可知，伽马分布、对数正态分布与真实分布的偏态方向一致，从而使得这两个分布在峰值位置上与真实分布差距不大，同时伽马分布在尾部数据拟合和分位点数值的接近程度上表现尤为突出。在所有的EDF拟合优度检验中，伽马分布的检验p值最大，说明样本来源于伽马分布最易被接受，由此判断该分布的拟合效果最好。通过以上不同角度的分析，发现对于上海2022年3月至6月的疫情中无症状感染者，其增长模型适合用伽马分布进行拟合，若某地发生类似疫情，可以利用伽马分布分析预测疫情的发展情况，从而制定合理的防疫政策。

NOTES

^*通讯作者。