# Selective Concept Models:
Permitting Stakeholder Customisation at Test-Time

Matthew Barker,1 Katherine M. Collins, 1 Krishnamurthy (Dj) Dvijotham, 2
Adrian Weller, 1 3 Umang Bhatt 1 3

###### Abstract

Concept-based models perform prediction using a set of concepts that are
interpretable to stakeholders. However, such models often involve a fixed,
large number of concepts, which may place a substantial cognitive load on
stakeholders. We propose Selective COncept Models (SCOMs) which make
predictions using only a subset of concepts and can be customised by
stakeholders at test-time according to their preferences. We show that SCOMs
only require a fraction of the total concepts to achieve optimal accuracy on
multiple real-world datasets. Further, we collect and release a new dataset,
CUB-Sel, consisting of human concept set selections for 900 bird images from
the popular CUB dataset. Using CUB-Sel, we show that humans have unique
individual preferences for the choice of concepts they prefer to reason about,
and struggle to identify the most theoretically informative concepts. The
customisation and concept selection provided by SCOM improves the efficiency
of interpretation and intervention for stakeholders.

## Introduction

Humans can reason about a limited number of concepts at once when making
decisions (Miller 1956; Luck and Vogel 1997; Cowan 2001). While concept-based
methods such as Concept Bottleneck Models (CBMs) (Koh et al. 2020) have been
proposed to support human interpretability and intervenability in machine
learning (ML) systems, such models typically involve dozens of concepts, well
beyond the number of concepts stakeholders can process at any given time
(Tenenbaum 1998; Ramaswamy et al. 2022).

To reduce the cognitive load of reasoning about many concepts, we propose
Selective COncept Models (SCOMs). SCOMs provide a streamlined extension of
CBMs by selecting the concepts that are most pertinent to any given task from
a larger set of available concepts. This enables a stakeholder to reason with
a reduced set of concepts without compromising task accuracy. Unlike CBMs
which require a fixed concept set, SCOMs make predictions using an arbitrary
concept subset which can be customised at inference-time  without retraining.
For example, one might wish to prohibit consideration of sensitive attributes
such as biological sex and subjective attractiveness during prediction. For
SCOMs, withdrawing certain concepts is trivial, whereas conventional CBMs make
such exclusion difficult.

SCOMs enable flexible customisation according to a stakeholder’s preferences
for the number of concepts to use and their personal trade-off between
cognitive load and predictive accuracy (Ramaswamy et al. 2022). On the task of
bird species recognition, SCOMs require only 6 out of 28 concepts to achieve
optimal prediction accuracy. Smaller concept sets decrease the human cost of
interventions, and increase the impact of each intervention. Thus this work is
complementary to research aimed at designing better intervention policies over
a given set of concepts, for example CooP (Chauhan et al. 2022). Since SCOMs
place no restrictions on the exact output network architecture, they provide a
simple extension to existing models used by practitioners.

Figure 1: In CBMs, human-interpretable concepts are predicted from the input
$X$ and used to infer the output $Y$. SCOMs extend CBMs by selecting the most
relevant concepts which maximise the mutual information between the selected
concepts ($C$) and the output ($Y$). The output model learns to make
predictions using an augmented concept vector, which contains the mask used to
select concepts. This allows the number of concepts ($k$) and the specific
concepts selected to be customised at inference time by the stakeholder
without retraining the model.

## Contributions

  * •

We show how feature selection can be applied to the problem of concept
selection, using mutual information (MI) maximisation to provide “useful”
selections, even in the presence of duplicated concepts.

  * •

We propose SCOMs; models which make predictions using an arbitrary subset of
concepts, allowing customisation by stakeholders without retraining.

  * •

We show that SCOMs provide comparable accuracy on real-world datasets CUB and
CelebA, while using only a fraction of the available concepts. We confirm that
SCOM selection algorithms outperform existing techniques for feature
selection.

  * •

We design a human subject experiment to validate the use of SCOMs, and gather
900 sets of concept selections, which we release as CUB-Sel. We identify that
people have varied preferences for the number of concepts they select, and
favour concepts which are less theoretically informative for the task. This
motivates the use of SCOMs to adaptively select lightweight concept sets.

## Definitions

We follow (Koh et al. 2020) in considering concepts as high-level, human-
interpretable – and human-specifiable – values. We define $\mathcal{Z}$ as the
complete set of available concepts, $\mathcal{C}\subset\mathcal{Z}$ as the
chosen set of concepts of size $k$, and $\mathcal{Y}$ as the target
predictions. SCOMs aim to identify the set $\mathcal{C}$ for a given $k$ such
that the performance for a mapping $\mathcal{C}\rightarrow\mathcal{Y}$ is
comparable to the mapping $\mathcal{Z}\rightarrow\mathcal{Y}$. We largely
focus on the setting where $\mathcal{Z}$ is obtained from a Sequential CBM
(Koh et al. 2020), which predicts the concepts based on inputs $\mathcal{X}$
(such as image pixels). However, SCOMs be applied to the inputs $\mathcal{X}$
directly in the case where these inputs are also human concepts. It is assumed
that $\mathcal{Z}$ is given for the SCOM. We use capitalised letters (e.g.
$X$) to denote random variables and lowercase letters to denote instances of a
random variable (e.g. $x$).

## Comparison to Prior Work

#### Concept-Based Models

There has been a proliferation of concept-based models which use concepts as
intermediaries for output predictions. There are several flavours of these
models: those that rely on datasets annotated with concept attributes (Koh et
al. 2020; Marconato, Passerini, and Teso 2022; Zarlenga et al. 2022), those
that learn concepts (Alvarez Melis and Jaakkola 2018), and those that
introduce concept features post-training (Yuksekgonul, Wang, and Zou 2022).
Most models emphasise the motivation of using human-understandable concepts,
yet unlike SCOMs, typically ignore the overall size of the concept set.

#### Model Intervenability

(Chauhan et al. 2022; Shin et al. 2022) give ways of making interventions in
CBMs more effective than random: their work can be used _in conjunction_ with
SCOM. SCOM reduces the size of the concept set used for prediction, increasing
the impact of each intervention, while still allowing intervention procedures
to select the best concept to intervene on. While in practice, stakeholders
may prefer to provide soft labels if they are unsure, (Collins et al. 2023)
show that existing models perform poorly with uncertain interventions.

#### Model Interpretability

There are a wide range of techniques which attempt to improve the
interpretability of concept-based models (Kim et al. 2018; Zhou et al. 2018;
Ghorbani et al. 2019; Yeh et al. 2020; Crabbé and van der Schaar 2022; Abid,
Yuksekgonul, and Zou 2022; Bai et al. 2022). These techniques often provide
metrics highlighting the importance of relevant concepts to the output
prediction. However, (Kim et al. 2022) show with human feedback experiments
that concept-based models using numbers for interpretability can be
overwhelming and lose relevance to the actual image. SCOMs are compatible with
these existing interpretability methods, while reducing the number of concepts
a human has to consider.

#### Feature Engineering and Selection

Feature engineering, the process of applying a set of hand-crafted transforms
to raw input data, is commonly used to increase model accuracy and robustness
of ML models (Zheng and Casari 2018), with significant success across domains
including text modelling (Lewis 1992) and computer vision (Nixon and Aguado
2019). By considering concepts as human-interpretable features, concept
selection in SCOMs is an application of feature engineering. Information-
theoretic methods, commonly used to maximise the statistical dependency
between features and the output, are well-suited to the problem of selecting
concept sets in SCOMs. (Fleuret 2004) have empirical success in approximating
the best global feature set using conditional MI.

Instead of selecting one feature set for the whole dataset, instance-level
techniques can be used to identify custom feature sets for each datapoint. For
example, Learning to Explain (L2X) (Chen et al. 2018) trains a neural network
to predict a variational approximation to MI, and uses a sampler to select the
feature set of a given size which maximises this quantity. Unlike SCOMs,
existing techniques either provide no control over the feature set size (e.g.
INVASE (Yoon, Jordon, and van der Schaar 2019)), or require the size to be
fixed during training (e.g. L2X). To the best of our knowledge, SCOMs are the
first models which allow customisation for instance-level selection _without
retraining_.

## Approach

In this section, we provide a high-level overview for each of the four key
stages in SCOM:

  1. 1.

Concept model (training). Train a model $x\rightarrow c$ which predicts the
concepts ($c$) present in each input ($x$):

$c=g(x)$ (1)

For some datasets, the concepts may be the input $x$ itself.

  2. 2.

Output model (training). Train a neural network which makes predictions using
any arbitrarily chosen concept set $\mathcal{C}\in\mathcal{Z}$. The original
concept vector is augmented to give $\tilde{c}$ according to the chosen
subset, which is then used for prediction:

$\hat{y}=f(\tilde{c})$ (2)

  3. 3.

Concept selection (test). Select the best concept set
$\mathcal{C}\subset\mathcal{Z}$ for each size $k$, using a greedy algorithm
which maximises MI between the concepts and the outputs:

$C_{s}=\operatorname*{arg\,max}_{\mathcal{C}\subset\mathcal{Z}}I(Y;C)\quad\text{subject
to}\quad|C|=k$ (3)

The MI is approximated using the trained network in Step 2.

  4. 4.

Inference (test): A value of $k$ is chosen and inference is performed using
the selected concept set of size $k$:

$y^{*}=f(\tilde{c}_{s})$ (4)

Stakeholders may customise the concept set depending on their individual
preferences.

SCOM addresses two main issues: _how to select a concept set for predictions_
, and _how to perform inference without retraining the model._ These two
problems are interlinked; MI estimation in Step 3 necessitates inference with
a concept set to avoid retraining, while test-time prediction requires concept
selection to find an efficient concept set for a given $k$. Since concept
selection is at the instance level, it occurs at test time. Further details on
each stage are provided next.

### Training the Concept Model

The concept model predicts the concepts present in each input. SCOMs impose no
restrictions on the details of this model, and only require concept
predictions for training further stages and inference. As a result, the
concept model may take on a variety of forms:

  * •

Sequential Bottleneck. Using human-annotated datapoints (input $x$, concepts
$c$), a sequential bottleneck (Koh et al. 2020) trains a model $x\rightarrow
c$ which predicts the concepts present in each input using a neural network.
This is a direct application of multi-label classification in the supervised
learning setting, which is well-studied in the literature (Zhang and Zhou
2013).

  * •

Data processing. In cases when it is not possible to obtain annotated concepts
for each datapoint, data-processing and augmentation techniques may be used to
predict concepts. This is common in the computer vision domain (Perez and Wang
2017), although it is widely used across other domains as well (Kotsiantis,
Kanellopoulos, and Pintelas 2006).

  * •

Direct Inputs. For certain datasets, such as tabular data, the inputs may
already be human-interpretable and thus considered as concepts directly. In
these cases, the concept model for SCOM will simply be the identity function,
i.e. $\hat{c}=x$.

The further stages below assume the concepts have already been obtained from
each datapoint.

### Training the Output Model

The output model predicts the output label using _any concept set_
$\mathcal{C}\subset\mathcal{Z}$. In comparison to the concept model, there is
a lack of prior work on models which can perform inference using arbitrary
predictions. Yoon, Jordon, and van der Schaar attempt to address this
challenge by “masking out” any unused inputs by replacing their values with
zeros. However, this approach prevents the network from distinguishing between
“masked zeros” and input zeros.

#### Concept Augmentation

Instead, SCOMs augment the concept vector with the binary mask vector used to
mask unused concepts, shown in Figure 1. Thus the network has enough
information to identify which zeros are masked, and which are genuine inputs.
The mask is applied to the concept vector $c$ to give the augmented vector:

$\tilde{c}=\begin{bmatrix}c\odot m\\\ m\end{bmatrix}$ (5)

Where $\odot$ denotes the elementwise Hadamard product. For example, consider
the case with 3 concepts, and concepts 1 and 3 are selected. In this case, the
concept vector, mask, and augmented concept vector are given below:

$c=\begin{bmatrix}{\color[rgb]{1,0,0}c_{1}}\\\ {\color[rgb]{1,0,0}c_{2}}\\\
{\color[rgb]{1,0,0}c_{3}}\end{bmatrix},\qquad
m=\begin{bmatrix}{\color[rgb]{0,0,1}1}\\\ {\color[rgb]{0,0,1}0}\\\
{\color[rgb]{0,0,1}1}\end{bmatrix},\qquad\tilde{c}=\begin{bmatrix}{\color[rgb]{1,0,0}c_{1}}\\\
{\color[rgb]{1,0,0}0}\\\ {\color[rgb]{1,0,0}c_{3}}\\\ {\color[rgb]{0,0,1}1}\\\
{\color[rgb]{0,0,1}0}\\\ {\color[rgb]{0,0,1}1}\end{bmatrix}$ (6)

For multi-dimensional concepts, each value in the concept is multiplied by the
corresponding mask value. Thus, SCOMs can select concepts of any dimension,
unlike many existing methods. Once the model has constructed $\tilde{c}$ from
the chosen subset, the output model makes predictions using $\tilde{c}$ as
input:

$\hat{y}=f(\tilde{c})$ (7)

The function $f$ is complex and can be approximated using a neural network,
the training of which is described below.

#### Training the Model

The output model learns a mapping from all possible subsets of $\mathcal{Z}$
to the output $\mathcal{Y}$, i.e.
$\mathcal{P}(\mathcal{Z})\rightarrow\mathcal{Y}$ where $\mathcal{P}$ denotes
the power set. To achieve this, a new mask $m$ is sampled for each batch,
which then augments the concept vector. The predictions are used to calculate
the chosen loss (e.g. cross-entropy), and the model parameters $\theta$ are
updated using Stochastic Gradient Descent (SGD). The training procedure is
detailed in Algorithm 1 below.

Algorithm 1 Training output model

1:for Batch $B$ in $[B_{1},...,B_{N}]$ do

2: $m\sim p(m)$ $\triangleright$ Sample a mask $m$

3: $c\leftarrow c\odot m$ $\triangleright$ Hadamard Product $\odot$

4: $\tilde{c}\leftarrow[c,m]$ $\triangleright$ Append mask to input vector

5: $\hat{y}=f(\tilde{c})$ $\triangleright$ Predictive function $f$

6: $\theta\leftarrow SGD(\hat{y},y,\theta)$ $\triangleright$ Update $\theta$
using SGD

7:end for

It is illustrative to compare the SCOM output model which learns
$\mathcal{P}(\mathcal{Z})\rightarrow\mathcal{Y}$ to a CBM which learns
$\mathcal{Z}\rightarrow\mathcal{Y}$. Since SCOM requires a more complex
mapping than CBM it requires greater model complexity in the form of
larger/additional hidden layers. In addition, the concept augmentation doubles
the size of the input layer in comparison to a CBM.

The probability distribution $p$ used to sample each mask needs to be
determined when training the model. To minimise the expected loss, the
distribution of the mask used during training should match the true mask
distribution. Consider rewriting $p(m)$ using the dependence on $k$:

$\displaystyle p(m)$ $\displaystyle=\sum_{k=1}^{n}p(m,k)$ (8)
$\displaystyle=\sum_{k=1}^{n}p(m|k)p(k)$ (9)

During training, $p(m|k)$ and $p(k)$ are unknown and for SCOMs we assume both
distributions are uniform. The two step sampling procedure is as follows:

  1. 1.

Sample $k\sim U(1,k)$

  2. 2.

Randomly set $k$ values of mask $m$ to 1, with equal probability.

Importantly, this ensures the output model learns to make predictions using
all possible concept set sizes equally. In cases where the distribution of $k$
is known prior to training, using this distribution when sampling $k$ may
achieve faster convergence and a lower expected loss. However, throughout this
work $k$ is assumed to be uniformly distributed.

### Concept Selection

We define _concept selection_ as the method of choosing the concept set in a
way that balances the number of concepts used for a given task and downstream
task performance. (Ghorbani et al. 2019) propose that chosen concepts should
satisfy the three desiderata of meaningfulness, coherency and importance. In
the case of SCOMs, concepts are provided as dataset annotations which are
chosen by humans to be meaningful and coherent. However, the importance of a
concept determines how necessary it is to make accurate predictions and is
task dependent.

We use MI between the concepts and the labels, $I(Y;C)$, as a measure of the
collective prediction accuracy of the concept set. Using an MI measure takes
into account the correlations between concepts, rather than ranking every
concept in isolation. For example, while concepts such as “belly colour” and
“wing colour” might both be important individually for predicting the bird
species, they are likely to be strongly correlated. Therefore, only one
concept may be required to make accurate predictions. Precisely, the task of
selecting $k$ concepts can be framed as one of maximising MI subject to a
cardinality constraint (Fleuret 2004):

$\operatorname*{arg\,max}_{\mathcal{C}\subset\mathcal{Z}}I(Y;C)\quad\text{subject
to}\quad|C|=k$ (10)

The naive approach of trying every subset of $\mathcal{Z}$ with size $k$
suffers from combinatorial explosion and is unfeasible. However, we can
achieve a _good_ solution in polynomial time using greedy algorithms
(Nemhauser, Wolsey, and Fisher 1978). There are two common greedy algorithms:
forward selection (FS) and backward elimination (BE). FS starts with an empty
set $\mathcal{C}=\emptyset$, and then sequentially adds individual concepts
$c\in\mathcal{Z}$, maximizing $I(Y;C\cup c_{i})$ at every stage. In
comparison, BE starts with the complete set $\mathcal{C}=\mathcal{Z}$, and
removes individual concepts $c\in\mathcal{C}$, maximizing $I(Y;C\setminus
c_{i})$ at every stage. The two procedures are detailed in Algorithms 2 and 3.

Algorithm 2 Forward selection procedure

1:$Z\leftarrow\\{c_{1},c_{2},...\\}$ $\triangleright$ Complete set of
available concepts

2:$C\leftarrow\emptyset$ $\triangleright$ Start with empty set

3:while $|C|\leq k$ do

4:
$c_{\text{next}}\leftarrow\operatorname*{arg\,max}_{c_{i}\in\mathcal{Z}}I(Y;C\cup
c_{i})$ $\triangleright$ Maximise MI

5: $C\leftarrow C\cup c_{\text{next}}$

6:end while

Algorithm 3 Backward elimination procedure

1:$Z\leftarrow\\{c_{1},c_{2},...\\}$ $\triangleright$ Complete set of
available concepts

2:$C\leftarrow Z$ $\triangleright$ Start with complete set

3:while $|C|\geq k$ do

4:
$c_{\text{next}}\leftarrow\operatorname*{arg\,max}_{c_{i}\in\mathcal{C}}I(Y;C\setminus
c_{i})$ $\triangleright$ Maximise MI

5: $C\leftarrow C\setminus c_{\text{next}}$

6:end while

Due to the incremental nature of greedy selection, every concept set of size
$\leq k$ (for FS) and $\geq k$ (for BE) is also obtained when finding the
concept size of size $k$. This allows concept sets of all sizes to be
calculated by running the selection procedure once for $k=|Z|$ (FS) and $k=0$
(BE). Since SCOMs place _no constraint on the dimensionality of concepts_ ,
calculating $I(Y;C)$ requires an approximation when concepts are continuous.
Existing methods for estimating MI between high-dimensional continuous
variables such as the kernel density estimator (Kolchinsky and Tracey 2017) or
k-nearest neighbours method (Kraskov, Stögbauer, and Grassberger 2004) are
sensitive to exact parameters and are not well suited to distributions of
concept logits, which tend to be bimodal for binary concepts. To avoid this
issue, rather than maximising MI directly, in Theorem 1 we show that
minimising the output entropy $H(\hat{Y}|C)$ instead gives the same optimal
concept set. Unlike MI, output entropy is simple to calculate from the model
predictions:

$H(\hat{Y}|C)=-\sum_{y\in\mathcal{Y}}p(\hat{y}|c)\log p(\hat{y}|c)$ (11)

Using entropy as a proxy for MI, the maximisation at each stage in FS
(Algorithm 2) can be solved using:

$c_{\text{next}}=\operatorname*{arg\,max}_{c_{i}\in\mathcal{Z}}I(Y;C\cup
c_{i})=\operatorname*{arg\,min}_{c_{i}\in\mathcal{Z}}H(Y|C\cup c_{i})$ (12)

Similarly, for BE (Algorithm 3) the optimisation at each stage is given by:

$c_{\text{next}}=\operatorname*{arg\,max}_{c_{i}\in\mathcal{C}}I(Y;C\setminus
c_{i})=\operatorname*{arg\,min}_{c_{i}\in\mathcal{C}}H(Y|C\setminus c_{i})$
(13)

### Inference

A core property of SCOMs is that predictions can be made from arbitrarily
chosen concept sets – without retraining the model. Once the output model has
been trained (Step 2) and the concept set has been selected (Step 3),
inference is simple and follows a three stage procedure:

  1. 1.

The concepts $c$ are predicted from the test datapoint $x^{*}$:

$c=g(x^{*})$ (14)

  2. 2.

A subset of concepts, $c_{s}$, is chosen according to the selection algorithm,
and optionally customised by a stakeholder.

  3. 3.

The concept subset is augmented to give $\tilde{c}_{s}$ and used to infer the
output:

$y^{*}=f(\tilde{c_{s}})$ (15)

## Computational Experiments

To validate the SCOM selection algorithms compared to existing methods, we
assess the performance of SCOMs on datasets for concept-based models popular
in the literature: the CUB dataset for bird species recognition (Welinder et
al. 2010) and the CelebFaces Attributes Dataset (CelebA) (Liu et al. 2015) for
facial identity recognition. For both datasets, a multi-layer perceptron with
two hidden layers, each containing 100 neurons, was trained for the output
model. Using all the concepts, a prediction accuracy of 75.3% and 63.8% was
achieved for CUB and CelebA respectively. Further details for each dataset are
provided below.

#### CUB

We follow (Koh et al. 2020) in using only 112 of the 312 original CUB binary
attributes (e.g. “has blue wings”, “has a pointed beak”). These attributes are
further categorised into the 28 multi-valued concepts (e.g. “wing colour”)
used in the original crowdsourced annotation by (Welinder et al. 2010). For
the purpose of intervention and interpretation, one multi-valued concept more
accurately represents a human “concept” than the binary attributes. Matching
(Koh et al. 2020), the concept labels are fixed at the class-level, so that
each bird from the same species is labelled with the same concept annotations.
This property, along with large number of concept groups, makes CUB
particularly well-suited for SCOMs. Further, the same predicted concept logits
for the Sequential bottleneck were used when training and evaluating SCOM
performance, which achieve a concept accuracy of 96.6%.

Since the CUB concepts are multi-valued, methods which select individual
values (e.g. L2X) aren’t suitable. The prediction accuracy of randomly
selected concepts is used as a baseline. Errors are reported over 3 random
seeds which determine the concept logit predictions, with the train/test split
given by (Wah et al. 2011).

#### CelebA

CelebA (Liu et al. 2015) is a dataset with 40 concept annotations for 200K
celebrity faces from 10K people, which has been used to evaluate concept
models (Ghorbani et al. 2019; Zarlenga et al. 2022). For computational
efficiency, we subsample CelebA by using random subsets containing 20
identities for the prediction task. Unlike the processed version of CUB,
attributes for CelebA are determined at the image level, meaning different
images of the same person may have different concept annotations. Although the
annotations are noisy, attributes such as “wearing a hat” cannot be sensibly
fixed at the class level. Thus CelebA is used to investigate the suitability
of SCOMs for instance-level concept annotations. A ResNet-50 (He et al. 2016)
was trained on all 10k identities to predict the attributes for each face,
with an accuracy of $0.916$.

Unlike CUB, the concepts are single-valued attributes, and thus L2X was used
as a benchmark. Since L2X cannot make predictions directly, concepts are
selected by L2X and then used as inputs to the trained SCOM output model.
Errors are reported over 5 random seeds which determine the samples taken from
the dataset.

### Results

#### Prediction Accuracy

Figure 2: Prediction accuracy of SCOMs on the CUB and CelebA datasets as a
function of $k$, compared to L2X and random baselines. On CUB (left), BE SCOM
performs significantly better than FS SCOM at the instance level. On CelebA
(right), BE SCOM performs much better than FS SCOM, and better than L2X for
low numbers of concepts. For all plots, errors are $\pm 1\sigma$, calculated
over random seeds which affects the concepts randomly selected, the concept
model chosen for CUB, and the dataset sample for CelebA.

Figure 3: Prediction accuracy for CUB and CelebA vs number of oracle
interventions, for certain values of $k$. Smaller concept sets have a greater
increase in accuracy after the same number of interventions. Errors are $\pm
1\sigma$, calculated over 10 random seeds which determines the concepts that
are intervened. Figure 4: Example interventions using selected concept sets
from CelebA (left, 8 concepts) and CUB (right, 6 concepts). The concept logit
predictions are rounded to give binary predictions, and only the positive
concepts are shown. Each intervention changes one concept which corrects the
original incorrect model prediction. When using selected concept sets, it is
much easier for experts to identify incorrect concepts and reduces the human
cost of interventions.

We first investigate how task accuracy compares between SCOMs, which permit
flexible test-time intervention without retraining against the L2X baseline,
which retrains a new model for each value of $k$. The results are shown in
Figure 2 and numerical values are provided in Tables 1 and 2.

For CUB and CelebA, the BE methods perform better than the greedier FS
methods. On CUB, BE SCOMs perform particularly well, achieving the optimal
accuracy of 76.1% when $k=6$. In addition, BE SCOM outperforms the L2X
baseline on CelebA for low numbers of concepts, with approximately equal
performance for large concept sets.

#### Intervention Accuracy

One of the key benefits of CBMs is that they allow interventions from a human
expert at the concept level to improve the accuracy of the model. To test this
property for SCOMs, concepts were replaced by predictions from an oracle,
chosen randomly for simplicity. We assume that human experts are able to
identify the true value of the concepts with perfect accuracy. The task of
intervening on concepts is separate from the concept selection procedure
performed in the human experiment.

The initial concept sets were selected using the BE variant of SCOMs, since
that gives the greatest prediction accuracy for both CUB and CelebA. For CUB,
the oracle provides binary values for each concept, which are determined by
the species that the bird belongs to. For CelebA, there is no clear class-
level oracle, and individual image annotations are too noisy to be considered
an oracle. Following work on human-derived soft annotations (Peterson et al.
2019; Collins, Bhatt, and Weller 2022), we assume that the CelebA oracle
intervenes on each concept with the “correct uncertainty”, estimated as the
mean value of that concept over all images of the same person. While we
acknowledge that this oracle is unrealistic, we aim to show the efficacy of
interventions assuming they can be obtained. Though beyond the scope of this
work, we refer to Collins et al. for an exposition of why such an oracle-
intervention assumption may not always be wise in practice.

Figure 3 highlights that the intervention accuracy increases more with each
intervention for smaller concept set sizes ($k$). After only a few
interventions, the accuracy gain for small $k$ overcomes the initial accuracy
penalty compared to the full concept set. As each concept in a small concept
set has a greater impact on the overall prediction, it is intuitive that each
intervention yields a larger effect. Examples of interventions are shown in
Figure 4, where a single oracle intervention corrects the initial prediction
made by the model. Since oracle interventions are more effective for smaller
concept sets, SCOMs exhibit greater intervenability compared to existing
concept models.

### Takeaways

On the real-world datasets CUB and CelebA, SCOMs perform better than existing
methods. BE methods are much more accurate than FS, and require a fraction of
the number of concepts for optimal performance.

## Human Experiment

Figure 5: Interface for the human experiment. Participants were shown an image
of a bird, and then asked to select the relevant concepts from the list of 28
total concepts. Figure 6: Left: Distribution of number of concepts selected
for each image (left) and the mean number of concepts selected by each
participant (center). Both distributions have a large variance, showing the
number of concepts humans prefer to reason varies. The standard deviation of
the number of concepts selected by each participant (right) is comparatively
large, motivating the instance-level selection SCOM provides.  Figure 7:
Example human selections of the most relevant concepts. Here, we depict cases
where three participants labelled the same image. Concepts with at least one
disagreement are shown in blue, and illustrate that people have varied
individual preferences. Figure 8: Bar chart comparing the frequency of
concepts selected by participants to the selections by the BE SCOM. Humans
prefer concepts which SCOMs rarely select, emphasising the importance of
customisation for the concept sets selected by SCOM. Errors bars of $\pm
1\sigma$ are provided for the human selections.

Our computational results demonstrate the methodological prowess of SCOMs.
However, we care principally about SCOMs to work with and support real humans.
One of the key motivations for SCOMs is that they allow stakeholder
customisation of concept sets _at inference time_. The need for customisation
rests on the assumption that humans have varied preferences for the number of
concepts they prefer to reason with ($k$). Further, the desire for algorithmic
selection of those concepts assumes that people may not always choose the most
informative concepts on their own. We validate these two assumptions through a
human subject experiment, wherein participants were asked to select the
concept groups for each image. We release the resulting annotations in a new
dataset, CUB-Sel.

#### Experimental Details

We recruited 30 US-based participants from the crowd-sourcing platform
Prolific (Palan and Schitter 2018). Participants were asked to select as many
concepts as they thought were relevant when predicting the bird species for
each image. Selection implicitly captures their preference for the number of
concepts to reason with ($k$). Each participant was shown 30 images, yielding
a total of 900 concept set selections (which compose CUB-Sel). The images were
randomly sampled from the CUB test dataset, with a total of 28 possible
concepts for each image. Participants were paid at a base rate of $8/hour,
with a bonus rate of up to $9/hour to encourage high-quality selections. The
study interface is shown in Figure 5. Further details on CUB are provided in
the Computational Experiments.

### Results

#### Preferred Concept Set Sizes

To motivate the customisation SCOMs provide, we first investigate the
distribution over $k$, i.e., the number of concepts people selected. Figure 6
shows the distribution of $k$ values both for each image and for each
participant. A participant’s $k$ was calculated as the mean $k$ selected over
the 30 images they were shown. Both distributions have a large variance,
showing that the number of concepts humans prefer to reason about varies
significantly – across people and individual images. There was a mean
difference of 9.6 concepts when several participants labelled the same image,
emphasising that different people prefer to reason with different concepts.
Examples of multiple concept selections for the same image are shown in Figure
7. Thus, the customisation SCOM provides is essential to cater for individual
stakeholder preferences.

#### Comparison Between Humans and SCOMs

SCOMs assume that humans struggle to select the relevant concepts from a
larger set. To evaluate the quality of human selections, we computed the
accuracy of the SCOM output model on CUB using human selections and define
this as the “human accuracy”. Importantly, although the concept sets are
chosen by humans, the SCOM model is still used for inference. Thus, the “human
accuracy” is a measure of how informative concepts are _for the model_ ; the
information concepts provide may differ if people were responsible for
classification instead.

We compare the human accuracy to SCOM and L2X instance-level selection methods
in Figure 9. Remarkably, the human accuracy is lower than random for all $k$
values and far lower than the BE SCOM algorithms. To explore human versus
model selection further, we compare the frequency of concepts selected by
humans to the selections by a SCOM 111We use BE as the SCOM baseline,
evaluated on the same images shown to participants (which represent a subset
of the test set). in Figure 8. There is a significant difference between the
distribution of concepts selected by humans vs SCOM. Concepts such as “Shape”
and “Bill length” are frequently selected by humans, but not favoured by SCOM.
Concepts which are not specific to bird species recognition are easier to
understand and thus treated preferentially by humans, even if they are not the
most informative concepts. Challenging concepts such as “Nape color” and
“Under tail color” are rarely selected by humans, and some participants noted
that they struggled to understand these concepts.

### Takeaways

The human subject experiment shows that humans have individual $k$ values
which vary significantly, motivating the test time customisation SCOMs
provide. Further, humans favour concepts which are less theoretically
informative, as measured by SCOM selections. Consequently, the accuracy of
model predictions when using human selection is significantly lower than SCOM
selection.

## Discussion

We investigated the performance of SCOMs compared to random and L2X baselines,
for two real world datasets. The selection algorithm uses MI to select the
best concept set of a given size at either the dataset or instance level, and
is robust to correlations between concepts. For CUB, prediction accuracy is
far superior to random and human selections. For CelebA, SCOMs outperform L2X,
especially for low numbers of concepts. After a few interventions, predictions
using the selected subset outperforms those using all concepts. Through our
human subject experiment, we show that humans have different individual
preferences for the number and type of concepts they prefer to reason about,
motivating the value of SCOMs. By avoiding retraining on different concept
sets, SCOMs allow stakeholders flexibility to customise the concepts used to
make predictions after the model has been trained. Next, we discuss
implications of the results and potential limitations of SCOMs.

### Choosing the Concept Set Size $k$ in Practice

When SCOMs are used by real stakeholders, there are several methods of
choosing the concept set size $k$. A practitioner could employ a validation
set and choose $k$ to either maximise the accuracy, or as the minimum value
which surpasses a certain accuracy threshold. Alternatively, $k$ may be chosen
according to the preferences of the experts who are intervening on the
concepts.

### Human vs Algorithm Concept Preferences

Our human experiment elucidated that humans and SCOM selection algorithms
disagree on which concepts should be selected; humans prefer generic concepts,
while SCOMs prefer specific concepts which provide more information. Although
it decreases the human cost, allowing humans to customise the concept sets
originally selected algorithmically is likely to _decrease_ the prediction
accuracy as well. Hence, there is a compromise between accuracy and cost of
intervention/interpretation which ought to be considered when using SCOMs in
practice. Note, we acknowledge that the participants who contributed to CUB-
Sel were not necessarily experts in the task. Experts may have different
preferences when selecting concept sets, potentially changing the overall
prediction accuracy as well.

### Limitations

While our approach holds promise for improving the usability of concept models
for stakeholders, it is worth noting conditions which permit maximal
effectiveness: SCOMs work best on datasets with a large number of concepts
(e.g., CUB) and rely on a well-trained output model. In cases where all the
concepts are required for good prediction, smaller concept sets selected by
SCOM will suffer an accuracy penalty. Our selection algorithms aim to minimise
output entropy which is only a good measure of MI when the output model is
_well-trained_. Since the output model learns a mapping from every concept
subset to the output, it may converge slower during training than models which
only consider one concept set. Moreover, so far our work has only studied
SCOMs in the classification setting; it is worth studying SCOMs in other
paradigms (e.g., regression) and over additional real-world concept datasets.

When removing concepts at test-time, we recognise that while this removes
explicit dependence from the model, it does not completely cleanse implicit
dependence on the removed concepts. Sensitive concepts may be highly
correlated with other concepts which are not removed, meaning some bias may
still be present in the model. Further, the output model was trained on all
the concepts, which each affect the parameters of the model even if they are
later removed from selection. For concept sets initially selected, it is
assumed that all concepts are available for selection. Removing a concept at
test time invalidates this assumption and reduces the efficacy of the concept
selection procedure.

### Future Work

One of the benefits of our prediction method, which uses arbitrary concept
subsets, is that it allows stakeholder customisation at inference-time without
retraining. The human experiment in this work was focused on evaluating human
concept selection. In reality, the concept sets are likely to be selected by
SCOM, and then provided to human experts for interpretation and intervention.
An exciting direction following this work is to validate these usability
benefits of SCOMs through further human experiments. By reducing the cost of
expert interventions, SCOMs facilitate end-to-end experiments involving
concept selection, stakeholder customisation and expert intervention. Such
experiments would test the use of concept models in real-world settings,
paving the way for the deployment of SCOMs in industry.

We encourage further investigation into discrepancies between human and
algorithm concept preferences – and potential for complementarity (Steyvers et
al. 2022) – through CUB-Sel and related datasets such as CUB-S (Collins et al.
2023). In addition, smarter intervention methods, such as CooP (Chauhan et al.
2022), could be used in conjunction with SCOMs to further increase the
efficacy of interventions.

## Conclusion

In this work, we have introduced SCOMs, an algorithmic framework which makes
predictions using a selected concept set for a given task. SCOMs provide a way
of improving the intervenability and interpretability of concept based models,
with no compromise on accuracy. On the task of bird species classification,
optimal accuracy (76.1%) was achieved using only 6 out of the available 28
concepts. Oracle interventions are more effective for the selected concept
sets, demonstrating greater intervenability for SCOMs compared to existing
concept models.

Using a human subject experiment we have shown that the concept sets humans
choose often differ from the concepts which are most useful to the model.
SCOMs account for this difference by allowing humans to customise the selected
concept sets at inference time _without retraining_ , according to their
individual preferences. SCOMs place no restrictions on the form of concepts
selected or the output model architecture, allowing for easy integration with
other techniques. We release the annotations from our human experiment in a
new dataset, CUB-Sel, to encourage further investigation into the utility of
SCOMs in concert with real humans.

## Acknowledgments

We thank the participants from Prolific who took part in our human experiment.
KMC acknowledges support from the Marshall Commission and the Cambridge Trust.
AW acknowledges support from a Turing AI Fellowship under grant EP/V025279/1,
The Alan Turing Institute, and the Leverhulme Trust via CFI. UB acknowledges
support from DeepMind and the Leverhulme Trust via the Leverhulme Centre for
the Future of Intelligence (CFI), and from the Mozilla Foundation.

## References

  * Abid, Yuksekgonul, and Zou (2022) Abid, A.; Yuksekgonul, M.; and Zou, J. 2022.  Meaningfully debugging model mistakes using conceptual counterfactual explanations.  In _International Conference on Machine Learning_ , 66–88. PMLR.
  * Alvarez Melis and Jaakkola (2018) Alvarez Melis, D.; and Jaakkola, T. 2018.  Towards robust interpretability with self-explaining neural networks.  _Advances in neural information processing systems_ , 31.
  * Bai et al. (2022) Bai, A.; Yeh, C.-K.; Ravikumar, P.; Lin, N. Y.; and Hsieh, C.-J. 2022.  Concept Gradient: Concept-based Interpretation Without Linear Assumption.  _arXiv preprint arXiv:2208.14966_.
  * Chauhan et al. (2022) Chauhan, K.; Tiwari, R.; Freyberg, J.; Shenoy, P.; and Dvijotham, K. 2022.  Interactive Concept Bottleneck Models.  _arXiv preprint arXiv:2212.07430_.
  * Chen et al. (2018) Chen, J.; Song, L.; Wainwright, M.; and Jordan, M. 2018.  Learning to explain: An information-theoretic perspective on model interpretation.  In _International Conference on Machine Learning_ , 883–892. PMLR.
  * Collins et al. (2023) Collins, K. M.; Barker, M.; Zarlenga, M. E.; Raman, N.; Bhatt, U.; Jamnik, M.; Sucholutsky, I.; Weller, A.; and Dvijotham, K. 2023.  Human Uncertainty in Concept-Based AI Systems.  _arXiv preprint arXiv:2303.12872_.
  * Collins, Bhatt, and Weller (2022) Collins, K. M.; Bhatt, U.; and Weller, A. 2022.  Eliciting and Learning with Soft Labels from Every Annotator.  _Proceedings of the AAAI Conference on Human Computation and Crowdsourcing_ , 10(1): 40–52.
  * Cowan (2001) Cowan, N. 2001.  The magical number 4 in short-term memory: A reconsideration of mental storage capacity.  _Behavioral and brain sciences_ , 24(1): 87–114.
  * Crabbé and van der Schaar (2022) Crabbé, J.; and van der Schaar, M. 2022.  Concept Activation Regions: A Generalized Framework For Concept-Based Explanations.  _arXiv preprint arXiv:2209.11222_.
  * Fleuret (2004) Fleuret, F. 2004.  Fast binary feature selection with conditional mutual information.  _Journal of Machine learning research_ , 5(9).
  * Ghorbani et al. (2019) Ghorbani, A.; Wexler, J.; Zou, J. Y.; and Kim, B. 2019.  Towards automatic concept-based explanations.  _Advances in Neural Information Processing Systems_ , 32.
  * He et al. (2016) He, K.; Zhang, X.; Ren, S.; and Sun, J. 2016.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 770–778.
  * Kim et al. (2018) Kim, B.; Wattenberg, M.; Gilmer, J.; Cai, C.; Wexler, J.; Viegas, F.; et al. 2018\.  Interpretability beyond feature attribution: Quantitative testing with concept activation vectors (tcav).  In _International conference on machine learning_ , 2668–2677. PMLR.
  * Kim et al. (2022) Kim, S. S.; Watkins, E. A.; Russakovsky, O.; Fong, R.; and Monroy-Hernández, A. 2022.  ” Help Me Help the AI”: Understanding How Explainability Can Support Human-AI Interaction.  _arXiv preprint arXiv:2210.03735_.
  * Koh et al. (2020) Koh, P. W.; Nguyen, T.; Tang, Y. S.; Mussmann, S.; Pierson, E.; Kim, B.; and Liang, P. 2020.  Concept bottleneck models.  In _International Conference on Machine Learning_ , 5338–5348. PMLR.
  * Kolchinsky and Tracey (2017) Kolchinsky, A.; and Tracey, B. D. 2017.  Estimating mixture entropy with pairwise distances.  _Entropy_ , 19(7): 361.
  * Kotsiantis, Kanellopoulos, and Pintelas (2006) Kotsiantis, S.; Kanellopoulos, D.; and Pintelas, P. 2006.  Data Preprocessing for Supervised Learning.  _International Journal of Computer Science_ , 1: 111–117.
  * Kraskov, Stögbauer, and Grassberger (2004) Kraskov, A.; Stögbauer, H.; and Grassberger, P. 2004.  Estimating mutual information.  _Physical review E_ , 69(6): 066138.
  * Lewis (1992) Lewis, D. D. 1992.  _Representation and learning in information retrieval_.  University of Massachusetts Amherst.
  * Liu et al. (2015) Liu, Z.; Luo, P.; Wang, X.; and Tang, X. 2015.  Deep Learning Face Attributes in the Wild.  In _Proceedings of International Conference on Computer Vision (ICCV)_.
  * Luck and Vogel (1997) Luck, S. J.; and Vogel, E. K. 1997.  The capacity of visual working memory for features and conjunctions.  _Nature_ , 390(6657): 279–281.
  * Marconato, Passerini, and Teso (2022) Marconato, E.; Passerini, A.; and Teso, S. 2022.  GlanceNets: Interpretabile, Leak-proof Concept-based Models.  _arXiv preprint arXiv:2205.15612_.
  * Miller (1956) Miller, G. A. 1956.  The magical number seven, plus or minus two: Some limits on our capacity for processing information.  _Psychological review_ , 63(2): 81.
  * Nemhauser, Wolsey, and Fisher (1978) Nemhauser, G. L.; Wolsey, L. A.; and Fisher, M. L. 1978.  An analysis of approximations for maximizing submodular set functions—I.  _Mathematical programming_ , 14(1): 265–294.
  * Nixon and Aguado (2019) Nixon, M.; and Aguado, A. 2019.  _Feature extraction and image processing for computer vision_.  Academic press.
  * Palan and Schitter (2018) Palan, S.; and Schitter, C. 2018.  Prolific. ac—A subject pool for online experiments.  _Journal of Behavioral and Experimental Finance_ , 17: 22–27.
  * Perez and Wang (2017) Perez, L.; and Wang, J. 2017.  The effectiveness of data augmentation in image classification using deep learning.  _arXiv preprint arXiv:1712.04621_.
  * Peterson et al. (2019) Peterson, J. C.; Battleday, R. M.; Griffiths, T. L.; and Russakovsky, O. 2019.  Human uncertainty makes classification more robust.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 9617–9626.
  * Ramaswamy et al. (2022) Ramaswamy, V. V.; Kim, S. S.; Fong, R.; and Russakovsky, O. 2022.  Overlooked factors in concept-based explanations: Dataset choice, concept salience, and human capability.  _arXiv preprint arXiv:2207.09615_.
  * Shin et al. (2022) Shin, S.; Jo, Y.; Ahn, S.; and Lee, N. 2022.  A Closer Look at the Intervention Procedure of Concept Bottleneck Models.  In _Workshop on Trustworthy and Socially Responsible Machine Learning, NeurIPS 2022_.
  * Steyvers et al. (2022) Steyvers, M.; Tejeda, H.; Kerrigan, G.; and Smyth, P. 2022.  Bayesian modeling of human–AI complementarity.  _Proceedings of the National Academy of Sciences_ , 119(11): e2111547119.
  * Tenenbaum (1998) Tenenbaum, J. 1998.  Bayesian modeling of human concept learning.  _Advances in neural information processing systems_ , 11.
  * Wah et al. (2011) Wah, C.; Branson, S.; Welinder, P.; Perona, P.; and Belongie, S. 2011.  The caltech-ucsd birds-200-2011 dataset.
  * Welinder et al. (2010) Welinder, P.; Branson, S.; Mita, T.; Wah, C.; Schroff, F.; Belongie, S.; and Perona, P. 2010.  Caltech-UCSD birds 200.
  * Yeh et al. (2020) Yeh, C.-K.; Kim, B.; Arik, S.; Li, C.-L.; Pfister, T.; and Ravikumar, P. 2020.  On completeness-aware concept-based explanations in deep neural networks.  _Advances in Neural Information Processing Systems_ , 33: 20554–20565.
  * Yoon, Jordon, and van der Schaar (2019) Yoon, J.; Jordon, J.; and van der Schaar, M. 2019.  INVASE: Instance-wise variable selection using neural networks.  In _International Conference on Learning Representations_.
  * Yuksekgonul, Wang, and Zou (2022) Yuksekgonul, M.; Wang, M.; and Zou, J. 2022.  Post-hoc Concept Bottleneck Models.  _arXiv preprint arXiv:2205.15480_.
  * Zarlenga et al. (2022) Zarlenga, M. E.; Barbiero, P.; Ciravegna, G.; Marra, G.; Giannini, F.; Diligenti, M.; Shams, Z.; Precioso, F.; Melacci, S.; Weller, A.; et al. 2022.  Concept embedding models.  _arXiv preprint arXiv:2209.09056_.
  * Zhang and Zhou (2013) Zhang, M.-L.; and Zhou, Z.-H. 2013.  A review on multi-label learning algorithms.  _IEEE transactions on knowledge and data engineering_ , 26(8): 1819–1837.
  * Zheng and Casari (2018) Zheng, A.; and Casari, A. 2018.  _Feature engineering for machine learning: principles and techniques for data scientists_.  ” O’Reilly Media, Inc.”.
  * Zhou et al. (2018) Zhou, B.; Sun, Y.; Bau, D.; and Torralba, A. 2018.  Interpretable basis decomposition for visual explanation.  In _Proceedings of the European Conference on Computer Vision (ECCV)_ , 119–134.

##  Appendix A Theorems

######  Theorem 1.

For a well-trained estimator, minimising the output entropy conditioned on the
concepts, $H(\hat{Y}|C)$, is equivalent to maximising the MI, $I(Y;C)$:

$\operatorname*{arg\,min}_{\mathcal{C}\subset\mathcal{Z}}H(\hat{Y}|C)=\operatorname*{arg\,max}_{\mathcal{C}\subset\mathcal{Z}}I(Y;C)$
(16)

###### Proof.

Using the independence of $H(Y)$ with respect to $C$:

$\operatorname*{arg\,min}_{\mathcal{C}\subset\mathcal{Z}}H(\hat{Y}|C)=\operatorname*{arg\,min}_{\mathcal{C}\subset\mathcal{Z}}H(\hat{Y}|C)-H(Y)$
(17)

For a well-trained estimator, $p(y|c)=p(\hat{y}|c)$, implying
$H(Y|C)=H(\hat{Y}|C)$. Using this assumption:

$\operatorname*{arg\,min}_{\mathcal{C}\subset\mathcal{Z}}H(\hat{Y}|C)=\operatorname*{arg\,min}_{\mathcal{C}\subset\mathcal{Z}}H(Y|C)-H(Y)$
(18)

Noting that the RHS is the negative MI:

$\operatorname*{arg\,min}_{\mathcal{C}\subset\mathcal{Z}}H(\hat{Y}|C)=\operatorname*{arg\,min}_{\mathcal{C}\subset\mathcal{Z}}-I(Y;C)$
(19)

Reversing the minus sign gives the required result:

$\operatorname*{arg\,min}_{\mathcal{C}\subset\mathcal{Z}}H(\hat{Y}|C)=\operatorname*{arg\,max}_{\mathcal{C}\subset\mathcal{Z}}I(Y;C)$
(20)

∎

##  Appendix B Further Results

Figure 9: Prediction accuracy on the CUB dataset as a function of $k$. Accuracy using human selections is even worse than random, motivating the need for SCOMs to provide algorithmic selections. Errors are $\pm 1\sigma$, calculated over random seeds which affects the concepts randomly selected and the concept model chosen. CUB | Method
---|---
$k$ | BE | FS | Human | Random
1 | $22.4\pm 1$ | $16.1\pm 1$ | $0.0\pm 0$ | $4.4\pm 4$
3 | $71.1\pm 0.7$ | $46.4\pm 1$ | $18.9\pm 6$ | $26.9\pm 7$
6 | $75.5\pm 1$ | $62.6\pm 1$ | $42.7\pm 2$ | $55.5\pm 6$
12 | $75.8\pm 2$ | $65.7\pm 0.5$ | $53.0\pm 0.7$ | $71.3\pm 2$
28 (All) | $75.3\pm 0.6$ | $75.3\pm 0.6$ | - | $75.3\pm 0.6$
Table 1: Prediction accuracy on the CUB dataset for SCOMs compared to the random and human benchmarks for selected values of $k$. Errors are $\pm 1\sigma$, calculated over 3 random seeds which determine the concept logit predictions and the initial parameters of the output model. CelebA | Method
---|---
$k$ | BE | FS | L2X | Random
1 | $35.8\pm 10$ | $30.0\pm 10$ | $15.9\pm 6$ | $8.8\pm 2$
5 | $60.8\pm 8$ | $42.5\pm 8$ | $45.8\pm 6$ | $22.0\pm 5$
10 | $63.3\pm 2$ | $45.8\pm 9$ | $56.2\pm 4$ | $35.5\pm 5$
20 | $65.0\pm 2$ | $48.3\pm 10$ | $62.9\pm 2$ | $51.2\pm 5$
40 (All) | $63.8\pm 3$ | $63.8\pm 3$ | $63.8\pm 3$ | $63.8\pm 3$
Table 2: Prediction accuracy on the CelebA dataset for SCOMs compared to the
random and L2X instance benchmarks for selected values of $k$. Errors are $\pm
1\sigma$, calculated over 5 random seeds which determine the train/test split
and the initial parameters of the output model.


# Measuring Accessibility using Gravity and Radiation Models

Duccio Piovani  Centre for Advanced Spatial Analysis (CASA), University
College London (UCL), 90 Tottenham Court Road , London, W1T 4TJ nam.r , 3 rue
de Steinkerque, 75018 Paris<EMAIL_ADDRESS>Elsa Arcaute  Centre
for Advanced Spatial Analysis (CASA), University College London (UCL), 90
Tottenham Court Road , London, W1T 4TJ Gabriela Uchoa  Centre for Advanced
Spatial Analysis (CASA), University College London (UCL), 90 Tottenham Court
Road , London, W1T 4TJ Prefeitura Municipal de Teresina, Praça Marechal
Deodoro da Fonseca 860, Teresina, Brasil 64000-160
Alan Wilson  Centre for Advanced Spatial Analysis (CASA), University College
London (UCL), 90 Tottenham Court Road , London, W1T 4TJ The Alan Turing
Institute, British Library, 96 Euston Road, London NW1  Michael Batty  Centre
for Advanced Spatial Analysis (CASA), University College London (UCL), 90
Tottenham Court Road , London, W1T 4TJ

###### Abstract

Since the presentation of the Radiation Model, much work has been done to
compare its findings with those obtained from Gravitational Models. These
comparisons always aim at measuring the accuracy with which the models
reproduce the mobility described by origin-destination matrices. This has been
done at different spatial scales using different datasets, and several
versions of the models have been proposed to adjust to various spatial
systems. However the models, to our knowledge, have never been compared with
respect to policy testing scenarios. For this reason, here we use the models
to analyze the impact of the introduction of a new transportation network, a
Bus Rapid Transport system, in the city of Teresina in Brazil. We do this by
measuring the estimated variation in the trip distribution, and formulate an
accessibility to employment indicator for the different zones of the city. By
comparing the results obtained with the two approaches, we are able, not only
to better assess the goodness of fit and the impact of this intervention, but
also to understand reasons for the systematic similarities and differences in
their predictions.

## Introduction

Assessing the impact of new infrastructural projects is a challenging and
demanding task that requires knowledge or estimates of the mobility of the
individuals living in the city. Many models have been developed to this
effect[1, 2, 3, 4], focusing on different scales of the urban system,
according to the quality of the data available. Traditionally, these models
allocate trips from one geographical zone to another, according to estimates
of where people live and work. Infrastructure projects are then assessed
following changes in accessibility which are computed from the model’s
predictions of the people living and working in these zones. Among the models
that have been proposed over the decades, the gravitational model [4] has been
one of the most widely adopted in various contexts (for example, [5, 6, 7, 8,
9, 10]). Depending on the information available on the demographics and
mobility of the individuals, this model exists in the following different
forms as an unconstrained, singly constrained, and doubly constrained model.
Naturally, the more information available to calibrate the model, the better
its performance can be assessed against observed data, notwithstanding the
fact that information is not always available for this purpose.

Recent years have seen a dramatic increase in available data on individuals
and their urban environments, allowing researchers to test these models more
effectively, thus providing more detail on the outstanding problems of human
mobility. This has prompted a surge in the literature, where new models have
been proposed, such as the _Radiation Model_ [11]. This model takes its
inspiration from the intervening opportunities model [1] where flows are
modelled without parameters and take only as input the population
distribution. It produces predictions with a high degree of accuracy at the
intra-county scale, hence introducing a new benchmark in the field of mobility
modelling. This has triggered the interest of many researchers, and many works
have appeared where its predictions are compared with those of the more
traditional gravity model [12, 13, 14, 15, 16, 17, 18]. These efforts have
focused on comparing the accuracy of the models in reproducing observed
origin-destination matrices. As shown in [15, 17, 18] the main limitations of
the Radiation Model is in its inability to produce adequate outcomes at
different spatial scales which is a direct consequence of its own virtue of
being parameter free. In order to overcome this limitation, several solutions
have been proposed: notably in [18] the authors introduce a _normalised_
version of the model in order to take into account of finite size effects,
while in [17] the authors propose an _Extended_ version of the Model
introducing a parameter that can be calibrated to the data. These studies show
how in general, the accuracy in reproducing the observed flows of the
_extended_ and the _normalised_ version are often comparable to those obtained
from the doubly-constrained gravitational model.

The origin-destination matrices used to compare the models in previous works,
are extracted both by conventional datasets, as for example mobility surveys
and census data, or unconventional datasets such as mobile phone or geo-
located social media data. Very often though the matrices are outdated,
incomplete or obtained by biased samples of the population [19]. We therefore
believe that the models need to be compared in a more practical manner where
their results can be read and interpreted more easily. For this reason, here
we want to take a different approach and exploit the two different models to
quantify and predict the impact of the introduction of a new Bus Rapid Transit
(BRT) system in the city of Teresina in Brazil. BRT systems are increasing in
popularity worldwide as an alternative cost-effective investment in comparison
to expensive urban rail transport projects [20]. Readers can see from
http://brtdata.org/ that there are more than 206 cities which have introduced
some kind of BRT system and with the number of new corridors under
construction increasing steadily. Indeed emerging economies have been seduced
by the publicized BRT success from cities as Curitiba and Bogotà [21, 22, 23,
24], cities that developed high performance with BRT systems, enhancing
mobility and sustainability at an affordable price. For this reason,
implementations of BRT systems are now being evaluated by many cities [25, 26]
and there are many studies of these proposals[27, 28, 29]. With this in mind,
the main goal of this paper is to use both the Gravitational and the Radiation
models to quantify the effects of such transportation interventions, measuring
the variation in accessibility, and comparing the two results in terms of an
analysis of their similarities and differences.

##  1 The Case Study: The BRT implementation in Teresina

As mentioned above, our case study will be the city of Teresina in Brazil, a
medium size metropolis which is currently implementing a BRT system. It is the
capital of the state of Piauì and its metropolitan region has almost
$1.2\text{x}10^{6}$ inhabitants according to the last census estimate made by
the ”Instituto Brasileiro de Geografia e Estatística”in 2015. The
administrative region, named ”Região Integrada de Desenvolvimento da Grande
Teresina - RIDE/Grand Teresina”, is composed of 15 municipalities. However,
only two municipalities are served by a connected urban transport network -
Teresina and Timon (see Fig. 1). Both cities together concentrate most of the
population in the region with just over $10^{6}$ inhabitants. The metropolitan
public transportation consists of a bus network, and a rail service connecting
the southeast zone to the city centre. The rail service operates sharing the
infrastructure with freight trains on a single track in both directions,
resulting in a sparse and low usage service. For this reason, it will not be
considered for cost comparisons in this study. The present structure of the
public transport system is thus non-hierarchical. The majority of the lines
form a radial scheme, departing from the suburbs towards the city centre with
a few services that directly connect zones in the suburbs. In recent years,
several exclusive bus lanes in the central area have been constructed which
aim to reduce travel time on congested roads.

Figure 1:  (Left Panel) A map of Teresina, zones 1 to 55, and Timon, zones 56
to 64, where zones are coloured according to their average income from white,
low income, to dark blue, high income. (Right panel) A map of the new BRT
stops and their corridors. The BRT stops are only found in Teresina but
mobility in Timon is clearly affected by the scheme.

In 2008, the municipality of Teresina approved its Transport and Mobility
Master Plan which proposes a new system of public transport in the city. It
relies on an Origin-Destination survey of trip making conducted in 2007 which
was developed to analyse travel patterns and predict future scenarios [30].
The original proposal suggests the implementation of a Bus Rapid Transit
system (BRT), splitting the existing routes and creating a hierarchical system
of feeder, inter-terminal and trunk services. This will be composed of eleven
terminals connected through express bus corridors. The proposal aims to
increase the effectiveness of public transportation in the city, improving the
accessibility to jobs, education and public services. In this paper, through
different measures of accessibilities, we evaluate the impact of such an
infrastructure project.

##  2 Methods

In this section we will present and briefly recount the details of the models
we have used to estimate the impact of the introduction of the BRT system in
Teresina and Timon, while also introducing the equations we have used to
measure the accessibilities. As we will see, we have taken into account the
journey to work distribution using both models and infrastructure to calculate
the accessibilities.

###  2.1 The Gravity Model

In the gravity model approach, the flow from zone $i$ to zone $j$ is
proportional to the opportunities,employment in this case, $E_{j}$ in
destination $j$, and to the demand in origin $i$ (which is represented by the
population) $P_{i}$, and weighted by the cost function $f(c_{ij})$. Given our
approach, we consider the cost $c_{ij}$ of going from $i$ to $j$ as the
expected time of travel using the public transportation network (see appendix
for details on the travel time calculations), and the function as being an
exponential decay of the form $f(c_{ij})=\text{e}^{-\beta c_{ij}}$, where
$\beta$ is a parameter that has to be calibrated on data. Traditionally to
model the journey to work trip distribution, the total number commuters is
constrained (outflow), as are the employees arriving at work (inflow). This
corresponds to the doubly-constrained model, where the flows are described by
the equation

$T^{\text{dbl}}_{ij}=A_{i}B_{j}P_{i}E_{j}\text{e}^{-\beta c_{ij}}$ (1)

where $A_{i}$ and $B_{j}$ are two normalization constants that one has to
solve iteratively. By imposing the constraints on the total outflow and on the
total inflow

$\sum_{j}T_{ij}=P_{i}\qquad\qquad\qquad\sum_{i}T_{ij}=E_{j}$ (2)

and following the procedure in [4], we get

$A_{i}=\frac{1}{\sum_{k}B_{k}E_{k}\text{e}^{-\beta c_{ik}}}\qquad\qquad
B_{j}=\frac{1}{\sum_{k}A_{k}P_{k}\text{e}^{-\beta c_{kj}}}$ (3)

In its single constrained version, where only the constrain on the outflow is
kept, the flow from origin $i$ and destination $j$ in this case has the form

$T^{\text{sng}}_{ij}=Z_{i}P_{i}E_{j}\text{e}^{-\beta c_{ij}}$ (4)

where $Z_{i}$ is the normalisation constant. Imposing the constraint on the
outflow $\sum_{j}T_{ij}=P_{i}$ leads to

$Z_{i}=\frac{1}{\sum_{k}E_{k}\text{e}^{-\beta c_{ik}}}$ (5)

In the next section we will use both forms to reproduce the origin-destination
matrix and to calculate the accessibilities.

###  2.2 The Radiation Model

The original radiation model [11] takes its inspiration from the intervening
opportunities model [1]. In this approach the probability of commuting between
two units $i$ and $j$ depends on the number of opportunities between the
origin and the destination, rather than on their distance. In its original
formulation the radiation model made use of the population in each zone, using
it also as a proxy for employment. The flow between zone $i$ and $j$ is
therefore quantified as

$T^{\text{rad}}_{ij}=T_{i}\frac{P_{i}P_{j}}{(P_{i}+P_{ij})(P_{i}+P_{j}+P_{ij})}$
(6)

where $P_{ij}$ is the population in zones included in a radius of distance
$d_{ij}$ (and excluding those of zones $P_{i}$ and $P_{j}$). These represent
the opportunities between them, and $T_{i}$ is the amount of commuters in $i$.
As presented, the model in eq.(6) was formulated to describe flows happening
on large scales and the absence of parameters to calibrate makes the model
hard to fit to smaller scales. For this reason, the form we have used is
slightly different and following [18] we have added a normalisation constant
that takes into account the _finite size_ of the system. Moreover we have used
the employment to characterize opportunities rather than the population. The
flows between zones $i$ and $j$ are now described by

$T^{\text{rad}}_{ij}=\frac{P_{i}}{(1-\frac{P_{i}}{P})}\frac{E_{i}E_{j}}{(E_{i}+E_{ij})(E_{i}+E_{j}+E_{ij})}$
(7)

where $P$ is the total amount of population in the system, and where $E_{ij}$
is the amount of employment between zones $i$ and $j$.

We have also used the Extended version of this model presented in [17] where a
parameter $\alpha$ is introduced, and whose calibration makes the model
adaptable to different spatial scales. The flow in this version is derived by
combining the original Radiation model with _survival analysis_[31] and in
this context the probability of commuting from $i$ to $j$ is described by the
equation

$P(1|E_{i},E_{j},E_{ij})=\frac{\left[(E_{i}+E_{j}+E_{ij})^{\alpha}-(E_{i}+E_{ij})^{\alpha})\right](E_{i}^{\alpha}+1)}{\left[(E_{i}+E_{ij})^{\alpha}+1\right]\left[(E_{i}+E_{j}+E_{ij})^{\alpha}+1\right]}$
(8)

The details of the calculations that lead to this form may be found in [17].
The flows in the extended version are the product eq.(8), the population in
the origin zone $i$, and normalisation term

$T^{\text{ext}}_{ij}=P_{i}\frac{P(1|E_{i},E_{j},E_{ij})}{\sum_{k}P(1|E_{i},E_{k},E_{ik})}$
(9)

One may notice that as per our construction, eq.(9) is constrained to meet the
outflow $\sum_{k}T_{ik}=P_{i}$, but not the inflow $\sum_{k}T_{ki}\neq E_{i}$
exactly as in the singly-constrained version of the gravity model.

###  2.3 Measures of Accessibility

As mentioned in the introduction, in order to quantify the impact of the new
infrastructure, we measure the accessibility predicted by the models before
and after the introduction of the BRT, which is done by using both
$c^{\text{old}}_{ij}$ and the updated cost matrix $c^{\text{brt}}_{ij}$.
Accessibility has become a central concept in physical planning in the past
decades [32, 33, 34, 35, 36], and many different definitions exist which
depend on the specific application. In general as stated in [35],
accessibility associates some measure of opportunity at a place with the cost
of actually realising that opportunity, or in other words as the cost of
getting to some place traded off against the benefits received once that place
is reached. In [35], two main types of accessibility are defined: type 1 takes
into account the locational behaviour described through models, with
infrastructure only implicitly considered; type 2, considers the physical
infrastructure and some generalised measure of the _distance_ from the zone of
interest to all others. For a zone $i$ we define accessibility of type 1 as

$A^{1}_{i}=\frac{\sum_{j}T_{ij}(1/c_{ij})}{\sum_{j}T_{ij}}$ (10)

where $c_{ij}$ is the cost of commuting from $i\rightarrow j$ (this is a
measured quantity which depends on the infrastructure), and $T_{ij}$ is the
predicted flow and depends on the model used to calculate it. The
accessibility in eq.(10) quantifies the inverse of the average cost of
commuting from the given area: high values of $A_{i}$ correspond to the flows
happening with low values of $c_{ij}$ and viceversa. We define, as a simple
measurement of the infrastructure, the accessibility of type 2 as

$A^{2}_{i}=\frac{1}{N_{d}}\sum_{j}\frac{E_{j}}{c_{ij}}$ (11)

where $E_{j}$ are the opportunities (employment) in zone $j$, and $N_{d}$ is
the total number of destinations. Once again eq.(11) is telling us the average
benefit-cost ratio for zone $i$. As we can see, there is no modelling involved
in this measure, and we have only exploited the cost matrix and the
opportunities distribution. Comparing eq.(10) and eq.(11) allows us to
understand the information benefit of adding a modelling layer to the
analysis, especially given the use of the different kinds of model we have
applied.

##  3 Results

In the first part of this section, we will produce the trip duration
distribution with the different models and calibrate it to the trip data
(details of which with respect to the origin-destination matrix data are to be
found in the appendix). As said, the reproduction of the observed origin
destination matrices is usually how the goodness of fit of such a model is
tested. From this analysis, it will be clear that the normalized radiation
model is not appropriate in fitting mobility problems at this scale (as it has
not been developed with intra city trips in mind), and will therefore, we will
not consider this model further. In the second part, we will compare in great
detail the results obtained when calculating the impact on the accessibility
of the new BRT bus network. We will do this only using the extended version of
the radiation and the gravitational models. Finally, to understand these
results, we study the characteristics of the zones whose predictions are
similar using the two approaches and compare them with the zones for which
they differ.

###  3.1 Trip Duration Distributions

To calibrate the parameters of the models we have used we used the maximum
likelihood method. As one can see in the appendix the origin-destination
survey counts only 138 households, and therefore does not contain information
on trips from every zone to every other. For each origin zone, we only have
information on trips to a limited number of destinations. This implies that we
cannot calibrate the models using the common part of commuters based on the
Sorensen’s index [37] as done in previous comparisons [15, 16, 17, 18]. We
therefore extract the trip duration distribution of the bus trips contained in
the dataset [30] , $p_{\text{obs}}(t)$ and for each value of the parameter
which we generically indicate as $p$, we measure the distribution of the
predicted duration of the trips $p_{\text{mod}}(t,p)$. We then find the values
of the parameters that maximise the likelihood expression

$L(p)=\sum_{t}p_{\text{obs}}(t)\log{(p_{\text{mod}}(t,p))}$ (12)

where $p=\alpha$ for the extended radiation and $p=\beta$ for the gravity
model.

Figure 2: In the top panels (a) and (b) we show the behaviour of the
likelihood for the gravity and the extended radiation models. If the gravity
models exhibit a clear maximum for $\beta_{sng}=0.022$, $\beta_{dbl}=0.026$,
in the extended radiation context this is reached asymptotically for
$\alpha\rightarrow 0$. In (c) we compare the trip duration distribution found
in the data with what comes out from the models.

In fig.(2)(a)-(b) we show the form of the likelihood for the doubly and singly
constrained gravity model, and for the extended radiation model. For both
versions of the gravity approach the likelihood presents a clear maximum with
$\beta_{single}=0.022$ and $\beta_{double}=0.026$, values we have used to
produce the results shown in this paper.

When repeating the calculation in eq.(12) for the extended radiation model, we
see how the maximum is asymptotically reached for $\alpha\rightarrow 0$. Given
that we are working at an intra-city scale this is not surprising, and in the
supplementary material of [17](in section 9) the authors have solved the
model’s form as eq.(9) for this limit with these scales in mind. The equation
that describes the flow from zone i to zone j for $\alpha\rightarrow 0$
becomes

$\lim_{\alpha\rightarrow
0}T^{\text{ext}}_{ij}=P_{i}\frac{n_{j}/(n_{i}+s_{ij})}{\sum_{k}n_{k}/(n_{i}+s_{ik})}$
(13)

In fig.(2)(c) we show the trip duration distribution obtained with the
different models. The black dashed curve indicates the trip duration
distribution found in the data $p_{\text{obs}}(t)$, the red and orange curves
represent the singly and doubly-constrained gravity models, the normalised and
extended radiation the cyan and blue curves.

It clear that with no calibration process, the normalised radiation is not in
agreement with what is seen in the data, and how the distribution obtained
with the extended version of the model, in the limit $\alpha\rightarrow 0$,
moves closer to what seen in the data. If on the one hand, this indicates that
indeed the extended version of the radiation model makes it a better fit for
simulating mobility at small scales, these results seem to be yet another
proof of superiority of the gravity approach when dealing with urban mobility
problems. We will therefore exclude the normalised radiation from the
following analysis, and concentrate on the radiation model in its extended
form.

###  3.2 Accessibility Variations after the Introduction of the BRT

Figure 3: In this figure we show the accessibility to employment with
$A_{\text{old}}$ in (a) and $\frac{A_{\text{brt}}}{A_{\text{old}}}$ in (b) on
the top left and top right respectively, for each zone of the city of Teresina
and Timon. As noted in the legends of the figures, the red and orange curves
show the results obtained from the gravity model, from the doubly and singly
constrained versions respectively, while the blue curve represents the results
obtained with the extended radiation model. The black curve in the
accessibility variation is the variation measured with respect accessibility
of type 2. In the insets, we present the $R^{2}$ values between the various
curves. The maps show the spatial distribution of $A^{\text{var}}$ obtained
with the doubly-constrained gravity model and with the extended radiation
model.

In order to quantify the impact of the introduction of the BRT, we have
measured, for all zones in Teresina, the quantities in eq.(10) and (11) using
the cost matrices before the intervention $c^{\text{old}}_{ij}$ and after
$c^{\text{brt}}_{ij}$. To explicitly study the variation introduced by the BRT
system, we have then analysed the ratio of the two quantities

$A_{\text{var}}=\frac{A_{\text{brt}}}{A_{\text{old}}}$ (14)

so that the zones with $A_{\text{var}}(i)>1$ are predicted to benefit from the
intervention and vice-versa. We repeat this process using the gravitational
and radiation approach. The results are summarized in fig(3), where we show in
detail all the quantities we have discussed and the spatial distribution of
the predicted impact on the city’s zoning system. In the top panel we have
shown $A_{\text{old}}$ and $A_{\text{var}}$ calculated with the various
models. Given the difference in the characteristic values between the
accessibilities of type 1 and 2, we only compare their predictions on the
$A_{\text{var}}$, where the variation of type 2 is represented by a black
curve. We will refer to $A(T^{\text{dbl}})$, $A(T^{\text{sng}})$ and
$A(T^{\text{ext}})$ for the accessibilities calculated using the double-
constrained gravity, single-constrained gravity and the extended model
respectively.

What immediately catches the eye is the substantial agreement between the
$A(T^{\text{ext}})$ (blue curve) and $A(T^{\text{dbl}})$ (red curve) both
before the BRT introduction, and then in predicting the impact of its
introduction. The $R^{2}$ values in the insets confirm this impression, with
higher values between the two models than that measured between
$A(T^{\text{dbl}})$ and $A(T^{\text{sng}})$. This is quite unexpected,
especially considering the fact that per construction, the extended radiation
model is a singly constrained model.

Figure 4: In the top figure (a) we show the distribution of $\Delta r$. The
blue zones highlight zones in category 1, the white area zones in category 2,
and the red in category 3. In the middle figures, we show the spatial
distribution of the 3 categories (b), and in (c) for category 1 and 3 we show
the average distance to employment. The shaded areas represent the standard
deviation of the distribution. Finally in the bottom figures, we present for
each category the distribution of population (d) and employment (e)

Looking at the accessibility variations $A_{\text{var}}$, it is clear that all
three models predict two main peaks, one for zones 17 and 16, and another for
zone 26. The first peak can be seen analysing both types of accessibility,
$A^{1}$ and $A^{2}$ (black curve), while the peak in zone 26 is not captured
by the measure of the infrastructure. The BRT map in fig.(1) shows that a new
BRT stop is planned near zone 26, so a spike in its accessibility is indeed
reasonable. This seems to suggest that the behavioural layer introduced
through the modelling of the flows adds more information to the simple
analysis of infrastructure. The variation of the accessibilities of type 1
predicted by the three models have a high $R^{2}$ value ($R^{2}>0.8$), while
$A^{1}(T^{\text{sng}})$ is the only one in good agreement with $A^{2}$.

To understand the spatial configuration of the results, we show a heat map of
the $A_{\text{var}}$ on the map of the city in fig.(3), only using the doubly
constrained and the extended radiation models. The maps show that the zones in
the city centre are expected to benefit from the new transport network, but
the main benefits are located in the south zones 16 and 17, and in the centre
north zones 26 and 27. We can see from the BRT map in fig.(1) that these zones
are close to where new bus stops will be positioned. The zones in the north of
the city, despite being close to the new stops and served by a BRT, are
expected to generate longer journeys to work, and therefore their
accessibility will be lower. The zones in the district of Timon are not
included in the BRT intervention, and it is therefore not surprising to find
that their employment accessibilities decrease according to all the three
models.

Despite producing different aggregated travel time distributions (fig.(2)),
the predictions obtained with the two models are in good agreement. This is
clear by looking at the high $R^{2}$ values obtained when comparing the
curves, and its especially true for the doubly constrained and the extended
radiation models. Moreover the results are reasonable bearing in mind the
spatial distribution of the BRT stops and corridors, and accessibility of type
1 seems a more appropriate measure for this study than the measure of type 2.

###  3.3 The Doubly-Constrained Gravity versus The Extended Radiation Model

As portrayed in the introduction, one of the main objectives of this work is
to analyse similarities and differences in the predictions made using the two
approaches. Despite the great agreement found between the doubly-constrained
and the extended radiation models, several zones show contradictory results.
We will now look into the properties of these specific zones, and check if
there are similarities among them. We do this by comparing the predictions on
the accessibility after the BRT introduction, $A_{\text{brt}}(T^{\text{ext}})$
and $A_{\text{brt}}(T^{\text{dbl}})$, by analyzing the difference in the rank
of each zone, namely

$\Delta r_{i}=r^{\text{dbl}}_{i}-r^{\text{ext}}_{i}$ (15)

where $r^{\text{dbl}}_{i}$ is the rank of zone $i$ using the gravity and
$r^{\text{ext}}_{i}$ is the rank obtained using the radiation model. If
$\Delta r_{i}>0$, this implies that $r^{\text{dbl}}_{i}>r^{\text{ext}}_{i}$,
which means that the accessibility of zone $i$ ranks higher if we use the
gravitational model to calculate the accessibilities, and vice versa. The
distribution of the quantity in eq.(15) calculated using the two models is
found in the top panel of fig.(4). The red and blue areas highlight the zones
for which the predictions vary considerably ($\|\Delta r\|>2$), which are
those whose characteristics we want to analyse. We have therefore divided the
zones into

  * •

Category 1: $\Delta r_{i}>>0$

  * •

Category 2: $\Delta r_{i}\simeq 0$

  * •

Category 3: $\Delta r_{i}<<0$

Furthermore we have studied several characteristics of the zones that belong
to each category: the population and employment distribution, their spatial
distribution and their distance to employment. As we can see from fig.(4), the
population and employment distributios are completely comparable in all three
categories (d)-(e), and no precise spatial pattern emerges when projecting
these categories on the city’s zoning system (b). But indeed an interesting
difference emerges when quantifying the distance to employment of the
different categories, as we can see from fig.(4)(c). What emerges is a clear
tendency of the zones of category 1 to be closer to the opportunities (blue
curves) than those of category 3 (red curves). By looking at the figure, it is
clear how the blue curve tends to increase faster than the red one. A zone
whose closest neighbouring zones are rich in employment tends to have a larger
predicted accessibility using the gravitational approach, while zones
surrounded by zones with low opportunities score higher accessibility using
the radiation approach. The reasons behind this result are probably to be
found in the details of the equations, and as to how the methods deal with the
concept of distance. We believe that this finding can be very useful in
understanding future applications of the two methods.

##  4 Discussion

In this work, we have used two different approaches to test the impact of a
new transportation facility. In doing so, we have achieved two results: the
first of course is to enable us to cross-check our findings, but the second is
comparing the gravity and radiation models models on real data and in
practice. We have seen how the observed trip duration distribution can be
reproduced quite well by the gravity models, in both their versions, while we
also see how the radiation model fails to provide as a good and as a
satisfactory description. The extended radiation model owing to the fact that
it is calibrated improves this model substantially. Interestingly the
agreement between the extended radiation and the double-constrained gravity
model, with different agglomerated trip duration distributions, is
outstanding. This may suggest that the detailed reproduction of the observed
distribution is not crucial in a city planning context with respect to
measuring accessibility variations at the scale we are working at. Furthermore
we have compared the estimated impact on the accessibilities, using the type 1
and 2 definitions. The type 1 accessibility, in this context, seems to be in
better agreement with what one may expect. By looking at fig.(3) we see that
the zones with the highest estimated improvement, zones 15, 16 and 26, have
low starting accessibilities and are positioned next to a BRT stop. The
modelling of flows does seem to add information that the policy maker cannot
extract simply by observing the differences in the infrastructure and the trip
durations. With this in mind the models that seem to better describe mobility
are the extended radiation and the doubly-constrained gravity model. We have
seen how, despite considerable agreement, the predictions obtained by the two
models differ for those zones with many opportunities around them and those
with very few. Indeed we have seen that the zones with a neighbourhood rich in
opportunities perform better with the gravity model, and vice-versa. This is
an interesting result, which is potentially useful in urban planing scenarios.
We may conclude by saying that these preliminary results show that the
radiation model in its extended version could be a valid alternative to study
urban mobility and test new transportation networks. More research on this
topic would help us better understand how the two models could support each
other.

##  5 Appendix

###  5.1 Data

The main data source is the STRANS/PMT (Teresina Transport Authority) database
created through an Origin-Destination survey conducted in 2007 including the
cities of Teresina and Timon [30] . The database includes household
conditions, personal socioeconomic information and travel diaries per person
on the day before the survey was taken. The main dataset contains trip data
(walking times, waiting times, travel times, mode, origin and destination
zones and activities, and trip costs) combined with disaggregated demographic
and socioeconomic information about the traveler (education level, income,
gender, employment, age) and traveler’s household data (traffic zone, comfort
and deprivation variables - number of families, bedrooms, bathrooms, sewage
system, access to water, energy consumption). The survey consists of 64
traffic zones in Teresina and Timon that coincide with the cities’ districts
which provide the opportunity to gather sociodemographic data from the
national census. In total, the dataset contains 5,177 journeys distributed
across 138 households. The dataset also contains information about the number
of employers, students and total population in each traffic zone. Geo-
referenced data is also available for bus routes and stops for all routes in
Teresina and Timon and the General Transit Feed Specification (GTFS) for the
public transportation in the city. Data about Timon’s buses routes were taken
from Moovit App. The Bus Journeys Dataset was also taken from STRANS/PMT and
contains 45,090 bus trips for a twenty-four hour interval for each bus line in
2006. The database describes single bus passengers journey through the
variables: Bus Route, Bus ID, Direction, Time at Origin, Origin Bus Stop and
Destination bus stop.

###  5.2 Cost Function Calculations

The effects of the new transportation network on trip distribution have been
measured by calculation of the generalized costs of travel, for each pair of
zones, $i\rightarrow j$. To do so we had to calculate two cost matrices, one
before $\\{C^{\text{old}}_{ij}\\}$ and one after $\\{C^{\text{brt}}_{ij}\\}$
the introduction of the BRT. Each element of the matrices indicate the time
necessary to travel from zone $i$ to zone $j$ using public transportation. We
have collected the travel times of the old transportation network from Google
Maps Directions, using the API, which provides a service to calculate
directions between requested locations considering available modes of
transportation. The API request takes as inputs origins and destinations,
which were taken from zone centroids coordinates, and gives an estimate of the
expected time of travel. The results were stored in the Origin Destination
cost matrix $C^{\text{old}}$. On the other hand the travel times after the BRT
introduction had to be predicted. To estimate new routes, we have built a
model using ArcGIS Spatial Analyst Extension for the existing bus network with
the current local bus GTFS (General Transit Feed Specification, given by the
Transport Authority of Teresina), and calibrated with real travel time data
and the predicted travel time between pairs of zones was calculated. Further
developments for the BRT systems were introduced in the network (BRT GTFS),
and the predicted travel time between pairs of zones was calculated. In this
case the results were stored in the origin destination cost matrix
$C^{\text{brt}}_{ij}$.

Data Availability: All the data used in this work can be found in this public
repository https://github.com/ducciopiovani/Data-on-Teresina-, where we have
stored the cost matrices and the origin and destination matrix in 3 separate
files.

Competing Interests:We declare we have no competing interests.

Authors’ Contributions: All authors conceived and designed the study, analyzed
and interpreted the results and drafted the manuscript. D.P. carried out the
analysis, written all the codes both for the simulations and the data
analysis, and coordinated the efforts. G.U. gathered and manipulated the data,
built the cost matrices and wrote the first draft.

Funding: All authors acknowledge the funding from the Engineering and Physical
Sciences Research Council, grant number: EP/M023583/1.

## References

  * [1] Samuel A Stouffer.  Intervening opportunities: a theory relating mobility and distance.  American sociological review, 5(6):845–867, 1940.
  * [2] Alex Anas.  Discrete choice theory, information theory and the multinomial logit and gravity models.  Transportation Research Part B: Methodological, 17(1):13–23, 1983\.
  * [3] Moshe E Ben-Akiva and Steven R Lerman.  Discrete choice analysis: theory and application to travel demand, volume 9.  MIT press, 1985.
  * [4] Alan Geoffrey Wilson.  The use of entropy maximising models, in the theory of trip distribution, mode split and route split.  Journal of Transport Economics and Policy, pages 108–126, 1969\.
  * [5] Woo-Sung Jung, Fengzhong Wang, and H Eugene Stanley.  Gravity model in the korean highway.  EPL (Europhysics Letters), 81(4):48005, 2008.
  * [6] James E Anderson.  The gravity model.  Annu. Rev. Econ., 3(1):133–160, 2011.
  * [7] Jameel Khadaroo and Boopen Seetanah.  The role of transport infrastructure in international tourism development: A gravity model approach.  Tourism management, 29(5):831–840, 2008.
  * [8] Gautier Krings, Francesco Calabrese, Carlo Ratti, and Vincent D Blondel.  Urban gravity: a model for inter-city telecommunication flows.  Journal of Statistical Mechanics: Theory and Experiment, 2009(07):L07003, 2009.
  * [9] Duygu Balcan, Vittoria Colizza, Bruno Gonçalves, Hao Hu, José J Ramasco, and Alessandro Vespignani.  Multiscale mobility networks and the spatial spreading of infectious diseases.  Proceedings of the National Academy of Sciences, 106(51):21484–21489, 2009.
  * [10] Alan Wilson.  Boltzmann, lotka and volterra and spatial structural evolution: an integrated methodology for some dynamical systems.  Journal of The Royal Society Interface, 5(25):865–871, 2008.
  * [11] Filippo Simini, Marta C González, Amos Maritan, and Albert-László Barabási.  A universal model for mobility and migration patterns.  Nature, 484(7392):96–100, 2012.
  * [12] Chaogui Kang, Yu Liu, Diansheng Guo, and Kun Qin.  A generalized radiation model for human mobility: Spatial scale, searching direction and trip constraint.  PloS one, 10(11):e0143500, 2015.
  * [13] Filippo Simini, Amos Maritan, and Zoltán Néda.  Human Mobility in a Continuum Approach.  PLoS ONE, 8(3), 2013.
  * [14] Xiao Liang, Jichang Zhao, Li Dong, and Ke Xu.  Unraveling the origin of exponential law in intra-urban human mobility.  arXiv preprint arXiv:1305.6364, 2013.
  * [15] Maxime Lenormand, Sylvie Huet, Floriana Gargiulo, and Guillaume Deffuant.  A Universal Model of Commuting Networks.  PLoS ONE, 7(10), 2012.
  * [16] Maxime Lenormand, Aleix Bassolas, and José J. Ramasco.  Systematic comparison of trip distribution laws and models.  Journal of Transport Geography, 51:158–169, 2016.
  * [17] Yingxiang Yang, Carlos Herrera, Nathan Eagle, and Marta C González.  Limits of predictability in commuting flows in the absence of data for calibration.  Scientific reports, 4, 2014.
  * [18] A. Paolo Masucci, Joan Serras, Anders Johansson, and Michael Batty.  Gravity versus radiation models: On the importance of scale and heterogeneity in commuting flows.  Physical Review E - Statistical, Nonlinear, and Soft Matter Physics, 88(2):1–8, 2013.
  * [19] Maxime Lenormand, Miguel Picornell, Oliva G Cantú-Ros, Antonia Tugores, Thomas Louail, Ricardo Herranz, Marc Barthelemy, Enrique Frias-Martinez, and Jose J Ramasco.  Cross-checking different sources of mobility information.  PLoS One, 9(8):e105184, 2014.
  * [20] Robert Cervero.  Bus rapid transit (brt): An efficient and competitive mode of public transport.  IURD Working Paper 2013-01, 2013.
  * [21] Daniel A Rodríguez and Felipe Targa.  Transport Reviews, 24(5):587–610, 2004.
  * [22] Daniel A Rodriguez and Carlos H Mojica.  Land value impacts of bus rapid transit: the case of bogotá?s transmilenio.  Land Lines, 2008.
  * [23] Luis Antonio Lindau, Dario Hidalgo, and Daniela Facchini.  Curitiba, the cradle of bus rapid transit.  Built Environment, 36(3):274–282, 2010.
  * [24] Herbert S Levinson, Samuel Zimmerman, Jennifer Clinger, and G Scott Rutherford.  Bus rapid transit: An overview.  Journal of Public Transportation, 5(2):1, 2002.
  * [25] Robert Cervero and Chang Deok Kang.  Bus rapid transit impacts on land uses and land values in seoul, korea.  Transport Policy, 18(1):102–116, 2011.
  * [26] Taotao Deng and John D Nelson.  Bus rapid transit implementation in beijing: An evaluation of performance and impacts.  Research in Transportation Economics, 39(1):108–113, 2013.
  * [27] David A Hensher and Thomas F Golob.  Bus rapid transit systems: a comparative assessment.  Transportation, 35(4):501–518, 2008.
  * [28] Graham Currie.  The demand performance of bus rapid transit.  Journal of Public Transportation, 8(1):3, 2005.
  * [29] Kristen Day, Lawrence Loh, Ryan Richard Ruff, Randi Rosenblum, Sean Fischer, and Karen K Lee.  Does bus rapid transit promote walking? an examination of new york city’s select bus service.  Journal of physical activity and health, 11(8):1512–1516, 2014\.
  * [30] Plano diretor de trasportes e mobilidade urbana de teresina.  http://semplan.teresina.pi.gov.br/wp-content/uploads/2017/03/PLANO-DIRETOR-DE-TRANSPORTE-E-MOBILIDADE-URBANA-DE-TERESINA.pdf.
  * [31] Rupert G Miller Jr.  Survival analysis, volume 66.  John Wiley & Sons, 2011.
  * [32] Luca Bertolini, Frank Le Clercq, and Loek Kapoen.  Sustainable accessibility: a conceptual framework to integrate transport and land use plan-making. two test-applications in the netherlands and a reflection on the way forward.  Transport policy, 12(3):207–220, 2005.
  * [33] Grégory Vandenbulcke, Thérèse Steenberghen, and Isabelle Thomas.  Mapping accessibility in belgium: a tool for land-use and transport planning?  Journal of Transport Geography, 17(1):39–53, 2009.
  * [34] Karst T Geurs and Bert Van Wee.  Accessibility evaluation of land-use and transport strategies: review and research directions.  Journal of Transport geography, 12(2):127–140, 2004.
  * [35] Michael Batty.  Accessibility: in search of a unified theory.  Environment and Planning B, 36, 191-194., 2009.
  * [36] SR Jones.  Accessibility measures: a literature review.  Technical report, 1981.
  * [37] Thorvald Sørensen.  A method of establishing groups of equal amplitude in plant sociology based on similarity of species and its application to analyses of the vegetation on danish commons.  Biol. Skr., 5:1–34, 1948.