Preprint 4 [gray]0.90

11institutetext: MercadoLibre Inc.
11email: {hernan.vazquez,jorge.sanchez,
<EMAIL_ADDRESS>
# Integrating Hyperparameter Search into Model-Free AutoML with Context-Free
Grammars

Hernán Ceferino Vázquez  Jorge Sanchez  Rafael Carrascosa

###### Abstract

Automated Machine Learning (AutoML) has become increasingly popular in recent
years due to its ability to reduce the amount of time and expertise required
to design and develop machine learning systems. This is very important for the
practice of machine learning, as it allows building strong baselines quickly,
improving the efficiency of the data scientists, and reducing the time to
production. However, despite the advantages of AutoML, it faces several
challenges, such as defining the solutions space and exploring it efficiently.
Recently, some approaches have been shown to be able to do it using tree-based
search algorithms and context-free grammars. In particular, GramML presents a
model-free reinforcement learning approach that leverages pipeline
configuration grammars and operates using Monte Carlo tree search. However,
one of the limitations of GramML is that it uses default hyperparameters,
limiting the search problem to finding optimal pipeline structures for the
available data preprocessors and models. In this work, we propose an extension
to GramML that supports larger search spaces including hyperparameter search.
We evaluated the approach using an OpenML benchmark and found significant
improvements compared to other state-of-the-art techniques.

###### Keywords:

Automated Machine Learning Hyperparameter Optimization Context-free Grammars
Monte Carlo Tree Search Reinforcement Learning Model-Free

##  1 Introduction

The practice of Machine Learning (ML) involves making decisions to find the
best solution for a problem among a myriad of candidates. Machine learning
practitioners often find themselves in a process of trial and error in order
to determine optimal processing pipelines (data preprocessing, model
architectures, etc.) and their configurations (hyperparameters). These
decisions go hand in hand with the many constraints imposed by the actual
application domain and have a great impact on the final system performance.
These constraints can be either intrinsic, such as working with heterogeneous
data, extrinsic, such as requiring using limited computational resources
and/or short development cycles, or both. Solving real-world problems by means
of ML is a time-consuming and difficult process.

In this context, the development of methods and techniques that allow for the
automated generation of machine learning solutions [28] is a topic of great
practical importance, as they help reduce the time and expertise required for
the design and evaluation of a broad range of ML-based products and
applications.

Besides the benefits it brings to the practice of machine learning, AutoML
poses several challenges. In our case, these challenges are how to define the
solution space (the set of possible pipeline configurations, i.e. pipeline
structure and its hyperparameters) in a way that is both comprehensive and
computationally feasible, and how to explore it efficiently both in terms of
time and resources.

Many techniques have been developed over the years to address these problems,
most of which rely on parametric models to guide the exploration of the search
space. However, such model-based solutions require a characterization of
configuration space and its possible solutions via meta-features, which are
not always easy to define and evolve over time, adding an extra (meta)
learning step as well as an additional layer of hyperparameters (those of the
meta-model) [26]. An alternative to this family of models corresponds to the
model-free counterpart. These approaches rely on different strategies to guide
the search stage, in which evolutionary [21, 12]and tree-based techniques [22,
19] appear as the most prominent ones.

A particularly interesting approach consists of using context-free grammars
for the search space specification and an efficient search strategy over the
production tree [30] derived from it [19, 8, 29]. A specific algorithm that
has been used in this case is Monte Carlo Tree Search (MCTS), which allows for
a more efficient exploration of the search space [22]. In particular, [29]
proposes GramML, a simple yet effective model-free reinforcement learning
approach based on an adaptation of the MCTS for trees and context-free
grammars that show superior performance compared to the state-of-the-art.
However, one of the limitations of their model is that it disregards
hyperparameters, leaving them to their default values. This effectively
reduces the search space pipeline architectures consisting of data
preprocessors and generic model configurations.

In this work, we propose an extension to GramML that supports larger search
spaces by incorporating hyperparameters into the model. We extend the basic
MCTS algorithm to account for this increase in complexity, by incorporating a
pruning strategy and the use of non-parametric selection policies. In
addition, we run experiments to evaluate the performance of our approach and
compare it with state-of-the-art techniques. The results show that our method
significantly improves performance and demonstrates the effectiveness of
incorporating hyperparameter search into grammar-based AutoML.

##  2 Related Work

This section provides a succinct overview of prior studies related to AutoML.
The main representative for solving this class of problems is AutoSklearn
[13], which is based on the Sequential Model-Based Algorithm Configuration
(SMAC) algorithm [16] and incorporates Bayesian optimization (BO) and
surrogate models to determine optimal machine learning pipelines. In
AutoSklearn hyperparameters and pipeline selection problem are tackled as a
single, structured, joint optimization problem named CASH (Combined Algorithm
Selection and Hyperparameter optimization [25]). Recent research has produced
model-based approaches that exhibit comparable performance to AutoSklearn. For
instance, MOSAIC [22] leverages a surrogate model formulation in conjunction
with a Monte Carlo Tree Search (MCTS) to search the hypothesis space. MOSAIC
tackle both pipeline selection and hyperparameter optimization problems using
a coupled hybrid strategy: MCTS is used for the first and Bayesian
optimization for the latter, where the coupling is enforced by using a
surrogate model. AlphaD3M [9] also builds on MCTS and uses a neural network as
surrogate model, considering only partial pipelines rather than whole
configurations. The difference with MOSAIC is that the pipeline selection is
decoupled from the hyperparameter search, it first looks for the best
pipelines and then performs hyperparameter optimization of the best pipelines
using SMAC. [8] studied the use of AlphaD3M in combination with a grammar,
relying on pre-trained models for guidance during the search. Another
approach, PIPER [19], utilizes a context-free grammar and a greedy heuristic,
demonstrating competitive results relative to other techniques. In PIPER, the
best pipeline structure is found and then its hyperparameters are searched
with a standard CASH optimizer. Other approaches include the use of genetic
programming [21, 23] and hierarchical planning [20, 17] techniques.

##  3 Background

A formal grammar is a collection of symbols and grammar rules that describe
how strings are formed in a specific formal language [24]. There are various
types of grammar, including context-free grammars [6]. In context-free
grammars, the grammar rules consist of two components: the left-hand side of a
grammar rule is always a single non-terminal symbol. These symbols never
appear in the strings belonging to the language. The right-hand side contains
a combination of terminal and non-terminal symbols that are equivalent to the
left-hand side. Therefore, a non-terminal symbol can be substituted with any
of the options specified on the right-hand side.

Once specified, the grammar provide a set of rules by which non-terminal
symbols can be substituted to form _strings_. These strings belong to the
language represented by the grammar. This generation process (production task)
is performed by starting from an initial non-terminal symbol and replacing it
recursively using the grammar rules available. An important property, in this
case, is that the set of rules can be expressed in the form of a tree and the
generation of any given production can be seen as a path from the root to a
particular leaf of this tree. The generation process can be thus seen as a
search problem over the tree induced by the grammar.

In the context of the AutoML problem, the grammar is used to encode the
different components and their combination within a desired pipeline. For
instance, we can set a grammar that accounts for different data preprocessing
strategies, like _PREPROCESSING := ”NumericalImputation” or SCALING :=
”MinMaxScaler” or ”StandardScaler”_ or the type of weight regularization used
during learning, as _PENALTY := ”l1” or ”l2”_. Finding optimal ML pipelines
corresponds to searching for a path from the root (the initial symbol of the
grammar) to a leaf of the tree (the last rule evaluated) that performs best in
the context of a given ML problem. Interestingly, we could also encode
hyperparameter decisions within the same grammar, _e.g_. by listing the
possible values (or ranges) that a particular parameter can take. By doing so,
we grow the tree size (configuration space of all possible ML pipelines and
their configurations) exponentially.

From the above, a key component of grammar-based approaches to AutoML is the
search algorithm, as it has to be efficient and scale gracefully with the size
of the tree. In the literature, one of the most reliable alternatives is the
Monte Carlo Tree Search (MCTS) algorithm [18]. MCTS is a search algorithm that
seeks to expand the search tree by using random sampling heuristics. It
consists of four steps: _selection_ , _expansion_ , _simulation_ and
_backpropagation_. During the selection step, the algorithm navigates the
search tree to reach a node that has not yet been explored. Selection always
starts from the root node and, at each level, selects the next node to visit
according to a given _selection policy_. The expansion, simulation, and
backpropagation steps ensure that the tree is expanded according to the rules
available at each node and, once a terminal node is reached, that the
information about the reliability of the results obtained at this terminal
state is propagated back to all the nodes along the path connecting this node
to the root. As a reference, Algorithm 1 illustrates the general form of the
MCTS algorithm.

Algorithm 1 General MCTS approach

1:function Search($s_{0}$)

2: create root node $a_{0}$ with state $s_{0}$

3: while within budget do

4: $a_{i}\leftarrow\text{Selection}(a_{0},SelectionPolicy)$

5: $a_{i}\leftarrow\text{Expansion}(a_{i})$

6: $\Delta\leftarrow\text{Simulation}(s(a_{i}))$

7: $\text{Backpropagation}(a_{i},\Delta)$

8:return $\text{BestChild}(a_{0},SelectionPolicy)$

In general, MCTS is a powerful method that allows to explore the search tree
in a systematic and efficient way and can be used for AutoML. However,
balancing the exploration and exploitation trade-off depends on find best
selection policy to use. Furthermore, when the tree grows up exponentially as
in the case of hyper-parameter search problem, finding good solutions in an
acceptable time can be challenging.

##  4 Hyperparameter Search in Grammar-Based AutoML

In grammar-based AutoML, the objective is to generate pipeline configurations
using the available production rules that encode the choices for a particular
non-terminal symbol. In our case, we extend the GramML111The Extended GramML
source code is available at: github.com/mercadolibre/fury_gramml-with-
hyperparams-search. formulation of [29] to account for the hyperparameters
associated with the different ML components that compose an ML pipeline. For
continuous numerical hyperparameters, we use the same sampling intervals as
AutoSklearn for a fairer comparison. In the case of categorical or discrete
parameters, we incorporate their values directly into the grammar.

Using MCTS the grammar is traversed generating productions and looking for the
best pipeline configurations. The algorithm can be stopped at any point, and
the current best pipeline configuration is returned. The longer the algorithm
runs, the more productions will be generated and evaluated, and the more
likely that the recommended pipeline is indeed the optimal one. The process of
generating grammar productions and evaluating pipeline configurations consists
of the traditional four steps of MCTS; however, it has some considerations
given the nature of the problem. In our setting, the only valid configurations
are leaf nodes, as they contain the complete configuration object. Any attempt
to execute a partial configuration would result in an error. Additionally,
each iteration of the algorithm (commonly referred to as an episode in RL)
ends in a leaf node or terminal state. These differences lead to some changes
in the original algorithm.

Algorithm 2 MCTS for AutoML

1:function Search($s_{0}$)

2: $a_{i}\leftarrow a_{0}$

3: while used budget $<$ budget $\And$ $a_{i}$ is not Terminal do

4: $a_{i}\leftarrow\text{Selection}(a_{0},SelectionPolicy)$

5: if $a_{i}$ is not Terminal then

6: $a_{i}\leftarrow\text{Expansion}(a_{i})$

7: $\Delta\leftarrow\text{Simulation}(s(a_{i}))$

8: else

9: $\Delta\leftarrow\text{GetReward}(s(a_{i}))$

10: $\text{Backpropagation}(a_{i},\Delta)$

11: $p_{i}^{*}\leftarrow\text{BestLeaf}(a_{0},SelectionPolicy)$

12: $\text{PruneTree}(p_{i}^{*})$

13: return $p_{i}^{*}$

Algorithm 2 illustrates an adaptation of the MCTS algorithm for AutoML. The
Search function is called by the main program until a budget is exhausted. The
budget can be specified in terms of the number of iterations or running time.
Each call to the Search function returns a pipeline configuration. This
function is executed until either a pipeline configuration is found or the
budget is exhausted.

In each step of the algorithm, if the selected node is not terminal, the
search function performs the expansion and simulation steps. In the expansion
step, the tree is navigated from the selected node to find a node that has not
been expanded. Once found, it is expanded, and the algorithm continues. In the
simulation step, the algorithm randomly selects a path from the selected node
to a leaf. Upon reaching a leaf, a reward function is computed and returned.
In the case of the selected node being terminal, this reward is used for
backpropagation. Finally, the selected leaf is pruned from the tree. If all
leaves of a branch are pruned, then the branch is also pruned.

Selection and backpropagation steps depend on the actual selection policy.
Selection policies are important as they guide the search, regulating the
trade-off between exploration and exploitation. In what follows, we present
three different alternatives to this combination that appear as particularly
suited to our problem.

#### Upper Confidence Bound applied to Trees (UTC).

UTC [18] is a standard algorithm that seeks to establish a balance between
exploration and exploitation to find the best option among a set of
candidates. It uses a constant C as the only parameter that regulates this
trade-off. UCT is based on the computation of two statistics, a cumulative sum
of the rewards obtained and the number of visits for each node. base on these
values, we compute a score for each node as:

$v(a_{i})=\frac{reward(a_{i})}{visits(a_{i})}+C\sqrt{\frac{\log{\left(visits(a_{i})\right)}}{visits(parent(a_{i})))}},$
(1)

The selection policy corresponds to choosing the node with the highest score.
Algorithm 3 shows the process by which rewards and visits are backpropagated
from a leaf up to the root.

Algorithm 3 Backpropagation for UCT

1:function Backpropagation($a$, $\Delta$)

2: $a_{i}\leftarrow a$ $\triangleright$ a is a leaf node

3: while $a_{i}$ is not Null do

4: $reward_{i}\leftarrow reward_{i}+\Delta$

5: $visits_{i}\leftarrow visits_{i}+1$

6: $a_{i}\leftarrow parent(a_{i})$ $\triangleright$ if $a_{i}$ is the root
node parent is Null

#### Bootstrap Thompson Sampling (BTS).

In Thompson Sampling (TS) the basic idea is to select actions stochastically
[1]. TS requires being able to compute the exact posterior [5]. In cases where
this is not possible, the posterior distribution can be approximated using
bootstrapping [11]. This method is known as Bootstrap Thompson Sampling (BTS)
[10]. In BTS the bootstrap distribution is a collection of bootstrap
replicates $j\in\\{1,...,J\\}$, where $J$ is a tunable hyperparameter that
balances the trade-off between exploration and exploitation. A smaller $J$
value results in a more greedy approach, while a larger $J$ increases
exploration, but has a higher computational cost [10]. Each replicate $j$ has
some parameters $\theta$, which are used to estimate $j$’s expected utility
given some prior distribution $P(\theta)$. At decision time, the bootstrap
distribution for each node is sampled, and the child with the highest expected
utility is selected. During backpropagation, the distribution parameters are
updated by simulating a coin flip for each replicate $j$. In our scenario, we
use a Normal distribution as prior and two parameters $\alpha$ y $\beta$ (by
consistency to other works [15].) If the coin flip comes up heads, the
$\alpha$ and $\beta$ parameters for $j$ are re-weighted by adding the observed
reward to $\alpha_{j}$ and the value $1$ to $\beta_{j}$. Algorithm 4
illustrates the backpropagation step for the BTS strategy.

Algorithm 4 Backpropagation for BTS

1:function Backpropagation($a$, $\Delta$)

2: $a_{i}\leftarrow a$ $\triangleright$ a is a leaf node

3: while $a_{i}$ is not Null do

4: for $j\in J$ do

5: sample $d_{j}$ from Bernoulli($\nicefrac{{1}}{{2}}$)

6: if $d_{j}=1$ then

7: $\alpha_{a_{i}}\leftarrow\alpha_{a_{i}}+\Delta$

8: $\beta_{a_{i}}\leftarrow\beta_{a_{i}}+1$

9: $a_{i}\leftarrow parent(a_{i})$ $\triangleright$ if $a_{i}$ is the root
node parent is Null

Regarding the selection policy, we estimate the value of the nodes as
$\alpha_{a_{i}}/\beta_{a_{i}}$, which represents the largest point estimate of
our prior (the mean). Afterwards, we select a child with a probability
proportional to that value.

#### Tree Parzen Estimator (TPE).

TPE is a widely used decision algorithm in hyperparameter optimization [2]
that can be easily adapted to work on trees. TPE defines two densities,
$l(a_{i})$ and $g(a_{i})$, where $l(a_{i})$ is computed by using the
observations $\\{a_{i}\\}$ such that corresponding reward $y=reward(a_{i})$ is
less than a threshold $y^{*}$, while $g(a_{i})$ is computed using the
remaining observations. The value of $y^{*}$ is chosen as the
$\gamma$-quantile of the observed $y$ values.

To compute the functions $g(a_{i})$ and $l(a_{i})$, it is necessary to keep
track of the rewards obtained from each node instead of accumulating them as
in the case of UCT. Algorithm 5 shows the overall process.

Algorithm 5 Backpropagation for TPE

1:function Backpropagation($a$, $\Delta$)

2: $a_{i}\leftarrow a$ $\triangleright$ a is a leaf node

3: while $a_{i}$ is not Null do

4: $[reward(a_{i})]\leftarrow[reward(a_{i})]+[\Delta]$

5: $a_{i}\leftarrow parent(a_{i})$ $\triangleright$ if $a_{i}$ is the root
node parent is Null

The tree-structured induced by the partitions used to compute $l$ and $g$
makes it easy to evaluate candidates based on the ratio
$g(a_{i})\textup{/}l(a_{i})$, similar to the selection policy used in BTS.

##  5 Experiments

To evaluate our work, we conducted experiments in two parts. In the first
part, we performed an ablation study of each of the non-parametric selection
policies proposed in the context of the proposed MCTS algorithm. In the second
part, we compare our approach with other from the state-of-the-art.

###  5.1 Experimental Setup

The complete grammar we use in our experiments is based on the same set of
components as in AutoSklearn. We employ the same hyperparameter ranges and
functions as AutoSklearn when sampling and incorporates them into the grammar.
For our experiments, we sample 3 values of each hyperparameter that, added to
the pipeline options, extend the space to more than 183 billions possible
combinations. This has to be compared with the $\sim$24K configurations of the
original GramML method [29].

For evaluation, we use the OpenML-CC18 benchmark [27], a collection of 72
binary and multi-class classification datasets carefully curated from the
thousands of datasets on OpenML [3]. We use standard train-test splits, as
provided by the benchmark suite. Also, to calculate validation scores (e.g.
rewards), we use AutoSklearn’s default strategy which consists on using 67% of
the data for training and 33% for validation. All experiments were performed
on an Amazon EC2 R5 spot instance (8 vCPU and 64GB of RAM).

###  5.2 Ablation Study

To conduct the ablation study, four tasks with low computational cost were
selected from the OpenML-CC18 benchmark222These tasks, identified by task IDs
11, 49, 146819, and 10093.. For each task, we run 100 iterations of the
algorithm for each policy configuration and report the mean and standard
deviation for the following metrics: time per iteration (Time Iter), number of
actions per iteration (Act/Iter), simulations’ repetition rate (Rep Ratio),
time spent by the first iteration (Time 1st), the number of actions for the
first iteration (1st Act), and the total time (Tot Time) and total number of
actions (Tot Act) at the end of the 100 iterations. Time measurements only
consider the time of the algorithm and do not take into account the fitting
time. The simulations’ repetition ratio, _i.e_. the number of times the
algorithm repeats the same paths during simulation, can be seen as a measure
of the exploration efficiency of the algorithm (the lower this value, the more
efficient the exploration of the search space).

Tables 1, 2 and 3 show results for the UTC, BTS and TPE strategies,
respectively. All three methods depend on a single hyperparameter to control
the overall behavior of the search algorithm, _i.e_. the exploration-
exploitation trade-off. An immediate effect of varying such parameters can be
observed on the time it takes for the algorithm to complete an iteration. This
has to do with the iteration ending when the algorithm reaches a leaf. If we
favor exploration (i.e. increase C and J for UCT and BTS, or decrease $\gamma$
for TPE), the algorithm explores the tree in width, performs more total
actions, and takes longer to reach a leaf. On the other hand, by favoring
exploitation, the algorithm performs fewer actions and reaches the terminal
nodes faster. The number of actions per iteration and the time per iteration
are strongly correlated.

Comparing the different strategies, we observe that UCT with $C=0.7$ is the
most efficient alternative (in terms of the repetition ratio), while UCT with
$C=0$ is the fastest, with the later being also the less efficient in terms of
exploration.

Finally, we consider the best version of each strategy prioritizing
exploration efficiency (low iteration rate) trying to keep the time per
iteration limited to $\sim$0.1. We select UCT with $C=0.7$, BTS with $J=1$,
and TPE with $\gamma=0.85$ and name them $\text{GramML}^{++}_{\text{UCT}}$,
$\text{GramML}^{++}_{\text{BTS}}$, and $\text{GramML}^{++}_{\text{TPE}}$,
respectively.

Table 1: Ablation results for the UCT strategy for different values of the parameter $C$. | C=0 | C=0.1 | C=0.7 | C=1
---|---|---|---|---
Time Iter | 0,04 (0,01) | 0,05 (0,01) | 0,10 (0,02) | 0,14 (0,04)
Time 1st | 0,76 (0,31) | 1,65 (0,49) | 7,16 (2,50) | 10,64 (3,99)
Tot Time | 3,94 (0,64) | 4,90 (0,54) | 9,71 (2,14) | 14,21 (3,84)
Act/Iter | 5,8 (0,3) | 7,1 (0,9) | 19,4 (4,7) | 28,1 (7,7)
1st Act | 165,5 (52,7) | 324,2 (95,9) | 1535,0 (545,8) | 2201,2 (810,5)
Tot Act | 682,5 (35,4) | 817,2 (95,7) | 2047,0 (475,9) | 2914,0 (776,5)
Rep Ratio | 0,65 (0,04) | 0,58 (0,07) | 0,40 (0,08) | 0,45 (0,10)
Table 2: Ablation results for the BTS strategy for different values of the parameter $J$. | J=1 | J=10 | J=100 | J=1000
---|---|---|---|---
Time Iter | 0,15 (0,05) | 0,15 (0,06) | 0,33 (0,23) | 0,28 (0,12)
Time 1st | 7,13 (2,19) | 7,42 (3,86) | 18,90 (15,42) | 14,56 (7,05)
Tot Time | 14,72 (4,51) | 14,98 (6,27) | 33,21 (23,20) | 28,18 (12,40)
Act/Iter | 24,6 (7,8) | 24,3 (11,5) | 46,0 (30,0) | 36,6 (15,4)
1st Act | 1307,5 (398,8) | 1333,2 (713,6) | 2744,0 (2075,2) | 2027,0 (927,2)
Tot Act | 2561,2 (783,8) | 2531,5 (1151,3) | 4706,7 (3008,7) | 3760,2 (1546,0)
Rep Ratio | 0,50 (0,13) | 0,51 (0,10) | 0,54 (0,09) | 0,59 (0,04)
Table 3: Ablation results for the TPE strategy for different values of the parameter $\gamma$. |  $\gamma$=50 |  $\gamma$=65 |  $\gamma$=75 |  $\gamma$=85
---|---|---|---|---
Time Iter | 0,33 (0,16) | 0,19 (0,03) | 0,18 (0,06) | 0,12 (0,02)
Time 1st | 11,16 (7,55) | 5,93 (1,47) | 4,02 (2,50) | 1,71 (0,57)
Tot Time | 33,31 (15,82) | 18,92 (3,35) | 17,86 (6,05) | 12,33 (1,76)
Act/Iter | 15,0 (5,5) | 9,1 (1,3) | 7,8 (1,8) | 5,9 (0,2)
1st Act | 740,2 (382,0) | 434,7 (81,8) | 286,5 (128,4) | 143,0 (25,7)
Tot Act | 1605,5 (555,8) | 1018,7 (129,8) | 882,5 (187,2) | 697,0 (23,0)
Rep Ratio | 0,47 (0,06) | 0,59 (0,06) | 0,60 (0,05) | 0,68 (0,01)

###  5.3 Comparison with other techniques

In this section, we compare our methods with others from the literature. For
each task, each method is run for an hour (time budget). We report performance
on the test set proposed in the benchmark suite. We compare the different
variations of our approach to AutoSklearn [13] and MOSAIC [22] as they allow
us to rely on the exact same set of basic ML components. In all cases, we set
a maximum fitting time of 300 seconds, as the time limit for a single call to
the machine learning model. If the algorithm runs beyond this time, the
fitting is terminated. In our case, we return a reward of zero for methods
exceeding this time limit. For each task, we ranked the performance of all
systems and reported the average ranking (lower is better) and average
performance score (higher is better) calculated as 1 - regret (the difference
between maximum testing performance found so far by each method and the true
maximum) at each time step.

Comparative results are shown in Figure 1. Figure 1(a) shows the average rank
across time. A marked difference of $\text{GramML}^{++}_{\text{BTS}}$ variant
can be seen over MOSAIC, AutoSklearn and GramML. In particular,
$\text{GramML}^{++}_{\text{BTS}}$ have the best ranking over time, although it
is surpassed in the first few minutes by MOSAIC and GramML. There also seems
to be a slight decreasing trend for $\text{GramML}^{++}_{\text{BTS}}$ which
would seem to indicate that the results could improve with more time.

In addition, Figure 1(b) shows the average score across time.
$\text{GramML}^{++}$ variants have a visible difference from MOSAIC,
AutoSKLearn and GramML. The best average score is also achieved by
$\text{GramML}^{++}_{\text{BTS}}$. An interesting observation is that the
previous version of GramML was eventually slightly outperformed in terms of
average score by AutoSKLearn. One reason for this was that GramML did not
utilize hyperparameters. Conversely, the variants of $\text{GramML}^{++}$
demonstrate a significant improvement over their predecessor and AutoSKLearn.
It should be noted that the average regret measure may not be very
representative if the accuracy varies greatly between tasks. However, we
believe that the graph provides important information together with the
average ranking for the comparison between techniques in the OpenML-CC18
benchmark.

(a) Average rank over 1 hour execution.

(b) Average score over 1 hour execution.

Figure 1: Results of the approaches on OpenML-CC18 benchmark.

In order to check the presence of statistically significant differences in the
average rank distributions, Figure 2 shows the results using critical
difference (CD) diagrams [7] at 15, 30, 45 and 60 minutes. We use a non-
parametric Friedman test at p $<$ 0.05 and a Nemenyi post-hoc test to find
which pairs differ [14]. If the difference is significant, the vertical lines
for the different approaches appear as well separated. If the difference is
not significant, they are joined by a thick horizontal line.

The diagrams show statistically significant differences at 15 minutes between
AutoSKLearn and the other variants. At 30 minutes both
$\text{GramML}^{++}_{\text{BTS}}$ and GramML have significant difference with
MOSAIC and AutoSKLearn. At minute 45, the difference between
$\text{GramML}^{++}_{\text{BTS}}$ and the other techniques remains, but there
is no significant difference between the previous version of GramML and
$\text{GramML}^{++}_{\text{TPE}}$. Finally, at 60 minutes,
$\text{GramML}^{++}{\text{BTS}}$ showed a significant difference compared to
MOSAIC and AutoSKLearn, while $\text{GramML}^{++}{\text{TPE}}$ showed a
significant difference only compared to MOSAIC. However,
$\text{GramML}^{++}_{\text{UCT}}$ did not show a statistically significant
difference compared to either MOSAIC or AutoSKLearn.

(a) 15 minutes

(b) 30 minutes

(c) 45 minutes

(d) 1 hour

Figure 2: CD plots with Nimenyi post-hoc test

##  6 Conclusions and Future Work

This article presents a model-free approach for grammar-based AutoML that
integrates hyperparameter search as an extension to GramML [29]. The approach
involves two steps: (1) incorporating hyperparameter values into the grammar
using grammar rules and (2) modifying the Monte Carlo Tree Search (MCTS)
algorithm to improve the collection of the best pipeline configuration,
implement pruning, and support different selection policies and
backpropagation functions. The evaluation involved an ablation study to assess
the efficiency of each selection policy and a comparison with state-of-the-art
techniques, which showed significant improvements, particularly with a variant
that uses bootstrapped Thompson Sampling [10]. Consecuently, this work
demonstrates the effectiveness of incorporating hyperparameter search into
grammar-based AutoML, and provides a promising approach for addressing the
challenges of larger search spaces.

Furthermore, this work opens up several avenues for future research. In model-
based AutoML, meta-learning can improve exploration efficiency [26]. However,
the application of meta-learning in model-free AutoML based on grammars
remains an open question. Additionally, incorporating resource information
into the AutoML objective function is also an important direction to follow
[28]. Finally, the speed of the search algorithms is an important factor that
can be improved through horizontal scalability [4]. Therefore, future research
should aim to parallelize the algorithm for more efficient resource
utilization.

## References

  * [1] Bai, A., Wu, F., Chen, X.: Bayesian mixture modelling and inference based thompson sampling in monte-carlo tree search. Advances in neural information processing systems 26 (2013)
  * [2] Bergstra, J., Komer, B., Eliasmith, C., Yamins, D., Cox, D.D.: Hyperopt: a python library for model selection and hyperparameter optimization. Computational Science & Discovery 8(1), 014008 (2015)
  * [3] Bischl, B., Casalicchio, G., Feurer, M., Hutter, F., Lang, M., Mantovani, R.G., van Rijn, J.N., Vanschoren, J.: Openml benchmarking suites. arXiv:1708.03731v2 [stat.ML] (2019)
  * [4] Bourki, A., Chaslot, G., Coulm, M., Danjean, V., Doghmen, H., Hoock, J.B., Hérault, T., Rimmel, A., Teytaud, F., Teytaud, O., et al.: Scalability and parallelization of monte-carlo tree search. In: Computers and Games: 7th International Conference, CG 2010, Kanazawa, Japan, September 24-26, 2010, Revised Selected Papers 7. pp. 48–58. Springer (2011)
  * [5] Chapelle, O., Li, L.: An empirical evaluation of thompson sampling. Advances in neural information processing systems 24 (2011)
  * [6] Chomsky, N.: Syntactic structures. In: Syntactic Structures. De Gruyter Mouton (2009)
  * [7] Demšar, J.: Statistical comparisons of classifiers over multiple data sets. The Journal of Machine learning research 7, 1–30 (2006)
  * [8] Drori, I., Krishnamurthy, Y., Lourenco, R., Rampin, R., Cho, K., Silva, C., Freire, J.: Automatic machine learning by pipeline synthesis using model-based reinforcement learning and a grammar. arXiv preprint arXiv:1905.10345 (2019)
  * [9] Drori, I., Krishnamurthy, Y., Rampin, R., Lourenco, R.d.P., Ono, J.P., Cho, K., Silva, C., Freire, J.: Alphad3m: Machine learning pipeline synthesis. arXiv preprint arXiv:2111.02508 (2021)
  * [10] Eckles, D., Kaptein, M.: Thompson sampling with the online bootstrap. arXiv preprint arXiv:1410.4009 (2014)
  * [11] Efron, B.: Bayesian inference and the parametric bootstrap. The annals of applied statistics 6(4), 1971 (2012)
  * [12] Evans, B., Xue, B., Zhang, M.: An adaptive and near parameter-free evolutionary computation approach towards true automation in automl. In: 2020 IEEE Congress on Evolutionary Computation (CEC). pp. 1–8. IEEE (2020)
  * [13] Feurer, M., Klein, A., Eggensperger, K., Springenberg, J., Blum, M., Hutter, F.: Efficient and robust automated machine learning. Advances in neural information processing systems 28 (2015)
  * [14] Gijsbers, P., Bueno, M.L., Coors, S., LeDell, E., Poirier, S., Thomas, J., Bischl, B., Vanschoren, J.: Amlb: an automl benchmark. arXiv preprint arXiv:2207.12560 (2022)
  * [15] Hayes, C.F., Reymond, M., Roijers, D.M., Howley, E., Mannion, P.: Distributional monte carlo tree search for risk-aware and multi-objective reinforcement learning. In: Proceedings of the 20th international conference on autonomous agents and multiagent systems. pp. 1530–1532 (2021)
  * [16] Hutter, F., Hoos, H.H., Leyton-Brown, K.: Sequential model-based optimization for general algorithm configuration. In: International conference on learning and intelligent optimization. pp. 507–523. Springer (2011)
  * [17] Katz, M., Ram, P., Sohrabi, S., Udrea, O.: Exploring context-free languages via planning: The case for automating machine learning. In: Proceedings of the International Conference on Automated Planning and Scheduling. vol. 30, pp. 403–411 (2020)
  * [18] Kocsis, L., Szepesvári, C.: Bandit based monte-carlo planning. In: European conference on machine learning. pp. 282–293. Springer (2006)
  * [19] Marinescu, R., Kishimoto, A., Ram, P., Rawat, A., Wistuba, M., Palmes, P.P., Botea, A.: Searching for machine learning pipelines using a context-free grammar. In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 35, pp. 8902–8911 (2021)
  * [20] Mohr, F., Wever, M., Hüllermeier, E.: Ml-plan: Automated machine learning via hierarchical planning. Machine Learning 107(8), 1495–1515 (2018)
  * [21] Olson, R.S., Moore, J.H.: TPOT: A tree-based pipeline optimization tool for automating machine learning. In: Workshop on automatic machine learning. pp. 66–74. PMLR (2016)
  * [22] Rakotoarison, H., Sebag, M.: Automl with monte carlo tree search. In: Workshop AutoML 2018@ ICML/IJCAI-ECAI (2018)
  * [23] de Sá, A.G., Pinto, W.J.G., Oliveira, L.O.V., Pappa, G.L.: Recipe: a grammar-based framework for automatically evolving classification pipelines. In: European Conference on Genetic Programming. pp. 246–261. Springer (2017)
  * [24] Segovia-Aguas, J., Jiménez, S., Jonsson, A.: Generating context-free grammars using classical planning. In: Proceedings of the Twenty-Sixth International Joint Conference on Artificial Intelligence (IJCAI-17); 2017 Aug 19-25; Melbourne, Australia. Melbourne: IJCAI; 2017. p. 4391-7. IJCAI (2017)
  * [25] Thornton, C., Hutter, F., Hoos, H.H., Leyton-Brown, K.: Auto-weka: Combined selection and hyperparameter optimization of classification algorithms. In: Proceedings of the 19th ACM SIGKDD international conference on Knowledge discovery and data mining. pp. 847–855 (2013)
  * [26] Vanschoren, J.: Meta-learning. Automated machine learning: methods, systems, challenges pp. 35–61 (2019)
  * [27] Vanschoren, J., Van Rijn, J.N., Bischl, B., Torgo, L.: Openml: networked science in machine learning. ACM SIGKDD Explorations Newsletter 15(2), 49–60 (2014)
  * [28] Vazquez, H.C.: A general recipe for automated machine learning in practice. In: Advances in Artificial Intelligence–IBERAMIA 2022: 17th Ibero-American Conference on AI, Cartagena de Indias, Colombia, November 23–25, 2022, Proceedings. pp. 243–254. Springer (2023)
  * [29] Vazquez, H.C., Sánchez, J., Carrascosa, R.: GramML: Exploring context-free grammars with model-free reinforcement learning. In: Sixth Workshop on Meta-Learning at the Conference on Neural Information Processing Systems (2022), https://openreview.net/forum?id=OpdayUqlTG
  * [30] Waddle, V.E.: Production trees: A compact representation of parsed programs. ACM Transactions on Programming Languages and Systems (TOPLAS) 12(1), 61–83 (1990)


# Discussion of the manuscript:
Spatial+: a novel approach to spatial confounding

Georgia Papadogeorgou
Department of Statistics, University of Florida

I congratulate Dupont, Wood and Augustin (DWA hereon) for providing an easy-
to-implement method for estimation in the presence of spatial confounding, and
for addressing some of the complicated aspects on the topic. The method
regresses the covariate of interest on spatial basis functions and uses the
residuals of this model in an outcome regression. The authors show that, if
the covariate is not completely spatial, this approach leads to consistent
estimation of the conditional association between the exposure and the
outcome.

Below I discuss conceptual and operational issues that are fundamental to
inference in spatial settings: (i) the target quantity and its
interpretability, (ii) the non-spatial aspect of covariates and their relative
spatial scales, and (iii) the impact of spatial smoothing. While DWA provide
some insights on these issues, I believe that the audience might benefit from
a deeper discussion. In what follows, I focus on the setting where a
researcher is interested in interpreting the relationship between a given
covariate and an outcome. I refer to the covariate of interest as the exposure
to differentiate it from the rest.

##  1 What are we estimating?

First, I discuss the target of estimation in the manuscript and in the
literature. This target is key in deciding what is the appropriate method for
estimation.

I begin by introducing some notation. Assume that $Z$ and $Y$ are the exposure
and outcome of interest, and covariates $\bm{X}$ include all the variables
that confound this relationship. I use the subscript $i$ to denote the value
of the variables at location $i=1,2,\dots,n$. In what follows, I assume that
$\bm{X}$ is comprised of four random variables:  (a) a measured variable $C$
that can be spatial or not, and unmeasured variables of the following three
types:  (b) $S_{1}$that is completely spatial,  (c) $S_{2}^{+}$that includes a
spatial and an independent component, and  (d) $U$that is not spatial.
Similarly to DWA, we can think of $S_{2}^{+}$ as a random variable that can be
written as $S_{2}^{+}=S_{2}+E$ for $S_{2}$ being a completely spatial random
variable and $E$ representing a component without spatial structure.

For simplicity of discussion, I also assume that the following additive model
holds:

$Y_{i}=\beta_{0}+\beta_{1}Z_{i}+\beta_{2}C_{i}+\beta_{3}S_{i1}+\beta_{4}S_{i2}^{+}+\beta_{5}U_{i}+\epsilon_{i},$
(1)

where $\epsilon_{i}\overset{iid}{\sim}N(0,\sigma^{2})$. This model implies
that any spatial structure in the outcome $Y$ is due to spatial dependencies
in its predictors, though this assumption could be relaxed.

###  1.1 The “unconditional” association

A researcher only has access to $(C_{i},Z_{i},Y_{i})$, for $i=1,2,\dots,n$,
while the remaining three random variables in $\bm{X}$ are not measured. With
these data available, a first-step analysis could regress the outcome, $Y$ on
the exposure and the sole observed covariate, $(Z,C)$. Even if not directly
stated, the target of estimation in this regression is the relationship
between the exposure and outcome of interest, while conditioning solely on the
measured variable:

$\displaystyle(Y,Z)\mid C.$ (Unconditional target)

Even though the remaining variables in $\bm{X}$ are not accounted for, this
relationship can be of explicit scientific interest, especially if the
measured variable $C$ is an interpretable characteristic of the population at
hand. For example, $Z,$ $Y$ and $C$ can represent a subject’s exposure to a
given pollutant, their blood pressure levels, and age, respectively. In this
case, the unconditional target represents the expected change in an
individual’s blood pressure for a unit increase in the pollutant, after
adjusting for age. I refer to this target as the unconditional association as
it does not condition on the spatial random variables. For this target of
estimation, the OLS estimator will be unbiased.

However, not conditioning on the remaining variables in $\bm{X}$ leads to
complications for inference for the estimated coefficient. If at least one of
the coefficients of the spatial predictors are non-zero ($\beta_{3}\neq 0$, or
$\beta_{4}\neq 0$), the residuals of this regression would be spatially
correlated. Since correlation in the model residuals that is not accounted for
would lead to incorrect inference on model coefficients, researchers often
include a spatial random effect with covariance structure that aims to mimic
the residual dependence. The goal of the random effect is to explain the
spatial dependence in the outcome that is not already explained by $(Z,C)$.
However, if the exposure $Z$ is, at least partially, spatial, and since the
additional model components are themselves spatial, they are often spuriously
associated due to structural dependencies (Hodges & Reich, 2010; Lee & Ogburn,
2020). Therefore, the random effect “involuntarily” captures part of the
spatial variability in the outcome that was originally explained by the
exposure, and the coefficient of the exposure in the model with the spatial
random effect no longer has the same interpretation as the coefficient of the
non-spatial model.

The restricted spatial regression (RSR) approach ensures that the random
effect is orthogonal to the exposure of interest. That means that the target
of estimation for RSR is the same as the one from the non-spatial model, and
RSR is truly useful, applicable, and unbiased for the relationship between the
exposure and the outcome conditional only on the measured variables.

###  1.2 The “spatially-conditional” association

The target of estimation for DWA is different. Instead of focusing on the
$(Y,Z)$ relationship when conditioning solely on $C$, DWA aim to also adjust
for all spatial confounders. Informally, the target is set to learning

$\displaystyle(Y,Z)\mid C,\text{ all spatial confounders},$ (Spatially-
conditional target)

though as we will see below, it truely achieves an alternative quantity.
Therefore, in the presence of unmeasured spatial confounders, the targets of
estimation for RSR and Spatial+ are different, the coefficient of the exposure
in the two models has a different interpretation, and the methods should not
be directly compared.

Since confounders are predictors of the outcome, an outcome model that does
not adjust for spatial confounders will exhibit residual spatial dependence.
At the same time, since they are predictors of the exposure, adding components
in the regression model that are orthogonal to the exposure (as in RSR) will
not capture them. However, since conditioning on any spatial variable that is
not a confounder does not change the relationship under study, the spatially-
conditional target above can be re-written as

$(Y,Z)\mid C,\text{ all spatial variables}.$

Writing the target of estimation explicitly as conditioning on any spatial
variable clarifies two aspects of the manuscript.  (1) If the exposure $Z$ is
fully spatial and it does not include any independent components, then the
spatially-conditional relationship cannot be learnt without additional
assumptions, since the exposure is perfectly collinear with the conditioning
set. 111To make progress in this setting, one might need to make assumptions
on the relative spatial scales of the spatial confounders and the exposure,
though I ignore this scenario here. This setting is investigated in Schnell &
Papadogeorgou (2020), and a related discussion in terms of the spatial scales
of the spatial confounders themselves can be found in Section 1.3. (2) If $Z$
includes an independent component, the spatially-conditional relationship
between the exposure and the outcome can be learnt while focusing on the non-
spatial aspect of the data.  The former is central in DWA’s manuscript. The
latter motivates three of the approaches discussed in the manuscript: the
spatial model that includes basis functions directly in the outcome model, the
gSEM model that controls for space both in the exposure and in the outcome
model, and the proposed Spatial+ model that regresses the exposure itself on
the basis functions and uses the exposure model residuals in the outcome
model. When no spatial smoothing is performed, all of these methods adjust for
any completely spatial structure in the data, and are therefore unbiased for
the spatially-conditional association, as expected, and as illustrated in the
simulations. The issues that arise with smoothing are discussed in Section 2.

The latter observation is also key in Schnell & Papadogeorgou (2020), where
the authors link this observation to the assumption of positivity in causal
inference, and discuss that estimation of the exposure coefficient can only be
accurately achieved if the exposure itself varies within levels of the spatial
confounder. In the setting of DWA where the spatial confounder can vary at any
spatial frequency, this observation by Schnell & Papadogeorgou (2020) is
translated to $Z$ including an independent component.

As discussed in Section 1.1, RSR studies the exposure-outcome relationship
when conditioning only on the measured variables. As discussed above, Spatial+
aims to extend the conditioning set to include spatial confounders. Even
though conditioning on a more complete set of confounders is generally
preferred, I illustrate that the estimated quantities from Spatial+ are not
always interpretable. Take for example the true outcome model in Equation 1,
and assume that both spatial variables are important $(\beta_{3}\neq
0,\beta_{4}\neq 0)$ while the independent variable is not ($\beta_{5}=0$).
Then, we can re-write this model as

$Y_{i}=\beta_{0}+\beta_{1}Z_{i}+\beta_{2}C_{i}+\beta_{3}S_{i1}+\beta_{4}(S_{i2}+E_{i})+\epsilon_{i},$

where $S_{2}$ and $E$ decompose $S_{2}^{+}$ in its spatial and independent
components. The target of estimation that conditions on measured variables and
all spatial confounders would be

$(Y,Z)\mid C,S_{1},S_{2}^{+}\quad\equiv\quad(Y,Z)\mid C,S_{1},S_{2},E.$

Therefore, this target would represent the exposure-outcome relationship when
conditioning on the variables $C,S_{1},S_{2}^{+}$. However, the $E$ part of
$S_{2}^{+}$ is not spatial and cannot be captured by any spatial basis
function. Therefore, even though the implicit target would condition on all
spatial variables, the method instead estimates

$(Y,Z)\mid C,S_{1},S_{2}.$ (Spatially-conditional quantity achieved)

Since the achieved target conditions solely on part of a confounder, the
estimated conditional exposure-outcome relationship has a complicated
interpretation. Whether this relationship is more or less interpretable than
the unconditional target of Section 1.1 is, I believe, up for debate!

In the previous blood pressure example, let’s think of $S_{1}$ as representing
outdoor temperature which is completely spatial, and $S_{2}^{+}$ as
representing eating habits (such as intake of salt, or fried food). Even
though eating habits are partially spatial due to local culture, they also
include an independent, subject-specific component. Then, the coefficient
estimated by Spatial+ would describe the expected change in blood pressure for
a unit increase in the pollutant when adjusting for age, temperature, and part
of someone’s eating habits.

###  1.3 Assumptions for interpretable spatially-conditional target

Therefore, the achieved target of estimation for the Spatial+ model often has
a complicated interpretation, largely due to variables that include both
spatial and independent components. So when can Spatial+ be used to acquire
interpretable conditional associations?

The first case corresponds to the scenario where there are no spatial
confounders of the exposure-outcome relationship that include an independent
component (Papadogeorgou et al., 2019). This assumption is implicitly made in
equation (1) of DWA. If all covariates in the analysis are either fully
spatial or not spatial at all (if, for example, $S_{2}^{+}=S_{2}$), then
Spatial+ estimates the exposure-outcome relationship conditional on all
spatial variables, which has a clear interpretation. In turn, if there exist
covariates that include both spatial and independent components such as
$S_{2}^{+}$, these variables must not confound the exposure-outcome
relationship when conditioning on fully spatial variables. This can be
satisfied if

$E\perp\\!\\!\\!\perp Y\mid Z,C,S_{1},S_{2},\quad\text{or}\quad
E\perp\\!\\!\\!\perp Z\mid C,S_{1},S_{2}$

where $\perp\\!\\!\\!\perp$ is used to denote independence. However, since all
of these spatial variables are unmeasured, assuming that none of the spatial
confounders includes an independent component remains, exactly that, an
assumption.

So what happens when we are not willing to make the assumption that all
spatial confounders are fully spatial? Progress can still be made in this
setting if we are willing to make an alternative assumption on the relative
spatial scales of spatial confounders, and update Spatial+ accordingly.
Remember that the issue arises because $S_{1}$ and $S_{2}^{+}$ are both
spatial and unmeasured, which forbids us from directly differentiating between
the two. To make progress, we could make an assumption which allows us to
separate the one spatial variable from the other. Such an assumption could be
the following: $S_{2}^{+}$ varies only in higher spatial frequencies than
$S_{1}$. If this assumption holds, then including spatial basis functions that
represent only low-frequency spatial components would capture $S_{1}$ fully,
while not capturing any aspect (spatial or independent) of $S_{2}^{+}$
(Antonelli et al., 2017; Keller & Szpiro, 2020). This could be acheived, for
example, using wavelet or fourier basis functions for which the frequency of
the basis functions used is directly controllable. Under this assumption and
the corresponding adjustment of Spatial+, the method would be estimating the
exposure-outcome relationship when conditioning on $C$ and $S_{1}$ only, which
would have a clear interpretation.

###  1.4 Spatial confounding when estimating causal effects

If the set of covariates $\bm{X}$ includes all the confounders of the
exposure-outcome relationship, then an outcome model which adjusts for the
exposure and all of the confounders and is correctly specified can be the base
for estimating causal effects (Robins, 1986; Snowden et al., 2011,
g-computation). I focus first on linear additive models such as the one in
Equation 1. In this case, the coefficient of the exposure, $\beta_{1}$, can be
interpreted causally as the effect of increasing the exposure $Z$ by one unit
on the expected outcome $Y$. However, the coefficient will not have the same,
causal interpretation if we do not adjust for the variables $S_{1},$ $S_{2}+$
and $U$, which are unmeasured, nor if we only adjust for a subset of them.
Since $U$ is a completely independent random variable, we cannot hope to
identify the causal effect of the exposure on the outcome in this manner.
What’s more, adjusting for some of the confounders, and not all, does not
necessarily mean that the bias for estimating causal effects is reduced.
Instead, it could increase, if confounding from spatial and independent
variables are in opposite directions! Therefore, if unmeasured independent
variables confound the exposure-outcome relationship, adjusting for spatial
unmeasured confounders does not necessarily bring us “closer” to estimating
causal effects.

However, if the true outcome model is non-linear, or includes exposure-
covariate interactions, the complications are not resolved, even if we assume
the absence of independent, unmeasured confounders (for example, if $U$ does
not exist, and $S_{2}+$ is fully spatial). In order to use the outcome model
for estimating causal effects, g-computation teaches us that we need to
integrate over the distribution of the confounders, and the coefficient of the
exposure is not necessarily causally interpretable (though log-linear models
are an exception for a specific type of causal contrast (Schnell &
Papadogeorgou, 2020)). When the spatial confounders are unmeasured, and their
distribution cannot be learnt from data, integrating over their distribution
cannot be achieved directly, and requires additional assumptions (see, for
example Christiansen et al., 2020). Recently, Reich et al. (2021) wrote a
review article regarding spatial causal inference, including a thorough
discussion on spatial confounding.

##  2 Spatial smoothing

Based on the discussion of Section 1, we have established that the non-spatial
and RSR models estimate a different quantity compared to the rest. We have
also established that regressing away all spatial aspect of the outcome or the
exposure (through, for example, unpenalized splines) would adjust fully for
the spatial component of all spatial confounders. Assuming that this quantity
has a useful interpretation, here I focus on investigating the following
questions: what happens with spatial smoothing, and how much should predictive
accuracy matter?

Since spatial smoothing excludes the higher frequency components, it implies
that the achieved target represents the exposure-outcome association
conditional on the spatial variables that vary at the included, lower
frequencies only. If all spatial confounders vary at lower frequencies, the
estimated quantity will truly align with the spatially-conditional
association. However, even though spatial smoothing is expected to lead to
lower prediction MSE, it is not immidiately clear whether it should be used
for coefficient estimation. In the presence of high frequency spatial
confounders, spatial smoothing might lead us to “miss” these confounders.

I investigate two scenarios where the Spatial and Spatial+ models might miss
strong, high frequency spatial confounders, and return biased estimates of the
spatially-conditional association. 222Remember that the Spatial model employs
spatial basis functions in the outcome model. These spatial variables can be
strong confounders in at least two different ways:  (1) they can be important
predictors of the exposure and weak predictors of the outcome, or  (2) they
can be weak predictors of the exposure, and important predictors of the
outcome.  In the former case, since the high-frequency spatial variables are
strong predictors of the exposure, one hopes that the choice of penalty in the
exposure model for Spatial+ will not penalize these high frequency components
heavily. The reverse is true for the Spatial model, for which we would expect
that it will penalize the high frequency spatial basis functions more heavily,
as these are weak predictors of the outcome. Therefore, in this scenario, we
would expect Spatial+ to perform better than the Spatial model for estimating
the spatially-conditional association. The reverse is true in the latter case.
If the spatial confounders are weak predictors of the exposure, it is expected
that the chosen penalty for the exposure model according to Spatial+ will
penalize the high frequency components heavily, returning biased estimates of
the spatially-conditional association. In contrast, we would expect the
Spatial model that includes spatial basis functions in the outcome model
directly to perform more accurately.

Even though it is prudent not to overly generalize the message from this
discussion, I believe that deeper investigation is warranted on how the
methods perform under different confounding scenarios. That is essential under
the light of similar observations from other parts of the literature that
focus on estimation rather than prediction, and encourage the use of both an
exposure and an outcome model simultaneously (e.g., Wang et al., 2012; Belloni
et al., 2014; Shortreed & Ertefaie, 2017; Antonelli et al., 2019).
Specificially in the setting of spatial confounding, the use of both models is
proposed in at least two approaches (gSEM, investigated in the manuscript;
Thaden & Kneib (2018) and Schnell & Papadogeorgou (2020)). If we view gSEM as
a combination of the Spatial and Spatial+ models in that it adjusts for
spatial effects in both the exposure and the outcome model, one might expect
gSEM to be more robust to high frequency spatial confounding than any of the
two methods separately.

Finally, the work referenced above has been crucial in illuminating that
prediction and estimation are two separate goals, and that methodology that is
appropriate for one can be inappropriate for the other. In light of this work,
even though the methods’ predictive accuracy can inform us of how well each
model fits the data, it should not be used as a model-choice criterion when we
are interested in interpreting estimated parameters. In fact, a model with
lower AIC might often have higher bias for the coefficient of interest.

## References

  * (1)
  * Antonelli et al. (2019) Antonelli, J., Parmigiani, G. & Dominici, F. (2019), ‘High-dimensional confounding adjustment using continuous spike and slab priors’, Bayesian analysis 14(3), 805.
  * Antonelli et al. (2017) Antonelli, J., Schwartz, J., Kloog, I. & Coull, B. A. (2017), ‘Spatial multiresolution analysis of the effect of pm2. 5 on birth weights’, The annals of applied statistics 11(2), 792.
  * Belloni et al. (2014) Belloni, A., Chernozhukov, V. & Hansen, C. (2014), ‘Inference on treatment effects after selection among high-dimensional controls’, The Review of Economic Studies 81(2), 608–650.
  * Christiansen et al. (2020) Christiansen, R., Baumann, M., Kuemmerle, T., Mahecha, M. D. & Peters, J. (2020), ‘Towards causal inference for spatio-temporal data: Conflict and forest loss in colombia’, arXiv preprint arXiv:2005.08639 .
  * Hodges & Reich (2010) Hodges, J. S. & Reich, B. J. (2010), ‘Adding spatially-correlated errors can mess up the fixed effect you love’, The American Statistician 64(4), 325–334.
  * Keller & Szpiro (2020) Keller, J. P. & Szpiro, A. A. (2020), ‘Selecting a scale for spatial confounding adjustment’, Journal of the Royal Statistical Society. Series A,(Statistics in Society) 183(3), 1121.
  * Lee & Ogburn (2020) Lee, Y. & Ogburn, E. L. (2020), ‘Network dependence can lead to spurious associations and invalid inference’, Journal of the American Statistical Association pp. 1–15.
  * Papadogeorgou et al. (2019) Papadogeorgou, G., Choirat, C. & Zigler, C. M. (2019), ‘Adjusting for unmeasured spatial confounding with distance adjusted propensity score matching’, Biostatistics 20(2), 256–272.
  * Reich et al. (2021) Reich, B. J., Yang, S., Guan, Y., Giffin, A. B., Miller, M. J. & Rappold, A. (2021), ‘A review of spatial causal inference methods for environmental and epidemiological applications’, International Statistical Review .
  * Robins (1986) Robins, J. (1986), ‘A new approach to causal inference in mortality studies with a sustained exposure period—application to control of the healthy worker survivor effect’, Mathematical modelling 7(9-12), 1393–1512.
  * Schnell & Papadogeorgou (2020) Schnell, P. M. & Papadogeorgou, G. (2020), ‘Mitigating unobserved spatial confounding when estimating the effect of supermarket access on cardiovascular disease deaths’, The Annals of Applied Statistics 14(4), 2069–2095.
  * Shortreed & Ertefaie (2017) Shortreed, S. M. & Ertefaie, A. (2017), ‘Outcome-adaptive lasso: variable selection for causal inference’, Biometrics 73(4), 1111–1122.
  * Snowden et al. (2011) Snowden, J. M., Rose, S. & Mortimer, K. M. (2011), ‘Implementation of g-computation on a simulated data set: demonstration of a causal inference technique’, American journal of epidemiology 173(7), 731–738.
  * Thaden & Kneib (2018) Thaden, H. & Kneib, T. (2018), ‘Structural equation models for dealing with spatial confounding’, The American Statistician 72(3), 239–252.
  * Wang et al. (2012) Wang, C., Parmigiani, G. & Dominici, F. (2012), ‘Bayesian effect estimation accounting for adjustment uncertainty’, Biometrics 68(3), 661–671.