Error installing FAS package from github

Following the instructions from fishR, I get the next message:

Error in function (type, msg, asError = TRUE) : couldn’t connect to host

Of course, my institution filters outgoing internet connection through a proxy, so here is some steps to overcome this issue:

## ------------------------------------------
install.packages("httr")
if (!require('devtools')) install.packages('devtools')
library(devtools)
library(httr)

## settings your proxy using 'httr' package
set_config(use_prox(url="http://proxyname.or.number",port=8080,username="youruser",password="yourpass"))
(set_config(verbose())) # check the configuration

## installing from github
devtools::install_github('FSA','droglenc')

reset_config()  # only if you want to revert the proxy configuration
## ------------------------------------------

 

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!!!!!

Using C++ directives to make ADMB report file

One way to get output from ADMB is making a suitable encode REPORT_SECTION, allowing custom output of the model result. However, when your model have hundreds of output, declare each one can be a repetitive task. Recently, I was exploring ADMB repositories and I got a nice option, provided by Steve Martell, that make use of C++ directives.

Consider this simple linear regression model (data gist):

DATA_SECTION
 init_int nobs
 init_matrix data(1,nobs,1,2)
 vector x(1,nobs)
 vector y(1,nobs)

PARAMETER_SECTION
 init_number a(1)
 init_number b(2)
 objective_function_value f
 LOCAL_CALCS
   x = column(data,2);
   y = column(data,1);
 END_CALCS

PROCEDURE_SECTION
 f=0;
 for (int i=1;i<=nobs;i++)
    {
     f+=square(y(i)-(a+b*x(i)));
    }

REPORT_SECTION
 report << "rss: " << f << endl;
 report << "data: " << data << endl;

Note that REPORT_SECTION have two lines with repetitive code. Using both #undef and #define C++ directives, you can save a lot of encode, principally if your model is multiparametric like most fisheries stock assessment models. All you need to do is add the following lines into the GLOBALS_SECTION and modify the REPORT_SECTION according to C++ directives:

GLOBALS_SECTION
 #undef myReport
 #define myReport(object) report << #object "\n" << object << endl;

DATA_SECTION
 init_int nobs
 init_matrix data(1,nobs,1,2)

PARAMETER_SECTION
 ....
 ....
 ....

REPORT_SECTION
 myReport(x);
 myReport(f);

Enjoy!!!

Otolith shape contour analysis

Otra vez tratando de ajustar los contornos de los otolitos de jurel. Aunque en R-project quizás sea más fácil, miren el libro Morphometrics with R para varios ejemplos utilizando elipticas de fourier,  ya llevo un tiempo realizandolo en Matlab©.  Les comparto parte de código…

%% Read data
clear
clc

path = ('home\jquiroz\Documents\Proyectos\FIP_2007_27 \analisis\Data Fourier\@COQUIMBO');
cd(path)
files = dir(fullfile('home\jquiroz\Documents\Proyectos\FIP_2007_27\analisis\Data Fourier\','\@COQUIMBO\*.txt'));

%% Fourier
for i = 1: length(files)
 outline = textread(char(files(i).name));
 iNoOfHarmonicsAnalyse = 100;
 bNormaliseSizeState = 0; % (T=1 - F=0)
 bNormaliseOrientationState = 0; % (T=1 - F=0)
 rFSDs(:,:,i) = fEfourier(outline, iNoOfHarmonicsAnalyse, bNormaliseSizeState, bNormaliseOrientationState);
end

%% Recontruc
for i = 1: length(files)
 iNoOfHarmonicsReconstruct = 50;
 iNoOfPointsReconstruct = 100;
 outln(:,:,i) = rEfourier(rFSDs(:,:,i), iNoOfHarmonicsReconstruct,iNoOfPointsReconstruct);
end

figure
for i=1:8
 subplot(3,3,i)
 outline = textread(char(files(i).name));
 plot(outline(:,1),outline(:,2),'-r'); hold on;
 plot(outln(:,1,i),outln(:,2,i),'--k')
 set(gca,'DataAspectRatio',[1 1 1],'PlotBoxAspectRatio',[1 1 1],...
 'XTickLabel',{'0.1';'0.5';'1'},'YTickLabel',{'0.1';'0.5';'1';'1.5'})
end
%% Plot 1

clear X XX Y YY
X = outln(:,1,1);
Y = outln(:,2,1);
for i = 1: length(files)-1
 X = [X outln(:,1,i+1)]; Y = [Y outln(:,2,i+1)];
end

XX = mean(X');
YY = mean(Y');
plot([XX XX(1)],[YY YY(1)],'-k','LineWidth',2);

%% discriminant analysis

nFDs = 50;
for i = 1: length(files)
 tmp = rFSDs(:,1:nFDs,i)';
 analy(:,i) = tmp(:);
end

%% PLot 2

n = [0.05 0.10 0.50 1]*iNoOfHarmonicsAnalyse;
color = {'r','k','m','g'};
p= 1; % otolito
outline_s = textread(char(files(p).name));

for i=1:length(n)
 iNoOfHarmonicsReconstruct = n(i);
 iNoOfPointsReconstruct = 100;
 outlnn = rEfourier(rFSDs(:,:,p), iNoOfHarmonicsReconstruct,iNoOfPointsReconstruct);
 hold on; plot(outline_s(:,1),outline_s(:,2),':k',outlnn(:,1),outlnn(:,2),'--rs',...
 'LineWidth',1,'MarkerEdgeColor','k','MarkerFaceColor',...
 char(color(i)),'MarkerSize',3);
 hold off
end