In KevinSee/QRFcapacity: Carrying Capacity with QRF
# load knitr for markdown
library(knitr)
knitr::opts_chunk$set(echo=FALSE,
warning=FALSE,
# error = FALSE,
message=FALSE)
#options(tinytex.verbose = TRUE)
options(knitr.kable.NA = '-')
options(knitr.table.format = "pandoc")
library(kableExtra)
# load needed packages
library(tidyverse)
library(sf)
library(QRFcapacity)
library(janitor)
library(ggpubr)
library(ggrepel)
library(scales)
library(magrittr)
library(here)
theme_set(theme_bw())
Introduction
We can validate Chinook summer parr estimates of capacity with spawner-recruit data from a number of watersheds. Morgan Bond at the Northwest Fisheries Science Center helped compile this data set with help from a variety of agencies and organizations including WDFW, ODFW, IDFG, USFS, BioAnalysts, Biomark, Nez Perce, etc. Most of the parr (recruits) data comes from rotary screw traps (with the exception of the Chiwawa, where snorkel surveys of parr were used to estimate parr directly). Estimates of parr-to-presmolt and parr-to-smolt survival were used to transform estimates of fall and spring emigrants to estimates of equivalent summer parr. Spawner estimates came from either redd counts or direct estimates of spawners.
From this dataset, we fit several spawner-recruit curves to each location, including Beverton-Holt, Ricker and hockey stick. The data at some locations led to some curves not fitting (no model convergence). Each of these spawner-recruit curves also generated an estimate of summer parr carrying capacity, which we can then compare to QRF estimates. We want to emphasize that the data for the spawner-recruit curves and the QRF estimates are completely independent, and therefore provide a good way to validate the QRF estimates of carrying capacity.
# read in trap polygons
trap_poly = read_sf(here('data/raw/spawn_rec/trap_regions.shp')) %>%
st_transform(crs = 5070) %>%
select(-Id)
# which QRF model to use?
mod_choice = c('juv_summer',
'juv_summer_dash',
'redds')[1]
load(here(paste0('output/modelFits/qrf_', mod_choice, '.rda')))
# define max capacity based on observed densities
cap_max = qrf_mod_df %>%
group_by(Species) %>%
summarise(quant95 = quantile(fish_dens, 0.95),
quant98 = quantile(fish_dens, 0.98),
quant99 = quantile(fish_dens, 0.99),
quant995 = quantile(fish_dens, 0.995),
quant100 = quantile(fish_dens, 1),
.groups = "drop") %>%
filter(Species == 'Chinook') %>%
pull(quant99)
# read in capacity estimates and create a shapefile
load(here(paste0('output/modelFits/extrap_200rch_', mod_choice, '.rda')))
in_covar_range_df = model_svy_df$pred_all_rchs[[1]] %>%
select(UniqueID, in_covar_range)
rm(mod_data_weights, model_svy_df, extrap_covars)
data("rch_200")
# shrink predictions over cap_max down to cap_max, make SE NA
all_preds %<>%
mutate(chnk_per_m_se = if_else(chnk_per_m > cap_max,
as.numeric(NA), chnk_per_m_se),
chnk_per_m = if_else(chnk_per_m > cap_max,
cap_max, chnk_per_m))
rch_200_cap = rch_200 %>%
filter(chnk) %>%
select(UniqueID, GNIS_Name, reach_leng:HUC8_code,
chnk, chnk_use, chnk_ESU_DPS:chnk_NWR_NAME) %>%
left_join(all_preds) %>%
mutate(chnk_tot = chnk_per_m * reach_leng,
chnk_tot_se = chnk_per_m_se * reach_leng) %>%
left_join(in_covar_range_df) %>%
# filter(in_covar_range) %>%
st_transform(st_crs(trap_poly)) %>%
select(everything(), geometry)
# rch_200_cap %>%
# st_drop_geometry() %>%
# xtabs(~ (chnk_per_m > cap_max) + in_covar_range, .)
rm(rch_200, all_preds)
qrf_caps = trap_poly %>%
st_join(rch_200_cap) %>%
mutate_at(vars(chnk_ESU_DPS, chnk_MPG, chnk_NWR_NAME),
list(fct_explicit_na)) %>%
group_by(chnk_ESU_DPS, chnk_MPG, trap_name) %>%
summarise(n_rch = n_distinct(UniqueID),
tot_lgth = sum(reach_leng, na.rm = T),
cap = sum(chnk_tot, na.rm = T),
cap_se = mean(chnk_per_m_se, na.rm = T) * tot_lgth) %>%
# cap_se = sqrt(sum(chnk_tot_se^2))) %>%
ungroup() %>%
mutate(cap_cv = cap_se / cap) %>%
rename(Trap = trap_name) %>%
mutate(Trap = recode(Trap,
'Minam R.' = 'Minam River',
'Entiat R.' = 'Entiat River',
'Crooked Fork Creek' = 'Crooked Fork',
'Upper Grande Ronde R.' ='Upper Grande Ronde River',
'Toucannon R.' = 'Tucannon River',
'Pahsimeroi' = 'Pahsimeroi River',
'Hood R.' = 'Hood River',
'John Day Middle Fork' = 'Middle Fork John Day',
'Upper Salmon River' = 'Upper Salmon River Chinook',
'John Day Upper Mainstem' = 'Upper Mainstem John Day',
'East Fork Salmon' = 'East Fork Salmon River',
'Lostine R.' = 'Lostine River'))
load(here(paste0('output/modelFits/extrap_200rch_RF_', mod_choice, '.rda')))
in_covar_range_df = model_rf_df$pred_all_rchs[[1]] %>%
select(UniqueID, in_covar_range)
rm(mod_data_weights, model_rf_df, extrap_covars)
rch_200_cap = st_read(here("output/gpkg/StreamNetwork_Capacity_Chnk.gpkg")) %>%
st_transform(st_crs(trap_poly)) %>%
filter(spp_domain) %>%
left_join(in_covar_range_df) %>%
mutate(sum_juv_tot = sum_juv_per_m * reach_leng,
sum_juv_tot_se = sum_juv_per_m_se * reach_leng)
qrf_caps = trap_poly %>%
st_join(rch_200_cap2) %>%
rename(chnk_ESU_DPS = ESU_DPS,
chnk_MPG = MPG,
chnk_NWR_NAME = NWR_NAME) %>%
mutate(across(c(chnk_ESU_DPS, chnk_MPG, chnk_NWR_NAME),
fct_explicit_na)) %>%
group_by(chnk_ESU_DPS, chnk_MPG, trap_name) %>%
summarise(n_rch = n_distinct(UniqueID),
tot_lgth = sum(reach_leng, na.rm = T),
cap = sum(sum_juv_tot, na.rm = T),
cap_se = mean(sum_juv_per_m_se, na.rm = T) * tot_lgth,
# cap_se = sqrt(sum(sum_juv_tot_se^2)),
.groups = "drop") %>%
mutate(cap_cv = cap_se / cap) %>%
rename(Trap = trap_name) %>%
mutate(Trap = recode(Trap,
'Minam R.' = 'Minam River',
'Entiat R.' = 'Entiat River',
'Crooked Fork Creek' = 'Crooked Fork',
'Upper Grande Ronde R.' ='Upper Grande Ronde River',
'Toucannon R.' = 'Tucannon River',
'Pahsimeroi' = 'Pahsimeroi River',
'Hood R.' = 'Hood River',
'John Day Middle Fork' = 'Middle Fork John Day',
'Upper Salmon River' = 'Upper Salmon River Chinook',
'John Day Upper Mainstem' = 'Upper Mainstem John Day',
'East Fork Salmon' = 'East Fork Salmon River',
'Lostine R.' = 'Lostine River'))
# read in capacity estimates and create a shapefile
load(here(paste0('output/modelFits/extrap_200rch_RF_', mod_choice, '.rda')))
in_covar_range_df = model_rf_df$pred_all_rchs[[1]] %>%
select(UniqueID, in_covar_range)
rm(mod_data_weights, model_rf_df, extrap_covars)
data("rch_200")
rch_200_cap = rch_200 %>%
filter(chnk) %>%
select(UniqueID, GNIS_Name, reach_leng:HUC8_code,
chnk, chnk_use, chnk_ESU_DPS:chnk_NWR_NAME) %>%
left_join(all_preds) %>%
mutate(chnk_tot = chnk_per_m * reach_leng,
chnk_tot_se = chnk_per_m_se * reach_leng) %>%
left_join(in_covar_range_df) %>%
# filter(in_covar_range) %>%
st_transform(st_crs(trap_poly)) %>%
select(everything(), geometry)
# rch_200_cap %>%
# st_drop_geometry() %>%
# xtabs(~ (chnk_per_m > cap_max) + in_covar_range, .)
#
# rch_200_cap %<>%
# mutate(chnk_per_m_se = if_else(chnk_per_m > cap_max,
# as.numeric(NA), chnk_per_m_se),
# chnk_per_m = if_else(chnk_per_m > cap_max,
# cap_max, chnk_per_m))
rm(rch_200, all_preds)
qrf_caps = trap_poly %>%
st_join(rch_200_cap) %>%
mutate_at(vars(chnk_ESU_DPS, chnk_MPG, chnk_NWR_NAME),
list(fct_explicit_na)) %>%
group_by(chnk_ESU_DPS, chnk_MPG, trap_name) %>%
summarise(n_rch = n_distinct(UniqueID),
tot_lgth = sum(reach_leng, na.rm = T),
cap = sum(chnk_tot, na.rm = T),
cap_se = mean(chnk_per_m_se, na.rm = T) * tot_lgth) %>%
# cap_se = sqrt(sum(chnk_tot_se^2))) %>%
ungroup() %>%
mutate(cap_cv = cap_se / cap) %>%
rename(Trap = trap_name) %>%
mutate(Trap = recode(Trap,
'Minam R.' = 'Minam River',
'Entiat R.' = 'Entiat River',
'Crooked Fork Creek' = 'Crooked Fork',
'Upper Grande Ronde R.' ='Upper Grande Ronde River',
'Toucannon R.' = 'Tucannon River',
'Pahsimeroi' = 'Pahsimeroi River',
'Hood R.' = 'Hood River',
'John Day Middle Fork' = 'Middle Fork John Day',
'Upper Salmon River' = 'Upper Salmon River Chinook',
'John Day Upper Mainstem' = 'Upper Mainstem John Day',
'East Fork Salmon' = 'East Fork Salmon River',
'Lostine R.' = 'Lostine River'))
# read in capacity estimates and create a shapefile
load(here(paste0('output/modelFits/extrap_mastPts_', mod_choice, '.rda')))
data("gaa")
data("chnk_domain")
# shrink predictions over cap_max down to cap_max, make SE NA
# all_preds %>%
# xtabs(~ chnk_per_m > cap_max, .)
all_preds %<>%
mutate(chnk_per_m_se = if_else(chnk_per_m > cap_max,
as.numeric(NA), chnk_per_m_se),
chnk_per_m = if_else(chnk_per_m > cap_max,
cap_max, chnk_per_m))
capacity_pts = all_preds %>%
left_join(gaa %>%
select(Site, Lon, Lat)) %>%
filter(!is.na(Lon)) %>%
st_as_sf(coords = c('Lon', 'Lat'),
crs = 4326) %>%
st_transform(st_crs(trap_poly))
qrf_caps2 = trap_poly %>%
st_join(capacity_pts) %>%
group_by(trap_name) %>%
summarise(n_pts = n_distinct(Site),
avg_cap = mean(chnk_per_m, na.rm = T),
avg_cap_se = mean(chnk_per_m_se, na.rm = T)) %>%
ungroup() %>%
left_join(trap_poly %>%
st_join(chnk_domain %>%
mutate(reach_leng = st_length(.))) %>%
mutate_at(vars(chnk_ESU_DPS = ESU_DPS,
chnk_MPG = MPG),
list(fct_explicit_na)) %>%
group_by(chnk_ESU_DPS, chnk_MPG, trap_name) %>%
summarise(tot_lgth = sum(reach_leng)) %>%
ungroup() %>%
st_drop_geometry()) %>%
group_by(chnk_ESU_DPS, chnk_MPG, trap_name, n_pts, tot_lgth) %>%
summarise(cap = avg_cap * as.numeric(tot_lgth),
cap_se = avg_cap_se * as.numeric(tot_lgth),
cap_cv = cap_se / cap) %>%
ungroup() %>%
rename(Trap = trap_name) %>%
mutate(Trap = recode(Trap,
'Minam R.' = 'Minam River',
'Entiat R.' = 'Entiat River',
'Crooked Fork Creek' = 'Crooked Fork',
'Upper Grande Ronde R.' ='Upper Grande Ronde River',
'Toucannon R.' = 'Tucannon River',
'Pahsimeroi' = 'Pahsimeroi River',
'Hood R.' = 'Hood River',
'John Day Middle Fork' = 'Middle Fork John Day',
'Upper Salmon River' = 'Upper Salmon River Chinook',
'John Day Upper Mainstem' = 'Upper Mainstem John Day',
'East Fork Salmon' = 'East Fork Salmon River',
'Lostine R.' = 'Lostine River')) %>%
select(starts_with('chnk'),
Trap,
everything(),
geometry) %>%
filter(n_pts > 1)
trp_rch = rch_200_cap %>%
st_join(trap_poly) %>%
filter(!is.na(trap_name))
trp_pts = capacity_pts %>%
st_join(trap_poly) %>%
filter(!is.na(trap_name)) %>%
st_join(chnk_domain %>%
st_buffer(dist = 500)) %>%
filter(!is.na(Species)) %>%
filter(trap_name %in% unique(trp_rch$trap_name))
comp_df = trp_pts %>%
group_by(Trap = trap_name) %>%
nest() %>%
rename(pt_data = data) %>%
left_join(trp_rch %>%
group_by(Trap = trap_name) %>%
nest() %>%
rename(rch_data = data)) %>%
left_join(chnk_domain %>%
st_join(trap_poly) %>%
filter(!is.na(trap_name)) %>%
mutate(reach_leng = st_length(.)) %>%
mutate_at(vars(reach_leng),
list(as.numeric)) %>%
group_by(Trap = trap_name) %>%
nest() %>%
rename(old_domain = data)) %>%
mutate(fig = map2(pt_data,
rch_data,
function(x, y) {
ggplot() +
geom_sf(data = x,
aes(color = chnk_per_m),
alpha = 0.6,
size = 3) +
geom_sf(data = y,
aes(color = chnk_per_m)) +
scale_color_viridis_c(direction = -1,
name = 'Chnk / m') +
theme(axis.text = element_blank(),
axis.title = element_blank())
}),
fig = map2(fig,
Trap,
function(x, y) {
x +
labs(title = y)
})) %>%
arrange(Trap)
all_p = ggpubr::ggarrange(plotlist = comp_df$fig,
ncol = 2,
nrow = 2)
# ggsave(here('output/figures/capacity_maps_all.pdf'),
# all_p,
# width = 8,
# height = 11)
# pdf(here('output/figures/capacity_maps_all.pdf'),
# width = 8,
# height = 11)
# for(i in 1:length(all_p)) print(all_p)
# dev.off()
# comp_df %>%
# mutate(rch_cap = map_dbl(rch_data,
# .f = function(x) {
# x %>%
# summarise_at(vars(chnk_tot),
# list(sum)) %>%
# pull(chnk_tot)
# })) %>%
# mutate(pts_cap = map2_dbl(pt_data,
# old_domain,
# .f = function(x, y) {
# avg_cap = x %>%
# summarise_at(vars(chnk_per_m),
# list(mean),
# na.rm = T) %>%
# pull(chnk_per_m)
# strm_lgth = y %>%
# summarise_at(vars(reach_leng),
# list(sum),
# na.rm = T) %>%
# pull(reach_leng)
# avg_cap * strm_lgth
# }))
# create a comparison dataframe
trp_comp = qrf_caps %>%
as_tibble() %>%
select(chnk_ESU_DPS:cap_se) %>%
mutate(method = "reach") %>%
bind_rows(qrf_caps2 %>%
as_tibble() %>%
select(chnk_ESU_DPS:cap_se) %>%
mutate(method = "point") %>%
mutate(tot_lgth = as.numeric(tot_lgth))) %>%
pivot_wider(id_cols = c(chnk_ESU_DPS:Trap),
names_from = "method",
values_from = c("tot_lgth",
"n_rch", "n_pts",
"cap",
"cap_se")) %>%
rename(n_pts = n_pts_point,
n_rch = n_rch_reach) %>%
select(-n_rch_point, -n_pts_reach) %>%
filter(!is.na(cap_reach),
!is.na(cap_point))
# recalculate point capacity, using reach lengths
trp_comp %<>%
mutate(cap_point = cap_point / tot_lgth_point * tot_lgth_reach)
# drop one outlier prediction
trp_comp %<>%
filter(cap_reach < 1e10)
# mark some SE as NA if too big
trp_comp %<>%
mutate_at(vars(starts_with("cap_se")),
list(~ if_else(. > 1e10,
NA_real_,
.)))
lng_p = trp_comp %>%
ggplot(aes(x = tot_lgth_point/1000,
y = tot_lgth_reach/1000,
color = chnk_MPG)) +
geom_abline(linetype = 2) +
geom_point(size = 5) +
geom_text_repel(aes(label = Trap)) +
theme(legend.position = "bottom") +
labs(x = "Master Sample Points",
y = "200 m Reaches",
color = 'MPG',
title = "Length of Chinook Domain (km)")
lng_rel_p = trp_comp %>%
mutate(lng_diff = tot_lgth_reach - tot_lgth_point,
lng_rel_diff = lng_diff / tot_lgth_point) %>%
ggplot(aes(x = Trap,
y = lng_rel_diff,
color = chnk_MPG)) +
geom_hline(yintercept = 0,
linetype = 2,
color = 'red') +
geom_point(size = 5) +
theme(axis.text.x = element_text(angle = 45,
hjust = 1),
legend.position = "bottom") +
labs(y = "Relative Difference",
color = 'MPG',
title = "Length of Chinook Domain (km)")
cap_p = trp_comp %>%
mutate_at(vars(starts_with("cap")),
list(~ . / 1000)) %>%
ggplot(aes(x = cap_point,
y = cap_reach,
color = chnk_MPG)) +
geom_abline(linetype = 2) +
geom_errorbar(aes(ymin = cap_reach + cap_se_reach * qnorm(0.025),
ymax = cap_reach + cap_se_reach * qnorm(0.975)),
width = 0) +
geom_errorbarh(aes(xmin = cap_point + cap_se_point * qnorm(0.025),
xmax = cap_point + cap_se_point * qnorm(0.975)),
height = 0) +
geom_point(size = 5) +
geom_text_repel(aes(label = Trap)) +
scale_x_continuous(trans = "log",
breaks = breaks_pretty(n = 8)) +
scale_y_continuous(trans = "log",
breaks = breaks_pretty(n = 8)) +
theme(legend.position = "bottom") +
labs(x = "Master Sample Points",
y = "200 m Reaches",
color = 'MPG',
title = "Capacity - Chinook Parr (x 1,000)")
trp_comp %>%
mutate(cap_diff = cap_reach - cap_point,
cap_rel_diff = cap_diff / cap_point,
cap_diff_se = sqrt(cap_se_reach^2 + cap_se_point^2),
incl_zero = if_else(cap_diff + qnorm(0.025) * cap_diff_se < 0 &
cap_diff + qnorm(0.975) * cap_diff_se > 0,
T, F)) %>%
ggplot(aes(x = Trap,
y = cap_diff,
color = chnk_MPG)) +
geom_hline(yintercept = 0,
linetype = 2,
color = 'red') +
geom_errorbar(aes(ymin = cap_diff + cap_diff_se * qnorm(0.025),
ymax = cap_diff + cap_diff_se * qnorm(0.975)),
width = 0,
position = position_dodge(width = 0.3)) +
geom_point(position = position_dodge(width = 0.3)) +
theme(axis.text.x = element_text(angle = 45,
hjust = 1),
legend.position = "bottom") +
labs(x = 'Trap',
y = 'Reach - Point',
color = "MPG")
trp_comp %>%
mutate(cap_diff = cap_reach - cap_point,
cap_rel_diff = cap_diff / cap_point,
cap_diff_se = sqrt(cap_se_reach^2 + cap_se_point^2),
incl_zero = if_else(cap_diff + qnorm(0.025) * cap_diff_se < 0 &
cap_diff + qnorm(0.975) * cap_diff_se > 0,
T, F)) %>%
ggplot(aes(x = Trap,
y = cap_rel_diff,
color = chnk_MPG)) +
geom_hline(yintercept = 0,
linetype = 2,
color = 'red') +
# geom_errorbar(aes(ymin = cap_diff + cap_diff_se * qnorm(0.025),
# ymax = cap_diff + cap_diff_se * qnorm(0.975)),
# width = 0,
# position = position_dodge(width = 0.3)) +
geom_point(position = position_dodge(width = 0.3)) +
theme(axis.text.x = element_text(angle = 45,
hjust = 1),
legend.position = "bottom") +
labs(x = 'Trap',
y = '(Reach - Point) / Point',
color = "MPG")
cap_dens_p = trp_comp %>%
mutate(cap_dens_point = cap_point/tot_lgth_point,
cap_dens_reach = cap_reach/tot_lgth_reach) %>%
# mutate(cap_dens_point = cap_point / tot_lgth_reach) %>%
ggplot(aes(x = cap_dens_point,
y = cap_dens_reach,
color = chnk_MPG)) +
geom_abline(linetype = 2) +
geom_point() +
geom_text_repel(aes(label = Trap)) +
theme(legend.position = "bottom") +
labs(x = "Master Sample Points",
y = "200 m Reaches",
color = 'MPG',
title = "Average Capacity Density (fish / m)")
trp_comp %>%
mutate(cap_dens_point = cap_point/tot_lgth_point,
cap_dens_reach = cap_reach/tot_lgth_reach) %>%
mutate(dens_diff = cap_dens_reach - cap_dens_point,
dens_rel_diff = dens_diff / cap_dens_point) %>%
ggplot(aes(x = Trap,
y = dens_rel_diff,
color = chnk_MPG)) +
geom_hline(yintercept = 0,
linetype = 2,
color = 'red') +
geom_point(size = 5) +
theme(axis.text.x = element_text(angle = 45,
hjust = 1),
legend.position = "bottom") +
labs(y = "Relative Difference",
color = 'MPG',
title = "Average Capacity Density (fish / m)")
# comp_df$rch_data[[which(comp_df$Trap == 'Upper Lemhi')]] %>%
# st_drop_geometry() %>%
# group_by(GNIS_Name) %>%
# summarise(n_rch = n_distinct(UniqueID),
# avg_per_m = mean(chnk_per_m),
# sd_per_m = sd(chnk_per_m, na.rm = T),
# avg_cap = mean(chnk_tot),
# sd_cap = sd(chnk_tot))
lem_strm_cap <- comp_df$pt_data[[which(comp_df$Trap == 'Upper Lemhi')]] %>%
group_by(StreamName) %>%
summarise_at(vars(chnk_per_m),
list(mean),
na.rm = T) %>%
st_drop_geometry()
lem_leng <- comp_df$old_domain[[which(comp_df$Trap == 'Upper Lemhi')]] %>%
group_by(StreamName) %>%
summarise_at(vars(reach_leng),
list(sum)) %>%
st_drop_geometry()
add_uplem_cap = lem_strm_cap %>%
full_join(lem_leng) %>%
filter(!is.na(reach_leng)) %>%
mutate(cap = chnk_per_m * reach_leng) %>%
filter(StreamName != "Lemhi River") %>%
summarise_at(vars(cap),
list(sum)) %>%
pull(cap)
qrf_caps$cap[qrf_caps$Trap == 'Upper Lemhi'] = qrf_caps$cap[qrf_caps$Trap == 'Upper Lemhi'] + add_uplem_cap
library(FSA)
bh1 <- srFuns('BevertonHolt',
param = 3)
data("spawn_recr_data")
# fit Beverton-Holt curve with capacity fixed at QRF capacity prediction
qrf_params = spawn_recr_data %>%
group_by(Population) %>%
nest() %>%
left_join(qrf_caps %>%
as_tibble() %>%
select(Trap, cap, cap_se),
by = c('Population' = 'Trap')) %>%
filter(!is.na(cap)) %>%
mutate(form = 'QRF',
inits = map2(.x = data,
.y = cap,
.f = function(x, y) {
inits = srStarts(Parr ~ Spawners,
data = x,
type = 'BevertonHolt',
param = 3)
inits$a = if_else(inits$a < 0,
2e-3,
inits$a)
inits$b = 1 / y
return(inits)
}),
mod_fit = map2(.x = data,
.y = inits,
.f = function(x, y) {
fit = suppressWarnings(try(nls(log(Parr) ~ log(bh1(Spawners, a, b)),
data = x,
start = y,
algorithm = 'port',
lower = c(0, y['b']),
upper = c(Inf, y['b']))))
}),
coefs = map(mod_fit,
.f = function(x) {
if(class(x) == 'try-error') {
return(as.numeric(NA))
} else coef(x)
}),
preds = map2(.x = data,
.y = mod_fit,
.f = function(x, y) {
if(class(y) != 'try-error') {
tibble(Spawners = 1:max(x$Spawners)) %>%
mutate(Parr = predict(y,
newdata = .),
Parr = exp(Parr))
} else return(tibble(Spawners = NA,
Parr = NA))
})) %>%
rename(est = cap,
se = cap_se) %>%
# select(-c(area:tot_length, tot_cap_cv)) %>%
mutate(cv = se / est)
data("spawn_rec_params")
spawn_rec_params %<>%
bind_rows(qrf_params)
# pull out estimates of capacity
cap_tbl = spawn_rec_params %>%
select(Population, form, est, se, cv) %>%
arrange(Population, form) %>%
mutate(form = fct_relevel(form,
"QRF",
after = Inf))
# make a few estimates NA (too big)
cap_tbl %<>%
mutate(across(c(est, se, cv),
~ if_else(Population %in% c('Entiat River',
'Marsh Creek',
'Minam River',
'Tucannon River') &
form == 'BevertonHolt',
NA_real_, .)))
# pull out fitted spawner recruit curves
curve_fits = spawn_rec_params %>%
select(Population, form, preds) %>%
unnest(cols = "preds") %>%
arrange(Population, form) %>%
mutate(form = fct_relevel(form,
"QRF",
after = Inf))
# define polygons for uncertainty shading
plot_max = Inf
plot_max = 5e5
cap_poly = cap_tbl %>%
mutate(lwrCI = est + se * qnorm(0.025),
uprCI = est + se * qnorm(0.975)) %>%
mutate(lwrCI = if_else(lwrCI < 0, 0, lwrCI),
uprCI = if_else(uprCI > plot_max, plot_max, uprCI)) %>%
left_join(spawn_rec_params %>%
select(Population, data) %>%
unnest(cols = "data") %>%
group_by(Population) %>%
summarise_at(vars(Spawners, Parr),
tibble::lst(min, max),
na.rm = T)) %>%
group_by(Population, form) %>%
summarise(coords = list(tibble(x = c(0, 0, Spawners_max, Spawners_max),
y = c(lwrCI, uprCI, uprCI, lwrCI)))) %>%
ungroup() %>%
unnest(cols = c(coords)) %>%
mutate(form = fct_relevel(form,
"QRF",
after = Inf))
# which watersheds to plot?
keep_wtsds = qrf_caps %>%
filter(!is.na(cap)) %>%
filter(!Trap %in% c('Methow R.')) %>%
pull(Trap)
# keep_wtsds = c('Catherine Creek',
# 'Chiwawa R.',
# 'Lostine River',
# 'Marsh Creek',
# 'Minam River',
# 'Pahsimeroi River',
# 'Upper Grande Ronde River',
# 'Upper Lemhi',
# 'Upper Mainstem John Day')
spawn_rec_comp_p = spawn_rec_params %>%
# filter out a couple populations with super poor fits
filter(Population %in% keep_wtsds) %>%
select(Population, data) %>%
unnest(cols = "data") %>%
distinct() %>%
ggplot(aes(x = Spawners,
y = Parr)) +
geom_polygon(data = cap_poly %>%
# filter out a couple populations with super poor fits
filter(Population %in% keep_wtsds),
aes(x = x,
y = y,
fill = form),
alpha = 0.2) +
geom_line(data = curve_fits %>%
# filter out a couple populations with super poor fits
filter(Population %in% keep_wtsds),
aes(color = form),
lwd = 1.5) +
geom_hline(data = cap_tbl %>%
# filter out a couple populations with super poor fits
filter(Population %in% keep_wtsds),
aes(yintercept = est,
color = form),
linetype = 2) +
geom_point(aes(shape = Spawner_type),
size = 3) +
facet_wrap(~ Population,
ncol = 3,
scales = 'free') +
scale_color_brewer(palette = 'Set1',
breaks = c('QRF',
'BevertonHolt',
'Ricker',
'Hockey'),
labels = c('QRF',
'Beverton Holt',
'Ricker',
'Hockey Stick'),
name = 'Model') +
scale_fill_brewer(palette = 'Set1',
breaks = c('QRF',
'BevertonHolt',
'Ricker',
'Hockey'),
labels = c('QRF',
'Beverton Holt',
'Ricker',
'Hockey Stick'),
name = 'Model') +
labs(shape = 'Spawner\nType') +
theme(legend.position = "bottom")
dodge_width = 0.5
sr_est_comp_p = cap_tbl %>%
filter(Population %in% keep_wtsds) %>%
mutate(lwrCI = est + qnorm(0.025) * se,
lwrCI = if_else(lwrCI <= 0, 1, lwrCI),
uprCI = est + qnorm(0.975) * se,
uprCI = if_else(uprCI > 1.3e6, 1.3e6, uprCI)) %>%
mutate(type = if_else(form == 'QRF',
'QRF', 'S-R')) %>%
ggplot(aes(x = Population,
y = est,
color = form,
group = form)) +
geom_errorbar(aes(ymin = lwrCI,
ymax = uprCI),
position = position_dodge(width = dodge_width),
width = 0) +
geom_point(aes(shape = type),
fill = 'white',
size = 2,
position = position_dodge(width = dodge_width)) +
scale_shape_manual(values = c('QRF' = 19,
'S-R' = 21),
name = 'Type') +
scale_color_brewer(palette = 'Set1',
breaks = c('QRF',
'BevertonHolt',
'Ricker',
'Hockey'),
labels = c('QRF',
'Beverton Holt',
'Ricker',
'Hockey Stick'),
name = 'Model') +
# scale_y_continuous(trans = 'log') +
# theme(axis.text.x = element_text(angle = 90),
# axis.text.y = element_blank()) +
theme(axis.text = element_blank(),
legend.position = "bottom") +
facet_wrap(~ Population,
scales = 'free') +
labs(y = 'Capacity Estimate',
x = NULL)
# make some figures for presentations
keep_wtsds = c('Catherine Creek',
'Chiwawa R.',
"Entiat River",
"Hayden Creek",
'Lostine River',
'Marsh Creek',
'Minam River',
"Pahsimeroi River",
'Upper Grande Ronde River',
# 'Upper Lemhi',
"Upper Mainstem John Day")
spawn_rec_comp_p2 = spawn_rec_params %>%
# filter out a couple populations with super poor fits
filter(Population %in% keep_wtsds) %>%
select(Population, data) %>%
unnest(cols = "data") %>%
distinct() %>%
ggplot(aes(x = Spawners,
y = Parr)) +
# geom_polygon(data = cap_poly %>%
# # filter out a couple populations with super poor fits
# filter(Population %in% keep_wtsds),
# aes(x = x,
# y = y,
# fill = form),
# alpha = 0.2) +
geom_line(data = curve_fits %>%
# filter out a couple populations with super poor fits
filter(Population %in% keep_wtsds),
aes(color = form),
lwd = 1.5) +
geom_hline(data = cap_tbl %>%
# filter out a couple populations with super poor fits
filter(Population %in% keep_wtsds),
aes(yintercept = est,
color = form),
linetype = 2) +
geom_point(aes(shape = Spawner_type),
size = 3) +
facet_wrap(~ Population,
# ncol = 3,
nrow = 2,
scales = 'free') +
scale_color_brewer(palette = 'Set1',
breaks = c('QRF',
'BevertonHolt',
'Ricker',
'Hockey'),
labels = c('QRF',
'Beverton Holt',
'Ricker',
'Hockey Stick'),
name = 'Model') +
scale_fill_brewer(palette = 'Set1',
breaks = c('QRF',
'BevertonHolt',
'Ricker',
'Hockey'),
labels = c('QRF',
'Beverton Holt',
'Ricker',
'Hockey Stick'),
name = 'Model') +
labs(shape = 'Spawner\nType') +
theme(legend.position = "bottom",
axis.text = element_text(size = 6))
dodge_width = 0.5
sr_est_comp_p2 = cap_tbl %>%
filter(Population %in% keep_wtsds) %>%
mutate(lwrCI = est + qnorm(0.025) * se,
lwrCI = if_else(lwrCI <= 0, 1, lwrCI),
uprCI = est + qnorm(0.975) * se,
uprCI = if_else(uprCI > 1.3e6, 1.3e6, uprCI)) %>%
mutate(type = if_else(form == 'QRF',
'QRF', 'S-R')) %>%
ggplot(aes(x = Population,
y = est,
color = form,
group = form)) +
geom_errorbar(aes(ymin = lwrCI,
ymax = uprCI),
position = position_dodge(width = dodge_width),
width = 0) +
geom_point(aes(shape = type),
fill = 'white',
size = 2,
position = position_dodge(width = dodge_width)) +
scale_shape_manual(values = c('QRF' = 19,
'S-R' = 21),
name = 'Type') +
scale_color_brewer(palette = 'Set1',
breaks = c('QRF',
'BevertonHolt',
'Ricker',
'Hockey'),
labels = c('QRF',
'Beverton Holt',
'Ricker',
'Hockey Stick'),
name = 'Model') +
# scale_y_continuous(trans = 'log') +
# theme(axis.text.x = element_text(angle = 90),
# axis.text.y = element_blank()) +
theme(axis.text = element_blank(),
legend.position = "bottom") +
facet_wrap(~ Population,
nrow = 2,
scales = 'free') +
labs(y = 'Capacity Estimate',
x = NULL)
ggsave(here("output/figures/spawn_rec_curves.pdf"),
spawn_rec_comp_p2,
width = 11.5,
height = 5)
ggsave(here("output/figures/spawn_rec_comp.pdf"),
sr_est_comp_p2,
width = 11.5,
height = 5)
sr_df = spawn_rec_params %>%
filter(Population == "Chiwawa R.") %>%
select(Population, data) %>%
unnest(cols = "data") %>%
distinct() %>%
mutate(across(Parr,
~ . / 1e5))
sr_curve = curve_fits %>%
mutate(across(Parr,
~ . / 1e5)) %>%
filter(Population == "Chiwawa R.") %>%
mutate(form = fct_relevel(form,
"QRF",
after = Inf))
sr_cap = cap_tbl %>%
mutate(across(c(est, se),
~ . / 1e5)) %>%
filter(Population == "Chiwawa R.") %>%
mutate(form = fct_relevel(form,
"QRF",
after = Inf))
sr_poly = cap_poly %>%
mutate(across(c(y),
~ . / 1e5)) %>%
filter(Population == "Chiwawa R.") %>%
mutate(form = fct_relevel(form,
"QRF",
after = Inf))
chiw_p1 = sr_df %>%
ggplot(aes(x = Spawners,
y = Parr)) +
geom_point(size = 3) +
labs(y = "Parr (x 100,000)") +
scale_y_continuous(limits = c(0, 5))
chiw_p2 = chiw_p1 +
geom_line(data = sr_curve %>%
filter(form != "QRF"),
aes(color = form),
lwd = 2) +
geom_hline(data = sr_cap %>%
filter(form != "QRF"),
aes(yintercept = est,
color = form),
linetype = 2,
lwd = 1.5) +
scale_color_brewer(palette = 'Set1',
breaks = c('QRF',
'BevertonHolt',
'Ricker',
'Hockey'),
labels = c('QRF',
'Beverton Holt',
'Ricker',
'Hockey Stick'),
name = 'Model') +
theme(legend.position = "bottom")
chiw_p3 = chiw_p1 +
geom_line(data = sr_curve,
aes(color = form),
lwd = 2) +
geom_hline(data = sr_cap,
aes(yintercept = est,
color = form),
linetype = 2,
lwd = 1.5) +
scale_color_brewer(palette = 'Set1',
breaks = c('QRF',
'BevertonHolt',
'Ricker',
'Hockey'),
labels = c('QRF',
'Beverton Holt',
'Ricker',
'Hockey Stick'),
name = 'Model') +
theme(legend.position = "bottom")
chiw_p4 = chiw_p3 +
geom_polygon(data = sr_poly,
aes(x = x,
y = y,
fill = form),
alpha = 0.2) +
scale_fill_brewer(palette = 'Set1',
breaks = c('QRF',
'BevertonHolt',
'Ricker',
'Hockey'),
labels = c('QRF',
'Beverton Holt',
'Ricker',
'Hockey Stick'),
name = 'Model')
chiw_p_list = list(chiw_p1,
chiw_p2,
chiw_p3,
chiw_p4)
for(i in 1:length(chiw_p_list)) {
ggsave(here("output/figures",
paste0("Chiwawa_SR_p", i, ".pdf")),
chiw_p_list[[i]],
width = 7,
height = 6)
}
cap_tbl_ci = cap_tbl %>%
mutate(lwrCI = est + qnorm(0.025) * se,
lwrCI = if_else(lwrCI <= 0, 1, lwrCI),
uprCI = est + qnorm(0.975) * se,
uprCI = if_else(uprCI > 1.3e6, 1.3e6, uprCI))
comp_df = cap_tbl_ci %>%
filter(form != 'QRF') %>%
filter(!is.na(est)) %>%
inner_join(cap_tbl_ci %>%
filter(form == 'QRF') %>%
select(Population, form, est) %>%
spread(form, est) %>%
rename(qrf_est = QRF) %>%
left_join(cap_tbl_ci %>%
filter(form == 'QRF') %>%
select(Population, form, lwrCI) %>%
spread(form, lwrCI) %>%
rename(qrf_lwrCI = QRF)) %>%
left_join(cap_tbl_ci %>%
filter(form == 'QRF') %>%
select(Population, form, uprCI) %>%
spread(form, uprCI) %>%
rename(qrf_uprCI = QRF))) %>%
mutate(in_qrf_ci = if_else(est <= qrf_uprCI & est >= qrf_lwrCI, T, F),
qrf_in_ci = if_else(qrf_est <= uprCI & qrf_est >= lwrCI, T, F),
overlap_ci = if_else(lwrCI <= qrf_uprCI & uprCI >= qrf_lwrCI, T, F))
comp_df %>%
group_by(form) %>%
summarise(n_pops = n_distinct(Population)) %>%
full_join(comp_df %>%
group_by(form) %>%
summarise_at(vars(ends_with("_ci")),
list(sum),
na.rm = T)) %>%
mutate_at(vars(perc_in_qrf_ci = in_qrf_ci,
perc_qrf_in_ci = qrf_in_ci,
perc_overlap_ci = overlap_ci),
list(~ . / n_pops))
comp_df %>%
filter(!overlap_ci)
We then plotted the data, as well as the fitted curves, including the one which uses QRF estimates of capacity as one of the parameters, on the same plot, to compare the various estimates.
spawn_rec_comp_p
sr_est_comp_p
KevinSee/QRFcapacity documentation built on Feb. 27, 2023, 3:57 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# load knitr for markdown library(knitr) knitr::opts_chunk$set(echo=FALSE, warning=FALSE, # error = FALSE, message=FALSE) #options(tinytex.verbose = TRUE) options(knitr.kable.NA = '-') options(knitr.table.format = "pandoc") library(kableExtra)
# load needed packages library(tidyverse) library(sf) library(QRFcapacity) library(janitor) library(ggpubr) library(ggrepel) library(scales) library(magrittr) library(here) theme_set(theme_bw())
Introduction
We can validate Chinook summer parr estimates of capacity with spawner-recruit data from a number of watersheds. Morgan Bond at the Northwest Fisheries Science Center helped compile this data set with help from a variety of agencies and organizations including WDFW, ODFW, IDFG, USFS, BioAnalysts, Biomark, Nez Perce, etc. Most of the parr (recruits) data comes from rotary screw traps (with the exception of the Chiwawa, where snorkel surveys of parr were used to estimate parr directly). Estimates of parr-to-presmolt and parr-to-smolt survival were used to transform estimates of fall and spring emigrants to estimates of equivalent summer parr. Spawner estimates came from either redd counts or direct estimates of spawners.
From this dataset, we fit several spawner-recruit curves to each location, including Beverton-Holt, Ricker and hockey stick. The data at some locations led to some curves not fitting (no model convergence). Each of these spawner-recruit curves also generated an estimate of summer parr carrying capacity, which we can then compare to QRF estimates. We want to emphasize that the data for the spawner-recruit curves and the QRF estimates are completely independent, and therefore provide a good way to validate the QRF estimates of carrying capacity.
# read in trap polygons trap_poly = read_sf(here('data/raw/spawn_rec/trap_regions.shp')) %>% st_transform(crs = 5070) %>% select(-Id) # which QRF model to use? mod_choice = c('juv_summer', 'juv_summer_dash', 'redds')[1]
load(here(paste0('output/modelFits/qrf_', mod_choice, '.rda'))) # define max capacity based on observed densities cap_max = qrf_mod_df %>% group_by(Species) %>% summarise(quant95 = quantile(fish_dens, 0.95), quant98 = quantile(fish_dens, 0.98), quant99 = quantile(fish_dens, 0.99), quant995 = quantile(fish_dens, 0.995), quant100 = quantile(fish_dens, 1), .groups = "drop") %>% filter(Species == 'Chinook') %>% pull(quant99)
# read in capacity estimates and create a shapefile load(here(paste0('output/modelFits/extrap_200rch_', mod_choice, '.rda'))) in_covar_range_df = model_svy_df$pred_all_rchs[[1]] %>% select(UniqueID, in_covar_range) rm(mod_data_weights, model_svy_df, extrap_covars) data("rch_200") # shrink predictions over cap_max down to cap_max, make SE NA all_preds %<>% mutate(chnk_per_m_se = if_else(chnk_per_m > cap_max, as.numeric(NA), chnk_per_m_se), chnk_per_m = if_else(chnk_per_m > cap_max, cap_max, chnk_per_m)) rch_200_cap = rch_200 %>% filter(chnk) %>% select(UniqueID, GNIS_Name, reach_leng:HUC8_code, chnk, chnk_use, chnk_ESU_DPS:chnk_NWR_NAME) %>% left_join(all_preds) %>% mutate(chnk_tot = chnk_per_m * reach_leng, chnk_tot_se = chnk_per_m_se * reach_leng) %>% left_join(in_covar_range_df) %>% # filter(in_covar_range) %>% st_transform(st_crs(trap_poly)) %>% select(everything(), geometry) # rch_200_cap %>% # st_drop_geometry() %>% # xtabs(~ (chnk_per_m > cap_max) + in_covar_range, .) rm(rch_200, all_preds) qrf_caps = trap_poly %>% st_join(rch_200_cap) %>% mutate_at(vars(chnk_ESU_DPS, chnk_MPG, chnk_NWR_NAME), list(fct_explicit_na)) %>% group_by(chnk_ESU_DPS, chnk_MPG, trap_name) %>% summarise(n_rch = n_distinct(UniqueID), tot_lgth = sum(reach_leng, na.rm = T), cap = sum(chnk_tot, na.rm = T), cap_se = mean(chnk_per_m_se, na.rm = T) * tot_lgth) %>% # cap_se = sqrt(sum(chnk_tot_se^2))) %>% ungroup() %>% mutate(cap_cv = cap_se / cap) %>% rename(Trap = trap_name) %>% mutate(Trap = recode(Trap, 'Minam R.' = 'Minam River', 'Entiat R.' = 'Entiat River', 'Crooked Fork Creek' = 'Crooked Fork', 'Upper Grande Ronde R.' ='Upper Grande Ronde River', 'Toucannon R.' = 'Tucannon River', 'Pahsimeroi' = 'Pahsimeroi River', 'Hood R.' = 'Hood River', 'John Day Middle Fork' = 'Middle Fork John Day', 'Upper Salmon River' = 'Upper Salmon River Chinook', 'John Day Upper Mainstem' = 'Upper Mainstem John Day', 'East Fork Salmon' = 'East Fork Salmon River', 'Lostine R.' = 'Lostine River'))
load(here(paste0('output/modelFits/extrap_200rch_RF_', mod_choice, '.rda'))) in_covar_range_df = model_rf_df$pred_all_rchs[[1]] %>% select(UniqueID, in_covar_range) rm(mod_data_weights, model_rf_df, extrap_covars) rch_200_cap = st_read(here("output/gpkg/StreamNetwork_Capacity_Chnk.gpkg")) %>% st_transform(st_crs(trap_poly)) %>% filter(spp_domain) %>% left_join(in_covar_range_df) %>% mutate(sum_juv_tot = sum_juv_per_m * reach_leng, sum_juv_tot_se = sum_juv_per_m_se * reach_leng) qrf_caps = trap_poly %>% st_join(rch_200_cap2) %>% rename(chnk_ESU_DPS = ESU_DPS, chnk_MPG = MPG, chnk_NWR_NAME = NWR_NAME) %>% mutate(across(c(chnk_ESU_DPS, chnk_MPG, chnk_NWR_NAME), fct_explicit_na)) %>% group_by(chnk_ESU_DPS, chnk_MPG, trap_name) %>% summarise(n_rch = n_distinct(UniqueID), tot_lgth = sum(reach_leng, na.rm = T), cap = sum(sum_juv_tot, na.rm = T), cap_se = mean(sum_juv_per_m_se, na.rm = T) * tot_lgth, # cap_se = sqrt(sum(sum_juv_tot_se^2)), .groups = "drop") %>% mutate(cap_cv = cap_se / cap) %>% rename(Trap = trap_name) %>% mutate(Trap = recode(Trap, 'Minam R.' = 'Minam River', 'Entiat R.' = 'Entiat River', 'Crooked Fork Creek' = 'Crooked Fork', 'Upper Grande Ronde R.' ='Upper Grande Ronde River', 'Toucannon R.' = 'Tucannon River', 'Pahsimeroi' = 'Pahsimeroi River', 'Hood R.' = 'Hood River', 'John Day Middle Fork' = 'Middle Fork John Day', 'Upper Salmon River' = 'Upper Salmon River Chinook', 'John Day Upper Mainstem' = 'Upper Mainstem John Day', 'East Fork Salmon' = 'East Fork Salmon River', 'Lostine R.' = 'Lostine River'))
# read in capacity estimates and create a shapefile load(here(paste0('output/modelFits/extrap_200rch_RF_', mod_choice, '.rda'))) in_covar_range_df = model_rf_df$pred_all_rchs[[1]] %>% select(UniqueID, in_covar_range) rm(mod_data_weights, model_rf_df, extrap_covars) data("rch_200") rch_200_cap = rch_200 %>% filter(chnk) %>% select(UniqueID, GNIS_Name, reach_leng:HUC8_code, chnk, chnk_use, chnk_ESU_DPS:chnk_NWR_NAME) %>% left_join(all_preds) %>% mutate(chnk_tot = chnk_per_m * reach_leng, chnk_tot_se = chnk_per_m_se * reach_leng) %>% left_join(in_covar_range_df) %>% # filter(in_covar_range) %>% st_transform(st_crs(trap_poly)) %>% select(everything(), geometry) # rch_200_cap %>% # st_drop_geometry() %>% # xtabs(~ (chnk_per_m > cap_max) + in_covar_range, .) # # rch_200_cap %<>% # mutate(chnk_per_m_se = if_else(chnk_per_m > cap_max, # as.numeric(NA), chnk_per_m_se), # chnk_per_m = if_else(chnk_per_m > cap_max, # cap_max, chnk_per_m)) rm(rch_200, all_preds) qrf_caps = trap_poly %>% st_join(rch_200_cap) %>% mutate_at(vars(chnk_ESU_DPS, chnk_MPG, chnk_NWR_NAME), list(fct_explicit_na)) %>% group_by(chnk_ESU_DPS, chnk_MPG, trap_name) %>% summarise(n_rch = n_distinct(UniqueID), tot_lgth = sum(reach_leng, na.rm = T), cap = sum(chnk_tot, na.rm = T), cap_se = mean(chnk_per_m_se, na.rm = T) * tot_lgth) %>% # cap_se = sqrt(sum(chnk_tot_se^2))) %>% ungroup() %>% mutate(cap_cv = cap_se / cap) %>% rename(Trap = trap_name) %>% mutate(Trap = recode(Trap, 'Minam R.' = 'Minam River', 'Entiat R.' = 'Entiat River', 'Crooked Fork Creek' = 'Crooked Fork', 'Upper Grande Ronde R.' ='Upper Grande Ronde River', 'Toucannon R.' = 'Tucannon River', 'Pahsimeroi' = 'Pahsimeroi River', 'Hood R.' = 'Hood River', 'John Day Middle Fork' = 'Middle Fork John Day', 'Upper Salmon River' = 'Upper Salmon River Chinook', 'John Day Upper Mainstem' = 'Upper Mainstem John Day', 'East Fork Salmon' = 'East Fork Salmon River', 'Lostine R.' = 'Lostine River'))
# read in capacity estimates and create a shapefile load(here(paste0('output/modelFits/extrap_mastPts_', mod_choice, '.rda'))) data("gaa") data("chnk_domain") # shrink predictions over cap_max down to cap_max, make SE NA # all_preds %>% # xtabs(~ chnk_per_m > cap_max, .) all_preds %<>% mutate(chnk_per_m_se = if_else(chnk_per_m > cap_max, as.numeric(NA), chnk_per_m_se), chnk_per_m = if_else(chnk_per_m > cap_max, cap_max, chnk_per_m)) capacity_pts = all_preds %>% left_join(gaa %>% select(Site, Lon, Lat)) %>% filter(!is.na(Lon)) %>% st_as_sf(coords = c('Lon', 'Lat'), crs = 4326) %>% st_transform(st_crs(trap_poly)) qrf_caps2 = trap_poly %>% st_join(capacity_pts) %>% group_by(trap_name) %>% summarise(n_pts = n_distinct(Site), avg_cap = mean(chnk_per_m, na.rm = T), avg_cap_se = mean(chnk_per_m_se, na.rm = T)) %>% ungroup() %>% left_join(trap_poly %>% st_join(chnk_domain %>% mutate(reach_leng = st_length(.))) %>% mutate_at(vars(chnk_ESU_DPS = ESU_DPS, chnk_MPG = MPG), list(fct_explicit_na)) %>% group_by(chnk_ESU_DPS, chnk_MPG, trap_name) %>% summarise(tot_lgth = sum(reach_leng)) %>% ungroup() %>% st_drop_geometry()) %>% group_by(chnk_ESU_DPS, chnk_MPG, trap_name, n_pts, tot_lgth) %>% summarise(cap = avg_cap * as.numeric(tot_lgth), cap_se = avg_cap_se * as.numeric(tot_lgth), cap_cv = cap_se / cap) %>% ungroup() %>% rename(Trap = trap_name) %>% mutate(Trap = recode(Trap, 'Minam R.' = 'Minam River', 'Entiat R.' = 'Entiat River', 'Crooked Fork Creek' = 'Crooked Fork', 'Upper Grande Ronde R.' ='Upper Grande Ronde River', 'Toucannon R.' = 'Tucannon River', 'Pahsimeroi' = 'Pahsimeroi River', 'Hood R.' = 'Hood River', 'John Day Middle Fork' = 'Middle Fork John Day', 'Upper Salmon River' = 'Upper Salmon River Chinook', 'John Day Upper Mainstem' = 'Upper Mainstem John Day', 'East Fork Salmon' = 'East Fork Salmon River', 'Lostine R.' = 'Lostine River')) %>% select(starts_with('chnk'), Trap, everything(), geometry) %>% filter(n_pts > 1)
trp_rch = rch_200_cap %>% st_join(trap_poly) %>% filter(!is.na(trap_name)) trp_pts = capacity_pts %>% st_join(trap_poly) %>% filter(!is.na(trap_name)) %>% st_join(chnk_domain %>% st_buffer(dist = 500)) %>% filter(!is.na(Species)) %>% filter(trap_name %in% unique(trp_rch$trap_name)) comp_df = trp_pts %>% group_by(Trap = trap_name) %>% nest() %>% rename(pt_data = data) %>% left_join(trp_rch %>% group_by(Trap = trap_name) %>% nest() %>% rename(rch_data = data)) %>% left_join(chnk_domain %>% st_join(trap_poly) %>% filter(!is.na(trap_name)) %>% mutate(reach_leng = st_length(.)) %>% mutate_at(vars(reach_leng), list(as.numeric)) %>% group_by(Trap = trap_name) %>% nest() %>% rename(old_domain = data)) %>% mutate(fig = map2(pt_data, rch_data, function(x, y) { ggplot() + geom_sf(data = x, aes(color = chnk_per_m), alpha = 0.6, size = 3) + geom_sf(data = y, aes(color = chnk_per_m)) + scale_color_viridis_c(direction = -1, name = 'Chnk / m') + theme(axis.text = element_blank(), axis.title = element_blank()) }), fig = map2(fig, Trap, function(x, y) { x + labs(title = y) })) %>% arrange(Trap) all_p = ggpubr::ggarrange(plotlist = comp_df$fig, ncol = 2, nrow = 2) # ggsave(here('output/figures/capacity_maps_all.pdf'), # all_p, # width = 8, # height = 11) # pdf(here('output/figures/capacity_maps_all.pdf'), # width = 8, # height = 11) # for(i in 1:length(all_p)) print(all_p) # dev.off() # comp_df %>% # mutate(rch_cap = map_dbl(rch_data, # .f = function(x) { # x %>% # summarise_at(vars(chnk_tot), # list(sum)) %>% # pull(chnk_tot) # })) %>% # mutate(pts_cap = map2_dbl(pt_data, # old_domain, # .f = function(x, y) { # avg_cap = x %>% # summarise_at(vars(chnk_per_m), # list(mean), # na.rm = T) %>% # pull(chnk_per_m) # strm_lgth = y %>% # summarise_at(vars(reach_leng), # list(sum), # na.rm = T) %>% # pull(reach_leng) # avg_cap * strm_lgth # }))
# create a comparison dataframe trp_comp = qrf_caps %>% as_tibble() %>% select(chnk_ESU_DPS:cap_se) %>% mutate(method = "reach") %>% bind_rows(qrf_caps2 %>% as_tibble() %>% select(chnk_ESU_DPS:cap_se) %>% mutate(method = "point") %>% mutate(tot_lgth = as.numeric(tot_lgth))) %>% pivot_wider(id_cols = c(chnk_ESU_DPS:Trap), names_from = "method", values_from = c("tot_lgth", "n_rch", "n_pts", "cap", "cap_se")) %>% rename(n_pts = n_pts_point, n_rch = n_rch_reach) %>% select(-n_rch_point, -n_pts_reach) %>% filter(!is.na(cap_reach), !is.na(cap_point)) # recalculate point capacity, using reach lengths trp_comp %<>% mutate(cap_point = cap_point / tot_lgth_point * tot_lgth_reach) # drop one outlier prediction trp_comp %<>% filter(cap_reach < 1e10) # mark some SE as NA if too big trp_comp %<>% mutate_at(vars(starts_with("cap_se")), list(~ if_else(. > 1e10, NA_real_, .))) lng_p = trp_comp %>% ggplot(aes(x = tot_lgth_point/1000, y = tot_lgth_reach/1000, color = chnk_MPG)) + geom_abline(linetype = 2) + geom_point(size = 5) + geom_text_repel(aes(label = Trap)) + theme(legend.position = "bottom") + labs(x = "Master Sample Points", y = "200 m Reaches", color = 'MPG', title = "Length of Chinook Domain (km)") lng_rel_p = trp_comp %>% mutate(lng_diff = tot_lgth_reach - tot_lgth_point, lng_rel_diff = lng_diff / tot_lgth_point) %>% ggplot(aes(x = Trap, y = lng_rel_diff, color = chnk_MPG)) + geom_hline(yintercept = 0, linetype = 2, color = 'red') + geom_point(size = 5) + theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "bottom") + labs(y = "Relative Difference", color = 'MPG', title = "Length of Chinook Domain (km)") cap_p = trp_comp %>% mutate_at(vars(starts_with("cap")), list(~ . / 1000)) %>% ggplot(aes(x = cap_point, y = cap_reach, color = chnk_MPG)) + geom_abline(linetype = 2) + geom_errorbar(aes(ymin = cap_reach + cap_se_reach * qnorm(0.025), ymax = cap_reach + cap_se_reach * qnorm(0.975)), width = 0) + geom_errorbarh(aes(xmin = cap_point + cap_se_point * qnorm(0.025), xmax = cap_point + cap_se_point * qnorm(0.975)), height = 0) + geom_point(size = 5) + geom_text_repel(aes(label = Trap)) + scale_x_continuous(trans = "log", breaks = breaks_pretty(n = 8)) + scale_y_continuous(trans = "log", breaks = breaks_pretty(n = 8)) + theme(legend.position = "bottom") + labs(x = "Master Sample Points", y = "200 m Reaches", color = 'MPG', title = "Capacity - Chinook Parr (x 1,000)") trp_comp %>% mutate(cap_diff = cap_reach - cap_point, cap_rel_diff = cap_diff / cap_point, cap_diff_se = sqrt(cap_se_reach^2 + cap_se_point^2), incl_zero = if_else(cap_diff + qnorm(0.025) * cap_diff_se < 0 & cap_diff + qnorm(0.975) * cap_diff_se > 0, T, F)) %>% ggplot(aes(x = Trap, y = cap_diff, color = chnk_MPG)) + geom_hline(yintercept = 0, linetype = 2, color = 'red') + geom_errorbar(aes(ymin = cap_diff + cap_diff_se * qnorm(0.025), ymax = cap_diff + cap_diff_se * qnorm(0.975)), width = 0, position = position_dodge(width = 0.3)) + geom_point(position = position_dodge(width = 0.3)) + theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "bottom") + labs(x = 'Trap', y = 'Reach - Point', color = "MPG") trp_comp %>% mutate(cap_diff = cap_reach - cap_point, cap_rel_diff = cap_diff / cap_point, cap_diff_se = sqrt(cap_se_reach^2 + cap_se_point^2), incl_zero = if_else(cap_diff + qnorm(0.025) * cap_diff_se < 0 & cap_diff + qnorm(0.975) * cap_diff_se > 0, T, F)) %>% ggplot(aes(x = Trap, y = cap_rel_diff, color = chnk_MPG)) + geom_hline(yintercept = 0, linetype = 2, color = 'red') + # geom_errorbar(aes(ymin = cap_diff + cap_diff_se * qnorm(0.025), # ymax = cap_diff + cap_diff_se * qnorm(0.975)), # width = 0, # position = position_dodge(width = 0.3)) + geom_point(position = position_dodge(width = 0.3)) + theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "bottom") + labs(x = 'Trap', y = '(Reach - Point) / Point', color = "MPG") cap_dens_p = trp_comp %>% mutate(cap_dens_point = cap_point/tot_lgth_point, cap_dens_reach = cap_reach/tot_lgth_reach) %>% # mutate(cap_dens_point = cap_point / tot_lgth_reach) %>% ggplot(aes(x = cap_dens_point, y = cap_dens_reach, color = chnk_MPG)) + geom_abline(linetype = 2) + geom_point() + geom_text_repel(aes(label = Trap)) + theme(legend.position = "bottom") + labs(x = "Master Sample Points", y = "200 m Reaches", color = 'MPG', title = "Average Capacity Density (fish / m)") trp_comp %>% mutate(cap_dens_point = cap_point/tot_lgth_point, cap_dens_reach = cap_reach/tot_lgth_reach) %>% mutate(dens_diff = cap_dens_reach - cap_dens_point, dens_rel_diff = dens_diff / cap_dens_point) %>% ggplot(aes(x = Trap, y = dens_rel_diff, color = chnk_MPG)) + geom_hline(yintercept = 0, linetype = 2, color = 'red') + geom_point(size = 5) + theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "bottom") + labs(y = "Relative Difference", color = 'MPG', title = "Average Capacity Density (fish / m)")
# comp_df$rch_data[[which(comp_df$Trap == 'Upper Lemhi')]] %>% # st_drop_geometry() %>% # group_by(GNIS_Name) %>% # summarise(n_rch = n_distinct(UniqueID), # avg_per_m = mean(chnk_per_m), # sd_per_m = sd(chnk_per_m, na.rm = T), # avg_cap = mean(chnk_tot), # sd_cap = sd(chnk_tot)) lem_strm_cap <- comp_df$pt_data[[which(comp_df$Trap == 'Upper Lemhi')]] %>% group_by(StreamName) %>% summarise_at(vars(chnk_per_m), list(mean), na.rm = T) %>% st_drop_geometry() lem_leng <- comp_df$old_domain[[which(comp_df$Trap == 'Upper Lemhi')]] %>% group_by(StreamName) %>% summarise_at(vars(reach_leng), list(sum)) %>% st_drop_geometry() add_uplem_cap = lem_strm_cap %>% full_join(lem_leng) %>% filter(!is.na(reach_leng)) %>% mutate(cap = chnk_per_m * reach_leng) %>% filter(StreamName != "Lemhi River") %>% summarise_at(vars(cap), list(sum)) %>% pull(cap) qrf_caps$cap[qrf_caps$Trap == 'Upper Lemhi'] = qrf_caps$cap[qrf_caps$Trap == 'Upper Lemhi'] + add_uplem_cap
library(FSA) bh1 <- srFuns('BevertonHolt', param = 3) data("spawn_recr_data") # fit Beverton-Holt curve with capacity fixed at QRF capacity prediction qrf_params = spawn_recr_data %>% group_by(Population) %>% nest() %>% left_join(qrf_caps %>% as_tibble() %>% select(Trap, cap, cap_se), by = c('Population' = 'Trap')) %>% filter(!is.na(cap)) %>% mutate(form = 'QRF', inits = map2(.x = data, .y = cap, .f = function(x, y) { inits = srStarts(Parr ~ Spawners, data = x, type = 'BevertonHolt', param = 3) inits$a = if_else(inits$a < 0, 2e-3, inits$a) inits$b = 1 / y return(inits) }), mod_fit = map2(.x = data, .y = inits, .f = function(x, y) { fit = suppressWarnings(try(nls(log(Parr) ~ log(bh1(Spawners, a, b)), data = x, start = y, algorithm = 'port', lower = c(0, y['b']), upper = c(Inf, y['b'])))) }), coefs = map(mod_fit, .f = function(x) { if(class(x) == 'try-error') { return(as.numeric(NA)) } else coef(x) }), preds = map2(.x = data, .y = mod_fit, .f = function(x, y) { if(class(y) != 'try-error') { tibble(Spawners = 1:max(x$Spawners)) %>% mutate(Parr = predict(y, newdata = .), Parr = exp(Parr)) } else return(tibble(Spawners = NA, Parr = NA)) })) %>% rename(est = cap, se = cap_se) %>% # select(-c(area:tot_length, tot_cap_cv)) %>% mutate(cv = se / est)
data("spawn_rec_params") spawn_rec_params %<>% bind_rows(qrf_params) # pull out estimates of capacity cap_tbl = spawn_rec_params %>% select(Population, form, est, se, cv) %>% arrange(Population, form) %>% mutate(form = fct_relevel(form, "QRF", after = Inf)) # make a few estimates NA (too big) cap_tbl %<>% mutate(across(c(est, se, cv), ~ if_else(Population %in% c('Entiat River', 'Marsh Creek', 'Minam River', 'Tucannon River') & form == 'BevertonHolt', NA_real_, .))) # pull out fitted spawner recruit curves curve_fits = spawn_rec_params %>% select(Population, form, preds) %>% unnest(cols = "preds") %>% arrange(Population, form) %>% mutate(form = fct_relevel(form, "QRF", after = Inf)) # define polygons for uncertainty shading plot_max = Inf plot_max = 5e5 cap_poly = cap_tbl %>% mutate(lwrCI = est + se * qnorm(0.025), uprCI = est + se * qnorm(0.975)) %>% mutate(lwrCI = if_else(lwrCI < 0, 0, lwrCI), uprCI = if_else(uprCI > plot_max, plot_max, uprCI)) %>% left_join(spawn_rec_params %>% select(Population, data) %>% unnest(cols = "data") %>% group_by(Population) %>% summarise_at(vars(Spawners, Parr), tibble::lst(min, max), na.rm = T)) %>% group_by(Population, form) %>% summarise(coords = list(tibble(x = c(0, 0, Spawners_max, Spawners_max), y = c(lwrCI, uprCI, uprCI, lwrCI)))) %>% ungroup() %>% unnest(cols = c(coords)) %>% mutate(form = fct_relevel(form, "QRF", after = Inf)) # which watersheds to plot? keep_wtsds = qrf_caps %>% filter(!is.na(cap)) %>% filter(!Trap %in% c('Methow R.')) %>% pull(Trap) # keep_wtsds = c('Catherine Creek', # 'Chiwawa R.', # 'Lostine River', # 'Marsh Creek', # 'Minam River', # 'Pahsimeroi River', # 'Upper Grande Ronde River', # 'Upper Lemhi', # 'Upper Mainstem John Day') spawn_rec_comp_p = spawn_rec_params %>% # filter out a couple populations with super poor fits filter(Population %in% keep_wtsds) %>% select(Population, data) %>% unnest(cols = "data") %>% distinct() %>% ggplot(aes(x = Spawners, y = Parr)) + geom_polygon(data = cap_poly %>% # filter out a couple populations with super poor fits filter(Population %in% keep_wtsds), aes(x = x, y = y, fill = form), alpha = 0.2) + geom_line(data = curve_fits %>% # filter out a couple populations with super poor fits filter(Population %in% keep_wtsds), aes(color = form), lwd = 1.5) + geom_hline(data = cap_tbl %>% # filter out a couple populations with super poor fits filter(Population %in% keep_wtsds), aes(yintercept = est, color = form), linetype = 2) + geom_point(aes(shape = Spawner_type), size = 3) + facet_wrap(~ Population, ncol = 3, scales = 'free') + scale_color_brewer(palette = 'Set1', breaks = c('QRF', 'BevertonHolt', 'Ricker', 'Hockey'), labels = c('QRF', 'Beverton Holt', 'Ricker', 'Hockey Stick'), name = 'Model') + scale_fill_brewer(palette = 'Set1', breaks = c('QRF', 'BevertonHolt', 'Ricker', 'Hockey'), labels = c('QRF', 'Beverton Holt', 'Ricker', 'Hockey Stick'), name = 'Model') + labs(shape = 'Spawner\nType') + theme(legend.position = "bottom") dodge_width = 0.5 sr_est_comp_p = cap_tbl %>% filter(Population %in% keep_wtsds) %>% mutate(lwrCI = est + qnorm(0.025) * se, lwrCI = if_else(lwrCI <= 0, 1, lwrCI), uprCI = est + qnorm(0.975) * se, uprCI = if_else(uprCI > 1.3e6, 1.3e6, uprCI)) %>% mutate(type = if_else(form == 'QRF', 'QRF', 'S-R')) %>% ggplot(aes(x = Population, y = est, color = form, group = form)) + geom_errorbar(aes(ymin = lwrCI, ymax = uprCI), position = position_dodge(width = dodge_width), width = 0) + geom_point(aes(shape = type), fill = 'white', size = 2, position = position_dodge(width = dodge_width)) + scale_shape_manual(values = c('QRF' = 19, 'S-R' = 21), name = 'Type') + scale_color_brewer(palette = 'Set1', breaks = c('QRF', 'BevertonHolt', 'Ricker', 'Hockey'), labels = c('QRF', 'Beverton Holt', 'Ricker', 'Hockey Stick'), name = 'Model') + # scale_y_continuous(trans = 'log') + # theme(axis.text.x = element_text(angle = 90), # axis.text.y = element_blank()) + theme(axis.text = element_blank(), legend.position = "bottom") + facet_wrap(~ Population, scales = 'free') + labs(y = 'Capacity Estimate', x = NULL)
# make some figures for presentations keep_wtsds = c('Catherine Creek', 'Chiwawa R.', "Entiat River", "Hayden Creek", 'Lostine River', 'Marsh Creek', 'Minam River', "Pahsimeroi River", 'Upper Grande Ronde River', # 'Upper Lemhi', "Upper Mainstem John Day") spawn_rec_comp_p2 = spawn_rec_params %>% # filter out a couple populations with super poor fits filter(Population %in% keep_wtsds) %>% select(Population, data) %>% unnest(cols = "data") %>% distinct() %>% ggplot(aes(x = Spawners, y = Parr)) + # geom_polygon(data = cap_poly %>% # # filter out a couple populations with super poor fits # filter(Population %in% keep_wtsds), # aes(x = x, # y = y, # fill = form), # alpha = 0.2) + geom_line(data = curve_fits %>% # filter out a couple populations with super poor fits filter(Population %in% keep_wtsds), aes(color = form), lwd = 1.5) + geom_hline(data = cap_tbl %>% # filter out a couple populations with super poor fits filter(Population %in% keep_wtsds), aes(yintercept = est, color = form), linetype = 2) + geom_point(aes(shape = Spawner_type), size = 3) + facet_wrap(~ Population, # ncol = 3, nrow = 2, scales = 'free') + scale_color_brewer(palette = 'Set1', breaks = c('QRF', 'BevertonHolt', 'Ricker', 'Hockey'), labels = c('QRF', 'Beverton Holt', 'Ricker', 'Hockey Stick'), name = 'Model') + scale_fill_brewer(palette = 'Set1', breaks = c('QRF', 'BevertonHolt', 'Ricker', 'Hockey'), labels = c('QRF', 'Beverton Holt', 'Ricker', 'Hockey Stick'), name = 'Model') + labs(shape = 'Spawner\nType') + theme(legend.position = "bottom", axis.text = element_text(size = 6)) dodge_width = 0.5 sr_est_comp_p2 = cap_tbl %>% filter(Population %in% keep_wtsds) %>% mutate(lwrCI = est + qnorm(0.025) * se, lwrCI = if_else(lwrCI <= 0, 1, lwrCI), uprCI = est + qnorm(0.975) * se, uprCI = if_else(uprCI > 1.3e6, 1.3e6, uprCI)) %>% mutate(type = if_else(form == 'QRF', 'QRF', 'S-R')) %>% ggplot(aes(x = Population, y = est, color = form, group = form)) + geom_errorbar(aes(ymin = lwrCI, ymax = uprCI), position = position_dodge(width = dodge_width), width = 0) + geom_point(aes(shape = type), fill = 'white', size = 2, position = position_dodge(width = dodge_width)) + scale_shape_manual(values = c('QRF' = 19, 'S-R' = 21), name = 'Type') + scale_color_brewer(palette = 'Set1', breaks = c('QRF', 'BevertonHolt', 'Ricker', 'Hockey'), labels = c('QRF', 'Beverton Holt', 'Ricker', 'Hockey Stick'), name = 'Model') + # scale_y_continuous(trans = 'log') + # theme(axis.text.x = element_text(angle = 90), # axis.text.y = element_blank()) + theme(axis.text = element_blank(), legend.position = "bottom") + facet_wrap(~ Population, nrow = 2, scales = 'free') + labs(y = 'Capacity Estimate', x = NULL) ggsave(here("output/figures/spawn_rec_curves.pdf"), spawn_rec_comp_p2, width = 11.5, height = 5) ggsave(here("output/figures/spawn_rec_comp.pdf"), sr_est_comp_p2, width = 11.5, height = 5)
sr_df = spawn_rec_params %>% filter(Population == "Chiwawa R.") %>% select(Population, data) %>% unnest(cols = "data") %>% distinct() %>% mutate(across(Parr, ~ . / 1e5)) sr_curve = curve_fits %>% mutate(across(Parr, ~ . / 1e5)) %>% filter(Population == "Chiwawa R.") %>% mutate(form = fct_relevel(form, "QRF", after = Inf)) sr_cap = cap_tbl %>% mutate(across(c(est, se), ~ . / 1e5)) %>% filter(Population == "Chiwawa R.") %>% mutate(form = fct_relevel(form, "QRF", after = Inf)) sr_poly = cap_poly %>% mutate(across(c(y), ~ . / 1e5)) %>% filter(Population == "Chiwawa R.") %>% mutate(form = fct_relevel(form, "QRF", after = Inf)) chiw_p1 = sr_df %>% ggplot(aes(x = Spawners, y = Parr)) + geom_point(size = 3) + labs(y = "Parr (x 100,000)") + scale_y_continuous(limits = c(0, 5)) chiw_p2 = chiw_p1 + geom_line(data = sr_curve %>% filter(form != "QRF"), aes(color = form), lwd = 2) + geom_hline(data = sr_cap %>% filter(form != "QRF"), aes(yintercept = est, color = form), linetype = 2, lwd = 1.5) + scale_color_brewer(palette = 'Set1', breaks = c('QRF', 'BevertonHolt', 'Ricker', 'Hockey'), labels = c('QRF', 'Beverton Holt', 'Ricker', 'Hockey Stick'), name = 'Model') + theme(legend.position = "bottom") chiw_p3 = chiw_p1 + geom_line(data = sr_curve, aes(color = form), lwd = 2) + geom_hline(data = sr_cap, aes(yintercept = est, color = form), linetype = 2, lwd = 1.5) + scale_color_brewer(palette = 'Set1', breaks = c('QRF', 'BevertonHolt', 'Ricker', 'Hockey'), labels = c('QRF', 'Beverton Holt', 'Ricker', 'Hockey Stick'), name = 'Model') + theme(legend.position = "bottom") chiw_p4 = chiw_p3 + geom_polygon(data = sr_poly, aes(x = x, y = y, fill = form), alpha = 0.2) + scale_fill_brewer(palette = 'Set1', breaks = c('QRF', 'BevertonHolt', 'Ricker', 'Hockey'), labels = c('QRF', 'Beverton Holt', 'Ricker', 'Hockey Stick'), name = 'Model') chiw_p_list = list(chiw_p1, chiw_p2, chiw_p3, chiw_p4) for(i in 1:length(chiw_p_list)) { ggsave(here("output/figures", paste0("Chiwawa_SR_p", i, ".pdf")), chiw_p_list[[i]], width = 7, height = 6) }
cap_tbl_ci = cap_tbl %>% mutate(lwrCI = est + qnorm(0.025) * se, lwrCI = if_else(lwrCI <= 0, 1, lwrCI), uprCI = est + qnorm(0.975) * se, uprCI = if_else(uprCI > 1.3e6, 1.3e6, uprCI)) comp_df = cap_tbl_ci %>% filter(form != 'QRF') %>% filter(!is.na(est)) %>% inner_join(cap_tbl_ci %>% filter(form == 'QRF') %>% select(Population, form, est) %>% spread(form, est) %>% rename(qrf_est = QRF) %>% left_join(cap_tbl_ci %>% filter(form == 'QRF') %>% select(Population, form, lwrCI) %>% spread(form, lwrCI) %>% rename(qrf_lwrCI = QRF)) %>% left_join(cap_tbl_ci %>% filter(form == 'QRF') %>% select(Population, form, uprCI) %>% spread(form, uprCI) %>% rename(qrf_uprCI = QRF))) %>% mutate(in_qrf_ci = if_else(est <= qrf_uprCI & est >= qrf_lwrCI, T, F), qrf_in_ci = if_else(qrf_est <= uprCI & qrf_est >= lwrCI, T, F), overlap_ci = if_else(lwrCI <= qrf_uprCI & uprCI >= qrf_lwrCI, T, F)) comp_df %>% group_by(form) %>% summarise(n_pops = n_distinct(Population)) %>% full_join(comp_df %>% group_by(form) %>% summarise_at(vars(ends_with("_ci")), list(sum), na.rm = T)) %>% mutate_at(vars(perc_in_qrf_ci = in_qrf_ci, perc_qrf_in_ci = qrf_in_ci, perc_overlap_ci = overlap_ci), list(~ . / n_pops)) comp_df %>% filter(!overlap_ci)
We then plotted the data, as well as the fitted curves, including the one which uses QRF estimates of capacity as one of the parameters, on the same plot, to compare the various estimates.
spawn_rec_comp_p
sr_est_comp_p
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.