Comparative analysis of practical identifiability methods for an SEIR model

1.   Introduction

Compartment models have been used extensively to study infectious diseases. Among them, susceptible-exposed-infectious-recovered (SEIR) models have been extensively used to study disease dynamics impacted by vertical transmission ^[1,2], vaccination strategies ^[3], delayed infectiousness ^[4], multistage infectiousness ^[5], and treatment strategies ^[6], and most recently, were applied to COVID-19 ^[7,8,9]. These models often include parameters for which numerical values are unknown a priori and cannot be directly measured. Since many parameters in a given model are not directly measurable, researchers often obtain parameters from outbreak data.

Parameter estimation relies on comparing empirical observations of the modeled system with the corresponding model output. Many computational techniques can be employed for parameter estimation. Most of these techniques rely on minimizing the difference between model output and observed data. However, before numerically estimating parameter values, it is important to address whether the parameters of the model are identifiable. Identifiability analysis determines to what extent and with what level of certainty the parameters of a model can be recovered from the available empirical data ^[10]. These relate to two general types of identifiability analysis: structural and practical identifiability. Structural identifiability is the theoretical possibility of determining the true values of parameters of a model from observations of its outputs and knowledge of its dynamic equations ^[11,12,13]. On the other hand, practical identifiability provides information on the accuracy with which parameters can be estimated from the available discrete and noise-corrupted measurements ^[14]. So, for an infinite amount of noise-free data, structural identifiability implies (the maximum) practical identifiability ^[15].

Three common methods for evaluating the practical identifiability of models are Monte Carlo simulations ^[16,17], correlation matrices ^[18,19,20], and the profile likelihood ^{[21,22,23,24,25]}. The Monte Carlo approach can be applied to evaluate disparate model structures each with different observation schema by implementing a random sampling algorithm ^[16]. While Monte Carlo simulations offer conceptual and algorithmic simplicity, their computational cost can be very high, as many samples are required to obtain a good approximation. The Correlation Matrix approach assesses the correlation between parameter estimates ^[18,19,20]. It is much less computationally intensive than the Monte Carlo method but only provides a pairwise analysis of model parameters. Likelihood profiling is a common way to perform a practical identifiability analysis when using a likelihood-based estimation procedure (e.g., approximate Bayesian computation or maximum likelihood estimation) because the profiles are used for quantifying uncertainty in parameter estimates (i.e., to approximate confidence intervals) ^[23].

Epidemiological data plays a critical role in infectious disease surveillance, as it describes how infectious diseases are distributed within populations and what factors contribute to transmission ^[26]. Several mathematical models have been established during the last decades to forecast disease progression using epidemiological data ^{[27,28,29,30,31,32]}. There are many technical challenges in implementing a standardized epidemiological data framework. The way data is reported and collected, and the type and frequency of data available through data-sharing institutions (public health departments, ministries of health, data collection, or aggregation services) differ significantly. For example, different data types (including prevalence, incidence, cumulative incidence, and mortality) are reported at differing time intervals (such as daily, weekly, monthly, or even more infrequently). As a result, understanding how these factors affect parameter identifiability and estimation is critical.

In this study, we systematically examined the influence of different data types (prevalence, incidence, and cumulative incidence) and sampling frequencies (daily, weekly, and monthly) on the practical identifiability of SEIR model parameters. We explored four scenarios with different peak infection times and data collection windows to assess the impact of different data collection strategies on parameter identifiability. Discrepancies between Monte Carlo and correlation matrix results highlight the importance of considering multiple criteria for identifiability. Overall, our results show that incidence data yielded the most identifiable parameters, prevalence data exhibited intermediate identifiability, and cumulative incidence data resulted in the least identifiable parameters. Varying data collection frequencies affected parameter identifiability, with longer time series and higher sampling rates generally improving identifiability.

In Section 2, we present the SEIR model and the structural identifiability results. Sections 3.1 and 3.2 outline the use of MC and CM methods to assess parameter identifiability, respectively. In Section 4, we present the results of our analysis, with Section 4.1 detailing the outcomes of the MC simulations, Section 4.2 illustrating the CM results, and Section 4.3 comparing the results from both MC and CM methods. In Section 5, we discuss our finding, their implications, and limitations. Finally, in Section 6, we provide a concise summary of our findings.

2.   The SEIR model

The overarching goal of our study is to examine how different data types and collection frequencies influence parameter identifiability under structural and practical identifiability paradigms. To this end, we chose a simple SEIR model framework containing a small number of parameters as the foundation for our study. Here, the SEIR model is a closed population with no vital dynamics (i.e., no births or deaths):

where $N(t) = S(t)+E(t)+I(t)+R(t)$ is the total population, the parameter $\beta$ is the transmission rate from susceptible to infected, $\gamma$ is the rate of transition from exposed to infectious, and $\alpha$ is the recovery rate from infected to recovered (Figure 1). The state variable $I(t)$ represents the prevalence of the infectious disease at time $t$. The auxiliary variable $C(t)$ tracks the cumulative number of infectious individuals at time $t$ since the start of the outbreak. $C(t)$ is not a state of the system but rather a class to track cumulative incidence, from which incidence can be derived. Thus, the instantaneous number of new infections (i.e., incidence) at time $t$ is $\dot{C}(t)$. In practice, the incidence of new cases is measured over a time interval $[t_{i-1}, t_i]$, and therefore can be computed using the difference of the cumulative incidence at these time points: $C(t_i)-C(t_{i-1})$. Initial conditions for the $S, E, I,$ and $R$ states will be noted by $S(0), E(0), I(0),$ and $R(0)$, respectively.

The SEIR model Eq (2.1) is well studied and often presented as an entry-level educational tool ^[33]. As individuals progress through the compartments, eventually the number of susceptible individuals may be insufficient to sustain transmission. This leads to the disease-free equilibrium, which is globally asymptotically stable when $\beta N /\alpha < 1$ and unstable when $\beta N /\alpha > 1$. This cutoff is the basic reproduction number $\mathcal{R}_0$, which is commonly used to characterize the infectiousness of a disease that is freshly introduced into a completely susceptible population. Importantly, this reproduction number $\mathcal{R}_0$, along with other epidemiological metrics, are calculated from model parameters. If a model is non-identifiable, then researchers risk reaching incorrect conclusions about disease persistence that follow from incorrect parameter values.

2.1.   Structural identifiability

Structural identifiability addresses the question of whether the parameters of a model can be uniquely determined from "perfect" input-output data, assuming the underlying model correctly represents the input-output data. Here, perfect data means the data is continuous in time and without noise. Structural identifiability serves as a prerequisite for practical identifiability analyses as it sets the foundation on whether the model parameters can be distinctly ascertained. Without first establishing local or global identifiability, the results attained in practical identifiability analysis can be challenging to interpret due to uncertainty in the numerical optimization. Consider a general compartment model:

where the $n$ state variables at time $t$ are represented by the state vector $\mathbf{x}(t)\in\mathbb{R}^n$, parameters are denoted by the vector $\mathbf{p}\in\mathbb{R}^k$, the initial condition vector is denoted as $\mathbf{x}_0\in\mathbb{R}^n$, and the $m$ observation variables are represented by $\mathbf{y}(t)\in\mathbb{R}^m$. The system Eq (2.2) is globally structurally identifiable for the parameter vector $\mathbf{p}_1$ if for every parameter vector $\mathbf{p}_2,$

We say that system Eq (2.2) is locally structurally identifiable for the parameter vector $\mathbf{p}_1$ if, for every $\mathbf{p}_2$,

Otherwise, system Eq (2.2) is unidentifiable. An equivalent characterization of local identifiability is if there are at least two or more finite parameter values of $\mathbf{p}$ that yield the same observations $\mathbf{y}$. Here, the state variables denoted by $\mathbf{x}$ are specifically $[x_1, x_2, x_3, x_4, x_5] = [S, E, I, R, C].$

There are several well-established methods for determining the structural identifiability of disease models such as the direct test ^[34,35] and the similarity transformation method ^[36,37,38], which can be applied to autonomous (no external input) systems. Alternatively, the Taylor series approach ^[39], differential algebra ^[40,41], the generating series approach ^[42], implicit functions approach ^[43], or differential geometry ^[44] are applicable to external input. Various software tools have been developed to perform structural identifiability analysis, such as the observability test ^[45], differential algebra for identifiability of systems (DAISY, ^[40]), exact arithmetic rank (EAR, ^[46]), COMBOS ^[47], Data2Dynamics ^[48], STRIKE-GOLDD ^[44], GenSSI2 ^[49], structural identifiability analyser (SIAN, ^[50]), and StructualIdentifiablity.jl^[51]. These toolboxes provide a broad range of features and implement methods based on differential algebra, semi-numerical differential algebra, generating series, local algebraic observability, and identifiability tables. However, some of these toolboxes have limited performance or can be time-consuming to use ^[52].

We establish the structural identifiability of the prevalence and cumulative incidence data types using the differential algebra approach with the aid of the software package $\texttt{StructuralIdentifiability.jl.}$ The details are given in the following propositions.

Proposition 1. The system Eq (2.1) is locally identifiable to prevalence data and globally identifiable if initial conditions are known.

Proof. Given the observation is prevalence data $I(t)$, that is, $y(t) = I(t)$, the input-output equation derived by $\texttt{StructuralIdentifiability.jl}$ is:

In this equation, $I'$ represents the derivative of the state variable $I(t)$. Using the definition of being globally identifiable, we set $y(t, \mathbf{p}_1) = y(t, \mathbf{p}_2)$, where $\mathbf{p}_1 = [\alpha_1, \beta_1, \gamma_1],$ and $\mathbf{p}_2 = [\alpha_2, \beta_2, \gamma_2]$ which yields the following equations:

Solving the system of equations above using the Groebner basis approach in terms of $\mathbf{p}_1$ in $\texttt{Python}$, we obtain two sets of solutions:

Hence, the system Eq (2.1) is locally identifable given prevalence data. Furthermore, the function $\texttt{find ioequations}$ within $\texttt{StructuralIdentifiability.jl}$ outputs all identifiable equations including the initial conditions which are listed below:

It is clear to see that, with the initial condition $S(0), E(0), \text{and } I(0)$ known, the identifiable equation $S(t)\gamma$ combined with $\alpha\gamma$ implies that both $\alpha$ and $\gamma$ become globally identifiable. □

Proposition 2. The system Eq (2.1) is locally identifiable to cumulative data and globally identifiable if initial conditions are known.

Proof. Given the observation being the prevalence data $C(t)$, that is, $y(t) = C(t)$, we can obtain the input-output equation using $\texttt{StructuralIdentifiability.jl}$. Following a similar procedure as in Proposition 2.1, we obtain the following set of equations:

We solve this system using the Groebner basis approach in terms of $\mathbf{p}_1$ in $\texttt{Python}$ and obtain two sets of solutions:

Thus, the system Eq (2.1) is locally identifiable given cumulative data.

Furthermore, to investigate the impact of initial conditions, we can find all identifiable equations including the initial conditions:

It is clear to see that, with the initial condition $S(0), E(0), \text{and } I(0)$ known, the identifiable equation $I(t)\beta, \beta\gamma$ combined with $\alpha+\gamma$ implies all three parameters become globally identifiable. □

As far as we know, all of these packages can only handle observations when explicitly expressed as functions of the state variables at any time. Consequently, we are unable to examine the structural identifiability of the SEIR model for incidence defined as a difference in time, so its analysis will be limited to practical identifiability. Table 1 summarizes the outcomes of these assessments. In cases with specified initial conditions, all parameters of the SEIR model exhibit structural identifiability concerning both prevalence and cumulative incidence outputs. However, when initial conditions are unspecified, only $\beta$ is structurally identifiable with respect to prevalence, while the remaining parameters are only locally identifiable. Cumulative incidence shows all parameters are only locally identifiable.

3.   Methods for practical identifiability

Compared to structural identifiability analysis, practical identifiability analysis accounts for the sampling rates and noisiness of experimental data ^[53]. For this reason, practical identifiability analyses typically involve fitting models to epidemic data. In general, a model is considered practically identifiable if a unique parameter vector for a given model can be consistently obtained. Researchers may use real or simulated epidemic data to assess practical identifiability ^{[17,54,55,56]}, and there are different approaches and criteria for assessing identifiability. In this study, we apply two approaches to analyze the practical identifiability of a SEIR model with hypothetical "true" parameters. We first use a Monte Carlo approach to consider the practical identifiability of our system of ODEs before conducting a similar analysis using a correlation matrix approach.

In general, practical identifiability can be defined for the dynamical model in Eq (2.2). In this study, we consider $\mathbf{x} = [S, E, I, R, C]$ and $f(\mathbf{x}(t), \mathbf{p})$ the right-hand side of Eq (2.1). We assume that $\mathbf{x}_0$ is known but that all parameters $\mathbf{p} = [\beta, \gamma, \alpha]$ must be recovered from data. For practical identifiability, observations $\mathbf{y}(t)$ at a finite set of time points $t_i$ are modeled by:

with random measurement noise $\epsilon(t_i)$ drawn from a normal distribution with zero mean and standard deviation $\sigma$, i.e., $\epsilon \sim \mathcal{N}(0, \sigma)$. The set of time points $t_i$ are sampled at different frequencies and time spans depending on the scenario. In this study, we denote $\mathbf{y}$ as prevalence, incidence, or cumulative incidence data and define $g\big(\mathbf{x}(t), \mathbf{p}\big)$ to select the appropriate state: $g(\mathbf{x}(t_i), \mathbf{p}) = x_3(t_i)$ for prevalence, $g\big(\mathbf{x}(t_i), \mathbf{p} \big) = x_5(t_i)$ for cumulative incidence, and $g\big(\mathbf{x}(t_i), \mathbf{p} \big) = x_5(t_i)-x_5(t_{i-1})$ for incidence. A model is practically identifiable if there is a unique parameter vector $\mathbf{p} = \hat{\mathbf{p}}$ that minimizes the difference between the model output $g$ and the data $\mathbf{y}$.

Data collection plays a key role when estimating parameters for an epidemic model. For example, case data may be collected on a daily, weekly, or monthly basis. To determine how various time series data impact the practical identifiability of the model parameters, we explored four scenarios with different peaks of the epidemic curve and lengths of data collection. See Supplemental Figure 3 for plots of each scenario.

Scenario 1: The peak for the epidemic curve occurred at day 109 with the time span of 365 days, where $\beta = 0.0001$, $\gamma = 0.2$, and $\alpha = 0.03$.

Scenario 2: The peak for the epidemic curve occurred at day 25 with the time span of 50 days, where $\beta = 0.001$, $\gamma = 0.2$, and $\alpha = 0.03$.

Scenario 3: The peak for the epidemic curve occurred at day 109 with the time span of 100 days.

Scenario 4: The peak for the epidemic curve occurred at day 25 with the time span of 20 days.

Note that the time spans for Scenarios 1 and 2 were chosen to include roughly twice the time it takes for incidence (new infections) to drop below one. We selected this cutoff as a balance between capturing the entire prevalence curve and avoiding many observations near the steady state for incidence and cumulative incidence. Scenarios 3 and 4 simulate the effect of only having access to data before the peak of the epidemic.

3.1.   Monte Carlo simulations

The history of Monte Carlo (MC) simulations can be traced back to the work of ^[57]. The MC method, as a sampling technique using random numbers and probability distributions, can be used to determine the practical identifiability of a model. This approach is due to its versatility and straightforward implementation. We perform MC simulations by generating $M =$10,000 synthetic data sets using the true parameter vector $\hat{\mathbf{p}}$ and adding noise to the data in increasing amounts. MC simulations are outlined in the following steps:

(1) Solve the SEIR ODEs (Eq 2.1) numerically with the true parameter vector $\hat{\mathbf{p}}$ to obtain the output vector $\mathbf{g}\big(\mathbf{x}(t), \hat{\mathbf{p}}\big)$ at the discrete data time points $\{t_i\}_{i = 1}^n$.

(2) Generate $M =$10,000 data sets with a given measurement error. We assume the error follows a normal distribution with mean $0$ and variance $\sigma^2(t)$; that is, the data are described by

where $\epsilon \sim \mathcal{N}(0, \sigma)$ at the discrete data time points $\{t_i\}_{i = 1}^n$ for all $j = 1, 2, ..., M$ data sets.

(3) Estimate the parameter vector $\mathbf{p}_j$ by fitting the dynamical model to each of the $M$ simulated data sets. This is achieved by minimizing the difference between model output and the data generated for the specific scenario:

This optimization problem is solved in MATLAB R2021a using the built-in function $\texttt{fminsearchbnd}$, which is part of the Optimization toolbox. Since $\texttt{fminsearchbnd}$ is a local solver, the optimized minimum value can be influenced by the starting point. To avoid issues related to the starting value, we use the true parameter values as the initial parameter starting point provided to $\texttt{fminsearchbnd}$. Furthermore, we used the optimization function $\texttt{fmincon}$ to test the consistency of the results in the methodology. Both functions produced the same qualitative results with similar AREs (the detailed ARE results are provided in the Supplemental document). For the remainder of the manuscript, the results shown refer to AREs from $\texttt{fminsearchbnd}$.

(4) Calculate the average relative estimation error (ARE) for each parameter in the set $\mathbf{p}$ following ^[16]:

where $p^{(k)}$ is the $k$-th parameter in the set $\mathbf{p}$, $\hat{p}^{(k)}$ is the $k$-th parameter in the true parameter vector $\hat{\mathbf{p}}$, and $p_j^{(k)}$ is the $k$-th element of $\mathbf{p}_j$.

(5) Repeat steps 1 through 4, increasing the level of noise ($\sigma = 0, 1, 5, 10, 20, 30\%)$.

Following the convention in ^[17], we assume that a given parameter is practically identifiable if the ARE of the parameter is less than or equal to the measurement error, $\sigma$.

3.2.   Correlation matrix method

Although the MC simulation approach is easy to understand and simple to implement, the associated computational cost is high due to a large number of repetitions and the associated optimization costs. One alternative is to utilize the sensitivity matrix of the model to compute the coefficient matrix for parameters of interest ^[20,58,59]. This requires much less computation and is relatively simple if measurement errors follow an identical and independent distribution. The correlation matrix (CM) method assesses the correlation between estimated parameters using a matrix of output sensitivities to model parameters. If estimated parameters are highly correlated, then they are considered practically unidentifiable.

CM assessments of identifiability are local to a particular parameter. The assessment also depends on the type of model observation, frequency of data collection, and assumed distribution of measurement noise. However, CM does not incorporate realizations of noisy observations, as in the MC approach. Instead, these considerations are incorporated in the following steps:

(1) Solve the SEIR model sensitivities numerically, for model observations $\mathbf{g}(\mathbf{x}(t), \mathbf{p})$ at the discrete collection times $\{t_i\}_{i = 1}^n$ with respect to each parameter in $\mathbf{p}$.

(a) For prevalence and cumulative incidence, this is obtained from the SEIR model sensitivity equations and the extracted solutions for $x_3$ ($I$) and $x_5$ ($C$):

(b) For incidence, this is obtained from numerically integrating the quantity

where $\dfrac{\partial I}{\partial \mathbf{p}}$ and $\dfrac{\partial S}{\partial \mathbf{p}}$ come from solutions to the sensitivity equations.

(2) Construct a sensitivity matrix where the $i, j^{\text{th}}$ component corresponds to the sensitivity of the model output $\mathbf{g}(\mathbf{x}(t), \mathbf{p})$ at time $t_i$ to the $j^{\text{th}}$ parameter in $\mathbf{p}$. This is denoted by

(3) Compute the inverse of the weighted Fisher information matrix, $IM = (F^TWF)^{-1}$, where $W$ is the diagonal weighted matrix for the least squares error. For the assumed distribution of error in Eq (3.1), the weights are $g(\mathbf{x}(t_i), \mathbf{p})$.

(4) Compute the correlation coefficients $\chi_{ij} = IM_{ij}/\sqrt{IM_{ii}IM_{jj}}$ for all $i, j = 1, 2, 3$ and $i\neq j$, corresponding to correlations between the three parameters in $\mathbf{p}$.

(5) If the correlation coefficients are below 0.9 for all parameter pairs, then the parameter is practically identifiable for the assumed model observations.

The choice of 0.9 for the identifiability threshold in step five is chosen to reflect that a value close to unity represents a non-identifiable parameter pair (that is, these parameters are perfectly correlated) ^[17]. However, one of the challenges of the CM approach is that the determination of identifiability is qualitative in nature, and to make it quantitative requires the user to choose a value for this threshold. We conduct the CM approach for all data types described at the start of Section 3, to assess the practical identifiability of the model for the true parameters. We also repeat this process for estimated parameters ${\mathbf{p}}_j$ obtained in Step 3 of the MC approach (Section 3.1). This allows us to assess the perceived identifiability of parameter estimates obtained from noisy data.

4.   Results

In this section, we will examine the output of the MC and CM methods described in Section 3 when applied to the SEIR model described in Section 2. For all of the estimated parameters, we chose the search bounds in the $\texttt{fminsearchbnd}$ algorithm to be $[0, 1]$. Although, in practice, we may not know the true bounds of the parameters, this interval encompasses biologically feasible values with respect to their epidemiological definitions. For this model, the upper search bound corresponds to a highest possible transmission rate of 1000 new infections per infected individual each day and the shortest possible exposure and recovery times of 1 day each. We are interested in how different factors impact practical identifiability results and whether these results are aligned for both methods (MC and CM). We do not expect to have a one-to-one relationship in each scenario; however, it is beneficial to examine when and why we obtain matching outcomes when they occur in the simulations.

4.1.   Monte Carlo (MC) results

Using the MC algorithm in Section 3.1, we tested the practical identifiability of the SEIR model parameters with respect to different data types, sampling frequencies, and time periods of available data. A model parameter is said to be practically identifiable if the ARE values given by Eq (3.2) are less than or equal to the noise percentage $\sigma_0$ for the given noise level. For example, if the ARE values for all levels of noise for $\beta$ are less than the corresponding noise level, we say $\beta$ is practically identifiable. If the ARE for $\beta$ is more than its corresponding noise level for any noise level, we would say this parameter is not practically identifiable. If all parameters are practically identifiable in a given scenario, then the model is practically identifiable. This is consistent with the definition of practical identifiability used in ^[17]. The identifiability results from the MC approach are summarized in Table 2, and the details of the ARE values can be found in the Supplemental Tables 5–16.

Before examining the ARE values, we generate violin plots for each scenario to observe the distribution of parameter values as noise is introduced in the MC procedure. Although we are not using the violin plots as a condition for practical identifiability, it does provide insight into how the values for each parameter change as we incorporate noise. Figure 2 represents the violin plots for prevalence, incidence, and cumulative incidence in Scenario 1 with respect to $\beta$. Predictably, as more noise is added to the data set, we observe a larger dispersion of parameter values for $\beta$, especially for 20% and 30% noise levels. The remaining plots for the other scenarios, data type, and parameters can be found in the Supplemental Figures 4, 5, 6, and 11.

As summarized in Table 2, a higher data density and longer temporal data availability tend to result in greater parameter identifiability overall. Moreover, the parameter $\gamma$ had the fewest identifiable individual scenarios, with the highest number of identifiable scenarios being $\beta$ and $\alpha$ depending on the disease data metric used. More scenarios of $\alpha$ were identifiable when using prevalence, while incidence resulted in more scenarios of $\beta$ being identifiable. Additionally, the incidence scenarios had the most cases where parameters were identifiable, while only $\beta$ is identifiable in some cases with cumulative incidence data. In particular, one special case where shorter temporal data yields better identifiability results occurred for parameter $\beta$ with cumulative incidence data. This may be attributed to the performance of the optimization algorithm. It is also worth noting that comparing the weekly (monthly) results from Scenarios 1 and 2 with the daily (weekly) results from Scenarios 3 and 4, one could conclude that knowing more about the whole epidemic curve is more important for obtaining an identifiable model than obtaining dense data for a shorter period of time. A detailed discussion and reasoning about these findings are presented in Section 5.

4.2.   Correlation matrix (CM) results

We next test the practical identifiability of the same scenarios as in the MC results, using the CM approach. Since CM reports identifiability of parameter pairs, instead of individual parameters, we holistically define a problem to be identifiable if all three parameter pairs have a correlation below 0.9. If any pair is more correlated, then the remaining estimates may be affected and result in an unidentifiable problem. In Supplemental Tables 29–31, we report these correlation results and, in Table 3, we summarize the resulting assessments of practical identifiability. We find that the most identifiable data type is incidence data, which leads to identifiable problems for any sampling rate in Scenario 1 and daily sampling rates in Scenarios 3 and 4. Estimating from prevalence data is an identifiable problem for daily and weekly sampling rates in Scenario 3, and problems using cumulative incidence are never identifiable.

In most cases, parameter correlation values do not significantly change with sampling frequency. Only six correlation values increase by more than 0.1 between daily and monthly data, and none of these changes cross the 0.9 threshold for identifiability. Four correlation values increase past the 0.9 threshold between daily and monthly data, changing the assessment of identifiability for two data types (monthly prevalence and incidence data in Scenario 3). We therefore find two cases (out of nine) where increased sampling frequency improves the CM assessment of identifiability. However, this is a limited number of occurrences and it does not suggest an overall relationship between sampling frequency and identifiability under the CM criteria. In all other cases, the CM results are not sensitive to the frequency of data collection. Contrary to expectations, four correlation values decrease between daily and monthly data, although the change is less than 0.1 and does not cross the 0.9 threshold for identifiability.

4.3.   Comparing MC and CM

In Supplemental Tables 32–34, we assess identifiability using the CM criteria for all estimates obtained in the MC process, under $0\%$ and $30\%$ noise levels, and then report the percentage of MC estimates that would be considered identifiable. At a $0\%$ noise level, the CM results for the MC estimates match the results for the true parameters, since the estimates are very close to the true parameters. However, at a $30\%$ noise level, the MC estimates may be far from the true parameter. For example, using daily prevalence data, the true parameter vector for Scenario 3 is identifiable by the CM approach. However, only $18\%$ of the MC parameter estimates at $30\%$ noise are identifiable by the CM criteria. Although $\beta$ and $\alpha$ are identifiable by the MC criteria and remain close to their true values, there is a higher error in the estimate for $\gamma$. Because the CM results are affected by these incorrect $\gamma$ values, the entire problem is considered unidentifiable. In other cases, incorrect MC estimates may be considered identifiable by the CM criteria. For example, using daily or weekly cumulative incidence data, the true parameter vector for Scenario 1 is not identifiable by the CM criteria. Under the MC criteria, $\beta$ is identifiable using daily data, but is very close to the cutoff, and all other parameters are not identifiable. Using weekly data, none of the parameters are identifiable. However, $66\%$ of the MC parameter estimates using daily data and $34\%$ of the MC estimates using weekly data are identifiable by the CM criteria.

We next consider parameter correlations in four different cases, which correspond to all possible assessments of identifiability compared across the MC and CM criteria. The cases include daily prevalence data from Scenario 1, weekly prevalence data from Scenario 3, daily incidence data from Scenario 1, and daily cumulative incidence data from Scenario 1. In Supplemental Figure 13, we plot the normalized error of MC estimates using data with $30\%$ noise, for estimated pairs of $\beta$:$\gamma$, $\beta$:$\alpha$, and $\alpha$:$\gamma$. We fit a straight line through the parameter error and report the slope, as a numerical assessment of correlation between MC parameter estimates. We consider two parameters correlated if the lines are sufficiently different from a strictly horizontal or vertical fit, which we define to be when the magnitude of the slope is between 0.5 and 2. This threshold is selected for symmetry around slopes with magnitude 1 and bounded away from horizontal and vertical lines, which would indicate no correlation. We select the specific bounds such that the two cases identifiable under the MC criteria also have uncorrelated parameter estimates. We compare these values to the CM results and present the correlations and identifiability in Table 4. We find that the sign of the correlation matches across the MC estimates and CM calculations in all but one case ($\alpha$:$\gamma$ for daily cumulative incidence data from Scenario 1). However, there is not a consistent match in the correlation magnitudes across MC estimates and CM calculations. The overall classification of "correlated" matches for nine (out of 12) pairs, but the thresholds for both criteria are somewhat arbitrary. For example, in Supplemental Tables 35 and 36, we assume wider and narrower ranges for MC slopes to be classified as "correlated", where the overall classification of "correlated" only matches in eight or seven cases.

5.   Discussion

In this study, we systematically investigated the impact of different data types (prevalence, incidence, and cumulative incidence) and sampling frequencies (daily, weekly, and monthly) on the practical identifiability of the SEIR model parameters, using two identifiability methods: the Monte Carlo (MC) method and the correlation matrix (CM) method. We found that incidence data, sampled at a higher frequency, resulted in the greatest degree of identifiability. While more data (as obtained at a higher sampling frequency) should intuitively lead to greater identifiability (although some exceptions are discussed in more detail below), it is less clear why incidence data should lead to greater identifiability compared to cumulative data. As such, we examine more closely the structure of these data types, and how they relate to model parameters. The disease metrics incidence and cumulative incidence are related to the flow into the exposed compartment, $E$, which depends on the transmission rate parameter $\beta$. On the other hand, prevalence is the result of integrating over the flows into and out of the infectious class, $I$, which are directly impacted by disease-progression rate $\gamma$ and recovery rate $\alpha$, and only indirectly by $\beta$. Because cumulative incidence "smooths out" the information encoded in incidence, it is possible that different combinations of disease progression and recovery rates can yield very similar cumulative incidence curves. For example, if a decrease in the incubation period causes an increased rate of flow into the infected class, this could be balanced out by an increase in the recovery rate out of the infected class and thus lead to similar values of the cumulative incidence at the current time. This hypothesis is further supported by the transmission rate being the only identifiable parameter for cumulative incidence.

A more significant contributor, however, to the unidentifiability of cumulative incidence data is how cumulative incidence is generated. We argue that the method for adding measurement error in the case of cumulative data is patently unrealistic. The error model presented in our work (motivated by other studies on identifiability in outbreak models) assumes that the mean of the measurement error is proportional to the solution curve $C(t)$, resulting in very large measurement errors late in the epidemic. A more realistic way to generate the bootstrapped cumulative data would be to add measurement error to the incidence data and accumulate the incidence to obtain the corresponding cumulative incidence data, which results in much smaller measurement errors later in the epidemic compared with the first approach described here–a result of incidence tapering off after the peak of the epidemic. The challenge with the second approach is that the accumulated incidence data produces measurement errors that are not independent across time.

There are identifiable parameters in each of the various levels of data collection frequency, but the same cannot be said for each scenario. That is, given the low enough noise in our data, we are able to identify at least one parameter for each frequency of data collection, but we find that supplying data for only a portion of an epidemic is a limiting factor for practical identifiability. This is evident in Scenarios 3 and 4, where truncated time series data is used in the estimation process, resulting in few, if any, practically identifiable parameters. Since these scenarios mimic what would occur during an outbreak, care should be taken when estimating parameters before the epidemic has reached its peak, as there is likely to be a great deal of uncertainty in the estimates, and therefore, in the trajectories that follow. Projections of hypothetical control problems should take this uncertainty into consideration.

There are more identifiable problems by the CM criteria than problems where all parameters are identifiable by the MC criteria. However, CM does not allow for partially identifiable problems, and there are more cases where the MC approach classifies some parameters as independently identifiable than fully identifiable problems under the CM approach. The MC and CM assessments of "completely identifiable" problems seldom match, except for Scenario 1 daily incidence data. The assessments do, however, agree on the "completely unidentifiable" problems for Scenarios 1 and 2 using weekly or monthly cumulative incidence data, Scenario 3 monthly data of all types, and Scenario 4 daily prevalence data. For both assessments, we find that problems using incidence data are most often identifiable and problems using cumulative incidence data are least often identifiable. The assessments do not appear to follow similar patterns across sampling rates or cutoff dates. The MC approach more often follows the pattern that increased sampling rates (daily versus monthly) or longer time series (Scenario 3 versus Scenario 1) lead to more identifiable problems, while CM is less sensitive to sampling rates and is inconsistent in its response to longer time series.

Overall, we find that more data does correspond to more identifiable problems for the MC approach or does not affect identifiable problems for the CM approach. However, contrary to expectations, we find some cases where recovering parameters using shorter time series (Scenarios 3 or 4) are identifiable problems, where problems using longer time series (Scenarios 1 or 2) were not identifiable. Under the MC criteria, when using weekly cumulative incidence data, $\beta$ is identifiable for Scenario 3 but not Scenario 1. This also occurs under the CM criteria using daily incidence data, where Scenario 4 is an identifiable problem but Scenario 2 is not. These results suggest that long tails without change in the time series may reduce identifiability. For the MC approach, the long tails may present an opportunity for the minimization to "fit to noise". This is supported by the difference across error levels, where the estimates have more error at lower noise levels for Scenario 3 (Supplemental Table 15) but as noise increases, the estimates for Scenario 1 have more error (Supplemental Table 13). Intuitively, this also makes sense for the CM approach, which relies on the model's sensitivity matrix; once the model solution reaches zero or levels off, there is less sensitivity to the model parameters. This explanation does not hold for prevalence data, which does not level off in our scenarios, but we find that under the CM criteria, when using daily or weekly prevalence data, Scenario 3 is an identifiable problem but Scenario 1 is not. We do not have an explanation for this outcome, except that Scenario 3 is borderline unidentifiable (0.89 correlation, in Supplemental Table 29) and that the 0.9 cutoff for identifiability may not be appropriate.

In Section 4.3, we further examine how MC and CM assessments of identifiability are related. In general, we found little overlap between the two, and by applying the CM criteria to MC estimates, we found that the CM criteria can lead to misleading results about parameter estimates. Most importantly, we found that parameter estimates that did not meet MC criteria for identifiability were sometimes identifiable by the CM criteria. This indicates that in practice, it is possible to estimate incorrect values for unidentifiable parameters but conclude by the CM criteria that the recovered parameters were identifiable. In other cases, the true parameter vector was identifiable by the CM criteria and the MC remained close to the true parameters, but the MC estimates were unidentifiable by the CM criteria. This discrepancy would also carry over to sensitivity-based confidence intervals, making it possible to have low confidence in a well-estimated parameter due to relatively small changes in the parameter space. However, it is important to note that our MC results represent a best case, in which numerical minimization starts from the true parameter estimate; in practice, the high sensitivities that affect the CM approach may affect estimates that start from other parts of the parameter space. In assessing the correlations of MC estimates, we found that the sign of correlation between parameters matched the results from the CM approach, but there was not a clear relationship between the magnitude of correlation for both approaches.

For the MC method, we tested two optimizers, $\texttt{fmincon}$ and $\texttt{fminsearchbnd}$, from MATLAB for numerically finding the optimal parameters. $\texttt{fmincon}$ is a function typically used for constrained optimization where the parameter space is limited through equality and inequality constraints and the objective function is at least second-order differentiable. It uses an interior point algorithm to find the optimal solution by default unless otherwise specified. On the other hand, $\texttt{fminsearchbnd}$ requires upper and lower bounds for the parameter spaces, but cannot handle any other constraints on the parameters. It uses the Nelder-Mead simplex direct search algorithm, which does not require differentiability of the objective function. We implemented these two functions not to directly compare the numerical results, but to determine if they provided similar identifiability outputs.

We computed the AREs six times with different optimization algorithms, number of iterations, and data generation processes for each of the 27 cases (9 cases for each of the three data types as shown in Table 2) and presented the results for the run with the most iterations using $\texttt{fminsearchbnd}$ and $\texttt{fmincon}$ in the Supplemental Material. Notably, these computations exhibit both qualitative and quantitative variations. The first two times, we generate 1000 simulated data for each case, employing the $\texttt{fmincon}$ and $\texttt{fminsearchbnd}$ optimization methods separately. These two optimizers give different identifiability for 5 cases. To mitigate the stochastic variability arising from data generation, we computed the AREs two additional times using the same 1000 simulated data for each case across both optimizers. Then, the number of different cases is only one, though a different case than the previous 5. Furthermore, to reduce inter-case stochastic variability resulting from the utilization of partial data from Scenarios 1 and 2 in Scenarios 3 and 4, we conduct the final two rounds of computations. In these rounds, we initially generated 10,000 simulated data sets for Scenarios 1 and 2 and subsequently extracted the required data sets for all other cases. We observe that both optimizers yielded consistent qualitative outcomes across all cases. Based on this observation, we recommend increasing the number of simulated data sets whenever feasible, in order to enhance the reliability of the results.

Quantitatively, there are cases where AREs fluctuate around the $30\%$ threshold over these six times. For example, the ARE corresponded to $30\%$ noise data for $\beta$ with prevalence as observation ranged from $28.3\%$ to $30.9\%$. It would be difficult to assert practical identifiability in such cases. There are arguments that practical identifiability can be assessed as long as these thresholds are relatively on the same level. However, we found there is one case where the AREs ranged from $16\%$ to $63\%$, while in another case, the AREs only vary from $16.1\%$ to $17.1\%$. This poses a difficulty in the certainty of the results from the MC method. Our intensive computations suggest that the typically used 1000 simulations for the MC method may be insufficient in practice.

In addition, two distinct computational behaviors emerge when computing AREs for the same data: 1) In some cases, when the parameters are not identifiable, $\texttt{fmincon}$ produces notably larger ARE than $\texttt{fminsearchbnd; }$ 2) There are occurrences where the ARE computed from $\texttt{fminsearchbnd}$ is nearly zero while the results from $\texttt{fmincon}$ distinctly deviate from zero. The first phenomenon is likely due to the utilization of an unbounded search region by the $\texttt{fmincon}$ optimizer, which arises from the absence of constraints on the parameter space. In contrast, when employing the $\texttt{fminsearchbnd}$ optimizer, we impose a range of $[0, 1]$ for all parameters, thereby obtaining less AREs by confining the search within a bounded parameter space. The second phenomenon could be attributed to the search algorithm employed by these two optimizers. Our investigation indicates that, in instances where a non-zero ARE is obtained for zero noise data, the $\texttt{fmincon}$ optimizer generates optimal parameter values that deviate from the accurate ones, even when the initial parameter values are set to the accurate values. Notably, the resulting parameters yield a larger error compared to the accurate values. This suggests that the algorithm encounters random deviations from the initial search point and is unable to converge back to the accurate parameters, potentially due to the intricate landscape of the error function.

6.   Conclusions

In this study, we found that, for a range of parameter values, sampling rates, and criteria for practical identifiability, parameters estimated from incidence data are most often identifiable and parameters estimated from cumulative incidence data are least often identifiable. In practice, estimates obtained from cumulative incidence data should be treated with caution. However, it is important to note that identifiability was sensitive to underlying parameters across different scenarios. Additionally, we found that assessments of identifiability seldom agreed across the MC and CM approaches, with each method relying on cutoffs for identifiability that may be arbitrary. A further complication is that parameter estimation and identifiability results were sensitive to the choice of minimization algorithm, even after implementing common safeguards through numerical tolerances. Taken together, sensitivities to parameter values, the definition of identifiability, and numerical algorithms demonstrate that a single assessment of identifiability may be misleading or incomplete. We recommend a range of tests for identifiability to ensure confidence in the results. Finally, we note that our results are obtained for a relatively simple epidemic model under a best-case scenario, in which the model perfectly matches the processes generating the data, and were generated with knowledge of the underlying parameter vector, allowing for numerical searches to begin from the desired solutions. These conditions are impossible to meet in practice, and we still found fewer identifiable problems than expected. Parameter estimation for more complex models must be treated with significant caution and may require more rigorous testing of identifiability than is currently practiced. Although challenging, work in this area would greatly benefit from concrete guidelines for assessing practical identifiability from an array of methods.

Author contributions

Omar Saucedo, Amanda Laubmeier, Tingting Tang, Benjamin Levy, Lale Asik, Tim Pollington, Olivia Prosper Feldman: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Validation, Visualization, Roles/writing-original draft, Writing-review & editing. All authors have read and approved the final version of the manuscript for publication.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This paper is based on work supported by the American Mathematical Society Mathematics Research Communities through National Science Foundation grants 1641020 and 1916439. OPF was supported in part by the National Science Foundation grant DMS-2045843.

Conflict of interest

The authors declare that they have no conflicts of interest.

Supplementary materials

Figures 7 through 10 plot normalized objective function value (fval) and normalized ARE for a given iteration versus other metrics of instance. Normalized ARE for a given scenario was obtained by dividing the given ARE value by the noise level for the scenario. Thus normalized ARE values below 1 are considered identifiable while values greater than 1 are unidentifiable. The normalized fval was calculated as

With this in mind, positive values indicate the estimated parameter values fit the data better than the true parameter values while negative values indicate the true parameters provided a lower objective function value.

Figure 7 shows ARE values plotted on a log scale against our normalized objective function value (fval). The data is further stratified by parameter type and noise level. The normalized objective function values tend to be above 0 indicating that the estimated parameters fit the noisy data better than the true parameter values. Scenarios with the lowest ARE and fval values with the smallest variation tend to correspond to scenarios with a larger number of data points. This indicates that the optimization procedure struggles to identify the true parameters when noise is added to a small number of data points. Increasing ARE values tend to be associated with scenarios that have fewer data points and fvals that take on a larger range of values.

We notice similar dynamics in Figure 8, where the number of data points on a log scale is plotted against the normalized fval for different noise levels and data types. The normalized fval tending to be above 0 indicates that the estimated parameter values provide a better fit relative to true parameter values. As the number of data points decreases, the fval shows a wider range, again demonstrating how adding noise to a small number of data points complicates the parameter estimation process.

Figures 9 and 10 plot the number of data points in a given scenario on the log scale versus ARE as broken down by parameter and noise level. ARE values below 1 are considered identifiable while values greater than 1 are unidentifiable. Beta clearly produced the lowest ARE values across all noise levels and data types. We see that as long as there are enough data points, beta is the most likely to be identifiable.

The following tables contain the correlation coefficients for the different scenarios and data types, rounded to two decimal places. If the correlation coefficient is above 0.9 for any pair of parameters, the parameter vector is "not identifiable" for the scenario and data type.

The following tables report the percent of identifiable parameters from MCMC iterates. In three scenarios, we omit parameter estimates, which result in a noninvertible Fischer information matrix. This happens for monthly Prevalence data under Scenario 3 (984 parameters omitted), weekly Cumulative Incidence data under Scenario 2 (4 parameters omitted), and monthly Cumulative Incidence data under Scenario 2 (641 parameters omitted). The corresponding results are indicated with an asterisk in the tables below.