In jacob-ogre/section7.quality: Section 7 Quality: Analyses for Evansen et al.

Data import

library(AICcmodavg)
library(dplyr)
library(ggplot2)
library(ggthemes)
library(glmulti)
library(knitr)
library(lm.beta)
library(lmtest)
library(MASS)
library(nnet)
library(ordinal)
library(readxl)
library(readr)
library(section7.quality)
library(viridis)

# Load the data
data(combo, package = "section7.quality")
data(formal, package = "section7.quality")
data(informal, package = "section7.quality")
data(no_prog, package = "section7.quality")
data(all_no_prog, package = "section7.quality")

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)
  if (is.null(layout)) {
    layout <- matrix(
      seq(1, cols * ceiling(numPlots/cols)),
      ncol = cols, nrow = ceiling(numPlots/cols)
    )
  }

  if (numPlots==1) {
    print(plots[[1]])
  } else {
    grid.newpage()
    pushViewport(
      viewport(
        layout = grid.layout(
          nrow(layout), 
          ncol(layout)
        )
      )
    )

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      matchidx <- as.data.frame(
        which(layout == i, arr.ind = TRUE)
      )
      print(
        plots[[i]], 
        vp = viewport(
          layout.pos.row = matchidx$row,
          layout.pos.col = matchidx$col
        )
      )
    }
  }
}

Analyses

We start broad and slowly drill down on the results:

1. overall quality for formal + informal
2. overall specificity of conservation actions for formals
3. components of specificity and quality for formals
4. components of informal quality

Overall quality

All consultations

all_overQ <- run_glm_combo(combo)
kable(aictab(all_overQ$mods))

bar <- ggplot(combo, aes(combo$Overall_Q, fill = Service)) + 
  geom_histogram(
    alpha = 0.8, 
    position = position_dodge(), 
    bins = 7) + 
  labs(x = "Overall Quality") +
  scale_fill_viridis(discrete = TRUE) +
  theme_hc() +
  theme(
    legend.position = c(.05, .95),
    legend.justification = c("left", "top"),
    legend.box.just = "right",
    legend.margin = margin(6, 6, 6, 6)
  )

box <- ggplot(combo, aes(y = combo$Overall_Q, x = Service)) +
  geom_boxplot(aes(fill = combo$Service), alpha = 0.8) + 
  labs(x = "Overall Quality") +
  scale_color_viridis(discrete = TRUE) +
  scale_fill_viridis(discrete = TRUE) +
  coord_flip() +
  labs(x = "", y = "") +
  theme_hc() + 
  theme(legend.position = "none",
        axis.text.x = element_blank(),
        axis.ticks = element_blank(),
        plot.margin = unit(c(2.5,0.3,0.1,1.1), "cm"))

multiplot(box, bar)

# pdf(file = "inst/figures/Fig1a.pdf", height = 7, width = 7)
# multiplot(box, bar)
# dev.off()

Models 9 (best) & 2 differ only by year...base inference on 9, but note year effect estimates. Now to calculate the confidence intervals for the odds ratios:

mod9_OR_CI <- get_ORs_CIs(all_overQ$mods[[9]])
kable(mod9_OR_CI)

Formal and informal consultations separate

formal_overQ <- run_glm_combo_formal(formal)
kable(aictab(formal_overQ$mods))
mod2_OR_CI <- get_ORs_CIs(formal_overQ$mods[[2]])
kable(mod2_OR_CI)

informal_overQ <- run_glm_combo_informal(informal)
kable(aictab(informal_overQ$mods))
mod2_OR_CI <- get_ORs_CIs(informal_overQ$mods[[2]])
kable(mod2_OR_CI)

#######################
# plot!
f_bar <- ggplot(formal, aes(formal$Overall_Q, fill = Service)) + 
  geom_histogram(
    alpha = 0.8, 
    position = position_dodge(), 
    bins = 7) + 
  labs(x = "Overall Quality", y = "# consultations") +
  scale_fill_viridis(discrete = TRUE) +
  theme_hc() +
  theme(
    legend.position = c(.05, .95),
    legend.justification = c("left", "top"),
    legend.box.just = "right",
    legend.margin = margin(6, 6, 6, 6)
  )

i_bar <- ggplot(informal, aes(informal$Overall_Q, fill = Service)) + 
  geom_histogram(
    alpha = 0.8, 
    position = position_dodge(), 
    bins = 7) + 
  labs(x = "Overall Quality") +
  scale_fill_viridis(discrete = TRUE) +
  theme_hc() +
  theme(
    axis.title.y = element_blank(),
    legend.position = "none",
    legend.justification = c("left", "top"),
    legend.box.just = "right",
    legend.margin = margin(6, 6, 6, 6)
  )

f_box <- ggplot(formal, aes(y = formal$Overall_Q, x = Service)) +
  geom_boxplot(aes(fill = formal$Service), alpha = 0.8) + 
  labs(x = "Overall Quality") +
  scale_color_viridis(discrete = TRUE) +
  scale_fill_viridis(discrete = TRUE) +
  coord_flip() +
  ggtitle("Formal") +
  labs(x = "", y = "") +
  theme_hc() + 
  theme(legend.position = "none",
        axis.text.x = element_blank(),
        axis.ticks = element_blank(),
        plot.margin = unit(c(2.5,0.2,0.1,0.1), "cm"),
        plot.title = element_text(hjust = 0.5))

i_box <- ggplot(informal, aes(y = informal$Overall_Q, x = Service)) +
  geom_boxplot(aes(fill = informal$Service), alpha = 0.8) + 
  labs(x = "Overall Quality") +
  scale_color_viridis(discrete = TRUE) +
  scale_fill_viridis(discrete = TRUE) +
  coord_flip() +
  ggtitle("Informal") +
  labs(x = "", y = "") +
  theme_hc() + 
  theme(legend.position = "none",
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.margin = unit(c(2.5,0.3,0.1,0.5), "cm"),
        plot.title = element_text(hjust = 0.5))

multiplot(f_box, f_bar, i_box, i_bar, cols = 2)

pdf(file = "inst/figures/Fig1b.pdf", height = 7, width = 7)
multiplot(f_box, f_bar, i_box, i_bar, cols = 2)
dev.off()

Quality components: formal consultations

formal$yr <- formal$Year - 2008
fws <- filter(formal, Service == "FWS")
fws_prog <- filter(formal, Service == "FWS" & Programmatic == 1)
fws_noprog <- filter(formal, Service == "FWS" & Programmatic == 0)

We next dive into the individuals consultation components to identify patterns of variation. We focus on just two predictors, Service and year in the models, which include a random effects term for the consultation. We exclude whether the consultation was programmatic because programmatics were only done by FWS, which confounds parameter estimation. Instead, we provide only a simple graphical analyses of programmatic effects.

Species status

statusQ_mods <- run_polrs(formal, formal$status_Q)
aictab(statusQ_mods$mods)
statusQ_mods$summaries[[2]] 
statusQ_mods$coef[[2]]
get_ORs_CIs_alt(statusQ_mods$mods[[2]])

sp_stat <- ggplot(formal, aes(status_Q)) +
  geom_bar(
    aes(fill = Service), 
    position = position_dodge(),
    alpha = 0.8) +
  scale_fill_viridis(discrete = TRUE) +
  ggtitle("Species Status") +
  labs(x = "", y = "# consultations") +
  theme_hc() + 
  theme(legend.position = c(0.2, 0.8),
        plot.title = element_text(hjust = 0.5))

ggplot(fws, aes(status_Q, fill = factor(Programmatic))) +
  geom_histogram(bins = 6,
                 position = "identity",
                 alpha = 0.3) +
  scale_fill_viridis(discrete = TRUE) +
  labs(x = "Species Status Score",
       y = "# of consultations",
       fill = "Programmatic?") +
  theme_hc() + 
  theme(legend.position = "top")

tapply(fws$status_Q, INDEX = fws$Programmatic, FUN = mean, na.rm = TRUE)

The three models all have very similar AICc values, but models 1 (Service) and 3 (year) are subsets of model 2. Even though its AICc is ~2 units higher, we use model 2 b/c all params are present and very similar. The main result is that neither Service nor year have strong signals, nor (per the plot above and the means between yes/no) does whether FWS consultations were programmatic.

Baseline conditions

baselineQ_mods <- run_polrs(formal, formal$baseline_Q)
aictab(baselineQ_mods$mods)
baselineQ_mods$summaries[[1]]
baselineQ_mods$coef[[1]]
get_ORs_CIs_alt(baselineQ_mods$mods[[1]])
tapply(formal$baseline_Q, INDEX = formal$Service, FUN = mean, na.rm = TRUE)

base_q <- ggplot(formal, aes(baseline_Q)) +
  geom_bar(
    aes(fill = Service), 
    position = position_dodge(preserve = "single"),           
    alpha = 0.8) +
  scale_fill_viridis(discrete = TRUE) +
  ggtitle("Env. Baseline") +
  labs(x = "Score", y = "# consultations") +
  theme_hc() + 
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5))

ggplot(fws, aes(baseline_Q, fill = factor(Programmatic))) +
  geom_histogram(bins = 4,
                 position = "identity",
                 alpha = 0.3) +
  scale_fill_viridis(discrete = TRUE) +
  labs(x = "Env. Baseline Score",
       y = "# of consultations",
       fill = "Programmatic?") +
  theme_hc() + 
  theme(legend.position = "top")

tapply(fws$baseline_Q, INDEX = fws$Programmatic, FUN = mean, na.rm = TRUE)

Unlike the status quality, Service stands out as the best predictor of the quality of the environmental baseline. Similarly, the figure and means indicate that the programmatic consultations were important to increasing the env. baseline score for FWS.

effects analysis quality

effQ_mods <- run_polrs(formal, formal$eff_Q)
aictab(effQ_mods$mods)
effQ_mods$summaries[[2]]
effQ_mods$coef[[2]]
get_ORs_CIs_alt(effQ_mods$mods[[2]])
tapply(formal$eff_Q, INDEX = formal$Service, FUN = mean, na.rm = TRUE)

eff_q <- ggplot(formal, aes(formal$eff_Q)) +
  geom_bar(
    aes(fill = formal$Service), 
    width = 0.9, 
    position = position_dodge(preserve = "single"),
    alpha = 0.8) +
  scale_fill_viridis(discrete = TRUE) +
  ggtitle("Effects Analysis") +
  labs(x = "", y = "") +
  theme_hc() + 
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5))

ggplot(fws, aes(eff_Q, fill = factor(Programmatic))) +
  geom_histogram(bins = 3,
                 position = "identity",
                 alpha = 0.3) +
  scale_fill_viridis(discrete = TRUE) +
  labs(x = "Effects Score",
       y = "# of consultations",
       fill = "Programmatic?") +
  theme_hc() + 
  theme(legend.position = "top")

tapply(fws$eff_Q, INDEX = fws$Programmatic, FUN = mean, na.rm = TRUE)

As with the status of the species section, the quality of the effects analysis did not differ between Services nor, within FWS, between (non-) programmatics.

cumulative effects

CEQ_mods <- run_polrs(formal, formal$CE_Q)
aictab(CEQ_mods$mods)
CEQ_mods$summaries[[1]]
CEQ_mods$coef[[1]]
get_ORs_CIs_alt(CEQ_mods$mods[[1]])
tapply(formal$CE_Q, INDEX = formal$Service, FUN = mean, na.rm = TRUE)

ce_q <- ggplot(formal, aes(CE_Q)) +
  geom_bar(
    aes(fill = Service), 
    position = position_dodge(preserve = "single"),
    alpha = 0.8) +
  scale_fill_viridis(discrete = TRUE) +
  ggtitle("Cum. Effects") +
  labs(x = "Score", y = "") +
  theme_hc() + 
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5))

ggplot(fws, aes(CE_Q, fill = factor(Programmatic))) +
  geom_histogram(bins = 3,
                 position = "identity",
                 alpha = 0.8) +
  scale_fill_viridis(discrete = TRUE) +
  labs(x = "Cum. Effects Score",
       y = "# of consultations",
       fill = "Programmatic?") +
  theme_hc() + 
  theme(legend.position = "top")

tapply(fws$CE_Q, INDEX = fws$Programmatic, FUN = mean, na.rm = TRUE)

As with the status of the species section, the quality of cumulative effects analysis did not differ between Services. There was a marginal difference - based on simple graphical analysis - between programmatics and non-programmatics.

pdf(file = "inst/figures/Fig2.pdf", height = 10, width = 10)
multiplot(sp_stat, base_q, eff_q, ce_q, cols = 2)
dev.off()

Conclusion The primary cause of the difference between FWS and NMFS overall scores was driven by differences in the environmental baseline analysis. NMFS's scores tended to be higher than FWS in the other quality components, and even though not statistically significant alone, also like contributed to the overall difference.

Quality components: informal consultations

Action mention

action_mods <- inf_binom(informal, informal$mention_action)
aictab(action_mods$mods)
action_mods$summaries[[1]] # nothing significant (or even close)

infmac <- filter(informal, !is.na(mention_action))
mact <- ggplot(infmac, aes(factor(mention_action))) +
  geom_histogram(
    aes(fill = Service), 
    position = position_dodge(preserve = "single"),
    stat = "count",
    alpha = 0.8) +
  scale_fill_viridis(discrete = TRUE) +
  ggtitle("Mention Action") +
  labs(x = "", y = "# consultations") +
  scale_x_discrete(labels = c(0, 1)) +
  theme_hc() + 
  theme(legend.position = c(0.1, 0.8),
        plot.title = element_text(hjust = 0.5))

There is no relationship between whether the action was mentioned and Service or consultation duration.

Action analysis

act_analysis_mods <- inf_binom(informal, informal$act_analysis)
aictab(act_analysis_mods$mods)
act_analysis_mods$summaries[[1]] 
get_ORs_CIs_alt(act_analysis_mods$mods[[1]])

infan <- filter(informal, !is.na(act_analysis))
analy <- ggplot(infan, aes(factor(act_analysis))) +
  geom_histogram(
    aes(fill = Service), 
    position = position_dodge(preserve = "single"),
    stat = "count",
    alpha = 0.8) +
  scale_fill_viridis(discrete = TRUE) +
  ggtitle("Action Analysis") +
  labs(x = "Score", y = "# consultations") +
  scale_x_discrete(labels = c(0, 1)) +
  theme_hc() + 
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5))

Only duration significant - a very slight decline with duration - and tiny effect size (1.02 OR).

Species analysis

sp_analysis_mods <- inf_binom(informal, informal$sp_analysis)
aictab(sp_analysis_mods$mods)
sp_analysis_mods$summaries[[2]]
get_ORs_CIs_alt(sp_analysis_mods$mods[[2]])

infsp <- filter(informal, !is.na(sp_analysis))
sp_an <- ggplot(infsp, aes(factor(sp_analysis))) +
  geom_histogram(aes(fill = Service), 
                 # position = "identity",
                 position = position_dodge(prespreserve = "single"),
                 stat = "count",
                 alpha = 0.8) +
  scale_fill_viridis(discrete = TRUE) +
  ggtitle("Species Analysis") +
  labs(x = "", y = "") +
  scale_x_discrete(labels = c(0, 1)) +
  theme_hc() + 
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5))

Duration X Service is significant, but the effect size is again tiny (OR = 1.004).

Reason informal

mention_reason_inf_mods <- inf_binom(informal, informal$mention_reason_inf)
aictab(mention_reason_inf_mods$mods)
mention_reason_inf_mods$summaries[[1]]
get_ORs_CIs_alt(mention_reason_inf_mods$mods[[1]])

infre <- filter(informal, !is.na(mention_reason_inf))
reason <- ggplot(infre, aes(factor(mention_reason_inf))) +
  geom_bar(
    aes(fill = Service), 
    position = position_dodge(preserve = "single"),
    stat = "count",
    alpha = 0.8) +
  scale_fill_viridis(discrete = TRUE) +
  ggtitle("Reason Informal") +
  labs(x = "Score", y = "") +
  scale_x_discrete(labels = c(0, 1)) +
  theme_hc() + 
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5))

Only duration significant - a very slight increase of whether the reason a consultation was informal with duration - and tiny effect size (1.03 OR).

Map included?

map_mods <- inf_binom(informal, informal$map)
aictab(map_mods$mods)
map_mods$summaries[[1]] # Service has OR = 5.08
get_ORs_CIs_alt(mention_reason_inf_mods$mods[[1]])

infmap <- filter(informal, !is.na(map))
wmap <- ggplot(infmap, aes(factor(map))) +
  geom_bar(
    aes(fill = Service), 
    bins = 2,
    stat = "count",
    position = position_dodge(prespreserve = "single"),
    alpha = 0.8) +
  scale_fill_viridis(discrete = TRUE) +
  ggtitle("Map?") +
  labs(x = "Score", y = "") +
  theme_hc() + 
  scale_x_discrete(labels = c(0, 1)) +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5))

pdf(file = "inst/figures/Fig3.pdf", height = 8, width = 12)
multiplot(mact, analy, sp_an, reason, wmap, cols = 3)
dev.off()

Additional Plots

This plots section is basically obsolete because figures have been added above.

Overall quality terms

Formal + informal

par(mfrow = c(1, 3))
boxplot(combo$Overall_Q ~ combo$Programmatic,
        main = "Programmatics",
        ylab = "% overall quality")
boxplot(combo$Overall_Q ~ combo$Service,
        main = "With programmatics",
        ylab = "")
boxplot(combo[combo$Programmatic == "0", ]$Overall_Q ~ combo[combo$Programmatic == "0", ]$Service,
        main = "Without programmatics",
        ylab = "")
par(mfrow = c(1, 1))

par(mfrow = c(2, 1))
ggplot(combo, aes(x = Service, y = Overall_Q)) +
  geom_boxplot() +
  coord_flip() +
  theme_hc()

Formal only

par(mfrow = c(1, 3))
boxplot(formal$Overall_Q ~ formal$Programmatic,
        main = "Programmatics",
        ylab = "% overall quality")
boxplot(formal$Overall_Q ~ formal$Service,
        main = "With programmatics",
        ylab = "")
boxplot(formal[formal$Programmatic == "0", ]$Overall_Q ~ formal[formal$Programmatic == "0", ]$Service,
        main = "Without programmatics",
        ylab = "")
par(mfrow = c(1, 1))

Year effects

par(mfrow = c(3, 1))
boxplot(combo$Overall_Q ~ factor(combo$Year),
        main = "All through time",
        ylab = "% overall quality")
boxplot(formal$Overall_Q ~ factor(formal$Year),
        main = "Formal through time",
        ylab = "% overall quality")
boxplot(informal$Overall_Q ~ factor(informal$Year),
        main = "Informals through time",
        ylab = "% overall quality")
par(mfrow = c(1, 1))

Service effect

boxplot(informal$Overall_Q ~ informal$Service,
        main = "Informals By Service",
        ylab = "% overall quality")

Duration effects

plot(combo$Overall_Q ~ combo$total_duration,
     xlab = "Duration (days)",
     ylab = "% overall quality",
     pch = c(17, 19)[as.numeric(factor(combo$Formal))],
     col = c("red", "blue")[as.numeric(factor(combo$Service))])
abline(lm(Overall_Q ~ total_duration, data = combo[combo$Service == "NMFS",]),
       col="red")
abline(lm(Overall_Q ~ total_duration, data = combo[combo$Service == "FWS",]),
       col="blue")

And some quick plots of quality components so I could visualize the data while troubleshooting the models

plt <- ggplot(data = formal, aes(status_Q)) +
       geom_bar() +
       facet_grid(. ~ Service + Programmatic) +
       theme_hc()
plt

plt <- ggplot(data = formal, aes(baseline_Q)) +
       geom_bar() +
       facet_grid(. ~ Service + Programmatic) +
       theme_hc()
plt

plt <- ggplot(data = formal, aes(eff_Q)) +
       geom_bar() +
       facet_grid(. ~ Service + Programmatic) +
       theme_hc()
plt

plt <- ggplot(data = formal, aes(CE_Q)) +
       geom_bar() +
       facet_grid(. ~ Service + Programmatic) +
       theme_hc()
plt

jacob-ogre/section7.quality documentation built on May 28, 2019, 9:56 p.m.