Specification testing for ordinary differential equation models with fixed design and applications to COVID-19 epidemic models.

Checking the models about the ongoing Coronavirus Disease 2019 (COVID-19) pandemic is an important issue. Some famous ordinary differential equation (ODE) models, such as the SIR and SEIR models have been used to describe and predict the epidemic trend. Still, in many cases, only part of the equations can be observed. A test is suggested to check possibly partially observed ODE models with a fixed design sampling scheme. The asymptotic properties of the test under the null, global and local alternative hypotheses are presented. Two new propositions about U-statistics with varying kernels based on independent but non-identical data are derived as essential tools. Some simulation studies are conducted to examine the performances of the test. Based on the available public data, it is found that the SEIR model, for modeling the data of COVID-19 infective cases in certain periods in Japan and Algeria, respectively, maybe not be appropriate by applying the proposed test.

Introduction

Since December 2019, a new infectious disease called the Coronavirus Disease 2019 (COVID-19) has spread worldwide, resulting in an ongoing pandemic. To help understand the characteristic of this epidemic and to act suitably to reduce its risk, numerous studies have been done by using several approaches to estimate the sizes of the COVID-19 infected populations in different countries and regions and further predict their trends (e.g., Lin et al. (2020); Altieri et al. (2021); Yang et al. (2020); Tian et al. (2020)). The modeling methods of ordinary differential equations (ODEs) are popularly used by using the mechanisms behind the data, including the Susceptible-Infectious-Removed (SIR) model, the Susceptible-Exposed-Infectious-Removed (SEIR) model, and several other variants. Note that we may only have part of the model equations observed from the public data. This is often the case in practice. Thus, in this paper, we particularly pay attention to such models.

Consider a general d-dimensional ODE model as with an initial condition . Here is a given parametric family of functions. The ODE system is commonly observed with noise as where are chosen from and . The measurement error satisfying has a nonsingular variance-covariance matrix , and is independent with for . If , the ODE system is fully observed, otherwise is a partially observed ODE system, where only a small subset of components can be measured (Dattner (2015)).

Using ODE models to describe dynamic processes can estimate parameters with particular relevance and conduct further statistical analysis. However, an important issue is whether the assumed ODE model rightly describes the measurement data. For instance, we may want to know whether an SEIR model is good enough to model real COVID-19 infective population data for a period of time. If not, we may have to be careful to use the results from the statistical analysis. Thus, a model checking for the assumed ODE model should be accompanied. Given a certain parametric family , the hypotheses we consider are formally stated as where is an unknown parameter vector.

The hypothesis testing problem (2) is a model specification testing problem. Some relevant tests, usually under random design, concerning regression model specification are built by using nonparametric estimations, and thus are called the local smoothing tests. Some literature, such as Fan and Li (2000), points out that local smoothing methods are sensitive to alternative models which are oscillating/highly frequent than global smoothing methods in general. Examples contain Härdle and Mammen (1993), Zheng (1996), Koul and Ni (2004) and Lavergne and Patilea (2012). Some others are based on residual-marked empirical processes and the average over an index set. Then they are called the global smoothing tests. Examples include Stute (1997), Stute et al. (1998a), Stute et al. (1998b), Zhu (2003) and Khmaladze and Koul (2004). These tests can better detect smooth alternative models. However, testing for ODE models specification still receives less attention. For testing the hypotheses in (2), Hooker (2009) proposed a goodness-of-fit test based on estimated forcing functions. It has an exact null distribution under the normality assumption on independent components and the homoscedasticity structure on error terms. It may have difficulty controlling the size when the above assumptions do not hold. As known from the qualitative theory of differential equations, an ODE model often has a periodic non-constant solution, which is called a limit cycle. Thus, the original function often oscillates. See Hirsch et al. (2013) as a good reference. This phenomenon suggests that local smoothing tests for model checking should be more appropriate than global smoothing tests. Liu et al. (2021) constructed a local smoothing test under more relaxed assumptions on error terms and gave its asymptotic properties with a random design. Their trajectory matching-based test, written as , takes advantage of the solution trajectory of the ODE system. We note that the good theoretical results of rely on a random design sampling scheme. The sample 's are independent and identically distributed (i.i.d.), thus the limiting properties can be derived by using the existing theory of U-statistics. But for ODE systems such as the epidemic model, the variable t usually represents time, and the response is measured at certain time points in practice. Therefore, it is suitable to treat under fixed design rather than i.i.d. random design. The test that is based on this non-i.i.d. sequence of 's needs to be investigated to see its usability.

In this paper, we give a modification of with a fixed design and study the respective properties under the null, local and global alternatives. Without notational confusion, we still write it as . We then check the SEIR model for two real COVID-19 data sets. Since the sample are still independent in the case of fixed design, we also consider U-statistics with varying kernels based on independent but not identically distributed random variables by giving two new propositions concerning degenerate and non-degenerate U-statistics.

The paper is organized as follows. In Section 2, we will present some general asymptotic properties of U-statistics. Section 3 will contain the construction of the test and its asymptotic properties. Section 4 will include some simulation results. In Section 5, we will use the proposed test to check whether a classic SEIR model is tenable for describing the tendency of COVID-19 infective cases from 15 January to 29 February, 2020 in Japan and from 25 February to 18 April, 2020 in Algeria respectively. A summary and a brief discussion will be given in Section 6. The technical proofs will be relegated to Appendix.

U-statistics with varying kernels

The theory of U-statistics with i.i.d. random variables has been well developed (see Hoeffding (1948); Lee (1990)). For U-statistics with varying kernels, Powell et al. (1989) established the asymptotic equivalence of a U-statistic and its projection when the U-statistic is non-degenerate. In degenerate case, Hall (1984) proved a central limit theorem, showing that a second order U-statistic has a limiting normal distribution under some regularity conditions. Fan and Li (1996) and Fan and Li (1999) extended Hall's result to more complicated cases. However, the theory of U-statistics with varying kernels based on independent but non-identically distributed data remains unknown. We then study the asymptotic properties of both degenerate and non-degenerate U-statistics in the non-i.i.d. scenarios.

Consider a second order U-statistic : with the kernel that varies with the sample size n. Recall that is an independent but non-identical sample. Assuming the expectation of exists for every and , we define is the first order projection of . If , the U-statistic is degenerate, otherwise is non-degenerate. The following proposition is for the non-degenerate case, which is exactly a generalization of Lemma 3.1 of Powell et al. (1989). Since and are asymptotically equivalent, we can use this proposition to establish the asymptotic normality of a non-degenerate U-statistic by using the central limit theorem.

Supposeis non-degenerate. Iffor everyand, then.

For a degenerate , we assume is centered without loss of generality. Denote . We state the following result that is a generalization of Theorem 1 of Hall (1984) to the non-i.i.d. case.

For every,and n, supposealmost surely and, if asthenis asymptotically standard normal.

The test with fixed design

Test statistic construction

We first define the fixed design sampling scheme. This condition is used to establish the asymptotic equivalence between the empirical distribution of and a given distribution , which is commonly used in fixed design sampling (Xue et al. (2010)). Hereafter we denote .

The sample pointsare fixed. There exists a distributionwith the corresponding density functionsuch thatwhereis the empirical distribution of. The density functionis bounded away from 0 and ∞. The first derivative of is bounded and continuous.

Denote as the unique solution trajectory of the ODE system . By using the trajectory , we can convert the hypothesis testing problem (2) to the problem of checking the vector function for some .

The following construction is a modification of the test proposed by Liu et al. (2021) that is for the problem with random design. Let with be the residual of the k-th component. Under , and leads to , while under , , and . At the sample level, we can use as an estimator of with the nonlinear least-squares estimator . Then, for each observed component k, we construct a test statistic that is the sample analogue of as follows: where h is the bandwidth parameter and is the kernel function. is in spirit similar to the test proposed in Zheng (1996). We compose all o tests to form a vector . The final test statistic is constructed by aggregating all the components of : with a normalized symmetric matrix : where and ⊙ means the element-wise product of two vectors.

Asymptotic properties

To investigate the asymptotic properties, we first give certain regularity conditions.

,satisfy. They have nonsingular variance-covariance matrix.is independent withfor every.

All partial derivatives ofup to 2 with respect to X and t are existent and continuous.

is a nonnegative, bounded, continuous and symmetric function with.

For all,is continuously differentiable and bounded by a measurable functionwith.

The parameter space Θ is a convex compact subset of.

is a Borel measurable real function onfor each θ and is twice continuously differentiable with respect to θ for each t.

For all, we have

.

The matrixis nonsingular.

Condition (1) is basic. Condition (2) assumes that the ODE system has a unique solution trajectory. The commonly used requirement for the kernel function is put in Condition (3). Conditions (4)-(9) are standard for the nonlinear least squares estimation, which is similar to the conditions given in White (1981), and Xue et al. (2010).

It can be shown that and are asymptotic U-statistics with non-i.i.d. observations. Thus we use Proposition 2.1, Proposition 2.2 in Section 2 to derive their asymptotic properties. The limiting null distribution of and its consistency under the global alternatives are stated in the following theorem.

Given Conditions (1)-(9), ifand, then

Under the null hypothesis,where, andis the chi-square distribution with o degrees of freedom.

Define anmatrixwith theelement, and. Then under the global alternative,whereis an o-dimensional vector whose k-th component is equal to.

We also study the power performance of under the following sequence of local alternative hypotheses: with the bounded multiple response function and the departure parameter . Then we give the following theorem which states the asymptotic power performance of under . The results show that the test has a typical rate existing local smoothing tests can achieve to detect local alternatives converging to the null model.

Given Conditions (1)-(9), ifand, then underwith,where μ is an o-dimensional vector with the i-th componentParticularly, if,whereis the noncentral chi-squared distribution where the noncentrality parameterwith.

Simulations

We conduct five simulation examples to evidence the performance of in finite sample scenarios. In Example 1, the null ODE system is linear. Examples 2 - 4 respectively use three nonlinear ODE models as the null models, which have been often used in economic, neuroscience and ecology. Example 5 reconsiders the above four models but with partially observed ODE systems. In this example, we assume that we only have the data of the second component. Furthermore, we consider both independent and highly correlated components of the error terms for each example.

For the ODE models used in the examples, the variable t can be designed time point taking values . With this fixed design, we apply the test for each cases. As a comparison, the test proposed by Hooker (2009) is computed as the competitor.

Given the ODE models, the trajectory is obtained by using the 4-stage Runge-Kutta algorithm such that we have the respective models . The error terms with independent components , independently follow the normal distribution , where is the identity matrix of size 2. Then the error terms with highly correlated components are made by , where . The initial values of ODE models are assumed to be known. We conduct 1000 experiments for each simulation case. The sample size is 300, and the significance level is 0.05.

To construct the test, we use the Matlab algorithm in OPTI Toolbox (Currie and Wilson (2012)) to get the nonlinear least-squares estimator . The trajectory is also obtained by utilizing the 4-stage Runge-Kutta algorithm. We choose the Epanechnikov kernel in the test construction. As for the choice of bandwidth h, we found a small h is helpful to ensure the asymptotic properties theoretically and makes the test has better performances in the simulations empirically. Thus we recommended a rate of , which is faster than the optimal rate of for the kernel estimation. By the rule of thumb, the coefficient can be chosen as , where is the length of the whole time interval of the ODE model and . In the following examples, the bandwidth is equal to . To compute Hooker (2009)'s test , the demonstration code given in the supplementary materials of that paper was revised to make it feasible to handle the cases with unknown parameters.

Three cases are considered in this example. In each case, the linear null ODE model is added with one kind of disturbance term to form alternative ODE models. From the above ODE models, it can be seen that and correspond to the null hypothesis, otherwise to the alternative hypotheses. To transform the arbitrary length of sample time interval to 1, we induce a timescale parameter τ. The true parameter is set to be and the initial values are . We choose and . The empirical sizes and powers are reported in Table 1 , where ϵ and denote error terms with independent and highly correlated components, respectively. The rows with and show the values of our test and Hooker's test with a fixed design.

Table 1

Empirical sizes and powers in Example 1.

Model	α	0	0.5	0	0.5	1	0	1
	β	0	0	0.5	0.5	0	1	1
H11 with ϵ	TMn	0.043	1.000	0.819	1.000	1.000	1.000	1.000
	TH	0.054	1.000	1.000	1.000	1.000	1.000	1.000
H12 with ϵ	TMn	0.053	1.000	0.990	1.000	1.000	1.000	1.000
	TH	0.059	1.000	1.000	1.000	1.000	1.000	1.000
H13 with ϵ	TMn	0.063	1.000	1.000	1.000	1.000	1.000	1.000
	TH	0.056	1.000	1.000	1.000	1.000	1.000	1.000
H11 with ϵ⁎	TMn	0.051	1.000	0.935	1.000	1.000	1.000	1.000
	TH	0.071	1.000	1.000	1.000	1.000	1.000	1.000
H12 with ϵ⁎	TMn	0.051	1.000	0.998	1.000	1.000	1.000	1.000
	TH	0.063	1.000	1.000	1.000	1.000	1.000	1.000
H13 with ϵ⁎	TMn	0.069	1.000	1.000	1.000	1.000	1.000	1.000
	TH	0.072	1.000	1.000	1.000	1.000	1.000	1.000

Data sets are generated from the following ODE models:

With the two kinds of error terms, the results show that the trajectory matching-based test maintains the significance level when both and . It also has very good powers with all of the alternative models. The performances of are also good with the independent components of error terms. has good powers with dependent components of error terms, yet its empirical sizes are relatively larger than the significance level consistently under the three null models.

Under the null hypothesis, the model is called the Solow growth model, which is popular to describe long-run economic growth. In this model, and represent the labor and capital respectively. We choose the true parameters , , , . We report the empirical sizes and powers in Table 2 .

Table 2

Empirical sizes and powers in Example 2.

Model	α	0	0.5	0	0.5	1	0	1
	β	0	0	0.5	0.5	0	1	1
H2 with ϵ	TMn	0.043	0.491	0.776	0.946	1.000	1.000	1.000
	TH	0.051	1.000	1.000	1.000	1.000	1.000	1.000
H2 with ϵ⁎	TMn	0.042	0.609	0.898	0.838	1.000	1.000	1.000
	TH	0.100	1.000	1.000	1.000	1.000	1.000	1.000

Data sets are generated from the following ODE model:

For this nonlinear ODE model, it can be seen that also has good powers and sizes with a fixed design. has superior performances when the components of the error are independent. It shows larger powers than . However, in the dependent error components case, it cannot maintain the significance level. This is not surprising because requires the condition of independent components on the derivation of the null distribution.

The null model is the famous FitzHugh-Nagumo ODE model (FitzHugh (1961); Nagumo et al. (1962)) in neuroscience, which describes the behavior of spike potentials in the giant axon of squid neurons. Following Ding and Wu (2014), the true parameter is , , , and the initial values are . The empirical sizes and powers are reported in Table 3 .

Table 3

Empirical sizes and powers in Example 3.

Model	α	0	0.5	0	0.5	1	0	1
	β	0	0	0.5	0.5	0	1	1
H3 with ϵ	TMn	0.049	0.844	0.301	1.000	1.000	0.999	1.000
	TH	0.052	1.000	1.000	1.000	1.000	1.000	1.000
H3 with ϵ⁎	TMn	0.047	0.645	0.263	0.995	1.000	0.988	1.000
	TH	0.108	1.000	1.000	1.000	1.000	1.000	1.000

The data sets are generated from the following ODE model:

still works very well in all cases. Similar to Example 2, performs well in the independent error components case, but fails to work in the dependent error components case.

The null ODE is called the Lotka-Volterra model in ecology, which is designed for modeling the evolution of prey-predator populations (Lotka (1910); Volterra (1928); Goel et al. (1971)). Following Brunel (2008), we set the true parameters , , and . The empirical sizes and powers are summarized in Table 4 . This model favors which can have much higher power than , where both of them can hold the significance level well.

Table 4

Empirical sizes and powers in Example 4.

Hypothesis	α	0	0.5	0	0.5	1	0	1
	β	0	0	0.5	0.5	0	1	1
H4 with ϵ	TMn	0.042	0.194	0.908	1.000	0.951	1.000	1.000
	TH	0.040	0.998	1.000	1.000	1.000	1.000	1.000
H4 with ϵ⁎	TMn	0.051	0.248	0.869	1.000	0.944	1.000	1.000
	TH	0.060	1.000	1.000	1.000	1.000	1.000	1.000

Data sets are generated from the following ODE models:

The models are the same as in the previous three examples, except that the data of the second component are measured.

The empirical sizes and powers are presented in Table 5 . Since we only use the second component data, the independent and dependent error components cases are similar. The results show that can well control the size under all the null hypotheses. It also shows good powers under most of the alternatives. The test still has powers when and , meaning that as the unobserved component is correlated with the observed component, the test could still be sensitive in a certain extent. For , and , the test does not show powers when and because the magnitude of disturbance is not sufficiently large. In fixed design we set , the powers rise up to 0.508 and 0.101 for and with the independent components error respectively. also performs well and it has larger powers than .

Table 5

Empirical sizes and powers in Example 5.

Hypothesis	α	0	0.5	0	0.5	1	0	1
	β	0	0	0.5	0.5	0	1	1
H11 with ϵ	TMn	0.042	0.045	0.443	0.395	0.054	1.000	1.000
	TH	0.057	0.041	1.000	0.994	0.276	1.000	1.000
H12 with ϵ	TMn	0.045	0.059	0.896	0.913	0.053	1.000	1.000
	TH	0.046	0.088	1.000	1.000	0.339	1.000	1.000
H13 with ϵ	TMn	0.050	0.078	0.543	0.074	0.214	0.999	0.196
	TH	0.047	0.693	1.000	0.689	0.989	1.000	0.989
H2 with ϵ	TMn	0.048	0.059	0.592	0.645	0.066	1.000	1.000
	TH	0.055	0.204	1.000	1.000	0.646	1.000	1.000
H3 with ϵ	TMn	0.062	0.290	0.173	0.808	0.993	0.965	1.000
	TH	0.038	1.000	1.000	1.000	1.000	1.000	1.000
H4 with ϵ	TMn	0.068	0.136	0.807	0.997	0.940	1.000	1.000
	TH	0.047	0.994	1.000	1.000	1.000	1.000	1.000
H11 with ϵ⁎	TMn	0.064	0.061	0.382	0.338	0.051	1.000	1.000
	TH	0.050	0.054	0.999	0.995	0.253	1.000	1.000
H12 with ϵ⁎	TMn	0.051	0.059	0.855	0.886	0.045	1.000	1.000
	TH	0.053	0.093	1.000	1.000	0.305	1.000	1.000
H13 with ϵ⁎	TMn	0.054	0.072	0.506	0.060	0.196	0.998	0.181
	TH	0.046	0.666	1.000	0.676	0.983	1.000	0.990
H2 with ϵ⁎	TMn	0.044	0.054	0.526	0.808	0.072	1.000	1.000
	TH	0.047	0.200	1.000	1.000	0.616	1.000	1.000
H3 with ϵ⁎	TMn	0.052	0.238	0.176	0.779	0.993	0.946	1.000
	TH	0.052	1.000	1.000	1.000	1.000	1.000	1.000
H4 with ϵ⁎	TMn	0.054	0.136	0.756	0.999	0.918	1.000	1.000
	TH	0.049	0.988	1.000	1.000	1.000	1.000	1.000

We conclude that has well-controlled sizes and good powers for ODE systems with fixed design. It is suitable for fully observed as well as partially observed ODE systems. has good performances with independent components error terms, showing higher powers than in some cases. However, as seen in Example 2, Example 3, cannot well maintain the significance level for some ODE models with dependent components error terms. Some unreported results also show its failure to maintain the significance level in heteroscedastic cases, while still works well in this aspect. This would be because needs more restrictive prerequisites to use than needs. Thus, is more robust against model settings than . Another advantage of is its computing stability. It avoids computing the determinants and eigenvalues of certain matrices the algorithm of computing needs, which could cause computational errors. This has also been shown in some unreported numerical studies. This advantage also means lower computing costs. For example, to complete 1000 experiments in case with ϵ and using MATLAB R2021a, computing spends 126.87 s, while computing only needs 17.10 s. Overall, these two tests could be complements to each other, and is of robustness property against model settings.

We have further comment on the test . The observation about the significance level maintenance problem of motivates using an adjusted value in the denominator where as such that is asymptotically equivalent to , and can bring down the size in finite sample cases. This is a frequently used approach in practice. However, we have tried this method and found that this strategy seems not to work for . In some cases with dependent components error terms or heteroscedasticity, we need a large to help maintain the significance level. In contrast, in some other cases, a large causes the test to be too conservative. Also, it is understandable that a large denominator also brings the test power down. It is hard to find a suitable denominator to use in all cases.

The application to COVID-19 epidemic models

The classic SEIR ODE model has been used for respectively describing the infected population in Japan and Algeria (see, Kuniya (2020), and Bentout et al. (2020)). We apply our test, because of its robustness against different models, to check whether the model is plausible.

The form of the mechanical ODE system is as follows: Here represent the susceptible, exposed, infective and removed populations respectively. We also add a timescale parameter τ to normalize the sample time interval. The ODE model (3) has been applied to fit reported infective cases data during a period starting from the date that the first case was reported. Specifically, we use the data from 15 January to 29 February in Japan and from 25 February to 18 April, 2020 in Algeria. Note that this is a partially observed ODE system since we only observed infective cases corresponding to . Following the instruction in Kuniya (2020) and Bentout et al. (2020), we fix and . Further, let , thus each population is actually the proportion to the total population. Also transfer the number of reported infective cases to the corresponding proportion by using , where is the total number of people and is the fraction of infective individuals that are identified by diagnosis. We then use the nonlinear least squares method to estimate . According to Kuniya (2020), almost does not affect the estimation. Here we suppose . As for , we have for Japan and for Algeria. The initial values are set to be . We use the Epanechnikov kernel and the bandwidth to ensure it is larger than the interval between two consecutive time points.

For the Japan epidemic data, the value of is 18.65 and the corresponding p-value is about 0. For the Algeria epidemic data, these values are 16.49 and 0 respectively. These results suggest that the SEIR model, which is the ODE model under the null, may not be plausible (Fig. 1, Fig. 2 ).

Fig. 1

The data of Japan infective cases: time course of responses and residuals.

Fig. 2

The data of Algeria infective cases: time course of responses and residuals.

Conclusion

This paper studies model checking for ordinary differential equations model with fixed design. The asymptotic properties of a test , under the null, global, and local alternative hypothesis, are presented. We give two new propositions about U-statistics with varying kernels for independent but non-identical data to derive these results. These propositions are useful to handle many test statistics that are U-statistics with non-i.i.d. data. As an important application, we use to check the SEIR models by using openly accessible data, finding that it may not be suitable for describing the COVID-19 infected size in Japan and Algeria in the period we described in the application.

As for checking epidemic ODE models, some critical issues are still unsolved. For instance, the parameters may be estimated by a Bayesian method instead of the frequentist method. This is because we may want to estimate them more accurately using the past data. This deserves a further study.