Validity of regression diagnostics for deterministic computer experiments

Question

The question

Is it possible to justify the use of regression diagnostics such as F-tests or the AIC for the analysis of deterministic computer experiments?

Background

One way to analyse computer experiments is by building surrogate models. The unknown function $f$ is approximated by a nice and explicitly known function $\hat f$, determined from some training data $(x_i)$ and the according responses $(f(x_i))$.

To make the question more concrete it suffices to focus on multivariate polynomials as surrogate functions. They are fit using standard regression tools (e.g. "lm" in R or "sklearn.linear_model.LinearRegression" in Python). Furthermore, since design matrices can get large very quick, methods relying on variable selection such as stepwise regression are often used.

Notice that in most cases neither of the two main regression assumptions are fulfilled: More often than not the unknown function $f$ will be ill-specified, i.e. not be a polynomial of any degree, and there will be no noisy observations, i.e. independent errors with a normal distribution.

My issue

The minimisation of residual sum of squares makes perfect sense for an approximation task. Hence, I do understand the use of regression tools in general and diagnostics such as $R^2$ even though the regression assumptions are violated.

But it seems to me that diagnostics such as F-tests and the AIC do rely crucially on the likelihood or other statistical assumptions. While tools calculate them readily (but quite formally it seems), I do not see any reason why stepwise regression or variable selection based on such diagnostics should work in the deterministic and ill-specified case.

That said, I have seen this approach used in practice and it is also mentioned in textbooks (e.g. Chapter 7.2.2 of [1]). So maybe there is something to it?

In other words: Can stepwise regression or variable selection using AIC (or similar statistics based on the likelihood or distributional assumptions of the regression model) be justified also in the case of deterministic computer experiments?

[1] "The Design and Analysis of Computer Experiments" by Thomas J. Santner, Brian J. Williams, William I. Notz

Answer 1

This question reminds me a bit of the question Could a mismatch between loss functions used for fitting vs. tuning parameter selection be justified? , where I answered that it can be useful to use a different loss function than the loss function which needs to be optimized. The reason is that statistical fluctuations might be best minimized by using a loss function that relates to the statistical distribution of the errors. For example: when we are estimating the location parameter of a Laplace distribution, then finding an estimator that minimizes the sum of absolute residuals will be better in optimizing the expected square error than an estimator that minimizes the sum of squared residuals.

In your example there is also in some way statistical variation. You have a deterministic function (the population), but you only take a sample $x_i$ where you evaluate $f(x_i)$, so you have a sub-selection of the population*, and that resembles the randomness of sampling.

I do not see any reason why stepwise regression or variable selection based on such diagnostics should work in the deterministic and ill-specified case.

The error distribution of this 'randomness' is not very clear, and to optimize AIC (based on an arbitrary likelihood) or F statistic will only make sense as a heuristic diagnostic. Yet they do work and that is because they introduce some way to penalize the number of terms. You can fit a high order polynomial and reach an $R^2 = 1$ but that is possibly introducing a larger prediction error.

Instead, to reduce this error, it might be better to use, for instance, leave one out cross validation, which relates to the mean squared error (Proof of LOOCV formula), and is asymptotically similar to optimizing a likelihood (Equivalence of AIC and LOOCV under mismatched loss functions)

---

^{*(If you would know the entire function $f$ then you could just optimize the entire function straight away, and if you have no limits to the number of terms then you get something like a perfectly fitting Taylor series or Fourier series).}

Answer 2

Statistical methods applied to deterministic phenomena: Statistical methods like regression analysis can be quite useful in situations where one is dealing with deterministic values, particularly in cases where the deterministic behaviour is complex enough that it is not helpful to express this behaviour in deterministic terms. In such cases, statistical methods can be employed, where we model deterministic data using a statistical model that includes an "error term" with a particular distribution. As with other applications, you can still use diagnostic tests to check if the (deterministic) residuals from the model appear to follow the assumed distribution, etc. Once it is determined that a statistical model is applicable to a set of data (deterministic or not) and diagnostic tests confirm plausibility of the model assumptions, the entire suite of tests and methods applicable to that model are okay to implement. In a regression context, this includes methods like F-tests, AIC computations, etc.

In regard to this issue, it is worth noting that simulations used in statistical practice are generated by pseudo-random number generators (PRNGs) that are deterministic in nature. (So long as you set the seed for these generators in a deterministic manner, the resulting series of numbers is deterministic.) Indeed, the entire field of constructing valid PRNGs and the entire field of simulation analysis using their values involves the application of statistical models and tests to deterministic phenomena.

An applied example: As an example of a regression analysis of this kind, you can have a look at O'Neill (2020) (pp. 6-9, 11-12). That paper examines the classical occupancy distribution and computes the accuracy of the normal approximation to that distribution over a large range of parameter values. Part of the analysis involves showing that the approximation tends to get more accurate as the parameters get large, and this is done by computing some summary quantities relating to the approximation error, and showing how they change as the parameter values increase. To this end, the paper includes a regression analysis looking at how these summary quantities for the approximation error reduce as the parameter values become large, using a model that has extremely high goodness-of-fit. It turns out that these quantities follow very closely (but not exactly) to a simple deterministic function and so the regression analysis is able to show that this simple relationship has a high goodness-of-fit. What is notable in this analysis is that the data used for the regression is purely deterministive --- it consists of computed RMSE values comparing two known distributions over a set of known parameter values. There is no randomness or real-world data in that analysis.

Philosophical considerations: When conducting this kind of analysis it is worth noting that it is possible to take an underlying philosophical position with respect to probability that does not assume the existence of (aleatory) randomness and so does not contradict determinism. Epistemic interpretations of probability generally consider probability distributions to describe our own uncertainty in values, even if those values are fixed. Statisticians who take an underlying epistemic interpretation of probability (which would include most Bayesians at a minimum) usually have no in-principle aversion to applying statistical methods to deterministic phenomena.

Answer 3

It sounds like your deterministic experiments are effectively evaluating an unknown function $f$ at many different input values $x$. If evaluating $f$ at new $x$ values is cheap, then you can easily get arbitrarily large samples; so hypothesis tests (which are partly a measure of sample size) aren't going to be a very useful way to choose a parsimonious $\hat f$. Similar arguments can be made for AIC.

Instead, I'd suggest you think of a measure of "How good does my approximation need to be to be useful for my purposes?" What are the units of $f$? Choose the smallest value of a maximum or average discrepancy between $f$ and $\hat f$ that would be acceptable for your needs, on the range of $x$ values that you care about. Then, evaluate $f$ at many $x$ values, so that you have more than enough data to get low variance in the coefficient estimates when fitting fairly complex $\hat f$ models. Choose the simplest of these $\hat f$ models that meets your criteria for a "good enough" approximation.

---

On the other hand, if evaluating $f$ for new $x$ values is expensive or slow and your sample size is essentially fixed, then the F test or AIC might help you choose the most-complex $\hat f$ you can afford to fit with this sample size. The fact that $f$ is deterministic isn't a problem (though it might affect your choice of which test to carry out; see below).

It often helps to think of significance tests (like the F test) as a measure of sample size. I will assume you're using the F test that compares two nested regression models. Then the F test is basically asking: "Is my sample large enough to trust the larger / more complex model? Or is there too much variability relative to this sample size, and it would be safer to use the smaller / less complex model at this sample size?"

Say the true unknown $f$ is pretty close to a quadratic function (but not quite exactly quadratic). And say you first try to compare a linear $\hat {f_1}$ vs a quadratic $\hat {f_2}$ fit to the pairs $(x, f(x))$. If you've evaluated $f$ at very few values of $x$, the F test might fail to reject $H_0$, telling you that you don't have enough data and your estimated $\hat {f_2}$ is too noisy. But if you've evaluated $f$ at enough different $x$ values, the F test might reject $H_0$, telling you that you've got enough data to safely trust your estimated $\hat {f_2}$ as a better approximation than $\hat {f_1}$.

But this issue of "Is my sample big enough?" is distinct from "Is $\hat f$ a good enough approximation to $f$?" The latter question requires subject matter expertise. Again, imagine that $f$ is nearly-but-not-quite quadratic. And let's say that the true differences between $f$ and its best quadratic approximation are negligible for your purposes; a cubic $\hat{f_3}$ could technically be an even better fit, but it would make almost no practical difference in whatever you're using $\hat f$ for. Even so, if you have a large enough sample size, the F test comparing quadratic vs cubic will reject $H_0$ and say you've got enough data to trust the cubic approximation $\hat{f_3}$.

My key point is that "negligible for your purposes" is something that you must decide and the F test cannot answer for you. I don't know the context of your computer experiments, but here's a physics example. Say you're trying to predict $f$ = boiling point of water from $x$ = atmospheric pressure. The true relationship is not linear, but if your purpose is to suggest approximate baking times in a recipe book, a linear $\hat f$ might be plenty good enough. On the other hand if your purpose is to design high-performance meteorological equipment, a linear $\hat f$ may be nowhere near good enough and you'll need a more complex model (and enough data to fit it).

---

Finally, if you do use hypothesis tests in this scenario, and your deterministic $f$ is smooth and your $x$ values are very close together... then the independent-residuals assumption isn't going to be appropriate.

Think of time-series data. Imagine an extreme example: $y$ is changing once a day (and is independent from day to day), but I record measurements every hour. Then I don't really have $n$ independent measurements -- my "effective sample size" is more like $n/24$. If I fit a regression model to $y$ vs $time$ and don't account for this autocorrelation, my F tests will wrongly think that I do have $n$ independent measurements, and my p-values will be way smaller than they really ought to be for the effective sample size.

In your case it might not be quite as extreme. But I'd still look into autocorrelation models to account for the fact that similar $x$ values might have extremely similar $f$ values.

Answer 4

In other words: Can stepwise regression or variable selection using AIC (or similar statistics based on the likelihood or distributional assumptions of the regression model) be justified also in the case of deterministic computer experiments?

Here's an argument for likelihood-based statistics. Let $Y = f(X)$, where $f$ is a deterministic function. Assume you observe a finite number of realizations s.t. $\mathbf{y} = f(\mathbf{X})$. If you knew both $f$ and $\mathbf{X}$, then you would know $\mathbf{y}$ with probability one, or very vaguely $P(\mathbf{y} = f(\mathbf{X}) | \mathbf{X}, f) = 1$. You don't know $f$, so you place a Gaussian process prior $f | \mathbf{X}\sim\mathcal{GP}(\cdot, \cdot)$. You integrate out $f$ to obtain the output marginal likelihood, which is MVN Gaussian: $\mathbf{y} | \mathbf{X}\sim\mathcal{MVN}(\cdot, \cdot)$. Even if the computer code output is deterministic and there's no "innate" variability on it, the output marginal likelihood is multivariate normal due to your uncertainty on $f$. So, statistics based on this marginal likelihood seem sensible to me.

More general, here are a few things that might help you decide how comfortable you feel about using such diagnostic techniques.

• Your computer model is deterministic, i.e., $Y = f(X)$. Since the computer model does not have a random error component, there's no aleatoric uncertainty. However, you don't know $f$, so there's epistemic uncertainty. See Aleatoric and epistemic.
• Even though the underlying function $f$ might be deterministic, the output value is unknown a priori (before running the simulator). This is sometimes called code uncertainty and is enough to justify the use of a stochastic process.
• Surrogates are most frequently set up in a Bayesian framework, and so the probabilistic statements are about your uncertainty about the function that you are approximating $f$ and not really about the inherent random noise in your computer code. This applies even when Empirical Bayes is used, which in many contexts can get quite close to simply "find MLE estimates but interpret like a Bayesian model".
• Bastos and O'Hagan (2009, sec. 3.1) discusses the challenges associated with diagnostics for linear models fitted to a deterministic function. In particular, the training set predictions will be perfect if you use an interpolator, thus the need to use cross-validation or a new data set. Out-of-sample validation introduces some uncertainty. However, they also note that there might be scoring functions for deterministic computer experiments such as the Mahalanobis distance. See 10.1198/TECH.2009.08019

An adjacent note.

• A linear model might not be commonly used for a surrogate. However, you can cast a linear model as a Gaussian process with linear kernel $\mathbf{XX}^\top$ (Rasmussen and Williams, 2006, sec. 4.2.2).

References

• Bastos, L. S., & O’Hagan, A. (2009). Diagnostics for Gaussian process emulators. Technometrics, 51(4), 425–438. https://doi.org/10/bw62bq
• Rasmussen, C. E., & Williams, C. K. I. (2006). Gaussian processes for machine learning. MIT Press.