# Testing Causality in Scientific Modelling Software

Andrew G. Clark
<EMAIL_ADDRESS>Michael Foster
<EMAIL_ADDRESS>Benedikt Prifling
<EMAIL_ADDRESS>Neil Walkinshaw
<EMAIL_ADDRESS>Robert M. Hierons
<EMAIL_ADDRESS>Volker Schmidt
<EMAIL_ADDRESS>Robert D. Turner
<EMAIL_ADDRESS>
###### Abstract

From simulating galaxy formation to viral transmission in a pandemic,
scientific models play a pivotal role in developing scientific theories and
supporting government policy decisions that affect us all. Given these
critical applications, a poor modelling assumption or bug could have far-
reaching consequences. However, scientific models possess several properties
that make them notoriously difficult to test, including a complex input space,
long execution times, and non-determinism, rendering existing testing
techniques impractical. In fields such as epidemiology, where researchers seek
answers to challenging causal questions, a statistical methodology known as
Causal Inference has addressed similar problems, enabling the inference of
causal conclusions from noisy, biased, and sparse data instead of costly
experiments. This paper introduces the Causal Testing Framework: a framework
that uses Causal Inference techniques to establish causal effects from
existing data, enabling users to conduct software testing activities
concerning the effect of a change, such as Metamorphic Testing and Sensitivity
Analysis, _a posteriori_. We present three case studies covering real-world
scientific models, demonstrating how the Causal Testing Framework can infer
test outcomes from reused, confounded test data to provide an efficient
solution for testing scientific modelling software.

##  1 Introduction

The use of scientific modelling software to model, simulate, and understand
complex phenomena has become commonplace. Such systems have played a pivotal
role in improving our scientific understanding across a wide range of
phenomena and disciplines, and are increasingly used outside of academia.
Governments, for example, make extensive use of scientific modelling software
to simulate and evaluate various policies and interventions [1]. Perhaps most
notably, this has included the use of epidemiological models to predict the
impact of a number of COVID-19 mitigation measures [2, 3].

Testing such models is particularly challenging [4]. They typically have vast
input spaces comprising hundreds of parameters, as well as complex output
spaces. Executing large numbers of tests is often impossible, because each
execution can require a significant amount of time and resource to execute.
Compounding this issue further, scientific models are often non-deterministic,
meaning developers must run each test case multiple times and observe the
distribution of outputs. Furthermore, these systems are often developed by
scientists with a limited amount of training as software engineers [5].

Collectively, these issues make it difficult (and sometimes impossible) to
determine whether the output of a test case or modelling scenario is correct
or not. This is referred to as the test oracle problem [6]. Instead, to
determine whether a software system is fit for purpose, a tester generally
corroborates evidence to investigate smaller, more specific relationships
between inputs and outputs. By making changes to particular input parameters
and observing changes to particular output variables, there is an implicit
assumption that the input parameters in question somehow influence the
computation (i.e. have a ‘causal’ effect) of the outputs.

In this paper we are specifically concerned with this intrinsic challenge: How
can we test the (implicitly causal) input-output relationships in a system
with a vast and complex input space, which may be non-deterministic and suffer
from the test oracle problem, without the ability to resort to large numbers
of test executions?

The challenge of analysing causal relationships in limited, noisy data instead
of running costly experiments is well-established in the statistical context.
In areas such as epidemiology, for example, causal questions that cannot be
answered experimentally due to ethical concerns, such as _Does smoking cause
lung cancer?_ [7], have been answered using a powerful statistical methodology
known as causal inference (CI). By incorporating domain knowledge about known
causal relationships between variables (or absence thereof), CI can produce
_estimands_ that isolate the causal quantity of interest. That is, ‘recipes’
for analysing data in a causally-valid way. Conventional statistical methods
can then be employed to quantify the presence (or absence) of specific causal
relationships, correcting for bias in the data, without the need for
experimental procedures.

This paper is motivated by the observation that CI and software testing share
a common goal in many cases: to establish precise and salient causal
relationships. Moreover, by viewing software testing through a causal lens, we
can leverage well-established CI techniques that conceptually address several
testing challenges presented by scientific models for causality-driven testing
activities, such as metamorphic testing and sensitivity analysis.

To this end, we introduce a testing framework that incorporates an explicit
model of causality into the testing process, facilitating the direct
application of CI techniques to a myriad of software testing problems. To
achieve this, we take a model-based testing (MBT) perspective [8], in which
testing is based on a model of the expected behaviour of the system-under-test
that typically either describes the allowed sequences of events or gives a
formal relation between the inputs and outputs [9, 10]. Traditionally, MBT has
focused on models expressed using state-based languages, such as finite state
machines [11] and labelled transition systems [12], or models that define the
allowed input-output relationships using languages, such as Z [13] and VDM
[14]. However, given the focus on causality in this paper, we require a model
that specifies the expected _causal relationships_ between system inputs and
outputs. Here, we assume that such causal information is represented by a
causal directed acyclic graph (DAG) [15, 16].

Our decision to incorporate causal DAGs into the testing process is motivated
by two main factors. First, testing can be viewed as a causal activity in
which the tester checks whether expected causal relationships hold; in order
to automate this process, we require the expected causal relationships to be
expressed. Second, the causal DAG is a lightweight and intuitive model that is
widely used by domain experts in areas such as epidemiology and social
sciences to make causal assumptions actionable and transparent [17, 18].

In this paper, we make three contributions. First, we introduce a conceptual
framework that approaches causality-driven software testing activities as CI
problems, and clarifies the components necessary to leverage state-of-the-art
CI techniques. Second, we provide a reference implementation of the framework
that can form the basis for future CI-driven tools for testing scientific
modelling software. Third, we conduct three case studies covering real-world
scientific models from different domains to demonstrate the advantages of a
CI-led approach to testing scientific modelling software.

The remainder of this paper is structured as follows. Section 2 provides a
motivating example and necessary background. Section 3 introduces our
conceptual framework that frames causality-driven testing activities as
problems of CI. Section 4 then introduces our reference implementation of this
framework, before demonstrating its application to three real-world scientific
models in Section 5 and discussing the main findings and threats to validity
in Section 6. Section 7 reviews related work, and Section 8 concludes the
paper.

##  2 Background and Preliminaries

This section defines the scope of the paper and introduces the main challenges
associated with testing scientific modelling software, as outlined in Kanewala
and Bieman’s survey on the same topic [4]. We present these challenges in the
context of a real-world, motivating example that is used as one of three case
studies in Section 5. We then provide a background on metamorphic testing
[19], a known solution to some of these challenges, and a brief introduction
to causal inference, the statistical methodology employed by the framework
presented in Section 3.

###  2.1 Black-Box Software Systems

In this paper, we view and test software from a black-box perspective [20],
focusing on the relationships between its inputs and outputs rather than its
inner-workings and source code. More formally, in this paper, we conceptualise
the system-under-test (SUT) as follows:

######  Definition 2.1.

A _system-under-test_ (SUT) is a software system comprising a set of input
variables, $I$, and output variables, $O$, such that $I\cap O=\oslash$. We
consider inputs to be parameters whose values are set prior to execution that
influence the resulting system behaviour. We consider outputs to be features
of the system that can be measured at any point during or after execution
without inspecting or modifying the source code.

Given our focus on causality in this paper, we provide an informal definition
of causality in Definition 2.2. This follows from Pearl’s characterisation of
causation, which states that “variables earn causal character through their
capacity to sense and respond to changes in other variables” [21].

######  Definition 2.2.

We say that a variable $X=x$ _causes_ a variable $Y$ if there exists some
value $x^{\prime}$ such that, had the value of $X$ been changed to
$x^{\prime}$, the value of $Y$ would change in response.

Furthermore, we are primarily interested in scientific modelling software.
Informally, we consider this to be any form of software that has a significant
computational component and simulates, models, or predicts the behaviour of
complex, uncertain phenomena to support policy and scientific decisions [4,
22]. We focus on this form of software as it typically possesses a number of
challenging characteristics that preclude the application of many conventional
testing techniques, but can be addressed by the framework introduced in
Section 3. In the following section, we introduce a motivating example to
familiarise the reader with these challenging properties.

###  2.2 Motivating Example: Covasim

Covasim [3, 23] is an epidemiological agent-based model that has been used to
inform COVID-19 policy decisions in several countries [24, 25, 26, 27]. Given
the critical applications of such scientific models, it is of paramount
importance that they are tested to the best of our abilities. However, Covasim
has a number of characteristics that make testing particularly challenging.

Covasim has a vast and complex input space, with 64 unique input parameters,
27 of which are complex objects characterised by further parameters.
Furthermore, the precise values for many of the inputs are unknown and are
instead described by a distribution, meaning that any given scenario can be
simulated using a potentially intractable number of input configurations.

Covasim also suffers from long execution times and high computational costs.
Non-trivial runs of Covasim can take hours and accumulate large amounts of
data. To compound this issue further, the model is also non-deterministic:
running the same simulation parameters multiple times (with a different seed)
will yield different results, meaning that each modelling scenario must be
simulated several times to observe a distribution of outcomes.

Additionally, Covasim encounters the oracle problem: for most modelling
scenarios, the precise expected output is unknown. This makes Covasim a
traditionally “untestable” [28] system as it is difficult to determine whether
the output of a given test is correct.

Despite these challenges, Covasim features a mixture of unit, integration, and
regression tests achieving 88% code coverage.111Code coverage obtained from
commit
https://github.com/InstituteforDiseaseModeling/covasim/commit/7da3bc46e2344fa8128dfa66c260cadf4213bea27da3bc4.
However, many of these tests lack a test oracle and appear to rely on the user
to determine correctness instead. For example, the vaccine intervention has
two tests [29] that instantiate and run the model with two different vaccines
and plot the resulting model outputs on a graph for manual inspection.

While the existing vaccination tests reveal the difference in outcome _caused_
by changing from one vaccine to another, the experimental approach employed
would not scale well if the tester wanted to test more general properties that
cover larger value ranges. For example, tests covering multiple versions of
vaccine (Pfizer, Moderna, etc.) and outcomes (infections, hospitalisations,
etc.). However, this is not a criticism of Covasim, but a statement that
conventional testing techniques are impractical for testing scientific
modelling software. Hence, there is a clear need for testing techniques more
sympathetic to their challenging characteristics.

###  2.3 Metamorphic Testing

One technique which has been advocated as a basis for testing scientific
software [4] is metamorphic testing [19]. The basic idea is to define
_metamorphic relations_ that describe the expected change in output in
response to a change in input. For example, to test an implementation of
$\sin$, we may assert that $\forall x.\,\sin(x)=\sin(\pi-x)$. These relations
provide a means of generating test cases and validating the observed behaviour
[30]. By stating the expected behaviour in terms of _changes_ to inputs and
outputs, we can test the system without knowing a precise expected outcome.

Statistical metamorphic testing (SMT) [31] generalises this to non-
deterministic systems, which produce different outputs when run repeatedly
under identical input configurations. Rather than comparing outputs directly,
the SUT is run multiple times for each input configuration and statistical
tests are performed on the corresponding distributions of outputs. However,
the high computational costs involved in this process are a major limitation
to the applicability of SMT to scientific models.

###  2.4 Causal Inference

The framework we present in Section 3 uses a family of statistical techniques,
known as causal inference (CI), designed to make claims about causal
relationships between variables [32]. Our goal is to use this family of
techniques to provide an efficient method for testing scientific software.
Here we provide a brief introduction to the essential notions of CI used in
this work. For a more comprehensive overview, we refer the reader to [33, 16].

####  2.4.1 Preliminaries

Causality is often presented in terms of the “ladder of causality” [34], which
groups different tasks into three ‘rungs’: Rung one is _observation and
association_ as per traditional statistical methods; Rung two is
_intervention_ , which imagines the effects of taking particular actions:
“What if I do…?”, and rung three is _counterfactual_ , which imagines the
effects of retrospective actions: “What if I had done…?”.

Traditional statistical approaches are limited to rung one. By simply
observing the association between variables (in our case input and output
variables), without systematically controlling the selection of values or
resorting to additional domain knowledge, it is impossible to answer
fundamentally _causal_ questions [33]. This problem is commonly captured by
the adage: “correlation does not imply causation”.

CI enables us to estimate and quantify causal effects in order to make claims
about causal relationships [32]. Informally, the causal effect of a treatment
$T$ on an outcome $Y$ is the change in $Y$ that is caused by a specific change
in $T$ [34]. In this context, a _treatment_ is a variable that represents a
particular action or intervention, such as changing a line of code, and an
_outcome_ is an observable feature or event, such as the occurrence of a
fault.

One of the main challenges underlying CI is the design of experiments or
statistical procedures that mitigate sources of bias to isolate and measure
causality (rungs two and three) as opposed to association (rung one). At a
high level, this challenge can be broken down into two tasks: identification
and estimation. Identification involves identifying sources of bias that must
be adjusted for (either experimentally or statistically) in order to obtain a
causal estimate. Estimation is the process of using statistical estimators,
adjusted for the identified biasing variables, to estimate the causal effect.

####  2.4.2 Metrics

Several metrics can be used to measure causal effects. Perhaps the most
desirable is the _individual treatment effect_ (ITE), which describes the
effect of a given treatment on a particular individual. In the majority of
cases, however, individual-level inferences are unattainable due to the
_fundamental problem of causal inference_ [35]; namely that, for a given
individual, it is usually only possible to observe the outcome of a single
version of treatment (e.g. an individual either takes an aspirin for their
headache or does not).

To address this, researchers typically turn to population-level causal
metrics, such as the _Average Treatment Effect_ (ATE):

$\text{ATE}=\sum_{z\in Z}\mathds{E}{[Y\mid X=x_{t},Z=z]}-\mathds{E}{[Y\mid
X=x_{c},Z=z]}$

The ATE quantifies the average additive change in outcome we expect to observe
in response to changing some treatment variable $X$ from the _control value_
$x_{c}$ to the _treatment value_ $x_{t}$, while adjusting for all biasing
variables $Z$. However, in some instances, it is desirable to refine our
inferences to specific sub-populations defined by some notable characteristic.
To this end, the conditional ATE (CATE) can be obtained by applying the ATE to
specific sub-populations of interest [36].

An alternative causal metric is the _Risk Ratio_ (RR):

$\text{RR}=\sum_{z\in Z}\frac{\mathds{E}{[Y\mid
X=x_{t},Z=z]}}{\mathds{E}{[Y\mid X=x_{c},Z=z]}}$

The RR captures the multiplicative change in an outcome $Y$ caused by changing
the treatment variable $X$ from the control value $x_{c}$ to the treatment
value $x_{t}$ while adjusting for all biasing variables $Z$.

Other effect metrics such as the _odds ratio_ (OR) and the _effect of
treatment on the treated_ (ATT) also exist but fall outside the scope of this
paper. Furthermore, to quantify uncertainty, effect measures are typically
accompanied by 95% confidence intervals that quantify the interval within
which we are 95% confident the true estimate lies [37].

###  2.5 Causal DAGs

CI generally depends on domain expertise and causal assumptions that cannot be
tested in practice [38]. Given that different domain experts may make
different assumptions about the same problem and that these may lead to
different results, it is essential that all assumptions are made transparent.
To this end, causal DAGs provide an intuitive graphical method for
communicating the causal assumptions necessary to solve CI problems [15].
Formally, a causal graph is defined as follows [16]:

######  Definition 2.3.

A causal graph $G$ is a directed acyclic graph (DAG) $G=(V,E)$ comprising a
set of nodes representing random variables, $V$, and a series of edges, $E$,
representing causality between these variables, where:

  1. 1.

The presence/absence of an edge $V_{i}\to V_{j}$ represents the
presence/absence of a direct causal effect of $V_{i}$ on $V_{j}$.

  2. 2.

All common causes of any pair of variables on the graph are themselves present
on the graph.

In Figure 1, $\mathtt{X}$, $\mathtt{Y}$, and $\mathtt{Z}$ are nodes
representing _random variables_ , which, in this context, are variables that
can take different values for different individuals (e.g. people or software
executions). We say that $\mathtt{X}$ is a _direct cause_ of $\mathtt{Y}$
because there is an edge from $\mathtt{X}$ directly into $\mathtt{Y}$. We
refer to $\mathtt{Y}$ as a _descendant_ of $\mathtt{Z}$ and $\mathtt{X}$
because there is a sequence of edges, known as a _path_ , such that, if you
follow the direction of those edges, you can reach $\mathtt{Y}$ from
$\mathtt{Z}$. That is, $\mathtt{Z}\to\mathtt{X}\to\mathtt{Y}$. Throughout this
paper, we will use the terms ‘causal graph’ and ‘causal DAG’ interchangeably.

$\mathtt{X}$$\mathtt{Y}$$\mathtt{Z}$ Figure 1: An example causal DAG for the
causal effect of $\mathtt{X}$ on $\mathtt{Y}$ confounded by $\mathtt{Z}$

As mentioned in the previous section, in order to estimate the causal effect
of $\mathtt{X}$ on $\mathtt{Y}$, we need to identify and adjust for all
variables that bias the relationship $\mathtt{X\to Y}$. Using a causal graph,
we can achieve this automatically by applying a pair of graphical tests, the
_back-door criterion_ and _d-separation_ , which are formally defined as
follows:

######  Definition 2.4.

A path $p$ is _blocked_ or _d-separated_ by a set of variables $Z$ if and only
if at least one of the following conditions hold [39]:

  1. 1.

$p$ contains a chain $i\to k\to j$ or a fork $i\leftarrow k\to j$ where $k\in
Z$

  2. 2.

$p$ contains a collider $i\to k\leftarrow j$ where $k\notin Z$ and for all
descendants $k^{\prime}$ of $k$, $k^{\prime}\notin Z$.

######  Definition 2.5.

A set of variables $Z$ is said to satisfy the _back-door criterion_ relative
to an ordered pair of variables $(X,Y)$ if both of the following conditions
hold [39]:

  1. 1.

No variable in $Z$ is a descendant of $X$.

  2. 2.

Z blocks every path between $X$ and $Y$ that contains an arrow into $X$.

A set of variables $Z$ is said to be a _sufficient adjustment set_ relative to
a pair of variables ($X$, $Y$) if adjusting for $Z$ blocks all back-door paths
between $X$ and $Y$. Conceptually, this corresponds to a set of variables
that, once adjusted for, mitigate all known sources of bias and that is
therefore capable of isolating the _causal effect_ of interest. For example,
in Figure 1, $\mathtt{Z}$ satisfies the back-door criterion relative to
$\mathtt{\left(X,Y\right)}$ because $\mathtt{Z}$ blocks every path between
$\mathtt{X}$ and $\mathtt{Y}$ with an arrow into $\mathtt{X}$. Therefore, we
can endow the ATE of $\mathtt{X}$ on $\mathtt{Y}$ with a causal interpretation
and estimate its value directly using the following closed-form statistical
expression:

$\sum_{\mathtt{z\in
Z}}\mathds{E}{[\mathtt{Y|X=1,Z=z}]}-\mathds{E}{[\mathtt{Y|X=0,Z=z}]}$

Overall, causal DAGs provide a principled and automated approach for designing
statistical ‘recipes’ capable of measuring causal relationships and endowing
statistical measures with causal interpretations. In the following section, we
introduce a framework that facilitates the application of this approach to the
testing of scientific modelling software. Furthermore, we opt to use graphical
CI over other CI frameworks, such as potential outcomes [40] or structural
equation modelling [41], as it provides a transparent and intuitive way to
both specify and test causal relationships, without requiring users to know
their precise functional form.

##  3 Causal Testing Framework

This section introduces the Causal Testing Framework (CTF): a conceptual
framework that approaches causality-driven testing activities as CI problems.
That is, testing activities that intend to establish the (inherently causal)
relationship between inputs and outputs, such as metamorphic testing. By
framing testing activities in this way, it is possible to leverage CI
techniques to make strong claims about causal relationships between inputs and
outputs, and to do so in an efficient manner by exploiting data from previous
test executions.

In the remainder of this section, we define four key components of our causal
testing framework: specifications ($S$), programs ($P$), tests ($T$), and
oracles ($O$) [42], giving an example using Covasim (see Section 2) for each
component. We also provide informal guidance for constructing causal DAGs and
examine the relationship between the CTF and metamorphic testing.

###  3.1 Causal Specification

In the CTF, our primary aim is to test scientific models in terms of the
effects of interventions. Given the diverse range of possible scenarios that a
typical scientific model can simulate, we further focus on testing individual
modelling scenarios. We define a modelling scenario as a series of constraints
placed over a subset of the SUT’s (see Definition 2.1) input variables that
characterise the scenario of interest. Therefore, in the causal testing
framework, the set of programs, $P$, are programs that implement modelling
scenarios $\mathcal{M}$ (Definition 3.1), denoted $\mathcal{P_{M}}$.

######  Definition 3.1.

A _modelling scenario_ $\mathcal{M}$ is a pair $(X,C)$ where $X$ is a non-
strict subset of the model’s input variables and $C$ is a set of constraints
over valuations of $X$, which may be empty.

The expected behaviour of scientific modelling software in a given scenario
depends on a series of underlying modelling assumptions. It is therefore
essential that such modelling assumptions are made transparent and readily
available, particularly for the purposes of testing. Indeed, past
investigations into modelling failures have highlighted the importance of
transparency and accountability [1]. In the same vein, causal testing requires
an explicit record of causal assumptions to enable the transparent and
reproducible application of graphical CI techniques. To this end, we use a
causal DAG that captures causality amongst a subset of the SUT’s input and
outputs. Therefore, we define a causal specification (Definition 3.2) as a
pair comprising a modelling scenario ($\mathcal{M}$) and a causal DAG
($\mathcal{G}$).

######  Definition 3.2.

A _causal specification_ is a pair $\mathcal{S}=(\mathcal{M},\mathcal{G})$
comprising a modelling scenario $\mathcal{M}$ and a causal DAG $\mathcal{G}$
capturing the causal relationships amongst the inputs and outputs of the SUT
that are central to the modelling scenario.

######  Example 3.1.

Consider a scenario in Covasim where we want to test the effect of
prioritising the elderly for vaccination $V$ on the total vaccine doses
administered $N_{D}$, total vaccinated agents $N_{V}$, maximum number of doses
per agent $M_{D}$, and cumulative infections $I$. Further, let us restrict our
simulation length to 50 days, the initial number of infected agents to 1000,
and the population size to 50,000. Our modelling scenario is then
characterised by the constraints
$\\{\mathtt{days}=50,\mathtt{pop\\_size}=50000,\mathtt{pop\\_infected}=1000\\}$,
and the causal DAG is the set of edges $V\to\\{N_{V},N_{D},I\\}$. Note the
absence of edge $V\to M_{D}$. Here we are asserting that $V$ may cause a
change in $N_{V}$, $N_{D}$, and $I$, but should cause no change to $M_{D}$.
This is because at most two doses of the vaccine are administered to each at
agent so changing the target population should not affect this.

###  3.2 Constructing Causal DAGs

In the testing context, causal DAGs offer a flexible, lightweight means by
which to capture potential causal relationships between inputs and outputs.
Here we present a set of guidelines for constructing the graph (informed by
our experience with the case studies).

We start by constructing a complete directed graph over the set of inputs and
output: $I\cup O$. Then, to simplify this structure, we apply the following
assumption:

######  Assumption 1.

Outputs must temporally succeed inputs and, therefore, cannot cause inputs.

Assumption 1 follows from temporal precedence (that a cause must precede its
effect) [43] and enables us to delete all edges from outputs to inputs.

Then, in many cases, we can also apply the following assumption to remove all
edges from inputs to inputs:

######  Assumption 2.

Inputs cannot cause changes to the values of other inputs and, therefore,
cannot share causal relationships.

As stated in Definition 2.1, in this paper, we assume that all inputs are
assigned their values prior to execution. Under this characterisation, changes
to the value of one input cannot _physically_ affect another input’s value
and, therefore, inputs cannot share causal relationships. Of course, there are
caveats to this; if a system has input validation, for example, the assignment
of one input’s value may _physically_ restrict which values can be selected
for a second input. Note that, in such cases, our framework is still
applicable, but the user would have to consider more edges manually to
construct their DAG.

This leaves us with the following forms of potential causal relationships to
consider: $I\to O$ and $O\to O$ (and $I\to I$ if Assumption 2 cannot be
applied). Output to output causality may occur in software where an earlier
output is used in the computation of a later output. For example, in a weather
forecasting model, a prediction of the weather in three days time is affected
by the weather predicted for one and two days time.

This is the point at which the tester’s domain knowledge is fed into the
model, by pruning edges where they are certain that there is no causal
relationship (see Definition 2.2 for an informal definition of causality). We
recommend following this approach of pruning edges from a complete directed
graph over adding edges to an initially empty graph, as the absence of an edge
carries a stronger assumption than the presence of one [18]. This follows from
the fact that the presence of an edge states that there exists _some_ causal
relationship, whereas the absence of an edge states that there is _precisely_
no causal relationship.

###  3.3 Causal Testing

Causal testing draws its main inspiration from CI, which focuses on the
effects of _interventions_ on _outcomes_. In this context, an intervention
manipulates an input configuration in a way that is expected to _cause_ a
specific outcome to change. Here, we refer to the pre-intervention input
configuration as a _control_ and the post-intervention input configuration as
a _treatment_. A causal test case then specifies the expected causal effect of
this intervention. When phrased this way, causal tests bear a remarkable
similarity to metamorphic tests, highlighting the fact that, at its core,
metamorphic testing can be viewed as an inherently a causal activity. We
explain this relationship further in Section 3.4.

######  Definition 3.3.

An _intervention_ $\Delta:\mathcal{X}\to\mathcal{X^{\prime}}$ is a function
which manipulates the values of a subset of input valuations.

######  Definition 3.4.

A _causal test case_ is a 4-tuple
$(\mathcal{M},\mathcal{X},\Delta,\mathcal{Y})$ that captures the expected
causal effect, $\mathcal{Y}$, of an intervention, $\Delta$, made to an input
valuation, $\mathcal{X}$, on some model outcome in the context of modelling
scenario $\mathcal{M}$.

######  Example 3.2.

Continuing with our vaccination example, our intervention
$\Delta(\mathcal{X})=\mathcal{X}[\mathtt{vaccine}:=\mathtt{Pfizer^{\prime}}]$
is a function that replaces the $\mathtt{Pfizer}$ vaccine in control input
configuration $\mathcal{X}=\mathtt{\\{vaccine=Pfizer\\}}$ with an age-
restricted version ($\mathtt{Pfizer^{\prime}}$), yielding the treatment input
configuration $\mathcal{X^{\prime}}=\mathtt{\\{vaccine=Pfizer^{\prime}\\}}$.
We then specify the expected outcomes for the causal test cases. For example,
the expected ATE of the intervention on $M_{D}$ is 0, since it should remain
unchanged: $\mathcal{Y}=(\mathtt{ATE_{M_{D}}=0})$.

Finally, we must consider the test oracle: the procedure used to determine
whether the outcome of a causal test case is correct. In the context of causal
testing, the oracle must ascertain the correctness of causal estimates
relative to a modelling scenario ($\mathcal{M}$). Therefore, we refer to our
oracle as a causal test oracle (Definition 3.3).

######  Definition 3.5.

A _causal test oracle_ $\mathcal{O}$ is an assertion placed over the causal
effect $\mathcal{Y}$ of an intervention $\Delta$ made to an input valuation
$\mathcal{X}$, to ascertain the validity of the outcome of a causal test case.

######  Example 3.3.

In our running Covasim vaccine example, the oracle to check whether the
maximum number of doses has remained unchanged corresponds to an assertion:
$\mathtt{ATE_{M_{D}}=0}$.

Any discrepancy between the test result and the expected outcome revealed by
the test oracle implies one of two problems: (i) the implementation contains a
bug or an error, or (ii) the causal DAG is incorrect. It follows that causal
testing lends itself to an iterative testing process [44], whereby the user
inspects the source code to explain any identified discrepancies and, if no
bugs are found, reviews the causal model to check if the underlying science is
correct.

Collectively, the components of the CTF enable the application of graphical CI
techniques to testing activities that concern the causal effect of some
intervention. In theory, the CTF should therefore provide the following
advantages over existing solutions:

  1. 1.

The ability to derive test outcomes _experimentally_ (by strategic model
executions) and _observationally_ (by applying CI techniques to past execution
data). In a testing context, this would enable the tester to not only run test
cases on the system directly, but also to infer test outcomes
_retrospectively_ from existing test data.

  2. 2.

The ability to identify and adjust for confounding bias in observational data
using a causal DAG. From a testing perspective, this effectively relaxes the
experimental conditions ordinarily required to reach causal conclusions.
Namely, the need for carefully controlled, unbiased test data.

  3. 3.

The ability to derive _counterfactual_ test outcomes using appropriate
statistical models. This would enable testers to infer how the model would
likely behave, had it been run under a different parameterisation. Therefore,
where practical constraints preclude further executions of the SUT,
counterfactual inference can offer a cost-effective alternative.

In Section 5, we apply the CTF to a series of real-world scientific models to
explore the extent to which these advantages are realised in a testing
context. Furthermore, we demonstrate how, collectively, the aforementioned
advantages can offer an efficient approach for exploratory testing [45, 46] of
scientific modelling software.

###  3.4 Relationship to Metamorphic Testing

At a high level, the CTF and metamorphic testing share the same objective: to
evaluate the _effect_ caused by making a change to some input.

Metamorphic testing provides a means of generating “follow-up test cases”
using metamorphic relations which should hold over a number of different
parameter values [6, 30]. In contrast to typical program invariants, which
must hold for every execution of a given program, metamorphic relations hold
between different executions. In other words, they investigate the effect of a
change (or _intervention_ in causal language) on an input. This is a key
similarity between causal testing and metamorphic testing.

In this sense, metamorphic tests can be thought of as experiments designed to
answer causal questions about the SUT. For example, a metamorphic test for our
property of the $\sin$ function in Section 2 that $\forall
x.\,\sin(x)=\sin(\pi-x)$ can be thought of as an experiment that confirms
whether changing the input from $X=x$ to $X=\pi-x$ causes no change to the
output. That is, there should be _no causal effect_. This synergism suggests
that metamorphic testing can be re-framed and solved as a problem of CI and,
therefore, benefit from its advantages. To this end, in Section 5, we
demonstrate how the CTF can conduct metamorphic testing using CI techniques.

One advantage of causal testing over traditional metamorphic testing is that
causal testing does not necessarily require dedicated test runs of the system
to be performed if sufficient log data already exists. Even (and especially)
if this data is biased, CI can account for this, meaning that testing can be
performed on systems which cannot be tested for reasons of practicality.
Furthermore, systems can be tested retroactively, enabling concerns about a
model’s correctness to be investigated even after the model has been run. This
is potentially advantageous in the context of scientific models, where their
integrity and correctness can be called into question years after policies
based on their output have already been made. In such situations, the DAG
makes clear the assumptions made about the functionality of the model so adds
weight to any conclusions made.

##  4 CTF Reference Implementation

This section provides an overview of our open-source Python reference
implementation of the Causal Testing Framework
(CTF),222https://github.com/CITCOM-project/CausalTestingFramework comprising
over 4000 lines of Python code, and outlines four stages of the CTF workflow:
Specification, Test Cases, Data Collection, and Testing.

###  4.1 Causal Specification

To begin causal testing, we form a causal specification (Definition 3.2),
comprising two components: a modelling scenario and a causal DAG. We form the
modelling scenario by specifying a set of constraints over the inputs that
characterise the scenario-under-test, such as $x_{1}<x_{2}$. Next, we specify
our causal DAG using the DOT language [47], in which graphs are expressed as a
series of edges, such as $x_{1}\to x_{2}$, following the guidelines outlined
in Section 3.2.

###  4.2 Causal Test Case

Now that we have a causal specification, we define a causal test case that
describes the intervention whose effect we wish to test. In our reference
implementation, a causal test case is an object that requires us to specify a
control input configuration, a treatment input configuration, and an expected
effect. In the following steps, this information will enable us to collect
appropriate test data (Data Collection), design experiments isolating the
causal effect of interest within this data, and define test oracles that
ascertain whether the expected causal effect is observed (Causal Testing).

###  4.3 Data Collection

After creating a causal specification and causal test case, the next step is
to collect data corresponding to the modelling scenario. We can achieve this
either experimentally (in situations where we are able to directly execute the
SUT) or observationally (in situations were we are not able to execute the
SUT, but are instead able to draw upon prior execution data).

####  4.3.1 Experimental Data Collection

Experimental data collection executes the model _directly_ under both the
control and treatment input configuration to isolate the causal effect of the
intervention. To this end, our reference implementation provides an abstract
experimental data collector class, requiring us to implement one method that
executes our model with a given input configuration. This method enables the
CTF to run the model under the experimental conditions necessary to isolate
causality directly.

####  4.3.2 Observational Data Collection

Since it is often infeasible to run models a statistically significant number
of times, we also provide the option to use observational, existing test data.
This data may not meet the experimental conditions necessary to isolate the
causal effect and thus may contain biases that lead purely statistical
techniques astray. However, by employing graphical CI techniques, the CTF can
identify and mitigate bias in the data, providing an efficient method for
testing scientific models _a posteriori_.

There are two caveats to this. First, the causal DAG must be correctly
specified. Second, the observational data must be consistent with the
constraints of the causal specification. To this end, our reference
implementation includes an observational data collector class that takes a CSV
file of existing test data as input and uses the Z3 theorem prover [48] to
identify and remove any model runs that violate constraints. Next we describe
how the CTF infers test outcomes from this data.

###  4.4 Causal Testing

Given a causal test case, testing is carried out in two stages: causal
inference (CI) and applying the test oracle.

####  4.4.1 Causal Inference

To infer the causal effect of interest, our reference implementation applies
the two steps of CI outlined in Section 2: identification and estimation. For
identification, the CTF algorithmically identifies an adjustment set (see
Section 2.4) for the causal effect of interest. Then, for estimation, we
design an appropriate estimator that adjusts for the identified adjustment
set, and apply the estimator to our data to estimate the desired causal metric
(e.g. ATE or RR, see Section 2). To this end, our reference implementation
provides customisable regression and causal forest [49] estimators. In
addition, the CTF includes an abstract estimator class that enables users to
define their own estimators. This step outputs a causal test result containing
the inferred causal estimate for the desired causal metric (e.g. ATE or RR,
see Section 2.4) and 95% confidence intervals.

####  4.4.2 Test Oracle

After applying CI, all that remains is the test oracle procedure. That is, to
check whether the causal test results match our expectations. For this
purpose, our reference implementation provides several test oracles that check
for positive, negative, zero, and exact effects. Alternatively, to handle more
complex outputs, a user can specify a custom oracle that ascertains whether a
causal test result should pass or fail.

Now that we have discussed the workflow of our CTF reference implementation,
in the following section, we demonstrate its application to three vastly
different real-world scientific models.

##  5 Case Studies

This section demonstrates the application of the Causal Testing Framework
(CTF) to three real-world scientific models from three different domains. Our
goal here is not to comprehensively test these models, but rather to
investigate the extent to which the CTF realises the theoretical advantages
outlined in Section 3. The code to execute the causal tests in the following
case studies can be found in the examples directory of the reference
implementation.333https://github.com/CITCOM-project/CausalTestingFramework

Each case study follows the same high-level structure. We start by providing a
brief overview of the testing activity (the broad testing objective) and the
process of acquiring data for the analysis. We then describe the application
of the CTF to analyse the generated data, and conclude by analysing the
outcomes and answering the relevant research questions.

###  5.1 Research Questions

Across our case studies, we corroborate evidence to collectively answer three
research questions. The contribution of each case study to the research
questions will be highlighted throughout the case studies and the collective
findings will be discussed in Section 6.

##### RQ1: Can we apply Causal Testing both experimentally and
observationally?

As stated in Section 3, one of the advantages of CI is the ability to draw
causal conclusions from both experimental and observational data. In a testing
context, this corresponds to the CTF having the ability to not only execute
test cases via controlled experiments (i.e. traditional test executions), but
also to infer test outcomes from existing test data.

We consider RQ1 to be answered in the affirmative if the following success
criterion is satisfied: The CTF must be able to measure the causal effect of
an intervention by (i) _directly executing_ the model under the conditions
necessary to isolate the causal effect of interest, and (ii) applying CI
techniques to _infer test outcomes from existing test data_.

##### RQ2: Can Causal Testing identify and mitigate confounding bias in
observational test data?

As stated in Section 2, the distinction between statistical and causal
estimates lies in the ability to identify and control for sources of
confounding bias. Therefore, in order to endow inferred test outcomes with a
causal interpretation, the CTF must also have the ability to identify and
adjust for confounding bias.

We consider RQ2 to be answered in the affirmative if the following success
criterion is satisfied: The CTF must be able to _identify confounding
variables_ and _adjust for their effect_ using statistical models.

##### RQ3: Can Causal Testing handle counterfactual test outcomes?

From a testing perspective, the ability to reason about counterfactual
outcomes would enable us to predict (causally) how a model will behave in new,
unobserved scenarios without having to run the SUT itself, based on our
knowledge of how it behaves in existing scenarios and how it is expected to
respond to changes (interventions).

We consider RQ3 to be answered affirmatively if the following success
criterion is satisfied: The CTF must be able to infer the outcome of causal
tests that _have not been observed in previous test data_ and _without futher
executions_ of the SUT.

###  5.2 Subject Systems

Our case studies apply the CTF to three different subject systems: a Poisson
tessellation model, a cardiac action potential model, and an epidemiological
agent-based model. Here, we provide a high-level description of these systems
and their practical importance.

##### Poisson Line Tessellation Model

Our first subject system is the Poisson Line Tessellation model (PLT). This
model uses a Poisson process to generate a series of lines that are positioned
and oriented at random within a given sampling window to form a tessellation.
While the behaviour of this model is predominantly random by design, it can be
configured using three numerical input parameters to produce tessellations
with predictable properties. In order to test these properties, we extract
four numerical outputs from the resulting tessellation.

We selected this model as a case study because it has been the subject of
prior research on statistical metamorphic testing [31] and comprises a number
of well-characterised input-output relationships. In addition, Poisson process
models are commonly used to model random processes for a range of
applications, including simulating road networks [50, 51] and modelling photon
arrival in 3D imaging [52]. It is the stochastic yet predictable behaviour of
Poisson process models that make them an interesting but difficult subject to
test.

##### Cardiac Action Potential Model

Our second subject system is the Luo-Rudy 1991 ventricular cardiac action
potential model [53] (LR91). This is a mathematical model comprising a system
of differential equations that describe the membrane action potential of a
mammalian ventricular cell. The behaviour of this model is controlled by 24
constants, 8 rate variables, 8 state variables, and 25 algebraic variables, to
produce a time-series describing how the action potential of the modelled
cardiac cell varies over time.

We selected LR91 as a case study as it follows a different modelling paradigm
to our other subject systems and has supported extensive and important
research into cardiovascular physiology. Furthermore, amongst its vast and
largely uncertain input space, LR91 has several well-characterised input-
output relationships suitable for causal analysis.

##### Epidemiological Agent-Based Model

Our third subject system is Covasim; the epidemiological agent-based model
introduced as a motivating example in Section 2. As a brief reminder, Covasim
is a complex, real-world scientific model that is primarily used to simulate
detailed COVID-19 scenarios in order to evaluate the impact of various
interventions, such as vaccination and contact tracing [3], in specific
demographics. These scenarios are configured via 64 input parameters and
described by 56 time-series outputs.

We selected Covasim as a case study as it has a number of challenging
characteristics that present a significant testing challenge (see Section 2).
Moreover, it has been used to inform a number of important policy decisions
across a range of countries, including the UK, US, and Australia [54, 24, 26,
55], demonstrating its practical importance.

We cover two testing scenarios using Covasim. In Section 5.5, we return to our
example from Section 2 and experimentally test the effect of vaccinating the
elderly on several vaccine-related outcomes via strategic executions of the
SUT, revealing an interesting bug in the process. Then, in Section 5.6 we test
the effect of increasing the $\beta$ parameter (infectiousness) on cumulative
infections using existing confounded test data. The second scenario is notable
because it deals with the realistic setting where the metamorphic property is
subject to a significant degree of uncertainty - a common issue when testing
scientific models [4].

###  5.3 Poisson Line Tessellation Model

In this case study, we use the CTF to conduct SMT on a Poisson Line
Tessellation model. We have two objectives here. First, we aim to show that
the CTF can test the same metamorphic relations as a conventional SMT
approach. Second, we aim to demonstrate how the CTF enables this to be
achieved without the need for the large numbers of carefully controlled test
executions that are ordinarily required for SMT. We now describe the behaviour
of the PLT model, referring to the example tessellation in Figure 2.

$-$1.0$-$0.50.00.51.0$-$1.0$-$0.50.00.51.01234

Figure 2: A tessellation generated by the PLT model with a width ($W$), height
($H$), and intensity ($I$) of 1. There are two lines which intersect the
sampling window ($L_{t}$, highlighted in grey). The intersection of these
lines forms four polygons in total ($P_{t}$).

The PLT model has three positive floating point input parameters: the width
$W$ and height $H$ of a sampling window (shaded in grey in Figure 2), and the
intensity $I$ of the Poisson process. Informally, the intensity parameter
controls the average rate at which lines are placed. Given these inputs, the
model generates a set of straight lines that intersect the origin-centred
sampling window by drawing from a Poisson process on
$[0,\infty)\times[0,2\pi)$444The interval $[0,\infty)$ corresponds to the
random distance of the lines to the origin, and the interval $[0,2\pi)$
corresponds to the random angle of the point on the line that is closest to
the origin. In the case of the orientation distribution, the upper interval
bound is $\pi$ since rotating a line by an angle of $\pi$ (i.e. 180 degrees)
leads to the same orientation., where the orientation is uniformly distributed
on $[0,\pi]$. The model then outputs the total number of lines intersecting
the sampling window, $L_{t}$, and the number of polygons formed by the
intersecting lines, $P_{t}$.

In Figure 2, for example, the inputs $W=H=I=1$ produce a tessellation in which
there are two lines intersecting the sampling window ($L_{t}=2$) that form
four polygons ($P_{t}=4$). Then, by dividing $L_{t}$ and $P_{t}$ by the
sampling window area (i.e. $W\times H$), we obtain two further outputs
corresponding to the number of lines and polygons per unit area ($L_{u}$ and
$P_{u}$, respectively). Since $W=H=1$ in Figure 2, it follows that
$L_{u}=L_{t}=2$ and $P_{u}=P_{t}=4$.

####  5.3.1 Testing Activity

In this case study, we replicate the SMT approach followed by Guderlei et al.
[31] to explore whether the CTF can achieve similar results to traditional SMT
approaches, and investigate whether it can do so without the need for a large
number of test cases (as is usually the case with SMT). To this end, we
implement the approach in the CTF using both the _experimental_ and
_observational_ mode of operation.

First, to conduct SMT experimentally, we directly run the SUT under the
conditions necessary to isolate and measure the causal effect of interest,
such as the effect of $I$ on $P_{u}$, numerous times. Then, to conduct SMT
observationally, we repeat the approach, but this time using graphical CI to
infer test outcomes from a smaller amount of existing, randomly generated
data.

In both the experimental and observational approach, we expect to observe that
two metamorphic properties from the original study will hold [31] :

  1. 1.

Doubling $I$ should cause $P_{u}$ to increase by a factor of 4.

  2. 2.

$P_{u}$ should be independent of $W$ and $H$.

####  5.3.2 Data Generation

We generated two sets of execution data. First, to replicate the SMT approach
followed in the original study [31], we sampled 50 input configurations, with
width and height incremented together over the interval $[1,10]$ (i.e.
$W=H=1,W=H=2,\ldots,W=H=10$), such that the sampling window is always square,
and the control and treatment values for intensity are powers of 2 up to 16.
We then executed each configuration $100$ times to account for non-
determinism, resulting in $5000$ model runs.

Second, to explore how the CTF enables us to reuse observational past
execution data to infer the outcome of metamorphic test cases, we generated
1000 random input configurations and executed each one a single time. To
generate these input configurations, we used Latin hypercube sampling [56, 57]
of the distributions $W,H\sim\mathbf{U}(0,10)$ and $I\sim\mathbf{U}(0,16)$ to
provide even coverage of the input space.

####  5.3.3 Causal Testing

To begin causal testing, we first specify our modelling scenario using
constraints $\\{0<W\leq 10,0<I\leq 16,W=H\\}$ and the causal DAG in Figure 3,
the relationships in which can be explained by theoretical results for $L_{t}$
and $P_{t}$ [58]. Specifically, for a given PLT run with $W=w$, $H=h$, and
$I=i$, we expect $L_{t}=2i(w+h)$ and $P_{t}=\pi i^{2}wh$. We include an edge
$L_{t}\to P_{t}$ because $P_{t}$ is determined by the intersection of lines.
We also include edges from $W$ and $H$ to $L_{u}$ and $P_{u}$ because these
quantities are calculated by dividing their respective total values by the
window area.

$W$$I$$H$$L_{t}$$P_{t}$$L_{u}$$P_{u}$ Figure 3: A causal DAG for the PLT
model.

Having created our causal specification, we now conduct a series of causal
tests to investigate the two properties mentioned above: (1) whether doubling
$I$ causes $P_{u}$ to increase by a factor of 4, and (2) whether the sample
window size has a causal effect on $P_{u}$.

##### Effect of $I$ on $P_{u}$

First, we test whether doubling $I$ causes $P_{u}$ to increase by a factor of
4 for $I\in\\{1,\ldots,16\\}$ and $W,H\in\\{1,\ldots,10\\}$. Since we are
interested in the multiplicative effect of $I$ on $P_{u}$, we use the _risk
ratio_ (RR, see Section 2), which quantifies the factor by which the
intervention (doubling $I$) causes the outcome change:

$\text{RR}=\frac{\mathds{E}{[P_{u}\mid I=i_{t}]}}{\mathds{E}{[P_{u}\mid
I=i_{c}]}}$

Using the SMT-style data, we calculate the RR directly by contrasting the
empirical $P_{u}$ means for the control and treatment values of $I$. On the
other hand, for the random dataset, we must use the CTF to perform
identification on the causal DAG in Figure 3. This reveals that there is no
confounding to adjust for in this scenario and that we can therefore use a
regression model of the form $P_{u}\sim x_{1}I+x_{2}I^{2}$ to derive a causal
estimate. We include an $I^{2}$ term here as we expect $P_{u}$ to vary
quadratically with $I$ (this decision is informed by domain expertise but can
be validated by varying $I$ and observing changes to $I$). Additionally,
notice that our regression model does not contain $W$ or $H$ as both the
theoretical results and the original study [31] indicate that $P_{u}$ should
be independent of window size.

We can then apply this regression model to the randomly sampled data to
predict the RR. Notice that, since the second dataset is randomly generated,
it is unlikely to contain the exact values for the control and treatment
necessary to directly compute the RR. Therefore, we must rely on our causally-
valid statistical estimator to infer _counterfactuals_ in this instance.

##### Effect of $W$ on $P_{u}$

Second, we test whether the sample window size has a causal effect on $P_{u}$.
We achieve this by observing whether changes to the width of the window $W$
_cause_ a change to $P_{u}$. Since we are only interested in whether there is
_some_ effect, we use the _average treatment effect_ (ATE, see Section 2),
which quantifies the additive change in outcome caused by the intervention
(increasing $W$):

$\text{ATE}=\mathds{E}{[P_{u}\mid W=w_{t}]}-\mathds{E}{[P_{u}\mid W=w_{c}]}$

Ordinarily, to investigate whether $W$ affects $P_{u}$ using SMT, we would
need to run another set of experiments, this time fixing the value of $I$ and
varying $W$. In the CTF, however, we can use our causally valid estimators to
infer the effect of $W$ on $P_{u}$ from the _same_ 1000 model runs. We achieve
this by modifying our regression model to include terms for $W$ and $W^{-1}$,
reflecting our hypotheses that $W$ does affect $P_{u}$ and that they share an
inverse relationship (this choice is informed by domain expertise but can be
validated by varying $W$ and observing $P_{u}$). We then apply this model to
the original data to estimate the above ATE. The effect of $H$ could be
investigated similarly, but we omit this due to space constraints.

####  5.3.4 Results

Table 1: RR of doubling $I$ under different values of $W$ and $H$. The bottom row gives the value estimated using regression. Bold values round to 3, violating the expected behaviour. $W$ | $H$ | $\frac{\mathds{E}{[P_{u}\mid I=2]}}{\mathds{E}{[P_{u}\mid I=1]}}$ | $\frac{\mathds{E}{[P_{u}\mid I=4]}}{\mathds{E}{[P_{u}\mid I=2]}}$ | $\frac{\mathds{E}{[P_{u}\mid I=8]}}{\mathds{E}{[P_{u}\mid I=4]}}$ | $\frac{\mathds{E}{[P_{u}\mid I=16]}}{\mathds{E}{[P_{u}\mid I=8]}}$
---|---|---|---|---|---
1 | 1 | 2.5888 | 3.4461 | 3.6178 | 3.6187
2 | 2 | 3.0359 | 3.5410 | 3.6003 | 3.7264
3 | 3 | 3.5025 | 3.5945 | 4.0191 | 3.6545
4 | 4 | 3.1138 | 3.5285 | 4.1562 | 3.7290
5 | 5 | 3.6686 | 3.7686 | 3.9408 | 3.8751
6 | 6 | 3.6933 | 3.6988 | 3.9219 | 3.9707
7 | 7 | 3.7127 | 3.6271 | 3.9862 | 3.9370
8 | 8 | 3.4957 | 3.8300 | 3.8861 | 4.0110
9 | 9 | 3.5633 | 4.0009 | 3.9342 | 3.9338
10 | 10 | 3.8275 | 3.7525 | 4.0128 | 4.0181
Estimated | 2.8280 | 3.1711 | 3.4772 | 3.6993

Table 1 shows the results for our investigation into the effect of $I$ on
$P_{u}$. The first 10 rows show the RRs associated with the experimental, SMT-
style data for various values of $W$ and $H$, and the final row shows the RRs
estimated using the observational, random dataset. These results show that
both approaches identify an inconsistency between the metamorphic relations
and implementation from the original study [31]: for lower values of $W$, $H$,
and $I$, the RR is closer to three than four, meaning our metamorphic relation
is not satisfied. This is a particularly interesting result since $P_{u}$
should be independent of $W$ and $H$. Furthermore, our random dataset only
contains 1000 data points whereas our SMT-style dataset contains 5000,
suggesting it is more economical to estimate counterfactuals than to design
and execute controlled experiments.

Table 2: ATE of increasing W from W${}_{\text{c}}$ to W${}_{\text{t}}$ on $P_{u}$ with $I=1$ in the PLT model. W${}_{\text{c}}$ | W${}_{\text{t}}$ | ATE | 95% CIs
---|---|---|---
1 | 2 | -7.3786 | [-13.9182, -0.8390]
2 | 3 | -2.7097 | [ -9.8029, 4.3836]
3 | 4 | -1.5424 | [-11.1209, 8.0361]
4 | 5 | -1.0755 | [-13.7084, 11.5574]
5 | 6 | -0.8421 | [-16.7413, 15.0572]
6 | 7 | -0.7087 | [-19.9729, 18.5556]
7 | 8 | -0.6253 | [-23.3084, 22.0578]
8 | 9 | -0.5697 | [-26.7043, 25.5649]
9 | 10 | -0.5308 | [-30.1383, 29.0767]

Table 2 shows the results for our investigation into the effect of $W$ on
$P_{u}$. Here, each row shows how $P_{u}$ changes when $W$ is increased from
$W_{c}$ to $W_{t}$ with the intensity fixed to $I=1$. According to the
original study [31], changes to $W$ should _not cause_ a change to $P_{u}$.
Our results show that this property holds for all but the first row because
these rows have confidence intervals that contain zero, meaning there is no
statistically significant causal effect. However, the confidence intervals for
the first row of Table 2 show that, when $W$ is increased from $W=1$ to $W=2$,
there is a statistically significant causal effect on $P_{u}$ of $-7.3786$
(the intervals do not contain zero).

This conflicting result indicates a problem with either the program or the
metamorphic property. In this case, we believe that the problem stems from
basic geometry: lines are less likely to intersect a smaller sample window. As
the sample window becomes larger, there is more area to average over so
$P_{u}$ becomes more reliable. Therefore, the metamorphic relations should
perhaps specify a minimum window size to which they apply.

Overall, this case study shows that not only can we conduct SMT using the CTF,
but we can do so _using previous execution data_ and _less data_ than required
by traditional SMT methods. Furthermore, we demonstrated how this approach
allowed us to refine our metamorphic relations and find faults _without
running the SUT further times_. In doing so, this case study has provided
insights into RQ1 and RQ3.

##### RQ1

The results in Table 1 demonstrate that the CTF is able to measure causal
effects using both experimental and observational data, and that the
associated results are similar. Furthermore, the results in Table 2 show how,
using the observational approach, we could identify an inconsistency in the
metamorphic relation proposed in the original study [31].

##### RQ3

This case study involved executing test cases using observational data which
did not contain the exact executions of interest. Instead, by leveraging the
causal assumptions in Figure 3 and appropriate statistical models, we borrowed
information from similar executions to infer the outcome of the counterfactual
test cases.

###  5.4 Cardiac Action Potential Model

Our goal in this case study is to demonstrate how the CTF can accomplish
sensitivity analysis of a scientific model in a straightforward and efficient
way. This technique is commonly used to validate and verify scientific models,
with a specific focus on identifying which inputs have the greatest impact on
model outputs [59, 60]. Here, we take a CI-led approach and measure the ATE of
several input parameters on one output, $APD_{90}$ (the time taken for the
action potential to depolarise by 90%), quantifying the extent to which this
output is affected by changes to the inputs. We compare our results to an
existing Gaussian Process approach, showing that we obtain comparable findings
with the same quantity of data while using a simpler statistical model.

####  5.4.1 Testing Activity

In this case study, we replicate part of an existing study [61] that conducts
uncertainty and sensitivity analysis on LR91 using a Gaussian Process Emulator
(GPE) [62] trained on runs of the model to quantify the uncertainty in outputs
with respect to different inputs. In short, the approach in [61] trained a GPE
on 200 runs of LR91, with input configurations sampled via Latin Hyper Cube
Sampling [63] from a series of normalised uniform design distributions to
ensure even coverage of the input space. The GPE was then used to calculate
the expectation of a given output, conditional on an input of interest, to
quantify the effect of varying each of the six inputs on the eight output
parameters, over the range of the design distribution.

From a CI perspective, we can obtain similar information by computing the ATE
of each input on each output over the range of the design distribution.
Specifically, we can set our control value to the mean value of the design
distribution and uniformly increment our treatment value from the minimum to
the maximum value of the design distribution. This yields a series of ATEs
that quantify the expected change in output caused by changing the input
parameters by specific amounts above and below their mean, revealing the
magnitude of each input’s effect on the outputs.

Due to space limitations, we limit our analysis to the effect of the six
inputs on one output, $APD_{90}$. We have selected this output because the
original paper uses it to exemplify the approach. Based on the results
reported in [61], we expect an increase in the parameters $G_{K}$, $G_{b}$,
and $G_{K1}$ to cause a decrease in $APD_{90}$, an increase in $G_{si}$ to
cause an increase in $APD_{90}$, and expect the remaining parameters ($G_{Na}$
and $G_{Kp}$) to have no significant effect on $APD_{90}$. We also anticipate
the following monotonic relationship over the (absolute) magnitude of the
inputs’ effects:

$\lvert APD_{90}^{G_{si}}\rvert>\lvert APD_{90}^{G_{K}}\rvert>\lvert
APD_{90}^{G_{b}}\rvert>\lvert APD_{90}^{G_{K1}}\rvert$

####  5.4.2 Data Generation

To gather data from LR91, we followed the same approach as [61], where the 200
input configurations were sampled from the design distributions using Latin
Hyper Cube sampling and then normalised. We then executed each of these input
configurations on an auto-generated Python implementation of LR91 from the
cellML modelling library [64]. We extended this implementation to enable us to
sample the input values via Latin Hyper Cube sampling and automatically
extract the outputs.555Our LR91 model is available at:
https://github.com/AndrewC19/LR91

####  5.4.3 Causal Testing

To approach sensitivity analysis as a CI problem, we first specify our
modelling scenario and causal DAG. For this set of tests, the modelling
scenario constrains each input to the range of its uniform design distribution
(as specified in the original paper [61]):

$\displaystyle\\{$ $\displaystyle 17.250\leq G_{Na}\leq 28.750,0.0675\leq
G_{si}\leq 0.1125,0.2115\leq G_{K}\leq 0.3525,$ $\displaystyle 0.4535\leq
G_{K1}\leq 0.7559,0.0137\leq G_{Kp}\leq 0.0229,0.0294\leq G_{b}\leq 0.0490\\}$

As in the original study, these input values were then normalised to the range
[0, 1].

We then specify the expected cause-effect relationships (and absence thereof)
as the causal DAG shown in Figure 4. For each relationship, we then create a
suite of causal test cases covering a series of interventions that
incrementally increase/decrease the value of the inputs over the range of the
design distribution. For each input, this is achieved by setting the control
value to 0.5 (the mean) and uniformly sampling 10 treatment values over the
range [0, 1]. This produces a total of 10 test cases per input that vary its
value from 0.5 to each of the treatment values: $[0,\,0.1,\,0.2,\,...\,1.0]$.
Using the CTF, we then perform identification and estimation. Here, the cause-
effect relationships are straightforward and there is no confounding to adjust
for, enabling us to fit a regression model $APD_{90}\sim x_{0}+x_{1}G_{z}$ for
each input $z\in\\{si,K,Na,Kp,K1,b\\}$. Using these models, we then predict
the ATE and 95% confidence intervals for each test.

$G_{si}$$G_{K}$$G_{Na}$$G_{Kp}$$G_{K1}$$G_{b}$$APD_{90}$ Figure 4: LR91
modelling scenario’s Causal DAG, where the sensitivity of $APD_{90}$ to each
conductance input is computed as the causal effect (ATE).

####  5.4.4 Results

Figure 5: Sensitivity of $APD_{90}$ in response to changes to the mean value
of input parameters in LR91.

The results, as summarised in Figure 5, show that all expected relationships
hold with statistical significance (95% confidence intervals do not contain 0)
and are visually similar to Figure 5 in the original study [61]. Specifically,
$G_{si}$ has a positive effect, $G_{K}$, $G_{K1}$, $G_{b}$ have negative
effects, and $G_{Na}$ and $G_{Kp}$ have no significant effect. The gradients
corresponding to these effects reveal that the effect sizes follow the
expected monotonic relationship: $\lvert APD_{90}^{G_{si}}\rvert>\lvert
APD_{90}^{G_{K}}\rvert>\lvert APD_{90}^{G_{b}}\rvert>\lvert
APD_{90}^{G_{K1}}\rvert$.

Overall, in this case study, we have shown that the CTF reaches the same
conclusions as the original study. This is however achieved by using a simpler
statistical model guided by causality instead of associations within the data.
Additionally, this case study suggests that, like SMT, sensitivity analysis
can be framed and solved as a CI problem. These results have provided insights
into RQ1 and RQ3.

##### RQ1

We collected observational data from LR91 and used our causal DAG and the CTF
to perform identification and, in turn, design a statistical estimator capable
of making causal estimates. Therefore, we have demonstrated the use of
observational data to test for sensitivity with respect to different inputs.

##### RQ3

In order to perform sensitivity analysis of LR91, for each cause-effect
relationship of interest, we executed a series of causal test cases that
varied (normalised) inputs over the range [0, 1] whilst observing the change
in the output $APD_{90}$. The observational data used for this purpose was
sampled at random (via Latin Hyper Cube sampling) and, therefore, did not
contain the exact executions required for this exercise. Instead, our
estimators borrowed information from similar executions to predict the outcome
of the counterfactual test cases.

###  5.5 Covasim: Experimental Casual Testing

In this case study, we demonstrate the ability of the CTF to conduct
statistical metamorphic testing (SMT) of Covasim [3] _experimentally_. That
is, isolating the causal effect of interest via strategic executions of the
SUT, rather than applying graphical CI to observational data. Our aim here is
to provide evidence to support our claim that metamorphic testing is a
fundamentally causal activity that can be framed and solved as a problem of
CI.

####  5.5.1 Testing Activity

Revisiting our example from Section 3, our aim is to experimentally determine
the effect of prioritising vaccination for the elderly on the following
outputs: cumulative number of infections, number of doses given, maximum
number of doses per agent, and number of agents vaccinated.

Our expectation here is that prioritising the elderly should lead to an
increase in infections. This is because, by prioritising the elderly, we are
less likely to vaccinate agents in the model with a greater propensity for
spreading the virus (e.g. younger individuals who attend a school or
workplace). We also expect the number of vaccines and doses administered to
decrease as there are fewer elderly agents in the model. In contrast, the
maximum number of doses should not change, as the vaccine is set to be
administered at most two times per agent.

####  5.5.2 Data Generation

We executed the model under two input configurations 30 times each using an
experimental data collector (see Section 4) for every test. For both input
configurations, we used the default Covasim parameters, but fixed the
simulation length to 50 days, initial infected agents to 1000, population size
to 50,000, and made the default Pfizer vaccine available from day seven.
However, for the second configuration, we also sub-targeted (prioritised)
vaccination to the elderly using the vaccinate_by_age method from the Covasim
vaccination
tutorial.666https://github.com/InstituteforDiseaseModeling/covasim/blob/master/examples/t05_vaccine_subtargeting.py

####  5.5.3 Causal Testing

We first form a causal specification using the constraints and DAG from
Example 3.1. We then use the CTF to perform identification and testing. Since
we are executing the causal tests directly (and are therefore able to
explicitly control for potential biases), there is no confounding to adjust
for in the resulting data. As a consequence, we can directly calculate the ATE
by contrasting the average cumulative infections produced by the control
(vaccinate everyone) and treatment executions (prioritise the elderly).

####  5.5.4 Results

As expected, prioritising the elderly causes the cumulative infections to
increase (ATE: 2399.7, 95% CIs: [2323.7, 2475.8]) and causes no change to the
maximum doses (ATE: $8.9\mathrm{\times 10}^{-16}$, 95% CIs:
[$3.7\mathrm{\times 10}^{-17}$, $4.1\mathrm{\times 10}^{-16}$]).

However, when we examine the number of doses given (which we would expect to
remain fixed), the tests in fact show that the SUT erroneously causes the
number of doses administered and the number of people vaccinated to increase
sharply by 481351 (95% CIs: [480550, 482152]) and 483506 (95% CIs: [482646,
484367]), respectively. This is an obvious and potentially problematic bug, as
it reveals that more agents have been vaccinated than there are agents in the
simulation (by a factor of 9.7).

We raised an
issue777https://github.com/InstituteforDiseaseModeling/covasim/issues/370 on
Covasim’s GitHub repository to report this bug in September 2021 and the
Covasim developers replied in November confirming that the bug had been fixed
for version 3.1. Although the developers did not explain the cause of the bug
nor how it was fixed. The change log for version 3.1 stated the following:
_Rescaling now does not reset vaccination status; previously, dynamic
rescaling erased it._

Overall, this testing scenario has provided evidence to support the
experimental aspect of RQ1 and uncovered an interesting vaccine-related bug.

##### RQ1

In this testing scenario, we demonstrated an experimental approach to causal
testing in the CTF. More specifically, in order to measure the effect of
prioritising the elderly for vaccination, we used the CTF to execute a series
of tests that experimentally determined the causal effect of sub-targeting the
elderly on a number of vaccine-related outcomes.

###  5.6 Covasim: Observational Causal Testing

We now consider the effect of doubling infectiousness ($\beta$) on cumulative
infections, but this time using confounded observational data. In this case
study, we have two goals. First, we test the intuitive property that
increasing infectiousness should cause an increase in cumulative infections.
Then we turn our attention to more precise properties for which we have less
confident expectations, and thus a more ambiguous test oracle. Specifically,
we explore how the magnitude of the effect on cumulative infections varies
across populations with a specific age and density.

This form of exploratory testing [45, 46] necessitates a degree of exploration
and refinement during the testing process that, in a conventional setting,
requires numerous potentially costly test executions. To address this issue,
we demonstrate how the CTF can infer test outcomes from _existing, confounded
test data_ to explore uncertain model behaviour _without further model
executions_ in the absence of a precise test oracle.

####  5.6.1 Data Generation

When reasoning about infectiousness and the spread of COVID-19 using Covasim,
there are several parameters that can affect this output. These include the
variant of the virus and population characteristics such as age and household
size, with older populations being more susceptible to infection and higher
household contacts leading to quicker viral spread. In Covasim these
population characteristics cannot be specified directly, but can be indirectly
altered by selecting a geographical location.

For this scenario, we used four locations (Gambia, Mozambique, Bangladesh, and
Oman) which are particularly distinct in terms of average age and household
contacts. We wanted a collection of locations with both overlapping and non-
overlapping age demographics and household sizes to enable comparison of the
metamorphic test outcomes in combinations of these variables i.e. (many
contacts, low age), (many contacts, high age), (few contacts, low age), (few
contacts, high age) to improve our understanding of the model.

We generate three sets of data here. First, we directly apply a conventional
SMT approach to Covasim in which we execute the model 100 times with
$\beta=0.016$ and $\beta=0.032$ for each location, before averaging and
contrasting their respective cumulative infections. These results are shown in
Table 3 and capture the true location-specific metamorphic test outcomes
achieved by executing Covasim directly, as well as the average age and
household contacts of each location.

Second, we simulate (uncontrolled) observational data. Instead of directly
manipulating the $\beta$ parameter, we assign a different dominant variant to
each location (Alpha, Beta, Delta, Gamma), each of which has its own specific
$\beta$ value
($\beta_{\alpha}=0.027,\beta_{\beta}=0.016,\beta_{\delta}=0.035,\beta_{\gamma}=0.032)$.
The challenge from a testing perspective is to try to establish the effect of
doubling $\beta$, when there are several other confounding factors (location,
age, contacts) that could also affect the number of infections. Following this
approach, we generate 10,000 input configurations by uniformly sampling a
location and, based on that location, assign its dominant variant with $75\%$
probability and uniformly sampled one of the remaining variants with $25\%$
probability. We then run each input configuration once and record the outputs.

Finally, to introduce counterfactuals, we remove $\beta=0.032$ from the
observational data, meaning that there are no observed data that explicitly
include a doubling of the $\beta$ value. Therefore, any attempts to answer
questions concerning the impact of doubling $\beta$ would need to reason about
unobserved events – i.e _counterfactuals_.

Table 3: Metamorphic test results for the effect of doubling infectiousness ($\beta$) on cumulative infections in different locations in Covasim. Location | Age | Contacts | ATE (%) | 95% CIs
---|---|---|---|---
Bangladesh | 30.17 | 3.47 | 13.05 | [12.90, 13.20]
Gambia | 22.22 | 7.23 | 10.89 | [10.72, 11.06]
Mozambique | 22.26 | 3.37 | 20.62 | [20.37, 20.86]
Oman | 29.98 | 7.02 | 5.88 | [5.81, 5.96]

####  5.6.2 Causal Testing

To begin causal testing, we form our causal specification by specifying a
modelling scenario and the causal DAG shown in Figure 6. Our modelling
scenario uses the default Covasim parameters apart from $\beta$ (the input
under study) and the location. We also fixed the duration, population size,
initial infected agents, and scaling parameter as follows:

$\\{\mathtt{days}=216,\mathtt{pop\\_size}=51633,\mathtt{pop\\_infected}=1000,\mathtt{pop\\_scale}=20\\}$

Here, the $\mathtt{pop\\_scale}$ parameter scales the simulated population to
a size of $20\times 51633=1032660$, improving the efficiency of the simulation
at the cost of model fidelity.

$L$$\beta$$I$$A$$C$ Figure 6: A causal DAG for the Covasim modelling scenario
where the causal effect of $\beta$ on $I$ is confounded.

Next, applying the CTF to the causal DAG in Figure 6, we identify the
following sufficient adjustment sets:

$\\{A,C\\},\\{A,C,L\\},\\{L\\},\\{L,A\\},\\{L,C\\}$

While any of the above adjustment sets are suitable for identifying the effect
of $\beta$ on $I$, to borrow data from executions with overlapping age and
household contacts, we select $\\{A,C\\}$. This yields the following closed-
form statistical expression that is capable of directly estimating the causal
effect (ATE) of interest:

$\sum_{a\in A,\;c\in
C}ATE=\mathds{E}[I\mid\beta=0.032,A=a,C=c]-\mathds{E}[I\mid\beta=0.016,A=a,C=c]$

Then, to estimate the value of this estimand, we implement a regression model
of the form:

$I\sim x_{0}+x_{1}\beta+x_{2}\beta^{2}+x_{3}A+x_{4}C$

This regression model includes the quadratic term $\beta^{2}$ because
cumulative infections vary quadratically with $\beta$. This choice is informed
by domain expertise but can also be easily verified by varying the $\beta$
parameter and observing the resulting change in cumulative infections.

At this point, we have specified a causally-valid statistical model that is
capable of directly estimating the causal effect (ATE) of $\beta$ on
cumulative infections (our first goal). We can therefore fit the above
regression model to the data, substitute the values $\beta=0.016$ and
$\beta=0.032$, and contrast the respective estimates for $I$ to obtain the ATE
over the entire data-set.

While this statistical procedure enables us to determine whether doubling
$\beta$ causes the cumulative infections to increase, it cannot reveal how the
magnitude of this effect varies across populations with different age
demographics and population densities (our second goal). Instead, to achieve
this, we can change our causal metric to the conditional ATE (CATE), as
defined in Section 2. This metric allows us to examine causality within
specific populations instead of simply averaging over all data and, therefore,
to test more advanced properties of the SUT. Furthermore, this change of
metric does not require any additional executions of the SUT; we simply apply
our estimator to subsets of the data with a specific age and density, thereby
providing an efficient method for exploring uncertain behaviour using existing
data.

####  5.6.3 Results

Table 4: Causal Testing results for increasingly specific causal questions that pertain to different metamorphic properties, with associational (Rung 1), causal (Rung 2), and counterfactual (Rung 3) estimates provided as answers to each question. Causal Question |  Adj. Set | Age | Contacts | Rung | Runs | ATE (%) | 95% CIs
---|---|---|---|---|---|---|---
(1) What is the effect of increasing $\beta$ from 0.016 to 0.032 on cumulative infections? | $\oslash$ | - | - | 1 | 10000 | 19.67 | [19.51, 19.84]
$\\{A,C\\}$ | 2 | 10000 | 14.65 | [14.54, 14.76]
$\\{A,C\\}$ | 3 | 7522 | 14.95 | [14.84, 15.06]
(2) What is the effect of increasing $\beta$ from 0.016 to 0.032 on cumulative infections in a population with a lower age demographic? | $\oslash$ | L | - | 1 | 4972 | 21.46 | [21.25, 21.66]
$\\{A,C\\}$ | 2 | 4972 | 16.63 | [16.54, 16.73]
$\\{A,C\\}$ | 3 | 4552 | 16.76 | [16.66, 16.86]
(3) What is the effect of increasing $\beta$ from 0.016 to 0.032 on cumulative infections in a population with an older age demographic? | $\oslash$ | H | - | 1 | 5028 | 9.38 | [9.21, 9.55]
$\\{A,C\\}$ | 2 | 5028 | 8.46 | [8.38, 8.55]
$\\{A,C\\}$ | 3 | 2970 | 8.56 | [8.48, 8.65]
(4) What is the effect of increasing $\beta$ from 0.016 to 0.032 on cumulative infections in a population with a lower age demographic and few household contacts? | $\oslash$ | L | L | 1 | 2536 | 20.38 | [20.33, 20.43]
$\\{A,C\\}$ | 2 | 2536 | 20.37 | [20.32, 20.42]
$\\{A,C\\}$ | 3 | 2327 | 20.65 | [20.61, 20.70]
(5) What is the effect of increasing $\beta$ from 0.016 to 0.032 on cumulative infections in a population with a lower age demographic and many household contacts? | $\oslash$ | L | H | 1 | 2436 | 10.15 | [10.12, 10.18]
$\\{A,C\\}$ | 2 | 2436 | 10.15 | [10.12, 10.18]
$\\{A,C\\}$ | 3 | 2225 | 10.17 | [10.14, 10.20]
(6) What is the effect of increasing $\beta$ from 0.016 to 0.032 on cumulative infections in a population with an older age demographic and few household contacts? | $\oslash$ | H | L | 1 | 2564 | 11.79 | [11.76, 11.82]
$\\{A,C\\}$ | 2 | 2564 | 11.79 | [11.77, 11.82]
$\\{A,C\\}$ | 3 | 2345 | 11.80 | [11.77, 11.83]
(7) What is the effect of increasing $\beta$ from 0.016 to 0.032 on cumulative infections in a population with an older age demographic and many household contacts? | $\oslash$ | H | H | 1 | 2464 | 5.61 | [5.59, 5.62]
$\\{A,C\\}$ | 2 | 2464 | 5.61 | [5.59, 5.62]
$\\{A,C\\}$ | 3 | 625 | 5.70 | [5.68, 5.71]

Table 4 summarises the results of applying the CTF to Covasim to understand
the effect of doubling infectiousness ($\beta$) on cumulative infections in
various populations. This table contains a series of increasingly specific
causal questions that pertain to increasingly precise metamorphic relations.
Each causal question has been answered using three approaches corresponding to
the rungs of the causal ladder (Section 2): association (simple regression
with no adjustment), intervention (regression model adjusted for identified
confounders), and counterfactual (executions with $\beta=0.032$ have been
removed from data).

The first row of Table 4 (Causal Question 1, Rung 2) shows that doubling
$\beta$ from 0.016 to 0.032 _causes_ an increase to infections of 14.65% with
95% confidence intervals of (14.54%, 14.76%). For comparison, we also include
a naïve estimate (Rung 1), calculated by simply regressing $I$ against
$\beta+\beta^{2}$. The causal estimate is almost 5% smaller than the
associational one, indicating that naïve regression leads to an overestimate
of the causal effect in this case. This reflects the difference between a
statistical prediction (Rung 1) and a causal estimate (Rungs 2 and 3); the
former is blind to the source of non-causal ‘noise’ while the latter
recognises and adjusts for this. In situations where a more precise oracle is
available (e.g. we expect to observe a specific change in outcome), this error
could result in misleading test outcomes. Additionally, our counterfactual
result (Rung 3), is similar to our causal result (Rung 2), indicating that the
CTF is able to operate in the counterfactual setting for this test case with
only a small loss of accuracy.

Table 4 also shows the CATEs computed for different age strata (Causal
Questions 2 and 3), and different combinations of age and household contacts
(Causal Questions 4-7). The results for questions 2 and 3 indicate that the
effect is more pronounced in lower age groups than in higher ones (16.67%
versus 8.46%). For question 2 in particular there is a notable difference of
4.83% between the estimate computed by adjusting for age and the associational
estimate. Results for questions 4-7 indicate that, within each age stratum,
populations with many household contacts are affected to a lesser extent than
populations with few household contacts. For these cases, the difference
between the causal and associational estimates (Rungs 1 and 2) is negligible.
This is because the choice of metric adjusts for age and household contacts by
design. As with question 1, our Rung 3 estimates are very close to those of
Rung 2, indicating that the CTF copes well in the counterfactual case for this
case study.

Furthermore, since Causal Questions 4-7 correspond to an individual location
(defined by a unique age-contact stratum), we can compare the inferred test
outcomes to the metamorphic test outcomes in Table 3. We achieve similar
results for most of the test cases, however our aim here is not to evaluate
the accuracy of our statistical models (which could potentially be improved by
more advanced statistical models, such as Causal Forest Estimators [65]), but
rather to demonstrate the potential of a CI-driven approach to testing
scientific modelling software by re-using test data.

Overall, this testing scenario has provided evidence to answer all research
questions in the affirmative.

##### RQ1

In this testing scenario, we demonstrated an observational approach to causal
testing. To understand the effect of increasing the infectiousness ($\beta$)
of a virus on cumulative infections, we collected observational data from
Covasim and used the CTF to infer test outcomes retrospectively via graphical
CI.

##### RQ2

In this case study, our observational data contained confounding bias, as
shown in Figure 6. Therefore, purely statistical techniques (rung one) would
be led astray in this situation (as shown in Table 4) to reach biased test
outcomes. However, using the CTF, we were able to identify age and household
contacts as a sufficient adjustment set (see Section 2.4) and adjust for their
biasing affect with an appropriate statistical model. This enabled us to
recover the causal effects of interest from the confounded observational data.

##### RQ3

To demonstrate the ability of the CTF to infer counterfactual test outcomes,
we repeated the analysis for each causal question after removing all
executions with $\beta=0.032$ from the data. As a result, these
(counterfactual) causal test cases had to borrow information from similar
executions to infer unobserved test outcomes (i.e. outcomes where
$\beta=0.032$). The results for these counterfactual inferences can be seen in
the third row of each Causal Question in Table 4 and demonstrate that, in this
case study, we could accurately infer counterfactual test outcomes.

##  6 Discussion

Section 5 outlined the research questions and success criteria for our case
studies, measuring the extent to which the Causal Testing Framework (CTF)
delivers the theoretical advantages stated in Section 3. We will now discuss
our findings for these research questions and a pair of interesting bugs found
in the process.

###  6.1 RQ1: Can we apply Causal Testing both experimentally and
observationally?

Across our case studies, we demonstrated both an experimental and
observational application of the CTF. In Section 5.5, for example, we designed
causal test cases to evaluate the effects of prioritising the elderly for
vaccination and then executed the model _directly_ to isolate the causal
effects of interest. In contrast, we then generated confounded observational
data in Section 5.6 and inferred from this the causal effect of doubling
infectiousness ($\beta$) on cumulative infections.

In practice, this functionality makes it possible for testers to establish and
reason about the salient causal relationships in software systems in two
distinct ways. First, as with many conventional testing methods, the CTF is
able to run the SUT _directly_ to examine causal behaviour. Second, when
practical constraints prevent further tests from being executed, the CTF is
able to infer test outcomes from existing test data. Hence, this observational
mode of testing is sympathetic to the challenging properties of scientific
modelling software, such as long execution times and non-determinism, that
make them costly and impractical to test experimentally.

The Causal Testing Framework can test properties using observational as well
as more conventional experimental data. Therefore, causal testing can be
carried out using existing test data instead of running the program.

An open problem related to RQ1 is to understand the properties and
requirements of observational data that are necessary to achieve reliable and
accurate inferences. Intuitively, we cannot infer a specific test outcome from
arbitrary data; it must have some _overlap_ with the target testing scenario.
Hence, a natural question is whether there is a minimal amount of overlap that
would be sufficient for the purposes of testing scientific software. If so,
future work could use this as a target for test generation or selection
techniques in order to produce test suites that, once executed, yield an
optimal set of data from which precise inferences can be made.

###  6.2 RQ2: Can Causal Testing mitigate confounding bias in observational
test data?

In Section 5.6, we inferred the causal effect of doubling infectiousness on
cumulative infections from _confounded data_. We achieved this by applying
graphical CI techniques to our causal DAG to identify and then adjust for age
and household contacts as confounders, revealing that naive statistical
prediction overestimated the effect by almost $5\%$.

This demonstrates the potential of the CTF to relax the conditions ordinarily
associated with testing (the need to run the model repeatedly under carefully
controlled conditions), while preserving the ability to draw causal
inferences. Hence, whereas a conventional approach requires carefully
controlled executions to draw causal conclusions, the CTF can achieve this
using existing, confounded data. Critically, this requires careful analysis of
the causal structure of the scenario that would not be possible using
statistical models alone, and thus represents an important contribution of the
CTF.

The Causal Testing Framework can identify and mitigate sources of confounding
bias in observational test data. This enables testing to be carried out
retroactively using biased execution data in situations where it is
impractical or even impossible to repeatedly run the SUT.

A caveat to our approach, as with any model-based testing technique, is that
we require the causal DAG to be produced by the domain expert. While this is
an intuitive and widely accepted practice in fields such as epidemiology and
social sciences [18], it still requires manual effort and domain expertise.
However, the field of causal discovery (CD) [66] provides a range of tools and
techniques for automatically learning causal structures, such as causal DAGs,
from observational data. Hence, to improve the usability of the CTF, future
work could investigate the potential of applying existing CD algorithms to
software execution data to automatically construct causal specifications.

###  6.3 RQ3: Can Causal Testing handle counterfactual test outcomes?

Throughout our case studies, we employed causally-valid statistical estimators
to infer counterfactual test outcomes. We achieved this by combining domain
expertise (for example, the decision to include quadratic terms in estimators)
with graphical CI techniques. This process enabled us to identify and adjust
for confounding and thus borrow information from similar test executions to
draw causal conclusions about the effect of unobserved interventions. In
Section 5.6, for example, we emphasised this point by deleting all executions
from the data containing $\beta=0.032$ before repeating our causal test cases
to measure the causal effect of increasing $\beta$ from $0.016$ to $0.032$.

From a testing perspective, our case studies have shown how counterfactual
inference can offer an efficient alternative to conventional testing methods
in situations where practical constraints exist. Specifically, by leveraging
causal assumptions and domain expertise in the form of a statistical model,
the tester can predict how the model would behave in unseen modelling
scenarios. Hence, where it is impractical to run the model directly for
testing purposes, counterfactual inference can be employed as an economical
alternative.

By leveraging the causal assumptions encoded in a causal DAG and appropriate
statistical estimators, the Causal Testing Framework can “borrow” data from
similar software executions to infer counterfactual test outcomes. This makes
is possible to learn test outcomes without possessing data that cover the
exact parameter values of interest.

Counterfactual inferences are only as valid as the functional form of the
specified statistical model. The process of specifying such a model typically
depends on domain expertise that, for systems exploring uncertain phenomena,
are inherently elusive. However, there exists a significant amount of work
into estimation models that help to relax these constraints, such as causal
forests [65] and doubly robust estimators [67]. Future work could investigate
the suitability of more advanced statistical estimators such as these for
handling counterfactual inference in the absence of sufficient domain
knowledge.

###  6.4 Additional Findings

Overall, our case studies have provided affirmative answers to RQ1, RQ2, and
RQ3, demonstrating how the CTF can realise the theoretical advantages of a CI-
led approach to software testing, as outlined in Section 3. We also identified
a number of directions for future work that would improve the practicality of
the CTF.

In addition, our case studies revealed two further findings. First, by
leveraging the advantages of CI in a testing context, we were able to
efficiently explore the uncertain behaviour of the subject systems. Second,
two of our case studies led to the discovery of interesting bugs. We briefly
discuss these findings below.

##### Exploratory Testing

Across our case studies, we applied the CTF to existing model execution data
to test additional properties without further model executions. In Section
5.3, for example, we used the CTF to conduct statistical metamorphic testing
using less data than a conventional approach. We then used the _same data_ to
test a second property without additional model runs.

Similarly, in Section 5.6, using a fixed set of model executions, we applied
the CTF to answer increasingly specific causal questions corresponding to more
precise metamorphic relations. In practice, this enabled us to employ an
economic approach to exploratory testing in order to better understand the
constituent impacts of age and household contacts on cumulative infections,
without running the model further times.

This finding suggests the CTF offers an efficient solution for testing
scientific modelling software through careful causal analysis of _existing
test data_ , rather than further costly model executions. It is this emphasis
on maximising what we can learn from existing data that makes the CTF suitable
for scientific modelling software.

##### Bugs Found

Our case studies also revealed two interesting, previously undiscovered bugs
in two of the studied scientific models: the Poisson Line Tessellation model
and Covasim.

First, in Section 5.3, we found that the relationship between intensity and
number of polygons per unit area described in [31] was more fragile at smaller
window sizes. This suggested that the window size (width and/or height) has a
causal effect on the number of polygons per unit area, while [31] stated that
these variables should be independent. We then designed a causal test case to
confirm that increasing the window width from 1 to 2 whilst holding intensity
constant _caused_ a significant change in he number of polygons per unit area.

Second, in Section 5.5, we found a bug in Covasim’s vaccine implementation
where, upon prioritising elderly for vaccination, the number of vaccinated
individuals grew to nearly ten times the number of individuals in the
simulation. While this does not appear to significantly affect the key outputs
of the model, it is not difficult to imagine how such a bug could lead to an
overestimation of the effects of particular interventions.

###  6.5 Threats to Validity

Our evaluation in Section 5 cannot be used, and does not claim, to make
general conclusions about the accuracy and scalability of the CTF. Instead, it
is intended to serve as a proof of concept that causality-focused software
testing activities can be framed and solved as problems of CI, and therefore
demonstrate the potential for graphical CI to be applied to software testing
problems. Nonetheless, there are some threats to validity worth considering
here.

####  6.5.1 External Validity

In this work, the main threat to external validity is that our case studies
only cover three subject systems involving a moderate number of input and
output variables. As graphical CI requires domain expertise for the data
generating mechanism in the form of a causal DAG, a significant amount of time
was spent familiarising ourselves with the subject systems and understanding
their constituent cause-effect relationships. As a result, this limited our
ability to systematically collect and analyse large numbers of varied subject
systems.

Furthermore, our subject systems were all implemented in Python. Therefore,
our findings do not necessarily generalise to scientific modelling software
implemented in other languages. However, the CTF only requires execution data
in CSV format to perform causal testing observationally and can thus be
applied, in theory, to tabular data produced by _any_ scientific model.

As a consequence of the aforementioned threats to external validity, we
acknowledge that our results may not generalise to _all_ forms of scientific
modelling software. However, we attempt to mitigate the aforementioned threats
to external validity by selecting models that differ in their complexity,
subject-matter, and modelling paradigm. In addition, as discussed in Section
5.2, the selected systems have important but vastly different applications
across a variety domains, and have all been the subject of prior research.

####  6.5.2 Internal Validity

In this paper, the main threat to internal validity is that we did not
optimise the estimators and configuration parameters thereof for our case
studies. While this avoids the problem of over-fitting, it means there may
exist statistical models that are more suitable for modelling and inferring
the behaviour of the input-output relationships under study.

In the same vein, when handling counterfactual estimation, we specified
regression equations that capture the expected functional form of various
input-output relationships. For example, when testing Covasim in Section 5.6,
we specified a regression model in which cumulative infections vary
quadratically with infectiousness. We called upon our experience with the
models and subject area to specify these equations. However, different domain
experts may have different opinions about the correct functional forms of the
input-out relationships and may therefore have specified these relationships
differently.

As a consequence of the above threats to internal validity, we acknowledge
that there may exist alternative statistical models that can achieve more
precise causal and counterfactual inferences for the subject systems. However,
we partially mitigate the above threats to internal validity by manually
inspecting the functional forms of the relationships between inputs and
outputs of interest in the SUT. We achieve this by varying one parameter at a
time and observing how the output in question changes in response (in a
similar way to our sensitivity analysis case study in Section 5.4).

##  7 Related Work

In this section, we provide a brief review of work related to the two main
topics concerning our paper: approaches for testing scientific software and
causality in software testing. Additionally, we summarise automatic approaches
to generating causal DAGs and highlight a number of open research challenges.

###  7.1 Testing Techniques for Scientific Software

As stated in Kanewala and Bieman’s survey [4], scientific models are seldom
tested using systematic approaches. Instead, techniques such as sensitivity
[68] and uncertainty analysis [69] are often employed to analyse and appraise
scientific models. However, these approaches generally require many costly
executions that make them prohibitively expensive at scale [70]. To address
this issue, modellers have turned to emulator approaches [62, 70], where a
surrogate model is developed to approximate the behaviour of the simulation
and provide an efficient way to validate behaviour [61, 71]. However, these
emulators are driven by statistical associations and are unable to draw
_causal_ inferences from existing test data.

Another issue that precludes the testing of scientific modelling software is
the oracle problem [6]; the lack of a mechanism that can be used to ascertain
whether the outcome of a test case is correct or not. Kanewala and Bieman’s
survey [4] identifies several approaches followed by scientific modellers to
overcome the oracle problem, including: pseudo oracles, comparison to
analytical solutions or experimental results, and expert judgement. In
addition to these solutions, modellers have also turned to metamorphic testing
(see Section 2) to overcome the lack of oracle. This approach relies on the
scientists being able to specify metamorphic relations capable of revealing
faults. However, these relationships are notoriously challenging to identify
[30].

To assist with the identification of metamorphic relations, Kanewala and
Bieman developed a machine learning approach for predicting metamorphic
relations for numerical software [72]. This is achieved by representing
numerical functions as a statement-level control flow graph and extracting
features from this graph to train a classifier. In recent years, several new
approaches for automatically predicting metamorphic relations for a specific
form of software have been proposed, including for cyber-physical systems [73,
74] and matrix calculation programs [75]. However, the generation of
metamorphic relations remains a difficult problem with automatic solutions
available for only a few specific forms of software.

###  7.2 Causality in Software Testing

In more conventional settings, CI techniques have been applied to the software
testing problem of fault localisation (FL). Informally, FL concerns
identifying locations of faults in a program [76] and often involves computing
a “suspiciousness metric” for software components, such as program statements.
However, these metrics are often confounded by other software components. To
address this, Baah et al. [77] translated FL to a CI problem, using program
dependence graphs as a model of causality to estimate the causal effects of
program statements on the occurrence of faults. Subsequent papers build on
this to handle additional sources of bias [78]; leverage more advanced
statistical models [78, 79]; and adapt to different software components [80,
81, 82, 79].

More recently, Lee et al. have introduced the Causal Program Dependence
Analysis Framework and applied it to FL. This is a CI-driven framework that
measures the strength of dependence between program elements by modelling
their causal structure [83]. Unlike the previous CI-based FL techniques, this
framework does not use static analysis to construct its underlying causal
model, and instead approximates the causal structure from existing,
observational data. Through a series of experiments, the framework has been
shown to outperform slicing-based and search-based FL techniques, and, due to
its focus on dependence relations instead of coverage, is less susceptible to
coincidental correctness (executions that pass but cover faulty components).

In a similar vein, software testing often involves understanding _why_ a
particular outcome occurs, such as a program failure. To this end, Johnson et
al. [84], developed a tool that explains the root cause of faulty software
behaviour. This tool creates “causal experiments” by mutating an existing test
to form a suite of minimally different tests that, contrary to the original,
are not fault-causing. The passing and failing tests can then be compared to
understand _why_ a fault occurred. Similarly, Chockler et al. [85] developed a
tool to _explain_ the decisions of deep neural network (DNN) image
classifiers. Following the actual causes framework [86], this tool offers
explanations in the form of minimal subsets of pixels sufficient for the DNN
to classify an image.

Another software testing technique concerning causality is cause-effect
graphing, a black-box approach adapted from hardware testing. Here, input-
output relationships are expressed in a variant of a combinatorial logic
network, known as a cause-effect graph, created by manually extracting causes,
effects, and boolean constraints from natural language specifications [87,
88]. Unlike the previous techniques, this approach does not use CI.

###  7.3 Automatic Generation of Causal DAGs

In this paper, we have assumed that all causal DAGs are specified manually by
a domain expert. While this is an intuitive and widely accepted approach for
conducting CI in fields such as epidemiology and social sciences, there are
two potential methods that could, in theory, (partially) automate this
process.

First, under certain strict assumptions and with large quantities of data, it
is possible to predict the structure of causal DAGs from observational data.
Where _model inference_ provides a source of models for traditional MBT
techniques [10], the field of _causal discovery_ (CD) [89] provides methods to
infer causal structures from data by exploiting asymmetries that distinguish
association from causation [66]. However, due to the need for large amounts of
data and their strict assumptions, we have had limited success in applying CD
algorithms to model execution data. We plan to investigate this route further
in future work.

Second, causal DAGs can be generated via static analysis of source code. DAGs
derived in this way have already been used for FL [79, 83]. However, this
approach relies on source code being openly available and produces a detailed,
low-level causal model of the SUT. While this level of granularity is ideal
for the purpose of FL, the resulting causal model would be less suitable for a
typical scientific modeller.

In addition to the aforementioned challenges, there is a fundamental barrier
to using automatically generated causal models for testing: inferred models
represent the implemented system rather than the true specification. In other
words, even if we could perfectly recover the DAG of the implementation, this
would contain any bugs the implementation may have. We would, in effect, be
testing the system against itself, so it would trivially look correct. Hence,
the correctness of any inferred DAGs must be verified by a domain expert.

##  8 Conclusion and Future Work

In this paper, we presented the Causal Testing Framework: a conceptual
framework that facilitates the application of causal inference techniques to
software testing problems. This framework follows a model-based testing
approach to incorporate an explicit model of causality into the software
testing process in the form of a causal DAG, enabling the direct application
of graphical causal inference methods to software testing activities. In
addition, we provided an open source reference implementation of the Causal
Testing Framework.

A key contribution of the Causal Testing Framework is the ability to infer
test outcomes from previous execution data, despite the presence of
confounding and counterfactuals, to provide an efficient method for testing
scientific models. To demonstrate this, we applied our framework to three
real-world scientific models of varying size and complexity, including a
Poisson line tessellation model, a cardiac action potential model, and an
epidemiological agent-based model. These case studies indicate that a causal
inference-led testing approach can reduce the amount of data necessary to test
a model and thus offers a more efficient solution for testing scientific
modelling software, particularly in the absence of a precise test oracle.
Furthermore, we identified a pair of notable bugs in two of the subject
systems: a fragile metamorphic relation and a vaccination-related error.

Software testing is an inherently causal process, and the field of causal
inference holds much-untapped potential. To this end, the Causal Testing
Framework lays the foundation for a new line of causality-driven software
testing techniques. In one line of future work, we plan to investigate the
overlap between causal inference and metamorphic testing with the aim of using
causal DAGs as a generative model for _causal_ metamorphic relations. This
could form the basis for novel test case generation and test adequacy metrics
based around causality. A separate direction of research would be to establish
a (semi-)automatic, reliable process for the discovery of causal DAGs of
software systems. Despite the issues surrounding automatic causal DAG
generation discussed in Section 7, it remains an attractive prospect, and
would lower the barrier to entry and improve the usability of the Causal
Testing Framework.

## Acknowledgements

Foster, Walkinshaw, and Hierons are funded by the EPSRC CITCoM grant
EP/T030526/1. For the purpose of open access, the author has applied a
Creative Commons Attribution (CC BY)888Where permitted by UKRI a CC-BY-ND
licence may be stated instead. licence to any Author Accepted Manuscript
version arising.

## References

  * [1] M. Oldfield and E. Haig, “Analytical modelling and uk government policy,” AI and Ethics, pp. 1–16, 2021.
  * [2] R. N. Thompson, “Epidemiological models are important tools for guiding covid-19 interventions,” BMC medicine, vol. 18, no. 1, pp. 1–4, 2020.
  * [3] C. C. Kerr, R. M. Stuart, D. Mistry, R. G. Abeysuriya, K. Rosenfeld, G. R. Hart, R. C. Núñez, J. A. Cohen, P. Selvaraj, B. Hagedorn, et al., “Covasim: an agent-based model of COVID-19 dynamics and interventions,” PLOS Computational Biology, vol. 17, no. 7, p. e1009149, 2021.
  * [4] U. Kanewala and J. M. Bieman, “Testing scientific software: A systematic literature review,” Information and software technology, vol. 56, no. 10, pp. 1219–1232, 2014.
  * [5] D. Kelly and R. Sanders, “The challenge of testing scientific software,” in Proceedings of the 3rd annual conference of the Association for Software Testing (CAST 2008: Beyond the Boundaries), pp. 30–36, Citeseer, 2008.
  * [6] E. T. Barr, M. Harman, P. McMinn, M. Shahbaz, and S. Yoo, “The oracle problem in software testing: A survey,” IEEE Transactions on Software Engineering, vol. 41, no. 5, pp. 507–525, 2015.
  * [7] J. Cornfield, W. Haenszel, E. C. Hammond, A. M. Lilienfeld, M. B. Shimkin, and E. L. Wynder, “Smoking and Lung Cancer: Recent Evidence and a Discussion of Some Questions,” JNCI: Journal of the National Cancer Institute, vol. 22, pp. 173–203, 01 1959.
  * [8] E. F. Moore, “Gedanken-experiments,” in Automata Studies (C. Shannon and J. McCarthy, eds.), Princeton University Press, 1956.
  * [9] R. M. Hierons, K. Bogdanov, J. P. Bowen, R. Cleaveland, J. Derrick, J. Dick, M. Gheorghe, M. Harman, K. Kapoor, P. Krause, G. Lüttgen, A. J. H. Simons, S. A. Vilkomir, M. R. Woodward, and H. Zedan, “Using formal specifications to support testing,” ACM Computing Surveys, vol. 41, no. 2, pp. 9:1–9:76, 2009.
  * [10] M. Utting, A. Pretschner, and B. Legeard, “A taxonomy of model-based testing approaches,” Software Testing, Verification and Reliability, vol. 22, no. 5, pp. 297–312, 2012.
  * [11] D. Lee and M. Yannakakis, “Principles and methods of testing finite-state machines - a survey,” Proceedings of the IEEE, vol. 84, no. 8, pp. 1089–1123, 1996.
  * [12] J. Tretmans, “Model based testing with labelled transition systems,” in Formal Methods and Testing, vol. 4949 of Lecture Notes in Computer Science, pp. 1–38, Springer, 2008.
  * [13] R. M. Hierons, “Testing from a Z specification,” The Journal of Software Testing, Verification and Reliability, vol. 7, no. 1, pp. 19–33, 1997\.
  * [14] J. Dick and A. Faivre, “Automating the generation and sequencing of test cases from model-based specifications,” in FME ’93, First International Symposium on Formal Methods in Europe, (Odense, Denmark), pp. 268–284, Springer-Verlag, Lecture Notes in Computer Science 670, 19-23 April 1993.
  * [15] J. Pearl, “Causal diagrams for empirical research,” Biometrika, vol. 82, pp. 669–688, 12 1995.
  * [16] M. A. Hernán and J. M. Robins, Causal Inference: what if.  Boca Raton: Chapman & Hall/CRC, 2020.
  * [17] S. Greenland, J. Pearl, and J. M. Robins, “Causal diagrams for epidemiologic research,” Epidemiology, pp. 37–48, 1999.
  * [18] P. W. Tennant, E. J. Murray, K. F. Arnold, L. Berrie, M. P. Fox, S. C. Gadd, W. J. Harrison, C. Keeble, L. R. Ranker, J. Textor, et al., “Use of directed acyclic graphs (DAGs) to identify confounders in applied health research: review and recommendations,” International Journal of Epidemiology, vol. 50, no. 2, pp. 620–632, 2021.
  * [19] T. Y. Chen, S. C. Cheung, and S. M. Yiu, “Metamorphic testing: A new approach for generating next test cases,” Tech. Rep. HKUST-CS98-01, The Hong Kong University of Science and Technology, 1998.
  * [20] S. Nidhra and J. Dondeti, “Black box and white box testing techniques-a literature review,” International Journal of Embedded Systems and Applications (IJESA), vol. 2, no. 2, pp. 29–50, 2012.
  * [21] J. Pearl, “Does obesity shorten life? or is it the soda? on non-manipulable causes,” Journal of Causal Inference, vol. 6, no. 2, 2018.
  * [22] K. Kreyman, D. L. Parnas, and S. Qiao, Inspection procedures for critical programs that model physical phenomena.  Citeseer, 1999.
  * [23] I. for Disease Modelling, “Covasim.” https://github.com/InstituteforDiseaseModeling/covasim, 2022.
  * [24] C. C. Kerr, D. Mistry, R. M. Stuart, K. Rosenfeld, G. R. Hart, R. C. Núñez, J. A. Cohen, P. Selvaraj, R. G. Abeysuriya, M. Jastrzębski, et al., “Controlling COVID-19 via test-trace-quarantine,” Nature Communications, vol. 12, no. 1, pp. 1–12, 2021.
  * [25] J. A. Cohen, D. Mistry, C. C. Kerr, and D. J. Klein, “Schools are not islands: Balancing COVID-19 risk and educational benefits using structural and temporal countermeasures,” medRxiv, 2020.
  * [26] J. Panovska-Griffiths, C. C. Kerr, R. M. Stuart, D. Mistry, D. J. Klein, R. M. Viner, and C. Bonell, “Determining the optimal strategy for reopening schools, the impact of test and trace interventions, and the risk of occurrence of a second COVID-19 epidemic wave in the uk: a modelling study,” The Lancet Child & Adolescent Health, vol. 4, no. 11, pp. 817–827, 2020.
  * [27] N. Scott, A. Palmer, D. Delport, R. Abeysuriya, R. Stuart, C. C. Kerr, D. Mistry, D. J. Klein, R. Sacks-Davis, K. Heath, et al., “Modelling the impact of reducing control measures on the COVID-19 pandemic in a low transmission setting,” Med J Aust, vol. 214, no. 2, pp. 79–83, 2020.
  * [28] E. Weyuker, “On testing non-testable programs,” Computer Journal, vol. 25, 11 1982.
  * [29] I. for Disease Modeling, “Covasim: Vaccine tests.” https://github.com/InstituteforDiseaseModeling/covasim/blob/master/tests/test_interventions.py, 2022\.
  * [30] S. Segura, G. Fraser, A. B. Sanchez, and A. Ruiz-Cortés, “A survey on metamorphic testing,” IEEE Transactions on software engineering, vol. 42, no. 9, pp. 805–824, 2016.
  * [31] R. Guderlei and J. Mayer, “Statistical metamorphic testing testing programs with random output by means of statistical hypothesis tests and metamorphic testing,” in Seventh International Conference on Quality Software (QSIC 2007), pp. 404–409, 2007.
  * [32] L. Keele, “The statistics of causal inference: A view from political methodology,” Political Analysis, pp. 313–335, 2015.
  * [33] J. Pearl, “Causal inference in statistics: An overview,” Statistics Surveys, vol. 3, pp. 96–146, 2009.
  * [34] J. Pearl and D. Mackenzie, The Book of Why.  Allen Lane, 2018.
  * [35] P. W. Holland, “Statistics and causal inference,” Journal of the American statistical Association, vol. 81, no. 396, pp. 945–960, 1986.
  * [36] J. Abrevaya, Y.-C. Hsu, and R. P. Lieli, “Estimating conditional average treatment effects,” Journal of Business & Economic Statistics, vol. 33, no. 4, pp. 485–505, 2015.
  * [37] S. F. O’Brien and Q. L. Yi, “How do i interpret a confidence interval?,” Transfusion, vol. 56, no. 7, pp. 1680–1683, 2016.
  * [38] K. J. Rothman and S. Greenland, “Causation and causal inference in epidemiology,” American journal of public health, vol. 95, no. S1, pp. S144–S150, 2005.
  * [39] J. Pearl, Causality.  Cambridge university press, 09 2009.
  * [40] D. B. Rubin, “Causal inference using potential outcomes: Design, modeling, decisions,” Journal of the American Statistical Association, vol. 100, no. 469, pp. 322–331, 2005.
  * [41] R. B. Kline, Principles and practice of structural equation modeling.  Guilford publications, 2015.
  * [42] M. Staats, M. W. Whalen, and M. P. Heimdahl, “Programs, tests, and oracles: the foundations of testing revisited,” in 2011 33rd international conference on software engineering (ICSE), pp. 391–400, IEEE, 2011.
  * [43] J. Pearl and T. S. Verma, “A theory of inferred causation,” in Studies in Logic and the Foundations of Mathematics, vol. 134, pp. 789–811, Elsevier, 1995.
  * [44] G. C. Murphy, D. Notkin, and K. Sullivan, “Software reflexion models: Bridging the gap between source and high-level models,” in Proceedings of the 3rd ACM SIGSOFT symposium on Foundations of software engineering, pp. 18–28, 1995.
  * [45] J. Itkonen and K. Rautiainen, “Exploratory testing: a multiple case study,” in 2005 International Symposium on Empirical Software Engineering, 2005., pp. 10–pp, IEEE, 2005.
  * [46] X. Lin, M. Simon, and N. Niu, “Exploratory metamorphic testing for scientific software,” Computing in science & engineering, vol. 22, no. 2, pp. 78–87, 2018.
  * [47] J. Ellson, E. Gansner, L. Koutsofios, S. C. North, and G. Woodhull, “Graphviz— open source graph drawing tools,” in Graph Drawing (P. Mutzel, M. Jünger, and S. Leipert, eds.), (Berlin, Heidelberg), pp. 483–484, Springer Berlin Heidelberg, 2002.
  * [48] L. de Moura and N. Bjørner, “Z3: An efficient SMT solver,” in Tools and Algorithms for the Construction and Analysis of Systems, pp. 337–340, Springer Berlin Heidelberg, 2008.
  * [49] S. Wager and S. Athey, “Estimation and inference of heterogeneous treatment effects using random forests,” Journal of the American Statistical Association, vol. 113, no. 523, pp. 1228–1242, 2018.
  * [50] V. V. Chetlur and H. S. Dhillon, “Coverage analysis of a vehicular network modeled as cox process driven by poisson line process,” IEEE Transactions on Wireless Communications, vol. 17, no. 7, pp. 4401–4416, 2018\.
  * [51] F. Morlot, “A population model based on a poisson line tessellation,” in 2012 10th International Symposium on Modeling and Optimization in Mobile, Ad Hoc and Wireless Networks (WiOpt), pp. 337–342, 2012.
  * [52] D. Shin, A. Kirmani, A. Colaço, and V. K. Goyal, “Parametric poisson process imaging,” in 2013 IEEE Global Conference on Signal and Information Processing, pp. 1053–1056, IEEE, 2013.
  * [53] C.-H. Luo and Y. Rudy, “A model of the ventricular cardiac action potential. Depolarization, repolarization, and their interaction.,” Circulation Research, vol. 68, no. 6, pp. 1501–1526, 1991.
  * [54] J. Panovska-Griffiths, C. C. Kerr, W. Waites, R. M. Stuart, D. Mistry, D. Foster, D. J. Klein, R. M. Viner, and C. Bonell, “The potential contribution of face coverings to the control of sars-cov-2 transmission in schools and broader society in the uk: a modelling study,” 2020.
  * [55] R. M. Stuart, R. G. Abeysuriya, C. C. Kerr, D. Mistry, D. J. Klein, R. Gray, M. Hellard, and N. Scott, “The role of masks in reducing the risk of new waves of covid-19 in low transmission settings: a modeling study,” MedRxiv, 2020.
  * [56] J. L. Deutsch and C. V. Deutsch, “Latin hypercube sampling with multidimensional uniformity,” Journal of Statistical Planning and Inference, vol. 142, no. 3, pp. 763–772, 2012.
  * [57] S. Moza, “sahilm89/lhsmdu: Latin Hypercube Sampling with Multi-Dimensional Uniformity (LHSMDU): Speed Boost minor compatibility fixes,” jul 2020.
  * [58] S. N. Chiu, D. Stoyan, W. S. Kendall, and J. Mecke, Stochastic Geometry and its Applications.  Chichester, West Sussex, United Kingdom: John Wiley & Sons Inc, 3rd ed., 2013.
  * [59] J. P. Kleijnen, “Verification and validation of simulation models,” European journal of operational research, vol. 82, no. 1, pp. 145–162, 1995\.
  * [60] F. Sarrazin, F. Pianosi, and T. Wagener, “Global sensitivity analysis of environmental models: Convergence and validation,” Environmental Modelling & Software, vol. 79, pp. 135–152, 2016.
  * [61] E. T. Chang, M. Strong, and R. H. Clayton, “Bayesian sensitivity analysis of a cardiac cell model using a Gaussian process emulator,” PloS one, vol. 10, no. 6, p. e0130252, 2015.
  * [62] C. E. Rasmussen, “Gaussian processes in machine learning,” in Summer school on machine learning, pp. 63–71, Springer, 2003.
  * [63] M. Stein, “Large sample properties of simulations using Latin hypercube sampling,” Technometrics, vol. 29, no. 2, pp. 143–151, 1987.
  * [64] cellML, “cellml: Luo-rudy 1991.” https://models.cellml.org/exposure/456b07d6a7a5b45ed71caad0ea2c0b9d, 2022\.
  * [65] S. Athey and S. Wager, “Estimating treatment effects with causal forests: An application,” Observational Studies, vol. 5, no. 2, pp. 37–51, 2019.
  * [66] C. Glymour, K. Zhang, and P. Spirtes, “Review of causal discovery methods based on graphical models,” Frontiers in genetics, vol. 10, p. 524, 2019\.
  * [67] M. J. Funk, D. Westreich, C. Wiesen, T. Stürmer, M. A. Brookhart, and M. Davidian, “Doubly robust estimation of causal effects,” American journal of epidemiology, vol. 173, no. 7, pp. 761–767, 2011.
  * [68] J. E. Oakley and A. O’Hagan, “Probabilistic sensitivity analysis of complex models: a bayesian approach,” Journal of the Royal Statistical Society: Series B (Statistical Methodology), vol. 66, no. 3, pp. 751–769, 2004.
  * [69] I. Farajpour and S. Atamturktur, “Error and uncertainty analysis of inexact and imprecise computer models,” Journal of Computing in Civil Engineering, vol. 27, no. 4, pp. 407–418, 2013.
  * [70] S. Conti and A. O’Hagan, “Bayesian emulation of complex multi-output and dynamic computer models,” Journal of statistical planning and inference, vol. 140, no. 3, pp. 640–651, 2010.
  * [71] I. Vernon, M. Goldstein, and R. Bower, “Galaxy formation: Bayesian history matching for the observable universe,” Statistical science, pp. 81–90, 2014.
  * [72] U. Kanewala and J. M. Bieman, “Using machine learning techniques to detect metamorphic relations for programs without test oracles,” in 2013 IEEE 24th International Symposium on Software Reliability Engineering (ISSRE), pp. 1–10, IEEE, 2013.
  * [73] J. Ayerdi, V. Terragni, A. Arrieta, P. Tonella, G. Sagardui, and M. Arratibel, “Generating metamorphic relations for cyber-physical systems with genetic programming: an industrial case study,” in Proceedings of the 29th ACM Joint Meeting on European Software Engineering Conference and Symposium on the Foundations of Software Engineering, pp. 1264–1274, 2021.
  * [74] J. Ayerdi, V. Terragni, A. Arrieta, P. Tonella, G. Sagardui, and M. Arratibel, “Evolutionary generation of metamorphic relations for cyber-physical systems,” in Proceedings of the Genetic and Evolutionary Computation Conference Companion, pp. 15–16, 2022.
  * [75] K. Rahman and U. Kanewala, “Predicting metamorphic relations for matrix calculation programs,” in 2018 IEEE/ACM 3rd International Workshop on Metamorphic Testing (MET), pp. 10–13, IEEE, 2018.
  * [76] W. E. Wong, R. Gao, Y. Li, R. Abreu, and F. Wotawa, “A survey on software fault localization,” IEEE Transactions on Software Engineering, vol. 42, no. 8, pp. 707–740, 2016.
  * [77] G. K. Baah, A. Podgurski, and M. J. Harrold, “Causal inference for statistical fault localization,” in Proceedings of the 19th international symposium on Software testing and analysis, pp. 73–84, 2010.
  * [78] G. K. Baah, A. Podgurski, and M. J. Harrold, “Mitigating the confounding effects of program dependences for effective fault localization,” in Proceedings of the 19th ACM SIGSOFT symposium and the 13th European conference on Foundations of software engineering, pp. 146–156, 2011.
  * [79] A. Podgurski and Y. Küçük, “Counterfault: Value-based fault localization by modeling and predicting counterfactual outcomes,” in 2020 IEEE International Conference on Software Maintenance and Evolution (ICSME), pp. 382–393, IEEE, 2020.
  * [80] G. Shu, B. Sun, A. Podgurski, and F. Cao, “Mfl: Method-level fault localization with causal inference,” in 2013 IEEE Sixth International Conference on Software Testing, Verification and Validation, pp. 124–133, IEEE, 2013.
  * [81] Z. Bai, G. Shu, and A. Podgurski, “Numfl: Localizing faults in numerical software using a value-based causal model,” in 2015 IEEE 8th International Conference on Software Testing, Verification and Validation (ICST), pp. 1–10, IEEE, 2015.
  * [82] R. Gore and P. F. Reynolds, “Reducing confounding bias in predicate-level statistical debugging metrics,” in 34th International Conference on Software Engineering (ICSE), pp. 463–473, IEEE, 2012.
  * [83] S. Lee, D. Binkley, R. Feldt, N. Gold, and S. Yoo, “Causal program dependence analysis,” arXiv preprint arXiv:2104.09107, 2021.