A. Example using ISO codes
A.1. Trade data preparation
Structure of trade data
This example uses simulated trade data for Northern American countries, consisting of a list of data frames with the required data. These data use country identifiers by ISO 3166-1 (alpha-2) codes. Trade data are arranged in three-months time periods.
The load_csv() function included in the
qPRAentry package can be used to import the data from CSV
files.
Total quantity of commodity imported from third countries.
This data frame must contain the columns: reporter (importing countries, in this case by ISO codes), partner (exporting countries), value (quantity of commodity), and time_period (time period of the trade activity).
Using the example data, we select imports from all third countries (column partner):
extra_total <- datatrade_NorthAm$extra_import
head(extra_total)
#> reporter partner time_period value
#> 1 BM CNTR_1 April-June 73.58
#> 2 BM CNTR_2 April-June 17.55
#> 3 BM CNTR_3 April-June 26.61
#> 4 BM CNTR_4 April-June 19349.36
#> 5 BM CNTR_5 April-June 8070.73
#> 6 CA CNTR_1 April-June 6628.23Quantity of commodity imported from third countries where the pest is present.
This data frame must contain the columns: reporter (importing countries, in this case by ISO codes), partner (exporting countries where the pest is present), value (quantity of commodity), and time_period (time period of the trade activity).
Here, we assume that the pest is present in countries “CNTR_1” and “CNTR_2”:
library(dplyr)
CNTR_pest <- c("CNTR_1", "CNTR_2")
extra_pest <- datatrade_NorthAm$extra_import %>% filter(partner%in%CNTR_pest)
head(extra_pest)
#> reporter partner time_period value
#> 1 BM CNTR_1 April-June 73.58
#> 2 BM CNTR_2 April-June 17.55
#> 3 CA CNTR_1 April-June 6628.23
#> 4 CA CNTR_2 April-June 910.00
#> 5 GL CNTR_1 April-June 1403.07
#> 6 GL CNTR_2 April-June 6783.81Quantity of commodity traded between countries of interest.
This data frame must contain the columns: reporter (importing countries, in this case by ISO codes), partner (exporting countries, in this case by ISO codes), value (quantity of commodity), and time_period (time period of the trade activity):
intra_trade <- datatrade_NorthAm$intra_trade
head(intra_trade)
#> reporter partner time_period value
#> 1 BM CA April-June 792.86
#> 2 BM GL April-June 1291.80
#> 3 BM PM April-June 830.42
#> 4 BM US April-June 11.57
#> 5 CA BM April-June 608.07
#> 6 CA GL April-June 6289.37Quantity of commodity produced internally in each country of interest.
This data frame must contain the columns: reporter (producing countries, in this case by ISO codes), value (quantity of commodity), and time_period (time period of production):
internal_production <- datatrade_NorthAm$internal_production
head(internal_production)
#> reporter time_period value
#> 1 BM January-March 119625.01
#> 2 CA January-March 55555.83
#> 3 GL January-March 17790.80
#> 4 PM January-March 70680.64
#> 5 US January-March 45125.31
#> 6 BM April-June 79240.12Generation of the TradeData object from the above data
frames
The trade_data() function assembles the trade data to
generate a TradeData object needed to subsequently
calculate the quantity of potentially infested/infected commodity
entering each country of interest.
Using the arguments filter_IDs and
filter_period we can select the countries and time periods
of interest, respectively. If nothing is specified in these arguments,
by default all countries and time periods included in the data will be
selected. For this example, the United States (US) and Canada (CA) are
selected, and for the time periods January-March and April-June.
trade_NorthAm <- trade_data(extra_total = extra_total,
extra_pest = extra_pest,
intra_trade = intra_trade,
internal_production = internal_production,
filter_IDs = c("US", "CA"),
filter_period = c("January-March", "April-June"))See total trade:
head(trade_NorthAm$total_trade)
#> country_IDs time_period extra_total extra_pest intra_import intra_export
#> 1 US January-March 7796.46 6888.32 5448.33 35.49
#> 2 US April-June 3567.79 1379.75 12.10 1035.70
#> 3 CA January-March 23129.97 470.94 35.49 5448.33
#> 4 CA April-June 14003.04 7538.23 1035.70 12.10
#> internal_production total_available export_prop
#> 1 45125.31 52921.77 1
#> 2 17928.84 21496.63 1
#> 3 55555.83 78685.80 1
#> 4 16840.86 30843.90 1See trade between countries:
head(trade_NorthAm$intra_trade)
#> reporter partner time_period value export_prop
#> 1 CA US April-June 1035.70 1
#> 2 CA US January-March 35.49 1
#> 3 US CA April-June 12.10 1
#> 4 US CA January-March 5448.33 1Below is an example of how to visualise data using ISO 3166-1
(alpha-2) country codes, displaying the total quantity of commodity
available in each country. The plot_countries() function
can be used to display other data organised by using ISO 3166-1
(alpha-2) country codes. This function allows to incorporate other
utilities of the ggplot2 package.
library(ggplot2)
plot_countries(data = trade_NorthAm$total_trade,
iso_col = "country_IDs",
values_col = "total_available",
title = "Total commodity available",
legend_title = "units") +
xlim(-180,-20) + ylim(0,90)
A.2. \(N_{trade}\) - Quantity of potentially infested imported commodity
Calculation of the \(N_{trade}\) for each time period
\(N_{trade}\) summary for the time periods
ntrade_NorthAm_summary <- ntrade(trade_data = trade_NorthAm,
summarise_result = c("mean", "sd",
"quantile(0.025)",
"median",
"quantile(0.975)"))
head(ntrade_NorthAm_summary)
#> country_IDs mean sd median q0.025 q0.975
#> 1 US 4116.27 3959.853 4116.27 1456.2333 6776.307
#> 2 CA 4022.35 5062.035 4022.35 621.9207 7422.779Plot the \(N_{trade}\) median for each country:
plot_countries(data = ntrade_NorthAm_summary,
iso_col = "country_IDs",
values_col = "median",
title = "Ntrade median",
legend_title = "units") +
xlim(-180,-20) + ylim(0,90)
A.3. \(N_{trade}\) redistribution from country level to principal subdivisions
The redist_iso() function requires an additional data
frame with values for each subdivision according to which the
redistribution is performed proportionally. The redistribution is shown
below using simulated commodity consumption data for each territorial
subdivision of the United States (US) and Canada (CA) using ISO 3166-2
codes.
# read data for redistribution and filter subdivisions of US and CA
redist_data <- datatrade_NorthAm$consumption_iso2 %>%
filter(substr(iso_3166_2, 1, 2) %in% c("US", "CA"))
data_redist <- redist_iso(data = ntrade_NorthAm_summary,
iso_col = "country_IDs",
values_col = "median",
redist_data = redist_data,
redist_iso_col = "iso_3166_2",
redist_values_col = "value")
head(data_redist)
#> # A tibble: 6 × 4
#> ISO_1 ISO_2 proportion median
#> <chr> <chr> <dbl> <dbl>
#> 1 US US-WA 0.000308 1.27
#> 2 CA CA-BC 0.0447 180.
#> 3 US US-ID 0.0485 199.
#> 4 US US-MT 0.00857 35.3
#> 5 CA CA-AB 0.217 873.
#> 6 CA CA-SK 0.0457 184.Note: The qPRAentry package currently does not include a
built-in function to plot data at the subdivision level using ISO 3166-2
codes, although it is available using NUTS codes (see B.3). However, users can easily combine the
output with other packages, such as rnaturalearth and
ggplot2, to create maps representing these data.
A.4. Pathway model
The number of potential founder populations (\(NPFP\)) of a pest entering a country or region can be estimated using a pathway model. This model combines the \(N_{trade}\) data with parameters that are relevant in the estimation of the entry of the pest under assessment. Each of these parameters must be assigned a suitable probability distribution. The following shows how the \(N_{trade}\) data obtained above at the country level are combined with other parameters to set up the pathway model and estimate the \(NPFP\).
First, the conceptual model is designed. Here, three parameters have been added in different ways as an illustrative demonstration: \[NPFP = N_{trade} \cdot (1/P1) \cdot ((P2 \cdot 1000) + P3)\]
A distribution is then assigned to each parameter and all relevant
information, along with the desired number of iterations, is
incorporated into the pathway_model() function. Note that
\(N_{trade}\) should not be included in
the model expression.
# pathway model (excluding ntrade)
model <- "(1/P1) * ((P2 * 1000) + P3)"
# parameter distributions
parameters_dist <- list(P1 = list(dist = "unif", min = 0, max = 1),
P2 = list(dist = "beta", shape1 = 1, shape2 = 5),
P3 = list(dist = "norm", mean = 0, sd = 1))
res_pathway <- pathway_model(ntrade_data = ntrade_NorthAm_summary,
IDs_col = "country_IDs",
values_col = "median",
expression = model,
parameters = parameters_dist,
niter = 100)
head(res_pathway)
#> # A tibble: 3 × 8
#> # Groups: country_IDs [3]
#> country_IDs Mean SD Min Q0.25 Median Q0.75 Max
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 US 4508059. 12139686. 31561. 430918. 1067629. 2206858. 79000791.
#> 2 CA 4405199. 11862696. 30841. 421085. 1043269. 2156504. 77198236.
#> 3 Total 8913259. 24002382. 62402. 852003. 2110898. 4363362. 156199027.The result also includes the total \(NPFP\) for the set of countries considered:
res_pathway[res_pathway$country_IDs == "Total",]
#> # A tibble: 1 × 8
#> # Groups: country_IDs [1]
#> country_IDs Mean SD Min Q0.25 Median Q0.75 Max
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 Total 8913259. 24002382. 62402. 852003. 2110898. 4363362. 156199027.Plot the \(NPFP\) median for each country:
plot_countries(data = res_pathway,
iso_col = "country_IDs",
values_col = "Median",
colors = c("#DCE319FF", "#55C667FF", "#33638DFF"),
title = "NPFP median",
legend_title = "NPFP") +
xlim(-180,-20) + ylim(0,90)
