Kajian Stok (Stock Assessment)
Ahmad Romdon
3/3/2024
Memuat Paket
Memuat paket tambahan yang digunakan.
library("TropFishR")Memuat dan Manipulasi Data
Memuat dan manipulasi data yang akan digunakan.
# memuat data
lfq3 <- read.csv("S. fasciatus.csv")
head(lfq3)## midlength X10.1.2019 X9.2.2019 X6.3.2019 X7.4.2019 X8.5.2019 X7.6.2019
## 1 11.2 0 0 13 2 1 0
## 2 12.7 0 3 32 7 15 0
## 3 14.2 0 4 39 6 10 0
## 4 15.7 3 7 33 11 10 4
## 5 17.2 6 9 31 16 7 4
## 6 18.7 6 4 25 8 13 8
## X4.7.2019 X5.8.2019 X10.9.2019 X3.10.2019 X6.11.2019 X8.12.2019
## 1 0 0 0 0 0 4
## 2 0 10 6 0 17 17
## 3 3 19 13 10 33 28
## 4 3 5 14 16 16 22
## 5 10 6 9 8 17 22
## 6 33 9 4 6 5 12lfq4 <- read.csv("S. fasciatus_catch.csv", header = TRUE)
head(lfq4)## Length Catch
## 1 11.2 20
## 2 12.7 107
## 3 14.2 165
## 4 15.7 144
## 5 17.2 145
## 6 18.7 133lwr <- read.csv("S. fasciatus_LWR.csv", header = TRUE)
head(lwr)## panjang_mm bobot_g logL logW
## 1 64 5.29 1.806 0.723
## 2 78 8.85 1.892 0.947
## 3 80 10.60 1.903 1.025
## 4 84 12.12 1.924 1.084
## 5 86 12.05 1.934 1.081
## 6 86 12.16 1.934 1.085# manipulasi data
dates <- colnames(lfq3) [-1]
dates <- strsplit(dates, "X")
dates <- unlist(lapply(dates, function(x) x[2]))
dates <- as.Date(dates, "%d.%m.%Y")
# data untuk analisis
mina <- list(dates = dates,
midLengths = lfq3$midlength,
catch = as.matrix(lfq3[,-1]))
class(mina) <- "lfq"
mina$catch[is.na(mina$catch)] <- 0Restrukturisasi Data
# plot frekuensi tangkapan (raw)
plot(mina, Fname = "catch")
# plot frekuensi (restructured)
mina <- lfqRestructure(mina, MA=3)
plot(mina, Fname = "rcounts")
Parameter Pertumbuhan
Pendugaan parameter pertumbuhan meliputi nilai panjang asimptotik ikan (Linf), koefisien pertumbuhan (K), dan umur teoritik ikan pada saat panjang ikan sama dengan nol (t0).
Linf dan K
Banyak cara dalam menentukan Linf dan K, penentuan metode dapat didasarkan pada literatur yang ada.
Contoh dalam perhitungan ini menggunakan response surface analysis dengan metode optimise (metode yang digunakan pada FiSAT).
# response surface analysis (varies Linf and K)
fit1 <- ELEFAN(
lfq = mina, MA = 3,
Linf_range = seq(28.1, 32, length.out = 100), #linF_range bisa didapatkan dari referensi
K_range = exp(seq(log(0.1),log(1.0), length.out = 100)), #length.out menentukan kerapatan
method = "optimise", #metode yang dapat digunakan optimise dan cross
cross.date = mina$dates[3],
cross.midLength = mina$midLengths[4],
contour = TRUE, add.values = FALSE,
hide.progressbar = FALSE #ubah menjadi 'TRUE' untuk mengikuti perkembangan algoritma
)
# memberi titik pada plot berdasarkan nilai Linf dan K
points(fit1$par["Linf"], fit1$par["K"], pch="*", cex=2, col=2)
# menampilkan hasil perhitungan parameter
unlist(fit1$par)
fit1$Rn_max
fit1$agemax
# plot frekuensi (restructured)
plot(fit1)
Selain itu, metode lain yang dapat digunakan adalah:
# Powell_Wetherall
PW <- powell_wetherall(mina, catch_columns = 1:3, reg_int = NULL)
PW$Linf_est
PW$confidenceInt_Linf
# K-Scan (varies K, Linf is fixed)
fit2 <- ELEFAN(
lfq = mina, method = "optimise", #metode yang dapat digunakan optimise dan cross
Linf_fix = PW$Linf_est,
K_range = round(exp(seq(from = log(0.1), to = log(1.0), length.out = 100)),2), #length.out menentukan kerapatan
cross.date = mina$dates[3], cross.midLength = mina$midLengths[4]
)
# bin with maximum score crossed
fit3 <- ELEFAN(
lfq = mina, method = "optimise", #metode yang dapat digunakan optimise dan cross
Linf_range = seq(from = 28.1, to = 34, length.out = 100), #length.out menentukan kerapatan
K_range = exp(seq(from = log(0.1), to = log(1.0), length.out = 100)),
cross.max = TRUE,
contour = TRUE
)
# default
fit4 <- ELEFAN(
lfq = mina, method = "optimise", #metode yang dapat digunakan optimise dan cross
Linf_range = seq(from = 28.1, to = 34, length.out =100), #length.out menentukan kerapatan
K_range = exp(seq(from = log(0.1), to = log(1.0), length.out = 100)),
contour = TRUE
)
# ELEFAN-SA (tanpa fluktuasi musiman)
set.seed(1234)
fit5 <- ELEFAN_SA(
lfq = mina,
seasonalised = FALSE, #TRUE, jika ada fluktuasi musiman
init_par = list(Linf=30.51, K=.6, t_anchor=0.41, ts=0.5, C=0.5),
low_par = list(Linf=28.32, K=.21, t_anchor=0.05, ts=0, C=0),
up_par = list(Linf=32.69, K=.99, t_anchor=0.78, ts=1, C=1),
SA_temp = 1e5,
SA_time = 60,
maxit = NULL,
MA = 3,
plot.score = TRUE,
verbose = TRUE
)
# ELEFAN-GA (tanpa fluktuasi musiman)
set.seed(1234)
fit6 <- ELEFAN_GA(
lfq = mina,
seasonalised = FALSE, #TRUE, jika ada fluktuasi musiman
low_par = list(Linf=28.32, K=.21, t_anchor=0.05, ts=0, C=0),
up_par = list(Linf=32.69, K=.99, t_anchor=0.78, ts=1, C=1),
popSize = 50,
pmutation = 0.1,
maxiter = 100,
run = 100,
MA = 3,
plot.score = TRUE,
monitor = FALSE,
parallel = FALSE
)t0
t0 <- -10^(-0.3922-(0.2752*log10(fit1$par$Linf))-(1.038*log10(fit1$par$K)))
print(t0)## [1] -0.4174318Pola Rekrutmen
mina$Linf <- fit1$par$Linf
mina$K <- fit1$par$K
mina$t0 <- t0
# masukan waktu pengambilan sampel
s_dates <- as.POSIXlt(mina$dates, format="%Y-%d-%m")
mina$dates## [1] "2019-01-10" "2019-02-09" "2019-03-06" "2019-04-07" "2019-05-08"
## [6] "2019-06-07" "2019-07-04" "2019-08-05" "2019-09-10" "2019-10-03"
## [11] "2019-11-06" "2019-12-08"rec <- recruitment(param = mina, tsample = s_dates$yday/365, plot = TRUE)
plot(rec$per_recruits, type = "l", ylab = "recruitment (%)", xlab = "bulan", xaxt = 'n', las = 1)
axis(1, 1:12, rec$months_abb)
Laju Mortalitas dan Eksploitasi
Mortalitas
Ms <- M_empirical(Linf = fit1$par$Linf, K_l = fit1$par$K, tmax=fit1$agemax, temp = 28.6, method = c("AlversonCarney","Hoenig","Pauly_Linf","Then_tmax","Then_growth"))
Ms## M
## Alverson and Carney (1975) 0.515
## Hoenig (1983) - Joint Equation 0.548
## Hoenig (1983) - Fish Equation 0.527
## Pauly (1980) - Length Equation 0.962
## Then (2015) - tmax 0.729
## Then (2015) - growth 0.665# memilih nilai M
bestM <- Ms[1,] #berdasrkan Alverson and Carney (1975)
# pemilihan nilai M dapat didasarkan literatur yang tersediaEksploitasi
mina2 <- synLFQ3
mina2$midLengths <- lfq4[, 1]
mina2$midLengths <- as.numeric(mina2$midLengths)
mina2$catch <- as.numeric(mina2$catch)
mina2$Linf <- fit1$par$Linf
mina2$K <- fit1$par$K
mina2$t0 <- t0
mina2$catch <- as.matrix(lfq4[, -1])
mina2$catch <- as.numeric(mina2$catch)
# estimasi nilai dengan kurva tangkapan (catch curve)
res_cc <- catchCurve(mina2, calc_ogive = TRUE, reg_int = NULL, auto = TRUE)
# masukan nilai pada data
mina2$Z <- res_cc$Z
mina2$M <- bestM
mina2$FM <- as.numeric(mina2$Z - bestM)
mina2$E <- mina2$FM/mina2$Z
# tampilkan hasil perhitungan
res <- res_cc
res$M <- bestM
res$FM <- mina2$FM
res$E <- mina2$E
res$t0 <- mina2$t0
res$midLengths <- NULL; res$catch <- NULL; res$t_midL <- NULL
res$lnC_dt <- NULL; res$reg_int <- NULL; res$linear_mod <- NULL
res$confidenceInt <- NULL; res$intercept <- NULL; res$linear_mod_sel <- NULL
res$Sobs <- NULL; res$ln_1_S_1 <- NULL; res$Sest <- NULL
unlist(res)## Linf K t0 Z se t50
## 31.33030303 0.39000000 -0.41743175 1.18350009 0.09386811 1.01192602
## t75 t95 L50 L75 L95 M
## 1.13932494 1.35337345 13.38844814 14.25811571 15.62541904 0.51500000
## FM E
## 0.66850009 0.56485005Hubungan Panjang-Berat
# model regresi
lm1 <- lm(logW ~ logL, data = lwr) #data yang digunakan dalam bentuk log
# plot
plot(logW ~ logL, data=lwr, pch=19, xlab="log(Panjag Total)",ylab="log(Bobot)")
# hasil regresi
summary(lm1)##
## Call:
## lm(formula = logW ~ logL, data = lwr)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.197633 -0.014832 0.000959 0.016196 0.299775
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.84447 0.02482 -195.2 <2e-16 ***
## logL 3.06472 0.01148 267.0 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03015 on 1406 degrees of freedom
## Multiple R-squared: 0.9807, Adjusted R-squared: 0.9806
## F-statistic: 7.128e+04 on 1 and 1406 DF, p-value: < 2.2e-16Model Yield Per Recruit (YpR)
# masukan nilai parameter
mina2$a <- 10^lm1$coefficients["(Intercept)"]
mina2$b <- lm1$coefficients["logL"]
mina2$Lmat <- 190 #panjang pertama kali matang gonad (mm)
mina2$Wmat <- mina2$a * mina2$Lmat ^ mina2$b #berat pertama kali matang gonad (g), diestimasi dari persamaan hubungan panjang-berat
mina2$tr <- mean(rec[["tS_frac"]]) #rata-rata umur rekrutmen
# parameter selektivitas
selectivity_list <- list(selecType = "trawl_ogive",
L50 = res_cc$L50, L75 = res_cc$L75)Model Thompson and Bell
# nilai F dengan rentang
TB1 <- predict_mod(mina2, type = "ThompBell",
FM_change = seq(from = 0, to = 2, by = 0.1), #rentang nilai dapat ditentukan berdasarkan hasil yang didapatkan atau literatur
stock_size_1 = 1,
curr.E = mina2$E,
curr.Lc = res_cc$L50,
s_list = selectivity_list,
plot = TRUE, hide.progressbar = TRUE)
# nilai F dan Lc dengan rentang
TB2 <- predict_mod(mina2, type = "ThompBell",
FM_change = seq(from = 0, to = 2, by = 0.1), #rentang nilai dapat ditentukan berdasarkan hasil yang didapatkan atau literatur
Lc_change = c(13,15,20),
stock_size_1 = 1,
curr.E = mina2$E,
curr.Lc = res_cc$L50,
s_list = selectivity_list,
hide.progressbar = TRUE, plot = TRUE)
# Plot
par(mfrow = c(2,1), mar = c(4,5,2,4.5), oma = c(1,0,0,0))
plot(TB1, mark = TRUE)
mtext("(a)", side = 3, at = -0.1, line = 0.6)
plot(TB2, type = "Isopleth", xaxis1 = "FM", mark = TRUE, contour = 6)
mtext("(b)", side = 3, at = -0.1, line = 0.6)
# hasil referensi biologi Biological reference levels
TB1$df_Es #biologi
TB2$df_Es
TB1$currents #biomasa
TB2$currentsModel Beverton and Holt’s
# nilai F dengan rentang
ypr1 <- predict_mod(mina2, type = "ypr",
FM_change = seq(from = 0, to = 1.5, by = 0.05),
Lc_change = seq(),
stock_size_1 = 1,
curr.E = mina2$E,
curr.Lc = res_cc$L50,
s_list = selectivity_list,
plot = TRUE, hide.progressbar = TRUE)
# nilai F dengan rentang dan Lc
ypr2 <- predict_mod(mina2, type = "ypr",
FM_change = seq(from = 0, to = 1.5, by = 0.05),
Lc_change = seq(from = 10, to = 18, length.out = 3),
stock_size_1 = 1,
curr.E = mina2$E,
curr.Lc = res_cc$L50,
s_list = selectivity_list,
plot = TRUE, hide.progressbar = TRUE)
# Plot
par(mfrow = c(2,1), mar = c(4,5,2,4.5), oma = c(1,0,0,0))
plot(ypr1, mark = TRUE)
mtext("(a)", side = 3, at = -0.1, line = 0.6)
plot(ypr2, type = "Isopleth", xaxis1 = "FM", mark = TRUE, contour = 6)
mtext("(b)", side = 3, at = -0.1, line = 0.6)
# hasil referensi biologi Biological reference levels
ypr1$df_Es #biologi## Lc tc F01 Fmax F05 E01 Emax E05
## 1 1 -0.3342562 0.6 0.6 0.2 0.5381166 0.5381166 0.2797203ypr2$df_Es## Lc tc F01 Fmax F05 E01 Emax E05
## 1 10 0.5683555 0.60 0.60 0.2 0.5381166 0.5381166 0.2797203
## 2 14 1.1008482 0.75 0.75 0.2 0.5928854 0.5928854 0.2797203
## 3 18 1.7737116 0.90 0.90 0.2 0.6360424 0.6360424 0.2797203ypr1$currents #biomasa## curr.Lc curr.tc curr.E curr.F curr.YPR curr.YPR.rel
## (Intercept) 13.38845 1.011926 0.5648501 0.6685001 0.04479589 0.05091599
## curr.BPR curr.BPR.rel
## (Intercept) 0.1355651 0.1540863ypr2$currents## curr.Lc curr.tc curr.E curr.F curr.YPR curr.YPR.rel
## (Intercept) 13.38845 1.011926 0.5648501 0.6685001 0.04479589 0.05091599
## curr.BPR curr.BPR.rel
## (Intercept) 0.1355651 0.1540863Session Info:
## R version 4.3.2 (2023-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: Asia/Jakarta
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] TropFishR_1.6.3 Matrix_1.6-5
##
## loaded via a namespace (and not attached):
## [1] expm_0.999-9 jsonlite_1.8.8 compiler_4.3.2 highr_0.10
## [5] crayon_1.5.2 Rcpp_1.0.12 GA_3.2.4 msm_1.7.1
## [9] stringr_1.5.1 parallel_4.3.2 jquerylib_0.1.4 splines_4.3.2
## [13] yaml_2.3.8 fastmap_1.1.1 lattice_0.21-9 R6_2.5.1
## [17] plyr_1.8.9 knitr_1.45 iterators_1.0.14 MASS_7.3-60
## [21] bslib_0.6.1 rlang_1.1.3 cachem_1.0.8 stringi_1.8.3
## [25] xfun_0.42 sass_0.4.8 doParallel_1.0.17 cli_3.6.2
## [29] magrittr_2.0.3 digest_0.6.34 foreach_1.5.2 grid_4.3.2
## [33] rstudioapi_0.15.0 mvtnorm_1.2-4 lifecycle_1.0.4 GenSA_1.1.14
## [37] evaluate_0.23 glue_1.7.0 codetools_0.2-19 survival_3.5-7
## [41] reshape2_1.4.4 rmarkdown_2.25 tools_4.3.2 htmltools_0.5.7