key: cord-0566415-b1mp3tn5
authors: Sung, Chih-Li; Hung, Ying
title: Efficient calibration for imperfect epidemic models with applications to the analysis of COVID-19
date: 2020-09-26
journal: nan
DOI: nan
sha: 9c2c125d9ff4dab6a520e0ff387070de4495e615
doc_id: 566415
cord_uid: b1mp3tn5

The estimation of unknown parameters in simulations, also known as calibration, is crucial for practical management of epidemics and prediction of pandemic risk. A simple yet widely used approach is to estimate the parameters by minimizing the sum of the squared distances between actual observations and simulation outputs. It is shown in this paper that this method is inefficient, particularly when the epidemic models are developed based on certain simplifications of reality, also known as imperfect models which are commonly used in practice. To address this issue, a new estimator is introduced that is asymptotically consistent, has a smaller estimation variance than the least squares estimator, and achieves the semiparametric efficiency. Numerical studies are performed to examine the finite sample performance. The proposed method is applied to the analysis of the COVID-19 pandemic for 20 countries based on the SEIR (Susceptible-Exposed-Infectious-Recovered) model with both deterministic and stochastic simulations. The estimation of the parameters, including the basic reproduction number and the average incubation period, reveal the risk of disease outbreaks in each country and provide insights to the design of public health interventions.

The estimation of unknown parameters in simulations, also known as calibration, is crucial for practical management of epidemics and prediction of pandemic risk. A simple yet widely used approach is to estimate the parameters by minimizing the sum of the squared distances between actual observations and simulation outputs. It is shown in this paper that this method is inefficient, particularly when the epidemic models are developed based on certain simplifications of reality, also known as imperfect models which are commonly used in practice. To address this issue, a new estimator is introduced that is asymptotically consistent, has a smaller estimation variance than the least squares esti- 

The coronavirus disease pandemic has shown profound impacts on public health and the economy worldwide. The development of efficient and effective public health interventions to prevent major outbreaks and contain the pandemic relies heavily on a quantitative understanding regarding the spread of the virus, such as the transmission rate and the average incubation period. A commonly used approach in epidemiology is to estimate these quantities of interest using epidemic mathematical models, such as the susceptible-infected recovered (SIR) model, with agentbased simulations which capture complex social networks and global scale into the models [1, 2, 3] .

To estimate the parameters of interest, a widely used frequentist approach is to minimize the sum of the squared distances between the observed data and the simulation outputs, which is often referred to as the least squares approach. See, for example, [4, 5, 6, 7, 8, 9, 10, 11] . This estimation approach is intuitive and easy to compute; however, it is shown in this paper that this method is inefficient, that is, its asymptotic variance is not theoretically minimal, particularly when the mathematical models associated with the simulators are built under certain assumptions or simplifications, which may not hold in reality. These simulators are called imperfect simulators in the computer experiment literature [12, 13, 14] . Imperfect simulators are common in epidemiology [2] , and therefore the parameter estimates based on the least squares approach are not efficient.

To improve the estimation efficiency, in this paper we propose a new estimation method for the parameters of interest in the epidemic models. In the computer experiment literature, these unknown parameters associated with the mathematical models are often called calibration parameters, and the process of estimating the parameters such that the model simulations agree with the observed data is called calibration [12, 15] . Although there are numerous developments on calibration, most of the work focus on continuous outputs while the discussions on non-Gaussian outputs, such as count data which are often observed in epidemiology, are scarce [16, 17] . In this paper, we propose a new estimation method for non-Gaussian outputs, particularly for count data for our applications in epidemiology, which minimizes the L 2 projection of the discrepancy between the true mean process and the simulation outputs. It can be shown that the proposed estimator is asymptotically consistent, and provides a smaller asymptotic variance than the least squares estimator. Furthermore, it can be shown that the proposed estimator achieves the semiparametric efficiency, even when the model simulations cannot match the reality due to certain assumptions or simplifications.

It is worth noting that there are extensive studies and applications of calibration by Bayesian procedures [18, 19, 20, 21] . However, without taking the model imperfection into account in the conventional Bayesian framework, the theoretical justification for the parameter estimation with imperfect simulators are not fully developed. On the other hand, Bayesian calibration of [12] takes into account the model imperfection through Gaussian process modeling, but it suffers from the unidentifiability issue when the parameter estimation is of interest [22, 23, 24, 25, 26] . Furthermore, most of the existing developments are based on continuous outputs with a Gaussian assumption, which is not valid for the count data in the epidemic models in our applications. Recent studies on addressing the unidentifiability issue can be found in [14] and [27] .

The remainder of the paper is organized as follows. In Section 2, the new estimation method based on L 2 projection is introduced. Theoretical properties of the proposed estimator are developed in Section 3. In Section 4, numerical studies are conducted to demonstrate the finite sample performance of the proposed estimator and the empirical comparison with the least squares estimator. In section 5, the estimation method is applied to the study of COVID-19 based on an SEIR (Susceptible-Exposed-Infectious-Recovered) compartmental model with both deterministic and stochastic simulations. Discussions and concluding remarks are given in Section 6. Computational details for the estimation are given in Appendix, and the mathematical proofs and the R [28] code for implementation are provided in Supporting Web Materials.

Compartmental models are widely used mathematical models for the study of infectious diseases, in which the population is assigned to compartments with labels such as S (Susceptible), I (Infectious), and R (Recovered and immune) [18] . For example, the most basic compartment model in epidemiology is the SIR model, in which each living individual is assumed to be in one of the above three compartments at any given time. The model specifies the rates at which individuals change their compartment, which include the progresses from S to I when infected, and from I to R upon recovery. Various extensions to SIR have been extensively developed in the literature of epidemiology [2] .

Let f (x, θ) denote the number of infected cases at time x ∈ Ω ⊆ + , where θ ∈ Θ ⊆ q is a set of unknown calibration parameters associated with the compartmental model. For instance, an SIR model contains two parameters, θ = (β , γ), where β is the infection rate controlling the transition from S to I , and γ is the recovery rate controlling the transition from I and R . Suppose that y i is the reported number of infected cases at time x i . Then, given the reported number of infected cases in n days, { (x i , y i ) } n i =1 , the commonly used approach to estimate the parameters is to minimize the sum of squared differences between actual numbers of infected cases and simulation outputs from compartmental models. The estimated parameters are denoted byθ LS n , where LS stands for least squares, and they are obtained byθ LS n = arg min

In addition to the least squares estimator, the maximum likelihood estimator (MLE) is also a commonly used estimation approach. Assume that y i ∼ Poi(f (x i , θ)), where i = 1, . . . , n, we obtain the MLE for the calibration

Despite the wide applications of the least squares approach and MLE, it can be shown that the least squares estimator does not achieve the semiparametric efficiency when the simulator f (x, θ) is imperfect, meaning that the simulation output cannot perfectly fit the response, even with the best fit of θ. The asymptotic variance can be reduced by the proposed estimator introduced in this subsection. It can also be shown that MLE is asymptotically inconsistent when the simulator f (x, θ) is imperfect. Theoretical justifications are provided in Section 3.

Assume that the number of cases y i follows a Poisson distribution: y i ∼ Poi(λ (x i )) for i = 1, . . . , n, and y i and y j are mutually independent for any i j , where λ (x i ) is the true mean function of y i . The function λ (x ) is often called the true process in the computer experiment literature [12, 13, 29] . Ideally, if the underlying mean function λ (x ) is known, the true parameter can be defined as the minimizer of the L 2 projection of the discrepancy between the true process and the simulation output, that is,

where

In reality, the underlying true process λ ( ·) is unknown that needs to be estimated by observed data. Therefore,

given the data { (x i , y i ) } n i =1 , we propose to estimate the true process by the kernel Poisson regression [30, 31] . Similar to the conventional Poisson regression [32] , we use the logarithm as the canonical link function, that is, logλ n ( ·) = ξ n ( ·), andξ n ( ·) is fitted byξ

where · 2 N Φ (Ω) is the norm of the reproducing kernel Hilbert space generated by a given positive definite reproducing kernel Φ, and κ n is a tuning parameter, which can be chosen by cross-validation methods. Thus, the proposed estimator of θ, which we call L 2 -estimator throughout this paper, is the minimizer of the L 2 projection as follows:

The optimal solution of (4) has the form ofξ n (x ) =b + n i =1â i Φ (x i , x ), whereb and {â i } n i =1 can be obtained by the iterative re-weighted least squares algorithm [33, 34, 35] . The detail of the algorithm is given in Appendix A. In practice, the calculation of the L 2 norm in (5) can be approximated by numerical integration methods, such as Monte Carlo integration [36] .

For some stochastic agent-based simulations in epidemiology that can be computationally intensive, it is infeasible to obtain f (x, θ) by conducting simulations for all possible combinations of the input parameters. In such cases, the simulator is often approximated by a emulator which is computationally efficient. There are extensive studies on the development of statistical emulators in the computer experiment literature [15] . Gaussian processes (GPs) are the most commonly used tools in the construction of emulators [37] . Based on computer experiments with sample size N , a statistical emulator is denoted byf N (x, θ). which produces a predictive distribution of f (x, θ) with any untried (x, θ) ∈ (Ω, Θ). Specifically, the distribution off N (x, θ) with any untried (x, θ) ∈ (Ω, Θ) is a normal distribution with the mean function, defined by m N (x, θ) , and the variance function, defined by v 2 N (x, θ). We refer more details to [37] . Thus, by Fubini's Theorem, the L 2 -estimator of (5) can be replaced bỹ

The applications of the proposed method with various existing emulators are demonstrated in Sections 4 and 5.

It is worth noting that, the Poisson regression, y i ∼ Poi(λ (x i )), may encounter overdispersion due to the presence of greater variability [32] . That is, the variance of the data is larger than the mean, which violates the assumption of Poisson distribution. The deviance goodness of fit test [32] can be used to assess the model assumption. To take into account the issue of overdispersion, a quasi-Poisson regression can be considered which assumes the variance of y i is φλ (x ), where φ > 1 is the overdispersion parameter. The overdispersion parameter can be estimated by the ratio of the deviance to the effective degree freedom. The details of the deviance goodness of fit test and the estimation of overdispersion parameter are provided in Appendix A.

Theoretical properties of the L 2 -estimator are discussed in this section, including the asymptotic consistency and the semiparametric efficiency. Theoretical comparisons with the least squares estimators are also provided by examining their asymptotic variances. The proofs are given in Supporting Web Materials.

The following theorem shows that the L 2 -estimatorθ n in (5) is asymptotically consistent and normally distributed.

Remark When the overdispersion parameter, φ, is present in the Poisson regression, the result in Theorem 1 can be rewritten as

By the delta method, the following corollary extends the result of Theorem 1 to a function of the L 2 -estimator, which we denote as g (θ n ).

For a function g satisfying the property that g (θ * ) exists and is non-zero valued, we have

as n → ∞.

Corollary 2 provides a theoretical support for the estimation and inference of some commonly used quantities of interest in epidemiology, such as the basic reproduction rate, which measures the transmission potential of a disease.

For instance, in an SIR model the basic reproduction rate is a ratio of two calibration parameters, that is, g (θ 1 , θ 2 ) = θ 1 /θ 2 , where θ 1 = β and θ 2 = γ. The result of Corollary 2 can then be applied to construct the confidence intervals for the basic reproduction rate.

When estimating the unknown parameters in compartmental models, the parameter of interest θ * in (3) is qdimensional, while the parameter space of the Poisson model, y i ∼ Poi(λ (x i )), contains an infinite dimensional function space that covers λ. Therefore, the calibration problem is regarded as a semiparametric problem. For these problems, the estimation method that can reach the highest estimation efficiency is called semiparametric efficient [38, 39] . Specifically, let Λ be an infinite dimensional parameter space whose true value is λ 0 . Suppose that T n is an estimator for θ * and √ n (T n − θ * ) is asymptotically normal. Let Λ 0 be an arbitrary finite dimensional space of Λ that satisfies λ 0 ∈ Λ 0 . Consider the same calibration problem but with the parameter space Λ 0 , then under this parametric assumption and some regularity conditions, an efficient estimator can be obtained by the maximum likelihood method, which is denoted by S Λ 0 n . Then, the estimator T n is called semiparametric efficient if there exists a Λ 0 such that S Λ 0 n has the same asymptotic variance as T n . More details regarding the semiparametric efficiency can be found in [13, 38, 39] .

It can be shown in the following theorem that the proposed L 2 -estimator is semiparametric efficient.

Under the regularity conditions in Theorem 1,θ n is semiparametric efficient.

When the simulator f is too costly to evaluate, as discussed in Section 2.2, an statistical emulator can be considered after conducting a computer experiment of size N on the simulator. Suppose that the emulator of f (x, θ), i.e., f N (x, θ), follows a normal distribution with the mean function, m N (x, θ), and the variance function, v 2 N (x, θ), i.e.,

and the L 2 -estimator is obtained by (6) asθ n . Then, the following theorem provides the asymptotic distribution ofθ n .

With the emulator (8) , it is of no surprise that the estimatorθ n is asymptotic inconsistent. However, when the size of the computer experiment, N , is sufficiently large, with an appropriate emulator (e.g., GP emulator) and under some regularity conditions, we have

, leading to θ N → θ * , which implies thatθ n is asymptotic inconsistent when N is sufficiently large.

In the next theorem, the asymptotic properties of the least squares estimator are developed and compared with those of the L 2 -estimator.

Similar to the L 2 -estimator, it is shown that the least squares estimator is asymptotically consistent and normally distributed. It can also be shown that W 2 (θ * ) ≥ W 0 (θ * ), which leads to

This implies that the asymptotic variance of the least squares estimatorθ LS n is greater or equal to that ofθ n . The equality in (9) holds if and only if

In the next theorem, the asymptotic properties of the MLE as in (2) are developed and compared with those of the L 2 -estimator.

Unlike the L 2 -estimator and the least squares estimator, the MLE is asymptotically inconsistent. For example, 

In this section, two artificial examples are conducted to examine the finite sample performance of the proposed method and compare the estimation performance with the least squares approach. These numerical studies are conducted on a desktop with 3.5 GHz CPU and 8GB of RAM, and 4 CPUs are available for parallel computing.

We consider an imperfect simulator adapted from [13] with one calibration parameter. The true process is assumed

, and it is illustrated in the left panel of Figure 1 as the solid line. The data are generated from equal-spaced inputs in [0, 2π ] with n = 50 and the outputs are generated from a Poisson distribution with the mean process λ (x i ) for i = 1, . . . , 50, which are shown as the solid dots in the left panel of Figure 1 .

We assume that the simulation output is

This simulator is imperfect because √ θ 2 − θ + 1 is always positive for any θ ∈ Θ. The true parameter can be analytically solved by minimizing (3), which gives that θ * = −0.1789. Plugging in the true calibration parameter, the simulator

is demonstrated as the dashed line, which is imperfect because, even with the true minimizer, the discrepancy between the simulation output and the true process is nonzero. The performance of the L 2 -estimator is compared with the least squares estimator and maximum likelihood estimator based on the mean squared errors (MSEs) obtained from 100 replicates, that is, 100 i =1 (θ i − θ * ) 2 /100, whereθ i is the estimate at the i -th replicate. Their MSEs are shown in the first three bars in the right panel of Figure 1 . It shows that the L 2 -estimator ("L2") yields a smaller MSE than the least squares estimator ("LSE") and maximum likelihood estimator ("MLE"). To quantify the uncertainty of the L 2 estimator, the 95% confidence intervals are constructed based on the asymptotic result in Theorem 1, where λ, θ * and ¾ are approximated byλ,θ n , and Monte-Carlo integration [36] , respectively. Out of the 100 replicates, the true parameter is contained by the confidence interval 96 times, which appears to be close to the nominal coverage 95%.

We further compare the estimation performance for the cases when the simulations are computationally demanding and therefore statistical emulators are built as surrogates. Before comparing the estimation performance, we first examine the emulation performance of two existing emulation methods that are applicable to count data, which are the multiresolution functional ANOVA emulation [41] and the heteroscedastic Gaussian process emulation [42] . Both methods have available packages in R [28] , which are MRFA [43] and hetGP [44] , respectively. These emulators are trained by conducting a computer experiment, which simulates the model outputs of f (x, θ) of size m, where the inputs are sampled from (x, θ) ∈ (Ω, Θ) ∈ 2 using a Latin hypercube design (LHD) [45] . For each input setting, simu-lations are conducted with a replicates. The emulation performance is examined by performing predictions on 10, 000 random untried input settings from (Ω, Θ). With four different combinations of m and a, the root mean squared prediction errors (RMSPEs) of the two emulators along with their computational time are reported in Table 1 of Appendix B. In this example, it appears that hetGP outperforms MRFA in terms of computational time and RMSPE. Thus, we select the emulator built by hetGP as the surrogate to the actual simulator in the following analysis.

Next, we compare the estimation performance with the hetGP emulator built by m = 100, a = 100 samples, leading to total sample size N = ma = 10, 000. The L 2 estimator is obtained by (6) with the emulator, and the least squares estimator is similarly obtained by minimizing

. For the MLE as in (2), the actual simulator f (x, θ) is replaced by the mean of the hetGP emulator, i.e., m N (x, θ). The MSEs are shown in the last three bars in the right panel of Figure 1 . Similar to the previous result without emulators, the L 2 -estimator provides a smaller MSE than the least squares estimator and MLE. By comparing the first three and last three bars, it is not surprising to see that the MSEs of "L2+emulator", "LSE+emulator", and "MLE+emulator" are larger than "L2", "LSE", and "MLE" due to the prediction uncertainty from emulation. It is also worth noting that the fourth bar is lower than the second and the third, which indicates that, by choosing the proposed estimator, the MSE can be smaller than the least squares estimator and MLE even if the simulations are approximated by an emulator. Similarly, we construct the 95% confidence intervals based on the asymptotic result in Theorem 4 for the L 2 estimator of (6), and out of the 100 replicates, the true parameter is contained by the confidence interval 95 times, which matches the nominal coverage of 95%.

We consider a more complex problem with three calibration parameters adapted from [14] . Assume that the true mean The results are shown in the first three bars in each plot of Figure 2 , in which the y -axis represents the MSEs. Similar to the previous example, it appears that the L 2 -estimator outperforms the least squares estimator and MLE for all of the three parameters.

In this example, we also examine the prediction performance of the two existing emulators, MRFA and hetGP.

A computer experiment is conducted to train the two emulators by running the simulation outputs of f (x, θ) at m unique sample locations with a replicates, in which the unique input locations are sampled from (x, θ) ∈ (χ, Θ) ⊆ 4 using an LHD. After the emulators are built, the RMSEs are computed based on the predictions of 10, 000 untried input locations, and the prediction results are summarized in Table 2 with different settings of m and a. Similar to the previous example, the hetGP method outperforms MRFA in terms of prediction accuracy and computational time.

With a larger a, i.e., more replicates, the prediction accuracy of hetGP appears to increase without much increase in computational time. Thus, we select hetGP as the emulator in the following analysis.

We now compare the estimation performance for the cases where emulators are constructed as surrogates to the actual simulations. The emulator is built by hetGP with m = 300, a = 100 and based on the emulator, the estimation performance is summarized by the last three bars in each of the three plots in Figure 2 . The results indicate that, either when the actual simulator is conducted or emulated, the L 2 -estimator provides smaller MSEs compared to other two estimators. 

The proposed method is applied to estimate unknown parameters in compartmental models for the analysis of COVID-19 pandemic. According to the epidemiology literature, the SEIR model, which consists of four compartments, Susceptible-Exposed-Infectious-Recovered, is widely recommended for COVID-19 simulations because it accounts for the incubation period through the exposed compartment [21, 46, 47, 48, 49] . We focus on two types of SEIR simulators: a deterministic simulator and a stochastic simulator, which are discussed in Sections 5.1 and 5.2, respectively. To estimate the unknown parameters in these two SEIR simulators, we collect the actual numbers of infected cases from Johns Hopkins University CCSE repository [50] through an R package covid19.analytics [51] . 

The SEIR model extends the idea of SIR to further consider an E compartment before entering the I compartment, which indicates a significant incubation period during which individuals have been infected but are not yet infectious themselves. Mathematically, a deterministic SEIR model can be written as:

where S , E , I and R represent the numbers of cases in the corresponding compartment, N = S + E + I + R is the total population, x is time, β is the contact rate that represents the average number of contacts per person per time in the susceptible compartment, γ is the recovery rate from the infectious compartment, and κ is the incubation rate which represents the rate of latent individuals becoming infectious, or equivalently, the average incubation period is 1/κ. These ordinary differential equations have no analytical solution but can be solved by numerical solvers, such as the ODEPACK [52] . The basic reproduction number, R 0 = β /γ, is the expected number of new infected cases from an infectious individual in a population where all subjects are susceptible. It is of great interest in epidemiology because it provides the insight of the virus transmission, which is essential to the development of effective strategies to contain the pandemic. An accurate and efficient estimation of R 0 is not only important for the public safety, but it also has significant impacts on global economy.

This model includes six unknown parameters: β , κ, γ and the initial numbers of infectious, exposed, and recovered cases (denoted by I (0), E (0), and R (0) respectively), so we have θ = (β , κ, γ, I (0), E (0), R (0)). The number of the daily infected cases based on the SEIR model can be obtained by f (x, θ) = κE (x ), where E (x ) is the solution of E in the ordinary differential equations of (12) . Before estimating the parameters, a deviance goodness of fit test is performed to examine the kernel Poisson regression as in (4), i.e., y i ∼ Poi(λ n (x i )). It appears that the p-values of the test are all smaller than 0.0001, which indicates that there is a lack-of-fit in the current model. Therefore, a more flexible model, the quasi-Poisson as described in 2.2, is applied to capture the potential overdispersion.

For each country, the L 2 -estimator of θ is obtained by minimizing (5), and the corresponding estimated reproduction number R 0 can be calculated by R 0 = β /γ. The point estimates of R 0 and their 95% confidence intervals, which are obtained by the result of Corollary 2, are summarized in Figure 3 for the 20 countries. It shows that, from March, 2020 to March, 2021, all of the 20 countries have the basic reproduction numbers greater than 1, which means that the COVID-19 outbreak still post threats to these countries. Note that, the recovery rate γ is in the denominator of R 0 , and therefore the variation of R 0 appears to be higher for the countries having smaller recovery rates. Figure 4 . Note that the confidence intervals are similarly constructed based on Corollary 2. That is, the variance of f (x,θ n ) can be approximated by

where θ is the partial derivative with respect to θ. In general, it appears that the simulation results can reasonably capture the overall trend observed from the actual numbers of infected cases, which are shown as the gray dots.

For Iran, Czechia, and Spain, the discrepancy between the simulation results and actual observations is relatively larger than the other countries. This is partly because SEIR is an imperfect simulator which is built based on some assumptions or simplifications, and these assumptions may have larger deviations from the reality for certain countries.

Another reason is that the intrinsic dynamics are neglected in the deterministic simulations. To take into account the dynamics, a stochastic simulator is considered in the next subsection. 

A stochastic SEIR simulation provides a more sophisticated and realistic framework to integrate infection dynamics in different compartments as continuous-time Markov chains [53, 54, 55] . To conduct these simulations, we implement an R package, SimInf [56] , in which the simulation results are obtained by the Gillespie stochastic algorithm [57] .

Stochastic SEIR simulations are computationally more demanding. For example, it takes more than 10 minutes to produce one simulation result for one country under a given parameter setting. It is computationally infeasible to perform simulations for all the possible combinations of the parameters; therefore, an emulator is constructed as an efficient surrogate to the actual simulation. In this study, we consider the hetGP emulator, which is built based on the simulations generated using a 60-run LHD for parameter settings with 20 equal-spaced time steps in x , which leads to the total sample size of m = 1200. For each parameter-input setting, 50 replicates are simulated, i.e., a = 50, so the total sample size of this computer experiment is N = ma = 60, 000. Based on this emulator, it takes less than two seconds to emulate the result for an untried parameter setting, which is significantly faster than the actual stochastic simulation.

With the hetGP emulator, which has the form of (8), the L 2 -estimators are obtained by minimizing (6) . The corresponding estimates of R 0 and their 95% confidence intervals are summarized in Figure 5 , where the variance is obtained based on the result of Theorem 4. It appears that Colombia, Brazil, and India have their basic reproduction numbers controlled below 1, which also show small basic reproduction numbers in the deterministic simulations (less than 1.05). We further report the estimated incubation period, 1/κ, for each country and the corresponding 95% confidence intervals in Figure 6 . The overall average incubation period is 5.15 as indicated by the red dashed line.

In Figure 7 , the actual numbers of infected cases are illustrated as the gray dots. By plugging in the L 2 -estimators, the simulation results for the top-12 countries with the highest R 0 are illustrated as the red curves, along with the 95% confidence intervals as the red dashed lines. Overall, the simulation results show a much better agreement with the actual observations compared to the deterministic ones in Section 5.1. In particular, by taking into account the intrinsic dynamics, the simulation discrepancy for Czechia is significantly reduced from the deterministic one shown in Figure   4 . Note that the confidence intervals are computed based on

which can be approximated by

using the result of Corollary 2. The first term, v 2 N (x,θ n ), distinguishes the variance from that of the deterministic SEIR (see (13) ), which contains the errors sourced from the emulation uncertainty and the variations under replication. As a result, the confidence intervals are generally wider than the ones of deterministic SEIR in Section 5.1.

A new calibration method is proposed to estimate the unknown parameters in epidemic models. The proposed estimator outperforms the least squares estimator by providing a smaller estimation variance and achieving the semiparametric efficiency. The proposed method is applied to the SEIR model for the analysis of COVID-19 pandemic. The estimates of the quantities of interest, such as the basic reproduction number and the average incubation period, and their confidence intervals are obtained based on the asymptotic results.

Apart from the frequentist approach studied in this paper, we are currently developing a Bayesian framework that extends the recent developments of Bayesian calibration to count data. For example, the orthogonal Gaussian process models [58] or the Bayesian projected calibration [27, 59] can be used to model the model discrepancy, which addresses the unidentifiability issue for continuous outputs, and it is conceivable to further extend the modeling to count data by incorporating the idea of the generalized calibration in [17] . It is also worth investigating the confidence set on the calibration parameters using the method of [60] for the application herein. Another interesting direction that deserves further studies is to relax the constant parameter assumption. Instead, the calibration parameters can be assumed to be functions of some factors, such as time or temperature, which not only increases the model flexibility but also can provide further insights to the time-course dynamics of the COVID-19 infection. 

Since the optimal solution has the form of ξ n (

, one can show that the penalized likelihood in (4) can be rewritten as

where a = (a 1 , . . . , a n ), ψ (x ) = (Φ (x, x 1 ) , . . . , Φ(x, x n )), and Φ = (Φ (x i , x j )) 1≤i ,j ≤n . The optimal solution of a and b can then be obtained by taking the first-order partial derivatives of the objective function with respect to a and b and setting them equal to zero, which can be solved by the iterative re-weighted least squares algorithm as follows.

where 1 n = [1, . . . , 1] T and 0 n = [0, . . . , 0] T , and denote W as an n × n diagonal matrix with diagonal elements

Then, in each step, one first solve for β :

with an initial guess of η, which is a vector of size n, and then update each element of η by

The estimateβ can then be obtained by continuing solving for β and η iteratively until some convergence criterion is met.

To examine the goodness-of-fit of the Poisson regression, the following deviance goodness of fit test is considered.

Since it can be shown that the deviance of the model follows a chi-square distribution asymptotically, that is

when n is sufficiently large, where the effective degree freedom, edf = trace(S), where

If the test indicates that overdispersion is present in the Poisson model, the overdispersion parameter φ can be estimated byφ = D /edf. , where m is the sample size of unique locations and a is the number of replicates. RMSPEs are reported for the two emulators based on 10, 000 random predictive locations.

Web Appendix A | ASYMPTOTIC RESULTS OFξ n OF (3)

We start with a result developed by [7] for general nonparametric regression (Lemma 11.4 and 11.5 in [7] ). Denote

Suppose F is the class of all regression functions equipped with the Sobolev norm · H m (Ω) , which is defined by

. Letξ n := arg min

for some κ n > 0. Then the convergence rate ofξ n is given in the following lemma.

Let ξ 0 ∈ F. Assume that there exists some nonnegative k 0 and k 1 so that

where z 0 = ξ 0 (x ) and x ∈ Ω. For κ −1 n = O (n 2m/(2m+1) ), we have

In fact, the norms of some reproducing kernel Hilbert spaces (RKHS) are equivalent to Sobolev norms. For instance, the RKHS generated by the Matérn kernel function, given by

where ν ≥ 1 and ρ ∈ + are tuning parameters and K ν is a Bessel function with parameter ν, is equal to the (fractional)

Sobolev space H ν+d /2 (Ω), and the corresponding norms · N Φ (Ω) and · H ν+d /2 (Ω) are equivalent [10, 6] . Therefore, as a consequence of Lemma 7, we have the following proposition forξ n obtained by (3) .

Suppose that ξ 0 ∈ F = N Φ (Ω), and N Φ (Ω) can be embedded into H m (Ω). Then, for κ −1 n = O (n 2m/(2m+d ) ), the estimatorξ n in (3) andλ n = exp{ξ n } satisfy

Proposition 8 suggests that one may choose κ n n −2m/(2m+d ) to obtain the best convergence rate λ n − λ 0 L 2 (Ω) = O p (n −m/(2m+d ) ), where a n b n denotes that the two positive sequences a n and b n have the same order of magnitude.

The regularity conditions for Theorem 1 are given below. For any θ ∈ Θ ⊂ q , write θ = (θ 1 , . . . , θ q ). Denote e i = y i − λ (x i ) and ξ ( ·) = log λ ( ·). C2: θ * is the unique solution to (2) and is an interior point of Θ.

C3: sup θ∈Θ f ( ·, θ) L 2 (Ω) < +∞.

for all θ ∈ U and all i , j = 1, . . . , q .

C6: W 0 (θ * ) = ¾ λ (X ) ∂f ∂θ (X , θ * ) ∂f ∂θ T (X , θ * ) is positive definite. C7: ξ ∈ N Φ (Ω) and N Φ (Ω, ρ) is Donsker for all ρ > 0.

The Donsker property is an important concept in the theoretical studies of empirical processes. The definition and detailed discussion are referred to [9] and [2] . [5] showed that if the conditions of Proposition 8 hold and m > d /2, then N Φ (Ω, ρ) is a Donsker. The authors also mentioned that under the condition C1 and ¾[exp{C |y i −λ (x i ) | }] < +∞ for some C > 0, the conditions of Proposition 8 are satisfied. Therefore, by choosing a suitable sequence of κ n , say κ n n −2m/(2m+d ) , one can show that condition C10 holds and C8 and C9 are ensured by Proposition 8.

Additional regularity conditions for Theorem 4 are given below, where we assume that the emulator,f N , is built based on a computer experiment with N samples, which follows the distribution as in (8) .

C11: θ N is the unique solution to (9) and is an interior point of Θ.

C12: sup θ∈Θ m N ( ·, θ) L 2 (Ω) < +∞.

for all θ ∈ U and all i , j = 1, . . . , q .

Additional regularity conditions for Theorem 5 are given below.

Suppose that there exists a neighborhood U of θ * such that f (x, ·) ∈ C 2,1 (U ) for all x ∈ Ω.

Additional regularity conditions for Theorem 6 are given below.

C18: θ is the unique solution to (12) and is an interior point of Θ.

Suppose that there exists a neighborhood U of θ such that f (x, ·) ∈ C 2,1 (U ) for all x ∈ Ω.

C20: Suppose that f (x, θ ) > 0 for all x ∈ Ω.

is positive definite.

We first show that under the regularity conditions C1-C10 in Web Appendix B,

The proof is developed along the lines described in Theorem 1 of [5] and Theorem 3.1 of [4] . We first prove the consistency,θ n p − → θ * . It suffices to prove that λ n ( ·) − f ( ·, θ) L 2 (Ω) converges to λ ( ·) − f ( ·, θ) L 2 (Ω) uniformly with respect to θ ∈ Θ in probability, which is ensured by

where the inequality follow from the Schwarz inequality and the triangle inequality.

Since it can be shown that g L 2 (Ω) ≤ Vol(Ω) g L∞ (Ω) , for any g ∈ L ∞ (Ω), where Vol(Ω) is the volume of Ω,

we have

By (S.4) and the conditions C3, C8, and C9, we have that (S.3) converges to 0 uniformly with respect to θ ∈ Θ, which proves the consistency ofθ n .

Sinceθ n minimizes (5), by invoking conditions C1, C2 and C5, we have

which implies that (3), we know thatξ n minimizes l over N Φ (Ω). Sincê θ n p − → θ * and by condition C5, ∂f ∂θ j ( ·,θ n ) ∈ N Φ (Ω) for j = 1, . . . , q with sufficiently large n. Then we have

We first consider C n . Let ρ) is Donsker. Thus, by Theorem 2.10.6 in [9] , F 1 = {exp{g } − exp{ξ } : g ∈ N Φ (Ω, ρ) } is also Donsker because the exponential function is a Lipschitz function. By condition C5, the class F 2 = { ∂f ∂θ j ( ·,θ n ), θ ∈ U } is Donsker. Since both F 1 and F 2 are uniformly bounded, by Example 2.10.8 in [9] the product class F 1 × F 2 is also Donsker. Thus, the asymptotic equicontinuity property holds, which implies that for any > 0 there exists a δ > 0 such that lim sup n→∞ Pr sup

. See Theorem 2.4 of [3] . This implies that for any > 0 there exists a δ > 0 such that

Suppose ε > 0 is a fixed value. Condition C9 implies that there exists ρ 0 > 0 such that P r ( ξ N Φ > ρ 0 ) ≤ ε/3. In addition, choose δ 0 to be a possible value of δ which satisfies (S.7) with = ε/3 and ρ = ρ 0 . Condition C8 implies that

Then, for sufficiently large n, we have

The first and third inequalities follow from the definition ofξ • n , and the last inequality follows from (S.7). Thus, this

implies that E 1n (ξ n , θ) tends to zero in probability, which gives

which implies

Then, by substituting (S.5) to (S.8) and using condition C2, Taylor expansion can be applied to (S.8) at θ * , which leads

whereθ n lies betweenθ n and θ * . By the consistency ofθ n , we then haveθ n p − → θ * , which implies that

Next, we consider D n . Define the empirical process

U } is a Donsker class, which ensures that E 2n ( ·) converges weakly in L ∞ (U ) to a tight Gaussian process, denoted by G ( ·). Without loss of generality, we assume G ( ·) has continuous sample paths. Then, by the continuous mapping theorem [8] and the consistency ofθ n , we have E 2n (θ n ) p − → G (θ * ). Because E 2n (θ * ) = 0 for all n, we have G (θ * ) = 0.

Then, we have E 2n (θ) p − → 0, which gives

Lastly, we consider E n . Applying conditions C5, C9, C10, we have

By combining (S.6), (S.9), (S.10) and (S.11), we havê

Thus, given condition C6 and

(S.12) implies the asymptotic result of Theorem 1 by the central limit theorem, that is,

If suffices to show thatθ n given in (4) has the same asymptotic variance as the estimator obtained by using maximum likelihood (ML) method. Consider the following q -dimensional parametric model indexed by γ, ξ γ ( ·) = ξ ( ·) + γ T ∂f ∂θ ( ·, θ * ), (S. 13) with γ ∈ q . By the fact that ξ ( ·) = log λ ( ·) and (S.13), it becomes a Poisson regression model with the coefficient γ.

The log-likelihood of γ is

and its first and second derivatives are

and

exp ξ ( ·) + γ T ∂f ∂θ ( ·, θ * ) ∂f ∂θ ( ·, θ * ) ∂f ∂θ T ( ·, θ * ) .

Since the true model is y i ∼ Poi(λ (x i )), the true value of γ is 0. Hence, under the regularity conditions of Theorem 1, the ML estimator of γ has the asymptotic expression

where W 0 is defined in (7). Then, a natural estimator for θ * in (2) iŝ

The asymptotic variance ofθγ n can be obtained by the delta method. First, define

for all t near zero, and let

Then, (S.15) implies that φ (θ (t ), t ) = 0 for all t near zero. Then, by the implicit function theorem, we have

By the delta method, we haveθ

and together with (S.14) and (S.16), it yieldŝ

Thus, the asymptotic variance of the ML estimator is 4V 0 (θ * ) −1 W 0 (θ * )V 0 (θ * ) −1 , which is the same as that ofθ n (see Theorem 1) . Therefore, the L 2 estimatorθ n in (4) is semiparametric efficient.

Similar to Web Appendix C, we first show that under the regularity conditions C1, C7-C15 in Web Appendix B,

We first prove the consistency,θ n p − → θ N . It suffices to prove that ¾ λ n ( ·) −f N ( ·, θ) 2 L 2 (Ω) converges to ¾ λ ( ·) − f N ( ·, θ) 2 L 2 (Ω) uniformly with respect to θ ∈ Θ in probability, which is ensured by

where the inequality follow from the Schwarz inequality and the triangle inequality. Then, by (S.4) and the conditions C8, C9, and C12, we have that (S.17) converges to 0 uniformly with respect to θ ∈ Θ, which proves thatθ n p − → θ N .

Sinceθ n minimizes (6), by invoking conditions C1, C11 and C14, we have

which implies that ( ·,θ n ) ∈ N Φ (Ω) for j = 1, . . . , q with sufficiently large n. Similar to the proof in Web Appendix C, we have 0 = ∂ ∂t l (ξ n ( ·) + t ∂m N ∂θ j ( ·,θ n )) | t =0 = 1 n n i =1 y i − exp{ξ n (x i ) } ∂m N ∂θ j (x i ,θ n ) + 2κ n <ξ n , ∂m N ∂θ j ( ·,θ n ) > N Φ (Ω) ∫ Ω ∂ 2 ∂θ i ∂θ j [m N (z ,θ n ) − λ (z ) ] 2 + v 2 N (z,θ n ) dz (θ n − θ N ) + o p (n −1/2 ), whereθ n lies betweenθ n and θ * . By the consistency ofθ n , we then haveθ n p − → θ N , which implies that

Thus, we have Thus, given condition C15 and

(S.12) implies the asymptotic result of Theorem 1 by the central limit theorem, that is,

Consider a parametric model, y i ind.

∼ N (f (x i , θ), 1), which is misspecified because the true model is y i ind.

∼ Poi(λ (x i )).

It is easy to see thatθ LS n is the maximum likelihood estimator (MLE) of θ under this misspecified model. Thus, under the regularity conditions C1-C4 and C16, the MLE converges to θ in probability which uniquely minimizes Kullback-Liebler divergence [1, 11] . That is,

where h 0 and h 1 are the density functions of y i under the misspecified model and the true model, respectively, and ¾ 0 denotes the expectation with respect to y i ∼ h 0 . Since the true model follows a Poisson distribution, we have

By the law of large numbers, we have

Therefore,θ LS n p → θ = arg min θ∈Θ ¾[λ (X ) − f (X , θ) ] 2 , which shows thatθ LS n is consistent because θ = θ * .

In addition, by [11] , the asymptotic normality ofθ LS n can be shown as follows,

= ¾ λ (X ) ∂f ∂θ (X , θ * ) ∂f ∂θ T (X , θ * ) + ¾ (λ (X ) − f (X , θ * )) 2 ∂f ∂θ (X , θ * ) ∂f ∂θ T (X , θ * ) = W 2 (θ * ). 

Suppose thatθ MLE n is the maximum likelihood estimator based the parametric model, y i ind.

∼ Poi(f (x i , θ)), which is misspecified because the true model is y i ind.

∼ Poi(λ (x i )). Thus, under the regularity conditions C1, C18-19, the MLE converges to θ in probability which uniquely minimizes Kullback-Liebler divergence [1, 11] . That is, 

The spread of awareness and its impact on epidemic outbreaks

Modeling infectious disease dynamics in the complex landscape of global health

Modelling to contain pandemics

SARS outbreaks in Ontario, Hong Kong and Singapore: the role of diagnosis and isolation as a control mechanism

Model parameters and outbreak control for SARS

Parameter estimation and uncertainty quantification for an epidemic model

Fitting dynamic models to epidemic outbreaks with quantified uncertainty: A primer for parameter uncertainty, identifiability, and forecasts

Data-based analysis, modelling and forecasting of the COVID-19 outbreak

Parameter estimation and prediction for coronavirus disease outbreak 2019 (COVID-19) in Algeria

Scenario analysis of non-pharmaceutical interventions on global COVID-19 transmissions. Covid Economics: Vetted and Real-Time Papers

Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy

Bayesian calibration of computer models

Efficient calibration for imperfect computer models

Bayesian calibration of inexact computer models

The Design and Analysis of Computer Experiments

A generalized Gaussian process model for computer experiments with binary time series

Generalized Computer Model Calibration for Radiation Transport Simulation

Mathematical Tools for Understanding Infectious Disease Dynamics

Bayesian emulation and calibration of a dynamic epidemic model for A/H1N1 influenza

An epidemiological forecast model and software assessing interventions on COVID-19 epidemic in China

Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: A modelling study

A framework for validation of computer models

Simultaneous determination of tuning and calibration parameters for computer experiments

Adding spatially-correlated errors can mess up the fixed effect you love

The importance of scale for spatial-confounding bias and precision of spatial regression estimators

Calibrating a large computer experiment simulating radiative shock hydrodynamics

Adjustments to Computer Models via Projected Kernel Calibration

R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing

A theoretical framework for calibration in computer models: parametrization, estimation and convergence properties

Empirical Processes in M-estimation

Kernel Poisson regression machine for stochastic claims reserving

Generalized linear models

Semi-parametric generalized linear models

Generalized Additive Models

Soft classification, a.k.a. risk estimation, via penalized log likelihood and smoothing spline analysis of variance

Monte Carlo and quasi-Monte Carlo methods

Surrogates: Gaussian Process Modeling, Design, and Optimization for the Applied Sciences

Efficient and Adaptive Estimation for Semiparametric Models

Introduction to Empirical Processes and Semiparametric Inference

On prediction properties of kriging: Uniform error bounds and robustness

Multiresolution functional ANOVA for large-scale, many-input computer experiments

Practical heteroscedastic Gaussian process modeling for large simulation experiments

MRFA: Fitting and Predicting Large-Scale Nonlinear Regression Problems using Multi-Resolution Functional ANOVA (MRFA) Approach

Heteroskedastic Gaussian Process Modeling and Design under Replication

Comparison of three methods for selecting values of input variables in the analysis of output from a computer code

A simulation of a COVID-19 epidemic based on a deterministic SEIR model

SEIR model for COVID-19 dynamics incorporating the environment and social distancing

SEIR modeling of the COVID-19 and its dynamics

Stability analysis and numerical simulation of SEIR model for pandemic COVID-19 spread in Indonesia

An interactive web-based dashboard to track COVID-19 in real time

Load and Analyze Live Data from the CoViD-19 Pandemic

ODEPACK, a systematized collection of ODE solvers

An introduction to stochastic epidemic models

Stochastic Epidemic Models and Their Statistical Analysis

A primer on stochastic epidemic models: Formulation, numerical simulation, and analysis

SimInf: An R Package for Data-Driven Stochastic Disease Spread Simulations

Exact Stochastic Simulation of Coupled Chemical Reactions

Calibrating functional parameters in the ion channel models of cardiac cells

Bayesian projected calibration of computer models

Computer model calibration with confidence and consistency

The behavior of maximum likelihood estimates under nonstandard conditions

Introduction to Empirical Processes and Semiparametric Inference

Penalized quasi-likelihood estimation in partial linear models

Calibration for computer experiments with binary responses and application to cell adhesion study

Efficient calibration for imperfect computer models

A theoretical framework for calibration in computer models: parametrization, estimation and convergence properties

Empirical Processes in M-estimation

Asymptotic Statistics (Cambridge Series in Statistical and Probabilistic Mathematics)

Weak Convergence and Empirical Processes: With Applications to Statistics

Scattered Data Approximation

Maximum likelihood estimation of misspecified models

This work was supported by NSF DMS 1660477 and NSF HDR TRIPODS award CCF 1934924.

Additional supporting information can be found online, including the mathematical proofs of Theorems 1, 3, 4, 5, and 6, and the R code for reproducing the results in the article.

where h 0 and h 1 are the density functions of y i under the misspecified model and the true model, respectively, and ¾ 0 denotes the expectation with respect to y i ∼ h 0 . Since the true model follows a Poisson distribution, we haveBy the law of large numbers, we haveTherefore,θ MLE, which proves the consistency.Similar to the proof of Theorem 5, the asymptotic normality ofθ LS n can be shown as follows,and