Try   HackMD

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

tags: 用R做 生態學

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

𝐿𝑡=𝐿[1𝑒𝐾(𝑡𝑡0)+(CK2𝜋)sin2𝜋(𝑡𝑡𝑠)]
其中
Lt
是年齡為
t
之體長;
L
為極限體長(asymptotic length),也是 VBGF 要逼近的目標;
t0
Lt=0
的理論年齡,一般在推論時,會假設
t0=0
K
為成長參數(growth coefficient);
e
是自然指數;
C
為季節性振蕩幅度參數(seasonal oscillation),與年溫差成正比,溫差 1 度相當於變化 0.1 單位,其值介於 0(無振盪)到 1(最大振盪)之間;
ts
為成長曲線起伏變化的起點,其值視生物棲息地的冬季低迷點(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 提出的方法:
    𝐵𝑖𝑛=RangerangeSD×4

    Bin=Range1+3.322logN

    其中
    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)
  1. 尋找極限體長
    使用兩種方法來尋求極限體長,需要調整極限體長範圍與 cross.datecross.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),又稱公代)體長資料作為輸入值,這步繪出來的生長曲線如下圖:

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

  1. 計算幼體入添比例
    使用 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

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

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

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

  1. 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)
  1. 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