# Relative Feature Importance
††thanks: This work is funded by the German Federal Ministry of Education and
Research (BMBF) under Grant No. 01IS18036A and supported by the Bavarian State
Ministry of Science and the Arts in the framework of the Centre
Digitisation.Bavaria (ZD.B). The authors of this work take full responsibility
for its content.

Gunnar König1,2, Christoph Molnar1, Bernd Bischl1, Moritz Grosse-Wentrup2,3,4
1Institute for Statistics, LMU Munich, 2Research Group Neuroinformatics,
University of Vienna,
3Research Platform Data Science @ Uni Vienna, 4Vienna Cognitive Science Hub

###### Abstract

Interpretable Machine Learning (IML) methods are used to gain insight into the
relevance of a feature of interest for the performance of a model. Commonly
used IML methods differ in whether they consider features of interest in
isolation, e.g., Permutation Feature Importance (PFI), or in relation to all
remaining feature variables, e.g., Conditional Feature Importance (CFI). As
such, the perturbation mechanisms inherent to PFI and CFI represent extreme
reference points. We introduce Relative Feature Importance (RFI), a
generalization of PFI and CFI that allows for a more nuanced feature
importance computation beyond the PFI versus CFI dichotomy. With RFI, the
importance of a feature relative to any other subset of features can be
assessed, including variables that were not available at training time. We
derive general interpretation rules for RFI based on a detailed theoretical
analysis of the implications of relative feature relevance, and demonstrate
the method’s usefulness on simulated examples.

###### Index Terms:

feature importance, interpretable machine learning, explainable artificial
intelligence, causality

##  I Introduction

Predictive modelling is increasingly deployed in high-stakes environments,
e.g., in the criminal justice system [11], loan approval [32], recruiting [9]
and medicine [27]. Due to legal regulations [29, 10] and ethical
considerations, ML methods need not only perform robustly in such environments
but also be able to justify their recommendations in a human-intelligible
fashion. This development has given rise to the field of interpretable machine
learning (IML) that involves studying methods that provide insight into the
relevance of features for model performance, referred to as feature
importance.
Prominent feature importance techniques include permutation feature importance
(PFI) [5, 12] and conditional feature importance (CFI) [25, 12, 19]. PFI is
based on replacing the feature of interest $X_{j}$ with a perturbed version
sampled from the marginal distribution $P(X_{j})$ while CFI perturbs $X_{j}$
such that the conditional distribution with respect to the set $R$ of
remaining features $P(X_{j}|X_{R})$ is preserved. The sampling strategy
defines the method’s reference point and therefore affects the method’s
implicit notion of relevance. While PFI quantifies the overall reliance of the
model on the feature of interest, CFI quantifies its unique contribution given
all remaining features.
While both PFI and CFI are useful, they fail to answer more nuanced questions
of feature importance. For instance, a stakeholder may be interested in the
importance of a feature relative to a subset of features. Also, the user may
want to know how important a feature is relative to variables that had not
been available at training time.
We suggest relative feature importance (RFI) as a generalization of PFI and
CFI that moves beyond the dichotomy between PFI, which breaks all dependencies
with features, and CFI, which preserves all dependencies with features. In
contrast to PFI and CFI, RFI is based on a perturbation that is restricted to
preserve the relationships with a set of variables $G$ that can be chosen
arbitrarily. We show that RFI is (1) semantically meaningful and (2)
practically useful.
We demonstrate the semantical meaning of RFI in Section IV. In particular, we
derive general interpretation rules that link nonzero RFI to (1) the
conditional dependence of the feature of interest with the target and non-
conditioned features $X_{\underline{R}}$ given the conditioned variables
$X_{G}$ in the data and (2) the conditional dependence of the input to the
feature of interest $X_{j}$ with the model’s prediction $\hat{Y}$ given fixed
inputs to the remaining features $X_{R}$ (Theorem 1). Furthermore, we show
that a nonzero difference between $\text{RFI}_{j}^{G}$ and
$\text{RFI}_{j}^{G\cup N}$, with $N$ being an arbitrary set disjunct with $G$,
implies the conditional dependence $X_{j}\not\Perp X_{N}|X_{G}$ (Theorem 2).
In Section V, we provide an implementation of RFI estimation that is based on
recent results from the related knockoff research field [7, 23]. Furthermore,
we translate the testing framework developed for conditional feature
importance [30] to RFI. We support our theoretical analysis and findings by
various simulation studies in Section VI. In particular, we show that RFI can
expose the indirect contribution of variables that are not directly used by
the model but provide information via dependent variables (Section VI-A).
Similarly, we show how RFI can be used to assess feature importance with
respect to variables not included at training time (Section VI-B).

###  I-A Contributions and Related Work

While conditioning on subsets of variables has been suggested before [25, 12],
the implications of this generalized variant of CFI have not yet been
rigorously analyzed. Some IML methods perturb or hide subsets of features,
e.g., in the context of multiple regression relative importance analysis is a
model-specific technique that averages over all importances of models trained
on feature subsets [6, 16]. Model-agnostic, local approximations to the
respective feature effect that avoid retraining and instead perturb subsets of
features have also been proposed [33, 17]. A very recent global, model-
agnostic feature importance proposal called SAGE quantifies feature importance
by perturbing multiple features [8].
While the aforementioned approaches are all based on removing several features
to provide more nuanced insight into the model, our proposal only modifies the
feature of interest. Our approach is model-agnostic and global, while most
aforementioned approaches are model-specific or local. The exception is the
global, model-agnostic SAGE [8], however the approaches are not only
computationally but also semantically different. E.g. our method assigns an
importance of zero for features that are not used by the model111A proof of
this property is given in Lemma 2., which is not the case for SAGE. While our
approach aims to provide nuanced insights into variable importance relative to
a specific set, SAGE aims to quantify the overall importance of variables for
the model.
Feature importance relative to variables that have not been included in the
training set has not been studied before. The indirect influence of variables
that the model does not computationally rely but statistically depend on has
been studied e.g. in [1].

##  II Background and Notation

###  II-A Notation

$\\{j\\}$$R$$\underline{R}$$G$$G^{*}$$\overline{G}$ Figure 1: Overview of our
notation.

We denote the target variable, i.e., the variable the model predicts, as $Y$
and feature variables by $X_{(.)}$. We refer to the variables as features to
emphasize when they were used in model training. Their observations are
denoted by $y$ and $x_{(.)}$. We use $D:=\\{1,\ldots,p\\}$ for the index set
of all features included in model training and $j$ for the index of our
feature of interest, $X_{j}$. The index set of the remaining variables is
denoted as $R:=D\backslash\\{j\\}$ (rest, remainder). The index set of
features, relative to which the importance of $X_{j}$ is considered, is
denoted as $G$. As $G$ can refer to any index set of variables, we denote its
intersection with $R$ as $\overline{G}=R\cap G$ and its complement as
$\underline{R}=R\backslash G$. We denote the index set of conditioning
variables that were not made available to the model during training as
$G^{*}=G\backslash R$.
In case we add new elements to the conditioning set $G$, we will denote this
set as $N$. The set may include variables within and outside $D$. The
respective components are denoted as $N^{*}=N\backslash R$ and as
$\overline{N}=R\cap N$. The remainder of $R$ without $G$ and $N$ is denoted as
$\underline{\underline{R}}=\underline{R}\backslash N$. We denote perturbed
variables of interest relative to $G$ as $\tilde{X}_{j}^{G}$. We refer to the
original and perturbed probability distribution of $X_{j}$ as the
observational and interventional distribution $P(X_{j},\dots)$ and
$P(\tilde{X}_{j}^{G},\dots)$. The inspected model is denoted as $f$, its
prediction as $\hat{Y}$. Independence of $Y$ and $X$ conditional on $Z$ is
denoted using $X\Perp Y|Z$, the respective conditional dependence as
$X\not\Perp Y|Z$.

###  II-B Feature Importance

Performance-based feature importance methods assess the relevance of a feature
of interest $X_{j}$ by assessing the impact of a perturbation of $X_{j}$ on
the model’s performance. Local feature importance methods focus on the
importance of features for specific data points, whereas global feature
importance methods assess the impact over the whole domain. In the following,
we focus on global methods.
Global feature importance is computed according to the following general
schemata:

$\text{FI}_{j}=\tilde{\mathcal{R}}^{j}-\mathcal{R}\text{ or
}\text{FI}_{j}=\frac{\tilde{\mathcal{R}}^{j}}{\mathcal{R}}$

where we denote the original risk of the model and the risk after perturbing
$X_{j}$ as $\mathcal{R}$ and $\tilde{\mathcal{R}}^{j}$, respectively. For
estimation, the true risk $\mathcal{R}$ is replaced with the empirical risk
$\mathcal{R}_{\text{emp}}$.
Feature importance methods furthermore differ in how they perturb and whether
they rely on retraining the model. While some methods retrain the model after
the perturbation (e.g. LOCO, [15]), others evaluate the impact of the
perturbation on the same original model (e.g. [5, 25]). In this work, we focus
on methods that avoid retraining.
For methods that avoid retraining, we observe a dichotomy between two general
perturbation approaches: resampling that preserves the marginal and resampling
that preserves the conditional distribution. Marginal resampling was
originally proposed to compute perturbed versions of $X_{j}$ by permuting the
observations $x_{j}^{(i)}$ within the sample [5]. The respective sample breaks
the dependence between $X_{j}$ and $(Y,X_{R})$ while preserving the marginal
distribution $P(X_{j})$. More recently, Model Reliance was proposed [12],
which takes the expectation over all possible permutations. Resampling from
the marginal distribution has been criticized to introduce bias, in particular
because it overestimates the importance of correlated variables [25],
resulting in incorrect feature rankings [26]. It also leads to extrapolation
under dependent features [14, 19], i.e. conclusions about the model are being
drawn using unrealistic data points on which the model was not trained. CFI,
on the other hand, samples from the conditional distribution $P(X_{j}|X_{R})$
[25, 28, 2, 7, 12, 14, 19]. A large variety of model-specific methods exist
[13, 31]. Conditional variants quantify the importance of a feature given the
information that all remaining features $R$ contain about $X_{j}$ [20],
thereby avoiding evaluation of the model on unrealistic datapoints [19].

##  III Relative Feature Importance

Relative Feature Importance is a general framework that assesses feature
importance relative to arbitrary variable sets $G$. The frameworks subsumes
PFI and CFI as two extreme special cases.
In PFI, $X_{j}$ is replaced with a perturbed version that preserves the
marginal distribution $P(X_{j})$ while breaking the dependencies with $Y$ and
all features. In CFI, a perturbed version of $X_{j}$ is used that preserves
the conditional distribution $P(X_{j}|X_{R})$, thereby only breaking
conditional dependence between $X_{j}$ and $Y$ given all features. As our
analysis in Section IV establishes, the replacement strategies of PFI and CFI
define extreme reference points. CFI quantifies the contribution relative to
all remaining features $R$, whereas PFI regards a feature in isolation.
We go beyond the PFI versus CFI dichotomy. We argue that it is (1) meaningful
(Section IV) and (2) practically useful (Section VI) to replace $X_{j}$ with
perturbed versions that preserve the conditional distribution $P(X_{j}|X_{G})$
with respect to arbitrary sets $G$ while requiring
$\tilde{X}_{j}^{G}\Perp(X_{\underline{R}},Y)|X_{G}$. $G$ can be a subset of
$R$, but can also include variables not available at training time such that
$G\backslash R\neq\emptyset$. We term the resulting method Relative Feature
Importance (RFI):

######  Definition 1 (Relative Feature Importance – RFI)

We define Relative Feature Importance with respect to a feature set $G$ with
$Y\notin G$ and a fixed model $f$ as

$\text{RFI}_{j}^{G}:=\tilde{\mathcal{R}}^{j|G}-\mathcal{R},$

where $\tilde{\mathcal{R}}^{j|G}:=\mathcal{R}(Y,f(X_{R},\tilde{X}_{j}^{G}))$
is the risk w.r.t. to a replacement variable $\tilde{X}_{j}^{G}$ and
$\mathcal{R}=\mathcal{R}(Y,f(X_{j},X_{R}))$ refers to the original risk. The
replacement variable has to satisfy

  * •

$\tilde{X}_{j}^{G}\sim P(X_{j}|X_{G})$ and

  * •

$\tilde{X}_{j}^{G}\Perp(X_{\underline{R}},Y)|X_{G}$.

In the following section, we discuss the semantic meaning of RFI. The
estimation of RFI is discussed in Section V.

##  IV Interpreting Relative Feature Importance

IML techniques aim to provide insight into the model and, possibly, into the
underlying data generating mechanism. However, IML techniques themselves are
subject to interpretation. The characterization of an IML method by its
mathematical definition is computationally precise, but has limited aid in
guiding users to make conclusions about the underlying model and data. In this
section we provide a (non-comprehensive) list of interpretation rules for RFI,
that characterize the method by how it behaves in its context. This context
includes both the model and the underlying data generating mechanism. More
specifically, we link RFI to (conditional) independence in the underlying data
set as well as to whether the model’s prediction $\hat{Y}$ is constant in the
argument $x_{j}$ for a fixed $x_{R}$. While RFI can be used for quantification
of feature importance, we focus our analysis on relevance as a binary property
and characterize relative feature relevance ($\text{RFI}\neq 0$). We show that
the implicit notion of relevance of RFI is defined by the choice of $G$. By
modifying the conditioning set $G$ beyond the PFI versus CFI dichotomy, we are
able to gain insight into more nuanced aspects of the model and the data
generating mechanism. The main results are given in Theorem 1 and Theorem 2.
Furthermore, we highlight limitations stemming from the choice of the loss
function $L$ and the model fit for the interpretation, which are, in our
humble opinion, underrepresented in the current discussion.
We structure our analysis by taking the user’s perspective and asking ”What
can we infer from relative feature relevance?”.

###  IV-A Implications of Relative Feature Relevance

In the following, we analyze the implications of RFI without further
assumptions about model and data. We thereby distinguish between two levels of
explanation. Relative feature relevance provides insight, both into model and
data.

###### Theorem 1

If $RFI_{j}^{G}\neq 0$ then

  * •

$X_{j}\not\Perp(Y,X_{\underline{R}})|X_{G}$ in the underlying distribution
(data level)

  * •

$\tilde{X}_{j}\not\Perp\hat{Y}|X_{R}$ w.r.t. the interventional distribution
$P(X_{j}|X_{G})P(X_{G},X_{\underline{R}})>0$ (model level)

We prove Theorem 1 in two steps. First, we assess the implications of the
respective independence for the underlying data set (Lemma 1). Then, we assess
the implications of the respective independence for the model (Lemma 2). The
contrapositions yield Theorem 1.

###### Lemma 1

If $X_{j}\Perp(Y,X_{\underline{R}})|X_{G}$ for any G with $Y\notin G$ then
$RFI_{j}^{G}=0$.

We base the proof of Lemma 1 on the insight that (because the model $f$ is
fixed) an equivalence in distribution implies an equivalence in risk
(Proposition 1). Therefore conditions under which the interventional
distribution $P(\tilde{X}_{j}^{G},X_{R},Y)$ coincides with the original
distribution $P(X_{j},X_{R},Y)$ are sufficient for $RFI=0$.

###### Proposition 1

If observational and interventional distribution coincide, then risks with and
without perturbation are equal:

$P(Y,X_{j},X_{R})=P(Y,\tilde{X}_{j}^{G},X_{R})\Rightarrow\mathcal{R}(f)=\tilde{\mathcal{R}}^{j|G}(f)$

###### Proof:

Given that $P(Y,X_{j},X_{R})=P(Y,\tilde{X}_{j},X_{R})$ we can write

$\displaystyle\mathcal{R}(f)$
$\displaystyle=\mathbb{E}_{Y,X_{j},X_{R}}[L(Y,f(X_{j},X_{R}))]$
$\displaystyle=\mathbb{E}_{Y,\tilde{X}_{j},X_{R}}[L(Y,f(\tilde{X}_{j},X_{R}))]=\tilde{\mathcal{R}}(f).$

∎

We show next that the conditional independence
$X_{j}\Perp(X_{\underline{R}},Y)|X_{G}$ is a sufficient condition for identity
of both distributions.

###### Proof:

It holds that

$\displaystyle P(Y,X_{j},X_{\underline{R}},X_{G})$ $\displaystyle=$
$\displaystyle P(X_{j}|Y,X_{\underline{R}},X_{G})P(Y,X_{\underline{R}},X_{G})$
$\displaystyle\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny$X_{j}\Perp(X_{\underline{R}},Y)|X_{G}$}}}{=}}$
$\displaystyle P(X_{j}|X_{G})P(Y,X_{\underline{R}},X_{G})$
$\displaystyle\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny(def)}}}{=}}$
$\displaystyle P(\tilde{X}_{j}^{G}|X_{G})P(Y,X_{\underline{R}},X_{G})$
$\displaystyle=$ $\displaystyle
P(\tilde{X}_{j}^{G},Y,X_{\underline{R}},X_{G}).$

Using Proposition 1 we can infer that $RFI_{j}^{G}=0$.
∎

So far, we have assessed implications for the underlying data generating
mechanism. Next, we assess implications for the inspected model $f$.

###### Lemma 2

If $\tilde{X}_{j}^{G}\Perp\hat{Y}|X_{R}$ w.r.t. the interventional
distribution $P(\tilde{X}_{j}^{G},X_{G},X_{\underline{R}})$ then
$\text{RFI}_{j}^{G}=0$ for any $G$.

###### Proof:

If the prediction for an observation $(x_{1},\dots,x_{p})$ is independent of
the value $x_{j}^{\prime}$ w.r.t. the interventional distribution, the
prediction is unaffected when replacing $x_{j}$ with any value
$x_{j}^{\prime}$ with
$P(x_{j}^{\prime}|X_{G}=x_{G})P(X_{G}=x_{G},X_{\underline{R}}=x_{\underline{R}})>0$.
Consequently, any sample from $\tilde{X}_{j}^{G}$ yields the same prediction.
Furthermore values $x_{j}^{\prime}$ with nonzero probability over the
interventional distribution also have nonzero probability over the
observational distribution. The interventional distribution can be rewritten
as

$\displaystyle P(\tilde{X}_{j}^{G},X_{G},X_{\underline{R}})$
$\displaystyle=P(\tilde{X}_{j}^{G}|X_{G},X_{\underline{R}})P(X_{G},X_{\underline{R}})$
$\displaystyle=P(\tilde{X}_{j}^{G}|X_{G})P(X_{G},X_{\underline{R}})$
$\displaystyle=P(X_{j}|X_{G})P(X_{G},X_{\underline{R}}).$

Similarly, the observational distribution can be factorized into
$P(X_{j}|X_{G},X_{\underline{R}})P(X_{G},X_{\underline{R}})$. As
$P(X_{j}|X_{G},X_{\underline{R}})>0\Rightarrow P(X_{j}|X_{G})>0$ (which can be
derived from, e.g., the law of total probability) it follows that
$P(\tilde{X}_{j}^{G},X_{G},X_{\underline{R}})>0\Rightarrow
P(X_{j},X_{G},X_{\underline{R}})>0$.
Consequently the prediction $\hat{y}$ for any value $x_{j}$ with positive
probability $P(X_{j}=x_{j}|X_{R}=x_{R})$ is identical given unchanged $x_{R}$.
As the conditional distributions of $X_{j}$ and $\tilde{X}_{j}^{G}$ overlap
and the distribution of $X_{R}$ is unaffected, the prediction $\hat{Y}$ is
identical with and without perturbation. Therefore
$\mathcal{R}=\tilde{\mathcal{R}}^{j|G}$ and $RFI_{j}^{G}=0$.
∎

To summarize, we have shown that independence on the dataset and on the model
level respectively imply $\text{RFI}_{j}^{G}=0$ and can thereby prove Theorem
1.

###### Proof:

The result follows from contraposition of Lemma 1 and contraposition of Lemma
2.
∎

Theorem 1 shows that nonzero $\text{RFI}_{j}^{G}$ implies dependencies between
sets of variables on the model level as well as on the data level. Which
dependencies are relevant for $\text{RFI}_{j}^{G}$ can be controlled with the
conditioning set $G$. Consequently, the conditioning set $G$ determines the
method’s implicit definition of relevance. I.e., on the data level, if
$X_{j}\Perp(X_{\underline{R}},Y)|X_{G}$ holds, $\text{RFI}_{j}^{G}$ is zero
irrespective of any other dependencies that may hold, e.g. with $X_{G}$ (Lemma
1). Nonzero RFI, a difference in performance on interventional and
observational distribution, can only be caused by dependencies that have been
destroyed in the interventional distribution, the dependencies with and via
$X_{G}$ are preserved by the replacement $\tilde{X}_{j}^{G}$ and can therefore
not be responsible for $\text{RFI}_{j}^{G}\neq 0$. Similarly, on the model
level, $\tilde{X}_{j}^{G}\Perp\hat{Y}|X_{R}$ over the interventional
distribution $P(X_{j}|X_{G})P(X_{G},X_{\underline{R}})$ yields zero RFI (Lemma
2). The behavior of the model outside the domain in which it is evaluated is
irrelevant for $RFI_{j}^{G}$. What domain the model is evaluated over depends
on the choice of $G$.
Because we can control RFI’s implicit definition of relevance with $G$, RFI
allows more nuanced insights into model and data than PFI or CFI alone. In
Theorem 1, we aim to make the implicit definition of relevance explicit. On
the data level, nonzero RFI implies the dependence of $X_{j}$ with the tuple
$(Y,X_{\underline{R}})$ given $X_{G}$
($X_{j}\not\Perp(Y,X_{\underline{R}})|X_{G}$). In order to understand the
aforementioned dependence, using the graphoid axioms contraction and weak
union [22], the equivalent formulation below can be adduced:

$(X_{j}\not\Perp Y|X_{G})\vee(X_{j}\not\Perp X_{\underline{R}}|X_{G},Y).$

At least one of the two conditional dependencies has to hold for nonzero
$\text{RFI}_{j}^{G}$. The first dependence can be rephrased as: $X_{j}$ is
informative of $Y$, even if we already know $X_{G}$. It is more difficult to
make sense of the second dependence. Under dependent features $(X_{j}\not\Perp
X_{\underline{R}}|X_{G},Y)$, the distribution of $X_{j}$ with
$X_{\underline{R}}$ is not preserved under perturbation $\tilde{X}_{j}^{G}$.
In the interventional distribution $P(\tilde{X}_{j}^{G},X_{\underline{R}})$
observations that are improbable or impossible w.r.t. the observational
distribution $P(X_{j},X_{\underline{R}})$ can be possible and probable (and
vice versa). Consequently, in the interventional distribution the feature
distribution differs from the observation feature distribution. Even if
$X_{j}\Perp Y|X_{G}$ holds, the model may perform suboptimally due to this
distribution shift and cause $\text{RFI}_{j}^{G}$ nonzero222Let e.g.
$X_{1},X_{2}$ be perfectly correlated and independent of $Y$. Then adding
$X_{1}-X_{2}$ does not alter its prediction performance, unless the dependence
between the variables is broken. Also see [14] for a discussion in PFI.. If
the conditioning set is a superset of $R$ ($G\supseteq R$), such that set of
remaining variables $X_{\underline{R}}$ is empty, it holds that $(X_{j}\Perp
X_{\underline{R}}|X_{G},Y)$. Therefore nonzero RFI must be attributed to
$(X_{j}\not\Perp Y|X_{G})$ for $G\supseteq R$.
On the model level, nonzero RFI implies that the model’s predictions are
conditionally dependent on $\tilde{X}_{j}^{G}$ given the remaining features
$R$ are fixed. E.g. for a linear model that has coefficient zero for all terms
involving $X_{j}$, this dependence would not be fulfilled, and
$\text{RFI}_{j}^{G}$ would be zero (Lemma 2). The model is evaluated over the
interventional distribution $P(X_{j}|X_{G})P(X_{G},X_{\underline{R}})>0$,
which varies depending on $G$. If $G$ contains a nearly perfect correlate of
$X_{j}$, $X_{j}$ can be reconstructed well. In contrast, if $G=\emptyset$, for
every possible $x_{R}$ the model is evaluated over the whole marginal
distribution of $X_{j}$. Although choosing a smaller set $G\subset R$ leads to
extrapolation under dependent features, it allows more insight into the
model’s mechanism. For interpretation purposes like safety, this is highly
desirable.
In the preceding paragraphs we have highlighted the importance of the
conditioning set $G$ for the method’s implicit notion of relevance and
illustrated the results from Theorem 1. We have argued that the conditioning
set controls which potential dependencies can be responsible for nonzero
$\text{RFI}_{j}^{G}$. The insights lead to a further, interesting application
of RFI. By assessing the difference $\Delta RFI_{j}^{G\to G\cup
N}=\text{RFI}_{j}^{G}-\text{RFI}_{j}^{G\cup N}$ when modifying the
conditioning set $G$ by adding new elements $N$, we are able to assess the
role of the dependencies with variables in $N$ relative to a baseline $G$.
While for $\text{RFI}_{j}^{G}$ only dependencies of $X_{j}$ with and via $G$
are preserved, for $\text{RFI}_{j}^{G\cup N}$ also dependencies with and via
$N$ are maintained. If $\Delta RFI_{j}^{G\to G\cup N}$ is nonzero, this change
has to be due to dependencies involving $N$, but not $G$. We substantiate this
claim with Theorem 2. In order for $\Delta RFI_{j}^{G\to G\cup N}$ to be
positive, the dependence $X_{j}\not\Perp X_{N}|X_{G}$ has to hold.

###### Theorem 2

If the difference $\Delta RFI_{j}^{G\to G\cup N}=RFI_{j}^{G}$ \-
$RFI_{j}^{G\cup N}\neq 0$, then $X_{j}\not\Perp X_{N}|X_{G}$.

###### Proof:

Under independence $X_{j}\Perp X_{n}|X_{G}$ it holds that

$\displaystyle
P(\tilde{X}_{j}^{G},Y,X_{\underline{\underline{R}}},X_{G},X_{N})$
$\displaystyle=P(\tilde{X}_{j}^{G}|Y,X_{\underline{\underline{R}}},X_{G},X_{N})P(Y,X_{\underline{\underline{R}}},X_{G},X_{N})$
$\displaystyle\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny(def
$\tilde{X}_{j}^{G}$)}}}{=}}P(X_{j}|X_{G})P(Y,X_{\underline{\underline{R}}},X_{G},X_{N})$
$\displaystyle\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny$X_{j}\Perp
X_{n}|X_{G}$}}}{=}}P(X_{j}|X_{G},X_{N})P(Y,X_{\underline{\underline{R}}},X_{G},X_{N})$
$\displaystyle\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny(def
$\tilde{X}_{j}^{G\cup N}$)}}}{=}}P(\tilde{X}_{j}^{G\cup
N}|X_{G},X_{N})P(Y,X_{\underline{\underline{R}}},X_{G},X_{N})$
$\displaystyle\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny(def
$\tilde{X}_{j}^{G\cup N}$)}}}{=}}P(\tilde{X}_{j}^{G\cup
N}|Y,X_{G},X_{N},X_{\underline{\underline{R}}})P(Y,X_{\underline{\underline{R}}},X_{G},X_{N})$
$\displaystyle=P(\tilde{X}_{j}^{G\cup
N},Y,X_{\underline{\underline{R}}},X_{G},X_{N})$

The equality
$P(\tilde{X}_{j}^{G},Y,X_{\underline{\underline{R}}},X_{G},X_{N})=P(\tilde{X}_{j}^{G\cup
N},Y,X_{\underline{\underline{R}}},X_{G},X_{N})$ implies
$P(\tilde{X}_{j}^{G},Y,X_{R})=(\tilde{X}_{j}^{G\cup N},Y,X_{R})$. Invoking
Proposition 1 it holds that the corresponding risks $\mathcal{R}^{j|G}$ and
$\mathcal{R}^{j|G\cup N}$ are equal. As
$\text{RFI}_{j}^{G}-\text{RFI}_{j}^{G\cup
N}=\mathcal{R}^{j|G}-\mathcal{R}^{j|G\cup N}$ it holds that $X_{j}\not\Perp
X_{n}|X_{G}\Rightarrow\Delta RFI_{j}^{G\to G\cup N}=0$. Contraposition proves
Theorem 2.
∎

While nonzero $\text{RFI}_{j}^{G}$ as well as nonzero $\Delta RFI_{j}^{G\to
G\cup N}$ have clear implications, interpreting zero $\text{RFI}_{j}^{G}$ or
zero $\Delta RFI_{j}^{G\to G\cup N}$ is difficult. For example, we may be
tempted to interpret $\text{RFI}_{j}^{G}=0$ as conditional independence in the
data. However, the general principle that absence of evidence is no evidence
for absence also applies in the context of RFI. A dependence in the data may
not be captured by the model when it has a poor fit and does not rely on the
respective variable. Similarly, although $f$ may be optimal, a dependence in
higher moments may simply not be modeled by $f$ or captured by the loss $L$.
As all aforementioned causes of nonzero RFI are potentially sufficient, but
not necessary, it is unclear which of the causes nonzero RFI can be attributed
to. Furthermore, the related problem of conditional independence testing is
provably hard [24].
The theoretical insights that we derive in this Section (Theorem 1 and 2) are
applied and illustrated in a simulation study in Section VI.

##  V Estimation and Testing

Estimating and sampling from the conditional distribution is in general
difficult, especially in high-dimensional continuous settings. Various
approaches for replacing $X_{j}$ with samples from its conditional
distribution exist, e.g., knockoff approaches [2, 7, 23], imputation and
weighting [12] or permutation within decision tree leaves [18]. We used
Model-X knockoffs [7] in this work, but note that the RFI approach is agnostic
to its algorithmic implementation.
Using (standard) empirical risk estimates, our RFI estimate is

$\hat{\text{RFI}}_{j}^{G}=\frac{1}{n}\sum_{i=1}^{n}L\left(y^{(i)},f(\tilde{x}^{(i)}_{j},x^{(i)}_{R})\right)-\frac{1}{n}\sum_{i=1}^{n}L\left(y^{(i)},f(x^{(i)}_{j},x^{(i)}_{R})\right)$

where $\tilde{x}^{(i)}_{j}$ is a sample from $\tilde{X}_{j}^{G}$. We can then
test for nonzero $\text{RFI}_{j}^{G}$ using procedures for conditional
independence tests, e.g., [30], thereby quantifying the uncertainty coming
from empirical risk minimization. Because of the central limit theorem, the
empirical risk converges (in probability) to a Gaussian distribution with
increasing number of observations. Therefore, one-sided, paired t-tests can be
used to infer tests and confidence intervals [30]. The test procedures
proposed in [30] are agnostic to the conditioning set for the perturbation
$\tilde{X}_{j}^{G}$. For smaller samples, the Exact Test by Fisher may be
used.
The t-test and Fisher Exact Test ignore uncertainty and bias of the estimation
procedures, i.e. the ML model and the knockoff-sampler are treated as “fixed”.
E.g. misspecified, suboptimal models may not capture dependencies. Or
dependencies are in higher moments that are not captured by the loss.
Consequently, without further assumptions, the framework does not provide a
test for conditional independence in the dataset.
The popular testing procedures for knockoffs proposed by [7] provide FDR over
all features, but does not test the significance of the importance of
individual features.

##  VI Simulation Studies

In the following, we demonstrate the usefulness of RFI on two simulation
studies. In the first example, we use RFI to expose indirect influence of
variables that are not computationally used by the model. In the second
example, we assess feature importance relative to a confounder that was
unavailable at training time. In both examples, we represent the underlying
data generating mechanism, that gives rise to the dependencies in the data,
with a causal directed acyclic graph (DAG). The code for the examples is
available online333Link to Code: https://github.com/gcskoenig/icpr2020-rfi.

###  VI-A Indirect Influence

A prominent application of interpretable machine learning is auditing models
regarding its reliance on protected attributes $A$ like age or sex. A reliance
on the respective attributes may result in unfair discrimination and requires
further inspection. With approaches like fairness through unawareness [3], the
model does not rely on protected attributes directly. However, by implicitly
reconstructing the sensitive attributes using seemingly harmless correlates,
the model can indirectly make use of the protected attribute resulting in
potentially harmful, unfair discrimination [3].
PFI and CFI cannot expose such indirect influence. As Lemma 2 proves,
$RFI_{A}^{G}$ is zero for a model that does not (directly) use the feature of
interest $A$ for the prediction for any conditioning set $G$. Furthermore,
from PFI and CFI alone, we cannot infer whether the importance of a variable
can be attributed to its dependence with an indirect influence. Using
$\text{RFI}_{j}^{G}$ with $G=A$ we preserve the influence of $A$ on the
prediction and can thereby restrict the attribution of importance to
contributions stemming from dependencies not involving $A$ (Theorem 1, Lemma
1). The difference to $\Delta RFI_{j}^{G\to G\cup N}$ with $G=\emptyset$ and
$N=A$ exposes the indirect influence.
Not every indirect influence from a sensitive attribute is considered
undesirable. Certain correlates of $A$ may indeed be valid criteria for a
decision (e.g. [4]). Importance stemming from dependencies with $A$ via such
resolving variables $Z$ would be considered acceptable. We can assess the
indirect influence beyond contributions stemming from dependence via $Z$ by
comparing to a baseline $G=Z$. In this baseline, contributions via $Z$ are
preserved and therefore irrelevant for RFI. Consequently, when setting $N=A$,
the difference $\Delta RFI_{j}^{G\to G\cup N}$ only quantifies indirect
influence that is not resolved by $Z$.
We demonstrate the usefulness of RFI to expose indirect influence in a
simulation study. The dataset is a sample drawn from the distribution induced
by a structural causal model (SCM) depicted in Figure 2. All relationships are
additive linear with coefficients $1$ and Gaussian noise terms
($\sigma_{1}=\sigma_{2}=\sigma_{4}=1$, $\sigma_{3}=0.3$ and $\sigma_{y}=0.5$).
An ordinary least squares linear regression model was fit to predict $Y$ from
$X_{1},\dots,X4$ (MSE = $0.25$,
$f(x_{1},x_{2},x_{3},x_{4})=0.00x_{1}-0.01x_{2}+1.01x_{3}+1.00x_{4}$). We
trained model-X knockoffs [7] on the training data and evaluated RFI on test
data. Sample size is $10^{5}$ with $10\%$ test data.
In order to quantify the direct influence of the features we compute PFI. As
we can see in Figure 3, $X_{1}$ and $X_{2}$ are considered irrelevant. In
order to expose their indirect influence, we additionally compute RFI with
respect to $G=\\{X_{1}\\}$ and $G=\\{X_{2}\\}$ respectively. For both
variables we observe a drop in importance of $X_{3}$ and $X_{4}$. Consequently
both $X_{1}$ and $X_{2}$ have an indirect influence on the target (Theorem 2).
Furthermore we are interested in whether the indirect influence of $X_{1}$ can
be resolved by $X_{2}$. We therefore compute $\text{RFI}_{j}^{G\cup N}$ with
$G=\\{X_{2}\\}$ and $N=\\{X_{1}\\}$. We see that for $X_{3}$ no change in
importance can be observed. This is due to the independence $X_{1}\Perp
X_{3}|X_{2}$444As faithfulness and causal markov condition hold,
$d$-separation in the graph and (conditional) independence coincide[21]. We
can therefore read the independence structures off Figures 2 and 4. (Theorem
2). The indirect influence is resolved. However, for $X_{4}$ the importance
decreases further and is therefore not resolved by $X_{2}$. This is in
alignment with the dependence $X_{1}\not\Perp X_{4}|X_{2}$ implied by the
graph (Figure 2).

$X_{1}$$X_{2}$$X_{3}$$X_{4}$$Y$ Figure 2: Variable $X_{1}$ influences $Y$ both
via the chain $X_{2}\rightarrow X_{3}$ and via $X_{4}$. $X_{1}$ may be some
undesired influence, and $X_{2}$ a variable resolving the undesired influence.
We find that the prediction can nevertheless be influenced via $X_{4}$ by
comparing $RFI_{4}^{X_{2}}$ with $RFI_{4}^{X_{2},X_{1}}$ (Figure 3). All
relationships are additive linear Gaussian with all coefficients being equal
to $1$ and $\sigma_{1}=\sigma_{2}=\sigma_{4}=1$, $\sigma_{3}=0.3$ and
$\sigma_{y}=0.5$. Figure 3: RFI’s for a linear regression model fitted on the
dataset illustrated in Figure 2. Feature importance values are averaged over
30 runs and rounded. Feature importance values are averaged over 30 runs and
rounded. We evaluated significance using a t-test for the first run. All
positive features were significant at $\alpha=0.01$, whereas for all zero RFI
values the null could not be rejected. For $X_{1}$ and $X_{2}$ all RFIs are
zero, whereas for $X_{3}$ and $X_{4}$ RFIs are positive. We see that $X_{1}$
and $X_{2}$ both have an indirect influence on $X_{3}$ and $X_{4}$, but that
$X_{2}$ can resolve the influence of $X_{1}$ on $X_{3}$.

###  VI-B Variables Outside Training Set

When designing a model $f$, a practitioner may have decided to exclude a
variable from the feature set, e.g., because it was then considered
irrelevant, it belongs to a different modality or would have required further
preprocessing. Furthermore, when auditing a machine learning model $f$,
variables that have not been available for the training of the model may be
accessible.
In this example, we demonstrate that variables outside the training set can be
included in the conditioning set for RFI. Consequently, importance of the
features relative to variables outside the training set and the indirect
influence of such variables can be assessed. More specifically, we simulate a
hypothetical situation where the influence of a previously unknown confounder
$C$ shall be evaluated. This variable $C$ is available for the model audit. In
particular, we wonder whether the features $X_{1}$, $X_{2}$ and $X_{3}$ are
only or partly important due to a dependence via $C$.
The dataset was sampled from a structural causal model (SCM) depicted in
Figure 4. Assuming faithfulness and the causal Markov condition, this DAG
implies the following (conditional) (in-)dependencies: $X_{1}$ is independent
of $C$, $X_{3}$ is independent of $Y$ conditional on $C$, and $X_{2}$ is
dependent on $Y$. Note that the dependence between $X_{2}$ and $Y$ is due to
the common cause $C$ as well as due to a direct effect of $X_{2}$ on $Y$. All
relationships are additive linear with coefficients $1$ and additive Gaussian
noise ($\sigma_{1}=\sigma_{2}=\sigma_{C}=1.0$ and
$\sigma_{3}=\sigma_{Y}=0.5$). We fit an ordinary least squares linear
regression model on $X_{1}$, $X_{2}$ and $X_{3}$ to predict $Y$ (MSE = $0.40$,
$f(x_{1},x_{2},x_{3})=1.0x_{1}+1.17x_{2}+0.67x_{3}$). $C$ was not available
for model training. We trained Model-X knockoffs [7] on training data and
sampled from $\tilde{X}_{j}^{G}$ on test data. Sample size is $10^{5}$ with
$10\%$ test data.
When computing $\text{RFI}_{j}^{C}$ ($G=\\{C\\}$) for each variable, the
different relationships with $C$ become apparent. The respective results are
depicted in Figure 5. For $X_{1}$ the feature importance relative to $C$
remains unchanged as the variables are pairwise independent (Theorem 2). For
$X_{3}$, that is only dependent with $Y$ via $C$, it completely vanishes
(Lemma 1). For $X_{2}$ the feature importance decreases but remains nonzero,
as $X_{2}$ is dependent with $Y$ directly and via $C$.
Consequently, using RFI, we can (1) identify variables that are important due
to a variable unavailable at training time and (2) distinguish between
variables that only depend on $Y$ via $C$ from those that do not. With PFI
($G=\emptyset$) or CFI ($G=R$) such a distinction is in general not possible.

$C$$X_{1}$$X_{2}$$X_{3}$$Y$$C$$X_{1}$$X_{2}$$X_{3}$$Y$ Figure 4: Left: We see
the causal graph $\mathcal{G}$ corresponding to the Structural Causal Model
that was used to generate the dataset used in Figure 5. All relationships are
additive linear Gaussian with all coefficients equal to $1$ and
$\sigma_{1}=\sigma_{2}=\sigma_{C}=1.0$ and $\sigma_{3}=\sigma_{Y}=0.5$. Right:
Pairwise dependencies after conditioning on $C$. Figure 5: Feature Importance
results corresponding to the dataset depicted in Figure 4. We averaged RFI
over 30 runs. RFI for $X_{1}$ is unaffected by changes in $G$, for $X_{2}$ RFI
drops with $C$ is added to $G$. For $X_{3}$ RFI vanishes relative to $C$. For
all except for $RFI_{X_{3}}^{C}$ the null can be rejected at $\alpha=0.01$ in
the first run.

##  VII Discussion

We proposed relative feature importance (RFI), a general conditional feature
importance framework which allows to condition on arbitrary sets of other
features, including features outside the training set. We underpin the method
with theoretical results allowing insight into both model and underlying
dataset. In a simulation study, the usefulness of the method for the exposure
of indirect influence is demonstrated.
Relative feature importance requires sampling from (unknown) conditional
distributions. For continuous variables and in high-dimensional settings this
task is challenging and an open area of research [7, 23]. Uncertainty stemming
from inaccurate sampling may affect the interpretation. The quality of insight
into the underlying dataset strongly depends on the training and evaluation of
the model. Dependencies in higher moments are usually not modeled and not
captured by standard loss functions and can therefore not be detected.
Especially the interpretation of zero RFI requires careful assessment of the
model specification. Further research is needed to assess necessary
assumptions for the interpretation of RFI. These challenges are not unique to
RFI, but apply more generally in the field of interpretable machine learning
[20].

## References

  * [1] Philip Adler, Casey Falk, Sorelle A. Friedler, Tionney Nix, Gabriel Rybeck, Carlos Scheidegger, Brandon Smith, and Suresh Venkatasubramanian.  Auditing black-box models for indirect influence.  Knowledge and Information Systems, 54(1):95–122, 2018.  arXiv: 1602.07043.
  * [2] Rina Foygel Barber, Emmanuel J Candès, and others.  Controlling the false discovery rate via knockoffs.  The Annals of Statistics, 43(5):2055–2085, 2015.  Publisher: Institute of Mathematical Statistics.
  * [3] Solon Barocas, Moritz Hardt, and Arvind Narayanan.  Fairness and Machine Learning.  fairmlbook.org, 2019.  http://www.fairmlbook.org.
  * [4] Vence L Bonham, Shawneequa L Callier, and Charmaine D Royal.  Will precision medicine move us beyond race?  The New England journal of medicine, 374(21):2003, 2016.
  * [5] Leo Breiman.  Random forests.  Machine Learning, pages 1–122, 2001.
  * [6] David V Budescu.  Dominance analysis: a new approach to the problem of relative importance of predictors in multiple regression.  Psychological bulletin, 114(3):542, 1993.
  * [7] Emmanuel Candès, Yingying Fan, Lucas Janson, and Jinchi Lv.  Panning for gold: ‘model-X’ knockoffs for high dimensional controlled variable selection.  Journal of the Royal Statistical Society. Series B: Statistical Methodology, 80(3):551–577, 2018.  arXiv: 1610.02351.
  * [8] Ian Covert, Scott Lundberg, and Su-In Lee.  Understanding Global Feature Contributions Through Additive Importance Measures.  arXiv preprint arXiv:2004.00668, 2020.
  * [9] Jeffrey (Reuters) Dastin.  Amazon scraps secret AI recruiting tool that showed bias against women.  Reuters, 2018.
  * [10] Lydia de la Torre.  A Guide to the California Consumer Privacy Act of 2018.  SSRN Electronic Journal, pages 1–17, 2018.
  * [11] Julia Dressel and Hany Farid.  The accuracy, fairness, and limits of predicting recidivism.  Science Advances, 4(1):1–6, 2018.
  * [12] Aaron Fisher, Cynthia Rudin, and Francesca Dominici.  All models are wrong, but many are useful: Learning a variable’s importance by studying an entire class of prediction models simultaneously.  Journal of Machine Learning Research, 20(177):1–81, 2019.
  * [13] Ulrike Grömping.  Variable importance assessment in regression: linear regression versus random forest.  The American Statistician, 63(4):308–319, 2009.  Publisher: Taylor & Francis.
  * [14] Giles Hooker and Lucas Mentch.  Please Stop Permuting Features: An Explanation and Alternatives.  arXiv preprint arXiv:1905.03151v, pages 1–15, 2019.  arXiv: 1905.03151v1.
  * [15] Jing Lei, Max G’Sell, Alessandro Rinaldo, Ryan J Tibshirani, and Larry Wasserman.  Distribution-free predictive inference for regression.  Journal of the American Statistical Association, 113(523):1094–1111, 2018.
  * [16] Stan Lipovetsky and Michael Conklin.  Analysis of regression in game theory approach.  Applied Stochastic Models in Business and Industry, 17(4):319–330, 2001.  Lipovetsky2001.
  * [17] Scott M. Lundberg and Su In Lee.  A unified approach to interpreting model predictions.  Advances in Neural Information Processing Systems, 2017-Decem(Section 2):4766–4775, 2017.  arXiv: 1705.07874.
  * [18] Christoph Molnar, Gunnar König, Bernd Bischl, and Giuseppe Casalicchio.  Model-agnostic feature importance and effects with dependent features–a conditional subgroup approach.  arXiv preprint arXiv:2006.04628, 2020.
  * [19] Christoph Molnar, Gunnar König, Bernd Bischl, and Giuseppe Casalicchio.  Model-agnostic feature importance and effects with dependent features-a conditional subgroup approach.  arXiv preprint arXiv:2006.04628, 2020.
  * [20] Christoph Molnar, Gunnar König, Julia Herbinger, Timo Freiesleben, Susanne Dandl, Christian A. Scholbeck, Giuseppe Casalicchio, Moritz Grosse-Wentrup, and Bernd Bischl.  Pitfalls to avoid when interpreting machine learning models.  arXiv preprint arXiv:2007.04131, 2020.
  * [21] Judea Pearl.  Causality.  Cambridge university press, 2009.
  * [22] Judea Pearl and Azaria Paz.  Graphoids: A graph-based logic for reasoning about relevance relations.  University of California (Los Angeles). Computer Science Department, 1985\.
  * [23] Yaniv Romano, Matteo Sesia, and Emmanuel Candès.  Deep knockoffs.  Journal of the American Statistical Association, pages 1–12, 2019\.  Publisher: Taylor & Francis.
  * [24] Rajen D Shah and Jonas Peters.  The hardness of conditional independence testing and the generalised covariance measure.  arXiv preprint arXiv:1804.07203, 2018.
  * [25] Carolin Strobl, Anne Laure Boulesteix, Thomas Kneib, Thomas Augustin, and Achim Zeileis.  Conditional variable importance for random forests.  BMC Bioinformatics, 9:1–11, 2008.
  * [26] Laura Toloşi and Thomas Lengauer.  Classification with correlated features: unreliability of feature ranking and solutions.  Bioinformatics, 27(14):1986–1994, 2011.  Publisher: Oxford University Press.
  * [27] Eric J Topol.  High performance medicine: the convergence of human and artificial intelligence.  Nature Medicine, 25(January), 2019.  Publisher: Springer US.
  * [28] Eugene Tuv, Alexander Borisov, George Runger, and Kari Torkkola.  Feature selection with ensembles, artificial variables, and redundancy elimination.  Journal of Machine Learning Research, 10(Jul):1341–1366, 2009.
  * [29] Paul Voigt and Axel dem Bussche.  The eu general data protection regulation (gdpr).  A Practical Guide, 1st Ed., Cham: Springer International Publishing, 2017.  Publisher: Springer.
  * [30] David S. Watson and Marvin N. Wright.  Testing Conditional Independence in Supervised Learning Algorithms.  arXiv preprint arXiv:1901.09917, 2019.  arXiv: 1901.09917.
  * [31] Pengfei Wei, Zhenzhou Lu, and Jingwen Song.  Variable importance analysis: a comprehensive review.  Reliability Engineering & System Safety, 142:399–432, 2015.  Publisher: Elsevier.
  * [32] Yufei Xia, Chuanzhe Liu, Yu Ying Li, and Nana Liu.  A boosted decision tree approach using Bayesian hyper-parameter optimization for credit scoring.  Expert Systems with Applications, 78:225–241, 2017.  Publisher: Elsevier Ltd.
  * [33] Erik Štrumbelj and Igor Kononenko.  Explaining prediction models and individual predictions with feature contributions.  Knowledge and information systems, 41(3):647–665, 2014.  Publisher: Springer.


# Learning on Health Fairness and Environmental Justice via Interactive
Visualization

Abdullah-Al-Raihan Nayeem1, Ignacio Segovia-Dominguez2,3, Huikyo Lee2, Dongyun
Han4,
Yuzhou Chen5, Zhiwei Zhen3, Yulia Gel3, and Isaac Cho1,4
1Department of Computer Science, University of North Carolina at Charlotte,
Charlotte, United States
2Jet Propulsion Laboratory, California Institute of Technology, Pasadena,
United States
3Department of Mathematical Sciences, The University of Texas at Dallas,
Dallas, United States
4Department of Computer Science, Utah State University, Logan, United States
5Department of Computer and Information Sciences, Temple University,
Philadelphia, United States

###### Abstract

This paper introduces an interactive visualization interface with a machine
learning consensus analysis that enables the researchers to explore the impact
of atmospheric and socioeconomic factors on COVID-19 clinical severity by
employing multiple Recurrent Graph Neural Networks. We designed and
implemented a visualization interface that leverages coordinated multi-views
to support exploratory and predictive analysis of hospitalizations and other
socio-geographic variables at multiple dimensions, simultaneously. By
harnessing the strength of geometric deep learning, we build a consensus
machine learning model to include knowledge from county-level records and
investigate the complex interrelationships between global infectious disease,
environment, and social justice. Additionally, we make use of unique NASA
satellite-based observations which are not broadly used in the context of
climate justice applications. Our current interactive interface focus on three
US states (California, Pennsylvania, and Texas) to demonstrate its scientific
value and presented three case studies to make qualitative evaluations.

###### Index Terms:

coronavirus severity, interactive visualization, multivariate visualization

##  I Introduction

Transmissibility of various respiratory infectious diseases can be affected by
weather conditions and atmospheric factors, which, in turn, have an impact on
hospitalization rates. Particularly, in the context of the coronavirus disease
2019 (COVID-19) pandemic and its spread, there is no consensus on the impact
of atmospheric factors on the disease dynamic and the type of relationship
among variables [1, 2, 3]. Modeling and understanding the interlink between
atmospheric variables and the disease spread is cumbersome while considering
other factors (e.g., socioeconomic variables) also involved during infection
spread. Using state-of-the-art instruments on NASA’s Earth Observing System
(EOS) satellites, we make the first step forward to fill in the
interdisciplinary gap among exceptional NASA’s satellite-based datasets,
COVID-19 hospitalization dynamics, and interactive visualization systems.

The COVID-19 pandemic provides a unique opportunity to study environmental
(in)justice from both global and local perspectives, including novel
dimensions of social problems that keep arising as new data and historic
information comes out [4]. COVID-19 has unveiled the racial disparities in
healthcare access and health outcomes and taught us that health fairness and
environmental (in)justice goes beyond local hazards and exposure to pollutants
[5]. Using modern tools from artificial intelligence (AI), we can potentially
find novel elements that could lead us to previously missed analysis on the
role of climate in contemporary environmental (in)justice.

In this paper, we use multiple Recurrent Graph Neural Networks (RGNNs) to
forecast COVID-19 clinical severity with varying atmospheric factors. By using
RGNNs, we rely on one of the most effective state-of-the-art methods to model
data in non-Euclidean spaces, which are part of the broader field of Geometric
Deep Learning (GDL) [6]. Given the irregular geographic structure of available
data for COVID-19, and the noisy official records, GDL using the key
spatiotemporal patterns in satellite observations as predictors appear to be
one of the most promising forecasting approaches for tracking the hidden
mechanisms behind spatiotemporal COVID-19 dynamics. Additionally, for large-
scale multivariate spatiotemporal data such as our case, interactive
visualizations are essential to rapidly identifying critical association
patterns. The following are the key contributions of this paper:

  * •

We introduce a voting system based on multiple RGNNs to account for
connectivity among US counties, which naturally introduces a transmission
network of the disease at each learning stage. Our experiments take into
account both the temporal and spatial dependencies simultaneously with an
emphasis on forecasting the number of hospitalizations due to COVID-19.

  * •

We provide interactive visualization that further enhances the finding of
relationships between COVID-19 clinical severity and three key atmospheric
factors (temperature, relative humidity, and argon oxygen decarburization
(AOD)) obtained through satellite observations.

  * •

We present three usage scenarios from the interactive coordinated
visualizations to discover underlying insights within the cross-domain
multivariate datasets.

##  II Related Work

###  II-A Deep Learning for Spatiotemporal Forecasting

A natural deep learning solution to sequential data modeling consists of
Recurrent Neural Networks (RNNs) and their successors [7]. GNNs effectively
incorporate spatial dependencies in the RNN framework, especially in traffic
prediction problems on road network graphs [8]. Diffusion Convolutional RNN
(DCRNN) [9] integrates diffusion graph convolutional layers into a GRU network
for long-term traffic forecasting. In turn, Guo et al. [10] further employ an
attention mechanism to learn spatial and temporal dependencies, while Wu et
al. [11] developed a self-adaptive adjacency matrix to perform graph
convolutions without a predefined graph structure. Space-Time-Separable Graph
Convolutional Network (STSGCN) [12] simultaneously captures the localized
spatial-temporal correlations by localized spatial-temporal graphs.

###  II-B Bio-surveillance on COVID-19

Since the outbreak of COVID-19 in December 2019, there have been countless
attempts to model the spatiotemporal dynamics with particular emphasis on
forecasting disease spread and mortality. Mathematical epidemiologists study
the efficacy of contemporary models [13], and compare different intervention
strategies [14], whilst other researchers focus on mixed data-driven and
statistical approaches such as Autoregressive integrated moving average
(ARIMA) models [15], neural networks [16], RNNs [17], and Fuzzy Logic in
combination with Particle Swarm Optimization [18]. However, modeling dynamics
of COVID-19 clinical severity through the perspective of recorded
hospitalization rates remains under-explored [19, 20]. Our visualization
framework is a novel research effort to gain an understanding of the non-
linear dynamics of COVID-19 clinical severity in the presence of atmospheric
and socioeconomic factors.

###  II-C Visualization Methods

Numerous visualization systems are developed for COVID-19 monitoring [21, 22,
23] that provide the confirmed cases, active cases, recoveries, and deaths.
For instance, CoronaViz [24] is a visualization system with an interactive
view and a data control panel to unfold patterns in the spatiotemporal
COVID-19 variables. The data control in CoronaViz provides location and
feature-based filtering as well as an option to normalize based on population.
Geo-OEDV [25] is a visual framework for the exploratory analysis of clinical,
socioeconomic, and environmental data for COVID-19 targeting the general
users, policymakers, and researchers. Another visual analytics environment is
proposed for public health officials to perform exploration and decision-
making tasks with COVID-19 clinical data [26]. For multivariate data analysis,
CrossVis [27] demonstrated an interactive parallel coordinates plot consisting
of categorical and numerical axes to explore large-scale heterogeneous
spatiotemporal data. Multifaceted parallel coordinate plots are utilized in
other research to illustrate the association of spatiotemporal variables [28,
29].

To our best knowledge, none of the existing visualization systems supports
incorporating the cross-domain data for analyzing the impact of COVID-19
clinical severity factors. Hence, our visualization that accommodates
multivariate cross-domain data can benefit researchers to include variables
from fragmented data sets in forecasting the clinical severity.

##  III Datasets

Our research utilizes three real-world heterogeneous datasets: the COVID-19
dataset, the NASA dataset, and the socioeconomic dataset.

COVID-19 dataset.  Since our approach produces COVID-19 hospitalization
forecasts at the county-level resolution, in accordance with data
availability, we focus our experiments on three US states: California (CA),
Pennsylvania (PA), and Texas (TX). Daily records on COVID-19 and
hospitalizations are taken from curated datasets at the CovidActNow
project[30]. We follow the disease progression and guidelines from the Models
of Infectious Disease Agent Study (MIDAS) for COVID-19 modeling research [31].
We use time series, with values obtained on a daily basis, from February 1 to
December 31, 2020, thus focusing on the period of time when hospitalization
numbers were publicly available for all three states, at a county level.

NASA dataset.  To avoid current limitations on ground-based climate
observations, in terms of resolution on covered areas across the US, we rely
on processed satellite-based observations of temperature, relative humidity,
and AOD. NASA’s Distributed Active Archive Centers (DAAC) servers provide
publicly available access to the original datasets through the products: 1)
Atmospheric Infrared Sounder (AIR)/Aqua L3 Daily Standard Physical Retrieval
(AIRS3STD) [32] for temperature and relative humidity, and 2) the Atmosphere
Daily Global Product from Moderate-resolution Imaging Spectroradiometers
(MODIS) on Terra (MOD08_D3) [33] for AOD. Our climatology datasets (i.e.
temperature and relative humidity variables) come from averaging 6209 days
from January 1, 2003 to December 31, 2019. Thus, we compress annual cycles to
build the whole annual cycle of 2020, using the level 3 product (L3) data from
regular Gaussian grids with $1^{\circ}\times 1^{\circ}$ resolution. Similarly,
our AOD dataset is calculated using AOD at 550 $nm$ wavelength from the 6939
days from January 1, 2001 to December 31, 2019; thus producing an average
annual cycle. We match each county-level time series with its corresponding
Federal Information Processing Standard Publication 6-4 (FIPS 6-4) code. Thus,
our datasets are easier to use along with other datasets in the same county-
level resolution; e.g., socioeconomic indicators, infectious diseases,
transportation, etc.

County-level connectivity.  We analyze the impact of two types of graph
connectivity: 1) border-based county connections, and 2) socioeconomic-based
county connections. Network connectivity, using shared borders, is built via
the official county adjacency file record layout from the US Census Bureau. We
generate a socioeconomic-based network for each state under study using 9
socioeconomic variables. A socioeconomic matrix distance, via the euclidean
metric, helps us select the strongest county-level connections which serve as
input to each GNN model. Five variables come from the Centers for Disease
Control and Prevention (CDC)/Agency for Toxic Substances and Disease Registry
(ATSDR) Social Vulnerability Index [34]: Socioeconomic Status, Household
Composition & Disability, Minority Status & Language, Housing Type &
Transportation, and Overall Vulnerability Index. The rest of the variables are
part of The COVID-19 Vaccine Coverage Index (CVAC) [35]: Historic
Undervaccination, Sociodemographic Barriers, Resource-Constrained Healthcare
System, and Healthcare Accessibility Barriers.

##  IV Design Method

To evaluate the overall performance of DL/GDL models on spatio-temporal
forecasting tasks, in our experiment, we use 8 DL and GDL models to forecast
the next $\omega$ steps data based on past $h$ steps historical data, i.e.,
$\\{X_{t+1},X_{t+2},\dots,X_{t+\omega}\\}=\mathcal{M}(X_{t-h+1},X_{t-h+2},\dots,X_{t})$,
where $\mathcal{M}(\cdot)$ denotes DL/GDL model and $X_{i}$ is the node
features of graph $G_{t}$ at timestamp $i$ ($i\in\\{0,1,\dots,\tau\\}$).

Figure 1: The system pipeline for GDL models and interactive visualization
interface for county-level geographical model consensus of COVID-19 clinical
severity risk.

Benchmark models. We benchmark two broad classes of neural networks (i) RNNs:
LSTM can forecast multivariate time series with LSTM hidden units; (ii)
Spatiotemporal Graph Convolutional Networks: spatiotemporal model with the
framework of graph convolutional network (GCN) exploit GCN and temporal
convolution to capture dynamic spatial and temporal patterns and correlations;
we report performances for 8 types of state-of-the-art methods on our
benchmark datasets: LSTM, DCRNN, LSTM R-GCN, Attention Temporal GCN (A3T-GCN),
Message Passing Neural Networks with LSTM (MPNN+LSTM), Evolving GCNs
(EvolveGCNO), and Gated Graph Neural Network (GG-NNs) for Dynamic Graphs
(DyGrEncoder).

Experimental setting for forecasting. We use daily data from February 1 to
December 31, 2020, and split the graph signals into a train set (the first 80%
of days) and a test set (the last 20% of days). The training step uses 5 lags
of daily reported values to produce 15 days ahead forecasting. We train each
model using a hidden layer dimension of 128, and output dimension of
$\\{55,60,251\\}$ for CA, PA, and TX, respectively; as well as 500 epochs, a
dropout rate of 0.5, and an AMSGrad optimizer using 0.02 as a learning rate.
We execute 10 runs for each GDL model on each state to report prediction
results with statistical significance at both state and county levels.

GDL models for spatiotemporal forecasting. We present the recurrent graph
convolution operation on spatiotemporal networks in Algorithm 1. Given the
time-varying networks (i.e.,
${\bf{\boldsymbol{X}}}=\\{X_{0},\dots,X_{t},\dots,X_{\tau}\\}$ and
${\bf{\boldsymbol{G}}}=\\{G_{0},\dots,G_{t},\dots,G_{\tau}\\}$), for each
timestamp $t$, we first operate graph convolution on each node $v$ and obtain
graph embedding $H_{t}$, then feed the output embedding into RNN layer (e.g.,
Gated Recurrent Units (GRU)) for temporal dynamics modeling. In Algorithm 1,
$f_{\textnormal{MLP}}$, $\textnormal{AGG}(\cdot)$,
$\textnormal{CONCAT}(\cdot)$, $\phi(\cdot)$ denote the multilayer perceptron,
aggregation function (e.g., average, sum, max), concatenation operation, and
nonlinear activation function (e.g., ReLU) respectively; $W_{g}$, $W_{z}$,
$W_{r}$, $W_{q}$, $b_{z}$, $b_{r}$, $b_{q}$ are trainable weights, and $\odot$
denotes the element-wise product.

Input: ${\bf{\boldsymbol{X}}}=\\{X_{0},\dots,X_{T}\\}$ and
${\bf{\boldsymbol{G}}}=\\{G_{0},\dots,G_{\tau}\\}$

Output:
${\bf{\boldsymbol{\tilde{Q}}}}_{t}=\\{\tilde{Q}_{1},\tilde{Q}_{2},\dots,\tilde{Q}_{t},\dots,\tilde{Q}_{\tau}\\}$

/* Graph convolution operation on $\\{G_{t},X_{t}\\}$ */

1 for  _the graph $\\{G_{t},X_{t}\\}$ at timestamp $t$ in mini-batch_ do

2 $h_{t}^{(0)}=f\textsubscript{MLP}(X_{t})$

3 for  _$\ell\leftarrow 1$ to _nlayers__ do

4 $h^{(\ell-1)}_{t,\mathcal{N}(v)}=\textnormal{\small
AGG}\left(h^{(\ell-1)}_{t,u},\forall u\in\mathcal{N}(v)\right)$

5 $h^{(\ell)}_{t,v}=\sigma\left(W_{g}\cdot\textnormal{\small
CONCAT}\left(h^{(\ell-1)}_{t,\mathcal{N}(v)},h^{(\ell-1)}_{t,v}\right)\right)$

6

7 end for

8

9 end for

$H_{t}=\textnormal{Readout}(\\{h^{(\ell)}_{t,v},\forall v\in\mathcal{V}\\})$
/* Embedding */

/* Recurrent neural network layer */

10 $Z_{t}=\phi(W_{z}[H_{t-1},H_{t}]+b_{z})$

11 $R_{t}=\phi(W_{r}[H_{t-1},H_{t}]+b_{r})$

12 $\hat{Q}_{t}=\tanh{(W_{q}[R_{t}\odot\hat{Q}_{t-1},H_{t}])+b_{q}}$

13 $\tilde{Q}_{t}=Z_{t}\odot\hat{Q}_{t-1}+(1-Z_{t})\odot\hat{Q}_{t}$

Algorithm 1 Recurrent graph convolution operation

Consensus analysis. We compute forecasting using 8 state-of-the-art models.
Our consensus analysis consists in adding votes, at a county level, whenever
an RNN/GNN model improves its results (when adding an atmospheric variable)
with respect to the baseline (removing the atmospheric variable). Each
improvement is validated through one-tailed hypothesis testing, at significant
level, on reduction of the root-mean-square error (RMSE). To the best of our
knowledge, there are no previous attempts to use consensus among multiple
RNN/GNN models to characterize the impact of atmospheric factors on COVID-19
clinical severity, along with an interactive visualization system to
complement the decision-making of experts.

Analysis Pipeline. Fig. 1 illustrates our data analysis workflow. First, we
collect and pre-process the spatiotemporal data, both COVID-19 hospitalization
numbers, and satellite-based observations of temperature, relative humidity,
and AOD. Next, we generate two types of county-level connectivity: 1)
geographic-based border county connections, and 2) socioeconomic-based county
connections. Using a chosen graph connectivity and dynamic measurements at
each node, we build a dynamic node-feature network
$\\{G_{0},...,G_{\tau-1},G_{\tau}\\}$ for the GNNs models. Our platform allows
us to analyze results using both types of networks by using these as input in
our GDL models. We train a set of models
$\\{M_{0},...,M_{\tau-1},M_{\tau}\\}$, as described in Section IV, to
forecasting hospitalization numbers with and without atmospheric factors.
Next, we statistically verify, and count, any model that produces an
improvement in forecasting results when adding atmospheric factors. This
consensus model produces geographical patterns which help us identify
communities with a higher risk of hospitalization due to COVID-19.

Note that this same framework can be used to study other infectious diseases,
such as Zika or Monkeypox, which in turn will impact government policies and
health decision-making.

Figure 2: Overview of the interactive visualization interface. A) Parallel
coordinates presenting a multivariate visualization of clinical, atmospheric,
and socioeconomic variables of the counties. B) A county network illustrates
the relationships among the counties based on socioeconomic similarities. C)
Choropleth map views presenting the geographical representation of the
counties with the summary of participating variables and the GDL voting
result. A control bar at the top supports reconfiguring these interactive
visualizations.

##  V Interactive Visualizations

The visual interface presents interactive and coordinated visualizations to
facilitate the user for exploring the clinical, atmospheric, and socioeconomic
variables as well as the voting results from GDL models. To support the
interactive analysis of both input and resulting variables, the interface
consists of four main views: A) Multivariate Visualization to explore the
input data variables (e.g., death rate, hospitalization rate, socioeconomic
status), B) County Network View to explore the connectivity among the counties
of a specific US state, C) Choropleth Map View to provide state choropleth
maps to explore the summarized data variables for each county by its
geographic location, and D) Voting Results View to explore the GDL results to
determine the impact of socioeconomic variables on COVID-19 clinical severity
factors. The interactive visualizations are served via a web-interface.111The
visual interface is available online: https://vizus.cs.usu.edu/app/gdlviz/

Multivariate Visualization. To examine the impact factors, we categorize the
input variables between two groups: primary and impact variables. The primary
variables are death and hospitalization rates. Since the county population has
a substantial difference among others, we normalized the primary data
variables to find the rate per 100K population. Moreover, the impact variables
have the five variables from the CDC/ATSDR social vulnerability index. To
support the investigation of correlation among the primary and impact data
variables, we leveraged parallel coordinates that are widely used
visualization for multivariate data [27].

The dimensions in the parallel coordinates combine the primary and impact
variables. Each axis in the parallel coordinates represents a dimension. Since
we are investigating the impact variables against the primary variables, the
axes for the primary variables are highlighted in orange color (Fig. 2A). Each
connecting line in the parallel coordinates represents a county of the
selected state. Hence, the total number of connecting lines in the parallel
coordinates matches the number of counties in the state. Each county line
connects the dimensions based on its corresponding value. The parallel
coordinates provide interactive data filtering to specify the window size for
each axis (Fig. 2G). The filtered county lines are highlighted in blue color
while other lines are in gray color.

An axis for each impact variable includes two mini scatter plots on the bottom
(Fig. 2H). The left and right scatter plots depict a glimpse of the
corresponding correlation with the death and hospitalization rates
respectively. The interface provides a detailed view of each scatter plot
including a trend line for the selected variables on clicking a mini scatter
plot.

County Network View. This view presents each county as a node and the
connections between the two counties are represented with an edge. The length
of the edge denotes the euclidean distance (i.e., socioeconomic similarity)
between the two connected counties. County nodes within the network can be
focused or selected with mouse interaction. Hovering a mouse cursor on a
county node shows all neighboring counties in the network as shown in Fig. 4A.
Moreover, we calculate network centrality to find the significant counties in
the network. Network centrality provides a tool to rate/rank the nodes in the
network based on the number of degrees (edges) [36].

The view has a control bar to change the distance threshold, network
centrality threshold, network strength, etc. The distance threshold filters
out the connections that are far from the selected value. Therefore, the
higher value of the threshold provides a dense-edge network as shown in Fig.
4F. The network centrality threshold provides control to the user to determine
the rank where a county is to be considered statistically significant.
Initially, all the counties are marked blue whereas the important counties
based on the network centrality threshold are marked green (Fig. 5).

Choropleth Map View. This view is designed to present a state’s county-level
information on independent and impact variables as well as to depict the
patterns based on geographical features. We utilized different color scales
associated with the data variables to illustrate the value for each county.
The view has two choropleth maps (Fig. 2C and D) to provide a comparative
capability to the user. These reconfigurable choropleth maps are labeled as
left and right choropleth where the user can select their variable of
interest.

The data presented in the choropleth map are categorical, regardless of
independent, impact, or resulting variables. For the continuous variables such
as hospitalization rate per 100k and relative humidity, we calculated the
quantiles to categorize the counties. For example, Fig. 4F suggests that
counties that lie in the top 80% intensity of relative humidity are located on
the east side of TX. Choropleth maps are also employed to visualize the
summary of the voting results discussed in Section IV. Since we used 8 GDL
models, each county can get 0 to 8 votes from the model output. When the user
hovers over a county, this view provides additional information on a tooltip
such as the county name with code, values for the corresponding atmospheric
and socioeconomic variables, and a breakdown of the forecasts predicted by the
individual deep-learning models.

A histogram is implemented below the choropleth map to summarize the
categorical data presented on the map. The bars in the histogram used the same
color scale used for the choropleth map. The histogram presents the number of
counties for each category (e.g. vote count, and quantile) as shown in Fig.
4D-F. Hovering a mouse cursor over a bar highlights the corresponding counties
on the choropleth map.

Figure 3: Pixel bar view provides a detailed illustration of the voting
results from GDL models. Interactive sorting and filtering strategies are
implemented to support the exploration and logical reasoning of counties’ vote
counts for individual states.

Pixel Bar View. A compound view that consists of two visualizations combined
with interactive features - a) a bar chart containing the vote count for the
counties (Fig. 3A) where each bar represents a county in the selected state,
b) followed by the pixel bar view that provides a breakdown of the voting
results (Fig. 3E). Both visualizations share the same x-axis and are
vertically aligned, combining 8 rows, each representing a GDL model. The pixel
bar view marks the county’s associated pixel with blue color if the
corresponding model results in a positive prediction.

This view provides interactive features to allow the user to reconfigure the
rendered visualizations based on data sorting and filtering. We implemented a
couple of sorting strategies to support the user in identifying patterns. The
sorting strategies are applied by the participating independent variable
(e.g., hospitalization rate/100k), impact variable (e.g., AOD), and vote count
for a county (Fig. 3). The sorting determines the order of the counties
presented in the pixel bar view. For example, Fig. 3B illustrates the pixel
bar sorted by the vote count where the counties with similar vote counts are
tied together. This feature provides an idea of how individual models
performed and whether there is any bias in the vote count. The user can also
filter the view while interacting with the bar chart at the top (Fig. 3F).
When the user selects a group of neighboring counties ordered in the bar
chart, the pixel bar view is reconfigured to present the selected counties.
Therefore, the user can focus on a specific group and identify the county
names interacting with the pixel bar (Fig. 3G). The selection works as a
sliding window where the view allows the user to shift left/right and change
the window size.

Figure 4: The interactive interface demonstrates exploratory user interactions
that coordinate multiple views facilitating visual sense-making.

User Interactions. We implemented user interactions based on data and view
manipulation capabilities to facilitate the human cognition process. The user
interactions are categorized based on their support for different visual
analysis tasks [37]. In our visual interface, we leveraged select, explore,
filter, reconfigure, and connect interactions.

The user can switch a state (TX, PA, or CA) from the control bar. In addition,
the user can specify by interacting with the scatterplots on the multivariate
parallel coordinates. For example, selecting a scatterplot presenting the
hospitalization rate and AOD (Fig. 2H) reconfigures the context for the
choropleth views. Selection is also utilized for specifying the independent
and impact variable while interacting with the scatterplots on the
multivariate parallel coordinates. For example, selecting a scatterplot
presenting the hospitalization rate and AOD (Fig. 2H) reconfigures the context
for the choropleth views. The user can switch between pre-computed model
results based on the socioeconomic distance and boundary-based output of the
GDL models. Moreover, hovering a mouse over a county on the vote-count
choropleth map view shows the GDL models that contribute to that count.
Interacting with a node in the county network also highlights the neighboring
counties (Fig. 4) based on the selected distance threshold. Filtering
interaction allows the user to focus on a specific subset or area of interest.
The parallel coordinates presenting multivariate data support filter
interaction for each data axis as shown in Fig. 2A. The county network view
also provides a group select option where the neighboring counties can be
identified as a cluster on the network (Fig. 2E). Previously, we also
discussed the filtering features for the pixel bar view and choropleth map
view interacting with the histogram.

The visual interface demonstrates a connected and coordinated response to the
user’s filter interactions. For example, in Fig. 2A, the filter applied on the
multivariate parallel coordinates essentially triggered a filtering effect on
the other visual components such as Fig. 2B-C-D. Similarly, interacting with a
node on the network or a bar chart imposes a coordinated effect on all other
views as shown in Fig. 4.

##  VI Usage Scenarios

We present three usage scenarios of our proposed visual interface. For each
usage scenario, we focus on presenting the coordinated visualization’s
capabilities to support exploratory analysis and provide logical reasoning to
the discovered knowledge to demonstrate the scalability and adaptability of
our interface.

Figure 5: Coordinated multi-views demonstrating the interactive features in
sense-making the GDL model result for TX and PA when associated with AOD and
relative humidity.

Improved forecasting in TX due to AOD:  In Fig. 5A, the choropleth map at the
county level of TX illustrates North Eastern and Great Houston regions are
better predicted, by multiple DL models at a $\alpha=0.1$ significant level,
when introducing measurements of particles of the atmosphere. This fact can be
associated with a higher amount of particles in the air of these counties when
contrasted with their average rates. As an example, the on-road air pollution
of Houston is 9 times greater than its metro-area counterparts [38].
Additionally, a view of the Overall Vulnerability Index [39] shows a
significant relationship between counties with the highest impact of air
quality on hospitalizations and those socioeconomically vulnerable counties as
shown in Fig. 5B. Further uncovering the structural inequalities and
disparities, in health access and health outcomes, of socioeconomically
disadvantaged communities (Fig. 5C).

Linked racial disparities with clinical severity in PA: Forecasting of
COVID-19 hospitalizations in Southeastern PA gets improved when associated
with daily rates of relative humidity. This further supports the health impact
of pollution on these communities when linked with the fact that PA ranks as
one of the 25 worst US metro areas for ozone and year-round particle pollution
[40]. There is a further significant improvement, in the total number of
votes, when adding dense socioeconomic-based counties’ connectivity (223)
against the boundary-based network (193). Even more, a visual comparison
between our consensus model and the Minority Status & Language index shows a
compelling relationship between racial disparities and COVID-19 clinical
severity. Thus showing the important role that environmental (in)justice is
playing to prevent COVID-19 clinical severity in PA, and across communities
around the world [4].

Figure 6: Comparing the GDL voting results where counties are grouped by vote
counts for TX, PA, and CA to learn the impact of AOD, relative humidity, and
temperature in COVID-19 clinical severity.

State level impact of atmospheric factors: The statistically verified votes in
our consensus model allow us to identify regions where specific atmospheric
factors could explain COVID-19 clinical severity. At a state level, we can
conclude that AOD plays a bigger role than relative humidity and temperature
when linked to forecasting hospitalization numbers, conditioned to multiple
socioeconomic indices. Fig. 6 presents the total number and distribution of
votes in the three states when adding each atmospheric variable AOD (TX: 1079,
PA: 236, CA: 73), relative humidity (TX: 409, PA: 145, CA: 23), and
temperature (TX: 421, PA: 226, CA: 21). These findings further support the
current trend of rethinking our engagement with the natural world, mainly
because of our role as generators of particles in the air, but also from a
wider environmental (in)justice perspective due to the impact in already
vulnerable communities [4].

###  VI-A Discussion and Limitations

The proposed visualization interface offers a number of multifold perspectives
for health analytics way beyond the focus of COVID-19 disease tracking. First,
the developed interface allows public health practitioners to perform
retrospective spatiotemporal analysis of potential atmospheric contributors to
clinical severity and mortality associated with respiratory diseases at
various time periods and scales, while systematically accounting for socio-
demographic information and quantifying uncertainties. As a result, such
analysis can shed light on, for example, the potential (dis)similarities
between the emerging and earlier recorded variants of a virus and the
associated sensitivity of these variants to weather and climate. Second,
coupled with weather forecasting, the visualized consensus among GDL models
can be used to identify areas of the highest vulnerability risks and potential
next hotspots requiring proactive medical emergency readiness. Third, more
generally, armed with the new visualization interface, public policymakers can
investigate such questions of critical societal importance as climate justice
and their impacts on various disease progression stages.

While the interactive visualization interface provides coordinated
visualizations and rich user interactions that support primary and resulting
data exploration and facilitate sense-making, it relies on pre-computed model
results. Hence, the user is responsible to run their own machine-learning
workflow to visually analyze data in the proposed interface. In our future
work, we will address this concern by leveraging a visual analytics framework
where the user is able to create workflows, prepare analytical scripts, and
execute from the interface which would essentially end up utilizing the
proposed interactive visualizations [41]. A visual analytics system will also
facilitate the user to interactively experiment with the parameters and
statistical significance. Moreover, we plan to conduct a user study with
domain experts to report a quantitative evaluation of our proposed
visualization system.

##  VII Conclusion

The COVID-19 pandemic provides a unique opportunity to study the interlinks
between environment, climate (in)justice, and biosurveillance. From a
contemporary AI approach, this paper aims to add up to the collective
knowledge of the lessons learned to identify and tackle structural
inequalities contributing to higher mortality in vulnerable communities, and
its relationships with atmospheric factors. We introduce an interactive
interface and consensus analysis with coordinated multi-views that supports
analyzing the impacts of atmospheric observations and socioeconomic factors on
COVID-19 clinical severity. We present three usage scenarios from the proposed
interactive interface to demonstrate the usability and its scientific value in
analyzing COVID-19 clinical severity. Our findings support the impact of
aerosols on the Earth’s radiation budget and air quality, and its consequences
for health fairness across communities of different socioeconomic backgrounds.
Since our machine learning and interactive framework leverages NASA’s
satellite-based observations, our consensus analysis can be extended to a
wider coverage over the entire globe.

## Acknowledgement

This project is partially supported by the grants NASA 21-AIST21_2-0059 and
NSF Data Infrastructure Building Blocks (DIBBs) Program (Award #1640818).

## References

  * [1] M. Rendana and W. M. R. Idris, “New COVID-19 variant (B.1.1.7): Forecasting the occasion of virus and the related meteorological factors,” _Journal of Infection and Public Health_ , vol. 14, no. 10, pp. 1320–1327, 2021, special Issue on COVID-19 – Vaccine, Variants and New Waves.
  * [2] O. Iloanusi and A. Ross, “Leveraging weather data for forecasting cases-to-mortality rates due to COVID-19,” _Chaos, Solitons & Fractals_, vol. 152, p. 111340, 2021.
  * [3] X. Zhang, V. Maggioni, P. Houser, Y. Xue, and Y. Mei, “The impact of weather condition and social activity on COVID-19 transmission in the United States,” _Journal of Environmental Management_ , vol. 302, p. 114085, 2022\.
  * [4] S. M. Wilson, R. Bullard, J. Patterson, and S. B. Thomas, “Roundtable on the Pandemics of Racism, Environmental Injustice, and COVID-19 in America,” _Environmental Justice_ , vol. 13, no. 3, pp. 56–64, 2020.
  * [5] D. Cooper and J. Nagel, “Lessons from the pandemic: Climate change and covid-19,” _International Journal of Sociology and Social Policy_ , vol. ahead-of-print, 01 2021.
  * [6] M. M. Bronstein, J. Bruna, Y. LeCun, A. Szlam, and P. Vandergheynst, “Geometric deep learning: going beyond euclidean data,” _IEEE Signal Processing Magazine_ , vol. 34, no. 4, pp. 18–42, 2017.
  * [7] Y. Yu, X. Si, C. Hu, and J. Zhang, “A review of recurrent neural networks: LSTM cells and network architectures,” _Neural Computation_ , vol. 31, no. 7, pp. 1235–1270, 2019.
  * [8] B. Yu, H. Yin, and Z. Zhu, “Spatio-temporal graph convolutional networks: A deep learning framework for traffic forecasting,” in _IJCAI_. IJCAI Organization, jul 2018.
  * [9] Y. Li, R. Yu, C. Shahabi, and Y. Liu, “Diffusion convolutional recurrent neural network: Data-driven traffic forecasting,” in _Proceedings of the ICLR_ , 2018.
  * [10] S. Guo, Y. Lin, N. Feng, C. Song, and H. Wan, “Attention based spatial-temporal graph convolutional networks for traffic flow forecasting,” in _AAAI_ , vol. 33, 2019, pp. 922–929.
  * [11] Z. Wu, S. Pan, G. Long, J. Jiang, and C. Zhang, “Graph WaveNet for Deep Spatial-Temporal Graph Modeling,” in _Proceedings of the 28th IJCAI_. IJCAI Organization, aug 2019.
  * [12] C. Song, Y. Lin, S. Guo, and H. Wan, “Spatial-temporal synchronous graph convolutional networks: A new framework for spatial-temporal network data forecasting,” _AAAI_ , vol. 34, no. 01, pp. 914–921, apr 2020.
  * [13] A. Vespignani, H. Tian, C. Dye, J. O. Lloyd-Smith, R. M. Eggo, M. Shrestha, S. V. Scarpino, B. Gutierrez, M. U. Kraemer, J. Wu _et al._ , “Modelling COVID-19,” _Nature Reviews Physics_ , vol. 2, no. 6, pp. 279–281, 2020.
  * [14] S. L. Chang, N. Harding, C. Zachreson, O. M. Cliff, and M. Prokopenko, “Modelling transmission and control of the COVID-19 pandemic in Australia,” _Nature Communications_ , vol. 11, no. 1, pp. 1–13, 2020.
  * [15] M. K. Nesa, M. R. Babu, and M. T. Mamun Khan, “Forecasting COVID-19 situation in Bangladesh,” _Biosafety and Health_ , 2021.
  * [16] K. C. de Carvalho, J. P. Vicente, and J. P. Teixeira, “COVID-19 Time Series Forecasting – Twenty Days Ahead,” _Procedia Computer Science_ , vol. 196, pp. 1021–1027, 2022.
  * [17] K. ArunKumar, D. V. Kalaga, C. M. S. Kumar, M. Kawaji, and T. M. Brenza, “Forecasting of COVID-19 using deep layer Recurrent Neural Networks (RNNs) with Gated Recurrent Units (GRUs) and Long Short-Term Memory (LSTM) cells,” _Chaos, Solitons & Fractals_, vol. 146, p. 110861, 2021.
  * [18] N. Kumar and S. Susan, “Particle swarm optimization of partitions and fuzzy order for fuzzy time series forecasting of COVID-19,” _Applied Soft Computing_ , vol. 110, p. 107611, 2021.
  * [19] I. Segovia-Dominguez, Z. Zhen, R. Wagh, H. Lee, and Y. R. Gel, “TLife-LSTM: Forecasting Future COVID-19 Progression with Topological Signatures of Atmospheric Conditions.” in _PAKDD_. Springer, 2021, pp. 201–212.
  * [20] I. Segovia Dominguez, H. Lee, Y. Chen, M. Garay, K. M. Gorski, and Y. R. Gel, “Does Air Quality Really Impact COVID-19 Clinical Severity: Coupling NASA Satellite Datasets with Geometric Deep Learning,” in _ACM SIGKDD_ , 2021, pp. 3540–3548.
  * [21] “Johns Hopkins University, coronavirus covid-19 global cases,” https://coronavirus.jhu.edu/, Accessed: 2022-09-03.
  * [22] “World Health Organization, novel coronavirus (covid-19) situation dashboard,” https://covid19.who.int/, Accessed: 2022-09-03.
  * [23] “University of Washington, novel coronavirus infection map,” https://hgis.uw.edu/virus/, Accessed: 2022-09-03.
  * [24] H. Samet, Y. Han, J. Kastner, and H. Wei, “Using animation to visualize spatio-temporal varying covid-19 data,” in _ACM SIGSPATIAL International Workshop on Modeling and Understanding the Spread of COVID-19_. New York, NY, USA: ACM, 2020, p. 53–62.
  * [25] A. Bernasconi and S. Grandi, “A Conceptual Model for Geo-Online Exploratory Data Visualization: The Case of the COVID-19 Pandemic,” _Information_ , vol. 12, no. 2, 2021.
  * [26] S. Afzal, S. Ghani, H. C. Jenkins-Smith, D. S. Ebert, M. Hadwiger, and I. Hoteit, “A Visual Analytics Based Decision Making Environment for COVID-19 Modeling and Visualization,” in _IEEE VIS_ , 2020, pp. 86–90.
  * [27] C. A. Steed, J. R. Goodall, J. Chae, and A. Trofimov, “CrossVis: A visual analytics system for exploring heterogeneous multivariate data with applications to materials and climate sciences,” _Graphics and Visual Computing_ , vol. 3, p. 200013, 2020.
  * [28] F. Wang, W. Li, S. Wang, and C. R. Johnson, “Association rules-based multivariate analysis and visualization of spatiotemporal climate data,” _ISPRS International Journal of Geo-Information_ , vol. 7, no. 7, 2018.
  * [29] Z. Li, Q. Huang, Y. Jiang, and F. Hu, “Sovas: a scalable online visual analytic system for big climate data analysis,” _International Journal of Geographical Information Science_ , vol. 34, no. 6, pp. 1188–1209, 2020.
  * [30] “U.S. COVID Tracker,” https://covidactnow.org/, Accessed: 2022-09-03.
  * [31] “MIDAS: Online Portal for COVID-19 Modeling Research,” https://midasnetwork.us/covid-19/, Accessed: 2022-09-03.
  * [32] “AIRS/Aqua L3 Daily Standard Physical Retrieval (AIRS-only), Goddard Earth Sciences Data and Information Services Center (GES DISC),” https://doi.org/10.5067/Aqua/AIRS/DATA303, Accessed: 2022-09-03.
  * [33] “MODIS Level 3 Atmosphere Gridded Product,” https://modis.gsfc.nasa.gov/data/dataprod/mod08.php, Accessed: 2022-09-03.
  * [34] CDC, “CDC SVI 2018 Documentation,” https://www.atsdr.cdc.gov/placeandhealth/svi/documentation/pdf/SVI2018Documentation-H.pdf, Accessed: 2022-09-03.
  * [35] S. Ventures, “The U.S. COVID-19 Vaccine Coverage Index (CVAC) Methodology,” https://cvi-data-output.s3.amazonaws.com/assets/CVAC_Methodology_Feb2021.pdf, Accessed: 2022-09-03.
  * [36] H. Ibarra, “Network centrality, power, and innovation involvement: Determinants of technical and administrative roles,” _Academy of Management journal_ , vol. 36, no. 3, pp. 471–501, 1993.
  * [37] Y. Lu, R. Garcia, B. Hansen, M. Gleicher, and R. Maciejewski, “The state-of-the-art in predictive visual analytics,” in _Computer Graphics Forum_ , vol. 36, no. 3. Wiley Online Library, 2017, pp. 539–562.
  * [38] ANT, “Air north texas,” https://www.airnorthtexas.org/ozone, Accessed: 2022-09-03.
  * [39] CDC-ATSDR, “Cdc/atsdr social vulnerability index,” https://www.atsdr.cdc.gov/placeandhealth/svi/documentation/pdf/SVI2018Documentation-H.pdf/, Accessed: 2022-09-03.
  * [40] ALA, “American Lung Association. https://www.lung.org/media/press-releases/air-quality-in-philadelphia,” 2022\.
  * [41] A.-A.-R. Nayeem, M. Elshambakey, T. Dobbs, H. Lee, D. Crichton, Y. Zhu, C. Chokwitthaya, W. J. Tolone, and I. Cho, “A visual analytics framework for distributed data analysis systems,” in _IEEE Big Data_ , 2021, pp. 229–240.