VPA(資源解析手法教科書)
前回は北大水産の実験で使用したVPAの計算式を用いましたが、今回はよりメジャーな「資源解析手法教科書」に載っている計算式を使用したVPAプログラムを載せます。
具体的にいうと、資源解析手法教科書 pp110‐112の4節「実用的なVPA」を参考としました。
一応、教科書の問題3の数値例(p107)と表2(p114)の数値を用いて動作確認を行いましたが、バグ等ございましたらご一報いただけると幸いです。
RでVPA
早速RでVPAしてみましょう。
VPAがどういうモノか知りたい方はこのページを参照してください。
私の自作したVPAプログラムはこちらからダウンロードしてください。右クリックして名前を付けてリンク先を保存するとよいです。
サンプルデータはこちらから落としてください。これは管理人の作成したシミュレーションデータであり、実データではありません。
では、VPAのプログラム「VPA_textbook.txt」とデータファイル「data.txt」が両方とも同じフォルダに格納されているという前提で話を進めていきます。
まずは、Rを起動した後コンソールのファイル→ディレクトリの変更をして、データや計算プログラムの入ったフォルダを指定してください。そのうえで、以下を実行します。
source(“VPA_textbook.txt”)
data <- read.table(“data.txt”, header=T)
一行目は計算プログラムを読み込むという指示で、2行目はデータの読み込みです。
今回はVPA.2という名前の関数を作成しました。私の作ったVPA.2という関数には引数として、データ、自然死亡係数、漁獲の影響力の初期パラメータの3つが必要となります。
これらを突っ込めば一瞬で海の中にいる魚の個体数と漁獲の影響力の経年変化を計算することが可能です。
やってみます。
VPA.2(data, M=0.2, F.first=0.5)
Mが自然死亡係数(自然死亡の大きさ。自然死亡率そのものではない)で、F.firstが漁獲係数(漁獲の影響力。漁獲割合そのものではないので注意)の初期値です。
初期値は名前の通り初期の値なので、データを使って計算しなおした結果を出力してくれるのですが、自然死亡の大きさに関しては入れた値をそのまま使用して計算されてしまいます。Mを変えると結果は結構変わります。これが魚の数を数えることが難しい理由の一つですね。
計算結果は長いのでこちらからダウンロードしてください。
馬場様
VPA勉強中の者です。
分かりやすいコードや記事でいつも参考にさせていただいております。
「VPA(資源解析手法教科書)」のコードについて質問させてください。
載せていただいたコードを手持ちのデータで実行すると、FやNの推定結果がマイナス値になってしまいます。
用意していただいたテストデータや、資源解析手法教科書に記載された表2の数値では問題なくNやFは推定されていました。
元データは2~7歳の34年分の漁獲尾数で、最小値と最大値は0.219~144.392(万尾)です。
Excelを使うと(同じ式で算出しています)プラス値になるので、Rで行っているFの最適化で何か問題が生じていると思い、別の最適化関数を使ってみるなどあれこれと試行錯誤したのですが、改善できませんでした。
この問題について何かお心当たりなどございましたら、教えていただけないでしょうか。
よろしくお願いします。
RT
RT様
コメントありがとうございます。
管理人の馬場です。
返信が遅れてすいません。
当方のコードは何も工夫をせずに最適化しているので、誤った結果が出る可能性はありますね。
Excelとの違いに関しては、データを見てみないとわかりかねるところがあります。
お役に立てず、すいません。
馬場様
お返事いただきありがとうございます。
現状では原因不明とのこと、承知いたしました。
地道に勉強しつつエラーと戦ってみようと思います。
RT