在R中使用ELEFAN做漁業資源分析

tags: 用R做 生態學

在漁業資源調查中，體長是最容易得到的測量值之一。透過統計體長分布隨時序的變化，我們可以估算魚群（或其他經濟性捕撈動物）的成長速度變化、幼體補充時間等等。奧地利生物學家 Ludwig von Bertalanffy 於 1930 年代發展出了魚類的成長方程式 von Bertalanffy growth function（VBGF）如下：

$$𝐿𝑡={𝐿}_{\mathrm{\infty }}[1-{𝑒}^{-𝐾(𝑡-{𝑡}_0)+(\frac{CK}{2𝜋})sin2𝜋(𝑡-{𝑡}_{𝑠})}]$$
其中

$L_t$ 是年齡為

$t$ 之體長；

$L_{\mathrm{\infty }}$ 為極限體長（asymptotic length），也是 VBGF 要逼近的目標；

$t_0$ 為

$L_t=0$ 的理論年齡，一般在推論時，會假設

$t_0=0$；

$K$ 為成長參數（growth coefficient）；

$e$ 是自然指數；

$C$ 為季節性振蕩幅度參數（seasonal oscillation），與年溫差成正比，溫差 1 度相當於變化 0.1 單位，其值介於 0（無振盪）到 1（最大振盪）之間；

$t_s$ 為成長曲線起伏變化的起點，其值視生物棲息地的冬季低迷點（winter point，WP）而定，把一年 12 個月換算為 0 到 1 之間，生物一年中成長最慢的時段即為其 WP 所在。在北半球一般為 2 月中旬，即 WP=0.2。

因為 VBGF 估算起來相當複雜，因此聯合國糧食及農業組織（Food and Agriculture Organization）與世界魚類研究中心（International Center for Living Aquatic Resources Management, ICLARM）聯合研發了一套電腦體長頻度分析系統 ELEFAN（Electronic Length Frequency Analysis）提供全球漁業資源估算者來使用，可以搭配一樣是 FAO 推出的軟體 FishStat 系列（現在在他們網頁上的是 FishStatJ）。類似的系統還有 MULTIFAN，搭配 MULTIFAN-CL 軟體使用。

當初為了估算計畫需要的貝類 VBGF，本來想用 MULTIFAN-CL 來做統計，沒想到註冊了之後遲遲沒收到回信，感覺可能沒人維護了，因此想找找在 R 裡有沒有替代的套件。結果還真的找到了一個 TropFishR，由 Taylor 與 Mildenberger 2017 年發表於 Fisheries Management and Ecology 期刊。

用法有教學，大致上分為幾步：

1. 輸入資料與決定組距
   這邊組距的決定參考葉樹藩在《試驗設計學》中引用 Snedecor 及 Sturges 提出的方法：
   
   $$𝐵𝑖𝑛=\frac{Range}{\frac{range}{SD}\times 4}$$
   
   $$Bin=\frac{Range}{1+3.322{\log }_N}$$
   其中
   $Range$ 為組距、
   $SD$ 為標準差、
   $N$ 為數量。
   可以將所得的 2 個組距先平均當作組距來使用，看看是否便於了解資料的分布情形。組距很重要，會影響到之後的分析結果。

D <- read.csv("資料.csv")
D$date <- as.Date(as.character(D[["Month"]]), format = "%Y%m%d")
lfq.o <- lfqCreate(D, Lname = "Width", Dname = "date")
plot(lfq.o, Fname = "catch")

range = range(D$Width)[2]-range(D$Width)[1]
bin = round(mean(c(range/((range/sd(D$Width))*4), range/(1+3.322*log(length(D$Width))))), digits=2)
bin

#先以估算出的bin來作圖看看
lfq <- lfqModify(lfq.o, bin_size = bin)
plot(lfq, Fname = "catch", main="Length distribution of __Species__")

#可以用之後估算的VBGF參數來畫看看
tmp <- lfqFitCurves(lfq, par = list(Linf=53.2944060, K=0.3532444, t_anchor=0.4425567, C=0.4548343, phiL=3.0014386, ts=0.6951600), draw = TRUE, col=2, lty=1)

2. 尋找極限體長
   使用兩種方法來尋求極限體長，需要調整極限體長範圍與 cross.date、cross.date 兩個參數，最終獲得最符合現狀的反應曲面圖（response surface），決定 VBGF 中的諸多參數。
   另外也可以透過 Powell Wetherall plot 來決定極限體長，所需要的參數比較不需要花費時間調整，就用之前估算組距繪出的 lfq（length-frequency）之中參數來計算即可，最後會輸出一個極限體長的範圍。其實也可以試看看參考這個方法估算出的極限體長，丟到 ELEFAN 裡面來估計看看。

#ELEFAN 找極限體長
length(lfq$midLengths)
lfq$midLengths

fit <- ELEFAN(lfq, method = "cross", MA = 5,
              Linf_range = seq(from = 40, to = 65, length.out = 20),
              K_range = seq(0.3,2,0.1), addl.sqrt = TRUE,
              cross.date = lfq$dates[12], cross.midLength = lfq$midLengths[9],
              contour = TRUE,
              hide.progressbar = F)

fit$Rn_max; unlist(fit$par)
points(fit$par["Linf"], fit$par["K"], pch="*", cex=2, col=2)

plot(fit, main="Growth curve of __Species__")
lfqFitCurves(fit, lty=1, col=2, par=fit$par, draw=TRUE)

#Powell Wetherall plot & Linf 找極限體長
PW <- powell_wetherall(param = lfq,
                       catch_columns = 1:ncol(lfq$catch),
                       reg_int = c(10, 60))
paste("Linf=", round(PW$Linf_est, digits = 2), "±", round(PW$se_Linf, digits = 2))

以 2020 至 2021 年的船形薄殼蛤（Laternula marilina (Reeve, 1860)，又稱公代）體長資料作為輸入值，這步繪出來的生長曲線如下圖：

3. 計算幼體入添比例
   使用 ELEFAN 估算出的 VBGF 參數來推斷每個月的幼體入添比例。這邊需要手動加入的參數是 lfq$t0，代表體長為 0 時的時間（將 1 月 1 日到 12 月 31 日化為 0 到 1 來表示），最後會輸出個月的幼體入添比例圖，另外也可以輸出表格，可以了解每個月估計有多少幼體會加入族群。

#計算入添比例
lfq$Linf <- fit2$par["Linf"]$Linf
lfq$K <- fit$par["K"]$K
lfq$t0 <- 0.3755
s_dates <- as.POSIXlt(lfq$dates, format="%d.%m.%Y")
rctmnt <- recruitment(lfq, tsample = s_dates$yday/365, plot = T)
#windowsFonts(GJ = windowsFont("Gen Jyuu Gothic P Medium"))
title(main = "Monthly recruitment of ___species___", cex.main = 1.6,
      #xlab = "月份", family = "GJ",
      cex.lab = 1.2)
rctmnt$per_recruits

一樣以上面的公代資料為例，最後繪出幼體逐月入添比例：

可以發現公代的幼貝入添比例在圖中呈現一個相當漂亮的單峰，代表每年只有一個入添高峰。

4. SA（simulated annealing）
   如果覺得經典 VBGF 模擬出的效果不理想，或者物種的成長速度有明顯的季節性（例如一些底棲生物），可以考慮帶有季節成長參數的模擬退火法（simulated annealing）可以獲得比較符合經驗的成長曲線。

set.seed(1111)
fit2 <- ELEFAN_SA(lfq, seasonalised = T,
  init_par = fit$par[1:5],
  low_par = list(Linf=PW$confidenceInt_Linf[1], K=0.3, t_anchor=0.7, ts=0, C=0),
  up_par = list(Linf=PW$confidenceInt_Linf[2], K=1.5, t_anchor=1, ts=1, C=1),
  SA_temp = 2e5,
  SA_time = 60,
  maxit = 400, MA = 3,
  plot.score = TRUE,
  verbose = F)
fit2$Rn_max; unlist(fit2$par)
plot(fit2, main="Growth curve of ___species___")
lfqFitCurves(lfq, par = fit2$par,
             draw = TRUE, col = "darkred", lty = 1, lwd=1.5)

5. GA（genetic algorithm）
   遺傳演算法（genetic algorithm）自 1975 年提出後，被應用到許多的函數最佳化問題，在 1990 年代也開始被用來尋找 VBGF 的最佳參數（對於遺傳演算法本身有興趣或可參考這篇）。

set.seed(1111)
fit3 <- ELEFAN_GA(lfq, seasonalised = T,
  low_par = list(Linf=PW$confidenceInt_Linf[1], K=0.3, t_anchor=0, ts=0, C=0),
  up_par = list(Linf=PW$confidenceInt_Linf[2], K=1.5, t_anchor=1, ts=1, C=1),
  popSize = 50, pmutation = 0.2,
  maxiter = 100, run = 20,
  MA = 5, plot.score = TRUE,
  monitor = F,
  parallel = FALSE)
fit3$Rn_max; unlist(fit3$par)
plot(fit3, main="Growth curve of ___species___")
lfqFitCurves(lfq, par = fit3$par,
             draw = TRUE, col = "darkred", lty = 1, lwd=1.5)

以上是在 R 中進行 ELEFAN 與其他漁業資源統計的簡介，一邊摸索一邊做，也還沒有很了解。若有需要修改或者討論的地方歡迎留言。
🐕‍🦺 2022.10.25