# Design choice and machine learning model performances

Rosa Arborettia, Riccardo Ceccatob, Luca Pegorarob and Luigi Salmasob,*
aDepartment of Civil, Environmental and Architectural Engineering, Università
degli Studi di Padova, Padua, Italy; bDepartment of Management and
Engineering, Università degli Studi di Padova, Vicenza, Italy *corresponding
author<EMAIL_ADDRESS>
###### Abstract

An increasing number of publications present the joint application of Design
of Experiments (DOE) and machine learning (ML) as a methodology to collect and
analyze data on a specific industrial phenomenon. However, the literature
shows that the choice of the design for data collection and model for data
analysis is often not driven by statistical or algorithmic advantages, thus
there is a lack of studies which provide guidelines on what designs and ML
models to jointly use for data collection and analysis. This article discusses
the choice of design in relation to the ML model performances. A study is
conducted that considers 12 experimental designs, 7 families of predictive
models, 7 test functions that emulate physical processes, and 8 noise
settings, both homoscedastic and heteroscedastic. The results of the research
can have an immediate impact on the work of practitioners, providing
guidelines for practical applications of DOE and ML.

###### keywords:

Simulation study , Physical experiments , Predictive analytics , Gaussian
Process , Artificial Neural Networks

##  1 Introduction

In many industries the development of mathematical models of physical
phenomena or processes is of primary interest. Objectives of such models may
include prediction of one or more quantities or identification of the optimal
configuration of parameters which leads to an optimization of some quantities
of interest, such as a quality criterion or a cost measure. Design of
Experiments (DOE) is a well-known statistical procedure that has been widely
applied for decades especially for optimization, but also for the construction
of surrogate models able to reliably predict one or more responses of
interest. Historically, the most widespread model employed for the analysis of
DOE data, which also stands at the core of the response surface methodology
(RSM) 1, is the second-order model 23.

Recently, with the advent of the fourth industrial revolution and the
consequential increase in data availability in the industrial environment,
machine learning (ML) techniques for data analysis have rapidly spread in
industry. In this context, a recent trend consists in the adoption of ML
algorithms for the analysis of experimental data, in lieu of the typical
second-order model, as proved by the increasing number of publications dealing
with these topics, as shown in Arboretti et al. 2. However, it has also been
shown that there is a lack of studies that try to provide some guidelines on
what experimental designs appear to be more appropriate when ML models are
selected as a tool for analysis 2. Instead, in the vast majority of papers
that jointly apply DOE and ML, the choice of one specific design appears to be
driven by incidental factors such as the analyst’s experience or diffusion in
one specific area 2. The aim of this article is to show in an empirical manner
what designs are more promising for collection of data that is then modeled
through ML techniques. As a collateral result, this will also lead to
indications being provided on what predictive analytics techniques to take
into consideration when analyzing DOE data. It should be pointed out here that
only predictive performance of the algorithms in terms of prediction error
will be taken into account, and other aspects, such as the ability to quantify
uncertainty of predictions and the possibility of investigating causality, are
features that should be considered when choosing an algorithm for prediction
or optimization 2. The authors encourage further research devoted to these
topics, but this is outside the scope of the present work.

The paper is organized as follows: section 2 presents the literature
background to the application of ML for the analysis of DOE data; section 3
presents the methodology and details of the simulation study conducted;
section 4 shows the results and section 5 discusses the findings in depth.
Finally, in section 6, a summary of the main contributions of the work and the
conclusions are presented.

##  2 Literature background

In their recent article, Arboretti et al. 2 provided a comprehensive overview
of the topics of DOE and ML applied in a product innovation (PI) setting.
Arboretti et al. 2 highlight both the advantages and challenges of the
application of ML in PI, discuss the implications of a joint adoption of DOE
and ML, and identify the most common experimental designs employed and ML
algorithms chosen for analysis. From their work, some clear advantages of the
adoption of ML for data analysis appear evident, including the ability of ML
to appropriately model data without requiring assumptions on the underlying
distribution, the ability to model complex non-linear relationships, and the
general capability of providing a better fit to the data, thus ensuring more
accurate predictions and the requirement for a less strict design structure,
to the point that undesigned data can be effectively modeled. Additionally, a
proper adoption of a DOE+ML framework can enhance further advantages,
including the minimization of the number of trials through the identification
of the most informative combinations of the factor levels, the ability to
explore experimental regions omitted by DOE alone, the possibility to automate
the experimentation procedure, and others 2.

An in-depth analysis of the type of DOE used in combination with ML modeling
was also conducted 2, revealing high fragmentation in the choice of designs.
The most recurring class of DOEs comes from the RSM literature including
Central Composite Designs 3 (CCDs), Box-Behnken Designs 4 (BBDs), and full and
fractional factorial designs (FFDs). Another relevant group of experimental
designs comes from the Robust Parameter Design (RPD) literature 5, 6, 7, while
space-filling designs are seldomly applied 8. Overall, the non-sequential
experimentation procedure is still predominant 2, despite the fact it has been
shown that several advantages are embedded in a sequential strategy 9, 10, 11.
The overall picture in regards to the choice of design in a DOE+ML framework
is not that surprising as another finding is that physical experiments are by
far more common than computer experiments 2, meaning that traditional designs
from RSM are favored. The main reason is that such designs have been
specifically developed for the optimization of physical systems and they can
count on a very well established literature set that proposes solutions for
many of the typical problems which characterize the industrial context, e.g.
restrictions to randomization with the split-plot designs 12, or the need for
an estimation of many effects while minimizing the experimental effort with
definitive screening designs 13. Other designs, such as space-filling designs,
have been developed for computer experiments in which the concepts of
randomization and replication, that are at the base of traditional DOE, no
longer count because experiments are typically deterministic and the frequent
change in factor levels is possible and does not imply a cost. Also, the
preference for non-sequential experimentation is likely a consequence of
physical experimentation, since companies may adopt strict procedures and
protocols that discourage the adoption of a sequential approach. For instance,
access to a system under investigation may be limited because experimentation
can cause a delay to production deadlines and management may be unwilling to
reserve the required resources for a prolonged period. Thus, in many
situations, execution of an entire DOE in one single batch of trials is
preferable from a production planning perspective.

Arboretti et al. 2 also provide a picture of the ML models most widely applied
to DOE data for product innovation, revealing extensive use of Artificial
Neural Networks (ANNs), and sporadic adoption of other prediction algorithms,
such as Gaussian Processes (GPs) and Support Vector Regression models (SVMs).

##  3 Design and model choice

In this section we detail the procedure and simulation study conducted to
assess the performances of different designs paired with different ML models
for data analysis. In section 3.1 we present the research objective and its
boundaries, in section 3.2 and 3.3 we briefly review the chosen designs and
the ML algorithms respectively. In sections 3.4 and 3.5 we describe the
selected test functions and different noise structures and levels, and in
section 3.6 we detail the ranking procedure adopted to obtain the final rank
of designs and algorithms.

###  3.1 Objective and boundaries

The objective of the simulation study is to provide, for the first time, an
empirical assessment of the suitability of different experimental designs for
use with ML models for data analysis and prediction. Consequently, the most
promising predictive models will also be indicated. This research is subject
to some constraints, that come mainly from the literature background as
presented in section 2:

  * 1.

The focus is on the designs, although several ML algorithms will be tested.
Some indications will also therefore be provided for the choice of ML model,
but solely predictive performance will be considered while crucial aspects
such as uncertainty quantification or model interpretability will not be
examined.

  * 2.

Only the non-sequential DOE case is investigated, despite the sequential
design approach being very promising. This is mainly for two reasons: (i) the
literature 2 shows that the non-sequential approach is still predominant; (ii)
in the sequential experimentation setting there is also the need for one
initial dataset to fit the algorithms in the first iteration. Since undesigned
data are often presently employed as the initial dataset 2, sequential
experimentation will also benefit from the results of this study.

  * 3.

The focus is on physical experiments. Consequently:

    * (a)

Non-linear test functions that emulate physical processes will mainly be
selected.

    * (b)

Small data will be considered.

    * (c)

Designs with a moderate number of factor levels will be favored.

    * (d)

Both the homoscedastic and heteroscedastic noise cases will be taken into
consideration.

To the best of our knowledge, this is the first study devoted to the analysis
of the impact of choice of experimental design when the method chosen for data
analysis is of the ML type. Similar studies have been published that focus on
the methods employed for data analysis, but they refer to the case of computer
experiments, thus select only space-filling designs for data collection and
choose GP models for prediction 14. In this study we take a larger number of
designs with heterogeneous characteristics, and also test a large number of
different models. A representation of the general structure of the simulation
study is reported in Figure 1, and further details are provided in the
following sections.

Figure 1: General framework of the simulation study.

###  3.2 Experimental designs

Table LABEL:DOEs_SS shows a summary of the experimental designs selected for
this study. The designs will be used to collect data on 7 different test
functions, all having 6 active dimensions (section 3.4). This is in line with
most industrial applications of DOE, as it has been shown that 75% of case
studies in engineering consider 6 or less factors 15.

Since the principal aim of this work is to provide guidelines to practitioners
regarding what experimental designs to choose in physical experimentation when
data is analyzed with predictive models, we decided to allocate for each
alternative design the same budget in terms of number of trials, to provide a
fair comparison between the competitive options. Often this is also the case
for practical applications of DOE in industry, where a given amount of
resources is allocated for a specific study and the issue is to choose one
experimental design to appropriately pursue the goals of the project. For this
reason, some powerful designs such as definitive screening designs 13 are not
considered in this study, as they would require a different number of trials
given the chosen number of experimental factors, making the comparison with
other designs unfair.

Furthermore, and as already pointed out, we suppose that the main objective of
the project is to obtain accurate predictions about a specific physical
phenomenon. Even though we recognize that other aspects (e.g. uncertainty
quantification, interpretability) are also important, we argue that if the
analyst selects an intricate predictive model (e.g. a ML algorithm) to analyze
a non-linear phenomenon, the main objective is to produce accurate
predictions. Some evidence from the literature 2 also confirms this, given
that in product innovation ML models are employed for the analysis of DOE data
mainly for their advantages in terms of predictive performance and ability to
model complex non-linear relationships.

The designs considered in this work have been grouped into three categories:
(i) “classical designs”, (ii) “optimal designs” and (iii) “space-filling
designs”. In the “classical designs” category we include the CCD, BBD and FFD.
Even though such designs have been around for decades, it has been shown that
they are still the most used for data collection in ML studies on product
innovation 2. The “optimal designs” included are the D-optimal and I-optimal
that have recently replaced classical designs in many applications. The
designs from the “space-filling” category have been included to provide a
benchmark for those practitioners that may be familiar with the literature on
computer experiments, and may therefore want to assess the potential of the
predictive models to provide accurate predictions when a maximum number of
levels is included in the design. This is because a large number of levels
provides more supporting points to the predictive models, and may thus favour
a better understanding of the non-linearities in the data. However, we
recognize that such designs are not practical in physical experimentation, and
that they are usually employed to analyze the behavior of deterministic
computer solvers, while here the main focus is on noisy physical experiments.
This is also why only a few designs in this category are selected, and other
widely used designs (e.g. maximin LHDs 16), are not considered here.
Nevertheless, one “hybrid” design from the space-filling literature with
characteristics that enable its application to physical experiments is also
considered.

In the following sections we provide details of the designs in each category.

####  3.2.1 Classical designs

Three classical designs have been selected from the RSM literature, namely
CCD, BBD and FFD. The CCD was built from a two-level factorial design of
resolution V with 4 center points added; an appropriate alpha parameter was
then chosen to obtain a rotatable design and axial (star) points were included
accordingly. The whole design was scaled to have the range determined by the
axial points equal to the range of the levels for each test function, as
detailed in section 3.4. In total, 52 runs were required to generate the
entire design, and the same budget in terms of runs was selected for all the
other designs. Similarly, the BBD includes 4 central points and has 52 runs.
Both the CCD and BBD have 3 levels for each factor 3, 4. The FFD was built as
a full factorial design including 6 levels for each factor. To fit the budget
of 52 runs, the design was then computer-optimized using the D-optimality
criterion. The rationale behind the choice of number of factor levels for the
FFD was that we expected a larger number of levels to provide more granularity
to the predictive models, thus enabling a more precise estimate of the
underlying non-linear function. Furthermore, we argue that this number of
levels is not prohibitive for application in physical experiments, and can
represent an appropriate trade-off between a better (expected) predictive
performance and a marginally more complicated execution of the experiment.

####  3.2.2 Optimal designs

In this paper we consider designs that satisfy the D- and I-optimality
criteria. The D-optimality criterion is the most commonly used, and seeks
designs that maximize the determinant of the Fisher information matrix, thus
ensuring a minimization of the volume of the confidence ellipsoid around the
model parameters. The main objective is to obtain a precise estimate of the
effects, therefore D-optimal designs can be considered estimation-oriented. If
the goal of the study is, as in this case, prediction-oriented, it may be
preferable to opt for I-optimal designs, since they minimize the average
prediction variance over the experimental region. Nowadays, this is the
default criterion in some commercial software 17 when the main goal of the DOE
study is to obtain a reliable prediction of the response (over all the design
space or in the vicinity of the optimal configuration), such as in RSM
applications. Furthermore, the evaluation of the prediction variance over the
design space is a popular quality criterion for the choice of a DOE 18.

Another well-known prediction-oriented criterion is G-optimality that ensures
a minimization of the maximum variance of prediction. However, it has been
shown that a minimization of the maximum variance of prediction implies a
greater prediction variance over a large part of the region of interest 19. In
this study the objective is to minimize the prediction error over the entire
design space, thus we discard the G-optimality criterion and prefer the
I-optimal case 20. Other estimation-oriented optimality criteria also exist
(e.g. A-optimality), but given the different objective of this work, we only
focus on the D-optimal case as it is the most prominent in the category. For a
more detailed discussion of the main optimality criteria, we refer readers to
the paper by Jones et al. 21, pages 370-372.

In the simulation study the optimal designs (D_opt, I_opt) were built levering
the capabilities of JMP®’s Custom Design platform - this employs the
coordinate exchange algorithm 22 over multiple initializations to generate the
preferable design given an optimality criterion. The analyst is required to
input the levels for each factor (6 equally-spaced levels), the number of runs
(52) and a list of effects that it is necessary or desirable to estimate 17.

Finally, the impact of replication is also assessed for the optimal designs.
To this end, D- and I-optimal designs were constructed with JMP®’s Custom
Design platform having a 50% level of replication (26 base runs replicated
twice, D_opt_50%repl and I_opt_50%repl respectively). A full replication of
the base design is often performed in DOE studies and it is one of the funding
principles of DOE for physical experimentation 23. Furthermore, on several
occasions it has also been shown to provide advantages for stochastic computer
experiments when predictive models (especially GPs) are used 24, 25, 26, 27.

####  3.2.3 Space-filling designs

Two designs have also been selected from the space-filling family, namely a
Random Latin Hypercube Design (LHD_rand) 28, 29 and a Maximum Projection
Design (MAXPRO) 30. The LHD_rand design is one of the first and most prominent
designs of the space-filling family 28 and is constructed to spread out the
experimental points in a way that ensures each of the factors has all
proportions of its distribution represented by input values 28, with the aim
of collecting as much information as possible on the design space. The MAXPRO
design is a more recent evolution of the LHD and has been shown to have
several advantages over the other space-filling designs, the main one being
optimization of the projection properties for all subspaces of the factors 30,
8. This is a fundamental concept for both computer and physical experiments
when only some of the factors are active, which is often the case due to the
sparsity-of-effects principle 31, 23. A consequence of the maximum projection
property 30 is that if some of the factors prove to be inactive during the
experiment, they may be dropped with no consequences to the space-filling
structure of the design, which is also advisable in a physical experiment.
However, both the LHD_rand and MAXPRO designs have a characteristic that
limits their applicability to physical experiments: they have as many factor
levels as there are runs, since they tend to maximize the spread of the points
in the design space. In practical applications in physical experimentation, a
large number of levels for each factor may be neither feasible nor advisable,
since the presence of noise may obscure some of the active effects. A recent
generalization of the Maximum Projection criterion overcomes this problem by
enabling MAXPRO designs to deal with multiple types of factors.

In their recent article, Joseph et al. 32 updated the MAXPRO designs allowing
both categorical and discrete numeric factors to be treated. In their work the
authors refer to the case of computer experiments in which some of the
variables of a computer simulator may be categorical (e.g. material type) or
discrete numeric (the continuous factor can only take some pre-specified
levels, e.g. number of flutes in a solid end milling process) 32. In this
article we mainly leverage one of the advantages of the updated MAXPRO
designs: if a limited number of levels is identified, the continuous factors
can be restricted to the discrete numeric case, meaning the MAXPRO designs are
also applicable in physical experimentation 33. In our view, and for the scope
of the present paper, this makes MAXPRO designs with discrete numeric factors
(MAXPRO_dis) hybrid designs: while they are constructed to be essentially
space-filling (although only on a limited number of levels), the limited
number of levels makes them competitive with RSM designs for physical
experimentation 33. To our knowledge, this is the only class of designs of the
space-filling type that can actually be applied to physical experimentation 33
by taking a limited number of levels for each continuous factor. The
MAXPRO_dis design selected for this study counts 6 equally-spaced levels for
each factor in the test functions. As such this design is directly comparable
to the FFD and optimal designs, having the same budget in terms of levels for
the factors.

Two replicated MAXPRO_dis designs have also been included in the simulation
study: one with a 50% level of replication (26 base runs replicated twice,
MAXPRO_dis_50%repl) and the other with a 25% level of replication
(MAXPRO_dis_25%repl), in which the configurations to replicate have been
selected randomly. In this case, given the flexibility of the MAXPRO_dis
designs, the authors wanted to evaluate whether a lower level of replication
could be advantageous, even though by doing so a proper separation of noise
from signal 27 could be attempted only at some randomly-selected locations and
not at all configurations (as is the case for a 50% level of replication).

Table 1: Summary of the experimental designs chosen for the simulation study.

ID | Description | Software | # levels | # factors | # runs | Replication
---|---|---|---|---|---|---
CCD | Central Composite Design | R57 | 3 | 6 | 52 | 0%
BBD | Box-Behnken Design | R57 | 3 | 6 | 52 | 0%
FFD | D-optimal full factorial design | Minitab® | 6 | 6 | 52 | 0%
D_opt | D-optimal design | JMP® | 6 | 6 | 52 | 0%
I_opt | I-optimal design | JMP® | 6 | 6 | 52 | 0%
LHD_rand | Random Latin Hypercube Design | R 58 | 52 | 6 | 52 | 0%
MAXPRO | MaxPro space-filling design | R 30 | 52 | 6 | 52 | 0%
MAXPRO_dis | MaxPro discrete numeric design | R 32 | 6 | 6 | 52 | 0%
D_opt_50%repl | D-optimal design | JMP® | 6 | 6 | 52 | 50%
I_opt_50%repl | I-optimal design | JMP® | 6 | 6 | 52 | 50%
MAXPRO_dis_50%repl | MaxPro discrete numeric design | R 32 | 6 | 6 | 52 | 50%
MAXPRO_dis_25%repl | MaxPro discrete numeric design | R 32 | 6 | 6 | 52 | 25%

###  3.3 Machine learning models

Arboretti et al. 2 showed that ANNs are by far the most adopted algorithm for
the analysis of data from experimental designs when the purpose is to predict
one or more responses. Other models adopted in the literature 2 are Support
Vector Regression models (SVMs)34, Gaussian Processes (GPs) 35 and linear
models (LMs) based on quadratic regression with interactions 23. In addition
to those methods, in this study we investigate the performance of Random
Forests (RFs)36 and the algorithms contained within the Automated Machine
Learning (aml) platform offered by H2O 37, that provides a comprehensive view
of the most used ML models.

Careful tuning was carried out to choose the best configuration of the
hyperparameters that minimizes the root mean squared error (RMSE) of each
algorithm, as detailed below:

  * 1.

ANN shallow (ANN_sh): this is the simplest implementation of ANNs, counting
one hidden layer. The number of neurons was searched for in the range
$[3-12]$, and to limit overfitting, weight decay was applied with a level of
$[0-0.5]$. After some initial tuning, linear activation function was selected
for the output node. All hyperparameters were optimized by means of fivefold
Cross Validation (5-fold CV) 38 applied on the training data from the
experimental designs as detailed in Table LABEL:DOEs_SS.

  * 2.

ANN deep (ANN_dp): in this case, two to four hidden layers were considered for
the ANN, with either 6 or 12 neurons per layer. Such configurations of the
multilayered ANN were identified as the most promising in a preliminary tuning
of the models, though a concrete risk of overfitting was identified,
especially for the deepest networks. Several activation functions were also
tested, including Tanh, Maxout and Rectified Linear 39. In order to limit
overfitting of the training data, dropout 40 and both $\ell_{1}$ (Lasso) and
$\ell_{2}$ (Ridge) regularization were employed 39, making the models more
stable and less susceptible to noise. 5-fold CV was used for the selection of
the best hyperparameter configuration, as it is the most balanced solution
both in terms of computational effort and reduction in variance of the error
estimation 38.

  * 3.

Random Forest (RF): 2000 trees were grown from the training data, constituting
the RF. This configuration was obtained after an accurate initial tuning, and
it has already been used in the literature for application on DOE data 41. The
number of candidate predictors selected for each split is chosen through
5-fold CV, as such obtaining decorrelation of the base learners and improved
predictions 38.

  * 4.

Support Vector Regression model (SVM): three different kernels were considered
for SVM, namely linear
$k(\mathbb{x},\mathbb{x^{\prime}})=<\mathbb{x},\mathbb{x}^{\prime}>$,
polynomial
$k(\mathbb{x},\mathbb{x}^{\prime})=(\texttt{scale}<\mathbb{x},\mathbb{x}^{\prime}>+\texttt{offset})$degree,
and Gaussian radial basis function (RBF)
$k(\mathbb{x},\mathbb{x}^{\prime})=\texttt{exp}(-\sigma||\mathbb{x}-\mathbb{x}^{\prime}||^{2})$
where $\mathbb{x}$ and $\mathbb{x}^{\prime}$ are two input vectors 42. 5-fold
CV was employed for tuning the hyperparameters (e.g. scale, offset and degree
in the polynomial kernel or $\sigma$ in the RBF kernel).

  * 5.

Gaussian Process (GP): it is usually the model of choice for the analysis of
data from computer experiments, typically fitting data from space-filling
designs 35. Consider a $d$-dimensional design $\mathbb{D}$ with $n$ runs
$\mathbb{D}=\\{\mathbb{x}_{1},...,\mathbb{x}_{n}\\}$, where
$\mathbb{x}=(x_{1},...,x_{d})$ denotes one of the input vectors and
$\mathbb{y}=(y(\mathbb{x}_{1}),...,y(\mathbb{x}_{n}))^{T}$ is the vector of
the outputs. The GP assumption states that $y(\mathbb{x})$ is a realization
of:

$Y(\mathbb{x})=\mu(\mathbb{x})+Z(\mathbb{x})$ (1)

where $\mu(\mathbb{x})$ is the trend component and $Z(\mathbb{x})$ is
multivariate normal with mean 0 and covariance function $k$ (kernel function):

$Z(\mathbb{x})\sim\mathcal{N}(0,\,\tau^{2}k(\cdot,\cdot))$ (2)

where $\tau^{2}$ is a multiplicative constant known as the scale parameter,
defining the amplitude of the function 35. The choice of kernel function
$k(\cdot,\cdot)$ is critical, and many different options have been proposed in
the literature. For an extensive discussion on the properties of the commonly
used covariance functions for GPs, please refer to the work by Rasmussen and
Williams 43, chapter 4.

Let ${x}$ and ${x}^{\prime}$ denote two components of the input vectors
$\mathbb{x}$ and $\mathbb{x}^{\prime}$ respectively; in this article we
consider the following kernel functions 44:

    * (a)

Gaussian:
$k({x},{x}^{\prime})=\texttt{exp}(-\frac{({x}-{x}^{\prime})^{2}}{2\theta^{2}})$

    * (b)

Exponential:
$k({x},{x}^{\prime})=\texttt{exp}(-\frac{|{x}-{x}^{\prime}|}{\theta})$

    * (c)

Power-Exponential:
$k({x},{x}^{\prime})=\texttt{exp}(-(\frac{|{x}-{x}^{\prime}|}{\theta})^{t})$,
with $0<t\leq 2$

    * (d)

Matérn $5/2$:
$k({x},{x}^{\prime})=(1+\frac{\sqrt{5}|{x}-{x}^{\prime}|}{\theta}+\frac{5({x}-{x}^{\prime})^{2}}{3\theta^{2}})\texttt{exp}(-\frac{\sqrt{5}|{x}-{x}^{\prime}|}{\theta})$

    * (e)

Matérn $3/2$:
$k({x},{x}^{\prime})=(1+\frac{\sqrt{3}|{x}-{x}^{\prime}|}{\theta})\texttt{exp}(-\frac{\sqrt{3}|{x}-{x}^{\prime}|}{\theta})$

Note that here the kernel functions are reported on the 1-$d$ case, but they
can be generalized to the $d$-dimensional case by taking
$K(\mathbb{x},\mathbb{x}^{\prime})=\prod\limits_{i=1}^{d}k({x}_{i},{x}_{i}^{\prime})$
44. $\theta$ is the lengthscale parameter which defines the rate of decay of
the spatial correlation among two data points; in this article, we consider
the vectorized version of $\bm{\theta}=(\theta_{1},...,\theta_{d})$, that
allows the strength of correlation to be modulated separately by distance in
each input coordinate 35. Both $\bm{\theta}$ and $\tau^{2}$ are estimated via
maximum likelihood estimation (MLE) 44. Another relevant parameter is the
“nugget” $g$, whose main practical role is to prevent problems of inversion of
the covariance matrix due to the presence of numerical instabilities 44, 35,
but can also provide protection from violations of the modeling assumptions
(e.g. inappropriate covariance function) or other practical occurrences such
as data sparsity attributable to high dimensionality and small data sizes 45.
The nugget effect is added to the diagonal terms of the covariance matrix, and
in this article $g=10\textsuperscript{-8}\texttt{var}(\mathbb{y})$ is chosen.

The choices for the trend component $\mu(\mathbb{x})$ are either:

    * (a)

Constant, with only the intercept $\beta_{0}$.

    * (b)

Quadratic with interactions and stepwise selection of significant terms.

  * 6.

Automated Machine Learning (aml): this is a user-friendly ML software offered
by H2O. Its main advantage is that it automates the ML workflow, which
includes automatic training and tuning of many models by 5-fold CV and within
a user-specified time limit 46. The default algorithms evaluated by aml
include: RFs, Extremely Randomized Trees, Generalized Linear Models, Gradient
Boosting Machines, Deep Neural Networks and Stacked Ensembles of all the base
learners 47, 48.

  * 7.

Second-order model 23 and stepwise selection of relevant effects (LM):

$y=\beta_{0}+\sum_{i=1}^{d}{\beta_{i}x_{i}}+\sum_{i=1}^{d}{\beta_{ii}x_{i}^{2}}+\mathop{\sum\sum}_{i<j}{\beta_{ij}x_{i}x_{j}}+\epsilon$
(3)

where $\beta_{0}$, $\beta_{i}$, $\beta_{ii}$ and $\beta_{ij}$ are the
regression coefficients and $x_{i}$, $x_{j}$ are the input variables with
$i,j=1,\dots,d$.

###  3.4 Test functions

The test functions were selected from the open literature 49, mainly emulating
physical processes. All the functions were restricted to 6 active dimensions
and are evaluated on the training data, i.e. the experimental designs (Tab.
LABEL:DOEs_SS), and on a noiseless random LHD with 1000 points that is used as
test data 14. In the study the dependent variables were standardized:

$y_{n}\textsuperscript{std}=\frac{y_{n}-\overline{y}}{\sigma_{y}}$ (4)

where $y_{n}\textsuperscript{std}$ is the standardized value and $y_{n}$ is
the observed value for the $n$-th observation, $\overline{y}$ and $\sigma_{y}$
are the mean and standard deviation of $y$ respectively. For each test
function, $\overline{y}$ and $\sigma_{y}$ were estimated from 100 random LHDs
each with $500000$ data points:
$\overline{y}=(\sum_{i=1}^{100}\overline{y}_{LHD_{i}})/100$ where
$\overline{y}_{LHD_{i}}$ is the mean response over the $i$-th design, and
$\sigma_{y}=\max\sigma_{yLHD_{i}}$ with $\sigma_{yLHD_{i}}$ the standard
deviation of the response calculated over the $i$-th design ($i=1,...,100$).
Independent variables are normalized to $0-1$. A summary of the test functions
is reported in Table LABEL:test_fnct, and further details are provided in
supplemental material.

Table 2: Summary of the test functions chosen for the simulation study.

ID | Description | # dimensions
---|---|---
Borehole | models water flow through a borehole | 6
OTL circuit | models an output transformerless push-pull circuit | 6
Piston | models the circular motion of a piston within a cylinder | 6
Piston Mod | a modification of Piston function, with increased non-linearity | 6
Robot arm | models the position of a robot arm which has 3 segments | 6
Rosenbrock Function | is a popular test problem for optimization algorithms | 6
Wing weight | models a light aircraft wing | 6

###  3.5 Noise levels

In this study we investigate both the homoscedastic and heteroscedastic noise
cases. For the homoscedastic case, we assume a random normal noise component
$\epsilon\sim\mathcal{N}(0,\sigma_{hom}^{2})$, with $\sigma_{hom}$ in the
range $[0,0.5\sigma_{y}]$ 50. For the heteroscedastic case, we assume a random
normal error $\epsilon\sim\mathcal{N}(0,\sigma_{het}^{2})$ that increases
linearly with the response. That is, at a specific design configuration
$\mathbb{x}$: $\sigma_{het}=a(f(\mathbb{x})+b)$, with
$a=(\texttt{m}\sigma_{y}-0.05\sigma_{y})/(\max y-\min y)$,
$b=(0.05\sigma_{y}-a\min y)/a$, where $\min y$ and $\max y$ are the minimum
and maximum $y$ registered over the 100 LHDs respectively (section 3.4) and m
is in $[0.5,5]$ depending on the magnitude of noise 26. A summary of the
different noise levels is reported in Table LABEL:noisetab, and the impact of
different noise magnitudes is visualized in Figure 2.

Table 3: Summary of the noise structures chosen for the simulation study.

Noise $\sigma$ | Type | Description
---|---|---
$\sigma_{hom}=0\sigma_{y}$ | – | 0% noise, deterministic function
$\sigma_{hom}=0.05\sigma_{y}$ | Homoscedastic | 5% noise
$\sigma_{hom}=0.125\sigma_{y}$ | Homoscedastic | 12.5% noise
$\sigma_{hom}=0.2\sigma_{y}$ | Homoscedastic | 20% noise
$\sigma_{hom}=0.5\sigma_{y}$ | Homoscedastic | 50% noise
$\sigma_{het,min}=0.05\sigma_{y}$, | Heteroscedastic | Moderate: 5% noise at $\min y$
$\sigma_{het,max}=0.5\sigma_{y}$ | and 50% noise at $\max y$ (low5_high50)
$\sigma_{het,min}=0.05\sigma_{y}$, | Heteroscedastic | Intermediate: 5% noise at $\min y$
$\sigma_{het,max}=1\sigma_{y}$ | and 100% noise at $\max y$ (low5_high100)
$\sigma_{het,min}=0.05\sigma_{y}$, | Heteroscedastic | Severe: 5% noise at $\min y$
$\sigma_{het,max}=5\sigma_{y}$ | and 500% noise at $\max y$ (low5_high500)

Figure 2: Visualization of the impact of different noise structures and
magnitudes, both homoscedastic (on the left) and heteroscedastic (on the
right) for the test function “Piston”. The residuals are computed as the
difference between the deterministic function and the noisy function on a
1000-point LHD. Similar results are obtained for the other test functions.

###  3.6 Ranking procedure

Given the large number of designs, models, test functions and noise structures
and magnitudes, it should already be clear that a definitive ranking of the
DOEs and models with respect to the prediction error obtained on a test
dataset would be quite difficult to achieve. As such, a robust ranking
procedure that can give strong and reliable insights on the results is
required, and for this article we select a procedure based on the execution of
nonparametric permutation tests. The main advantages of this procedure include
high flexibility and the possibility of applying it without requiring any
strict assumption on the size or underlying distribution of data, while
nonetheless providing inferential information 51. Other ranking approaches
have been proposed in the literature 52, mainly parametric. However, in
comparison to the approach considered here those methods are more demanding
from the point of view of statistical assumptions (e.g. about data
distribution) that cannot always be verified in practical applications.
Furthermore, other ranking and selection procedures need to specify the
probability of correct classification and how far two ordered populations
should be from each other in the ranking 52, requirements that are not needed
in the approach adopted here 53, 54. Additionally, the considered
nonparametric ranking procedure also proved its efficacy from the point of
view of the reliability of results as it has been extensively validated by
means of simulation studies 53 and practical applications considering both
observational and experimental data in several fields 54, including medicine,
new product development and marketing studies. Recently it has also been
included as a core component in a methodology that performs variable selection
in near-infrared spectroscopy using ML models55.

In this paper the permutation tests are applied using the difference in means
as test statistics and assuming independent or paired data, depending on the
specific situation. Considering $G_{i}$ and $G_{j}$ with $i,j=1,...,C,i\neq j$
two different groups of data to be compared (e.g. the different experimental
designs), the permutation testing framework is employed to test the
directional alternative hypothesis
$\texttt{RMSE}_{G_{i}}>\texttt{RMSE}_{G_{j}}$, where RMSE is the prediction
error calculated on the test data, i.e. a noiseless random LHD with 1000 data
points 14. 2000 permutations are considered and the relevant p-values are
computed by means of the Conditional Monte Carlo (CMC) procedure described in
Pesarin and Salmaso 51. In this study we refer to the one-way MANOVA layout
(although in a nonparametric framework), meaning that interest lies in the
detection of differences between $C$ groups on $p$ dependent variables (e.g.
the test error obtained with different predictive models). Therefore, a
multivariate ranking procedure is needed and we opt for the permutation
approach for ranking multivariate populations described by Arboretti et al.
53. After the execution of all possible comparisons between the $C$ groups
over the $p$ dependent variables and the application of Fisher’s combining
function as suggested in 51, 53 for the combination of the p-values over the
$p$ variables, matrix $\mathbb{P}$ with dimension $C\times C$, that includes
the combined p-values related to the comparisons between the $C$ groups, is
obtained.

The ranking procedure 53 is as follows:

  1. 1.

Generate the matrix $\mathbb{S}$ where $S_{ij}=1$ if ${P}_{ij}\leq\alpha/2$,
$S_{ij}=0$ otherwise. ${P}_{ij}$ is an entry of $\mathbb{P}$.

  2. 2.

Compute the vector $\mathbb{r}^{D}$ of downward rank estimates
$r_{j}^{D}=1+\sum_{i=1}^{C}S_{ij}$, $j=1,...,C$.

  3. 3.

Compute the vector $\mathbb{r}^{U}$ of upward rank estimates
$r_{i}^{U}=1+\\{\\#(C-\sum_{j=1}^{C}S_{ij})>(C-\sum_{j=1}^{C}S_{i^{\prime}j}),i^{\prime}=1,...,C,i^{\prime}\neq
i\\}$, $i=1,...,C$.

  4. 4.

Compute the vector $\mathbb{r}$ of ranking estimates
$r_{i}=1+\\{\\#(r_{i}^{D}+r_{i}^{U})/2>(r_{j}^{D}+r_{j}^{U})/2,j=1,...,C,i\neq
j\\},i=1,...,C$.

In the next section we will present the final rankings for both the DOEs and
ML models considered in this work. In order to achieve a final ranking of the
different groups, two applications of the ranking procedure described above
will be required.

Consider the case in which the interest lies in ranking the experimental
designs. In the first application of the ranking procedure, both the test
function and noise setting will be fixed, meaning that a ranking of the
designs for each specific test function in a particular noise setting will be
obtained. In this case the $C$ groups refer to the experimental designs (Tab.
LABEL:DOEs_SS), while the $p$ dependent variables are the $p$ vectors of test
errors obtained with the $p=7$ predictive models (section 3.3). The groups $G$
are assumed to be independent. For each of the $p$ variables, each group $G$
will have 10 observations, that are 10 values of RMSE on the test data
obtained with 10 repetitions of the simulation in which the 5-fold CV
procedure and the predictive models are initialized with a different random
seed.

At this point, a second application of the ranking procedure is performed,
taking advantage of the results obtained in the first iteration. In this case,
the noise levels will remain fixed, and the ranking of the designs will be
performed over the different test functions, meaning that the $C$ groups will
still refer to the designs, but now each $G$ will have 7 observations that are
the ranks of the designs from the first iteration obtained on the 7 test
functions given a specific noise setting. In this case $p=1$, meaning that the
univariate version of the method in Arboretti et al. 53 is used. The final
result consists in an overall rank of each group $G$ for each specific noise
setting.

When interest lies in achieving a rank of the predictive models, the procedure
is similar: the $C$ groups will refer to the 7 predictive models (section
3.3), and the $p$ dependent variables are the $p$ vectors of test errors
obtained by fitting the algorithms on data from the $p=12$ experimental
designs (Tab. LABEL:DOEs_SS). In this case the groups $G$ cannot be assumed
independent since for a given experimental design all the models are fitted to
the same data. In all hypothesis testing procedures, $\alpha$ is assumed to be
$5\%$.

##  4 Results

This section presents the results of the simulation study and provides the
intermediate and final ranks for both the experimental designs and predictive
algorithms; sections 4.1 and 4.2 respectively.

###  4.1 Results for the experimental designs

Figure 3 shows the results in terms of the test error obtained by the models
trained on the data collected through the different experimental designs. For
the sake of brevity, the figure focuses on the results of the test function
“Piston”, considering only one level of error for both the homoscedastic and
heteroscedastic cases. Similar results are obtained for the other noise levels
and test functions, while synthetic information will be reported in the
following sections for all test functions and noise levels. It should be noted
that for replicated designs, it was not possible to fit the LM model, since
too few unique data points were present in the training data to allow the
estimation of parameters.

Figure 3: Results of the 10 iterations of the simulation study for the
function “Piston”, considering the test error (RMSE) obtained by the different
models trained on the different experimental designs for an intermediate level
of homoscedastic (left) and heteroscedastic (right) noise. The mean value for
each group is also indicated by a black star. Similar results are obtained for
the remaining test functions and noise settings.

####  4.1.1 First ranking

At this stage the first application of the multivariate ranking procedure
described in section 3.6 is conducted. Figure 1 (Appendix A) shows the results
for the function “Piston” considering all the noise settings. The ranking is
conducted for each test function, and a summary of the results is displayed in
Figures 2 and 3, with each plot referring to a specific noise setting.

####  4.1.2 Final ranking

The following step consists in the execution of an additional iteration of the
ranking procedure as described in section 3.6. The final rank provides a
synthetic result for the intermediate rankings displayed in each panel of
Figures 2 and 3, essentially carrying out a ranking of the ranks from the
first iteration. Results are displayed in Tables LABEL:tab_rank_DOE_hom and
LABEL:tab_rank_DOE_het for the homoscedastic and heteroscedastic cases
respectively. Each table should be read column-wise: the relative ranks are
computed for each noise setting.

###  4.2 Results for the predictive models

The procedure for ranking the predictive models is more or less equivalent to
the one followed for the experimental designs, with slight modifications as
described in section 3.6. Figures 1 and 2 (Appendix B) show a summary of the
intermediate ranks obtained for each test function and noise setting, both
homoscedastic and heteroscedastic. The results of the final ranks are shown in
Tables LABEL:tab_rank_MODEL_hom and LABEL:tab_rank_MODEL_het.

Table 4: Final rank of the designs for the homoscedastic noise cases.

Design | 0% | 5% | 12.5% | 20% | 50% |
---|---|---|---|---|---|---
CCD | 8 | 9 | 11 | 10 | 12 |
BBD | 6 | 6 | 4 | 4 | 4 |
FFD | 5 | 1 | 1 | 1 | 1 |
D_opt | 9 | 8 | 8 | 7 | 4 |
I_opt | 3 | 1 | 4 | 2 | 1 |
MAXPRO_dis | 1 | 1 | 1 | 2 | 4 |
MAXPRO | 2 | 1 | 1 | 4 | 1 |
LHD_rand | 3 | 1 | 6 | 6 | 8 |
D_opt_50%repl | 12 | 12 | 12 | 12 | 11 |
I_opt_50%repl | 11 | 11 | 10 | 10 | 4 |
MAXPRO_dis_50%repl | 10 | 9 | 9 | 9 | 8 |
MAXPRO_dis_25%repl | 7 | 7 | 7 | 8 | 8 |

Table 5: Final rank of the designs for the heteroscedastic noise cases.

Design | low5_high50 | low5_high100 | low5_high500 |
---|---|---|---|---
CCD | 11 | 11 | 11 |
BBD | 3 | 4 | 1 |
FFD | 5 | 7 | 8 |
D_opt | 8 | 4 | 8 |
I_opt | 3 | 1 | 6 |
MAXPRO_dis | 2 | 1 | 1 |
MAXPRO | 1 | 3 | 6 |
LHD_rand | 6 | 7 | 8 |
D_opt_50%repl | 11 | 11 | 12 |
I_opt_50%repl | 10 | 9 | 1 |
MAXPRO_dis_50%repl | 9 | 9 | 1 |
MAXPRO_dis_25%repl | 6 | 4 | 1 |

Table 6: Final rank of the models for the homoscedastic noise cases.

Model | 0% | 5% | 12.5% | 20% | 50% |
---|---|---|---|---|---|---
LM | 2 | 2 | 2 | 2 | 3 |
GP | 1 | 1 | 1 | 1 | 1 |
SVM | 4 | 3 | 2 | 3 | 2 |
RF | 5 | 5 | 5 | 5 | 5 |
ANN_sh | 3 | 3 | 4 | 4 | 3 |
ANN_dp | 7 | 7 | 5 | 7 | 7 |
aml | 6 | 5 | 5 | 5 | 6 |

Table 7: Final rank of the models for the heteroscedastic noise cases.

Model | low5_high50 | low5_high100 | low5_high500 |
---|---|---|---|---
LM | 3 | 4 | 4 |
GP | 1 | 1 | 3 |
SVM | 1 | 2 | 1 |
RF | 5 | 5 | 2 |
ANN_sh | 4 | 3 | 6 |
ANN_dp | 6 | 7 | 6 |
aml | 6 | 6 | 5 |

##  5 Discussion

In this section we discuss the results of the simulation study, focusing both
on the experimental designs (section 5.1) and predictive models (section 5.2).

###  5.1 Discussion: experimental designs

In general, and as expected, the results of the study show that the choice of
experimental design impacts on the final outcome of the analysis. From the
unprocessed results (e.g. Fig. 3) it can already be seen that not all the
designs perform equally, and some trends can already be detected (e.g. for
intermediate noise magnitudes the replicated designs do not offer an
advantage). However, the raw results clearly indicate the need for the
application of a robust ranking procedure to provide valuable synthetic
information.

The first application of the ranking procedure provides some preliminary
indications about the experimental designs. Given a test function (e.g. Fig
1), the presence of noise largely affects the efficacy of the designs: in
general, for modest levels of noise, the designs with many factor levels
appear to be favored, while as noise increases, the performance of designs
with a moderate number of levels and the addition of replicates increases.
Comparing the ranks of the designs for the different test functions appears to
confirm this trend (Fig. 2, 3). Furthermore, an interaction effect seems to
exist between the test functions and the experimental designs because, in
general, notwithstanding the level of noise, some designs tend to perform
better for some test functions than others: take for example the functions
“Robot arm” and “Rosenbrock” for which the designs with many factor levels
tend to consistently perform better than designs with only three factor
levels. We put this down to the different degrees of non-linearity of each
function. However, as the magnitude of noise increases, the situation
progressively becomes less clear and for severe levels of noise, whether
homoscedastic or heteroscedastic, it becomes burdensome for the analyst to
choose by simply considering the results of the first rank.

The results of the second rank provide clearer indications, and confirm that
not all the designs perform equally as well. Focusing on the homoscedastic
case (Tab. LABEL:tab_rank_DOE_hom), FFD, MAXPRO_dis and MAXPRO are the overall
best performers, followed closely by the I-optimal design. BBD and LHD_rand
follow, with the BBD being the best of the designs with only three levels for
each factor. The D-optimal design and the CCD are among the worst performers.
In general, the presence of replicates makes the predictions worse, but a
lower level of replication is preferable to a larger one.

The situation differs when focus turns to the heteroscedastic noise setting
(Tab. LABEL:tab_rank_DOE_het). In this case, the MAXPRO_dis is the
unquestionable winner, ranking first for the intermediate and severe noise
situations, and second in the case of moderate heteroscedasticity. The BBD,
I_opt, MAXPRO and MAXPRO_dis_25%repl follow, all behaving similarly. The
MAXPRO_dis_50%repl, D_opt, FFD, I_opt_50%repl and LHD_rand perform less well,
and the ranking is completed by the CCD and D_opt_50%repl designs.

From the analysis it emerges that the MAXPRO_dis is the overall best performer
if both the homoscedastic and heteroscedastic settings are considered.
Apparently, the space-filling property of the design can propose combinations
of the factor settings that maximize the ability of several different
predictive models to capture the non-linearity inherent in the phenomenon. At
the same time, the reduced number of factor levels makes the design robust to
the presence of noise and applicable for physical experimentation. The
MAXPRO_dis significantly outperforms all the other designs with 6 levels for
the factors, namely FFD, D_opt and I_opt. Interestingly, while such designs
perform similarly (except for D_opt) in the homoscedastic setting, in the
heteroscedastic scenario the MAXPRO_dis significantly outperforms its
competitors despite having the same number of factor levels. Particularly
relevant is the fact that not only the FFD and estimation-oriented D-optimal
designs are outperformed, but also the I-optimal design that is specifically
built to favor accurate predictions of the unknown underlying function within
the design space. The specific reasons for the outstanding performance of the
MAXPRO_dis are still unknown, but the authors believe that they are related to
the intrinsic space-filling criterion adopted by the designs. Confirmation of
the merit of such designs is that the traditional space-filling version of the
MAXPRO also outperforms its direct competitor, i.e. the LHD_rand design, both
in the homoscedastic and heteroscedastic scenarios. While the LHD_rand design
achieves good results only for rather small noise levels, the MAXPRO design
proves to be quite effective in intermediate and large noise situations as
well. However, it should again be stressed that traditional space-filling
designs are not a viable solution in many practical applications of physical
experimentation, and are reported here only as a benchmark. Furthermore, as
expected from the optimal designs, the prediction-oriented I_opt design
consistently outperforms the estimation-oriented D_opt option. This implies
that even though the considered optimal designs have historically been
developed to be used together with parametric second-order models (LM), their
properties also appear to be retained when the model for data analysis is from
the ML family.

Another surprise is the difference between the CCD and BBD. These classic
designs both have three levels for each factor, but in general the BBD
performs better, especially in settings with large noise in which it is among
the overall best performers. Apparently, the edge points of BBDs are
preferable to the corner and axial points included in CCDs. One reason might
be that for the selected test functions, which model different physical
processes, the non-linearity is better captured if points are distributed on
the edges of the design space rather than on the axes built from the center of
the design space. The predictive performance obtained using the BBD makes such
a design comparable with other experimental designs that have twice the number
of levels, especially in the heteroscedastic case, which is rather surprising.

Another finding is that, as expected, the addition of replicates provides a
valuable contribution only as the noise increases, especially in the
heteroscedastic case. However, in general, there is no proof that replicated
designs should be adopted over non-replicated designs. Apparently, the
exploration of less unique input configurations in replicated designs hinders
the ability of the predictive models to appropriately learn the behavior of
the underlying test functions, and this effect is more critical than the
possibility of separating noise from signal that is provided by the presence
of replicates. From a practical standpoint, it appears that replicated designs
should only be chosen if the underlying phenomenon presents severe
heteroscedasticity. Again, the designs from the MAXPRO_dis family outperform
both the competitors (D_opt_50%repl and I_opt_50%repl), although sometimes
only marginally.

###  5.2 Discussion: predictive models

The unprocessed results shown in Figure 3 already convey the idea that the
choice of model for data analysis is critical to the value of final
predictions: given an experimental design, a “bad model” can easily provide a
test error which is twice as large as the one achieved with a “good model”.

Considering the results of the first rank, a very strong, clear pattern
emerges indicating that the GP model is consistently better than all its
competitors for all the test functions, especially in low and intermediate
homoscedastic noise settings (Fig. 1). The situation is more uncertain if we
have to decide what algorithms follow in the rank.

For moderate and intermediate levels of heteroscedasticity, the GP is
confirmed as the best or one of the best algorithms, but as the inhomogeneity
in noise increases, some algorithms, namely SVM and RF, significantly
outperform the GP (Fig. 2).

The results from the second application of the ranking procedure (Tab.
LABEL:tab_rank_MODEL_hom) corroborate our first impressions; in the
homoscedastic case, the GP is ranked first for all the noise levels
considered. The LM, SVM and ANN_sh follow, while the RF, ANN_dp and aml behave
poorly in all the analyzed scenarios (Tab. LABEL:tab_rank_MODEL_hom).

A similar situation is shown in the final results for the heteroscedastic
setting (Tab. LABEL:tab_rank_MODEL_het), but this time the SVM performs
slightly better than the GP. Moreover, as the heteroscedasticity becomes
severe, both the SVM and RF outperform the GP model. Behind the GP and SVM,
the LM, RF, and ANN_sh achieve acceptable results, while, as for the
homoscedastic setting, the aml and ANN_dp perform poorly.

A general comment that can be made is that, within the scope of the present
study, the GP appears to be the preferable choice for data analysis. This
justifies the widespread adoption of the model for the analysis of computer
experiments. However, in this case the algorithm proves its merits also when
fitting designs with a moderate number of factor levels, mimicking the
situation of physical experimentation. Nonetheless, attention should be paid
to the training data because in the case of severe heteroscedastic noise,
other options appear superior.

The SVM and LM also accomplish acceptable results, with the SVM the best
performer in the heteroscedastic setting and, as expected, the LM performing
well only in the homoscedastic case. Another finding worthy of note is the
different performances of RF for homogeneous and inhomogeneous noise, as the
algorithm increases its performance as the impact of heteroscedasticity
increases, becoming one of the best options.

The results for the ANN algorithms, that we recall being by far the most
adopted ML algorithms for analysis of DOE data 2, are also interesting. The
ANN_sh is never among the best performers, neither in the homoscedastic nor
heteroscedastic setting, while the ANN_dp is the overall poorest performer
among the models considered. The aml also performs badly, however it should be
pointed out that both the algorithms contained in the aml platform and ANN_dp
are tools that take advantage of big data, while in this study we only
consider small data for training the algorithms.

Finally, we remind the reader that in this study we only considered the
prediction error as a valuable metric to rank the models, while other
indicators, such as the possibility of providing a quantification of
prediction uncertainty and model interpretability are crucial aspects that
should also be considered 41, 2. To this end, we mention that GPs natively
provide a quantification of uncertainty and are recognized as being fairly
interpretable in the ML literature 43, 56.

##  6 Conclusions

In this paper we introduce a procedure and a simulation study to provide an
empirical assessment of the most promising experimental designs and predictive
models to be used jointly for data collection and analysis. The focus is on
physical experimentation therefore we mainly consider test functions that
emulate physical processes, use small data and favor designs with a moderate
number of factor levels. Several designs, predictive models, test functions
and noise settings are considered, thus the study covers circumstances which
are encountered in practical industrial applications, making the analysis
highly relevant for practitioners.

The results show the significant potential of designs from the MAXPRO family
30, 32, especially when the factors are restricted to the discrete numeric
case, reducing the number of levels to fit the budget of the experimenter.
Overall, these designs outperform their competitors from the literature even
when the same number of factor levels is considered, meaning that the MAXPRO
criterion can maximize the ability of the predictive models to capture the
non-linearity in the underlying phenomenon. Other good designs are the
I-optimals, as they provide only marginally worse performances than the MAXPRO
discrete numeric designs in both the homoscedastic and heteroscedastic
situations. Considering ML models, the GP appears to be the overall best
performer, but other methods provide better predictions in the case of severe
heteroscedasticity.

Nevertheless, it should be remembered that the choice of experimental design
is more critical than the choice of model for data analysis because the first
must be chosen before the project starts, while the second can be tuned after
the experimental data has been collected. Thus, several methods can be
compared to find the best performer for the specific set of data.

##  7 Acknowledgments

The authors gratefully thank Fondazione Cariparo for partially supporting this
research.

##  Appendix A Appendix A

figureIn this section we report the results of the first rank for the
experimental designs.

Figure 1: Summary of the ranks of the experimental designs for the function
“Piston”. The homoscedastic cases are reported on the left and the
heteroscedastic settings are reported on the right. Figure 2: Summary of the
ranks of the experimental designs for each test function and all homoscedastic
noise settings. Figure 3: Summary of the ranks of the experimental designs for
each test function and all heteroscedastic noise settings.

##  Appendix B Appendix B

figureIn this section we report the results of the first rank for the
predictive models.

Figure 1: Summary of the ranks of the predictive models for each test function
and all homoscedastic noise settings. Figure 2: Summary of the ranks of the
predictive models for each test function and all heteroscedastic noise
settings.

## References

  * Myers et al. [2004] R. H. Myers, D. C. Montgomery, G. G. Vining, C. M. Borror, S. M. Kowalski,  Response surface methodology: A retrospective and literature survey,  Journal of Quality Technology 36 (2004) 53–77. doi:10.1080/00224065.2004.11980252.
  * Arboretti et al. [2021] R. Arboretti, R. Ceccato, L. Pegoraro, L. Salmaso,  Design of experiments and machine learning for product innovation: a systematic literature review,  Quality and Reliability Engineering International (2021) 1–26. doi:10.1002/qre.3025.
  * Box and Wilson [1951] G. E. P. Box, K. B. Wilson,  On the experimental attainment of optimum conditions,  Journal of the Royal Statistical Society. Series B (Methodological) 13 (1951) 1–45. URL: http://www.jstor.org/stable/2983966.
  * Box and Behnken [1960] G. E. Box, D. W. Behnken,  Some new three level designs for the study of quantitative variables,  Technometrics 2 (1960) 455–475.
  * Taguchi and Wu [1979] G. Taguchi, Y. Wu, Introduction to off-line quality control, Central Japan Quality Control Assoc., 1979\.
  * Taguchi [1986] G. Taguchi, Introduction to quality engineering: designing quality into products and processes, Technical Report, 1986.
  * Taguchi [1987] G. Taguchi, System of experimental design; engineering methods to optimize quality and minimize costs, Technical Report, 1987.
  * Joseph [2016] V. R. Joseph,  Space-filling designs for computer experiments: A review,  Quality Engineering 28 (2016) 28–35. doi:10.1080/08982112.2015.1100447.
  * Box [1999] G. E. Box,  Statistics as a catalyst to learning by scientific method part ii—a discussion,  Journal of Quality Technology 31 (1999) 16–29.
  * Vining [2011] G. Vining,  Technical advice: Design of experiments, response surface methodology, and sequential experimentation,  Quality Engineering 23 (2011) 217–220. doi:10.1080/15226514.2011.555280.
  * Jensen [2018] W. A. Jensen,  Open problems and issues in optimal design,  Quality Engineering 30 (2018) 583–593.
  * Cortes et al. [2018] L. A. Cortes, J. R. Simpson, P. A. Parker,  Response surface split-plot designs: A literature review,  Quality and Reliability Engineering International 34 (2018) 1374–1389.
  * Jones and Nachtsheim [2011] B. Jones, C. J. Nachtsheim,  A class of three-level designs for definitive screening in the presence of second-order effects,  Journal of Quality Technology 43 (2011) 1–15.
  * Chen et al. [2016] H. Chen, J. L. Loeppky, J. Sacks, W. J. Welch, et al.,  Analysis methods for computer experiments: How to assess and what counts?,  Statistical science 31 (2016) 40–60.
  * Ilzarbe et al. [2008] L. Ilzarbe, M. J. Álvarez, E. Viles, M. Tanco,  Practical applications of design of experiments in the field of engineering: a bibliographical review,  Quality and Reliability Engineering International 24 (2008) 417–428.
  * Morris and Mitchell [1995] M. D. Morris, T. J. Mitchell,  Exploratory designs for computational experiments,  Journal of statistical planning and inference 43 (1995) 381–402.
  * Inc. [2019] S. I. Inc., JMP® 14 Design of Experiments Guide, SAS Institute Inc., 2019\.
  * Zahran et al. [2003] A. Zahran, C. M. Anderson-Cook, R. H. Myers,  Fraction of design space to assess prediction capability of response surface designs,  Journal of Quality Technology 35 (2003) 377–386. URL: https://doi.org/10.1080/00224065.2003.11980235. doi:10.1080/00224065.2003.11980235. arXiv:https://doi.org/10.1080/00224065.2003.11980235.
  * Rodriguez et al. [2010] M. Rodriguez, B. Jones, C. M. Borror, D. C. Montgomery,  Generating and assessing exact g-optimal designs,  Journal of quality technology 42 (2010) 3–20.
  * Goos et al. [2016] P. Goos, B. Jones, U. Syafitri,  I-optimal design of mixture experiments,  Journal of the American Statistical Association 111 (2016) 899–911.
  * Jones et al. [2021] B. Jones, K. Allen-Moyer, P. Goos,  A-optimal versus d-optimal design of screening experiments,  Journal of Quality Technology 53 (2021) 369–382. URL: https://doi.org/10.1080/00224065.2020.1757391. doi:10.1080/00224065.2020.1757391. arXiv:https://doi.org/10.1080/00224065.2020.1757391.
  * Meyer and Nachtsheim [1995] R. K. Meyer, C. J. Nachtsheim,  The coordinate-exchange algorithm for constructing exact optimal experimental designs,  Technometrics 37 (1995) 60–69. URL: https://www.tandfonline.com/doi/abs/10.1080/00401706.1995.10485889. doi:10.1080/00401706.1995.10485889. arXiv:https://www.tandfonline.com/doi/pdf/10.1080/00401706.1995.10485889.
  * Montgomery [2017] D. C. Montgomery, Design and analysis of experiments, John wiley & sons, 2017.
  * Ankenman et al. [2010] B. Ankenman, B. L. Nelson, J. Staum,  Stochastic kriging for simulation metamodeling,  Operations Research 58 (2010) 371–382. URL: http://www.jstor.org/stable/40605923.
  * Boukouvalas et al. [2014] A. Boukouvalas, D. Cornford, M. Stehlík,  Optimal design for correlated processes with input-dependent noise,  Computational Statistics & Data Analysis 71 (2014) 1088–1102. URL: https://www.sciencedirect.com/science/article/pii/S0167947313003484. doi:https://doi.org/10.1016/j.csda.2013.09.024.
  * Jalali et al. [2017] H. Jalali, I. Van Nieuwenhuyse, V. Picheny,  Comparison of kriging-based algorithms for simulation optimization with heterogeneous noise,  European Journal of Operational Research 261 (2017) 279–301.
  * Binois et al. [2019] M. Binois, J. Huang, R. B. Gramacy, M. Ludkovski,  Replication or exploration? sequential design for stochastic simulation experiments,  Technometrics 61 (2019) 7–23.
  * McKay et al. [1979] M. D. McKay, R. J. Beckman, W. J. Conover,  A comparison of three methods for selecting values of input variables in the analysis of output from a computer code,  Technometrics 21 (1979) 239–245. URL: http://www.jstor.org/stable/1268522.
  * Stein [1987] M. Stein,  Large sample properties of simulations using latin hypercube sampling,  Technometrics 29 (1987) 143–151.
  * Joseph et al. [2015] V. R. Joseph, E. Gul, S. Ba,  Maximum projection designs for computer experiments,  Biometrika 102 (2015) 371–380.
  * Li et al. [2006] X. Li, N. Sudarsanam, D. D. Frey,  Regularities in data from factorial experiments,  Complexity 11 (2006) 32–45.
  * Joseph et al. [2020] V. R. Joseph, E. Gul, S. Ba,  Designing computer experiments with multiple types of factors: The maxpro approach,  Journal of Quality Technology 52 (2020) 343–354.
  * Krishna et al. [2021] A. Krishna, V. R. Joseph, S. Ba, W. A. Brenneman, W. R. Myers,  Robust experimental designs for model calibration,  Journal of Quality Technology 0 (2021) 1–12. URL: https://doi.org/10.1080/00224065.2021.1930618. doi:10.1080/00224065.2021.1930618. arXiv:https://doi.org/10.1080/00224065.2021.1930618.
  * Hastie et al. [2009] T. Hastie, R. Tibshirani, J. Friedman, The elements of statistical learning: data mining, inference, and prediction, Springer Science & Business Media, 2009\.
  * Gramacy [2020] R. B. Gramacy, Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences, Chapman Hall/CRC, Boca Raton, Florida, 2020. http://bobby.gramacy.com/surrogates/.
  * Breiman [2001] L. Breiman,  Random forests,  Machine learning 45 (2001) 5–32.
  * LeDell and Poirier [2020] E. LeDell, S. Poirier,  H2o automl: Scalable automatic machine learning,  in: Proceedings of the AutoML Workshop at ICML, volume 2020, 2020.
  * James et al. [2013] G. James, D. Witten, T. Hastie, R. Tibshirani, An introduction to statistical learning, volume 112, Springer, 2013.
  * Candel and LeDell [2021] A. Candel, E. LeDell, Deep learning with H2O, 2021. URL: http://docs.h2o.ai/h2o/resources/, sixth Edition.
  * Srivastava et al. [2014] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, R. Salakhutdinov,  Dropout: a simple way to prevent neural networks from overfitting,  The journal of machine learning research 15 (2014) 1929–1958.
  * Arboretti et al. [2021] R. Arboretti, R. Ceccato, L. Pegoraro, L. Salmaso, C. Housmekerides, L. Spadoni, E. Pierangelo, S. Quaggia, C. Tveit, S. Vianello,  Machine learning and design of experiments with an application to product innovation in the chemical industry,  Journal of Applied Statistics 0 (2021) 1–26. URL: https://doi.org/10.1080/02664763.2021.1907840. doi:10.1080/02664763.2021.1907840.
  * Karatzoglou et al. [2004] A. Karatzoglou, A. Smola, K. Hornik, A. Zeileis,  kernlab – an S4 package for kernel methods in R,  Journal of Statistical Software 11 (2004) 1–20. URL: http://www.jstatsoft.org/v11/i09/.
  * Rasmussen and Williams [2006] C. E. Rasmussen, C. K. I. Williams, Gaussian Processes for Machine Learning, The MIT Press, 2006.
  * Roustant et al. [2012] O. Roustant, D. Ginsbourger, Y. Deville,  Dicekriging, diceoptim: Two r packages for the analysis of computer experiments by kriging-based metamodeling and optimization,  Journal of Statistical Software, Articles 51 (2012) 1–55. URL: https://www.jstatsoft.org/v051/i01. doi:10.18637/jss.v051.i01.
  * Gramacy and Lee [2012] R. B. Gramacy, H. K. Lee,  Cases for the nugget in modeling computer experiments,  Statistics and Computing 22 (2012) 713–722.
  * H2O [2021] H2O, Automl: Automatic machine learning, 2021\. URL: https://docs.h2o.ai/h2o/latest-stable/h2o-docs/automl.html.
  * Cook [2016] D. Cook, Practical machine learning with H2O: powerful, scalable techniques for deep learning and AI, ” O’Reilly Media, Inc.”, 2016.
  * Landry [2020] M. Landry, Machine Learning with R and H2O, 2020\. URL: http://h2o-release.s3.amazonaws.com/h2o/rel-zeno/2/docs-website/h2o-docs/booklets/RBooklet.pdf, seventh Edition.
  * Surjanovic and Bingham [2021] S. Surjanovic, D. Bingham, Virtual library of simulation experiments: Test functions and datasets, 2021. URL: http://www.sfu.ca/~ssurjano.
  * Picheny et al. [2013] V. Picheny, T. Wagner, D. Ginsbourger,  A benchmark of kriging-based infill criteria for noisy optimization,  Structural and Multidisciplinary Optimization 48 (2013) 607–626.
  * Pesarin and Salmaso [2010] F. Pesarin, L. Salmaso, Permutation tests for complex data: theory, applications and software, John Wiley & Sons, 2010.
  * Gupta and Panchapakesan [2002] S. S. Gupta, S. Panchapakesan, Multiple Decision Procedures, Society for Industrial and Applied Mathematics, 2002. URL: https://epubs.siam.org/doi/abs/10.1137/1.9780898719161. doi:10.1137/1.9780898719161. arXiv:https://epubs.siam.org/doi/pdf/10.1137/1.9780898719161.
  * Arboretti et al. [2014] R. Arboretti, S. Bonnini, L. Corain, L. Salmaso,  A permutation approach for ranking of multivariate populations,  Journal of Multivariate Analysis 132 (2014) 39–57.
  * Corain et al. [2016] L. Corain, R. Arboretti, S. Bonnini, Ranking of Multivariate Populations: A Permutation Approach with Applications, Chapman Hall/CRC, Boca Raton, Florida, 2016. doi:10.1201/b19673.
  * Arboretti et al. [2021] R. Arboretti, R. Ceccato, L. Pegoraro, L. Salmaso,  Interval selection: A case-study-based approach,  Applied Stochastic Models in Business and Industry (2021).
  * Liu et al. [2020] H. Liu, Y.-S. Ong, X. Shen, J. Cai,  When gaussian process meets big data: A review of scalable gps,  IEEE transactions on neural networks and learning systems 31 (2020) 4405–4423.
  * Groemping [2020] U. Groemping, DoE.wrapper: Wrapper Package for Design of Experiments Functionality, 2020. URL: https://CRAN.R-project.org/package=DoE.wrapper, r package version 0.11.
  * Dupuy et al. [2015] D. Dupuy, C. Helbert, J. Franco,  DiceDesign and DiceEval: Two R packages for design and analysis of computer experiments,  Journal of Statistical Software 65 (2015) 1–38. URL: http://www.jstatsoft.org/v65/i11/.