In ThomasDelSantoONeill/nash: An Efficient Algorithm for Nash Equilibria when Harvesting Interacting Living Resources
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
Disclaimer
The following vignette serves two main purposes. It is (i) a step-by-step introduction to the nash package, and (ii) a reproducible script where the code for the original article https://doi.org/10.1111/faf.12817 is archived.
Benchmark Models
In the original Fish and Fisheries article, the nash package was applied to three models of distinct complexity with regards to the number of species ($S$) within the stock-complex for which $\mathbf{F}_\textbf{Nash}$ is computed. These are, (i) @Taylor2005 modified two-species competitive Lotka-Volterra (LV) model; and (ii) @ICES2017 Baltic Sea and (iii) @ICES2016 North Sea Ecopath with Ecosim [@Christensen2004;@Steenbeek2016] models, both of which were ported into R using the Rpath https://github.com/NOAA-EDAB/Rpath [@Rpath2016;@Lucey2020].
Complexity increases with respect to the number of interacting species and trophic links modelled; from $2$ species and $2$ links for the modified LV model to $22$ and $69$ compartments and $82$ and $968$ links for the Baltic Sea and North Sea models, respectively. The following snapshots represent these latter two EwE models showcasing as non-grey nodes the $S$ species with commercial value whose harvesting rates we optimised, whilst capturing its effect on the entire food web.
``` {r, echo=FALSE, warning = FALSE, message=FALSE, out.width="100%"}
Plotting routine for food web visualisation of Rpath models (Figure 1
in F&F manuscript).
load("BalticSeaModel.RData")
Libraries
library(data.table)
library(igraph)
library(ggraph)
library(tidygraph)
library(extrafont)
library(cowplot)
loadfonts(device = "all")
Script
NODES
nodes <- data.table(spnum = 1:(Rpath.model$NUM_GROUPS-10),
spname = Rpath.model$Group[1:22])
nodes$sptype <- rep(0,22)
nodes$sptype[4:10] <- 1
nodes[sptype == 0, type.label := factor('Background')]
nodes[sptype == 1, type.label := 'Commercial']
nodes[, biomass := Rpath.model$Biomass[1:22]]
nodes[, weight := abs(biomass - mean(biomass)) / sd(biomass)]
LINKS
DietComp <- Rpath.model$DC[1:22,]
predprey <- data.table(from = row(DietComp)[DietComp > 0],
to = col(DietComp)[DietComp > 0],
weight = DietComp[DietComp > 0])
# # FISHERY
# # Combine landings and discards
catch <- Rpath.model$Landings + Rpath.model$Discards
fishery <- data.table(from = row(catch)[catch > 0],
to = col(catch)[catch > 0],
weight = catch[catch > 0])
# Need to fix these weights as they are not equal between predprey and fleet
# Transform fishery column number to their assigned number in spnum
fleet.num <- nodes[sptype == 3, spnum]
for(i in 1:ncol(catch)) fishery[to == i, to := fleet.num[i]][]
linksbs <- predprey
NETWORK AESTHETIC PROPERTIES
net <- graph_from_data_frame(d=predprey, vertices=nodes, directed=T)
graph <- as_tbl_graph(net)
Create layout
lay <- create_layout(net, 'grid')
Set y value to Trophic Level
lay$y <- Rpath.model$TL[1:22]
Recenter on x-axis
Trophic Level 1
ppnum <- which(lay$y < 2)
n <- length(ppnum)
lay$x[ppnum] <- order(lay$x[ppnum]) * 1/n - 0.5/n
Trophic Levels 2+
TL.groups <- round((max(lay$y) - 2) / 0.5)+1
for(i in 1:TL.groups){
if(i == 1) TL <- 2
TLnum <- which(lay$y >= TL & lay$y < (TL + 0.5))
n <- length(TLnum)
lay$x[TLnum] <- order(lay$x[TLnum]) * 1/n - 0.5/n
TL <- TL + 0.5
}
Defining a theme
foodweb_theme <- theme(
axis.line.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.x = element_blank(),
axis.line.y = element_line(),
axis.ticks.y = element_line(),
axis.title.y = element_text(angle = 360,
margin = margin(
t = 0,
r = 10,
b = 0,
l = 0
)),
axis.text.y.left = element_text(),
panel.grid.major.y = element_line(linetype = 2,
color = rgb(0, 0, 0, alpha = 0.5)),
panel.grid.minor.y = element_blank(),
panel.background = element_blank(),
legend.key = element_blank(),
legend.background = element_blank(),
aspect.ratio = 1,
plot.margin = unit(c(0, 3, 0, 0), "cm"),
text = element_text(family = "Consolas")
)
Defining node colours (c(consumers, detritus, fleets, producers))
node_colours <- c("grey75", "#009688")
Nice labels
nicelabbs <- nodes$spname
Plotting
p1 <- ggraph(lay) +
geom_edge_link0(
aes(edge_colour = linksbs$weight),
edge_width = linksbs$weight * 5,
lineend = "round"
) +
scale_edge_color_gradient(
name = "Link Strength",
low = "#6495ED85",
high = "#C1403D85",
limits = c(0, 1.01),
guide = guide_edge_colorbar(
title = "Link Strength",
title.position = "top",
title.hjust = 0.5,
ticks = TRUE,
label = TRUE,
frame.colour = "black",
barwidth = 7.5,
direction = "horizontal",
position = "top"
)
) +
geom_node_point(
aes(fill = type.label),
shape = 21,
size = 5,
stroke = 1,
color = "#000000",
alpha = 0.85,
show.legend = FALSE
) +
geom_node_label(
aes(label = nicelabbs),
repel = TRUE,
label.padding = 0.75,
label.size = NA,
fill = NA,
size = 2,
family = "Consolas",
fontface = "bold",
max.overlaps = Inf
) +
labs(y = "Trophic\nLevel") +
foodweb_theme +
scale_fill_manual(
values = node_colours,
guide = guide_legend(
direction = "vertical",
position = "top",
title = element_blank()
)
)
p1 + theme(legend.position="none") + labs(subtitle = "Baltic Sea Ecosystem")
```r
load("NorthSeaModel.RData")
## Script
# NODES
nodes <- data.table(spnum = 1:(Rpath.model$NUM_GROUPS-11),
spname = Rpath.model$Group[1:69])
nodes$sptype <- rep(0,69)
nodes$sptype[c(13:20,23,28:29,33:34,38)] <- 1
nodes[sptype == 0, type.label := factor('Background')]
nodes[sptype == 1, type.label := 'Commercial']
nodes[, biomass := Rpath.model$Biomass[1:69]]
nodes[, weight := abs(biomass - mean(biomass)) / sd(biomass)]
# LINKS
DietComp <- Rpath.model$DC[1:69,]
predprey <- data.table(from = row(DietComp)[DietComp > 0],
to = col(DietComp)[DietComp > 0],
weight = DietComp[DietComp > 0])
# # # FISHERY
# # # Combine landings and discards
# catch <- Rpath.model$Landings + Rpath.model$Discards
# fishery <- data.table(from = row(catch)[catch > 0],
# to = col(catch)[catch > 0],
# weight = catch[catch > 0])
# # Need to fix these weights as they are not equal between predprey and fleet
# # Transform fishery column number to their assigned number in spnum
# fleet.num <- nodes[sptype == 3, spnum]
# for(i in 1:ncol(catch)) fishery[to == i, to := fleet.num[i]][]
links <- predprey
# NETWORK AESTHETIC PROPERTIES
net <- graph_from_data_frame(d=predprey, vertices=nodes, directed=T)
graph <- as_tbl_graph(net)
# Create layout
lay <- create_layout(net, 'grid')
# Set y value to Trophic Level
lay$y <- Rpath.model$TL[1:69]
# Recenter on x-axis
# Trophic Level 1
ppnum <- which(lay$y < 2)
n <- length(ppnum)
lay$x[ppnum] <- order(lay$x[ppnum]) * 1/n - 0.5/n
# Trophic Levels 2+
TL.groups <- round((max(lay$y) - 2) / 0.5)
for(i in 1:TL.groups){
if(i == 1) TL <- 2
TLnum <- which(lay$y >= TL & lay$y < (TL + 0.5))
n <- length(TLnum)
lay$x[TLnum] <- order(lay$x[TLnum]) * 1/n - 0.5/n
TL <- TL + 0.5
}
# Defining node colours (c(consumers, detritus, fleets, producers))
node_colours <- c("grey75", "#009688")
# Nice labels
nicelab <- nodes$spname
# Plotting
p2 <- ggraph(lay) +
geom_edge_link0(
aes(edge_colour = links$weight),
edge_width = links$weight * 5,
lineend = "round"
) +
scale_edge_color_gradient(
name = "Link Strength",
low = "#6495ED85",
high = "#C1403D85",
limits = c(0, 1.01),
guide = guide_edge_colorbar(
title = "Link Strength",
title.position = "top",
title.hjust = 0.5,
ticks = TRUE,
label = TRUE,
frame.colour = "black",
barwidth = 7.5,
direction = "horizontal",
position = "top"
)
) +
geom_node_point(
aes(fill = type.label),
shape = 21,
size = 5,
stroke = 1,
color = "#000000",
alpha = 0.85,
show.legend = FALSE
) +
geom_node_label(
aes(label = nicelab),
repel = TRUE,
label.padding = 0.75,
label.size = NA,
fill = NA,
size = 2,
family = "Consolas",
fontface = "bold",
max.overlaps = Inf
) +
labs(y = "Trophic\nLevel") +
foodweb_theme +
scale_fill_manual(
values = node_colours,
guide = guide_legend(
direction = "vertical",
position = "top",
title = element_blank()
)
)
p2 + theme(legend.position="none") + labs(subtitle = "North Sea Ecosystem")
It is important to note that both the Baltic and North Sea models are used as ‘key-run’ parameterisations peer-reviewed by the Working Group on Multispecies Assessment Methods (WGSAM) to inform management advice provided by the International Council for the Exploration of the Sea (ICES).
Running nash with Rpath models
To support users employing Rpath as their operating model, a built-in function fn_rpath has been included within the nash package to provide the input function fn. The function is capable of transparently handling models in which all or some of the exploited compartments are stage structured
(e.g. all species in the Baltic Sea model and Cod, Whiting, Haddock, Saithe and Herring in the North Sea model).
### NOT EVALUATING THIS CHUNK GIVEN THE RUNNING TIME NEEDED TO FIND NE FOR ###
### COMPLEX ECOLOGICAL MODELS. IT IS POSSIBLE TO REPRODUCE THE RESULTS BY ###
### COPY/PASTING IN YOUR LOCAL MACHINE. ###
# Load libraries
library(nash)
library(Rpath)
# Load the BS model
load("BalticSeaModel.RData")
# Commercial stock-complex
spp = c("AdCod", "AdHerring", "AdSprat", "AdFlounder")
# Initialising search with last year of data Fs
par <- as.numeric(tail(Rsim.model$fishing$ForcedFRate[, spp], n = 1))
# Running nash via LV
nash.eq.LV.BS <- nash(
par = par,
fn = fn_rpath,
aged.str = TRUE,
data.years = 10,
rsim.mod = Rsim.model,
IDnames = c(
"JuvCod",
"AdCod",
"JuvHerring",
"AdHerring",
"JuvSprat",
"AdSprat",
"JuvFlounder",
"AdFlounder"
),
method = "LV",
yield.curves = TRUE,
rpath.params = Rpath.parameters,
avg.window = 250,
simul.years = 500,
integration.method = "AB",
conv.criterion = 0.01,
F.increase = 0.01
)
# Load the NS model
load("NorthSeaModel.RData")
# Commercial stock-complex
spp = c(
"AduCod",
"AduWhiting",
"AduHaddock",
"AduSaithe",
"AduHerring",
"NorwayPout",
"Sandeels",
"Plaice",
"Sole"
)
# Initialising search with last year of data Fs
par <- Rsim.model$fishing$ForcedFRate[sim.years, spp]
# Running nash via LV
nash.eq.LV.NS <- nash(
par = par,
fn = fn_rpath,
aged.str = TRUE,
data.years = 23,
rsim.mod = Rsim.model,
IDnames = c(
"JuvCod",
"AduCod",
"JuvWhiting",
"AduWhiting",
"JuvHaddock",
"AduHaddock",
"JuvSaithe",
"AduSaithe",
"JuvHerring",
"AduHerring",
"NorwayPout",
"Sandeels",
"Plaice",
"Sole"
),
method = "LV",
yield.curves = TRUE,
rpath.params = Rpath.parameters,
avg.window = 500,
simul.years = 1000,
integration.method = "AB",
conv.criterion = 0.01
)
Once executed, it is possible to test whether NE-MSY has been achieved by varying the fishing mortality on each stock in turn while keeping the others fixed at $\mathbf{F}_\textbf{Nash}$. This is precisely what the following figures show for all $S$ species and across models, where nash correctly computed the Nash equilibrium fishing mortality rates (vertical dashed lines), such that, for no stock the long-term yeld can be increased by choosing a different fishing mortality.
library(ggplot2)
library(extrafont)
library(reshape2)
loadfonts()
# Personalised theme
theme_tjdso <- theme(
text = element_text(family = "Consolas", size = 20),
panel.background = element_blank(),
panel.border = element_rect(fill = FALSE),
panel.grid.major = element_line(
linetype = 3,
colour = "grey",
size = 0.25
),
panel.grid.minor = element_line(
linetype = 3,
colour = "grey",
size = 0.1
),
strip.background = element_blank(),
strip.text = element_text(size = 20),
strip.placement = "outside",
panel.spacing.x = unit(5, "mm"),
axis.ticks.length = unit(-2, "mm"),
axis.text.x.top = element_blank(),
axis.text.y.right = element_blank(),
axis.title.x.top = element_blank(),
axis.title.y.right = element_blank(),
axis.title = element_text(size = 18),
axis.text.x.bottom = element_text(margin = margin(4, 0, 0, 0, "mm")),
axis.text.y.left = element_text(margin = margin(0, 4, 0, 0, "mm"))
)
nash.eq.LV.BS <- readRDS("BS-NE-Fmsy-LV.rds")
sppname <- c("Cod", "Herring", "Sprat", "Flounder")
# Transform data into something ggplot2 reads
Fvec <- c()
Yvec <- c()
sppvec <- c()
for (i in 1:length(nash.eq.LV.BS$YieldEQ)) {
Fvec <- append(Fvec, nash.eq.LV.BS$YieldEQ[[i]][, 1])
Yvec <- append(Yvec, nash.eq.LV.BS$YieldEQ[[i]][, 2])
sppvec <-
append(sppvec, rep(sppname[i], nrow(nash.eq.LV.BS$YieldEQ[[i]])))
}
Yeq <- data.frame("Spp" = sppvec,
"Fval" = Fvec,
"Yield" = Yvec)
Fnash <- data.frame("Spp" = sppname,
"Fnash" = as.numeric(nash.eq.LV.BS$par))
lab <- round(Fnash$Fnash, digits = 3)
# round(Fnash$Fnash, digits = 3)
lab <- c(expression(paste("F"["Nash"], "=", 0.205)),
expression(paste("F"["Nash"], "=", 0.277)),
expression(paste("F"["Nash"], "=", 0.606)),
expression(paste("F"["Nash"], "=", 0.309)))
Fnash$Labs <- as.character(lab)
Fnash$MSYnash <- nash.eq.LV.BS$value
# Plot
ggplot(data = Yeq, aes(x = Fval, y = Yield)) +
geom_line(size = 1.5) +
geom_vline(
data = Fnash,
aes(xintercept = Fnash),
linetype = "dashed",
size = 1.25
) +
facet_wrap( ~ Spp, scales = "free") +
scale_x_continuous(labels = scales::number_format(accuracy = 0.01),
sec.axis = dup_axis()) +
scale_y_continuous(labels = scales::number_format(accuracy = 0.01),
sec.axis = dup_axis()) +
theme_tjdso +
labs(y = bquote(paste("Yield (Kg x ", 10 ^ 3, " x ", Km ^ -2, ")")),
x = bquote(paste("Fishing Mortality (", yr ^ -1, ")"))) +
geom_text(
data = Fnash,
aes(
x = Fnash,
y = c(0.1, 0.5, 0.1, 0.05),
label = Labs
),
family = "Consolas",
parse = TRUE,
size = 7,
col = 2
)
nash.eq.LV.NS <- readRDS("NS-NE-Fmsy-LV.rds")
sppname <-
c(
"Cod",
"Whiting",
"Haddock",
"Saithe",
"Herring",
"NorwayPout",
"Sandeels",
"Plaice",
"Sole"
)
# Transform data into something ggplot2 reads
Fvec <- c()
Yvec <- c()
sppvec <- c()
for (i in 1:length(nash.eq.LV.NS$YieldEQ)) {
Fvec <- append(Fvec, nash.eq.LV.NS$YieldEQ[[i]][, 1])
Yvec <- append(Yvec, nash.eq.LV.NS$YieldEQ[[i]][, 2])
sppvec <-
append(sppvec, rep(sppname[i], nrow(nash.eq.LV.NS$YieldEQ[[i]])))
}
Yeq <- data.frame("Spp" = sppvec,
"Fval" = Fvec,
"Yield" = Yvec)
Fnash <- data.frame("Spp" = sppname,
"Fnash" = as.numeric(nash.eq.LV.NS$par))
lab <- round(Fnash$Fnash, digits = 3)
lab <- c(
expression(paste("F"["Nash"], " = ", 0.332)),
expression(paste("F"["Nash"], " = ", 0.310)),
expression(paste("F"["Nash"], " = ", 0.203)),
expression(paste("F"["Nash"], " = ", 0.175)),
expression(paste("F"["Nash"], " = ", 0.246)),
expression(paste("F"["Nash"], " = ", 0.663)),
expression(paste("F"["Nash"], " = ", 0.606)),
expression(paste("F"["Nash"], " = ", 0.283)),
expression(paste("F"["Nash"], " = ", 0.171))
)
Fnash$Labs <- as.character(lab)
# Plot
ggplot(data = Yeq, aes(x = Fval, y = Yield)) +
geom_line(size = 1.5) +
geom_vline(
data = Fnash,
aes(xintercept = Fnash),
linetype = "dashed",
size = 1.25
) +
scale_x_continuous(labels = scales::number_format(accuracy = 0.01),
sec.axis = dup_axis()) +
scale_y_continuous(labels = scales::number_format(accuracy = 0.01),
sec.axis = dup_axis()) +
theme_tjdso +
labs(y = bquote(paste("Yield (Kg x ", 10 ^ 3, " x ", Km ^ -2, ")")),
x = bquote(paste("Fishing Mortality (", yr ^ -1, ")"))) +
geom_text(
data = Fnash,
aes(
x = Fnash,
y = c(0.1, 0.02, 0.02, 0.05, 0.2, 0.1, 0.2, 0.2, 0.02),
label = Labs
),
family = "Consolas",
parse = TRUE,
size = 4,
col = 2
) +
facet_wrap( ~ Spp, scales = "free") + theme(text = element_text(size = 10))
Applying nash with Conservation Constraints
It is conceivable that in some cases harvesting at the Nash equilibrium might drive one or more species to extinction. Even if this does not happen, some species might be placed under unacceptable risk levels of stock collapse. Therefore, the nash package includes an option that allows the user to explicitly incorporate conservation constraints [@Matsuda2006], specifying a conservation biomass threshold Bcons below which stock sizes are not permitted to fall.
We showcase this functionality with the Baltic Sea ecosystem by adding a conservation constraint for the biomass of Cod; which we chose arbitrarily by holding Cod $1 \times 10^5 \text{ Tn}$ above the original $B_{\text{Nash,Cod}}$ (i.e. $\approx 84 \times 10^4 \text{ Tn}$).
# Load the BS model and reference Nash equilibrium results
load("Data/BalticSeaModel.RData")
nash.eq.LV.BS <- readRDS("BS-NE-Fmsy-LV.rds")
# Commercial stock-complex
spp = c("AdCod", "AdHerring", "AdSprat", "AdFlounder")
# Initialising search with last year of data Fs
par <- as.numeric(tail(Rsim.model$fishing$ForcedFRate[, spp], n = 1))
# Baltic Sea Modelled Area
model.area <- 240669
# Target for Cod biomass = 1e5 (kg^3)
blim <-
c((
as.numeric(nash.eq.LV.BS$value / nash.eq.LV.BS$par)[1] + (100000 / model.area)
), 0, 0, 0)
# Running nash via LV
nash.eq.LV.BS.CC <- nash(
par = par,
fn = fn_rpath,
aged.str = TRUE,
data.years = 10,
rsim.mod = Rsim.model,
IDnames = c(
"JuvCod",
"AdCod",
"JuvHerring",
"AdHerring",
"JuvSprat",
"AdSprat",
"JuvFlounder",
"AdFlounder"
),
method = "LV",
yield.curves = FALSE,
rpath.params = Rpath.parameters,
avg.window = 250,
simul.years = sim.years,
integration.method = "RK4",
Bcons = blim,
track = TRUE,
conv.criterion = 0.01
)
This results in a new set of $\mathbf{F_{Nash}}$ values in which, as expected, the fishing pressure on Cod needs to be tapered from $0.204 \text{ year}^{−1}$ to $0.150 \text{ year}^{−1}$ with its consequent reduction in yield from $\approx 15 \times 10^4 \text{ Tn year}^{−1}$ to $\approx 13 \times 10^4 \text{ Tn year}^{−1}$. Application of the constraint increased $F_{\text{Nash,Herring}}$ from $0.277 \text{ year}^{−1}$ to $0.392 \text{ year}^{−1}$ generating a surplus of $\approx 17 \times 10^4 \text{ Tn year}^{−1}$ but had little influence on the NE-MSY harvesting rates for Flounder. For Sprat, by contrast, the conservation constraint on Cod led to a higher $F_{\text{Nash,Sprat}}$ of $0.660 \text{ year}^{−1}$ with a marginal reduction in yield of $\approx 1 \times 10^3 \text{ Tn year}^{−1}$.
``` {r, echo = FALSE, warning = FALSE, message=FALSE, out.width="100%"}
Load data for plotting (Figure 2 in EAP manuscript)
nash.eq.LV.BS <- readRDS("BS-NE-Fmsy-LV.rds")
BSCCdataY2 <- readRDS("CCYield_Spp_2.rds")
BSCCdataY3 <- readRDS("CCYield_Spp_3.rds")
BSCCdataY4 <- readRDS("CCYield_Spp_4.rds")
nash.eq.LV.BS.CC <- readRDS("BS-CCNE-Fmsy.rds")
Plotting -----------------------------------------------------------------
library(ggplot2)
library(extrafont)
library(reshape2)
loadfonts()
Baltic Sea
sppname <- c("Cod", "Herring", "Sprat", "Flounder")
Transform data into something ggplot2 reads
Fvec <- c()
Yvec <- c()
sppvec <- c()
for (i in 1:length(nash.eq.LV.BS$YieldEQ)) {
Fvec <- append(Fvec, nash.eq.LV.BS$YieldEQ[[i]][, 1])
Yvec <- append(Yvec, nash.eq.LV.BS$YieldEQ[[i]][, 2])
sppvec <-
append(sppvec, rep(sppname[i], nrow(nash.eq.LV.BS$YieldEQ[[i]])))
}
Yeq <- data.frame("Spp" = sppvec,
"Fval" = Fvec,
"Yield" = Yvec)
Fnash <- data.frame("Spp" = sppname,
"Fnash" = as.numeric(t(tail(
na.omit(nash.eq.LV.BS$par), n = 1
))))
Code is kept at the conservation constraint biomass so populate a matrix with
NA for plotting purposes.
BSCCdataY2[[1]] <- matrix(rep(NA, 100), ncol = 2, nrow = 50)
BSCCdataY2[[3]] <- BSCCdataY3[[3]]
BSCCdataY2[[4]] <- BSCCdataY4[[4]]
Fvec.cons5 <- c()
Yvec.cons5 <- c()
sppvec.cons5 <- c()
for (i in 1:length(BSCCdataY2)) {
Fvec.cons5 <- append(Fvec.cons5, BSCCdataY2[[i]][, 1])
Yvec.cons5 <- append(Yvec.cons5, BSCCdataY2[[i]][, 2])
sppvec.cons5 <-
append(sppvec.cons5, rep(sppname[i], nrow(BSCCdataY2[[i]])))
}
Yeq.cons5 <- data.frame("Spp" = sppvec.cons5,
"Fval" = Fvec.cons5,
"Yield" = Yvec.cons5)
Fnash.cons5 <- data.frame("Spp" = sppname,
"Fnash.cons5" = as.numeric(t(tail(
na.omit(nash.eq.LV.BS.CC$par), n = 1
))))
Fnash.cons5$Yield <- c(nash.eq.LV.BS.CC$value[1], rep(NA, 3))
lab <- round(Fnash.cons5$Fnash, digits = 3)
lab <- c(expression(paste("CC-F"["Nash"], " = ", 0.150)),
expression(paste("CC-F"["Nash"], " = ", 0.392)),
expression(paste("CC-F"["Nash"], " = ", 0.660)),
expression(paste("CC-F"["Nash"], " = ", 0.309)))
Fnash.cons5$Labs <- as.character(lab)
lab <- round(Fnash$Fnash, digits = 3)
lab <- c(expression(paste("F"["Nash"], " = ", 0.205)),
expression(paste("F"["Nash"], " = ", 0.277)),
expression(paste("F"["Nash"], " = ", 0.606)),
expression(paste("F"["Nash"], " = ", 0.309)))
Fnash$Labs <- as.character(lab)
Personalised theme
theme_tjdso <-
theme(
text = element_text(family = "Consolas", size = 20),
panel.background = element_blank(),
panel.border = element_rect(fill = FALSE),
panel.grid.major = element_line(
linetype = 3,
colour = "grey",
size = 0.25
),
panel.grid.minor = element_line(
linetype = 3,
colour = "grey",
size = 0.1
),
strip.background = element_blank(),
strip.text = element_text(size = 20),
strip.placement = "outside",
panel.spacing.x = unit(5, "mm"),
axis.ticks.length = unit(-2, "mm"),
axis.text.x.top = element_blank(),
axis.text.y.right = element_blank(),
axis.title.x.top = element_blank(),
axis.title.y.right = element_blank(),
axis.title = element_text(size = 18),
axis.text.x.bottom = element_text(margin = margin(4, 0,
0, 0,
"mm")),
axis.text.y.left = element_text(margin = margin(0, 4,
0, 0,
"mm"))
)
Plot
ggplot(data = Yeq, aes(x = Fval, y = Yield)) +
geom_line(size = 1.5, col = 2) +
geom_vline(
data = Fnash,
aes(xintercept = Fnash),
linetype = "dashed",
size = 1.25,
col = 2
) +
geom_text(
data = Fnash,
aes(
x = Fnash,
y = c(0.1, 0.5, 0.1, 0.05),
label = Labs
),
family = "Consolas",
parse = TRUE,
size = 5,
col = 2
) +
facet_wrap( ~ Spp, scales = "free") +
scale_x_continuous(labels = scales::number_format(accuracy = 0.01),
sec.axis = dup_axis()) +
scale_y_continuous(labels = scales::number_format(accuracy = 0.01),
sec.axis = dup_axis()) +
theme_tjdso +
labs(y = bquote(paste("Yield (Kg x ", 10 ^ 3, " x ", Km ^ -2, ")")),
x = bquote(paste("Fishing Mortality (", yr ^ -1, ")"))) +
# geom_line(data = Yeq.cons, aes(x = Fval, y = Yield), size = 1.5, col = "blue") +
geom_line(
data = Yeq.cons5,
aes(x = Fval, y = Yield),
size = 1.5,
col = "steelblue"
) +
# geom_vline(data = Fnash.cons, aes(xintercept = Fnash.cons),
# linetype = "dashed", size = 1.25, col = "blue") +
geom_vline(
data = Fnash.cons5,
aes(xintercept = Fnash.cons5),
linetype = "dashed",
size = 1.25,
col = "steelblue"
) +
geom_point(
data = Fnash.cons5,
aes(x = Fnash.cons5, y = Yield),
col = "steelblue",
size = 4
) +
geom_text(
data = Fnash.cons5,
aes(
x = Fnash.cons5,
y = c(0.25, 1.5, 0.3, 0.15),
label = Labs
),
family = "Consolas",
parse = TRUE,
size = 5,
col = "steelblue"
)
# `nash`'s `LV` algorithm *vs* `round-robin`
A key advantage of the `LV` method over the simple `round-robin` method, that is, sequential optimisation of each $F_i$ until convergence is reached, is that it requires much less computation time, especially for complex ecological communities. As a
platform-independent metric of performance, we compared the number of objective function evaluations between the two methods. For this analysis, we applied both methods across the three tested models with $n=77$ different initial harvesting rate values.
```r
SMFmat <- readRDS("SMFmat.rds")
BSFmat <- readRDS("BSFmat.rds")
NSFmat <- readRDS("NSFmat.rds")
# Theme
theme_tjdso <-
theme(
text = element_text(family = "Consolas", size = 20),
panel.background = element_blank(),
panel.border = element_rect(fill = FALSE),
panel.grid.major = element_line(
linetype = 3,
colour = "grey",
size = 0.25
),
panel.grid.minor = element_line(
linetype = 3,
colour = "grey",
size = 0.1
),
strip.background = element_blank(),
strip.text = element_text(size = 20),
strip.placement = "outside",
panel.spacing.x = unit(5, "mm"),
axis.ticks.length = unit(-2, "mm"),
axis.text.x.top = element_blank(),
axis.text.y.right = element_blank(),
axis.title.x.top = element_blank(),
axis.title.y.right = element_blank(),
axis.title = element_text(size = 18),
axis.text.x.bottom = element_text(margin = margin(4, 0,
0, 0,
"mm")),
axis.text.y.left = element_text(margin = margin(0, 4,
0, 0,
"mm"))
)
# SM
SMF <- reshape2::melt(SMFmat)[, -1]
SMF$Var2 <- rep(c("Spp.1", "Spp.2"), each = 29)
colnames(SMF) <- c("Species", "Fini")
SMF <- as.data.frame(SMF)
SMref <- data.frame(Species = c("Spp.1", "Spp.2"),
Ref = c(0.2, 0.3))
ggplot(data = SMF, aes(x = Species, y = Fini)) +
geom_jitter(col = "grey75", alpha = 0.75) +
geom_boxplot(size = 1.2,
col = "black",
alpha = 0.001) +
labs(x = "",
y = "") +
theme_tjdso +
geom_point(data = SMref, aes(x = Species, y = Ref), col = 2) +
labs(x = "",
y = bquote(paste("Fishing Mortality (", yr ^ -1, ")")))
# BS
colnames(BSFmat) <- c("Cod", "Herring", "Sprat", "Flounder")
BSF <- reshape2::melt(BSFmat)[, -1]
colnames(BSF) <- c("Species", "Fini")
BSF <- as.data.frame(BSF)
BSref <-
data.frame(
Species = c("Cod", "Herring", "Sprat", "Flounder"),
Ref = c(0.4342084, 0.06663042, 0.1139571, 0.3030173)
)
ggplot(data = BSF, aes(x = Species, y = Fini)) +
geom_jitter(col = "grey75", alpha = 0.75) +
geom_boxplot(size = 1.2,
col = "black",
alpha = 0.001) +
labs(x = "",
y = bquote(paste("Fishing Mortality (", yr ^ -1, ")"))) +
theme_tjdso +
geom_point(data = BSref, aes(x = Species, y = Ref), col = 2)
# NS
colnames(NSFmat) <-
c(
"Cod",
"Whiting",
"Haddock",
"Saithe",
"Herring",
"NorwayPout",
"Sandeels",
"Plaice",
"Sole"
)
NSF <- reshape2::melt(NSFmat)[, -1]
colnames(NSF) <- c("Species", "Fini")
NSF <- as.data.frame(NSF)
NSref <-
data.frame(
Species = c(
"Cod",
"Whiting",
"Haddock",
"Saithe",
"Herring",
"NorwayPout",
"Sandeels",
"Plaice",
"Sole"
),
Ref = c(
0.2433686,
0.3393892,
0.1778537,
0.1137075,
0.16818,
0.108255,
0.1468737,
0.2017529,
0.1575716
)
)
ggplot(data = NSF, aes(x = Species, y = Fini)) +
geom_jitter(col = "grey75", alpha = 0.75) +
geom_boxplot(size = 1.2,
col = "black",
alpha = 0.001) +
labs(x = "",
y = "") +
theme_tjdso +
scale_x_discrete(guide = guide_axis(angle = 45)) +
geom_point(data = NSref, aes(x = Species, y = Ref), col = 2) +
labs(x = "",
y = bquote(paste("Fishing Mortality (", yr ^ -1, ")")))
In all cases, the LV method outperformed the round-robin method for the same convergence threshold (set to a value of $0.01$ for this analysis). Noteworthy is, in particular, that for the North Sea model,
the LV method requires only $1.88$ times more function evaluations than for the Baltic Sea model, even though the dimension of the $\mathbf{F}_\textbf{Nash}$ vector increases by $2.25$ from $4$ to $9$. Using the round-robin method, by comparison, the number of function evaluations increases $7$-fold.
SMdata <- readRDS("SM_fn.counts_data.rds")
BSdata <- readRDS("BS_fn.counts_data.rds")
NSdata <- readRDS("NS_fn.counts_data.rds")
performance <- data.frame(Simple.Model = c(colMeans(SMdata)),
BalticSea.Model = c(colMeans(BSdata)),
NorthSea.Model = c(colMeans(NSdata, na.rm = TRUE)))
colnames(performance) <- c("Simple Model", "Baltic Sea Model",
"North Sea Model")
rownames(performance) <- c("LV", "round-robin")
performance <- as.matrix(performance)
barplot(performance, beside = TRUE, ylim = c(0,2000),
legend = rownames(performance),
args.legend = list(x = "topleft", bty = "n",
title = "Method",
title.adj = 0.1),
ylab = "Function Evaluations (no.)", cex.lab = 1, cex.names = 1,
cex.axis = 1, border = TRUE, col = c("gray45","gray85"))
# width.error <- 0.025
segments(x0 = 1.5, y0 = mean(na.omit(SMdata[,1]))-sd(na.omit(SMdata[,1])),
x1 = 1.5, y1 = mean(na.omit(SMdata[,1]))+sd(na.omit(SMdata[,1])), lwd = 2)
# segments(x0 = 1.5-width.error, y0 = mean(SMdata[,1])+sd(SMdata[,1]),
# x1 = 1.5+width.error, y1 = mean(SMdata[,1])+sd(SMdata[,1]), lwd = 2)
# segments(x0 = 1.5-width.error, y0 = mean(SMdata[,1])-sd(SMdata[,1]),
# x1 = 1.5+width.error, y1 = mean(SMdata[,1])-sd(SMdata[,1]), lwd = 2)
segments(x0 = 2.5, y0 = mean(na.omit(SMdata[,2]))-sd(na.omit(SMdata[,2])),
x1 = 2.5, y1 = mean(na.omit(SMdata[,2]))+sd(na.omit(SMdata[,2])), lwd = 2)
# segments(x0 = 2.5-width.error, y0 = mean(SMdata[,2])+sd(SMdata[,2]),
# x1 = 2.5+width.error, y1 = mean(SMdata[,2])+sd(SMdata[,2]), lwd = 2)
# segments(x0 = 2.5-width.error, y0 = mean(SMdata[,2])-sd(SMdata[,2]),
# x1 = 2.5+width.error, y1 = mean(SMdata[,2])-sd(SMdata[,2]), lwd = 2)
segments(x0 = 4.5, y0 = mean(na.omit(BSdata[,1]))-sd(na.omit(BSdata[,1])),
x1 = 4.5, y1 = mean(na.omit(BSdata[,1]))+sd(na.omit(BSdata[,1])),
lwd = 2)
# segments(x0 = 4.5-width.error, y0 = mean(na.omit(BSdata[,1]))+sd(na.omit(BSdata[,1])),
# x1 = 4.5+width.error, y1 = mean(na.omit(BSdata[,1]))+sd(na.omit(BSdata[,1])),
# lwd = 2)
# segments(x0 = 4.5-width.error, y0 = mean(na.omit(BSdata[,1]))-sd(na.omit(BSdata[,1])),
# x1 = 4.5+width.error, y1 = mean(na.omit(BSdata[,1]))-sd(na.omit(BSdata[,1])),
# lwd = 2)
segments(x0 = 5.5, y0 = mean(na.omit(BSdata[,2]))-sd(na.omit(BSdata[,2])),
x1 = 5.5, y1 = mean(na.omit(BSdata[,2]))+sd(na.omit(BSdata[,2])),
lwd = 2)
# segments(x0 = 5.5-width.error, y0 = mean(na.omit(BSdata[,2]))+sd(na.omit(BSdata[,2])),
# x1 = 5.5+width.error, y1 = mean(na.omit(BSdata[,2]))+sd(na.omit(BSdata[,2])),
# lwd = 2)
# segments(x0 = 5.5-width.error, y0 = mean(na.omit(BSdata[,2]))-sd(na.omit(BSdata[,2])),
# x1 = 5.5+width.error, y1 = mean(na.omit(BSdata[,2]))-sd(na.omit(BSdata[,2])),
# lwd = 2)
segments(x0 = 7.5, y0 = mean(na.omit(NSdata[,1]))-sd(na.omit(NSdata[,1])),
x1 = 7.5, y1 = mean(na.omit(NSdata[,1]))+sd(na.omit(NSdata[,1])),
lwd = 2)
# segments(x0 = 7.5-width.error, y0 = mean(na.omit(NSdata[,1]))+sd(na.omit(NSdata[,1])),
# x1 = 7.5+width.error, y1 = mean(na.omit(NSdata[,1]))+sd(na.omit(NSdata[,1])),
# lwd = 2)
# segments(x0 = 7.5-width.error, y0 = mean(na.omit(NSdata[,1]))-sd(na.omit(NSdata[,1])),
# x1 = 7.5+width.error, y1 = mean(na.omit(NSdata[,1]))-sd(na.omit(NSdata[,1])),
# lwd = 2)
segments(x0 = 8.5, y0 = mean(na.omit(NSdata[,2]))-sd(na.omit(NSdata[,2])),
x1 = 8.5, y1 = mean(na.omit(NSdata[,2]))+sd(na.omit(NSdata[,2])),
lwd = 2)
grid()
box()
References
ThomasDelSantoONeill/nash documentation built on Aug. 10, 2024, 1:49 a.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.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
Disclaimer
The following vignette serves two main purposes. It is (i) a step-by-step introduction to the nash package, and (ii) a reproducible script where the code for the original article https://doi.org/10.1111/faf.12817 is archived.
Benchmark Models
In the original Fish and Fisheries article, the nash package was applied to three models of distinct complexity with regards to the number of species ($S$) within the stock-complex for which $\mathbf{F}_\textbf{Nash}$ is computed. These are, (i) @Taylor2005 modified two-species competitive Lotka-Volterra (LV) model; and (ii) @ICES2017 Baltic Sea and (iii) @ICES2016 North Sea Ecopath with Ecosim [@Christensen2004;@Steenbeek2016] models, both of which were ported into R using the Rpath https://github.com/NOAA-EDAB/Rpath [@Rpath2016;@Lucey2020].
Complexity increases with respect to the number of interacting species and trophic links modelled; from $2$ species and $2$ links for the modified LV model to $22$ and $69$ compartments and $82$ and $968$ links for the Baltic Sea and North Sea models, respectively. The following snapshots represent these latter two EwE models showcasing as non-grey nodes the $S$ species with commercial value whose harvesting rates we optimised, whilst capturing its effect on the entire food web.
``` {r, echo=FALSE, warning = FALSE, message=FALSE, out.width="100%"}
Plotting routine for food web visualisation of Rpath models (Figure 1
in F&F manuscript).
load("BalticSeaModel.RData")
Libraries
library(data.table) library(igraph) library(ggraph) library(tidygraph) library(extrafont) library(cowplot) loadfonts(device = "all")
Script
NODES
nodes <- data.table(spnum = 1:(Rpath.model$NUM_GROUPS-10), spname = Rpath.model$Group[1:22]) nodes$sptype <- rep(0,22) nodes$sptype[4:10] <- 1 nodes[sptype == 0, type.label := factor('Background')] nodes[sptype == 1, type.label := 'Commercial'] nodes[, biomass := Rpath.model$Biomass[1:22]] nodes[, weight := abs(biomass - mean(biomass)) / sd(biomass)]
LINKS
DietComp <- Rpath.model$DC[1:22,] predprey <- data.table(from = row(DietComp)[DietComp > 0], to = col(DietComp)[DietComp > 0], weight = DietComp[DietComp > 0])
# # FISHERY
# # Combine landings and discards
catch <- Rpath.model$Landings + Rpath.model$Discards
fishery <- data.table(from = row(catch)[catch > 0],
to = col(catch)[catch > 0],
weight = catch[catch > 0])
# Need to fix these weights as they are not equal between predprey and fleet
# Transform fishery column number to their assigned number in spnum
fleet.num <- nodes[sptype == 3, spnum]
for(i in 1:ncol(catch)) fishery[to == i, to := fleet.num[i]][]
linksbs <- predprey
NETWORK AESTHETIC PROPERTIES
net <- graph_from_data_frame(d=predprey, vertices=nodes, directed=T) graph <- as_tbl_graph(net)
Create layout
lay <- create_layout(net, 'grid')
Set y value to Trophic Level
lay$y <- Rpath.model$TL[1:22]
Recenter on x-axis
Trophic Level 1
ppnum <- which(lay$y < 2) n <- length(ppnum) lay$x[ppnum] <- order(lay$x[ppnum]) * 1/n - 0.5/n
Trophic Levels 2+
TL.groups <- round((max(lay$y) - 2) / 0.5)+1 for(i in 1:TL.groups){ if(i == 1) TL <- 2 TLnum <- which(lay$y >= TL & lay$y < (TL + 0.5)) n <- length(TLnum) lay$x[TLnum] <- order(lay$x[TLnum]) * 1/n - 0.5/n TL <- TL + 0.5 }
Defining a theme
foodweb_theme <- theme( axis.line.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.x = element_blank(), panel.grid.minor.x = element_blank(), panel.grid.major.x = element_blank(), axis.line.y = element_line(), axis.ticks.y = element_line(), axis.title.y = element_text(angle = 360, margin = margin( t = 0, r = 10, b = 0, l = 0 )), axis.text.y.left = element_text(), panel.grid.major.y = element_line(linetype = 2, color = rgb(0, 0, 0, alpha = 0.5)), panel.grid.minor.y = element_blank(), panel.background = element_blank(), legend.key = element_blank(), legend.background = element_blank(), aspect.ratio = 1, plot.margin = unit(c(0, 3, 0, 0), "cm"), text = element_text(family = "Consolas") )
Defining node colours (c(consumers, detritus, fleets, producers))
node_colours <- c("grey75", "#009688")
Nice labels
nicelabbs <- nodes$spname
Plotting
p1 <- ggraph(lay) + geom_edge_link0( aes(edge_colour = linksbs$weight), edge_width = linksbs$weight * 5, lineend = "round" ) + scale_edge_color_gradient( name = "Link Strength", low = "#6495ED85", high = "#C1403D85", limits = c(0, 1.01), guide = guide_edge_colorbar( title = "Link Strength", title.position = "top", title.hjust = 0.5, ticks = TRUE, label = TRUE, frame.colour = "black", barwidth = 7.5, direction = "horizontal", position = "top" ) ) + geom_node_point( aes(fill = type.label), shape = 21, size = 5, stroke = 1, color = "#000000", alpha = 0.85, show.legend = FALSE ) + geom_node_label( aes(label = nicelabbs), repel = TRUE, label.padding = 0.75, label.size = NA, fill = NA, size = 2, family = "Consolas", fontface = "bold", max.overlaps = Inf ) + labs(y = "Trophic\nLevel") + foodweb_theme + scale_fill_manual( values = node_colours, guide = guide_legend( direction = "vertical", position = "top", title = element_blank() ) ) p1 + theme(legend.position="none") + labs(subtitle = "Baltic Sea Ecosystem")
```r load("NorthSeaModel.RData") ## Script # NODES nodes <- data.table(spnum = 1:(Rpath.model$NUM_GROUPS-11), spname = Rpath.model$Group[1:69]) nodes$sptype <- rep(0,69) nodes$sptype[c(13:20,23,28:29,33:34,38)] <- 1 nodes[sptype == 0, type.label := factor('Background')] nodes[sptype == 1, type.label := 'Commercial'] nodes[, biomass := Rpath.model$Biomass[1:69]] nodes[, weight := abs(biomass - mean(biomass)) / sd(biomass)] # LINKS DietComp <- Rpath.model$DC[1:69,] predprey <- data.table(from = row(DietComp)[DietComp > 0], to = col(DietComp)[DietComp > 0], weight = DietComp[DietComp > 0]) # # # FISHERY # # # Combine landings and discards # catch <- Rpath.model$Landings + Rpath.model$Discards # fishery <- data.table(from = row(catch)[catch > 0], # to = col(catch)[catch > 0], # weight = catch[catch > 0]) # # Need to fix these weights as they are not equal between predprey and fleet # # Transform fishery column number to their assigned number in spnum # fleet.num <- nodes[sptype == 3, spnum] # for(i in 1:ncol(catch)) fishery[to == i, to := fleet.num[i]][] links <- predprey # NETWORK AESTHETIC PROPERTIES net <- graph_from_data_frame(d=predprey, vertices=nodes, directed=T) graph <- as_tbl_graph(net) # Create layout lay <- create_layout(net, 'grid') # Set y value to Trophic Level lay$y <- Rpath.model$TL[1:69] # Recenter on x-axis # Trophic Level 1 ppnum <- which(lay$y < 2) n <- length(ppnum) lay$x[ppnum] <- order(lay$x[ppnum]) * 1/n - 0.5/n # Trophic Levels 2+ TL.groups <- round((max(lay$y) - 2) / 0.5) for(i in 1:TL.groups){ if(i == 1) TL <- 2 TLnum <- which(lay$y >= TL & lay$y < (TL + 0.5)) n <- length(TLnum) lay$x[TLnum] <- order(lay$x[TLnum]) * 1/n - 0.5/n TL <- TL + 0.5 } # Defining node colours (c(consumers, detritus, fleets, producers)) node_colours <- c("grey75", "#009688") # Nice labels nicelab <- nodes$spname # Plotting p2 <- ggraph(lay) + geom_edge_link0( aes(edge_colour = links$weight), edge_width = links$weight * 5, lineend = "round" ) + scale_edge_color_gradient( name = "Link Strength", low = "#6495ED85", high = "#C1403D85", limits = c(0, 1.01), guide = guide_edge_colorbar( title = "Link Strength", title.position = "top", title.hjust = 0.5, ticks = TRUE, label = TRUE, frame.colour = "black", barwidth = 7.5, direction = "horizontal", position = "top" ) ) + geom_node_point( aes(fill = type.label), shape = 21, size = 5, stroke = 1, color = "#000000", alpha = 0.85, show.legend = FALSE ) + geom_node_label( aes(label = nicelab), repel = TRUE, label.padding = 0.75, label.size = NA, fill = NA, size = 2, family = "Consolas", fontface = "bold", max.overlaps = Inf ) + labs(y = "Trophic\nLevel") + foodweb_theme + scale_fill_manual( values = node_colours, guide = guide_legend( direction = "vertical", position = "top", title = element_blank() ) ) p2 + theme(legend.position="none") + labs(subtitle = "North Sea Ecosystem")
It is important to note that both the Baltic and North Sea models are used as ‘key-run’ parameterisations peer-reviewed by the Working Group on Multispecies Assessment Methods (WGSAM) to inform management advice provided by the International Council for the Exploration of the Sea (ICES).
Running nash with Rpath models
To support users employing Rpath as their operating model, a built-in function fn_rpath has been included within the nash package to provide the input function fn. The function is capable of transparently handling models in which all or some of the exploited compartments are stage structured
(e.g. all species in the Baltic Sea model and Cod, Whiting, Haddock, Saithe and Herring in the North Sea model).
### NOT EVALUATING THIS CHUNK GIVEN THE RUNNING TIME NEEDED TO FIND NE FOR ### ### COMPLEX ECOLOGICAL MODELS. IT IS POSSIBLE TO REPRODUCE THE RESULTS BY ### ### COPY/PASTING IN YOUR LOCAL MACHINE. ### # Load libraries library(nash) library(Rpath) # Load the BS model load("BalticSeaModel.RData") # Commercial stock-complex spp = c("AdCod", "AdHerring", "AdSprat", "AdFlounder") # Initialising search with last year of data Fs par <- as.numeric(tail(Rsim.model$fishing$ForcedFRate[, spp], n = 1)) # Running nash via LV nash.eq.LV.BS <- nash( par = par, fn = fn_rpath, aged.str = TRUE, data.years = 10, rsim.mod = Rsim.model, IDnames = c( "JuvCod", "AdCod", "JuvHerring", "AdHerring", "JuvSprat", "AdSprat", "JuvFlounder", "AdFlounder" ), method = "LV", yield.curves = TRUE, rpath.params = Rpath.parameters, avg.window = 250, simul.years = 500, integration.method = "AB", conv.criterion = 0.01, F.increase = 0.01 ) # Load the NS model load("NorthSeaModel.RData") # Commercial stock-complex spp = c( "AduCod", "AduWhiting", "AduHaddock", "AduSaithe", "AduHerring", "NorwayPout", "Sandeels", "Plaice", "Sole" ) # Initialising search with last year of data Fs par <- Rsim.model$fishing$ForcedFRate[sim.years, spp] # Running nash via LV nash.eq.LV.NS <- nash( par = par, fn = fn_rpath, aged.str = TRUE, data.years = 23, rsim.mod = Rsim.model, IDnames = c( "JuvCod", "AduCod", "JuvWhiting", "AduWhiting", "JuvHaddock", "AduHaddock", "JuvSaithe", "AduSaithe", "JuvHerring", "AduHerring", "NorwayPout", "Sandeels", "Plaice", "Sole" ), method = "LV", yield.curves = TRUE, rpath.params = Rpath.parameters, avg.window = 500, simul.years = 1000, integration.method = "AB", conv.criterion = 0.01 )
Once executed, it is possible to test whether NE-MSY has been achieved by varying the fishing mortality on each stock in turn while keeping the others fixed at $\mathbf{F}_\textbf{Nash}$. This is precisely what the following figures show for all $S$ species and across models, where nash correctly computed the Nash equilibrium fishing mortality rates (vertical dashed lines), such that, for no stock the long-term yeld can be increased by choosing a different fishing mortality.
library(ggplot2) library(extrafont) library(reshape2) loadfonts() # Personalised theme theme_tjdso <- theme( text = element_text(family = "Consolas", size = 20), panel.background = element_blank(), panel.border = element_rect(fill = FALSE), panel.grid.major = element_line( linetype = 3, colour = "grey", size = 0.25 ), panel.grid.minor = element_line( linetype = 3, colour = "grey", size = 0.1 ), strip.background = element_blank(), strip.text = element_text(size = 20), strip.placement = "outside", panel.spacing.x = unit(5, "mm"), axis.ticks.length = unit(-2, "mm"), axis.text.x.top = element_blank(), axis.text.y.right = element_blank(), axis.title.x.top = element_blank(), axis.title.y.right = element_blank(), axis.title = element_text(size = 18), axis.text.x.bottom = element_text(margin = margin(4, 0, 0, 0, "mm")), axis.text.y.left = element_text(margin = margin(0, 4, 0, 0, "mm")) ) nash.eq.LV.BS <- readRDS("BS-NE-Fmsy-LV.rds") sppname <- c("Cod", "Herring", "Sprat", "Flounder") # Transform data into something ggplot2 reads Fvec <- c() Yvec <- c() sppvec <- c() for (i in 1:length(nash.eq.LV.BS$YieldEQ)) { Fvec <- append(Fvec, nash.eq.LV.BS$YieldEQ[[i]][, 1]) Yvec <- append(Yvec, nash.eq.LV.BS$YieldEQ[[i]][, 2]) sppvec <- append(sppvec, rep(sppname[i], nrow(nash.eq.LV.BS$YieldEQ[[i]]))) } Yeq <- data.frame("Spp" = sppvec, "Fval" = Fvec, "Yield" = Yvec) Fnash <- data.frame("Spp" = sppname, "Fnash" = as.numeric(nash.eq.LV.BS$par)) lab <- round(Fnash$Fnash, digits = 3) # round(Fnash$Fnash, digits = 3) lab <- c(expression(paste("F"["Nash"], "=", 0.205)), expression(paste("F"["Nash"], "=", 0.277)), expression(paste("F"["Nash"], "=", 0.606)), expression(paste("F"["Nash"], "=", 0.309))) Fnash$Labs <- as.character(lab) Fnash$MSYnash <- nash.eq.LV.BS$value # Plot ggplot(data = Yeq, aes(x = Fval, y = Yield)) + geom_line(size = 1.5) + geom_vline( data = Fnash, aes(xintercept = Fnash), linetype = "dashed", size = 1.25 ) + facet_wrap( ~ Spp, scales = "free") + scale_x_continuous(labels = scales::number_format(accuracy = 0.01), sec.axis = dup_axis()) + scale_y_continuous(labels = scales::number_format(accuracy = 0.01), sec.axis = dup_axis()) + theme_tjdso + labs(y = bquote(paste("Yield (Kg x ", 10 ^ 3, " x ", Km ^ -2, ")")), x = bquote(paste("Fishing Mortality (", yr ^ -1, ")"))) + geom_text( data = Fnash, aes( x = Fnash, y = c(0.1, 0.5, 0.1, 0.05), label = Labs ), family = "Consolas", parse = TRUE, size = 7, col = 2 ) nash.eq.LV.NS <- readRDS("NS-NE-Fmsy-LV.rds") sppname <- c( "Cod", "Whiting", "Haddock", "Saithe", "Herring", "NorwayPout", "Sandeels", "Plaice", "Sole" ) # Transform data into something ggplot2 reads Fvec <- c() Yvec <- c() sppvec <- c() for (i in 1:length(nash.eq.LV.NS$YieldEQ)) { Fvec <- append(Fvec, nash.eq.LV.NS$YieldEQ[[i]][, 1]) Yvec <- append(Yvec, nash.eq.LV.NS$YieldEQ[[i]][, 2]) sppvec <- append(sppvec, rep(sppname[i], nrow(nash.eq.LV.NS$YieldEQ[[i]]))) } Yeq <- data.frame("Spp" = sppvec, "Fval" = Fvec, "Yield" = Yvec) Fnash <- data.frame("Spp" = sppname, "Fnash" = as.numeric(nash.eq.LV.NS$par)) lab <- round(Fnash$Fnash, digits = 3) lab <- c( expression(paste("F"["Nash"], " = ", 0.332)), expression(paste("F"["Nash"], " = ", 0.310)), expression(paste("F"["Nash"], " = ", 0.203)), expression(paste("F"["Nash"], " = ", 0.175)), expression(paste("F"["Nash"], " = ", 0.246)), expression(paste("F"["Nash"], " = ", 0.663)), expression(paste("F"["Nash"], " = ", 0.606)), expression(paste("F"["Nash"], " = ", 0.283)), expression(paste("F"["Nash"], " = ", 0.171)) ) Fnash$Labs <- as.character(lab) # Plot ggplot(data = Yeq, aes(x = Fval, y = Yield)) + geom_line(size = 1.5) + geom_vline( data = Fnash, aes(xintercept = Fnash), linetype = "dashed", size = 1.25 ) + scale_x_continuous(labels = scales::number_format(accuracy = 0.01), sec.axis = dup_axis()) + scale_y_continuous(labels = scales::number_format(accuracy = 0.01), sec.axis = dup_axis()) + theme_tjdso + labs(y = bquote(paste("Yield (Kg x ", 10 ^ 3, " x ", Km ^ -2, ")")), x = bquote(paste("Fishing Mortality (", yr ^ -1, ")"))) + geom_text( data = Fnash, aes( x = Fnash, y = c(0.1, 0.02, 0.02, 0.05, 0.2, 0.1, 0.2, 0.2, 0.02), label = Labs ), family = "Consolas", parse = TRUE, size = 4, col = 2 ) + facet_wrap( ~ Spp, scales = "free") + theme(text = element_text(size = 10))
Applying nash with Conservation Constraints
It is conceivable that in some cases harvesting at the Nash equilibrium might drive one or more species to extinction. Even if this does not happen, some species might be placed under unacceptable risk levels of stock collapse. Therefore, the nash package includes an option that allows the user to explicitly incorporate conservation constraints [@Matsuda2006], specifying a conservation biomass threshold Bcons below which stock sizes are not permitted to fall.
We showcase this functionality with the Baltic Sea ecosystem by adding a conservation constraint for the biomass of Cod; which we chose arbitrarily by holding Cod $1 \times 10^5 \text{ Tn}$ above the original $B_{\text{Nash,Cod}}$ (i.e. $\approx 84 \times 10^4 \text{ Tn}$).
# Load the BS model and reference Nash equilibrium results load("Data/BalticSeaModel.RData") nash.eq.LV.BS <- readRDS("BS-NE-Fmsy-LV.rds") # Commercial stock-complex spp = c("AdCod", "AdHerring", "AdSprat", "AdFlounder") # Initialising search with last year of data Fs par <- as.numeric(tail(Rsim.model$fishing$ForcedFRate[, spp], n = 1)) # Baltic Sea Modelled Area model.area <- 240669 # Target for Cod biomass = 1e5 (kg^3) blim <- c(( as.numeric(nash.eq.LV.BS$value / nash.eq.LV.BS$par)[1] + (100000 / model.area) ), 0, 0, 0) # Running nash via LV nash.eq.LV.BS.CC <- nash( par = par, fn = fn_rpath, aged.str = TRUE, data.years = 10, rsim.mod = Rsim.model, IDnames = c( "JuvCod", "AdCod", "JuvHerring", "AdHerring", "JuvSprat", "AdSprat", "JuvFlounder", "AdFlounder" ), method = "LV", yield.curves = FALSE, rpath.params = Rpath.parameters, avg.window = 250, simul.years = sim.years, integration.method = "RK4", Bcons = blim, track = TRUE, conv.criterion = 0.01 )
This results in a new set of $\mathbf{F_{Nash}}$ values in which, as expected, the fishing pressure on Cod needs to be tapered from $0.204 \text{ year}^{−1}$ to $0.150 \text{ year}^{−1}$ with its consequent reduction in yield from $\approx 15 \times 10^4 \text{ Tn year}^{−1}$ to $\approx 13 \times 10^4 \text{ Tn year}^{−1}$. Application of the constraint increased $F_{\text{Nash,Herring}}$ from $0.277 \text{ year}^{−1}$ to $0.392 \text{ year}^{−1}$ generating a surplus of $\approx 17 \times 10^4 \text{ Tn year}^{−1}$ but had little influence on the NE-MSY harvesting rates for Flounder. For Sprat, by contrast, the conservation constraint on Cod led to a higher $F_{\text{Nash,Sprat}}$ of $0.660 \text{ year}^{−1}$ with a marginal reduction in yield of $\approx 1 \times 10^3 \text{ Tn year}^{−1}$.
``` {r, echo = FALSE, warning = FALSE, message=FALSE, out.width="100%"}
Load data for plotting (Figure 2 in EAP manuscript)
nash.eq.LV.BS <- readRDS("BS-NE-Fmsy-LV.rds") BSCCdataY2 <- readRDS("CCYield_Spp_2.rds") BSCCdataY3 <- readRDS("CCYield_Spp_3.rds") BSCCdataY4 <- readRDS("CCYield_Spp_4.rds") nash.eq.LV.BS.CC <- readRDS("BS-CCNE-Fmsy.rds")
Plotting -----------------------------------------------------------------
library(ggplot2) library(extrafont) library(reshape2) loadfonts()
Baltic Sea
sppname <- c("Cod", "Herring", "Sprat", "Flounder")
Transform data into something ggplot2 reads
Fvec <- c() Yvec <- c() sppvec <- c() for (i in 1:length(nash.eq.LV.BS$YieldEQ)) { Fvec <- append(Fvec, nash.eq.LV.BS$YieldEQ[[i]][, 1]) Yvec <- append(Yvec, nash.eq.LV.BS$YieldEQ[[i]][, 2]) sppvec <- append(sppvec, rep(sppname[i], nrow(nash.eq.LV.BS$YieldEQ[[i]]))) } Yeq <- data.frame("Spp" = sppvec, "Fval" = Fvec, "Yield" = Yvec) Fnash <- data.frame("Spp" = sppname, "Fnash" = as.numeric(t(tail( na.omit(nash.eq.LV.BS$par), n = 1 ))))
Code is kept at the conservation constraint biomass so populate a matrix with
NA for plotting purposes.
BSCCdataY2[[1]] <- matrix(rep(NA, 100), ncol = 2, nrow = 50) BSCCdataY2[[3]] <- BSCCdataY3[[3]] BSCCdataY2[[4]] <- BSCCdataY4[[4]]
Fvec.cons5 <- c() Yvec.cons5 <- c() sppvec.cons5 <- c() for (i in 1:length(BSCCdataY2)) { Fvec.cons5 <- append(Fvec.cons5, BSCCdataY2[[i]][, 1]) Yvec.cons5 <- append(Yvec.cons5, BSCCdataY2[[i]][, 2]) sppvec.cons5 <- append(sppvec.cons5, rep(sppname[i], nrow(BSCCdataY2[[i]]))) } Yeq.cons5 <- data.frame("Spp" = sppvec.cons5, "Fval" = Fvec.cons5, "Yield" = Yvec.cons5) Fnash.cons5 <- data.frame("Spp" = sppname, "Fnash.cons5" = as.numeric(t(tail( na.omit(nash.eq.LV.BS.CC$par), n = 1 )))) Fnash.cons5$Yield <- c(nash.eq.LV.BS.CC$value[1], rep(NA, 3)) lab <- round(Fnash.cons5$Fnash, digits = 3) lab <- c(expression(paste("CC-F"["Nash"], " = ", 0.150)), expression(paste("CC-F"["Nash"], " = ", 0.392)), expression(paste("CC-F"["Nash"], " = ", 0.660)), expression(paste("CC-F"["Nash"], " = ", 0.309))) Fnash.cons5$Labs <- as.character(lab)
lab <- round(Fnash$Fnash, digits = 3) lab <- c(expression(paste("F"["Nash"], " = ", 0.205)), expression(paste("F"["Nash"], " = ", 0.277)), expression(paste("F"["Nash"], " = ", 0.606)), expression(paste("F"["Nash"], " = ", 0.309))) Fnash$Labs <- as.character(lab)
Personalised theme
theme_tjdso <- theme( text = element_text(family = "Consolas", size = 20), panel.background = element_blank(), panel.border = element_rect(fill = FALSE), panel.grid.major = element_line( linetype = 3, colour = "grey", size = 0.25 ), panel.grid.minor = element_line( linetype = 3, colour = "grey", size = 0.1 ), strip.background = element_blank(), strip.text = element_text(size = 20), strip.placement = "outside", panel.spacing.x = unit(5, "mm"), axis.ticks.length = unit(-2, "mm"), axis.text.x.top = element_blank(), axis.text.y.right = element_blank(), axis.title.x.top = element_blank(), axis.title.y.right = element_blank(), axis.title = element_text(size = 18), axis.text.x.bottom = element_text(margin = margin(4, 0, 0, 0, "mm")), axis.text.y.left = element_text(margin = margin(0, 4, 0, 0, "mm")) )
Plot
ggplot(data = Yeq, aes(x = Fval, y = Yield)) + geom_line(size = 1.5, col = 2) + geom_vline( data = Fnash, aes(xintercept = Fnash), linetype = "dashed", size = 1.25, col = 2 ) + geom_text( data = Fnash, aes( x = Fnash, y = c(0.1, 0.5, 0.1, 0.05), label = Labs ), family = "Consolas", parse = TRUE, size = 5, col = 2 ) + facet_wrap( ~ Spp, scales = "free") + scale_x_continuous(labels = scales::number_format(accuracy = 0.01), sec.axis = dup_axis()) + scale_y_continuous(labels = scales::number_format(accuracy = 0.01), sec.axis = dup_axis()) + theme_tjdso + labs(y = bquote(paste("Yield (Kg x ", 10 ^ 3, " x ", Km ^ -2, ")")), x = bquote(paste("Fishing Mortality (", yr ^ -1, ")"))) + # geom_line(data = Yeq.cons, aes(x = Fval, y = Yield), size = 1.5, col = "blue") + geom_line( data = Yeq.cons5, aes(x = Fval, y = Yield), size = 1.5, col = "steelblue" ) + # geom_vline(data = Fnash.cons, aes(xintercept = Fnash.cons), # linetype = "dashed", size = 1.25, col = "blue") + geom_vline( data = Fnash.cons5, aes(xintercept = Fnash.cons5), linetype = "dashed", size = 1.25, col = "steelblue" ) + geom_point( data = Fnash.cons5, aes(x = Fnash.cons5, y = Yield), col = "steelblue", size = 4 ) + geom_text( data = Fnash.cons5, aes( x = Fnash.cons5, y = c(0.25, 1.5, 0.3, 0.15), label = Labs ), family = "Consolas", parse = TRUE, size = 5, col = "steelblue" )
# `nash`'s `LV` algorithm *vs* `round-robin` A key advantage of the `LV` method over the simple `round-robin` method, that is, sequential optimisation of each $F_i$ until convergence is reached, is that it requires much less computation time, especially for complex ecological communities. As a platform-independent metric of performance, we compared the number of objective function evaluations between the two methods. For this analysis, we applied both methods across the three tested models with $n=77$ different initial harvesting rate values. ```r SMFmat <- readRDS("SMFmat.rds") BSFmat <- readRDS("BSFmat.rds") NSFmat <- readRDS("NSFmat.rds") # Theme theme_tjdso <- theme( text = element_text(family = "Consolas", size = 20), panel.background = element_blank(), panel.border = element_rect(fill = FALSE), panel.grid.major = element_line( linetype = 3, colour = "grey", size = 0.25 ), panel.grid.minor = element_line( linetype = 3, colour = "grey", size = 0.1 ), strip.background = element_blank(), strip.text = element_text(size = 20), strip.placement = "outside", panel.spacing.x = unit(5, "mm"), axis.ticks.length = unit(-2, "mm"), axis.text.x.top = element_blank(), axis.text.y.right = element_blank(), axis.title.x.top = element_blank(), axis.title.y.right = element_blank(), axis.title = element_text(size = 18), axis.text.x.bottom = element_text(margin = margin(4, 0, 0, 0, "mm")), axis.text.y.left = element_text(margin = margin(0, 4, 0, 0, "mm")) ) # SM SMF <- reshape2::melt(SMFmat)[, -1] SMF$Var2 <- rep(c("Spp.1", "Spp.2"), each = 29) colnames(SMF) <- c("Species", "Fini") SMF <- as.data.frame(SMF) SMref <- data.frame(Species = c("Spp.1", "Spp.2"), Ref = c(0.2, 0.3)) ggplot(data = SMF, aes(x = Species, y = Fini)) + geom_jitter(col = "grey75", alpha = 0.75) + geom_boxplot(size = 1.2, col = "black", alpha = 0.001) + labs(x = "", y = "") + theme_tjdso + geom_point(data = SMref, aes(x = Species, y = Ref), col = 2) + labs(x = "", y = bquote(paste("Fishing Mortality (", yr ^ -1, ")"))) # BS colnames(BSFmat) <- c("Cod", "Herring", "Sprat", "Flounder") BSF <- reshape2::melt(BSFmat)[, -1] colnames(BSF) <- c("Species", "Fini") BSF <- as.data.frame(BSF) BSref <- data.frame( Species = c("Cod", "Herring", "Sprat", "Flounder"), Ref = c(0.4342084, 0.06663042, 0.1139571, 0.3030173) ) ggplot(data = BSF, aes(x = Species, y = Fini)) + geom_jitter(col = "grey75", alpha = 0.75) + geom_boxplot(size = 1.2, col = "black", alpha = 0.001) + labs(x = "", y = bquote(paste("Fishing Mortality (", yr ^ -1, ")"))) + theme_tjdso + geom_point(data = BSref, aes(x = Species, y = Ref), col = 2) # NS colnames(NSFmat) <- c( "Cod", "Whiting", "Haddock", "Saithe", "Herring", "NorwayPout", "Sandeels", "Plaice", "Sole" ) NSF <- reshape2::melt(NSFmat)[, -1] colnames(NSF) <- c("Species", "Fini") NSF <- as.data.frame(NSF) NSref <- data.frame( Species = c( "Cod", "Whiting", "Haddock", "Saithe", "Herring", "NorwayPout", "Sandeels", "Plaice", "Sole" ), Ref = c( 0.2433686, 0.3393892, 0.1778537, 0.1137075, 0.16818, 0.108255, 0.1468737, 0.2017529, 0.1575716 ) ) ggplot(data = NSF, aes(x = Species, y = Fini)) + geom_jitter(col = "grey75", alpha = 0.75) + geom_boxplot(size = 1.2, col = "black", alpha = 0.001) + labs(x = "", y = "") + theme_tjdso + scale_x_discrete(guide = guide_axis(angle = 45)) + geom_point(data = NSref, aes(x = Species, y = Ref), col = 2) + labs(x = "", y = bquote(paste("Fishing Mortality (", yr ^ -1, ")")))
In all cases, the LV method outperformed the round-robin method for the same convergence threshold (set to a value of $0.01$ for this analysis). Noteworthy is, in particular, that for the North Sea model,
the LV method requires only $1.88$ times more function evaluations than for the Baltic Sea model, even though the dimension of the $\mathbf{F}_\textbf{Nash}$ vector increases by $2.25$ from $4$ to $9$. Using the round-robin method, by comparison, the number of function evaluations increases $7$-fold.
SMdata <- readRDS("SM_fn.counts_data.rds") BSdata <- readRDS("BS_fn.counts_data.rds") NSdata <- readRDS("NS_fn.counts_data.rds") performance <- data.frame(Simple.Model = c(colMeans(SMdata)), BalticSea.Model = c(colMeans(BSdata)), NorthSea.Model = c(colMeans(NSdata, na.rm = TRUE))) colnames(performance) <- c("Simple Model", "Baltic Sea Model", "North Sea Model") rownames(performance) <- c("LV", "round-robin") performance <- as.matrix(performance) barplot(performance, beside = TRUE, ylim = c(0,2000), legend = rownames(performance), args.legend = list(x = "topleft", bty = "n", title = "Method", title.adj = 0.1), ylab = "Function Evaluations (no.)", cex.lab = 1, cex.names = 1, cex.axis = 1, border = TRUE, col = c("gray45","gray85")) # width.error <- 0.025 segments(x0 = 1.5, y0 = mean(na.omit(SMdata[,1]))-sd(na.omit(SMdata[,1])), x1 = 1.5, y1 = mean(na.omit(SMdata[,1]))+sd(na.omit(SMdata[,1])), lwd = 2) # segments(x0 = 1.5-width.error, y0 = mean(SMdata[,1])+sd(SMdata[,1]), # x1 = 1.5+width.error, y1 = mean(SMdata[,1])+sd(SMdata[,1]), lwd = 2) # segments(x0 = 1.5-width.error, y0 = mean(SMdata[,1])-sd(SMdata[,1]), # x1 = 1.5+width.error, y1 = mean(SMdata[,1])-sd(SMdata[,1]), lwd = 2) segments(x0 = 2.5, y0 = mean(na.omit(SMdata[,2]))-sd(na.omit(SMdata[,2])), x1 = 2.5, y1 = mean(na.omit(SMdata[,2]))+sd(na.omit(SMdata[,2])), lwd = 2) # segments(x0 = 2.5-width.error, y0 = mean(SMdata[,2])+sd(SMdata[,2]), # x1 = 2.5+width.error, y1 = mean(SMdata[,2])+sd(SMdata[,2]), lwd = 2) # segments(x0 = 2.5-width.error, y0 = mean(SMdata[,2])-sd(SMdata[,2]), # x1 = 2.5+width.error, y1 = mean(SMdata[,2])-sd(SMdata[,2]), lwd = 2) segments(x0 = 4.5, y0 = mean(na.omit(BSdata[,1]))-sd(na.omit(BSdata[,1])), x1 = 4.5, y1 = mean(na.omit(BSdata[,1]))+sd(na.omit(BSdata[,1])), lwd = 2) # segments(x0 = 4.5-width.error, y0 = mean(na.omit(BSdata[,1]))+sd(na.omit(BSdata[,1])), # x1 = 4.5+width.error, y1 = mean(na.omit(BSdata[,1]))+sd(na.omit(BSdata[,1])), # lwd = 2) # segments(x0 = 4.5-width.error, y0 = mean(na.omit(BSdata[,1]))-sd(na.omit(BSdata[,1])), # x1 = 4.5+width.error, y1 = mean(na.omit(BSdata[,1]))-sd(na.omit(BSdata[,1])), # lwd = 2) segments(x0 = 5.5, y0 = mean(na.omit(BSdata[,2]))-sd(na.omit(BSdata[,2])), x1 = 5.5, y1 = mean(na.omit(BSdata[,2]))+sd(na.omit(BSdata[,2])), lwd = 2) # segments(x0 = 5.5-width.error, y0 = mean(na.omit(BSdata[,2]))+sd(na.omit(BSdata[,2])), # x1 = 5.5+width.error, y1 = mean(na.omit(BSdata[,2]))+sd(na.omit(BSdata[,2])), # lwd = 2) # segments(x0 = 5.5-width.error, y0 = mean(na.omit(BSdata[,2]))-sd(na.omit(BSdata[,2])), # x1 = 5.5+width.error, y1 = mean(na.omit(BSdata[,2]))-sd(na.omit(BSdata[,2])), # lwd = 2) segments(x0 = 7.5, y0 = mean(na.omit(NSdata[,1]))-sd(na.omit(NSdata[,1])), x1 = 7.5, y1 = mean(na.omit(NSdata[,1]))+sd(na.omit(NSdata[,1])), lwd = 2) # segments(x0 = 7.5-width.error, y0 = mean(na.omit(NSdata[,1]))+sd(na.omit(NSdata[,1])), # x1 = 7.5+width.error, y1 = mean(na.omit(NSdata[,1]))+sd(na.omit(NSdata[,1])), # lwd = 2) # segments(x0 = 7.5-width.error, y0 = mean(na.omit(NSdata[,1]))-sd(na.omit(NSdata[,1])), # x1 = 7.5+width.error, y1 = mean(na.omit(NSdata[,1]))-sd(na.omit(NSdata[,1])), # lwd = 2) segments(x0 = 8.5, y0 = mean(na.omit(NSdata[,2]))-sd(na.omit(NSdata[,2])), x1 = 8.5, y1 = mean(na.omit(NSdata[,2]))+sd(na.omit(NSdata[,2])), lwd = 2) grid() box()
References
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.