In Defenders-ESC/threatdemog: Analysis of Threat and Demography Scoring Manuscript

library(ggplot2)
library(ggthemes)
library(threatdemog)

Are FWS's current status metrics sufficient for accurately reflecting the conservation status of ESA-listed species?

Multinomial models

First, a naive model that assumes the response variable is the threat or demography score, even though this implies causality flows the opposite direction.

# Threats and demography by status change
amod <- lm(change$Threat_allev ~ change$status_cat + 0)
summary(amod)

bmod <- lm(change$Improve_demo ~ change$status_cat + 0)
summary(bmod)

cmod <- lm(change$combined_components ~ change$status_cat + 0)
summary(cmod)

But the response is actually FWS's recommended status change, which should be a function of threat and demographic status. This requires using multinomial models. Because we are interested in the relationship between threat and demography scores individually and combined, we consider four models:

change$status_cat <- relevel(factor(change$status_cat), ref="No Change")
combo_mod <- nnet::multinom(status_cat ~ 0 + combined_components, 
                            data = change,
                            trace = FALSE)
combo_p <- get_multinom_p(combo_mod)
thr_mod <- nnet::multinom(status_cat ~ 0 + Threat_allev, 
                          data = change,
                          trace = FALSE)
thr_p <- get_multinom_p(thr_mod)
demo_mod <- nnet::multinom(status_cat ~ 0 + Improve_demo, 
                           data = change,
                           trace = FALSE)
demo_p <- get_multinom_p(demo_mod)
thr_demo_mod <- nnet::multinom(status_cat ~ 0 + Threat_allev + Improve_demo, 
                               data = change,
                               trace = FALSE)
thr_demo_p <- get_multinom_p(thr_demo_mod)

AICcmodavg::aictab(list(combo_mod, thr_mod, demo_mod, thr_demo_mod))

Model 1 (threat and demography scores added for a single predictor) is the most parsimonious model given the data. The details of the model are:

summary(combo_mod)

# cat(paste0("p-values: ", thr_p[,1]))
cat("p-values:\n")
combo_p

While more parsimonious, model 4 (threat and demography scores as separate predictors) is nearly as good a fit with a deltaAICc of 3.61. Model 4 has K = 4 vs. K = 2 for model 1; with 2 AIC units added per parameter, this means the model fit is essentially the same excepting the additional parameters. Because model 1 does not allow administrators to determine threat vs. demographic status independently, we think model 4 is better from a useability standpoint if not strictly statistically better. The summary for model 4 is:

summary(thr_demo_mod)

# cat(paste0("p-values: ", thr_p[,1]))
cat("p-values:\n")
thr_demo_p

We can plot the data by FWS recommendation and threat and demography scores to 'see' how things shake out:

change$tmp_stat_cat <- ifelse(change$status_cat == "Degrade",
                                   "Uplist",
                                   ifelse(change$status_cat == "Improve",
                                          "Down- or de-list",
                                          "No change"))
change$tmp_stat_cat <- as.factor(change$tmp_stat_cat)
change$tmp_stat_cat <- relevel(change$tmp_stat_cat, ref="No change")

aplot <- ggplot(change, aes(tmp_stat_cat, Threat_allev)) +
         geom_boxplot(outlier.shape=NA, alpha=0.5) +
         geom_jitter(alpha=0.3, size=4,
                     position=position_jitter(height=0.05, width=0.5)) +
         labs(x="\nFWS status change recommendation",
              y="",
              title="Threats\n") + 
         theme_minimal(base_size=14) +
         theme(axis.text.x=element_text(vjust=1.5),
               axis.ticks=element_blank(),
               text=element_text(size=14, family="Open Sans"),
               legend.position="none")

bplot <- ggplot(change, aes(tmp_stat_cat, Improve_demo)) +
         geom_boxplot(outlier.shape=NA, alpha=0.5) +
         geom_jitter(alpha=0.3, size=4,
                     position=position_jitter(height=0.05, width=0.4)) +
         labs(x="\n",
              y="",
              title="Demography\n") + 
         theme_minimal(base_size=14) +
         theme(axis.text.x=element_text(vjust=1.5),
               axis.ticks=element_blank(),
               text=element_text(size=14, family="Open Sans"),
               legend.position="none")

cplot <- ggplot(change, aes(tmp_stat_cat, combined_components)) +
         geom_boxplot(outlier.shape=NA, alpha=0.5) +
         geom_jitter(alpha=0.3, size=4,
                     position=position_jitter(height=0.05, width=0.4)) +
         labs(x="\n",
              y="Score",
              title="Threats + Demography\n") + 
         theme_minimal(base_size=14) +
         theme(axis.text.x=element_text(vjust=1.5),
               axis.ticks=element_blank(),
               text=element_text(size=14, family="Open Sans"),
               legend.position="none")
cplot
bplot
aplot

# multiplot(cplot, aplot, bplot, cols=3)

This basically comports with the multinomial models: there is some predictability of FWS status change recommendations, but those recommendations don't always seem to be consistent with threats and demography individually.

---

Discriminant function analysis

The second evaluation of whether FWS's status change recommendations are consistent with threat and demography change scores is done with discriminant function analysis (or linear discriminant analysis). We consider the same four model structures as with the multinomial models.

dfa_combo <- MASS::lda(status_cat ~ combined_components,
                       data = change,
                       na.action = "na.omit",
                       CV=TRUE)
combo_comps <- get_lda_comps(change, dfa_combo)

dfa_both <- MASS::lda(status_cat ~ Threat_allev + Improve_demo,
                      data = change,
                      na.action = "na.omit",
                      CV=TRUE)
both_comps <- get_lda_comps(change, dfa_both)

dfa_threat <- MASS::lda(status_cat ~ Threat_allev,
                        data = change,
                        na.action = "na.omit",
                        CV=TRUE)
thr_comps <- get_lda_comps(change, dfa_threat)

dfa_demog <- MASS::lda(status_cat ~ Improve_demo,
                       data = change,
                       na.action = "na.omit",
                       CV=TRUE)
demog_comps <- get_lda_comps(change, dfa_demog)

consists <- c(combo = combo_comps$consistency,
              both = both_comps$consistency,
              threat = thr_comps$consistency,
              demography = demog_comps$consistency)

consists

The LDA model using threat and demography scores as separate predictors has the highest consistency. To see the (mis-)classifications, with FWS's assignment in the columns and the assignment expected given the scores in rows, we have:

t(both_comps$marg_tab)

The species FWS has not recommended for a status change are the species most likely to have threat or demographic change scores that indicate a change is happening. We can plot the cross-classifications to get another view:

lda_class$Model <- paste0("Model ", lda_class$Model)
lda_class$Pct_cross <- paste0(lda_class$Pct_cross)

LDA_plot2 <- ggplot(data=lda_class, aes(x=classify, y=N_cross, fill=factor(col))) +
             geom_bar(stat="identity") +
             geom_text(y=26, aes(label=Pct_cross), size=3, colour="gray40") +
             scale_fill_manual(values=c('#0A4783', 'darkolivegreen4', '#f49831'),
                               labels=c("overprotect", "consistent", "underprotect"),
                               name="") +
             facet_grid(Model ~ FWS_rec) +
             ylim(0, 27) +
             labs(x = "\nLDA classification",
                  y = "# species classified\nper FWS x LDA category\n",
                  title = "FWS status change recommendation") +
             theme_bw() +
             theme(text=element_text(size=12, family="Open Sans"),
                   panel.background=element_rect(fill="white"))
LDA_plot2

Defenders-ESC/threatdemog documentation built on May 6, 2019, 1:58 p.m.