######相対密度(CPUE、モニタリングにおける発見頻度など)&防除努力量→ ###### 自然増加率、平均的な相対密度の動態、捕獲率と努力量の関係推定 ###### (個体数は推定しない) ###尤度関数:自然増加率もデータから推定(estim.pd関数で内部的に使用) ##引数 #par:パラメータ #ra:観測された相対密度 #eff:防除努力量 #occasion:時間ステップラベル #site:空間ユニットラベル #noc:時間ステップ数 #nsite:空間ユニット数 # データはそれぞれ時系列順に並べたものを使う ##返り値:対数尤度 loglf.pd<-function(par,nobs,surveyeff,removeeff,occasion,site,noc,nsite){ gr<-exp(par[1]) #自然増加率 cap<-exp(par[2]) #単位努力量当たり捕獲率 obssize<-exp(par[3]) racoef<-exp(par[(1:nsite)+3]) #初年度の平均相対密度 ny<-length(nobs) #年数 mu.ra<-eff.mat<-nobs.mat<-matrix(NA,nsite,noc) #平均相対密度、相対密度格納用行列 removeeff.mat<-surveyeff.mat<-matrix(0,nsite,noc) #努力量格納用行列 removeeff.mat[cbind(site,occasion)]<-removeeff surveyeff.mat[cbind(site,occasion)]<-surveyeff nobs.mat[cbind(site,occasion)]<-nobs mu.ra[1,]<-racoef for(i in 1:nsite){ for(j in 2:noc){ mu.ra[i,j]<-mu.ra[i,j-1]*gr*exp(-cap*(removeeff.mat[i,j-1])) #パラメータに基づく平均相対密度計算 } } print(mu.ra) cat("\n") mu.nobs<-mu.ra*surveyeff.mat loglf.samp<-dnbinom(c(nobs.mat),mu=c(mu.nobs),size=obssize,log=T) #対数尤度の計算 loglf<-sum(loglf.samp[!is.na(c(nobs.mat))]) #対数尤度の合計 return(loglf) } ###尤度関数:自然増加率は決め打ち(estim.pd関数で内部的に使用) ##引数 #par:パラメータ #nobs:観測された個体数 #surveyeff:観測努力量 #removeeff:防除努力量 # データはそれぞれ時系列順に並べたものを使う #fixr:事前に与える自然増加率 #occasion:時間ステップラベル #site:空間ユニットラベル #noc:時間ステップ数 #nsite:空間ユニット数 ##返り値:対数尤度 loglf.pd.fixr<-function(par,nobs,surveyeff,removeeff,occasion,site,noc,nsite,fixr){ cap<-exp(par[1]) #単位努力量当たり捕獲率 obssize<-exp(par[2]) racoef<-exp(par[(1:nsite)+2]) #初年度の平均相対密度 ny<-length(nobs) #年数 mu.ra<-eff.mat<-nobs.mat<-matrix(NA,nsite,noc) #平均相対密度、相対密度格納用行列 removeeff.mat<-surveyeff.mat<-matrix(0,nsite,noc) #努力量格納用行列 removeeff.mat[cbind(site,occasion)]<-removeeff surveyeff.mat[cbind(site,occasion)]<-surveyeff nobs.mat[cbind(site,occasion)]<-nobs mu.ra[1,]<-racoef for(i in 1:nsite){ for(j in 2:noc){ mu.ra[i,j]<-mu.ra[i,j-1]*fixr*exp(-cap*removeeff.mat[i,j-1]) #パラメータに基づく平均相対密度計算 } } print(mu.ra) cat("\n") mu.nobs<-mu.ra*surveyeff.mat loglf.samp<-dnbinom(c(nobs.mat),mu=c(mu.nobs),size=obssize,log=T) #対数尤度の計算 loglf<-sum(loglf.samp[!is.na(c(nobs.mat))]) #対数尤度の合計 return(loglf) } ###パラメータ推定関数 ##引数 #data:入力データ(必須) #noc:時間ステップ数 #nsite:空間ユニット数 #fixr:自然増加率を決め打ちにする場合にのみそれを入力(オプション) #init:パラメータ推定の初期値を入力(オプション、なければ自動的に決定) #maxit:推定計算の反復回数(オプション、収束しないときは大きな値にする) #...:control.optimの引数(オプション、上級者向き) ##返り値 #リスト┬推定結果(関数optimの返り値) # ├データ #   ├時間ステップ数 # ├空間ユニット数 # └決め打ち自然増加率(データから推定した場合はNULL値) estim.pd<-function(data,noc,nsite,fixr=NULL,init=NULL,maxit=10000,method="BFGS",...){ occasion<-data[,1] #時間のラベル site<-data[,2] #空間ユニットのラベル nobs<-data[,3] #発見数 surveyeff<-data[,4] #調査努力量 removeeff<-data[,5] #防除努力量 if(is.null(init)){ #推定計算の初期値を決めていない場合 ra1<-(nobs/surveyeff)[occasion==1] init.ln.racoef<-ifelse(is.null(nobs),0,log(ra1)) #各サイトのln(平均初期相対密度)の初期値 init.ln.cap<-log(0.01/max(removeeff)) #ln(単位努力量当たり捕獲率)の初期値 init.ln.obssize<-log(1) #ln(観測誤差の尺度母数)の初期値 if(is.null(fixr)){ #自然増加率も推定する場合 init.ln.gr<-0 #自然増加率の初期値 init<-c(init.ln.gr,init.ln.cap,init.ln.obssize,init.ln.racoef) }else{ #自然増加率を事前に与える場合 init<-c(init.ln.cap,init.ln.obssize,init.ln.racoef) } } res<-vector("list",5); names(res)<-c("result","data","noc","nsite","fixr") res[[2]]<-data res[[3]]<-noc res[[4]]<-nsite #対数尤度の最大化 if(is.null(fixr)){ res[[1]]<-optim(init,loglf.pd,control=list(fnscale=-1,maxit=maxit,...),hessian=T,nobs=nobs,surveyeff=surveyeff,removeeff=removeeff,occasion=occasion,site=site,noc=noc,nsite=nsite,method=method) }else{ res[[1]]<-optim(init,loglf.pd.fixr,control=list(fnscale=-1,maxit=maxit),hessian=T,nobs=nobs,surveyeff=surveyeff,removeeff=removeeff,fixr=fixr,occasion=occasion,site=site,noc=noc,nsite=nsite,method=method) res[[5]]<-fixr } return(res) } ###要約統計量 ##引数:estim.pdの結果 ##返り値: #ln(平均初期相対密度u1)、ln(自然増加率r)、 #ln(単位努力量当たり捕獲率c)、ln(観測誤差の尺度母数theta)の最尤推定値、標準誤差、95%CI summary.pd<-function(res){ if(res$result$convergence!=0){ #もし推定が収束していない場合は中断 stop(message="パラメータが収束していないようです。設定を変えて計算しなおしてください。") } estim<-res$result$par #最尤推定値 se.pd<-sqrt(diag((-1)*solve(res$result$hessian))) #SE lci=estim-se.pd*1.96 uci=estim+se.pd*1.96 res.summ<-cbind(MLE=estim,SE=se.pd,CI2.5=lci,CI97.5=uci) #行名 len.estim<-length(estim) u1names<-paste("ln(u1",1:res$nsite,")",sep="") if(is.null(res$fixr)){ rownames(res.summ)<-c("ln(r)","ln(c)","ln(size)",u1names) cat(paste("自然増加率 = ",signif(exp(estim[1]),4),"; 単位努力量当たり捕獲率 = ",signif(exp(estim[2]),4),"\n",sep="")) }else{ rownames(res.summ)<-c("ln(c)","ln(size)",u1names) cat(paste("自然増加率 = ",signif(res$fixr,4),"; 単位努力量当たり捕獲率 = ",signif(exp(estim[1]),4),"\n",sep="")) } return(res.summ) } ###努力量と除去率の関係 ##引数:estim.pdの結果 ##返り値: #防除努力量に対する除去率の予測値(平均&95%CI)と作図 removalrate.pd<-function(res){ is.fixr<-!is.null(res$fixr) x<-seq(0,max(res$data[,5]),length=100) estim.c<-res$result$par[ifelse(is.fixr,1,2)] se.c<-sqrt(diag((-1)*solve(res$result$hessian)))[ifelse(is.fixr,1,2)] #SE lci.c=estim.c-se.c*1.96 uci.c=estim.c+se.c*1.96 p.remove<-1-exp(-exp(estim.c)*x) p.remove.lci<-1-exp(-exp(lci.c)*x) p.remove.uci<-1-exp(-exp(uci.c)*x) plot(x,p.remove,type="l",ylim=c(0,1),yaxs="i",xaxs="i",xlab="Removal Effort",ylab="Removal Rate") points(x,p.remove.lci,type="l",lty=2) points(x,p.remove.uci,type="l",lty=2) res<-data.frame(effort=x,removalrate=p.remove,ci0.25=p.remove.lci,ci97.5=p.remove.uci) return(res) } ###防除シミュレーション(初期個体数指定) #初期個体数、自然増加率、防除係数、毎年の防除努力量を指定して、モンテカルロシミュレーションによる将来予測 #個体群動態の不確実性としてポアソン分布を仮定 ##引数: #r:自然増加率 #c:防除係数 #effort:防除努力量のベクトル #nyrs:シミュレーション年数 #nind.init:初期個体群サイズ(ベクトルで与えれば、その回数シミュレーションを繰り返す。デフォルトでは初期個体数を50〜100の一様乱数で決めたシミュレーションを100回繰り返す) ##返り値: #反復×年数のシミュレーション個体数が出力された行列 pred.nind<-function(r,c,effort,nyrs,nind.init=round(runif(100,50,100))){ niter<-length(nind.init) res.mat<-matrix(NA,niter,nyrs) #シミュレーション結果を入れる箱 res.mat[,1]<-nind.init for(i in 2:nyrs){ mu.nind<-res.mat[,i-1]*r*exp(-c*effort[i-1]) #平均個体数の計算 res.mat[,i]<-rpois(niter,mu.nind) #ポアソン乱数による個体数の発生 } return(res.mat) } ###防除シミュレーション(相対個体数の動態) #自然増加率、防除係数、毎年の防除努力量のみを指定して、平均的な相対密度の変化を計算する #個体群動態の不確実性は考慮しない ##引数: #r:自然増加率 #c:防除係数 #effort:防除努力量のベクトル #nyrs:シミュレーション年数 ##返り値: #初期状態を1としたときの相対平均個体密度を計算 pred.rdens<-function(r,c,effort,nyrs){ res<-rep(NA,nyrs) #シミュレーション結果を入れる箱 res[1]<-1 for(i in 2:nyrs){ res[i]<-res[i-1]*r*exp(-c*effort[i-1]) #平均個体数の計算 } return(res) }