再生産関係チェックリスト
再生産関係の選択・診断を行うためのチェックリストです。ここで、将来予測・管理基準値計算チュートリアル (../docs/future-doc-abc.html) の例を使って説明します。
1. 再生産関係型の比較
Hockey-Stick型/Beverton-Holt型/Ricker型の再生産関係を比較します。
SRdata <- get.SRdata(res.pma) #res.pma: vpaの結果
resHS <- fit.SR(SRdata,SR="HS",method="L2",AR=0)
resBH <- fit.SR(SRdata,SR="BH",method="L2",AR=0)
resRI <- fit.SR(SRdata,SR="RI",method="L2",AR=0)
plot(SRdata$R ~ SRdata$SSB, cex=2, type = "b",xlab="SSB",ylab="R",
main="HS vs. BH vs. RI",ylim=c(0,max(SRdata$R)*1.3),xlim=c(0,max(SRdata$SSB)*1.1))
points(rev(SRdata$SSB)[1],rev(SRdata$R)[1],col=1,type="p",lwd=3,pch=16,cex=2)
points(resHS$pred$SSB,resHS$pred$R,col=2,type="l",lwd=3)
points(resBH$pred$SSB,resBH$pred$R,col=3,type="l",lwd=3,lty=2)
points(resRI$pred$SSB,resRI$pred$R,col=4,type="l",lwd=3,lty=3)
legend("topleft",
legend=c(sprintf("HS %5.2f",resHS$AICc),sprintf("BH %5.2f",resBH$AICc),sprintf("RI %5.2f",resRI$AICc)),
lty=1:3,col=2:4,lwd=2,title="AICc",ncol=3)
resSR <- resHS #HSを選択選択した再生産関係：Hockey-Stick
理由（特にBHやRIの場合は詳しく）： BHとRIでは再生産関係がほぼ直線になるため、HSを選択した。AICcもHSがやや小さい。
2. 最小絶対値法や残差の自己相関の検討
サンプル数が少ない場合や残差が正規分布に従っていない場合（3節参照）、外れ値に対して頑健な推定方法である最小絶対値法（中央値推定）が有効なオプションとして考えられます。 また、加入の残差が環境影響などにより時間的なトレンドをもつ場合には（4節参照）、残差の自己相関を考慮する方法が考えられます。 ここでは、これらの手法により再生産関係がどの程度変わるのかをチェックします。
resAR1 <- fit.SR(SRdata,SR="HS",method="L2",AR=1)
resL1 <- fit.SR(SRdata,SR="HS",method="L1",AR=0)
plot(SRdata$R ~ SRdata$SSB, cex=2, type = "b",xlab="SSB",ylab="R",
main="Effects of autocorrelation and L1",ylim=c(0,max(SRdata$R)*1.3),xlim=c(0,max(SRdata$SSB)*1.1))
points(rev(SRdata$SSB)[1],rev(SRdata$R)[1],col=1,type="p",lwd=3,pch=16,cex=2)
points(resSR$pred$SSB,resSR$pred$R,col=2,type="l",lwd=3)
points(resAR1$pred$SSB,resAR1$pred$R,col=3,type="l",lwd=3,lty=2)
points(resL1$pred$SSB,resL1$pred$R,col=4,type="l",lwd=3,lty=3)
legend("topleft",
legend=c(sprintf("L2&AR0 %5.2f",resSR$AICc),sprintf("L2&AR1 %5.2f",resAR1$AICc),sprintf("L1&AR0 %5.2f",resL1$AICc)),
lty=1:3,col=2:4,lwd=2,title="AICc",ncol=3)
resSR <- resL1 #L1 normを採用選択した再生産関係：L1&AR0
理由： 自己相関を推定した場合、自己相関係数は低く(rho=0.05)、AICcは高くなった。残差が正規分布に従っていなかったため（下記参照）、再生産関係は大きく変わらないものの、最小絶対値法による中央値推定を採用した。
3. 正規性のチェック
再生産関係から予測される加入量と観測値（資源評価値）の残差が正規分布に従っているかをチェックします。 Shapiro-Wilk検定とKolmogorov-Smirnov 検定を行い、「残差が正規分布に従っている」という帰無仮説を検定します。 また、QQ plotを描き、理論予測値 (y=x) から大きく逸脱していないかをチェックします。 これらの結果、正規性が疑われる場合には最小絶対値法を検討することが望ましいです。
check1 <- shapiro.test(resSR$resid)
check2 <- ks.test(resSR$resid,y="pnorm")
par(mfrow=c(1,2),mar=c(4,4,2,2))
hist(resSR$resid,xlab="Residuals",main="Normality test",freq=FALSE)
X <- seq(min(resSR$resid)*1.3,max(resSR$resid)*1.3,length=200)
points(X,dnorm(X,0,resSR$pars$sd),col=2,lwd=3,type="l")
mtext(text=" P value",adj=1,line=-1,lwd=2,font=2)
mtext(text=sprintf(" SW: %1.3f",check1$p.value),adj=1,line=-2)
mtext(text=sprintf(" KS: %1.3f",check2$p.value),adj=1,line=-3)
qqnorm(resSR$resid2,cex=2)
qqline(resSR$resid2,lwd=3)
診断結果：Kolmogorov-Smirnov 検定では有意となり、正規分布に従ってないことが示唆された。QQ plotもやや直線とはずれていた。
4. 残差のトレンドと自己相関係数
残差の時間的なトレンドをチェックします。 トレンドが見られたり、自己相関係数が有意である場合には、残差の自己相関を考慮した再生産関係を検討することが望ましいです。
par(mfrow=c(1,2),mar=c(4,4,2,2))
plot(SRdata$year, resSR$resid2,pch=16,main="",xlab="Year",ylab="Residual")
abline(0,0,lty=2)
par(new=T)
scatter.smooth(SRdata$year, resSR$resid2, lpars=list(col="red", lwd=2),ann=F,axes=FALSE)
ac.res <- acf(resSR$resid2,plot=FALSE)
plot(ac.res,main="",lwd=3)
診断結果：残差に時間的なトレンドは見られるが、自己相関係数は有意でなかった。
5. 残差ブートストラップ
パラメータ推定の信頼性をチェックするために、残差ブートストラップを行います。信頼区間が広い場合や、ブートストラップの中央値と点推定値の乖離が大きい場合には、パラメータ推定の信頼性が低いことになります。
boot.res <- boot.SR(resSR)
par(mfrow=c(2,2),mar=c(4,4,2,2))
hist(sapply(1:length(boot.res), function(i) boot.res[[i]]$pars$a),xlab="",ylab="",main="a")
abline(v=resSR$pars$a,col=2,lwd=3)
abline(v=median(sapply(1:length(boot.res), function(i) boot.res[[i]]$pars$a)),col=3,lwd=3,lty=2)
arrows(quantile(sapply(1:length(boot.res), function(i) boot.res[[i]]$pars$a),0.1),0,
quantile(sapply(1:length(boot.res), function(i) boot.res[[i]]$pars$a),0.9),0,
col=4,lwd=3,code=3)
legend("topright",
legend=c("Estimate","Median","CI(0.8)"),lty=1:2,col=2:4,lwd=2,ncol=1,cex=1)
hist(sapply(1:length(boot.res), function(i) boot.res[[i]]$pars$b),xlab="",ylab="",main="b")
abline(v=resSR$pars$b,col=2,lwd=3)
abline(v=median(sapply(1:length(boot.res), function(i) boot.res[[i]]$pars$b)),col=3,lwd=3,lty=2)
arrows(quantile(sapply(1:length(boot.res), function(i) boot.res[[i]]$pars$b),0.1),0,
quantile(sapply(1:length(boot.res), function(i) boot.res[[i]]$pars$b),0.9),0,
col=4,lwd=3,code=3)
hist(sapply(1:length(boot.res), function(i) boot.res[[i]]$pars$sd),xlab="",ylab="",main="sd")
abline(v=resSR$pars$sd,col=2,lwd=3)
abline(v=median(sapply(1:length(boot.res), function(i) boot.res[[i]]$pars$sd)),col=3,lwd=3,lty=2)
arrows(quantile(sapply(1:length(boot.res), function(i) boot.res[[i]]$pars$sd),0.1),0,
quantile(sapply(1:length(boot.res), function(i) boot.res[[i]]$pars$sd),0.9),0,
col=4,lwd=3,code=3)
if (resSR$input$AR==1) {
hist(sapply(1:length(boot.res), function(i) boot.res[[i]]$pars$rho),xlab="",ylab="",main="rho")
abline(v=resSR$pars$rho,col=2,lwd=3)
abline(v=median(sapply(1:length(boot.res), function(i) boot.res[[i]]$pars$rho)),col=3,lwd=3,lty=2)
arrows(quantile(sapply(1:length(boot.res), function(i) boot.res[[i]]$pars$rho),0.1),0,
quantile(sapply(1:length(boot.res), function(i) boot.res[[i]]$pars$rho),0.9),0,
col=4,lwd=3,code=3)
}
par(mfrow=c(1,1))
plot(SRdata$R ~ SRdata$SSB, cex=2, type = "b",xlab="SSB",ylab="R",
main="Residual bootstrap",ylim=c(0,max(SRdata$R)*1.3))
points(rev(SRdata$SSB)[1],rev(SRdata$R)[1],col=1,type="p",lwd=3,pch=16,cex=2)
for (i in 1:length(boot.res)) {
points(boot.res[[i]]$pred$SSB,boot.res[[i]]$pred$R,type="l",lwd=2,col=rgb(0,0,1,alpha=0.1))
}
points(resSR$pred$SSB,resSR$pred$R,col=2,type="l",lwd=3)
診断結果：パラメータaとbは、ブートストラップ推定値の中央値と点推定値がほぼ一致した。sigmaはややずれていた。 再生産関係は比較的ロバストであった。
6. ジャックナイフ推定
パラメータ推定の頑健性を調べるために、一点ずつ除いてデータをジャックナイフ解析を行います。これにより、どの年のデータの影響が大きいかが明らかになります。
jack.res <- lapply(1:length(SRdata$year), function(i){
jack <- resSR
jack$input$w[i] <- 0
do.call(fit.SR,jack$input)
})
par(mfrow=c(2,2),mar=c(4,4,2,2))
plot(SRdata$year,sapply(1:length(SRdata$year), function(i) jack.res[[i]]$pars$a),type="b",
xlab="Removed year",ylab="",main="a",pch=19)
abline(resSR$pars$a,0,lwd=3,col=2)
plot(SRdata$year,sapply(1:length(SRdata$year), function(i) jack.res[[i]]$pars$b),type="b",
xlab="Removed year",ylab="",main="b",pch=19)
abline(resSR$pars$b,0,lwd=3,col=2)
plot(SRdata$year,sapply(1:length(SRdata$year), function(i) jack.res[[i]]$pars$sd),type="b",
xlab="Removed year",ylab="",main="sd",pch=19)
abline(resSR$pars$sd,0,lwd=3,col=2)
if (resSR$input$AR==1){
plot(SRdata$year,sapply(1:length(SRdata$year), function(i) jack.res[[i]]$pars$rho),type="b",
xlab="Removed year",ylab="",main="rho",pch=19)
abline(resSR$pars$rho,0,lwd=3,col=2)
}
par(mfrow=c(1,1)) 
plot(SRdata$R ~ SRdata$SSB, cex=2, type = "b",xlab="SSB",ylab="R",
main="Jackknife estimate",ylim=c(0,max(SRdata$R)*1.3))
points(rev(SRdata$SSB)[1],rev(SRdata$R)[1],col=1,type="p",lwd=3,pch=16,cex=2)
for (i in 1:length(jack.res)) {
points(jack.res[[i]]$pred$SSB,jack.res[[i]]$pred$R,type="l",lwd=3,col=rgb(0,0,1,alpha=0.1))
}
points(resSR$pred$SSB,resSR$pred$R,col=2,type="l",lwd=3)
診断結果：1992年・2000年のデータを除くとbの推定値が大きくなる。 そのほかのパラメータは各データの除去に対して比較的頑健である。
7. プロファイル尤度
パラメータa, bの値を変えたときに尤度がどの程度変化するかを解析します。ここでは、sdやrhoは最尤推定値に固定しています。この結果とブートストラップ信頼区間およびジャックナイフ解析の結果も同時に図示します。
ngrid <- 100
a.grid <- seq(resSR$pars$a*0.5,resSR$pars$a*1.5,length=ngrid)
b.grid <- seq(min(SRdata$SSB),max(SRdata$SSB),length=ngrid)
ba.grids <- expand.grid(b.grid,a.grid)
prof.lik.res <- sapply(1:nrow(ba.grids),function(i) prof.lik(resSR,a=as.numeric(ba.grids[i,2]),b=as.numeric(ba.grids[i,1])))
image(b.grid,a.grid,matrix(prof.lik.res,nrow=ngrid),ann=F,col=cm.colors(12),
ylim=c(resSR$pars$a*0.5,resSR$pars$a*1.5),xlim=c(min(SRdata$SSB),max(SRdata$SSB)))
par(new=T, xaxs="i",yaxs="i")
contour(b.grid,a.grid,matrix(prof.lik.res,nrow=ngrid),
ylim=c(resSR$pars$a*0.5,resSR$pars$a*1.5),xlim=c(min(SRdata$SSB),max(SRdata$SSB)),
xlab="b",ylab="a",main="Profile likelihood")
for(i in 1:length(jack.res)) points(jack.res[[i]]$pars$b,jack.res[[i]]$pars$a,lwd=1,col=1)
lines(y=as.numeric(quantile(sapply(1:length(boot.res),function(i)boot.res[[i]]$pars$a),c(0.1,0.9))),
x=rep(resSR$pars$b,2),col=4,lwd=2)
lines(x=as.numeric(quantile(sapply(1:length(boot.res),function(i)boot.res[[i]]$pars$b),c(0.1,0.9))),
y=rep(resSR$pars$a,2),col=4,lwd=2)
legend("bottomleft",c("Bootstrap CI(0.8)","Jackknife"),lty=1:0,pch=c("","○"),col=c(4,1),lwd=2:1)
診断結果：パラメータbの値が大きくなっても尤度の変化は比較的小さい