epiflows: an R package for risk assessment of travel-related spread of disease

As international travel increases worldwide, new surveillance tools are needed to help identify locations where diseases are most likely to be spread and prevention measures need to be implemented. In this paper we present epiflows, an R package for risk assessment of travel-related spread of disease. epiflows produces estimates of the expected number of symptomatic and/or asymptomatic infections that could be introduced to other locations from the source of infection. Estimates (average and confidence intervals) of the number of infections introduced elsewhere are obtained by integrating data on the cumulative number of cases reported, population movement, length of stay and information on the distributions of the incubation and infectious periods of the disease. The package also provides tools for geocoding and visualization. We illustrate the use of epiflows by assessing the risk of travel-related spread of yellow fever cases in Southeast Brazil in December 2016 to May 2017.


Introduction
Infectious disease outbreaks cause significant suffering and mortality in the affected populations, and damage the health, social and economic well-being of the families affected by diseases as well as producing significant economic costs for local and national governments. As we have seen with Ebola and SARS, disease outbreaks can spread beyond national borders 1 . Travelers can acquire a disease while staying in a foreign country, and then seed new outbreaks in their home country after their return. As international travel increases worldwide, new surveillance tools are needed to help identify locations where diseases are most likely to be spread and prevention measures need to be implemented. This is essential to limit the global spread of local outbreaks.
Recently, Dorigatti et al. 2 developed a method to assess the risk of travel-related international spread of disease by integrating epidemiological and travel (by air, land and water) volume. The model developed by Dorigatti et al. 2 estimates the expected number of infections introduced elsewhere by taking into account population flows, lengths of stay, as well as the variability of the disease incubation and infectious periods. The method was applied to quantify the risk of spread of a recent outbreak of yellow fever in Southeast Brazil in December 2016 to May 2017, and was able to identify the countries that could have received travel-related disease cases capable of seeding local transmission.
In this paper we present epiflows, an R package that implements the method presented by Dorigatti et al. 2 for risk assessment of travel-related spread of disease. Using data on population movement between the location that is source of the infection and other locations, lengths of stay, as well as information about the disease incubation and infectious period distributions, the package allows the estimation of the number of (symptomatic and/or asymptomatic) infections that could be spread to other locations together with uncertainty measures. The package also provides tools for geocoding and visualization of population flows.
The remainder of the paper is organized as follows. First, we briefly describe the modelling framework that is implemented in the epiflows package. Second, we introduce the main components of epiflows including instructions for installation and main functions. Third, we illustrate the use of the package via the assessment of the risk of travel-related spread of yellow fever cases due to population flows between Southeast Brazil and other countries in December 2016 to May 2017. Specifically, we discuss the data required and show how to perform the statistical analyses, how to interpret the results, and the visualization options. Finally, the conclusions are presented.

Model
In this Section we explain the modelling framework presented in 2 for estimating the expected number of infections departing from one infectious location during the incubation or infectious periods. These cases comprise exportations and importations. Exportations refer to the infected residents of the infectious location (i.e. location with sustained disease transmission) that travel to other locations. Importations (also referred to as returning travelers) are people that are infected during a temporary stay in the infectious location and then return to their home location. The following Sections describe how to model exportations and importations to produce the total number of expected cases that could be spread to other locations together with uncertainty measures.

Exportations
Let C S,W denote the cumulative number of infections in location S in time window W. Here, W denotes the temporal window between the first and the last disease case in location S. Note that Dorigatti et al. 2 calculated C S,W by multiplying the number of confirmed and reported yellow fever cases by 10 to account for underreporting of asymptomatic and mild yellow fever cases.

Amendments from Version 2
The previous version of the manuscript had a small error. In Section "Arguments of the estimate_risk_spread() function" we wrote `num_sim` instead of `n_sim`. This has been corrected in the new version.
Any further responses from the reviewers can be found at the end of the article REVISED Let pop S be the resident population of the infectious location S, and W S,D T the number of residents of location S travelling to location D in time window W. The per capita probability that a resident from the infectious location travelled to other location D during the time window W is given by We assume that the incubation period (D E ) and the infectious period (D I ) are random variables, with associated probability distributions that are disease-specific. Using these, we can calculate the probability p i that an infection incubated or is infectious in time window W as minimum , 1 Finally, the number of residents of the infectious location S that are infected and travel abroad during their incubation or infectious period during the time window W can be calculated as That is, E S,D is a product of the cumulative number of infections in location S in time window W, the per capita probability that a resident of S travels to location D, and the probability that an infection incubated or is infectious in time window W.
Note here that if travel data are expressed annually ( ) The probability of returning to the home location while incubating or infectious is given by Finally, the expected number of travelers infected during their stay in the infectious location and returning to their home location before the end of the infectious period can be calculated as the product of the number of travelers, the per capita risk of infection and the probability of returning home while incubating or infectious,

Implementation
The R package epiflows [21] is hosted in the Comprehensive R Archive Network (CRAN) which is the main repository for R packages: http://CRAN.R-project.org/package=epiflows. Users can install epiflows in R by executing the following code: install.packages("epiflows") There is also a development version from GitHub which can be accessed at https://github.com/reconhub/ epiflows. This version of the package may contain new features which are not incorporated in the version on CRAN yet but may be useful for some users. GitHub also includes issue tracking where users can note problems or suggestions for improvements. This development version from GitHub can be installed by using the install_github() function from the R package devtools 3 : install.packages("devtools") library("devtools") install_github("reconhub/epiflows") When installing epiflows, other R packages which epiflows depends on are also automatically installed. These packages include sp 4 for manipulating spatial objects; geosphere 5 for calculating distances between locations; and leaflet 6 for visualization.

Operation
The main function of the package is estimate_risk_spread() which calculates the mean and 95% confidence intervals of the number of cases spread to different locations from an infectious location. It is also possible to use this function to produce a data frame with all simulations (not just the mean and 95% confidence intervals that is computed from the simulations). This permits the user to aggregate the estimates and calculate confidence intervals with different levels using single simulations. To execute this function the following information is needed: • population of the infectious location, • number of infections in the infectious location, and the first and last dates of reported cases, • number of travelers between the infectious location and other locations, • average length of stay of travelers from other locations visiting the infectious location, • distributions of the incubation and infectious periods, • number of simulations to be drawn from the incubation and infectious period distributions, • logical value indicating whether the returned object should be a data frame with all simulations, or a data frame with the mean and lower and upper limits of a 95% confidence interval of the number of infections spread to each location.
Other useful functions are plot() which produces visualizations of population flows between locations, and add_coordinates() which finds the coordinates of the locations.

Use cases
In this Section we provide an example on how to use epiflows to calculate the number of yellow fever cases spreading from south-east Brazil to other countries due to human movement. We show how to define the arguments of the estimate_risk_spread() function, interpret the results, and make visualizations with the population flows.

Data
We use the data YF_flows and YF_locations which are contained in the epiflows package as data (YF_flows) and data(YF_locations), respectively. These data contain the population size, the assumed number of yellow fever infections, dates of first and last case reporting, number of travelers and length of stay for the states of Espirito Santo, Minas Gerais, Rio de Janeiro, Sao Paulo, and for the whole region of Southeast Brazil (which comprises the four states of Espirito Santo, Minas Gerais, Rio de Janeiro and Sao Paulo) in the period December 2016 to May 2017 [19], [20].
Following Dorigatti et al. 2 , the total number of yellow fever infections in each of the Brazilian states was calculated by multiplying the cumulative number of confirmed yellow fever cases reported in 7 by 10 to account for underreporting of asymptomatic and mild yellow fever cases. The dates of first and last case reported in each state were derived as described by Dorigatti et al. 2 . Population data were obtained from the Brazilian Institute of Geography and Statistics website 8 . These data also contain the number of travelers in the specified time window between the states of Espirito Santo, Minas Gerais, Rio de Janeiro, Sao Paulo (and the whole Southeast Brazilian region) and other countries. These estimates were obtained from World Tourism Organization data on the volume of air, land and water border crossings for Brazil for the year 2015 9 , having assumed that travelers were distributed across the Brazilian states according to the relative population density and having accounted for information on the monthly distribution of tourism and on the average duration of stay of international visitors to Brazil 10 , as detailed in 2.

The epiflows object
To aid in data organization between flows and metadata, we have implemented the epiflows object. This inherits the epicontacts object from the epicontacts package 11 , storing three elements: 1. flows -a data frame defining the number of cases flowing from one location to another 2. locations -a data frame listing the locations present in flows and relevant metadata.
3. vars -a dictionary mapping column names in locations to known global variables defined in global_vars(). These global variables are used as default values in estimate_risk_spread().
Because a flow of cases from one location to another can be thought of as a contact with a wider scope, the epiflows object inherits from the epicontacts object, where locations are stored in the linelist element and flows are stored in the contacts element (though the user does not need to interact with these elements by name). By building on the epicontacts object, we ensure that all the methods for sub-setting an object of class epicontacts also applies to epiflows, reducing the maintenance effort.
An epiflows object can be created with the make_epiflows() function by providing a data frame flows with the number of travelers between locations, a data frame locations with information about the locations, and the names of the columns of data frame locations indicating the name of each variable.
In the data frame flows each row represents the number of travelers from one location to the next. flows has at least three columns: columns from and to indicating where the flow starts and ends, respectively, and column n indicating the number of travelers that are in the flow. Data frame YF_flows contains the population flows of the Brazil data. Espirito Santo Spain 3270.500 In data frame locations each row represents a location, and columns specify useful information about the locations such as ID, population, number of cases, dates and length of stay. locations must contain at least one column specifying the location ID used in the flows data frame. YF_locations contains, for each Brazilian state considered in our example, the code (location_code), the population (location_ population), the number of assumed infections in the time window (num_cases_time_window), and the dates of the first and last case reported (first_date_cases and last_date_cases, respectively). It also contains length_of_stay which are the lengths of stay (in days) of the travelers visiting Brazil from other countries. Then, we can create an epiflows object called Brazil_epiflows as follows.

Arguments of the estimate_risk_spread() function
The arguments that need to be specified in estimate_risk_spread() to calculate the cases or infections introduced in other countries are as follows. The first argument is an epiflows object containing the number of travelers between locations, the population size, the number of cases, and the first and last dates of reporting in the infectious location, and the average length of stay in days of travelers from other locations visiting the infectious location.
The second argument of estimate_risk_spread() is location_code which is a character string denoting the infectious location code. We also need to specify the incubation and infectious period distributions. Specifically, we need to provide functions with a single argument n that generate n random incubation and infectious periods. To do this, we can use random generation functions of distributions that are implemented in R including Normal rnorm(), LogNormal rlnorm(), Gamma rgamma(), Weibull rweibull(), and Exponential rexp(). Details about the meaning and arguments of these functions can be obtained by typing ? and the function name (e.g., ?rnorm). We should consider the literature carefully before deciding on appropriate distributions. Examples of systematic reviews of incubation period distributions are 12 and 13. In this example, we use the specified distributions and parameterisation following 14 and 15. Thus, we assume that the incubation period D E is log-normally distributed with mean 4.6 days and variance 2.7 days, and that the infectious period D I is normally distributed with mean 4.5 days and variance 0.6 days. We can define functions incubation() and infectious() as Argument n_sim is the number of simulations to be drawn from the incubation and infectious period distributions. It is recommended to use at least 1,000 simulations. The last argument of estimate_risk_spread() is return_all_simulations. This is a logical value indicating whether the returned object should be a data frame with all simulations (return_all_simulations= TRUE), or a data frame with the mean and lower and upper limits of a 95% confidence interval of the number of infections spread to each location (return_all_simulations = FALSE).

Execution of the estimate_risk_spread() function
Once we have constructed the objects needed to call estimate_risk_spread() we can execute the function and obtain the estimated mean number of cases spread to each country and the 95% confidence intervals. The code to calculate the cases spread from Espirito Santo is the following: set.seed(2018-07-25) res <-estimate_risk_spread(Brazil_epiflows, location_code = "Espirito Santo", r_incubation = incubation, r_infectious = infectious, n_sim = 1e5 ) The results returned by estimate_risk_spread() are stored in the res object. This is a data frame with the columns mean_cases indicating the mean number of cases spread to each location, and lower_limit_95CI and upper_limit_95CI indicating the lower and upper limits of 95% confidence intervals. The result object is shown below. We can plot the results with ggplot() as follows ( Figure 1).

Figure 1. Mean number of yellow fever cases and 95% CI spread from Espirito Santo to other locations.
Note that if we set return_all_simulations equal to TRUE, the result object res will be a data frame with all simulations.

Visualize population flows
We can visualize flows of people travelling between locations using plot() and passing as first parameter an epiflows object containing the population flows, and as second parameter the type of plot we wish to produce. Population flows can be displayed on an interactive map, as a network or as a grid between origins and destinations as described in the following sections.

Flows displayed on an interactive map
We can visualize population flows on an interactive map using plot() with the parameter type="map".
For this option to work, the epiflows object needs to include the longitude and latitude of the locations in decimal degree format. If coordinates are known, they can be added to the epiflows object using the add_coordinates() function from the epiflows package. In our example, the longitude and latitude data are in the data frame YF_coordinates. They can be added as follows.

Brazil_epiflows <-add_coordinates(Brazil_epiflows, coordinates = YF_coordinates[, -1])
If coordinates are unknown, we may resort to one of the freely available tools for geocoding. For example, we can use the geocode() function from the ggmap package 16 . This function finds the latitude and longitude of a given location using either the Data Science Toolkit or Google Maps. We can also use add_coordinates() which uses geocode() to find the coordinates and directly add them to the epiflows object as follows.

Brazil_epiflows <-add_coordinates(Brazil_epiflows, overwrite = TRUE)
Once we have assigned coordinates to the epiflows object, we can use plot() with type="map" to visualize the population flows between locations in an interactive map ( Figure 2).

plot(Brazil_epiflows, type = "map")
The produced map can be zoomed and permits an easy examination of flows. plot() uses the gcIntermediate() function from the geosphere package 5 to obtain the great circle arcs between locations, and then uses the leaflet package 6 to create an interactive map with the connection lines. The connection lines are coloured according to flow volume, and as the mouse passes over the lines, lines highlight and information about connections is shown. We can also include parameters to specify a title, the center of the map or a color palette. An interactive version of this visualization is shown here: https://www.repidemicsconsortium.org/ epiflows/articles/introduction.html#introduction-epiflows-map.

Flows displayed as a network
Population flows can also be displayed as a dynamic network using plot() with type = "network" (Figure 3).

plot(Brazil_epiflows, type = "network")
This option uses the package visNetwork 17 to show the locations as nodes of a network and connections between them representing population flows. This plot is interactive and it is possible to highlight a given location and examine its population flows, as well as its population, number of cases, dates and length of stay. This type of plot can be used when coordinates of locations are missing. An interactive version of this plot can be viewed here: https://www.repidemicsconsortium.org/epiflows/articles/introduction.html#introduction-epiflows-vis.

Flows displayed as a grid between origins and destinations
Finally, population flows can also be shown as a grid between locations with the option type="grid" (Figure 4).

plot(Brazil_epiflows, type = "grid")
This plot shows flows between locations as points by positioning origins and destination in y and x axes, respectively. When using this option, additional arguments can be passed to set the size, color or shape of the points as in function geom_point() of package ggplot2 18 . As the network plot, the grid plot can be used when coordinates of locations are missing.  International travel has an important role in the spread of infectious diseases across national borders. We think the epiflows package represents a useful tool for disease surveillance that can help public health officials identify locations where diseases are most likely to spread and prevention measures are most needed.

Noam Ross
EcoHealth Alliance, New York City, NY, USA The authors present and provide a tutorial to epiflows, an R package for calculating the risk of travel-related disease export from an epidemic area. It is a useful implementation of an algorithm, with associated visualization tools. It is technically sound though the scaffolding around the core algorithm is somewhat over-engineered.
Both the package design and paper description imply this package is designed for rapid risk assessment. My comments are primarily in regard to the clarity of the description and the usability of the package API in this context.
The `global_vars()` function is a thin wrapper around R's `options()` mechanism that obfuscates what is actually happening, and the name itself is somewhat confusing. (There are many different uses and abuses of global variables in R). It is difficult to see how this function improves over simply telling the user that default variables are defined by `epiflows.vars` in R's options mechanism (epiflows.varnames might be clearer).
While it is mentioned that the `epiflows` object inherits from `epicontacts`, it is not clear what this means and how it is relevant to the user. Given the epiflows object has different contents than thè epicontacts` object, it should be explained. I see in the package vignette that subsetting methods are inherited, but given that the object contents are different, this should be demonstrated in the paper. Otherwise it is hard to see what the advantage of this object is at all, over simply passing data frames to the algorithm function.
The paper (as well as the code demonstration and vignette in the package), spends considerable time on the conversion of the partially processed data in `YF_Brazil` to data appropriately structured to be used in the `make_epiflows()` and `estimate_risk_spread()` function. This is confusing and not very useful -most users will not have data in exactly the format of `YF_Brazil`, and it distracts from the description of the core package functionality. It would be far clearer to introduce and demonstrate the functions first using data in its ready-to-analyze form. If the authors wish to provide an example of an actual full workflow, they should start with actual raw source data from the the supplementary data. A useful addition would also be to describe possible sources of the data needed.
The distributions of incubation and infectious periods is important and glossed over rather quickly here. First, on page 3, the text reads, "The method assumes that the incubation period TE and the infectious period TI follow specific probability distributions." This is unclear, I think "disease-specific" would be clearer. Moreover, in the code tutorial, these distributions are simply assumed. It would make more sense to describe why such distributions are selected and the source of the parameters, e.g., "For yellow fever we choose these distributions based on clinical literature describing the disease course as X to Y days of incubation and V to Z days infectious period (Citation 1, Citation 2). Lognormal or Gamma distributions are typically used for these distributions..." maintenance effort. We have clarified this in Section "The epiflows object" of the manuscript.
3. In order to present an example that is as clear as possible for readers, we have modified the yellow fever example and now we do not describe datasets that do not have ready-to-analyze form. Now we start by describing the YF_flows and YF_locations objects which can be directly passed to the make_epiflows() function to create the Brazil_epiflows object.
4. Thank you for pointing out the importance of the distributions of incubation and infectious periods. In Section "Exportations" we have replaced "the method assumes that the incubation period T_E and the infectious period T_I follow specific probability distributions." by "We assume that the incubation period (D_E) and the infectious period (D_I) are random variables, with associated probability distributions that are disease-specific." In Section "Arguments of the estimate_risk_spread() function", we have cited the papers we used to choose the incubation and infectious period distributions of the yellow fever example. We decided not to provide data for incubation and infectious periods of other pathogens since we expect users to have some knowledge of the life history of the disease they want to apply it to. We have also mentioned that users should consider the literature carefully before deciding on appropriate distributions, and we reference two review papers that may be useful.
5. Thank you for the notation review. Now we have changed the notation of the incubation period T_E by D_E, and infectious period T_I by D_I, where letter D stands for duration. Now letter T is only used for travelers.
No competing interests were disclosed.

Competing Interests:
The benefits of publishing with F1000Research: Your article is published within days, with no editorial bias You can publish traditional articles, null/negative results, case reports, data notes and more The peer review process is transparent and collaborative Your article is indexed in PubMed after passing peer review Dedicated customer support at every stage For pre-submission enquiries, contact research@f1000.com