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)
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())
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)
stf=blicc_fit(dl, ntarget=100, nwarmup=200, nchain=1)
Back to Top