A modified Susceptible-Infected-Recovered model for observed under-reported incidence data

Fitting Susceptible-Infected-Recovered (SIR) models to incidence data is problematic when not all infected individuals are reported. Assuming an underlying SIR model with general but known distribution for the time to recovery, this paper derives the implied differential-integral equations for observed incidence data when a fixed fraction of newly infected individuals are not observed. The parameters of the resulting system of differential equations are identifiable. Using these differential equations, we develop a stochastic model for the conditional distribution of current disease incidence given the entire past history of reported cases. We estimate the model parameters using Bayesian Markov Chain Monte-Carlo sampling of the posterior distribution. We use our model to estimate the transmission rate and fraction of asymptomatic individuals for the current Coronavirus 2019 outbreak in eight American Countries: the United States of America, Brazil, Mexico, Argentina, Chile, Colombia, Peru, and Panama, from January 2020 to May 2021. Our analysis reveals that the fraction of reported cases varies across all countries. For example, the reported incidence fraction for the United States of America varies from 0.3 to 0.6, while for Brazil it varies from 0.2 to 0.4.


Introduction
Susceptible-Infected-Recovered (SIR) models, introduced by Kermack and McKendrick and further developed by Wilson and Worcester [1,2], have been extensively used to describe the temporal dynamics of infectious disease outbreaks [3][4][5]. They have also been widely used to estimate the disease transmission rate by fitting the models to observed incidence data [6][7][8], such as time series of daily or weekly reported number of new cases provided by [9][10][11][12], for example. Implicit in all these model fittings is the assumption that all the infected individuals have been observed. Yet that assumption is problematic when disease incidences are underreported. Under-reporting of incidence is prevalent in health surveillance of emerging diseases [13,14], and also occurs when a disease presents a large fraction of asymptomatic carriers, e.g., Typhoid fever, Hepatitis B, Epstein-Barr virus [15] and Zika [16]. Lack

Model development
Our epidemic model is developed in three steps. First, we extend a generalized SIR model to describe the dynamics of the observed (under-counted) infections. Second, we introduce a local version of that SIR model to describel the evolution of the epidemic in a series of observational time windows given the past time serie of observed incidences. This more flexible model is used to compute the conditional expectation of current observed incidence given the past history. Third, we develop a computationally tractable approximation for the conditional expectation to speed up Monte-Carlo Markov Chain (MCMC) inferences of our model parameters.
Generalized SIR model. Classical mass-action epidemic models, such as the SIR models, are simple yet useful mathematical descriptions of the temporal dynamics of disease outbreak [3][4][5]. These models describe the temporal evolution of the number of susceptible S(t), infected I(t) and recovered R(t) individuals in a population of fixed size N = S(t) + I(t) + R(t). We model their dynamics through the set of integral-differential equations [34]: IðtÞ ¼ RðtÞ ¼ with initial conditions S(0), I(0), R(0). The parameter β measures the transmission rate (also called infection rate [35,36]) and the function F(t) is the cumulative distribution of the time from infection to recovery. When F(t) = 1 − e −γt , the exponential distribution with mean γ −1 , our model reduces to the standard SIR model (see Murray [37] for example). For completeness, the proof of existence and uniqueness of the solution of System (1)-(3) is provided in the appendix. An alternative proof can be found in [34]. The model parameters β and F(t) are epidemiologically relevant and provide insights into the outbreak. For example, the basic reproductive number as defined by Lotka [38,39]: The term S(t)/N in Eq (4) is the fraction of susceptible individuals in the population that can be infected and g À 1 ¼ R 1 0 ð1 À FðtÞÞdt is the average recovery time. R 0 is arguably the most widely used measure of the severity of an outbreak [40,41], at least in the absence of interventions to control it. It measures the expected number of secondary infections attributed to the index case in a naïve population. Other quantities of interest, such as the maximum number of infected individuals and the total number of infections, can be expressed in terms of the reproductive number R 0 , e.g. Weiss [36].
For many diseases, it is reasonable to assume that the disease progression from infection to recovery is known, either because the disease is well characterized, or because the date of onset of symptoms, hospital admissions, and discharge data are available [42]. Thus we will assume throughout this paper that we know the distribution of the recovery period F(t) and we will focus on estimating the transmission rate β.
Modeling the observed disease incidence. LetSðtÞ,ĨðtÞ andRðtÞ denote the observed number of susceptible, infected and recovered individuals as a function of time. We make the following modeling assumptions: (A1) The true underlying dynamics follows the SIR dynamics described by Eqs (1)-(3) with known fixed population size N and time-to-recovery distribution F(t).
(A3) The recovery distribution is the same for observed and unobserved infected individuals.
Under these assumptions, the observed number of infected individuals at time t is and similarly,RðtÞ ¼ pRðtÞ. The number of observed susceptible individuals is Eq (6) follows by solving the differential equationS 0 ðtÞ ¼ pS 0 ðtÞ and using the identitỹ Sð0Þ ¼ Nð1 À pÞ þ pSð0Þ, which results from (A2) and N = S(0) + I(0) + R(0). These equations capture the intuitive idea that under-reported incidence results in a larger number of observed susceptible and fewer infected and recovered individuals through the epidemic evolution. Consider the ratio, which yields from Assumption (A2) and Eqs (5) and (6): For a standard SIR model with p = 1, the ratio v(t) is unity. However, for the observed process, the ratio v(t) starts at one and then monotonically decreases over time. It follows that fitting an SIR model to observed incidence data, neglecting the under-reporting, will produce a nearly unbiased, but possibly noisy estimate for β early in the outbreak when v(t) � 1. As more data becomes available and v(t) decreases, the estimated transmission rate β will under-estimate the true value. As a consequence, one might at later times in an outbreak underestimate the severity of the outbreak and call the epidemic under control prematurely.
The following theorem describes the dynamics of the observed number of susceptible, infected and recovered individuals when only a fraction p of the infected individuals are observed.
Theorem 1 Under assumptions (A1), (A2), and (A3), the process of the observed individuals evolves according to the following set of integral-differential equations: RðtÞ ¼ The conclusion of the theorem follows from algebraic manipulations of Eqs (1) to (6). The addition of the positive term ðð1 À pÞb=pÞĨðtÞ toS 0 ðtÞ implies a slower depletion rate of the observed susceptible population than would be expected under the standard SIR model. Note that this positive term must be small enough such that ÀS 0 ðtÞ � 0, for all t � 0 and all p > 0, condition imposed from Assumption (A2). Assumption (A2) of observing the same fraction p of initial infected and recovered individuals was established only for the technical mathematical proofs of Eq (5) andRðtÞ ¼ pRðtÞ. This mathematical assumption will be relaxed in the following local dynamics definition.
A stochastic model for the observed incidence. Observed incidences of disease are typically reported at regular time intervals. Precisely, let 0 = t 0 < t 1 < t 2 < . . . < t n denote the boundaries of the observation windows. For simplicity, we assume that t k = kΔ, and we denote by Y k the number of new cases of the disease observed in the interval (t k−1 , t k ], k = 1, 2, . . ., n. We also assume that the new cases Y k depend on the actual observed past history of incidences H kÀ 1 ¼ fY 1 ; Y 2 ; . . . ; Y kÀ 1 g. As a result, our model takes into account the impact of fluctuations in the reports. Indeed, imagine that the reported cases Y k are much larger than what is predicted by Model (8)-(10). That excess of cases will alter the observed dynamics of the outbreak, making it progress faster. Similarly, smaller numbers of incidences will slow down the outbreak. The following model takes into account past fluctuations in the incidence to model locally the dynamics of the process at each time interval given the past history.
Definition 1 Let Y 1 , Y 2 , . . ., Y k be the sequence of observed incidences and assume that the cumulative probability distribution F for the time to recovery is continuous. We model the local dynamics of the observed number of susceptibleS k ðtÞ and infectedĨ k ðtÞ individuals at time t in the interval (t k−1 , t k ] through the set of differential-integral equations NpS k ðtÞĨ k ðtÞ þ bð1 À pÞ pĨ k ðtÞ ð11Þ with initial conditionsS k ðt kÀ 1 Þ ¼Sð0Þ À P kÀ 1 j¼1 Y j and with the convention that P 0 j¼1 Y j ¼ 0, where bothS k ðt kÀ 1 Þ > 0 and ÀS 0 k ðtÞ � 0 for all t � 0 and p > 0. For this model, the conditional expectation of incidence given the past history is for all k = 1, 2, . . ., n.

Remark 1
Continuity of the cumulative distribution of the time to recovery F implies that I k ðt k Þ is left continuous. Furthermore, if F has a probability density, thenĨ k ðtÞ admits a righthand derivative at t k−1 .
The local model described in Definition 1 has the same infection dynamics, Eq (11), as the global model. What differs is the evolution of the number of infected individuals, and how it relates to the history of past incidences. The following heuristic serves to motivate Eq (12) in Definition 1. Decompose the integral in Eq (9) for the number of infected individuals into a sum over each observed window (t j−1 , t j ] to writẽ To getĨ k ðtÞ, replaceS 0 ðtÞ by its local instantiationS 0 k ðtÞ on (t k−1 , t k ] andS 0 ðtÞ by Y j /Δ, the empirical rate of new infections, on the interval (t j−1 , t j ]. This is interpreted as assuming that the Y j new infections in the interval (t j−1 , t j ] occur uniformly in that interval. This allows us to take into account the actual number of observed incidence in each time interval instead of using modeled derived quantities, which provides the needs flexibility for our local epidemic model to better track more complex epidemic dynamics than is possible using a global generalized SIR model.
We use the expression for the conditional expectation of incidences in the interval (t k−1 , t k ] given the time series of past observed incidences in Definition 1 to model the conditional dis- A negative binomial counts the number of success in a sequence of identically and independently Bernoulli with probability of success p = μ k /(μ k + r) before r failures (with probability 1p) occur, μ k is the conditional expectation defined in Eq (13). With this parametrization, the conditional expectation and variance are respectively. The shape parameter r controls the amount of over dispersion when compared to a Poisson distribution for which V½Y k jH kÀ 1 � ¼ E½Y k jH kÀ 1 �. In particular, as the shape parameter r grows to infinity, the negative binomial model converges to a Poisson distribution with rate μ k . Thus, the negative binomial distribution allows us to account for the extra-Poisson variability that arises in the data. Other distributions are possible, such as beta negative binomial distribution [43] or the Conway-Maxwell-Poisson distribution [44].
With repeated application of the chain rule, we combine the set of conditional distributions for Y k |Y 1 , Y 2 , . . ., Y k−1 into a joint likelihood for the model parameters Gðy k þ rÞ where Γ denotes the gamma function and μ k depends only on β and p. Since, in the model formulation, the distribution of Y 1 does not contain any information about the transmission rate and the fraction of observed cases, the term P½Y 1 � is dropped from the likelihood. Approximation of the conditional expectation. To reduce the computational burden required to numerically solve the set of differential-integral Eqs (11) and (12), and the ensuing integration in Eq (13) to evaluate the conditional expectation, we propose to approximate the conditional expectation μ k by linearizing bothS k ðuÞ andĨ k ðuÞ around t k−1 in Eq (13), and integrate the result explicitly. The following lemma encapsulates the resulting approximation.
Lemma 1 Assume that the cumulative probability distribution F for the time to recovery has a probability density f. The conditional expectation μ k can be approximated by whenĨ kÀ 1 6 ¼ 0 andS kÀ 1 =N � ð1 À pÞ, and μ k = 0 otherwise. Here, for all k = 1, 2, . . ., n. Proof of Lemma 1. Eqs (19)- (22) follow directly from the definition ofSðt kÀ 1 Þ and the evaluation of Eqs (11) and (12) at t k−1 . To prove Eq (22), we first take the derivative ofĨ k ðtÞ, Eq (12), with respect to t and simplify it as follows: Then, we evaluate at t k−1 and simplify the resulting equation, using the definition of each From the definition, bothS kÀ 1 andĨ kÀ 1 are non-negative quantities, and the hypothesis S kÀ 1 =N � ð1 À pÞ implies that ÀS 0 kÀ 1 � 0, for all p > 0. Therefore, all these equations are well defined. In the proof of Eq (18), the linear approximation of bothS k ðuÞ andĨ k ðuÞ around t k−1 are:S k ðuÞ �S k ðt kÀ 1 Þ þ ðu À t kÀ 1 ÞS 0 k ðt kÀ 1 Þ ð23Þ Substituting these equations in the integrand of Eq (13) and solving it yields: (20) and (22) kÀ 1 in the previous equation and simplifying it yields: where the conclusion of Eq (18) follows. Remark 2 Better approximations for μ k can be obtained using higher order Taylor expansions forS k ðuÞ andĨ k ðuÞ. This requires the distribution F of time to recovery to have higher order derivatives.
Identifiability. It is known that the measured growth rates in early SIR outbreaks are insensitive to under-reporting. Indeed, in early outbreaks, S(t) � N and hence I 0 (t) � (β − γ)I (t). Under Assumption (A2), we have thatĨðtÞ ¼ pIðtÞ andĨ 0 ðtÞ ¼ pI 0 ðtÞ, which imply that It follows that the disease incidence grows exponentially with rate β − γ, irrespective on the fraction p of observed incidence. Hence the transmission rate β can be estimated if the recovery rate γ is known, but the fraction p cannot be estimated at that early stage of the outbreak. As the outbreak matures and moves away from its early exponential growth phase, it becomes possible to estimate both the transmission rate β and the fraction p of observed cases. The following theorem provides verifiable conditions for both these parameters to be identifiable.
The proof of Theorem 2 is found in the appendix. Remark 3 As we note earlier, β and and p are not identifiable in the early stages of an outbreak. This is also evident in Theorem 2: In the early stages, we have thatS k ðuÞ � N, so that the vectors (U 1 , . . ., U k ) and (V 1 , . . ., V k ) are essentially co-linear. Later in the outbreak, as S k (u) is no longer close to N, both parameters become identifiable.

Bayesian parameter estimation
We use the Metropolis-Hastings algorithm to draw Monte-Carlo Markov chain (MCMC) [45] samples from the posterior distribution of the model parameters given the epidemic outbreak data. Our implementation transforms the original parameters Θ = (β, p, r) intõ Y ¼ ðx; Z; rÞ 2 R 3 , where ξ = log(β), η = log(p/(1−p)) and ρ = log(r), and selects proposals from a multivariate Gaussian distribution N ð�jỸ m Þ ¼ NðỸ m ; SÞ; with meanỸ m and diagonal covariance matrix S with entries 0.001, 0.01, and 0.01. The results presented in the next section are from 40,000 MCMC samples gathered after 40,000 burn-in iterations when starting from Θ 0 = (0.5, 0.5, 25). Our implementation used the approximation for μ k presented in Lemma 1.
Following [46,47], we model the distribution of time to recovery from COVID-19 as the convolution of a lognormal distribution (with mean = 5.2 and sdlog = 0.662) with a Weibull distribution (with mean = 5 and sd = 1.9). The mean and standard error of the resulting recovery time distribution are 10.27 and 4.32, respectively. We refer the interested reader to [30] for a detailed description of additional disease progression parameters of SARS-CoV-2 infection.
Separate chains were run for the time series of incidence data from each country, using all the data from the date of the first confirmed COVID-19 cases to May 18th, 2021 (see Table 1).
The assumption of a constant transmission rate does not hold, as each country implemented various mitigation and control strategies, from national lockdown orders to closing of public meeting places (see Table 1 which shows the date on first implementation of mitigation as reported in [48]). To avoid having to model the change in the transmission rate resulting from the implementation of mitigations, our parameter estimation starts on the first day of intervention as reported in Table 1. We still use the whole time series from the time of first confirmed incidence to estimate the number of infected individuals as defined by Eq (12).
To reduce the impact of weekly reporting patterns (e.g. fewer cases are reported over the weekend) we apply a moving average of seven days to the raw incidence counts before executing the MCMC algorithm. Finally, the initial conditionsS 0 ¼ N À pIð0Þ À pRð0Þ, R(0) = 0, I 0 ¼ pIð0Þ are set using the reported national population counts and number of initial cases as reported in Table 1.

Analysis of COVID-19 incidence data
We performed separate Bayesian inferences for eight American Countries: the United States of America (USA), Brazil, Mexico, Argentina, Chile, Colombia, Peru, and Panama. Fig 1 shows histograms of the marginal posterior distribution of the transmission rate β after the start of mitigation, the fraction observed p, and the negative binomial shape parameter r for each country. The median and 95% credible intervals of these posterior distributions are presented in Table 2. Even though each country used different mitigation strategies, with various level of enforcement, the credible intervals for the transmission parameter of each of the eight countries overlap, with the exception of Peru. There are several hypothesis for why this may be the case: the effectiveness of the various mitigation strategies is compromised by having a small fraction of non-compliant individuals, or most of the benefits of the mitigation strategy are achieved by wearing face masks and moderate social distancing. A third hypothesis is that the estimated transmission rate in our model is a time average of the instantaneous transmission rates, and that averaging lessens the differences in transmission rates.
Similarly, the posterior distributions for the fraction of observed incidence are similar across most of the analyzed countries. The two exceptions are Peru and Mexico, with the under-reporting in Mexico being particularly acute. This is consistent with the observation that Mexico has one of lowest numbers of tests performed per reported case [49]. While an under-reporting factor of about 15 is very large, we believe this effect is real because of how well the model fits the data (see the appendix) and narrowness of the posterior distribution.
Related analyses of COVID-19 data in Mexico have used values for the fraction of reported cases of p = 0.2 or p = 0.4 to analyze and forecast the evolution of the COVID-19 pandemic and hospital demands [50,51]. These values are closer to the values that we found for the other Latin American countries. However, these values were not derived from the data. It would be interesting to use our model to investigate the under-reporting in Mexico at a county level to see how the results would differ from local to national levels.
Excess deaths [52] provide an alternative measure of the true impact of COVID-19. Using that measure, [53] reports that COVID-19 deaths in Mexico are under-reported by a factor of 3, whereas we show a factor of 15 for under-reported incidence. This difference may be due differential testing rates of deceased and infected individuals which may arise from the standard of care of severely ill patients admitted to intensive care units that requires COVID-19 testing [50,54].
Our analysis flags Peru as being different from the other countries both in term of having a higher transmission rate, and a lower reported fraction. Our analysis does not reveal why this is the case, and further analysis incorporating country level explanatory variables to predict transmission rates and under-reporting is needed to uncover the reasons why Peru is different from the other countries in America we studied.
Finally, the estimate of the shape parameter r of the negative binomial distribution shows that the relative inflation of the Poisson variance ranges from 2%-5%. That effect is statistically significant. Again, the distributions across the eight countries are commensurate, with the United States and Peru exhibiting more extra Poisson variability than the other countries.

Under-estimation of the transmission rate
In light of Eq (7), we suggested in the introduction that failing to account for under-reporting leads to underestimating the transmission rate β. Here, we numerically demonstrate this effect by fitting an SIR-type model directly to raw incidence data, deliberately neglecting to model under-reporting. The median and 95% credible intervals of the posterior distribution for the transmission rate when modeling under-reporting, and when not are displayed in Table 3.
The parameter β p coincides with the values from Table 2, while β 1 refers to estimates when the fraction observed is p = 1. The posterior distribution for shape parameter r were similar for p unknown and p fixed.
Observe that in all cases, the 95% credible intervals for the transmission rate do not overlap. This shows that knowledge of the fraction of reported incidence is statistically important.

Variation on the fraction of reported cases
In this section, we consider modeling and estimating a time dependent fraction p(t) of reported incidence, which can arise from uneven availability of COVID-19 tests [31,50,55]. To this end, we model the reported fraction p(t) with a piece-wise constant function: for all 0 < t � t n � ξ M , where I ½x kÀ 1 ;x k Þ denotes the indicator function for each interval [ξ k−1 , ξ k ). We regularize the sequence of reported fractions p 1 , p 2 , . . ., p M by adding the penalty to the loglikelihood. We assume that the variation between reported fraction, p k − p k−1 , are identically and independently normally distributed with mean zero and variance 1/λ. Similarly as in the previous section, we performed separate Bayesian inferences to estimate the posterior distributions of β, p 1 , p 2 , . . ., p M , r and λ for each analyzed country: the United States of America, Brazil, Mexico, Argentina, Chile, Colombia, Peru, and Panama. We transform p k and λ into η k = log(p k /(1 − p k )) and l = log(λ) and use uniform improper priors on the transformed parameters. We defined p(t) with constant pieces of length modulo 90 days and we use the equations from Lemma 2 to compute the expected incidence μ k . These equations generalize the equations from Lemma 1 when p = p k for all k = 1, 2, . . ., M, see the appendix for further details. The median and 95% credible intervals of the posterior distributions of β, r, and λ are presented in Table 4. For clarity in the presented results for λ values, we decided to round them to the nearest integer values. The analogous results for the posterior distributions for each reported fraction, p 1 , p 2 , p 3 , p 4  Additionally, in the first panel of Figs 2-4 and S1-S5 Figs we show the credibility of the model estimates for the daily COVID-19 incidence for each country, where the expected median of reported cases, μ k , are plotted in red lines, the upper and lower credible intervals are plotted in blue lines, while the expected incidences lie in the blue-shadow area with probability of 95%. The negative binomial distribution function, Eq (14), was used to build the credible intervals. To estimate the expected cases, the parameter values for β and r were set equal to the values provided in Table 4 and the p k values were set to the estimated median values of p(t) as shown in the second panel of each figure and for each country, respectively.
From Table 4, the marginal posterior distributions for the parameter λ overlap for all analyzed country. All these marginal posterior distributions skewed to the right with large values. For most countries, the credible intervals for p(t) include a constant function. That is, statistically, we do not have enough evidence to reject the hypothesis that the reported fraction p(t) for each country is not a constant function during the entire analyzed data set. And for countries that have a small variance λ −1 for the increment p k − p k−1 , we have further evidence that p (t) is nearly constant. The one country for which a constant p(t) is not retained is Mexico (see S1 Fig).
The second panels of Figs 2-4 shows that there are some variations across all p k credible intervals for the United States of America, Brazil, and Peru. Interestingly, the credible intervals of p k for each country are all contained in a wider credible interval than those obtained when assuming a constant fraction p of observed cases as reported in Table 2. Comparing Tables 2  and 4, we see that there are not significant changes on the posterior distributions for β and r when we assume p constant and p variable. In general, we observe more variation for the observed proportion p across the countries than within a country. The latter result is not surprising, as countries implemented different testing policies which may affect the way the incidence data were reported [49,50,55].

Strengths and weaknesses of the proposed local SIR model
Our model locally exploits the SIR dynamics, using past observations to set the initial conditions. This results in a flexible model that can fit complex patterns, such as multiple waves that typically require a time varying transmission rate, with a single parameter. This flexibility comes at a cost: our single estimated transmission rate is a time average of the true time varying one. And while we show that our model empirically fits the data well within the credible intervals, we over-estimate the expected incidence in the valleys and under-estimate near the peaks. It follows that the derived estimates for the reproductive number near a local bottom of an outbreak will have a positive bias, leading to a more conservative view of the effect of mitigation.
Our formulation can be generalized to build epidemic models having non-parametric transmission rates. Such models will alleviate the weakness discussed above, and can be used to identify model-based uncertainties in models. These extensions will be presented in a forthcoming paper. We are also planning to extend the model by incorporating the exposed class, which will provide a more realistic model to study COVID-19 pandemic. As COVID-19 disease progression depends on both the length of time an individual remains in the exposed and infectious classes [30]. This model extension would help us to analyze the effect of different

PLOS ONE
infectious period distributions that could change at the early outbreak due to interventions such as testing, isolation or contact tracing. Finally, our model has a limited ability to estimate time-varying under-count fractions. Numerical experiments have shown that adding more flexibility to how the latter varies over time degrades our ability to estimate the transmission rate.

Conclusion
We present a new extension of the standard SIR epidemiological models to study the underreported incidence of infectious diseases. The new model reveals that fitting a SIR model type directly to raw incidence data will under-estimate the true infectious rate when neglecting under-reported cases. Using the epidemic model we also present a Bayesian methodology to estimate the transmission rate and fraction of under-reported incidence with credible intervals that result directly from incidence data. We also argue that our statistical model can properly track and estimate complex incidence reports, where the resulted estimates update as more data are incorporated.
Using our methodology on the COVID-19 example, we found that the credible intervals for the transmission rates overlap across the eighth analyzed American countries: the United States of America, Brazil, Argentina, Chile, Colombia, Peru, and Panama. In all the cases, the median transmission rates are above 0.105 and below 0.122 (see Tables 2-4). And, for most countries, the credible intervals for the time dependent fraction of reported cases p(t) include a constant function, and they also provide a range values for the fraction of reported cases per each country. In average, from January 03, 2020 to May 18, 2021: the reported incidence fraction for the United States of America and Panama varies from 0.3 to 0.6; the reported  which is obtained by taking the derivative with respect to t of Eqs (2) and (3) and using Property 1. Therefore, it is enough to prove existence and uniqueness of solutions of System (29)- (31). It follows that the function G: U ! R 3 defined by GðS; I; R; tÞ ¼ ðS 0 ðtÞ; I 0 ðtÞ; R 0 ðtÞÞ ð32Þ is continuously differentiable in U, see for example [56, pp. 32]. Since @G @S , @G @I , @G @R , and @G @t exist and are continuous in U, then G is continuously differentiable in U. Therefore, the solution of System (29)-(31) exists for the initial condition S(0), I(0), R(0) and is unique in K.
Proof of Theorem 2 Set α 1 = β/p and α 2 = β(1 − p)/p. Since identifiability of α 1 and α 2 implies identifiability of β and p. We can estimate α 1 and α 2 by minimizing the sum of squares The two parameters are identifiable if and only if the vectors (U 1 , . . ., U m ) and (V 1 , . . ., V m ) are not co-linear.

Modeling the time dependence fraction of reported incidence
The following definition describes the dynamics of the observed susceptible and infected individuals when constant fractions p k of infected individuals are observed at each interval (t k−1 , t k ], i.e.,S 0 k ðtÞ ¼ p k S 0 ðtÞ for all t in that interval of time. This hypothesis allows us to study the case when the parameter p is a piece-wise time dependent function, as it is defined in Eq (27).
Definition 2 Let Y 1 , Y 2 , . . ., Y k be the sequence of observed incidences and assume that the cumulative probability distribution F for the time to recovery is continuous. We model the local dynamics of the observed number of susceptibleS k ðtÞ and infectedĨ k ðtÞ individuals at time t in the interval (t k−1 , t k ] through the set of differential-integral equations: Np kS k ðtÞĨ k ðtÞ þ bĨ k ðtÞ