1-1 Construct a table of the number of fish in each age and 1 inch TL category

載入需要套件

library(dplyr)
library(tidyverse)
library(tibble)
library(graphics)

匯入Stripe Bass資料，查看資料長相與散布圖

SB <- read.table("C:\\Users\\user\\Desktop\\台灣大學\\族群生態與漁業資源永續\\W4\\StripedBass3.txt", header = T, sep = ",")
str(SB)
## 'data.frame':    1201 obs. of  2 variables:
##  $ tl : int  18 18 18 18 18 18 18 18 18 18 ...
##  $ age: int  2 2 2 2 3 3 4 4 5 7 ...
plot(SB[2:1], xlab = "Age", ylab = "Length", pch = 20, main = "Age-Length of Stripe Bass")

1-2 Construct an observed Age-Length Key Table

去除NA值，計算年紀對應體長的count

SB <- SB %>% na.omit() %>% group_by(tl, age) %>% summarise(count = n())

生成Age-length key table，Row為Length，Column為Age

Age_Length <- SB %>% pivot_wider(names_from = age, values_from = count, values_fill = 0) %>% as.data.frame()
rownames(Age_Length) <- Age_Length$tl
Age_Length <- Age_Length[-1]
Age_Length

1-3 Proportion of fish in the 30 in TL category should be assigned age 10

計算體長為30，對應到的10歲的比例，由於體長30對應10歲的魚並無資料，因此比例為0。

Por <- matrix(nrow = dim(Age_Length)[1], ncol = dim(Age_Length)[2])
for (i in 1:30) {
  Por[i, 1:17] <- round(as.numeric(Age_Length[i, ] / sum(Age_Length[i, ])), digits = 1)
}
row.names(Por) <- rownames(Age_Length)
colnames(Por) <- colnames(Age_Length)

Por[30, ]
##  2  3  4  5  7  6  8  9 10 11 13 12 14 15 18 19 17 
##  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0

2-1 Estimate L∞ and K using the Ford-Walford method

匯入Pacific Bluefin Tuna資料，查看資料長相與散布圖

setwd("C:\\Users\\user\\Desktop\\台灣大學\\族群生態與漁業資源永續\\W4")
PBF <- read.csv("PBF_Len_Age.csv")
str(PBF)
## 'data.frame':    2381 obs. of  2 variables:
##  $ Age   : int  16 7 9 16 24 16 19 8 8 12 ...
##  $ Length: int  235 204 203 235 252 230 251 190 217 242 ...
plot(PBF, pch = 20, main = "Age-Length of Pacific Bluefin Tuna")

計算Age對應AveRage Length，建立表格

Avg_PBF <-aggregate(Length ~ Age, data = PBF, FUN = mean)
Avg_PBF

計算Length t與Length t+1的Intercept & slope

FW_table <- data.frame(Age = rbind(0, Avg_PBF[1]),
                         Length = rbind(0, Avg_PBF[2]),
                         Age_1 = rbind(Avg_PBF[1], 0),
                         Length_1 = rbind(Avg_PBF[2], 0))
FW_table <- FW_table[-c(1, 26), ]

FW_linear <- lm(Length.1 ~ Length, data = FW_table)
#summary(FW_linear)
coef(FW_linear)
## (Intercept)      Length 
##  43.1504350   0.8253364

估計Linf與K，分別為246.9971與0.1920083

#Linf = Intercept / 1-Slope
LW_Linf <- 43.1504 / (1 - 0.8253)
print(LW_Linf)
## [1] 246.9971
#K = -ln(Slope)
LW_K <- -log(0.8253)
print(LW_K)
## [1] 0.1920083
plot(FW_table[c(2, 4)], pch = 20, ylab = "Length+1", main = "Length of Pacific Bluefin Tuna")
abline(FW_linear, col = "red", lwd = 2)
text(200,230, "Lt+1 = 0.8253*Lt + 43.1504", col = "red")

2-2 Estimate t0 with information: 50 cm FL at age 1

將Age為1時，Length為50cm之資訊帶入，獲得t0為-0.1780091

LW_t0 <- 1 + (1/LW_K)*(log((LW_Linf - 50) / LW_Linf))
print(LW_t0)
## [1] -0.1780091

3-1 Fit a Schnute growth function to the length-at-age data using nonlinear fitting in R.

t1為sample中minimum Age : 4
t2為sample中maximum Age : 28
L1為sample中minimum Age的對應Average length : 172.6667
t1為sample中maximum Age的對應Average length : 246.5

L1 <- Avg_PBF[1, 2]
L2 <- Avg_PBF[25, 2]
t1 <- Avg_PBF[1, 1]
t2 <- Avg_PBF[25, 1]

帶入Schnute growth function，估計K值為0.1741895

Schnute <- function(K, Age){
  Length <- L1 + (L2 - L1) * (1 - exp(-K * (Age - t1))) / (1 - exp(-K * (t2 - t1)))
  return(Length)
}

fit <- nls(Length ~ Schnute(K, Age), 
            start = list(K = 0.2), 
            data = Avg_PBF, 
            trace = F)
#summary(fit)
K <- coef(fit)
print(K)
##         K 
## 0.1741895

3-2 Estimated growth coefficients (L∞, K, and t0 )

將估值K(0.1741895)，已知t1, t2, L1與L2帶入估計Linf與t0
Linf與t0分別為247.6465與-2.859104

Linf <- (L2 - L1 * exp(-K * (t2 - t1))) / (1 - exp(-K * (t2 - t1)))
print(Linf)
##        K 
## 247.6465
t0 <- t1 + (1 / K) * log((L2 - L1) / (L2 - L1 * exp(-K * (t2 - t1))))
print(t0)
##         K 
## -2.859104

3-3 Scatter plot of length(y)-age(x) including the fitted growth curve

將利用Schnute growth function估計出的K, Linf, t0帶回VBGF，計算出fitted growth curve

plot(Avg_PBF, pch = 20, main = "Schnute estimated growth curve")
lines(Avg_PBF$Age, predict(fit, Avg_PBF$Age), lwd = 2, col = "blue")

4-1 Based on (3), fix 𝑡0 at zero and re-estimate the growth parameters L∞ and Kc

建立VBGF,並將t0固定為0

t0 <- 0
VBGF <- function(Linf, K, Age) {
  Length <- Linf*(1 - exp(-K*(Age - t0)))
  return(Length)
}

重新估計Linf與K值，分別為244.0035130與0.2694651

fit2 <- nls(Length ~ VBGF(Linf, K, Age), start = list(Linf = Linf, K = K), data = Avg_PBF, trace = F)
#summary(fit2)
coef(fit2)
##        Linf           K 
## 244.0035130   0.2694651

4-2 Compare the differences in the growth curves of (3) & (4)

藍線 : Schnute method將估計資料鎖定於資料Age範圍內，因此不會估計到4歲以下的Length數值。
紅線 : 若將t0固定為0，估計值會延伸到t0，造成轉折點有高估Length之現象(納入錯估之Age-Length資料)。

plot(Avg_PBF, pch = 20, main = "Schnute estimated growth curve at estimated t0 and fixed t0")
lines(Avg_PBF$Age, predict(fit, Avg_PBF$Age), lwd = 2, col = "blue")
lines(Avg_PBF$Age, predict(fit2, Avg_PBF$Age), lwd = 2, col = "red")
legend(20, 220, legend = c("Estimated t0 line", "Fixed t0 line"),
       col = c("blue", "red"), lty = 1:1,cex = 1)