This vignette accompanies a manuscript we will post to PeerJ Preprints and submit to Conservation Letters. The preprint--which will be updated in the near future--can be found here.
Recovery planning is an essential part of implementing the U.S. Endangered Species Act (ESA), but conservationists and government agencies recognize there are problems with the process. A common perception is that too many species lack recovery plans, these plans take too long to write, and plans are rarely updated to include new information. Using data from all U.S. domestic and transboundary ESA-listed species--those species we expect should have recovery plans--we quantify basic characteristics of ESA recovery planning. We show that ~1/4 of these listed species lack recovery plans; the average recovery plan has taken >5 years to be developed after listing; and half of recovery plans are 20 or more years old. These results are not unexpected given dwindling budgets and more species to protect, but they underscore the need for systematic improvements to recovery planning. We discuss how many of the shortcomings we identify here can be ameliorated, at least in part, by transitioning to modern, web-based recovery plans.
# library(ecosscraper) library(dplyr) library(extrafont) library(ggplot2) library(ggthemes) library(knitr) library(lubridate) library(plotly) library(recovery.plan.overview) library(skimr) library(stringr) library(tibble) library(viridis)
# have <- readRDS("data/fin_have_plans_2018-01-08.rds") # miss <- readRDS("data/missing_plans_2018-01-08.rds") # tecp <- readRDS("data/TECP_2018-01-08.rds") have <- readRDS("../data/fin_have_plans_2018-01-08.rds") miss <- readRDS("../data/missing_plans_2018-01-08.rds") tecp <- readRDS("../data/TECP_2018-01-08.rds") # manual fix...unclear how I missed this: miss <- filter(miss, common_name != "Sperm whale") miss <- filter(miss, common_name != "North Atlantic Right Whale") spwh <- list(as.Date("2010-12-01"), "Final Recovery Plan for the Sperm Whale", "Final", "Physeter macrocephalus", "Sperm whale", as.Date("1970-06-02"), "Mammals", "NMFS", "Endangered", "Wherever found", "https://ecos.fws.gov/ecp0/profile/speciesProfile?sId=8363", 14792, 40.52, "Physeter macrocephalusWherever found") narw <- list( as.Date("2005-05-26"), "RECOVERY PLAN FOR THE NORTH ATLANTIC RIGHT WHALE (EUBALAENA GLACIALIS) REVISION", "Final Revision 2", "Eubalaena glacialis", "North Atlantic right whale", as.Date("1970-06-02"), "Mammals", "NMFS", "Endangered", "Wherever found", "https://ecos.fws.gov/ecp0/profile/speciesProfile?sId=8363", 12777, 35.01, "Eubalaena glacialisWherever found") have <- rbind(have, spwh) have <- rbind(have, narw) final_have <- filter(have, plan_status == "Final") tecp_elig <- filter( tecp, as.numeric(as.Date("2018-01-08") - first_listed) > 365*2.5 ) have$year_plan <- year(have$date) have$year_list <- year(have$first_listed) miss$year_list <- year(miss$first_listed) tecp$year_list <- year(tecp$first_listed) official <- filter(have, grepl(have$plan_status, pattern = "Final|Revision")) drafts <- dplyr::setdiff(have, official) %>% filter(plan_status != "Exempt") official$plan_age <- as.numeric(as.Date("2018-01-08") - official$date) by_plan <- distinct(have, title, .keep_all = TRUE)
# Have to use a combination of distinct w/ population or title because of # quirkiness with join on species and how FWS presents their data tables on ECOS w_final <- filter(have, plan_status == "Final") %>% distinct(species, where_listed, .keep_all = TRUE) str_c("# species with final plans: ", dim(w_final)[1], "\n\n") %>% cat() revs <- c("Final Revision 1", "Final Revision 2", "Final Revision 3", "Draft Revision 1", "Draft Revision 2") w_revision <- filter(have, plan_status %in% revs) %>% distinct(species, where_listed, .keep_all = TRUE) str_c("# species with revised plans: ", dim(w_revision)[1], "\n\n") %>% cat() w_plan <- rbind(w_final, w_revision) str_c("# species with final or revised plans: ", dim(w_plan)[1], "\n\n") %>% cat() w_draft <- filter(have, plan_status == "Draft") %>% distinct(title, species, .keep_all = TRUE) str_c("# species with draft plans: ", dim(w_draft)[1], "\n\n") %>% cat() unofficial <- c("Conservation Strategy", "Draft", "Outline") unoffic <- filter(have, plan_status %in% unofficial) %>% distinct(title, species, .keep_all = TRUE) str_c("# species with draft or outline plans: ", dim(unoffic)[1], "\n\n") %>% cat()
loc_summary <- function(x) { c(min = round(min(x, na.rm = TRUE), 1), median = round(median(x, na.rm = TRUE), 1), mean = round(mean(x, na.rm = TRUE), 1), max = round(max(x, na.rm = TRUE))) } loc_sum_2 <- function(x) { c(min = as.character(min(x, na.rm = TRUE)), median = as.character(median(x, na.rm = TRUE)), mean = as.character(mean(x, na.rm = TRUE)), max = as.character(max(x, na.rm = TRUE))) } delta_tmp <- loc_summary(w_final$elapsed_years) listd_tmp <- loc_sum_2(w_final$first_listed) pland_tmp <- loc_sum_2(w_final$date) rev_listd_tmp <- loc_sum_2(w_revision$first_listed) rev_pland_tmp <- loc_sum_2(w_revision$date) draft_tmp <- loc_sum_2(w_draft$date) draft_del <- loc_summary(w_draft$elapsed_years) placehold <- c(" ", " ", " ", " ") to_show_df <- rbind(placehold, listd_tmp, pland_tmp, delta_tmp, placehold, rev_listd_tmp, rev_pland_tmp, placehold, draft_tmp, draft_del) %>% as_data_frame() guide <- c("Spp. w/ Final Plans", "Listed Date", "Plan Date", "Years Elapsed", "Spp. w/ Revised Plans", "Listed Date", "Plan Date", "Spp. w/ Draft Plans", "Draft Date", "Years Elapsed") to_show_df <- cbind(guide, to_show_df) kable(to_show_df, caption = "Summary statistics of time between listing and final recovery plan. Min, median, mean, and max are given in years.", digits = 1, align = "c") # write.table(data.frame(to_show_df), # file = "vignettes/summary_stats.tsv", # row.names = FALSE, # sep = "\t", # quote = FALSE)
reg_elig <- table(tecp_elig$lead_region) %>% as_data_frame() plan_elig <- table(w_plan$lead_region) %>% as_data_frame() res <- left_join(plan_elig, reg_elig, by = "Var1") names(res) <- c("region", "with_plan", "n_eligible") res$prop_plan <- round(res$with_plan / res$n_eligible, 3) * 100 row.names(res) <- NULL kable(res, caption = "Distribution of plans among FWS regions and NMFS.", align = "c", digits = 3) # write.table(data.frame(res), # file = "vignettes/summary_by_region.tsv", # row.names = FALSE, # sep = "\t", # quote = FALSE)
age_by_spp <- median(official$plan_age, na.rm = TRUE) cat("Age by species: ", age_by_spp, "\n") qq <- distinct(official, title, .keep_all = TRUE) plan_age <- median(qq$plan_age, na.rm = TRUE) cat("Age by plan: ", plan_age, "\n")
Per Cons. Lett. editor suggestion, adding a couple of analyses of multispecies plans.
multisp <- table(have$title) %>% sort(decreasing = TRUE) %>% as.data.frame(stringsAsFactors = FALSE) %>% rename(title = Var1, n_spp = Freq) plan_yrs <- data_frame(title = have$title, year = have$year_plan) %>% distinct(.keep_all = TRUE) multisp <- left_join(multisp, plan_yrs, by = "title") singles <- filter(multisp, n_spp == 1) multis <- filter(multisp, n_spp > 1) have$multispp <- ifelse( have$title %in% multis$title, TRUE, FALSE ) qplot(data = have, x = multispp, y = elapsed, geom = "boxplot") + labs(x = "Multispecies plan?", y = "Time-to-plan (days)") + theme_hc() yrtab_byspp <- table(have$year_plan, have$multispp) %>% as.data.frame.matrix() %>% rownames_to_column("year") %>% mutate(pct_multi = `TRUE` / (`TRUE` + `FALSE`)) yrtab_byspp$year <- as.numeric(yrtab_byspp$year) ggplot(data = yrtab_byspp, aes(x = year, y = pct_multi)) + geom_bar(stat = "identity") + labs(x = "Year", y = "Percent of species covered by multispecies plan") + theme_hc() + scale_x_continuous(breaks = seq(1970, 2017, 5)) single_yr <- table(singles$year) %>% as.data.frame(stringsAsFactors = FALSE) %>% rename(year = Var1, n_singles_plans = Freq) multi_yr <- table(multis$year) %>% as.data.frame(stringsAsFactors = FALSE) %>% rename(year = Var1, n_multi_plans = Freq) by_year_multi <- full_join(single_yr, multi_yr, by = "year") by_year_multi$n_multi_plans <- ifelse( is.na(by_year_multi$n_multi_plans), 0, by_year_multi$n_multi_plans ) by_year_multi$pct_single <- with( data = by_year_multi, n_singles_plans / (n_singles_plans + n_multi_plans) ) by_year_multi$pct_multi <- with( data = by_year_multi, 1 - pct_single ) by_year_multi$year <- as.numeric(by_year_multi$year) ggplot(data = by_year_multi, aes(x = year, y = pct_multi)) + geom_bar(stat = "identity") + labs(x = "Year", y = "Proportion multi-species plans") + theme_hc() + scale_x_continuous(breaks = seq(1970, 2017, 5)) mod1 <- lm(have$elapsed ~ have$year_plan + have$multispp) summary(mod1) hist(resid(mod1)) qqnorm(resid(mod1))
Across all species and all plans, it looks like multispecies plans have shorter time-to-plan by about 734 days (median):
tapply(have$elapsed, have$multispp, median)
Now, basic models for partitioning variation among regions, taxonomic groups:
mod <- lm(final_have$elapsed_years ~ final_have$lead_region) summary(mod) anova(mod) hist(resid(mod)) final_have$plan_age <- as.numeric(as.Date("2018-01-08") - final_have$date ) mod <- lm(final_have$plan_age ~ final_have$lead_region) summary(mod) anova(mod) hist(resid(mod)) cor.test(final_have$plan_age, final_have$elapsed)
taxo_elig <- table(tecp_elig$species_group) %>% as_data_frame() intermed <- distinct(w_plan, species, title, .keep_all = TRUE) plan_taxo <- table(intermed$species_group) %>% as_data_frame() res <- left_join(plan_taxo, taxo_elig, by = "Var1") names(res) <- c("taxon", "with_plan", "n_eligible") res$prop_plan <- round(res$with_plan / res$n_eligible, 3) * 100 row.names(res) <- NULL kable(res, caption = "Distribution of recovery plans among taxonomic groups.", align = "c", digits = 3) # write.table(data.frame(res), # file = "vignettes/summary_by_taxo.tsv", # row.names = FALSE, # sep = "\t", # quote = FALSE)
Time-to-plan and plan age by taxonomic group
qplot(data = final_have, x = species_group, y = elapsed, geom = "boxplot") + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) + theme_hc() + labs(x = "", y = "Time-to-plan (days)") mod2 <- lm(elapsed ~ species_group, data = final_have) summary(mod2) final_have$age <- as.numeric(as.Date("2018-01-08") - final_have$date) qplot(data = final_have, x = species_group, y = age, geom = "boxplot") + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) + theme_hc() + labs(x = "", y = "Plan age (days)") mod3 <- lm(age ~ species_group, data = final_have) summary(mod3) tapply(final_have$age, final_have$species_group, median) have_final_mult <- filter(have, plan_status %in% c("Final", "Final Revision 1", "Final Revision 2")) tapply(have_final_mult$elapsed, have_final_mult$multispp, median)
year <- seq(min(tecp$year_list, na.rm = TRUE), max(tecp$year_list, na.rm = TRUE)) n_listed <- lapply(year, FUN = function(x) { sum(tecp$year_list <= x, na.rm = TRUE) }) %>% unlist() n_official <- lapply(year, function(x) { sum(official$year_plan <= x) }) %>% unlist() n_unofficial <- lapply(year, function(x) { sum(unoffic$year_plan <= x) }) %>% unlist() n_list_df <- data_frame(year, n_listed, n_official, n_unofficial) n_list_df$n_no_plan <- with(n_list_df, n_listed - n_official - n_unofficial) plt <- ggplot(n_list_df, aes(y = n_listed, x = year)) + geom_line(size = 1) + geom_line(aes(y = n_no_plan, x = year), colour = "orange", size = 1) + geom_line(aes(y = n_official, x = year), colour = "blue", linetype = "dashed") + geom_line(aes(y = n_unofficial, x = year), colour = "#998ec3", linetype = "dotted") + labs(x = "", y = "Count") + theme_hc(base_size = 18) + theme(text = element_text(family = "Arial")) ggplotly(plt) # ggsave("inst/figures/Fig1A.pdf", plot = plt)
n_listed <- lapply( year, function(x) sum(tecp$year_list == x, na.rm = TRUE) ) %>% unlist() n_off_spp <- lapply( year, function(x) sum(official$year_plan == x, na.rm = TRUE) ) %>% unlist() n_draft_spp <- lapply( year, function(x) sum(unoffic$year_plan == x, na.rm = TRUE) ) %>% unlist() n_plans <- lapply( year, function(x) sum(by_plan$year_plan == x, na.rm = TRUE) ) %>% unlist() n_list_df <- data_frame(year, n_listed, n_off_spp, n_draft_spp, n_plans) %>% mutate(pct_spp_w_plan = n_off_spp / n_listed) # n_list_df$pct_w_plan <- n_list_df$n_w_plans / n_list_df$n_listed plt <- ggplot(n_list_df, aes(y = n_plans, x = year)) + geom_line(colour = "gray70") + geom_line(aes(y = n_off_spp, x = year)) + geom_line(aes(y = n_draft_spp, x = year), linetype = "dotted", color = "gray70") + labs(x = "Year", y = "Count") + theme_hc() + theme(text = element_text(family = "Arial")) ggplotly(plt) # ggsave("inst/figures/Fig1B.pdf", plot = plt)
pct_w_plan <- lapply( year, function(x) { ls_n <- sum(tecp$year_list == x, na.rm = TRUE) of_n <- sum(official$year_list == x, na.rm = TRUE) ifelse(ls_n == 0, NA, of_n / ls_n) } ) %>% unlist() n_list_df$pct_w_plan <- pct_w_plan plt <- ggplot(n_list_df, aes(y = pct_w_plan, x = year)) + geom_point(size = 3, alpha = 0.5) + geom_smooth(level = 0.9, alpha = 0, colour = "black") + scale_y_continuous(limits = c(0,1)) + labs(x = "Year ESA-listed", y = "Proportion of species with official plans") + theme_hc() + theme(text = element_text(family = "Arial")) ggplotly(plt) # ggsave("inst/figures/Fig2.pdf", plot = plt)
plt <- qplot(official$elapsed_years, geom = "histogram") + labs(x = "Time between date listed and recovery plan date (years)", y = "Number of species") + theme_hc() + theme(plot.title=element_text(hjust = 0, size = 18), text = element_text(family = "Arial")) plt # ggsave("inst/figures/Fig3a.pdf", plot = plt)
time_to_plan <- sort(final_have$elapsed_years, decreasing = FALSE) order <- seq_along(time_to_plan) percentile <- 100 * order / max(order) pctile <- data_frame(time_to_plan, order, percentile) plt <- ggplot(pctile, aes(y = time_to_plan, x = percentile)) + geom_line() + labs(x = "Percentile", y = "Time-to-recovery plan (years)\n") + theme_hc() + theme(text = element_text(family = "Arial")) ggplotly(plt) %>% layout(font = list(family = "Arial")) # ggsave("inst/figures/Fig3b.pdf", plot = plt)
tmp <- final_have tmp$year_listed <- year(tmp$first_listed) tmp$year_plan <- year(tmp$date) years <- seq(min(tmp$year_listed, na.rm = TRUE), max(tmp$year_listed, na.rm = TRUE)) cur_fx <- function(df, x, fn) { sub_df <- filter(df, year_listed == x) if(dim(sub_df)[1] == 0) return(NA) switch(fn, mean = mean(sub_df$elapsed_years, na.rm = TRUE), min = min(sub_df$elapsed_years, na.rm = TRUE), max = max(sub_df$elapsed_years, na.rm = TRUE) ) } # # This is to calculate a potential "minimum time-to-plan" for species # # that do not have final plans or final plans in revision. However, # # the results aren't particularly informative and are not included in # # the paper or the figure below. # get_potential <- function(x, df) { # sub_df <- filter(df, Year_Listed == x) # if(dim(sub_df)[1] == 0) return(NA) # mean = mean(sub_df$Age_Min_Age, na.rm = TRUE) # return(mean / 365) # } n_plans <- lapply(years, function(x) {dim(filter(tmp, year_plan == x))[1]} ) %>% unlist() tm_mean <- unlist(lapply(years, cur_fx, df = tmp, fn = "mean")) tm_min <- unlist(lapply(years, cur_fx, df = tmp, fn = "min")) tm_max <- unlist(lapply(years, cur_fx, df = tmp, fn = "max")) # potential_age <- unlist(lapply(years, get_potential, df = listed)) dat <- data_frame(years, tm_mean, tm_min, tm_max, # potential_age, n_plans) names(dat) <- c("Year", "Mean", "Min", "Max", #"Potential", "N_plans") plt <- ggplot(dat, aes(y = Mean, x = Year)) + geom_line() + geom_line(aes(y = Min, x = Year), colour = "gray70") + geom_line(aes(y = Max, x = Year), colour = "gray70") + # geom_line(aes(y = Potential, x = Year), colour = "red") + labs(x = "Year of Listing", y = "Time-to-plan (years)\n") + theme_hc() + theme(text = element_text(family = "Arial")) ggplotly(plt) %>% layout(font = list(family = "Arial")) # ggsave("inst/figures/Fig4.pdf", plot = plt)
Including one analysis here because the data frame is built in the chunk above:
# Subset to > 1978 because plans were not required before then... mod <- lm(elapsed_years ~ year_listed, data = filter(tmp, year_listed > 1978)) summary(mod) aov(mod) hist(resid(mod))
official$plan_age <- as.numeric(as.Date("2018-01-08") - official$date) / 365 skimr::skim(official$plan_age) plt <- ggplot(official, aes(x = plan_age)) + geom_histogram(color = "white", bins = 40) + labs(x = "Plan age (years)", y = "Number of species") + theme_hc(base_size = 18) + theme(plot.title=element_text(hjust = 0, size = 20), text = element_text(family = "Arial")) plt # ggsave("inst/figures/Fig5a.pdf", plot = plt)
age <- sort(official$plan_age, decreasing = FALSE) order <- seq(1:length(age)) percentile <- 100 * order / max(order) pctile <- data_frame(age, order, percentile) gg <- ggplot(pctile, aes(y = age, x = percentile)) + geom_line() + labs(x = "Percentile", y = "Age of species' plan (years)\n") + theme_hc() + theme(text = element_text(family = "Arial")) ggplotly(gg) %>% layout(font = list(family = "Arial")) # ggsave("inst/figures/Fig5b.pdf", plot = gg)
tmp <- distinct(official, title, .keep_all = TRUE) age <- sort(tmp$plan_age, decreasing = FALSE) order <- seq_along(age) percentile <- 100 * (order / max(order)) pctile <- data_frame(age, order, percentile) gg <- ggplot(pctile, aes(y = age, x = percentile)) + geom_line() + labs(x = "Percentile", y = "Plan age (years)\n") + theme_hc() + theme(text = element_text(family = "Arial")) ggplotly(gg) %>% layout(font = list(family = "Arial")) # ggsave("inst/figures/Fig5c.pdf", plot = gg) skim(tmp$plan_age) skim(official$plan_age)
plt <- ggplot(official, aes(y = elapsed_years, x = lead_region)) + geom_violin(fill = "gray", alpha = 0.8, colour = NA) + # geom_jitter(height = 0, alpha = 0.2) + geom_boxplot(colour = "gray20", fill = NA) + labs(x = "FWS Region / NMFS", y = "Time-to-Plan (years)") + theme_hc() + theme(text = element_text(family = "Arial")) plt # ggsave("inst/figures/Fig6a.pdf", plot = plt)
plt <- ggplot(official, aes(y = plan_age, x = lead_region)) + geom_violin(fill = "gray", alpha = 0.7, colour = NA) + # geom_jitter(height = 0, alpha = 0.2) + geom_boxplot(colour = "gray20", fill = NA) + labs(x = "FWS Region / NMFS", y = "Plan Age (years)") + theme_hc() + theme(text = element_text(family = "Arial")) plt # ggsave("inst/figures/Fig6b.pdf", plot = plt) lm(plan_age ~ lead_region, data = official) %>% summary()
And one quick final analysis here, just to show the negative time-to-plan, plan age correlation:
ttp <- tapply(official$elapsed_years, INDEX = official$lead_region, FUN = median, na.rm = TRUE) age <- tapply(official$plan_age, INDEX = official$lead_region, FUN = median, na.rm = TRUE) cor.test(age, ttp)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.