A data-driven analysis on the mediation effect of compartment models between control measures and COVID-19 epidemics

To assess the influence of control measures on the transmission dynamics of COVID-19 in a quantitative way, here we collect 16 different types of control measures and their respective implementation strength marked by discrete grades from 127 countries/territories (see data file 1 in Supporting Information), which have been promulgated and taken by each government during the first wave of COVID-19 epidemics before July 7, 2020. Although, strictly speaking, the last category – medical resources is not a measure for prevention or containment, we still take it into consideration, because it characterizes the capacity of a country to provide professional medical care of those infected people and thus would influence the cure and mortality rates to a great extent. On the other hand, three groups of official epidemic data, including numbers of daily new cases, cured cases and death cases, are collected from the WHO’s reports [14] (see data file 2 in Supporting Information). For consistency, we redefine the outbreak date of COVID-19 epidemic in each country as the one when the number of accumulative cases in this country reaches 100.

The mediation indicates that the effect of an independent variable on a dependent variable is transmitted through a third variable called mediator, which is frequently utilized in research of applied disciplines, such as psychology, education and social science [21]. To go beyond single variable, one can refer to the Structural Equation Modeling (SEM) [22], which enables simultaneously analyzing multiple independent variables and multiple mediator variables, as well as dealing with manifest and latent variables at the same time.

For our purpose, here we illustrate the basic idea of structural equation modeling for mediation analysis with multiple mediators operating in parallel. Without loss of generality, let $X_{1}$ and $X_{2}$ be two independent variables, $Y$ the dependent variable, and $M_{1}$ and $M_{2}$ the mediators independent of each other. The corresponding linear regression models read:

Here the first line regresses the dependent variable solely on independent variables, the second and third lines regress the mediators on independent variables, while the last line regresses the dependent variable on both the mediators and independent variables. $\tau$’s and $\eta$’s are linear regression coefficients and $\varepsilon$’s denote errors.

Based on above formulas, the differences of coefficients

are recognized as the total mediation effect of $X_{1}$ (and $X_{2}$, respectively) on $Y$, where the constitutive term $\mu_{i}\eta_{ij}$ represents the corresponding individual mediation effect from $X_{j}$ to $Y$ through $M_{i}$ ($i,j=1,2$).

To conduct mediation analysis, the first step is to centralize the data by subtracting mean values from the independent variables, dependent variables, and mediators, respectively. Subsequently, we employ the bootstrap technique [22] to generate a resampled and sufficiently large dataset. Finally, the SEM is used to estimate the undetermined coefficients. If the goodness of fit is acceptable, the significance of the mediation effect is judged based on its confidence interval obtained by the bias-corrected bootstrap. The criteria is, for the mediation effect under testing, if its confidence interval does not contain zero, the corresponding mediation effect is considered to be significant. Otherwise, if the fit of the SEM is unacceptable, the analysis will be stopped.

With respect to the compartment model in this study, the four categories of policies are considered as independent variables $(X_{1},\cdots,X_{4})$, whose indicators are 16 control measures in Table 1. The logarithm of the five coefficients in the SEIR-QD model are used as the mediator variables, i.e. $(M_{1},\cdots,M_{5})=(\log(\alpha),\log(\beta),\log(\gamma),\log(\delta),\log(\kappa))$. Finally, the logarithms of six dynamical features characterizing the spreading of COVID-19 epidemics are taken as the dependent variable one by one, meaning $Y\in\{\log(Q_{100}/N),\log(R_{100}/N),\log(D_{100}/N),\log(t_{1/2}),\log(t_{lag}),\log(k_{app})\}$. The “lavaan” package in R language [23] is used to perform our mediation analysis with structural equation modeling. The corresponding results can be found in data file 3 in Supporting Information.

As a newly emerging highly infectious disease, COVID-19 has some unique characteristics in comparison to other infectious diseases, e.g., a relatively high ratio of asymptomatic patients. Besides, during the first wave of COVID-19 pandemics, public health emergency responses, such as locking down the domestic and international public transportation, quarantining infected patients and so on, had been implemented globally. These facts make the dynamic behaviors of COVID-19 epidemics dramatically different from previous ones.

By taking the specificity of COVID-19 into consideration, the SEIR-QD model (Eqs. (3a)-(3g)) has demonstrated its capability in modeling the transmission dynamics of COVID-19 epidemics and attracted extensive attentions and followers [24, 25]. Compared to the classic SEIR model, three additional compartments have been introduced into this model, including a compartment representing the protected individuals due to the strengthened prevention measures, a compartment for quarantined patients which cannot infect susceptible individuals as well as a separated recovery compartment from the one for death cases.

In summary, the SEIR-QD model contains seven compartments – susceptible cases $S(t)$, protected and insusceptible cases $P(t)$, exposed cases $E(t)$, infected cases $I(t)$, quarantined cases $Q(t)$, recovered cases $R(t)$, and death cases $D(t)$. The transition rates between two connecting compartments are given through six parameters: $\alpha$ characterizing how fast the susceptible cases convert into insusceptible cases, $\beta$ for the rate of susceptible cases getting infected upon contacting with infectious cases, $\gamma$ being the rate of exposed cases convert into infectious cases, $\delta$ as the rate of infectious cases into quarantined, while $\lambda$ and $\kappa$ representing the rate of quarantined cases being recovered or dead. Moreover, $N$ denotes the total population number.

By utilizing the simulated annealing algorithm in MATLAB 2021b [26], the parameters in the SEIR-QD model are fitted with respect to the real epidemic data we collected. In this way, optimal values for the free model parameters are estimated automatically (see data file 4 in Supporting Information).

Due to the significant diversity in the circumstances of each country/territory, it is neccessary to analyze the data for each class, which requires clustering the countries/territories before subsequent analysis. Here the K-means method is employed with the hyperparameter $4$ for the number of data clusters and Euclidean distance as the similarity metric for data samples.

To test the significance of the difference among clusters, the Kruskal-Wallis rank sum test is carried out. The basic idea of Kruskal-Wallis test is to replace the original observation value with its rank for one-way analysis of variance. Suppose there are $M$ samples which are divided into $k$ groups. The sample size in the $i$-th group is $m_{i}$, and $M=\Sigma^{k}_{i=1}m_{i}$. Order all data samples from small to large and record the order number as their rank. In particular, if there are several data at the same position, take their average order number as the common rank. Once the value of Kruskal-Wallis $H$ test statistic

is larger than the significance level, where $R_{i}$ denotes the rank sums, one can conclude the sample groups are statistically significantly different.

The Pearson’s correlation coefficients are used to measure the linear correlation between random variables. Supposing two random variables $A$ and $B$ have $n$ scalar observed values, the Pearson’s correlation coefficient of $A$ and $B$ is defined by

where $\mu_{A}$ and $\sigma_{A}$ are the mean and standard deviation of random variable $A$, $\mu_{B}$ and $\sigma_{B}$ are those of $B$ respectively. The calculation of Pearson’s correlation coefficients can be directly extended to multiple random variables. Here we use the “Corrcoef” function of Matlab to perform the calculation and use the “heatmap” function for visualization.

To provide a relatively comprehensive description of the epidemic dynamics, we refer to six dynamic features, which are constituted by the cumulative numbers of confirmed/cured/dead cases on the 100th day normalized by the total population of the country ($Q_{100}/N$, $R_{100}/N$ and $D_{100}/N$), the half time $t_{1/2}$ (days required for reaching one half of the cumulative number of quarantined cases on the 100th day), the lag time $t_{lag}$ (days required for reaching $5\%$ of the cumulative number of quarantined cases on the 100th day), and the apparent spreading rate $k_{app}$ (the slope of the curve for the cumulative number of confirmed cases at the half time).

The influence of control measures on the transmission dynamics of COVID-19 epidemics is a key issue considered in the current study, and is analyzed by both linear and nonlinear regression models. Here, the multiple linear regression is performed mainly by the least squares method and is supplemented by the stepwise regression and the principal component regression by using the functions “regress”, “stepwise” and “pcacov” in the Statistics and Machine Learning Toolbox of MATLAB [26].

The nonlinear regression analysis is realized through a neural network model – a multilayer perceptron (MLP) constituted by one input layer, one hidden layer with 100 nodes and the ReLU activation function $\sigma(z)=\max(0,z)$, and one output layer using the Sigmoid activation function $\sigma(z)=(1+e^{-z})^{-1}$. The inputs are normalized grades of 16 control measures, while the outputs are values of six dynamical features for the epidemic. The $L_{2}$ norm of the deviations between true values and the predicted ones gives the loss function. The network is trained for 10000 epochs with the Adam optimizer and the learning rate is set to $10^{-3}$. In addition, reasonably adjusting the batch size of the whole dataset will accelerate convergence, which in this case is set as $10$.

As the 16 different types of control measures for 127 countries/territories involved in this study (see Methods and Table 1) are not mutually independent, e.g. strict restrictions on public gathering would be closely related to the measures of school and workplace closure, it is necessary to analyze the correlations among them.

To this end, heatmap is used for descriptive analysis of the control measures (Fig. 1), according to which, it is clear that the sixteen control measures can be aggregated into four categories (Table 1). The first category includes restrictions on public transportation and closure of public places, which reduce the susceptible population and contacts among people. The second category includes policies on contact tracking and testing, which directly influence the extent and efficacy of case detection and documentation. The third category refers to the promotion of testing and the implementation of quarantine measures on suspected cases, reflecting the efficacy of quarantine protocols. The last category includes the status of medical resources.

On the other hand, it has been widely known that there is a dramatic distinction among the control measures and their implementation intensity taken by each nation during the first wave of COVID-19 pandemic. By using the K-means method, these countries are also clustered into 4 groups based on their own implementation strength of control measures. As visualized in Fig. 2, Cluster 1 includes 3 Asian countries, China, South Korea, and Singapore, which takes strict control measures and keeps the pandemic well contained during the early 2020. Cluster 2 includes 53 countries/territories, more than half of which are European countries. Most countries in this group have rich medical resources and take mediate-level measures of prevention and containment. Cluster 3 comprises 61 nations/territories, predominantly encompassing developing countries of the Third World. These countries face inadequate medical resources, and the implementation of COVID-19 control measures within this cluster lacks the necessary stringency to yield desired outcomes. Cluster 4 includes the remaining 10 countries. This cluster is of greater diversity than the preceding three ones. Similar results are obtained by the agglomerative hierarchical clustering method and the density-based spatial clustering and application with noise method, (see Figs. S1, S2 in Supporting Information), which confirm the robustness and reliability of our clustering.

To provide a comprehensive description of the transmission dynamics of the first wave of COVID-19 epidemics, six dynamical features are introduced, that is $Q_{100}/N$, $R_{100}/N$, $D_{100}/N$, $t_{1/2}$, $t_{lag}$, and $k_{app}$. Their definitions are given in the section of Methods, and their values have been extracted from the WHO data (see data file 4 in Supporting Information).

As illustrated in Fig. 3A, the half time $t_{1/2}$ and the apparent spreading rate $k_{app}$ show a clear monotonic dependence on the control strength over the countries in the first three clusters. The non-monotonic dependence for the cumulative confirmed, cured and death cases on the 100th day could be attributed to the fact the epidemics in those less developed countries/territories are still out of control on the 100th day since its outbreak. The extraordinarily long lag-time for countries in Cluster 1 shows that the spreading of COVID-19 virus has been suppressed at a very low level due to the extremely strict containing measures taken by these three countries.

A more quantitative analysis is carried out by adopting multiple linear regression. Based on p-values and R-square values summarized in Figs. 3B-3G, some general conclusions are reached, i.e.

• •

  Except for the lag time, all other five dynamical features are statistically significantly correlated with the overall implementation intensity of control measures in each country, where the p-values for correlation coefficients are all less than $0.01$.

• •

  The impacts of the control measures in Category 2 on the cumulative confirmed/cured/death cases on the 100th day $Q_{100}/N,R_{100}/N,D_{100}/N$ and the spreading rate $k_{app}$ are statistically significant. These findings reveal the importance of keeping track of suspected cases, who have had close contact with confirmed cases, for maintaining the epidemic under control.

• •

  The p values for the linear correlations between the half time/lag time and measures in Category 4, the half-time and measures in Category 3 are less than $0.05$, meanwhile their corresponding R-square values are relatively small too. This fact reveals the intrinsic diversity and complexity of the spreading dynamics of COVID-19 epidemics. Solely raising up the control strength of single (or several) measure may not be as effective as expected.

• •

  With respect to nations in Cluster 2, statistically significant correlations are observed between $Q_{100}/N$ and measures in Category 1, $D_{100}/N$ and measures in Category 1 and Category 4. Meanwhile, as a contrary, correlations are absent for nations in Cluster 3 (with p-values much larger than $0.05$). This result hints that strengthening the public control, like school and workplace closures, can be effective only when the epidemic is still under control. So is the increase of medical resources for reducing the death cases.

In the last section, we have shown that the control measures exhibit statistically significant impacts on the spreading dynamics of COVID-19 epidemics. However, to gain a better understanding on their intrinsic correlations, especially to aid public healthy policy decision, we need to refer to concrete predictive models. Here we make use of the SEIR-QD model in Eqs. (3a)-(3g), which has been shown to be a suitable model for studying COVID-19 epidemics [24, 25].

Based on the mediation analysis with structural equation modeling (see Methods) [22, 23], the SEIR-QD model is found to play a crucial mediation effect between control measures and the epidemic dynamics. To be concrete,

• •

  The total mediation effect of the SEIR-QD model is statistically significant for four epidemic dynamical features, i.e. the normalized cumulative numbers of confirmed/cured/death cases on the 100th day $Q_{100}/N,R_{100}/N,D_{100}/N$ and the spreading rate $k_{app}$.

• •

  The infection rate $\beta$ plays a significant role in the mediation effect of the SEIR-QD model. In particular, the mediation effect along the pathway from control measures in Category 2 to the four dynamical features through the infection rate $\beta$ is the most prominent, as illustrated through red arrows in Fig. 4. This fact can be clearly understood, since more frequent nucleic acid tests and continued tracking of suspected cases allow for a shorter time to identify these infected cases and isolate them from normal people, and thus effectively reduce the infection rate $\beta$.

• •

  The mediation effect along the pathway from control measures in Category 2 to $Q_{100}/N$ through the protection rate $\alpha$ is also statistically significant. It shows that the inclusion of the protection rate $\alpha$ into the classical SEIR model is crucial for proper characterization of the impacts of control measures on the COVID-19 epidemics.

Our above findings agree with the report by Zhou et al. [13], in which they also highlighted the key role of tracing and testing during the outbreaks of Delta variant in Xi’an, Yangzhou and Guangzhou in China.

The SEIR-QD model provides a quantitative way to estimate the impacts of various control measures and their implementation intensities on the transmission dynamics of COVID-19 epidemics. To begin with, we argue that there is a close relation between the control measures and parameters in the SEIR-QD model.

• •

  Control measures in Category 1 mainly affects the protection rate $\alpha$ and infection rate $\beta$ in the SEIR-QD model. Increasing the control strength of policies in Category 1 will lead to a larger $\alpha$ and a smaller $\beta$.

• •

  In addition to the infection rate $\beta$, the control measures of Category 2 adjust the quarantine rate $\delta$ too. More frequent nuclei acid testing and continuous tracking of suspected cases allow a shorter time to identify those infected cases and isolate them from normal people. This statement is directly supported by the linear correlation analysis, which shows that the correlations for both $\beta$ and $\delta$ are statistically significant, as well as the mediation analysis, which gives an affirmative answer to the coefficient $\beta$.

• •

  The linear correlations between control measures in Category 3 & 4 and the mortality rate $\kappa$ are statistically significant. Richer medical resources mean patients can get more professional care and medical treatment in general, resulting in a decrease in the mortablity rate $\kappa$.

Now based on above statements, the parameters corresponding to control measures in each category are adjusted by $\pm 20\%$ with respect to their default values, while the rest parameters in the SEIR-QD model are kept unchanged to reflect the different policy impacts. As depicted in Fig. 5 through three representative countries – South Korea for Cluster 1, Germany for Cluster 2, and Mexico for Cluster 3, the COVID-19 epidemics in countries belonging to different clusters exhibit significantly distinguishable reliance on control measures. For instance, a twenty percentage relaxation on the implementation intensity of measures in either Category 1 or Category 2 results in limited perturbations on the daily new quarantined cases in South Korea, whereas in Mexico, the impact becomes exceptionally large. Germany falls somewhere in between. The influence of control measures in Category 2 is weaker compared to that of measures in Category 1. These observations accurately reflect the varying situations in these countries. Interested readers are referred to the Supporting Information for data on other representative countries for further comparison and validation.

In addition, a global sensitivity analysis is performed and summarized in Fig. 6. As measured by the normalized cumulative numbers of confirmed/cured/death cases on the 100th day, $Q_{100}/N,R_{100}/N,D_{100}/N$, tightening measures in Category 1 seems to be more effective than those in Category 2. However, the cost for shutting down public transportations and public places is much higher too. Furthermore, by increasing the control strength of measures in Category 1, the fold changes in $Q_{100}/N$ among four country clusters show a monotonic behavior, meaning its average impact on Cluster 1 is the strongest while on Cluster 3 and 4 are the weakest. Similar conclusions could be reached for $R_{100}/N$ and $D_{100}/N$ too. Contrarily, non-monotonic behaviors are observed for the impacts of control measures in Category 2 among four country clusters. The average influence of nuclei acid testing and keeping tracing suspected cases for countries in Cluster 3 is more apparent than expected. Especially, control measures in Category 1 and 2 exhibit opposite impacts on the half time. This fact shows that performing nuclei acid testing and keeping tracing suspected cases only suppress the outbreak of COVID-19 epidemics and thus prolong the time to reach its peak. In summary, the policy impacts on the transmission dynamics of COVID-19 are significantly different among four country clusters.

Along with above compartment-model-based analysis, the multiple linear regression and nonlinear regression by MLP are carried out too. The results by MLP show an excellent fitting to the real epidemic data, while the performance of a multiple linear regression model is moderate. However, neither linear nor nonlinear regression models offer consistent results during the sensitivity analysis (Figs. S4, S10, S11 in Supporting Information), and thus are not suitable for evaluating the policy impacts on the epidemic dynamics.

The online version contains supplementary material available at https://doi.org/xxxxxxxxxxx.

Not applicable.

Not applicable.

The datasets generated during and/or analyzed in this study are available from the corresponding authors on reasonable request.

The authors declare that they have no confict of interest.

This work was supported by the National Natural Science Foundation of China (11801020, 12205135, 21877070), Guangdong Basic and Applied Basic Research Foundation (2023A1515010157), the Natural Science Foundation of Fujian Province of China (2020J05172), Startup Research Funding of Minjiang University (mjy19033), Start-up Research Fund of Yanqi Lake Beijing Institute of Mathematical Sciences and Applications, Special Pre-research Project of Beijing University of Technology for Fighting the Outbreak of Epidemics.

LH and CZ designed and implemented the study, collected the data, analyzed the data, writing–original draft, writing–review & edit. DZ, WW and WY analyzed the data, writing–original draft, writing–review & edit. LP writing–original draft, writing–review & edit. All authors have read and approved the content of the manuscript. All authors read and approved the final manuscript.

Not applicable.