Global costs of biological invasions R codes and scripts
INVACOST team
30/01/2021
Foreword
This page contains all the scripts and analyses used for the paper on the global costs of biological invasions. This code also includes exploratory analyses which were not used in the paper.
The code is based on the invacost R package which is available on CRAN. To read a step-by-step tutorial, see here: https://github.com/Farewe/invacost
Note that the InvaCost database is a dynamic resource, and its content is continuously updated. We provide here in the Data loading step a fixed snapshot of the database used for the paper, for full reproducibility.
For more information, or questions regarding the analyses, please write to leroy.boris@gmail.com, christophe.diagne@u-psud.fr and franck.courchamp@u-psud.fr
Data loading
library(invacost)
# Snapshot of the InvaCost version used in the manuscript
invacost <- getInvaCostVersion("2.1")
invacost <- invacost[-grep("NE", invacost$Cost_ID),]
invacost <- invacost[-grep("SC", invacost$Cost_ID),]
invacost <- invacost[-which(is.na(invacost$Cost_estimate_per_year_2017_USD_exchange_rate)), ]
length(which(invacost$Implementation == "Potential"))
invacost <- invacost[invacost$Implementation == "Observed", ]
length(which(invacost$Method_reliability == "Low"))
invacost <- invacost[which(invacost$Method_reliability == "High"), ]
nrow(invacost)
# Expanding and formating the database
db.over.time <- expandYearlyCosts(invacost,
startcolumn = "Probable_starting_year_adjusted",
endcolumn = "Probable_ending_year_adjusted")Figure A. Global costs
# Raw average ---------------
global.raw <- summarizeCosts(db.over.time,
minimum.year = 1970,
maximum.year = 2017)
# Cost trend over time ---------------
db.over.time$Publication_lag <- db.over.time$Publication_year - db.over.time$Impact_year
quantiles <- quantile(db.over.time$Publication_lag, probs = c(.25, .5, .75))
#
# # Assigning weights
# # Below 25% the weight does not matter because years will be removed
# year_weights[names(year_weights) >= (2017 - quantiles["25%"])] <- 0
# # Between 25 and 50%, assigning 0.25 weight
# year_weights[names(year_weights) >= (2017 - quantiles["50%"]) &
# names(year_weights) < (2017 - quantiles["25%"])] <- .25
# # Between 50 and 75%, assigning 0.5 weight
# year_weights[names(year_weights) >= (2017 - quantiles["75%"]) &
# names(year_weights) < (2017 - quantiles["50%"])] <- .5
global.trend <- modelCosts(
db.over.time, # The EXPANDED database
minimum.year = 1970,
maximum.year = 2017,
# Some years are so incomplete that we eliminate with our 25% threshold (see above)
incomplete.year.threshold = 2017 - quantiles["25%"],
incomplete.year.weights = NULL,
mars.nprune = 10)
global.raw$average.cost.per.period$middle.years <- global.raw$average.cost.per.period$initial_year +
(global.raw$average.cost.per.period$final_year -
global.raw$average.cost.per.period$initial_year) / 2
global.trend$cost.data$Calibration <- factor(global.trend$cost.data$Calibration,
levels = c("Included", "Excluded"))
p1 <- ggplot(global.raw$average.cost.per.period) +
ylab(paste0("Annual cost in US$ millions")) +
# Points
geom_point(aes_string(x = "middle.years",
y = "annual_cost")) +
# Lines between points
geom_line(aes_string(x = "middle.years",
y = "annual_cost"),
linetype = 2) +
# Horizontal bars (year span)
geom_segment(aes_string(x = "initial_year",
xend = "final_year",
y = "annual_cost",
yend = "annual_cost")) +
geom_point(data = global.trend$cost.data,
aes(x = Year, y = Annual.cost,
shape = Calibration),
alpha = .6)
p1 <- p1 +
xlab("Year") +
scale_x_continuous(breaks = global.raw$year.breaks) +
theme_bw() +
scale_y_log10(breaks = 10^(-15:15),
labels = scales::comma) +
annotation_logticks(sides = "l")
p1 <- p1 +
# stat_smooth(data = global.trend$estimated.annual.costs[-grep("Quantile", global.trend$estimated.annual.costs$model), ],
# aes(x = Year, y = fit)) + # To see the shape of model average
geom_line(data = global.trend$estimated.annual.costs[global.trend$estimated.annual.costs$model == "OLS regression" &
global.trend$estimated.annual.costs$Details == "Linear", ],
aes_string(x = "Year",
y = "fit"),
size = 1, alpha = .6) +
geom_text(data = global.raw$average.cost.per.period,
aes(x = initial_year,
y = annual_cost + 0.25 * annual_cost,
label = scales::comma(annual_cost)),
hjust = 0, size = 4) +
theme(panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14),
panel.grid.minor = element_blank(),
legend.position = c(.15, .85),
legend.text = element_text(size = 12),
legend.title = element_text(size = 12))
pdf("./figure1.pdf", height = 5, width = 8)
p1
dev.off()## png
## 2We found that biological invasions costed in US$ billion 1,288.09 (2017 US$) to human societies since 1970 (N = 1319 cost data; see Methods for details on the data processing). On average, invasions costed US$ billion 26.84 per year (Figure 1); however, the average cost during recent years was higher with US$ billion 83.31 per year between 2000 and 2009.
Raw cost for 2010-2017: 29.20
However, these figures do not account for the unsaturated increase in invasions for the past 50 years and their effects on costs. When we modelled the temporal trends of costs (Appendix global.1, global.2), we found a more-than-3-fold increase in annual costs every 10 years since 1970. (Examples for linear regressions: OLS linear regression 3.53-fold increase every 10 year; Robust linear regression, 3.07-fold increase every 10 year).
Our regressions over time predicted that the average annual cost of invasions in 2017 ranged between US$ billion 47 and 163 depending on the model, with an average value at US$ billion 93.
Table Global costs
Here are the predicted values of all models for key years:
key.years <- global.trend$estimated.annual.costs
key.years <- key.years[which(key.years$Year %in% c(1970, 1980, 1990, 2000, 2010, 2017)), ]
# key.years <- key.years[-grep("Quantile", key.years$model), ] # Omit quantile because they don't estimate average trend
key.years <- key.years[-which(key.years$model == "OLS regression" &
key.years$Details == "Quadratic"), ] # Omit OLS quad because not significant
key.years <- tidyr::spread(key.years[, 1:4], key = Year, value = fit)
key.years[, c("1970", "1980", "1990", "2000", "2010", "2017")] <- data.frame(lapply(key.years[c("1970", "1980", "1990", "2000", "2010", "2017")], scales::comma))
# Ordering models
key.years <- key.years[c(3, 7, 8, 2, 1, 4, 5, 6), ]
key.years| model | Details | 1970 | 1980 | 1990 | 2000 | 2010 | 2017 | |
|---|---|---|---|---|---|---|---|---|
| 3 | OLS regression | Linear | 245.4 | 866.0 | 3,056 | 10,783 | 38,051 | 91,980 |
| 7 | Robust regression | Linear | 239.4 | 735.6 | 2,260 | 6,945 | 21,340 | 46,823 |
| 8 | Robust regression | Quadratic | 372.6 | 678.9 | 1,674 | 5,582 | 25,182 | 86,535 |
| 2 | MARS | 725.6 | 725.6 | 1,025 | 32,548 | 18,322 | 162,719 | |
| 1 | GAM | 414.8 | 673.7 | 1,720 | 7,771 | 29,968 | 78,795 | |
| 4 | Quantile regression | Quantile 0.1 | 150.4 | 407.1 | 1,102 | 2,984 | 8,079 | 16,224 |
| 5 | Quantile regression | Quantile 0.5 | 215.4 | 707.1 | 2,321 | 7,620 | 25,017 | 57,493 |
| 6 | Quantile regression | Quantile 0.9 | 695.9 | 3,872.7 | 21,553 | 119,950 | 667,562 | 2,219,916 |
Our results also showed that there is a strong inter-annual variation in cost data, with a limited number of years having a very high economic cost, and the majority of years having below-average economic costs (Figure 1, Appendix global.2). This inter-annual heterogeneity can either reflect an actual pattern of boom and bust impacts, where invasive species have high economic impacts on specific years followed by a lower impact afterwards (i.e., systems collapse or adapt), or a lack of data for most years.
Figure B. Costs per taxonomic group
# year_weights -> tmp
year_weights <- NULL
# year_weights <- tmp
library(purrr)
library(cowplot)
Plants <- list(
Plants = db.over.time[db.over.time$Kingdom %in% "Plantae", ],
Liliopsida = db.over.time[db.over.time$Class %in% c("Liliosida",
"Liliopsida"), ],
Magnoliopsida = db.over.time[db.over.time$Class %in% "Magnoliopsida", ],
Pinopsida = db.over.time[db.over.time$Class %in% "Pinopsida", ],
# poly = db.over.time[db.over.time$Class %in% "Polypodiopsida", ], # Cost really small
Ulvophyceae = db.over.time[db.over.time$Class %in% "Ulvophyceae", ]
)
Invertebrates <-list(
Invertebrates = db.over.time[db.over.time$Phylum %in%
c("Arthropoda",
"Ctenophora",
"Mollusca",
"Mollusca/Mollusca",
"Nematoda"), ],
Insecta = db.over.time[db.over.time$Class %in% "Insecta", ],
Bivalvia = db.over.time[db.over.time$Class %in% "Bivalvia", ],
Gastropoda = db.over.time[db.over.time$Class %in% "Gastropoda", ],
Secernentea = db.over.time[db.over.time$Class %in% "Secernentea", ]
)
Vertebrates <- list(
Vertebrates = db.over.time[db.over.time$Phylum %in% "Chordata", ],
Amphibia = db.over.time[db.over.time$Class %in% "Amphibia", ],
Reptilia = db.over.time[db.over.time$Class %in% "Reptilia", ],
Aves = db.over.time[db.over.time$Class %in% "Aves", ],
Mammalia = db.over.time[db.over.time$Class %in% "Mammalia", ]
)
Mixed <- list(
Mixed = db.over.time[which(
db.over.time$Kingdom %in% c("Animalia/Plantae", "Diverse/Unspecified") |
!(db.over.time$Kingdom %in% "Plantae") &
db.over.time$Phylum %in% c("Diverse/Unspecified")), ]
)
obs.impact <- list()
pred.impact <- list()
groups <- c("Plants", "Invertebrates", "Vertebrates", "Mixed")
summary.costs <- data.frame()
for (group in groups)
{
obs.impact[[group]] <- map(get(group),
summarizeCosts,
minimum.year = 1970,
maximum.year = 2017)
pred.impact[[group]] <-
map(names(get(group)),
function(x, db, raw, ...)
if(length(which(is.na(raw[[x]]$average.cost.per.period$annual_cost))) <= 1) # Only fit models for taxa for which we have data for at least 4 decades
{
modelCosts(db[[x]],
...)
} else
{
NA
},
db = get(group),
raw = obs.impact[[group]],
minimum.year = 1970,
maximum.year = 2017,
final.year = 2017,
# Some years are so incomplete that we eliminate with our 25% threshold (see above)
incomplete.year.threshold = 2017 - quantiles["25%"],
# For the other incomplete years we apply the vector of weights that we defined above
incomplete.year.weights = year_weights
)
names(pred.impact[[group]]) <- names(get(group))
for (taxon in names(get(group)))
{
summary.costs <- rbind.data.frame(
summary.costs,
data.frame(
Group = group,
Taxon = taxon,
Type.cost = "Observed",
Cost.period = paste0(obs.impact[[group]][[taxon]]$average.total.cost$initial_year,
" - ", obs.impact[[group]][[taxon]]$average.total.cost$final_year),
Average.annual.cost = obs.impact[[group]][[taxon]]$average.total.cost$annual_cost,
Signif = NA,
Cumulated.cost = obs.impact[[group]][[taxon]]$average.total.cost$total_cost,
lwr = NA,
upr = NA))
summary.costs <- rbind.data.frame(
summary.costs,
data.frame(
Group = group,
Taxon = taxon,
Type.cost = "Observed",
Cost.period = paste0(obs.impact[[group]][[taxon]]$average.cost.per.period$initial_year[
obs.impact[[group]][[taxon]]$average.cost.per.period$initial_year %in% c(2000, 2010)
],
" - ", obs.impact[[group]][[taxon]]$average.cost.per.period$final_year[
obs.impact[[group]][[taxon]]$average.cost.per.period$initial_year %in% c(2000, 2010)
]),
Average.annual.cost = obs.impact[[group]][[taxon]]$average.cost.per.period$annual_cost[
obs.impact[[group]][[taxon]]$average.cost.per.period$initial_year %in% c(2000, 2010)
],
Signif = NA,
Cumulated.cost = obs.impact[[group]][[taxon]]$average.cost.per.period$total_cost[
obs.impact[[group]][[taxon]]$average.cost.per.period$initial_year %in% c(2000, 2010)
],
lwr = NA,
upr = NA))
if(all(is.na(pred.impact[[group]][[taxon]])))
{
estimated.costs <- data.frame(Group = group,
Taxon = taxon,
Type.cost = "Estimated",
Cost.period = NA,
Average.annual.cost = NA,
Signif = NA,
Cumulated.cost = NA,
lwr = NA,
upr = NA)
} else
{
estimated.costs <- data.frame(
Group = group,
Taxon = taxon,
Type.cost = "Estimated",
Cost.period = as.character(pred.impact[[group]][[taxon]]$parameters$final.year),
Average.annual.cost = pred.impact[[group]][[taxon]]$final.year.cost["robust.linear"],
Signif = pred.impact[[group]][[taxon]]$model.summary$robust.linear$coefficients["Year", "Pr(>|t|)"],
Cumulated.cost = NA,
lwr = pred.impact[[group]][[taxon]]$estimated.annual.costs$lwr[
pred.impact[[group]][[taxon]]$estimated.annual.costs$Year ==
pred.impact[[group]][[taxon]]$parameters$final.year &
pred.impact[[group]][[taxon]]$estimated.annual.costs$model == "Robust regression" &
pred.impact[[group]][[taxon]]$estimated.annual.costs$Details == "Linear"
],
upr = pred.impact[[group]][[taxon]]$estimated.annual.costs$upr[
pred.impact[[group]][[taxon]]$estimated.annual.costs$Year ==
pred.impact[[group]][[taxon]]$parameters$final.year &
pred.impact[[group]][[taxon]]$estimated.annual.costs$model == "Robust regression" &
pred.impact[[group]][[taxon]]$estimated.annual.costs$Details == "Linear"
])
}
summary.costs <- rbind.data.frame(
summary.costs,
estimated.costs)
}
}
# Removing estimates for which model parameters are not significant
summary.costs[which(summary.costs$Signif > 0.05), c("Average.annual.cost", "lwr", "upr")] <- NA
# summary.costs[which(summary.costs$Signif <= 0.05), ]
# Workaround for error bars
summary.costs$lwr <- ifelse(is.na(summary.costs$lwr),
summary.costs$Average.annual.cost,
summary.costs$lwr)
summary.costs$upr <- ifelse(is.na(summary.costs$upr),
summary.costs$Average.annual.cost,
summary.costs$upr)
# Removing gasteropods and mixed data because of low data quality
summary.costs$Average.annual.cost[
which(summary.costs$Taxon == "Gastropoda" &
summary.costs$Type.cost == "Estimated")
] <- NA
summary.costs$Average.annual.cost[
which(summary.costs$Taxon == "Bivalvia" &
summary.costs$Type.cost == "Estimated")
] <- NA
summary.costs$Average.annual.cost[
which(summary.costs$Taxon == "Mixed" &
summary.costs$Type.cost == "Estimated")
] <- NA
summary.costs$Taxon[summary.costs$Taxon == "Mixed"] <- "Multi-taxa\nstudies"
summary.costs$Group[summary.costs$Group == "Mixed"] <- "Mixed\ntaxa"
# Removing multi-group studies
summary.costs <- summary.costs[-which(summary.costs$Group %in% "Mixed\ntaxa"), ]
# Changing order of facets
summary.costs$Group <- factor(as.character(summary.costs$Group),
levels = c("Invertebrates", "Vertebrates", "Plants"))
# Adding letters to facets for the plot
summary.costs$Group.a <- summary.costs$Group.b <- summary.costs$Group
levels(summary.costs$Group.a) <- c("a. invertebrates", "b. vertebrates", "c. plants")
levels(summary.costs$Group.b) <- c("d. invertebrates", "e. vertebrates", "f. plants")
# Preparing x-axis positions
summary.costs$Taxon <- factor(summary.costs$Taxon, levels = unique(summary.costs$Taxon))
summary.costs$x.pos <- as.numeric(summary.costs$Taxon)
summary.costs$x.pos[which(!(summary.costs$x.pos %in% c(1, 6, 11)))] <- summary.costs$x.pos[which(!(summary.costs$x.pos %in% c(1, 6, 11)))] + .5
summary.costs$x.pos[which(summary.costs$x.pos >= 11)] <- summary.costs$x.pos[which(summary.costs$x.pos >= 11)] + 2
summary.costs$x.pos[which(summary.costs$x.pos >= 6)] <- summary.costs$x.pos[which(summary.costs$x.pos >= 6)] + 2
levels(summary.costs$Taxon)[which(levels(summary.costs$Taxon) %in% c("Plants", "Invertebrates", "Vertebrates"))] <- c("All plants", "All\ninvertebrates", "All vertebrates")
# summary.costs$Cost.period <- as.factor(summary.costs$Cost.period)
# Removing the 2010-2017 period because it does not provide a lot of information
pavg <- ggplot(summary.costs[which(summary.costs$Cost.period != "2010 - 2017"), ]) +
geom_pointrange(aes(x = x.pos, y = Average.annual.cost,
col = Cost.period,
shape = Type.cost,
ymin = lwr,
ymax = upr),
position = position_dodge(width = .5),
size = .5) +
scale_x_continuous(breaks = unique(summary.costs$x.pos),
minor_breaks = NULL,
labels = unique(summary.costs$Taxon)) +
theme_bw() +
scale_y_log10(breaks = 10^(-15:15),
labels = scales::comma,
limits = c(0.09, NA)) +
annotation_logticks(sides = "l") +
facet_grid(~ Group.b, scales = "free_x", space = "free_x") +
ylab("Average annual\ncost in\nUs$ millions") +
scale_color_manual(na.translate=FALSE,
values = c(`1970 - 2017` = rgb(86, 180, 233, maxColorValue = 255), # Colours for colourblind people
`2000 - 2009` = rgb(213, 94, 0, maxColorValue = 255),
`2017` = rgb(204, 121, 167, maxColorValue = 255))) +
scale_shape_discrete(labels = c("Observed average cost",
"Estimated cost ± 95%IC")) + # Signe +-
guides(col = guide_legend(title = "Time interval",
override.aes = list(linetype = c(0, 0, 0),
shape = c(15, 15, 15),
size = 1.5),
ncol = 1, order = 1),
shape = guide_legend(title = "Cost nature",
override.aes = list(shape = c(17, 16),
linetype = c(0, 1)),
ncol = 1), order = 2) +
theme(legend.position = c(-.15, 1.47),
axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.25,
size = 13, face = c("bold", rep("plain", 4))),
axis.text.y = element_text(size = 13),
plot.margin = margin(t = 0, .2, 0, 3, "lines"),
legend.box = "vertical",
plot.title = element_text(hjust = 0),
panel.grid.major.x = element_blank(),
panel.grid.minor.y = element_blank(),
axis.title.x = element_text(size = 17),
axis.title.y = element_text(size = 17),
strip.text.x = element_text(size = 13)) +
# ggtitle("Average annual cost") +
xlab("Taxon")
# Trick to get colours & bars stacked correctly
summary.costs2 <- summary.costs[nrow(summary.costs):1, ]
# summary.costs2$Cost.period <-
# factor(summary.costs2$Cost.period,
# levels = rev(levels(summary.costs2$Cost.period)))
ptot <- ggplot(summary.costs[which(summary.costs$Cost.period != "2010 - 2017" &
!is.na(summary.costs$Cumulated.cost)), ]) +
geom_col(aes(x = x.pos, y = Cumulated.cost,
fill = Cost.period),
width = c(.9, .75),
position = position_dodge(width = 0.1,
preserve = "total")) +
theme_bw() +
scale_fill_manual(na.translate=FALSE,
values = c(rgb(86, 180, 233, maxColorValue = 255), # Same colours as above
rgb(213, 94, 0, maxColorValue = 255))) +
scale_y_continuous(breaks = seq(0, 2e6, by = 1e5),
labels = scales::comma) +
scale_x_continuous(breaks = unique(summary.costs2$x.pos),
minor_breaks = NULL,
labels = unique(summary.costs2$Taxon)) +
# scale_y_log10(breaks = 10^(-15:15),
# labels = scales::comma,
# limits = c(1, NA)) +
# annotation_logticks(sides = "l") +
facet_grid(~ Group.a, scales = "free_x", space = "free_x") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0,
size = 13,
face = c(rep("plain", 4), "bold")), # We reverse here because rows are reversed!
axis.text.y = element_text(size = 13),
plot.margin = margin(t = 0, 0.2, 0, 3, "lines"),
plot.title = element_text(hjust = 0),
panel.grid.major.x = element_blank(),
panel.grid.minor.y = element_blank(),
axis.title.x = element_text(size = 17),
axis.title.y = element_text(size = 17),
strip.text.x = element_text(size = 13)) +
ylab("Cumulative\ncost in\nUS$ millions") +
guides(fill = FALSE) +
# ggtitle("Cumulative cost") +
xlab("")
# png("./Figure3.png", units = "in", height = 8.5, width = 10,
# res = 300)
pdf("./figure3.pdf", height = 8.5, width = 10)
plot_grid(ptot, pavg, ncol = 1, align = "v")
dev.off()Methodological details
This figure represents (a) the average annual costs as observed in the database (1970 - 2017 and 2000-2009) and as predicted by the linear regression over time for 2017 for those taxons for which we had enough data points (i.e., more than 30 years of data); and (b) the cumulated cost over time for 1970-2017 and 2000-2009. Note that these values include only estimates that could be derived for a single taxonomic group; hence, a significant proportion of the economic impacts of invasive species could not be split per taxonomic group (see Appendix Taxo.1). The 1970-2017 period was chosen as the range of time periods in the database for which we have relatively complete data. The 2000-2009 time period was chosen as the 10-year interval for which we have the most complete data and the highest economic impacts of IAS.
The average annual costs for 1970-2017 and 2000-2009 are represented without error bars, following Weissgerber et al. (2014) for two reasons. First, there are not enough data points for error bars to be meaningful; second, the distribution of data points is skewed, with most years having a lower-than-average economic cost, and only a few years having a high economic cost. Therefore, we refer readers to Appendix Taxo.2 where we illustrate for each individual taxon the average economic costs over the full distribution of economic impacts per year. In spite of the skewed distribution of data chose to illustrate average annual costs rather than median annual costs, because we assume that the skewness of data is caused by the strong incompleteness of economic data for most years (see our rationale in …[main text or appendix?]); therefore, we deemed that the average annual cost is probably closer to the actual annual cost than the median annual cost. Regardless, we provide the median anual cost for each taxon in Appendix XX by the means of 0.5 quantile regressions.
Results
Out of the global cumulated cost of invasive species, US$ billion 698 (54.2%) were derived from estimates combining multiple taxonomic groups, and US$ billion 591 (45.8%) were derived from estimates that could be assigned to a single taxonomic group (Appendix Taxo.1). Update values if database change! Because they are derived from the appendix below!
Among these single-groupe estimates, we found that the single costliest group was invertebrates, almost exclusively because of insects, with a cumulated economic cost of US$ billion 416 (Figure 2a). Vertebrates had the second highest impact with US$ billion 166; the majority of this impact was caused by mammals. Plants had the lowest reported impact with US$ billion 9.
Likewise, the observed and estimated annual costs of invasive invertebrates were an order of magnitude higher than those of invertebrates and plants (Figure 2b). Worryingly, the predicted average annual cost for 2017 was almost always higher than the observed average annual cost for previous periods of time. In other words, for all sufficiently-known groups, the annual cost of invasions is continuously increasing (Figure 2b, Appendix Taxo.2). The only exception to this pattern is vertebrates, for which extremely high costs were reported in the period 2000-2009, which resulted in an average annual cost for this period much higher than predicted by a linear trend over 1970-2017.
Discuss whether or not we go into details for subgroups
Figure C. Geographic distribution of costs
Figure C.1 Cost per REGION
db.over.time$Geographic_region[which(db.over.time$Geographic_region == "Europa")] <- "Europe"
library(purrr)
library(tmap)
Regions <- list(
Africa = db.over.time[db.over.time$Geographic_region %in% "Africa", ],
Asia = db.over.time[db.over.time$Geographic_region %in% "Asia", ],
Europe = db.over.time[db.over.time$Geographic_region %in% "Europe", ],
`Central America` = db.over.time[db.over.time$Geographic_region %in% "Central America", ],
`North America` = db.over.time[db.over.time$Geographic_region %in% "North America", ],
`Oceania` = db.over.time[db.over.time$Geographic_region %in% "Oceania", ],
`South America` = db.over.time[db.over.time$Geographic_region %in% "South America", ]
)
obs.regions <- map(Regions,
summarizeCosts,
minimum.year = 1970,
maximum.year = 2017)
library(sf)
library(sp)
library(rgdal)
library(rgeos)
# Download the data here from natural earth data, 10m resolution
# https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_admin_0_map_units.zip
# Map units, admin 0: https://www.naturalearthdata.com/downloads/10m-cultural-vectors/
invacost.world <- readOGR("d:/r/projects/invacost_nature/data/ne_10m_admin_0_map_units.shp")## OGR data source with driver: ESRI Shapefile
## Source: "D:\R\Projects\invacost_nature\data\ne_10m_admin_0_map_units.shp", layer: "ne_10m_admin_0_map_units"
## with 295 features
## It has 94 fields
## Integer64 fields read as strings: NE_IDinvacost.world$invacost.region <- as.character(invacost.world$CONTINENT)
invacost.world$invacost.region[invacost.world$SUBREGION %in% c("Central America",
"Caribbean")] <- "Central America"
invacost.world <- invacost.world[-which(invacost.world$invacost.region %in% c("Seven seas (open ocean)", "Antarctica")), ]
invacost.world2 <- do.call(raster::bind,
sapply(levels(as.factor(invacost.world$invacost.region)),
function(x) gUnaryUnion(invacost.world[which(invacost.world$invacost.region == x), ]),
USE.NAMES = FALSE))
invacost.world <- SpatialPolygonsDataFrame(invacost.world2,
data = data.frame(Region = levels(as.factor(invacost.world$invacost.region))),
match.ID = FALSE)
rm(invacost.world2)
invacost.world <- st_as_sf(invacost.world)
invacost.world$cumulated.cost <-
invacost.world$total.cost.2000.2009 <-
invacost.world$average.cost.2000.2009 <-
invacost.world$estimated.annual.cost.2017 <- NA
for(region in names(Regions))
{
invacost.world$cumulated.cost[invacost.world$Region == region] <-
obs.regions[[region]]$average.total.cost$total_cost
invacost.world$total.cost.2000.2009[invacost.world$Region == region] <-
obs.regions[[region]]$average.cost.per.period[
which(obs.regions[[region]]$average.cost.per.period$initial_year == 2000), "total_cost"]
invacost.world$average.cost.2000.2009[invacost.world$Region == region] <-
obs.regions[[region]]$average.cost.per.period[
which(obs.regions[[region]]$average.cost.per.period$initial_year == 2000), "annual_cost"]
# invacost.world$estimated.annual.cost.2017[invacost.world$Region == region] <-
# pred.regions[[region]]$final.year.cost["linear"]
}
tm_shape(invacost.world,
projection = "+proj=eck4") +
tm_polygons("cumulated.cost",
palette = "OrRd",
title = "Cost in US$ millions",
style="pretty") +
tm_layout(# title = "Average annual cost from 2000 to 2009",
# title.size = .9,
legend.bg.color = "white",
# title.position=c("center", "top"),
inner.margins=c(.04, .01, .1, .01),
bg.color="#AEDFE5",
outer.bg.color="white",
earth.boundary=c(-180, 180, -70, 90),
earth.boundary.color="white",
earth.boundary.lwd=.4,
space.color="white",
attr.outside=T,
attr.color="grey20",
frame = FALSE,
legend.title.size = 1.5,
legend.text.size = 1)
As expected, the average annual cost between 2000 and 2009 was highest in North America (Figure 3), with a total of US$ millions48,093 (figure 3, appendix region.2). This pattern probably reflects a bias in the literature. North America had the second highest number of unique cost estimates (304 - the first is Oceania-Pacific islands with 490 unique estimates - Appendix region.1) and the highest combined total number of years over which values were estimated (1062). Second, Asia had the second largest average annual cost of IAS, with US$ millions20,862. Third, the average annual cost was between US$ million 370 and 7,991. This pattern is similar when looking at the distribution of cumulated costs for the entire period (Appendix region.3).
Figure C.2 Cost per country
# SAME AS ABOVE, Download the data here from natural earth data, 10m resolution
# https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_admin_0_map_units.zip
# Map units, admin 0: https://www.naturalearthdata.com/downloads/10m-cultural-vectors/
countries <- readOGR("d:/r/projects/invacost_nature/data/ne_10m_admin_0_map_units.shp")## OGR data source with driver: ESRI Shapefile
## Source: "D:\R\Projects\invacost_nature\data\ne_10m_admin_0_map_units.shp", layer: "ne_10m_admin_0_map_units"
## with 295 features
## It has 94 fields
## Integer64 fields read as strings: NE_IDinvacost.countries <- data.frame(original.name = unique(invacost$Official_country))
invacost.countries$valid <- FALSE
invacost.countries$valid[which(invacost.countries$original.name %in% countries$ADMIN)] <- TRUE
invacost.countries$renamed <- as.character(invacost.countries$original.name)
# Corrections to country names
invacost.countries$original.name[!invacost.countries$valid]## [1] "USA"
## [2] "Scotland"
## [3] "Diverse/Unspecified"
## [4] "Great Britain"
## [5] "Antigua"
## [6] "Argentina/Chile"
## [7] "Malaysia/Myanmar/Thailand"
## [8] "Philippines/Thailand/Vietnam"
## [9] "Australia/Brazil/Venezuela"
## [10] "England"
## [11] "Northern Ireland"
## [12] "Austria/Denmark/Finland/Ireland/Latvia/Lithuania/Netherlands/Spain/Sweden/UK"
## [13] "Mexico/USA"
## [14] "Belize/Guatemala/Mexico/USA"
## [15] "Australia/Papua New Guinea"
## [16] "India/Pakistan"
## [17] "Guatemala/Mexico/USA"
## [18] "Tanzania"
## [19] "Benin/Guinea Conakry/Togo"
## [20] "Canada/USA"
## [21] "Germany/Unspecified"
## [22] "Portugal/Spain"
## [23] "USA/Unspecified"
## [24] "Timor-Leste"
## [25] "Brazil/Guatemala/Panama/Salvador/Venezuela"
## [26] "Cambodia/Malaysia/Thailand"
## [27] "Salvador"
## [28] "Bangladesh/India/Pakistan/Sri Lanka"invacost.countries$renamed[invacost.countries$original.name %in% "USA"] <- "United States of America"
invacost.countries$renamed[invacost.countries$original.name %in% c("Scotland", "Great Britain", "England", "Northern Ireland")] <- "United Kingdom"
invacost.countries$renamed[invacost.countries$original.name %in% "Antigua"] <- "Antigua and Barbuda"
invacost.countries$renamed[invacost.countries$original.name %in% "Bahamas"] <- "The Bahamas"
invacost.countries$renamed[invacost.countries$original.name %in% "Tanzania"] <- "United Republic of Tanzania"
invacost.countries$renamed[invacost.countries$original.name %in% "Timor-Leste"] <- "East Timor"
invacost.countries$renamed[invacost.countries$original.name %in% "Salvador"] <- "El Salvador"
# Omitted countries
# invacost.countries$renamed[!invacost.countries$renamed %in% countries$ADMIN]
db.over.time$country2 <- invacost.countries$renamed[match(db.over.time$Official_country, invacost.countries$original.name)]
country_summary <- data.frame()
for(country in unique(db.over.time$country2))
{
subdb <- db.over.time[which(db.over.time$country2 %in% country), ]
curestimate <- summarizeCosts(subdb,
minimum.year = 1970,
maximum.year = 2017)
country_summary <- rbind.data.frame(country_summary,
data.frame(country = country,
nb.estimates = curestimate$average.total.cost$number_estimates,
cumulated.cost = curestimate$average.total.cost$total_cost,
average.annual.cost = curestimate$average.total.cost$annual_cost,
average.annual.cost.2000.2009 = curestimate$average.cost.per.period$annual_cost[
which(curestimate$average.cost.per.period$initial_year == 2000)
]))
}
countries2 <- countries[countries$ADMIN %in% country_summary$country, ]
countries2$cumulated.cost <- country_summary$cumulated.cost[match(countries2$ADMIN,
country_summary$country)]
countries2$nb.estimates <- country_summary$nb.estimates[match(countries2$ADMIN,
country_summary$country)]
countries2sf <- st_as_sf(countries2)
# Countries with multiple polygons have to be corrected
countries2sf$ADMIN[which(duplicated(countries2sf$ADMIN))]## [1] "France" "Belgium"
## [3] "South Korea" "United Kingdom"
## [5] "United Kingdom" "United Kingdom"
## [7] "Belgium" "Netherlands"
## [9] "France" "France"
## [11] "Antigua and Barbuda" "France"
## [13] "France" "United Republic of Tanzania"
## [15] "Portugal" "Portugal"
## [17] "China" "New Zealand"countries2sf$BRK_NAME[which(duplicated(countries2sf$ADMIN))]## [1] "France" "Walloon" "Republic of Korea"
## [4] "Wales" "England" "Scotland"
## [7] "Brussels" "Caribbean Netherlands" "Martinique"
## [10] "Guadeloupe" "Barbuda" "Réunion"
## [13] "Mayotte" "Zanzibar" "Madeira"
## [16] "Azores" "Paracel Is." "Tokelau"countries2sf$BRK_NAME[which(countries2sf$ADMIN == "New Zealand")]## [1] "New Zealand" "Tokelau"# countries2sf$nb.estimates[which(duplicated(countries2sf$ADMIN))] <- NA
# Correction to center the UK point on England
countries2sf$nb.estimates[which(countries2sf$ADMIN == "United Kingdom")][c(1, 2, 4)] <- NA
countries2sf$nb.estimates[which(countries2sf$ADMIN == "France")][c(1, 3, 4, 5, 6)] <- NA
countries2sf$nb.estimates[which(countries2sf$ADMIN == "South Korea")][1] <- NA
countries2sf$nb.estimates[which(countries2sf$ADMIN == "Belgium")][c(1, 3)] <- NA
countries2sf$nb.estimates[which(countries2sf$ADMIN == "Netherlands")][2] <- NA
countries2sf$nb.estimates[which(countries2sf$ADMIN == "United Republic of Tanzania")][2] <- NA
countries2sf$nb.estimates[which(countries2sf$ADMIN == "Portugal")][c(2, 3)] <- NA
countries2sf$nb.estimates[which(countries2sf$ADMIN == "China")][2] <- NA
countries2sf$nb.estimates[which(countries2sf$ADMIN == "New Zealand")][2] <- NA
fig4 <- tm_shape(invacost.world,
projection = "+proj=eck4") +
tm_polygons() +
tm_shape(countries2sf,
projection="eck4") +
tm_polygons("cumulated.cost",
palette = "OrRd",
title = "Cost in US$ millions",
style = "fixed",
breaks = c(0, 15000, 50000, 150000, 350000, 550000)) +
tm_layout(legend.outside = TRUE) +
tm_symbols(size = "nb.estimates",
title.size = "Number of cost estimates\nper country",
style = "fixed",
sizes.legend = c(10, 50, 100, 150, 250, 350),
perceptual = TRUE,
scale = 4,
alpha = .5,
col = rgb(0, 114, 178, maxColorValue = 255)
) +
tm_layout(# title = "Average annual cost from 2000 to 2009",
# title.size = .9,
# legend.bg.color = "white",
# title.position=c("center", "top"),
inner.margins=c(.01, .01, .05, .01),
outer.margins = c(.01, .01, .01, .01) ,
between.margin = 0.01,
# bg.color="#AEDFE5",
# outer.bg.color="white",
earth.boundary=c(-180, 180, -70, 90),
earth.boundary.color="white",
earth.boundary.lwd=.4,
# space.color="white",
attr.outside=T,
attr.color="grey20",
frame = FALSE,
legend.title.size = 1.5,
legend.text.size = 1,
legend.width = 0.1,
legend.outside.size = 0.2)
cairo_pdf("figure4.pdf", width = 16, height = 6.5)
fig4
dev.off()## png
## 2fig4
Figure D. Types of costs
Figure Da. Cumulative costs over time
library(purrr)
plyr::count(invacost$Type_of_cost_merged)
costtype <- list(
Damage = db.over.time[db.over.time$Type_of_cost_merged %in% "Damage_costs", ],
Management = db.over.time[db.over.time$Type_of_cost_merged %in% "Management_costs", ]
)
obs.costtype <- map(costtype,
summarizeCosts,
minimum.year = 1970,
maximum.year = 2017)
cost.data <- rbind(
data.frame(obs.costtype$Damage$cost.data,
type = "Damage"),
data.frame(obs.costtype$Management$cost.data,
type = "Management")
)
costperiod <- rbind(
data.frame(obs.costtype$Damage$average.cost.per.period,
type = "Damage"),
data.frame(obs.costtype$Management$average.cost.per.period,
type = "Management")
)
costperiod$middle.years <- costperiod$initial_year +
(costperiod$final_year -
costperiod$initial_year) / 2
plot.breaks = 10^(-15:15)
yeargroups.damage <- dplyr::group_by(obs.costtype$Damage$cost.data,
get(obs.costtype$Damage$parameters$year.column))
yeargroups.management <- dplyr::group_by(obs.costtype$Management$cost.data,
get(obs.costtype$Management$parameters$year.column))
yearly.cost <- rbind.data.frame(
setNames(data.frame(dplyr::summarise(yeargroups.damage,
Annual.cost = sum(get(obs.costtype$Damage$parameters$cost.column))),
"Damage"), c("Year", "Annual.cost", "type")),
setNames(data.frame(dplyr::summarise(yeargroups.management,
Annual.cost = sum(get(obs.costtype$Management$parameters$cost.column))),
"Management"), c("Year", "Annual.cost", "type"))
)
ggplot(costperiod) +
ylab("Average annual cost per period in US$ millions") +
xlab("Year") +
scale_x_continuous(breaks = obs.costtype$Damage$year.breaks) +
theme_bw() +
scale_y_log10(breaks = plot.breaks,
labels = scales::comma) +
annotation_logticks() +
# Points
geom_point(aes_string(x = "middle.years",
y = "annual_cost",
col = "type"),
size = 3) +
# Lines between points
geom_line(aes_string(x = "middle.years",
y = "annual_cost",
col = "type"),
linetype = 2,
size = 1) +
# Horizontal bars (year span)
geom_segment(aes_string(x = "initial_year",
xend = "final_year",
y = "annual_cost",
yend = "annual_cost",
col = "type"),
size = 1) +
# Points
geom_point(data = yearly.cost,
mapping = aes_string(x = "Year",
y = "Annual.cost",
col = "type"),
size = 2,
alpha = .8) +
guides(col = guide_legend(title = "Category\nof cost")) +
theme(panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.text = element_text(size = 13),
axis.title = element_text(size = 18),
panel.grid.minor = element_blank(),
legend.text = element_text(size = 13),
legend.title = element_text(size = 14),
legend.position = c(.15, .9))
We found different temporal patterns between damage and management costs (Figure 4a). Damage costs have strongly and continuously increased since 1970, reaching an annual average of US$ millions 56,212 in 2000-2009. However, management costs only increased during the 1990s, and are 10-times smaller than damage costs, with ran annual average of US$ millions 3,019 in 2000-2009. Interestingly, the trend over time is extremely similar between damage and management costs, except that management costs started to increase approximately 15 years after damage costs (Figure 4.b ou appendix au choix), especially when we estimate trends exclusively on the basis of the the most complete years (i.e., years that have at least 75% completeness, all years before 2006).
Figure Db. Trends over time based on linear robust regressions
pred.costtype.25 <- map(costtype,
modelCosts,
minimum.year = 1970,
maximum.year = 2017,
final.year = 2017,
# Some years are so incomplete that we eliminate with our 25% threshold (see above)
incomplete.year.threshold = 2017 - quantiles["25%"],
# For the other incomplete years we apply the vector of weights that we defined above
incomplete.year.weights = year_weights)
cost.data <- rbind(
data.frame(pred.costtype.25$Damage$cost.data,
type = "Damage"),
data.frame(pred.costtype.25$Management$cost.data,
type = "Management")
)
cost.data$Calibration <- factor(cost.data$Calibration,
levels = c("Included", "Excluded"))
model.preds <- rbind(
data.frame(pred.costtype.25$Damage$estimated.annual.costs,
type = "Damage",
calib = 2017 - quantiles["25%"]),
data.frame(pred.costtype.25$Management$estimated.annual.costs,
type = "Management",
calib = 2017 - quantiles["25%"])
)
model.preds$calib <- as.factor(model.preds$calib)
# model.preds <- model.preds[which(model.preds$model == "GAM"), ]
model.preds <- model.preds[which(model.preds$model == "Robust regression" &
model.preds$Details == "Linear"), ]
fig2 <- ggplot() +
ylab("Annual cost in US$ millions") +
xlab("Year") +
theme_bw() +
scale_y_log10(breaks = plot.breaks,
labels = scales::comma) +
annotation_logticks() +
geom_point(data = cost.data,
aes_string(x = "Year",
y = "Annual.cost",
col = "type",
shape = "Calibration"),
size = 2, alpha = .8) +
geom_line(data = model.preds[which(model.preds$model == "Robust regression" &
model.preds$Details == "Linear"), ],
aes_string(x = "Year",
y = "fit",
col = "type"),
size = 1.1,
alpha = .8) +
geom_ribbon(data = model.preds[which(model.preds$model == "Robust regression" &
model.preds$Details == "Linear"), ],
aes_string(x = "Year",
ymin = "lwr",
ymax = "upr",
fill = "type"),
alpha = .1,
linetype = 0) +
scale_color_discrete(name = "Type of cost") +
theme() +
guides(col = guide_legend(title = "Category of cost"),
fill = FALSE) +
theme(panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.text = element_text(size = 13),
axis.title = element_text(size = 18),
panel.grid.minor = element_blank(),
legend.text = element_text(size = 12,
margin = margin(t = 5, b = 5, unit = "pt")),
legend.title = element_text(size = 12),
legend.position = c(.2, .75))
pdf("./figure2.pdf", height = 5, width = 8)
fig2
dev.off()
fig2
Increases in costs per 10 year periods:
Damage: 5.95
Management: 1.75
Figure Dc. Attempt to merge figure 4a & b in a single one
gglines <- rbind.data.frame(data.frame(Year = model.preds$Year,
Annual_cost = model.preds$fit,
Data.type = "Predicted trend over time (Robust regression)",
Cost.category = model.preds$type,
upr = model.preds$upr,
lwr = model.preds$lwr),
data.frame(Year = costperiod$middle.years,
Annual_cost = costperiod$annual_cost,
Data.type = "Average cost among 10-year periods",
Cost.category = costperiod$type,
upr = NA,
lwr = NA))
gglines$Data.type <- factor(gglines$Data.type,
levels = c("Predicted trend over time (Robust regression)",
"Average cost among 10-year periods"))
ggplot(costperiod) +
ylab("Average annual cost per period in US$ millions") +
xlab("Year") +
scale_x_continuous(breaks = obs.costtype$Damage$year.breaks) +
theme_bw() +
scale_y_log10(breaks = plot.breaks,
labels = scales::comma) +
annotation_logticks(sides = "l") +
geom_ribbon(data = gglines[-which(is.na(gglines$upr)), ], aes_string(x = "Year",
ymin = "lwr",
ymax = "upr",
fill = "Cost.category"),
alpha = .1, linetype = 0) +
# Lines between points
geom_line(data = gglines,
aes_string(x = "Year",
y = "Annual_cost",
col = "Cost.category",
linetype = "Data.type"),
size = 1) +
# Horizontal bars (year span) +
geom_segment(aes_string(x = "initial_year",
xend = "final_year",
y = "annual_cost",
yend = "annual_cost",
col = "type"),
linetype = 2,
size = 1) +
# Points
geom_point(data = yearly.cost,
mapping = aes_string(x = "Year",
y = "Annual.cost",
col = "type"),
size = 2,
alpha = .8) +
guides(col = guide_legend(title = "Category of cost"),
linetype = guide_legend(title = "Line type",
reverse = TRUE),
fill = FALSE) +
theme(panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.text = element_text(size = 13),
axis.title = element_text(size = 18),
panel.grid.minor = element_blank(),
legend.text = element_text(size = 13),
legend.title = element_text(size = 14),
legend.position = c(.04, .999),
legend.justification = c(0, 1))
Figure Dd. Alluvial diagram of management vs damage costs
library(ggalluvial)
ggplot(data = yearly.cost,
aes(x = Year, y = Annual.cost, alluvium = type)) +
geom_alluvium(aes(fill = type, colour = type),
alpha = .75, decreasing = NA) +
theme_bw() +
theme(axis.text.x = element_text(angle = -30, hjust = 0)) +
scale_fill_brewer(type = "qual", palette = "Set3") +
scale_color_brewer(type = "qual", palette = "Set3") +
theme_bw() +
scale_y_log10(breaks = plot.breaks,
labels = scales::comma,
limits = c(1, 1e9)) +
annotation_logticks(sides = "l") 
Figure E Top 10 species
library(png)
library(RCurl)
library(ggtext)
# Lumping together all species that are frequently found together or confused
db.over.time$sp.list <- as.character(db.over.time$Species)
db.over.time$genus.list <- as.character(db.over.time$Genus)
db.over.time$sp.list[which(db.over.time$Genus == "Aedes")] <- "Aedes spp."
db.over.time$sp.list[which(db.over.time$Genus == "Rattus" | db.over.time$Genus == "Mus")] <- "Rattus spp./Mus spp."
db.over.time$sp.list[which(db.over.time$Species %in% c("Rattus sp./Mus sp.",
"Mus musculus/Rattus rattus",
"Mus musculus/Rattus norvegicus",
"Mus sp./Rattus sp."))] <- "Rattus spp./Mus spp."
db.over.time$genus.list[which(db.over.time$sp.list == "Rattus spp./Mus spp.")] <- "Rattus/Mus"
db.over.time$sp.list[which(db.over.time$Genus == "Felis/Rattus")] <- "Felis catus/Rattus spp."
db.over.time$sp.list[which(db.over.time$Genus == "Oryctolagus/Rattus")] <- "Oryctolagus spp./Rattus spp."
db.over.time$sp.list[which(db.over.time$Genus == "Canis")] <- "Canis lupus spp."
# Unique identifier
db.over.time$unique.sp.id <- do.call("paste", db.over.time[, c("Kingdom", "Phylum", "Class", "Family", "genus.list", "sp.list")])
db.over.time <- db.over.time[-which(db.over.time$Impact_year < 1970), ]
species.summary <- data.frame()
for(sp in unique(db.over.time$unique.sp.id))
{
# Db for current species
cur.db <- db.over.time[which(db.over.time$unique.sp.id %in% sp), ]
cur.raw <- summarizeCosts(cur.db, minimum.year = 1970,
maximum.year = 2017)
# Management db
cur.db.management <- cur.db[which(cur.db$Type_of_cost_merged == "Management_costs"), ]
if(all(cur.db.management$Impact_year < 1970) | all(cur.db.management$Impact_year > 2017))
{
cur.db.management <- cur.db.management[-0, ]
}
if(nrow(cur.db.management))
{
cur.management <- summarizeCosts(cur.db.management, minimum.year = 1970,
maximum.year = 2017)
} else
{
cur.management <- list(average.total.cost = data.frame(annual_cost = 0,
total_cost = 0,
number_estimates = 0,
number_year_values = 0))
}
# Damage db
cur.db.damage <- cur.db[which(cur.db$Type_of_cost_merged == "Damage_costs"), ]
if(all(cur.db.damage$Impact_year < 1970) | all(cur.db.damage$Impact_year > 2017))
{
cur.db.management <- cur.db.management[-0, ]
}
if(nrow(cur.db.damage))
{
cur.damage <- summarizeCosts(cur.db.damage, minimum.year = 1970,
maximum.year = 2017)
} else
{
cur.damage <- list(average.total.cost = data.frame(annual_cost = 0,
total_cost = 0,
number_estimates = 0,
number_year_values = 0))
}
# Mixed db
cur.db.mixed <- cur.db[which(cur.db$Type_of_cost_merged == "Mixed_costs"), ]
if(all(cur.db.mixed$Impact_year < 1970) | all(cur.db.mixed$Impact_year > 2017))
{
cur.db.management <- cur.db.management[-0, ]
}
if(nrow(cur.db.mixed))
{
cur.mixed <- summarizeCosts(cur.db.mixed, minimum.year = 1970,
maximum.year = 2017)
} else
{
cur.mixed <- list(average.total.cost = data.frame(annual_cost = 0,
total_cost = 0,
number_estimates = 0,
number_year_values = 0))
}
species.summary <- rbind.data.frame(species.summary,
data.frame(
Kingdom = cur.db$Kingdom[1],
Phylum = cur.db$Phylum[1],
Class = cur.db$Class[1],
Family = cur.db$Family[1],
Genus = cur.db$Genus[1],
Species = cur.db$sp.list[1],
Cost.type = c("Total", "Management", "Damage", "Mixed"),
Average.annual.cost = c(cur.raw$average.total.cost$annual_cost,
cur.management$average.total.cost$annual_cost,
cur.damage$average.total.cost$annual_cost,
cur.mixed$average.total.cost$annual_cost),
Cumulated.cost = c(cur.raw$average.total.cost$total_cost,
cur.management$average.total.cost$total_cost,
cur.damage$average.total.cost$total_cost,
cur.mixed$average.total.cost$total_cost),
Number.estimates = c(cur.raw$average.total.cost$number_estimates,
cur.management$average.total.cost$number_estimates,
cur.damage$average.total.cost$number_estimates,
cur.mixed$average.total.cost$number_estimates),
Number.year.values = c(cur.raw$average.total.cost$number_year_values,
cur.management$average.total.cost$number_year_values,
cur.damage$average.total.cost$number_year_values,
cur.mixed$average.total.cost$number_year_values),
Unique.sp.id = cur.db$unique.sp.id[1]
))
}
# Adjust the file path to your needs
write.csv(species.summary, "d:/r/Projects/invacost_nature/data/species_summary.csv")
# Table containing only total costs (for identification of top species)
reduced.sp.summary <- species.summary[which(species.summary$Cost.type == "Total"), ]
reduced.sp.summary <- reduced.sp.summary[order(reduced.sp.summary$Cumulated.cost, decreasing = TRUE), ]
# Ordering full table according to cumulated cost
species.summary$Unique.sp.id <- factor(species.summary$Unique.sp.id,
levels = reduced.sp.summary$Unique.sp.id)
species.summary <- dplyr::arrange(species.summary,
Unique.sp.id)
# Filtering top 10 species
top.10 <- reduced.sp.summary$Species[reduced.sp.summary$Species != "Diverse/Unspecified"][1:10]
# csv file available here: www.borisleroy.com/permanent/images.csv
images <- read.csv("d:/r/Projects/invacost_nature/data/images.csv", sep = ";")
images <- images[images$Taxon %in% top.10, ]
images <- images[match(top.10, images$Taxon), ]
# Preparing silhouettes
cat(paste0("annotation_custom(grid::rasterGrob(readPNG(getURLContent('",
images$URL, "')), \n interpolate = TRUE),\n xmin = ",
1:10 - .5,
",\n xmax = ",
1:10 + .5,
",\n ymin = ",
round(reduced.sp.summary[reduced.sp.summary$Species %in% top.10, "Cumulated.cost"]),
" + y1,\n ymax = ",
round(reduced.sp.summary[reduced.sp.summary$Species %in% top.10, "Cumulated.cost"]),
" + y2)", collapse = " +\n "))
# Preparing plot
y1 <- 5.000
y2 <- 20.000
ggsummary <- species.summary[which(species.summary$Species %in% top.10 &
species.summary$Cost.type != "Total"), ]
ggsummary$Species <- factor(as.character(ggsummary$Species),
levels = as.character(top.10))
# Ordering species names
reduced.top10 <- reduced.sp.summary[reduced.sp.summary$Species %in% top.10, ]
# Tidying species names
reduced.top10$Species <- gsub(" ", "\\\n", as.character(reduced.top10$Species))
levels(ggsummary$Species) <- gsub(" ", "\\\n", levels(ggsummary$Species))
# Reducing a long name
# reduced.top10$Species[grep("Bursa", reduced.top10$Species)] <- "Bursaphe-\nlenchus\nxylophilus"
# levels(ggsummary$Species)[grep("Bursa", levels(ggsummary$Species))] <- "Bursaphe-\nlenchus\nxylophilus"
# Adding Africanized honey bee subspecies name
reduced.top10$Species[reduced.top10$Species == "Apis\nmellifera"] <- "Apis\nmellifera\nscutellata"
levels(ggsummary$Species)[levels(ggsummary$Species) == "Apis\nmellifera"] <- "Apis\nmellifera\nscutellata"
ggsummary$SpLabels <- ggsummary$Species
levels(ggsummary$SpLabels) <- c("*Aedes*<br>spp.",
"*Rattus*<br>spp./*Mus*<br>spp.",
"*Felis<br>catus*",
"*Coptotermes<br>formosanus*",
"*Solenopsis<br>invicta*",
"*Boiga<br>irregularis*",
"*Anthonomus<br>grandis*",
"*Lymantria<br>dispar*",
"*Apis<br>mellifera<br>scutellata*",
"*Cochliomyia<br>hominivorax*")
# Moving to billions
ggsummary$Cumulated.cost <- ggsummary$Cumulated.cost / 1000
# reduced.top10$Cumulated.cost <- reduced.top10$Cumulated.cost / 1000
reduced.top10$ynew <- 5
ptop10 <- ggplot(ggsummary) +
geom_col(aes(x = SpLabels, y = Cumulated.cost, fill = Cost.type)) +
theme_minimal() +
ylab("Cumulative cost for 1970-2017 in 2017 US$ billions") +
xlab("Taxon") +
scale_y_continuous(#labels = scales::comma,
limits = c(0, 175),
breaks = c(0, 10, 20, 50, 75, 100, 150),
minor_breaks = c(5, 15),
) +
scale_fill_manual(values = c(Management = rgb(red = 0, 158, 115, maxColorValue = 255), # Yellow -> green
Damage = rgb(red = 213, 94, 0, maxColorValue = 255), # Blue -> brown
Mixed = rgb(red = 86, 180, 233, maxColorValue = 255)), # Pink -> blue
name = "Cost type") +
theme(panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.text.y = element_text(size = 13),
axis.text.x = element_markdown(size = 13),
axis.title = element_text(size = 16),
panel.grid.minor.x = element_blank(),
panel.grid.major.x = element_blank(),
legend.text = element_text(size = 14),
legend.title = element_text(size = 16),
legend.position = c(.75, .65)) +
annotation_custom(grid::rasterGrob(readPNG(getURLContent('http://phylopic.org/assets/images/submissions/f538aa99-5c08-4f96-97d9-2e094ef5d84f.thumb.png')),
interpolate = TRUE),
xmin = 0.5,
xmax = 1.5,
ymin = 148.669 + y1,
ymax = 148.669 + y2) +
annotation_custom(grid::rasterGrob(readPNG(getURLContent('http://phylopic.org/assets/images/submissions/53f8f372-adc0-4247-bc66-0075e241fc95.thumb.png')),
interpolate = TRUE),
xmin = 1.5,
xmax = 2.5,
ymin = 67.106 + y1,
ymax = 67.106 + y2) +
annotation_custom(grid::rasterGrob(readPNG(getURLContent('http://phylopic.org/assets/images/submissions/23cd6aa4-9587-4a2e-8e26-de42885004c9.thumb.png')),
interpolate = TRUE),
xmin = 2.5,
xmax = 3.5,
ymin = 51.556 + y1,
ymax = 51.556 + y2) +
annotation_custom(grid::rasterGrob(readPNG(getURLContent('http://phylopic.org/assets/images/submissions/188d20e7-a6da-4d7c-99e1-52935c2b0c6d.thumb.png')),
interpolate = TRUE),
xmin = 3.5,
xmax = 4.5,
ymin = 19.003 + y1,
ymax = 19.003 + y2) +
annotation_custom(grid::rasterGrob(readPNG(getURLContent('http://phylopic.org/assets/images/submissions/028577c7-8fb7-4c0e-83ee-5f8035c6e3df.thumb.png')),
interpolate = TRUE),
xmin = 4.5,
xmax = 5.5,
ymin = 17.304 + y1,
ymax = 17.304 + y2) +
annotation_custom(grid::rasterGrob(readPNG(getURLContent('http://phylopic.org/assets/images/submissions/6045a4fc-d4f9-485b-9b4e-dd0ef32aeed5.thumb.png')),
interpolate = TRUE),
xmin = 5.5,
xmax = 6.5,
ymin = 15.014 + y1,
ymax = 15.014 + y2) +
annotation_custom(grid::rasterGrob(readPNG(getURLContent('http://phylopic.org/assets/images/submissions/d01de196-29cf-42d0-b52f-8ae433de8c30.thumb.png')),
interpolate = TRUE),
xmin = 6.5,
xmax = 7.5,
ymin = 7.432 + y1,
ymax = 7.432 + y2) +
annotation_custom(grid::rasterGrob(readPNG(getURLContent('http://phylopic.org/assets/images/submissions/c224abfd-ee39-4923-98e5-c2606dcc56cb.thumb.png')),
interpolate = TRUE),
xmin = 7.5,
xmax = 8.5,
ymin = 6.340 + y1,
ymax = 6.340 + y2) +
annotation_custom(grid::rasterGrob(readPNG(getURLContent('http://phylopic.org/assets/images/submissions/b199a5f5-20c4-4cc9-9c54-1b51578c2487.thumb.png')),
interpolate = TRUE),
xmin = 8.5,
xmax = 9.5,
ymin = 5.676 + y1,
ymax = 5.676 + y2) +
annotation_custom(grid::rasterGrob(readPNG(getURLContent('http://phylopic.org/assets/images/submissions/7252c46a-6bf1-42dd-ac5a-1018f404dfc8.thumb.png')),
interpolate = TRUE),
xmin = 9.5,
xmax = 10.5,
ymin = 4.057 + y1,
ymax = 4.057 + y2)
# +
# annotation_custom(
# grob = textGrob(label = df$n[i], hjust = 0, gp = gpar(cex = 1.5)),
# ymin = df$y[i], # Vertical position of the textGrob
# ymax = df$y[i],
# xmin = 14.3, # Note: The grobs are positioned outside the plot area
# xmax = 14.3)
# }
library(grid)
for (i in 1:nrow(reduced.top10)) {
ptop10 <- ptop10 + annotation_custom(
grob = textGrob(label = reduced.top10$Number.estimates[i], hjust = .5, gp = gpar(cex = 1)),
ymin = -3, # Vertical position of the textGrob
ymax = -3,
xmin = i, # Note: The grobs are positioned outside the plot area
xmax = i)
}
ptop10
cairo_pdf("figure5.pdf", width = 12.5, height = 8)
ptop10
dev.off()General Appendix:
Number of unspecified periods in the database
nrow(invacost[which(is.na(invacost$Probable_starting_year) |
is.na(invacost$Probable_ending_year)), ])## [1] 700nrow(invacost)## [1] 1319Number of estimates
nb.estimates <- rbind(Global = c(global.raw$parameters$number.of.estimates,
global.raw$average.cost.per.period$number_estimates),
Plants = c(obs.impact$Plants$Plants$parameters$number.of.estimates,
obs.impact$Plants$Plants$average.cost.per.period$number_estimates),
Invertebrates = c(obs.impact$Invertebrates$Invertebrates$parameters$number.of.estimates,
obs.impact$Invertebrates$Invertebrates$average.cost.per.period$number_estimates),
Vertebrates = c(obs.impact$Vertebrates$Vertebrates$parameters$number.of.estimates,
obs.impact$Vertebrates$Vertebrates$average.cost.per.period$number_estimates),
Africa = c(obs.regions$Africa$parameters$number.of.estimates,
obs.regions$Africa$average.cost.per.period$number_estimates),
Asia = c(obs.regions$Asia$parameters$number.of.estimates,
obs.regions$Asia$average.cost.per.period$number_estimates),
Europe = c(obs.regions$Europe$parameters$number.of.estimates,
obs.regions$Europe$average.cost.per.period$number_estimates),
`Central America` = c(obs.regions$`Central America`$parameters$number.of.estimates,
obs.regions$`Central America`$average.cost.per.period$number_estimates),
`North America` = c(obs.regions$`North America`$parameters$number.of.estimates,
obs.regions$`North America`$average.cost.per.period$number_estimates),
`Oceania-Pacific islands` = c(obs.regions$`Oceania-Pacific islands`$parameters$number.of.estimates,
obs.regions$`Oceania-Pacific islands`$average.cost.per.period$number_estimates),
`South America` = c(obs.regions$`South America`$parameters$number.of.estimates,
obs.regions$`South America`$average.cost.per.period$number_estimates),
Damage = c(obs.costtype$Damage$parameters$number.of.estimates,
obs.costtype$Damage$average.cost.per.period$number_estimates),
Management = c(obs.costtype$Management$parameters$number.of.estimates,
obs.costtype$Management$average.cost.per.period$number_estimates)
)
colnames(nb.estimates) <- c("1970 - 2017",
paste0(global.raw$average.cost.per.period$initial_year,
" - ",
global.raw$average.cost.per.period$final_year))
data.frame(nb.estimates)| X1970…2017 | X1970…1979 | X1980…1989 | X1990…1999 | X2000…2009 | X2010…2017 | |
|---|---|---|---|---|---|---|
| Global | 1312 | 25 | 88 | 276 | 783 | 336 |
| Plants | 221 | 3 | 5 | 63 | 168 | 18 |
| Invertebrates | 469 | 16 | 43 | 105 | 241 | 164 |
| Vertebrates | 526 | 6 | 39 | 84 | 314 | 131 |
| Africa | 130 | 2 | 4 | 36 | 76 | 46 |
| Asia | 105 | 1 | 8 | 11 | 68 | 44 |
| Europe | 151 | 3 | 2 | 12 | 101 | 52 |
| Central America | 47 | 3 | 3 | 11 | 26 | 10 |
| North America | 310 | 13 | 30 | 91 | 163 | 73 |
| South America | 68 | 0 | 1 | 11 | 32 | 29 |
| Damage | 402 | 8 | 23 | 66 | 209 | 148 |
| Management | 873 | 15 | 61 | 201 | 552 | 177 |
# adjust file path
write.table(nb.estimates,
"d:/r/projects/invacost_nature/data/nb_estimates.csv", sep = ";")Observed values
Average.values <- rbind(Global = c(global.raw$average.total.cost$annual_cost,
global.raw$average.cost.per.period$annual_cost),
Plants = c(obs.impact$Plants$Plants$average.total.cost$annual_cost,
obs.impact$Plants$Plants$average.cost.per.period$annual_cost),
Invertebrates = c(obs.impact$Invertebrates$Invertebrates$average.total.cost$annual_cost,
obs.impact$Invertebrates$Invertebrates$average.cost.per.period$annual_cost),
Vertebrates = c(obs.impact$Vertebrates$Vertebrates$average.total.cost$annual_cost,
obs.impact$Vertebrates$Vertebrates$average.cost.per.period$annual_cost),
Africa = c(obs.regions$Africa$average.total.cost$annual_cost,
obs.regions$Africa$average.cost.per.period$annual_cost),
Asia = c(obs.regions$Asia$average.total.cost$annual_cost,
obs.regions$Asia$average.cost.per.period$annual_cost),
Europe = c(obs.regions$Europe$average.total.cost$annual_cost,
obs.regions$Europe$average.cost.per.period$annual_cost),
`Central America` = c(obs.regions$`Central America`$average.total.cost$annual_cost,
obs.regions$`Central America`$average.cost.per.period$annual_cost),
`North America` = c(obs.regions$`North America`$average.total.cost$annual_cost,
obs.regions$`North America`$average.cost.per.period$annual_cost),
`Oceania-Pacific islands` = c(obs.regions$`Oceania-Pacific islands`$average.total.cost$annual_cost,
obs.regions$`Oceania-Pacific islands`$average.cost.per.period$annual_cost),
`South America` = c(obs.regions$`South America`$average.total.cost$annual_cost,
obs.regions$`South America`$average.cost.per.period$annual_cost),
Damage = c(obs.costtype$Damage$average.total.cost$annual_cost,
obs.costtype$Damage$average.cost.per.period$annual_cost),
Management = c(obs.costtype$Management$average.total.cost$annual_cost,
obs.costtype$Management$average.cost.per.period$annual_cost)
)
Cumulated.values <- rbind(Global = c(global.raw$average.total.cost$total_cost,
global.raw$average.cost.per.period$total_cost),
Plants = c(obs.impact$Plants$Plants$average.total.cost$total_cost,
obs.impact$Plants$Plants$average.cost.per.period$total_cost),
Invertebrates = c(obs.impact$Invertebrates$Invertebrates$average.total.cost$total_cost,
obs.impact$Invertebrates$Invertebrates$average.cost.per.period$total_cost),
Vertebrates = c(obs.impact$Vertebrates$Vertebrates$average.total.cost$total_cost,
obs.impact$Vertebrates$Vertebrates$average.cost.per.period$total_cost),
Africa = c(obs.regions$Africa$average.total.cost$total_cost,
obs.regions$Africa$average.cost.per.period$total_cost),
Asia = c(obs.regions$Asia$average.total.cost$total_cost,
obs.regions$Asia$average.cost.per.period$total_cost),
Europe = c(obs.regions$Europe$average.total.cost$total_cost,
obs.regions$Europe$average.cost.per.period$total_cost),
`Central America` = c(obs.regions$`Central America`$average.total.cost$total_cost,
obs.regions$`Central America`$average.cost.per.period$total_cost),
`North America` = c(obs.regions$`North America`$average.total.cost$total_cost,
obs.regions$`North America`$average.cost.per.period$total_cost),
`Oceania-Pacific islands` = c(obs.regions$`Oceania-Pacific islands`$average.total.cost$total_cost,
obs.regions$`Oceania-Pacific islands`$average.cost.per.period$total_cost),
`South America` = c(obs.regions$`South America`$average.total.cost$total_cost,
obs.regions$`South America`$average.cost.per.period$total_cost),
Damage = c(obs.costtype$Damage$average.total.cost$total_cost,
obs.costtype$Damage$average.cost.per.period$total_cost),
Management = c(obs.costtype$Management$average.total.cost$total_cost,
obs.costtype$Management$average.cost.per.period$total_cost)
)
colnames(Cumulated.values) <- colnames(Average.values) <- c("1970 - 2017",
paste0(global.raw$average.cost.per.period$initial_year,
" - ",
global.raw$average.cost.per.period$final_year))
data.frame(Cumulated.values)| X1970…2017 | X1970…1979 | X1980…1989 | X1990…1999 | X2000…2009 | X2010…2017 | |
|---|---|---|---|---|---|---|
| Global | 1288087.867 | 6786.650730 | 11887.065382 | 202749.6227 | 833082.633 | 233581.8950 |
| Plants | 8902.476 | 5.750304 | 119.138716 | 3305.5941 | 2715.327 | 2756.6657 |
| Invertebrates | 415696.898 | 4279.977292 | 11641.175266 | 24631.3002 | 227826.743 | 147317.7016 |
| Vertebrates | 165949.140 | 2500.923134 | 121.068662 | 529.5412 | 143040.446 | 19757.1609 |
| Africa | 10190.337 | 6.368962 | 26.444103 | 2983.5722 | 3696.007 | 3477.9448 |
| Asia | 269192.518 | 17.554213 | 1621.171151 | 573.5950 | 208620.189 | 58360.0088 |
| Europe | 28798.739 | 192.488046 | 6.240115 | 240.6275 | 10640.305 | 17719.0790 |
| Central America | 5758.668 | 43.628069 | 277.671895 | 127.0589 | 4668.727 | 641.5825 |
| North America | 528738.476 | 6395.066605 | 9101.970316 | 30726.1629 | 480928.954 | 1586.3222 |
| South America | 84940.845 | NA | 3.266089 | 117.3816 | 79914.538 | 4905.6590 |
| Damage | 892179.562 | 846.726336 | 5123.951925 | 186202.3532 | 562118.209 | 137888.3213 |
| Management | 66013.198 | 3649.054056 | 3427.519751 | 9487.2330 | 30189.543 | 19259.8486 |
# adjust file path
write.table(Average.values,
"d:/r/projects/invacost_nature/data/average_costs.csv", sep = ";")
# adjust file path
write.table(Cumulated.values,
"d:/r/projects/invacost_nature/data/cumulated_costs.csv", sep = ";")Predicted values
Predicted.values <- format(data.frame(Global = global.trend$final.year.cost,
Plants = pred.impact$Plants$Plants$final.year.cost,
Invertebrates = pred.impact$Invertebrates$Invertebrates$final.year.cost,
Vertebrates = pred.impact$Vertebrates$Vertebrates$final.year.cost,
# Africa = pred.regions$Africa$final.year.cost,
# Asia = pred.regions$Asia$final.year.cost,
# Europe = pred.regions$Europe$final.year.cost,
# `Central America` = pred.regions$`Central America`$final.year.cost,
# `North America` = pred.regions$`North America`$final.year.cost,
# `Oceania-Pacific islands` = pred.regions$`Oceania-Pacific islands`$final.year.cost,
# `South America` = pred.regions$`Oceania-Pacific islands`$final.year.cost,
Damage.25 = pred.costtype.25$Damage$final.year.cost,
Management.25 = pred.costtype.25$Management$final.year.cost
), scientific = FALSE)
data.frame(Predicted.values)| Global | Plants | Invertebrates | Vertebrates | Damage.25 | Management.25 | |
|---|---|---|---|---|---|---|
| ols.linear | 91979.90 | 490.9546987 | 33359.412 | 2382.2035 | 84972.961 | 2002.0848 |
| ols.quadratic | 102620.61 | 11.1742742 | 45830.139 | 17140.9960 | 16753.597 | 1181.3576 |
| robust.linear | 46822.77 | 3842.3084162 | 23757.670 | 1355.2810 | 68716.060 | 2029.3408 |
| robust.quadratic | 86535.24 | 733.9635103 | 59444.789 | 13390.8015 | 25209.512 | 570.6399 |
| mars | 162718.71 | 0.1758477 | 11513.989 | 5794.5769 | 16691.964 | 356.4238 |
| gam | 78794.54 | 23.0378000 | 32008.330 | 2795.7431 | 11959.276 | 555.4946 |
| quantile.0.1 | 16223.81 | 0.7585002 | 6490.447 | 401.6485 | 4957.238 | 523.3445 |
| quantile.0.5 | 57492.87 | 4225.7159107 | 17970.458 | 1009.6185 | 84498.995 | 2156.3986 |
| quantile.0.9 | 2219915.57 | 17235.9374532 | 132228.772 | 595338.4192 | 6713500.332 | 19828.9286 |
# adjust file path
write.table(Predicted.values,
"d:/r/projects/invacost_nature/data/predicted_costs.csv", sep = ";")Appendix to figure A
Appendix global.1: time lag in publications
Time lag
# Analysis of time lag
lag_plot <- ggplot(db.over.time,
aes(y = Publication_lag)) +
geom_boxplot(outlier.alpha = .2) +
ylab("Publication lag (in years)") +
theme_minimal() +
theme(panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank(),
panel.grid.minor.x = element_blank(),
axis.text.y = element_blank(),
axis.text = element_text(size = 16),
axis.title = element_text(size = 18)) +
scale_y_continuous(breaks = c(-25, 0, 5, 10, 15, 25, 50, 75, 100)) +
xlab("") +
coord_flip()
cairo_pdf("extend_data1.pdf", height = 4, width = 12)
lag_plot
dev.off()## png
## 2lag_plot
Relationship between number of estimates and cost
# Relationship between number of estimates & cost value
ggplot(global.raw$cost.per.year,
aes(x = number_estimates,
y = cost,
col = year)) +
geom_point() +
geom_smooth() +
ylab("Cost per year in US$ million") +
xlab("Number of estimates") +
theme_minimal() +
scale_y_log10(breaks = 10^(-15:15),
labels = scales::comma) +
annotation_logticks(sides = "l") +
scale_color_gradient(low = "yellow", high = "red") +
guides(col = guide_legend(override.aes = list(size = 2)),
shape = guide_legend(override.aes = list(size = 2))) +
theme(panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14),
panel.grid.minor = element_blank(),
legend.text = element_text(size = 14),
legend.title = element_text(size = 18),
strip.text = element_text(size = 14))
# Number of estimates per year & cost value
ggplot(global.raw$cost.per.year,
aes(x = year,
y = number_estimates,
size = cost)) +
geom_point() +
geom_smooth() +
ylab("Number of estimates") +
xlab("Year") +
theme_minimal() +
guides(col = guide_legend(override.aes = list(size = 2)),
shape = guide_legend(override.aes = list(size = 2))) +
theme(panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14),
panel.grid.minor = element_blank(),
legend.text = element_text(size = 14),
legend.title = element_text(size = 18),
strip.text = element_text(size = 14))
Appendix global.2: Model fitting and evaluation
[This section will be completed in the Appendix in the main paper]
The 0.1 and 0.9 quantile regressions suggested increasing and extreme differences (0.1 quantile estimate in 2017: US$ billion 16.22; 0.9 quantile estimate in 2017: US$ billion 2219.92) between the years with lowest impacts and years with highest impacts. However, the incompleteness of data for recent years makes these estimations highly speculative for periods later than 2004.
Years below / above average trend
lin.years <- global.trend$estimated.annual.costs[which(global.trend$estimated.annual.costs$Year < global.trend$parameters$maximum.year &
global.trend$estimated.annual.costs$Details == "Linear"), ]
lin.years$observed <- global.trend$cost.data$Annual.cost[match(lin.years$Year,
global.trend$cost.data$Year)]
lin.years$diff <- lin.years$observed - lin.years$fit
lin.years$diff.upr <- lin.years$observed - lin.years$upr
lin.years$diff.lwr <- lin.years$observed - lin.years$lwr
length(which(lin.years$diff.upr > 0))## [1] 21length(which(lin.years$diff.lwr < 0))## [1] 31Model summaries
global.trend$model.summary##
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Summary of model fits ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
##
##
## ______________________________ Ordinary Least Square regression models _______________________________
##
##
## >>>>>>>> Linear regression
##
## R squared: 0.6742191 - Adjusted R squared: 0.6742191 Estimate Std. Error t value Pr(>|t|)
## (Intercept) -105.4890449 13.999669712 -7.535110 2.179750e-09
## Year 0.0547609 0.007061762 7.754566 1.058425e-09
## attr(,"class")
## [1] "coeftest"
## attr(,"method")
## [1] "t test of coefficients"
## attr(,"df")
## [1] 43
## attr(,"nobs")
## [1] 45
## attr(,"logLik")
## 'log Lik.' -32.15056 (df=3)
## ------------------------------------------------------------------------------------------------------------
##
##
##
## >>>>>>>> Quadratic regression
##
## R squared: 0.6745479 - Adjusted R squared: 0.6745479 Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.078936e+02 2.216600e+03 0.1389035 0.8901903
## Year -3.602996e-01 2.224609e+00 -0.1619608 0.8721134
## I(Year^2) 1.041819e-04 5.581298e-04 0.1866624 0.8528241
## attr(,"class")
## [1] "coeftest"
## attr(,"method")
## [1] "t test of coefficients"
## attr(,"df")
## [1] 42
## attr(,"nobs")
## [1] 45
## attr(,"logLik")
## 'log Lik.' -32.12784 (df=4)
## ------------------------------------------------------------------------------------------------------------
##
## ______________________________ Robust regression models _______________________________
##
##
## >>>>>>>> Linear regression
##
##
## Call:
## robustbase::lmrob(formula = transf.cost ~ Year, data = yearly.cost.calibration,
## weights = incomplete.year.weights)
## \--> method = "MM"
## Residuals:
## Min 1Q Median 3Q Max
## -0.41349 -0.17224 -0.01014 0.34127 1.53296
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -93.661247 8.990966 -10.42 2.46e-13 ***
## Year 0.048751 0.004537 10.74 9.31e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Robust residual standard error: 0.3216
## Multiple R-squared: 0.7892, Adjusted R-squared: 0.7843
## Convergence in 12 IRWLS iterations
##
## Robustness weights:
## observation 35 is an outlier with |weight| = 0 ( < 0.0022);
## 5 weights are ~= 1. The remaining 39 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.009949 0.825200 0.924800 0.836000 0.986500 0.998700
## Algorithmic parameters:
## tuning.chi bb tuning.psi refine.tol
## 1.548e+00 5.000e-01 4.685e+00 1.000e-07
## rel.tol scale.tol solve.tol eps.outlier
## 1.000e-07 1.000e-10 1.000e-07 2.222e-03
## eps.x warn.limit.reject warn.limit.meanrw
## 3.663e-09 5.000e-01 5.000e-01
## nResample max.it best.r.s k.fast.s k.max
## 500 50 2 1 200
## maxit.scale trace.lev mts compute.rd fast.s.large.n
## 200 0 1000 0 2000
## psi subsampling cov
## "bisquare" "nonsingular" ".vcov.avar1"
## compute.outlier.stats
## "SM"
## seed : int(0)
## ------------------------------------------------------------------------------------------------------------
##
##
##
## >>>>>>>> Quadratic regression
##
##
## Call:
## robustbase::lmrob(formula = transf.cost ~ Year + I(Year^2), data = yearly.cost.calibration,
## weights = incomplete.year.weights, cov = ".vcov.w")
## \--> method = "MM"
## Residuals:
## Min 1Q Median 3Q Max
## -0.43887 -0.17229 -0.02735 0.30656 1.57689
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.511e+03 1.255e+03 2.001 0.0519 .
## Year -2.566e+00 1.260e+00 -2.037 0.0480 *
## I(Year^2) 6.562e-04 3.161e-04 2.076 0.0441 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Robust residual standard error: 0.2752
## Multiple R-squared: 0.8342, Adjusted R-squared: 0.8263
## Convergence in 10 IRWLS iterations
##
## Robustness weights:
## 4 observations c(30,31,32,35) are outliers with |weight| = 0 ( < 0.0022);
## 3 weights are ~= 1. The remaining 38 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.3812 0.8298 0.9624 0.8867 0.9884 0.9989
## Algorithmic parameters:
## tuning.chi bb tuning.psi refine.tol
## 1.548e+00 5.000e-01 4.685e+00 1.000e-07
## rel.tol scale.tol solve.tol eps.outlier
## 1.000e-07 1.000e-10 1.000e-07 2.222e-03
## eps.x warn.limit.reject warn.limit.meanrw
## 7.378e-06 5.000e-01 5.000e-01
## nResample max.it best.r.s k.fast.s k.max
## 500 50 2 1 200
## maxit.scale trace.lev mts compute.rd fast.s.large.n
## 200 0 1000 0 2000
## psi subsampling cov
## "bisquare" "nonsingular" ".vcov.w"
## compute.outlier.stats
## "SM"
## seed : int(0)
## ------------------------------------------------------------------------------------------------------------
##
## ______________________________ Generalized Additive Models _______________________________
##
##
##
## Family: gaulss
## Link function: identity logb
##
## Formula:
## transf.cost ~ s(Year, k = gam.k)
## <environment: 0x0000000015f16120>
## ~s(Year, k = gam.k)
## <environment: 0x0000000015f16120>
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.48116 0.05636 61.762 <2e-16 ***
## (Intercept).1 -1.09501 0.11617 -9.426 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(Year) 2.529 3.005 235.03 < 2e-16 ***
## s.1(Year) 4.427 5.248 23.51 0.000416 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Deviance explained = 88.1%
## -REML = 32.439 Scale est. = 1 n = 45
## ------------------------------------------------------------------------------------------------------------
##
## ______________________________ Multiple Adaptive Regression Splines _______________________________
##
##
## Call: earth(formula=transf.cost~Year, data=yearly.cost.calibration,
## weights=incomplete.year.weights, pmethod="backward",
## nprune=mars.nprune, nfold=5, ncross=30, varmod.method="lm")
##
## coefficients
## (Intercept) 2.8607147
## h(Year-1989) 0.1501641
## h(Year-2004) -0.5690581
## h(Year-2007) 0.5543884
##
## Selected 4 of 6 terms, and 1 of 1 predictors (nprune=10)
## Termination condition: RSq changed by less than 0.001 at 6 terms
## Importance: Year
## Number of terms at each degree of interaction: 1 3 (additive model)
## GCV 0.2061979 RSS 6.616662 GRSq 0.7372182 RSq 0.8039996 CVRSq 0.59976
##
## Note: the cross-validation sd's below are standard deviations across folds
##
## Cross validation: nterms 3.51 sd 0.82 nvars 1.00 sd 0.00
##
## CVRSq sd MaxErr sd
## 0.6 0.201 1.62 0.857
##
## varmod: method "lm" min.sd 0.0482 iter.rsq 0.161
##
## stddev of predictions:
## coefficients iter.stderr iter.stderr%
## (Intercept) -0.1348638 0.196846 146
## transf.cost 0.1715754 0.0596978 35
##
## mean smallest largest ratio
## 95% prediction interval 1.888986 1.395355 2.910277 2.08569
##
## 68% 80% 90% 95%
## response values in prediction interval 82 91 98 98
## ------------------------------------------------------------------------------------------------------------
##
## ______________________________ Quantile regressions _______________________________
##
##
## >>>>>>>> 0.1 quantile
##
## $call
## quantreg::rq(formula = transf.cost ~ Year, tau = 0.1, data = yearly.cost.calibration,
## weights = incomplete.year.weights)
##
## $terms
## transf.cost ~ Year
## attr(,"variables")
## list(transf.cost, Year)
## attr(,"factors")
## Year
## transf.cost 0
## Year 1
## attr(,"term.labels")
## [1] "Year"
## attr(,"order")
## [1] 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: 0x0000000015f16120>
## attr(,"predvars")
## list(transf.cost, Year)
## attr(,"dataClasses")
## transf.cost Year
## "numeric" "numeric"
##
## $coefficients
## coefficients lower bd upper bd
## (Intercept) -83.03585596 -93.40128046 -52.12008831
## Year 0.04325533 0.02758737 0.04842056
##
## $residuals
## 1 2 3 4 5
## 8.153171e-01 3.435940e-01 3.015330e-01 2.570833e-01 2.138280e-01
## 6 7 8 9 10
## 1.979126e-01 1.553110e-01 1.865275e-01 7.371635e-02 8.643924e-01
## 11 12 13 14 15
## 2.520190e-01 1.009602e+00 1.775473e-01 1.396046e-01 1.528939e-01
## 16 17 18 19 20
## 3.944301e-02 6.041518e-03 5.042964e-04 -3.361350e-03 -6.060688e-04
## 21 22 23 24 25
## 4.729454e-01 -7.070252e-02 3.171465e-01 2.081584e-01 8.293227e-01
## 26 27 28 29 30
## 5.168462e-01 -6.860498e-02 1.199041e-14 3.411796e-01 1.790943e+00
## 31 32 33 34 35
## 1.756788e+00 1.635481e+00 8.879341e-01 4.912530e-01 1.921811e+00
## 36 37 38 39 40
## 1.125952e+00 2.776900e-01 4.657592e-01 4.440892e-16 5.751826e-02
## 41 42 43 44 45
## 6.809377e-02 9.352929e-01 6.540670e-01 3.826862e-01 7.850811e-01
##
## $rdf
## [1] 43
##
## $tau
## [1] 0.1
##
## attr(,"class")
## [1] "summary.rq"
## ------------------------------------------------------------------------------------------------------------
##
## >>>>>>>> 0.5 quantile
##
## $call
## quantreg::rq(formula = transf.cost ~ Year, tau = 0.5, data = yearly.cost.calibration,
## weights = incomplete.year.weights)
##
## $terms
## transf.cost ~ Year
## attr(,"variables")
## list(transf.cost, Year)
## attr(,"factors")
## Year
## transf.cost 0
## Year 1
## attr(,"term.labels")
## [1] "Year"
## attr(,"order")
## [1] 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: 0x0000000015f16120>
## attr(,"predvars")
## list(transf.cost, Year)
## attr(,"dataClasses")
## transf.cost Year
## "numeric" "numeric"
##
## $coefficients
## coefficients lower bd upper bd
## (Intercept) -99.36909811 -110.99860376 -89.51713784
## Year 0.05162554 0.04665242 0.05741113
##
## $residuals
## 1 2 3 4 5
## 6.592555e-01 1.791622e-01 1.287310e-01 7.591114e-02 2.428560e-02
## 6 7 8 9 10
## 4.440892e-15 -5.097181e-02 -2.812547e-02 -1.493069e-01 6.329990e-01
## 11 12 13 14 15
## 1.225534e-02 7.614686e-01 -7.895669e-02 -1.252696e-01 -1.203506e-01
## 16 17 18 19 20
## -2.421716e-01 -2.839433e-01 -2.978508e-01 -3.100866e-01 -3.157015e-01
## 21 22 23 24 25
## 1.494797e-01 -4.025384e-01 -2.305962e-02 -1.404179e-01 4.723762e-01
## 26 27 28 29 30
## 1.515295e-01 -4.422919e-01 -3.820571e-01 -4.924773e-02 1.392145e+00
## 31 32 33 34 35
## 1.349620e+00 1.219943e+00 4.640260e-01 5.897470e-02 1.481162e+00
## 36 37 38 39 40
## 6.769332e-01 -1.796990e-01 -2.664535e-15 -4.741294e-01 -4.249813e-01
## 41 42 43 44 45
## -4.227760e-01 4.360529e-01 1.464568e-01 -1.332942e-01 2.607305e-01
##
## $rdf
## [1] 43
##
## $tau
## [1] 0.5
##
## attr(,"class")
## [1] "summary.rq"
## ------------------------------------------------------------------------------------------------------------
##
## >>>>>>>> 0.9 quantile
##
## $call
## quantreg::rq(formula = transf.cost ~ Year, tau = 0.9, data = yearly.cost.calibration,
## weights = incomplete.year.weights)
##
## $terms
## transf.cost ~ Year
## attr(,"variables")
## list(transf.cost, Year)
## attr(,"factors")
## Year
## transf.cost 0
## Year 1
## attr(,"term.labels")
## [1] "Year"
## attr(,"order")
## [1] 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: 0x0000000015f16120>
## attr(,"predvars")
## list(transf.cost, Year)
## attr(,"dataClasses")
## transf.cost Year
## "numeric" "numeric"
##
## $coefficients
## coefficients lower bd upper bd
## (Intercept) -144.01951778 -173.73182018 -47.46778706
## Year 0.07454926 0.03056296 0.08959365
##
## $residuals
## 1 2 3 4 5
## 1.499479e-01 -3.530692e-01 -4.264241e-01 -5.021677e-01 -5.767169e-01
## 6 7 8 9 10
## -6.239263e-01 -6.978218e-01 -6.978992e-01 -8.420043e-01 -8.262215e-02
## 11 12 13 14 15
## -7.262895e-01 -1.820766e-14 -8.633490e-01 -9.325856e-01 -9.505903e-01
## 16 17 18 19 20
## -1.095335e+00 -1.160031e+00 -1.196862e+00 -1.232021e+00 -1.260560e+00
## 21 22 23 24 25
## -8.183023e-01 -1.393244e+00 -1.036689e+00 -1.176971e+00 -5.871007e-01
## 26 27 28 29 30
## -9.308712e-01 -1.547616e+00 -1.510305e+00 -1.200420e+00 2.180500e-01
## 31 32 33 34 35
## 1.526006e-01 4.440892e-15 -7.788407e-01 -1.206816e+00 1.924481e-01
## 36 37 38 39 40
## -6.347047e-01 -1.514261e+00 -1.357485e+00 -1.854538e+00 -1.828314e+00
## 41 42 43 44 45
## -1.849032e+00 -1.013127e+00 -1.325647e+00 -1.628322e+00 -1.257221e+00
##
## $rdf
## [1] 43
##
## $tau
## [1] 0.9
##
## attr(,"class")
## [1] "summary.rq"
## ------------------------------------------------------------------------------------------------------------Summary of model summaries
summary.models <- prettySummary(global.trend)
summary.models| 1 | 2 | 3 | 4 | 5 | 6 |
|---|---|---|---|---|---|
| Ordinary Least Square regression - Linear | |||||
| Estimate | Standard error | t value | p-value | ||
| Intercept | -105.489044881804 | 13.9996697121894 | -7.53510954547412 | 2.17975036848317e-09 | |
| Year | 0.0547609012543007 | 0.00706176190481712 | 7.75456635219407 | 1.05842544867867e-09 | |
| Adjusted R² | R² | ||||
| 0.666642835480016 | 0.674219134673652 | ||||
| _________________ | _________________ | _________________ | _________________ | _________________ | _________________ |
| Ordinary Least Square regression - Quadratic | |||||
| Estimate | Standard error | t value | p-value | ||
| Intercept | 307.89364298252 | 2216.60033161386 | 0.138903544581873 | 0.890190272481419 | |
| Year | -0.360299600623382 | 2.22460943198462 | -0.16196083476188 | 0.87211336123343 | |
| Adjusted R² | R² | ||||
| 0.65905020720394 | 0.674547925058306 | ||||
| _________________ | _________________ | _________________ | _________________ | _________________ | _________________ |
| Robust regression - Linear | |||||
| Estimate | Standard error | t value | p-value | ||
| Intercept | -93.6612474720254 | 8.99096629470269 | -10.4172615492073 | 2.45992981202492e-13 | |
| Year | 0.0487514648565585 | 0.00453702659918959 | 10.7452455458971 | 9.31494549229777e-14 | |
| Adjusted R² | R² | ||||
| 0.784258712932898 | 0.789161924002605 | ||||
| Summary of model weights | |||||
| Min | 25% | 50% | 75% | Max | |
| 0 | 0.855040806645515 | 0.947822254937868 | 0.992814333923444 | 0.999978599527095 | |
| Number of outliers | |||||
| 1 | |||||
| _________________ | _________________ | _________________ | _________________ | _________________ | _________________ |
| Robust regression - Quadratic | |||||
| Estimate | Standard error | t value | p-value | ||
| Intercept | 2510.60936235261 | 1254.59974160552 | 2.00112376807903 | 0.0518685148548383 | |
| Year | -2.56573294679132 | 1.25956772803921 | -2.03699482741229 | 0.0479852705373734 | |
| Adjusted R² | R² | ||||
| 0.826345216114038 | 0.834238615381582 | ||||
| Summary of model weights | |||||
| Min | 25% | 50% | 75% | Max | |
| 0 | 0.79787443115614 | 0.961674968633241 | 0.989012173589972 | 0.999785429572526 | |
| Number of outliers | |||||
| 4 | |||||
| _________________ | _________________ | _________________ | _________________ | _________________ | _________________ |
| Multivariate Adaptive Regression Splines | |||||
| log10(cost) | |||||
| (Intercept) | 2.86071468835477 | ||||
| h(Year-2004) | -0.569058144886949 | ||||
| h(Year-1989) | 0.150164093225849 | ||||
| h(Year-2007) | 0.554388408884719 | ||||
| Generalized R² | R² | Generalized Cross-Validation | Root Sum of Squares | ||
| 0.737218238108215 | 0.803999553630301 | 0.20619790917489 | 6.61666179663426 | ||
| Variance model | |||||
| Estimate | Standard error (last iteration) | Standard error/coefficient % | |||
| Intercept | -0.141609942882201 | 0.197711891600464 | 139.617238434262 | ||
| Intercept | 0.172483087955956 | 0.0600351142715953 | 34.8063772414171 | ||
| R² for last iteration | |||||
| 0.161046654947648 | |||||
| _________________ | _________________ | _________________ | _________________ | _________________ | _________________ |
| Generalized Additive Models | |||||
| Parametric coefficients | |||||
| Estimate | Standard error | z value | p-value | ||
| Intercept (mean) | 3.48115794132384 | 0.0563637125252753 | 61.7623961473932 | 0 | |
| Intercept (sd) | -1.09501266879915 | 0.116168267180074 | -9.42609109509875 | 4.25656733225113e-21 | |
| Smooth terms | |||||
| Estimated degree of freedom | Residual degree of freedom | Chi² | p-value | ||
| smooth (mean) | 2.52867546861836 | 3.00513116096814 | 235.032869469843 | 0 | |
| smooth (sd) | 4.42701336746575 | 5.24828108145192 | 23.5104232437925 | 0.000415848386082274 | |
| Explained deviance (%) | |||||
| 88.1453910642218 | |||||
| _________________ | _________________ | _________________ | _________________ | _________________ | _________________ |
| Quantile regression | |||||
| Coefficients quantile 0.1 | Coefficients quantile 0.5 | Coefficients quantile 0.9 | |||
| Intercept | -83.0358559569772 | -99.3690981149183 | -144.019517780825 | ||
| Year | 0.0432553340244475 | 0.051625538965215 | 0.0745492584223919 | ||
| _________________ | _________________ | _________________ | _________________ | _________________ | _________________ |
# adjust file path
write.csv(summary.models,
"d:/r/Projects/invacost_nature/data/globaltrend_modelsummary.csv")Model figures
plot(global.trend,
models = c("ols.linear", "robust.linear", "robust.quadratic", "gam",
"mars", "quantile"),
graphical.parameters =
ylab(paste0("Annual cost in US$ millions"))) +
theme_minimal() + xlab("Year") +
scale_y_log10(breaks = 10^(-15:15),
labels = scales::comma) +
annotation_logticks(sides = "l") +
guides(col = guide_legend(override.aes = list(size = 2)),
shape = guide_legend(override.aes = list(size = 2))) +
theme(panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14),
panel.grid.minor = element_blank(),
legend.text = element_text(size = 14),
legend.title = element_text(size = 18),
strip.text = element_text(size = 14))
Values
global.raw##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 1312 (number of individual year values: 3573)
## - Cost values in US$ millions:
## o Total cost over the entire period 1,288,087.87
## o Average annual cost over the entire period 26,835.16
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 6,786.65 678.67 25
## 2 1980 1989 10 11,887.07 1,188.71 88
## 3 1990 1999 10 202,749.62 20,274.96 276
## 4 2000 2009 10 833,082.63 83,308.26 783
## 5 2010 2017 8 233,581.89 29,197.74 336
## number_year_values middle.years
## 1 84 1974.5
## 2 268 1984.5
## 3 781 1994.5
## 4 1970 2004.5
## 5 470 2013.5global.trend##
## Estimation of annual cost values of invasive alien species over time
##
## - Temporal interval of data: [1970, 2017]
## - Temporal interval used for model calibration: [1970, 2014]
## - Cost transformation: log10
## - Values transformed in US$ million: Yes
## - Estimated average annual cost of invasive alien species in 2017:
##
## o Linear regression:
## . Linear: US$ million 91,979.90
## . Quadratic: US$ million 102,620.61
## o Robust regression:
## . Linear: US$ million 46,822.77
## . Quadratic: US$ million 86,535.24
## o Multiple Adapative Regression Splines: US$ million 162,718.71
## o Generalized Additive Model: US$ million 78,794.54
## o Quantile regression:
## . Quantile 0.1: US$ million 16,223.81
## . Quantile 0.5: US$ million 57,492.87
## . Quantile 0.9: US$ million 2,219,915.57
##
## You can inspect the summary of each fitted model with summary(object)Appendix global.3
global.trend$estimated.annual.costs| model | Year | Details | fit | lwr | upr | |
|---|---|---|---|---|---|---|
| 1970 | OLS regression | 1970 | Linear | 245.4317 | 136.7201 | 5.822818e+02 |
| 1971 | OLS regression | 1971 | Linear | 278.4143 | 158.5962 | 5.956299e+02 |
| 1972 | OLS regression | 1972 | Linear | 315.8293 | 183.7766 | 6.149900e+02 |
| 1973 | OLS regression | 1973 | Linear | 358.2724 | 212.7004 | 6.419215e+02 |
| 1974 | OLS regression | 1974 | Linear | 406.4192 | 245.8468 | 6.787726e+02 |
| 1975 | OLS regression | 1975 | Linear | 461.0363 | 283.7342 | 7.289600e+02 |
| 1976 | OLS regression | 1976 | Linear | 522.9931 | 326.9176 | 7.971108e+02 |
| 1977 | OLS regression | 1977 | Linear | 593.2762 | 375.9849 | 8.887843e+02 |
| 1978 | OLS regression | 1978 | Linear | 673.0042 | 431.5534 | 1.009713e+03 |
| 1979 | OLS regression | 1979 | Linear | 763.4467 | 494.2652 | 1.165108e+03 |
| 1980 | OLS regression | 1980 | Linear | 866.0433 | 564.7859 | 1.359636e+03 |
| 1981 | OLS regression | 1981 | Linear | 982.4275 | 643.8046 | 1.597923e+03 |
| 1982 | OLS regression | 1982 | Linear | 1114.4521 | 732.0392 | 1.885044e+03 |
| 1983 | OLS regression | 1983 | Linear | 1264.2190 | 830.2450 | 2.226778e+03 |
| 1984 | OLS regression | 1984 | Linear | 1434.1125 | 939.2290 | 2.629682e+03 |
| 1985 | OLS regression | 1985 | Linear | 1626.8373 | 1059.8671 | 3.101089e+03 |
| 1986 | OLS regression | 1986 | Linear | 1845.4616 | 1193.1228 | 3.649076e+03 |
| 1987 | OLS regression | 1987 | Linear | 2093.4660 | 1340.0674 | 4.282441e+03 |
| 1988 | OLS regression | 1988 | Linear | 2374.7988 | 1501.8981 | 5.010688e+03 |
| 1989 | OLS regression | 1989 | Linear | 2693.9387 | 1679.9557 | 5.844031e+03 |
| 1990 | OLS regression | 1990 | Linear | 3055.9667 | 1875.7403 | 6.793421e+03 |
| 1991 | OLS regression | 1991 | Linear | 3466.6462 | 2090.9261 | 7.870613e+03 |
| 1992 | OLS regression | 1992 | Linear | 3932.5153 | 2327.3766 | 9.088287e+03 |
| 1993 | OLS regression | 1993 | Linear | 4460.9907 | 2587.1597 | 1.046024e+04 |
| 1994 | OLS regression | 1994 | Linear | 5060.4859 | 2872.5644 | 1.200168e+04 |
| 1995 | OLS regression | 1995 | Linear | 5740.5449 | 3186.1190 | 1.372969e+04 |
| 1996 | OLS regression | 1996 | Linear | 6511.9945 | 3530.6118 | 1.566384e+04 |
| 1997 | OLS regression | 1997 | Linear | 7387.1161 | 3909.1141 | 1.782711e+04 |
| 1998 | OLS regression | 1998 | Linear | 8379.8419 | 4325.0056 | 2.024712e+04 |
| 1999 | OLS regression | 1999 | Linear | 9505.9763 | 4782.0032 | 2.295787e+04 |
| 2000 | OLS regression | 2000 | Linear | 10783.4475 | 5284.1931 | 2.600208e+04 |
| 2001 | OLS regression | 2001 | Linear | 12232.5931 | 5836.0655 | 2.943429e+04 |
| 2002 | OLS regression | 2002 | Linear | 13876.4837 | 6442.5544 | 3.332500e+04 |
| 2003 | OLS regression | 2003 | Linear | 15741.2904 | 7109.0797 | 3.776617e+04 |
| 2004 | OLS regression | 2004 | Linear | 17856.7012 | 7841.5957 | 4.287830e+04 |
| 2005 | OLS regression | 2005 | Linear | 20256.3939 | 8646.6423 | 4.881959e+04 |
| 2006 | OLS regression | 2006 | Linear | 22978.5719 | 9531.4031 | 5.579766e+04 |
| 2007 | OLS regression | 2007 | Linear | 26066.5729 | 10503.7682 | 6.408442e+04 |
| 2008 | OLS regression | 2008 | Linear | 29569.5583 | 11572.4037 | 7.403508e+04 |
| 2009 | OLS regression | 2009 | Linear | 33543.2963 | 12746.8280 | 8.611251e+04 |
| 2010 | OLS regression | 2010 | Linear | 38051.0495 | 14037.4952 | 1.009189e+05 |
| 2011 | OLS regression | 2011 | Linear | 43164.5821 | 15455.8874 | 1.192378e+05 |
| 2012 | OLS regression | 2012 | Linear | 48965.3025 | 17014.6157 | 1.420902e+05 |
| 2013 | OLS regression | 2013 | Linear | 55545.5592 | 18727.5307 | 1.708110e+05 |
| 2014 | OLS regression | 2014 | Linear | 63010.1110 | 20609.8450 | 2.071525e+05 |
| 2015 | OLS regression | 2015 | Linear | 71477.7948 | 22678.2661 | 2.534267e+05 |
| 2016 | OLS regression | 2016 | Linear | 81083.4178 | 24951.1440 | 3.126997e+05 |
| 2017 | OLS regression | 2017 | Linear | 91979.9032 | 27448.6323 | 3.890589e+05 |
| 19701 | OLS regression | 1970 | Quadratic | 264.7175 | 111.5784 | 7.321475e+02 |
| 19711 | OLS regression | 1971 | Quadratic | 297.2102 | 138.9244 | 7.526587e+02 |
| 19721 | OLS regression | 1972 | Quadratic | 333.8513 | 171.4500 | 7.793883e+02 |
| 19731 | OLS regression | 1973 | Quadratic | 375.1897 | 209.4027 | 8.132671e+02 |
| 19741 | OLS regression | 1974 | Quadratic | 421.8491 | 252.5847 | 8.555056e+02 |
| 19751 | OLS regression | 1975 | Quadratic | 474.5387 | 300.1256 | 9.076508e+02 |
| 19761 | OLS regression | 1976 | Quadratic | 534.0655 | 350.4063 | 9.716403e+02 |
| 19771 | OLS regression | 1977 | Quadratic | 601.3479 | 401.4083 | 1.049840e+03 |
| 19781 | OLS regression | 1978 | Quadratic | 677.4316 | 451.5286 | 1.145044e+03 |
| 19791 | OLS regression | 1979 | Quadratic | 763.5077 | 500.2949 | 1.260435e+03 |
| 19801 | OLS regression | 1980 | Quadratic | 860.9338 | 548.3866 | 1.399490e+03 |
| 19811 | OLS regression | 1981 | Quadratic | 971.2578 | 597.1441 | 1.565889e+03 |
| 19821 | OLS regression | 1982 | Quadratic | 1096.2449 | 648.1082 | 1.763444e+03 |
| 19831 | OLS regression | 1983 | Quadratic | 1237.9099 | 702.8044 | 1.996071e+03 |
| 19841 | OLS regression | 1984 | Quadratic | 1398.5527 | 762.7089 | 2.267825e+03 |
| 19851 | OLS regression | 1985 | Quadratic | 1580.8003 | 829.2909 | 2.582953e+03 |
| 19861 | OLS regression | 1986 | Quadratic | 1787.6544 | 904.0776 | 2.945950e+03 |
| 19871 | OLS regression | 1987 | Quadratic | 2022.5462 | 988.7192 | 3.361620e+03 |
| 19881 | OLS regression | 1988 | Quadratic | 2289.4002 | 1085.0535 | 3.835120e+03 |
| 19891 | OLS regression | 1989 | Quadratic | 2592.7064 | 1195.1669 | 4.372023e+03 |
| 19901 | OLS regression | 1990 | Quadratic | 2937.6046 | 1321.4582 | 4.978381e+03 |
| 19911 | OLS regression | 1991 | Quadratic | 3329.9805 | 1466.7046 | 5.660834e+03 |
| 19921 | OLS regression | 1992 | Quadratic | 3776.5776 | 1634.1307 | 6.426764e+03 |
| 19931 | OLS regression | 1993 | Quadratic | 4285.1250 | 1827.4824 | 7.284541e+03 |
| 19941 | OLS regression | 1994 | Quadratic | 4864.4858 | 2051.1007 | 8.243881e+03 |
| 19951 | OLS regression | 1995 | Quadratic | 5524.8279 | 2309.9949 | 9.316393e+03 |
| 19961 | OLS regression | 1996 | Quadratic | 6277.8210 | 2609.9044 | 1.051638e+04 |
| 19971 | OLS regression | 1997 | Quadratic | 7136.8647 | 2957.3418 | 1.186199e+04 |
| 19981 | OLS regression | 1998 | Quadratic | 8117.3517 | 3359.5958 | 1.337692e+04 |
| 19991 | OLS regression | 1999 | Quadratic | 9236.9721 | 3824.6767 | 1.509278e+04 |
| 20001 | OLS regression | 2000 | Quadratic | 10516.0649 | 4361.1670 | 1.705234e+04 |
| 20011 | OLS regression | 2001 | Quadratic | 11978.0260 | 4977.9468 | 1.931400e+04 |
| 20021 | OLS regression | 2002 | Quadratic | 13649.7786 | 5683.7491 | 2.195733e+04 |
| 20031 | OLS regression | 2003 | Quadratic | 15562.3196 | 6486.5192 | 2.509002e+04 |
| 20041 | OLS regression | 2004 | Quadratic | 17751.3511 | 7392.5638 | 2.885568e+04 |
| 20051 | OLS regression | 2005 | Quadratic | 20258.0137 | 8405.5246 | 3.344289e+04 |
| 20061 | OLS regression | 2006 | Quadratic | 23129.7358 | 9525.2795 | 3.909557e+04 |
| 20071 | OLS regression | 2007 | Quadratic | 26421.2187 | 10746.9274 | 4.612645e+04 |
| 20081 | OLS regression | 2008 | Quadratic | 30195.5806 | 12060.0926 | 5.493566e+04 |
| 20091 | OLS regression | 2009 | Quadratic | 34525.6837 | 13448.7448 | 6.603732e+04 |
| 20101 | OLS regression | 2010 | Quadratic | 39495.6762 | 14891.6768 | 8.009683e+04 |
| 20111 | OLS regression | 2011 | Quadratic | 45202.7840 | 16363.5946 | 9.798226e+04 |
| 20121 | OLS regression | 2012 | Quadratic | 51759.3933 | 17836.6525 | 1.208342e+05 |
| 20131 | OLS regression | 2013 | Quadratic | 59295.4725 | 19282.1294 | 1.501606e+05 |
| 20141 | OLS regression | 2014 | Quadratic | 67961.3904 | 20671.9899 | 1.879655e+05 |
| 20151 | OLS regression | 2015 | Quadratic | 77931.1957 | 21980.1199 | 2.369258e+05 |
| 20161 | OLS regression | 2016 | Quadratic | 89406.4371 | 23183.1979 | 3.006327e+05 |
| 20171 | OLS regression | 2017 | Quadratic | 102620.6120 | 24261.1968 | 3.839266e+05 |
| 19702 | Robust regression | 1970 | Linear | 239.4078 | 159.2535 | 3.599048e+02 |
| 19712 | Robust regression | 1971 | Linear | 267.8488 | 180.3775 | 3.977381e+02 |
| 19722 | Robust regression | 1972 | Linear | 299.6686 | 204.1476 | 4.398841e+02 |
| 19732 | Robust regression | 1973 | Linear | 335.2685 | 230.8575 | 4.869018e+02 |
| 19742 | Robust regression | 1974 | Linear | 375.0975 | 260.8253 | 5.394343e+02 |
| 19752 | Robust regression | 1975 | Linear | 419.6581 | 294.3938 | 5.982223e+02 |
| 19762 | Robust regression | 1976 | Linear | 469.5124 | 331.9313 | 6.641191e+02 |
| 19772 | Robust regression | 1977 | Linear | 525.2893 | 373.8325 | 7.381082e+02 |
| 19782 | Robust regression | 1978 | Linear | 587.6923 | 420.5193 | 8.213231e+02 |
| 19792 | Robust regression | 1979 | Linear | 657.5087 | 472.4430 | 9.150684e+02 |
| 19802 | Robust regression | 1980 | Linear | 735.6190 | 530.0862 | 1.020844e+03 |
| 19812 | Robust regression | 1981 | Linear | 823.0087 | 593.9669 | 1.140372e+03 |
| 19822 | Robust regression | 1982 | Linear | 920.7800 | 664.6431 | 1.275626e+03 |
| 19832 | Robust regression | 1983 | Linear | 1030.1663 | 742.7188 | 1.428862e+03 |
| 19842 | Robust regression | 1984 | Linear | 1152.5474 | 828.8506 | 1.602660e+03 |
| 19852 | Robust regression | 1985 | Linear | 1289.4671 | 923.7562 | 1.799961e+03 |
| 19862 | Robust regression | 1986 | Linear | 1442.6525 | 1028.2227 | 2.024120e+03 |
| 19872 | Robust regression | 1987 | Linear | 1614.0359 | 1143.1158 | 2.278957e+03 |
| 19882 | Robust regression | 1988 | Linear | 1805.7793 | 1269.3895 | 2.568824e+03 |
| 19892 | Robust regression | 1989 | Linear | 2020.3012 | 1408.0956 | 2.898679e+03 |
| 19902 | Robust regression | 1990 | Linear | 2260.3078 | 1560.3941 | 3.274167e+03 |
| 19912 | Robust regression | 1991 | Linear | 2528.8266 | 1727.5638 | 3.701724e+03 |
| 19922 | Robust regression | 1992 | Linear | 2829.2448 | 1911.0137 | 4.188681e+03 |
| 19932 | Robust regression | 1993 | Linear | 3165.3518 | 2112.2946 | 4.743397e+03 |
| 19942 | Robust regression | 1994 | Linear | 3541.3875 | 2333.1126 | 5.375405e+03 |
| 19952 | Robust regression | 1995 | Linear | 3962.0953 | 2575.3428 | 6.095576e+03 |
| 19962 | Robust regression | 1996 | Linear | 4432.7821 | 2841.0454 | 6.916312e+03 |
| 19972 | Robust regression | 1997 | Linear | 4959.3853 | 3132.4821 | 7.851762e+03 |
| 19982 | Robust regression | 1998 | Linear | 5548.5475 | 3452.1351 | 8.918069e+03 |
| 19992 | Robust regression | 1999 | Linear | 6207.7008 | 3802.7279 | 1.013366e+04 |
| 20002 | Robust regression | 2000 | Linear | 6945.1598 | 4187.2478 | 1.151956e+04 |
| 20012 | Robust regression | 2001 | Linear | 7770.2270 | 4608.9704 | 1.309976e+04 |
| 20022 | Robust regression | 2002 | Linear | 8693.3101 | 5071.4877 | 1.490167e+04 |
| 20032 | Robust regression | 2003 | Linear | 9726.0531 | 5578.7377 | 1.695654e+04 |
| 20042 | Robust regression | 2004 | Linear | 10881.4833 | 6135.0373 | 1.930007e+04 |
| 20052 | Robust regression | 2005 | Linear | 12174.1757 | 6745.1187 | 2.197301e+04 |
| 20062 | Robust regression | 2006 | Linear | 13620.4366 | 7414.1691 | 2.502186e+04 |
| 20072 | Robust regression | 2007 | Linear | 15238.5097 | 8147.8740 | 2.849973e+04 |
| 20082 | Robust regression | 2008 | Linear | 17048.8057 | 8952.4652 | 3.246723e+04 |
| 20092 | Robust regression | 2009 | Linear | 19074.1602 | 9834.7727 | 3.699359e+04 |
| 20102 | Robust regression | 2010 | Linear | 21340.1216 | 10802.2826 | 4.215783e+04 |
| 20112 | Robust regression | 2011 | Linear | 23875.2735 | 11863.1994 | 4.805016e+04 |
| 20122 | Robust regression | 2012 | Linear | 26711.5949 | 13026.5152 | 5.477361e+04 |
| 20132 | Robust regression | 2013 | Linear | 29884.8640 | 14302.0850 | 6.244580e+04 |
| 20142 | Robust regression | 2014 | Linear | 33435.1095 | 15700.7096 | 7.120102e+04 |
| 20152 | Robust regression | 2015 | Linear | 37407.1149 | 17234.2259 | 8.119264e+04 |
| 20162 | Robust regression | 2016 | Linear | 41850.9845 | 18915.6065 | 9.259576e+04 |
| 20172 | Robust regression | 2017 | Linear | 46822.7744 | 20759.0685 | 1.056103e+05 |
| 19703 | Robust regression | 1970 | Quadratic | 372.5585 | 198.7284 | 6.984401e+02 |
| 19713 | Robust regression | 1971 | Quadratic | 390.2582 | 220.0625 | 6.920825e+02 |
| 19723 | Robust regression | 1972 | Quadratic | 410.0358 | 243.3672 | 6.908463e+02 |
| 19733 | Robust regression | 1973 | Quadratic | 432.1194 | 268.7127 | 6.948954e+02 |
| 19743 | Robust regression | 1974 | Quadratic | 456.7706 | 296.1345 | 7.045427e+02 |
| 19753 | Robust regression | 1975 | Quadratic | 484.2892 | 325.6267 | 7.202607e+02 |
| 19763 | Robust regression | 1976 | Quadratic | 515.0196 | 357.1433 | 7.426857e+02 |
| 19773 | Robust regression | 1977 | Quadratic | 549.3575 | 390.6153 | 7.726109e+02 |
| 19783 | Robust regression | 1978 | Quadratic | 587.7581 | 425.9885 | 8.109598e+02 |
| 19793 | Robust regression | 1979 | Quadratic | 630.7460 | 463.2818 | 8.587440e+02 |
| 19803 | Robust regression | 1980 | Quadratic | 678.9263 | 502.6532 | 9.170160e+02 |
| 19813 | Robust regression | 1981 | Quadratic | 732.9986 | 544.4510 | 9.868417e+02 |
| 19823 | Robust regression | 1982 | Quadratic | 793.7722 | 589.2365 | 1.069306e+03 |
| 19833 | Robust regression | 1983 | Quadratic | 862.1860 | 637.7765 | 1.165557e+03 |
| 19843 | Robust regression | 1984 | Quadratic | 939.3304 | 691.0235 | 1.276862e+03 |
| 19853 | Robust regression | 1985 | Quadratic | 1026.4742 | 750.0991 | 1.404680e+03 |
| 19863 | Robust regression | 1986 | Quadratic | 1125.0971 | 816.2940 | 1.550720e+03 |
| 19873 | Robust regression | 1987 | Quadratic | 1236.9277 | 891.0853 | 1.716996e+03 |
| 19883 | Robust regression | 1988 | Quadratic | 1363.9890 | 976.1695 | 1.905884e+03 |
| 19893 | Robust regression | 1989 | Quadratic | 1508.6544 | 1073.5098 | 2.120184e+03 |
| 19903 | Robust regression | 1990 | Quadratic | 1673.7129 | 1185.3956 | 2.363190e+03 |
| 19913 | Robust regression | 1991 | Quadratic | 1862.4492 | 1314.5124 | 2.638786e+03 |
| 19923 | Robust regression | 1992 | Quadratic | 2078.7403 | 1464.0239 | 2.951565e+03 |
| 19933 | Robust regression | 1993 | Quadratic | 2327.1711 | 1637.6644 | 3.306981e+03 |
| 19943 | Robust regression | 1994 | Quadratic | 2613.1763 | 1839.8417 | 3.711564e+03 |
| 19953 | Robust regression | 1995 | Quadratic | 2943.2111 | 2075.7475 | 4.173191e+03 |
| 19963 | Robust regression | 1996 | Quadratic | 3324.9599 | 2351.4700 | 4.701467e+03 |
| 19973 | Robust regression | 1997 | Quadratic | 3767.5906 | 2674.0971 | 5.308236e+03 |
| 19983 | Robust regression | 1998 | Quadratic | 4282.0655 | 3051.7961 | 6.008293e+03 |
| 19993 | Robust regression | 1999 | Quadratic | 4881.5215 | 3493.8459 | 6.820350e+03 |
| 20003 | Robust regression | 2000 | Quadratic | 5581.7375 | 4010.5946 | 7.768373e+03 |
| 20013 | Robust regression | 2001 | Quadratic | 6401.7088 | 4613.3177 | 8.883384e+03 |
| 20023 | Robust regression | 2002 | Quadratic | 7364.3550 | 5313.9701 | 1.020588e+04 |
| 20033 | Robust regression | 2003 | Quadratic | 8497.3954 | 6124.8695 | 1.178894e+04 |
| 20043 | Robust regression | 2004 | Quadratic | 9834.4309 | 7058.4003 | 1.370226e+04 |
| 20053 | Robust regression | 2005 | Quadratic | 11416.2888 | 8126.8906 | 1.603709e+04 |
| 20063 | Robust regression | 2006 | Quadratic | 13292.6925 | 9342.7900 | 1.891252e+04 |
| 20073 | Robust regression | 2007 | Quadratic | 15524.3445 | 10719.1991 | 2.248351e+04 |
| 20083 | Robust regression | 2008 | Quadratic | 18185.5269 | 12270.6416 | 2.695160e+04 |
| 20093 | Robust regression | 2009 | Quadratic | 21367.3569 | 14013.8870 | 3.257939e+04 |
| 20103 | Robust regression | 2010 | Quadratic | 25181.8727 | 15968.6552 | 3.971071e+04 |
| 20113 | Robust regression | 2011 | Quadratic | 29767.1698 | 18158.1331 | 4.879821e+04 |
| 20123 | Robust regression | 2012 | Quadratic | 35293.8769 | 20609.3419 | 6.044141e+04 |
| 20133 | Robust regression | 2013 | Quadratic | 41973.3362 | 23353.4329 | 7.543906e+04 |
| 20143 | Robust regression | 2014 | Quadratic | 50067.9622 | 26425.9896 | 9.486119e+04 |
| 20153 | Robust regression | 2015 | Quadratic | 59904.3890 | 29867.3896 | 1.201490e+05 |
| 20163 | Robust regression | 2016 | Quadratic | 71890.1966 | 33723.2459 | 1.532533e+05 |
| 20173 | Robust regression | 2017 | Quadratic | 86535.2399 | 38044.9517 | 1.968289e+05 |
| 19704 | MARS | 1970 | 725.6291 | 145.5584 | 3.617362e+03 | |
| 19714 | MARS | 1971 | 725.6291 | 145.5584 | 3.617362e+03 | |
| 19724 | MARS | 1972 | 725.6291 | 145.5584 | 3.617362e+03 | |
| 19734 | MARS | 1973 | 725.6291 | 145.5584 | 3.617362e+03 | |
| 19744 | MARS | 1974 | 725.6291 | 145.5584 | 3.617362e+03 | |
| 19754 | MARS | 1975 | 725.6291 | 145.5584 | 3.617362e+03 | |
| 19764 | MARS | 1976 | 725.6291 | 145.5584 | 3.617362e+03 | |
| 19774 | MARS | 1977 | 725.6291 | 145.5584 | 3.617362e+03 | |
| 19784 | MARS | 1978 | 725.6291 | 145.5584 | 3.617362e+03 | |
| 19794 | MARS | 1979 | 725.6291 | 145.5584 | 3.617362e+03 | |
| 19804 | MARS | 1980 | 725.6291 | 145.5584 | 3.617362e+03 | |
| 19814 | MARS | 1981 | 725.6291 | 145.5584 | 3.617362e+03 | |
| 19824 | MARS | 1982 | 725.6291 | 145.5584 | 3.617362e+03 | |
| 19834 | MARS | 1983 | 725.6291 | 145.5584 | 3.617362e+03 | |
| 19844 | MARS | 1984 | 725.6291 | 145.5584 | 3.617362e+03 | |
| 19854 | MARS | 1985 | 725.6291 | 145.5584 | 3.617362e+03 | |
| 19864 | MARS | 1986 | 725.6291 | 145.5584 | 3.617362e+03 | |
| 19874 | MARS | 1987 | 725.6291 | 145.5584 | 3.617362e+03 | |
| 19884 | MARS | 1988 | 725.6291 | 145.5584 | 3.617362e+03 | |
| 19894 | MARS | 1989 | 725.6291 | 145.5584 | 3.617362e+03 | |
| 19904 | MARS | 1990 | 1025.3657 | 183.1066 | 5.741871e+03 | |
| 19914 | MARS | 1991 | 1448.9149 | 230.3408 | 9.114125e+03 | |
| 19924 | MARS | 1992 | 2047.4201 | 289.7594 | 1.446693e+04 | |
| 19934 | MARS | 1993 | 2893.1507 | 364.5056 | 2.296349e+04 | |
| 19944 | MARS | 1994 | 4088.2284 | 458.5333 | 3.645016e+04 | |
| 19954 | MARS | 1995 | 5776.9585 | 576.8163 | 5.785767e+04 | |
| 19964 | MARS | 1996 | 8163.2546 | 725.6116 | 9.183801e+04 | |
| 19974 | MARS | 1997 | 11535.2612 | 912.7900 | 1.457753e+05 | |
| 19984 | MARS | 1998 | 16300.1472 | 1148.2530 | 2.313905e+05 | |
| 19994 | MARS | 1999 | 23033.2712 | 1444.4558 | 3.672882e+05 | |
| 20004 | MARS | 2000 | 32547.6557 | 1817.0670 | 5.829999e+05 | |
| 20014 | MARS | 2001 | 45992.1601 | 2285.7968 | 9.254011e+05 | |
| 20024 | MARS | 2002 | 64990.2040 | 2875.4400 | 1.468898e+06 | |
| 20034 | MARS | 2003 | 91835.7958 | 3617.1872 | 2.331594e+06 | |
| 20044 | MARS | 2004 | 129770.5325 | 4550.2753 | 3.700961e+06 | |
| 20054 | MARS | 2005 | 49463.1801 | 2398.8879 | 1.019892e+06 | |
| 20064 | MARS | 2006 | 18853.3263 | 1264.6846 | 2.810566e+05 | |
| 20074 | MARS | 2007 | 7186.1112 | 666.7369 | 7.745213e+04 | |
| 20084 | MARS | 2008 | 9817.2147 | 820.1336 | 1.175146e+05 | |
| 20094 | MARS | 2009 | 13411.6634 | 1008.8224 | 1.782997e+05 | |
| 20104 | MARS | 2010 | 18322.1740 | 1240.9230 | 2.705261e+05 | |
| 20114 | MARS | 2011 | 25030.6058 | 1526.4231 | 4.104571e+05 | |
| 20124 | MARS | 2012 | 34195.2449 | 1877.6085 | 6.227681e+05 | |
| 20134 | MARS | 2013 | 46715.4004 | 2309.5914 | 9.448981e+05 | |
| 20144 | MARS | 2014 | 63819.6521 | 2840.9610 | 1.433652e+06 | |
| 20154 | MARS | 2015 | 87186.4088 | 3494.5831 | 2.175215e+06 | |
| 20164 | MARS | 2016 | 119108.6073 | 4298.5847 | 3.300356e+06 | |
| 20174 | MARS | 2017 | 162718.7143 | 5287.5635 | 5.007482e+06 | |
| 19705 | GAM | 1970 | 414.8239 | 237.4243 | 7.247735e+02 | |
| 19715 | GAM | 1971 | 427.5848 | 264.4702 | 6.913020e+02 | |
| 19725 | GAM | 1972 | 441.4140 | 290.7849 | 6.700705e+02 | |
| 19735 | GAM | 1973 | 457.2298 | 315.0350 | 6.636058e+02 | |
| 19745 | GAM | 1974 | 475.9957 | 337.0230 | 6.722744e+02 | |
| 19755 | GAM | 1975 | 498.5025 | 357.7193 | 6.946921e+02 | |
| 19765 | GAM | 1976 | 525.2333 | 378.3766 | 7.290883e+02 | |
| 19775 | GAM | 1977 | 556.3217 | 399.9247 | 7.738802e+02 | |
| 19785 | GAM | 1978 | 591.6096 | 423.0964 | 8.272390e+02 | |
| 19795 | GAM | 1979 | 630.8037 | 448.6460 | 8.869205e+02 | |
| 19805 | GAM | 1980 | 673.7168 | 477.2428 | 9.510764e+02 | |
| 19815 | GAM | 1981 | 720.5464 | 509.3203 | 1.019372e+03 | |
| 19825 | GAM | 1982 | 772.1292 | 545.2740 | 1.093365e+03 | |
| 19835 | GAM | 1983 | 830.1110 | 585.9597 | 1.175993e+03 | |
| 19845 | GAM | 1984 | 897.0064 | 633.0339 | 1.271054e+03 | |
| 19855 | GAM | 1985 | 976.1635 | 688.7847 | 1.383444e+03 | |
| 19865 | GAM | 1986 | 1071.6753 | 755.6616 | 1.519844e+03 | |
| 19875 | GAM | 1987 | 1188.2777 | 836.0164 | 1.688967e+03 | |
| 19885 | GAM | 1988 | 1331.2647 | 932.2403 | 1.901082e+03 | |
| 19895 | GAM | 1989 | 1506.4280 | 1046.9158 | 2.167629e+03 | |
| 19905 | GAM | 1990 | 1720.0295 | 1182.5398 | 2.501820e+03 | |
| 19915 | GAM | 1991 | 1978.8127 | 1341.0386 | 2.919901e+03 | |
| 19925 | GAM | 1992 | 2290.0631 | 1523.8433 | 3.441554e+03 | |
| 19935 | GAM | 1993 | 2661.7123 | 1732.7998 | 4.088592e+03 | |
| 19945 | GAM | 1994 | 3102.4586 | 1971.1106 | 4.883161e+03 | |
| 19955 | GAM | 1995 | 3621.8395 | 2243.3868 | 5.847285e+03 | |
| 19965 | GAM | 1996 | 4230.1808 | 2555.0557 | 7.003538e+03 | |
| 19975 | GAM | 1997 | 4938.3689 | 2912.2289 | 8.374166e+03 | |
| 19985 | GAM | 1998 | 5757.4611 | 3322.5139 | 9.976891e+03 | |
| 19995 | GAM | 1999 | 6698.2527 | 3795.9214 | 1.181968e+04 | |
| 20005 | GAM | 2000 | 7771.0285 | 4344.7404 | 1.389931e+04 | |
| 20015 | GAM | 2001 | 8985.7809 | 4982.6175 | 1.620519e+04 | |
| 20025 | GAM | 2002 | 10353.1456 | 5724.5558 | 1.872418e+04 | |
| 20035 | GAM | 2003 | 11886.1528 | 6589.1124 | 2.144153e+04 | |
| 20045 | GAM | 2004 | 13602.6642 | 7601.5706 | 2.434135e+04 | |
| 20055 | GAM | 2005 | 15528.1342 | 8794.8982 | 2.741623e+04 | |
| 20065 | GAM | 2006 | 17698.1930 | 10206.7954 | 3.068799e+04 | |
| 20075 | GAM | 2007 | 20160.5461 | 11875.1018 | 3.422687e+04 | |
| 20085 | GAM | 2008 | 22975.7956 | 13834.8296 | 3.815639e+04 | |
| 20095 | GAM | 2009 | 26216.9492 | 16113.1915 | 4.265626e+04 | |
| 20105 | GAM | 2010 | 29967.5734 | 18708.6613 | 4.800212e+04 | |
| 20115 | GAM | 2011 | 34318.8017 | 21546.9466 | 5.466112e+04 | |
| 20125 | GAM | 2012 | 39365.8787 | 24450.5825 | 6.337977e+04 | |
| 20135 | GAM | 2013 | 45205.7855 | 27196.2110 | 7.514146e+04 | |
| 20145 | GAM | 2014 | 51938.9112 | 29645.4218 | 9.099720e+04 | |
| 20155 | GAM | 2015 | 59679.7875 | 31798.8566 | 1.120065e+05 | |
| 20165 | GAM | 2016 | 68574.3493 | 33719.6420 | 1.394570e+05 | |
| 20175 | GAM | 2017 | 78794.5396 | 35470.8925 | 1.750331e+05 | |
| 19706 | Quantile regression | 1970 | Quantile 0.1 | 150.3668 | 119.0731 | 1.898849e+02 |
| 19716 | Quantile regression | 1971 | Quantile 0.1 | 166.1144 | 132.5982 | 2.081024e+02 |
| 19726 | Quantile regression | 1972 | Quantile 0.1 | 183.5113 | 147.6423 | 2.280944e+02 |
| 19736 | Quantile regression | 1973 | Quantile 0.1 | 202.7300 | 164.3718 | 2.500396e+02 |
| 19746 | Quantile regression | 1974 | Quantile 0.1 | 223.9615 | 182.9703 | 2.741361e+02 |
| 19756 | Quantile regression | 1975 | Quantile 0.1 | 247.4165 | 203.6402 | 3.006035e+02 |
| 19766 | Quantile regression | 1976 | Quantile 0.1 | 273.3280 | 226.6040 | 3.296861e+02 |
| 19776 | Quantile regression | 1977 | Quantile 0.1 | 301.9530 | 252.1059 | 3.616560e+02 |
| 19786 | Quantile regression | 1978 | Quantile 0.1 | 333.5760 | 280.4139 | 3.968167e+02 |
| 19796 | Quantile regression | 1979 | Quantile 0.1 | 368.5107 | 311.8205 | 4.355073e+02 |
| 19806 | Quantile regression | 1980 | Quantile 0.1 | 407.1040 | 346.6450 | 4.781079e+02 |
| 19816 | Quantile regression | 1981 | Quantile 0.1 | 449.7392 | 385.2344 | 5.250449e+02 |
| 19826 | Quantile regression | 1982 | Quantile 0.1 | 496.8395 | 427.9651 | 5.767981e+02 |
| 19836 | Quantile regression | 1983 | Quantile 0.1 | 548.8724 | 475.2437 | 6.339083e+02 |
| 19846 | Quantile regression | 1984 | Quantile 0.1 | 606.3547 | 527.5089 | 6.969855e+02 |
| 19856 | Quantile regression | 1985 | Quantile 0.1 | 669.8570 | 585.2318 | 7.667190e+02 |
| 19866 | Quantile regression | 1986 | Quantile 0.1 | 740.0097 | 648.9183 | 8.438880e+02 |
| 19876 | Quantile regression | 1987 | Quantile 0.1 | 817.5094 | 719.1104 | 9.293727e+02 |
| 19886 | Quantile regression | 1988 | Quantile 0.1 | 903.1255 | 796.3897 | 1.024166e+03 |
| 19896 | Quantile regression | 1989 | Quantile 0.1 | 997.7079 | 881.3808 | 1.129388e+03 |
| 19906 | Quantile regression | 1990 | Quantile 0.1 | 1102.1958 | 974.7578 | 1.246295e+03 |
| 19916 | Quantile regression | 1991 | Quantile 0.1 | 1217.6265 | 1077.2505 | 1.376295e+03 |
| 19926 | Quantile regression | 1992 | Quantile 0.1 | 1345.1460 | 1189.6540 | 1.520961e+03 |
| 19936 | Quantile regression | 1993 | Quantile 0.1 | 1486.0203 | 1312.8385 | 1.682047e+03 |
| 19946 | Quantile regression | 1994 | Quantile 0.1 | 1641.6482 | 1447.7601 | 1.861502e+03 |
| 19956 | Quantile regression | 1995 | Quantile 0.1 | 1813.5746 | 1595.4722 | 2.061492e+03 |
| 19966 | Quantile regression | 1996 | Quantile 0.1 | 2003.5065 | 1757.1368 | 2.284420e+03 |
| 19976 | Quantile regression | 1997 | Quantile 0.1 | 2213.3296 | 1934.0360 | 2.532956e+03 |
| 19986 | Quantile regression | 1998 | Quantile 0.1 | 2445.1270 | 2127.5834 | 2.810064e+03 |
| 19996 | Quantile regression | 1999 | Quantile 0.1 | 2701.2001 | 2339.3356 | 3.119040e+03 |
| 20006 | Quantile regression | 2000 | Quantile 0.1 | 2984.0912 | 2571.0048 | 3.463549e+03 |
| 20016 | Quantile regression | 2001 | Quantile 0.1 | 3296.6089 | 2824.4723 | 3.847667e+03 |
| 20026 | Quantile regression | 2002 | Quantile 0.1 | 3641.8559 | 3101.8024 | 4.275938e+03 |
| 20036 | Quantile regression | 2003 | Quantile 0.1 | 4023.2599 | 3405.2589 | 4.753418e+03 |
| 20046 | Quantile regression | 2004 | Quantile 0.1 | 4444.6076 | 3737.3230 | 5.285745e+03 |
| 20056 | Quantile regression | 2005 | Quantile 0.1 | 4910.0822 | 4100.7121 | 5.879200e+03 |
| 20066 | Quantile regression | 2006 | Quantile 0.1 | 5424.3049 | 4498.4023 | 6.540785e+03 |
| 20076 | Quantile regression | 2007 | Quantile 0.1 | 5992.3812 | 4933.6519 | 7.278307e+03 |
| 20086 | Quantile regression | 2008 | Quantile 0.1 | 6619.9508 | 5410.0272 | 8.100467e+03 |
| 20096 | Quantile regression | 2009 | Quantile 0.1 | 7313.2446 | 5931.4322 | 9.016970e+03 |
| 20106 | Quantile regression | 2010 | Quantile 0.1 | 8079.1455 | 6502.1396 | 1.003863e+04 |
| 20116 | Quantile regression | 2011 | Quantile 0.1 | 8925.2577 | 7126.8258 | 1.117752e+04 |
| 20126 | Quantile regression | 2012 | Quantile 0.1 | 9859.9815 | 7810.6089 | 1.244708e+04 |
| 20136 | Quantile regression | 2013 | Quantile 0.1 | 10892.5970 | 8559.0903 | 1.386230e+04 |
| 20146 | Quantile regression | 2014 | Quantile 0.1 | 12033.3561 | 9378.4004 | 1.543991e+04 |
| 20156 | Quantile regression | 2015 | Quantile 0.1 | 13293.5846 | 10275.2488 | 1.719855e+04 |
| 20166 | Quantile regression | 2016 | Quantile 0.1 | 14685.7941 | 11256.9786 | 1.915901e+04 |
| 20176 | Quantile regression | 2017 | Quantile 0.1 | 16223.8070 | 12331.6266 | 2.134446e+04 |
| 19707 | Quantile regression | 1970 | Quantile 0.5 | 215.3841 | 132.2347 | 3.508181e+02 |
| 19717 | Quantile regression | 1971 | Quantile 0.5 | 242.5712 | 152.2634 | 3.864408e+02 |
| 19727 | Quantile regression | 1972 | Quantile 0.5 | 273.1900 | 175.1637 | 4.260743e+02 |
| 19737 | Quantile regression | 1973 | Quantile 0.5 | 307.6736 | 201.2941 | 4.702724e+02 |
| 19747 | Quantile regression | 1974 | Quantile 0.5 | 346.5100 | 231.0400 | 5.196900e+02 |
| 19757 | Quantile regression | 1975 | Quantile 0.5 | 390.2486 | 264.8103 | 5.751059e+02 |
| 19767 | Quantile regression | 1976 | Quantile 0.5 | 439.5081 | 303.0315 | 6.374498e+02 |
| 19777 | Quantile regression | 1977 | Quantile 0.5 | 494.9854 | 346.1409 | 7.078348e+02 |
| 19787 | Quantile regression | 1978 | Quantile 0.5 | 557.4655 | 394.5787 | 7.875939e+02 |
| 19797 | Quantile regression | 1979 | Quantile 0.5 | 627.8321 | 448.7807 | 8.783201e+02 |
| 19807 | Quantile regression | 1980 | Quantile 0.5 | 707.0808 | 509.1743 | 9.819097e+02 |
| 19817 | Quantile regression | 1981 | Quantile 0.5 | 796.3327 | 576.1790 | 1.100606e+03 |
| 19827 | Quantile regression | 1982 | Quantile 0.5 | 896.8506 | 650.2139 | 1.237041e+03 |
| 19837 | Quantile regression | 1983 | Quantile 0.5 | 1010.0565 | 731.7128 | 1.394282e+03 |
| 19847 | Quantile regression | 1984 | Quantile 0.5 | 1137.5518 | 821.1444 | 1.575879e+03 |
| 19857 | Quantile regression | 1985 | Quantile 0.5 | 1281.1404 | 919.0339 | 1.785920e+03 |
| 19867 | Quantile regression | 1986 | Quantile 0.5 | 1442.8536 | 1025.9846 | 2.029101e+03 |
| 19877 | Quantile regression | 1987 | Quantile 0.5 | 1624.9792 | 1142.6941 | 2.310817e+03 |
| 19887 | Quantile regression | 1988 | Quantile 0.5 | 1830.0938 | 1269.9664 | 2.637269e+03 |
| 19897 | Quantile regression | 1989 | Quantile 0.5 | 2061.0992 | 1408.7194 | 3.015597e+03 |
| 19907 | Quantile regression | 1990 | Quantile 0.5 | 2321.2634 | 1559.9903 | 3.454037e+03 |
| 19917 | Quantile regression | 1991 | Quantile 0.5 | 2614.2671 | 1724.9394 | 3.962106e+03 |
| 19927 | Quantile regression | 1992 | Quantile 0.5 | 2944.2556 | 1904.8552 | 4.550813e+03 |
| 19937 | Quantile regression | 1993 | Quantile 0.5 | 3315.8971 | 2101.1601 | 5.232906e+03 |
| 19947 | Quantile regression | 1994 | Quantile 0.5 | 3734.4494 | 2315.4181 | 6.023151e+03 |
| 19957 | Quantile regression | 1995 | Quantile 0.5 | 4205.8339 | 2549.3443 | 6.938662e+03 |
| 19967 | Quantile regression | 1996 | Quantile 0.5 | 4736.7194 | 2804.8162 | 7.999280e+03 |
| 19977 | Quantile regression | 1997 | Quantile 0.5 | 5334.6164 | 3083.8876 | 9.228006e+03 |
| 19987 | Quantile regression | 1998 | Quantile 0.5 | 6007.9836 | 3388.8030 | 1.065151e+04 |
| 19997 | Quantile regression | 1999 | Quantile 0.5 | 6766.3472 | 3722.0155 | 1.230072e+04 |
| 20007 | Quantile regression | 2000 | Quantile 0.5 | 7620.4359 | 4086.2050 | 1.421149e+04 |
| 20017 | Quantile regression | 2001 | Quantile 0.5 | 8582.3329 | 4484.2999 | 1.642540e+04 |
| 20027 | Quantile regression | 2002 | Quantile 0.5 | 9665.6462 | 4919.4998 | 1.899069e+04 |
| 20037 | Quantile regression | 2003 | Quantile 0.5 | 10885.7018 | 5395.3009 | 2.196328e+04 |
| 20047 | Quantile regression | 2004 | Quantile 0.5 | 12259.7601 | 5915.5236 | 2.540802e+04 |
| 20057 | Quantile regression | 2005 | Quantile 0.5 | 13807.2602 | 6484.3432 | 2.940011e+04 |
| 20067 | Quantile regression | 2006 | Quantile 0.5 | 15550.0951 | 7106.3229 | 3.402680e+04 |
| 20077 | Quantile regression | 2007 | Quantile 0.5 | 17512.9210 | 7786.4498 | 3.938925e+04 |
| 20087 | Quantile regression | 2008 | Quantile 0.5 | 19723.5065 | 8530.1752 | 4.560477e+04 |
| 20097 | Quantile regression | 2009 | Quantile 0.5 | 22213.1253 | 9343.4574 | 5.280946e+04 |
| 20107 | Quantile regression | 2010 | Quantile 0.5 | 25016.9987 | 10232.8094 | 6.116113e+04 |
| 20117 | Quantile regression | 2011 | Quantile 0.5 | 28174.7937 | 11205.3507 | 7.084285e+04 |
| 20127 | Quantile regression | 2012 | Quantile 0.5 | 31731.1844 | 12268.8636 | 8.206694e+04 |
| 20137 | Quantile regression | 2013 | Quantile 0.5 | 35736.4841 | 13431.8554 | 9.507966e+04 |
| 20147 | Quantile regression | 2014 | Quantile 0.5 | 40247.3565 | 14703.6260 | 1.101667e+05 |
| 20157 | Quantile regression | 2015 | Quantile 0.5 | 45327.6183 | 16094.3414 | 1.276593e+05 |
| 20167 | Quantile regression | 2016 | Quantile 0.5 | 51049.1410 | 17615.1152 | 1.479420e+05 |
| 20177 | Quantile regression | 2017 | Quantile 0.5 | 57492.8686 | 19278.0962 | 1.714604e+05 |
| 19708 | Quantile regression | 1970 | Quantile 0.9 | 695.8591 | 193.9546 | 2.496563e+03 |
| 19718 | Quantile regression | 1971 | Quantile 0.9 | 826.1722 | 241.9242 | 2.821382e+03 |
| 19728 | Quantile regression | 1972 | Quantile 0.9 | 980.8889 | 301.3893 | 3.192360e+03 |
| 19738 | Quantile regression | 1973 | Quantile 0.9 | 1164.5794 | 374.9539 | 3.617099e+03 |
| 19748 | Quantile regression | 1974 | Quantile 0.9 | 1382.6694 | 465.7494 | 4.104728e+03 |
| 19758 | Quantile regression | 1975 | Quantile 0.9 | 1641.6010 | 577.5155 | 4.666288e+03 |
| 19768 | Quantile regression | 1976 | Quantile 0.9 | 1949.0225 | 714.6815 | 5.315219e+03 |
| 19778 | Quantile regression | 1977 | Quantile 0.9 | 2314.0147 | 882.4455 | 6.067983e+03 |
| 19788 | Quantile regression | 1978 | Quantile 0.9 | 2747.3588 | 1086.8414 | 6.944878e+03 |
| 19798 | Quantile regression | 1979 | Quantile 0.9 | 3261.8549 | 1334.7877 | 7.971078e+03 |
| 19808 | Quantile regression | 1980 | Quantile 0.9 | 3872.7004 | 1634.1090 | 9.177973e+03 |
| 19818 | Quantile regression | 1981 | Quantile 0.9 | 4597.9385 | 1993.5217 | 1.060487e+04 |
| 19828 | Quantile regression | 1982 | Quantile 0.9 | 5458.9914 | 2422.5834 | 1.230116e+04 |
| 19838 | Quantile regression | 1983 | Quantile 0.9 | 6481.2932 | 2931.6102 | 1.432904e+04 |
| 19848 | Quantile regression | 1984 | Quantile 0.9 | 7695.0408 | 3531.5799 | 1.676690e+04 |
| 19858 | Quantile regression | 1985 | Quantile 0.9 | 9136.0862 | 4234.0476 | 1.971354e+04 |
| 19868 | Quantile regression | 1986 | Quantile 0.9 | 10846.9952 | 5051.1098 | 2.329336e+04 |
| 19878 | Quantile regression | 1987 | Quantile 0.9 | 12878.3049 | 5995.4488 | 2.766277e+04 |
| 19888 | Quantile regression | 1988 | Quantile 0.9 | 15290.0168 | 7080.4793 | 3.301819e+04 |
| 19898 | Quantile regression | 1989 | Quantile 0.9 | 18153.3684 | 8320.5977 | 3.960590e+04 |
| 19908 | Quantile regression | 1990 | Quantile 0.9 | 21552.9380 | 9731.5163 | 4.773451e+04 |
| 19918 | Quantile regression | 1991 | Quantile 0.9 | 25589.1428 | 11330.6513 | 5.779052e+04 |
| 19928 | Quantile regression | 1992 | Quantile 0.9 | 30381.2051 | 13137.5344 | 7.025806e+04 |
| 19938 | Quantile regression | 1993 | Quantile 0.9 | 36070.6738 | 15174.2239 | 8.574366e+04 |
| 19948 | Quantile regression | 1994 | Quantile 0.9 | 42825.6056 | 17465.7066 | 1.050076e+05 |
| 19958 | Quantile regression | 1995 | Quantile 0.9 | 50845.5292 | 20040.2924 | 1.290035e+05 |
| 19968 | Quantile regression | 1996 | Quantile 0.9 | 60367.3387 | 22930.0114 | 1.589278e+05 |
| 19978 | Quantile regression | 1997 | Quantile 0.9 | 71672.2913 | 26171.0284 | 1.962826e+05 |
| 19988 | Quantile regression | 1998 | Quantile 0.9 | 85094.3151 | 29804.0882 | 2.429547e+05 |
| 19998 | Quantile regression | 1999 | Quantile 0.9 | 101029.8726 | 33875.0066 | 3.013146e+05 |
| 20008 | Quantile regression | 2000 | Quantile 0.9 | 119949.6718 | 38435.2172 | 3.743422e+05 |
| 20018 | Quantile regression | 2001 | Quantile 0.9 | 142412.5695 | 43542.3862 | 4.657838e+05 |
| 20028 | Quantile regression | 2002 | Quantile 0.9 | 169082.0796 | 49261.1036 | 5.803514e+05 |
| 20038 | Quantile regression | 2003 | Quantile 0.9 | 200745.9717 | 55663.6614 | 7.239722e+05 |
| 20048 | Quantile regression | 2004 | Quantile 0.9 | 238339.5405 | 62830.9275 | 9.041047e+05 |
| 20058 | Quantile regression | 2005 | Quantile 0.9 | 282973.2326 | 70853.3276 | 1.130135e+06 |
| 20068 | Quantile regression | 2006 | Quantile 0.9 | 335965.4475 | 79831.9450 | 1.413880e+06 |
| 20078 | Quantile regression | 2007 | Quantile 0.9 | 398881.4803 | 89879.7530 | 1.770214e+06 |
| 20088 | Quantile regression | 2008 | Quantile 0.9 | 473579.7581 | 101122.9932 | 2.217871e+06 |
| 20098 | Quantile regression | 2009 | Quantile 0.9 | 562266.7342 | 113702.7174 | 2.780443e+06 |
| 20108 | Quantile regression | 2010 | Quantile 0.9 | 667562.0631 | 127776.5105 | 3.487645e+06 |
| 20118 | Quantile regression | 2011 | Quantile 0.9 | 792575.9804 | 143520.4157 | 4.376915e+06 |
| 20128 | Quantile regression | 2012 | Quantile 0.9 | 941001.1735 | 161131.0856 | 5.495421e+06 |
| 20138 | Quantile regression | 2013 | Quantile 0.9 | 1117221.8568 | 180828.1830 | 6.902600e+06 |
| 20148 | Quantile regression | 2014 | Quantile 0.9 | 1326443.2739 | 202857.0643 | 8.673357e+06 |
| 20158 | Quantile regression | 2015 | Quantile 0.9 | 1574845.4510 | 227491.7735 | 1.090210e+07 |
| 20168 | Quantile regression | 2016 | Quantile 0.9 | 1869765.7436 | 255038.3863 | 1.370783e+07 |
| 20178 | Quantile regression | 2017 | Quantile 0.9 | 2219915.5694 | 285838.7435 | 1.724058e+07 |
Appendix to figure B
Appendix Taxo.1: single-group estimates versus multi-group estimates
Number of estimates per group:
plyr::count(db.over.time$Kingdom)| x | freq |
|---|---|
| Animalia | 2752 |
| Animalia/Plantae | 11 |
| Chromista | 1 |
| Diverse/Unspecified | 172 |
| Plantae | 670 |
plyr::count(db.over.time$Phylum)| x | freq |
|---|---|
| Arthropoda | 1567 |
| Arthropoda/Chordata/Tracheophyta | 1 |
| Chlorophyta | 7 |
| Chordata | 1028 |
| Ctenophora | 1 |
| Diverse/Unspecified | 276 |
| Mollusca | 104 |
| Nematoda | 39 |
| Ochrophyta | 1 |
| Tracheophyta | 582 |
plyr::count(db.over.time$Class)| x | freq |
|---|---|
| Actinopterygii | 2 |
| Actinopterygii/Magnoliopsida/Malacostraca/Mammalia | 1 |
| Amphibia | 120 |
| Amphibia/Reptilia | 1 |
| Arachnida | 2 |
| Arachnida/Insecta | 3 |
| Aves | 56 |
| Aves/Mammalia | 1 |
| Bacillariophyceae | 1 |
| Bivalvia | 73 |
| Cephalaspidomorphi | 2 |
| Diverse/Unspecified | 278 |
| Gastropoda | 31 |
| Insecta | 1558 |
| Liliopsida | 99 |
| Magnoliopsida | 461 |
| Malacostraca | 3 |
| Mammalia | 817 |
| Pinopsida | 15 |
| Polypodiopsida | 7 |
| Reptilia | 29 |
| Secernentea | 39 |
| Ulvophyceae | 7 |
single.group <- summarizeCosts(db.over.time[which(
db.over.time$Kingdom %in% "Plantae" |
db.over.time$Phylum %in% c("Arthropoda",
"Chordata",
"Ctenophora",
"Mollusca",
"Mollusca/Mollusca",
"Nematoda")), ],
minimum.year = 1970,
maximum.year = 2017)
multi.group <- summarizeCosts(db.over.time[which(
db.over.time$Kingdom %in% c("Animalia/Plantae", "Diverse/Unspecified") |
!(db.over.time$Kingdom %in% "Plantae") &
db.over.time$Phylum %in% c("Diverse/Unspecified")), ],
minimum.year = 1970,
maximum.year = 2017)In total, there are 1319 economic estimates in the database after we applied our filters. Out of these, there were 1223 economic estimates that could be attributed to a single major taxonomic group (Plantae, Arthropoda or Chordata), i.e. 92.72%. However, these single-group estimates accounted in total for US$ billion 591, i.e. 45.85 % of the global cumulated cost of invasive species. Multi-group estimates accounted for US$ billion 698, i.e. 54.15 % of the global cumulated cost.
p1 <- plot(single.group, graphical.parameters =
ylab(paste0("Annual cost in US$ millions"))) +
theme_minimal() + xlab("Year") +
scale_x_continuous(breaks = global.raw$year.breaks) +
scale_y_log10(breaks = 10^(-15:15),
labels = scales::comma) +
annotation_logticks(sides = "l") +
theme(panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14),
panel.grid.minor = element_blank(),
legend.text = element_text(size = 12)) +
ggtitle("a. Single-group estimates")
p2 <- plot(multi.group, graphical.parameters =
ylab(paste0("Annual cost in US$ millions"))) +
theme_minimal() + xlab("Year") +
scale_x_continuous(breaks = global.raw$year.breaks) +
scale_y_log10(breaks = 10^(-15:15),
labels = scales::comma) +
annotation_logticks(sides = "l") +
theme(panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14),
panel.grid.minor = element_blank(),
legend.text = element_text(size = 12)) +
ggtitle("b. Multi-group estimates")
plot_grid(p1, p2, ncol = 1, align = "v")
Appendix Taxo.2: values
Observed
obs.impact## $Plants
## $Plants$Plants
##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 221 (number of individual year values: 670)
## - Cost values in US$ millions:
## o Total cost over the entire period 8,902.48
## o Average annual cost over the entire period 185.47
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 5.75 0.58 3
## 2 1980 1989 10 119.14 11.91 5
## 3 1990 1999 10 3,305.59 330.56 63
## 4 2000 2009 10 2,715.33 271.53 168
## 5 2010 2017 8 2,756.67 344.58 18
## number_year_values
## 1 12
## 2 28
## 3 198
## 4 401
## 5 31
##
## $Plants$Liliopsida
##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 33 (number of individual year values: 99)
## - Cost values in US$ millions:
## o Total cost over the entire period 1,792.75
## o Average annual cost over the entire period 37.35
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 <NA> <NA> 0
## 2 1980 1989 10 106.01 10.60 2
## 3 1990 1999 10 119.31 11.93 6
## 4 2000 2009 10 187.13 18.71 24
## 5 2010 2017 8 1,380.30 172.54 6
## number_year_values
## 1 0
## 2 11
## 3 22
## 4 51
## 5 15
##
## $Plants$Magnoliopsida
##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 142 (number of individual year values: 461)
## - Cost values in US$ millions:
## o Total cost over the entire period 3,694.23
## o Average annual cost over the entire period 76.96
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 5.75 0.58 3
## 2 1980 1989 10 13.13 1.31 3
## 3 1990 1999 10 2,668.80 266.88 48
## 4 2000 2009 10 724.50 72.45 106
## 5 2010 2017 8 282.05 35.26 9
## number_year_values
## 1 12
## 2 17
## 3 149
## 4 270
## 5 13
##
## $Plants$Pinopsida
##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 2 (number of individual year values: 15)
## - Cost values in US$ millions:
## o Total cost over the entire period 27.16
## o Average annual cost over the entire period 0.57
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 <NA> <NA> 0
## 2 1980 1989 10 <NA> <NA> 0
## 3 1990 1999 10 9.03 0.90 1
## 4 2000 2009 10 18.13 1.81 2
## 5 2010 2017 8 <NA> <NA> 0
## number_year_values
## 1 0
## 2 0
## 3 5
## 4 10
## 5 0
##
## $Plants$Ulvophyceae
##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 2 (number of individual year values: 7)
## - Cost values in US$ millions:
## o Total cost over the entire period 22.79
## o Average annual cost over the entire period 0.47
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 <NA> <NA> 0
## 2 1980 1989 10 <NA> <NA> 0
## 3 1990 1999 10 <NA> <NA> 0
## 4 2000 2009 10 22.79 2.28 2
## 5 2010 2017 8 <NA> <NA> 0
## number_year_values
## 1 0
## 2 0
## 3 0
## 4 7
## 5 0
##
##
## $Invertebrates
## $Invertebrates$Invertebrates
##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 469 (number of individual year values: 1711)
## - Cost values in US$ millions:
## o Total cost over the entire period 415,696.90
## o Average annual cost over the entire period 8,660.35
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 4,279.98 428.00 16
## 2 1980 1989 10 11,641.18 1,164.12 43
## 3 1990 1999 10 24,631.30 2,463.13 105
## 4 2000 2009 10 227,826.74 22,782.67 241
## 5 2010 2017 8 147,317.70 18,414.71 164
## number_year_values
## 1 56
## 2 168
## 3 286
## 4 975
## 5 226
##
## $Invertebrates$Insecta
##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 436 (number of individual year values: 1558)
## - Cost values in US$ millions:
## o Total cost over the entire period 378,411.47
## o Average annual cost over the entire period 7,883.57
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 4,279.98 428.00 16
## 2 1980 1989 10 9,767.90 976.79 37
## 3 1990 1999 10 21,869.88 2,186.99 93
## 4 2000 2009 10 196,789.81 19,678.98 217
## 5 2010 2017 8 145,703.90 18,212.99 160
## number_year_values
## 1 56
## 2 141
## 3 220
## 4 920
## 5 221
##
## $Invertebrates$Bivalvia
##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 12 (number of individual year values: 73)
## - Cost values in US$ millions:
## o Total cost over the entire period 3,133.57
## o Average annual cost over the entire period 65.28
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 <NA> <NA> 0
## 2 1980 1989 10 331.55 33.15 3
## 3 1990 1999 10 1,803.90 180.39 4
## 4 2000 2009 10 997.83 99.78 11
## 5 2010 2017 8 0.30 0.04 1
## number_year_values
## 1 0
## 2 5
## 3 35
## 4 32
## 5 1
##
## $Invertebrates$Gastropoda
##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 8 (number of individual year values: 31)
## - Cost values in US$ millions:
## o Total cost over the entire period 3,365.01
## o Average annual cost over the entire period 70.10
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 <NA> <NA> 0
## 2 1980 1989 10 1,386.33 138.63 2
## 3 1990 1999 10 276.80 27.68 4
## 4 2000 2009 10 89.13 8.91 3
## 5 2010 2017 8 1,612.75 201.59 2
## number_year_values
## 1 0
## 2 14
## 3 10
## 4 4
## 5 3
##
## $Invertebrates$Secernentea
##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 3 (number of individual year values: 39)
## - Cost values in US$ millions:
## o Total cost over the entire period 440.82
## o Average annual cost over the entire period 9.18
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 <NA> <NA> 0
## 2 1980 1989 10 155.40 15.54 1
## 3 1990 1999 10 211.56 21.16 2
## 4 2000 2009 10 73.87 7.39 3
## 5 2010 2017 8 <NA> <NA> 0
## number_year_values
## 1 0
## 2 8
## 3 19
## 4 12
## 5 0
##
##
## $Vertebrates
## $Vertebrates$Vertebrates
##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 526 (number of individual year values: 995)
## - Cost values in US$ millions:
## o Total cost over the entire period 165,949.14
## o Average annual cost over the entire period 3,457.27
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 2,500.92 250.09 6
## 2 1980 1989 10 121.07 12.11 39
## 3 1990 1999 10 529.54 52.95 84
## 4 2000 2009 10 143,040.45 14,304.04 314
## 5 2010 2017 8 19,757.16 2,469.65 131
## number_year_values
## 1 16
## 2 69
## 3 233
## 4 495
## 5 182
##
## $Vertebrates$Amphibia
##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 18 (number of individual year values: 120)
## - Cost values in US$ millions:
## o Total cost over the entire period 39.07
## o Average annual cost over the entire period 0.81
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 <NA> <NA> 0
## 2 1980 1989 10 3.48 0.35 4
## 3 1990 1999 10 15.02 1.50 6
## 4 2000 2009 10 19.82 1.98 15
## 5 2010 2017 8 0.74 0.09 3
## number_year_values
## 1 0
## 2 13
## 3 55
## 4 49
## 5 3
##
## $Vertebrates$Reptilia
##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 18 (number of individual year values: 29)
## - Cost values in US$ millions:
## o Total cost over the entire period 15,015.68
## o Average annual cost over the entire period 312.83
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 <NA> <NA> 0
## 2 1980 1989 10 <NA> <NA> 0
## 3 1990 1999 10 13.76 1.38 4
## 4 2000 2009 10 14,999.09 1,499.91 11
## 5 2010 2017 8 2.84 0.35 4
## number_year_values
## 1 0
## 2 0
## 3 11
## 4 13
## 5 5
##
## $Vertebrates$Aves
##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 29 (number of individual year values: 56)
## - Cost values in US$ millions:
## o Total cost over the entire period 6,230.96
## o Average annual cost over the entire period 129.81
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 4.68 0.47 1
## 2 1980 1989 10 9.89 0.99 2
## 3 1990 1999 10 40.74 4.07 3
## 4 2000 2009 10 6,010.98 601.10 19
## 5 2010 2017 8 164.67 20.58 6
## number_year_values
## 1 1
## 2 2
## 3 10
## 4 31
## 5 12
##
## $Vertebrates$Mammalia
##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 455 (number of individual year values: 784)
## - Cost values in US$ millions:
## o Total cost over the entire period 144,618.01
## o Average annual cost over the entire period 3,012.88
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 2,496.25 249.62 5
## 2 1980 1989 10 107.70 10.77 33
## 3 1990 1999 10 452.90 45.29 69
## 4 2000 2009 10 121,992.95 12,199.29 266
## 5 2010 2017 8 19,568.22 2,446.03 117
## number_year_values
## 1 15
## 2 54
## 3 155
## 4 399
## 5 161
##
##
## $Mixed
## $Mixed$Mixed
##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 95 (number of individual year values: 196)
## - Cost values in US$ millions:
## o Total cost over the entire period 697,531.31
## o Average annual cost over the entire period 14,531.90
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 <NA> <NA> 0
## 2 1980 1989 10 5.68 0.57 1
## 3 1990 1999 10 174,283.19 17,428.32 24
## 4 2000 2009 10 459,492.07 45,949.21 59
## 5 2010 2017 8 63,750.37 7,968.80 23
## number_year_values
## 1 0
## 2 3
## 3 64
## 4 98
## 5 31Predicted
pred.impact## $Plants
## $Plants$Plants
##
## Estimation of annual cost values of invasive alien species over time
##
## - Temporal interval of data: [1970, 2017]
## - Temporal interval used for model calibration: [1970, 2014]
## - Cost transformation: log10
## - Values transformed in US$ million: Yes
## - Estimated average annual cost of invasive alien species in 2017:
##
## o Linear regression:
## . Linear: US$ million 490.95
## . Quadratic: US$ million 11.17
## o Robust regression:
## . Linear: US$ million 3,842.31
## . Quadratic: US$ million 733.96
## o Multiple Adapative Regression Splines: US$ million 0.18
## o Generalized Additive Model: US$ million 23.04
## o Quantile regression:
## . Quantile 0.1: US$ million 0.76
## . Quantile 0.5: US$ million 4,225.72
## . Quantile 0.9: US$ million 17,235.94
##
## You can inspect the summary of each fitted model with summary(object)
##
## $Plants$Liliopsida
##
## Estimation of annual cost values of invasive alien species over time
##
## - Temporal interval of data: [1980, 2017]
## - Temporal interval used for model calibration: [1980, 2014]
## - Cost transformation: log10
## - Values transformed in US$ million: Yes
## - Estimated average annual cost of invasive alien species in 2017:
##
## o Linear regression:
## . Linear: US$ million 7.03
## . Quadratic: US$ million 1.99
## o Robust regression:
## . Linear: US$ million 10.19
## . Quadratic: US$ million 20.65
## o Multiple Adapative Regression Splines: US$ million 1.12
## o Generalized Additive Model: US$ million 12.43
## o Quantile regression:
## . Quantile 0.1: US$ million 0.15
## . Quantile 0.5: US$ million 8.97
## . Quantile 0.9: US$ million 151.83
##
## You can inspect the summary of each fitted model with summary(object)
##
## $Plants$Magnoliopsida
##
## Estimation of annual cost values of invasive alien species over time
##
## - Temporal interval of data: [1970, 2014]
## - Temporal interval used for model calibration: [1970, 2014]
## - Cost transformation: log10
## - Values transformed in US$ million: Yes
## - Estimated average annual cost of invasive alien species in 2017:
##
## o Linear regression:
## . Linear: US$ million 49.93
## . Quadratic: US$ million 1.80
## o Robust regression:
## . Linear: US$ million 410.64
## . Quadratic: US$ million 818.27
## o Multiple Adapative Regression Splines: US$ million 0.00
## o Generalized Additive Model: US$ million 5.49
## o Quantile regression:
## . Quantile 0.1: US$ million 0.42
## . Quantile 0.5: US$ million 434.64
## . Quantile 0.9: US$ million 1,692.73
##
## You can inspect the summary of each fitted model with summary(object)
##
## $Plants$Pinopsida
## [1] NA
##
## $Plants$Ulvophyceae
## [1] NA
##
##
## $Invertebrates
## $Invertebrates$Invertebrates
##
## Estimation of annual cost values of invasive alien species over time
##
## - Temporal interval of data: [1970, 2017]
## - Temporal interval used for model calibration: [1970, 2014]
## - Cost transformation: log10
## - Values transformed in US$ million: Yes
## - Estimated average annual cost of invasive alien species in 2017:
##
## o Linear regression:
## . Linear: US$ million 33,359.41
## . Quadratic: US$ million 45,830.14
## o Robust regression:
## . Linear: US$ million 23,757.67
## . Quadratic: US$ million 59,444.79
## o Multiple Adapative Regression Splines: US$ million 11,513.99
## o Generalized Additive Model: US$ million 32,008.33
## o Quantile regression:
## . Quantile 0.1: US$ million 6,490.45
## . Quantile 0.5: US$ million 17,970.46
## . Quantile 0.9: US$ million 132,228.77
##
## You can inspect the summary of each fitted model with summary(object)
##
## $Invertebrates$Insecta
##
## Estimation of annual cost values of invasive alien species over time
##
## - Temporal interval of data: [1970, 2017]
## - Temporal interval used for model calibration: [1970, 2014]
## - Cost transformation: log10
## - Values transformed in US$ million: Yes
## - Estimated average annual cost of invasive alien species in 2017:
##
## o Linear regression:
## . Linear: US$ million 29,991.52
## . Quadratic: US$ million 54,582.91
## o Robust regression:
## . Linear: US$ million 25,240.60
## . Quadratic: US$ million 66,994.30
## o Multiple Adapative Regression Splines: US$ million 42,149.22
## o Generalized Additive Model: US$ million 27,967.40
## o Quantile regression:
## . Quantile 0.1: US$ million 6,875.69
## . Quantile 0.5: US$ million 15,592.20
## . Quantile 0.9: US$ million 131,052.80
##
## You can inspect the summary of each fitted model with summary(object)
##
## $Invertebrates$Bivalvia
##
## Estimation of annual cost values of invasive alien species over time
##
## - Temporal interval of data: [1988, 2011]
## - Temporal interval used for model calibration: [1988, 2011]
## - Cost transformation: log10
## - Values transformed in US$ million: Yes
## - Estimated average annual cost of invasive alien species in 2017:
##
## o Linear regression:
## . Linear: US$ million 10.26
## . Quadratic: US$ million 0.09
## o Robust regression:
## . Linear: US$ million 65.30
## . Quadratic: US$ million 7.40
## o Multiple Adapative Regression Splines: US$ million 0.38
## o Generalized Additive Model: US$ million 15.10
## o Quantile regression:
## . Quantile 0.1: US$ million 2.57
## . Quantile 0.5: US$ million 43.75
## . Quantile 0.9: US$ million 183.46
##
## You can inspect the summary of each fitted model with summary(object)
##
## $Invertebrates$Gastropoda
##
## Estimation of annual cost values of invasive alien species over time
##
## - Temporal interval of data: [1980, 2011]
## - Temporal interval used for model calibration: [1980, 2011]
## - Cost transformation: log10
## - Values transformed in US$ million: Yes
## - Estimated average annual cost of invasive alien species in 2017:
##
## o Linear regression:
## . Linear: US$ million 0.07
## . Quadratic: US$ million 13.66
## o Robust regression:
## . Linear: US$ million 147.40
## . Quadratic: US$ million 0.00
## o Multiple Adapative Regression Splines: US$ million 5.01
## o Generalized Additive Model: US$ million 16.59
## o Quantile regression:
## . Quantile 0.1: US$ million 0.00
## . Quantile 0.5: US$ million 0.00
## . Quantile 0.9: US$ million 2,634.76
##
## You can inspect the summary of each fitted model with summary(object)
##
## $Invertebrates$Secernentea
## [1] NA
##
##
## $Vertebrates
## $Vertebrates$Vertebrates
##
## Estimation of annual cost values of invasive alien species over time
##
## - Temporal interval of data: [1970, 2017]
## - Temporal interval used for model calibration: [1970, 2014]
## - Cost transformation: log10
## - Values transformed in US$ million: Yes
## - Estimated average annual cost of invasive alien species in 2017:
##
## o Linear regression:
## . Linear: US$ million 2,382.20
## . Quadratic: US$ million 17,141.00
## o Robust regression:
## . Linear: US$ million 1,355.28
## . Quadratic: US$ million 13,390.80
## o Multiple Adapative Regression Splines: US$ million 5,794.58
## o Generalized Additive Model: US$ million 2,795.74
## o Quantile regression:
## . Quantile 0.1: US$ million 401.65
## . Quantile 0.5: US$ million 1,009.62
## . Quantile 0.9: US$ million 595,338.42
##
## You can inspect the summary of each fitted model with summary(object)
##
## $Vertebrates$Amphibia
##
## Estimation of annual cost values of invasive alien species over time
##
## - Temporal interval of data: [1986, 2011]
## - Temporal interval used for model calibration: [1986, 2011]
## - Cost transformation: log10
## - Values transformed in US$ million: Yes
## - Estimated average annual cost of invasive alien species in 2017:
##
## o Linear regression:
## . Linear: US$ million 0.65
## . Quadratic: US$ million 0.03
## o Robust regression:
## . Linear: US$ million 1.68
## . Quadratic: US$ million 0.41
## o Multiple Adapative Regression Splines: US$ million 0.07
## o Generalized Additive Model: US$ million 1.80
## o Quantile regression:
## . Quantile 0.1: US$ million 0.70
## . Quantile 0.5: US$ million 1.30
## . Quantile 0.9: US$ million 8.15
##
## You can inspect the summary of each fitted model with summary(object)
##
## $Vertebrates$Reptilia
## [1] NA
##
## $Vertebrates$Aves
##
## Estimation of annual cost values of invasive alien species over time
##
## - Temporal interval of data: [1978, 2015]
## - Temporal interval used for model calibration: [1978, 2014]
## - Cost transformation: log10
## - Values transformed in US$ million: Yes
## - Estimated average annual cost of invasive alien species in 2017:
##
## o Linear regression:
## . Linear: US$ million 3.75
## . Quadratic: US$ million 0.59
## o Robust regression:
## . Linear: US$ million 0.37
## . Quadratic: US$ million 0.07
## o Multiple Adapative Regression Splines: US$ million 6.37
## o Generalized Additive Model: US$ million 5.53
## o Quantile regression:
## . Quantile 0.1: US$ million 0.11
## . Quantile 0.5: US$ million 2.44
## . Quantile 0.9: US$ million 7,914.02
##
## You can inspect the summary of each fitted model with summary(object)
##
## $Vertebrates$Mammalia
##
## Estimation of annual cost values of invasive alien species over time
##
## - Temporal interval of data: [1970, 2017]
## - Temporal interval used for model calibration: [1970, 2014]
## - Cost transformation: log10
## - Values transformed in US$ million: Yes
## - Estimated average annual cost of invasive alien species in 2017:
##
## o Linear regression:
## . Linear: US$ million 1,518.98
## . Quadratic: US$ million 14,629.40
## o Robust regression:
## . Linear: US$ million 931.09
## . Quadratic: US$ million 11,088.94
## o Multiple Adapative Regression Splines: US$ million 9,793.10
## o Generalized Additive Model: US$ million 2,220.76
## o Quantile regression:
## . Quantile 0.1: US$ million 348.37
## . Quantile 0.5: US$ million 190.68
## . Quantile 0.9: US$ million 306,590.62
##
## You can inspect the summary of each fitted model with summary(object)
##
##
## $Mixed
## $Mixed$Mixed
##
## Estimation of annual cost values of invasive alien species over time
##
## - Temporal interval of data: [1987, 2014]
## - Temporal interval used for model calibration: [1987, 2014]
## - Cost transformation: log10
## - Values transformed in US$ million: Yes
## - Estimated average annual cost of invasive alien species in 2017:
##
## o Linear regression:
## . Linear: US$ million 2,551.72
## . Quadratic: US$ million 17.19
## o Robust regression:
## . Linear: US$ million 3,367.12
## . Quadratic: US$ million 20.92
## o Multiple Adapative Regression Splines: US$ million 805.63
## o Generalized Additive Model: US$ million 216.21
## o Quantile regression:
## . Quantile 0.1: US$ million 12.75
## . Quantile 0.5: US$ million 2,527.62
## . Quantile 0.9: US$ million 35,402,556.79
##
## You can inspect the summary of each fitted model with summary(object)Appendix Taxo.3: Linear model summary for each taxon
Plants
"Plants --------------------------------"## [1] "Plants --------------------------------"pred.impact$Plants$Plants$model.summary$robust.linear##
## Call:
## robustbase::lmrob(formula = transf.cost ~ Year, data = yearly.cost.calibration,
## weights = incomplete.year.weights)
## \--> method = "MM"
## Residuals:
## Min 1Q Median 3Q Max
## -3.97827 -0.41335 0.00961 0.28378 1.38496
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.670e+02 1.168e+01 -14.30 <2e-16 ***
## Year 8.459e-02 5.878e-03 14.39 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Robust residual standard error: 0.4949
## Multiple R-squared: 0.8338, Adjusted R-squared: 0.8299
## Convergence in 9 IRWLS iterations
##
## Robustness weights:
## 3 observations c(42,44,45) are outliers with |weight| = 0 ( < 0.0022);
## 4 weights are ~= 1. The remaining 38 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.4138 0.8900 0.9591 0.9170 0.9888 0.9989
## Algorithmic parameters:
## tuning.chi bb tuning.psi refine.tol
## 1.548e+00 5.000e-01 4.685e+00 1.000e-07
## rel.tol scale.tol solve.tol eps.outlier
## 1.000e-07 1.000e-10 1.000e-07 2.222e-03
## eps.x warn.limit.reject warn.limit.meanrw
## 3.663e-09 5.000e-01 5.000e-01
## nResample max.it best.r.s k.fast.s k.max
## 500 50 2 1 200
## maxit.scale trace.lev mts compute.rd fast.s.large.n
## 200 0 1000 0 2000
## psi subsampling cov
## "bisquare" "nonsingular" ".vcov.avar1"
## compute.outlier.stats
## "SM"
## seed : int(0)"Liliopsida --------------------------------"## [1] "Liliopsida --------------------------------"pred.impact$Plants$Liliopsida$model.summary$robust.linear##
## Call:
## robustbase::lmrob(formula = transf.cost ~ Year, data = yearly.cost.calibration,
## weights = incomplete.year.weights)
## \--> method = "MM"
## Residuals:
## Min 1Q Median 3Q Max
## -1.668199 -0.022225 -0.005185 0.052346 2.138822
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.648921 6.712568 -0.395 0.696
## Year 0.001813 0.003377 0.537 0.595
##
## Robust residual standard error: 0.05897
## Multiple R-squared: 0.07383, Adjusted R-squared: 0.04576
## Convergence in 21 IRWLS iterations
##
## Robustness weights:
## 10 observations c(10,14,24,26,30,31,32,33,34,35)
## are outliers with |weight| = 0 ( < 0.0029);
## 4 weights are ~= 1. The remaining 21 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.5822 0.9271 0.9495 0.9172 0.9916 0.9987
## Algorithmic parameters:
## tuning.chi bb tuning.psi refine.tol
## 1.548e+00 5.000e-01 4.685e+00 1.000e-07
## rel.tol scale.tol solve.tol eps.outlier
## 1.000e-07 1.000e-10 1.000e-07 2.857e-03
## eps.x warn.limit.reject warn.limit.meanrw
## 3.663e-09 5.000e-01 5.000e-01
## nResample max.it best.r.s k.fast.s k.max
## 500 50 2 1 200
## maxit.scale trace.lev mts compute.rd fast.s.large.n
## 200 0 1000 0 2000
## psi subsampling cov
## "bisquare" "nonsingular" ".vcov.avar1"
## compute.outlier.stats
## "SM"
## seed : int(0)"Magnoliopsida --------------------------------"## [1] "Magnoliopsida --------------------------------"pred.impact$Plants$Magnoliopsida$model.summary$robust.linear##
## Call:
## robustbase::lmrob(formula = transf.cost ~ Year, data = yearly.cost.calibration,
## weights = incomplete.year.weights)
## \--> method = "MM"
## Residuals:
## Min 1Q Median 3Q Max
## -4.43535 -0.37124 0.02545 0.27098 2.09029
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.425e+02 1.006e+01 -14.16 <2e-16 ***
## Year 7.193e-02 5.054e-03 14.23 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Robust residual standard error: 0.4712
## Multiple R-squared: 0.8083, Adjusted R-squared: 0.8038
## Convergence in 8 IRWLS iterations
##
## Robustness weights:
## 3 observations c(42,44,45) are outliers with |weight| = 0 ( < 0.0022);
## 5 weights are ~= 1. The remaining 37 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.01068 0.90850 0.95480 0.88250 0.98820 0.99750
## Algorithmic parameters:
## tuning.chi bb tuning.psi refine.tol
## 1.548e+00 5.000e-01 4.685e+00 1.000e-07
## rel.tol scale.tol solve.tol eps.outlier
## 1.000e-07 1.000e-10 1.000e-07 2.222e-03
## eps.x warn.limit.reject warn.limit.meanrw
## 3.663e-09 5.000e-01 5.000e-01
## nResample max.it best.r.s k.fast.s k.max
## 500 50 2 1 200
## maxit.scale trace.lev mts compute.rd fast.s.large.n
## 200 0 1000 0 2000
## psi subsampling cov
## "bisquare" "nonsingular" ".vcov.avar1"
## compute.outlier.stats
## "SM"
## seed : int(0)Invertebrates
"Invertebrates --------------------------------"## [1] "Invertebrates --------------------------------"pred.impact$Invertebrates$Invertebrates$model.summary$robust.linear##
## Call:
## robustbase::lmrob(formula = transf.cost ~ Year, data = yearly.cost.calibration,
## weights = incomplete.year.weights)
## \--> method = "MM"
## Residuals:
## Min 1Q Median 3Q Max
## -0.506716 -0.175411 0.001158 0.264475 1.189871
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -82.81547 11.80075 -7.018 1.21e-08 ***
## Year 0.04323 0.00596 7.253 5.55e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Robust residual standard error: 0.2354
## Multiple R-squared: 0.8018, Adjusted R-squared: 0.7972
## Convergence in 20 IRWLS iterations
##
## Robustness weights:
## observation 32 is an outlier with |weight| = 0 ( < 0.0022);
## 6 weights are ~= 1. The remaining 38 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.05022 0.82770 0.93150 0.81940 0.97500 0.99820
## Algorithmic parameters:
## tuning.chi bb tuning.psi refine.tol
## 1.548e+00 5.000e-01 4.685e+00 1.000e-07
## rel.tol scale.tol solve.tol eps.outlier
## 1.000e-07 1.000e-10 1.000e-07 2.222e-03
## eps.x warn.limit.reject warn.limit.meanrw
## 3.663e-09 5.000e-01 5.000e-01
## nResample max.it best.r.s k.fast.s k.max
## 500 50 2 1 200
## maxit.scale trace.lev mts compute.rd fast.s.large.n
## 200 0 1000 0 2000
## psi subsampling cov
## "bisquare" "nonsingular" ".vcov.avar1"
## compute.outlier.stats
## "SM"
## seed : int(0)"Insecta --------------------------------"## [1] "Insecta --------------------------------"pred.impact$Invertebrates$Insecta$model.summary$robust.linear##
## Call:
## robustbase::lmrob(formula = transf.cost ~ Year, data = yearly.cost.calibration,
## weights = incomplete.year.weights)
## \--> method = "MM"
## Residuals:
## Min 1Q Median 3Q Max
## -0.8425001 -0.2201341 -0.0006057 0.2209433 0.9976475
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -85.722033 10.396421 -8.245 2.14e-10 ***
## Year 0.044682 0.005247 8.516 8.92e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Robust residual standard error: 0.2968
## Multiple R-squared: 0.74, Adjusted R-squared: 0.7339
## Convergence in 16 IRWLS iterations
##
## Robustness weights:
## 3 weights are ~= 1. The remaining 42 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.2355 0.8264 0.9436 0.8511 0.9862 0.9980
## Algorithmic parameters:
## tuning.chi bb tuning.psi refine.tol
## 1.548e+00 5.000e-01 4.685e+00 1.000e-07
## rel.tol scale.tol solve.tol eps.outlier
## 1.000e-07 1.000e-10 1.000e-07 2.222e-03
## eps.x warn.limit.reject warn.limit.meanrw
## 3.663e-09 5.000e-01 5.000e-01
## nResample max.it best.r.s k.fast.s k.max
## 500 50 2 1 200
## maxit.scale trace.lev mts compute.rd fast.s.large.n
## 200 0 1000 0 2000
## psi subsampling cov
## "bisquare" "nonsingular" ".vcov.avar1"
## compute.outlier.stats
## "SM"
## seed : int(0)Vertebrates
"Vertebrates --------------------------------"## [1] "Vertebrates --------------------------------"pred.impact$Vertebrates$Vertebrates$model.summary$robust.linear##
## Call:
## robustbase::lmrob(formula = transf.cost ~ Year, data = yearly.cost.calibration,
## weights = incomplete.year.weights)
## \--> method = "MM"
## Residuals:
## Min 1Q Median 3Q Max
## -1.4417 -0.6745 0.0880 0.5764 2.5767
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -108.55548 26.92426 -4.032 0.000222 ***
## Year 0.05537 0.01355 4.088 0.000187 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Robust residual standard error: 0.7975
## Multiple R-squared: 0.3926, Adjusted R-squared: 0.3784
## Convergence in 12 IRWLS iterations
##
## Robustness weights:
## one weight is ~= 1. The remaining 44 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.2750 0.8406 0.9417 0.8745 0.9762 0.9989
## Algorithmic parameters:
## tuning.chi bb tuning.psi refine.tol
## 1.548e+00 5.000e-01 4.685e+00 1.000e-07
## rel.tol scale.tol solve.tol eps.outlier
## 1.000e-07 1.000e-10 1.000e-07 2.222e-03
## eps.x warn.limit.reject warn.limit.meanrw
## 3.663e-09 5.000e-01 5.000e-01
## nResample max.it best.r.s k.fast.s k.max
## 500 50 2 1 200
## maxit.scale trace.lev mts compute.rd fast.s.large.n
## 200 0 1000 0 2000
## psi subsampling cov
## "bisquare" "nonsingular" ".vcov.avar1"
## compute.outlier.stats
## "SM"
## seed : int(0)"Amphibia --------------------------------"## [1] "Amphibia --------------------------------"pred.impact$Vertebrates$Amphibia$model.summary$robust.linear##
## Call:
## robustbase::lmrob(formula = transf.cost ~ Year, data = yearly.cost.calibration,
## weights = incomplete.year.weights)
## \--> method = "MM"
## Residuals:
## Min 1Q Median 3Q Max
## -2.440980 -0.119922 0.001734 0.121201 0.600395
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.818783 12.128027 -0.892 0.381
## Year 0.005475 0.006082 0.900 0.377
##
## Robust residual standard error: 0.157
## Multiple R-squared: 0.06106, Adjusted R-squared: 0.02194
## Convergence in 17 IRWLS iterations
##
## Robustness weights:
## observation 26 is an outlier with |weight| = 0 ( < 0.0038);
## 4 weights are ~= 1. The remaining 21 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1115 0.9280 0.9446 0.8765 0.9949 0.9988
## Algorithmic parameters:
## tuning.chi bb tuning.psi refine.tol
## 1.548e+00 5.000e-01 4.685e+00 1.000e-07
## rel.tol scale.tol solve.tol eps.outlier
## 1.000e-07 1.000e-10 1.000e-07 3.846e-03
## eps.x warn.limit.reject warn.limit.meanrw
## 3.658e-09 5.000e-01 5.000e-01
## nResample max.it best.r.s k.fast.s k.max
## 500 50 2 1 200
## maxit.scale trace.lev mts compute.rd fast.s.large.n
## 200 0 1000 0 2000
## psi subsampling cov
## "bisquare" "nonsingular" ".vcov.avar1"
## compute.outlier.stats
## "SM"
## seed : int(0)"Aves --------------------------------"## [1] "Aves --------------------------------"pred.impact$Vertebrates$Aves$model.summary$robust.linear##
## Call:
## robustbase::lmrob(formula = transf.cost ~ Year, data = yearly.cost.calibration,
## weights = incomplete.year.weights)
## \--> method = "MM"
## Residuals:
## Min 1Q Median 3Q Max
## -0.59541 -0.03606 0.08311 0.61295 3.33148
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 87.38702 26.91996 3.246 0.00343 **
## Year -0.04354 0.01348 -3.229 0.00358 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Robust residual standard error: 0.4908
## Multiple R-squared: 0.5133, Adjusted R-squared: 0.493
## Convergence in 11 IRWLS iterations
##
## Robustness weights:
## 5 observations c(12,13,17,21,23) are outliers with |weight| = 0 ( < 0.0038);
## 5 weights are ~= 1. The remaining 16 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.8434 0.8797 0.9436 0.9345 0.9944 0.9986
## Algorithmic parameters:
## tuning.chi bb tuning.psi refine.tol
## 1.548e+00 5.000e-01 4.685e+00 1.000e-07
## rel.tol scale.tol solve.tol eps.outlier
## 1.000e-07 1.000e-10 1.000e-07 3.846e-03
## eps.x warn.limit.reject warn.limit.meanrw
## 3.663e-09 5.000e-01 5.000e-01
## nResample max.it best.r.s k.fast.s k.max
## 500 50 2 1 200
## maxit.scale trace.lev mts compute.rd fast.s.large.n
## 200 0 1000 0 2000
## psi subsampling cov
## "bisquare" "nonsingular" ".vcov.avar1"
## compute.outlier.stats
## "SM"
## seed : int(0)"Mammalia --------------------------------"## [1] "Mammalia --------------------------------"pred.impact$Vertebrates$Mammalia$model.summary$robust.linear##
## Call:
## robustbase::lmrob(formula = transf.cost ~ Year, data = yearly.cost.calibration,
## weights = incomplete.year.weights)
## \--> method = "MM"
## Residuals:
## Min 1Q Median 3Q Max
## -1.42841 -0.67240 0.07446 0.61802 2.69244
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -103.07196 25.71508 -4.008 0.000239 ***
## Year 0.05257 0.01294 4.064 0.000201 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Robust residual standard error: 0.904
## Multiple R-squared: 0.3299, Adjusted R-squared: 0.3143
## Convergence in 11 IRWLS iterations
##
## Robustness weights:
## 2 weights are ~= 1. The remaining 43 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.3550 0.8830 0.9502 0.8886 0.9696 0.9975
## Algorithmic parameters:
## tuning.chi bb tuning.psi refine.tol
## 1.548e+00 5.000e-01 4.685e+00 1.000e-07
## rel.tol scale.tol solve.tol eps.outlier
## 1.000e-07 1.000e-10 1.000e-07 2.222e-03
## eps.x warn.limit.reject warn.limit.meanrw
## 3.663e-09 5.000e-01 5.000e-01
## nResample max.it best.r.s k.fast.s k.max
## 500 50 2 1 200
## maxit.scale trace.lev mts compute.rd fast.s.large.n
## 200 0 1000 0 2000
## psi subsampling cov
## "bisquare" "nonsingular" ".vcov.avar1"
## compute.outlier.stats
## "SM"
## seed : int(0)Appendix Taxo.3: individual plots for each taxon
Plants
plots <- map(names(obs.impact[["Plants"]]),
function(x, preds)
plot(preds[[x]], graphical.parameters =
ylab(paste0("Annual cost in US$ millions"))) +
theme_minimal() + xlab("Year") +
scale_x_continuous(breaks = global.raw$year.breaks) +
scale_y_log10(breaks = 10^(-15:15),
labels = scales::comma) +
annotation_logticks(sides = "l") +
guides(col = guide_legend(override.aes = list(size = 2)),
shape = guide_legend(override.aes = list(size = 2))) +
theme(panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14),
panel.grid.minor = element_blank(),
legend.text = element_text(size = 14),
legend.title = element_text(size = 18)) + ggtitle(ifelse(x %in% c("Plants", "Invertebrates", "Vertebrates"),
paste0("All ", x), x)),
obs.impact[["Plants"]])
if(any(sapply(plots, is.null)))
{
if(any(sapply(plots, is.null))) { plots <- plots[-which(sapply(plots, is.null))] }
}
plots[[length(plots) + 1]] <- cowplot::get_legend(plots[[1]])
plots[1:(length(plots) - 1)] <- lapply(plots[1:(length(plots) - 1)],
function(x)
{x + theme(legend.position='none')}
)
plot_grid(plotlist = plots, ncol = 2, align = "hv", labels = "auto")
plots <- map(names(pred.impact[["Plants"]]),
function(x, preds)
if(!all(is.na(preds[[x]])))
{
plot(preds[[x]], models = "robust.linear",
graphical.parameters =
ylab(paste0("Annual cost in US$ millions"))) +
theme_minimal() + xlab("Year") +
scale_x_continuous(breaks = global.raw$year.breaks) +
scale_y_log10(breaks = 10^(-15:15),
labels = scales::comma) +
annotation_logticks(sides = "l") +
guides(col = guide_legend(override.aes = list(size = 2)),
shape = guide_legend(override.aes = list(size = 2))) +
theme(panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14),
panel.grid.minor = element_blank(),
legend.text = element_text(size = 14),
legend.title = element_text(size = 18),
strip.text = element_text(size = 14)) + ggtitle(ifelse(x %in% c("Plants", "Invertebrates", "Vertebrates"), paste0("All ", x), x))
},
pred.impact[["Plants"]])
if(any(sapply(plots, is.null)))
{
plots <- plots[-which(sapply(plots, is.null))]
}
plots[[length(plots) + 1]] <- cowplot::get_legend(plots[[1]])
plots[1:(length(plots) - 1)] <- lapply(plots[1:(length(plots) - 1)],
function(x)
{x + theme(legend.position='none')}
)
plot_grid(plotlist = plots, ncol = 2, align = "hv", labels = "auto")
Invertebrates
plots <- map(names(obs.impact[["Invertebrates"]]),
function(x, preds)
plot(preds[[x]], graphical.parameters =
ylab(paste0("Annual cost in US$ millions"))) +
theme_minimal() + xlab("Year") +
scale_x_continuous(breaks = global.raw$year.breaks) +
scale_y_log10(breaks = 10^(-15:15),
labels = scales::comma) +
annotation_logticks(sides = "l") +
guides(col = guide_legend(override.aes = list(size = 2)),
shape = guide_legend(override.aes = list(size = 2))) +
theme(panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14),
panel.grid.minor = element_blank(),
legend.text = element_text(size = 14),
legend.title = element_text(size = 18)) + ggtitle(ifelse(x %in% c("Plants", "Invertebrates", "Vertebrates"), paste0("All ", x), x)),
obs.impact[["Invertebrates"]])
if(any(sapply(plots, is.null))) { plots <- plots[-which(sapply(plots, is.null))] }
plots[[length(plots) + 1]] <- cowplot::get_legend(plots[[1]])
plots[1:(length(plots) - 1)] <- lapply(plots[1:(length(plots) - 1)],
function(x)
{x + theme(legend.position='none')}
)
plot_grid(plotlist = plots, ncol = 2, align = "hv", labels = "auto")
# Removing Gastropods & Bivalves
pred.impact[["Invertebrates"]]$Gastropoda <- NA
pred.impact[["Invertebrates"]]$Bivalvia <- NA
plots <- map(names(pred.impact[["Invertebrates"]]),
function(x, preds)
if(!all(is.na(preds[[x]])))
{
plot(preds[[x]],
models = "robust.linear",
graphical.parameters =
ylab(paste0("Annual cost in US$ millions"))) +
theme_minimal() + xlab("Year") +
scale_x_continuous(breaks = global.raw$year.breaks) +
scale_y_log10(breaks = 10^(-15:15),
labels = scales::comma) +
annotation_logticks(sides = "l") +
guides(col = guide_legend(override.aes = list(size = 2)),
shape = guide_legend(override.aes = list(size = 2))) +
theme(panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14),
panel.grid.minor = element_blank(),
legend.text = element_text(size = 14),
legend.title = element_text(size = 18),
strip.text = element_text(size = 14)) + ggtitle(ifelse(x %in% c("Plants", "Invertebrates", "Vertebrates"), paste0("All ", x), x))
},
pred.impact[["Invertebrates"]])
if(any(sapply(plots, is.null))) { plots <- plots[-which(sapply(plots, is.null))] }
plots[[length(plots) + 1]] <- cowplot::get_legend(plots[[1]])
plots[1:(length(plots) - 1)] <- lapply(plots[1:(length(plots) - 1)],
function(x)
{x + theme(legend.position='none')}
)
plot_grid(plotlist = plots, ncol = 2, align = "hv", labels = "auto")
Vertebrates
plots <- map(names(obs.impact[["Vertebrates"]]),
function(x, preds)
plot(preds[[x]],
graphical.parameters =
ylab(paste0("Annual cost in US$ millions"))) +
theme_minimal() + xlab("Year") +
scale_x_continuous(breaks = global.raw$year.breaks) +
scale_y_log10(breaks = 10^(-15:15),
labels = scales::comma) +
annotation_logticks(sides = "l") +
guides(col = guide_legend(override.aes = list(size = 2)),
shape = guide_legend(override.aes = list(size = 2))) +
theme(panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14),
panel.grid.minor = element_blank(),
legend.text = element_text(size = 14),
legend.title = element_text(size = 18)) + ggtitle(ifelse(x %in% c("Plants", "Invertebrates", "Vertebrates"), paste0("All ", x), x)),
obs.impact[["Vertebrates"]])
if(any(sapply(plots, is.null))) { plots <- plots[-which(sapply(plots, is.null))] }
plots[[length(plots) + 1]] <- cowplot::get_legend(plots[[1]])
plots[1:(length(plots) - 1)] <- lapply(plots[1:(length(plots) - 1)],
function(x)
{x + theme(legend.position='none')}
)
plot_grid(plotlist = plots, ncol = 2, align = "hv", labels = "auto")
plots <- map(names(pred.impact[["Vertebrates"]]),
function(x, preds)
if(!all(is.na(preds[[x]])))
{
plot(preds[[x]],
models = "robust.linear",
graphical.parameters =
ylab(paste0("Annual cost in US$ millions"))) +
theme_minimal() + xlab("Year") +
scale_x_continuous(breaks = global.raw$year.breaks) +
scale_y_log10(breaks = 10^(-15:15),
labels = scales::comma) +
annotation_logticks(sides = "l") +
guides(col = guide_legend(override.aes = list(size = 2)),
shape = guide_legend(override.aes = list(size = 2))) +
theme(panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14),
panel.grid.minor = element_blank(),
legend.text = element_text(size = 14),
legend.title = element_text(size = 18),
strip.text = element_text(size = 14)) + ggtitle(ifelse(x %in% c("Plants", "Invertebrates", "Vertebrates"), paste0("All ", x), x))
},
pred.impact[["Vertebrates"]])
if(any(sapply(plots, is.null))) { plots <- plots[-which(sapply(plots, is.null))] }
plots[[length(plots) + 1]] <- cowplot::get_legend(plots[[1]])
plots[1:(length(plots) - 1)] <- lapply(plots[1:(length(plots) - 1)],
function(x)
{x + theme(legend.position='none')}
)
plot_grid(plotlist = plots, ncol = 2, align = "hv", labels = "auto")
Appendix to figure C
Appendix Region.1: Number and percentage of estimates per region
Number of estimates per region:
plyr::count(db.over.time$Geographic_region)| x | freq |
|---|---|
| Africa | 537 |
| Africa/Asia/Europe | 1 |
| Asia | 420 |
| Central America | 164 |
| Central America/North America | 43 |
| Central America/North America/South America | 22 |
| Central America/Oceania-Pacific islands | 1 |
| Central America/South America | 5 |
| Diverse/Unspecified | 9 |
| Europe | 339 |
| Europe/ Asia | 1 |
| Europe/North America | 13 |
| North America | 947 |
| Oceania | 791 |
| Oceania/Pacific Islands | 96 |
| South America | 217 |
single.region <- summarizeCosts(db.over.time[which(
db.over.time$Geographic_region %in% names(Regions)), ],
minimum.year = 1970,
maximum.year = 2017)
multi.group <- summarizeCosts(db.over.time[which(
!(db.over.time$Geographic_region %in% names(Regions))), ],
minimum.year = 1970,
maximum.year = 2017)There were 1228 economic estimates that could be attributed to a single region, i.e. 93.1% of our filtered database However, these single region estimates accounted in total for US$ billion 959, i.e. 74.44 % of the global cumulated cost of invasive species. Multi-region estimates accounted for US$ billion 329, i.e. 25.56 % of the global cumulated cost.
Appendix Region.2: Observed and estimated costs for each region
The lack of data for multiple regions prevent a reasonable comparison of linear regression estimates for 2017 between regions.
obs.regions## $Africa
##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 130 (number of individual year values: 537)
## - Cost values in US$ millions:
## o Total cost over the entire period 10,190.34
## o Average annual cost over the entire period 212.30
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 6.37 0.64 2
## 2 1980 1989 10 26.44 2.64 4
## 3 1990 1999 10 2,983.57 298.36 36
## 4 2000 2009 10 3,696.01 369.60 76
## 5 2010 2017 8 3,477.94 434.74 46
## number_year_values
## 1 11
## 2 25
## 3 142
## 4 296
## 5 63
##
## $Asia
##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 105 (number of individual year values: 420)
## - Cost values in US$ millions:
## o Total cost over the entire period 269,192.52
## o Average annual cost over the entire period 5,608.18
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 17.55 1.76 1
## 2 1980 1989 10 1,621.17 162.12 8
## 3 1990 1999 10 573.59 57.36 11
## 4 2000 2009 10 208,620.19 20,862.02 68
## 5 2010 2017 8 58,360.01 7,295.00 44
## number_year_values
## 1 1
## 2 24
## 3 40
## 4 297
## 5 58
##
## $Europe
##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 151 (number of individual year values: 339)
## - Cost values in US$ millions:
## o Total cost over the entire period 28,798.74
## o Average annual cost over the entire period 599.97
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 192.49 19.25 3
## 2 1980 1989 10 6.24 0.62 2
## 3 1990 1999 10 240.63 24.06 12
## 4 2000 2009 10 10,640.30 1,064.03 101
## 5 2010 2017 8 17,719.08 2,214.88 52
## number_year_values
## 1 12
## 2 10
## 3 56
## 4 189
## 5 72
##
## $`Central America`
##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 47 (number of individual year values: 164)
## - Cost values in US$ millions:
## o Total cost over the entire period 5,758.67
## o Average annual cost over the entire period 119.97
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 43.63 4.36 3
## 2 1980 1989 10 277.67 27.77 3
## 3 1990 1999 10 127.06 12.71 11
## 4 2000 2009 10 4,668.73 466.87 26
## 5 2010 2017 8 641.58 80.20 10
## number_year_values
## 1 3
## 2 3
## 3 14
## 4 133
## 5 11
##
## $`North America`
##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 310 (number of individual year values: 947)
## - Cost values in US$ millions:
## o Total cost over the entire period 528,738.48
## o Average annual cost over the entire period 11,015.38
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 6,395.07 639.51 13
## 2 1980 1989 10 9,101.97 910.20 30
## 3 1990 1999 10 30,726.16 3,072.62 91
## 4 2000 2009 10 480,928.95 48,092.90 163
## 5 2010 2017 8 1,586.32 198.29 73
## number_year_values
## 1 49
## 2 125
## 3 292
## 4 380
## 5 101
##
## $Oceania
##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 410 (number of individual year values: 758)
## - Cost values in US$ millions:
## o Total cost over the entire period 31,284.89
## o Average annual cost over the entire period 651.77
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 12.50 1.25 2
## 2 1980 1989 10 105.32 10.53 35
## 3 1990 1999 10 918.29 91.83 73
## 4 2000 2009 10 14,189.14 1,418.91 278
## 5 2010 2017 8 16,059.65 2,007.46 55
## number_year_values
## 1 3
## 2 58
## 3 171
## 4 445
## 5 81
##
## $`South America`
##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 68 (number of individual year values: 217)
## - Cost values in US$ millions:
## o Total cost over the entire period 84,940.84
## o Average annual cost over the entire period 1,769.60
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 <NA> <NA> 0
## 2 1980 1989 10 3.27 0.33 1
## 3 1990 1999 10 117.38 11.74 11
## 4 2000 2009 10 79,914.54 7,991.45 32
## 5 2010 2017 8 4,905.66 613.21 29
## number_year_values
## 1 0
## 2 6
## 3 17
## 4 148
## 5 46plots <- map(names(obs.regions),
function(x, preds)
if(!all(is.na(preds[[x]])))
{
plot(preds[[x]], graphical.parameters =
ylab(paste0("Annual cost in US$ millions"))) +
theme_minimal() + xlab("Year") +
scale_x_continuous(breaks = global.raw$year.breaks) +
scale_y_log10(breaks = 10^(-15:15),
labels = scales::comma) +
annotation_logticks(sides = "l") +
guides(col = guide_legend(override.aes = list(size = 2)),
shape = guide_legend(override.aes = list(size = 2))) +
theme(panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14),
panel.grid.minor = element_blank(),
legend.text = element_text(size = 14),
legend.title = element_text(size = 18)) + ggtitle(x)
},
obs.regions)
if(any(sapply(plots, is.null))) { plots <- plots[-which(sapply(plots, is.null))] }
plots[[length(plots) + 1]] <- cowplot::get_legend(plots[[1]])
plots[1:(length(plots) - 1)] <- lapply(plots[1:(length(plots) - 1)],
function(x)
{x + theme(legend.position='none')}
)
plot_grid(plotlist = plots, ncol = 2, align = "hv", labels = "auto")
Appendix Region.3: Maps of cumulated costs for all regions
tm_shape(invacost.world,
projection="+proj=eck4") +
tm_polygons("cumulated.cost",
palette = "OrRd",
title = "Cost in US$ millions",
style="pretty") +
tm_layout(title = "Average annual cost from 2000 to 2009",
title.size = .9,
legend.bg.color = "white",
title.position=c("center", "top"),
inner.margins=c(.04, .01, .1, .01),
bg.color="#AEDFE5",
outer.bg.color="white",
earth.boundary=c(-180, 180, -70, 90),
earth.boundary.color="white",
earth.boundary.lwd=.4,
space.color="white",
attr.outside=T,
attr.color="grey20",
frame = FALSE)
Appendix to figure D
Appendix costtype.1: Number and percentage of estimates per category of cost
Number of estimates per region:
invacost[is.na(invacost$Type_of_cost_merged), ]| Cost_ID | Repository | Reference_ID | Reference_title | Authors | Publication_year | Type_of_material | Previous_materials | Availability | Kingdom | Phylum | Class | Order | Family | Genus | Species | Common_name | Environment | Geographic_region | Official_country | State.Province | Location | Spatial_scale | Period_of_estimation | Time_range | Probable_starting_year | Probable_ending_year | Occurrence | Raw_cost_estimate_local_currency | Min_Raw_cost_estimate_local_currency | Max_Raw_cost_estimate_local_currency | Raw_cost_estimate_2017_USD_exchange_rate | Raw_cost_estimate_2017_USD_PPP | Cost_estimate_per_year_local_currency | Cost_estimate_per_year_2017_USD_exchange_rate | Cost_estimate_per_year_2017_USD_PPP | Currency | Applicable_year | Type_of_applicable_year | Implementation | Acquisition_method | Impacted_sector | Type_of_cost | Method_reliability | Details | Benefit_value.s. | Contributors | Island | verbatimHabitat | Habitat | protectedArea | Abstract | Language | Probable_starting_year_adjusted | Probable_ending_year_adjusted | Type_of_cost_merged | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 915 | 915 | TC | 561 | The Economic cost of invasive and non-native species in Ireland and Northern Ireland | Kelly et al | 2013 | Official report | Marbuah, G., Gren, I. M., & McKie, B. (2014). Economics of Harmful Invasive Species: A Review. Diversity, 6(3), 500–523. https://doi.org/10.3390/d6030500 | Yes | Animalia/Plantae | Diverse/Unspecified | Diverse/Unspecified | Diverse/Unspecified | Diverse/Unspecified | Diverse/Unspecified | Diverse/Unspecified | Diverse (animals/plants) | Aquatic/Terrestrial | Europe | Ireland | NA | NA | Country | 2012 | Year | 2012 | NA | Potentially ongoing | 202894406 | NA | NA | 278319891.4 | 204829893.2 | 202894406 | 2.783199e+08 | 2.048299e+08 | EUR | 2012 | Effective | Observed | Report/Estimation | NA | NA | High | /! there are a lot of details in this paper, here I only report the general sum (including plants) | no | C.D, C.A., L.N. | NA | NA | NA | NA | NA | EN | 2012 | 2012 | NA |
| 916 | 916 | TC | 561 | The Economic cost of invasive and non-native species in Ireland and Northern Ireland | Kelly et al | 2013 | Official report | Marbuah, G., Gren, I. M., & McKie, B. (2014). Economics of Harmful Invasive Species: A Review. Diversity, 6(3), 500–523. https://doi.org/10.3390/d6030500 | Yes | Animalia/Plantae | Diverse/Unspecified | Diverse/Unspecified | Diverse/Unspecified | Diverse/Unspecified | Diverse/Unspecified | Diverse/Unspecified | Diverse (animals/plants) | Aquatic/Terrestrial | Europe | Northern Ireland | NA | NA | Country | 2012 | Year | 2012 | NA | Potentially ongoing | 58623034 | NA | NA | 80415999.51 | 56469149.35 | 58623034 | 8.041600e+07 | 5.646915e+07 | EUR | 2012 | Effective | Observed | Report/Estimation | NA | NA | High | /! there are a lot of details in this paper, here I only report the general sum (including plants) | no | C.D, C.A., L.N. | NA | NA | NA | NA | NA | EN | 2012 | 2012 | NA |
| 927 | 927 | TC | 381 | Invasive Species : The Search for Solutions | Dybas | 2004 | Peer-reviewed article | none | Yes | Diverse/Unspecified | Diverse/Unspecified | Diverse/Unspecified | Diverse/Unspecified | Diverse/Unspecified | Diverse/Unspecified | Diverse/Unspecified | Invasive species | Aquatic/Terrestrial | North America | USA | NA | NA | Country | 2004 | Year | 2004 | NA | Potentially ongoing | 1.37E+11 | NA | NA | 1.77789E+11 | 1.77789E+11 | 137000000000 | 1.777890e+11 | 1.777890e+11 | USD | 2004 | Publication year | Observed | Report/Estimation | NA | NA | High | /!reference !! from Pimentel ? | no | C.D, C.A., L.N. | NA | NA | NA | NA | NA | EN | 2004 | 2004 | NA |
| 1076 | 1076 | Go | 90 | The Actual and Potential Economic Impact of Invasive Species on the Adirondack Park : A Preliminary Assessment | Yellow Wood Associates Inc. | 2014 | Official report | none | Yes | Diverse/Unspecified | Diverse/Unspecified | Diverse/Unspecified | Diverse/Unspecified | Diverse/Unspecified | Diverse/Unspecified | Diverse/Unspecified | NA | Aquatic/Terrestrial | North America | USA | NA | New York (Adirondacks) | Site | 2013 | Year | 2013 | NA | Potentially ongoing | 708000 | NA | NA | 744964.1046 | 744964.1046 | 708000 | 7.449641e+05 | 7.449641e+05 | USD | 2013 | Effective | Observed | Report/Estimation | NA | NA | High | “In addition to these monetary investments, survey respondents reported coordinating more than 12,000 volunteer hours in 2013 at an approximate value of $708,000”. Defined based on accounts from individual respondents + The value of volunteer hours is calculated by multiplying the number of volunteer hours reported by survey respondents by $28.73 — the value of one hour of volunteer time in New York State. Source: “Economic Impact: 36 Adirondack Nonprofits.” 2013. https://www.generousact.org/leading/economic-impact-study Accessed 06/09/14 | no | C.D, C.A., L.N. | NA | NA | NA | NA | NA | EN | 2013 | 2013 | NA |
| 1241 | 1243 | TC | 54 | Technical support to EU strategy on invasive alien species (IAS) - Assessment of the impacts of IAS in Europe and the EU | Kettunen et al | 2009 | Official report | none | Yes | Animalia | Arthropoda | Insecta | Diptera | Culicidae | Aedes | Aedes albopictus | Asian tiger mosquito | Terrestrial | Europe | Diverse/Unspecified | NA | NA | Continental | NA | Year | NA | NA | Potentially ongoing | 1340000 | NA | NA | 2126878.112 | NA | 1340000 | 2.126878e+06 | NA | EUR | 2009 | Publication year | Observed | Report/Estimation | NA | NA | High | [Unknown] The reference is replaced by references on the Asian longhorned beetle in Kettunen et al. 2009 | no | C.D, C.A., L.N. | NA | NA | NA | NA | NA | EN | 2009 | 2009 | NA |
| 1269 | 1271 | TC | 339 | Future Challenges and Opportunities for Agricultural R&D in the Semi-Arid Tropics | Ryan and Spencer | 2001 | Official report | none | Yes | Animalia | Arthropoda | Insecta | Lepidoptera | Noctuidae | Helicoverpa | Helicoverpa armigera | Cotton bollworm | Terrestrial | Diverse/Unspecified | Diverse/Unspecified | NA | NA | Global | 1992 | Year | 1992 | NA | Potentially ongoing | 317000000 | NA | NA | 438834193.8 | NA | 317000000 | 4.388342e+08 | NA | USD | 2001 | Publication year | Observed | Report/Estimation | Agriculture | NA | High | Cited refs (Shanower et al. 1999, Abate et al. 2000) not retrieved. ICRISAT (International Crops Research Institute for the Semi-Arid Tropics) contacted, no answer | no | C.D, C.A., L.N. | NA | NA | NA | NA | NA | EN | 1992 | 1992 | NA |
| 1620 | 1622 | Go | 46 | Economic Impacts of Invasive Species in the Pacific Northwest Economic Region | PNWER Invasive Species Working Group | 2012 | Official report | none | Yes | Animalia | Arthropoda | Insecta | Coleoptera | Curculionidae | Anthonomus | Anthonomus grandis | Boll weevil | Terrestrial | North America | USA | NA | NA | Country | Since 1890s (Until 2012 = publication year) | Period | 1890 | NA | Potentially ongoing | 406504065 | NA | NA | 433992801.3 | 433992801.3 | 3304911 | 3.528397e+06 | 3.528397e+06 | USD | 2012 | Publication year | Observed | Report/Estimation | Agriculture | NA | High | Costs to U.S. cotton industry since 1890s. No specifications on the type of cost. No reference | no | C.D, C.A., L.N. | NA | NA | NA | NA | NA | EN | 1890 | 2012 | NA |
| 1622 | 1624 | Go | 46 | Economic Impacts of Invasive Species in the Pacific Northwest Economic Region | PNWER Invasive Species Working Group | 2012 | Official report | none | Yes | Animalia | Arthropoda | Insecta | Lepidoptera | Erebidae | Lymantria | Lymantria dispar | Gypsy moths | Terrestrial | North America | USA | NA | Eastern U.S. forests | Site | 1981 | Year | 1981 | NA | Potentially ongoing | 764000000 | NA | NA | 815663430.4 | 815663430.4 | 764000000 | 8.156634e+08 | 8.156634e+08 | USD | 2012 | Publication year | Observed | Report/Estimation | Forestry | NA | High | Costs to Eastern U.S. forests in 1981. No specifications on the type of cost. No reference | no | C.D, C.A., L.N. | NA | NA | NA | NA | NA | EN | 1981 | 1981 | NA |
plyr::count(db.over.time$Type_of_cost_merged)| x | freq |
|---|---|
| Damage_costs | 1109 |
| Management_costs | 2332 |
| Mixed_costs | 115 |
| NA | 50 |
single.type <- summarizeCosts(db.over.time[which(
db.over.time$Type_of_cost_merged %in% paste0(names(costtype), "_costs")), ],
minimum.year = 1970,
maximum.year = 2017)
multi.type <- summarizeCosts(db.over.time[which(
!(db.over.time$Type_of_cost_merged %in% paste0(names(costtype), "_costs"))), ],
minimum.year = 1970,
maximum.year = 2017)There were 1282 economic estimates that could be attributed to a single cost type, i.e. 97.19% of our filtered database However, these single category estimates accounted in total for US$ billion 958, i.e. 74.39 % of the global cumulated cost of invasive species. Estimates based on both damage and management costs accounted for US$ billion 330, i.e. 25.61 % of the global cumulated cost.
Appendix costtype.2: Observed and estimated costs for each cost type
p1 <- plot(obs.costtype$Damage, graphical.parameters =
ylab(paste0("Annual cost in US$ millions"))) +
theme_minimal() + xlab("Year") +
scale_x_continuous(breaks = global.raw$year.breaks) +
scale_y_log10(breaks = 10^(-15:15),
labels = scales::comma) +
annotation_logticks(sides = "l") +
guides(col = guide_legend(override.aes = list(size = 2)),
shape = guide_legend(override.aes = list(size = 2))) +
theme(panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14),
panel.grid.minor = element_blank(),
legend.text = element_text(size = 14),
legend.title = element_text(size = 18)) + ggtitle("Damage")
p2 <- plot(obs.costtype$Management, graphical.parameters =
ylab(paste0("Annual cost in US$ millions"))) +
theme_minimal() + xlab("Year") +
scale_x_continuous(breaks = global.raw$year.breaks) +
scale_y_log10(breaks = 10^(-15:15),
labels = scales::comma) +
annotation_logticks(sides = "l") +
guides(col = guide_legend(override.aes = list(size = 2)),
shape = guide_legend(override.aes = list(size = 2))) +
theme(panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14),
panel.grid.minor = element_blank(),
legend.text = element_text(size = 14),
legend.title = element_text(size = 18)) + ggtitle("Management")
plot_grid(p1, p2, ncol = 1, align = "hv", labels = "auto")
p1 <- plot(pred.costtype.25$Damage, graphical.parameters =
ylab(paste0("Annual cost in US$ millions"))) +
theme_minimal() + xlab("Year") +
scale_x_continuous(breaks = global.raw$year.breaks) +
scale_y_log10(breaks = 10^(-15:15),
labels = scales::comma) +
annotation_logticks(sides = "l") +
guides(col = guide_legend(override.aes = list(size = 2)),
shape = guide_legend(override.aes = list(size = 2))) +
theme(panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14),
panel.grid.minor = element_blank(),
legend.text = element_text(size = 14),
legend.title = element_text(size = 18)) + ggtitle("Damage")
p2 <- plot(pred.costtype.25$Management, graphical.parameters =
ylab(paste0("Annual cost in US$ millions"))) +
theme_minimal() + xlab("Year") +
scale_x_continuous(breaks = global.raw$year.breaks) +
scale_y_log10(breaks = 10^(-15:15),
labels = scales::comma) +
annotation_logticks(sides = "l") +
guides(col = guide_legend(override.aes = list(size = 2)),
shape = guide_legend(override.aes = list(size = 2))) +
theme(panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14),
panel.grid.minor = element_blank(),
legend.text = element_text(size = 14),
legend.title = element_text(size = 18)) + ggtitle("Management")
pleg <- cowplot::get_legend(p1)
p1 <- p1 + theme(legend.position='none')
p2 <- p2 + theme(legend.position='none')
plot_grid(p1, p2, pleg, ncol = 2, align = "hv", labels = "auto")
Model summaries
"________________________________ Damage ________________________________"## [1] "________________________________ Damage ________________________________"pred.costtype.25$Damage$model.summary##
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Summary of model fits ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
##
##
## ______________________________ Ordinary Least Square regression models _______________________________
##
##
## >>>>>>>> Linear regression
##
## R squared: 0.6756504 - Adjusted R squared: 0.6756504 Estimate Std. Error t value Pr(>|t|)
## (Intercept) -149.83731071 22.10040128 -6.779846 2.683344e-08
## Year 0.07673108 0.01113779 6.889255 1.861863e-08
## attr(,"class")
## [1] "coeftest"
## attr(,"method")
## [1] "t test of coefficients"
## attr(,"df")
## [1] 43
## attr(,"nobs")
## [1] 45
## attr(,"logLik")
## 'log Lik.' -47.18365 (df=3)
## ------------------------------------------------------------------------------------------------------------
##
##
##
## >>>>>>>> Quadratic regression
##
## R squared: 0.712572 - Adjusted R squared: 0.712572 Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.281433e+03 2.736211e+03 -2.295668 0.02675462
## Year 6.233213e+00 2.746192e+00 2.269766 0.02841689
## I(Year^2) -1.545302e-03 6.890137e-04 -2.242774 0.03024662
## attr(,"class")
## [1] "coeftest"
## attr(,"method")
## [1] "t test of coefficients"
## attr(,"df")
## [1] 42
## attr(,"nobs")
## [1] 45
## attr(,"logLik")
## 'log Lik.' -44.46453 (df=4)
## ------------------------------------------------------------------------------------------------------------
##
## ______________________________ Robust regression models _______________________________
##
##
## >>>>>>>> Linear regression
##
##
## Call:
## robustbase::lmrob(formula = transf.cost ~ Year, data = yearly.cost.calibration,
## weights = incomplete.year.weights)
## \--> method = "MM"
## Residuals:
## Min 1Q Median 3Q Max
## -1.71297 -0.23587 0.01655 0.35356 1.77031
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.513e+02 1.548e+01 -9.776 1.71e-12 ***
## Year 7.742e-02 7.813e-03 9.908 1.14e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Robust residual standard error: 0.4449
## Multiple R-squared: 0.799, Adjusted R-squared: 0.7943
## Convergence in 19 IRWLS iterations
##
## Robustness weights:
## 7 weights are ~= 1. The remaining 38 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.07769 0.74800 0.92020 0.79390 0.98840 0.99840
## Algorithmic parameters:
## tuning.chi bb tuning.psi refine.tol
## 1.548e+00 5.000e-01 4.685e+00 1.000e-07
## rel.tol scale.tol solve.tol eps.outlier
## 1.000e-07 1.000e-10 1.000e-07 2.222e-03
## eps.x warn.limit.reject warn.limit.meanrw
## 3.663e-09 5.000e-01 5.000e-01
## nResample max.it best.r.s k.fast.s k.max
## 500 50 2 1 200
## maxit.scale trace.lev mts compute.rd fast.s.large.n
## 200 0 1000 0 2000
## psi subsampling cov
## "bisquare" "nonsingular" ".vcov.avar1"
## compute.outlier.stats
## "SM"
## seed : int(0)
## ------------------------------------------------------------------------------------------------------------
##
##
##
## >>>>>>>> Quadratic regression
##
##
## Call:
## robustbase::lmrob(formula = transf.cost ~ Year + I(Year^2), data = yearly.cost.calibration,
## weights = incomplete.year.weights, cov = ".vcov.w")
## \--> method = "MM"
## Residuals:
## Min 1Q Median 3Q Max
## -1.40449 -0.24684 -0.03941 0.40101 1.67601
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.561e+03 2.452e+03 -1.452 0.154
## Year 3.503e+00 2.463e+00 1.422 0.162
## I(Year^2) -8.602e-04 6.182e-04 -1.392 0.171
##
## Robust residual standard error: 0.4529
## Multiple R-squared: 0.7845, Adjusted R-squared: 0.7742
## Convergence in 20 IRWLS iterations
##
## Robustness weights:
## 5 weights are ~= 1. The remaining 40 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1414 0.7543 0.9398 0.8194 0.9852 0.9988
## Algorithmic parameters:
## tuning.chi bb tuning.psi refine.tol
## 1.548e+00 5.000e-01 4.685e+00 1.000e-07
## rel.tol scale.tol solve.tol eps.outlier
## 1.000e-07 1.000e-10 1.000e-07 2.222e-03
## eps.x warn.limit.reject warn.limit.meanrw
## 7.378e-06 5.000e-01 5.000e-01
## nResample max.it best.r.s k.fast.s k.max
## 500 50 2 1 200
## maxit.scale trace.lev mts compute.rd fast.s.large.n
## 200 0 1000 0 2000
## psi subsampling cov
## "bisquare" "nonsingular" ".vcov.w"
## compute.outlier.stats
## "SM"
## seed : int(0)
## ------------------------------------------------------------------------------------------------------------
##
## ______________________________ Generalized Additive Models _______________________________
##
##
##
## Family: gaulss
## Link function: identity logb
##
## Formula:
## transf.cost ~ s(Year, k = gam.k)
## <environment: 0x0000000035052510>
## ~s(Year, k = gam.k)
## <environment: 0x0000000035052510>
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.85560 0.07429 38.437 < 2e-16 ***
## (Intercept).1 -1.66065 0.31477 -5.276 1.32e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(Year) 4.628 5.199 2815.94 <2e-16 ***
## s.1(Year) 7.290 8.058 59.07 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Deviance explained = 100%
## -REML = 35.891 Scale est. = 1 n = 45
## ------------------------------------------------------------------------------------------------------------
##
## ______________________________ Multiple Adaptive Regression Splines _______________________________
##
##
## Call: earth(formula=transf.cost~Year, data=yearly.cost.calibration,
## weights=incomplete.year.weights, pmethod="backward",
## nprune=mars.nprune, nfold=5, ncross=30, varmod.method="lm")
##
## coefficients
## (Intercept) 4.22250745
## h(2004-Year) -0.09162633
##
## Selected 2 of 7 terms, and 1 of 1 predictors
## Termination condition: RSq changed by less than 0.001 at 7 terms
## Importance: Year
## Number of terms at each degree of interaction: 1 1 (additive model)
## GCV 0.4412562 RSS 17.29724 GRSq 0.712974 RSq 0.7384742 CVRSq 0.5901739
##
## Note: the cross-validation sd's below are standard deviations across folds
##
## Cross validation: nterms 2.34 sd 0.72 nvars 1.00 sd 0.00
##
## CVRSq sd MaxErr sd
## 0.59 0.322 -1.73 1.26
##
## varmod: method "lm" min.sd 0.0643 iter.rsq 0.255
##
## stddev of predictions:
## coefficients iter.stderr iter.stderr%
## (Intercept) -0.04730925 0.133685 283
## transf.cost 0.22918869 0.0597492 26
##
## mean smallest largest ratio
## 95% prediction interval 2.519646 0.809274 3.608065 4.458397
##
## 68% 80% 90% 95%
## response values in prediction interval 73 82 93 96
## ------------------------------------------------------------------------------------------------------------
##
## ______________________________ Quantile regressions _______________________________
##
##
## >>>>>>>> 0.1 quantile
##
## $call
## quantreg::rq(formula = transf.cost ~ Year, tau = 0.1, data = yearly.cost.calibration,
## weights = incomplete.year.weights)
##
## $terms
## transf.cost ~ Year
## attr(,"variables")
## list(transf.cost, Year)
## attr(,"factors")
## Year
## transf.cost 0
## Year 1
## attr(,"term.labels")
## [1] "Year"
## attr(,"order")
## [1] 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: 0x0000000035052510>
## attr(,"predvars")
## list(transf.cost, Year)
## attr(,"dataClasses")
## transf.cost Year
## "numeric" "numeric"
##
## $coefficients
## coefficients lower bd upper bd
## (Intercept) -115.5050298 -134.91245755 -89.89413149
## Year 0.0590978 0.04626538 0.06892761
##
## $residuals
## 1 2 3 4 5
## 3.545868e-01 2.954890e-01 2.591971e-01 1.772934e-01 1.181956e-01
## 6 7 8 9 10
## 5.909780e-02 5.329071e-15 4.278275e-01 -2.132964e-02 1.352480e+00
## 11 12 13 14 15
## 6.330803e-01 1.893290e+00 5.718542e-01 5.279675e-01 4.536586e-01
## 16 17 18 19 20
## 3.945608e-01 5.792805e-01 4.993907e-01 6.400636e-01 6.564130e-01
## 21 22 23 24 25
## 7.459335e-01 5.315078e-01 1.037385e+00 4.590570e-01 1.687707e+00
## 26 27 28 29 30
## 1.319079e+00 -3.865689e-02 -5.773160e-15 8.617145e-01 2.582371e+00
## 31 32 33 34 35
## 2.479201e+00 2.382560e+00 1.353244e+00 9.018915e-01 2.349320e+00
## 36 37 38 39 40
## 1.651996e+00 7.204197e-01 1.024658e+00 2.385391e-01 3.897013e-01
## 41 42 43 44 45
## 4.133517e-01 1.543682e+00 1.218814e+00 5.946491e-01 -6.261149e-01
##
## $rdf
## [1] 43
##
## $tau
## [1] 0.1
##
## attr(,"class")
## [1] "summary.rq"
## ------------------------------------------------------------------------------------------------------------
##
## >>>>>>>> 0.5 quantile
##
## $call
## quantreg::rq(formula = transf.cost ~ Year, tau = 0.5, data = yearly.cost.calibration,
## weights = incomplete.year.weights)
##
## $terms
## transf.cost ~ Year
## attr(,"variables")
## list(transf.cost, Year)
## attr(,"factors")
## Year
## transf.cost 0
## Year 1
## attr(,"term.labels")
## [1] "Year"
## attr(,"order")
## [1] 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: 0x0000000035052510>
## attr(,"predvars")
## list(transf.cost, Year)
## attr(,"dataClasses")
## transf.cost Year
## "numeric" "numeric"
##
## $coefficients
## coefficients lower bd upper bd
## (Intercept) -156.01607921 -174.67543767 -138.2589696
## Year 0.07979322 0.07082382 0.0892516
##
## $residuals
## 1 2 3 4 5
## 9.565977e-02 1.586655e-02 -4.112083e-02 -1.437199e-01 -2.235131e-01
## 6 7 8 9 10
## -3.033063e-01 -3.830996e-01 2.403250e-02 -4.458200e-01 9.072942e-01
## 11 12 13 14 15
## 1.671990e-01 1.406713e+00 6.458215e-02 3.108624e-15 -9.500429e-02
## 16 17 18 19 20
## -1.747975e-01 -1.077327e-02 -1.113585e-01 8.619006e-03 4.273027e-03
## 21 22 23 24 25
## 7.309809e-02 -1.620231e-01 3.231586e-01 -2.758647e-01 9.320900e-01
## 26 27 28 29 30
## 5.427660e-01 -8.356648e-01 -8.177034e-01 2.331572e-02 1.723277e+00
## 31 32 33 34 35
## 1.599412e+00 1.482075e+00 4.320631e-01 -3.998439e-02 1.386749e+00
## 36 37 38 39 40
## 6.687290e-01 -2.835425e-01 -1.065814e-14 -8.068139e-01 -6.763471e-01
## 41 42 43 44 45
## -6.733921e-01 4.362431e-01 9.067902e-02 -5.541810e-01 -1.795640e+00
##
## $rdf
## [1] 43
##
## $tau
## [1] 0.5
##
## attr(,"class")
## [1] "summary.rq"
## ------------------------------------------------------------------------------------------------------------
##
## >>>>>>>> 0.9 quantile
##
## $call
## quantreg::rq(formula = transf.cost ~ Year, tau = 0.9, data = yearly.cost.calibration,
## weights = incomplete.year.weights)
##
## $terms
## transf.cost ~ Year
## attr(,"variables")
## list(transf.cost, Year)
## attr(,"factors")
## Year
## transf.cost 0
## Year 1
## attr(,"term.labels")
## [1] "Year"
## attr(,"order")
## [1] 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: 0x0000000035052510>
## attr(,"predvars")
## list(transf.cost, Year)
## attr(,"dataClasses")
## transf.cost Year
## "numeric" "numeric"
##
## $coefficients
## coefficients lower bd upper bd
## (Intercept) -206.8129353 -260.0281937 -84.0504546
## Year 0.1059196 0.0452133 0.1326794
##
## $residuals
## 1 2 3 4 5
## -5.764968e-01 -6.824164e-01 -7.655302e-01 -8.942557e-01 -1.000175e+00
## 6 7 8 9 10
## -1.106095e+00 -1.212015e+00 -8.310089e-01 -1.326988e+00 4.884981e-14
## 11 12 13 14 15
## -7.662216e-01 4.471664e-01 -9.210912e-01 -1.011800e+00 -1.132930e+00
## 16 17 18 19 20
## -1.238850e+00 -1.100952e+00 -1.227664e+00 -1.133813e+00 -1.164285e+00
## 21 22 23 24 25
## -1.121587e+00 -1.382834e+00 -9.237789e-01 -1.548929e+00 -3.671002e-01
## 26 27 28 29 30
## -7.825506e-01 -2.187108e+00 -2.195273e+00 -1.380380e+00 2.934547e-01
## 31 32 33 34 35
## 1.434631e-01 4.973799e-14 -1.076138e+00 -1.574312e+00 -1.737057e-01
## 36 37 38 39 40
## -9.178516e-01 -1.896250e+00 -1.638833e+00 -2.471774e+00 -2.367433e+00
## 41 42 43 44 45
## -2.390605e+00 -1.307096e+00 -1.678786e+00 -2.349773e+00 -3.617359e+00
##
## $rdf
## [1] 43
##
## $tau
## [1] 0.9
##
## attr(,"class")
## [1] "summary.rq"
## ------------------------------------------------------------------------------------------------------------"______________________________ Management ______________________________"## [1] "______________________________ Management ______________________________"pred.costtype.25$Management$model.summary##
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Summary of model fits ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
##
##
## ______________________________ Ordinary Least Square regression models _______________________________
##
##
## >>>>>>>> Linear regression
##
## R squared: 0.3069115 - Adjusted R squared: 0.3069115 Estimate Std. Error t value Pr(>|t|)
## (Intercept) -41.86991092 12.38904116 -3.379593 0.0015532259
## Year 0.02239534 0.00623703 3.590705 0.0008408074
## attr(,"class")
## [1] "coeftest"
## attr(,"method")
## [1] "t test of coefficients"
## attr(,"df")
## [1] 43
## attr(,"nobs")
## [1] 45
## attr(,"logLik")
## 'log Lik.' -26.60838 (df=3)
## ------------------------------------------------------------------------------------------------------------
##
##
##
## >>>>>>>> Quadratic regression
##
## R squared: 0.3276923 - Adjusted R squared: 0.3276923 Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.033943e+03 2.263221e+03 -0.8986936 0.3739394
## Year 2.022553e+00 2.275155e+00 0.8889740 0.3790805
## I(Year^2) -5.020477e-04 5.717734e-04 -0.8780536 0.3849102
## attr(,"class")
## [1] "coeftest"
## attr(,"method")
## [1] "t test of coefficients"
## attr(,"df")
## [1] 42
## attr(,"nobs")
## [1] 45
## attr(,"logLik")
## 'log Lik.' -25.92345 (df=4)
## ------------------------------------------------------------------------------------------------------------
##
## ______________________________ Robust regression models _______________________________
##
##
## >>>>>>>> Linear regression
##
##
## Call:
## robustbase::lmrob(formula = transf.cost ~ Year, data = yearly.cost.calibration,
## weights = incomplete.year.weights)
## \--> method = "MM"
## Residuals:
## Min 1Q Median 3Q Max
## -1.22183 -0.17972 -0.09111 0.28722 1.17495
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -45.465164 15.142476 -3.002 0.00445 **
## Year 0.024181 0.007625 3.171 0.00280 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Robust residual standard error: 0.3083
## Multiple R-squared: 0.4335, Adjusted R-squared: 0.4203
## Convergence in 23 IRWLS iterations
##
## Robustness weights:
## 3 weights are ~= 1. The remaining 42 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0810 0.8575 0.9582 0.8410 0.9803 0.9954
## Algorithmic parameters:
## tuning.chi bb tuning.psi refine.tol
## 1.548e+00 5.000e-01 4.685e+00 1.000e-07
## rel.tol scale.tol solve.tol eps.outlier
## 1.000e-07 1.000e-10 1.000e-07 2.222e-03
## eps.x warn.limit.reject warn.limit.meanrw
## 3.663e-09 5.000e-01 5.000e-01
## nResample max.it best.r.s k.fast.s k.max
## 500 50 2 1 200
## maxit.scale trace.lev mts compute.rd fast.s.large.n
## 200 0 1000 0 2000
## psi subsampling cov
## "bisquare" "nonsingular" ".vcov.avar1"
## compute.outlier.stats
## "SM"
## seed : int(0)
## ------------------------------------------------------------------------------------------------------------
##
##
##
## >>>>>>>> Quadratic regression
##
##
## Call:
## robustbase::lmrob(formula = transf.cost ~ Year + I(Year^2), data = yearly.cost.calibration,
## weights = incomplete.year.weights, cov = ".vcov.w")
## \--> method = "MM"
## Residuals:
## Min 1Q Median 3Q Max
## -0.8042717 -0.2154520 -0.0004554 0.3062549 1.3413288
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.385e+03 1.569e+03 -2.157 0.0368 *
## Year 3.382e+00 1.576e+00 2.146 0.0377 *
## I(Year^2) -8.440e-04 3.957e-04 -2.133 0.0388 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Robust residual standard error: 0.3184
## Multiple R-squared: 0.4144, Adjusted R-squared: 0.3865
## Convergence in 27 IRWLS iterations
##
## Robustness weights:
## 2 weights are ~= 1. The remaining 43 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.03662 0.89850 0.93780 0.85570 0.98950 0.99890
## Algorithmic parameters:
## tuning.chi bb tuning.psi refine.tol
## 1.548e+00 5.000e-01 4.685e+00 1.000e-07
## rel.tol scale.tol solve.tol eps.outlier
## 1.000e-07 1.000e-10 1.000e-07 2.222e-03
## eps.x warn.limit.reject warn.limit.meanrw
## 7.378e-06 5.000e-01 5.000e-01
## nResample max.it best.r.s k.fast.s k.max
## 500 50 2 1 200
## maxit.scale trace.lev mts compute.rd fast.s.large.n
## 200 0 1000 0 2000
## psi subsampling cov
## "bisquare" "nonsingular" ".vcov.w"
## compute.outlier.stats
## "SM"
## seed : int(0)
## ------------------------------------------------------------------------------------------------------------
##
## ______________________________ Generalized Additive Models _______________________________
##
##
##
## Family: gaulss
## Link function: identity logb
##
## Formula:
## transf.cost ~ s(Year, k = gam.k)
## <environment: 0x0000000033bfd140>
## ~s(Year, k = gam.k)
## <environment: 0x0000000033bfd140>
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.74465 0.05539 49.55 <2e-16 ***
## (Intercept).1 -1.21813 0.11596 -10.51 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(Year) 4.772 5.752 54.96 < 2e-16 ***
## s.1(Year) 2.628 3.211 14.90 0.00248 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Deviance explained = 60.1%
## -REML = 26.933 Scale est. = 1 n = 45
## ------------------------------------------------------------------------------------------------------------
##
## ______________________________ Multiple Adaptive Regression Splines _______________________________
##
##
## Call: earth(formula=transf.cost~Year, data=yearly.cost.calibration,
## weights=incomplete.year.weights, pmethod="backward",
## nprune=mars.nprune, nfold=5, ncross=30, varmod.method="lm")
##
## coefficients
## (Intercept) 2.44011883
## h(Year-1989) 0.06339453
## h(Year-2004) -0.12793838
##
## Selected 3 of 6 terms, and 1 of 1 predictors
## Termination condition: RSq changed by less than 0.001 at 6 terms
## Importance: Year
## Number of terms at each degree of interaction: 1 2 (additive model)
## GCV 0.2141771 RSS 7.615185 GRSq 0.2571143 RSq 0.3860449 CVRSq -0.07600937
##
## Note: the cross-validation sd's below are standard deviations across folds
##
## Cross validation: nterms 2.65 sd 0.52 nvars 1.00 sd 0.00
##
## CVRSq sd MaxErr sd
## -0.076 0.689 2.17 1.02
##
## varmod: method "lm" min.sd 0.0471 iter.rsq 0.077
##
## stddev of predictions:
## coefficients iter.stderr iter.stderr%
## (Intercept) -0.4942044 0.487398 99
## transf.cost 0.3522279 0.18574 53
##
## mean smallest largest ratio
## 95% prediction interval 1.848102 1.431846 2.744785 1.916956
##
## 68% 80% 90% 95%
## response values in prediction interval 71 93 96 96
## ------------------------------------------------------------------------------------------------------------
##
## ______________________________ Quantile regressions _______________________________
##
##
## >>>>>>>> 0.1 quantile
##
## $call
## quantreg::rq(formula = transf.cost ~ Year, tau = 0.1, data = yearly.cost.calibration,
## weights = incomplete.year.weights)
##
## $terms
## transf.cost ~ Year
## attr(,"variables")
## list(transf.cost, Year)
## attr(,"factors")
## Year
## transf.cost 0
## Year 1
## attr(,"term.labels")
## [1] "Year"
## attr(,"order")
## [1] 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: 0x0000000033bfd140>
## attr(,"predvars")
## list(transf.cost, Year)
## attr(,"dataClasses")
## transf.cost Year
## "numeric" "numeric"
##
## $coefficients
## coefficients lower bd upper bd
## (Intercept) -26.45240297 -31.438273755 -15.99774158
## Year 0.01446266 0.009163603 0.01696573
##
## $residuals
## 1 2 3 4 5
## 8.272480e-01 2.070768e-02 6.245018e-03 -8.217645e-03 -2.268031e-02
## 6 7 8 9 10
## 4.229089e-02 2.961687e-02 1.087368e-01 -7.105427e-15 1.094145e+00
## 11 12 13 14 15
## 3.722119e-01 6.230207e-01 3.445156e-01 3.343320e-01 1.495442e-01
## 16 17 18 19 20
## 2.836583e-01 1.174490e-01 2.273977e-01 9.647195e-02 1.168717e-01
## 21 22 23 24 25
## 1.926413e-01 1.561823e-01 1.537637e-01 8.394847e-01 7.859189e-02
## 26 27 28 29 30
## 3.433523e-01 6.717108e-01 7.900643e-01 7.444635e-01 7.512429e-01
## 31 32 33 34 35
## 1.369307e+00 7.244661e-01 5.724714e-01 7.070359e-01 4.459824e-01
## 36 37 38 39 40
## 1.646898e+00 3.844214e-01 1.106602e-01 8.467218e-02 -4.440892e-15
## 41 42 43 44 45
## 3.497003e-01 -1.919879e-01 7.561563e-01 1.514574e+00 -6.624120e-01
##
## $rdf
## [1] 43
##
## $tau
## [1] 0.1
##
## attr(,"class")
## [1] "summary.rq"
## ------------------------------------------------------------------------------------------------------------
##
## >>>>>>>> 0.5 quantile
##
## $call
## quantreg::rq(formula = transf.cost ~ Year, tau = 0.5, data = yearly.cost.calibration,
## weights = incomplete.year.weights)
##
## $terms
## transf.cost ~ Year
## attr(,"variables")
## list(transf.cost, Year)
## attr(,"factors")
## Year
## transf.cost 0
## Year 1
## attr(,"term.labels")
## [1] "Year"
## attr(,"order")
## [1] 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: 0x0000000033bfd140>
## attr(,"predvars")
## list(transf.cost, Year)
## attr(,"dataClasses")
## transf.cost Year
## "numeric" "numeric"
##
## $coefficients
## coefficients lower bd upper bd
## (Intercept) -51.89331847 -65.19504863 -34.84728592
## Year 0.02738079 0.01887511 0.03411257
##
## $residuals
## 1 2 3 4 5
## 8.194585e-01 -5.329071e-15 -2.738079e-02 -5.476157e-02 -8.214236e-02
## 6 7 8 9 10
## -3.008929e-02 -5.568144e-02 1.052041e-02 -1.111346e-01 9.700928e-01
## 11 12 13 14 15
## 2.352411e-01 4.731318e-01 1.817085e-01 1.586068e-01 -3.909908e-02
## 16 17 18 19 20
## 8.209686e-02 -9.703053e-02 -2.220446e-15 -1.438438e-01 -1.363622e-01
## 21 22 23 24 25
## -7.351075e-02 -1.228878e-01 -1.382246e-01 5.345783e-01 -2.392327e-01
## 26 27 28 29 30
## 1.260967e-02 3.280500e-01 4.334854e-01 3.749665e-01 3.688277e-01
## 31 32 33 34 35
## 9.739739e-01 3.162147e-01 1.513019e-01 2.729483e-01 -1.023356e-03
## 36 37 38 39 40
## 1.186974e+00 -8.842065e-02 -3.751000e-01 -4.140061e-01 -5.115964e-01
## 41 42 43 44 45
## -1.748142e-01 -7.294205e-01 2.058055e-01 9.513054e-01 -1.238599e+00
##
## $rdf
## [1] 43
##
## $tau
## [1] 0.5
##
## attr(,"class")
## [1] "summary.rq"
## ------------------------------------------------------------------------------------------------------------
##
## >>>>>>>> 0.9 quantile
##
## $call
## quantreg::rq(formula = transf.cost ~ Year, tau = 0.9, data = yearly.cost.calibration,
## weights = incomplete.year.weights)
##
## $terms
## transf.cost ~ Year
## attr(,"variables")
## list(transf.cost, Year)
## attr(,"factors")
## Year
## transf.cost 0
## Year 1
## attr(,"term.labels")
## [1] "Year"
## attr(,"order")
## [1] 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: 0x0000000033bfd140>
## attr(,"predvars")
## list(transf.cost, Year)
## attr(,"dataClasses")
## transf.cost Year
## "numeric" "numeric"
##
## $coefficients
## coefficients lower bd upper bd
## (Intercept) -57.11428758 -1.180635e+02 -13.23632039
## Year 0.03044699 8.324736e-03 0.06100015
##
## $residuals
## 1 2 3 4 5
## -7.993606e-15 -8.225247e-01 -8.529717e-01 -8.834187e-01 -9.138657e-01
## 6 7 8 9 10
## -8.648788e-01 -8.935372e-01 -8.304015e-01 -9.551227e-01 1.230384e-01
## 11 12 13 14 15
## -6.148794e-01 -3.800550e-01 -6.745444e-01 -7.007124e-01 -9.014845e-01
## 16 17 18 19 20
## -7.833547e-01 -9.655483e-01 -8.715840e-01 -1.018494e+00 -1.014079e+00
## 21 22 23 24 25
## -9.542934e-01 -1.006737e+00 -1.025140e+00 -3.554029e-01 -1.132280e+00
## 26 27 28 29 30
## -8.835040e-01 -5.711298e-01 -4.687607e-01 -5.303458e-01 -5.395508e-01
## 31 32 33 34 35
## 6.252918e-02 -5.982962e-01 -7.662752e-01 -6.476950e-01 -9.247329e-01
## 36 37 38 39 40
## 2.601986e-01 -1.018263e+00 -1.308008e+00 -1.349980e+00 -1.450637e+00
## 41 42 43 44 45
## -1.116921e+00 -1.674593e+00 -7.424336e-01 -4.440892e-15 -2.192971e+00
##
## $rdf
## [1] 43
##
## $tau
## [1] 0.9
##
## attr(,"class")
## [1] "summary.rq"
## ------------------------------------------------------------------------------------------------------------Summary of model summaries
"________________________________ Damage ________________________________"## [1] "________________________________ Damage ________________________________"summary.damage <- prettySummary(pred.costtype.25$Damage)
# adjust file path
write.csv(summary.models,
"d:/r/Projects/invacost_nature/data/damage_modelsummary.csv")
summary.damage| 1 | 2 | 3 | 4 | 5 | 6 |
|---|---|---|---|---|---|
| Ordinary Least Square regression - Linear | |||||
| Estimate | Standard error | t value | p-value | ||
| Intercept | -149.837310710193 | 22.1004012779335 | -6.77984570623161 | 2.68334391659469e-08 | |
| Year | 0.0767310815387082 | 0.0111377909767825 | 6.88925494280327 | 1.86186336818208e-08 | |
| Adjusted R² | R² | ||||
| 0.668107353641688 | 0.675650368331649 | ||||
| _________________ | _________________ | _________________ | _________________ | _________________ | _________________ |
| Ordinary Least Square regression - Quadratic | |||||
| Estimate | Standard error | t value | p-value | ||
| Intercept | -6281.43290852757 | 2736.21141585214 | -2.29566797073366 | 0.0267546188413239 | |
| Year | 6.23321324892053 | 2.74619189797767 | 2.26976609082226 | 0.0284168884325929 | |
| Adjusted R² | R² | ||||
| 0.698884964819263 | 0.712572011872933 | ||||
| _________________ | _________________ | _________________ | _________________ | _________________ | _________________ |
| Robust regression - Linear | |||||
| Estimate | Standard error | t value | p-value | ||
| Intercept | -151.314155155946 | 15.4788032649375 | -9.77557195902232 | 1.71140342425717e-12 | |
| Year | 0.0774175574641614 | 0.00781326942360505 | 9.90847150749357 | 1.14029639579754e-12 | |
| Adjusted R² | R² | ||||
| 0.794343404174341 | 0.799017417715833 | ||||
| Summary of model weights | |||||
| Min | 25% | 50% | 75% | Max | |
| 0.0776888550407571 | 0.769744516452497 | 0.97298572329036 | 0.997494846222581 | 0.999993965602862 | |
| Number of outliers | |||||
| 0 | |||||
| _________________ | _________________ | _________________ | _________________ | _________________ | _________________ |
| Robust regression - Quadratic | |||||
| Estimate | Standard error | t value | p-value | ||
| Intercept | -3560.88783985051 | 2452.29676382075 | -1.45206236552812 | 0.153914165073996 | |
| Year | 3.50273954229397 | 2.46255393148877 | 1.42240114927203 | 0.162294674525059 | |
| Adjusted R² | R² | ||||
| 0.774240389257175 | 0.784502189745485 | ||||
| Summary of model weights | |||||
| Min | 25% | 50% | 75% | Max | |
| 0.141391661701479 | 0.78928760317949 | 0.964283590543045 | 0.993084946806837 | 0.999970205263011 | |
| Number of outliers | |||||
| 0 | |||||
| _________________ | _________________ | _________________ | _________________ | _________________ | _________________ |
| Multivariate Adaptive Regression Splines | |||||
| log10(cost) | |||||
| (Intercept) | 4.22250744506421 | ||||
| h(2004-Year) | -0.0916263325323822 | ||||
| Generalized R² | R² | Generalized Cross-Validation | Root Sum of Squares | ||
| 0.712973973132936 | 0.738474219321539 | 0.441256217256249 | 17.2972437164449 | ||
| Variance model | |||||
| Estimate | Standard error (last iteration) | Standard error/coefficient % | |||
| Intercept | -0.0473092502386516 | 0.133685470586444 | 282.577867778641 | ||
| Intercept | 0.229188694668325 | 0.0597492101934328 | 26.0698767362413 | ||
| R² for last iteration | |||||
| 0.254942998393425 | |||||
| _________________ | _________________ | _________________ | _________________ | _________________ | _________________ |
| Generalized Additive Models | |||||
| Parametric coefficients | |||||
| Estimate | Standard error | z value | p-value | ||
| Intercept (mean) | 2.85559708691651 | 0.0742923420982404 | 38.4373006189577 | 0 | |
| Intercept (sd) | -1.66064999943417 | 0.314773507683936 | -5.27569810958052 | 1.32251651108549e-07 | |
| Smooth terms | |||||
| Estimated degree of freedom | Residual degree of freedom | Chi² | p-value | ||
| smooth (mean) | 4.62753439068499 | 5.19886066584045 | 2815.94314093181 | 0 | |
| smooth (sd) | 7.29008628336001 | 8.05750830846531 | 59.0660301399824 | 0 | |
| Explained deviance (%) | |||||
| 99.9602856254404 | |||||
| _________________ | _________________ | _________________ | _________________ | _________________ | _________________ |
| Quantile regression | |||||
| Coefficients quantile 0.1 | Coefficients quantile 0.5 | Coefficients quantile 0.9 | |||
| Intercept | -115.505029838761 | -156.016079212433 | -206.812935309598 | ||
| Year | 0.0590978034771958 | 0.0797932229814547 | 0.105919625346723 | ||
| _________________ | _________________ | _________________ | _________________ | _________________ | _________________ |
"______________________________ Management ______________________________"## [1] "______________________________ Management ______________________________"summary.management <- prettySummary(pred.costtype.25$Management)
# adjust file path
write.csv(summary.management,
"d:/r/Projects/invacost_nature/data/management_modelsummary.csv")
summary.management| 1 | 2 | 3 | 4 | 5 | 6 |
|---|---|---|---|---|---|
| Ordinary Least Square regression - Linear | |||||
| Estimate | Standard error | t value | p-value | ||
| Intercept | -41.8699109182845 | 12.389041160915 | -3.37959252652868 | 0.00155322591144634 | |
| Year | 0.0223953363369852 | 0.00623702994903359 | 3.59070527478473 | 0.000840807438103559 | |
| Adjusted R² | R² | ||||
| 0.290793126632316 | 0.3069114646634 | ||||
| _________________ | _________________ | _________________ | _________________ | _________________ | _________________ |
| Ordinary Least Square regression - Quadratic | |||||
| Estimate | Standard error | t value | p-value | ||
| Intercept | -2033.94255084315 | 2263.22135175241 | -0.898693602933832 | 0.373939391447905 | |
| Year | 2.02255328677198 | 2.27515469158555 | 0.888973964826307 | 0.379080482061596 | |
| Adjusted R² | R² | ||||
| 0.295677643256896 | 0.327692295836128 | ||||
| _________________ | _________________ | _________________ | _________________ | _________________ | _________________ |
| Robust regression - Linear | |||||
| Estimate | Standard error | t value | p-value | ||
| Intercept | -45.465163518976 | 15.1424756951875 | -3.00249209139728 | 0.00444820485523086 | |
| Year | 0.0241807231113921 | 0.0076245831540984 | 3.17141575121971 | 0.00279720002911644 | |
| Adjusted R² | R² | ||||
| 0.420325806018833 | 0.433500219518405 | ||||
| Summary of model weights | |||||
| Min | 25% | 50% | 75% | Max | |
| 0.0810045951398613 | 0.863461592531261 | 0.962184914099026 | 0.982982105741942 | 0.999964700941953 | |
| Number of outliers | |||||
| 0 | |||||
| _________________ | _________________ | _________________ | _________________ | _________________ | _________________ |
| Robust regression - Quadratic | |||||
| Estimate | Standard error | t value | p-value | ||
| Intercept | -3384.81306645313 | 1569.28726697328 | -2.15691106254975 | 0.0367830531785806 | |
| Year | 3.38184884354167 | 1.57605964928108 | 2.14576196090313 | 0.0377173960385713 | |
| Adjusted R² | R² | ||||
| 0.386501025158076 | 0.414387342196346 | ||||
| Summary of model weights | |||||
| Min | 25% | 50% | 75% | Max | |
| 0.0366208634124731 | 0.900695907099135 | 0.950316793432793 | 0.990424183716609 | 0.999999813618279 | |
| Number of outliers | |||||
| 0 | |||||
| _________________ | _________________ | _________________ | _________________ | _________________ | _________________ |
| Multivariate Adaptive Regression Splines | |||||
| log10(cost) | |||||
| (Intercept) | 2.44011882793484 | ||||
| h(Year-2004) | -0.12793838124346 | ||||
| h(Year-1989) | 0.0633945281103967 | ||||
| Generalized R² | R² | Generalized Cross-Validation | Root Sum of Squares | ||
| 0.257114312200837 | 0.386044886116394 | 0.214177080668753 | 7.61518509044455 | ||
| Variance model | |||||
| Estimate | Standard error (last iteration) | Standard error/coefficient % | |||
| Intercept | -0.494204389961856 | 0.487398149714687 | 98.6227883876761 | ||
| Intercept | 0.352227865922781 | 0.185739821768437 | 52.7328583960352 | ||
| R² for last iteration | |||||
| 0.0771769074685341 | |||||
| _________________ | _________________ | _________________ | _________________ | _________________ | _________________ |
| Generalized Additive Models | |||||
| Parametric coefficients | |||||
| Estimate | Standard error | z value | p-value | ||
| Intercept (mean) | 2.74464921986815 | 0.0553898608416214 | 49.5514734676089 | 0 | |
| Intercept (sd) | -1.21812534975943 | 0.115958602494781 | -10.5048295128794 | 8.20709583343242e-26 | |
| Smooth terms | |||||
| Estimated degree of freedom | Residual degree of freedom | Chi² | p-value | ||
| smooth (mean) | 4.7718290528099 | 5.75217717433698 | 54.9642402692054 | 0 | |
| smooth (sd) | 2.62778040023576 | 3.21141480314004 | 14.8956606996087 | 0.0024799566325846 | |
| Explained deviance (%) | |||||
| 60.1304608804298 | |||||
| _________________ | _________________ | _________________ | _________________ | _________________ | _________________ |
| Quantile regression | |||||
| Coefficients quantile 0.1 | Coefficients quantile 0.5 | Coefficients quantile 0.9 | |||
| Intercept | -26.4524029726533 | -51.8933184748162 | -57.1142875796834 | ||
| Year | 0.0144626626733758 | 0.0273807870684787 | 0.030446993965438 | ||
| _________________ | _________________ | _________________ | _________________ | _________________ | _________________ |
Model RMSE
"________________________________ Damage ________________________________"## [1] "________________________________ Damage ________________________________"pred.costtype.25$Damage$RMSE## RMSE.calibration RMSE.alldata
## ols.linear 0.6904491 0.8699918
## ols.quadratic 0.6499642 0.7495925
## robust.linear 0.6991168 0.8621573
## robust.quadratic 0.6665810 0.7813917
## mars 0.6199864 0.7255249
## gam 0.6519571 0.7325923
## qt0.1 1.0762641 1.0755920
## qt0.5 0.6960889 0.8732677
## qt0.9 1.4087360 1.7022995"______________________________ Management ______________________________"## [1] "______________________________ Management ______________________________"pred.costtype.25$Management$RMSE## RMSE.calibration RMSE.alldata
## ols.linear 0.4370791 0.5139826
## ols.quadratic 0.4304768 0.4819343
## robust.linear 0.4394067 0.5163892
## robust.quadratic 0.4434196 0.4621144
## mars 0.4113713 0.4208931
## gam 0.3867199 0.4121480
## qt0.1 0.5910992 0.5921148
## qt0.5 0.4514046 0.5291704
## qt0.9 0.9128215 1.0341218Observed cost trends over time
Damage
obs.costtype$Damage##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 402 (number of individual year values: 1109)
## - Cost values in US$ millions:
## o Total cost over the entire period 892,179.56
## o Average annual cost over the entire period 18,587.07
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 846.73 84.67 8
## 2 1980 1989 10 5,123.95 512.40 23
## 3 1990 1999 10 186,202.35 18,620.24 66
## 4 2000 2009 10 562,118.21 56,211.82 209
## 5 2010 2017 8 137,888.32 17,236.04 148
## number_year_values
## 1 26
## 2 53
## 3 188
## 4 656
## 5 186Management
obs.costtype$Management##
## Average annual cost of invasive species over time periods
##
## - Temporal interval of data : [1970, 2017]
## - Values transformed in US$ million: Yes
## - Number of cost estimates: 873 (number of individual year values: 2299)
## - Cost values in US$ millions:
## o Total cost over the entire period 66,013.20
## o Average annual cost over the entire period 1,375.27
## o Average annual cost over each period
##
## initial_year final_year time_span total_cost annual_cost number_estimates
## 1 1970 1979 10 3,649.05 364.91 15
## 2 1980 1989 10 3,427.52 342.75 61
## 3 1990 1999 10 9,487.23 948.72 201
## 4 2000 2009 10 30,189.54 3,018.95 552
## 5 2010 2017 8 19,259.85 2,407.48 177
## number_year_values
## 1 38
## 2 193
## 3 558
## 4 1239
## 5 271Estimated cost trends over time
Damage
pred.costtype.25$Damage##
## Estimation of annual cost values of invasive alien species over time
##
## - Temporal interval of data: [1970, 2017]
## - Temporal interval used for model calibration: [1970, 2014]
## - Cost transformation: log10
## - Values transformed in US$ million: Yes
## - Estimated average annual cost of invasive alien species in 2017:
##
## o Linear regression:
## . Linear: US$ million 84,972.96
## . Quadratic: US$ million 16,753.60
## o Robust regression:
## . Linear: US$ million 68,716.06
## . Quadratic: US$ million 25,209.51
## o Multiple Adapative Regression Splines: US$ million 16,691.96
## o Generalized Additive Model: US$ million 11,959.28
## o Quantile regression:
## . Quantile 0.1: US$ million 4,957.24
## . Quantile 0.5: US$ million 84,498.99
## . Quantile 0.9: US$ million 6,713,500.33
##
## You can inspect the summary of each fitted model with summary(object)Management
pred.costtype.25$Management##
## Estimation of annual cost values of invasive alien species over time
##
## - Temporal interval of data: [1970, 2017]
## - Temporal interval used for model calibration: [1970, 2014]
## - Cost transformation: log10
## - Values transformed in US$ million: Yes
## - Estimated average annual cost of invasive alien species in 2017:
##
## o Linear regression:
## . Linear: US$ million 2,002.08
## . Quadratic: US$ million 1,181.36
## o Robust regression:
## . Linear: US$ million 2,029.34
## . Quadratic: US$ million 570.64
## o Multiple Adapative Regression Splines: US$ million 356.42
## o Generalized Additive Model: US$ million 555.49
## o Quantile regression:
## . Quantile 0.1: US$ million 523.34
## . Quantile 0.5: US$ million 2,156.40
## . Quantile 0.9: US$ million 19,828.93
##
## You can inspect the summary of each fitted model with summary(object)Global appendix : RMSEs
ATTENTION Supprimer la valeur RMSE du modèle quadratique OLS global, car il est non significatif donc non évalué !
RMSEs <- data.frame(Global = global.trend$RMSE[, 1],
Plants = pred.impact$Plants$Plants$RMSE[, 1],
Invertebrates = pred.impact$Invertebrates$Invertebrates$RMSE[, 1],
Vertebrates = pred.impact$Vertebrates$Vertebrates$RMSE[, 1],
# Africa = pred.regions$Africa$RMSE[, 1],
# Asia = pred.regions$Asia$RMSE[, 1],
# Europe = pred.regions$Europe$RMSE[, 1],
# `Central America` = pred.regions$`Central America`$RMSE[, 1],
# `North America` = pred.regions$`North America`$RMSE[, 1],
# `Oceania-Pacific islands` = pred.regions$`Oceania-Pacific islands`$RMSE[, 1],
# `South America` = pred.regions$`Oceania-Pacific islands`$RMSE[, 1],
Damage.25 = pred.costtype.25$Damage$RMSE[, 1],
Management.25 = pred.costtype.25$Management$RMSE[, 1]
)
RMSEs <- RMSEs[-grep("qt", rownames(RMSEs)), ]
data.frame(RMSEs)| Global | Plants | Invertebrates | Vertebrates | Damage.25 | Management.25 | |
|---|---|---|---|---|---|---|
| ols.linear | 0.4943647 | 0.8827656 | 0.3696940 | 0.9771180 | 0.6904491 | 0.4370791 |
| ols.quadratic | 0.4941152 | 0.6962406 | 0.3668752 | 0.9352001 | 0.6499642 | 0.4304768 |
| robust.linear | 0.5205157 | 0.9744516 | 0.3813940 | 0.9855533 | 0.6991168 | 0.4394067 |
| robust.quadratic | 0.5343973 | 0.8540153 | 0.3788097 | 0.9585592 | 0.6665810 | 0.4434196 |
| mars | 0.3834539 | 0.5877908 | 0.3085550 | 0.9064081 | 0.6199864 | 0.4113713 |
| gam | 0.4965620 | 0.6417153 | 0.3467522 | 0.9502886 | 0.6519571 | 0.3867199 |
# adjust file path
write.table(RMSEs,
"d:/r/projects/invacost_nature/data/RMSE.csv", sep = ";")High margin with relaxed assumptions in a vein similar to Pimentel et al.
# Expanding and formating the database
# db.over.time.high <- expandYearlyCosts(invacost,
# startcolumn = "Probable_starting_year_high_margin",
# endcolumn = "Probable_ending_year_high_margin")
# # global.raw <- calculateRawAvgCosts(db.over.time.high, minimum.year = 1970)
# global.pred <- costTrendOverTime(db.over.time.high, minimum.year = 1970)
# plot(global.pred, graphical.parameters =
# ylab(paste0("Annual cost in US$ millions"))) +
# theme_minimal() + xlab("Year") +
# scale_x_continuous(breaks = global.raw$year.breaks) +
# scale_y_log10(breaks = 10^(-15:15),
# labels = scales::comma) +
# annotation_logticks(sides = "l") +
# guides(col = guide_legend(override.aes = list(size = 2)),
# shape = guide_legend(override.aes = list(size = 2))) +
# theme(panel.border = element_blank(),
# axis.line = element_line(color = "black"),
# axis.text = element_text(size = 12),
# axis.title = element_text(size = 14),
# panel.grid.minor = element_blank(),
# legend.text = element_text(size = 14),
# legend.title = element_text(size = 18),
# strip.text = element_text(size = 14))
# global.pred
# global.pred$estimated.annual.costs[which(global.pred$estimated.annual.costs$Year == 1998), ]Global Appendix: Analyses for the revised manuscript
Difference between models with and without incompleteness weights
# Creating the vector of weights
year_weights <- rep(1, length(1960:2017))
names(year_weights) <- 1960:2017
# # Below 25% the weight does not matter because years will be removed
year_weights[names(year_weights) >= (2017 - quantiles["25%"])] <- 0
# Between 25 and 50%, assigning 0.25 weight
year_weights[names(year_weights) >= (2017 - quantiles["50%"]) &
names(year_weights) < (2017 - quantiles["25%"])] <- .25
# Between 50 and 75%, assigning 0.5 weight
year_weights[names(year_weights) >= (2017 - quantiles["75%"]) &
names(year_weights) < (2017 - quantiles["50%"])] <- .5
global.trend_weights <- modelCosts(
db.over.time, # The EXPANDED database
minimum.year = 1970,
maximum.year = 2017,
final.year = 2017,
# Some years are so incomplete that we eliminate with our 25% threshold (see above)
incomplete.year.threshold = 2017 - quantiles["25%"],
# For the other incomplete years we apply the vector of weights that we defined above
incomplete.year.weights = year_weights)
pg1 <- plot(global.trend, graphical.parameters =
ylab(paste0("Annual cost in US$ millions"))) +
theme_minimal() + xlab("Year") +
scale_x_continuous(breaks = global.raw$year.breaks) +
scale_y_log10(breaks = 10^(-15:15),
labels = scales::comma) +
annotation_logticks(sides = "l") +
guides(col = guide_legend(override.aes = list(size = 2)),
shape = guide_legend(override.aes = list(size = 2))) +
theme(panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14),
panel.grid.minor = element_blank(),
legend.text = element_text(size = 14),
legend.title = element_text(size = 18),
strip.text = element_text(size = 14))
pg2 <- plot(global.trend_weights, graphical.parameters =
ylab(paste0("Annual cost in US$ millions"))) +
theme_minimal() + xlab("Year") +
scale_x_continuous(breaks = global.raw$year.breaks) +
scale_y_log10(breaks = 10^(-15:15),
labels = scales::comma) +
annotation_logticks(sides = "l") +
guides(col = guide_legend(override.aes = list(size = 2)),
shape = guide_legend(override.aes = list(size = 2))) +
theme(panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14),
panel.grid.minor = element_blank(),
legend.text = element_text(size = 14),
legend.title = element_text(size = 18),
strip.text = element_text(size = 14))
pleg <- cowplot::get_legend(pg1)
pg1 <- pg1 + theme(legend.position='none')
pg2 <- pg2 + theme(legend.position='none')
plot_grid(pg1, pg2, pleg, ncol = 2, align = "hv", labels = "auto")
matrix(data = scales::comma(c(global.trend$final.year.cost,
global.trend_weights$final.year.cost)),
ncol = 2,
nrow = length(global.trend$final.year.cost),
dimnames = list(names(global.trend$final.year.cost),
c("Without.weights", "With.weights")))## Without.weights With.weights
## ols.linear "91,980" "123,981"
## ols.quadratic "102,621" "255,141"
## robust.linear "46,823" "38,086"
## robust.quadratic "86,535" "91,994"
## mars "162,719" "5,639"
## gam "78,795" "78,795"
## quantile.0.1 "16,224" "9,452"
## quantile.0.5 "57,493" "50,425"
## quantile.0.9 "2,219,916" "4,319,907"