# Use of Simulation Models for the Development of a Statistical Production
Framework for Mobile Network Data with the simutils Package111Presented at the
Conference Use of R in Official Statistics 2021, 24-26 November 2021,
Bucharest (Romania).

B. Oancea  Dept. Applied Economics and Quantitative Analysis, University of
Bucharest
90, Panduri Street
Sector 5, 050663 Bucharest, Romania  D. Salgado  Dept. Methodology and
Development of Statistical Production, Statistics Spain (INE)
Avda. Manoteras, 52
28050 Madrid, Spain  Dept. Statistics and Operations Research, Complutense
University of Madrid
Pza. Ciencias, 3 - Ciudad Universitaria
28040 Madrid, Spain  S. Barragán  Dept. Methodology and Development of
Statistical Production, Statistics Spain (INE)
Avda. Manoteras, 52
28050 Madrid, Spain  M. Necula  Dept. Innovative Tools in Official Statistics,
National Institute of Statistics (INS)
16 Libertatii Bvd.
Sector 5, 050706 Bucharest, Romania

###### Abstract

We propose to use agent-based simulation models for the development of
statistical methods in Official Statistics, especially in relation with the
new digital data sources. We present a mobile network data simulator which is
managed through the simutils R package which provides geospatial
representations of the simulated data. While the synthetic data are produced
by an external tool, our simutils package allows an R user to parameterize and
run this external simulation tool, to build geospatial data structures from
the simulation output or to compute several aggregates. The geospatial data
structures were designed with the purpose of using them in a visualization
package too. Useful simulation models require the incorporation of real
metadata from mobile telecommunication networks driving us to the inclusion of
functionalities allowing the user to specify and validate them. All metadata
are specified using XML file whose structure are defined in corresponding XSD
files. Our R package includes example data sets and we show here how validate
the metadata, how to run a simulation and how build the geospatial data
structures and how to compute different aggregates.

##  1 Introduction

###  1.1 Background

During the last fifteen years the production of official statistics in
national and international statistical offices is facing a double challenge,
namely the modernization and industrialization of the production process and
the incorporation of new digital data sources. The first challenge was clearly
identified in the first decade of the present century (HLG-MOS, 2011) and gave
rise to the advent of international production models such as the GSBPM
(UNECE, 2019b), the GSIM (UNECE, 2019c), and the GAMSO (UNECE, 2019a), to name
a few. These standards are being gradually adopted by most statistical offices
thus conforming a common working space in the international community (UNECE,
2021). The second challenge was already identified eight years ago (see e.g.
DGINS, 2013) and is taking efforts in the form of complementary institutional
initiatives (see e.g. Eurostat, 2014, 2021; UNECE, 2021). This challenge “will
require amendments to the statistical business architecture, processes,
production models, IT infrastructures, methodological and quality frameworks,
and the corresponding governance structures […]” (DGINS, 2018), i.e. basically
to modify all aspects of the entire production framework.

Multiple analyses have been conducted regarding the use of new digital data
for the production of official statistics (see Kitchin, 2015; Hand, 2018;
Salgado and Oancea, 2020, and multiple references therein). Most of them are
still far from being solved. Apart from the identification of concrete issues
such as the data access, the new methodological and quality assurance
frameworks, and the technological environment, more entagled implications
arise because of their interaction. More specifically, these new data sources
require new statistical methods in the production of official statistics
partially abandoning and hopefully integrating the use of sampling designs
(Särndal et al., 1992) with statistical models (Valliant et al., 2000),
especially in the fields of statistical learning (Hastie et al., 2009; Murphy,
2013), Bayesian modelling (Gelman et al., 2013) or geostatistics
(Schabenberger and Gotway, 2005). This new official statistical methodology is
developing too slow, among other things, because the access to this sort of
data is proving to be extremely hard and new legislation will need to be put
into force (Baldacci et al., 2021). Thus, strategies to boost the development
of new statistical production frameworks in this context are needed. As we
shall explain in the following sections, the use of simulation models and
their analysis with R brings an excellent scenario to face this situation.

###  1.2 The contribution of this paper

A basic tenet of our present contribution is that simulation models producing
synthetic data incorporating real metadata must become a fundamental tool to
develop these new methodologies and to assess their potentiality. Simulated
data must help develop the methodological and quality frameworks as well as
assess the technological environment to device a new statistical production
chain. This should be combined with the use of real data in such a way that
synthetic data are to be finally substituted by them to constitute the ever
evolving production cycle. This idea is illustrated in figure 1. Even when
having access to real data, simulated data provides an analytical and testing
playground of different aspects of any kind of statistical model used to
produce the official statistics.

Figure 1: Schematic representation of the strategy to develop and evolve a
statistical production framework.

As we shall argue below, the use of R stands as an excellent statistical
software development tool to implement this strategy. Following Chambers
(2009), this versatility arises mainly because R provides outstanding
characteristics for the development and evolution of statistical production
frameworks: (i) it is an _interface_ to computational procedures of many
kinds, (ii) it is _interactive_ , (iii) it incorporates _functional_
programming facilities, (iv) it incorporates _object-oriented_ programming
facilities, (v) it is _modular_ , and (vi) it is _collaborative_. The need to
interact with other computing languages and platforms is beyond doubt a need
for any modern computing system. The interactivity allows the user to check in
a fast and seamless way the results of intermediate computations. Functional
and object-oriented paradigms are essential ingredients in a computing system
to process statistical data and to produce reliable statistics. Modularity is
the key feature for a scalable system to cope with complexity. The
collaboration of the international community in an open environment is
fundamental to reach widely adopted production standards. There already exists
an important initiative to share open source statistical software for the
production of official statistics (ten Bosch et al., 2020).

In this contribution we focus on the use of mobile network data as a novel
promising data source for the production of official statistics. This data
source is indeed facing a shortage of data access for Official Statistics due
to intricate reasons, thus the use of synthetic data will allow us to move
forward in the development of statistical methodology.

The work is organised as follows. Firstly, we describe from a statistical
point of view the complex data ecosystem of mobile telecommunication networks
for the production of statistics. Next, we provide generic information about
the use of simulation models, especially in connection with mobile network
data. Then, we describe the general lines of the statistical software
development with the creation of several R packages to implement a modular
end-to-end process to produce population counts and origin-destination
matrices from this data source. We focus on the package simutils, which
provides a set of functionalities to manage simulated mobile network data. We
conclude with a brief summary and discussion.

##  2 Ecosystem of mobile network data for statistical production

Mobile telecommunication networks are extraordinarily complex systems
providing communication between moving agents through electromagnetic
interactions (see e.g. Miao et al., 2016). These networks daily produce a huge
amount of digital data to make this communication feasible and to allow
engineers monitor the situation to provide an increasingly higher-quality
communication service. As it is widely known and perceived, this kind of
communication involves personal data to some extent (in order to establish the
desired contact) and, thus, its use for producing official statistics raises
more often than not some concern. In this line of thought, it is mandatory to
explain as clear as possible what data are to be used to produce what
statistics.

Firstly, data to be used are generated in the network systems, never ever in
the mobile devices. That is to say, the goal is to reuse the digital trace
left behind in the mobile telecommunication networks by the mobile devices,
i.e. by individuals of a given target population of analysis. Hence, the more
appropriate expression _mobile network data_ instead of _mobile phone data_.

Secondly, what should be understood by “digital trace”? Nowadays two main
usages of mobile network data are considered, namely geolocation of mobile
devices (see Ahas et al., 2007; Salgado et al., 2021, for a pioneering work
and an extensive list of references, respectively) and Internet usage (Ucar et
al., 2021) both for sociodemographic analyses. In the former case,
telecommunication variables driving us to an estimation of the position of
devices through some statistical model are needed (Munoz et al., 2009). In the
latter case, network variables showing mobile digital usage are needed (Ucar
et al., 2021). As we shall argue, analyses with synthetic data may allow us to
identify those variables needed to reach the accuracy determined by official
statistics quality standards.

Thirdly, should we use individual or aggregated data? For official statistics
purposes, this is a misleading question. Aggregated mobile network data, if
any, arises from a computational procedure upon individual mobile network
data. This computational procedure is of statistical nature, especially
regarding the geolocation of mobile devices and more intensively when referred
to administrative territorial units (city districts, municipalities, …). This
procedure is executed by someone in a given computing platform following a
previously designed algorithm. Thus, the starting point is always individual
mobile network data since the network provides service at an individual level.
The question should be: are official statistics to be based on some form of
raw digital data or on some form of preprocessed statistics (i.e. of mobile
network statistics, we should say)? This question is far from having a
definitive answer and touches the role of Official Statistics in society. In
this work we concentrate on the proposal of executing in-situ production
processes hopefully designed jointly by statistical officers and
telecommunication experts (Salgado et al., 2021).

In this line of reasoning we shall work with network event data, i.e. data
generated due to the interaction between a mobile device and the network. We
shall denote this data by $\mathbf{E}_{kt}$, where $k$ denotes the mobile
device ID and $t$ denotes the timestamp generated by the interaction (with a
precision of seconds). Complementarily, we shall also make use of
telecommunication variables about the configuration of the network such as the
position of each antenna, their emission power, their path loss exponent,
…(Salgado et al., 2021). We shall denote this data by $\mathbf{N}_{a}$, where
$a$ stands for the antenna222We assume for simplicity’s sake that the
configuration of the network does not change during the period of analysis..
We used the simplified term antenna for what mobile network engineers cal a
Base Transceiver Station. These data allow us to compute an estimation of the
signal strength from each antenna or a similar measure in each tile of a
territorial grid of analysis through a so-called radio propagation model (Miao
et al., 2016). Notice that the information contained in $\mathbf{E}_{kt}$ will
be sensitive for users whereas the information contained in $\mathbf{N}_{a}$
will be sensitive for mobile network operators (hence the intricate complexity
of access to this data source).

Network event data can comprise multiple telecommunication variables, some of
them with no utility for our statistical purposes and in constant evolution
with the technology. For statistical purposes regarding geolocation, we shall
simulate datasets collecting the following information from each network
event: (i) the timestamp of the network event, (ii) a constant identification
variable for the device during the time period of analysis, (iii) an
identification variable for the radio cell where the event is produced, (iv) a
code for the type of event, (v) the network technology type (3G, 4G, etc.) and
(vi) the _Timing Advance_ of the signal (basically a code denoting the time
for the signal to reach the antenna from a mobile device). We complement these
variables with the simulated real position of the device where the network
event is produced (coordinates and grid tile identification). This synthetic
ground truth will allow us to assess the performance of any statistical model
of geolocation and of estimation of population counts and mobility patterns
(see figure 2).

Figure 2: Example of simulated network event data.

Regarding the network configuration parameter datasets, we collect variables
regarding (i) antenna (radio cell) identification, (ii) the mobile network
operator identification, (iii) the mobile network operator name, (iv) the
maximum number of connections per radio cell, (v) the emission power (in W),
(vi) the path loss exponent, (vii) the type of radio cell (omnidirectional,
120∘-directional), (viii) the minimum signal strength to establish a
connection (in dBm), (ix) the minimum signal dominance (as defined in Tennekes
and Gootzen (2021)) to establish a connection, (x) the midpoint parameter for
the computation of the signal dominance (in dBm), (xi) the steepness parameter
for the computation of the signal dominance, (xii-xviii) diverse orientation
parameters for the antenna (tilt, elevation, etc.), and (xix) the grid tile
identification of the antenna position (see figure 3).

Figure 3: Example of simulated network parameter configuration.

Both network event and network parameter variables are generated through a
network event data simulator developed in a modular form for this purpose, as
explained below.

From the perspective of the statistical analysis to be conducted, this wealth
of complex and technological data $\mathbf{E}_{kt}$, $\mathbf{N}_{a}$ should
not overwhelm the statistician since the ultimate goal is to build a
statistical model to compute the location probability
$\mathbb{P}\left(\textrm{tile}|\mathbf{E}_{kt},\mathbf{N}_{a}\right)$ for
every single mobile device $k$ at each time instant $t$ of analysis. These
location probabilities allow us to detach the analysis from the extremely
technological substratum so that deduplication of devices carried by the same
individual, identification of target individuals (tourists, commuters, etc.),
aggregation to produce probability distributions of population counts, and
inference to the whole target population can be undertaken on a statistical
basis (Tennekes et al., 2020; Tennekes and Gootzen, 2021; Salgado et al.,
2021; Ogulenko et al., 2021a). More sophisticated models involving another
telco variables can also be possibly devised, hence the hopeful collaboration
between statistical offices and mobile network operators. Synthetic data and
comparison with the simulated ground truth will allow the statistician to
assess the model performance and the quality indicators (accuracy, bias,
etc.).

##  3 Data simulation for the development of statistical production frameworks

A scientific model can be defined as a “simplified representation of a system
under study, which can be used to explore, to understand better or to predict
the behaviour of the system it represents” (O’Sullivan and Perry, 2013). A key
word in this definition is “simplified”, which is a conceptual need arising
from the scientific exercise itself of understanding and describing reality.
This is also present in the use of any type of data, either simulated or real.
As a matter of fact, despite the huge deluge of data inviting to some kind of
objective apprehension of reality thereof, this is not so: some sort of
conceptual, mathematical, phsyical, or scientific modelling exercise must be
undertaken on these data, all during generation, collection, processing, and
inference. Both real and simulated data are _a sociotechnical artefact_
(Kitchin, 2017). Hence, the well-known motto “Essentially, all models are
wrong…but some are useful” by Box (1979) (see also Wit et al., 2012, and
articles in the same monographic journal volume number). Notice, however, that
we have multiple ways of expressing the goodness of fit of a model (how wrong
it is), but we still lack the formal idea of a model utility (how useful it
is), despite being a fundamental question (see e.g. Field, 2015, for a
concrete discussion related to earthquake risk modelling). In our context,
simulation models are cleary oriented towards the development of statistical
methods driving us to the identification and selection of telecommunication
variables, geolocation methods, and statistical procedures in general
producing official statistics of interest for the public good. Official
Statistics is not about real data, but about an accurate description of
reality for policy-making and decision-taking.

This concept of model utility, however it is defined, is related to how
scientists use models, which O’Sullivan and Perry (2013) illustrate in terms
of two critical aspects, namely data availability and phenomenon understanding
(see figure 4). When we lack data and have a poor understanding of a
phenomenon, simulation models can be used to gain learning. On the other hand,
when both data and understanding abound, we can use models to make accurate
predictions. Currently, the production system of official statistics based on
mobile network data is clearly in a learning stage, hence the natural use of
simulation models.

Figure 4: Data versus understanding (from O’Sullivan and Perry (2013))

In this context, agent-based models (ABMs) stand as an excellent tool to
explore and test different proposals to create a statistical process based on
mobile network data. They are simulation models that describe the behaviour of
individual organisms. In the past, the complexity of these kinds of models
caused a lack of use. However, in recent years due to the computational
progress they have undergone an expansion. These models are used in a wide
range of fields such as Ecology (DeAngelis, 2018), Economics (Tesfatsion,
2002), Social Science (Gilbert and Troitzsch, 2005) and much more. The
importance of ABMs is obvious since they allow to represent complex systems
and analyse them in multiple and different scenarios. Even more in the case of
the development of statistical methodology, these models are essential to
assess the quality of the process and make comparisons among different
methods.

Simulated populations have already been traditionally used in Official
Statistics (e.g. in survey sampling to assess the performance of finite-
population design-based estimators). However, ABMs are not an extended tool,
among other things due to their complexity, to computational issues, and to
interpretation challenges. This is not exclusive in Official Statistics and
many other fields also face the same problem. With the aim of facilitating the
use and understanding of these models a standard protocol was proposed some
years ago (Grimm et al., 2006). This protocol basically conveys the key
elements of an ABM by making explicit seven items grouped in three blocks: (i)
overview (purpose, state variables and scales, process overview and
scheduling), (ii) design concepts, (iii) details (initialization, input,
submodels). The idea is to present first the general considerations, then the
strategy and finally the technical details (see Grimm et al., 2006, for
details). If ABMs are to be used in Official Statistics, this protocol or an
adaptation thereof should be agreed and followed.

There exists a methodology to create simulation models (see e.g. Law, 2008,
and multiple references therein) but, in our view, there exists a fundamental
aspect in their use for Official Statistics both from a technical and a
strategic point of view. Despite being necessarily a simplification of the
final production system, simulation models must incorporate _real_ metadata
into their construction so that simulated agents and their interacting context
resemble as much as possible real-life conditions. This provides a minimal
guarantee of reproducing some aspects of reality. However, the goal is to
develop versatile and robust statistical methods capable of reliably
estimating the synthetic ground truth, simulated with multiple
parameterizations of the metadata. If a new ground truth is simulated, the
statistical model should be able to estimate it. This is also paramount for
the assessment of the quality framework for the devised statistical process.

This incorporation of metadata, in our view, should be an essential element
for the collaboration between statistical offices and mobile network
operators, both at the development and the execution stages of the
incorporation of this data source into the production of official statistics.

##  4 simutils: an R package to manage simulated mobile network data

In this paper we show the importance of trustworthy software in the
development of statistical production frameworks. As stated by Chambers
(2008), “users of the analysis have no option but to trust the analysis, and
by extension the software that produced it. Both the data analyst and the
software provider therefore have a strong responsibility to produce a result
that is trustworthy, and, if possible, one that can be shown to be trustworthy
[…]”. Moreover, regarding software prototyping, Davis (1995) said that "The
best way to assess what users really need is to give them a “working” system
[…]".

Bearing these ideas in mind, we have chosen the R language due to the numerous
advantages that make it a trustworthy software, for instance: (i) it is _open
source_ , (ii) it is _user-friendly_ for statisticians and data scientists,
(iii) it is _agile enough_ to allow a reasonable speed _in implementing new
prototypes_ , (iv) there is a _wide range of packages and functions_ well
developed and related to statistical analysis, especially the latest
improvements on geocomputation functionalities (Lovelace et al., 2022), apart
from the essential feature of taking into account the functional and object-
oriented paradigms.

###  4.1 A modular approach: the software stack

Modularity is a key element in any production process, thus also in
statistical production (Baldwin and Clark, 2000; Esteban et al., 2018). The
idea of modularity is to build the whole process as a set of fully independent
modules where any building block follows the generic structure of a
statistical production step (see figure 5). The modularity allows us to
update, modify, and adapt each module according to the statistical domain,
statistical needs, and data ecosystem we face in each case. The modularity
also allows us to approach the issue about data access with a higher degree of
flexibility.

Figure 5: Building block for the statistical process (adapted from van der Loo
(2019))

The R language is without any doubt a good choice for doing implementation of
a modular process as it is presented by van der Loo (2021).

The big picture of the whole process to analyse mobile network data goes from
the raw telecommunication data (real or simulated) to the final estimation of
population counts and origin-destination matrices along a fully modular
process. The methodology framework of this process has been developed in the
European project called ESSnet on Big Data (Eurostat, 2021). The
implementation has been done mostly in R with some preliminary prototypes for
each process step. In particular, five big modules were proposed, namely, (i)
the geolocation of network events, (ii) the device multiplicity deduplication,
(iii) the statistical filtering of target individuals, (iv) the aggregation of
mobile devices, and (v) the inference over the target population (see figure
6).

Figure 6: The layered structure of the software implementation (the
statistical filtering module is missing because we are estimating the whole
population, i.e. all individuals are target)

In relation to this project the developed software is described by Oancea et
al. (2020), Oancea et al. (2021) (see also ESSnet on Big Data, 2021). There is
a lot of work in progress regarding a general improvement of all the packages
presented as prototypes, leaving for future work the issue of scalability for
the application in real production conditions.

###  4.2 Mobile network data simulator

In the context of the European project ESSnet on Big Data, a mobile network
data simulator has been developed in C++ (Oancea et al., 2019; ESSnet on Big
Data, 2021), whose constant evolution is available as open source software in
github (Oancea, 2021). This simulator needs a set of input files in XML format
and a geographical map or polygon alike in WKT format. These files contain
data about (i) the simulation (time span information, mobility patterns,
probabilities for individuals of having one or two mobile devices), (ii)
persons (how many individuals comprise the population, how fast they walk or
move by car, how long they spend at home, at work, etc.), (iii) antennas
(their position, emission power, path loss exponent, and other configuration
parameters), and (iv) prior specifications for the computation of location
probabilities (see figure 7).

Once the simulator is executed with all these parameters, the following
outputs can be obtained: (i) the grid of analysis in which the geographical
territory is divided, (ii) the ground truth about the persons at each time
instant (their positions and their devices), (iii) the parameters of each
antenna, (iv) the network events arising from the connections of the devices
to each antenna at each time instant, (v) the coverage area for each antenna,
and (vi) the signal strength (or the signal dominance) of each antenna in the
center of each grid tile. In figure 7 the whole structure of the simulator’s
inputs and outputs is represented.

The structure of each file is specified in a corresponding dictionary in XML,
whose contents, in turn, are specified also in a corresponding XSD file. These
file structures allow the user to represent the metadata by using a standard
format in XML as well as to conduct their validation with the dictionaries in
XSD format. This validation is done using a Java program integrated in the
simutils package through the rJava package (Urbanek, 2021) allowing for
compatibility with some features introduced in XML ver. 1.1 (features not yet
in R packages dealing with XML processing).

Figure 7: Inputs and Outputs of the simulator

###  4.3 The simutils package

The execution pipeline to deal with all the information related to the
simulator is shown in figure 8. Although the simulator has been implemented in
C++, all the information related to it can be managed with the R language by
using the simutils package, as it is shown in this section. This package is
available in github at https://github.com/bogdanoancea/simutils and is
currently under constant evolution. With this tool, it is possible to execute
the simulator as well as to manage, specify and change their inputs and
outputs. More, the outputs of the simutils package can be sent to simviz,
another R package that we created to build different visualizations of the
individual or aggregated data regarding the mobile devices and persons or have
a visual representation of different aspects of the mobile network and its
characteristics.

Figure 8: Execution pipeline MND

Next, we show how to execute step by step the pipeline to obtain simulated
datasets as well as some preliminary computations with the ground truth data
(total population counts and origin-destination matrices). All files used in
the following examples are available together with the simutils package.

The first step is to specify and validate the input files for the simulator.
This is seamlessly accomplished using the XML/XSD file structures mentioned
above:

rootPath <\- file.path(system.file(package = "simutils"), ’extdata’)# Validate
the simulation.xml file with its corresponding xsd filexml_fn <\-
file.path(rootPath, "input_files", "simulation.xml")xsd_fn <\-
file.path(rootPath,"metadata/input_files/schema_definition","simulation_dict.xsd")validate_xml(xsd_fn,
xml_fn)# Validate the antennas.xml file with its corresponding xsd filexml_fn
<\- file.path(rootPath, "input_files", "antennas.xml")xsd_fn <\-
file.path(rootPath,"metadata/input_files/schema_definition","antennas_dict.xsd")validate_xml(xsd_fn,
xml_fn)

Next, any modification of the simulator input files can be possibly undertaken
with the update_simInput() function:

# rootPath is taken from the preceding code snippetxmlSimInput <\-
file.path(rootPath, ’extdata/input_files’, ’simulation.xml’)newParam.lst <\-
list(end_time = 11,movement_pattern = structure(list(manhattan_grid =
structure(list(x_step = 40, y_step = 40, x_origin = 0, y_origin = 0),type =
’home_work_manhattan’))))xsdName <\-
file.path(rootPath,’extdata/metadata/input_files/schema_definition’,’simulation_dict.xsd’)#
New simulation input file to be created in working directorynewSimInputFile
<\- file.path(getwd(), ’newSimulation.xml’)update_SimInput(xmlSimInput =
xmlSimInput,newParam = newParam.lst,xsdName = xsdName,newFileName =
newSimInputFile)

The next step is to run the simulation, but firstly we need the proper
configuration of the simulator on your computer. The simulator has to be
installed separately. It is distributed as source code but we also provide an
installation kit. The source code can be downloaded from
https://github.com/bogdanoancea/simulator/ and the executable file can be
built using a C++ compiler. Detailed instructions on how to build the
simulator for Windows, MacOs or Linux operating systems are provided on the
above mentioned github page. The easiest way to install it under Windows is to
download the installation kit from the following link:
https://github.com/bogdanoancea/simulator/releases/tag/1.2.0-kit. Assuming the
default installation path, once installed, it is possible to run the
simulation from R with the runExeFile() function as follows:

# rootPath is taken from the preceding code snippetrun_ExeFile(path_to_exe =
"C:/Program Files/Simulator",simulator_version = "1.2.0",input_folder =
file.path(rootPath, "extdata/input_files"),simulationCFGFile =
"simulation.xml",mapFile = "map.wkt",personsCFGFile =
"persons.xml",antennasCFGFile = "antennas.xml")

The simulation tool is also distributed as a docker image. If one wants to use
the docker image instead of the executable file, the runDockerImage() function
should be run:

# rootPath is taken from the preceding code snippetinput_folder <\-
file.path(rootPath, "extdata/input_files")output_folder <\-
file.path(Sys.getenv(’HOME’), ’example_docker’)run_DockerImage(input_folder =
input_folder,simulationCFGFile = ’simulation.xml’,mapFile =
’map.wkt’,personsCFGFile = ’persons.xml’,antennasCFGFile =
’antennas.xml’,output_folder = output_folder)

This function automatically pulls the latest version of the simulator’s docker
image from the central docker hub.

Once the simulation has been produced, the resulting files can be read with
the read_simData() function as follows. We mention that the map.xml file used
here augments the initial map of the simulation by adding subdivisions
(territorial units) of the territory under consideration.

filename_map <\- c(xml= system.file("extdata/input_files", "map.xml", package
= "simutils"),xsd= ’’)filename_network <\- c(csv=
system.file("extdata/output_files/antennas.csv",package = "simutils"),xml=
system.file("extdata/metadata/output_files/antennas_dict.xml",package =
"simutils"))filename_signal <\- c(csv=
system.file("extdata/output_files/SignalMeasure_MNO1.csv",package =
"simutils"),xml=
system.file("extdata/metadata/output_files/SignalMeasure_dict.xml",package =
"simutils"))filename_coverage <\- c(csv= system.file("extdata/output_files",
"AntennaCells_MNO1.csv",package = "simutils"),xml=
system.file("extdata/metadata/output_files/AntennaCells_dict.xml",package =
"simutils"))filename_events <\- c(csv=
system.file("extdata/output_files/AntennaInfo_MNO_MNO1.csv",package =
"simutils"),xml=
system.file("extdata/metadata/output_files/events_dict.xml",package =
"simutils"))filename_grid <\- c(csv=
system.file("extdata/output_files/grid.csv",package = "simutils"),xml=
system.file("extdata/metadata/output_files/grid_dict.xml",package =
"simutils"))filename_individ <\- c(csv=
system.file("extdata/output_files/persons.csv",package = "simutils"),xml=
system.file("extdata/metadata/output_files/persons_dict.xml",package =
"simutils"))filenames <\- list(map = filename_map,network_parameters =
filename_network,signal = filename_signal,events =
filename_events,coverage_cells = filename_coverage,grid =
filename_grid,individuals = filename_individ)simData <\-
read_simData(filenames, crs = 2062)

The output of this function is the workhorse of the package. It is a list of
spatial objects from the stars and sf packages (Pebesma, 2021, 2018) used to
represent the territory in the map, the network information, the coverage of
each antenna, the network events variables and the signal strength in the grid
as well as the ground truth about the individuals. The choice of these
packages to represent geospatial objects is complemented with the choice of
the data.table package (Dowle and Srinivasan, 2021) to allow for an efficient,
fast and high-volume data management. The structure of this object can be
observed in figure 9.

Figure 9: Objects in simData

Another functionality of this package is to compute some aggregates of the
ground truth (total population counts and origin-destination matrices) that
can be computed with the compute_total() and compute_odmatrix() functions. The
code to obtain each aggregate is included in the following lines:

  * •

Counting individuals by territorial unit (Subregion_long in this example) and
time is undertaken as follows:

compute_odmatrix(individuals.sf = simData$individuals,what = ’individuals’,by
= c(’t’, ’Subregion_long’))

  * •

Counting devices by territorial unit (Subregion_long also in this example) and
time:

compute_odmatrix(individuals.sf = simData$individuals,what = ’devices’,by =
c(’t’, ’Subregion_long’))

  * •

Counting individuals by territorial unit (Subregion_long also in this example)
and time with 0 devices:

compute_odmatrix(individuals.sf = simData$individuals,what =
’individuals_dev0’,by = c(’t’, ’Subregion_long’))

  * •

Counting individuals by territorial unit (Subregion_long also in this example)
and time with 1 device:

compute_odmatrix(individuals.sf = simData$individuals,what =
’individuals_dev1’,by = c(’t’, ’Subregion_long’))

  * •

Counting individuals by territorial unit (Subregion_long also in this example)
and time with 2 devices:

compute_total(individuals.sf = simData$individuals,what =
’individuals_dev2’,by = c(’t’, ’Subregion_long’))

  * •

Counting multiple totals by territorial unit (Subregion_long also in this
example) and time:

totals <\- c(’individuals’, ’individuals_dev0’,
’devices’)compute_odmatrix(individuals.sf = simData$individuals,what =
totals,by = c(’t’, ’Subregion_long’))

The functions described above are the primary functions, but there is a bunch
of secondary functions which are needed to deal with the information. In
figure 10 there is a graph with the dependency relations among all the
functions in the simutils package. Primary functions are represented in green.

Figure 10: Function dependencies in the simUtils package (drawn with the
DependenciesGraphs package (Robert, 2016)).

In figure 11 there is a graph representing the packages dependencies.

Figure 11: Packages dependencies (drawn with the depgraph package (Aust,
2021))

The outputs of simutils can be visualized with the simviz package, which is
under early construction. Firstly, we have prioritised the choice and
implementation of the representation of these simulated network data (the
simulator’s inputs and outputs) into R geospatial objects to facilitate next
both the statistical analyses and the visualization with standard tools
already offered in the R ecosystem. The visualization of these objects will be
undertaken in the simviz package making use of the grammar of graphics
functionalities providing combined visualizations of the map and their
territorial units, the grid, the antenna coverage, the network events and the
individuals movements. The computation and statistical analysis (e.g. the
geolocation of mobile devices) will be undertaken in other packages (see e.g.
Sanguiao et al., 2020), which we are refactoring to include the aforementioned
representation of geospatial objects.

##  5 Summary and discussion

The advent of new data sources in Official Statistics is also bringing the
need for new statistical methods to infer about target populations with the
traditional high-quality standards. This new data sources are mainly digital
data with the characteristics of Big Data, thus rendering their access,
management and processing remarkably complex. In this line, the development of
new statistical production frameworks becomes thus a challenge, since the use
of real data to research stands as very costly (in a wide sense) or even
impossible.

To face this situation, we propose to use agent-based simulation models, which
can provide a wealth of data to experiment and to test. Computer simulations
are not free of controversy (see Edmonds and Meyer, 2013; Saam, 2017, and
references therein) and one can rightfully argue that a mobile
telecommunication network and the behaviours of its subscribers and the human
population, in general, is too complex to be modelled in an agent-based
simulation. This is basically the distinction between _verification_ and
_validation_ (Winsberg, 2009) and it does not invalidate the utility of
computer simulations.

A strong point for simulations is the provision of a synthetic ground truth
_for each parameterization_ , thus allowing the analyst to investigate and
test robust methods valid for a wide range of potential situations (never
known in real conditions). A necessary element is to provide scenarios which
are realistic enough, which can be undertaken by incorporating real metadata
into the simulations (e.g. realistic network configurations and mobility
patterns). This points once more at a necessary collaboration with data
holders not in terms of data transfers or data transmissions at this stage,
but to share critical information for the creating of an industrialised
statistical process for the public good.

The use of simulation models under these conditions, in particular, with
mobile network data allows us also to tackle methodological questions in a
(semi-)empirical way which otherwise would consume many resources. Take for
example the recent study by Ogulenko et al. (2021b) empirically showing the
fallacy of the closest antenna to estimate the geolocation of mobile devices,
invalidating thus many proposals based on Voronoi tessellations of the
geographical territory of analysis (with implications also for the disclosure
of information and the reidentification of individuals). Shortcomings from
this simple use of Voronoi tessellations was already identified as early as
2015 by Ricciato et al. (2015). This kind of situation can be replicated and
thus analysed with detail using the simulator presented here.

Therefore, for the development of statistical methodology, we see immediate
advantages of using simulation models: (i) to have actual ground truth figures
allowing us to conduct a thorough performance assessment of methods and
parameters and a better understanding by comparison between actual population
counts and their estimates, (ii) to identify different concrete aspects of the
problem by configuring different scenarios in order to propose specific
elements in the methodology to deal with them, (iii) to avoid the issue about
the access to real data (see above) and its consequences (lack of data,
confidentiality and privacy risks, legal concerns,…), and (iv) to provide a
body of technical knowledge to reach informed partnership agreements between
statistical offices and data holders. Real data cannot provide these
conditions for research.

In general terms, the use of agent-based simulation models enjoys attractive
features for the research on novel statistical methodologies, which we should
be aware about and which are shared in other scientific disciplines.

These advantages are complemented with the use of R to deploy in a fast and
trustworthy manner those software tools needed to test the novel methods. R
can deal with many different features regarding these new data sources,
namely, the geospatial representation of information, the visualization of
combined information, the use of different statistical methods (statistical
learning models, state space models, …) and reporting.

In this line, we have developed the simutils package to manage all the
geospatial information around these specific data simulations. Our package
provides means to parameterize and run simulations using and external tool
specially designed for this purpose, to manage the output data of the
simulation building geospatial objects needed to carry on specific analyzes
and to compute some aggregates. We showed concrete examples of the main
functionalities using data sets included in the package.

## References

  * Ahas et al. (2007) R. Ahas, A. Aasa, S. Silm, and M. Tiru.  Mobile Positioning Data in Tourism Studies and Monitoring: Case Study in Tartu, Estonia.  In M. Sigala, L. Mich, and J. Murphy, editors, _Information and Communication Technologies in Tourism 2007_ , pages 119–128, Vienna, 2007. Springer.
  * Aust (2021) F. Aust.  _depgraph: Plot package dependency graph_ , 2021.  URL https://github.com/crsh/depgraph.  R package version 0.1.0.
  * Baldacci et al. (2021) E. Baldacci, F. Ricciato, and A. Wirthmann.  A Reflection on The Re(Use) of New Data Sources for Official Statistics.  _Revista Índice_ , 83:8–11, 2021.
  * Baldwin and Clark (2000) C. Baldwin and K. Clark.  _Design Rules (vol.1): The power of modularity_.  MIT Press, 2000.
  * Box (1979) G. Box.  Some problems of statistics and everyday life.  _Journal of the American Statistical Association_ , 74:1–4, 1979.
  * Chambers (2008) J. Chambers.  _Software for data analysis_.  Springer, 2008.
  * Chambers (2009) J. Chambers.  Facets of R.  _The R Journal_ , 1(1):5–8, 2009.
  * Davis (1995) A. M. Davis.  Software prototyping.  In _Advances in computers_ , volume 40, pages 39–63. Elsevier, 1995\.
  * DeAngelis (2018) D. L. DeAngelis.  _Individual-based models and approaches in ecology: populations, communities and ecosystems_.  CRC Press, 2018.
  * DGINS (2013) DGINS.  The Scheveningen Memorandum.  Technical report, European Union, 2013.  https://ec.europa.eu/eurostat/cros/news/scheveningen-memorandum-big-data-and-official-statistics-adopted-essc_en.
  * DGINS (2018) DGINS.  The Bucharest Memorandum.  Technical report, European Union, 2018.  https://ec.europa.eu/eurostat/web/european-statistical-system/-/dgins2018-bucharest-memorandum-adopted.
  * Dowle and Srinivasan (2021) M. Dowle and A. Srinivasan.  _data.table: Extension of ‘data.frame‘_ , 2021.  URL https://CRAN.R-project.org/package=data.table.  R package version 1.14.2.
  * Edmonds and Meyer (2013) B. Edmonds and R. Meyer, editors.  _Simulating Social Complexity: A Handbook_.  Springer, 2013.
  * ESSnet on Big Data (2021) ESSnet on Big Data.  Software for Mobile Network Data in Official Statistics, 2021.  https://github.com/MobilePhoneESSnetBigData/.
  * Esteban et al. (2018) E. Esteban, M. Novás, S. S. na, D. Salgado, and L. Sanguiao.  Data organisation and process design based on functional modularity for a standard production process.  _Journal of Official Statistics_ , 34(4):811–833, 2018.
  * Eurostat (2014) Eurostat.  Vision 2020 Implementation Portfolio, 2014.  URL http://ec.europa.eu/eurostat/web/ess/about-us/ess-vision-2020/implementation-portfolio.
  * Eurostat (2021) Eurostat.  ESSnet on Big Data, 2021.  https://ec.europa.eu/eurostat/cros/content/essnet-big-data-1_en.
  * Eurostat (2021) Eurostat.  Big Data, 2021.  URL https://ec.europa.eu/eurostat/cros/content/big-data_en.  Collaboration in Research and Methodology for Official Statistics.
  * Field (2015) E. Field.  “All Models Are Wrong, but Some Are Useful”.  _Seismological Research Letters_ , 86(2A):291–293, 2015.
  * Gelman et al. (2013) A. Gelman, B. Carlin, H. Stern, D. Dunson, A. Vehtari, and D. Rubin.  _Bayesian data analysis_.  CRC Press, 2013.
  * Gilbert and Troitzsch (2005) N. Gilbert and K. Troitzsch.  _Simulation for the social scientist_.  McGraw-Hill Education (UK), 2005.
  * Grimm et al. (2006) V. Grimm, U. Berger, F. Bastiansen, S. Eliassen, V. Ginot, J. Giske, J. Goss-Custard, T. Grand, S. K. Heinz, G. Huse, et al.  A standard protocol for describing individual-based and agent-based models.  _Ecological modelling_ , 198(1-2):115–126, 2006\.
  * Hand (2018) D. Hand.  Statistical challenges of administrative and transaction data.  _Journal of the Royal Statistical Society A_ , 181:555–605, 2018.
  * Hastie et al. (2009) T. Hastie, R. Tibshirani, and J. Friedman.  _The Elements of Statistical Learning_.  Springer, New York, 2009.
  * HLG-MOS (2011) HLG-MOS.  Strategic vision of the High-Level Group for strategic developments in business architecture in Statistics.  _Conference of European Statisticians, Geneva, 14-16 June_ , 2011\.
  * Kitchin (2015) R. Kitchin.  Big Data and Official Statistics: Opportunities, challenges and risks.  _Statistical Journal of the IAOS_ , 31:471–481, 2015.
  * Kitchin (2017) R. Kitchin.  Thinking critically about and researching algorithms.  _Information, Communication & Society_, 20:14–29, 2017\.
  * Law (2008) A. Law.  How To Build Valid and Credible Simulation Models.  In S. Mason, R. Hill, L.Mönch, O. Rose, T. Jefferson, and J. Fowler, editors, _Proceedings of the 2008 Winter Simulation Conference_ , pages 39–47, 2008.
  * Lovelace et al. (2022) R. Lovelace, J. Nowosad, and J. Muenchow.  _Geocomputation with R_.  CRC Press, 2022.  https://geocompr.robinlovelace.net/.
  * Miao et al. (2016) G. Miao, J. Zander, K. Sung, and S. Slimane.  _Fundamental of Mobile Data Networks_.  Cambridge University Press, 2016.
  * Munoz et al. (2009) D. Munoz, F. Bouchereau, C. Vargas, and R. Enríquez-Caldera.  _Position Location Techniques and Applications_.  Academic Press, Burlington, 2009.
  * Murphy (2013) K. Murphy.  _Machine learning: a probabilistic perspective_.  MIT Press, 2013.
  * Oancea (2021) B. Oancea.  Mobile Network Data Simulator, 2021.  https://github.com/bogdanoancea/simulator.
  * Oancea et al. (2019) B. Oancea, M. Necula, D. Salgado, L. Sanguiao, and S. Barragán.  Deliverable i.2 (data simulator) a simulator for network event data.  Technical report, Eurostat, 2019.  URL https://ec.europa.eu/eurostat/cros/sites/default/files/WPI_Deliverable_I2_Data_Simulator_-_A_simulator_for_network_event_data.pdf.
  * Oancea et al. (2020) B. Oancea, S. Barragán, L. Sanguiao, and D. Salgado.  Deliverable i.4 (information techonologies) some it tools for the production of official statistics with mobile network data.  Technical report, Eurostat, 2020.  URL https://ec.europa.eu/eurostat/cros/system/files/wpi_deliverable_i3_a_proposed_production_framework_with_mobile_network_data_2020_11_26_final.pdf.
  * Oancea et al. (2021) B. Oancea, D. Salgado, L. Sanguiao, and S. Barragán.  A Set of R Packages to Estimate Population Counts from Mobile Phone Data.  _Romanian Statistical Review_ , (1):17–38, 2021.
  * Ogulenko et al. (2021a) A. Ogulenko, I. Benenson, I. Omer, and B. Alon.  Probabilistic positioning in mobile phone network and its consequences for the privacy of mobility data.  _Computers, Environment and Urban Systems_ , 85:101550, 2021a.  ISSN 0198-9715.  doi: https://doi.org/10.1016/j.compenvurbsys.2020.101550.  URL https://www.sciencedirect.com/science/article/pii/S0198971520302830.
  * Ogulenko et al. (2021b) A. Ogulenko, I. Benenson, M. Toger, J. Östh, and A. Siretskiy.  The fallacy of the closest antenna: Towards an adequate view of device location in the mobile network.  _arXiv preprint arXiv:2109.02154_ , 2021b.
  * O’Sullivan and Perry (2013) D. O’Sullivan and G. L. Perry.  _Spatial simulation: exploring pattern and process_.  John Wiley & Sons, 2013.
  * Pebesma (2018) E. Pebesma.  Simple Features for R: Standardized Support for Spatial Vector Data.  _The R Journal_ , 10(1):439–446, 2018.  doi: 10.32614/RJ-2018-009.  URL https://doi.org/10.32614/RJ-2018-009.
  * Pebesma (2021) E. Pebesma.  _stars: Spatiotemporal Arrays, Raster and Vector Data Cubes_ , 2021\.  URL https://CRAN.R-project.org/package=stars.  R package version 0.5-3.
  * Ricciato et al. (2015) F. Ricciato, P. Widhalm, M. Craglia, and F. Pantisano.  Estimating population density distribution from network-based mobile phone data.  Technical report, Joint Research Centre, 2015.  Technical Report JRC96568. Available at https://ec.europa.eu/jrc/file/document/26733.
  * Robert (2016) T. Robert.  _DependenciesGraphs: Dependencies visualization between functions and environments_ , 2016.  URL https://github.com/datastorm-open/DependenciesGraphs.  R package version 0.3.
  * Saam (2017) N. Saam.  What is a Computer Simulation? A Review of a Passionate Debate.  _Journal of General Philosophy of Science_ , 48:293––309, 2017.  doi: 10.1007/s10838-016-9354.  URL https://doi.org/10.1007/s10838-016-9354.
  * Salgado and Oancea (2020) D. Salgado and B. Oancea.  On new data sources for the production of official statistics, 2020.  arXiv:2003.06797v1.
  * Salgado et al. (2021) D. Salgado, L. Sanguiao, B. Oancea, S. Barragán, and M. Necula.  An end-to-end statistical process with mobile network data for official statistics.  _EPJ Data Science_ , 10:20, 2021.  doi: 10.1140/epjds/s13688-021-00275-w.  URL https://epjdatascience.springeropen.com/articles/10.1140/epjds/s13688-021-00275-w.
  * Sanguiao et al. (2020) L. Sanguiao, D. Salgado, and B. Oancea.  _destim: R package for mobile devices position estimation_ , 2020\.  URL https://github.com/Luis-Sanguiao/destim.  R package version 0.1.0.
  * Särndal et al. (1992) C.-E. Särndal, B. Swensson, and J. Wretman.  _Model assisted survey sampling_.  Springer, New York, 1992.
  * Schabenberger and Gotway (2005) O. Schabenberger and C. Gotway.  _Statistical Methods for Spatial Data Analysis_.  Texts in Statistical Science. Chapman & Hall/CRC, Boca Raton, 2005.
  * ten Bosch et al. (2020) O. ten Bosch, M. van der Loo, and A. Kowarik.  The awesome list of official statistics software: 100 …and counting.  The Use of R in Official Statistics Conference, 2-4 December 2020. See also https://github.com/SNStatComp/awesome-official-statistics-software, 2020\.
  * Tennekes and Gootzen (2021) M. Tennekes and Y. Gootzen.  A Bayesian approach to location estimation of mobile devices from mobile network operator data, 2021.  arXiv:2110.00439v1.
  * Tennekes et al. (2020) M. Tennekes, Y. Gootzen, and S. Shah.  A Bayesian approach to location estimation of mobile devices from mobile network operator data.  Technical report, CBS/Statistics Netherlands, 2020.  Working Paper 06-20. Available at https://www.cbs.nl/-/media/_pdf/2020/22/cbds_working_paper_location_estimation.pdf.
  * Tesfatsion (2002) L. Tesfatsion.  Agent-based computational economics: Growing economies from the bottom up.  _Artificial life_ , 8(1):55–82, 2002.
  * Ucar et al. (2021) I. Ucar, M. Gramaglia, M. Fiore, Z. Smoreda, and E. Moro.  News or social media? Socio-economic divide of mobile service consumption.  _Journal of The Royal Society Interface_ , 18:20210350, 2021\.  doi: 10.1098/rsif.2021.0350.  URL https://royalsocietypublishing.org/doi/10.1098/rsif.2021.0350.
  * UNECE (2019a) UNECE.  Generic activity model for statistical organizations v1.2, 2019a.  https://statswiki.unece.org/display/GAMSO/GAMSO+v1.2.
  * UNECE (2019b) UNECE.  Generic Statistical Business Process Model v5.1, 2019b.  https://statswiki.unece.org/display/GSBPM/GSBPM+v5.1.
  * UNECE (2019c) UNECE.  Generic statistical information model v1.2, 2019c.  https://statswiki.unece.org/display/gsim.
  * UNECE (2021) UNECE.  High-level group for the modernisation of statistical production and services, 2021.  https://unece.org/statistics/networks-of-experts/high-level-group-modernisation-statistical-production-and-services.
  * UNECE (2021) UNECE.  Big Data, 2021.  URL https://unece.org/statistics/ces/big-data.
  * Urbanek (2021) S. Urbanek.  _rJava: Low-Level R to Java Interface_ , 2021.  URL https://CRAN.R-project.org/package=rJava.  R package version 1.0-5.
  * Valliant et al. (2000) R. Valliant, A. Dorfmann, and R. Royall.  _Finite population sampling and inference. A prediction approach_.  Wiley, New York, 2000.
  * van der Loo (2019) M. P. van der Loo.  _Systematic approaches to data validation and data cleaning for statistical production_ , 2019.  Seminar at Statistics Spain (INE).
  * van der Loo (2021) M. P. van der Loo.  A Method for Deriving Information from Running R Code.  _The R Journal_ , 13(1):42–52, 2021.  doi: 10.32614/RJ-2021-056.  URL https://doi.org/10.32614/RJ-2021-056.
  * Winsberg (2009) E. Winsberg.  Computer Simulation and the Philosophy of Science.  _Philosophy Compass_ , 4/5:835––845, 2009.  doi: 10.1111/j.1747-9991.2009.00236.x.  URL https://doi.org/10.1111/j.1747-9991.2009.00236.x.
  * Wit et al. (2012) E. Wit, E. van den Heuvel, and J.-W. Romeijn.  ‘All models are wrong…’: an introduction to model uncertainty.  _Statistica Neerlandica_ , 66(3):217–236, 2012\.


# Combining astrophysical datasets with CRUMB

Fiona A. M. Porter & Anna M. M. Scaife
Jodrell Bank Centre for Astrophysics
University of Manchester
{fiona.porter<EMAIL_ADDRESS>
The Alan Turing Institute, 96 Euston Rd, London, UK<EMAIL_ADDRESS>
###### Abstract

At present, the field of astronomical machine learning lacks widely-used
benchmarking datasets; most research employs custom-made datasets which are
often not publicly released, making comparisons between models difficult. In
this paper we present CRUMB, a publicly-available image dataset of Fanaroff-
Riley galaxies constructed from four “parent” datasets extant in the
literature. In addition to providing the largest image dataset of these
galaxies, CRUMB uses a two-tier labelling system: a “basic” label for
classification and a “complete” label which provides the original class labels
used in the four parent datasets, allowing for disagreements in an image’s
class between different datasets to be preserved and selective access to
sources from any desired combination of the parent datasets.

##  1 Introduction

The field of astronomy is entering the era of Big Data astrophysical surveys,
with peta- or exascale volumes of data anticipated with the advent of
instruments such as the Square Kilometre Array [1]. The traditional method of
classifying sources – visual inspection by astronomers – will clearly be
impractical for these surveys, and a natural solution is to turn instead to
machine learning classification. However, astronomical machine learning
research does not at present have standard benchmarking datasets; instead, it
is common for researchers to develop their own datasets for their
publications, most of which are not made publicly available, making it
difficult to reproduce results or compare the performance of different
techniques. While this is being addressed in part by the release of large
labelled datasets by e.g. [2], an additional measure that can be taken is to
build combined datasets by identifying compatible datasets in the existing
literature and integrating them together. This work will discuss the
considerations needed when combining astronomical datasets, the methods used
to construct a merged dataset and the applications of the resulting dataset.

##  2 Dataset construction

For this work we chose to focus upon datasets of Fanaroff-Riley (FR) galaxies
[3], which are active galactic nuclei (AGN) which emit very brightly at radio
wavelengths. These galaxies are classified as either FRI, FRII or hybrid based
on their morphology (see Figure 1), with said morphologies being informative
of the properties of the host AGN engines and surrounding environments.
Despite a number of publications using this population for machine learning
[e.g. 4, 5, 6, 7], existing publicly-available datasets of FR galaxies are
comparatively small; the largest at the time this work began, the MiraBest
dataset, contains only 1256 sources [8]. Creating a combined dataset using
MiraBest and any other suitable datasets of FR galaxies was hence seen to have
two benefits: it would serve as an excellent test case for dataset merging and
also provide an appreciably larger population of FR galaxies within a single
dataset for the astronomical community. We chose to dub this combined dataset
Collected Radiogalaxies Using MiraBest
(CRUMB)111https://doi.org/10.5281/zenodo.7948346.

###  2.1 Selecting datasets to merge

When attempting to merge astronomical machine learning datasets, a key
consideration is whether they draw their image data from surveys with
sufficiently similar properties. Astrophysical sources can have markedly
different morphological characteristics at different wavelengths, so for
consistency in feature detection all image data should be drawn from surveys
of the same wavelength. As well as this, the appearance of sources naturally
varies when observed by instruments with different resolution; a survey with
higher angular resolution can reveal extended structure in what would appear
to be point sources when observed by a survey with lower angular resolution,
but may “resolve out” (and fail to detect) faint extended emission that is
visible in the lower-resolution survey. Again, this can result in
inconsistencies in the appearance of sources of the same class, as they are
effectively being shown at different spatial scales.

A natural solution to these concerns is to simply make use of datasets which
were constructed using the same survey. In radio astronomy, the VLA FIRST
survey [9] is a relatively common choice for dataset construction, being large
and readily accessible; it provides coverage of a large area of the northern
sky at the commonly-used frequency of 1.4 GHz, at a comparably high resolution
of 1.8 arcseconds per pixel, and its data are publicly available. Limiting our
search to extant image datasets of FR galaxies using FIRST data, two datasets
were identified as being suitable for combination with the MiraBest dataset
(hereafter MB) [8] and its supplemental hybrid dataset (hereafter MB-Hyb): the
FR-DEEP dataset [10] and the unnamed dataset of [4] (hereafter AT17).

(a) (b) (c)

Figure 1: Examples of FRI (a), FRII (b) and hybrid (c) morphology from CRUMB.
(a) is found in FR-DEEP and AT17; (b) is found in MB and FRDEEP; (c) is found
in MB-Hyb only.

###  2.2 Cross-matching sources

Astronomical datasets are often derive their source lists from catalogues of
the desired source classes, which provide class labels and coordinates for
each source’s centre; however, different catalogues may disagree on the exact
location of the centre, leading to multiple sets of coordinates being used for
the same source. Since the four “parent” datasets used to create CRUMB are
derived from a total of six different catalogues, it was necessary to cross-
match their sources to remove any duplicates.

To cross-match, the coordinates of all sources in the four datasets were
compared to one another, with pairs of sources being flagged if they were less
than 270 arcseconds (the width of one image) apart. Flagged sources were
visually inspected to determine whether the two sets of coordinates
corresponded to the same source; if so, we retained the best match to the
visual centre as the “preferred” coordinates and discarded any others. Using
this method, a total of 2100 unique sources were identified from the 2731 sets
of coordinates in the parent datasets. Seven of these sources were
additionally found to be poorly centred and were realigned accordingly.

###  2.3 Data processing

Image data for each source were collected from the FIRST survey [9] using
SkyView Virtual Telescope [11], which was used to produce a $300\times 300$
pixel (540") image centred upon the source in FITS format. Radio noise was
removed by using the `astropy` package’s `sigma_clipped_stats` function [12]
to identify all pixels with value $<3\sigma$ and setting them to zero, at
which point the images were cropped to their final size of $150\times 150$
pixels. Next, a circular mask of diameter 150 pixels was applied to remove any
bright background sources present around the edge of the images. Finally, the
source data were normalised by determining the minimum and maximum pixel flux
values for each image and scaling each image pixel as follows:

$Normalised\,pixel\,value=255\times\frac{Pixel\,value-
minimum\,flux}{Maximum\,flux-minimum\,flux}.$ (1)

Here, the factor of 255 allows for maximum dynamic range for PNG format, which
all images were saved as. While keeping the images in FITS format would have
allowed for a higher dynamic range to be preserved, this format is not common
outside of astronomy and the resulting file sizes were found to be two orders
of magnitude larger than the same source in PNG format; PNG was hence
preferred both for accessibility to the broader AI community and reduced file
size.

###  2.4 Creating class labels

Each of the datasets used in the creation of CRUMB makes use of a different
set of class labels for its sources (see Table 1); MB is the most complex,
having a total of ten classes to mark both morphological subclass and
certainty in classification, while AT17 has three classes and both FR-DEEP and
MB-Hyb have two. Determining a unified labelling system is hence somewhat
complex, particularly as the labelling systems of MB and AT17 are somewhat
incompatible, with bent-tailed sources being considered a subclass of FRIs by
MB but an entirely separate class by AT17. Additionally, of the 541 sources
which were present in more than one dataset, 64 were found to have
disagreements in label between the different datasets and 15 had labels which
were entirely contradictory; clearly, there is disagreement event amongst
experts about how some sources should be labelled, and selecting a single
“correct” label is not trivial.

Table 1: Class labels used in each of the parent datasets. “WAT” and “HT” are bent-tail morphologies treated as FRIs by MB, equivalent to “bent” sources in AT17; “DD” are double-double FRIIs. Label | MB | FR-DEEP | AT-17 | MB-Hyb
---|---|---|---|---
0 | Confident FRI | FRI | FRI | Confident hybrid
1 | Confident WAT | FRII | FRII | Uncertain hybrid
2 | Confident HT |  | Bent |
3 | Uncertain FRI |  |  |
4 | Uncertain WAT |  |  |
5 | Confident FRI |  |  |
6 | Confident DD |  |  |
7 | Uncertain FRII |  |  |
8 | Confident hybrid |  |  |
9 | Uncertain hybrid |  |  |

Rather than attempting to merge the disparate labelling systems, we make use
of the alternative approach of a two-tier labelling system: each image is
given both a “basic” label of either FRI, FRII or hybrid and a “complete”
label which contains its class labels in all of the parent datasets. This
compromise allows for CRUMB to offer both an easy-to-understand labelling
system for users who simply wish to train models using astronomical data and a
secondary, more complex system for users who wish to reproduce results using a
specific parent dataset or limit the sources they use based on requirements
such as e.g. removing any sources with label ambiguity.

CRUMB uses three classes for its basic labels - FRI (0), FRII (1) or hybrid
(2) - and contains a total of 1006 FRIs, 997 FRIIs and 97 hybrids. All sources
labelled as “bent” by AT17 were folded into the FRI class to align with the
primary class labels of MB; if users wish to instead treat bent-tailed sources
as a distinct class, this can be done by applying logical filtering to the
complete labels (see Section 3.1). Each complete label is a vector with four
entries and provides the label for the source in each of the parent datasets
per their respective labelling schemes, in the order MB, FR-DEEP, AT17, MB-
Hyb. If a source is not present in a particular dataset, it is denoted “-1”.
As an example, a source with MB label “0” (FRI) and AT17 label “2” (bent
source) would have a complete label of “[0, -1, 2, -1]”. In this way, multiple
contradictory labels can be preserved.

###  2.5 Building the batched dataset

CRUMB was designed as a batched dataset [13], capable of loading each data
batch into memory sequentially rather than loading the entire dataset at once,
making it practical to use on machines with limited memory such as personal
laptops. It was decided to split the dataset into seven batches of 300 images
each, with one batch being reserved as a test set.

There is a significant class imbalance in CRUMB, with both the FRI and FRII
classes containing around ten times as many sources as the hybrid class.
Randomly assigning images to each batch might hence result in significant
variation in the number of hybrid sources per batch, and it was deemed that an
underpopulation of hybrids in the test set in particular might negatively
affect model evaluation. To ensure this did not occur, hybrid sources were
evenly distributed between all batches before FRIs and FRIIs were added,
resulting in a minimum of 13 hybrids in each batch. With this done, the image
data and both sets of labels were collected for each batch and used to build
the final dataset, which was then made publicly accessible on Zenodo [14].

##  3 Using CRUMB

In addition to being (to our knowledge) the largest publicly-available machine
learning dataset of FR galaxies, CRUMB was designed with the intent of being
adaptable to its users’ needs and capable of being built on further should
more suitable data be found. This section will discuss how to access various
features of CRUMB and the possibility of its expansion.

###  3.1 Accessing labels and subclasses

CRUMB is designed for use with Python machine learning packages such as
Pytorch [15] and Keras [16], and the dataset is provided with a Python class
using a structure inherited from Pytorch’s data.Dataset class.

As well as the standard targets method providing the class for each source,
CRUMB offers the complete_labels method to access the complete class label for
each of the parent surveys (see Section 2.4). Using complete_labels, it is
possible to logically filter CRUMB to include only sources which match user
criteria for either training or testing. This may be done by loading the full
dataset, identifying desired source properties (e.g. labelled as "bent" by
AT17), applying logical operations to complete_labels to identify these
sources (e.g. complete_labels[2] == 2), and passing the resulting indices to a
subloader to use only those images; this method is demonstrated on CRUMB’s
github page222https://github.com/fmporter/CRUMB.

Additionally, a number of subclasses have been created for subsets that we
expect users might find useful. These include loaders for each of the parent
datasets, combined loaders to access MB and MBHyb simultaneously or include no
sources in MB or MBHyb, and a four-class loader which treats bent sources as a
separate class from FRIs. All parent dataset loaders include a flag to default
to either the basic label or the dataset’s original labelling scheme, and all
combination loaders include a flag to specify which dataset’s labels should be
used in the case of a contradiction. Additional subclasses may be easily
constructed using these loaders as templates, making CRUMB flexible and
customisable to its users’ needs.

###  3.2 Building on CRUMB

All sources in CRUMB are fully traceable to both their parent datasets via the
complete_labels method and the J2000 coordinates of the pictured source using
the filenames method. This means that it would be possible to process
additional image data in the same manner as CRUMB (as detailed both here and
in [8]), cross-match sources and append an additional column to the complete
label accordingly to expand upon CRUMB if desired. Similarly, because of the
sources’ traceability it would be possible to produce a dataset of the same
sources using survey data at a different frequency or resolution if desired
for e.g. transfer learning applications [10]. We intend to add more sources to
CRUMB if other compatible datasets are created in the future, maintaining it
as a single dataset from which all previous versions can be obtained by
filtering the complete labels.

##  4 Conclusions

Astronomical machine learning datasets are, at present, often not published
along with the papers they were built for, making it more difficult to draw
comparisons between techniques from different publications. To improve
reproducibility of our research and offer a resource in the form of a dataset
for both astronomers who want to investigate machine learning and for computer
scientists considering using astronomical data, we have presented CRUMB, a
dataset of FR galaxies designed to combine four existing astronomical machine
learning datasets in a flexible, accessible and reproducible manner. CRUMB was
purposefully designed to be straightforward to expand as new data becomes
available, so the authors anticipate its expansion in the future to ensure it
remains a useful resources when comparing studies of FR galaxies.

## Acknowledgments and Disclosure of Funding

FMP and AMS gratefully acknowledge support from Alan Turing Institute AI
Fellowship EP/V030302/1. FMP would also like to thank Hongming Tang and
Kshitij Thorat for sharing the source lists for their datasets.

## References

  * [1] AMM Scaife.  Big telescope, big data: towards exascale with the square kilometre array.  Philosophical Transactions of the Royal Society A, 378(2166):20190060, 2020.
  * [2] Mike Walmsley, Chris Lintott, Tobias Géron, et al.  Galaxy Zoo DECaLS: Detailed visual morphology measurements from volunteers and deep learning for 314 000 galaxies.  MNRAS, 509(3):3966–3988, January 2022.
  * [3] B. L. Fanaroff and J. M. Riley.  The morphology of extragalactic radio sources of high and low luminosity.  Monthly Notices of the Royal Astronomical Society, 167(1):31P–36P, 1974.
  * [4] A. K. Aniyan and K. Thorat.  Classifying Radio Galaxies with the Convolutional Neural Network.  ApJS, 230(2):20, June 2017.
  * [5] Zhixian Ma, Haiguang Xu, Jie Zhu, et al.  A Machine Learning Based Morphological Classification of 14,245 Radio AGNs Selected from the Best-Heckman Sample.  ApJ, 240(2):34, February 2019.
  * [6] Anna M. M. Scaife and Fiona Porter.  Fanaroff-Riley classification of radio galaxies using group-equivariant convolutional neural networks.  MNRAS, 503(2):2369–2379, May 2021.
  * [7] Devina Mohan, Anna M. M. Scaife, Fiona Porter, et al.  Quantifying uncertainty in deep learning approaches to radio galaxy classification.  MNRAS, 511(3):3722–3740, April 2022.
  * [8] Fiona A. M. Porter and Anna M. M. Scaife.  MiraBest: a data set of morphologically classified radio galaxies for machine learning.  RAS Techniques and Instruments, 2(1):293–306, January 2023.
  * [9] Robert H. Becker, Richard L. White, and David J. Helfand.  The VLA’s FIRST Survey.  In D. R. Crabtree, R. J. Hanisch, and J. Barnes, editors, Astronomical Data Analysis Software and Systems III, volume 61 of Astronomical Society of the Pacific Conference Series, page 165, January 1994\.
  * [10] H. Tang, A. M. M. Scaife, and J. P. Leahy.  Transfer learning for radio galaxy classification.  MNRAS, 488(3):3358–3375, September 2019.
  * [11] T. A. McGlynn, N. E. White, and K. Scollick.  SkyView: The Digital Multi-wavelength Sky on the Internet.  In American Astronomical Society Meeting Abstracts #184, volume 184 of American Astronomical Society Meeting Abstracts, page 27.08, May 1994.
  * [12] Astropy Collaboration.  The Astropy Project: Sustaining and Growing a Community-oriented Open-source Project and the Latest Major Release (v5.0) of the Core Package.  ApJ, 935(2):167, August 2022.
  * [13] Dominic Masters and Carlo Luschi.  Revisiting small batch training for deep neural networks.  arXiv preprint arXiv:1804.07612, 2018.
  * [14] Fiona A. M. Porter.  CRUMB: the Collected Radiogalaxies Using MiraBest dataset, March 2023\.  Access via Zenodo, https://doi.org/10.5281/zenodo.7948346.
  * [15] Adam Paszke, Sam Gross, Francisco Massa, et al.  PyTorch: An Imperative Style, High-Performance Deep Learning Library.  arXiv e-prints, page arXiv:1912.01703, December 2019.
  * [16] Francois Chollet et al.  Keras, 2015.