Kobe plot updated

Some helpful comment on my last post (in spanish) helped me to improve the kobe-plot, I hope this will help to somebody. Thanks go to Dennis Murphy at ggplot-groups for your extended response to my question about the setting-layer using ggplot.

Download the toy data example from github, and run the next script:

rm(list=ls())
library(ggplot2)
library(reshape2)
library(KernSmooth)

# ---------------- DATA -----------------------------------------

data <- dget('kobe_data.txt')
dta <- data[[1]] # fishing mortality and spawning biomas
v <- data[[2]] # Reference points {F20%, F40% and SBo}
F20 <- v[1]
F40 <- v[2]
F40.F20 <- F40/F20
So <- v[3]
yymax <- 1.4*max(dta$Ft,F40.F20)

x.mcmc <- dta$SB[33]/v[3]*rlnorm(350,0,0.2)
y.mcmc <- dta$Ft[33]/v[1]*rlnorm(350,0,0.12)
bw = 15
levs = seq(0,1,0.1)
bwx = (max(x.mcmc) - min(x.mcmc))/bw
bwy = (max(y.mcmc) - min(y.mcmc))/bw
est <- bkde2D(cbind(x.mcmc, y.mcmc), bandwidth = c(bwx, bwy), gridsize = c(51,51))

est$fhat = est$fhat/max(est$fhat)
tmp.1 <- expand.grid(x = est$x1, y = est$x2)
tmp.2 <- melt(est$fhat)
tmp.3 <- data.frame(x=tmp.1$x, y=tmp.1$y, z=tmp.2$value)

# ---------------- KOBE - PLOT -------------------------------

kobe <- ggplot(dta) +
theme_bw() + scale_x_continuous(limits=c(0,1)) +
scale_y_continuous(limits=c(0,yymax),
breaks = round(sort(c(0.5,1.0,F40.F20,1+(F40.F20-1)/2,yymax-(yymax-F40.F20)/2)),digits = 3),
labels = c(' 0.5','1.0', paste(round(1+(F40.F20-1)/2,digits=2)),
expression(F[20]/F[40]),paste(round(yymax-(yymax-F40.F20)/2,digits=2))))

kobe <- kobe +
geom_rect(xmin = 0.0, xmax = 1.0, ymin = 0.0, ymax = yymax, fill = 'yellow', alpha = 0.65) +
geom_rect(xmin = 0.0, xmax = 0.4, ymin = 1.0, ymax = yymax, fill = '#FF0000FF', alpha = 0.35) +
geom_rect(xmin = 0.4, xmax = 1.0, ymin = 1.0, ymax = 0.0, fill = 'green', alpha = 0.55) +
geom_rect(xmin = 0.0, xmax = 0.2, ymin = F40.F20, ymax = yymax, fill = 'red', alpha = 0.85) +
geom_contour(data = tmp.3, aes(x = x, y = y, z = z, colour = ..level..), bins = 9, size = 0.5) +
scale_colour_gradient(low = 'yellow', high = 'blue') +
geom_path(aes(x = SB/v[3], y = Ft/v[1]), linetype = 2, size = 0.4) +
geom_point(aes(x = SB/v[3], y = Ft/v[1], size = yr, fill = yr), shape = 21) +
scale_size(range=c(0,4)) + labs(y = 'Mortalidad por pesca, relativa a F40%',
x = 'Fracción de la Biomasa Virginal',colour = 'Prob', size = 'Año', fill = 'Año') +
scale_fill_gradient(limits=c(dta$yr[1], dta$yr[length(dta$yr)]), low="black", high='#0092FFFF') + # colour = yr
geom_text(data = head(dta,1), aes(x = SB/v[3], y = Ft/v[1]), label = dta$yr[1], colour = 'black', vjust = 2, size = 3) +
geom_text(data = tail(dta,1), aes(x = SB/v[3], y = Ft/v[1]), label = dta$yr[length(dta$yr)], colour = 'black', vjust = -2, size = 3) +
geom_text(data = dta, aes(0.1,yymax-(yymax-F40.F20)/2), label = 'Colapso', size = 4, colour = 'white', alpha = 0.2) +
geom_text(data = dta, aes(0.2,0.5), label = 'Sobrepescado \n BD < 0.4*BDo', size = 3, colour = 'grey50', alpha = 0.2) +
geom_text(data = dta, aes(0.7,F40.F20), label = 'Sobrepesca \n Ft > F40%', size = 3, colour = 'grey50', alpha = 0.2) +
geom_text(data = dta, aes(0.7,0.5), label = 'Sin Riesgo', size = 4, colour = 'black', alpha = 0.2)

print(kobe)
ggsave(kobe, file='kobe_plot_esc1.png')

Kobe plot using ggplot

Es recurrente utilizar gráficas de fases, tambien conocidas como gráficas Kobe,  para evaluar la condición de una población explotada con base en la mortalidad por pesca (F) y biomasa (en este caso, biomasa desovante SB) asociadas con puntos biológicos de referencias (ver Figura abajo). Sin embargo, la construcción de estas gráficas no es una tarea trivial, más aún cuando es necesario evaluar un número importante de escenarios.

ggplot es una librería de R-project altamente versatil que me ha facilitado la construcción de las gráficas Kobe, si les interesa pueden bajar el script ( como una incipiente función) y la trayectoria para los datos (uno de los tantos escenarios) de merluza del sur. Ejecuten en la consola de R las siguientes instrucciones:

rm(list=ls())
library(ggplot2)
library(reshape)
library(reshape2)

data <- dget('kobe_data.txt')
dta <- data[[1]] # fishing mortality and spawning biomas
v <- data[[2]] # Reference points {F20%, F40% and SBo}
kobePlot(dta, v) 

Si la mortalidad por pesca actual (Ft) está por encima del punto biológico de referencia (F40% en el caso de la Figura), se considera que está ocurriendo sobrepesca; si la biomasa actual (B; o alguna medida de producción del desove) está por debajo del 40% del nivel que existió previo a la explotación pesquera, se considera que la población está sobrepescada.

La idea de este gráfico es representar la trayectoria de la población a lo largo del tiempo para permitir ver el estatus histórico de la población. Típicamente, la población comienza en el cuadro superior derecho, a medida que se desarrolla la pesquería se traslada al cuadrado inferior izquierdo, donde la población ingresa a un estado sobreexplotación. Teóricamente, al aplicar una ordenación apropiada, la trayectoria deberia girar alrededor del centro de la gráfica.

Los comentarios y mejoras son más que bienvenidos!!!!!