成長曲線

最終更新:2015年11月22日

魚は、大きくなって価値が高くなってから取るのが ベスト。小さいものは獲らず、成長するまで待つことが肝心です。

でも、いったいどれくらい待てばいいの?
雌雄によって待つ時間は変えたほうがいいの? 一緒でいいの?

この問いに答えるために、魚の成長の特徴を調べ る。これが成長曲線の推定と言う代物です。

サンプルデータとRパッケージ

今回もこまかい理屈は置いておいて、Rを使ってvon Bertalanffyの成長曲線を推定します。さらに、雌雄によって成長曲線が異なるとみなせるかを検定するところまで行います。

今回はパッケージFSAを用います。

コードをコピペする際はこちらのソースコード一式をご利用されると便利です。

FSAは「パッケージのインストール」からは落とせません。FSAをインストールするには、以下のコードを実行します。ただし、やや時間がかかります。

基本的にはこのページに書かれていることだけで事足りますが、細かいことは
http://www.rforge.net/FSA/Installation.html
も参考にしてください(英語です)。

パッケージのインストールの手順を説明します。
まずは、Rを右クリックして「管理者として実行」します。そのうえで以下のコードを実行。

source(“http://www.rforge.net/FSA/InstallFSA.R”)

上記のコードはVPAのページでも使ったsource関数を使っています。ウェブサイトにUPされているRファイルをよみこんで実行してね、という指示です。

もしもうまくいかなかった場合はhttp://www.rforge.net/FSA/InstallFSA.Rから直接スクリプトをコピペしてRのエディタに張り付けてから実行してください。

さらに、データを読み込むために以下を実行する必要もあります。

install.packages(“FSAdata”,repos=”http://www.rforge.net/”,type=”source”)

パッケージをインストールする際に国の名前が表示された場合は、Japanならどこでもよいので3つくらいあるJapanのなかから好きなものを選んでダブルクリックしてください。そうしたら、パッケージのインストールが完了します。

パッケージのインストールが完了したら、以下のコードを実行して、パッケージを使用する下準備をします。

library(FSA)
library(FSAdata)
library(nlstools)

これでようやっと計算に映ることができます。
FSAパッケージは、私はあまり使ったことはありませんが水産資源解析関連のいろいろな関数があるようなので、興味のある方は調べてみてください。

パッケージFSAdataの中に入っているCroaker2と言うデータを用います。
どんなデータか見てみます。

 # データの読み込みとデータの表示
data(Croaker2)
head(Croaker2)

# 性別ごとの体長データの要約
tapply(Croaker2$tl,Croaker2$sex,summary)

# 性別ごとに分けた年齢体長関係
plot(tl~age,pch=c(21,19)[Croaker2 $sex],col=c(1,2)[Croaker2$sex],
data=Croaker2,main=”雌雄別年齢体長プロット”)
legend(“topleft”,legend=c(“F”,”M”),pch=c(21,19),col=c(1,2))

結果はこちら

> head(Croaker2)
age tl sex
1 1 243 F
2 1 247 F
3 1 248 M
4 2 330 F
5 2 320 F
6 2 285 F
> tapply(Croaker2$tl,Croaker2$sex,summary)
$F
Min. 1st Qu. Median Mean 3rd Qu. Max.
200.0 326.5 352.0 358.9 410.0 496.0

$M
Min. 1st Qu. Median Mean 3rd Qu. Max.
210.0 292.0 325.0 319.6 343.8 462.0

雌雄別プロット

tapplyの結果や、グラフを見ると、雌の方が若干大きそうな感じがします。

と言うことで、雄雌一緒にしたデータと、雌雄別のデータの両方を用いて、成長曲線がどう変わるか比較していきましょう。雌雄で分けるためのプログラムはこちら。

Croaker.M <- subset(Croaker2,sex==”M”)
Croaker.F <- subset(Croaker2,sex==”F”)

データの準備ができたので、次からは実際に成長曲線を推定していきましょう。

 

von Bertalanffy 曲線の推定

上のデータを用いて、代表的な成長曲線の一つである von Bertalanffy 曲線を推定してみます。

# 初期値の設定 FSAパッケージより
start.full <- vbStarts(tl~age,data=Croaker2,type=”typical”)
start.full

# 雌雄まとめたモデル
Bertalanffy.full<-nls(tl~Linf*(1-exp(-K*(age-t0))),data=Croaker2,start=start.full)
overview(Bertalanffy.full)

#雌雄別モデル
Bertalanffy.M<-nls(tl~Linf*(1-exp(-K*(age-t0))),data=Croaker.M,start=start.full)
Bertalanffy.F<-nls(tl~Linf*(1-exp(-K*(age-t0))),data=Croaker.F,start=start.full)

# プロット
plot(tl~age,xlab=”年齢”,ylab= “全長 (mm)”,pch=16,col=c(4, 2)[Croaker2$sex],
data=Croaker2,main=”雌雄別年齢 体長プロット”)
lines(c(1,2,3,4,5,6,7,8,9,10),predict(Bertalanffy.full,
new=data.frame(age=c (1,2,3,4,5,6,7,8,9,10))),lwd=2)
lines(c(1,2,3,4,5,6,7,8,9,10),col=2,predict(Bertalanffy.M,
new=data.frame(age=c (1,2,3,4,5,6,7,8,9,10))),lwd=2)
lines(c(1,2,3,4,5,6,7,8,9,10),col=4,predict(Bertalanffy.F,
new=data.frame(age=c (1,2,3,4,5,6,7,8,9,10))),lwd=2)
legend(“bottomright”,legend= c(“合わせた”,”雄”,”雌”),lwd=2,col=c(1,2,4))

文献3を見ると、重み付き最小二乗法を使った方がよさそうな気がしましたが、今回は重みを設定していません。重み付き最小二乗法を するならば、nls()の中にweightsを指定してやる必要があります。

成長曲線は、非線形のモデルに最小二乗法を当てはめて推定するわけですが、毎回初期値をどうしようかと言うことで悩みます。そこで、vbStarts(tl~age,data=Croaker2,type=”typical”)を使いました。

これは、ワルフォードの定差図を使ってパラメータを推定する関数です。ワルフォードの定差図による推定は統計的な誤りがあるので本来は使ってはなりませんが、非線形モデルを推定するための初期値にこれを使ったわけです。FSAの説明書(http://www.rforge.net/doc/packages/FSA/vbStarts.html
にもこの初期値は余り信用しすぎるなと書かれてあります。この値だけを信用すべきと言う確固たる理由は特にないようです。でも、正確な値と割と近そうな初期値を選べるということで、今回も採用しました。

今回は雌雄を差別しないモデル以外では省略しましたが、overview (モデル)でモデルの詳細を見ることができます。こちらはパッケージnlstoolsの関数で、パラメータの95%信頼区間な ども表示してくれます。

overviewの結果はこんな感じ。

> overview(Bertalanffy.full)

——
Formula: tl ~ Linf * (1 – exp(-K * (age – t0)))

Parameters:
Estimate Std. Error t value Pr(>|t|)
Linf 416.25682 19.34549 21.517 < 2e-16 ***
K 0.24180 0.06242 3.874 0.00013 ***
t0 -2.17074 0.85673 -2.534 0.01177 *

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 45.33 on 315 degrees of freedom

Number of iterations to convergence: 7
Achieved convergence tolerance: 5.744e-06

——
Residual sum of squares: 647000

——
Asymptotic confidence interval:
2.5% 97.5%
Linf 378.1941202 454.3195283
K 0.1189993 0.3646082
t0 -3.8563825 -0.4851029

——
Correlation matrix:
Linf K t0
Linf 1.0000000 -0.9637788 -0.8779838
K -0.9637788 1.0000000 0.9674683
t0 -0.8779838 0.9674683 1.0000000

最後のプログラムで、以下のような図が描けます。

成長曲線

 

雌雄で成長曲線に差があるか、検定

F検定で、成長曲線の雌雄差の有無を検定します。

Smf<-deviance(Bertalanffy.full)
Sf<-deviance(Bertalanffy.M)
Sm<-deviance(Bertalanffy.F)

F<- (Smf-Sm-Sf)/3 / ( (Sm+Sf)/(length(Croaker2$age)-3*2))

 

> F
[1] 20.24402

> 1-pf(F,3,length(Croaker2$age)-3*2)
[1] 5.166978e-12

RではExcelと異なり、推定されたモデルのパラメータをoverview関数で細かく見れます。また、F検定も比較的楽にできる のではないかと思います。

尚、重み付き最小二乗法で推定した場合には、F検定で無くχ2検定を行う必要があるようです。

 

文献

赤嶺達郎:水産資源解析の基礎, 恒星社厚生閣, 2007年

FSAのページ(FSA – Fisheries stock assessment methods.
http://www.rforge.net/FSA/index.html

 

 

成長曲線” に対して6件のコメントがあります。

  1. 吉倉 桂三 より:

    お世話になります。私は、富山県立滑川高等学校海洋科 教諭吉倉桂三と申します。
    当方、Rを使用しようと思っている超初心者です。
    飼育するサクラマスの成長曲線(von Bertalanffy)を作成しようとこちらのサイトを見つけました。
    Rは先日インストールしたVer.3.1.1です。

    現在、Rの使用もままならない段階で、無謀なほど自身の知識も稚拙です。
    こちらのサイトを使用し練習をと思ったのですが、最初から躓いています。
    こちらの案内のようにすすめ、最終的にグラフまで行ければ手ごたえもあるかと思っていました。

    (本文)
    パッケージのインストールの手順を説明します。
    まずは、Rを右クリックして「管理者として実行」します。そのうえで以下のコードを実行。

    source(“http://www.rforge.net/FSA/InstallFSA.R”)

    上記のコードはVPAのページでも使ったsource関数を使っています。ウェブサイトにUPされているRファイルをよみこんで実行してね、という指示です。

    もしもうまくいかなかった場合はhttp://www.rforge.net/FSA/InstallFSA.Rから直接スクリプトをコピペしてRのエディタに張り付けてから実行してください。
    (ここまで)

    ですでになにも起きません。

    Rコンソール画面に > のあとにsource…で行いました。
    するとすぐ下に > が現れたので、

    install.packages(“FSAdata”,repos=”http://www.rforge.net/”,type=”source”)

    も入れましたが、また > です。国名も聞いてきません。時間がかかるとのことでしたが、何も起きません。
    そもそもしていない操作があるのでしょうか?

    このサイトで手ごたえを得られれば、手持ちのデータでできればと考えていました。
    ご面倒をおかけし申し訳ございませんが、
    この状態で前に進める方法を教えてください。

    1. 馬場真哉 より:

      吉倉 桂三さま
      管理人の馬場です。

      サンプルコードが回らないということで確認したところ、確かに動きませんでした。
      原因はダブルクォーテーションマークのにあるようです。

      ……私の打ち間違いかと思ったら、訂正した後またもとに戻ってしまいました。
      ワードプレスの仕様なのかもしれません。

      ダブルクォーテーションマーク(”)を自分で打ち込んでいただければ、おそらく動くのではないかと思います。
      半角にしてからShift + 2をすれば正しいダブルクォーテーションが出てくるはずです。

      こちらの手違いで回らなくなってしまい、申し訳ございませんでした。
      お詫びいたします。

      また何か問題があればコメントを頂けると幸いです。

  2. 馬場真哉 より:

    確認したところ、全てのコードにおいてダブルクォーテーションマークがおかしなことになってしまっていました。

    サンプルコードをテキストファイルとして本文にリンクを張っておきましたので、こちらをご参照いただけると幸いです。
    こちらはそのままコピペして大丈夫なはずです。

    ご迷惑をおかけして、誠に申し訳ありませんでした。
    貴重なご意見、ありがとうございました。

  3. 駒澤一朗 より:

    馬場様。
    いつも著書とサイトで勉強させていただいております。
    また、質問させていただきたいのですが、手持ちのデータで雌雄別に成長曲線を求めたところ、Smfが10170.37、Sfが4374.778、Smが5134.005となりました。また、このデータのageのデータ数は1380でした。
    これらからFを求めると31.86614となり、さらにp値を求めると0と計算されるのですが、この結果は雌雄間で成長曲線に差があるということでいいのでしょうか。
    お忙しい中、申し訳ありませんが、ご教授いただけると助かります。

    1. 馬場真哉 より:

      駒澤様

      管理人の馬場です。
      すいません、コメントに気が付かず、返信が大変遅れました。
      p値がほぼ0ということであれば、有意水準を0.05とおいた場合、これを下回っていることになりますので、成長曲線に有意な差があると主張できることになります。

      1. 駒澤一朗 より:

        馬場様。
        お答えいただき、ありがとうございます。
        今後ともよろしくお願いいたします。

コメントを残す

メールアドレスが公開されることはありません。 が付いている欄は必須項目です

このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください