---
layout: 12pt,review
linenumbers: true
title: "Avoiding tipping points in fisheries management through Gaussian Process Dynamic Programming"
author:
- name: Carl Boettiger
affiliation: cstar
email: cboettig(at)gmail.com
footnote: Corresponding author
- name: Marc Mangel
affiliation: cstar
- name: Stephan Munch
affiliation: noaa
address:
- code: cstar
address: |
Center for Stock Assessment Research,
Department of Applied Math and Statistics,
University of California, Mail Stop SOE-2,
Santa Cruz, CA 95064, USA
- code: marc
address: |
Center for Stock Assessment Research, Department of Applied Math
and Statistics, University of California,
Mail Stop SOE-2, Santa Cruz, CA 95064, USA and Department of
Biology, University of Bergen, Bergen, Norway 9020
- code: noaa
address: |
Southwest Fisheries Science Center,
National Oceanic and Atmospheric Administration,
110 Shaffer Road, Santa Cruz, CA 95060, USA
abstract: |
Model uncertainty and limited data are fundamental challenges to
robust management of human intervention in a natural system.
These challenges are acutely highlighted by concerns that
many ecological systems may contain tipping points, such as
Allee population sizes. Before a collapse, we do not know
where the tipping points lie, if they exist at all. Hence,
we know neither a complete model of the system dynamics nor
do we have access to data in some large region of state-space
where such a tipping point might exist. We illustrate how
a Bayesian Non-Parametric (BNP) approach using a Gaussian
Process (GP) prior provides a flexible representation of this
inherent uncertainty. We embed GPs in a Stochastic Dynamic
Programming (SDP) framework in order to make robust management
predictions with both model uncertainty and limited data. We
use simulations to evaluate this approach as compared with the
standard approach of using model selection to choose from a set
of candidate models. We find that model selection erroneously
favors models without tipping points -- leading to harvest
policies that guarantee extinction. The GPDP performs nearly as
well as the true model and significantly outperforms standard
approaches. We illustrate this using examples of simulated
single-species dynamics, where the standard model selection
approach should be most effective, and find that it still
fails to account for uncertainty appropriately and leads to
population crashes, while management based on the GPDP does
not, since it does not underestimate the uncertainty outside
of the observed data.
keywords:
- Bayesian
- Structural Uncertainty
- Nonparametric
- Optimal Control
- Decision Theory
- Gaussian Processes
- Fisheries Management
bibliography: components/references.bib
csl: components/ecology.csl
documentclass: components/elsarticle
## rmarkdown render options
output:
pdf_document:
fig_caption: true
template: components/elsarticle.latex
keep_tex: true
---
```{r caching, include=FALSE}
library("methods")
library("knitr")
basename <- "manuscript"
opts_chunk$set(fig.path = paste("components/figure/", basename, "-", sep=""),
cache.path = paste("components/cache/", basename, "/", sep=""))
opts_chunk$set(cache = 2)
opts_chunk$set(tidy=FALSE, warning=FALSE, message=FALSE,
comment = NA, verbose = TRUE, echo=FALSE)
# PDF-based figures
opts_chunk$set(dev='pdf')
```
```{r plot-options, message=FALSE, warning=FALSE, include=FALSE, echo=FALSE}
cbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73",
"#F0E442", "#0072B2", "#D55E00", "#CC79A7")
opts_chunk$set(dev='pdf', dev.args=c(version="1.3"))
```
```{r libraries}
library("nonparametricbayes")
library("pdgControl")
library("ggplot2")
library("reshape2")
library("plyr")
library("data.table")
library("R2jags")
library("emdbook") # for as.mcmc.bugs (?)
library("coda") # for as.mcmc
library("modeest")
library("MASS")
library("pander")
library("cboettigR")
opts_chunk$set(fig.width=4, fig.height=3, echo=FALSE)
theme_set(theme_bw(base_size=12))
toggle = "hide" # results argument
```
```{r posterior-mode}
library(modeest)
posterior.mode <- function(x) {
mlv(x, method="shorth")$M
}
```
```{r stateeq}
f <- RickerAllee
p <- c(2, 8, 5)
K <- 10 # approx, a li'l' less
allee <- 5 # approx, a li'l' less
```
```{r sdp-pars, dependson="stateeq"}
sigma_g <- 0.05
sigma_m <- 0.0
z_g <- function() rlnorm(1, 0, sigma_g)
z_m <- function() 1
x_grid <- seq(0, 1.5 * K, length=50)
h_grid <- x_grid
profit <- function(x,h) pmin(x, h)
delta <- 0.01
OptTime <- 50 # stationarity with unstable models is tricky thing
reward <- 0
xT <- 0
Xo <- allee+.5# observations start from
x0 <- K # simulation under policy starts from
Tobs <- 40
MaxT = 1000 # timeout for value iteration convergence
```
```{r obs, dependson="sdp-pars", fig.keep='none'}
set.seed(1234)
#harvest <- sort(rep(seq(0, .5, length=7), 5))
x <- numeric(Tobs)
x[1] <- Xo
nz <- 1
for(t in 1:(Tobs-1))
x[t+1] = z_g() * f(x[t], h=0, p=p)
obs <- data.frame(x = c(rep(0,nz),
pmax(rep(0,Tobs-1), x[1:(Tobs-1)])),
y = c(rep(0,nz),
x[2:Tobs]))
raw_plot <- ggplot(data.frame(time = 1:Tobs, x=x), aes(time,x)) + geom_line()
```
```{r mle, dependson="obs"}
set.seed(12345)
estf <- function(p){
mu <- f(obs$x,0,p)
-sum(dlnorm(obs$y, log(mu), p[4]), log=TRUE)
}
par <- c(3,
9,
6,
0.06)
o <- optim(par, estf, method="L", lower=c(1e-5,1e-5,1e-5,1e-5))
f_alt <- f
p_alt <- c(as.numeric(o$par[1]), as.numeric(o$par[2]), as.numeric(o$par[3]))
sigma_g_alt <- as.numeric(o$par[4])
est <- list(f = f_alt, p = p_alt, sigma_g = sigma_g_alt, mloglik=o$value)
```
```{r mle-output, dependson="mle", results=toggle}
true_means <- sapply(x_grid, f, 0, p)
est_means <- sapply(x_grid, est$f, 0, est$p)
```
```{r gp-priors}
s2.p <- c(5,5)
d.p = c(10, 1/0.1)
```
```{r gp, dependson=c("gp-priors", "obs")}
gp <- gp_mcmc(obs$x, y=obs$y, n=1e5, s2.p = s2.p, d.p = d.p)
gp_dat <- gp_predict(gp, x_grid, burnin=1e4, thin=300)
```
```{r gp_traces_densities, dependson="gp"}
gp_assessment_plots <- summary_gp_mcmc(gp, burnin=1e4, thin=300)
```
```{r gp-output, dependson="gp", results=toggle}
tgp_dat <-
data.frame( x = x_grid,
y = gp_dat$E_Ef,
ymin = gp_dat$E_Ef - 2 * sqrt(gp_dat$E_Vf),
ymax = gp_dat$E_Ef + 2 * sqrt(gp_dat$E_Vf) )
```
```{r jags-setup}
y <- x
N <- length(x);
jags.data <- list("N","y")
n.chains <- 6
n.iter <- 1e6
n.burnin <- floor(10000)
n.thin <- max(1, floor(n.chains * (n.iter - n.burnin)/1000))
n.update <- 10
```
```{r common-priors}
stdQ_prior_p <- c(1e-6, 100)
stdR_prior_p <- c(1e-6, .1)
stdQ_prior <- function(x) dunif(x, stdQ_prior_p[1], stdQ_prior_p[2])
stdR_prior <- function(x) dunif(x, stdR_prior_p[1], stdR_prior_p[2])
```
```{r allen-model}
K_prior_p <- c(0.01, 20.0)
r0_prior_p <- c(0.01, 6.0)
theta_prior_p <- c(0.01, 20.0)
bugs.model <-
paste(sprintf(
"model{
K ~ dunif(%s, %s)
r0 ~ dunif(%s, %s)
theta ~ dunif(%s, %s)
stdQ ~ dunif(%s, %s)",
K_prior_p[1], K_prior_p[2],
r0_prior_p[1], r0_prior_p[2],
theta_prior_p[1], theta_prior_p[2],
stdQ_prior_p[1], stdQ_prior_p[2]),
"
iQ <- 1 / (stdQ * stdQ);
y[1] ~ dunif(0, 10)
for(t in 1:(N-1)){
mu[t] <- log(y[t]) + r0 * (1 - y[t]/K)* (y[t] - theta) / K
y[t+1] ~ dlnorm(mu[t], iQ)
}
}")
writeLines(bugs.model, "allen_process.bugs")
```
```{r allen-priors, dependson="common-priors"}
K_prior <- function(x) dunif(x, K_prior_p[1], K_prior_p[2])
r0_prior <- function(x) dunif(x, r0_prior_p[1], r0_prior_p[2])
theta_prior <- function(x) dunif(x, theta_prior_p[1], theta_prior_p[2])
par_priors <- list(K = K_prior, deviance = function(x) 0 * x,
r0 = r0_prior, theta = theta_prior,
stdQ = stdQ_prior)
```
```{r allen_priors_xtable, include=TRUE, echo=FALSE}
allen_priors_xtable <- data.frame(parameter = c("$r$", "$K$", "$X_C$", "$\\sigma$"),
"lower bound" = c(r0_prior_p[1], K_prior_p[1], theta_prior_p[1], stdQ_prior_p[1]),
"upper bound" = c(r0_prior_p[2], K_prior_p[2], theta_prior_p[2], stdQ_prior_p[2]))
```
```{r allen-mcmc, dependson=c("allen-model", "allen-pars", "jags-setup"), results=toggle}
jags.params=c("K","r0","theta","stdQ") # be sensible about the order here
jags.inits <- function(){
list("K"= 10 * rlnorm(1,0, 0.1),
"r0"= 1 * rlnorm(1,0, 0.1) ,
"theta"= 5 * rlnorm(1,0, 0.1) ,
"stdQ"= abs( 0.1 * rlnorm(1,0, 0.1)),
.RNG.name="base::Wichmann-Hill", .RNG.seed=123)
}
set.seed(1234)
# parallel refuses to take variables as arguments (e.g. n.iter = 1e5 works, but n.iter = n doesn't)
allen_jags <- do.call(jags, list(data=jags.data, inits=jags.inits,
jags.params, n.chains=n.chains,
n.iter=n.iter, n.thin=n.thin,
n.burnin=n.burnin,
model.file="allen_process.bugs"))
# Run again iteratively if we haven't met the Gelman-Rubin convergence criterion
recompile(allen_jags) # required for parallel
allen_jags <- do.call(autojags,
list(object=allen_jags, n.update=n.update,
n.iter=n.iter, n.thin = n.thin))
```
```{r allen-traces, dependson="allen-mcmc"}
tmp <- lapply(as.mcmc(allen_jags), as.matrix) # strip classes the hard way...
allen_posteriors <- melt(tmp, id = colnames(tmp[[1]]))
names(allen_posteriors) = c("index", "variable", "value", "chain")
plot_allen_traces <- ggplot(allen_posteriors) + geom_line(aes(index, value)) +
facet_wrap(~ variable, scale="free", ncol=1)
```
```{r allen-posteriors, dependson=c("allen-traces", "allen-priors")}
allen_priors <- ddply(allen_posteriors, "variable", function(dd){
grid <- seq(min(dd$value), max(dd$value), length = 100)
data.frame(value = grid, density = par_priors[[dd$variable[1]]](grid))
})
plot_allen_posteriors <- ggplot(allen_posteriors, aes(value)) +
stat_density(geom="path", position="identity") +
# geom_line(data=allen_priors, aes(x=value, y=density), col="red") +
facet_wrap(~ variable, scale="free", ncol=3)
```
```{r allen-output, dependson=c("posterior-mode", "allen-traces"), results=toggle}
A <- allen_posteriors
A$index <- A$index + A$chain * max(A$index) # Combine samples across chains by renumbering index
pardist <- acast(A, index ~ variable)
bayes_coef <- apply(pardist,2, posterior.mode)
bayes_pars <- unname(c(bayes_coef["r0"], bayes_coef["K"], bayes_coef["theta"])) # parameters formatted for f
allen_f <- function(x,h,p) unname(RickerAllee(x,h, unname(p[c("r0", "K", "theta")])))
allen_means <- sapply(x_grid, f, 0, bayes_pars)
bayes_pars
head(pardist)
```
```{r ricker-model}
K_prior_p <- c(0.01, 40.0)
r0_prior_p <- c(0.01, 20.0)
bugs.model <-
paste(sprintf(
"model{
K ~ dunif(%s, %s)
r0 ~ dunif(%s, %s)
stdQ ~ dunif(%s, %s)",
K_prior_p[1], K_prior_p[2],
r0_prior_p[1], r0_prior_p[2],
stdQ_prior_p[1], stdQ_prior_p[2]),
"
iQ <- 1 / (stdQ * stdQ);
y[1] ~ dunif(0, 10)
for(t in 1:(N-1)){
mu[t] <- log(y[t]) + r0 * (1 - y[t]/K)
y[t+1] ~ dlnorm(mu[t], iQ)
}
}")
writeLines(bugs.model, "ricker_process.bugs")
```
```{r ricker-priors, dependson="common-priors"}
K_prior <- function(x) dunif(x, K_prior_p[1], K_prior_p[2])
r0_prior <- function(x) dunif(x, r0_prior_p[1], r0_prior_p[2])
par_priors <- list(K = K_prior, deviance = function(x) 0 * x,
r0 = r0_prior, stdQ = stdQ_prior)
```
```{r ricker_priors_xtable, include=TRUE, echo=FALSE}
ricker_priors_xtable <- data.frame(
parameter = c("$r$", "$K$", "$\\sigma$"),
"lower bound" = c(r0_prior_p[1], K_prior_p[1], stdQ_prior_p[1]),
"upper bound" = c(r0_prior_p[2], K_prior_p[2], stdQ_prior_p[2]))
```
```{r ricker-mcmc, dependson="ricker-model", results=toggle}
jags.params=c("K","r0", "stdQ")
jags.inits <- function(){
list("K"= 10 * rlnorm(1,0,.5),
"r0"= rlnorm(1,0,.5),
"stdQ"=sqrt(0.05) * rlnorm(1,0,.5),
.RNG.name="base::Wichmann-Hill", .RNG.seed=123)
}
set.seed(12345)
ricker_jags <- do.call(jags,
list(data=jags.data, inits=jags.inits,
jags.params, n.chains=n.chains,
n.iter=n.iter, n.thin=n.thin, n.burnin=n.burnin,
model.file="ricker_process.bugs"))
recompile(ricker_jags)
ricker_jags <- do.call(autojags,
list(object=ricker_jags, n.update=n.update,
n.iter=n.iter, n.thin = n.thin,
progress.bar="none"))
```
```{r ricker-traces, dependson="ricker-mcmc"}
tmp <- lapply(as.mcmc(ricker_jags), as.matrix) # strip classes the hard way...
ricker_posteriors <- melt(tmp, id = colnames(tmp[[1]]))
names(ricker_posteriors) = c("index", "variable", "value", "chain")
plot_ricker_traces <- ggplot(ricker_posteriors) + geom_line(aes(index, value)) +
facet_wrap(~ variable, scale="free", ncol=1)
```
```{r ricker-posteriors, dependson=c("ricker-traces", "ricker-priors")}
ricker_priors <- ddply(ricker_posteriors, "variable", function(dd){
grid <- seq(min(dd$value), max(dd$value), length = 100)
data.frame(value = grid, density = par_priors[[dd$variable[1]]](grid))
})
# plot posterior distributions
plot_ricker_posteriors <- ggplot(ricker_posteriors, aes(value)) +
stat_density(geom="path", position="identity") +
# geom_line(data=ricker_priors, aes(x=value, y=density), col="red") + # don't plot priors
facet_wrap(~ variable, scale="free", ncol=2)
```
```{r ricker-output, dependson=c("posterior-mode", "ricker-traces"), results=toggle}
A <- ricker_posteriors
A$index <- A$index + A$chain * max(A$index) # Combine samples across chains by renumbering index
ricker_pardist <- acast(A, index ~ variable)
bayes_coef <- apply(ricker_pardist,2, posterior.mode)
ricker_bayes_pars <- unname(c(bayes_coef["r0"], bayes_coef["K"]))
ricker_f <- function(x,h,p){
sapply(x, function(x){
x <- pmax(0, x-h)
pmax(0, x * exp(p["r0"] * (1 - x / p["K"] )) )
})
}
ricker_means <- sapply(x_grid, Ricker, 0, ricker_bayes_pars[c(1,2)])
head(ricker_pardist)
ricker_bayes_pars
```
```{r myers-model}
r0_prior_p <- c(.0001, 10.0)
theta_prior_p <- c(.0001, 10.0)
K_prior_p <- c(.0001, 40.0)
bugs.model <-
paste(sprintf(
"model{
r0 ~ dunif(%s, %s)
theta ~ dunif(%s, %s)
K ~ dunif(%s, %s)
stdQ ~ dunif(%s, %s)",
r0_prior_p[1], r0_prior_p[2],
theta_prior_p[1], theta_prior_p[2],
K_prior_p[1], K_prior_p[2],
stdQ_prior_p[1], stdQ_prior_p[2]),
"
iQ <- 1 / (stdQ * stdQ);
y[1] ~ dunif(0, 10)
for(t in 1:(N-1)){
mu[t] <- log(r0) + theta * log(y[t]) - log(1 + pow(abs(y[t]), theta) / K)
y[t+1] ~ dlnorm(mu[t], iQ)
}
}")
writeLines(bugs.model, "myers_process.bugs")
```
```{r myers-priors}
K_prior <- function(x) dunif(x, K_prior_p[1], K_prior_p[2])
r_prior <- function(x) dunif(x, r0_prior_p[1], r0_prior_p[2])
theta_prior <- function(x) dunif(x, theta_prior_p[1], theta_prior_p[2])
par_priors <- list( deviance = function(x) 0 * x, K = K_prior,
r0 = r_prior, theta = theta_prior,
stdQ = stdQ_prior)
```
```{r myers_prior_xtable, include=TRUE, echo=FALSE}
myers_priors_xtable <- data.frame(parameter = c("$r$", "$K$", "$\\theta$", "$\\sigma$"),
"lower bound" = c(r0_prior_p[1], K_prior_p[1], theta_prior_p[1], stdQ_prior_p[1]),
"upper bound" = c(r0_prior_p[2], K_prior_p[2], theta_prior_p[2], stdQ_prior_p[2]))
```
```{r myers-mcmc, dependson="myers-model", results=toggle}
jags.params=c("r0", "theta", "K", "stdQ")
jags.inits <- function(){
list("r0"= 1 * rlnorm(1,0,.1),
"K"= 10 * rlnorm(1,0,.1),
"theta" = 1 * rlnorm(1,0,.1),
"stdQ"= sqrt(0.2) * rlnorm(1,0,.1),
.RNG.name="base::Wichmann-Hill", .RNG.seed=123)
}
set.seed(12345)
myers_jags <- do.call(jags,
list(data=jags.data, inits=jags.inits,
jags.params, n.chains=n.chains,
n.iter=n.iter, n.thin=n.thin,
n.burnin=n.burnin,
model.file="myers_process.bugs"))
recompile(myers_jags)
myers_jags <- do.call(autojags,
list(myers_jags, n.update=n.update,
n.iter=n.iter, n.thin = n.thin,
progress.bar="none"))
```
```{r myers-traces, dependson="myers-mcmc"}
tmp <- lapply(as.mcmc(myers_jags), as.matrix) # strip classes
myers_posteriors <- melt(tmp, id = colnames(tmp[[1]]))
names(myers_posteriors) = c("index", "variable", "value", "chain")
plot_myers_traces <- ggplot(myers_posteriors) +
geom_line(aes(index, value)) + # priors, need to fix order though
facet_wrap(~ variable, scale="free", ncol=1)
```
```{r myers-posteriors, dependson="myers-traces"}
par_prior_curves <- ddply(myers_posteriors, "variable", function(dd){
grid <- seq(min(dd$value), max(dd$value), length = 100)
data.frame(value = grid, density = par_priors[[dd$variable[1]]](grid))
})
plot_myers_posteriors <- ggplot(myers_posteriors, aes(value)) +
stat_density(geom="path", position="identity") +
# geom_line(data=par_prior_curves, aes(x=value, y=density), col="red") + # Whoops, these are misaligned. see table instead
facet_wrap(~ variable, scale="free", ncol=3)
```
```{r myers-output, dependson=c("posterior-mode", "myers-traces"), results=toggle}
A <- myers_posteriors
A$index <- A$index + A$chain * max(A$index) # Combine samples across chains by renumbering index
myers_pardist <- acast(A, index ~ variable)
bayes_coef <- apply(myers_pardist,2, posterior.mode) # much better estimates
myers_bayes_pars <- unname(c(bayes_coef["r0"], bayes_coef["theta"], bayes_coef["K"]))
myers_means <- sapply(x_grid, Myer_harvest, 0, myers_bayes_pars)
myers_f <- function(x,h,p) Myer_harvest(x, h, p[c("r0", "theta", "K")])
head(myers_pardist)
myers_bayes_pars
```
```{r assemble-models, dependson=c("myers-output", "ricker-output", "allen-output", "gp-output", "mle-output")}
models <- data.frame(x=x_grid,
GP=tgp_dat$y,
True=true_means,
MLE=est_means,
Ricker=ricker_means,
Allen = allen_means,
Myers = myers_means)
models <- melt(models, id="x")
# some labels
names(models) <- c("x", "method", "value")
# labels for the colorkey too
model_names = c("GP", "True", "MLE", "Ricker", "Allen", "Myers")
colorkey=cbPalette
names(colorkey) = model_names
step_ahead <- function(x, f, p){
h = 0
x_predict <- sapply(x, f, h, p)
n <- length(x_predict) - 1
y <- c(x[1], x_predict[1:n])
y
}
```
```{r gp-opt, dependson="gp-output"}
matrices_gp <- gp_transition_matrix(gp_dat$Ef_posterior, gp_dat$Vf_posterior, x_grid, h_grid)
opt_gp <- value_iteration(matrices_gp, x_grid, h_grid, MaxT, xT, profit, delta, reward)
```
```{r mle-opt, dependson="mle-output"}
matrices_true <- f_transition_matrix(f, p, x_grid, h_grid, sigma_g)
opt_true <- value_iteration(matrices_true, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta)
matrices_estimated <- f_transition_matrix(est$f, est$p, x_grid, h_grid, est$sigma_g)
opt_estimated <- value_iteration(matrices_estimated, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta)
```
```{r allen-opt, dependson="allen-output"}
matrices_allen <- parameter_uncertainty_SDP(allen_f, x_grid, h_grid, pardist, 4)
opt_allen <- value_iteration(matrices_allen, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta)
```
```{r ricker-opt, dependson="ricker-output"}
matrices_ricker <- parameter_uncertainty_SDP(ricker_f, x_grid, h_grid, as.matrix(ricker_pardist), 3)
opt_ricker <- value_iteration(matrices_ricker, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta)
```
```{r myers-opt, dependson="myers-output"}
matrices_myers <- parameter_uncertainty_SDP(myers_f, x_grid, h_grid, as.matrix(myers_pardist), 4)
myers_alt <- value_iteration(matrices_myers, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta)
```
```{r assemble-opt, dependson=c("gp-opt", "mle-opt", "allen-opt", "ricker-opt", "myers-opt")}
OPT = data.frame(GP = opt_gp$D, True = opt_true$D, MLE = opt_estimated$D, Ricker = opt_ricker$D, Allen = opt_allen$D, Myers = myers_alt$D)
colorkey=cbPalette
names(colorkey) = names(OPT)
```
```{r sims, dependson="assemble-opt"}
sims <- lapply(OPT, function(D){
set.seed(1)
lapply(1:100, function(i)
ForwardSimulate(f, p, x_grid, h_grid, x0, D, z_g, profit=profit, OptTime=OptTime)
)
})
# turn the list into a data.frame
keep.as.columns <- names(sims[[1]][[1]])
dat <- melt(sims, id=keep.as.columns)
sims_data <- data.table(dat)
setnames(sims_data, c("L1", "L2"), c("method", "reps"))
# Legend in original ordering please, not alphabetical:
sims_data$method = factor(sims_data$method, ordered=TRUE, levels=names(OPT))
```
```{r profits, dependson="sims", results=toggle}
Profit <- sims_data[, sum(profit), by=c("reps", "method")]
tmp <- dcast(Profit, reps ~ method)
#tmp$Allen <- tmp[,"Allen"] + rnorm(dim(tmp)[1], 0, 1) # jitter for plotting
tmp <- tmp / tmp[,"True"]
tmp <- melt(tmp[2:dim(tmp)[2]])
actual_over_optimal <-subset(tmp, variable != "True")
```
```{r dic_calc, dependson=c("posterior-mode", "myers-output", "ricker-output", "allen-output", "gp-output", "mle-output"), include=FALSE, echo=FALSE}
dic.dic <- function (x) sum(x$deviance) + sum(x[[2]])
recompile(allen_jags)
allen_dic <- dic.dic(dic.samples(allen_jags$model, n.iter=1000, type="popt"))
recompile(ricker_jags)
ricker_dic <- dic.dic(dic.samples(ricker_jags$model, n.iter=1000, type="popt"))
recompile(myers_jags)
myers_dic <- dic.dic(dic.samples(myers_jags$model, n.iter=1000, type="popt"))
dictable <- data.frame(Allen = allen_dic + 2*length(bayes_pars),
Ricker = ricker_dic + 2*length(ricker_bayes_pars),
Myers = myers_dic + 2*length(myers_bayes_pars),
row.names = c("DIC"))
```
```{r sensitivity-trends}
source("components/sensitivity-trends.R") # load the function for looping
sigmas <- c(0.01, 0.05, 0.1)
allee <- c(1,2,3,4,5)
vary_sigma <- lapply(sigmas, function(s) try(sensitivity_trends(4, sigma = s, seed=c(111,222,333))))
vary_allee <- lapply(allee, function(a) try(sensitivity_trends(a, sigma = sigma_g, seed=c(111,222,333))))
vary_sigma <- melt(vary_sigma, id=names(vary_sigma[[1]]))
vary_allee <- melt(vary_allee, id=names(vary_allee[[1]]))
```
```{r refresh-plot-options, message=FALSE, warning=FALSE, include=FALSE, echo=FALSE, cache=FALSE}
cbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73",
"#F0E442", "#0072B2", "#D55E00", "#CC79A7")
opts_chunk$set(dev='pdf', dev.args=c(version="1.7"))
opts_chunk$set(fig.width=4, fig.height=3, echo=FALSE)
theme_set(theme_bw(base_size=12))
```
Introduction
============
Decision making under uncertainty is a ubiquitous challenge in the
management of human intervention in natural resources and conservation.
Decision-theoretic approaches provide a framework to determine the
best sequence of actions in face of uncertainty, but only when that
uncertainty can be meaningfully quantified [@Fischer2009]. Over the last
four decades (beginning with
@Clark1976, @Clark2009 and @Walters1978)
dynamic optimization methods, particularly Stochastic Dynamic Programming
(SDP), have become increasingly important as a means of
understanding how to manage human intervention into natural systems.
Simultaneously, there has been increasing recognition of the importance of
multiple steady states or `tipping points' [@Scheffer2001; . @Scheffer2009;
@Polasky2011] in ecological systems.
We develop a novel approach to address these concerns in the context of
fisheries; although the challenges and methods are germane to
other problems of conservation or natural resource exploitation. Economic value and ecological concern have made
marine fisheries the crucible for much of the founding work on management
under uncertainty [@Gordon1954;
@Clark1976; -@Clark2009; @Reed1979; @May1979; @Ludwig1982]
Even if we know the proper deterministic description of the biological system, there is intrinsic stochasticity in biological dynamics, measurements, and
implementation of policy [_e.g._ @Reed1979; @Clark1986; @Roughgarden1996;
@Sethi2005]. We may also lack knowledge about the parameters of the biological dynamics [parametric uncertainty, _e.g._ @Ludwig1982; @Hilborn1997;
@McAllister1998; @Schapaugh2013], or even not know which model is proper description of the system
[structural uncertainty, _e.g._ @Williams2001; @Cressie2009; @Athanassoglou2012]. Of these,
the latter is generally the hardest to quantify. Typical approaches
confront the data with a collection of models, assuming that the true dynamics (or reasonable approximation) is among the collection and then use model choice or model averaging to arrive at a conclusion
[@Williams2001; @Athanassoglou2012; @Cressie2009]. Even setting aside
other concerns (see @Cressie2009), these
approaches are unable to describe uncertainty
outside the observed data range.
Structural uncertainty is particularly insidious when we try to predict
outside of the range of observed data [@Mangel2001] because we are extrapolating into unknown regions. In management applications, this extrapolation uncertainty is particularly important since (a) management involves considering actions that may move the
system outside the range of observed behavior, and (b) the decision tools (optimal control theory, SDP) rely on both reasonable estimates of the expected
outcomes and on the weights given to those outcomes
[_e.g._ @Weitzman2013]. Thus characterizing uncertainty is as
important as characterizing the expected outcome.
Tipping points in ecological dynamics [@Scheffer2001;
@Polasky2011] highlight this problem because
precise models are not available and data are
limited such as around high
stock levels or an otherwise desirable state. With perfect information
, one would know just how far a system could
be pushed before crossing the tipping point, and management would be
simple. But we face imperfect models and limited data and, with tipping points,
even small errors can have
very large consequences, as we shall illustrate later. Because intervention is often (but not always, see
@Hughes2013) too late after a tipping point has been crossed, management
is often concerned with avoiding
tipping points before any data about them is available or the
transition that would more clearly reveal a tipping point occurs
[e.g. @Bestelmeyer2012].
The dual concerns of model uncertainty and incomplete data
create a substantial challenge to existing decision-theoretic approaches
[@Brozovic2011]. We illustrate how Stochastic Dynamic Programming (SDP)
[@Mangel1988; @Marescot2013] can be implemented using a
Bayesian Non-Parametric (BNP) model of population dynamics [@Munch2005a].
The BNP method has two distinct advantages
. First, using a BNP model sidesteps the need for an accurate
model-based description of the system dynamics. Second, a BNP model
better reflects uncertainty when extrapolating beyond the observed
data. This is crucial to providing robust
decision-making when the correct model is not known (as is almost always true).
[We use _robust_ to characterize approaches that provide nearly optimal solutions
without being sensitive to the choice of the (unknown) underlying model.]
This paper is the first ecological application of the
SDP without an _a priori_ model of the
underlying dynamics. Unlike parametric approaches that can only
reflect uncertainty in parameter estimates, the BNP method provides a
broader representation of uncertainty, including uncertainty beyond the observed data.
We will show that Gaussian Process Dynamic Programming
(GPDP) allows us to find robust management solutions in face of limited
data without knowing the correct model structure.
For comparisons, we consider the performance of management based on GPDP against
management policies derived under several alternative parametric models [@Reed1979; @Ludwig1982; @Mangel1988]. Rather than compare models in terms
of best fit to data, we compare model performance in the concrete terms
of the decision-maker's objectives.
Approach and Methods
====================
We first describe the requirements of dynamic optimization for the
management of human intervention in natural resource systems. After that
we describe three parametric models for population dynamics and the Gaussian
Process (GP)[^1] description of population dynamics.
[^1]: We abbreviate Gaussian Process as GP, which refers to
the statistical model we use to approximate the population dynamics
, and we use the
term Gaussian Process Dynamic Programming [GPDP,], to refer to the _use_ of a
GP as the underlying process model when solving a Dynamic Programming equation
. Hence we will refer to the models as: GP, Ricker, Allen, etc, and the novel method we put forward here as GPDP.
### Requirements of dynamic optimization
Dynamic optimization requires characterizing the dynamics of a state
variable (or variables), a control action, and a value function.
For simplicity, we consider only a single state variable, which allows us to
focus on one-dimensional models of population dynamics
for the comparison of parametric models and the GP.
This is a best-case scenario for the parametric models because we
simulate underlying dynamics from one of the three parametric models
, whereas in the natural world we never know the "true" mode
. In addition, by choosing one-dimensional models with just a few parameters, we
limit the chance that poor performance will be due to inability to
estimate parameters accurately, something that becomes a more severe problem
for higher-dimensional parametric models. Finally, the parametric models
we consider are those most commonly used in modeling stock-recruitment
dynamics or to model sudden transitions between alternative stable states.
### Parametric models
We let $X(t)$ denote the size (numbers or biomass) of the focal
population at time $t$ and assume that in the absence of take its
dynamics are.
.
\begin{equation}
X(t+1) = Z(t) f(X(t), p) \label{eq1}
\end{equation}
Where Z(t) is log-normally distributed process stochasticity [@Reed1979]
and $p$ is a vector of parameters to be estimated from the data.
We describe the three choices for $f(X(t),p)$ in the next section.
The control action is a harvest or take, $h(t)$, measured
in the same units as $X$, at time $t$. Thus, in the presence of take,
the population size on the right hand side of Eqn 1 is replaced by
$S(t) = X(t) - h(t)$.
To construct the value function, we consider a return when $X(t) = x(t)$ and
harvest $h(t) = h$ denoted as the reward, $R(x(t), h)$. For example, if the
return is the harvest at time $t$, then $R(x(t), h(t)) = \min(x(t), h(t))$. We
assume that future harvests are discounted relative to current ones at
a constant rate of discount $\delta$ and ask for the harvest policy
that maximizes total discounted harvest between the current time $t$
and a final time $T$. That is, we seek to maximize over choices of
harvest $E_{X(t+1)} [ \sum_{t = 0}^{T} R(X(t), h(t), t) \delta^t]$,
where the state dynamics are given by Eqn 1 and $E$ denotes the
expectation over the process stochasticity in the future population state.
In order to find that policy, we introduce the value function $V(x(t), t)$
representing the total discounted catch from time $t$ onwards
given that $X(t) = x(t)$. This value function satisfies an equation of
SDP [@Mangel1988; @Clark2000; @Clark2009 ;@Mangel2014],
\begin{equation}
V(x(t), t) = \max_{h(t)} \lbrace R(h(t), x(t)) + \delta \cdot \mathbf{\mathrm{E}}_{X(t+1)} \left[ V(X(t+1), t+1) | x(t), h(t) \right] \rbrace
\end{equation}
where expectation is taken over all possible values of the next state,
$X(t+1)$, and maximized over all possible choices of harvest, $h(t)$.
That is, at time $t$, when population size is $x(t)$ and harvest $h(t)$ is
applied, the immediate return is $R(h(t), x(t))$. When
the sole source of uncertainty is the process stochasticity term, $Z$,
and thus the expectation in Eqn 2 is equivalent to taking expectations over
$Z(t)$. That is
\begin{equation}
\mathbf{\mathrm{E}}_{X(t+1)} \left[ V( X(t+1),t+1) | x(t), h(t) \right] = \mathbf{\mathrm{E}}_{Z(t)} \left[V( Z(t) f(x(t) - h(t))|p), t+1 | x(t), h(t) \right]
\end{equation}
where the population size after the take is $x(t) - h(t)$, which is then
translated into $X(t+1)$ by Eqn 1 (that is, we replace $X(t+1)$ by $Z(t) f(x(t)-h(t)|p)$).
When the parameters governing the dynamics
are also uncertain (the case of parametric uncertainty), the expectation in Eqn 2 involves averaging over the
posterior distribution for the parameters. That is,
\begin{equation}
\mathbf{\mathrm{E}}_{X(t+1)} \left[ V(X(t+1), t+1) | x(t), h(t) \right] = \mathbf{\mathrm{E}}_{p|\mathrm{data}} \{ \mathbf{\mathrm{E}}_{Z(t) | p, \mathrm{data}} \left[ V(Z(t) f(x(t) - h(t)|p), t) \right] \}
\end{equation}
When the underlying population dynamics are unknown (the case of structural uncertainty), the function $f$ itself is uncertain and the expectation for the next state includes
uncertainty in $f$ as well. That is
\begin{equation}
\mathbf{\mathrm{E}}_{X(t+1)} \left[ V(X(t+1), t+1) | x(t), h(t) \right] = \mathbf{\mathrm{E}}_{p|\mathrm{data}} \{ \mathbf{\mathrm{E}}_{f, Z(t) | p, \mathrm{data}} \left[ V( Z(t) f(x(t) - h(t)|p), t) \right] \}
\end{equation}
We consider the finite time problem with $T=$ `r MaxT`, which we
solve using the standard value iteration algorithm [see @Mangel1988;
@Clark2000].
\subsubsection{Parametric models}\label{parametric-models}
We consider three candidate parametric models for the population dynamics: The Ricker model, the Allen model [@Allen2005], and the Myers
model [@Myers1995], Eqns \eqref{ricker}-\eqref{myers}. In all three, we let $K$ denote the carrying capacity and $r$ the maximum per capita growth rate. The Ricker model has two parameters and the right hand side of Eqn 1 is
\begin{equation}
f(S(t)|r,K) = S(t) e^{r \left(1 - \frac{S(t)}{K} \right) } \label{ricker}
\end{equation}
The Allen model has three parameters
\begin{equation}
f(S(t)|r, K, X_c) = S(t) e^{r \left(1 - \frac{S(t)}{K}\right)\left(S(t) - X_c\right)} \label{allen}
\end{equation}
where $X_c$ denotes the location of the unstable steady state (i.e., the tipping point).
The Myers model also has three parameters
\begin{equation}
f(S(t) | r, K, \theta) = \frac{r S(t)^{\theta}}{1 + \frac{S(t)^\theta}{K}} \label{myers}
\end{equation}
where $\theta = 1$ corresponds to Beverton-Holt dynamics and $\theta$ > 2 leads to Allee effects and multiple stable states.
The Ricker model does not lead to multiple steady states.
The Allen model resembles the Ricker
dynamics with an added Allee effect parameter [@Courchamp2008], below
which the population cannot persist. The Myers model also has three
parameters and contains an Allee threshold, but unlike the Ricker model
saturates at high population size. The multiplicative log-normal stochasticity in Eqn 1
introduces one additional parameter $\sigma$ that must be estimated.
Because of our interest in management performance in the presence of
tipping points, all of our simulations are based on the Allen model.
The Allen model is thus the state of nature and is expected to provide
the best-case scenario. The Ricker model is a reasonable approximation of these dynamics far from the Allee threshold (but lacks threshold dynamics), while the Myers model shares
the essential feature of a threshold but differs in structure from the Allen model.
Throughout, we refer to the "True" model when the underlying parameters
_are known without error_, and refer to the "Allen" model
when these parameters have been estimated from the sample data.
We consider a period of `r Tobs` in which data are obtained to estimate the parameters or the GP. This is long
enough that the estimates do not depend on the particular realization,
and longer times are not likely to provide substantial improvement.
We fit each of the models is fit to the same data (Figure 1).
We inferred posterior distributions for the parameters of each model by Gibbs sampling (@Gelman2003 implemented in R
[@RTeam] using `jags`, [@R2jags]). We choose uniform priors
for all parameters (See appendix Tables S1-S3; R code provided). We
show one-step-ahead predictions of these model fits in Figure 1.
We tested each chain for Gelman-Rubin convergence and results were
robust to longer runs. For each simulation we also applied several
commonly used model selection criteria (AIC, BIC, DIC, see @Burnham2002)
to identify the best fitting model.
### The Gaussian Process model
The core difference for our purpose between the GP and the parametric
models is that the GP model is defined explicitly by the observed data. That is, it cannot be specified by the value of parameters alone.
The use of GPs to characterize dynamical systems is relatively new
[@Kocijan2005], and was first introduced in the context ecological
modeling and fisheries management in @Munch2005. GP models have
subsequently been used to test for the presence of Allee effects
[@Sugeno2013a], estimate the maximum reproductive rate [@Sugeno2013b],
determine temporal variation in food availability [@Sigourney2012], and
provide a basis for identifying model-misspecification [@Thorson2014].
An accessible and thorough introduction to the formulation and use of
GPs can be found in @Rasmussen2006.
A GP is a
stochastic process for which any realization of for
a finite sample of points $n$ is a multivariate normal distribution
of dimension $n$. Thus, to characterize the GP we need a mean function
and a covariance function. We proceed as follows.
As before, we assume that the data $X(t)$ are observed with process stochasticity around a mean function $ g(X(t))$
\begin{equation}
X(t+1) = g(X(t)) + \varepsilon,
\end{equation}
where $\varepsilon$ are IID normal random variables with zero-mean and variance $\sigma_n^2$. Note that we have chosen to
assume additive stochasticity. While we could consider log-normal stochasticity
as in the parametric models, we make this choice to emphasize that the
Gaussian process approach need not have structurally correct stochasticity to be effective.
Now consider the situation where we have a set of previous observations of the _timeseries_, $X_{\textrm{obs}}(t)$,
and wish to _predict_ future values $X_{\textrm{p}}(t)$. We consider the observations as pairs of points
$X_{\textrm{obs}}(t), X_{\textrm{obs}}(t+1)$ for $t=1, 2, …$ up through final observation time $T$. As described above, in the GP, a finite sample of these predictions at $n$ points into the future, $X_{\textrm p}(\tau)$ for $\tau \in [T+1, T+2, ..., T+n]$, is a multivariate normal ($\mathcal{N}$):
\begin{equation}
g(X_p|X_{\textrm{obs}}(t), X_{\textrm{obs}}(t+1)) \sim \mathcal{N}(E,C)
\end{equation}
with mean and covariance matrix given by (see Eqn 2.22 of @Rasmussen2006 for a more detailed derivation)
\begin{equation}
E = K(X_{\textrm p}(\tau), X_{\textrm{obs}}(t)) \left(K(X_{\textrm{obs}}(t),X_{\textrm{obs}}(t)) - \sigma \mathbb{I} \right)^{-1} X_{\textrm{obs}}(t+1)
\end{equation}
\begin{equation}
C = K(X_{\textrm p}(\tau), X_{\textrm p}(\tau)) - K(X_{\textrm p}(\tau), X_{\textrm{obs}}(t)) \left(K(X_{\textrm{obs}}(t),X_{\textrm{obs}}(t)) - \sigma \mathbb{I} \right)^{-1} K(X_{\textrm{obs}}(t), X_{\textrm p}(\tau))
\end{equation}
Where $K$, called the covariance kernel, measures how much influence one observation has on another. In this paper, for $K$ use a matrix whose $(i,j)$th element is given by:
\begin{equation}
K_{i,j}(x_i, y_j) = \exp\left( \frac{ -(x_i - y_j)^2}{2 \ell^2} \right)
\end{equation}
so that $\ell$ gives the characteristic length-scale over which correlation between
two observations decays.
We use inverse Gamma priors on both the length-scale $\ell$ and $\sigma$, thus for example
\begin{equation}
f(\ell; \alpha, \beta) = \frac{\beta^\alpha}{\Gamma(\alpha)} \ell^{-\alpha - 1}\exp\left(-\frac{\beta}{\ell}\right)
\end{equation}
For the prior on $\ell$, $\alpha = $ `r d.p[1]` and $\beta = $ `r d.p[2]`. The prior on $\sigma$, $\alpha = $ `r s2.p[1]` and $\beta = $ `r s2.p[2]`.
We use a Metropolis-Hastings Markov Chain Monte Carlo (@Gelman2003) to
infer posterior distributions of the parameters of the GP (Figure
S4, code in appendix). Since the posterior distributions differ
substantially from the priors (Figure S4), most of the information in the posterior comes from the data rather than the
prior belief.
The method of Gaussian Process Dynamic Programming (GPDP)
--------------------------------------------------
We derive the harvest policy from the
estimated GP by inserting it into a SDP algorithm. Given the GP posteriors, we construct the transition
matrix representing the probability of going to each state $X(t+1)$
given any current state $x(t)$ and any harvest $h(t)$ (See the function
`gp_transition_matrix()` in the provided R package). Given this transition
matrix, we use the same value iteration algorithm as in the parametric
case to determine the optimal policy.
In doing so,
the uncertainty in the future state under the GP, $X(t+1)$,
includes both process uncertainty (based on the estimation of $\sigma$)
and structural uncertainty of the posterior collection of curves.
Unlike the case of a parametric model, which is a single curve given the parameters,
the GP remains a distribution of curves.
Results
=======
Parametric and GP models for population dynamics
------------------------------------------------
To ensure our results are robust to the choice of parameters, we
consider 96 different scenarios, described in detail below. To help better
understand the process, we first describe in detail the
results of a single scenario.
```{r figure_1, dependson=c("assemble-models", "par-fns", "plot-options"), fig.cap="Points show the training data of stock-size over time. Curves show the posterior step-ahead predictions based on each of the estimated models. Observe that all models are fitting the data reasonably well. In the GP panel, black lines are means of the GP for random samples from the posterior distributions, while grey bands around them reflect the standard deviations of the GP distribution. (Recall that while given a particular set of parameters, a parametric model gives just one curve, the GP is always a distribution, whose mean we plot in black.)", fig.width=8, fig.height=6, cache=FALSE}
step_ahead_posteriors <- function(x){
gp_f_at_obs <- gp_predict(gp, x, burnin=1e4, thin=300)
df_post <- melt(lapply(sample(100, 30),
function(i){
data.frame(time = 1:length(x), stock = x,
GP = mvrnorm(1, gp_f_at_obs$Ef_posterior[,i], gp_f_at_obs$Cf_posterior[[i]]),
True = step_ahead(x,f,p),
MLE = step_ahead(x,f,est$p),
Allen = step_ahead(x, allen_f, pardist[i,]),
Ricker = step_ahead(x, ricker_f, ricker_pardist[i,]),
Myers = step_ahead(x, myers_f, myers_pardist[i,]))
}), id=c("time", "stock"))
}
df_post <- step_ahead_posteriors(x)
figure1b_posteriors <- ggplot(df_post) + geom_point(aes(time, stock)) +
geom_line(aes(time, value, col=variable, group=interaction(L1,variable)), alpha=0.2) +
facet_wrap(~variable) +
scale_colour_manual(values=colorkey) +
theme(legend.position="none")
figure1b_posteriors
write.csv(df_post, "components/data/figure1.csv")
```
All of the models fit the observed data rather closely and with relatively small uncertainty. In Figure 1, we show the posterior predictive curves. The training data of stock sizes observed over time are points, overlaid with the step-ahead predictions of each estimated model using the parameters sampled from their posterior distributions. Compared to the true model most estimates appear to over-fit, predicting patterns that are actually due purely to stochasticity. Model selection criteria (Table 1) penalize more complex models and show a preference for the simpler Ricker model over the models with alternative stable states (Allen and Myers). Details on MCMC estimates for each model, traces, and posterior distributions can be found in the appendix.
```{r deviances, dependson=c("dic_calc"), include=FALSE, echo=FALSE}
allen_deviance <- - posterior.mode(pardist[,'deviance'])
ricker_deviance <- - posterior.mode(ricker_pardist[,'deviance'])
myers_deviance <- - posterior.mode(myers_pardist[,'deviance'])
true_deviance <- 2*estf(c(p, sigma_g))
mle_deviance <- 2*estf(c(est$p, est$sigma_g))
aictable <- data.frame(Allen = allen_deviance + 2*(1+length(bayes_pars)), # +1 for noise parameter
Ricker = ricker_deviance + 2*(1+length(ricker_bayes_pars)),
Myers = myers_deviance + 2*(1+length(myers_bayes_pars)),
row.names = c("AIC"))
bictable <- data.frame(Allen = allen_deviance + log(length(x))*(1+length(bayes_pars)),
Ricker = ricker_deviance + log(length(x))*(1+length(ricker_bayes_pars)),
Myers = myers_deviance + log(length(x))*(1+length(myers_bayes_pars)),
row.names = c("BIC"))
pandoc.table(rbind(dictable, aictable, bictable))
```
```{r Table1, dependson="deviances", include=TRUE, results="asis", echo=FALSE}
pandoc.table(rbind(dictable, aictable, bictable), caption="Model selection scores for several common criteria (DIC: Deviance Information Criterion, AIC: Akaike Information Criterion, BIC: Bayesian Information Criterion) all select the wrong model. As the true (Allen) model is not distinguishable from the simpler (Ricker) model in the region of the observed data, this error cannot be avoided regardless of the model choice criterion. This highlights the danger of model choice when the selected model will be used outside of the observed range of the data.")
```
```{r figure_2, dependson=c("assemble-models", "plot-options"), fig.cap="The inferred Gaussian process compared to the true process and maximum-likelihood estimated process. We show the expected value for the function $f$ under each model. Two standard deviations from the estimated Gaussian process covariance with (light grey) and without (darker grey) measurement error are also shown. The training data are also shown as black points. The GP is conditioned on (0,0), shown as a pseudo-data point, as explained in the text", fig.width=8, fig.height=6, cache=FALSE}
x_grid_short <- x_grid[1:40]
gp_short <- gp_predict(gp, x_grid_short, burnin=1e4, thin=300)
models_posteriors <-
melt(lapply(sample(100, 50),
function(i){
sample_gp <- mvrnorm(1,
gp_short$Ef_posterior[,i],
gp_short$Cf_posterior[[i]])
data.frame(stock = x_grid_short,
GP = sample_gp,
y = sample_gp,
ymin = sample_gp - 2 * sqrt(gp_short$E_Vf),
ymax = sample_gp + 2 * sqrt(gp_short$E_Vf),
True = sapply(x_grid_short,f,0, p),
MLE = sapply(x_grid_short,f,0, est$p),
Allen = sapply(x_grid_short, allen_f, 0, pardist[i,]),
Ricker = sapply(x_grid_short, ricker_f, 0, ricker_pardist[i,]),
Myers = sapply(x_grid_short, myers_f, 0, myers_pardist[i,]))
}),
id=c("stock", "y", "ymin", "ymax"))
ggplot(models_posteriors) +
geom_ribbon(aes(x=stock, y=y, ymin=ymin, ymax=ymax, group=L1),
fill = "gray80",
data=subset(models_posteriors, variable == "GP")) +
geom_line(aes(stock, value, col = variable,
group=interaction(L1,variable)), alpha=0.2) +
geom_point(data = obs, aes(x,y)) +
xlab(expression(X[t])) + ylab(expression(X[t+1])) +
facet_wrap(~variable) +
scale_colour_manual(values=colorkey) +
theme(legend.position="none")
write.csv(models_posteriors, "components/data/figure2.csv")
write.csv(obs, "components/data/training_data.csv")
```
We show the mean inferred population dynamics of each model relative to the
true model used to generate the data in Figure 2, predicting
the relationship between observed population size (x-axis) to the population
size after recruitment the following year.
In addition to the raw data, the GP is conditioned
on going through the point 0,0 without error. This conditioning on (0,0) is equivalent to making the assumption that the population is closed, so that once it hits 0 it stays at 0, despite the lack of any data in the observed sequence to justify this. This assumption illustrates how the GP can capture common-sense biology without having to assume more explicit details about the dynamics at low population numbers that have never been observed. If the population were not closed, one could repeat the entire analysis without this assumption. Unlike parametric models, the GP corresponds to a
distribution of curves - the gray band with mean shown in black. Uncertainty in the parameters of the GP (not shown) further widens the band of possible population sizes.
In Figure S1 (see supplement), we show the performance of the models outside the observed training data.
```{r figure_3, dependson = c("assemble-opt"), fig.cap="The steady-state optimal policy (infinite boundary) calculated under each model. Policies are shown in terms of target escapement, $S(t)$, as under models such as this a constant escapement policy is expected to be optimal [@Reed1979].", fig.width=8, fig.height=6, cache=FALSE}
policies <- melt(data.frame(stock=x_grid, sapply(OPT, function(x) x_grid[x])), id="stock")
names(policies) <- c("stock", "method", "value")
ggplot(policies, aes(stock, stock - value, color=method)) +
geom_line(lwd=1.2) +
facet_wrap(~method) +
xlab("stock size, x(t)") +
ylab("escapement, S(t)") +
scale_colour_manual(values=colorkey, guide=FALSE)
write.csv(policies, "components/data/figure3.csv")
```
Despite the similarities in model fits to the observed data, the
policies inferred under each model differ widely (Figure 3).
Policies are shown in terms of target escapement, $S(t) = x_t - h$. Under models
such as this a constant escapement policy is expected to be optimal
[@Reed1979], whereby population levels below a certain size $S$ are
unharvested, while above that size the harvest strategy aims to return the
population to $S$. Whenever a model predicts that the population will not persist below a certain threshold, the optimal solution is to harvest the entire population immediately, resulting in an escapement $S=0$, as seen in the true (correct form, exact parameters) model, the Allen model (correct form, estimated parameters) and the GP.
Only the structurally correct model (Allen model) and the GP produce
policies close to the true optimum policy.
```{r figure_4, dependson=c("sim"), fig.cap="In the management context, GPDP outperforms approaches based on parametric models. We show 100 replicate simulations of the stock dynamics (Eqn 1) under the policies derived from each of the estimated models, as well as the policy based on the exact underlying model.", fig.width=6, fig.height=6, cache=FALSE}
ggplot(sims_data) +
geom_line(aes(time, fishstock, group=interaction(reps,method), color=method), alpha=0.1) +
scale_colour_manual(values=colorkey, guide = FALSE) +
xlab("stock size, x(t)") +
facet_wrap(~method) +
guides(legend.position="none", color=FALSE)
write.csv(sims_data, "components/data/figure4.csv")
```
In Figure 4, we show the consequences of managing 100 replicate
realizations of the simulated fishery under policies derived from each model
. The structurally
correct model under-harvests, leaving the stock to vary around its
unfished optimum. The structurally incorrect Ricker model over-harvests
the population past the tipping point consistently, resulting in the
immediate crash of the stock and thus leads to minimal long-term catch.
The results across replicate stochastic simulations are most
easily compared by using the relative differences in net present
value realized by each of the model (Figure 5). Although not perfect, the GPDP
consistently realizes a value close to the optimal solution,
and avoids ever driving the system across the tipping point,
which results in the near-zero value cases in the parametric models.
```{r figure_5, dependson=c("profits"), fig.cap="Histograms of the realized net present value of the fishery over a range of simulated data and resulting parameter estimates. For each data set, the three models are estimated as described above.", fig.width=6, fig.height=6, cache=FALSE}
ggplot(actual_over_optimal, aes(value)) + geom_histogram(aes(fill=variable)) +
facet_wrap(~variable, ncol=2) +
guides(legend.position = "none") +
xlab("Net present value") +
scale_fill_manual(values=colorkey, guide=FALSE) # density plots fail when delta fn
write.csv(actual_over_optimal, "components/data/figure5.csv")
```
Sensitivity Analysis
--------------------
These results are not sensitive to the modeling details of the simulation.
The GPDP estimate remains very close to the optimal solution (obtained
by knowing the true model) across changes to the training simulation,
scale of stochasticity, parameters or structure of the underlying model. In the Supplement, we consider
both a Latin hypercube approach and a more focused investigation of the effects of the relative
distance to the Allee threshold and the variance of process stochasticity.
Changing the intensity of the stochasticity or the distance between stable and unstable steady states
does not impact the performance of the GP relative to the optimal solution
obtained from the true model and true parameters. The parametric models are
more sensitive to this difference. Large values of $\sigma$ relative to the distance between the stable and unstable
point increases the chance of a stochastic transition. More precisely,
if we let $L$ is the distance between the stable and unstable steady states, then the probability that fluctuations drive the population across the unstable steady state scales as
$\exp\left(\frac{L^2}{\sigma^2}\right)$
(see @Gardiner2009 or @Mangel2006 for the derivation).
Thus, the impact of using a model that underestimates the risk of harvesting
beyond the critical point is considerable, since this such a situation
occurs more often. Conversely, with large enough distance between the
optimal escapement and unstable steady state relative to $\sigma$, the chance
of a transition becomes vanishingly small and all models can be estimated
near-optimally. Models that underestimate the cost incurred by population
sizes fluctuating significantly below the optimal escapement level will
not perform poorly as long as those fluctuations are sufficiently small.
The GPDP is only weakly influenced by increasing stochasticity or increasing
Allee effects over much of the range (Figure 6). Larger $\sigma$ or higher Allee
levels make even the optimal solution without any model or parameter uncertainty
unable to harvest the population effectively (e.g. the stochasticity is large enough to violate the
self-sustaining criterion of @Reed1979).
```{r figure_6, fig.width=6, fig.height=3, dependson=c("sensitivity-trends"), fig.cap="The effect of increasing noise or decreasing Allee threshold levels on the net present value of the fishery when managed under the GPDP, relative to managing under the true model (with known parameters). Other than the focal parameter (stochasticity, Allee threshold), other parameters are held fixed as above to illustrate this effect.", cache=FALSE}
plot_sigmas <- ggplot(vary_sigma, aes(noise, value)) +
geom_point() + ylab("Net present value") +
xlab("Level of growth stochasticity") + theme_bw()
plot_allees <- ggplot(vary_allee, aes(allee, value)) +
geom_point()+ ylab("Net present value") +
xlab("Allee threshold stock size") + theme_bw()
multiplot(plot_sigmas, plot_allees, cols=2)
write.csv(vary_sigma, "components/data/figure6a.csv")
write.csv(vary_allee, "components/data/figure6b.csv")
```
Discussion
==========
Simple, mechanistically motivated models offer the
potential to increase our basic understanding of ecological processes
[@Cuddington2013; @Geritz2012]. But such models can be both inaccurate and
misleading when used in a decision making framework.
In this paper we tackled two aspects of uncertainty that are
common to many ecological decision-making problems and fundamentally
challenging to existing approaches that largely rely on parametric
models:
1. We do not know the correct models for ecological systems.
2. We have limited data from which to estimate the model.
We have illustrated how the use of non-parametric methods provides
more reliable solutions in the sequential decision-making problem.
### Traditional model-choice approaches can be positively misleading.
Our results illustrate that model-choice approaches can be absolutely
misleading – by providing support to models that cannot capture tipping point
dynamics because they have fewer parameters and the data are far from the tipping point. That is, when the data
come from around the stable steady state,
all the parametric models are approximately linear and approximately
identical. Thus, it is intuitive that all model selection methods choose the simplest model.
In a complex world, the result is suboptimal. But in a world that might
contain tipping points, the result is disastrous.
Many managers both in fisheries and beyond face a similar situation: they
have not observed the population dynamics at all possible densities. A
lack of comprehensive data at all population sizes, combined with the
inability to formulate accurate population models for low population sizes
in the absence of data, makes this situation the rule more than the exception.
Relying on parametric models and model choice processes that favor simplicity
ignores this basic reality. For a long time, Carl Walters [e.g. @Walters1978] has argued that if we began
by fishing any newly exploited population down to very low levels and
then let it recover, we would be much better at estimating population
dynamics and thus predicting the optimal harvest levels. While certainly true,
this presents a rather risky policy in the face of potential tipping points.
The GPDP offers an alternative way to acknowledge this uncertainty.
### GPDP population dynamics capture larger uncertainty in regions where the data are poor
Parametric models perform most poorly when we seek a management strategy
outside the range of the observed data. The GPDP,
in contrast, leads to a predictive model that expresses a great deal of uncertainty
about the probable dynamics _outside_ the range of the observed data, while retaining very
good predictive accuracy _inside_ the range. The management policy
based on by the GPDP balances uncertainty outside the range of the observed data against
the immediate value of the harvest, and acts to stabilize the population
dynamics in a region of state space in which the predictions are
reliably reflected by the data.
Such problems are
ubiquitous across ecological decision-making
and conservation where the greatest concerns involve decisions that lead to population sizes that have never been observed and for which we do not know the response
-- whether this is the collapse of a fishery, the spread
of an invasive, or the loss of habitat.
### The role of the prior
Outside of the observed range of the data, the GP reverts to
the prior, and consequently the choice of the prior can also play a
significant role in determining the optimal policy inferred by the
SDP. In the examples shown here we have selected a prior that is both
relatively uninformative (due to the broad priors placed on its parameters
$\ell$ and $\sigma$ and simple (the choice of our covariance function, Eqns 12 and 13 ).
In practice, these
should be chosen to confer particular biological properties.
In principle, this may allow a manager to improve the performance of the
GPDP by adding detail as is justified.
For instance, it would be possible to use a linear or a Ricker-shaped
mean in the prior without making the much stronger assumption that the
Ricker is the structurally correct model [@Sugeno2013a]. One fruitful avenue of future research is identifying criteria to ensure the prior and the value function
are chosen appropriately for the problem at hand.
Acknowledgments
===============
This work was partially supported by NOAA-IAM grant to SM and Alec McCall and administered
through the Center for Stock Assessment Research, a partnership between
the University of California Santa Cruz and the Fisheries Ecology
Division, Southwest Fisheries Science Center, Santa Cruz, CA and by NSF
grant EF-0924195 to MM and NSF grant DBI-1306697 to CB.