BACKGROUND
a key tool in systems biology is mathematical modeling
<cit> .
modeling allows for a more complete analysis of the true relationship between experimental data and possible mechanistic explanations, compared to what is feasible using only classical biochemical reasoning.
nevertheless, because the data are limited and the systems are highly complex, and because many of the model parameters have to be estimated and cannot be uniquely determined, drawing mechanistic conclusions from modeling is challenging.
for instance, it is hard to produce validated models, or to find core predictions, i.e.
model predictions with low uncertainty
<cit> .
while model validation, in the strict sense, is not possible
<cit> , model rejection and hypothesis testing are possible and highly useful applications in modeling, also for biological research
<cit> .
formally, model rejection methods evaluate whether the null hypothesis,
h <dig>  that a specific model has generated some given data can be rejected or not.
one common way to do this is to test whether the residuals, i.e.
the differences between the simulated and measured data points, are too big.
this can be checked using the χ2-test statistic.
alternatively, one might also wish to check whether the residuals are too correlated.
this can be done using the whiteness test or the durbin-watson  test
<cit> .
however, there is a problem.
these tests are dependent on analytical derivations for the distributions of the test statistic under
h <dig>  but these derivations are based on a number of assumptions, which might not be fulfilled
<cit> .
for instance, some commonly used assumptions are that the experimental noise is normally or log-normally distributed, that the parameter estimates have converged, and that the parameters appear linearly in the model
<cit> .
because many of these assumptions are unfulfilled in systems biology problems, it is problematic to use these analytical expression.
some of the reasons why the assumptions often are unfulfilled include that the availability of data in systems biology examples often is severely limiting, that the signal-to-noise ratio is poor, that the number of parameters that appears non-linearly and/or are unidentifiable often are high, and, for model comparison approaches, such as the likelihood ratio test, that the tested models are not nested
<cit> .
for more information on these assumptions and limitations, we refer the reader to our previous paper
<cit> .
to help overcome the problem of unfulfilled assumptions, one may try to replace the analytical expressions with empirical distributions of the test statistics.
one way to derive the empirical distributions is to use bootstrap approaches.
in general, bootstrap samples are artificially generated data sets, where the distribution of the bootstrap samples should reflect the variability of the data.
although most mathematical proofs for bootstrap approaches usually also are derived under asymptotic conditions, an almost asymptotic setting is often achieved already for moderate sample sizes.
there are two types of bootstrap approaches: parametric and non-parametric
<cit> .
non-parametric bootstrap samples are generated from the original data set by drawing with replacement.
parametric bootstrap samples are generated from a specific model, e.g.
an error model, that also incorporates some null-hypothesis about the underlying system.
there is a rich literature for both parametric and non-parametric methods and their applications to statistical testing in biology
<cit> .
a specific but seemingly unexplored advantage of using a bootstrap setting is that it allows for the natural combination of different test statistics .
this advantage comes because, using bootstrapping, such combined statistical distributions can be calculated empirically, whereas the combination of such distributions analytically largely remains an unresolved problem.
there is a field that deals with the combination of information
<cit> , but this field primarily deals with combinations off different data sources, as in meta-analysis.
for the combination of different statistical tests, one approach that has been considered is to combine the p-values
<cit> .
there are some straightforward simplistic ways in which you could do these combinations.
for instance, given two tests statistics,
ta and
tb, for a specific model and data set, one could look at the maximum or minimum of the respective p-values etc.
<cit> .
this corresponds to the principle of rejecting only if both tests reject, or if at least one of them rejects, respectively.
however, there is a need to evaluate such naive combinations, in general and in the context of systems biology, and to provide more refined alternatives.
in this paper we examine how one could, should, and should not combine test statistics using parametric bootstrapping.
the analysis shows that such above simplistic strategies to combine p-values are unsound, i.e.
the tests are overly conservative or liberal, and in order to avoid this issue, one can instead use our proposed joint two-dimensional distribution of these test statistics.
analysis on simple examples where the truth is known demonstrates the potential gain in power obtained by including an extra dimension.
furthermore, it is found that a beneficial combination may be to combine two χ <dig> statistics, where the second one comes from the ability of a second model to describe data from the first, tested, model.
this leads to a new and more general motivation for the log-likelihood ratio test , in the more general sense of non-nested nonlinear models.
importantly, our 2d approach allows for an easy illustration of when and why a combination of tests is advantageous, and when and why an additional model is helpful.
finally, our results and new approaches are also repeated and illustrated on a previously published example from insulin signaling, with real data, and with a biological question that now is resolved statistically for the first time.
methods
theoretical setting
the herein presented bootstrap approach is applicable to any predictor-based model, i.e.
for any model that can return a vector of predicted outputs
y^, given a vector of parameters θ. importantly, this includes both linear and nonlinear, as well as static and dynamic, models.
nevertheless, since most models in systems biology are based on nonlinear ordinary differential equations
<cit> , we here present the method in an ode-based framework.
in the results section, examples of both linear, nonlinear, static, and dynamic models are included.
let the states in the model be denoted x, and let the time derivatives,
x˙, of these states be governed by a nonlinear smooth function, f. the states, x, usually correspond to the amounts or concentrations of particular molecular compounds.
the function f is usually given by summing up kinetic rate expressions of the involved compounds, assuming mass action laws, or in some cases, more detailed rate equation expressions such as michaelis-menten dynamics
<cit> .
let the function f depend on the states, the parameters, and some input signals, u. let the initial values for the states be denoted x <dig>  and note that these most often are defined as part of the parameter vector θ. finally, let the model outputs,
y^, be determined by a smooth nonlinear function g, which, just like f, may depend on the states, the parameters, and the inputs.
with these notations, the state-space description of the model may be written as:
x˙=f
x=x <dig>
y^=g
the noise, v, is assumed to enter only additively, and only in the measurement equations.
hence, with the measured output denoted as y, the assumption is that
y=y^+ν,ν∼d
for all t, and where ν follows the distribution d. a model,
m, is defined by the specific choice of the functions f and g, i.e.
by the model structure
, and some parameters, θ.
model implementation
all models have been implemented and analyzed using matlab®; r2011b
<cit> .
static models were fitted using standard linear regression methods, such as polyfit.
ode models were implemented using the systems biology toolbox
<cit>  and the add-on package sbao.
parameters of ode models were estimated using the global optimization algorithm simannealingsbao, available in sbtb, by minimizing the χ2-test statistic.
bootstrap setting
a bootstrap sample, b, is an artificial vector of observations.
a set of such bootstrap samples,
bi, is generated with the intent of representing the natural variation of the experimental data set, according to some specific procedures and assumptions.
here we consider parametric bootstrap samples, i.e.
samples that have been generated from a specific model structure, denoted
mi, whose parameters have been fitted to agree with the experimental data.
if nothing else is noted, the default set size, also called cloud size, used in this paper is  <dig> samples per set.
the bootstrap samples are generated by adding noise, drawn from the assumed distribution d, to a simulated output of a given model.
in this paper, the assumed noise distribution is gaussian with a mean of  <dig>  and a standard deviation of  <dig>  and  <dig>  for the static and dynamic case, respectively.
these noise level were chosen to be in the order of 5-10% of the average model output.
conversely, the assumed noise level for the insulin signaling test case corresponds at each time point to the noise in the experimental data.
however, for the first and second time point, where the signal has been normalized, and the noise therefore is zero, an average noise level for the data set is assumed.
each bootstrap sample corresponds to a new realization of the noise with the same measurement signals and time points as in the observed data.
the empirical distribution of any given test statistic, such as the χ <dig> or dw, is obtained by fitting the model of interest to all bootstrap samples, and then for each fit, calculating the appropriate test statistic.
empirical testings and conceptual basis behind the methods
a statistical test is a formal procedure for checking if a null hypothesis, here denoted
h <dig>  can be rejected.
in practice, a test maps a given set of observations, denoted
, to a test statistic,
t. a p-value for a given test statistic is the cumulative probability of that value and all other values that are even more extreme, given the distribution under
h <dig>  where
h <dig> typically corresponds to the hypothesis that the model you are testing is true.
in a bootstrapping environment we construct these distributions empirically, as described in the section bootstrap setting, rather than using analytical expressions for them .
, a model structure,
m <dig>  to be investigated, and two test statistics.
first the model is fitted to the experimental data and then the fitted model is used to generate bootstrap samples.
next, the model is fitted to all bootstrap samples.
the resulting two-dimensional cloud is used to generate a density estimation.
the cutoff plane is defined as the equidensity contour that goes through the coordinates of the experimental data .
the corresponding p-value is equal to the portion of the distribution below the plane, i.e.
outside the corresponding density contour line.
the p-value is then used for decision on whether or not to reject either model.
more specifically, if the null hypothesis that the model is true is correct, then bootstrap samples, generated from the fitted model, actually represent real samples from the true data generating process .
thus, under
h <dig>  the joint distribution of any number of goodness of fit  statistics represents the distribution that we would expect to see if we were able to repeatedly sample from the true dgp, and each time were to fit our model and calculate a corresponding vector of gof statistics.
we can therefore calculate the probability, under
h <dig>  of seeing a vector of gof statistics at least as extreme as the original observed vector of gof statistics, i.e.
a p-value.
these p-values are calculated by estimating the densities of these gof-distributions and then, as described for each test below, evaluate these at the coordinates of the observed data
to obtain the desired p-value.
the null hypothesis is then rejected if the observed vector of test statistics is very unlikely under
h <dig>  usually this amounts to specifying a significance level α and checking whether the observed value is more extreme than the corresponding threshold.
in this paper we consider a seemingly unexplored potential with bootstrap approaches: that they, unlike analytical approaches, allow for an easy calculation of the joint distribution of several test statistics.
consider a vector of k such test statistics, .
given the null hypothesis that the tested model is true, one can then generate data that is assumed to come from the true dgp, and for each data series calculate corresponding values for all these k test statistics.
these values then form vectors of values, and each vector constitute a point in a k-dimensional space.
these points together form an approximation of the true k-dimensional distribution if the null hypothesis is true.
finally, the experimental data series
also corresponds to such a point, and we can use a multi-dimensional density estimation to see whether or not it is realistic to assume that the experimental point lies within the obtained distribution.
this idea can in principle be used for any number of combined test statistics, but the computational cost of approximating a multi-dimensional distribution grows quickly with the dimension.
therefore, we here limit ourselves to combinations of two test statistics, i.e.
to k =  <dig> and to 2d distributions.
one-dimensional tests of a single model
the bootstrapped χ2-test
the χ2-test evaluates the size of the residuals, which are defined as the differences between the measured and predicted outputs:
ri≔yi-y^i
the test statistic,
tχ <dig>  is given by the residual sum of squares
tχ2=∑i,jyi-y^iσi <dig>
where the summation runs over all time points, tj, and all measurement signals, yi.
an empirical distribution of
tχ <dig> is obtained by generating bootstrap samples from a model and fitting this model to the samples, as described above.
the resulting distribution is evaluated using matlab and the empirical cumulative distribution function, ecdf.
a p-value,
pχ <dig>  under
h <dig>  is obtained by considering a right-tailed  distribution and interpolating the value of the cumulative empirical distribution at the coordinate corresponding to the χ2-value of the original data set.
the bootstrapped durbin-watson test
the dw test can be used to test whether the residuals in equation  are too correlated or anti-correlated.
the test statistic,
tdw, is given by
tdw=∑i,j>2ri-ri2∑i,jri <dig> ≤tdw≤ <dig>
the numerator is a measure of the covariance of adjacent residuals, the denominator for the respective variance.
for perfectly uncorrelated residuals the test statistic is equal to two.
a value close to zero indicates a positive correlation, whereas a value close to four indicates a negative correlation.
in this paper we have chosen to look only at correlation, and not at anti-correlation and therefore chosen a left-tailed distribution.
an empirical distribution of
tdw is obtained by generating bootstrap samples from a model and fitting this model to the samples, as described above.
the resulting distribution is evaluated using matlab and the empirical cumulative distribution function, ecdf.
a p-value, pdw, under
h <dig>  is obtained by considering a left-tailed distribution and interpolating the value of the cumulative empirical distribution at the coordinate corresponding to the dw value of the original data set.
simplistic combinations of bootstrapped tests
in this paper, p-values obtained from an empirical 1d χ2-distribution and an empirical 1d dw distribution are combined in various ways.
four of these ways are referred to as simplistic.
these combination tests are defined as follows:
pmin:=minpχ <dig> pdwpmax:=maxpχ <dig> pdwpmean:=pχ2+pdw/2pprod:=pχ2∗pdw
and the interpretations and motivations behind the combinations are described in the corresponding section in results .
two-dimensional tests of a single model
bootstrapped two-dimensional χ2 vs dw test
now follows a description of the novel algorithm for a joint two-dimensional χ2 vs dw test.
although described as a combination of these two specific tests, the method is generalizable to any combination of two test statistics, by simply replacing one or more of the test statistics.
our proposed algorithm consists of the following steps .
algorithm:
given a model,
m1; an experimental data set,
; two test statistics,
tχ <dig> and
tdw; and a significance level α:
<dig>  fit the model to the original data set
and let
θ^m <dig> denote the estimated parameter vector.
calculate the statistics
tχ2m <dig> and
tdwm <dig> according to  and  respectively.
<dig>  use
m <dig> to generate a set of bootstrap samples.
this set is denoted
b <dig>
<dig>  fit
m <dig> to each bootstrap sample, b, in
b <dig> and calculate the corresponding test statistics for each fit.
this results in one set of χ2-values and one set of dw-values, which together form a two-dimensional cloud:
•
c <dig>  cloud  <dig>  consisting of
tχ2m <dig> and
tdwm <dig>
<dig>  estimate  the two-dimensional empirical distribution
ρ <dig>  let
ρ1:=ρ <dig> tdwm1) denote the obtained density at the coordinate corresponding to the χ2 vs dw values of the original data set
.
for the given distribution, we define the cutoff plane as the equidensity contour that goes through
ρ <dig>
<dig>  using the two-dimensional distribution, calculate the p-value for the given model
m <dig>
pχ2-dwm1:=∫ρ1<ρ1ρ1tχ2m <dig> tdwm1dtχ2m1dtdwm <dig>
if
pχ2-dwm1<α, then
m <dig> should be rejected.
two-dimensional density estimation
the two-dimensional density of a cloud is estimated continuously with a smooth gaussian kernel
<cit> , and evaluated over a grid, i.e.
a 2d histogram.
the integral in equation  <dig> is then approximated by summing over all bins.
the total volume is then normalized to one.
two-dimensional p-value calculation
the calculations of p-values in 2d introduces some new considerations, and a few comments are in order.
consider figure
<dig>  and the property a, considered as a 1d distribution.
then, the probably most common way of calculating the p-value is p = 1 - p).
turning to 2d distributions, this formula can no longer be used, since there now are two properties, a and b. instead a more general formula is needed.
one such option is to use some formula based on the probability density function, ρ. then the corresponding formula is p = 1 - p) = p) .
in general, the p-value should give the probability that the obtained value, or an even more extreme one, is found, under the given null hypothesis, and this is ensured by both the 1d and 2d formulas.
note, however, that the 2d formula, using ρ, includes all regions of low density, even the ones where the model is surprisingly good, similar to a two-tailed test in one dimension.
a more detailed discussion on these issues is found in the discussion, and in additional file
1: figure s <dig>
tests involving two models
bootstrapped two-dimensional χ2 vs χ <dig> test
our proposed method for a two-dimensional χ2 vs χ <dig> test is similar to the two-dimensional χ2 vs dw test, where the dw test statistic has been replaced by the χ2-statistic of a second model.
the detailed steps are explained in the additional file
<dig> and in additional file
1: figure s <dig>
bootstrapped log-likelihood ratio test
given some data
, and two models
m <dig> and
m <dig>  an empirical distribution of the lhr,
tlhr, is obtained by generating bootstrap samples from either model  and fitting both models to the samples, as described above.
the resulting distribution of log-likelihoods  are evaluated using matlab and the empirical cumulative distribution function, ecdf.
a p-value, plhr, under
h <dig>  is obtained by considering a two-tailed distribution and interpolating the value of the cumulative empirical distribution at the coordinate corresponding to the lhr value of the original data set.
these steps are explained in detail in the additional file
<dig>
test cases
static models
two static models are considered;
ms <dig>  a straight line, and
ms <dig>  an exponential curve .
ms1̲:f=θs11x+θs12=y^
ms2̲:f=θs21ex+θs22=y^
ms <dig>  and the green and blue lines are corresponding fits from the two model structures.
the two dynamic model structures, involving a decay of the substrate x <dig>  which occurs via mass action or michaelis-menten kinetics, respectively.
vertical lines are data points  from a bootstrap generated by
md <dig>  and the blue and green lines are corresponding model fits by the two model structures.
equations and details regarding fits etc.
are found in materials and methods.
dynamic models
two dynamic non-nested examples are considered;
md <dig>  mass action kinetics, and
md <dig>  michaelis-menten kinetics with one free parameter .
md1̲:x˙1=-θd11x <dig> x1= <dig> y^=x1
md2̲:x˙1=-θd21x <dig> +x <dig> x1= <dig> y^=x1
analyses of methods
the receiver operator characteristic
the power of a statistical test is often determined by the relationship between the false positive rate  and the true positive rate
<cit> .
a false positive is the rejection of a true model, whereas a true positive is the rejection of a false model.
the dependency of the tpr on the fpr is called a receiver operator characteristic  curve.
the more concave the curve, i.e.
the larger the area under the curve , the better the discrimination between true and false models.
here, roc curves are constructed by considering a large number of artificially generated data sets, on which two hypotheses are tested, of which one is the true underlying model.
the obtained p-values for each hypothesis and data set are calculated and for any given fpr  the tpr is obtained.
type i error rate
for a given significance level α, it is expected that  <dig> ·α % of all true values would be rejected.
if the observed fpr is higher than the expected fpr, the test is prone to making type i errors, and is considered liberal.
in contrast, if the observed fpr is lower than the expected fpr, the test is considered conservative.
this method property is evaluated by considering a large number of artificially generated data sets, where the true model is known, and where the calculated p-values thus can be compared to the underlying truth.
any given significance level, i.e.
stated fpr, can thus be compared to the observed fpr, and the resulting relationship can be plotted in a graph .
ideally, the expected fpr should coincide with the observed fpr.
a convex plot would indicate a conservative test, whereas a concave plot would indicate a liberal test.
RESULTS
test cases
in this paper we have employed an ensemble of various bootstrap methods on both static  and dynamic  test case models.
we tested these approaches on static linear examples mainly for the following two reasons: firstly static models are common in science and our methods should therefore aim to be applicable to these kind of problems; secondly the solutions to the corresponding optimization problems are for these static linear examples unique and analytically attainable.
in contrast, we also performed our analysis on two test cases in an ode-framework.
parameters in ode models usually have to be estimated, often by solving non-convex optimization problems, and it is then not guaranteed that the optimal solution will be found.
in both settings, the number of parameters in the test case models were kept low, in order to more reliably evaluate the methods.
in each setting, static and dynamic, we let both models serve as the underlying truth to create  <dig> artificial data sets each.
for each such data set both models served as
h <dig>  and were consecutively fitted to the data, and the goodness of fit  was evaluated using various bootstrap approaches, starting with the simplistic combinations described earlier.
thus, for each setting and bootstrap approach, this resulted in  <dig> true positives and  <dig> true negatives, which were used to construct roc curves to evaluate the power of the tests.
combining χ <dig> and dw statistics
the χ2-test is used to check whether the residuals, i.e.
the differences between the model output and the data points, are too big.
conversely, the dw test is used the check if these residuals are too correlated.
both tests are useful in a model rejection scenario, and in the below analysis, looking at how these tests can be combined, the two tests are used as in their bootstrapped form .
simplistic combinations are unsound
the first part of the analysis concerns the simplistic combinations: pmin, pmax, pmean, and pprod .
although simple, these tests are not without interpretation, and several of them are what at first might seem like the obvious idea
<cit> .
the options min and max corresponds to rejecting if either or if both individual tests reject, respectively.
the mean could be thought of as a balancing between the two extremes, and prod, the product, could be thought of as the joint probability.
all these four simplistic combinations can be discarded based solely on an analysis of their observed and expected type i error rate.
these rates are plotted for all four tests in figure
<dig>  and the interpretations of these curves is as follows.
if the plotted lines lie away from the identity line, the expected false positive rate does not coincide with the observed false positive rate, and if this deviation from the identity line is large we call the test unsound.
a large deviation means one of two things: either the test is liberal , or the test is conservative .
a liberal method is generally regarded as unacceptable, since one wants to be able to trust rejections, but a small level of conservativeness could be accepted, so long as the test is useful.
in both the static  and the dynamic  case, the tested combinations are unsound.
the min  and prod  approaches are strikingly liberal, the max approach is highly conservative , and the mean  switches from below to above.
these plots should be compared to the single tests: χ <dig>  and dw  which lie along the identity line.
this difference between the single tests and the simplistic combinations clearly illustrates that the deviations from the identity line are big.
since the results are essentially the same for both the static and dynamic cases, the results were deemed sufficiently convincing to be able to reject all of the tested simplistic approaches as unsound.
a two-dimensional approach is both sound and informative
the second part of the analysis considers the 2d analysis, combining the χ <dig> and the dw tests .
although a precursor to the 2d methods presented herein has been mentioned and outlined in previous papers
<cit> , this is the first time all the implementation details have been solved, and the performance of the method is tried on examples.
one of the novel details concerns the density estimation.
in contrast to the simplistic combinations, this 2d approach is sound, or only slightly conservative, for both the static  and the dynamic  case.
the conservativeness is tolerated, since the test is informative, as can be gathered from the receiver operator characteristic  curves in figure
5a and c. these roc curves are to be interpreted as follows.
on the x-axis, the rate of erroneous rejections are plotted; this value is therefore to be as small as possible.
on the y-axis, the rate of correct rejections are plotted; this value is therefore to be as high as possible.
taken together, this means that the auc should be as big as possible, especially for the region of interest where the fpr is below  <dig> .
this region of interest was chosen because in biology α is rarely higher than  <dig> .
from figure
5a,c it is clear that the new 2d approach  outperforms both χ <dig>  and dw  considered as individual tests.
<dig>  however, the combination used here  does lie close to the identity line, together with the single test statistics.
the upper plots, , show the results for the static example, and  for the dynamic example.
as can be seen, the 2d χ2 vs dw test is consistent  and superior to both of its constituents tests, which is evident from the greater auc in both the static and the dynamic case .
introducing a second model
replacing the dw test statistic with the χ2-test statistic of a second model
the above 2d approach  can be generalized to combinations of other tests as well.
an important aspect of this is that one can consider a combination of two χ <dig> tests, where the second χ <dig> test comes from a second model,
m <dig> .
it is therefore intuitively sensible to test whether such a usage of two models is an advantageous usage of this 2d approach.
this property of one model’s ability to imitate the behavior of a second model is known as model mimicry, and the idea of using this in a model selection setting has been utilized by e.g.
<cit> .
this second model can in itself be an uninteresting model, i.e.
we are then not necessarily interested in the second model as an explanation of the original data, but only in how it interacts with the model being tested.
such a model is called a help model.
alternatively, the second model could be a competing model and its ability to describe the original data is then of equal importance as that of the first model.
if this latter situation is the case, one would typically perform the analysis with both models serving as
h <dig>  generate bootstrap samples from each model, and so on .
this version of the 2d test then becomes a form of model comparison, even though there are important differences.
for instance, this 2d analysis, unlike model discrimination tests like the conventional non-bootstrapped lhr, can result in all four cases of rejecting none, either, or both of the models.
in contrast, a conventional non-bootstrapped lhr can only result in the rejection of the simpler of the two models, or none of them.
furthermore, in this new setting, the two models do not have to be nested, i.e.
one of the models does not have to be a special case of the other, and the models can be general nonlinear odes.
the results on the same dynamic and static examples as before are plotted in figure
<dig>  as can be seen from figure
6b and d, this 2d χ2 vs χ <dig> method  also has some slight issues with conservativeness, but from figure
6a and c, it is clear that this slight conservativeness should be tolerated: the new 2d version outperforms the previous 2d method .
5:  are roc curves,  are type i error curves,  are for the static, and  for the dynamic example.
as can be seen, all methods are consistent, or weakly conservative.
the roc curves show that the 2d χ2 vs χ <dig> test  is better than the 2d χ2 vs dw test , and that the bootstrapped lhr , is best of all.
in contrast lhp  is completely uninformative, as it lies along the identity line.
the bootstrapped lhr test is the best approach in the case of a good help model
the final test included in this comparison is a bootstrapped version of the lhr .
this method has no issues with conservativeness , and outperforms all the other methods in terms of a roc analysis .
new motivation for the lhr test in the more general case of bootstrapping non-nested models
it is now clear that there are examples where the lhr is advantageous to the 2d χ2 vs χ <dig> analysis; let us now understand why this is the case.
at a first look, it seems like it should be the other way around: that the 2d χ2 vs χ <dig> analysis should be superior, almost by definition.
the reason for this is that the information in the lhr is contained in the 2d χ2 vs χ <dig> analysis.
the close relationship between the two methods can be seen by comparing a cloud from the analysis plotted in the χ2 vs χ <dig> plane , with the same cloud plotted in the lhr vs log-likelihood product  plane .
as can be seen, the shape of the cloud and its relation to the red dot is identical, only tilted  <dig> degrees.
this relation also follows from simple linear algebra.
a further inspection of the cloud in figure
7a, which is taken from the static example, reveals that the cloud is of a character that makes a 2d approach superior to the individual tests: just as the cloud in figure
<dig>  the cloud in figure
7a lies away from the axes, and the red dot can be distinguished better in a 2d analysis, than by looking along the individual axes.
however, when the cloud has been tilted to the lhr vs lhp plane , the red dot can be separated from the cloud when considering only one of the directions: the lhr direction.
that this preservation of the information is preserved when projecting the 2d cloud to the lhr line is also corroborated by the fact that the lhp is, for this example, completely non-informative .
in other words, the 1d lhr has extracted virtually all the relevant information of the 2d cloud.
all this means that if there would not be any price to pay for doing a 2d analysis, the lhr and 2d χ2 vs χ <dig> analysis would be equally good.
however, there is a price to pay for moving to a 2d bootstrapping analysis, and this is the reason why the lhr bootstrapping analysis is superior.
there are several components to this price.
first, the estimation of a distribution is more cumbersome.
second, and more importantly, this 2d estimation converges slower than a corresponding 1d density estimation.
this has to do with the curse of dimensionality, which simply means that the same number of points quickly become more scarcely packed as the dimensions increase, and that a corresponding density estimation will be based on fewer and fewer neighboring points.
this reduction in convergence speed can also be seen in additional file
1: figure s <dig>  where the lhr has converged already with cloud sizes of ∼  <dig> data points, but where the 2d χ2 vs χ <dig> analysis requires at least two more orders of magnitude for its convergence.
finally, there are also other prices of moving to a 2d analysis, such as the inability to clearly define a one-sided or two-sided test .
the existence of such additional factors is also reflected by the fact that the 2d test does not converge to the same roc curve as the lhr test .
choosing the second model
having established that the inclusion of a second help model may improve the statistical power of tests evaluating the first model, leads to the natural question of whether all help models would do.
the answer to this is “no”: if the help model is too simple or too flexible, the advantage is lost, and the resulting model comparison tests - lhr or the 2d χ2 vs χ <dig> analysis - perform worse than the other sound tests presented herein.
as an example of this, consider the completely flexible model, which simply goes through all data points .
if this model,
mflex, is used as the help model instead of the suggested competing model in the static example, the 2d cloud collapses to a line: since the cost of the help model is always zero .
therefore, there is no 2d advantage to make use of, and the lhr distribution will simply be zero minus the χ <dig> distribution , and lhr thus performs as bad as the two-tailed χ <dig> test .
in the additional file
<dig>  another simple help model is considered: a constant model that simply approximates a data-series with its mean value.
here, the picture is a little bit more mixed.
for the static example, the model is too simple, and the two-model tests are in-advantageous .
for the dynamic test case, the constant model does provide some additional information: the 2d χ2 vs χ <dig> analysis performs slightly better, and the lhr test slightly worse, than the single χ2-test .
finally, for all of the above examples with too simple or too flexible models, the 2d χ2 vs χ <dig> analysis is superior to the lhr test, showing that the lhr is more sensitive to the situation of having chosen a bad help model.
application to insulin signaling
as a real modeling example, we used data and models from a previous work
<cit> .
in that paper we analyzed experimental data from insulin signaling in primary human adipocytes.
some of the experimental data are shown in figure
8b.
the data consist of a time series which displays an overshoot: the response increases rapidly from zero, and reaches a maximal value around  <dig> min, and then decreases to an intermediate steady state value.
the observed response is caused by the addition of  <dig> nm insulin at time zero to previously unstimulated fat cells, and the measurements are performed using sds-page and immunoblotting to determine the degree of auto-phosphorylation of the insulin receptor .
the data are normalized such that the first point is zero, and the maximal value is  <dig>  for further details regarding the data, we refer to
<cit> .
using mathematical modeling, we were able to reject several model structures aspiring to explain these data, and we used model selection tools, such as the akaike information criterion
<cit> , on surviving competing hypotheses.
mi,c and the chosen help model
mi,b.
this example was analyzed previously
<cit>  and is of interest, since the rejection of
mi,c would suggest that the recycling is a necessary mechanism to explain the data.
depicted in panel  is the experimental data z , and fits of
mi,b  and
mi,c .
the measured data represent the increased auto-phosphorylation of the insulin receptor in response to  <dig> nm insulin in isolated primary human adipocytes, as measured by sds-page and immunoblotting.
panel  shows the bootstrapping cloud in the χ2 vs dw plane, when the bootstraps have been generated by
mi,c .
as can be seen, the cloud lies along the axes, and there is no benefit of using a 2d analysis.
panel  shows the χ2 vs χ <dig> scatter plot of the
mi,c cloud  generated after fitting both models to bootstrap sets from
mi,c.
the corresponding χ2-values for the experimental data  is also plotted .
as can be seen, the cloud lies away from the axis, and the experimental data point explores the obtained direction.
panel  summarizes the results.
as the clouds have indicated, the χ2 vs dw combination does not improve upon the individual tests, but still lies on the border of rejection.
the χ2 vs χ <dig> tests on the other hand perform better than the individual tests, and the lhr is best of all.
here, we re-examined one of the models,
mi,c, from
<cit> , that could previously not be unambiguously rejected.
the model structure of
mi,c and the chosen help model,
mi,b, are included also in this paper for convenience .
mi,b, the smaller model, contains only three reactions: insulin binding with auto-phosphorylation of the receptor, internalization with dephosphorylation, and recycling back to the plasma membrane.
mi,b fits to the data with a distinct overshoot, reaching an intermediate quasi-steady state after a couple of minutes .
mi,c, on the other hand, is considerably more detailed in terms of the binding, auto-phosphorylation, and internalization, but it does not contain any recycling.
mi,c fits to the data in a qualitatively very different way .
mi,c has its maximum almost immediately after stimulation, but thereafter never reaches any quasi-steady state as
mi,b does.
instead the signal declines almost linearly during the observed time period.
this example is of biological interest, since the rejection of
mi,c would suggest that the recycling is a necessary mechanism to explain the data.
this conclusion would follow because
mi,c is a more complex, and a more realistic interpretation of the underlying biological system in all other aspects except recycling.
in
<cit> , an aic analysis and a χ2-square test were performed on the various competing models and although
mi,c was discarded based on its lack of agreement with data and the fact that inclusion of a recycling reaction yielded a better model, neither of these results were statistically convincing.
indeed, although the results pointed towards the rejection of
mi,c, it was pointed out already in
<cit>  that a more accurate analysis would require a more generally applicable method such as bootstrapping.
in the re-analysis of this insulin signaling example, the analysis was done using all bootstrapped tests that have been found sound on the simpler test cases, where the truth was known.
unlike in those examples, here the truth is not known, but subsequent experimental analysis of the system has revealed that
mi,c indeed should be rejected.
the results are summarized in figure
8e.
as can be seen, a bootstrapped χ2-test yielded p =  <dig> , and a bootstrapped  dw test yielded p =  <dig> , both on the border of rejection.
a 2d χ2 vs dw analysis did in this case not improve this value because the resulting cloud and data point  did not lie in an informative configuration , p =  <dig> .
conversely, for the 2d χ2 vs χ <dig> the cloud and data point did lie in an informative configuration .
as a result, the two best methods on the test cases, the 2d χ2 vs χ <dig> and the bootstrapped lhr showed improved performance as compared to the above tests, yielding p-values that were clearly below the threshold , p = 8∗10- <dig> and p ∼  <dig> respectively.
also, the internal order of these methods was preserved: lhr was better than the 2d χ2 vs χ <dig>  these findings are the same as in all tested examples with a sufficiently good help model, and speaks for their generality.
discussion
in this paper we have extended the traditional parametric bootstrap approach to also look at combinations of different test statistics, here exemplified by the χ2-test statistic and the dw test statistic.
we have shown how simplistic combinations, as considered in the literature, of these two statistics are unsound, but that a two-dimensional empirical distribution, as enabled by the bootstrap setting, is both sound and informative.
we have further shown that it is even more informative to replace the dw statistic with the χ2-statistic from a second model, resulting in a two-dimensional χ2 vs χ <dig> test.
however, the one-dimensional bootstrapped lhr is found to be even better, and an analysis of this has led to a new motivation and understanding of the lhr, in the more general case of nonlinear and non-nested models.
a 2d bootstrap approach may be superior to a corresponding 1d approach, but there is in fact a trade-off between negative and positive aspects.
one positive aspect is of course that a 2d analysis contains more information than a 1d analysis, since e.g.
the insights contained in two statistics can be combined.
second, even though one in principle can combine test statistics in 1d as well, figure
<dig> shows that the simplistic combinations give inconsistent results, and therefore are unacceptable.
conversely, a 2d-combination of the same test statistics is sound .
a third positive aspect of 2d combinations is that they have a higher power than either of the 1d tests including only one of the test statistics .
all these positive aspects are due to the benefit illustrated in figure
<dig>  where it is clear that only the 2d combination of the test statistics reveal that the indicated point lies outside the empirical distributions.
there are, however, also several negative complications inherent in the addition of an extra dimension, and these complications are solvable to different extents.
the first complication concerns the more complicated density estimation that is required in 2d, but this has been resolved in the herein presented 2d approaches.
the second complication stems from the relative scarcity of neighboring points in two dimensions and its effect on the density approximation.
if the same number of points is used to estimate a 1d distribution and a corresponding 2d distribution, the end result will be a worse approximation for the 2d distribution, simply because of the curse of dimensionality.
this second complication can probably be completely overcome by adding more points to the estimation of the empirical distribution, but the addition of more points comes at the price of a higher computational cost.
the third complication is that a right-tailed, left-tailed, or two-tailed tolerance interval cannot be clearly defined in two dimensions.
this issue is described more below.
finally, the positive aspects of the additional dimension only appears if the resulting cloud does not lie along one of the axis, but as in figure
<dig>  in a direction that is not parallel to the axis.
all in all, this means that the advantages of a switch to a 2d approach are not guaranteed, while the negative consequences are guaranteed.
for this reason the choice between 1d and 2d is a trade-off.
similarly to the issue of the second dimension, our results seem to indicate that the addition of a second model provides an advantage, compared to analyzing a single model based only on its own residuals, but there is in fact also here a trade-off between positive and negative aspects.
the positive aspects are clear from the examples where we use a help model that is roughly equally good as the tested model: then the roc curves have a higher area under the curve , and because the model that should be rejected in the insulin signaling example gives a lower p-value .
these positive sides mean that the additional information provided by a second model provides valuable knowledge regarding the quality of the first model.
more specifically, this arguably means that the second dimension from the second models in the examples implies 2d clouds that do not lie along any of the axes, but, as in figure
<dig>  in a direction that is not parallel to the axes.
a final positive aspect of this bootstrapped two-model approach is that it goes beyond the capacity of a normal model comparison test, e.g.
lhr, since the resulting conclusion can be that none, either, or both models are rejected.
a standard model comparison test can never reject both models.
the negative aspects of adding a second model are of two types: i) those that have to do with a 2d approach, and which are described above, and ii) those that are associated to the fact that not all additional models provide an improvement.
as is clear from e.g.
figure
7f, it is clear that a poorly chosen model yields a worse test compared to a mere χ2-test.
the second negative aspect is that it is for the non-nested cases not possible to know in advance when a model is good or poor help model.
here it should be re-stated that the 2d χ2 vs χ <dig> test is more robust towards bad help models than the lhr test in all examples tested herein.
in summary, a help model should not be too flexible or too inflexible, and one can see whether the right complexity of the help model has been struck from the shape of the cloud: if it lies along one of the axes it is too flexible or too inflexible.
one of the negative aspects mentioned above needs a little further clarification: the definition of the empirical tolerance intervals, which are used to calculate the empirical p-values.
first, the 1d versions of the methods that we are using either operate in a one-sided way , or a in a two-sided way .
there is no obvious translation of sides and tails in a 1d distribution, to a 2d equivalent.
we here adopt the definition of the 2d tolerance region as the region with highest probability density .
in practice this is similar to a two-sided cutoff since one may reject a model because it is unrealistically good at describing the data, compared to the given noise level.
however, there are differences, such as the possibility to have several fragmented regions instead of a single joint one.
therefore, when comparing our method with a 1d-version of the same test, one could consider defining the 1d tolerance interval in a likewise manner , since this more closely mimics the cut-off we do in 2d.
however, all comparisons here are done with the one-sided or two-sided 1d-versions of the methods, since it is those that are used in practice, and those that our method should out-compete.
a key validation step, and a demonstration of the usefulness of the results herein, is that they are corroborated on a real-world example, which now has been resolved statistically for the first time: we can now reject
mi,c with a p-value that lies clearly below the threshold.
we have confidence that this rejection of
mi,c is correct, because in more recent works we have experimentally shown, by blocking internalization and measuring a downstream component, that recycling of the receptor does play a major role in insulin signaling
<cit> .
we have also measured the amount of internalized insulin receptor and shown that there is far too little internalized ir to explain the observed overshoot.
however, even though we in light of these new data have rejected this model, it is interesting to observe that it was possible, with our new method, to reject
mi,c based only on the data available at the time of
<cit> .
there are some limitations when interpreting the results that should be mentioned.
first, the results are only shown for three example cases, and there is no guarantee that they hold for all other examples.
nevertheless, the results are essentially the same for all these three examples: 2d is better than 1d for the static and dynamic examples, and for all examples the tests with a non-extreme help model are better than the single rejection tests, and lhr is best of all.
therefore, since the examples include static, dynamic, linear, nonlinear, and real-world aspects, these overall trends probably have some generality.
second, the generality is also limited by the fact that we do not have analytical proofs for the results.
this, however, is normal for bootstrap approaches.
third, another limitation with this approach is that it only considers the usage of a single help model or help statistic.
however, this is not a major limitation, since we anyway only advocate the usage of these methods in special cases, where power and accuracy, rather than computational speed, are desired.
in other words, our suggestion is to use this approach only in cases where you have a specific model that requires a more detailed analysis.
also, it is unlikely that a generalization of this approach to 3d would be beneficial, since then the price of estimating density in a high-dimensional space, and the corresponding slower convergence due to the curse of dimensionality, would be even more severe.
it is important to put our results in relation to the existing literature in related fields, such as statistics
<cit> , systems biology
<cit> , econometrics
<cit> , mathematical psychology
<cit> , phylogenetics
<cit>  etc.
first, our method is advantageous only in cases where you have a particularly important and challenging rejection case, where computational time is not a big issue.
this stands in contrast to the typical situation of aic and bayesian information criterion , where a big number of models can be sorted through a simple criterion
<cit> .
similarly, the herein presented methods are not suitable to the sub-field of optimal experimental design for the purpose of improved model rejection, since such studies requires an optimization over different experimental designs, which in turn mandates less computationally heavy approaches
<cit> .
second, the perhaps most commonly used method for model rejection, the χ2-test, has a problem - that the degrees of freedom in the χ <dig> distribution usually is unknown
<cit>  - but this problem is overcome by using the methods considered herein.
however, this is not a new result, but is true for all bootstrap approaches, characterizing the distribution empirically.
third, there are a number of commonly used test statistics that we have not considered
<cit> .
for instance, as an alternative to the dw test to measure correlation among the residuals, the whiteness and the run test may be used.
it is still unknown whether our results for how to combine test statistics holds also for these other tests.
the final, and now following, two comparisons with literature have to do with the lhr and with bayesian approaches.
the herein presented analysis presents a new way of showing why and when the lhr is advantageous compared to an individual χ2-test, for the more general setting of nonlinear and non-nested models.
since lhr has been both extensively used and studied in the literature, it is important to relate this claim to previous results.
lhr was first advocated by neyman and pearson in 1928
<cit> , and the basic χ <dig> distribution relation for nested linear models was known already in the 30’s
<cit> .
these results were generalized to the important case of nonnested models by cox in
<cit>  and to the case of neither of the competing models being true by vuong
<cit> .
however, these results are of limited use, since they rely on analytical derivations of mean and variance terms
<cit> , and the results by vuong do not even apply to time-series models
<cit> .
also note that there are important cases where the traditional likelihood ratio test is not applicable, e.g.
for usage in stochastic models based on the chemical master equation.
all of these limitations can be avoided, by adopting a bootstrap approach.
this approach basically only relies on the ability to draw bootstrap samples in a way that approximates the true data gathering process.
the simulation based type of bootstrap approaches studied herein was originally proposed by williams et al.
<cit> .
the williams approach has all the essential qualities of how we implement the bootstrapped lhr herein: both models are fitted to the data, and the fitted parameters are used to generate bootstrap samples that explicitly incorporates the null hypothesis that the used model is true, and finally both models are fitted to all bootstraps and corresponding distributions of lhr values are calculated.
this approach has also been widely used using minor modifications
<cit> , including an approach where the bootstrap samples are generated using drawing with replacement of the residuals
<cit> .
there are also some papers where theoretical properties of the williams approach have been investigated.
for instance
<cit> , shows that the empirical distribution of williams asymptotically converges to the correct distribution under certain conditions.
however, none of those papers use a 2d approach such as ours to achieve an intuitive understanding for why the lhr may be advantageous: that it incorporates the potential added value of the 2d approach compared to the individual χ2-tests, without paying the price of a 2d density estimation.
the most important and novel part herein is perhaps that it allows the user to quickly check whether and why the bootstrapped lhr is advantageous or disadvantageous to use compared to the individual χ2-test: it depends on whether the second χ2-test yields a cloud that lies away from being parallel to the axes, which in turn requires that the help model is neither too simple, nor too complex .
the final important comparison with literature concerns that with bayesian approaches.
bayesian approaches are centered around the combination of a prior distribution or belief with experimental data to obtain a posterior distribution.
although bayesian calculations in practice can be done using simple calculations like the bic, the perhaps most common approach involves markov chain monte carlo
<cit> , and such calculations have big similarities to bootstrap approaches.
one important such mcmc-based approach, which has been used in systems biology, is the bayes factor
<cit> .
bf can be viewed as a generalization of the lhr to a bayesian setting.
in particular, this means that the likelihoods are integrated over the prior distributions of the parameters, to obtain the ratio of the marginalized distributions.
methods to do these marginalizations have been investigated in e.g.
<cit> , and simplified calculations, not requiring the likelihood, using approximate bayesian computations , are considered e.g.
in
<cit> .
this inclusion of the parameter uncertainties is important, because in systems biology the parameters are often undetermined
<cit> , and an important continuation of this work will therefore be to compare the lhr with future extensions of the herein presented frequentist approaches to also include parameter uncertainty.
on this note, it should be mentioned that we have done a small analysis to see the effect of the inclusion of such parameter uncertainties on a specific cloud by exploiting the profile likelihood
<cit> .
this small scale analysis indicates that although the results may change upon such an inclusion, the change is not big compared to other uncertainties within the method.
another way to take parameter uncertainty into account is by introducing an initial step of non-parametric bootstrapping into the bootstrap sample generation, as done e.g.
in
<cit> .
once parameter uncertainty is taken into account in the bootstrap setting in this way, there is a striking similarity to the bayesian posterior predictive  checks
<cit> .
in bpp, the predicted distribution of future experiments is compared to the observed data.
this is done by generating new data by sampling and simulating from the posterior distribution, and then comparing the resulting distribution of goodness-of-fit  with the gof from the experimental data
<cit> .
with all these similarities pointed out, it should also be recalled that bayesian approaches are never identical to frequentist approaches, since frequentist approaches do not require a prior.
CONCLUSIONS
in a bootstrap setting, it is possible to obtain joint distributions for combinations of test statistics in a more straightforward way than is possible in an analytical setting, but this possibility has previously been little explored.
we here show that such combinations often do provide additional knowledge not contained in the individual tests, but that the considered simplistic combinations, like max and min, yield inconsistent, i.e.
overly conservative or liberal, results .
a new 2d approach , on the other hand, is only mildly conservative , and is superior to the individual tests .
these results were obtained on both a static and dynamic case, where the truth is known .
on the same examples, a 2d χ2 vs χ <dig> test is superior to a 2d χ2 vs dw test , where the additional χ2-value comes from the ability of a second model to describe bootstrap samples from the tested model .
the 2d χ2 vs χ <dig> test is, in turn, outperformed by the 1d bootstrapped lhr .
these results are also confirmed on a previously published rejection example from insulin signaling in human fat cells, which has now been statistically resolved for the first time .
further analysis of these results show that whether or not a 2d combination is advantageous depends on a balancing between positive and negative aspects.
the positive sides are found if the cloud as in figure
<dig> lies in a direction non-parallel to either of the axes, and the price to exploit this is e.g.
that density estimation in 2d converges more slowly, and that one cannot define one-sided or two-sided cutoffs for the tolerance regions.
similarly, the additional model only provides a benefit if it is of a rightly balanced ability to describe the data; otherwise using the additional model worsens the performance.
it is because of these balancing acts between positive and negative aspects that lhr may be the better choice: if the additional model is of appropriate complexity, lhr extracts all the useful information of the χ2 vs χ <dig> plot with a one-dimensional analysis, which thus avoids the estimation of a 2d density .
this analysis thus provides a new motivation for the lhr test, which is valid in the case of non-linear and non-nested models.
in summary, these results provide useful insights into the important systems biology problem of model rejection: when to use, and when not to use, 2d approaches and additional models.
these methods are to be exploited in challenging and important cases, when accuracy and power rather than computational speed are prioritized.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
rj developed, implemented, and analyzed the methods.
rj also constructed the test cases, and drafted the manuscript.
ps contributed with biological models, experimental data, and interpretations.
gc was the main designer of the study and co-drafted the manuscript.
all authors read and approved the final manuscript.
supplementary material
additional file 1
supplementary material.
the pdf-file “rj et al.
- supplementary material.pdf” contains additional elaboration on methods used, some additional results and analysis, and discussion on some of the issues covered in this manuscript, as well as all supplementary figures referred to.
click here for file
acknowledgements
this work was supported by the european commission through the biosim project, by the lions sweden research funding, by the swedish research council, and by the german federal ministry of education and research virtual liver grant 0313074d.
the authors also wish to thank clemens kreutz for crucial input on statistical analysis techniques, input regarding the development of the method, and for comments on a previous version of the manuscript.
the authors also want to thank jens timmer for input on the general study design and for financial support.

 a key tool in systems biology is mathematical modeling. modeling allows for a more complete analysis of the true relationship between experimental data and possible mechanistic explanations, compared to what is feasible using only classical biochemical reasoning. nevertheless, because the data are limited and the systems are highly complex, and because many of the model parameters have to be estimated and cannot be uniquely determined, drawing mechanistic conclusions from modeling is challenging. for instance, it is hard to produce validated models, or to find core predictions, i.e. model predictions with low uncertainty. while model validation, in the strict sense, is not possible, model rejection and hypothesis testing are possible and highly useful applications in modeling, also for biological research. formally, model rejection methods evaluate whether the null hypothesis, h that a specific model has generated some given data can be rejected or not. one common way to do this is to test whether the residuals, i.e. the differences between the simulated and measured data points, are too big. this can be checked using the χ2-test statistic. alternatively, one might also wish to check whether the residuals are too correlated. this can be done using the whiteness test or the durbin-watson test. however, there is a problem. these tests are dependent on analytical derivations for the distributions of the test statistic under h but these derivations are based on a number of assumptions, which might not be fulfilled. for instance, some commonly used assumptions are that the experimental noise is normally or log-normally distributed, that the parameter estimates have converged, and that the parameters appear linearly in the model. because many of these assumptions are unfulfilled in systems biology problems, it is problematic to use these analytical expression. some of the reasons why the assumptions often are unfulfilled include that the availability of data in systems biology examples often is severely limiting, that the signal-to-noise ratio is poor, that the number of parameters that appears non-linearly and/or are unidentifiable often are high, and, for model comparison approaches, such as the likelihood ratio test, that the tested models are not nested. for more information on these assumptions and limitations, we refer the reader to our previous paper. to help overcome the problem of unfulfilled assumptions, one may try to replace the analytical expressions with empirical distributions of the test statistics. one way to derive the empirical distributions is to use bootstrap approaches. in general, bootstrap samples are artificially generated data sets, where the distribution of the bootstrap samples should reflect the variability of the data. although most mathematical proofs for bootstrap approaches usually also are derived under asymptotic conditions, an almost asymptotic setting is often achieved already for moderate sample sizes. there are two types of bootstrap approaches: parametric and non-parametric. non-parametric bootstrap samples are generated from the original data set by drawing with replacement. parametric bootstrap samples are generated from a specific model, e.g. an error model, that also incorporates some null-hypothesis about the underlying system. there is a rich literature for both parametric and non-parametric methods and their applications to statistical testing in biology. a specific but seemingly unexplored advantage of using a bootstrap setting is that it allows for the natural combination of different test statistics. this advantage comes because, using bootstrapping, such combined statistical distributions can be calculated empirically, whereas the combination of such distributions analytically largely remains an unresolved problem. there is a field that deals with the combination of information, but this field primarily deals with combinations off different data sources, as in meta-analysis. for the combination of different statistical tests, one approach that has been considered is to combine the p-values. there are some straightforward simplistic ways in which you could do these combinations. for instance, given two tests statistics, ta and tb, for a specific model and data set, one could look at the maximum or minimum of the respective p-values etc.. this corresponds to the principle of rejecting only if both tests reject, or if at least one of them rejects, respectively. however, there is a need to evaluate such naive combinations, in general and in the context of systems biology, and to provide more refined alternatives. in this paper we examine how one could, should, and should not combine test statistics using parametric bootstrapping. the analysis shows that such above simplistic strategies to combine p-values are unsound, i.e. the tests are overly conservative or liberal, and in order to avoid this issue, one can instead use our proposed joint two-dimensional distribution of these test statistics. analysis on simple examples where the truth is known demonstrates the potential gain in power obtained by including an extra dimension. furthermore, it is found that a beneficial combination may be to combine two χ statistics, where the second one comes from the ability of a second model to describe data from the first, tested, model. this leads to a new and more general motivation for the log-likelihood ratio test, in the more general sense of non-nested nonlinear models. importantly, our 2d approach allows for an easy illustration of when and why a combination of tests is advantageous, and when and why an additional model is helpful. finally, our results and new approaches are also repeated and illustrated on a previously published example from insulin signaling, with real data, and with a biological question that now is resolved statistically for the first time. methods theoretical setting the herein presented bootstrap approach is applicable to any predictor-based model, i.e. for any model that can return a vector of predicted outputs y^, given a vector of parameters θ. importantly, this includes both linear and nonlinear, as well as static and dynamic, models. nevertheless, since most models in systems biology are based on nonlinear ordinary differential equations, we here present the method in an ode-based framework. in the results section, examples of both linear, nonlinear, static, and dynamic models are included. let the states in the model be denoted x, and let the time derivatives, x˙, of these states be governed by a nonlinear smooth function, f. the states, x, usually correspond to the amounts or concentrations of particular molecular compounds. the function f is usually given by summing up kinetic rate expressions of the involved compounds, assuming mass action laws, or in some cases, more detailed rate equation expressions such as michaelis-menten dynamics. let the function f depend on the states, the parameters, and some input signals, u. let the initial values for the states be denoted x and note that these most often are defined as part of the parameter vector θ. finally, let the model outputs, y^, be determined by a smooth nonlinear function g, which, just like f, may depend on the states, the parameters, and the inputs. with these notations, the state-space description of the model may be written as: x˙=f x=x y^=g the noise, v, is assumed to enter only additively, and only in the measurement equations. hence, with the measured output denoted as y, the assumption is that y=y^+ν,ν∼d for all t, and where ν follows the distribution d. a model, m, is defined by the specific choice of the functions f and g, i.e. by the model structure, and some parameters, θ. model implementation all models have been implemented and analyzed using matlab®; r2011b. static models were fitted using standard linear regression methods, such as polyfit. ode models were implemented using the systems biology toolbox and the add-on package sbao. parameters of ode models were estimated using the global optimization algorithm simannealingsbao, available in sbtb, by minimizing the χ2-test statistic. bootstrap setting a bootstrap sample, b, is an artificial vector of observations. a set of such bootstrap samples, bi, is generated with the intent of representing the natural variation of the experimental data set, according to some specific procedures and assumptions. here we consider parametric bootstrap samples, i.e. samples that have been generated from a specific model structure, denoted mi, whose parameters have been fitted to agree with the experimental data. if nothing else is noted, the default set size, also called cloud size, used in this paper is samples per set. the bootstrap samples are generated by adding noise, drawn from the assumed distribution d, to a simulated output of a given model. in this paper, the assumed noise distribution is gaussian with a mean of and a standard deviation of and for the static and dynamic case, respectively. these noise level were chosen to be in the order of 5-10% of the average model output. conversely, the assumed noise level for the insulin signaling test case corresponds at each time point to the noise in the experimental data. however, for the first and second time point, where the signal has been normalized, and the noise therefore is zero, an average noise level for the data set is assumed. each bootstrap sample corresponds to a new realization of the noise with the same measurement signals and time points as in the observed data. the empirical distribution of any given test statistic, such as the χ or dw, is obtained by fitting the model of interest to all bootstrap samples, and then for each fit, calculating the appropriate test statistic. empirical testings and conceptual basis behind the methods a statistical test is a formal procedure for checking if a null hypothesis, here denoted h can be rejected. in practice, a test maps a given set of observations, denoted, to a test statistic, t. a p-value for a given test statistic is the cumulative probability of that value and all other values that are even more extreme, given the distribution under h where h typically corresponds to the hypothesis that the model you are testing is true. in a bootstrapping environment we construct these distributions empirically, as described in the section bootstrap setting, rather than using analytical expressions for them., a model structure, m to be investigated, and two test statistics. first the model is fitted to the experimental data and then the fitted model is used to generate bootstrap samples. next, the model is fitted to all bootstrap samples. the resulting two-dimensional cloud is used to generate a density estimation. the cutoff plane is defined as the equidensity contour that goes through the coordinates of the experimental data. the corresponding p-value is equal to the portion of the distribution below the plane, i.e. outside the corresponding density contour line. the p-value is then used for decision on whether or not to reject either model. more specifically, if the null hypothesis that the model is true is correct, then bootstrap samples, generated from the fitted model, actually represent real samples from the true data generating process. thus, under h the joint distribution of any number of goodness of fit statistics represents the distribution that we would expect to see if we were able to repeatedly sample from the true dgp, and each time were to fit our model and calculate a corresponding vector of gof statistics. we can therefore calculate the probability, under h of seeing a vector of gof statistics at least as extreme as the original observed vector of gof statistics, i.e. a p-value. these p-values are calculated by estimating the densities of these gof-distributions and then, as described for each test below, evaluate these at the coordinates of the observed data to obtain the desired p-value. the null hypothesis is then rejected if the observed vector of test statistics is very unlikely under h usually this amounts to specifying a significance level α and checking whether the observed value is more extreme than the corresponding threshold. in this paper we consider a seemingly unexplored potential with bootstrap approaches: that they, unlike analytical approaches, allow for an easy calculation of the joint distribution of several test statistics. consider a vector of k such test statistics,. given the null hypothesis that the tested model is true, one can then generate data that is assumed to come from the true dgp, and for each data series calculate corresponding values for all these k test statistics. these values then form vectors of values, and each vector constitute a point in a k-dimensional space. these points together form an approximation of the true k-dimensional distribution if the null hypothesis is true. finally, the experimental data series also corresponds to such a point, and we can use a multi-dimensional density estimation to see whether or not it is realistic to assume that the experimental point lies within the obtained distribution. this idea can in principle be used for any number of combined test statistics, but the computational cost of approximating a multi-dimensional distribution grows quickly with the dimension. therefore, we here limit ourselves to combinations of two test statistics, i.e. to k =  and to 2d distributions. one-dimensional tests of a single model the bootstrapped χ2-test the χ2-test evaluates the size of the residuals, which are defined as the differences between the measured and predicted outputs: ri≔yi-y^i the test statistic, tχ is given by the residual sum of squares tχ2=∑i,jyi-y^iσi where the summation runs over all time points, tj, and all measurement signals, yi. an empirical distribution of tχ is obtained by generating bootstrap samples from a model and fitting this model to the samples, as described above. the resulting distribution is evaluated using matlab and the empirical cumulative distribution function, ecdf. a p-value, pχ under h is obtained by considering a right-tailed distribution and interpolating the value of the cumulative empirical distribution at the coordinate corresponding to the χ2-value of the original data set. the bootstrapped durbin-watson test the dw test can be used to test whether the residuals in equation are too correlated or anti-correlated. the test statistic, tdw, is given by tdw=∑i,j>2ri-ri2∑i,jri ≤tdw≤ the numerator is a measure of the covariance of adjacent residuals, the denominator for the respective variance. for perfectly uncorrelated residuals the test statistic is equal to two. a value close to zero indicates a positive correlation, whereas a value close to four indicates a negative correlation. in this paper we have chosen to look only at correlation, and not at anti-correlation and therefore chosen a left-tailed distribution. an empirical distribution of tdw is obtained by generating bootstrap samples from a model and fitting this model to the samples, as described above. the resulting distribution is evaluated using matlab and the empirical cumulative distribution function, ecdf. a p-value, pdw, under h is obtained by considering a left-tailed distribution and interpolating the value of the cumulative empirical distribution at the coordinate corresponding to the dw value of the original data set. simplistic combinations of bootstrapped tests in this paper, p-values obtained from an empirical 1d χ2-distribution and an empirical 1d dw distribution are combined in various ways. four of these ways are referred to as simplistic. these combination tests are defined as follows: pmin:=minpχ pdwpmax:=maxpχ pdwpmean:=pχ2+pdw/2pprod:=pχ2∗pdw and the interpretations and motivations behind the combinations are described in the corresponding section in results. two-dimensional tests of a single model bootstrapped two-dimensional χ2 vs dw test now follows a description of the novel algorithm for a joint two-dimensional χ2 vs dw test. although described as a combination of these two specific tests, the method is generalizable to any combination of two test statistics, by simply replacing one or more of the test statistics. our proposed algorithm consists of the following steps. algorithm: given a model, m1; an experimental data set,; two test statistics, tχ and tdw; and a significance level α: fit the model to the original data set and let θ^m denote the estimated parameter vector. calculate the statistics tχ2m and tdwm according to and respectively. use m to generate a set of bootstrap samples. this set is denoted b fit m to each bootstrap sample, b, in b and calculate the corresponding test statistics for each fit. this results in one set of χ2-values and one set of dw-values, which together form a two-dimensional cloud:  • c cloud consisting of tχ2m and tdwm estimate the two-dimensional empirical distribution ρ let ρ1:=ρ tdwm1) denote the obtained density at the coordinate corresponding to the χ2 vs dw values of the original data set. for the given distribution, we define the cutoff plane as the equidensity contour that goes through ρ using the two-dimensional distribution, calculate the p-value for the given model m pχ2-dwm1:=∫ρ1<ρ1ρ1tχ2m tdwm1dtχ2m1dtdwm if pχ2-dwm1<α, then m should be rejected. two-dimensional density estimation the two-dimensional density of a cloud is estimated continuously with a smooth gaussian kernel, and evaluated over a grid, i.e. a 2d histogram. the integral in equation is then approximated by summing over all bins. the total volume is then normalized to one. two-dimensional p-value calculation the calculations of p-values in 2d introduces some new considerations, and a few comments are in order. consider figure and the property a, considered as a 1d distribution. then, the probably most common way of calculating the p-value is p = 1 - p). turning to 2d distributions, this formula can no longer be used, since there now are two properties, a and b. instead a more general formula is needed. one such option is to use some formula based on the probability density function, ρ. then the corresponding formula is p = 1 - p) = p). in general, the p-value should give the probability that the obtained value, or an even more extreme one, is found, under the given null hypothesis, and this is ensured by both the 1d and 2d formulas. note, however, that the 2d formula, using ρ, includes all regions of low density, even the ones where the model is surprisingly good, similar to a two-tailed test in one dimension. a more detailed discussion on these issues is found in the discussion, and in additional file 1: figure s tests involving two models bootstrapped two-dimensional χ2 vs χ test our proposed method for a two-dimensional χ2 vs χ test is similar to the two-dimensional χ2 vs dw test, where the dw test statistic has been replaced by the χ2-statistic of a second model. the detailed steps are explained in the additional file and in additional file 1: figure s bootstrapped log-likelihood ratio test given some data, and two models m and m an empirical distribution of the lhr, tlhr, is obtained by generating bootstrap samples from either model and fitting both models to the samples, as described above. the resulting distribution of log-likelihoods are evaluated using matlab and the empirical cumulative distribution function, ecdf. a p-value, plhr, under h is obtained by considering a two-tailed distribution and interpolating the value of the cumulative empirical distribution at the coordinate corresponding to the lhr value of the original data set. these steps are explained in detail in the additional file test cases static models two static models are considered; ms a straight line, and ms an exponential curve. ms1̲:f=θs11x+θs12=y^ ms2̲:f=θs21ex+θs22=y^ ms and the green and blue lines are corresponding fits from the two model structures. the two dynamic model structures, involving a decay of the substrate x which occurs via mass action or michaelis-menten kinetics, respectively. vertical lines are data points from a bootstrap generated by md and the blue and green lines are corresponding model fits by the two model structures. equations and details regarding fits etc. are found in materials and methods. dynamic models two dynamic non-nested examples are considered; md mass action kinetics, and md michaelis-menten kinetics with one free parameter. md1̲:x˙1=-θd11x x1= y^=x1 md2̲:x˙1=-θd21x +x x1= y^=x1 analyses of methods the receiver operator characteristic the power of a statistical test is often determined by the relationship between the false positive rate and the true positive rate. a false positive is the rejection of a true model, whereas a true positive is the rejection of a false model. the dependency of the tpr on the fpr is called a receiver operator characteristic curve. the more concave the curve, i.e. the larger the area under the curve, the better the discrimination between true and false models. here, roc curves are constructed by considering a large number of artificially generated data sets, on which two hypotheses are tested, of which one is the true underlying model. the obtained p-values for each hypothesis and data set are calculated and for any given fpr the tpr is obtained. type i error rate for a given significance level α, it is expected that ·α % of all true values would be rejected. if the observed fpr is higher than the expected fpr, the test is prone to making type i errors, and is considered liberal. in contrast, if the observed fpr is lower than the expected fpr, the test is considered conservative. this method property is evaluated by considering a large number of artificially generated data sets, where the true model is known, and where the calculated p-values thus can be compared to the underlying truth. any given significance level, i.e. stated fpr, can thus be compared to the observed fpr, and the resulting relationship can be plotted in a graph. ideally, the expected fpr should coincide with the observed fpr. a convex plot would indicate a conservative test, whereas a concave plot would indicate a liberal test.  test cases in this paper we have employed an ensemble of various bootstrap methods on both static and dynamic test case models. we tested these approaches on static linear examples mainly for the following two reasons: firstly static models are common in science and our methods should therefore aim to be applicable to these kind of problems; secondly the solutions to the corresponding optimization problems are for these static linear examples unique and analytically attainable. in contrast, we also performed our analysis on two test cases in an ode-framework. parameters in ode models usually have to be estimated, often by solving non-convex optimization problems, and it is then not guaranteed that the optimal solution will be found. in both settings, the number of parameters in the test case models were kept low, in order to more reliably evaluate the methods. in each setting, static and dynamic, we let both models serve as the underlying truth to create artificial data sets each. for each such data set both models served as h and were consecutively fitted to the data, and the goodness of fit was evaluated using various bootstrap approaches, starting with the simplistic combinations described earlier. thus, for each setting and bootstrap approach, this resulted in true positives and true negatives, which were used to construct roc curves to evaluate the power of the tests. combining χ and dw statistics the χ2-test is used to check whether the residuals, i.e. the differences between the model output and the data points, are too big. conversely, the dw test is used the check if these residuals are too correlated. both tests are useful in a model rejection scenario, and in the below analysis, looking at how these tests can be combined, the two tests are used as in their bootstrapped form. simplistic combinations are unsound the first part of the analysis concerns the simplistic combinations: pmin, pmax, pmean, and pprod. although simple, these tests are not without interpretation, and several of them are what at first might seem like the obvious idea. the options min and max corresponds to rejecting if either or if both individual tests reject, respectively. the mean could be thought of as a balancing between the two extremes, and prod, the product, could be thought of as the joint probability. all these four simplistic combinations can be discarded based solely on an analysis of their observed and expected type i error rate. these rates are plotted for all four tests in figure and the interpretations of these curves is as follows. if the plotted lines lie away from the identity line, the expected false positive rate does not coincide with the observed false positive rate, and if this deviation from the identity line is large we call the test unsound. a large deviation means one of two things: either the test is liberal, or the test is conservative. a liberal method is generally regarded as unacceptable, since one wants to be able to trust rejections, but a small level of conservativeness could be accepted, so long as the test is useful. in both the static and the dynamic case, the tested combinations are unsound. the min and prod approaches are strikingly liberal, the max approach is highly conservative, and the mean switches from below to above. these plots should be compared to the single tests: χ and dw which lie along the identity line. this difference between the single tests and the simplistic combinations clearly illustrates that the deviations from the identity line are big. since the results are essentially the same for both the static and dynamic cases, the results were deemed sufficiently convincing to be able to reject all of the tested simplistic approaches as unsound. a two-dimensional approach is both sound and informative the second part of the analysis considers the 2d analysis, combining the χ and the dw tests. although a precursor to the 2d methods presented herein has been mentioned and outlined in previous papers, this is the first time all the implementation details have been solved, and the performance of the method is tried on examples. one of the novel details concerns the density estimation. in contrast to the simplistic combinations, this 2d approach is sound, or only slightly conservative, for both the static and the dynamic case. the conservativeness is tolerated, since the test is informative, as can be gathered from the receiver operator characteristic curves in figure 5a and c. these roc curves are to be interpreted as follows. on the x-axis, the rate of erroneous rejections are plotted; this value is therefore to be as small as possible. on the y-axis, the rate of correct rejections are plotted; this value is therefore to be as high as possible. taken together, this means that the auc should be as big as possible, especially for the region of interest where the fpr is below. this region of interest was chosen because in biology α is rarely higher than. from figure 5a,c it is clear that the new 2d approach outperforms both χ and dw considered as individual tests. however, the combination used here does lie close to the identity line, together with the single test statistics. the upper plots,, show the results for the static example, and for the dynamic example. as can be seen, the 2d χ2 vs dw test is consistent and superior to both of its constituents tests, which is evident from the greater auc in both the static and the dynamic case. introducing a second model replacing the dw test statistic with the χ2-test statistic of a second model the above 2d approach can be generalized to combinations of other tests as well. an important aspect of this is that one can consider a combination of two χ tests, where the second χ test comes from a second model, m. it is therefore intuitively sensible to test whether such a usage of two models is an advantageous usage of this 2d approach. this property of one model’s ability to imitate the behavior of a second model is known as model mimicry, and the idea of using this in a model selection setting has been utilized by e.g.. this second model can in itself be an uninteresting model, i.e. we are then not necessarily interested in the second model as an explanation of the original data, but only in how it interacts with the model being tested. such a model is called a help model. alternatively, the second model could be a competing model and its ability to describe the original data is then of equal importance as that of the first model. if this latter situation is the case, one would typically perform the analysis with both models serving as h generate bootstrap samples from each model, and so on. this version of the 2d test then becomes a form of model comparison, even though there are important differences. for instance, this 2d analysis, unlike model discrimination tests like the conventional non-bootstrapped lhr, can result in all four cases of rejecting none, either, or both of the models. in contrast, a conventional non-bootstrapped lhr can only result in the rejection of the simpler of the two models, or none of them. furthermore, in this new setting, the two models do not have to be nested, i.e. one of the models does not have to be a special case of the other, and the models can be general nonlinear odes. the results on the same dynamic and static examples as before are plotted in figure as can be seen from figure 6b and d, this 2d χ2 vs χ method also has some slight issues with conservativeness, but from figure 6a and c, it is clear that this slight conservativeness should be tolerated: the new 2d version outperforms the previous 2d method. 5: are roc curves, are type i error curves, are for the static, and for the dynamic example. as can be seen, all methods are consistent, or weakly conservative. the roc curves show that the 2d χ2 vs χ test is better than the 2d χ2 vs dw test, and that the bootstrapped lhr, is best of all. in contrast lhp is completely uninformative, as it lies along the identity line. the bootstrapped lhr test is the best approach in the case of a good help model the final test included in this comparison is a bootstrapped version of the lhr. this method has no issues with conservativeness, and outperforms all the other methods in terms of a roc analysis. new motivation for the lhr test in the more general case of bootstrapping non-nested models it is now clear that there are examples where the lhr is advantageous to the 2d χ2 vs χ analysis; let us now understand why this is the case. at a first look, it seems like it should be the other way around: that the 2d χ2 vs χ analysis should be superior, almost by definition. the reason for this is that the information in the lhr is contained in the 2d χ2 vs χ analysis. the close relationship between the two methods can be seen by comparing a cloud from the analysis plotted in the χ2 vs χ plane, with the same cloud plotted in the lhr vs log-likelihood product plane. as can be seen, the shape of the cloud and its relation to the red dot is identical, only tilted degrees. this relation also follows from simple linear algebra. a further inspection of the cloud in figure 7a, which is taken from the static example, reveals that the cloud is of a character that makes a 2d approach superior to the individual tests: just as the cloud in figure the cloud in figure 7a lies away from the axes, and the red dot can be distinguished better in a 2d analysis, than by looking along the individual axes. however, when the cloud has been tilted to the lhr vs lhp plane, the red dot can be separated from the cloud when considering only one of the directions: the lhr direction. that this preservation of the information is preserved when projecting the 2d cloud to the lhr line is also corroborated by the fact that the lhp is, for this example, completely non-informative. in other words, the 1d lhr has extracted virtually all the relevant information of the 2d cloud. all this means that if there would not be any price to pay for doing a 2d analysis, the lhr and 2d χ2 vs χ analysis would be equally good. however, there is a price to pay for moving to a 2d bootstrapping analysis, and this is the reason why the lhr bootstrapping analysis is superior. there are several components to this price. first, the estimation of a distribution is more cumbersome. second, and more importantly, this 2d estimation converges slower than a corresponding 1d density estimation. this has to do with the curse of dimensionality, which simply means that the same number of points quickly become more scarcely packed as the dimensions increase, and that a corresponding density estimation will be based on fewer and fewer neighboring points. this reduction in convergence speed can also be seen in additional file 1: figure s where the lhr has converged already with cloud sizes of ∼ data points, but where the 2d χ2 vs χ analysis requires at least two more orders of magnitude for its convergence. finally, there are also other prices of moving to a 2d analysis, such as the inability to clearly define a one-sided or two-sided test. the existence of such additional factors is also reflected by the fact that the 2d test does not converge to the same roc curve as the lhr test. choosing the second model having established that the inclusion of a second help model may improve the statistical power of tests evaluating the first model, leads to the natural question of whether all help models would do. the answer to this is “no”: if the help model is too simple or too flexible, the advantage is lost, and the resulting model comparison tests - lhr or the 2d χ2 vs χ analysis - perform worse than the other sound tests presented herein. as an example of this, consider the completely flexible model, which simply goes through all data points. if this model, mflex, is used as the help model instead of the suggested competing model in the static example, the 2d cloud collapses to a line: since the cost of the help model is always zero. therefore, there is no 2d advantage to make use of, and the lhr distribution will simply be zero minus the χ distribution, and lhr thus performs as bad as the two-tailed χ test. in the additional file another simple help model is considered: a constant model that simply approximates a data-series with its mean value. here, the picture is a little bit more mixed. for the static example, the model is too simple, and the two-model tests are in-advantageous. for the dynamic test case, the constant model does provide some additional information: the 2d χ2 vs χ analysis performs slightly better, and the lhr test slightly worse, than the single χ2-test. finally, for all of the above examples with too simple or too flexible models, the 2d χ2 vs χ analysis is superior to the lhr test, showing that the lhr is more sensitive to the situation of having chosen a bad help model. application to insulin signaling as a real modeling example, we used data and models from a previous work. in that paper we analyzed experimental data from insulin signaling in primary human adipocytes. some of the experimental data are shown in figure 8b. the data consist of a time series which displays an overshoot: the response increases rapidly from zero, and reaches a maximal value around min, and then decreases to an intermediate steady state value. the observed response is caused by the addition of nm insulin at time zero to previously unstimulated fat cells, and the measurements are performed using sds-page and immunoblotting to determine the degree of auto-phosphorylation of the insulin receptor. the data are normalized such that the first point is zero, and the maximal value is for further details regarding the data, we refer to. using mathematical modeling, we were able to reject several model structures aspiring to explain these data, and we used model selection tools, such as the akaike information criterion, on surviving competing hypotheses. mi,c and the chosen help model mi,b. this example was analyzed previously and is of interest, since the rejection of mi,c would suggest that the recycling is a necessary mechanism to explain the data. depicted in panel is the experimental data z, and fits of mi,b and mi,c. the measured data represent the increased auto-phosphorylation of the insulin receptor in response to nm insulin in isolated primary human adipocytes, as measured by sds-page and immunoblotting. panel shows the bootstrapping cloud in the χ2 vs dw plane, when the bootstraps have been generated by mi,c. as can be seen, the cloud lies along the axes, and there is no benefit of using a 2d analysis. panel shows the χ2 vs χ scatter plot of the mi,c cloud generated after fitting both models to bootstrap sets from mi,c. the corresponding χ2-values for the experimental data is also plotted. as can be seen, the cloud lies away from the axis, and the experimental data point explores the obtained direction. panel summarizes the results. as the clouds have indicated, the χ2 vs dw combination does not improve upon the individual tests, but still lies on the border of rejection. the χ2 vs χ tests on the other hand perform better than the individual tests, and the lhr is best of all. here, we re-examined one of the models, mi,c, from, that could previously not be unambiguously rejected. the model structure of mi,c and the chosen help model, mi,b, are included also in this paper for convenience. mi,b, the smaller model, contains only three reactions: insulin binding with auto-phosphorylation of the receptor, internalization with dephosphorylation, and recycling back to the plasma membrane. mi,b fits to the data with a distinct overshoot, reaching an intermediate quasi-steady state after a couple of minutes. mi,c, on the other hand, is considerably more detailed in terms of the binding, auto-phosphorylation, and internalization, but it does not contain any recycling. mi,c fits to the data in a qualitatively very different way. mi,c has its maximum almost immediately after stimulation, but thereafter never reaches any quasi-steady state as mi,b does. instead the signal declines almost linearly during the observed time period. this example is of biological interest, since the rejection of mi,c would suggest that the recycling is a necessary mechanism to explain the data. this conclusion would follow because mi,c is a more complex, and a more realistic interpretation of the underlying biological system in all other aspects except recycling. in, an aic analysis and a χ2-square test were performed on the various competing models and although mi,c was discarded based on its lack of agreement with data and the fact that inclusion of a recycling reaction yielded a better model, neither of these results were statistically convincing. indeed, although the results pointed towards the rejection of mi,c, it was pointed out already in that a more accurate analysis would require a more generally applicable method such as bootstrapping. in the re-analysis of this insulin signaling example, the analysis was done using all bootstrapped tests that have been found sound on the simpler test cases, where the truth was known. unlike in those examples, here the truth is not known, but subsequent experimental analysis of the system has revealed that mi,c indeed should be rejected. the results are summarized in figure 8e. as can be seen, a bootstrapped χ2-test yielded p =, and a bootstrapped dw test yielded p =, both on the border of rejection. a 2d χ2 vs dw analysis did in this case not improve this value because the resulting cloud and data point did not lie in an informative configuration, p =. conversely, for the 2d χ2 vs χ the cloud and data point did lie in an informative configuration. as a result, the two best methods on the test cases, the 2d χ2 vs χ and the bootstrapped lhr showed improved performance as compared to the above tests, yielding p-values that were clearly below the threshold, p = 8∗10- and p ∼ respectively. also, the internal order of these methods was preserved: lhr was better than the 2d χ2 vs χ these findings are the same as in all tested examples with a sufficiently good help model, and speaks for their generality. discussion in this paper we have extended the traditional parametric bootstrap approach to also look at combinations of different test statistics, here exemplified by the χ2-test statistic and the dw test statistic. we have shown how simplistic combinations, as considered in the literature, of these two statistics are unsound, but that a two-dimensional empirical distribution, as enabled by the bootstrap setting, is both sound and informative. we have further shown that it is even more informative to replace the dw statistic with the χ2-statistic from a second model, resulting in a two-dimensional χ2 vs χ test. however, the one-dimensional bootstrapped lhr is found to be even better, and an analysis of this has led to a new motivation and understanding of the lhr, in the more general case of nonlinear and non-nested models. a 2d bootstrap approach may be superior to a corresponding 1d approach, but there is in fact a trade-off between negative and positive aspects. one positive aspect is of course that a 2d analysis contains more information than a 1d analysis, since e.g. the insights contained in two statistics can be combined. second, even though one in principle can combine test statistics in 1d as well, figure shows that the simplistic combinations give inconsistent results, and therefore are unacceptable. conversely, a 2d-combination of the same test statistics is sound. a third positive aspect of 2d combinations is that they have a higher power than either of the 1d tests including only one of the test statistics. all these positive aspects are due to the benefit illustrated in figure where it is clear that only the 2d combination of the test statistics reveal that the indicated point lies outside the empirical distributions. there are, however, also several negative complications inherent in the addition of an extra dimension, and these complications are solvable to different extents. the first complication concerns the more complicated density estimation that is required in 2d, but this has been resolved in the herein presented 2d approaches. the second complication stems from the relative scarcity of neighboring points in two dimensions and its effect on the density approximation. if the same number of points is used to estimate a 1d distribution and a corresponding 2d distribution, the end result will be a worse approximation for the 2d distribution, simply because of the curse of dimensionality. this second complication can probably be completely overcome by adding more points to the estimation of the empirical distribution, but the addition of more points comes at the price of a higher computational cost. the third complication is that a right-tailed, left-tailed, or two-tailed tolerance interval cannot be clearly defined in two dimensions. this issue is described more below. finally, the positive aspects of the additional dimension only appears if the resulting cloud does not lie along one of the axis, but as in figure in a direction that is not parallel to the axis. all in all, this means that the advantages of a switch to a 2d approach are not guaranteed, while the negative consequences are guaranteed. for this reason the choice between 1d and 2d is a trade-off. similarly to the issue of the second dimension, our results seem to indicate that the addition of a second model provides an advantage, compared to analyzing a single model based only on its own residuals, but there is in fact also here a trade-off between positive and negative aspects. the positive aspects are clear from the examples where we use a help model that is roughly equally good as the tested model: then the roc curves have a higher area under the curve, and because the model that should be rejected in the insulin signaling example gives a lower p-value. these positive sides mean that the additional information provided by a second model provides valuable knowledge regarding the quality of the first model. more specifically, this arguably means that the second dimension from the second models in the examples implies 2d clouds that do not lie along any of the axes, but, as in figure in a direction that is not parallel to the axes. a final positive aspect of this bootstrapped two-model approach is that it goes beyond the capacity of a normal model comparison test, e.g. lhr, since the resulting conclusion can be that none, either, or both models are rejected. a standard model comparison test can never reject both models. the negative aspects of adding a second model are of two types: i) those that have to do with a 2d approach, and which are described above, and ii) those that are associated to the fact that not all additional models provide an improvement. as is clear from e.g. figure 7f, it is clear that a poorly chosen model yields a worse test compared to a mere χ2-test. the second negative aspect is that it is for the non-nested cases not possible to know in advance when a model is good or poor help model. here it should be re-stated that the 2d χ2 vs χ test is more robust towards bad help models than the lhr test in all examples tested herein. in summary, a help model should not be too flexible or too inflexible, and one can see whether the right complexity of the help model has been struck from the shape of the cloud: if it lies along one of the axes it is too flexible or too inflexible. one of the negative aspects mentioned above needs a little further clarification: the definition of the empirical tolerance intervals, which are used to calculate the empirical p-values. first, the 1d versions of the methods that we are using either operate in a one-sided way, or a in a two-sided way. there is no obvious translation of sides and tails in a 1d distribution, to a 2d equivalent. we here adopt the definition of the 2d tolerance region as the region with highest probability density. in practice this is similar to a two-sided cutoff since one may reject a model because it is unrealistically good at describing the data, compared to the given noise level. however, there are differences, such as the possibility to have several fragmented regions instead of a single joint one. therefore, when comparing our method with a 1d-version of the same test, one could consider defining the 1d tolerance interval in a likewise manner, since this more closely mimics the cut-off we do in 2d. however, all comparisons here are done with the one-sided or two-sided 1d-versions of the methods, since it is those that are used in practice, and those that our method should out-compete. a key validation step, and a demonstration of the usefulness of the results herein, is that they are corroborated on a real-world example, which now has been resolved statistically for the first time: we can now reject mi,c with a p-value that lies clearly below the threshold. we have confidence that this rejection of mi,c is correct, because in more recent works we have experimentally shown, by blocking internalization and measuring a downstream component, that recycling of the receptor does play a major role in insulin signaling. we have also measured the amount of internalized insulin receptor and shown that there is far too little internalized ir to explain the observed overshoot. however, even though we in light of these new data have rejected this model, it is interesting to observe that it was possible, with our new method, to reject mi,c based only on the data available at the time of. there are some limitations when interpreting the results that should be mentioned. first, the results are only shown for three example cases, and there is no guarantee that they hold for all other examples. nevertheless, the results are essentially the same for all these three examples: 2d is better than 1d for the static and dynamic examples, and for all examples the tests with a non-extreme help model are better than the single rejection tests, and lhr is best of all. therefore, since the examples include static, dynamic, linear, nonlinear, and real-world aspects, these overall trends probably have some generality. second, the generality is also limited by the fact that we do not have analytical proofs for the results. this, however, is normal for bootstrap approaches. third, another limitation with this approach is that it only considers the usage of a single help model or help statistic. however, this is not a major limitation, since we anyway only advocate the usage of these methods in special cases, where power and accuracy, rather than computational speed, are desired. in other words, our suggestion is to use this approach only in cases where you have a specific model that requires a more detailed analysis. also, it is unlikely that a generalization of this approach to 3d would be beneficial, since then the price of estimating density in a high-dimensional space, and the corresponding slower convergence due to the curse of dimensionality, would be even more severe. it is important to put our results in relation to the existing literature in related fields, such as statistics, systems biology, econometrics, mathematical psychology, phylogenetics etc. first, our method is advantageous only in cases where you have a particularly important and challenging rejection case, where computational time is not a big issue. this stands in contrast to the typical situation of aic and bayesian information criterion, where a big number of models can be sorted through a simple criterion. similarly, the herein presented methods are not suitable to the sub-field of optimal experimental design for the purpose of improved model rejection, since such studies requires an optimization over different experimental designs, which in turn mandates less computationally heavy approaches. second, the perhaps most commonly used method for model rejection, the χ2-test, has a problem - that the degrees of freedom in the χ distribution usually is unknown - but this problem is overcome by using the methods considered herein. however, this is not a new result, but is true for all bootstrap approaches, characterizing the distribution empirically. third, there are a number of commonly used test statistics that we have not considered. for instance, as an alternative to the dw test to measure correlation among the residuals, the whiteness and the run test may be used. it is still unknown whether our results for how to combine test statistics holds also for these other tests. the final, and now following, two comparisons with literature have to do with the lhr and with bayesian approaches. the herein presented analysis presents a new way of showing why and when the lhr is advantageous compared to an individual χ2-test, for the more general setting of nonlinear and non-nested models. since lhr has been both extensively used and studied in the literature, it is important to relate this claim to previous results. lhr was first advocated by neyman and pearson in 1928, and the basic χ distribution relation for nested linear models was known already in the 30’s. these results were generalized to the important case of nonnested models by cox in and to the case of neither of the competing models being true by vuong. however, these results are of limited use, since they rely on analytical derivations of mean and variance terms, and the results by vuong do not even apply to time-series models. also note that there are important cases where the traditional likelihood ratio test is not applicable, e.g. for usage in stochastic models based on the chemical master equation. all of these limitations can be avoided, by adopting a bootstrap approach. this approach basically only relies on the ability to draw bootstrap samples in a way that approximates the true data gathering process. the simulation based type of bootstrap approaches studied herein was originally proposed by williams et al.. the williams approach has all the essential qualities of how we implement the bootstrapped lhr herein: both models are fitted to the data, and the fitted parameters are used to generate bootstrap samples that explicitly incorporates the null hypothesis that the used model is true, and finally both models are fitted to all bootstraps and corresponding distributions of lhr values are calculated. this approach has also been widely used using minor modifications, including an approach where the bootstrap samples are generated using drawing with replacement of the residuals. there are also some papers where theoretical properties of the williams approach have been investigated. for instance, shows that the empirical distribution of williams asymptotically converges to the correct distribution under certain conditions. however, none of those papers use a 2d approach such as ours to achieve an intuitive understanding for why the lhr may be advantageous: that it incorporates the potential added value of the 2d approach compared to the individual χ2-tests, without paying the price of a 2d density estimation. the most important and novel part herein is perhaps that it allows the user to quickly check whether and why the bootstrapped lhr is advantageous or disadvantageous to use compared to the individual χ2-test: it depends on whether the second χ2-test yields a cloud that lies away from being parallel to the axes, which in turn requires that the help model is neither too simple, nor too complex. the final important comparison with literature concerns that with bayesian approaches. bayesian approaches are centered around the combination of a prior distribution or belief with experimental data to obtain a posterior distribution. although bayesian calculations in practice can be done using simple calculations like the bic, the perhaps most common approach involves markov chain monte carlo, and such calculations have big similarities to bootstrap approaches. one important such mcmc-based approach, which has been used in systems biology, is the bayes factor. bf can be viewed as a generalization of the lhr to a bayesian setting. in particular, this means that the likelihoods are integrated over the prior distributions of the parameters, to obtain the ratio of the marginalized distributions. methods to do these marginalizations have been investigated in e.g., and simplified calculations, not requiring the likelihood, using approximate bayesian computations, are considered e.g. in. this inclusion of the parameter uncertainties is important, because in systems biology the parameters are often undetermined, and an important continuation of this work will therefore be to compare the lhr with future extensions of the herein presented frequentist approaches to also include parameter uncertainty. on this note, it should be mentioned that we have done a small analysis to see the effect of the inclusion of such parameter uncertainties on a specific cloud by exploiting the profile likelihood. this small scale analysis indicates that although the results may change upon such an inclusion, the change is not big compared to other uncertainties within the method. another way to take parameter uncertainty into account is by introducing an initial step of non-parametric bootstrapping into the bootstrap sample generation, as done e.g. in. once parameter uncertainty is taken into account in the bootstrap setting in this way, there is a striking similarity to the bayesian posterior predictive checks. in bpp, the predicted distribution of future experiments is compared to the observed data. this is done by generating new data by sampling and simulating from the posterior distribution, and then comparing the resulting distribution of goodness-of-fit with the gof from the experimental data. with all these similarities pointed out, it should also be recalled that bayesian approaches are never identical to frequentist approaches, since frequentist approaches do not require a prior.  in a bootstrap setting, it is possible to obtain joint distributions for combinations of test statistics in a more straightforward way than is possible in an analytical setting, but this possibility has previously been little explored. we here show that such combinations often do provide additional knowledge not contained in the individual tests, but that the considered simplistic combinations, like max and min, yield inconsistent, i.e. overly conservative or liberal, results. a new 2d approach, on the other hand, is only mildly conservative, and is superior to the individual tests. these results were obtained on both a static and dynamic case, where the truth is known. on the same examples, a 2d χ2 vs χ test is superior to a 2d χ2 vs dw test, where the additional χ2-value comes from the ability of a second model to describe bootstrap samples from the tested model. the 2d χ2 vs χ test is, in turn, outperformed by the 1d bootstrapped lhr. these results are also confirmed on a previously published rejection example from insulin signaling in human fat cells, which has now been statistically resolved for the first time. further analysis of these results show that whether or not a 2d combination is advantageous depends on a balancing between positive and negative aspects. the positive sides are found if the cloud as in figure lies in a direction non-parallel to either of the axes, and the price to exploit this is e.g. that density estimation in 2d converges more slowly, and that one cannot define one-sided or two-sided cutoffs for the tolerance regions. similarly, the additional model only provides a benefit if it is of a rightly balanced ability to describe the data; otherwise using the additional model worsens the performance. it is because of these balancing acts between positive and negative aspects that lhr may be the better choice: if the additional model is of appropriate complexity, lhr extracts all the useful information of the χ2 vs χ plot with a one-dimensional analysis, which thus avoids the estimation of a 2d density. this analysis thus provides a new motivation for the lhr test, which is valid in the case of non-linear and non-nested models. in summary, these results provide useful insights into the important systems biology problem of model rejection: when to use, and when not to use, 2d approaches and additional models. these methods are to be exploited in challenging and important cases, when accuracy and power rather than computational speed are prioritized. competing interests the authors declare that they have no competing interests. authors’ contributions rj developed, implemented, and analyzed the methods. rj also constructed the test cases, and drafted the manuscript. ps contributed with biological models, experimental data, and interpretations. gc was the main designer of the study and co-drafted the manuscript. all authors read and approved the final manuscript. supplementary material additional file 1 supplementary material. the pdf-file “rj et al. - supplementary material.pdf” contains additional elaboration on methods used, some additional results and analysis, and discussion on some of the issues covered in this manuscript, as well as all supplementary figures referred to. click here for file acknowledgements this work was supported by the european commission through the biosim project, by the lions sweden research funding, by the swedish research council, and by the german federal ministry of education and research virtual liver grant 0313074d. the authors also wish to thank clemens kreutz for crucial input on statistical analysis techniques, input regarding the development of the method, and for comments on a previous version of the manuscript. the authors also want to thank jens timmer for input on the general study design and for financial support.
a new two-dimensional approach based on parametric bootstrapping, on the other hand, is found both consistent and with a higher power than the individual tests, when tested on static and dynamic examples where the truth is known. in the same examples, the most superior test is a 2d χ2vsχ where the second χ2-value comes from an additional help model, and its ability to describe bootstraps from the tested model. this superiority is lost if the help model is too simple, or too flexible. if a useful help model is found, the most powerful approach is the bootstrapped log-likelihood ratio. this problem is partly alleviated by the usage of bootstrapping, a computationally heavy approach to calculate an empirical distribution. bootstrapping also allows for a natural extension to estimation of joint distributions, but this feature has so far been little exploited. we have shown how to, and how not to, combine tests in a bootstrap setting, when the combination is advantageous, and when it is advantageous to include a second model. these results also provide a deeper insight into the original motivation for formulating the lhr, for the more general setting of nonlinear and non-nested models. analytical formulas for the corresponding distributions rely on assumptions that typically are not fulfilled. model rejections lie at the heart of systems biology, since they provide conclusive statements: that the corresponding mechanistic assumptions do not serve as valid explanations for the experimental data. we herein show that simplistic combinations of bootstrapped tests, like the max or min of the individual p-values, give inconsistent, i.e. these insights are valuable in cases when accuracy and power, rather than computational speed, are prioritized. the chi-square test or the durbin-watson test. these approaches statistically resolve a previously published rejection example for the first time.
 model rejections lie at the heart of systems biology, since they provide conclusive statements: that the corresponding mechanistic assumptions do not serve as valid explanations for the experimental data. rejections are usually done using e.g. the chi-square test or the durbin-watson test. analytical formulas for the corresponding distributions rely on assumptions that typically are not fulfilled. this problem is partly alleviated by the usage of bootstrapping, a computationally heavy approach to calculate an empirical distribution. bootstrapping also allows for a natural extension to estimation of joint distributions, but this feature has so far been little exploited.  we herein show that simplistic combinations of bootstrapped tests, like the max or min of the individual p-values, give inconsistent, i.e. overly conservative or liberal, results. a new two-dimensional approach based on parametric bootstrapping, on the other hand, is found both consistent and with a higher power than the individual tests, when tested on static and dynamic examples where the truth is known. in the same examples, the most superior test is a 2d χ2vsχ where the second χ2-value comes from an additional help model, and its ability to describe bootstraps from the tested model. this superiority is lost if the help model is too simple, or too flexible. if a useful help model is found, the most powerful approach is the bootstrapped log-likelihood ratio. we show that this is because the lhr is one-dimensional, because the second dimension comes at a cost, and because lhr has retained most of the crucial information in the 2d distribution. these approaches statistically resolve a previously published rejection example for the first time.  we have shown how to, and how not to, combine tests in a bootstrap setting, when the combination is advantageous, and when it is advantageous to include a second model. these results also provide a deeper insight into the original motivation for formulating the lhr, for the more general setting of nonlinear and non-nested models. these insights are valuable in cases when accuracy and power, rather than computational speed, are prioritized. model rejectionbootstrappingcombining information2dinsulin signalingmodel mimicrylikelihood ratio 