# Analysing ecological dynamics with relational event models: the case of
biological invasions

Ruta Juozaitiene
VU Institute of Data Science and Digital Technologies, Vilnius, Lithuania;
Vytautas Magnus University, Kaunas, Lithuania
Hanno Seebens
Senckenberg Biodiversity and Climate Research Centre, Frankfurt, Germany
Guillaume Latombe
Institute of Evolutionary Biology, University of Edinburgh, UK
Franz Essl
Bioinvasions. Global Change. Macroecology Group,
Department of Botany and Biodiversity Research, University of Vienna, Austria
Ernst C. Wit111To whom correspondence should be addressed<EMAIL_ADDRESS>
Institute of Computing, Università della Svizzera italiana, Lugano,
Switzerland

###### Abstract

Aims
Spatio-temporal processes play a key role in ecology, from genes to large-
scale macroecological and biogeographical processes. Existing methods studying
such spatio-temporally structured data either simplify the dynamic structure
or the complex interactions of ecological drivers. The aim of this paper is to
present a generic method for ecological research that allows analysing spatio-
temporal patterns of biological processes at large spatial scales by including
the time-varying variables that drive these dynamics.

Location
Global analysis at the level of 272 regions.

Methods
We introduce a method called relational event modelling (REM). REM relies on
temporal interaction dynamics, that encode sequences of relational events
connecting a sender node to a recipient node at a specific point in time. We
apply REM to the spread of alien species around the globe between 1880 and
2005, following accidental or deliberate introductions into geographical
regions outside of their native range. In this context, a relational event
represents the new occurrence of an alien species given its former
distribution.

Results
The application of relational event models to the first reported invasions of
4835 established alien species outside of their native ranges from four major
taxonomic groups enables us to unravel the main drivers of the dynamics of the
spread of invasive alien species. Combining the alien species first records
data with other spatio-temporal information enables us to discover which
factors have been responsible for the spread of species across the globe.
Besides the usual drivers of species invasions, such as trade, land use and
climatic conditions, we also find evidence for species-interconnectedness in
alien species spread.

Conclusions
Relational event models offer the capacity to account for the temporal
sequences of ecological events such as biological invasions and to investigate
how relationships between these events and potential drivers change over time.

Keywords: _Ecological dynamics; First records database; Relational event
model; Species invasions._

##  1 Introduction

The number of alien species has increased continuously over recent centuries
(Seebens et al., 2017) with far-reaching consequences on the biogeography of
species, ecosystem functioning and native species (Pyšek et al., 2020). The
introduction of species beyond their native ranges is attributable to human
mediated activities that have, deliberately or accidentally, enabled them to
overcome fundamental geographical barriers. Rapid range expansion of an
invasive species is a primary threat to global biodiversity and ecosystems. A
better understanding of biological invasions may help identify potential
strategies for preventing and mitigating invasion impacts. However, the
dynamics of such inter-regional spread can be complex, varying in time, space
and among species. In addition, dynamics are usually driven by multiple
factors simultaneously, involving exogenous ecological and socio-economic
factors, such as transport hubs (Floerl et al., 2009), habitat suitability,
climatic conditions or land use (Pyšek et al., 2010; Dyer et al., 2017; Essl
et al., 2019). Given the complex interplay of natural and anthropogenic
spatio-temporal factors, it is crucial to study the effects of multiple
factors affecting alien species spread simultaneously.

Existing statistical approaches aiming to predict the occurrence and spread of
alien species, ranging from correlative (e.g. Elith and Leathwick, 2009) to
dynamic modelling tools (Dullinger et al., 2009; Chapman et al., 2016; Paini
and Yemshanov, 2012), are usually constrained to individual species and
disregard the timing of an invasion in a particular region, ignoring potential
temporal changes in the magnitude and even direction of drivers on invasions.
In particular, species distribution models (SDMs) are widely used in ecology
to estimate species geographic ranges and to understand the drivers of species
distributions. However, these models are criticized for lacking important
mechanisms (such as species dispersal, establishment, and biotic interactions)
and independent validation, among other things (Araújo and Peterson, 2012). It
is argued that the ability of SDMs to predict independent data is higher when
considering single-species versus multi-species studies. A detailed literature
review revealed that 77% of studies focusing on a single species found that
these models successfully predicted occurrence, whereas 23% of studies
incorporating more than one species found unambiguous support for the SDMs
prediction abilities (Lee-Yaw et al., 2022). Moreover, of the vast majority
previous research focused only on static SDMs, which do not account for
temporal variation of predictor and response variables (Bellard et al., 2016;
Ni and Hulme, 2021). Other approaches suggest analysing residuals of a model
with a single predictor to remove first the confounding effect, such as
sampling effort (Sofaer and Jarnevich, 2017; Bonnamour et al., 2021) or region
area (Brown et al., 2021). Several attempts to account for temporal variation
in the drivers of alien species introduction performed separate analyses for
each period of interest. For instance, multiple linear regression models were
developed to compare relationships between the number of alien species and
explanatory variables across different time periods (Dyer et al., 2017; Brown
et al., 2021). However, each subset consisted of aggregated data representing
the overall period for each region. In a recent study, global temporal
dynamics of alien species invasions were analysed using the annual data of
trade openness and the number of invasions per year (Bonnamour et al., 2021).
In this case, the analysis of relationship between these two time-varying
quantities was based on Pearson’s product moment correlation. Thus, this
approach does not allow to simultaneously quantify the effect size of multiple
drivers and ignores possible confounding factors, that are very common when
performing univariate modelling.

This paper aims to describe a generative event history model that captures the
time-varying dynamics of alien species introductions and spread by integrating
potentially relevant factors. So-called relational event model was first
developed in the field of social network analysis (Butts, 2008; Perry and
Wolfe, 2013; Juozaitienė and Wit, 2022), before being applied in multiple
other fields, for example to interactions among cattle (Patison et al., 2015),
to transfers of patients between hospitals (Vu et al., 2017) and to political
interactions (DuBois and Smyth, 2010). A relational event describes a time-
stamped interaction, such as the appearance of a species in a new location. A
relational event model thus analyses sequences of relational events following
their observed temporal order. Ecological invasion processes have crucial
relational event characteristics: an invasion is a temporal event whereby a
species becomes associated with a certain new region. The methodological
advantage of using relational event models is the fact that they are able to
describe the underlying general interaction patterns, can deal with time-
varying variables and time-varying effects, and have a suit of diagnostic
tools to test the suitability of the models on actual data.

First, we describe the data of invasion events of four taxonomic groups (i.e.
mammals, birds, insects, vascular plants) from the Alien Species First Record
Database (FirstRecords), the most exhaustive source of first records of
species in regions of the world (Seebens et al., 2017). Then, we introduce the
relational event model (REM) for studying the factors driving the spread of
alien species across large spatial (i.e. globally) and temporal (i.e.
1880-2005) scales. We analyse effects of various drivers on the timing of
first records of alien species. Based on these results, we provide a synthesis
of the advantages and caveats of using REMs for studying ecological invasions.
We conclude our work with general ideas for REM application of analysing
temporal sequences of ecological relational events in various ecological
contexts.

##  2 Methods

In this section, we first describe the data that are used in this study. Then
we give an overview of the relational event model for ecological data and
describe how it can be applied to the alien species invasion process.

###  2.1 Data

This study is based on the Alien Species First Records database (Seebens et
al., 2017), which contains the years of the first records of established alien
species across taxonomic groups in regions worldwide (i.e. countries, and sub-
national regions such as islands, archipelagos and federal states).
Established alien species refer to species introduced, accidentally or
intentionally, outside of their natural geographic range and that have
established a viable population there (Blackburn et al., 2011). The data
contains more than 47,000 records of 16,922 established alien species across
275 regions worldwide. Most data are from plants (48%), followed by insects
(26%) and birds (6%), whereas mammals are present in 3% of the records.
Regions usually correspond to countries or islands having a particularly large
number of samples. In this manuscript, we use first records of four major
land-based taxonomic groups, namely vascular plants, insects, birds and
mammals, during the period of from 1880 until 2005, which represent those
taxonomic groups and times with most available information. The reason for the
cut-off year of 2005 is data availability of first records. In previous
studies (Seebens et al., 2017, 2018), it has been found that the first records
after 2005 are distinctly affected by a lag in reporting. It takes time until
a species has been recognized in a new environment, assessed as being an alien
species with established (i.e., reproducing) populations and the new record
has been published. The estimate from previous studies is that this process
takes on average 10-15 years, which is the reason why more recent first
records cannot be reliably interpreted at a global scale. The selected
taxonomic groups cover 81% of the data contained in the Alien Species First
Record Database.

We consider a range of explanatory variables to characterise the spatial
locations, climatic conditions, political relationships such as former
colonial ties among countries and various socio-economic variables during that
period (see Table 1). These variables have been shown to be of importance for
explaining the current distribution of alien species (e.g. Pyšek et al., 2010;
Dawson et al., 2017; Essl et al., 2019).

Land cover types have been shown to be associated with different levels of
alien species presence (Chytrỳ et al., 2009; Polce et al., 2011; Sánchez-Ortiz
et al., 2020). Here, information on land cover is used to characterise the
potential suitability of regions for the establishment of alien species and
the degree of anthropogenic disturbance, which facilitates establishments.
Land cover is described by the proportion of cropland, pasture and urban area;
data were obtained from the Harmonized Global Land Use data set, V1 (Hurtt et
al., 2011). The grid-based information has been averaged for the regions and
years (1971-2012) considered in this study. As land cover data were provided
on a decadal basis, while the analysis was run on an annual scale, land cover
variables were interpolated from a 10-year resolution to obtain an annual
resolution.

It is expected that favourable climate conditions facilitate the spread of
alien species and unsuitable environments with stressful conditions often lead
to a lower establishment success (Finch et al., 2021). To account for habitat
suitability we use the air temperature at 2 meters above the ground as the
near-surface air temperature. This variable is used to represent climate
conditions, using simulated annual values covering the period from 1871 to
2012 (Watanabe et al., 2011). Annual near-surface air temperature values were
averaged for each region and year between 1871 and 2012.

International trade has been identified as a major driver of alien species
introductions (via increased propagule and colonization pressures), and
therefore increases the probability of species to become established (Westphal
et al., 2008; Chapman et al., 2017; Blackburn et al., 2020; Bagnara et al.,
2022). International trade values were extracted from the Correlates of War
project (Barbieri et al., 2009). The dataset includes information on trade
flows between countries from 1871 to 2009. The dyadic trade dataset values
represent imports from one country to another in current U.S. dollars. The
dataset has missing values because trade flow is unknown for some countries
within certain time periods and some countries are not included in the
dataset. To address this issue, missing values that appear at the beginning of
the observational period were replaced with zeros in those time series that
start at zero value. In other cases, i.e., for time series for which the first
recorded trade value between countries $i$ and $j$ is not equal to zero, we
extrapolated missing values using the log-linear model,
$\log(\mbox{trade}(i,j,t)+1)=\alpha_{ij}+\beta_{ij}t$. This model assumes a
linear relationship between the logarithm of trade volume between countries
$i$ and $j$ at time $t$ and the time point $t$. In principle, any
extrapolation model could be used for this purpose. However, given the nature
of economic growth, an exponential growth model is appropriate and therefore
was used to determine the missing values.

A substantial part of the most important trade routes has formed since
colonial times. It is expected that former colonial ties facilitated the
spread and establishment of alien species (Yang et al., 2021; Lenzner et al.,
2022). Colonial ties were therefore included in this analysis using the
Colonial Dates Dataset (COLDAT) (Becker, 2019), which provides information on
the reach and duration of European colonial empires. Each country is uniquely
characterised by the colonial empire it belonged to or classified as
independent, in case there were no colonial ties.

Previous studies revealed that observed invasion dynamics are highly
correlated with variations in sampling effort and, thus, it is important to
address this bias before reaching conclusions regarding temporal trends of
biological invasions (Bonnamour et al., 2021; Sofaer and Jarnevich, 2017). To
account for variations in sampling effort, we used the number of species
present in a particular region in 1880 as a proxy. The number of present
species was calculated by summing the number of First Records up until 1880
and the number of native species. This variable was used as an approximate
measure of sampling effort under the assumption that there was no saturation
of the pool of native species yet to be recorded in each country at that time
and that sampling effort is constant over time. To acquire information on
native ranges for all species considered in our analyses, we combined data
from various sources. For mammals and birds, native ranges were extracted as
range maps from the IUCN range maps (https://www.iucnredlist.org, accessed
08.07.2016) and BirdLife (http://datazone.birdlife.org/species/requestdis,
accessed 01.11.2016), respectively. A species was then assigned to a region
listed in the first record database if the range map intersected with the
region. For vascular plants and insects, native ranges were extracted from
(van Kleunen et al., 2018) and the CABI Invasive Species Compendium
(https://www.cabi.org, accessed 15.07.2016), respectively. Region names were
harmonised to match the regions provided in the first record database.

To confirm independence between explanatory variables, we used a correlation
heatmap indicating correlations among all variables (Figure 2). Coefficient
values were small with the highest observed correlation equal to 0.42.
Therefore we conclude that multicollinearity does not cause a problem in the
fitting of the model.

For each species, the final data set used for analysis contained the native
and alien ranges with the respective first records. The final dataset included
16403 invasion events recorded for 4835 species (615 birds, 186 mammals, 3920
plants, and 114 insect species). The average number of invasion events per
taxonomic group was $\sim 5$ records per insect species, $\sim 4$ per bird and
mammal species, and $\sim 3$ recorded events per plant species. On average,
$\sim 46$ invasion events were recorded per region, ranging from 1 to 1685
events (Figure 3). 157 regions had less than 15 observations, with 33% of
these regions located in Africa; 36 regions had more than 100 records. The
region with the highest number of invasion events was Australia. Figure 4
illustrates the spread patterns of four species from different taxonomic
groups: the muskrat (mammal), the common ragweed (plant), the cassava mealybug
(insect) and the monk parakeet (bird). These plots give important indications
regarding the existence of natural environmental barriers affecting the spread
and establishment of alien species. For instance, the cassava mealybug is
native to South America, but it started spreading throughout sub-Saharan
Africa since its inadvertent introduction into the continent in the early
1970s (Parsa et al., 2012). The data supporting the findings of this study is
publicly available at https://doi.org/10.6084/m9.figshare.20017691.v1. We also
provide a full source code which is available at
https://github.com/Juozaitiene/REM-for-species-invasion.

###  2.2 Relational event model for ecological data

A relational event model (REM) is a probabilistic model that describes the
temporal interactions between a (possibly dynamic) set of sender nodes $S(t)$
and a set of receiver nodes $R(t)$ (Butts, 2008). Sender and receiver nodes
can either represent elements of the same group, so-called _one-mode networks_
, e.g., grooming relationships between a group of primates, or between
elements of different groups, so-called _bipartite networks_ , e.g., predation
relationships between predators and prey. In contrast to traditional network
analysis, the relational event model explicitly considers these networks to be
embedded in time. There are two main ways in which ecological relationships
can be embedded in time. First, the duration of an ecological relationship can
be an important and interesting aspect of the relationship itself. For
example, in certain animal species, there are sequential periods of monogamy,
whose durations may vary. Second, the temporal positioning of the interaction,
rather than its duration, can be the primary focus of interest, as for
grooming events. The relational event model is intended for modelling
interactions of the second type.

In this section we introduce the relational event model. First, we focus on a
generic model formulation for any type of ecological data. The aim is to
explain the idea behind the model and how it can be estimated from the data.
Thereafter, we describe a specific implementation of the relational event
model designed to model species invasions.

####  2.2.1 Generic model formulation

Relational event models are so-called _generative probabilistic models_. This
means that they aim to describe the underlying processes that give rise to the
observed interaction events between sender and receiver nodes. This is to
distinguish these models from _purely phenomenological, statistical models_ ,
which aim to describe the correlations between some observed variables, such
as the number of grooming events involving a particular primate in a
particular period and its age. In fact, rather than modelling some aggregate
quantity, relational event models aim to model each individual event
separately. In what follows, we aim to give a description of relational event
models, as well as the techniques that are used to estimate them from actual
data. An example of all concepts introduced is shown in Figure 1.

A relational event model is considered on a certain time interval
$[t_{b},t_{e}]$, where $t_{b}$ is the beginning and $t_{e}$ the end of the
observation period of the ecological process. During this period interaction
events occur between the group of sender nodes and the group of receiver nodes
at time instances $t\in[t_{e},t_{b}]$. The set of available senders and
receivers can change over time and does not have to be static. This is an
important advantage over other network models that typically require the
interacting elements to stay the same over the observation period.

##### Risk set.

As an example, let’s consider a water hole, where groups of prey and predators
regularly come to drink. At the beginning of the observation process
$S(t_{b})$ represents the set of predators and $R(t_{b})$ the group of prey.
The predators sometimes attack and kill some prey, while they are gathered to
drink. The first attack can be represented as a particular predator $s\in
S(t_{b})$ that kills a particular prey $r\in R(t_{b})$ at a certain time
$T_{sr}>t_{b}$. After that attack, i.e., $t>T_{sr}$, the prey set diminishes,
i.e., $R(t)=R(t_{b})\backslash\\{r\\}$, whereas the predator set stays the
same, $S(t)=S(t_{b})$. However, if a new predator $s_{\mbox{\scriptsize new}}$
joins the group at time $t$, then $S(t)=S(t_{b})\cup\\{s_{\mbox{\scriptsize
new}}\\}$

At any time $t$, the possible interactions $(s,r)$ that can occur is an
element of the so-called _risk set_

$\mathcal{R}(t)\subset S(t)\times R(t).$

For each element $(s,r)\in\mathcal{R}(t)$, there exists, in principle, the
event time $T_{sr}>t$, defined as the first interaction time after $t$ between
$s$ and $r$. However, only the first of these event times

$T={\arg\min}_{(s,r)\in\mathcal{R}(t)}T_{sr}$

is the one that is actually observed. After that, as we have seen above, the
risk set may change and certain event times may be censored or new ones may be
introduced.

##### History.

There is another important way that the situation potentially changes after an
event occurs. If some predator $s$ has killed prey $r$ at time $T=t$, then we
can hypothesize the system to change in two ways: possibly the predator will
be satisfied for a while, and the time $T_{sr^{\prime}}$ that he will attack
another prey $r^{\prime}\neq r$ will be delayed; moreover, the other prey may
be spooked and more at their guard, thereby delaying the attack times
$T_{s^{\prime}r^{\prime}}$ even from other predators $s^{\prime}\neq s$. To
describe the collection of previous events, we introduce the concept of the
_history of the process_ ,

$\mathcal{H}_{t}=\\{\mbox{all time-stamped interaction events
$(s_{i},r_{i},t_{i})$ by time $t$, i.e., $t_{i}\leq t$}\\}$

Related to this is the set $\mathcal{H}_{t-}$, which stands for the history of
all events strictly before time $t$.

##### Rate function.

With the above definitions of risk set and history in place, it is now
possible to define a generic and generative model for the interaction event
sequences. Although this is formally done through the definition of an
underlying counting process (Perry and Wolfe, 2013), we use a more intuitive
description via the sequence of events
$\\{(s_{i},r_{i},t_{i})~{}|~{}i=1,\ldots,n\\}$. The dependence of event times
on the constantly changing risk set and history can be expressed by assuming
that each event time $T_{sr}$ has a generalized exponential distribution with
(non-constant) _rate function_ $\lambda_{sr}(t|\mathcal{H}_{t-})$. The rate
function should be chosen in such a way as to express the scientific
hypotheses about the ecological interaction process that the ecologist would
like to test. Given that the rate function is strictly positive, it is common
to write the function in an exponential form (Cox, 1972),

$\lambda_{sr}(t|\mathcal{H}_{t-})=\lambda_{0}(t)e^{x_{sr}(t)^{\prime}\beta},$
(1)

where $x_{sr}(t)\in\mathbb{R}^{p}$ are $p$ known _dyadic_ covariates at time
$t$ that are allowed to depend on $\mathcal{H}_{t-}$, the history of the
process right before time $t$. The parameter $\beta\in\mathbb{R}^{p}$ express
the effect size of each of the above covariates. The aim of an empirical
ecological study is typically to use the data to test the significance, sign
and size of the effect sizes of $\beta$ in order to test one or more
ecological hypotheses of interest. The _baseline hazard_ is an unknown
function of time, which aims to capture residual effects not strictly related
to the sender and receivers.

For example, in the predator-prey example, we hypothesized that the rate of
the kill process might depend on the appetite of the predator and the
alertness of the prey. If we define $x_{s}(t)$ as the last time before $t$
that predator $s$ killed a prey, then we could hypothesize the following rate
function,

$\lambda_{sr}(t|\mathcal{H}_{t-})=\lambda_{0}(t)e^{\beta_{1}\left(t-x_{s}(t)\right)+\beta_{2}\left(1-\max_{s^{\prime}}x_{s^{\prime}}(t)\right)},$

where $\beta_{1}$ describes the saturation process of the predator and
$\beta_{2}$ the alertness of the prey. Testing whether $\beta_{1}=0$ or
$\beta_{2}=0$ would allow the scientist to decide whether any of these two
ecological hypotheses are real. The baseline hazard $\lambda_{0}(t)$ may
describe a diurnal cycle, capturing the fact that possibly predators are more
active in certain times of the day than others.

##### Effect estimation from data.

Using the standard properties of the exponential distribution, it can be shown
that each inter-arrival time $t_{i}-t_{i-1}$ has an exponential density with
rate
$\int_{t_{i-1}}^{t_{i}}\sum_{(s,r)\in\mathcal{R}(t)}\lambda_{sr}(t|\mathcal{H}_{t-})~{}dt$,
and, more importantly, that the conditional probability of event
$(s_{i},r_{i})$ at time $t_{i}$ is given by
$\frac{\lambda_{s_{i}r_{i}}(t_{i}|\mathcal{H}_{t_{i}-})}{\sum_{(s,r)\in\mathcal{R}(t_{i})}\lambda_{sr}(t_{i}|\mathcal{H}_{t_{i}-})}$.
This ratio is the result of the relative size of the hazard rate of actual
event $(s_{i},r_{i})$, in the numerator, and the hazards rates from all other
possible events $(s,r)$ that could happen at time $t_{i}$, in the denominator.
Using the exponential hazard formulation from (1), the baseline hazard
$\lambda_{0}(t_{i})$ cancels, resulting in

$P\left((s_{i},r_{i})\mbox{ happens
}|~{}t_{i}\right)=\frac{e^{x_{s_{i}r_{i}}(t_{i})^{\prime}\beta}}{\sum_{(s,r)\in\mathcal{R}(t_{i})}e^{x_{sr}(t_{i})^{\prime}\beta}}.$

The effect $\beta$ can be estimated by maximizing the probability of all the
occurred events, i.e.,

$\hat{\beta}={\arg\max}_{\beta}\prod_{i=1}^{n}\frac{e^{x_{s_{i}r_{i}}(t_{i})^{\prime}\beta}}{\sum_{(s,r)\in\mathcal{R}(t_{i})}e^{x_{sr}(t_{i})^{\prime}\beta}}.$

This $\hat{\beta}$ is typically referred to as the partial likelihood
estimate, which possesses many desirable properties (Cox, 1975). The baseline
hazard can be estimated posthoc according to Breslow (1975).

####  2.2.2 Relational event model for first invasion events

In this section, we will define the relational event model for first invasion
events of alien species into new regions. Here, the sender modes are species
and the receiver nodes are regions, and an invasion event describes the first
record of an alien species in a region as provided in the FirstRecords
database. In principle, it is also possible to consider senders and receivers
as geographical units, but as information about the region, where the species
came from (i.e., the sender), is lacking, we could not consider this approach.
The main aim of applying REM to biological invasions was to disentangle the
generative process of the “first invasion” relationship between species and
regions over time. We analyse four groups of species,
$L_{i}=\\{s_{1},\ldots,s_{n}\\}$ ($i$=mammals, birds, plants, insects),
separately across a fixed time-frame $[1880,2005]$, in which the invasion
process of the species is recorded relative to geographic set
$C=\\{c_{1},\ldots,c_{272}\\}$ of the 272 pre-specified regions. If at time
$t$ a species emerged in a new region $c$, then this invasion event can be
described by a triple $e=(s,c,t)$, whereby species $s$ spread to a region $c$
at time $t$. The relational event model describes a general framework for
modelling and inferring the stochastic species invasion process
$\\{(s,c,t)\\}_{t\in[1880,2005]}$ in order to identify the reasons why
particular species invade particular regions and why they did so at the time
they did. Although we study many species of a given taxonomic group together,
this does not mean that the relational event model assumes that the species
are physically interacting. In fact, the joint mixed effect analysis is able
to identify global drivers that are common to the spread of species as well as
capture various aspects of individual species dynamics.

We focus on several possible drivers of invasions that can be distinguished
into four groups (see Table 2). The first two groups encompass two types of
dynamic processes. The first group consists of time-varying factors (i.e.
drivers whose values change over time) with constant effects over time. This
is expressed by constant coefficients $\beta$ in the hazard function. This
group comprises temperature, colonial ties and the number of prior invasions.
The second group focuses on, again time-varying, factors whose effects change
over time. For example, not only can the volume of international trade vary
over time, but the effect of trade volume on the rate of invasion events may
also vary over time. This is expressed by time-varying coefficients in the
REM. This group encompasses distance, land cover and trade. Whether an effect
should belong to the first or the second group depends on whether or not the
evidence in the data suggests that a constant effect is sufficient. The third
group covers random effects of a monadic type, which are variables measured at
the level of a region or species (species invasiveness and region
invasibility), whereas the final group of factors considers dyadic random
effects, which denote variables measured at the level of pairs of species,
capturing the inter-dependence in alien species spread.

For each taxonomic group separately (plants, birds, insects, mammals), we fit
a relational event model representing global invasion dynamics over the period
from 1880 to 2005. The overall time-dependent hazard function of species
invasion dynamics considered for each taxonomic group is given as,

$\displaystyle\log\lambda_{sc}(t)$ $\displaystyle=$
$\displaystyle\log\lambda_{0}(t)+d_{sc}(t)\beta_{1}(t)+tr_{sc}(t)\beta_{2}(t)+dt_{sc}(t)\beta_{3}+l_{c}(t)\beta_{4}(t)$
(2)
$\displaystyle+u_{c}(t)\beta_{5}+k_{sc}(t)\beta_{6}+pi_{c}(t)\beta_{7}+se_{c}\beta_{8}+b_{s}+b_{c}+b_{ss_{c}(t)},$

where

  * •

The baseline hazard $\lambda_{0}(t)$ captures the underlying changing
intensity of the overall invasion process.

  * •

$d_{sc}(t)$ denotes the distance to the nearest region $c$ in which species
$s$ is present by time $t$. The proposed definition accounts for the distances
between invaded countries as well as native ranges. Long-distance invasion
events through natural dispersal are typically rare. Therefore, the distance
between two regions that species have to overcome could be an important factor
responsible for alien species invasions. The distance between two regions is
defined as the distance between their closest borders, and thus the distance
between neighbouring regions is zero. The proposed distance covariate is time-
varying, as the minimum distance between a region and another invaded region
will change over time as new regions get invaded. An example is shown in
Figure 1, which also explains the estimation procedure.

  * •

The time-varying quantity $tr_{sc}(t)$ is the logarithm of the sum of annual
trade flows between region $c$ and other regions invaded by species $s$ by
time $t$ – the log-transformation is used to capture saturation.

  * •

The variable $dt_{sc}(t)$ represents the minimum absolute difference in
temperature between region $c$ and other regions invaded by species $s$ by
time $t$. It is a time-dependent variable since the list of regions invaded by
species $s$ is changing over time.

  * •

$l_{c}(t)$ is the time-varying variable that describes the proportion of
agricultural land, defined as the sum of cropland and pasture proportions in
the region $c$ at time $t$.

  * •

$u_{c}(t)$ is the proportion of urban area in region $c$ at time $t$.

  * •

$k_{sc}(t)$ indicates the presence of species $s$ at time $t$ in a region
within the area of a colonial power.

  * •

The variable $pi_{c}(t)$ is a weighted number of prior invasions in region
$c$. We use a standard inverse exponential downweighing of recent invasions to
capture the current invasibility of certain regions for alien species
invasions, i.e, $pi_{c}(t)=\sum_{s\in
L_{i}}1_{\\{t_{sc}<t\\}}{0.95}^{t-t_{sc}}$, where $t_{sc}$ is the moment that
species $s$ invaded region $c$. This variable summarizes current invasibility
by giving more weight to more recent events.

  * •

The variable $se_{c}$ aims to correct for confounding effects by capturing the
differential sampling effort within each country. It is defined as the number
of species that were recorded in each country by 1880.

  * •

The relational event model also includes random effects representing species
invasiveness ($b_{s}$) and region invasibility ($b_{c}$). It is expected that
different species may vary in their spread behaviour, whereas similarly
certain regions may be more attractive destinations than others.

  * •

The last group of model explanatory variables considers dyadic random effects
that aim to capture interdependence in alien species spread. The effect
$b_{ss^{\prime}}$ describes the affinity of species $s$ and $s^{\prime}$ in
their spread patterns. The species indicator $s_{c}(t)$ is defined as the last
species to invade $c$ before $t$. We consider both a general formulation in
which the order of invasion by the species is accounted for, as well as the
more parsimonious symmetric formulation that constrains
$b_{ss^{\prime}}=b_{s^{\prime}s}$. It turns out that although the symmetric
formulation is often sufficient, for some taxonomic groups the general non-
symmetric random effect offers additional explanatory power. However, some
species combinations are very rare or never appear in the data. To robustly
estimate the effect of species interdependence in spread, which is not
affected by low invasion counts, we focus on species with the highest number
of recorded invasion events. Thus, this specific effect is estimated using
only the information on the top thirty most widespread alien species in each
taxonomic group.

Time-varying effects were modelled by estimating a piecewise constant
coefficient for each of five equally sized periods (1880-1905, 1906-1930,
1931-1955, 1956-1980, and 1981-2005). Analyses were conducted using the
coxme() function from the coxme package version 2.2-16 within the R
statistical software to fit a Cox proportional hazard model with mixed effects
(R Core Team, 2022; Therneau, 2020). Chi–square tests were used to compare
different models. To assess whether a fitted Cox regression models adequately
described the data we performed model diagnostics. The proportional hazards
assumption was checked using statistical tests and graphical diagnostics based
on the scaled Schoenfeld residuals. Thus, we used the cox.zph() function from
the survival package version 3.2-13 (Therneau, 2022) to test for independence
between residuals and time, the ggcoxzph() function from the survminer package
version 0.4.9 (Kassambara et al., 2021) to perform a graphical diagnostics and
the ggcoxdiagnostics() function from the survminer package version 0.4.9 to
test influential observations or outliers. The diagnostic tests suggested an
overall adequate fit of the model, for more details see Supplementary
Materials (Section S2). A more detailed description of the inference and model
diagnostics is provided in the online supplementary materials.

##  3 Results

In this section, we present the results of the alien species invasion analysis
using the relational event model. The final REMs provide a parsimonious
description of the data as summarized in Table 4. The models explained between
12.39 and 26.1% of the variance, depending on the taxonomic groups as shown in
Table 4 (all p-values $<10^{-4}$). Note that the provided $R^{2}$ values
should not be confused with the $R^{2}$ for linear regression. In event
history models, such as the REM, the $R^{2}$ coefficient has a fundamental
upper limit of roughly 56%, if there is no additional over-dispersion and all
the fundamental causes and their effect sizes were known perfectly (Zheng and
Agresti, 2000). In this light, the reported $R^{2}$ values suggest that the
models provide significant insight in the underlying generative process. To
confirm the hypothesis that the magnitude of the effect of some factors may
change over time, we tested for the benefit of including time-varying effects.
Thus, we analysed an extensive list of different models, where the smallest
model included only static coefficients and the largest included all time-
varying coefficients. To select the best model for each taxonomic group we
performed Wald tests for simple effects and chi-squared tests for time-varying
effects. Goodness-of-fit tests presented in Table 3 show that the proposed
models with several time-varying coefficients are a significant improvement
over the model including only static coefficients as well as the model
including only time-varying coefficients. More detailed results are provided
in the Supplementary Materials (Tables S1-S2). The performed analysis reveals
that models that include distance, land cover and trade are highly
significant, with the sole exception of cropland-pasture effects for insects.

##### Trade flow.

It is well-established that trade is crucial for alien species spread across
the globe (Seebens et al., 2015; Dawson et al., 2017; Essl et al., 2011).
However, this picture seems to be more complicated for invasion events than
might be expected. Although trade was certainly of major importance in the
early 20th century, the effect seems to have largely disappeared more
recently. During the first two periods (1880-1930), the effect of trade was
positive, i.e., the spread rate increased with the amount of trade between a
region (Figure 5). However, the effect has become negative since. This finding
may seem counter-intuitive, but can partially be explained by the fact that
the number of introduced species did not keep pace with the increase in trade
volumes due to a saturating relationship between imported commodities and the
number of alien species (Seebens et al., 2017). As a consequence, the relative
importance of trade as a predictor for alien species invasions declined over
time. Further, the proportion of goods with known high relevance for
introducing alien species such as agricultural and silvicultural goods (e.g.
wood, grain) has declined severely in the last decades (Luppold, 1988). There
has also been a structural shift in commodity composition, e.g., bulk
commodities have been substituted by processed consumer-ready products. The
share of bulk commodities declined from 50% in 1980 to 32% in 1995, while the
share of non-bulk products increased from 50% to 68%. Moreover, the share of
products with low risks of introducing alien species (e.g. electronics, goods
from other technology-intensive industries, and labour-intensive products,
particularly clothing) have grown most rapidly in world exports during the
period 1980-2000 (Mayer et al., 2003). A broader analysis of the absolute
values of traded goods would give a deeper understanding of the development of
international trade flows. However, it is very difficult to establish to what
extent the absolute trade (in units of weight such as tons) of raw goods has
changed over the past century. Even the value of bulk trade is not a good
proxy and, even then, is only known for a few instances.

##### Distance.

Distance has a negative effect on the rate of invasion events, which means
that long-distance invasion events are typically rare. Moreover, for most
taxonomic groups the effect size grew over time, i.e., the importance of
distance on species invasions increased in recent years (Figure 6a). We argue
that this can possibly be explained by the development of national and
international legal instruments focusing on reducing the risk of introducing
alien species through transportation (Turbelin et al., 2017). The effect of
distance has become especially strong after 1980 for mammals and plants, which
coincides with the fact that since the 1990s there has been a significant
increase in legislation regarding invasive alien species, especially in
countries with high rates of species invasion. The opposite trend has been
inferred for insects, i.e., the effect of distance decreased over time but has
remained negative.

##### Land cover variables.

The effect of the proportion of urban area in region is only significant for
vascular plants. This effect has been consistent over the whole time span we
analysed (Table 5). Increases in the fraction of urban areas in a region are
negatively correlated with spread rate for plants, as shown by the negative
coefficients. More specifically, for every increase in urbanization by 1%, the
rate of invasion events in that region decreases by 16.8%. On the other hand,
the effect of cropland and pasture on the rate of invasion events is more
complex and has changed over time in a non-linear fashion. It is significant
for plants, mammals and birds and its time-varying coefficients is negative
and has been growing for most of the period under consideration. Recently,
this effect for birds seems to be decreasing again (Figure 6b). This means
that larger fractions of cropland and pasture in a region resulted in lower
spread rates. Modern agricultural techniques may be responsible for these
particular changes.

##### Sampling effort.

The proxy for sampling effort has a significant effect on the rate of species
invasions for plants, insects and mammals, confirming previous studies by
Bonnamour et al. (2021). For our purposes, sampling effort has shown to be an
important confounding effect that should be accounted for in the model, in
order to correctly interpret other structural effects.

##### The number of alien species found in a particular region.

Both plants and insects are significantly affected by the number of species
that have already invaded a certain region. Both plants and insects were
recorded less frequently in regions that had already been invaded by many
other species. The effect for plants was nonetheless very small but
statistically significant, i.e., decreasing the hazard by 0.5% per invasion.
For insects the observed effect was higher previous invasions decreased the
hazard by 10.4%.

##### Thermal similarity.

For all taxonomic groups, rates of invasion events are significantly affected
by the climatic conditions of the receiving region. The temperature
differential effect for all taxonomic groups is negative (Table 5). The
parameter values associated with $dt_{sc}(t)$ mean that for each degree of
mean annual temperature difference to the nearest invaded region, the spread
rate has decreased by 15.8%, 8.5%, 6.2% and 5.1% for plants, insects, mammals
and birds, respectively, keeping everything else constant. The negative effect
of temperature for all taxonomic groups shows that biological invasions are
strongly affected by climatic conditions. This finding supports the importance
of climate matching between regions for alien species spread and establishment
(Abellán et al., 2017; Lovell et al., 2021).

##### Colonial ties.

The effect of colonial ties is not statistically significant for any taxonomic
group.

##### Species invasiveness.

Analysing the random effects, we can identify the species with the highest
invasiveness propensity (i.e. rate of spread across regions after accounting
for the influences of covariates). After accounting for the other effects
described above, the random invasiveness rate captures the tendency of
different species to spread. Table S3 in Supplementary Material presents the
top five alien species in each taxonomic group. Siberian Chipmunk (Tamias
sibiricus), Raccoon dogs (Nyctereutes procyonoides), and American mink
(Neovison vison) were the mammals with the highest invasiveness rate. These
results are in line with previous studies (Tedeschi et al., 2021b, a) stating
that these species are the most widespread in Europe, having invaded at least
27 countries each. The prevalence of these species may be related to the fact
that many game species and fur-bearing animals have been moved, introduced, or
have even escaped from farms into the wild.

For plants, water plants (Hydrocotyle ranunculoides) and conifers such as
Caribbean pine (Pinus caribaea) and Douglas fir (Pseudotsuga menziesii) ranked
first. It is argued that conifer species escaping from plantations and
becoming invasive are often those that have been cultivated the most widely
and for the longest time (Křivánek et al., 2006), while small seed mass, short
juvenile period, and short intervals between large seed crops have been shown
to contribute to the invasiveness of conifer trees (Richardson and Rejmánek,
2004).

The top three birds species with the highest invasiveness rate were the Rock
pigeon (Columba livia), followed by the Alexandrine parakeet (Psittacula
eupatria) and the Egyptian goose (Alopochen aegyptiaca). Indeed, Rock pigeons
have an extensive distribution in Eurasia and North Africa (Goodwin, 1960),
whereas Alexandrine parakeets and Egyptian geese are considered among the most
widespread alien species in Europe (Gyimesi and Lensink, 2012; Şahin and
Arslangündoğdu, 2019).

The highest ranking insects were pests, i.e. a fruit fly (Bactrocera invadens)
which was first introduced in Kenya and has rapidly expanded in 28 African
countries (Khamis et al., 2012); pink hibiscus mealybug (Maconellicoccus
hirsutus), which is considered a highly invasive species that has been
introduced to 75 countries in all over the world (Milonas and Partsinevelos,
2017), and western flower thrips (Frankliniella occidentalis), which is
indicated as a major worldwide crop pest (Kirk and Terry, 2003).

##### Region invasibility.

Figure 7 represents the random region effects reflecting region invasibility
for birds. Plots for other taxonomic groups are provided in the supplementary
materials, Fig. S15. Several countries that were part of the former British
Empire (New Zealand, South Africa, United Kingdom, USA) are on top of the list
of the country invasion propensity list for all four taxonomic groups. This
finding is in line with previous findings on hot spots of global invasion
levels across eight taxonomic groups, where these countries all show up
prominently (Dawson et al., 2017). There are likely several reasons that
jointly explain this result. Interestingly, most of these countries are large
and have a long history of an advanced ecological-scientific infrastructure
that has kept accurate records of biological data including species invasions.
Further, these countries have been part of one of the largest European Empires
and the integration into a near-global political entity and higher levels of
trade and people movements among these regions may have fostered biological
invasions as well (Lenzner et al., 2022). On the other hand, several of these
countries such as New Zealand and Australia have implemented strict
biosecurity policies to stop the introduction of alien species (Turbelin et
al., 2017). However, these policies are mostly of rather recent times and may
be masked by the time frame (1880-2005) under observation here.

##### Interdependence of alien species spread.

Standard deviations of the random effects show that the symmetric formulation
of the interdependence in spread propensity of pairs of species provides
similar explanatory power as the non-symmetric formulation in most cases
(Table 5). Only for mammals did the non-symmetric random effect offer
additional explanatory power. In this case, the temporally symmetric effect
for species interactions between mammals is not significant, whereas the
temporally asymmetric effect for mammals is highly significant. This clearly
shows that temporal ordering is an important aspect in explaining species
invasions. Models that would not distinguish between temporal order would not
be able to identify this effect. The estimated random effects for insects are
represented in Figure 8, whereas the others are included in the Supplementary
Materials.

##### Baseline hazard.

The baseline hazard $\lambda_{0}(t)$ is a remainder term of the relational
event model. It represents the underlying rate of spread over time that is
unexplained by the drivers included in the model. Figure 9 shows the
cumulative baseline hazard estimates for the four taxonomic groups covered.
The linear shape of the baseline hazard for both plants and birds suggests
that they have an approximately constant baseline hazard
$\lambda_{0}(t)=\frac{d\Lambda_{0}(t)}{dt}$, indicating that the spread rate
as explained by the model was mostly constant for these species. For mammals,
the cumulative baseline hazard has flattened since the 1960s, suggesting that
their spread rate has been decreasing. Earlier, the spread rate of mammals
outpaced that of the other taxonomic groups by a factor of two to three.
Interestingly, the two big waves of the mammal curve correspond to the two
World Wars, suggesting a real variation in species invasions as a result of
the World War (less international trade and travel, except military vehicles),
and it is related to a recording artifact since the number of available first
records likely decreased because of lower sampling. Instead, we observe a
rather constant hazard for birds and plants, suggesting no general increase in
the spread rate. However, the picture for insects is quite different, i.e.,
insects have seen a marked increase in their baseline hazard since the 1980s.
This implies that the spread rate has been accelerating for insects since then
in a way that cannot be explained by the current variables in the model, which
is in line with previous findings by (Bonnamour et al., 2021). These results
partly suggest that international biosecurity policies and legislation have
varying levels of success for different taxonomic groups. We hypothesize that
preventing and controlling early-stage mammal and bird invasions is much
easier than stopping insect invasions. Most species of mammals and birds have
been intentionally introduced for sport hunting, commercialization and
domestication as livestock, pets, or for pest control (Carpio et al., 2020).
Some other species have been semi-intentionally introduced (intentional import
but accidental escape), while most insect invasions are accidental. This is
because most insects are small and often go unnoticed until they have become
established (Eschen et al., 2015).

##  4 Discussion

In this section, we will first discuss the complexities and limitations of the
application of the relational event model to the alien species invasion
process. Furthermore, we discuss the opportunities as well as the difficulties
of applying this model more broadly in ecology.

Relational event models offer the capacity to account for the dynamic nature
of ecological relational events, such as biological invasions, and to
investigate how relationships between these events and potential drivers
change over time. These models have in principle the ability to identify
causal drivers of the relational process by accounting for confounding effects
that almost certainly are present in any complex ecological scenario. However,
REMs require comprehensive and accurate temporal data on potential drivers and
events, which are often scarce or difficult to come by. Our analyses should be
seen against the backdrop of possible recording delays of the first
introduction of a species in a region. Other sources of bias include uneven
sampling over space and time as well as species-specific detectability issues.
The detectability may be influenced by various factors, such as species
richness, body size, and environmental factors, such as vegetation structure
or successional change (Isaac et al., 2011). The uneven sampling of alien
species across regions mirrors the uneven distribution of human population and
socio-economic activity. Some regions have a long history of an advanced
ecological-scientific infrastructure where recording of biological data
including species invasions is much better than in other regions. These
uncertainties might explain the comparatively low degree of explained
variation. Information on trade, land use, and economic variables was not
easily available, and therefore proxies had to be employed.

Moreover, alien species spread is a much more complex process than the one
described in the presented REM, which is based on several simplifying
assumptions. For instance, distance can affect post-introduction natural
dispersal and also human-mediated dispersal (Wilson et al., 2009). In the case
of natural dispersal, the effect of distance on alien species spread is
species-specific (because of dispersal abilities) and, in the case of human-
mediated dispersal, the effect of distance on spread depends on pathways
(e.g., intentional or accidental), vectors (e.g., car or boat, or plane) and
also species traits (e.g., attaching capacities, egg laying behaviour,
preferred substrate) (Meurisse et al., 2019). Some species might be poor
dispersers naturally but easily transported through human-mediated dispersal
and other species can spread widely by their own following the initial
introduction (e.g., Vespa velutina) (Verdasca et al., 2021). This means that
although our relational event models give a more subtle picture of the time-
dependent invasion process of species, the parameter estimates have to be
interpreted in the light of these limitations.

This paper introduced a flexible and powerful REM modelling approach for the
analysis of ecological data and demonstrated its practicality and
effectiveness by applying the approach to the process of alien species
invasion. These models are highly generic and merely require information about
the temporal sequence of interactions. Analysis software is available in any
statistical package, which means that there are very few practical obstacles
for using these models more widely in ecology.

The REM can be generalized to any ecological interaction process where a
number of actors interact together over time. Biological interactions, which
we more generally refer to as ’relational events’, are the fundamental process
shaping the diversity of life (Del-Claro and Torezan-Silingardi, 2021a). These
relational events capture important basic relationships in ecology (Del-Claro
and Torezan-Silingardi, 2021b). They can be intraspecific, i.e., such as male-
male interference competition (Dijkstra and Border, 2018) or mating,
interspecific, such as predation or pollination, or form some more generalized
interaction, e.g., between a species and its habitat or the environment. Often
it is possible to associate particular temporal events with these
interactions, where both the event type and the timing of the event are
relevant. For example, the relational event model is well suited for data
structures common in animal behaviour studies: behavioural data collected via
proximity loggers provide reliable, highly granular information about social
interactions between animals. The REM has been used to model animal social
networks in cattle (Patison et al., 2015; Tranmer et al., 2015). This model
can also be applied to investigate scenarios where social connectivity is
related to population outcomes that are linked to social factors such as
cooperation (Clutton-Brock, 2009), disease transmission (Silk et al., 2017) or
social learning (Hoppitt and Laland, 2011).

Originally, relational event models were designed for small-scale networks as
their computational costs are quite extensive. In this study, we have shown
how they can nonetheless be applied to more complex networks based on high-
quality standardized databases of ecological data. Further developments of
novel, computationally efficient implementations and user-friendly software
will help quantitative ecologists to make the standard application and
diffusion of the methodology possible and more widely investigate time-
dependent ecological phenomena.

## Supplementary Materials

Supplementary material is available at
http://biostatistics.oxfordjournals.org.

Conflict of Interest: None declared.

## References

  * Abellán et al. (2017) P. Abellán, J. L. Tella, M. Carrete, L. Cardador, and J. D. Anadón.  Climate matching drives spread rate but not establishment success in recent unintentional bird introductions.  _Proceedings of the National Academy of Sciences_ , 114(35):9385–9390, 2017.  doi: https://doi.org/10.1073/pnas.1704815114.
  * Araújo and Peterson (2012) M. B. Araújo and A. T. Peterson.  Uses and misuses of bioclimatic envelope modeling.  _Ecology_ , 93(7):1527–1539, 2012.  doi: https://doi.org/10.1890/11-1930.1.
  * Bagnara et al. (2022) M. Bagnara, L. Nowak, H. J. Boehmer, F. Schöll, F. M. Schurr, and H. Seebens.  Simulating the spread and establishment of alien species along aquatic and terrestrial transport networks–a multi-pathway and high-resolution approach.  _Journal of Applied Ecology_ , 2022.
  * Barbieri et al. (2009) K. Barbieri, O. Keshk, and B. Pollins.  Trading data: Evaluating our assumptions and coding rules.  _Faculty Publications_ , 26, 11 2009.  doi: https://doi.org/10.1177/0738894209343887.
  * Becker (2019) B. Becker.  Colonial Dates Dataset (COLDAT).  Technical Report SFB1342, Socium, Bremen, 2019.
  * Bellard et al. (2016) C. Bellard, B. Leroy, W. Thuiller, J.-F. Rysman, and F. Courchamp.  Major drivers of invasion risks throughout the world.  _Ecosphere_ , 7(3):e01241, 2016.  doi: https://doi.org/10.1002/ecs2.1241.
  * Blackburn et al. (2020) T. Blackburn, G. Edwards, S. Tulaczyk, M. Scudder, G. Piccione, B. Hallet, N. McLean, J. Zachos, B. Cheney, and J. Babbe.  Ice retreat in Wilkes Basin of East Antarctica during a warm interglacial.  _Nature_ , 583:554–559, 07 2020.  doi: https://doi.org/10.1038/s41586-020-2484-5.
  * Blackburn et al. (2011) T. M. Blackburn, P. Pyšek, S. Bacher, J. T. Carlton, R. P. Duncan, V. Jarošík, J. R. Wilson, and D. M. Richardson.  A proposed unified framework for biological invasions.  _Trends in Ecology & Evolution_, 26(7):333–339, 2011.  doi: https://doi.org/10.1016/j.tree.2011.03.023.
  * Bonnamour et al. (2021) A. Bonnamour, J. M. Gippet, and C. Bertelsmeier.  Insect and plant invasions follow two waves of globalisation.  _Ecology letters_ , 24(11):2418–2426, 2021.  doi: https://doi.org/10.1111/ele.13863.
  * Breslow (1975) N. E. Breslow.  Covariance analysis of censored survival data.  _Biometrics_ , 30(1):89–99, 1975.
  * Brown et al. (2021) M. E. Brown, R. Oviedo Prieto, J. D. Corbin, J. H. Ness, R. Borroto-Páez, T. S. McCay, and M. S Farnsworth.  Plant pirates of the caribbean: is cuba sheltered by its revolutionary economy?  _Frontiers in Ecology and the Environment_ , 19(4):208–215, 2021.  doi: https://doi.org/10.1002/fee.2311.
  * Butts (2008) C. T. Butts.  A relational event framework for social action.  _Sociological Methodology_ , 38(1):155–200, 2008\.  doi: https://doi.org/10.1111/j.1467-9531.2008.00203.x.
  * Carpio et al. (2020) A. J. Carpio, Y. Álvarez, J. Oteros, F. León, and F. S. Tortosa.  Intentional introduction pathways of alien birds and mammals in Latin America.  _Global Ecology and Conservation_ , 22:e00949, 2020.  doi: https://doi.org/10.1016/j.gecco.2020.e00949.
  * Chapman et al. (2017) D. Chapman, B. V. Purse, H. E. Roy, and J. M. Bullock.  Global trade networks determine the distribution of invasive non-native species.  _Global Ecology and Biogeography_ , 26(8):907–917, 2017.  doi: https://doi.org/10.1111/geb.12599.
  * Chapman et al. (2016) D. S. Chapman, L. Makra, R. Albertini, M. Bonini, A. Páldy, V. Rodinkova, B. Šikoparija, E. Weryszko-Chmielewska, and J. M. Bullock.  Modelling the introduction and spread of non-native species: International trade and climate change drive ragweed invasion.  _Global Change Biology_ , 22(9):3067–3079, 2016\.  doi: https://doi.org/10.1111/gcb.13220.
  * Chytrỳ et al. (2009) M. Chytrỳ, P. Pyšek, J. Wild, J. Pino, L. C. Maskell, and M. Vilà.  European map of alien plant invasions based on the quantitative assessment across habitats.  _Diversity and Distributions_ , 15(1):98–107, 2009.  doi: https://doi.org/10.1111/j.1472-4642.2008.00515.x.
  * Clutton-Brock (2009) T. Clutton-Brock.  Cooperation between non-kin in animal societies.  _Nature_ , 462(7269):51–57, 2009.  doi: https://doi.org/10.1038/nature08366.
  * Cox (1972) D. R. Cox.  Regression models and life-tables.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 34(2):187–202, 1972.
  * Cox (1975) D. R. Cox.  Partial likelihood.  _Biometrika_ , 62(2):269–276, 1975.
  * Dawson et al. (2017) W. Dawson, D. Moser, M. Van Kleunen, H. Kreft, J. Pergl, P. Pyšek, P. Weigelt, M. Winter, B. Lenzner, T. M. Blackburn, et al.  Global hotspots and correlates of alien species richness across taxonomic groups.  _Nature Ecology & Evolution_, 1(7):1–7, 2017\.  doi: https://doi.org/10.1038/s41559-017-0186.
  * Del-Claro and Torezan-Silingardi (2021a) K. Del-Claro and H. Torezan-Silingardi.  _Plant-Animal Interactions Source of Biodiversity_.  Springer Nature, 05 2021a.  ISBN 978-3-030-66876-1.  doi: https://doi.org/10.1007/978-3-030-66877-8.
  * Del-Claro and Torezan-Silingardi (2021b) K. Del-Claro and H. M. Torezan-Silingardi.  An evolutionary perspective on plant-animal interactions.  _Plant-Animal Interactions: Source of Biodiversity_ , page 1, 2021b.  doi: https://doi.org/10.1007/978-3-030-66877-8˙1.
  * Dijkstra and Border (2018) P. D. Dijkstra and S. E. Border.  How does male–male competition generate negative frequency-dependent selection and disruptive selection during speciation?  _Current Zoology_ , 64(1):89–99, 2018.  doi: https://doi.org/10.1093/cz/zox079.
  * DuBois and Smyth (2010) C. DuBois and P. Smyth.  Modeling relational events via latent classes.  In _Proceedings of the 16th ACM SIGKDD international conference on Knowledge discovery and data mining_ , pages 803–812, 2010.  doi: https://doi.org/10.1145/1835804.1835906.
  * Dullinger et al. (2009) S. Dullinger, I. Kleinbauer, J. Peterseil, M. Smolik, and F. Essl.  Niche based distribution modelling of an invasive alien plant: effects of population status, propagule pressure and invasion history.  _Biological Invasions_ , 11(10):2401–2414, 2009\.  doi: https://doi.org/10.1007/s10530-009-9424-5.
  * Dyer et al. (2017) E. E. Dyer, P. Cassey, D. W. Redding, B. Collen, V. Franks, K. J. Gaston, K. E. Jones, S. Kark, C. D. L. Orme, and T. M. Blackburn.  The global distribution and drivers of alien bird species richness.  _PLoS Biology_ , 15(1):e2000942, 2017.  doi: https://doi.org/10.1371/journal.pbio.2000942.
  * Elith and Leathwick (2009) J. Elith and J. R. Leathwick.  Species distribution models: ecological explanation and prediction across space and time.  _Annual Review of Ecology, Evolution, and Systematics_ , 40:677–697, 2009.  doi: https://doi.org/10.1146/annurev.ecolsys.110308.120159.
  * Eschen et al. (2015) R. Eschen, A. Roques, and A. Santini.  Taxonomic dissimilarity in patterns of interception and establishment of alien arthropods, nematodes and pathogens affecting woody plants in e urope.  _Diversity and Distributions_ , 21(1):36–45, 2015\.  doi: https://doi.org/10.1111/ddi.12267.
  * Essl et al. (2011) F. Essl, S. Dullinger, W. Rabitsch, P. E. Hulme, K. Hülber, V. Jarošík, I. Kleinbauer, F. Krausmann, I. Kühn, W. Nentwig, et al.  Socioeconomic legacy yields an invasion debt.  _Proceedings of the National Academy of Sciences_ , 108(1):203–207, 2011.  doi: https://doi.org/10.1073/pnas.1011728108.
  * Essl et al. (2019) F. Essl, W. Dawson, H. Kreft, J. Pergl, P. Pyšek, M. Van Kleunen, P. Weigelt, T. Mang, S. Dullinger, B. Lenzner, D. Moser, N. Maurel, H. Seebens, A. Stein, E. Weber, C. Chatelain, Inderjit, P. Genovesi, J. Kartesz, O. Morozova, M. Nishino, P. M. Nowak, S. Pagad, W.-S. Shu, and M. Winter.  Drivers of the relative richness of naturalized and invasive plant species on Earth.  _AoB PLANTS_ , 11(5), 09 2019.  ISSN 2041-2851.  doi: https://doi.org/10.1093/aobpla/plz051.  URL https://doi.org/10.1093/aobpla/plz051.  plz051.
  * Finch et al. (2021) D. M. Finch, J. L. Butler, J. B. Runyon, C. J. Fettig, F. F. Kilkenny, S. Jose, S. J. Frankel, S. A. Cushman, R. C. Cobb, J. S. Dukes, , J. A. Hicke, and S. K. Amelon.  Effects of climate change on invasive species.  In T. M. Poland, T. Patel-Weynand, D. M. Finch, C. F. Miniat, D. C. Hayes, and V. M. Lopez, editors, _Invasive Species in Forests and Rangelands of the United States_ , pages 57–83. Springer, 2021.
  * Floerl et al. (2009) O. Floerl, G. Inglis, K. Dey, and A. Smith.  The importance of transport hubs in stepping-stone invasions.  _Journal of Applied Ecology_ , 46(1):37–45, 2009\.  doi: https://doi.org/10.1111/j.1365-2664.2008.01540.x.
  * Goodwin (1960) D. Goodwin.  Comparative ecology of pigeons in inner London.  _British Birds_ , 53(5):201–212, 1960.
  * Gyimesi and Lensink (2012) A. Gyimesi and R. Lensink.  Egyptian goose Alopochen aegyptiaca: An introduced species spreading in and from the Netherlands.  _Wildfowl_ , 62:128–145, 12 2012.
  * Hijmans et al. (2017) R. J. Hijmans, E. Williams, and C. Vennes.  geosphere: Spherical trigonometry. r package version 1.5-7, 2017.
  * Hoppitt and Laland (2011) W. Hoppitt and K. N. Laland.  Detecting social learning using networks: a users guide.  _American Journal of Primatology_ , 73(8):834–844, 2011.  doi: https://doi.org/10.1002/ajp.20920.
  * Hurtt et al. (2011) G. Hurtt, L. Chini, S. Frolking, R. Betts, J. Feddema, G. Fischer, J. Fisk, K. Hibbard, R. Houghton, A. Janetos, C. Jones, G. Kindermann, T. Kinoshita, K. K. Goldewijk, K. Riahi, and E.  Harmonization of land-use scenarios for the period 1500–2100: 600 years of global gridded annual land-use transitions, wood harvest, and resulting secondary lands.  _Climatic Change_ , 109(1):117–161, November 2011\.  doi: https://doi.org/10.1007/s10584-011-0153-2.  URL https://ideas.repec.org/a/spr/climat/v109y2011i1p117-161.html.
  * Isaac et al. (2011) N. J. Isaac, K. L. Cruickshanks, A. M. Weddle, J. Marcus Rowcliffe, T. M. Brereton, R. L. Dennis, D. M. Shuker, and C. D. Thomas.  Distance sampling and the challenge of monitoring butterfly populations.  _Methods in Ecology and Evolution_ , 2(6):585–594, 2011.  doi: https://doi.org/10.1111/j.2041-210X.2011.00109.x.
  * Juozaitienė and Wit (2022) R. Juozaitienė and E. C. Wit.  Nodal heterogeneity may induce ghost triadic effects in relational event models, 2022.  URL https://arxiv.org/abs/2203.16386.
  * Kassambara et al. (2021) A. Kassambara, M. Kosinski, and P. Biecek.  _survminer: Drawing Survival Curves using ’ggplot2’_ , 2021.  URL https://CRAN.R-project.org/package=survminer.  R package version 0.4.9.
  * Khamis et al. (2012) F. M. Khamis, D. K. Masiga, S. A. Mohamed, D. Salifu, M. de Meyer, and S. Ekesi.  Taxonomic identity of the invasive fruit fly pest, Bactrocera invadens: concordance in morphometry and dna barcoding.  _PLoS ONE_ , 7(9), 2012.  doi: https://doi.org/10.1371/journal.pone.0044862.
  * Kirk and Terry (2003) W. D. Kirk and L. I. Terry.  The spread of the western flower thrips Frankliniella occidentalis (pergande).  _Agricultural and Forest Entomology_ , 5(4):301–310, 2003.  doi: https://doi.org/10.1046/j.1461-9563.2003.00192.x.
  * Křivánek et al. (2006) M. Křivánek, P. Pyšek, and V. Jarošík.  Planting history and propagule pressure as predictors of invasion by woody species in a temperate region.  _Conservation Biology_ , 20(5):1487–1498, 2006\.  doi: https://doi.org/10.1111/j.1523-1739.2006.00477.x.
  * Lee-Yaw et al. (2022) J. Lee-Yaw, J. L. McCune, S. Pironon, and S. N. Sheth.  Species distribution models rarely predict the biology of real populations.  _Ecography_ , 2022(6):e05877, 2022.  doi: https://doi.org/10.1111/ecog.05877.
  * Lenzner et al. (2022) B. Lenzner, G. Latombe, A. Schertler, H. Seebens, Q. Yang, M. Winter, P. Weigelt, M. van Kleunen, P. Pyšek, J. Pergl, et al.  Naturalized alien floras still carry the legacy of european colonialism.  _Nature Ecology & Evolution_, 6(11):1723–1732, 2022.  doi: https://doi.org/10.1038/s41559-022-01865-1.
  * Lovell et al. (2021) R. S. Lovell, T. M. Blackburn, E. E. Dyer, and A. L. Pigot.  Environmental resistance predicts the spread of alien species.  _Nature Ecology & Evolution_, 5(3):322–329, 2021.  doi: https://doi.org/10.1038/s41559-020-01376-x.
  * Luppold (1988) W. G. Luppold.  _Hardwood import trends_ , volume 619.  US Department of Agriculture, Forest Service, Northeastern Forest Experiment, 1988.  doi: https://doi.org/10.2737/NE-RP-619.
  * Mayer et al. (2003) J. Mayer, A. Butkevicius, A. Kadri, and J. Pizarro.  Dynamic products in world exports.  _Review of World Economics_ , 139(4):762–795, 2003.  doi: https://doi.org/10.1007/BF02653112.
  * Meurisse et al. (2019) N. Meurisse, D. Rassati, B. P. Hurley, E. G. Brockerhoff, and R. A. Haack.  Common pathways by which non-native forest insects move internationally and domestically.  _Journal of Pest Science_ , 92:13–27, 2019.  doi: https://doi.org/10.1007/s10340-018-0990-0.
  * Milonas and Partsinevelos (2017) P. Milonas and G. Partsinevelos.  The pink hibiscus mealybug Maconellicoccus hirsutus (green)(hemiptera: Pseudococcidae) in Greece.  _Hellenic Plant Protection Journal_ , 10(2):80–83, 2017.  doi: https://doi.org/10.1515/hppj-2017-0009.
  * Ni and Hulme (2021) M. Ni and P. E. Hulme.  Botanic gardens play key roles in the regional distribution of first records of alien plants in china.  _Global Ecology and Biogeography_ , 30(8):1572–1582, 2021.  doi: https://doi.org/10.1111/geb.13319.
  * Paini and Yemshanov (2012) D. R. Paini and D. Yemshanov.  Modelling the arrival of invasive organisms via the international marine shipping network: a Khapra beetle study.  _PloS One_ , 7(9):e44589, 2012.  doi: https://doi.org/10.1371/journal.pone.0044589.
  * Parsa et al. (2012) S. Parsa, T. Kondo, and A. Winotai.  The cassava mealybug (Phenacoccus manihoti) in Asia: First records, potential distribution, and an identification key.  _PloS One_ , 7(10):e47675, 2012.  doi: 10.1371/journal.pone.0047675.
  * Patison et al. (2015) K. Patison, E. Quintane, D. Swain, G. Robins, and P. Pattison.  Time is of the essence: an application of a relational event model for animal social networks.  _Behavioral Ecology and Sociobiology_ , 69(5):841–855, 2015.  doi: https://doi.org/10.1007/s00265-015-1883-3.
  * Perry and Wolfe (2013) P. O. Perry and P. J. Wolfe.  Point process modelling for directed interaction networks.  _Journal of the Royal Statistical Society: SERIES B: Statistical Methodology_ , pages 821–849, 2013.  doi: https://doi.org/10.1111/rssb.12013.
  * Polce et al. (2011) C. Polce, W. E. Kunin, J. C. Biesmeijer, J. Dauber, O. L. Phillips, and A. F. S. Network.  Alien and native plants show contrasting responses to climate and land use in Europe.  _Global Ecology and Biogeography_ , 20(3):367–379, 2011.  doi: https://doi.org/10.1111/j.1466-8238.2010.00619.x.
  * Pyšek et al. (2010) P. Pyšek, V. Jarošík, P. E. Hulme, I. Kühn, J. Wild, M. Arianoutsou, S. Bacher, F. Chiron, V. Didžiulis, F. Essl, et al.  Disentangling the role of environmental and human pressures on biological invasions across Europe.  _Proceedings of the National Academy of Sciences_ , 107(27):12157–12162, 2010.  doi: https://doi.org/10.1073/pnas.1002314107.
  * Pyšek et al. (2020) P. Pyšek, P. E. Hulme, D. Simberloff, S. Bacher, T. M. Blackburn, J. T. Carlton, W. Dawson, F. Essl, L. C. Foxcroft, P. Genovesi, et al.  Scientists’ warning on invasive alien species.  _Biological Reviews_ , 95(6):1511–1534, 2020\.
  * R Core Team (2022) R Core Team.  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria, 2022.  URL https://www.R-project.org/.
  * Richardson and Rejmánek (2004) D. M. Richardson and M. Rejmánek.  Conifers as invasive aliens: a global survey and predictive framework.  _Diversity and Distributions_ , 10(5-6):321–331, 2004.  doi: https://doi.org/10.1111/j.1366-9516.2004.00096.x.
  * Şahin and Arslangündoğdu (2019) D. Şahin and Z. Arslangündoğdu.  Breeding status and nest characteristics of rose-ringed (Psittacula krameri) and alexandrine parakeets (Psittacula eupatria) in Istanbul’s city parks.  _Applied Ecology and Environmental Research_ , 17(2):2461–2471, 2019.  doi: http://dx.doi.org/10.15666/aeer/1702˙24612471.
  * Sánchez-Ortiz et al. (2020) K. Sánchez-Ortiz, K. J. Taylor, A. De Palma, F. Essl, W. Dawson, H. Kreft, J. Pergl, P. Pyšek, M. van Kleunen, P. Weigelt, et al.  Effects of land-use change and related pressures on alien and native subsets of island communities.  _PloS One_ , 15(12):e0227169, 2020.  doi: https://doi.org/10.1371/journal.pone.0227169.
  * Seebens et al. (2015) H. Seebens, F. Essl, W. Dawson, N. Fuentes, D. Moser, J. Pergl, P. Pyšek, M. van Kleunen, E. Weber, M. Winter, et al.  Global trade will accelerate plant invasions in emerging economies under climate change.  _Global Change Biology_ , 21(11):4128–4140, 2015\.  doi: https://doi.org/10.1111/gcb.13021.
  * Seebens et al. (2017) H. Seebens, T. Blackburn, E. Dyer, P. Genovesi, P. Hulme, J. Jeschke, S. Pagad, P. Pyšek, M. Winter, M. Arianoutsou, S. Bacher, B. Blasius, G. Brundu, C. Capinha, L. Celesti-Grapow, W. Dawson, S. Dullinger, N. Fuentes, H. Jäger, and F. Essl.  No saturation in the accumulation of alien species worldwide.  _Nature Communications_ , 8, 02 2017.  doi: https://doi.org/10.1038/ncomms14435.
  * Seebens et al. (2018) H. Seebens, T. Blackburn, E. Dyer, P. Genovesi, P. Hulme, J. Jeschke, S. Pagad, P. Pyšek, M. van Kleunen, M. Winter, M. Ansong, M. Arianoutsou, S. Bacher, B. Blasius, E. Brockerhoff, G. Brundu, C. Capinha, C. Causton, L. Celesti-Grapow, and F. Essl.  Global rise in emerging alien species results from increased accessibility of new source pools.  _Proceedings of the National Academy of Sciences_ , 115:201719429, 02 2018.  doi: https://doi.org/10.1073/pnas.1719429115.
  * Silk et al. (2017) M. J. Silk, D. P. Croft, R. J. Delahay, D. J. Hodgson, N. Weber, M. Boots, and R. A. McDonald.  The application of statistical network models in disease research.  _Methods in Ecology and Evolution_ , 8(9):1026–1041, 2017.  doi: https://doi.org/10.1111/2041-210X.12770.
  * Sofaer and Jarnevich (2017) H. R. Sofaer and C. S. Jarnevich.  Accounting for sampling patterns reverses the relative importance of trade and climate for the global sharing of exotic plants.  _Global Ecology and Biogeography_ , 26(6):669–678, 2017.  doi: https://doi.org/10.1111/geb.12577.
  * Tedeschi et al. (2021a) L. Tedeschi, D. Biancolini, C. Capinha, C. Rondinini, and F. Essl.  Invasive alien mammals of European union concern.  _bioRxiv_ , 2021a.  doi: https://doi.org/10.1101/2021.04.21.440832.  URL https://www.biorxiv.org/content/early/2021/04/22/2021.04.21.440832.
  * Tedeschi et al. (2021b) L. Tedeschi, D. Biancolini, C. Capinha, C. Rondinini, and F. Essl.  Introduction, spread, and impacts of invasive alien mammal species in Europe.  _Mammal Review_ , 2021b.  doi: https://doi.org/10.1111/mam.12277.
  * Therneau (2020) T. M. Therneau.  _coxme: Mixed Effects Cox Models_ , 2020.  URL https://CRAN.R-project.org/package=coxme.  R package version 2.2-16.
  * Therneau (2022) T. M. Therneau.  _A Package for Survival Analysis in R_ , 2022.  URL https://CRAN.R-project.org/package=survival.  R package version 3.3-1.
  * Tranmer et al. (2015) M. Tranmer, C. S. Marcum, F. B. Morton, D. P. Croft, and S. R. de Kort.  Using the relational event model (rem) to investigate the temporal dynamics of animal social networks.  _Animal Behaviour_ , 101:99–105, 2015.
  * Turbelin et al. (2017) A. J. Turbelin, B. D. Malamud, and R. A. Francis.  Mapping the global state of invasive alien species: patterns of invasion and policy responses.  _Global Ecology and Biogeography_ , 26(1):78–92, 2017.  doi: https://doi.org/10.1111/geb.12517.
  * van Kleunen et al. (2018) M. van Kleunen, P. Pyšek, W. Dawson, F. Essl, H. Kreft, J. Pergl, P. Weigelt, A. Stein, S. Dullinger, C. König, B. Lenzner, N. Maurel, D. Moser, H. Seebens, J. Kartesz, M. Nishino, A. Aleksanyan, M. Ansong, L. Antonova, and M. Winter.  The global naturalized alien flora (glonaf) database.  _Ecology_ , 100, 10 2018.  doi: https://doi.org/10.1002/ecy.2542.
  * Verdasca et al. (2021) M. J. Verdasca, H. Rebelo, L. Carvalheiro, and R. Rebelo.  Invasive hornets on the road: motorway-driven dispersal must be considered in management plans of Vespa velutina.  _NeoBiota_ , 69:177–198, 11 2021.  doi: 10.3897/neobiota.69.71352.
  * Vu et al. (2017) D. Vu, A. Lomi, D. Mascia, and F. Pallotti.  Relational event models for longitudinal network data with an application to interhospital patient transfers.  _Statistics in Medicine_ , 36(14):2265–2287, 2017\.  doi: https://doi.org/10.1002/sim.7247.
  * Watanabe et al. (2011) S. Watanabe, T. Hajima, K. Sudo, T. Nagashima, T. Takemura, H. Okajima, T. Nozawa, H. Kawase, M. Abe, T. Yokohata, T. Ise, H. Sato, E. Kato, K. Takata, S. Emori, and M. Kawamiya.  Miroc-esm 2010: model description and basic results of cmip5-20c3m experiments.  _Geoscientific Model Development_ , 4(4):845–872, 2011.  doi: https://doi.org/10.5194/gmd-4-845-2011.  URL https://gmd.copernicus.org/articles/4/845/2011/.
  * Westphal et al. (2008) C. Westphal, R. Bommarco, G. Carré, E. Lamborn, N. Morison, T. Petanidou, S. G. Potts, S. P. Roberts, H. Szentgyörgyi, T. Tscheulin, et al.  Measuring bee diversity in different European habitats and biogeographical regions.  _Ecological Monographs_ , 78(4):653–671, 2008\.  doi: https://doi.org/10.1890/07-1292.1.
  * Wilson et al. (2009) J. R. Wilson, E. E. Dormontt, P. J. Prentis, A. J. Lowe, and D. M. Richardson.  Something in the way you move: dispersal pathways affect invasion success.  _Trends in ecology & evolution_, 24(3):136–144, 2009.  doi: https://doi.org/10.1016/j.tree.2008.10.007.
  * Yang et al. (2021) Q. Yang, P. Weigelt, T. S. Fristoe, Z. Zhang, H. Kreft, A. Stein, H. Seebens, W. Dawson, F. Essl, C. König, et al.  The global loss of floristic uniqueness.  _Nature Communications_ , 12(1):1–10, 2021.  doi: https://doi.org/10.1038/s41467-021-27603-y.
  * Zheng and Agresti (2000) B. Zheng and A. Agresti.  Summarizing the predictive power of a generalized linear model.  _Statistics in medicine_ , 19(13):1771–1781, 2000\.  doi: https://doi.org/10.1002/1097-0258(20000715)19:13¡1771::aid-sim485¿3.0.co;2-p.

Variable | Description | Data source
---|---|---
$d_{sc}(t)$ | Closest distance between regions | (Hijmans et al., 2017)
$tr_{sc}(t)$ | International trade | (Barbieri et al., 2009)
$dt_{sc}(t)$ | Near-surface air temperature | (Watanabe et al., 2011)
$l_{c}(t)$ | Sum of cropland and pasture proportions | (Hurtt et al., 2011)
$u_{c}(t)$ | Proportion of urban area | (Hurtt et al., 2011)
$k_{sc}(t)$ | Colonial ties | (Becker, 2019)
$pi_{c}(t)$ | The Alien Species First Records database | (Seebens et al., 2018)
Table 1:  List of variables and data sources used in analysis. Although each of the various data sources contained a wider temporal range, our analysis focused on the period from 1880 until 2005 and we subsequently focused on this temporal range for each of the databases.  Effect |  Definition |  Type |  Justification
---|---|---|---
$d_{sc}(t)$ |  The distance to the nearest region $c$ in which species $s$ is present by time $t$. |  Time-varying variable with effect changing over time |  Geographic distance may influence migratory behaviours and strategies
$tr_{sc}(t)$ |  Logarithm of the sum of annual trade flows between region $c$ and other regions invaded by species $s$ by time $t$ |  Time-varying variable with effect changing over time |  International trade is a key pathway for the global spread of alien species
$dt_{sc}(t)$ |  Minimum absolute difference in temperature between region $c$ and other regions invaded by species $s$ by time $t$ |  Time-varying variable with constant effect |  Climatically suitable ranges ease establishment of alien species
$l_{c}(t)$ |  Sum of cropland and pasture proportions in the region $c$ at time $t$ |  Time-varying variable with constant effect |  Land use changes provide opportunities for the spread and establishment of alien species
$u_{c}(t)$ |  Proportion of urban area in region $c$ at time $t$ |  Time-varying variable with constant effect |  Land use changes provide opportunities for the spread and establishment of alien species
$k_{sc}(t)$ |  Presence of species $s$ at time $t$ in a region within the area of a colonial power |  Time-varying variable with constant effect |  In the Colonial era alien species have been intentionally introduced
$pi_{c}(t)$ |  Weighted number of prior invasions in region $c$ |  Time-varying variable with constant effect |  Temporal variable representing the changing importance of region $c$ over time
$se_{c}$ |  Sampling effort |  Constant variable with constant effect |  The number of native species per region is used as a proxy for sampling effort
$b_{s}$ |  Species invasiveness |  Monadic random effect |  Different species may vary in their spread behaviour
$b_{c}$ |  Region invasibility |  Monadic random effect |  Certain regions may be more attractive destinations than others
$b_{ss^{\prime}}$ |  Affinity of species $s$ and $s^{\prime}$ in their spread patterns |  Dyadic random effect |  Alien species may take similar invasion pathways
Table 2: List of explanatory variables used for model fitting.

$1$$2$Distance$t_{2}~{}-$ABC$\mathcal{R}(t_{2})=\\{(\mbox{\scriptsize\PHtunny},C),~{}(\mbox{\scriptsize\PHdove},A)~{},(\mbox{\scriptsize\PHdove},B)\\}$$d_{sc}(t_{2})=\\{~{}~{}~{}~{}{\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}2}~{}~{}~{}~{},~{}~{}~{}~{}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}3}~{}~{}~{}~{}~{},~{}~{}~{}~{}{\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}2}~{}~{}~{}~{}\\}$$\mbox{Prob}~{}=\\{~{}~{}0.42~{}~{},~{}~{}~{}0.16~{}~{}~{},~{}~{}0.42~{}~{}\\}$$L_{2}(\beta)=\frac{e^{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}3}\beta}}{2e^{{\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}2}\beta}+e^{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}3}\beta}}$$t_{1}~{}-$ABC$\mathcal{R}(t_{1})=~{}\\{~{}(\mbox{\scriptsize\PHtunny},B)~{},(\mbox{\scriptsize\PHtunny},C),(\mbox{\scriptsize\PHdove},A),(\mbox{\scriptsize\PHdove},B)\\}$$d_{sc}(t_{1})=\\{~{}~{}~{}~{}~{}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1}~{}~{}~{}~{}~{},~{}~{}~{}~{}~{}{\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}3}~{}~{}~{}~{},~{}~{}~{}~{}{\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}3}~{}~{}~{}~{},~{}~{}~{}~{}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2}~{}~{}~{}~{}\\}$$\mbox{Prob}~{}=~{}\\{~{}~{}~{}0.61~{}~{}~{},~{}~{}0.08~{}~{},~{}~{}~{}0.08~{}~{},~{}~{}0.23~{}~{}\\}$$L_{1}(\beta)=\frac{e^{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1}\beta}}{e^{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1}\beta}+2e^{{\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}3}\beta}+e^{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2}\beta}}$$t_{0}~{}-$ABC

Figure 1: An example of an invasion process of two species $\\{\mbox{\scriptsize\PHtunny, \PHdove}\\}$ in three countries $\\{A,B,C\\}$ according to the hypothetical hazard function $\lambda_{sc}(t)=\lambda_{0}(t)e^{-1\times d_{sc}(t)}$ assuming that the rate of invasion event depends on the distances to the nearest invaded country $d_{sc}(t)$. The effect size of the distance is negative, i.e., $\beta_{0}=-1$, meaning that long distance invasions are less likely. At time $t_{0}$ species  is present in country A and species  is present in country C. At time $t_{1}$, species  invades country B and at time $t_{2}$ species  invades country A. $\mathcal{R}(t)$ is the risk consisting of all possible invasion events at time $t$ and $d_{sc}(t)$ is the distance to the nearest invaded region at time $t$ of the corresponding invasion event. The probability of each invasion event is calculated following the formula of the event indicator (see section 2.2.1). In the figure, observed events are shown by a box around the element in the risk set. To estimate the effect of distance (i.e., parameter $\beta$) we use a partial maximum likelihood approach. The partial likelihood of $i$th invasion event is denoted as $L_{i}(\beta)$. After observing two invasion events $(\mbox{\scriptsize\PHtunny},B)$ and $(\mbox{\scriptsize\PHdove},A)$ the partial log-likelihood, $\log(L_{1}(\beta))+\log(L_{2}(\beta))$, is maximized at $\hat{\beta}=-0.64$, not far from the true effect size of distance $\beta_{0}=-1$.  Figure 2: A heatmap representing correlations between the explanatory variables in the model. As can be seen, multicollinearity between the variables is low, suggesting that the estimation is stable. Figure 3: Global map of the number of invasion events per country. The total number of analysed invasion events is $n=16403$. Australia is the country with most first records ($n=1685$).  |
---|---
|
Figure 4: Spread trajectories of the four species from different taxonomic groups: (top left) muskrat, (top right) common ragweed, (bottom left) cassava mealybug, (bottom right) monk parakeet. Countries are graded from the country invaded earliest (darker) to the latest (lighter). Black colour represents area originally occupied by the species in 1880. Countries without the presence of the species are shown in white. |  | Model loglik |  $\chi^{2}$ stat | df | p-value
---|---|---|---|---|---
Birds | Only constant effect | -16746 |  |  |
Proposed model | -16597 | 298.4 | 12 | $p<10^{-4}$
Only time-varying effects | -16592.6 | 8.8 | 4 | 0.07
Mammals | Only constant effect | -3513.6 |  |  |
Proposed model | -3299.9 | 427.37 | 12 | $p<10^{-4}$
Only time-varying effects | -3295.7 | 8.4 | 4 | 0.08
Plants | Only constant effect | -88361.0 |  |  |
Proposed model | -85782.0 | 5156.7 | 12 | $p<10^{-4}$
Only time-varying effects | -85772.0 | 20.1 | 12 | 0.07
Insects | Only constant effect | -3163.1 |  |  |
Proposed model | -3136.5 | 53.3 | 8 | $p<10^{-4}$
Only time-varying effects | -3132.2 | 8.5 | 8 | 0.39
Table 3: Summary table of the ANOVA results, showing the preference for the proposed model including several time-varying effects for all the taxonomic groups.  Taxonomic group | $R^{2}$ | Null loglik | Model loglik | AIC | BIC |  $\chi^{2}$ stat | df | p-value
---|---|---|---|---|---|---|---|---
Birds | 18.9 | -20463.1 | -16597.2 | 7695.9 | 7597.1 | 7731.9 | 18 | $p<10^{-4}$
Mammals | 25.6 | -4437.2 | -3299.9 | 2236.6 | 2159.4 | 2274.6 | 19 | $p<10^{-4}$
Plants | 26.1 | -116051.9 | -85782.3 | 60497.2 | 60348.7 | 60539.2 | 21 | $p<10^{-4}$
Insects | 12.3 | -3577.9 | -3136.5 | 852.8 | 794.3 | 882.8 | 15 | $p<10^{-4}$
Table 4: Summary table of the final invasion models selected, showing an overall fit between 12.3% and 26.1% for all taxonomic groups.  effect | plants | insects | mammals | birds
---|---|---|---|---
$pi_{c}(t)$: prior invasions | -0.005* | -0.11* | n.s. | n.s.
$dt_{sc}(t)$: temperature difference | -0.172* | -0.089* | -0.064* | -0.052*
$k_{sc}(t)$: colonial ties | n.s. | n.s. | n.s. | n.s.
$u_{c}(t)$: urban landscape | -18.42* | n.s. | n.s. | n.s.
$se_{c}$: sampling effort | 1.8* | 0.65* | 0.99* | n.s.
$\sigma_{\mbox{\scriptsize spc}}$: species invasiveness | 0.64* | 0.85* | 0.98* | 0.92*
$\sigma_{\mbox{\scriptsize cnt}}$: region invasibility | 3.59* | 1.22* | 1.36* | 1.49*
$\sigma_{\mbox{\scriptsize int}}$: interdependence in spread | n.s. | 0.71* |  0.22* | 0.41*
$\sigma_{\mbox{\scriptsize int}}^{\mbox{\scriptsize sym}}$: symmetric interdependence in spread | n.s. |  0.82* | n.s. |  0.34*
Table 5: Parameter estimates for the non-temporal fixed effects and random effect standard deviations across the four taxonomic groups (* p-value $<0.001$, n.s.: p-value $>0.05$). The bold values indicate the preferred formulation of interdependence in spread effect in each case, i.e., the random effect with the highest trade-off between explanatory power and the number of parameters needed. Figure 5: Effect of trade on the process of alien species spread. After an initial positive effect, in more recent times the effect of trade seems to have become negative. The grey shade area represents the 68% confidence interval. For all taxonomic groups the null-hypotheses of no time-varying trade effect can be rejected, as all the associated p-values $<0.001$. |
---|---
(a) | (b)
Figure 6:  Effect of a) distance and b) cropland and pasture proportion on
rates of invasion. a). Negative effect of distance on the rate of invasion
events has increased over time for mammals and plants. For all taxonomic
groups p-value $<0.001$. b). Crop land and pasture proportions have negative
effects on the rate of invasion events. The grey shade area represents the 68%
confidence interval. For all taxonomic groups the null-hypotheses of no time-
varying distance and land cover effects can be rejected, as all the associated
p-values $<0.001$. Figure 7: Global map of region invasibility, indicating
regions susceptibility to alien species invasions. Countries are graded from
the regions with the highest invasion propensity (black) to the smallest
(light grey), while correcting for all fixed effects in the REM. United
States, South Africa, Australia, New Zealand are geographic hot spots of
invasion for birds. Figure 8: Symmetric insect interdependence in spread
effects. Each node represents a species and an edge depicts a interdependence
in spread propensity. Green edges indicate that the spread of two species
tends to follow each other, while red edges suggest that two species tend to
avoid each other in their spatio-temporal spread. Moreover, the thickness of
the edge represents the strength of the measured interdependence. The provided
graph include only the subset of edges having the highest weight values, i.e.
increasing or decreasing the invasion hazard by a factor of 2. Smaller
interdependence strengths have been removed to make the graph easier to
interpret. Figure 9: Baseline hazard for the underlying invasion speeds of the
different taxonomic groups. Although mammals have shown the quickest pace, it
is tapering off recently, whereas the spread of insects is speeding up since
1980.