FLR packages
library(FLCore)
library(FLBRP)
library(FLasher)
library(FLife)
library(ggplotFL)
library(FLCandy)
library(plyr)
library(dplyr)
library(reshape)
library(ggplot2)
theme_set(theme_bw())
library(ggpubr)
library(ggcorrplot)
library(GGally)
library(fishblicc)
par = lhPar(FLPar(linf = 100))
eq = lhEql(par)
om = as(eq, "FLStock")
om = window(fwd(om, fbar = m(om)[4, -1], sr = eq), start = 19)
[1] "maxfbar has been changed to accomodate new plusgroup"
dimnames(om)$year = an(dimnames(om)$year) + 1980
range(om)["plusgroup"] = range(om)["maxage"]
om = fwd(om, fbar = window(fbar(om), start = 2021) %*% FLQuant(seq(1, 10, length.out = 61)),
sr = eq)
dim(window(om, end = 2080))
[1] 41 82 1 1 1 1
ak = invALK(par, cv = 0.1, age = an(dimnames(om)$age), bin = 1)
lfd = lenSamples(catch.n(om), ak, n = 5000)
units(lfd) = "cm"
# plotLengths(group(lfd, sum, year=year-year%%10))
dat = subset(transform(as.data.frame(lfd), lustrum = 5 * (year%/%5)), data > 0)
dt2 = ddply(dat, .(lustrum), with, {
lmodeFn <- function(len, n, bins = 25) {
dat = data.frame(bin = cut(len, breaks = bins), n = n)
res = ddply(dat, .(bin), with, data.frame(freq = sum(n)))
res = as.character(subset(res, freq == max(freq))[1, "bin"])
unbin(res)$mid
}
unbin = function(x) {
left = as.numeric(substr(x, 2, unlist(gregexpr(",", x)) - 1))
right = as.numeric(substr(x, unlist(gregexpr(",", x)) + 1, nchar(x) - 1))
mid = (left + right)/2
data.frame(left = left, right = right, mid = mid, n = x)
}
lc = lmodeFn(len, data) * 0.5
data.frame(lc = lc, lmean = weighted.mean(len, data * (len >= lc)))
})
gghistogram(dat, x = "len", weight = "data", bins = 120) + coord_flip() + scale_x_reverse() +
facet_grid(. ~ lustrum, scale = "free") + xlab("Length (cm)") + geom_vline(aes(xintercept = lc),
data = dt2, col = "red") + geom_vline(aes(xintercept = lmean), data = dt2, col = "blue") +
# scale_x_continuous(limits=c(0,60))+
theme_bw(16) + theme(legend.position = "none", axis.title.x = element_blank(), axis.text.x = element_blank(),
axis.ticks.x = element_blank())
Seasonal LFD
nseason = 4
nunit = 2
## Set up OM
par = lhPar(FLPar(linf = 40))
## Adjust parameters
dmns = dimnames(par)
dmns[["season"]] = seq(nseason)
dmns[["unit"]] = seq(nunit)
dmns = dmns[c(1, 4, 3, 2)]
par42 = FLPar(array(c(par), dim = laply(dmns, length), dimnames = dmns))
## Units
par42[, 2] = lhPar(FLPar(linf = 30))
## Seasons
par42["t0"] = par42["t0"] - seq(0, dim(par42)[3] - 1)/dim(par42)[3]
eq = lhEql(par)
om = as(eq, "FLStock")
om = window(fwd(om, fbar = m(om)[4, -1], sr = eq), start = 19)
[1] "maxfbar has been changed to accomodate new plusgroup"
dimnames(om)$year = an(dimnames(om)$year) + 1980
range(om)["plusgroup"] = range(om)["maxage"]
om = fwd(om, fbar = window(fbar(om), start = 2021) %*% FLQuant(seq(1, 10, length.out = 61)),
sr = eq)
dim(window(om, end = 2080))
[1] 41 82 1 1 1 1
## Make it seasonal with units
source("~/Desktop/active/FLCandy/R/seasonalise.R")
om42 = seasonalise(om)
om42 = FLCore:::expand(om, unit = seq(nunit), season = seq(nseason))
## ALKs
res = list()
ak42 = list()
for (i in seq(dim(par42)[2])) {
for (j in seq(dim(par42)[3])) res[[j]] = invALK(par42[, i, j], cv = 0.1, age = an(dimnames(om)$age),
bin = 0.1)
ak42[[i]] = FLPars(res)
}
## LFDs
res = list()
lfds = list()
for (i in seq(dim(par42)[2])) {
for (j in seq(dim(par42)[3])) res[[j]] = lenSamples(catch.n(om42)[, , i, j],
ak42[[i]][[j]], n = 500)
lfds[[i]] = FLQuants(res)
}
lfd42 = FLCore:::expand(lfds[[1]][[2]], season = seq(nseason), unit = seq(nunit))
for (iSeason in seq(dim(par42)[3])) for (iUnit in seq(dim(par42)[2])) lfd42[dimnames(lfds[[iUnit]][[iSeason]])[[1]],
, iUnit, iSeason] = lfds[[iUnit]][[iSeason]]
sv(lfd42)
ld()
dat = subset(as.data.frame(lfd42, drop = T), year %in% c(2000, 2080))
gghistogram(dat, x = "len", weight = "data", fill = "unit", position = "stack", bins = 60) +
coord_flip() + scale_x_reverse() + facet_grid(factor(season, labels = c("Season 1",
"2", "3", "4")) ~ factor(year, labels = c("Bmsy", "Over Fished")), scale = "free") +
xlab("Length") + #geom_vline(aes(xintercept=lc), data=dt2,col='red')+ geom_vline(aes(xintercept=lmean), data=dt2,col='blue')+ xlab("Length")
xlab("Length") + #geom_vline(aes(xintercept=lc), data=dt2,col='red')+ geom_vline(aes(xintercept=lmean), data=dt2,col='blue')+ +
xlab("Length") + #geom_vline(aes(xintercept=lc), data=dt2,col='red')+ geom_vline(aes(xintercept=lmean), data=dt2,col='blue')+ #geom_vline(aes(xintercept=lc),
xlab("Length") + #geom_vline(aes(xintercept=lc), data=dt2,col='red')+ geom_vline(aes(xintercept=lmean), data=dt2,col='blue')+ data=dt2,col='red')+
xlab("Length") + #geom_vline(aes(xintercept=lc), data=dt2,col='red')+ geom_vline(aes(xintercept=lmean), data=dt2,col='blue')+ geom_vline(aes(xintercept=lmean),
xlab("Length") + #geom_vline(aes(xintercept=lc), data=dt2,col='red')+ geom_vline(aes(xintercept=lmean), data=dt2,col='blue')+ data=dt2,col='blue')+
scale_x_continuous(limits = c(0, 60)) + scale_fill_manual("Component", values = c("red",
"blue")) + theme_bw(16) + theme(legend.position = "bottom", axis.title.y = element_blank(),
axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.text.x = element_blank(),
axis.ticks.x = element_blank())
library(remotes)
install_github("PaulAHMedley/fishblicc")
library(fishblicc)
## Prepare some data in the required format
frq=subset(lfd,year==2000&iter==1&len%in%10:101)$data
frq=c(34,42,45,68,65,56,91,123,126,126,113,123,107,111,113,103,91,85,74,53,43,
48,35,38,26,31,21,25,20,17,22,16,19,18,19,12,11,11,20,11,7,11,7,13,15,6,9,
7,8,4,7,6,8,3,3,5,6,7,5,3,7,7,2,5,1,3,2,2,0,0,1,2,0,1,3,2,3,2,2,2,0,1,2,2,
1,0,1,2,0,0,0)
dl=blicc_dat(
#model_name ="Base: Mk ~ Length-inverse",
LLB = 10:100, # Lower boundaries of length bins
fq = list(`Test`=frq), # frequency data
Linf = c(100,10),
sel_fun = c("logistic"), # Selectivity functions
Catch = c(1) # Relative catch for each gear
#gear_names = c("test"),
#Mk = 1.5,# Natural mortality prior mean (at reference length)
#ref_length =50, # Reference length for the length-inverse natural mortality model
#a =0.3, # Length-weight scale parameter (optional)
#b =3, # Length-weight exponent
#L50 =52 # Length at 50% maturity
)
## Fit the model to these data
slim <- blicc_mpd(dl)
rp_res <- blicc_ref_pts(slim, dl)
plot_expected_frequency(rp_res, gear=1:2)
stf=blicc_fit(dl, ntarget=100, nwarmup=200, nchain=1)
Back to Top