Diagnostics for Stock Recruitment Relationships - ShotaNishijima/frasyr GitHub Wiki

• 事前準備
• 通常の再生産関係（レジームなし）
  • 再生産関係の推定とプロット
  • 残差の分布のチェック
  • 残差の自己相関のチェック
  • ブートストラップ
  • ジャックナイフ
  • プロファイル尤度
• レジームありの再生産関係
  • 再生産関係の推定とプロット
  • 残差の分布のチェック
  • 残差の自己相関のチェック
  • ブートストラップ
  • ジャックナイフ
  • プロファイル尤度

再生産関係の診断の手順の一例をマイワシ対馬暖流系群を使って紹介します
診断の手法としては

• 残差の確率分布のチェック - 残差の自己相関のチェック - ブートストラップ
• ジャックナイフ - プロファイル尤度
  の5種類を実装しています
  通常の再生産関係にも、レジーム分けをした再生産関係にも適用できます

frasyrのインストール・読み込み

devtools::install_github("ichimomo/frasyr@dev")
library(frasyr)
library(stringr)

rvpaデータの読み込み

data(res_vpa)
SRdata = get.SRdata(res_vpa)

まず、レジームを考慮しない通常の再生産関係についての診断について説明します

ここではHS再生産関係④の①L2で自己相関ナシ、②L1で自己相関あり、③L1で自己相関を外側推定、④L2で自己相関を内側推定の4つのモデルを検討します

resL2 = fit.SR(SRdata, SR = "HS", method = "L2", out.AR = FALSE, AR = 0, rep.opt = TRUE)
resL1 = fit.SR(SRdata, SR = "HS", method = "L1", out.AR = FALSE, AR = 0, rep.opt = TRUE)
resL1outer = fit.SR(SRdata, SR = "HS", method = "L1", out.AR = TRUE, AR = 1, rep.opt = TRUE)
resL2inner = fit.SR(SRdata, SR = "HS", method = "L2", out.AR = FALSE, AR = 1, rep.opt = TRUE)

例として④をプロットします：

(g1 = SRplot_gg(resL2inner))

• check.SRdistという関数を使って、残差の分布が仮定にあっているかをチェックできます
• 図が3つ出力され、左が標準化された（標準偏差で割られた）残差のヒストグラムで赤い線が標準正規分布を表します。右上のSWはShapiro-Wilk検定により、残差が正規分布に従っているかを検定したときのP値です。右上のKSはKolmogorov-Smirnov検定により、残差が正規分布に従っているかを検定したときのP値です。どちらも、帰無仮説は「正規分布に従っている」です
• 真ん中は残差の累積確率密度のヒストグラムを示しており、理論的には0~1の一様分布に従うので、このヒストグラムが大きく凸凹しているようだと仮定した確率分布と合致していないことになります
• 右も累積確率密度を0~1の一様分布と仮定したときのQQプロットで、線は理論値を表しており、この線から外れていないほうがよいということになります。注意すべき点は、このQQプロットは一様分布を仮定しているので、端は線の載りやすく、中心付近が外れやすいことです。つまり、qqnorm等で出力される正規分布を仮定したQQプロットとは単純な比較はできません。これも、L1と比較するためにこういう処理をしています

check.SRdist(resL2)

output=TRUEとすれば、pngファイルが出力され、保存されます

check.SRdist(resL2, output = TRUE, filename = "ResidDistCheck_L2")

• L1を使う場合も同様の図が描かれますが、左の図に青い点線がプロットされます。これはラプラス分布を表しており、L1の場合の尤度はラプラス分布を使って計算されるため示してあります。右上のKS (Kolmogorov-Smirnov検定) のP値もラプラス分布との一致性を検定しています。SWの方は正規分布にのみ適用できる検定なので、正規分布との一致性を検定しています。L1推定の場合は正規分布を仮定しているわけではないので、この値は参考程度でいいと思います
• 真ん中と右の図は、ラプラス分布を仮定したときの累積確率密度のヒストグラムとQQプロットになります

check.SRdist(resL1)

自己相関を外側推定した場合にはファイルが2つ出力され、1つ目が再生産関係の推定、2つ目が自己相関の推定における残差です。しがたって、1枚目の出力は上と同じになります。自己相関の推定はmethodがL1/L2によらず正規分布を仮定しています

check.SRdist(resL1outer)

自己相関を内側で推定した場合には再生産関係と自己相関係数を同時に推定するので、自己相関を除いた残差に対する図が出力されます

check.SRdist(resL2inner)

• 残差の自己相関を事後的に推定するための関数calc.residARを作りました。fit.SRのout.ARオプションとやっていることは同じですが、AICの差だけではなく尤度の値やAICc等の値も取り出せます
• 後述しますが、fit.SRregimeでも外側で自己相関を推定できます
  ※このチェックは必須ではないです

outer1 = calc.residAR(resL1)
outer1 = calc.residAR(resL1)
rbind(outer1$pars, resL1outer$pars)

           a        b       sd        rho
1 0.02786598 53372.54 0.259883 0.02080394
2 0.02786598 53372.54 0.259883 0.02080394

c(diff(outer1$AIC), diff(resL1outer$AIC.ar))

   AR(1)        1 
1.986945 1.986945 

outer1$AICc

   AR(0)    AR(1) 
6.441233 8.729765 

• plot.autocorは残差の自己相関に関するプロットを3つ出力します。左は再生産関係との残差の時系列を表しており、赤い線は平滑化された曲線です。この線が0から外れていると残差がホワイトノイズではなく、トレンドがあると言えます
• 真ん中はラグを1,2,3,…年と増やした場合の自己相関係数を表しています。青い線は95％信頼区間を示しています。1年のラグで有意なときは自己相関の推定を検討したほうがよいと考えられます
• 右側はLjung-Box検定におけるP値です。Ljung-Box検定の帰無仮説は「ラグが1~mまでの全ての自己相関が0である」であり、対立仮説は「ラグ1から mまでの自己相関のうち、少なくとも一つが0でない」となります (参考: https://to-kei.net/time-series-analysis/hypothesis-testing/ )。この場合、全ての点で有意（赤点で表示される）となっているのは、ラグ1の相関係数が高いからだと考えられます

plot.autocor(resL1outer)

上は再生産関係との残差についてチェックをしましたが、自己相関を推定した後での残差についても、use.resid=2とすれば同様の図を描くことができます (デフォルトはuse.resid=1)。

plot.autocor(resL1outer, use.resid = 2)

output=TRUEとすればpngファイルが出力されます

plot.autocor(resL1outer, output = TRUE, filename = "devianceAR")
plot.autocor(resL1outer, use.resid = 2, output = TRUE, filename = "residAR")

• boot.SRには3種類のブートストラップ方法を実装しており、うち2つが残差のブートストラップで、残り一つがデータのブートストラップです。残差ブートストラップには、パラメトリックとノンパラメトリックの2つあります
• 残差のパラメトリックブートストラップは、fit.SRのpars$sdを標準偏差とする正規分布からランダムに乱数を発生させ、予測値からのずれを加えて加入量のブートストラップデータを生成し、再推定しています
• 残差のパラメトリックブートストラップはmethod="p"で実行可能で、図のプロットにはplot.bootSRを使用する
• ここでは100回でやっていますが、最終的に使用する予定の再生産関係には数を増やしたほうがよいと思われます（500回とか1000回）
• 自己相関を推定していない場合は最後のrhoの図は表示されません

boot.res1 = boot.SR(resL1outer, n = 100, method = "p")
plot.bootSR(boot.res1)

• 残差のノンパラメトリックブートストラップは残差の確率分布を仮定せず、残差を重複ありでリサンプリングして、加入量のブートストラップデータを生成します
• 確率分布を仮定しないので、L1とL2を比較する場合にはパラメトリックよりもこちらの方がよいかもしれません
• ノンパラメトリックブートストラップはmethod="n"で行います

boot.res2 = boot.SR(resL1outer, n = 100, method = "n")
plot.bootSR(boot.res2)

• データのブートストラップはデータを重複ありでリサンプリングしたデータを使用して、再生産関係の再推定を行います
• 親魚量データもリサンプリングにより変化するため、親魚量の不確実性も考慮されることになります
• 親魚量データに偏りがあったり、データ数が少なかったり、あるデータ点に推定値が大きく依存している場合はバイアスや不確実性が大きくなりやすいと思われます
• データブートストラップはmethod="d"で行います
• この場合データブートストラップの方が残差ブートストラップよりもバイアスが小さいのは親魚量データのバイアスは小さく、加入変動が大きいからでしょうかね…

boot.res3 = boot.SR(resL1outer, n = 100, method = "d")
plot.bootSR(boot.res3)

図を画像ファイルとして保存する場合はoutput=TRUEとします

plot.bootSR(boot.res1, filename = "para_boot", output = TRUE)
plot.bootSR(boot.res2, filename = "nonpara_boot", output = TRUE)
plot.bootSR(boot.res3, filename = "data_boot", output = TRUE)

データを一点ずつ除いて再推定し、各データが与える影響を評価します

jack1 = jackknife.SR(resL1outer, is.plot = TRUE)

図を画像ファイルとして保存する場合：

jack1 = jackknife.SR(resL1outer, output = TRUE, is.plot = TRUE)

• 再生産関係のパラメータa,bを変化させたときの尤度を計算し、点推定値以外に、尤度が高くなっている山が無いかをチェックします
• 自己相関を内部推定している場合はa,bを固定して尤度が最も高くなる自己相関係数rhoを推定しています 　
• HSで自己相関を内部推定すると複数の山が現れることが知られていますので、HSや自己相関の内部推定を検討する場合は重要な診断となります

prof1 = prof.likSR(resL1outer)

図を保存する場合:

prof1 = prof.likSR(resL2inner, output = TRUE)

次に、レジームで再生産関係を分けて推定した場合について説明します

• 推定にはfit.SRregimeを使います。ここではHS, L1で、2005年にレジームシフトが生じた場合を候補とします
• 仮に1995年と2005年に2回レジームシフトが起こったとして、レジームがA->B->Aと変化する場合はregime.key = c(0,1,0)、A->B->Cと変化する場合はregime.key = 0:2と指定します
• regime.year = c(1995,2005)と指定してください
• regime.parはレジームごとに異なるパラメータを設定します。今はすべてのパラメータが異なる場合 (resR1) とaは共通でbとsdが異なる場合 (resR2) を使います
• AICc等の結果はfit.SRと比較できます
• 詳しくは?fit.SRregimeをご覧ください

resR1 <- fit.SRregime(SRdata, SR = "HS", method = "L1", regime.year = c(2005), regime.par = c("a", 
    "b", "sd")[1:3], use.fit.SR = TRUE, regime.key = c(0, 1))
resR2 <- fit.SRregime(SRdata, SR = "HS", method = "L1", regime.year = c(2005), regime.par = c("a", 
    "b", "sd")[2:3], use.fit.SR = TRUE, regime.key = c(0, 1))
c(resL2$AICc, resL1$AICc, resL2inner$AICc, resR1$AICc, resR2$AICc)

[1] 10.79385 14.07974 13.40053 11.29241 19.24803

(g2 = SRregime_plot(resR1, regime.name = c("Low", "High")))

残差の確率分布のプロットはfit.SRと同様に出力できます

check.SRdist(resR1)

• レジームを推定してから自己相関を事後的に推定することができます
• calc.residARを使い、per_regime = TRUEでレジームごとに自己相関を推定、per_regime = FALSEでレジーム間で共通の自己相関係数を推定できます
• ここで自己相関係数が大きかったりAICcが大きく下がったからと言って将来予測および管理基準値計算に自己相関を考慮すべきかはわかりません。個人的には、自己相関はレジームの代替手段として考慮するものですし、自己相関とレジームを両方考慮した再生産関係を見たことがないので、自己相関は考慮しないでいいのではないかと思います

outer_R1 = calc.residAR(resR1, per_regime = TRUE)
outer_R1$regime_pars

# A tibble: 2 x 5
  regime      a      b    sd     rho
   <dbl>  <dbl>  <dbl> <dbl>   <dbl>
1      0 0.0330 50504. 0.266 -0.141 
2      1 0.0255 43166. 0.192 -0.0686

outer_R1$AICc

   AR(0)    AR(1) 
2.071602 6.841064 

outer_R0 = calc.residAR(resR1, per_regime = FALSE)
outer_R0$regime_pars

# A tibble: 2 x 5
  regime      a      b    sd    rho
   <dbl>  <dbl>  <dbl> <dbl>  <dbl>
1      0 0.0330 50504. 0.267 -0.112
2      1 0.0255 43166. 0.192 -0.112

outer_R0$AICc

   AR(0)    AR(1) 
2.071602 4.199656 

もしどうしてもレジームありで自己相関を内部推定したい場合には、fit.SRの引数のwを使用することで、レジームを完全に分けての推定はできます

w = as.numeric(!(SRdata$year %in% 1976:1987))
resR3 <- fit.SR(SRdata, SR = "HS", method = "L2", AR = 1, out.AR = FALSE, w = w)
resR3$pars

           a        b        sd       rho
1 0.02953365 47322.39 0.2577441 0.1068873

自己相関に関するプロットはplot.autocorで出力されますが、標準偏差が年代によって変わっているので標準化残差を使っています

plot.autocor(resR1)

• fit.SRregimeでもboot.SRを使ってfit.SRと同様に、①残差のパラメトリックブートストラップ、②残差のノンパラメトリックブートストラップ、③データのブートストラップを行うことができます
• 残差の標準偏差がレジームによって変わるので、標準化残差に直したうえで残差ブートストラップをやり、元の標準偏差の大きさに戻しています
• ここでは例として残差のノンパラメトリックブートストラップを示します。レジーム毎に各パラメータのヒストグラムが表示されます

boot.res4 = boot.SR(resR1, n = 100, method = "n")
plot.bootSR(boot.res4)

レジーム間で共通のパラメータは各レジームで同じ図となります

boot.res5 = boot.SR(resR2, n = 100, method = "n")
plot.bootSR(boot.res5)

• ジャックナイフも同様にレジームごとに結果が表示されます
• レジームシフトの位置を青点線で示しています
• レジーム間で共通のパラメータがある場合はレジーム間で同じ図が表示されます

jack2 = jackknife.SR(resR1, is.plot = TRUE)

jack3 = jackknife.SR(resR2, is.plot = TRUE)

• プロファイル尤度もfit.SRのときと同様の関数prof.likSRで計算・プロットできます
• レジームごとにa,bのパラメータを動かし、別のレジームで推定が必要なパラメータは尤度が最大になるように推定しています（そのため計算に時間がかかります）
• 図はレジームごとに出力されます

prof2 = prof.likSR(resR1)