A novel statistical method for long-term coronavirus modelling

Introduction

Statistical aspects of coronavirus disease 2019 (COVID-19) and other similar recent epidemics were receiving much attention in the modern research community.¹ Generally, it is quite challenging to calculate realistic biological system reliability factors and outbreak probabilities under actual epidemic conditions by using conventional theoretical statistical methods.^1–14 The latter is usually due to many degrees of system freedom and random variables governing dynamic biological systems, spread over extensive terrain. In principle, the reliability of a complex biological system may be accurately estimated straightforwardly by having enough measurements or by direct Monte Carlo simulations.¹⁰ For COVID-19, however, the only available observation numbers are limited as the observations are only available from the beginning of the year 2020 up to now. Motivated by the latter argument, the authors have introduced a novel reliability method for biological and health systems to predict and manage epidemic outbreaks more accurately, this study was focused on COVID-19 epidemics in northern China, with focus on cross-correlations between different provinces within same climatic zone. For other studies related to statistical variations per country see,¹⁵ where spatial lags were addressed. China was chosen to test the methodology proposed in this study, because of its COVID-19 origin and extensive health observations and related research available online.^15–29

Statistical modelling of lifetime data or extreme value theory (EVT) is widespread in medicine and engineering. For example, Gumbel used EVT to estimate the demographic of various populations in.³⁰ Similarly, in Ref. 31, the author builds on this EVT theory to determine the fitness effect using a Beta-Burr distribution. While in Ref. 32, the author discusses using a bivariate logistic regression model, which was then used to access multiple sclerosis patients with walking disabilities and in a cognitive experiment for visual recognition. Finally, Ref. 33 is a paper of relevance, which used EVT to estimate the probability of an influenza outbreak in China. The author demonstrated a forecasting prediction potential amid the epidemic in this paper. While in Ref. 34 similarly used EVT to predict and detect anomalies of influenza epidemics.

In this paper epidemic outbreak is viewed as unexpected incident that may occur at any province of a given country at any time, therefore spatial spread is accounted for. Moreover, a specific non-dimensional factor $\lambda$ is introduced to predict the latter epidemic risk at any time and any place.

Biological systems are subjected to ergodic environmental influences. The other alternative is to view the process as being dependent on specific environmental parameters whose variation in time may be modelled as an ergodic process on its own. The incidence data of COVID-19 in all provinces of the People's Republic of China (PRC) from February 2020 until today were retrieved from the official public PRC health website, for simplicity only northern provinces were selected for this study. As this dataset is organized by province (more than 30 provinces in China), the biological system under consideration can be regarded as a multi-degree of freedom (MDOF) dynamic system with highly inter-correlated regional components/dimensions. Some recent studies have already used statistical tools to predict COVID-19 development, for linear log model see Ref. 9, these studies however did not address fully dynamic space-time dynamic bio-system as this study does.

Note that while this study aims at reducing risk of future epidemic outbreaks by predicting them, it is solely focused on daily registered patient numbers and not on symptoms themselves. For long-lasting COVID-19 symptoms, the so-called “long COVID”, and its risk factors and whether it is possible to predict a protracted course early in the disease, see Ref. 35, for mortality research see Ref. 36.

Methods

This section presents novel bio-system reliability method, presently invented by authors and previously not existing. Advocated methodology used generalized extreme value (GEV) theory as a basis for its buildup. For numerical part authors used commercial software MATLAB, (Mathworks, V 8.6), namely its optimization routines, otherwise authors have used extrapolation code available from ACER. Only the code available from ACER was used to complete all sections of the methods – both numerical part as well as final extrapolation.

The novel method introduces MDOF (multi-degree of freedom) health response vector process $\mathit{R}\left(t\right)=\left(X\left(t\right),Y\left(t\right),Z\left(t\right),\dots \right)$ that was measured over a sufficiently long period of time$\left(0,T\right)$. Unidimensional global maxima over the entire time span $\left(0,T\right)$ denoted as $X_T^{\max }=\underset{0\le t\le T}{max}X\left(t\right)$, $Y_T^{\max }=\underset{0\le t\le T}{max}Y\left(t\right)$, $Z_T^{\max }=\underset{0\le t\le T}{max}Z\left(t\right),\dots$.

Let $X_1,\dots ,X_{N_X}$ be time local maxima of the process $X\left(t\right)$ consequent in time, recorded at discrete time instants $t_1^X<\dots <t_{N_X}^X$ that are monotonously increasing in$\left(0,T\right)$. A similar definition follows¹ for other MDOF response components $Y\left(t\right),Z\left(t\right),\dots$ with $Y_1,\dots ,Y_{N_Y};$ $Z_1,\dots ,Z_{N_Z}$ and so on. For simplicity, all $\mathit{R}\left(t\right)$ components, and therefore its maxima are assumed to be non-negative.

The target is to estimate system failure probability, in other words the probability of exceedance

More specifically, the moment when either $X\left(t\right)$ exceeds ${\eta }_X$, or $Y\left(t\right)$ exceeds ${\eta }_Y$, or $Z\left(t\right)$ exceeds ${\eta }_Z$, and so on, the system is regarded as immediately failed. Fixed failure levels ${\eta }_X$, ${\eta }_Y$, ${\eta }_Z$,…are, of course, individual for each unidimensional response component of $\mathit{R}\left(t\right)$. $X_{N_X}^{\max }=max\left\{X_j;j=1,\dots ,N_X\right\}=X_T^{\max }$, $Y_{N_Y}^{\max }=max\left\{Y_j;j=1,\dots ,N_Y\right\}=Y_T^{\max }$, $Z_{N_z}^{\max }=max\left\{Z_j;j=1,\dots ,N_Z\right\}=Z_T^{\max }$, and so on.

Now, the local maxima time instants $\left[t_1^X<\cdots <t_{N_X}^X;t_1^Y<\cdots <t_{N_Y}^Y;t_1^Z<\cdots <t_{N_Z}^Z\right]$ are sorted in monotonously non-decreasing order into one single merged time vector $t_1\le \cdots \le t_N$. Note that $t_N=max\left\{t_{N_X}^X,t_{N_Y}^Y,t_{N_Z}^Z,\cdots \right\}$, $N=N_X+N_Y+N_Z+\dots$. In this case $t_j$ represents local maxima of one of MDOF structural response components either $X\left(t\right)$ or $Y\left(t\right)$, or $Z\left(t\right)$ and so on. That means that having $\mathit{R}\left(t\right)$ time record, one just needs continuously and simultaneously screen for unidimensional response component local maxima and record its exceedance of MDOF limit vector $\left({\eta }_X,{\eta }_Y,{\eta }_Z,\dots \right)$ in any of its components $X,Y,Z,\dots$ Local unidimensional response component maxima are merged into one temporal non-decreasing vector $\overrightarrow{R}=\left(R_1,R_2,\dots ,R_N\right)$ following the merged time vector $t_1\le \cdots \le t_N$, see Figure 1. That is to say, each local maxima $R_j$ is, in fact, actual encountered local maxima corresponding to either $X\left(t\right)$ or $Y\left(t\right)$, or $Z\left(t\right)$ and so on. Finally, the unified limit vector $\left({\eta }_1,\dots ,{\eta }_N\right)$ is introduced with each component ${\eta }_j$ is either ${\eta }_X$, ${\eta }_Y$ or ${\eta }_Z$ and so on, depending on which of $X\left(t\right)$ or $Y\left(t\right)$, or $Z\left(t\right)$ etc., corresponds to the current local maxima with the running index$j$.

Figure 1. Illustration on how two exemplary processes X and Y are combined into new synthetic vector R.

Now, the scaling parameter$0<\lambda \le 1$ is introduced to, artificially, simultaneously decrease limit values for all response components, namely the new MDOF limit vector $\left({\eta }_X^{\lambda },{\eta }_Y^{\lambda },{\eta }_z^{\lambda },\dots \right)$ with ${\eta }_X^{\lambda }\equiv \lambda \bullet {\eta }_X$, $\equiv \lambda \bullet {\eta }_Y$, ${\eta }_z^{\lambda }\equiv \lambda \bullet {\eta }_Z$, … is introduced, see Ref. 37. The unified limit vector $\left({\eta }_1^{\lambda },\dots ,{\eta }_N^{\lambda }\right)$ is introduced with each component ${\eta }_j^{\lambda }$ is either ${\eta }_X^{\lambda }$, ${\eta }_Y^{\lambda }$or ${\eta }_z^{\lambda }$ and so on. The latter automatically defines probability $P\left(\lambda \right)$ as a function of $\lambda$, note that $P\equiv P\left(1\right)$ from Equation (1). Non-exceedance probability $P\left(\lambda \right)$ can be introduced as follows

Next, cascade of approximations which is based on conditioning is briefly outlined. In practice,^11–14,38 the dependence between the neighboring $R_j$ is not negligible; thus, the following one-step (will be called conditioning level $k=1$) memory approximation is introduced

where $3\le j\le N$ (will be called conditioning level $k=3$), and so on. The idea is to monitor each independent failure that happened locally first in time, thus avoiding cascading local inter-correlated exceedances.

Equation (4) exhibits subsequent refinements with respect to the statistical independence assumption. These approximations increasingly accurately approximate statistical dependence between neighboring maxima. Since the original MDOF process $\mathit{R}\left(t\right)$ has been assumed ergodic and thus stationary, the probability $p_k\left(\lambda \right)≔\text{Prob}\left\{R_j>{\eta }_j^{\lambda }\right|R_{j-1}\le {\eta }_{j-1}^{\lambda },R_{j-k+1}\le {\eta }_{j-k+1}^{\lambda }\}$ for $j\ge k$ is independent of $j$, however dependent on the conditioning level $k$. Thus, non-exceedance probability may be approximated as in the average conditional exceedance rate method, see Refs. 30, 37 for more details on exceedance probability

Note that Equation (5) follows from Equation (1) if neglecting$\text{Prob}\left(R_1\le {\eta }_1^{\lambda }\right)\approx 1$, as design failure probability must be epsilon order o(1), with $N\gg k$. Equation (5) is analogous to the well-known mean up-crossing rate equation for the stochastic process probability of exceedance.^10,37 There is convergence with respect to $k$, called here conditioning level

Note that Equation (5) for $k=1$ is equivalent to a well-known non-exceedance probability relationship with the mean up-crossing rate function

where ${\nu }^{+}\left(\lambda \right)$ denotes the mean up-crossing rate of the response level $\lambda$ for the above assembled non-dimensional vector $R\left(t\right)$ assembled from scaled MDOF system response $\left(\frac{X}{{\eta }_X},\frac{Y}{{\eta }_Y},\frac{Z}{{\eta }_Z},\dots \right)$. The mean up-crossing rate is given by the Rice's formula given in Equation (7) with $p_{R\dot{R}}$ being joint probability density for $\left(R,\dot{R}\right)$ with $\dot{R}$ being time derivative$R^{\prime }\left(t\right)$. Equation (7) relies on the Poisson assumption that is up-crossing events of high $\lambda$ levels (in this paper, it is $\lambda \ge 1$) can be assumed to be independent. The latter may not be the case for narrowband responses and higher-level dynamical systems that exhibit cascading failures in different dimensions, subsequent in time, caused by intrinsic inter-dependency between extreme events, manifesting itself in the appearance of highly correlated local maxima clusters within the assembled vector $\overrightarrow{R}=\left(R_1,R_2,\dots ,R_N\right)$.

In the above, the stationarity assumption has been used. However, the proposed methodology can also treat the nonstationary case. For nonstationary case, the scattered diagram of $m=1,..,M$ seasonal epidemic conditions, each short-term seasonal state has the probability$q_m$, so that${\sum }_{m=1}^M{q}_m=1$. Next, let one introduce the long-term equation

with $p_k\left(\lambda ,m\right)$ being the same function as in Equation (6) but corresponding to a specific short-term seasonal epidemic state with the number $m$.

Note that the accuracy of the suggested approach for a large variety of one-dimensional dynamic systems was successfully verified by authors in previous years.^10,37

Use case

Methods described in this paper are novel and state of art. Prediction of influenza type epidemics has long been the focus of attention in mathematical biology and epidemiology. It is known that public health dynamics is a seasonally and spatially varying dynamic system that is always challenging to analyse. This section illustrates the efficiency of the above-described methodology using the new method applied to the real-life COVID-19 data sets, presented as a new daily recorded infected patient time series, spread over different regions.

COVID-19 and influenza are contagious diseases with high transmissibility and ability to spread. Seasonal influenza epidemics caused by influenza A and B viruses typically occur annually during winter, presenting a burden on worldwide public health, resulting in around 3–5 million cases of severe illness and 250,000–500,000 deaths worldwide each year, according to the World Health Organization (WHO).³

This section presents a real-life application of the above-described method. The statistical data in the present section are taken from the official website of the National Health Commission of the people's Republic of China. The website provides the number of newly diagnosed cases every day in 13 administrative regions in northern China from 22 January 2020 to 6 April 2022. Patient numbers from thirteen different Chinese administrative regions were chosen as components $X,Y,Z,\dots$ thus constituting an example of a thirteen dimensional (13D) dynamic biological system.

In order to unify all three measured time series $X,Y,Z$, the following scaling was performed as follows

Next, all local maxima from three measured time series were merged into one single time series by keeping them in time non-decreasing order: $\overrightarrow{R}=\left(max\left\{X_1,Y_1,Z_1,\dots \right\},\dots ,max\left\{X_N,Y_N,Z_N,\dots \right\}\right)$ with the whole $\overrightarrow{R}$vector being sorted according to non-decreasing times of occurrence of these local maxima.

Figure 2 presents new daily recorded patients number plotted as a time-space 2D surface using MATLAB. Figure 3 presents the number of new daily recorded patients as a 13D vector $\overrightarrow{R}$, consisting of assembled regional new daily patient numbers. Note that vector $\overrightarrow{R}$ does not have physical meaning on its own, as it is assembled of different regional components with different epidemic backgrounds. Index $j$ is just a running index of local maxima encountered in a non-decreasing time sequence.

Figure 2. New daily recorded patients number plotted as a 2D time-space surface, data was according to http://www.nhc.gov.cn.

Figure 3. Number of new daily recorded patients as 13D vector $\overrightarrow{R}$.

Left: as it is, Right: scaled by Equation (10).

Figure 4 presents 100 years return level extrapolation according to Equation (9) towards epidemic outbreak with 100 year return period, indicated by the horizontal dotted line, and somewhat beyond, $\lambda =0.1$ cut-on value was used. Dotted lines indicate extrapolated 95% confidence interval according to Equation (10). According to Equation (5) $p\left(\lambda \right)$ is directly related to the target failure probability $1-P$ from Equation (1). Therefore, in agreement with Equation (5), system failure probability $1-P\approx 1-P_k\left(1\right)$ can be estimated. Note that in Equation (5), $N$ corresponds to the total number of local maxima in the unified response vector$\overrightarrow{R}$.Conditioning parameter$k=3$ was found to be sufficient due to occurrence of convergence with respect to $k$, see Equation (6). Figure 4 exhibits reasonably narrow 95% CI. The latter is an advantage of the proposed method.

Figure 4. 100 years return level (horizontal dotted line) extrapolation of $p_k\left(\lambda \right)$ towards critical (system failure) level (indicated by star) and beyond.

Extrapolated 95% CI indicated by dotted lines, using https://github.com/cran/acer.

Note that while being novel, the above-described methodology has a clear advantage of utilizing available measured data set quite efficiently due to its ability to treat health system multi-dimensionality and perform accurate extrapolation based on quite limited data set. Note that, predicted non-dimensional $\lambda$ level, indicated by star in Figure 4, represents probability of epidemic outbreak at any northern province in China in the years to come.

Discussion

Traditional health bio-systems reliability methods dealing with observed/measured time series do not have the ability to efficiently deal with high dimensionality and cross-correlation between different system responses. The main advantage of the introduced methodology is its ability to study high dimensional non-linear dynamic systems reliability.

Despite the simplicity, the present study successfully offers a novel multidimensional modelling strategy and a methodological avenue to implement the forecasting of an epidemic during its course, if it is assumed to be stationary in time. Proper setting of epidemiological alarm limits (failure limits) per province has been discussed, see Section Use case.

This paper studied recorded COVID-19 patient numbers from thirteen different Chinese northern provinces, constituting an example of a thirteen dimensional (13D) and ten-dimensional (10D) dynamic biological system respectively observed in 2020-2022. The novel reliability method was applied to new daily patient numbers as a multidimensional system in real-time.³⁰

The main conclusion is that if the public health system under local environmental and epidemiologic conditions in northern China is well managed. Predicted 100 year return period risk level $\lambda$ of epidemic outbreak is reasonably low. However, there is an ultra-low risk of a future epidemic outbreak in both countries, at least in 100 years horizon.

This study outlines a general-purpose, robust and straightforward multidimensional bio-system reliability method. The method introduced in this study has been previously successfully validated by application to a wide range of engineering models,^30,40–49 but only for one and two-dimensional system responses and, in general, very accurate predictions were obtained. Both measured and simulated time series responses can be analysed using the proposed method. It is shown that the method produced an acceptable 95% confidence interval, see Figure 4. Thus, the suggested methodology may be used as a tool in various non-linear dynamic biological systems reliability studies. The presented COVID-19 example does not limit potential areas of new method applicability by any means.

Data availability

The datasets analyzed during the current study are publicly available from the daily recorded patients dataset in China during 2020-2022 years, are available at http://www.nhc.gov.cn/. Daily recorded patients data was structured per province and per calendar day, namely it was straightforward to extract and systematize joint statistical distribution as a function of both space and time.

Software availability

Source code along with demo and user manual and examples used for extrapolation available from: ACER this third-party software is under license GPL-3.