Profile likelihood.

Confidence intervals for point estimates are traditionally calculated from the Fisher information matrix (FIM), which is derived from second-order derivatives of the likelihood at the MLE and thus represents the local curvature of the likelihood [20, 21]. This yields exact results for linear models but does not account for non-zero higher-order derivatives of the likelihood with respect to the parameters. In contrast, the profile likelihood approach represents a generalisation for nonlinear systems. To calculate the profile likelihood (10) for one parameter of interest θ_i that is scanned along its axis around the MLE while re-optimizing all remaining parameters θ_i≠j. Since this approach directly evaluates the likelihood function, it enables the calculation of reliable confidence intervals of the estimated parameters also for nonlinear models [9, 22–24].

In analogy to FIM-based standard errors motivated from the analysis of linear models, the uncertainty of parameter estimates of nonlinear ODE models can be assessed by the profile likelihood-based confidence interval (11) for a confidence level 1 − α. It represents the set of acceptable parameter values, characterized by a profile likelihood value PL(θ_i) below the threshold Δ_α compared to the log-likelihood value at the MLE .

This confidence interval can be interpreted as a continuous likelihood-ratio statistic (12) for which the null hypothesis cannot be rejected if D < Δ_α, following the discussed case with a simplification of nested models with one fixed parameter. Using Wilks’ theorem, the critical value of the test statistic T_α can be analogously defined. In this line, the threshold Δ_α is calculated under the asymptotic assumption N_data → ∞ from the (1 − α)-quantile, i.e. using the -distribution with one degree of freedom.

It should be noted, that in addition to the presented so-called frequentist approaches, also Bayesian methods might be considered as an alternative as they employ a distinct notion of parameter uncertainty. Bayesian methods involve incorporating prior knowledge and updating probabilities based on new data, resulting in posterior probabilities. Indeed, Bayesian methods offer an effective way to handle finite dataset effects, as the posterior does not encounter issues with finite sample sizes and does not require asymptotic assumptions. Such approaches e.g. by utilizing extensive Bayesian Markov chain Monte Carlo (MCMC) sampling have been successfully used for analogous mathematical modelling projects in systems biology and have been compared to the frequentist profile likelihood approach [25]. In general, MCMC sampling results have been shown to reliably propagate the uncertainty of parameter estimates to model predictions, but it has also been discussed that MCMC sampling results might be misleading in particular in the presence of non-identifiability [26]. In addition, defining adequate prior probabilities in practical applications for such approaches can be challenging, even with so-called uninformative priors, as they depend on upper and lower bounds, and differ when applying transformations to the parameters, e.g. logarithmic transformation. In the remainder of this work, we apply the discussed frequentist approach for profile likelihood-based confidence intervals as it was used also in many of the presented benchmark model examples. An additional advantage of the frequentist framework is that does not require the specification of priors since it can be performed by maximizing the likelihood calculated from the experimental data.

In addition to the definition of confidence intervals, the profile likelihood approach also allows for a data-based identifiability analysis of the estimated parameters. Based on the shape of the likelihood profile a classification scheme is available, that again depends on the the statistical threshold in terms of the log-likelihood [9, 27]. A parameter is considered as practically identifiable if its likelihood profile exceeds the statistical threshold given by for values of θ_i smaller and larger than the MLE , i.e. yielding finite profile likelihood-based confidence intervals, as depicted in Fig 1A. If the amount or quality of the available data is not sufficient to indicate a confidence bound in both directions of θ_i, the parameter is termed as practically non-identifiable. In this case, a unique point estimate may be obtained, but the profile does not exceed the threshold Δ_α in at least one direction and thus does not have a confidence interval with finite bounds, c.f. Fig 1B. A constant profile likelihood for the whole parameter domain as shown in Fig 1C indicates a structurally non-identifiable parameter. Such a case arises if the estimated parameters are not unique for the given data set, for example if the model structure contains a redundant parametrization which cannot be resolved by additional data.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Profile likelihood for identifiability analysis.

(A) Illustrative profile likelihood as blue solid line for an identifiable parameter. For analyses of the identifiability status, the intersection with statistical threshold value (red line) at is crucial. Furthermore, profile likelihood-based confidence intervals can be constructed from the projections of the intersections on the parameter axis, defining the lower θ_lb and upper bound θ_ub of the confidence interval, if these intersections exist. The MLE is indicated by the black asterisk. (B) Profile likelihood of a practically non-identifiable parameter that is open to the left resulting in an unbounded (C) Likelihood profile of a structural non-identifiability indicated by a flat line with unbounded confidence interval in both directions.

https://doi.org/10.1371/journal.pcbi.1011417.g001

Bartlett correction.

An improved approximation of such a test statistic by a χ²-distribution has been proposed by Bartlett [28]. It relies on a problem-specific Barlett correction factor C = 1/(1 + df/q), where in the original formulation df is the difference between the dimensions of the parameter spaces of the model under the alternative hypothesis and the null hypothesis and q is a constant which needs to be determined for the individual application. It has been shown that the modified likelihood ratio statistic (13) is closer to the -distribution in applications than the unmodified statistic, i.e. the first three cumulants agree up to an error of [29]. Depending on the actual setting and model formulation, the appropriate Bartlett correction factor is difficult to obtain or may be impossible to compute in realistic applications [18]. The idea of the Bartlett correction was later generalized by Lawley [30] with a generalized expression for the correction factor C as a function of the moments of the log-likelihood derivatives. The analytical derivation of Lawley’s approach, however is cumbersome, as it requires joint cumulants of log-likelihood derivatives up to the fourth order. While there are adaptions for certain model classes that allow a more efficient calculation of the respective Bartlett correction factor [31–34], no methods are available for nonlinear ODE models.

Chebyshev’s and Cantelli’s inequality.

Critical values of the test statistic are essential as statistical thresholds for the discussed likelihood-ratio test related methods. In addition to their derivation based on Wilk’s theorem and calculation from the χ²-distribution, we propose an alternative and more general attempt of defining such thresholds using Chebyshev’s and Cantelli’s inequality.

For a simple illustration, let’s first consider a random variable generated from a normal distribution with mean μ and standard deviation σ. Using the so-called empirical rule, for example for the 3σ neighborhood, the probability P(|X − μ| < 3σ) is known to be 99.7%. A popular implication from this rule is that almost all of the variables of a respective data sample are contained in the 3σ interval, if the numbers are normally distributed [35]. Moreover, it is approximately valid if the underlying distribution is at least bell-shaped with mean μ and standard deviation σ. Assuming that the distribution of a two-tailed test statistic would follow such a (approximative) normal distribution, statistical threshold could be constructed from the empirical rule according to a significance level α = 1 − P. For such a two-tailed test statistic a lower and upper statistical threshold could be used as lower and upper critical values for the two rejection regions at the outer parts of the distribution’s tails.

Under much weaker assumptions on the distribution, Chebyshev’s inequality can be utilized [36]. If the empirical mean μ and standard deviation σ of the sample are finite and known, Chebyshev’s inequality guarantees (14) i.e. in average at least (1 − 1/k²) × 100% of random numbers X drawn from the distribution to be contained in the interval [μ − kσ, μ + kσ]. This is completely independent from the underlying probability distribution and is thus also suitable even in cases when the form of the probability distribution is unknown.

It can be concluded that a realisation of a random number X of any kind of distribution lies beyond k = 4.47 standard deviations of its mean in a maximum of 5% of the cases, i.e. with α ≤ 5%. Again, for the rejection region of an hypothesis test an upper bound on the left tail of the distribution can be defined, as well as a lower bound for the rejection region on the right tail. Chebyshev’s inequality can be also applied to a positive one-tailed distribution, if , which is valid for the -distribution for .

For random variables drawn from single-tailed distributions, a one-sided threshold can be similarly derived from Cantelli’s inequality [37]. For a realisation of random numbers X in a sample with empirical mean μ and standard deviation σ, the probability that a single random number exceeds the sample’s mean by any real number c > 0 is (15) regardless of the underlying probability distribution. Motivated from the conclusion of Wilk’s theorem that the likelihood ratios should be at least asymptotically -distributed, Cantelli’s inequality can be used to derive a value for the statistical threshold in the finite-sample case, i.e. when a certain deviation from the asymptotic distribution is expected. For this, assume an in detail unknown but -like distribution with mean μ = n and standard deviation . It then follows from Eq (15) that a realisation of a random number X from such a distribution exceeds a critical value of T_ξ = μ + c with probability α ≤ 2n/(c² + 2n). The same rationale then applies when assuming an in detail unknown test statistic, which has -like mean μ = n and standard deviation . The critical value threshold T_ξ of this test statistic can then be calculated from Cantelli’s inequality for a significance level of maximal α ≤ 2n/(c² + 2n). In conclusion and for n = 1 degrees of freedom, this yields a critical value T_ξ = 1 + 6.16 = 7.16 for the statistical threshold for a maximal 5% significance level under the assumption of the in detail unknown, but -like test statistic. Of note: the threshold of a -distribution for a significance level α = 0.05 would be .

Parametric bootstrapping.

Bootstrapping methods are a frequently applied to estimate the sampling distribution of an arbitrary statistic [38]. An extension to these resampling methods based on the data sample alone is parametric bootstrapping, which is utilized in the following. Therefore, model trajectories are calculated based on the MLE , which is assumed as underlying truth. The model trajectories are then used for generation of realisations of artificial data samples from the distribution of the measurement errors. For the nonlinear ODE models analyzed here, model trajectories are simulated only from the observables according to the original experimental design using the MLE of the original data set. For each model, N_real = 500 realisations of artificial data sets are generated from simulations and by adding i.i.d. noise from the original error model.

It should be noted, that these simulated data sets follow the original experimental design, but their generation relies entirely on the original data and on the published model structure. The simulated datasets are generated based on the model trajectories and the distribution of measurement errors given by the error model that both were determined by estimating the model parameters according to the original experimental data in the first place. While the simulated data realizations cannot provide new insights into the true biological system, they are useful for investigating uncertainties and assessing reliability in terms of coverage. In principle, also parameter confidence intervals might be calculated from such bootstrap data samples, albeit in applications they are considered as computationally much more demanding than the preferred profile likelihood-based confidence intervals [39, 40].

In this study, artificial data sets produced according to the described parametric bootstrapping procedure serve to generate empirical likelihood ratios and to investigate their distribution. To this end, the model parameters are re-estimated for a bootstrap data realization , yielding the value of the log-likelihood function termed as . In the same way, the log-likelihood value termed as is determined for a nested model with one fixed parameter using the same data realization . This is done by fixing each parameter of interest θ_j to the value obtained from the fit to the original data, while re-fitting the model by re-optimizing all remaining parameters.

Since the parameter of interest is fixed to its true value , i.e. to the value which generates the model trajectories from which the artificial data sets are sampled, the model with L(θ|θ ∈ Θ₀) with nesting is a valid simplification of the full model with likelihood function L(θ|θ ∈ Θ₁). Thus, the likelihood ratio or the difference in terms of log-likelihood according to Eq (8) is constructed from the log-likelihood L of both fits, i.e. the fit of all parameters and the fit with fixed j-th parameter θ_j. This yields the empirical likelihood ratios (16) for each parameter θ_j and data realization from the bootstrapping approach. According to Wilks’ theorem and due to the difference of free model parameters of one (17) should asymptotically hold with df = 1 degrees of freedom of the chi-square distribution.

The described parametric bootstrapping procedure as summarized in pseudocode in S1 Fig is repeated for each of the N_real = 500 data realisations and for each estimated parameter θ_j of the m-th model, yielding the empirical realization of the likelihood ratio λ_ijm for the i-th data realization. Whether the asymptotic assumption from Eq (17) holds in realistic applications, is the central question in the analysis of the empirical likelihood ratios from Eq (16) in the remainder of this study. By calculating the empirical cumulative density function (ECDF) of λ_ijm over all simulated data samples, the empirical distribution can be related to the theoretically expected cumulative density function (CDF), for example to the expected -distribution in order to check the appropriateness of the asymptotic approximation. For this, empirical thresholds of λ for different significance levels α can be compared to the expected quantiles .

Probability-probability plots (pp-plots).

To compare the ECDF of the empirical likelihood ratio samples against the theoretical CDF of the asymptotically expected reference distribution, probability-probability-plots (pp-plots) are utilized [41]. If the sample is distributed according to the theoretically expected distribution, the approximately identical cumulative probabilities for the respective quantiles in both the theoretical CDF and in the ECDF will be observed. For visual inspection, the pp-plot is constructed by plotting the ECDF of the likelihood ratio sample, i.e. the cumulative probability values on the y-axis against the cumulative probability values of the theoretical CDF on the x-axis of the pp-plot, cf. Fig 2A. A data sample distributed accordingly to the theoretical CDF would lie on the diagonal, as shown by the red line in Fig 2A. This outcome is expected for the asymptotic case when the empirical likelihood ratio sample is indeed distributed as .

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Empirical likelihood ratio classification.

A: Probability-probability plot (pp-plot) with four possible graph scenarios. B: -samples show a noticeable deviation from the diagonal in the pp-plot. The pp-plot in panel B shows graphs from 1000 samples each with N_real = 500 random numbers drawn from a -distribution. C: Classifications plot with classification regions and classification scheme for mixed cases: the white graph is classified as conservative, although its upper part lies in the perfect consensus area. Likewise, the light gray graph is identified with the anti-conservative case, while the black graph represents an alternating case which populates all three areas.

https://doi.org/10.1371/journal.pcbi.1011417.g002

If however, the cumulative probability of the ECDF is larger than the expected theoretical CDF, the pp-plot graph lies in the upper left triangle. This case is termed conservative, since the ECDF implies a smaller empirical threshold value for the significance level α, compared to the expected threshold from the theoretical CDF. That is, the same likelihood ratio corresponds to a larger cumulative probability in the ECDF than for the theoretical CDF. In other words, the empirical sample contains more smaller values of the likelihood ratios than expected, i.e. more than the expected 1 − α% of the likelihood ratio values are smaller than the expected of the reference distribution. The consequence of applying such a non-adequate threshold derived from the asymptotic theory to this sample would be for example a too conservative confidence interval. A respective profile likelihood-based confidence interval with a threshold would be too large and the true value of a parameter would lie more often than in 1 − α% of the cases within the confidence interval. The pp-plot graph of such a conservative case is shown in yellow in Fig 2A. It indicates that the ECDF has a larger cumulative probability value compared to the value of the theoretical CDF at the chosen confidence level, cf. the example for a confidence level 1 − α = 0.9 in Fig 2A. For the shown example, more than 90% of the data sample are in fact smaller than the expected .

In contrast, the opposite case of a too small threshold is termed anti-conservative and is identified by a pp-plot graph lying entirely in the lower right triangle, c.f. blue line in Fig 2A. Here, a profile likelihood-based confidence interval with would be too small, since in this example only 71% of the likelihood ratios are smaller than T^theo, in contrast to the expected 90%.

Taking into account the course of the whole pp-plot graph, also a mixture of both cases can occur, as shown by the purple line in Fig 2A. In this alternating case, the graph lies for example in the conservative triangle for small values of the CDF, but crosses the diagonal in the middle of the cumulative probability domain of the CDF and remains below the diagonal for larger p-values in the anti-conservative regime.

Classification of realistic samples.

For the further analysis, a classification scheme for the pp-plots of the empirical likelihood ratios according to the discussed cases is implemented. This is necessary since realistic samples are less smooth than the graphs depicted in Fig 2A and even realisations from a -distributed random variable show a noticeable deviation from the expected straight diagonal in the pp-plot (light green line in Fig 2B). For illustration, 1000 random samples each containing N_real = 500 i.i.d. random variables drawn from the expected -distribution are generated, resulting in pp-plot graphs depicted in shades of gray in Fig 2B. In order to be able to identify empirical likelihood ratios which in fact correspond to the asymptotic theory, but show a similar deviation from the diagonal, a tolerance region around the diagonal is constructed for the classification scheme. This area is constructed from fitting a higher-order polynomial around the 3σ-deviation of the samples around the diagonal This polynomial is bounded at p = 0 and 1 and shows the largest deviation from the baseline at p = 0.5. The resulting tolerance band contains approximately 95% of the -sample graphs. In the context of the bootstrapping procedure for empirical likelihood ratios, a graph which remains entirely within this area of the pp-plot is termed as a perfect consensus with the asymptotic theory and is classified accordingly.

Graphs lying entirely above or below the red perfect consensus-area are identified either with the conservative case (yellow) or anti-conservative case (blue), respectively. Despite these clear cases, also mixed cases needs to be considered. Often, the conservative or anti-conservative graphs do not lie exclusively in their respective areas, but their tails overlap with the perfect consensus-area. Such a pp-plot graph with a deviation from perfect consensus-area in only one direction is assigned to the respective class. For example the white graph in Fig 2C is be assigned to the conservative case, even if parts of it lie in the perfect consensus-area, but not in the anti-conservative area. Likewise, the gray graph is identified as an anti-conservative case.

Contrarily, graphs that have parts in both, the conservative and the anti-conservative part are assigned to the fourth class that is termed alternating graph in the following. This scenario is illustrated by the black graph in Fig 2C. This characterizes a situation where a general statement about the conservativeness of statistical thresholds for likelihood ratios of the respective model parameter is not possible.

In fact, the appropriateness of the asymptotic approximation depends on the specific significance level for which e.g. a likelihood ratio threshold is derived. For this reason, the further analysis of empirical likelihood ratios λ from the parametric bootstrapping procedure covers both aspects: first, a classification of the whole pp-plots graphs in the presented four classes. And in a second analysis, an evaluation of the pp-plot graphs at certain relevant quantiles, i.e. for values of statistical thresholds for specific significance levels α. Such an example for 1 − α = 0.9 is depicted in Fig 2A, using the same classification areas as indicated in Fig 2C for the whole pp-plots graph analysis.

Analysis of whole pp-plots.

In a first step, the whole pp-plots of all parameters are compared to asymptotic theory by assigning their graphs according to the discussed scheme to the four classes. The results of the classification for all parameters and separately for each parameter type are summarized in Fig 4. Fig 5 shows the pp-plots for all 784 parameters, grouped by the models and separated in identifiable parameters in the upper panel in Fig 5A and practically non-identifiable parameters in the lower panel Fig 5B.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Classification results of whole pp-plots graphs compared to the -distribution.

Results are sorted by the percentage of non-problematic, i.e. conservative and perfect consensus cases. The anti-conservative and alternating fraction of cases is indicated by a negative sign. (A) Overall results for all model parameters and (B-E) results separately for each model parameter type.

https://doi.org/10.1371/journal.pcbi.1011417.g004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. All pp-plots graphs of the estimated model parameters, separated by their identifiability status.

(A) identifiable parameters and (B) practically non-identifiable based on their profile likelihoods. All graphs are classified in perfect (red), conservative (yellow), anti-conservative (blue) or alternating (purple) cases.

https://doi.org/10.1371/journal.pcbi.1011417.g005

The analysis reveals a heterogeneous classification outcome throughout the different models, covering models with almost all parameters being in the conservative regime or showing a perfect consensus with theory(8 Hass), but also models with approximately 100% of anti-conservative cases(17 Tönsing-Flu). A surprisingly low overall fraction of only 16.6% of the parameters in the perfect consensus class shows that only few empirical likelihood ratios are indeed -distributed. Only two models (12 Kurzhunov and 18 Tönsing—Zika) have more than 50% of parameters with a perfect overlap with asymptotic theory. The majority of analysed model parameters shows a distribution that is not in accordance with asymptotic theory. However, the non-problematic cases with at least conservative pp-plot graphs represent 51.8% of all analyzed parameters, shown as positive value bars in Fig 4A. Most models have a considerable fraction of anti-conservative or alternating pp-plot graphs. These problematic cases are grouped together with negative percentage in Fig 4. In such cases it can not be guaranteed that the statistical thresholds for the likelihood ratio, e.g. for confidence intervals, are at least conservative. With 6.8% of all parameters, only a small fraction shows an alternating profile, i.e. with parts of the pp-plot graph in the conservative and in the anti-conservative part above or below the diagonal in the pp-plot.

An impression of the effect size of the pp-plot graph deviation from asymptotic theory is given in Fig 6, which shows the average distance between the actual pp-plot graphs and the diagonal. For this, the vertical distance from to the graph to the diagonal is averaged over all 500 empirical likelihood ratios, i.e. data realisations for one parameter and normalized to the maximal possible deviation. The results are ordered with respect to the classification ranking for all parameters from Fig 4A. It shows a similar trend of the mean of the distances per model and the classification ranking, although the variance within the models is rather high. However, most pp-plot graphs show a relatively low distance to the diagonal, as indicated by the gray histogram on the right hand side in Fig 6. If graphs show a large deviation from the diagonal, they tend to be more likely in the upper conservative region than in the lower anti-conservative region. In some cases in the conservative part, however, they reach the maximal value of the distance. Contrarily, in the more problematic case of pp-plot graphs in the anti-conservative region, they stay rather close to the diagonal in the perfect consensus region.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Average distance of pp-plot graphs from the diagonal.

Distances are normalized by the maximal possible distance value, i.e. a graph lying on the upper x-axis at p_emp = 1 for all p_emp in the pp-plot. Models are ordered according to the ranking of the classification results from Fig 4A. The histogram on the right shows the distribution of average distances from the diagonal summarized for all parameters.

https://doi.org/10.1371/journal.pcbi.1011417.g006

When comparing the classification results for all parameters in Fig 4A with the results for only dynamic parameters in Fig 4B, it appears that the partitioning of pp-plot graph classifications changes slightly, but would result in a similar ranking of models based on the fraction of non-problematic cases. Dynamic parameters dominate the overall classification of the models, but their non-problematic fraction is slightly higher. Initial value and observational parameters in Fig 4C and 4D show a similar, yet less pronounced correlation with the overall results from Fig 4A. A comprehensive statement concerning the appropriateness of the asymptotic assumption for these parameter types remains vague, since many models do not have either initial value or observational parameters (cf. Fig 3A).

Error parameters are in most cases classified as anti-conservative, meaning that they fall in the problematic case where a statistical threshold for the likelihood-ratio test from asymptotic theory following Wilks’ theorem is too small compared to the empirical likelihood ratios from the bootstrapping procedure. It follows that, e.g. profile likelihood-based confidence intervals would be too small, regardless of the chosen significance level α. Thus, the uncertainty of the estimated error parameters would be underestimated and consequently the error model predicts the uncertainty to be smaller than it actually is. Interestingly, error parameters share the preferable property of being rather often identifiable which is desirable, but at the same time their pp-plots are mostly anti-conservative, which represents a rather problematic case.

A similar observation can be made concerning the relation between identifiability status and pp-plot classification, when observing at the actual pp-plots in Fig 5 where identifiable and practically non-identifiable parameters are shown separately. Within the identifiable parameters in in the upper panel Fig 5A, all four classes are present, but only few models (2 Bachmann, 9 Isensee, 10 Lucarelli and 16 Swameye) show a considerable amount of more than five conservative cases (yellow lines). Contrarily, in the lower panel Fig 5B showing the non-identifiable parameters almost all pp-plots indicate conservative pp-plot graphs. In contrast, anti-conservative parameters occur only in rare occasions for practically non-identifiable parameters (blue lines in Fig 5B), while almost every model exhibits at least one identifiable, yet anti-conservative parameter in Fig 5A.

Perfect consensus graphs are merely visible in Fig 5B ecause these graphs are concentrated on a rather small region around the diagonal in the pp-plots. Absolute numbers are given for each class on top of the plots. Conservative and anti-conservative graphs are clearly visible, since they are spread over the whole plotting area. More precisely, it can be noted that the blue anti-conservative graphs touch the lower or right hand sided axes of the plots only in few cases. In contrast, the yellow graphs of conservative and non-identifiable parameters in Fig 5B predominantly show a profile which instantaneously increases for p = 0 of the theoretical CDF indicating a rather large number of data realisations with very small empirical likelihood ratios close to zero.

Analysis of relevant quantiles.

In practice, not all quantiles of the test statistic are relevant, since usually study results do not report the whole confidence distribution [61], but rather confidence intervals for a specific significance level α or confidence level 1 − α. Thus in this second analysis, the empirical likelihood ratios are examined only at specific thresholds values of the test statistic for typical assumed confidence levels. In contrast to the analysis of the whole pp-plot graph, here, only the specified confidence levels are employed in the pp-plot graphs, cf. Fig 2A. For this, the three classification regions from Fig 2C are used, but only at the specific significance levels, corresponding to the respective quantiles of the test statistic. Consequently, the class of alternating graphs becomes obsolete, since at a specific level it is either in one of the non-problematic (perfect consensus or conservative) or in the anti-conservative region.

Fig 7 shows the results for the frequently applied threshold of the test statistic for a significance level of α = 0.05. For this, the pp-plots of the empirical likelihood ratios are classified only at the confidence level of 1 − α = 0.95. Specifically, it is assessed whether the calculated statistical threshold of the test statistic’s distribution under the asymptotic assumptions holds by checking if at least 95% of the empirical likelihood ratios from the bootstrapping procedure are smaller than 3.84. This is done by checking, if the pp-plot graph of the empirical CDF at p_theo = 0.95 lies within the perfect consensus region, below or above it. Compared to the results from the whole pp-plot analysis (cf. Fig 4), the analysis at the 0.95-quantile in Fig 7 shows a similar ranking of the models, although considerably more parameters reveal the preferable perfect consensus case. However, the overall amount of anti-conservative cases remains at a similar level and most models again show a heterogeneous pattern, with in the majority of models at least one third of as anti-conservative classified parameters.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Classification results of pp-plots graphs at the 95% confidence level.

Results are sorted by the percentage of non-problematic, i.e. conservative and perfect consensus cases. The model ranking remains similar to the whole pp-plot graph analysis (cf. Fig 4).

https://doi.org/10.1371/journal.pcbi.1011417.g007

The same kind of analysis is repeated for a set of specific confidence levels, in order to check if the appropriateness of the asymptotic assumption improves when considering the whole range of possible confidence levels, including the prominent 68% and 99% levels. Since the ideal situation of mostly perfect consensus cases is not reachable anyhow, the following analyses examine non-problematic cases, i.e. perfect consensus and conservative cases versus the problematic anti-conservative parameters. For this, the conservativeness ratio (18) is defined, which serves as a measure for the appropriateness of the assumed test statistic distribution. In other words, the CR is the fraction of parameters with a pp-plot graph not in the anti-conservative region for a given confidence level. A CR-value close to 1 indicates a good asymptotic approximation of the likelihood ratio statistic at the respective confidence level.

Fig 8A shows the conservativeness ratio CR for a range of confidence levels for all analyzed parameters of the benchmark models, separated in their parameter type, their identifiability status and for the individual models. In the heatmap in Fig 8A it can be observed that the CR indeed varies with the confidence level. As a general trend, it can be concluded that lower confidence levels show higher values of the CR, meaning that they are more often in accordance with asymptotic theory than higher confidence levels.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Appropriateness of the -distribution for empirical likelihood ratios at specific confidence levels.

(A) Heatmap of conservativeness ratio CR, i.e. fraction of non-problematic cases in the respective parameter group or model. CR-values close to 1 indicate a good agreement with asymptotic theory, while lower values occur in situations where the asymptotic approximation might not be appropriate. (B) Parameter groups or models with CR larger than 95% at the confidence level are indicated by black tiles in the upper panel. (C) A less strict criterion checks if the CR value is larger than 1 − α% of the respective confidence level and parameter groups or model.

https://doi.org/10.1371/journal.pcbi.1011417.g008

For a better overview of the results and to assess whether an individual tile of the heatmap in Fig 8A Corresponds to an acceptable CR or if the asymptotic approximation is inappropriate for this combination of confidence level and model or parameter group, two threshold criteria for the CR are constructed.

The first criterion is the stricter one and requires that the asymptotic approximation is valid in at least 95% of the pp-plots at the, regardless of the considered confidence level. The heatmap in Fig 8B indicates that this requirement of CR ≥ 95% is fulfilled for almost no model at any confidence level. Precisely, only for a very small confidence level of 10%, at least 95% of the model parameters show a sufficient accordance with asymptotic theory. Of course, the formulation of such low confidence levels is troublesome, as it corresponds to a not desirable type I error of e.g. 90%. Two exceptions are model ‘6 Bruno’, which does not reach a CR ≥ 0.95 for any of the analyzed confidence levels and model ‘8 Hass’, which exhibits the opposite case of a sufficient CR ≥ 0.95 for all analyzed confidence levels.

The second, yet weaker criterion for a sufficient CR is motivated from having the same type I error at all stages of the analysis, i.e. for the likelihood ratio statistic, the respective confidence intervals, and for the appropriateness of the choice of the -distribution. This translates effectively into a lower threshold for most tiles: a CR ≥ 1 − α would be a sufficient condition for a valid asymptotic assumption per heatmap tile. For example, at the 68%-confidence level, only 68% of the empirical likelihood ratios need to be non-problematic. Fig 8C shows the results for this criterion applied on the heatmap from panel A. It appears that considerable more tiles pass this weaker criterion. However, already for a confidence level of 1 − α ≥ 0.5, many parameter groups or models do not show an acceptable CR and thus, neither for the 95% criterion nor for the weaker criterion can the general validity of the asymptotic approximation from Wilks’ theorem be assumed.

These results confirm what was already concluded from the previous analysis of the pp-plots: The validity of the asymptotic assumption in Wilks’ theorem applied to realistic ODE models cannot be guaranteed, neither for the whole pp-plot, nor for typically utilized quantiles of the distribution, such as e.g. 0.68, 0.9, 0.95 or even 0.99. By comparing the CR results for the error model parameters with the average CR or other parameter types, it appears that they show considerably less agreement with the χ²-distribution. More precisely, error model parameters rather exhibit the anti-conservative case, in which confidence intervals tend to be too small. Furthermore, also the group of identifiable parameters are more often problematic, when compared to the CR of practically non-identifiable parameters. Interestingly, practically non-identifiable parameters rather exhibit the conservative case, i.e. confidence intervals tend to be too large (see also Fig 5). Looking at the individual models, the analysis reveals a large heterogeneity in the CR-values. There are models that show a good agreement of the asymptotic assumption (e.g. model 8 Hass) over all quantiles, or at least for quantiles up to 0.68 or 0.8. But it can be observed that there are also models where the asymptotic assumptions even for a confidence level of 1 − α = 0.5 holds at least conservatively in only 50% of the cases. It should be further noted as conclusion of the analysis at the relevant quantiles, that approximately 50% of the parameters are classified as anti-conservative at the commonly utilized confidence level of 0.95, cf. Fig 7.

Impact of increased data sample size.

After observing the deviation from asymptotic theory in a considerable fraction of cases and having formulated a geometrical explanation of its origin, three approaches handling the finite-sample case are discussed. First of all, the transition from the finite-sample case towards the asymptotic limit is inspected, showing on the one hand that indeed the amount of information contained in the data is the reason for the deviation and on the other hand highlighting that increasing the number of data points does cure the issue for the selected models and experimental designs for the presented models and experimental designs.

To this end, the impact on the pp-plots of an altered experimental design with additional data points in the benchmark models is investigated. Using a rather straightforward approach, the original experimental design is extended by simply simulating additional data data points for time courses of observables which were already contained in the original design. This is done for example by doubling the temporal sampling of all observables. Fig 10 shows the results of two experimental designs with increased temporal sampling for three models and compares the results to the original design.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Impact of experimental designs with increased amounts of data.

Whole pp-plot graph classification results and average normalized distances from the pp-plot graphs to the diagonal for three illustrative benchmark models. The outcome for the originally published design , and two artificial designs and with increased temporal sampling is shown. The respective number of estimated parameters, total number of data points and their ratio is given below.

https://doi.org/10.1371/journal.pcbi.1011417.g010

By shifting the design towards the asymptotic case, i.e. by increasing the number of data points, it can be observed that the classification of the whole pp-plot indeed tends to a higher number of parameters with non-problematic cases for the three models. For the model ‘17 Tönsing-Flu’, the classification even changes from completely anti-conservative to desirable 100% perfect consensus. It can be further observed for the models ‘3 Becker’ and ‘13 Raia’ that the average distances from the diagonal change only slightly, while the classification result indicates a clear improvement of the appropriateness of the asymptotic assumptions. However, a certain amount of conservative and anti-conservative cases persists, even for large amounts of data, compared to the original design. In these cases, the asymptotic regime is not yet reached by this very basic approach of increasing the temporal sampling. Presumably, even very large amounts of additional data points would not be sufficient to reach the asymptotic limit for all parameters, while the incorporation of comparatively few data points from additional observables might carry more relevant information about the systems and would yield a more suitable asymptotic experimental design. Other models did not show such a large effect in the pp-plots by using the increased temporal sampling design. In this line, it should be made clear that in general it is not guaranteed that a simple increase of data points of an experimental design, the asymptotic regime is rapidly reached.

The question how an improved design with additional data and additional observables that drives the setting closer to the asymptotic limit could be constructed, remains open in general, but might be guided by adapting existing optimal experimental design principles cf. [66, 67]. Another practical limitation to the experimental designs arises in typical systems biology related projects, where the increase of the sample size is not a viable option, since data generation in molecular biology is often laborious and costly. Changing the experimental design by including additional observables might also be impossible for example due to missing antibodies for certain compounds. Therefore, we discuss in the following two different approaches to cope with the finite-sample case when analyzing likelihood ratios.

Bartlett correction by extensive search.

In order to correct for the deviation of the likelihood ratio distribution from the asymptotic limit, the Bartlett correction can be used. Although it is not possible to calculate the appropriate correction factor in advance from the model properties in realistic applications of nonlinear ODE models investigated here, an approximate can be derived based on the results from the bootstrapping procedure. For this, the previous analysis of pp-plot classes at the 95%-quantile is repeated but with altered empirical likelihood ratios (21) using a set of 600 Bartlett-like correction factors C_ξ ∈ [0.1, …, 1, …, 600], cf. Eqs (13) and (16). Corresponding pp-plots are calculated, evaluated at the 95%-quantile and assigned to either the perfect consensus, conservative or anti-conservative case, according to the previously discussed scheme. Despite the large number of results from the 600 correction factors for each of the 784 parameters, this enables for a relatively cheap exhaustive search for an optimal Bartlett-like correction factor .

For all except one parameter, at least one individual Bartlett-like correction factor C can be found within the set of C_ξ ∈ [0.1, …, 1, …600], yielding a perfect consensus pp-plot at the 95%-quantile. Only in one special case, a parameter from the model ‘18 Tönsing-Zika’, which implements prior knowledge information, is assigned to the conservative case for all applied corrections factors even in an extended range of . As expected from the previous analysis (cf. Fig 7), the procedure yields approximately 41.2% of the parameters revealing a perfect consensus for uncorrected likelihood ratios with factor C = 1. For the remaining parameters, the individual optimal correction factor is defined as the C_ξ ≠ 1 closest to 1, that yields a perfect consensus pp-plot at 95%.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 11. Bartlett correction results of pp-plot classifications at the 95% confidence level.

Heatmap colors according to pp-plot classification for each Bartlett correction factor as indicated on the y-axis. Black tiles show optimal Bartlett correction factor closest to 1 of the individual parameter in each column. Since columns are ordered by this correction factor, the x-axis corresponds to the approximative conservativeness ratio aCR. White line and black bold numbers indicate the uncorrected outcome for C = 1, whereas green line and numbers show overall optimal Bartlett correction for a CR ≈ 95%. Deviations from 0.95 on x-axis originate from binning issues. (A) Results for all 768 parameters, (B) same results grouped by parameter type and (C) identifiability status. (D) Outcome for individual models reveal that a comprehensive overall optimal correction factor is difficult to determine. (E) Overall optimal Bartlett correction factors C* from exhaustive search with resulting thresholds for the likelihood ratio statistic.

https://doi.org/10.1371/journal.pcbi.1011417.g011

Fig 11A shows the 95%-pp-plot classification result as coloured tile for all correction factors C_ξ as indicated on the y-axis and for all analyzed parameters shown on the x-axis. These parameters are sorted by their individual optimal correction factor as indicated by the black line in the heatmap. Thus, each column shows the effect on the pp-plot classification for the specific model parameter with varying correction factors. It appears that the Bartlett correction affects the classification results in a non-monotonic way: the non-problematic perfect consensus and conservative cases fluctuate for smaller correction factors, whereas there is a sharp separation between the problematic cases and the anti-conservative classification. When sorting the parameters with respect to their optimal correction factors, an approximative conservativeness ratio (aCR) for a specific correction factor appears from the percentage of parameters on the right hand side of the intersection between the specified correction factor and the black line of optimal correction factors. For example, the white horizontal line in Fig 11A indicates the uncorrected scenario with correction factor C = 1. It intersects with the black line of the parameters’ optimal correction factors and a conservativeness ratio CR of 52% can be read off the x-axis, since this fraction of the parameters reveal a perfect consensus or a conservative classification, whereas most of the parameters at the right hand side of the intersection are classified as anti-conservative. Exceptions from the anti-conservative classification on the right hand side occur when the classification is non-monotonic with respect to the correction factor, i.e. when a parameter is classified as conservative for the respective correction factor (C = 1), but the optimal correction factor for a perfect consensus is larger. For this reason, only an approximative conservativeness ratio aCR is available in this visualisation from the projection of the intersection on the x-axis of the heatmap. It slightly differs from the previously calculated conservativeness ratio CR as for example in Fig 7.

The heatmaps in Fig 11 Allow to determine an overall optimal Bartlett-like correction factor C* for all analyzed parameters, but also separate overall correction factors C* of the groups of the specific model parameter types, based on their identifiability status or individually for each benchmark model. This overall optimal correction factor C* is defined such that an appropriate correction yields a conservativeness ratio CR = 0.95 or close to this value. This corresponds to 95% of the parameters in the heatmap classified preferably as perfect consensus case or at least conservative. The overall optimal Bartlett-like correction factor C* can be visually determined in the heatmaps by intersecting the vertical green line at aCR = 0.95 or its closest value with the black curve of individual optimal correction factors. This procedure implies that all parameters on the left-hand side of the intersection are classified as at least conservative for the respective correction factor.

In fact, for all examples in Fig 11, the approximative conservativeness ratio aCR matches the exact conservativeness ratio CR, since the non-monotonic classifications occur only for the upper (right hand sided) 5% of the columns. The green numbers on the y-axis of the heatmaps indicate the C* of the respective parameter subset.

The respective statistical threshold for the likelihood-ratio test can be determined from the Bartlett correction factor when assuming that the corrected likelihood ratio with the respective overall optimal correction factor, is now -distributed, cf. Eq (21). Fig 11E summarizes the overall optimal correction factors for each parameter group and model with the resulting alternative statistical thresholds which could be used in applications instead of the typically assumed T_{α = 0.05} = 3.84 from Wilks’ threorem. The analysis of the individual benchmark models shows a large variability with corrected alternative threshold values ranging from approximately 3.9 to 18.5.

It appears that a generally applicable correction factor C*, i.e. for all parameter types as well as for the group of only dynamic parameters would be C* = 0.52. For the parameter groups of only observational and initial value parameters are closer to the uncorrected scenario with C* = 0.66 and C* = 0.69, respectively. Error parameters, however, require a rather strong correction of the likelihood ratios, indicated by a Bartlett correction factor of C* = 0.21, resulting in a rather high value of the threshold .

The overall analysis analysis for all investigated parameters yields for an alternative Bartlett-corrected statistical threshold. However, the analysis of parameter groups and individual models reveals a large variability of the optimal correction factors, depending on the specific application, i.e. model and experimental design setting. Thus, despite the detailed and extensive search approach it remains difficult to identify a reliable overall Bartlett-like correction factor and threshold that is generally applicable in most situations.

Adapted thresholds using Chebyshev’s and Cantelli’s inequality.

In order to enable a more appropriate description of the likelihood ratio statistic, another approach for an alternative approximation of their unknown distribution is employed that might be used in realistic application scenarios.

While the empirical rule allows for approximative boundaries for e.g. confidence intervals in the case of Gaussian or approximatively bell-shaped distributions, Chebyshev’s and Cantelli’s inequality can be used for any kind of distribution under mild conditions. Here, they are specifically used to draw valid conclusions about the statistical thresholds also in the finite-sample case. For the empirical likelihood ratios, it was shown in the previous analyses that a -distribution in principle is a reasonable choice to describe the likelihood-ratio test statistic. Thus, the idea is to use both inequalities with an in detail unknown, but -like distribution. Only the mean μ = n = 1 and the variance σ² = 2n = 2 of such a -like distribution with n = 1 degrees of freedom need to be plugged into both inequalities for the following conclusions. This reduces the required distributional assumptions drastically without discarding the in principle appropriate -assumption entirely.

It follows from Chebyshev’s inequality in the kσ-neighbourhood for k = 4.47 around the mean that the probability of a random variable drawn from this in detail unknown distribution is P(X ∉ [0, 7.32]) < 0.05. Using k = 4.47 this leads to an alternative threshold of for a 5% significance level or by using k = 3.16 a threshold for a 10% significance level. Likewise, the typically assumed asymptotic threshold for a significance level α = 0.05 following Wilk’s theorem corresponds to an alternative Chebyshev-adapted threshold with significance level α = 0.2477. This corresponds to an actual confidence level of only 75.23% instead of the desired 95%.

In contrast to Chebyshev’s inequality, Cantelli’s inequality is restricted to a positive one-tailed distribution. It yields similar, yet slightly smaller values for the thresholds and . The asymptotic theory threshold value for α = 0.05 yields an alternative Cantelli-adapted threshold for significance level α = 0.1987. Again the resulting confidence level under these assumptions is only 80.13% instead of the expected confidence level of 95% from asymptotic theory.

Using these alternative thresholds, the conservativeness ratio CR is evaluated for the empirical likelihood ratios from the bootstrapping procedure and visualized in Fig 12 for the previously discussed parameter groups and also for the individual models. In terms of the pp-plot it is checked whether the graph at p_theo = 0.95, implying a threshold of T = 3.84 in the theoretical χ²-distribution corresponds to an empirical cumulative probability of at least p_emp > 0.8013 = 1 − α for the Cantelli-adapted threshold . In this scenario, the empirical sample has at least p_emp > 80.13% of its values below the threshold T = 3.84. The pp-plots are analyzed and classified as described earlier, so that a conservativeness ratio CR can be derived for the respective parameters. Fig 12C illustrates the results for the threshold value T = 3.84 with the Cantelli interpretation of the confidence level. Likewise, it is shown how often the non-problematic classifications occur for the empirical likelihood ratios of the respective parameter groups or models for T = 3.84, but with the probability interpretation from asymptotic theory in Fig 12A and and from Chebyshev’s inequality in Fig 12F. For larger values of the alternative thresholds, e.g. in Figs 9E, 9G, 9H and 12D the same pp-plots are evaluated at the corresponding larger p_theo-values of the theoretical CDF of the -distribution and checked against the probabilities p_emp = 0.9 or p_emp = 0.95 of the ECDF of the empirical likelihood ratios.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 12. Conservativeness ratio CR at asymptotic and adapted thresholds.

(A) Histogram of conservativeness ratio CR for the typically utilized 95%-threshold from asymptotic theory for all analyzed parameters, separately for parameter groups and individual models as a reference. (B) Table summarizes the CR for all parameters for alternative confidence levels using asymptotic theory and for adapted thresholds or adapted confidence levels using Chebyshev’s or Cantelli’s inequality. (C, F) All blue histograms including panel (A) share the same threshold value T = 3.84, but differ in the interpretation of the confidence level. (D, G) The green histograms correspond to a confidence level of 90% with slightly different thresholds using Chebyshev’s or Cantelli’s inequality. (E, H) Likewise, purple histograms correspond to a confidence level of 95%. Vertical black solid lines in the histograms indicate the CR > 1 − α and CR > 0.95 limit, which coincides for the asymptotic case.

https://doi.org/10.1371/journal.pcbi.1011417.g012

Overall, these adapted thresholds show a preferable conservativeness ratio CR for most parameter groups and models. In contrast, the asymptotic theory threshold in Fig 12A does by far not reach either the CR > 0.95, nor the CR > 1 − α limit, except for one single model especially error parameters are far off from a desirable CR. The largest confidence level for which at least CR > 1 − α holds on average for all investigated parameters using asymptotic theory is 1 − α = 0.6, cf. Fig 12B. Using the adapted thresholds, the weaker limit of CR > 1 − α is reached for both, the Chebyshev- and the Cantelli-based re-interpretation of the asymptotic threshold value T = 3.84 in all parameter groups and models. Both criteria are met for the overall CR for all parameters and for most models. Error parameters, which were often classified as problematic during this study, now show a conservativeness ratio, which is sufficient for the CR > 1 − α criterion in Fig 12C and 12F. When increasing the confidence level, however, error parameters remain problematic, also when using the Chebyshev- or Cantelli-adaptation.

In general, the adapted threshold for a confidence level of 95% shows only a minor deviation in the CR from the required 95%-limit on average for all parameters, cf. Fig 12B, 12E and 12H. Thus, could serve as an appropriate alternative threshold for likelihood ratios in the finite-sample case for the 95%-confidence level in the future. For existing studies, the typically utilized threshold value of T = 3.84 for an originally assumed confidence level of 1 − α = 95% could be re-interpreted in this line for the finite-sample case using the alternative confidence level 1 − α ≥ 80.13% from the Cantelli-adaptation.