はじめに

本日は、非線形モデルの一種であるvon Bertalanffy(フォン・ベルタランフィ)の成長曲線1を扱っていく。
このモデルは、簡潔に述べると、年齢ごとの魚の体重や体長の変化を表している。
von Bertalanffyの体長の成長曲線の公式は、以下の下記の通りである。

\[ L_t=L_∞(1-e^{-K(t-t0)})\tag{1} \]

式で使用している記号が表すものは、次のとおりである。
Ltは「ｔ時点(年齢時)の推定体長」、L∞は「最大到達体長」、
Kは「成長係数」、tは「年齢」、t0は「体長が０となる時の年齢」である。

今回は、乱数を生成させ、von-Bertalanffyの体長の成長曲線を導出する。

データ生成

初めに、年齢と体長の乱数を生成する。 今回は、1歳魚から6歳魚を採取し、下記のデータを得られたと仮定する。2

#yo1 <- data.frame(age=rep(1,1),tl=rnorm(1,246,8.6))
#yo2 <- data.frame(age=rep(2,5),tl=rnorm(5,511,2))
#yo3 <- data.frame(age=rep(3,10),tl=rnorm(10,584,1.4))
#yo4 <- data.frame(age=rep(4,160),tl=rnorm(160,687,7.5))
#yo5 <- data.frame(age=rep(5,19),tl=rnorm(19,770,2.2))
#yo6 <- data.frame(age=rep(6,5),tl=rnorm(5,802,1.5))

yoは年齢ごと体長の乱数を生成させ、データフレームにしたものを格納したものである。
yoの横の数字は、年齢である。例えば、yo1は「1歳魚(1年魚)」である。
yoには、“age”という「年齢」と“tl”という「体長」の2つの変数が入っている。
1つ紹介すると、yo3では体長の平均が584(mm)、分散が1.4の乱数を10個生成している。
更に、yo3は3年魚であるため、rep関数で10個“3”を生成している。

各年齢ごとのデータフレームを1つにまとめると以下の通りになる。
乱数だと出力した時に値が変わるため、csvに書き出して読み込ませる。 

#dat <- rbind(yo1,yo2,yo3,yo4,yo5,yo6)
#write.csv(dat,"vb.csv")
library(tidyverse)
dat2 <- read_csv("vb.csv")
library(DT)
DT::datatable(dat2)

datというデータフレームには、以上のように、2つの変数が200個入っている。 

summary(dat2)
      age             tl       
 Min.   :1.00   Min.   :262.7  
 1st Qu.:4.00   1st Qu.:680.7  
 Median :4.00   Median :687.0  
 Mean   :4.03   Mean   :685.9  
 3rd Qu.:4.00   3rd Qu.:694.6  
 Max.   :6.00   Max.   :803.5  
big <-max(dat2$tl)
big
[1] 803.4783
plot(dat2$age,dat2$tl,xlab = "年齢",ylab = "体長")

hist(dat2$tl,xlab ="体長")
上記の図より、600㎜台が多く採取することができた。

分析

初めに、von Bertalanffyの成長曲線を再提示する。 \[ L_t=L_∞(1-e^{-K(t-t0)})\tag{1} \]

この数式では、L∞、K、t0の3つのバラメータがある。
よって、L∞をa、Kをb、t0をcと置いて分析した。 

In = nls(tl ~ a*(1-exp(-b*(age - c))),data = dat2,
         start = list(a = big ,b = 0.3,c = 0.3))
summary(In)

Formula: tl ~ a * (1 - exp(-b * (age - c)))

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
a 994.13098   21.80000  45.602  < 2e-16 ***
b   0.26290    0.01590  16.531  < 2e-16 ***
c  -0.47810    0.09134  -5.234 4.23e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 10.7 on 197 degrees of freedom

Number of iterations to convergence: 4 
Achieved convergence tolerance: 3.368e-07

a=994.13091、b=0.26290、C=-0.47810とパラメータが推定されている。 したがって、(1)に代入すると、

\[ L_t=994.13091(1-e^{-0.26290(t+0.47810)})\tag{2} \] となる。

次に、年齢ごとに推定された体長を取り出し、データフレームに格納する。

predict(In,newdata = dat2)
  [1] 320.0989 475.9220 475.9220 475.9220 475.9220 475.9220 595.7218 595.7218
  [9] 595.7218 595.7218 595.7218 595.7218 595.7218 595.7218 595.7218 595.7218
 [17] 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262
 [25] 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262
 [33] 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262
 [41] 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262
 [49] 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262
 [57] 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262
 [65] 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262
 [73] 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262
 [81] 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262
 [89] 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262
 [97] 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262
[105] 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262
[113] 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262
[121] 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262
[129] 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262
[137] 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262
[145] 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262
[153] 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262
[161] 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262
[169] 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262 687.8262
[177] 758.6379 758.6379 758.6379 758.6379 758.6379 758.6379 758.6379 758.6379
[185] 758.6379 758.6379 758.6379 758.6379 758.6379 758.6379 758.6379 758.6379
[193] 758.6379 758.6379 758.6379 813.0793 813.0793 813.0793 813.0793 813.0793
d <- data.frame(predict(In))
dat3 <- data.frame(dat2,d)
DT::datatable(dat3)
age <- c(1,2,3,4,5,6)
tl_mm <- c(320.0989,475.9220,595.7218,687.8262,758.6379,813.0793)  
hi  <- cbind(age,tl_mm)
datatable(hi)

上記の表では、年齢と年齢ごとの体長の推定値を示している。
ageは年齢で、tl_mmは体長の推定値である。
体長は、1歳魚が320.0989mm、2歳魚が475.9220mm、3歳魚が595.7218mm、
4歳魚が687.8262mm、5歳魚が758.6379mm、6歳魚が813.0793mmと推定された。

この体長の推定値を図示する。 

plot(dat3$age,dat3$predict.In.,type="l",col="green",lwd=3
     ,xlab = "年齢",ylab = "体長(mm)",ylim=c(100,800))

x軸は年齢で、y軸を推定された体長とした。 少しであるが、弧を描いていることがわかる。

これを観測値をプロットした図に当てはめて比較する。 

plot(dat2$age,dat2$tl,xlab = "年齢",ylab = " ",type="b",ylim=c(100,800)
     ,lwd=4)
lines(dat2$age,fitted(In),col=2,lty=2,lwd=10,ylim=c(100,800))
par(new=T)
plot(dat3$age,dat3$predict.In.,type="l",col="green",lwd=3
     ,xlab = "年齢",ylab = "体長(mm)",ylim=c(100,800))
黒の実線が観測値、赤の点線がvon Bertalanffyの式の値、 青の実線が推定された体長の値である。
赤の点線と青の実線は同じ曲線になることは自明であるが、プロットしてみた。

このように、von Bertalanffyの成長曲線で想定されていた非線形のモデルを再現することができた。

おわりに

このように、1歳魚が約320mm、2歳魚が約476mm、3歳魚が約596mm、4歳魚が約688mm、5歳魚が約759mm、6歳魚が約813mmと体長が推定できた。友人によると、1歳魚はなかなか調査できないため不明だが、そのほかの年齢別体長は推定値に近いと言われた。
今回は、乱数を生成させ、仮想の魚の採取データを分析している。したがって分散やサンプルサイズなどは、実際の調査とかけ離れている可能性が高いことを最後にお伝えしたいと思う。

今回は、von Bertalanffyの成長曲線を扱ったが、非線形モデルは他にも存在するため、他のモデル統計ソフトで分析していきたい。

脚注

  1. “VBモデル”や“VBLM”と表記されていることもある。↩︎

  2. サンプルサイズ、平均、分散は、魚に詳しい友人に教えていただいた。↩︎