Section 6.4.1 Wildebeest
Logistic growth model, matrix model with dependencies of model parameters upon a number of rainfall as well as independent estimates of survival. Note the model with constant adult survival, constant lambda and calf survival dependent upon previous year’s rainfall has been commented out to save computing time.
Wildebeest data
yr=c(61,63,65,67,71,72,77,78,82,84,86,89) Nhat=c(263,357,439,483,693,773,1444,1249,1209,1338,1146,1686)/1000 se=c(NA,NA,NA,NA,28.8,76.7,200,355,272,138,133,176)/1000 cvhat=se/Nhat names(Nhat)=names(se)=names(cvhat)=yr # randomly select from exisitng CVs to create CVs fsro first four years set.seed=10293845 nas=which(is.na(cvhat)) cvhat[nas]=sample(na.omit(cvhat),length(nas),replace=TRUE) sehat=cvhat*Nhat allys=60:89 ymin=min(allys) y=yr-ymin ys=allys-ymin # Rainfall: rains=c(100,38,100,104,167,167,165,79,91,77,134,192,235,159,211,257,204,300,187,84,99,163,97,228,208,83, 44,112,191,202)/100 names(rains)=allys rain=c(38,104,167,79,192,235,300,187,97,208,44,202)/100 # just in years with N estimates names(rain)=yr # Catch: Ca=c(rep(0,length(allys))) Ca[allys>=77]=40/1000 # Hilborn & Mangel's best estimate of poaching: 40,000 poached all years from 1977 names(Ca)=allys # Calf survival estimates: phi1hat=c(.5,.25,.3,.26,.36,.32,.35,.36) # calf survival estimates yrc=c(64,65,66,67,68,70,71,72) # calf survival years yc=yrc-ymin # shift origin to same as for other years # Adult survival estimates: phi2hat=c(1-0.017,1-0.014,1-0.008,1-0.005,1-0.027,1-0.021)^4 # adult survival estimates yra=c(68,69,71,72,82,83) # adult survival years ya=yra-ymin # shift origin to same as for other years lower=Nhat-1.96*sehat upper=Nhat+1.96*sehat
Supplemental functions for wildebeest analysis
# ================================ Logistic Model Fitting ============================================ Ny=function(N0,r,k,y,ys,Ca){ nys=length(ys) Ns=rep(N0,nys) names(Ns)=names(Ca) for(i in 2:nys){ Ns[i]=Ns[i-1]+r*Ns[i-1]*(1-Ns[i-1]/k) - Ca[i-1] } N=Ns[which(is.element(ys,y))] return(N) } Ny.rain=function(N0,r,k,y,ys,Ca){ nys=length(ys) Ns=rep(N0,nys) names(Ns)=names(Ca) for(i in 2:nys){ Ns[i]=Ns[i-1]+r[i-1]*Ns[i-1]*(1-Ns[i-1]/k) - Ca[i-1] } N=Ns[which(is.element(ys,y))] return(N) } Ny.rain.k=function(N0,r,k,y,ys,Ca){ nys=length(ys) Ns=rep(N0,nys) names(Ns)=names(Ca) for(i in 2:nys){ Ns[i]=Ns[i-1]+r[i-1]*Ns[i-1]*(1-Ns[i-1]/k[i-1]) - Ca[i-1] } N=Ns[which(is.element(ys,y))] return(N) } wild.negllik=function(pars,y,Nhat,sehat,ys,Ca){ nll=sum(log(sehat)+(Nhat-Ny(exp(pars[1]),exp(pars[2]),exp(pars[3]),y,ys,Ca))^2/(2*sehat^2)) return(nll) } wild.rain.negllik=function(pars,y,Nhat,sehat,rains,ys,Ca){ r=exp(pars[2]+pars[3]*rains) nll=sum(log(sehat)+(Nhat-Ny.rain(exp(pars[1]),r,exp(pars[4]),y,ys,Ca))^2/(2*sehat^2)) return(nll) } wild.rain.k.negllik=function(pars,y,Nhat,sehat,rains,ys,Ca){ r=exp(pars[2]+pars[3]*rains) k=exp(pars[4]+pars[5]*rains) nll=sum(log(sehat)+(Nhat-Ny.rain.k(exp(pars[1]),r,k,y,ys,Ca))^2/(2*sehat^2)) return(nll) } wild.fit=function(N0,r,k,y,Nhat,sehat,ys,Ca,control=list(trace=0),hessian=FALSE) { pars=log(c(N0,r,k)) fit=optim(par=pars,fn=wild.negllik,y=y,Nhat=Nhat,sehat=sehat,ys=ys,Ca=Ca,control=control,hessian=hessian) N0=exp(fit$par[1]) r=exp(fit$par[2]) k=exp(fit$par[3]) est=list(N0=N0,r=r,k=k,N=Ny(N0,r,k,y,ys,Ca)) est$AIC=fit$value*2+2*length(fit$par) est$fit=fit est$type="logistic" return(est) } wild.rain.fit=function(N0,a,b,k,y,Nhat,sehat,rains,ys,Ca,control=list(trace=0),hessian=FALSE) { pars=c(log(N0),a,b,log(k)) fit=optim(par=pars,fn=wild.rain.negllik,y=y,Nhat=Nhat,sehat=sehat,rains=rains,ys=ys,Ca=Ca, control=control,hessian=hessian) N0=exp(fit$par[1]) a=fit$par[2] b=fit$par[3] k=exp(fit$par[4]) r=exp(a+b*rains) est=list(N0=N0,a=a,b=b,k=k,r=r,N=Ny.rain(N0,r,k,y,ys,Ca)) est$AIC=fit$value*2+2*length(fit$par) est$fit=fit est$type="logistic.rain" return(est) } wild.rain.k.fit=function(N0,a,b,k1,k2,y,Nhat,sehat,rains,ys,Ca,control=list(trace=0),hessian=FALSE) { pars=c(log(N0),a,b,k1,k2) fit=optim(par=pars,fn=wild.rain.k.negllik,y=y,Nhat=Nhat,sehat=sehat,rains=rains,ys=ys,Ca=Ca, control=control,hessian=hessian) N0=exp(fit$par[1]) a=fit$par[2] b=fit$par[3] k1=fit$par[4] k2=fit$par[5] k=exp(k1+k2*rains) r=exp(a+b*rains) est=list(N0=N0,a=a,b=b,k1=k1,k2=k2,k=k,r=r,N=Ny.rain.k(N0,r,k,y,ys,Ca)) est$AIC=fit$value*2+2*length(fit$par) est$fit=fit est$type="logistic.rain.k" return(est) } wild.rain.negllik=function(pars,y,Nhat,sehat,rains,ys,Ca){ r=exp(pars[2]+pars[3]*rains) nll=sum(log(sehat)+(Nhat-Ny.rain(exp(pars[1]),r,exp(pars[4]),y,ys,Ca))^2/(2*sehat^2)) return(nll) } wild.rain.k.negllik=function(pars,y,Nhat,sehat,rains,ys,Ca){ r=exp(pars[2]+pars[3]*rains) k=exp(pars[4]+pars[5]*rains) nll=sum(log(sehat)+(Nhat-Ny.rain.k(exp(pars[1]),r,k,y,ys,Ca))^2/(2*sehat^2)) return(nll) } wild.rain.fit=function(N0,a,b,k,y,Nhat,sehat,rains,ys,Ca,control=list(trace=0),hessian=FALSE) { pars=c(log(N0),a,b,log(k)) fit=optim(par=pars,fn=wild.rain.negllik,y=y,Nhat=Nhat,sehat=sehat,rains=rains,ys=ys,Ca=Ca, control=control,hessian=hessian) N0=exp(fit$par[1]) a=fit$par[2] b=fit$par[3] k=exp(fit$par[4]) r=exp(a+b*rains) est=list(N0=N0,a=a,b=b,k=k,r=r,N=Ny.rain(N0,r,k,y,ys,Ca)) est$AIC=fit$value*2+2*length(fit$par) est$fit=fit est$type="logistic.rain" return(est) } wild.rain.k.fit=function(N0,a,b,k1,k2,y,Nhat,sehat,rains,ys,Ca,control=list(trace=0),hessian=FALSE) { pars=c(log(N0),a,b,k1,k2) fit=optim(par=pars,fn=wild.rain.k.negllik,y=y,Nhat=Nhat,sehat=sehat,rains=rains,ys=ys,Ca=Ca, control=control,hessian=hessian) N0=exp(fit$par[1]) a=fit$par[2] b=fit$par[3] k1=fit$par[4] k2=fit$par[5] k=exp(k1+k2*rains) r=exp(a+b*rains) N=Ny.rain.k(N0,r,k,y,ys,Ca) all.N=Ny.rain.k(N0,r,k,ys,ys,Ca) est=list(N0=N0,a=a,b=b,k1=k1,k2=k2,k=k,r=r,N=N,ys=ys,all.N=all.N) est$AIC=fit$value*2+2*length(fit$par) est$fit=fit est$type="logistic.rain.k" return(est) } # ================================ Matrix Model Fitting ============================================ wildbst.MM.fit=function(parobj,ys,y,Nhat,sehat,Ca,rains, yc,phi1hat,ya,phi2hat,addata=FALSE, control=list(trace=0),hessian=FALSE) { tfmobj=tfm.MM(parobj) tpars=tfmobj$tpars # pars=c(log(N0),logit(b1[1]),log(b1[2]),logit(b2[1]),log(b2[2])) fit=optim(par=tpars,fn=wildbst.MM.negllik,ys=ys,y=y,Nhat=Nhat,sehat=sehat,Ca=Ca,rains=rains,tfmobj=tfmobj, yc=yc,phi1hat=phi1hat,ya=ya,phi2hat=phi2hat,addata=addata, control=control,hessian=hessian) parobj=untfm.MM(tfmobj) startvals=parobj$pars tfmobj$tpars=fit$par est.pars=untfm.MM(tfmobj)$pars parobj=untfm.MM(tfmobj) N=Ny.MM(parobj,y,ys,Ca,rains) # abundances at years with Nhats # phis=calc.phis(parobj,yc,ys,Ca,rains) all.N=Ny.MM(parobj,ys,ys,Ca,rains) # abundances at all years all.phis=calc.phis(parobj,ys,ys,Ca,rains) est=list(pars=est.pars,startvals=startvals,fixedpars=parobj$fixed,start.tfmvals=tpars,N=N,phis=all.phis,all.N=all.N,ys=ys) est$AIC=fit$value*2+2*length(fit$par) est$fit=fit est$type="MMadd" return(est) } wildbst.MM.negllik=function(pars,ys,y,Nhat,sehat,Ca,rains,tfmobj,yc,phi1hat,ya,phi2hat,addata=FALSE){ tfmobj$tpars=pars parobj=untfm.MM(tfmobj) Ny=Ny.MM(parobj,y,ys,Ca,rains) # abundances at years with Nhats nll.Nhat=sum(log(sehat)+(Nhat-Ny)^2/(2*sehat^2)) if(addata){ # additional estimates of calf and adult survival (with CVs=0.3) sigma2.phi1hat=log(0.3^2+1) # assuming CV(phi1hat)=0.3 since CV(x)=sqrt(exp(sigma^2)-1) sigma2.phi2hat=log(0.3^2+1) # assuming CV(phi2hat)=0.3, since CV(x)=sqrt(exp(sigma^2)-1) phi1s=calc.phis(parobj,yc,ys,Ca,rains)$phi1 # phi1 in years with phi1 estimates phi2s=calc.phis(parobj,ya,ys,Ca,rains)$phi2 # phi2 in years with phi2 estimates nll.phi1=sum((log(phi1hat)-phi1s)^2/(2*sigma2.phi1hat)) nll.phi2=sum((log(phi2hat)-phi2s)^2/(2*sigma2.phi2hat)) nll=nll.Nhat+nll.phi1+nll.phi2 } else { nll=nll.Nhat } return(nll) } Ny.MM=function(parobj,y,ys,Ca,rains){ fpar=unpack.pars(parobj) nys=length(ys) Ns=rep(fpar["N0"],nys) names(Ns)=ys for(i in 2:nys){ phi1=phi(c(fpar["phi1"],fpar["rain1"]),rains[i-1],Ns[i-1]) phi2=phi(c(fpar["phi2"],fpar["rain2"]),rains[i-1],Ns[i-1]) Ns[i]=phi1*lambda*Ns[i-1]+phi2*Ns[i-1]-Ca[i-1] } N=Ns[which(is.element(ys,y))] return(N) } phi=function(b,Rt1,Nt) { return(1.25*b[1]*Rt1/(1.25*Rt1+b[2]*Nt/1)) # A=1 in millions } calc.phis=function(parobj,y,ys,Ca,rains){ fpar=unpack.pars(parobj) nys=length(ys) Ns=rep(fpar["N0"],nys) phi1=phi2=rep(0,nys) names(Ns)=names(phi1)=names(phi2)=ys for(i in 2:nys){ phi1[i]=phi(c(fpar["phi1"],fpar["rain1"]),rains[i-1],Ns[i-1]) phi2[i]=phi(c(fpar["phi2"],fpar["rain2"]),rains[i-1],Ns[i-1]) Ns[i]=phi1[i]*lambda*Ns[i-1]+phi2[i]*Ns[i-1]-Ca[i-1] } phi1=phi1[which(is.element(ys,y))] phi2=phi2[which(is.element(ys,y))] return(list(phi1=phi1,phi2=phi2)) } #tfmobj=tfm.MM(N0,phi1,rain1,phi2,rain2,lambda,fixed) #untfm.MM(tfmobj) #' @description Transforms wildebeest matrix model parameters onto appropriate scales for numerical optimization. tfm.MM=function(parobj){ pars=parobj$pars fixed=names(parobj$fixed) if(length(fixed)==0) fixed="" if(!is.element("rain1",fixed)) { # estimate rain1 if(!is.element("rain2",fixed)) { # estimate rain2 if(!is.element("lambda",fixed)) { # estimate all tpars=c(log(pars["N0"]),logit(pars["phi1"]),log(pars["rain1"]),logit(pars["phi2"]),log(pars["rain2"]),logit(pars["lambda"])) names(tpars)=c("log(N0)","logit(phi1)","log(rain1)","logit(phi2)","log(rain2)","logit(lambda)") } else { # estimate rain1, rain2 tpars=c(log(pars["N0"]),logit(pars["phi1"]),log(pars["rain1"]),logit(pars["phi2"]),log(pars["rain2"])) names(tpars)=c("log(N0)","logit(phi1)","log(rain1)","logit(phi2)","log(rain2)") } } else { # don't estimate rain2: if(!is.element("lambda",fixed)) { # estimate rain1, lambda tpars=c(log(pars["N0"]),logit(pars["phi1"]),log(pars["rain1"]),logit(pars["phi2"]),logit(pars["lambda"])) names(tpars)=c("log(N0)","logit(phi1)","log(rain1)","logit(phi2)","logit(lambda)") } else { # estimate rain1 tpars=c(log(pars["N0"]),logit(pars["phi1"]),log(pars["rain1"]),logit(pars["phi2"])) names(tpars)=c("log(N0)","logit(phi1)","log(rain1)","logit(phi2)") } } } else { # don't estimate rain1 if(!is.element("rain2",fixed)) { # estimate rain2 if(!is.element("lambda",fixed)) { # estimate rain2, lambda tpars=c(log(pars["N0"]),logit(pars["phi1"]),logit(pars["phi2"]),log(pars["rain2"]),logit(pars["lambda"])) names(tpars)=c("log(N0)","logit(phi1)","logit(phi2)","log(rain2)","logit(lambda)") } else { # estimate rain2 tpars=c(log(pars["N0"]),logit(pars["phi1"]),logit(pars["phi2"]),log(pars["rain2"])) names(tpars)=c("log(N0)","logit(phi1)","logit(phi2)","log(rain2)") } } else { # don't estimate rain2: if(!is.element("lambda",fixed)) { # estimate lambda tpars=c(log(pars["N0"]),logit(pars["phi1"]),logit(pars["phi2"]),logit(pars["lambda"])) names(tpars)=c("log(N0)","logit(phi1)","logit(phi2)","logit(lambda)") } else { # estimate none tpars=c(log(pars["N0"]),logit(pars["phi1"]),logit(pars["phi2"])) names(tpars)=c("log(N0)","logit(phi1)","logit(phi2)") } } } tpars[tpars==-Inf]=-99999999999999999 # replace negativeinfinity with large negative number return(list(tpars=tpars,fixed=parobj$fixed)) } get.tfm.names=function(tfmobj){ np=length(tfmobj$tpars) nam=names(tfmobj$tpars) tname=rep("",np) for(i in 1:np){ tname[i]=strsplit(strsplit(nam[i],"(",fixed=TRUE)[[1]][2],")",fixed=TRUE)[[1]] } return(tname) } untfm.MM=function(tfmobj){ fixedval=tfmobj$fixed fixed=names(fixedval) tpars=tfmobj$tpars # set indices pnames=get.tfm.names(tfmobj) N0=which(pnames=="N0") phi1=which(pnames=="phi1") rain1=which(pnames=="rain1") phi2=which(pnames=="phi2") rain2=which(pnames=="rain2") lambda=which(pnames=="lambda") if(!is.element("rain1",fixed)) { # estimate rain1 if(!is.element("rain2",fixed)) { # estimate rain2 if(!is.element("lambda",fixed)) { # estimate all pars=c(exp(tpars[N0]),inv.logit(tpars[phi1]),exp(tpars[rain1]),inv.logit(tpars[phi2]),exp(tpars[rain2]),inv.logit(tpars[lambda])) names(pars)=c("N0","phi1","rain1","phi2","rain2","lambda") } else { # estimate rain1, rain2 pars=c(exp(tpars[N0]),inv.logit(tpars[phi1]),exp(tpars[rain1]),inv.logit(tpars[phi2]),exp(tpars[rain2])) names(pars)=c("N0","phi1","rain1","phi2","rain2") } } else { # don't estimate rain2: if(!is.element("lambda",fixed)) { # estimate rain1, lambda pars=c(exp(tpars[N0]),inv.logit(tpars[phi1]),exp(tpars[rain1]),inv.logit(tpars[phi2]),inv.logit(tpars[lambda])) names(pars)=c("N0","phi1","rain1","phi2","lambda") } else { # estimate rain1 pars=c(exp(tpars[N0]),inv.logit(tpars[phi1]),exp(tpars[rain1]),inv.logit(tpars[phi2])) names(pars)=c("N0","phi1","rain1","phi2") } } } else { # don't estimate rain1 if(!is.element("rain2",fixed)) { # estimate rain2 if(!is.element("lambda",fixed)) { # estimate rain2, lambda pars=c(exp(tpars[N0]),inv.logit(tpars[phi1]),inv.logit(tpars[phi2]),exp(tpars[rain2]),inv.logit(tpars[lambda])) names(pars)=c("N0","phi1","phi2","rain2","lambda") } else { # estimate rain2 pars=c(exp(tpars[N0]),inv.logit(tpars[phi1]),inv.logit(tpars[phi2]),exp(tpars[rain2])) names(pars)=c("N0","phi1","phi2","rain2") } } else { # don't estimate rain2: if(!is.element("lambda",fixed)) { # estimate lambda pars=c(exp(tpars[N0]),inv.logit(tpars[phi1]),inv.logit(tpars[phi2]),inv.logit(tpars[lambda])) names(pars)=c("N0","phi1","phi2","lambda") } else { # estimate none pars=c(exp(tpars[N0]),inv.logit(tpars[phi1]),inv.logit(tpars[phi2])) names(pars)=c("N0","phi1","phi2") } } } return(list(pars=pars,fixed=fixedval)) } unpack.pars=function(parobj){ if(parobj$fixed[1]=="") { pars=parobj$pars } else { pars=c(parobj$pars,parobj$fixed) } pars=pars[order(names(pars))] return(pars) } pack.pars=function(N0,phi1,rain1,phi2,rain2,lambda,fixed="") { pars=c(N0,phi1,rain1,phi2,rain2,lambda) names(pars)=c("N0","phi1","rain1","phi2","rain2","lambda") fixeds=which(is.element(names(pars),fixed)) if(length(fixeds)==0) { return(list(pars=pars,fixed="")) } else { fixedpars=pars[fixeds] pars=pars[-fixeds] if(length(fixedpars)>0) fixedpars=fixedpars[order(names(fixedpars))] if(length(pars)>0) pars=pars[order(names(pars))] return(list(pars=pars,fixed=fixedpars)) } } wildbst.MM.bs.trajectory=function(est,ys,y,Nhat,sehat,Ca,rains,B,seed){ fixedpars=est$fixedpars mu=est$fit$par np=length(mu) sigma=solve(est$fit$hessian) nys=length(ys) N=matrix(rep(-1,B*nys),nrow=B) tparmat=parmat=matrix(rep(0,B*np),nrow=B) colnames(parmat)=names(est$pars) colnames(tparmat)=names(est$fit$par) if(!is.null(seed)) set.seed(seed) # calculate trajectories: b=bad.bs=0 while(b < B){ tfmpars=rmvnorm(1,mu,sigma) # B samples from multivariate normal parameters # transform parameters: names(tfmpars)=names(est$fit$par) tfmobj=list(tpars=tfmpars,fixed=fixedpars) parobj=untfm.MM(tfmobj) try.N=Ny.MM(parobj,ys,ys,Ca,rains) # abundances at all years if(try.N[1]<=try.N[length(try.N)]) { # remove stupid trajectories b=b+1 N[b,]=try.N tparmat[b,]=tfmobj$tpars parmat[b,]=parobj$pars } else { bad.bs=bad.bs+1 } } return(list(N=N,pars=parmat,tfmpars=tparmat,bad.bs=bad.bs)) }
Wildebeest analysis
library(boot) # needed for logit() and inv.logit() library(mvtnorm) # needed for parameteric bootstrap control=list(trace=0,maxit=3000) # controls optim() #====================== Logistic Model Estimation ======================== # estimate basic logistic model r=0.13 k=1.7 N0=Nhat[1] N=Ny(N0,r,k,y,ys,Ca) est.logistic=wild.fit(N0,r,k,y,Nhat,sehat,ys,Ca,control=control) est.logistic$AIC plot(y+60,Nhat,ylim=range(lower,upper),xlab="Year",ylab=expression(paste(hat(italic(N))," (millions)"),cex=0.5)) segments(x0=y+60,y0=lower,y1=upper) lines(y+60,est.logistic$N,lty=1) title("Wildebeest logistic",cex.main=0.9) # estimate logistic model with rain affecting r a=log(0.1) b=0.2 r=exp(a+b*rains) k=1.7 N0=Nhat[1] N=Ny.rain(N0,r,k,y,ys,Ca) est.logistic.rain=wild.rain.fit(N0,a,b,k,y,Nhat,sehat,rains,ys,Ca,control=control) est.logistic.rain$AIC plot(y+60,Nhat,ylim=range(lower,upper),xlab="Year",ylab=expression(paste(hat(italic(N))," (millions)"),cex=0.5)) segments(x0=y+60,y0=lower,y1=upper) lines(y+60,est.logistic$N,lty=1) lines(y+60,est.logistic.rain$N,lty=1,col="red") title("Wildebeest logistic (black) with r(rain) (red)",cex.main=0.9) # estimate logistic model with rain affecting r and k a=log(0.1) b=0.2 r=exp(a+b*rains) k1=log(1.25) k2=0.1 k=exp(k1+k2*rains) N0=Nhat[1] N=Ny.rain.k(N0,r,k,y,ys,Ca) est.logistic.rain.k=wild.rain.k.fit(N0,a,b,k1,k2,y,Nhat,sehat,rains,ys,Ca,control=control) est.logistic.rain.k$AIC plot(y+60,Nhat,ylim=range(lower,upper),xlab="Year",ylab=expression(paste(hat(italic(N))," (millions)"),cex=0.5)) segments(x0=y+60,y0=lower,y1=upper) lines(est.logistic.rain.k$ys[-1]+60,est.logistic.rain.k$all.N[-1],lty=1,col="blue") lines(y+60,est.logistic.rain$N,lty=1,col="red") title("Wildebeest logistic with r(rain) (red) & k(rain) (blue)",cex.main=0.9) # Look at logistic AICs logistic.AICs=c(est.logistic$AIC,est.logistic.rain$AIC,est.logistic.rain.k$AIC) names(logistic.AICs)=c("logistic","logistic.rain","logistic.rain.K") sort(logistic.AICs)-min(logistic.AICs) # Delta AIC #======================== Matrix Model Estimation ======================== # Fit with different combinations of fixed parameters # FIRST NOT using addional estimates of calf and adult survival; N0=0.27 phi1=0.65 rain1=0.5 phi2=0.8 rain2=0.1 lambda=0.4 parobj.l=pack.pars(N0,phi1,rain1,phi2,rain2,lambda,fixed="lambda") # constant lambda, calf & adult survival depending on last year's rainfall est.MM.l=wildbst.MM.fit(parobj.l,ys,y,Nhat,sehat,Ca,rains,yc,phi1hat,ya,phi2hat,addata=FALSE,control=control) est.MM.l$fit$convergence est.MM.l[1:3] est.MM.l$AIC plot(y+60,Nhat,ylim=range(lower,upper),xlab="Year",ylab=expression(paste(hat(italic(N))," (millions)"))) segments(x0=y+60,y0=lower,y1=upper,lwd=2) lines(est.MM.l$ys+60,est.MM.l$all.N,lty=4,lwd=2) title("Wildebeest MM: lambda fixed",cex.main=0.9) N0=0.27 phi1=0.65 rain1=0 phi2=0.9 rain2=0.01 lambda=0.4 parobj.lr1=pack.pars(N0,phi1,0,phi2,rain2,lambda,fixed=c("lambda","rain1")) # constant calf survival, constant lambda, adult survival depending on last year's rainfall est.MM.lr1=wildbst.MM.fit(parobj.lr1,ys,y,Nhat,sehat,Ca,rains,yc,phi1hat,ya,phi2hat,addata=FALSE,control=control) est.MM.lr1$fit$convergence est.MM.lr1[1:3] est.MM.lr1$AIC plot(y+60,Nhat,ylim=range(lower,upper),xlab="Year",ylab=expression(paste(hat(italic(N))," (millions)"))) segments(x0=y+60,y0=lower,y1=upper,lwd=2) lines(est.MM.lr1$ys+60,est.MM.lr1$all.N,lty=4,lwd=2) title("Wildebeest MM: lambda fixed, const calf survival",cex.main=0.9) #N0=0.24 #phi1=0.6 #rain1=1 #phi2=0.95 #rain2=0 #lambda=0.4 #parobj.lr2=pack.pars(N0,phi1,rain1,phi2,0,lambda,fixed=c("lambda","rain2")) # constant adult survival, constant lambda, calf survival depending on last year's rainfall #est.MM.lr2=wildbst.MM.fit(parobj.lr2,ys,y,Nhat,sehat,Ca,rains,yc,phi1hat,ya,phi2hat,addata=FALSE,control=control) #est.MM.lr2$fit$convergence #est.MM.lr2[1:3] #est.MM.lr2$AIC #plot(y+60,Nhat,ylim=range(lower,upper),xlab="Year",ylab=expression(paste(hat(italic(N))," (millions)"))) #segments(x0=y+60,y0=lower,y1=upper,lwd=2) #lines(est.MM.lr2$ys+60,est.MM.lr2$all.N,lty=4,lwd=2) #title("Wildebeest MM: lambda fixed, const adult survival",cex.main=0.9) N0=0.27 phi1=0.6 rain1=0 phi2=0.99 rain2=0 lambda=0.4 parobj.lr1r2=pack.pars(N0,phi1,rain1,phi2,0,lambda,fixed=c("lambda","rain1","rain2")) # constant adult survival, constant lambda, calf survival depending on last year's rainfall est.MM.lr1r2=wildbst.MM.fit(parobj.lr1r2,ys,y,Nhat,sehat,Ca,rains,yc,phi1hat,ya,phi2hat,addata=FALSE,control=control) est.MM.lr1r2$fit$convergence est.MM.lr1r2[1:3] est.MM.lr1r2$AIC plot(y+60,Nhat,ylim=range(lower,upper),xlab="Year",ylab=expression(paste(hat(italic(N))," (millions)"))) segments(x0=y+60,y0=lower,y1=upper,lwd=2) lines(est.MM.lr1r2$ys+60,est.MM.lr1r2$all.N,lty=4,lwd=2) title("Wildebeest MM: lambda fixed, const calf & adult survival",cex.main=0.9) # Compare Matrix model AICs MM.AICs=c(est.MM.l$AIC,est.MM.lr1$AIC, #est.MM.lr2$AIC, est.MM.lr1r2$AIC) names(MM.AICs)=c("MM.fixed l","MM.fixed l&r1", #"MM.fixed l&r2", "MM.fixed l&r1&r2") sort(MM.AICs)-min(MM.AICs) # Delta AIC Sys.sleep(3) # Compare AICs from both kinds of model: AICs=c(MM.AICs,logistic.AICs) sort(AICs)-min(AICs) Sys.sleep(3) # constant lambda, calf & adult survival depending on last year's rainfall N0=0.26 phi1=0.2 rain1=0 phi2=0.95 rain2=0 lambda=0.4 parobj.lr1r2.add=pack.pars(N0,phi1,rain1,phi2,rain2,lambda,fixed=c("lambda","rain1","rain2")) # constant lambda, calf & adult survival depending on last year's rainfall est.MM.lr1r2.add=wildbst.MM.fit(parobj.lr1r2.add,ys,y,Nhat,sehat,Ca,rains,yc,phi1hat,ya,phi2hat,addata=TRUE,control=control,hessian=TRUE) est.MM.lr1r2.add$fit$convergence est.MM.lr1r2.add[1:3] est.MM.lr1r2.add$AIC plot(y+60,Nhat,ylim=range(lower,upper),xlab="Year",ylab=expression(paste(hat(italic(N))," (millions)"))) segments(x0=y+60,y0=lower,y1=upper,lwd=2) lines(est.MM.lr1r2.add$ys+60,est.MM.lr1r2.add$all.N,lty=4,lwd=2) title(main="Wildebeest MM: lambda fixed, const calf & adult survival\n Using independent adult & calf survival estimates",cex.main=0.9) # ====================== Plots for Chapter 6 ============================ par(cex.lab=0.8) plot(ys+60,rains,type="l",xlab="Year",ylab="Rainfall",lty=2) title("Annual dry-season rainfall") est=est.MM.lr1r2.add # best model with additional survival estimates bs=wildbst.MM.bs.trajectory(est,ys,y,Nhat,sehat,Ca,rains,B=100,seed=1958) # bootstrap best MM model with additonal estimates par(mar=c(5,5,4,2)) plot(y+60,Nhat,ylim=range(lower,upper),xlab="Year",ylab=expression(paste(hat(italic(N))," (millions)"))) for(i in 1:dim(bs$N)[1]){ lines(est$ys+60,bs$N[i,],col="gray") } points(y+60,Nhat) segments(x0=y+60,y0=lower,y1=upper,lwd=2) lines(est$ys+60,est$all.N,lty=1) # best model with additional survival estimates #lines(ys+60,est.MM.lr2$all.N,lty=2) # best model without additional survival estimates lines(y+60,est.logistic$N,lty=3) # best logistic model title("Wildebeest plot for Chapter 6")
Section 6.4.2 Gray whales
Deterministic logistic population growth model for gray whale survey data along with age-structured deterministic matrix model and trajectories from parametric bootstrap resamples of estimated matrix model parameters. Note number of bootstrap replicate trajectories has been limited to 100 for illustration here.
Graywhale data
# ----- Data from Punt & Wade (2010) ------------------------------- y=c(68:80,85,86,88,93,94,96,98) # survey years y=y-min(y) # start years at zero # abundance estimates and CVs Nhat=c(13426,14548,14553,12771,11079,17365,17375,15290,17564,18377,19538,15384,19763,23499,22921,26916,15762, 20103,20944,21135) cvhat=c(0.094,0.080,0.083,0.081,0.092,0.079,0.082,0.084,0.086,0.080,0.088,0.080,0.083,0.089,0.081,0.058,0.067, 0.055,0.061,0.068) sehat=cvhat*Nhat # convert to thousands Nhat=Nhat/1000 sehat=sehat/1000 # Estimated normal confidence bounds for N for plotting: lower=Nhat-1.96*sehat upper=Nhat+1.96*sehat # Catch ys=67:98 # catch years # male catches: Cm=c(151,92,93,70,62,66,98,94,58,69,86,94,57,53,36,56,46,59,55,46,47,43,61,67,69,0,0,21,48,18,48,64) # female catches: Cf=c(223,109,121,81,91,116,80,90,113,96,101,90,126,129,100,112,125,110,115,125,112,108,119,95,100,0,0,23,43,25, 31,61) # convert to thousands Cm=Cm/1000 Cf=Cf/1000 Ca=Cm+Cf # total catches
Graywhale functions
# Logistic fit # ------------ gray.fit.logistic=function(N0,r,k,sigma.pe,y,Nhat,sehat,Ca,control=list(trace=0),hessian=FALSE) { # fits logistic model # transform parameters: pars=log(c(N0,r,k,sigma.pe)) fit=optim(par=pars,fn=gray.negllik.pe,y=y,Nhat=Nhat,sehat=sehat,Ca=Ca,control=control,hessian=hessian) # untransform parameters: N0=exp(fit$par[1]) r=exp(fit$par[2]) k=exp(fit$par[3]) pe=exp(fit$par[4]) # pack stuff in a list and return it est=list(N0=N0,r=r,k=k,sigma.pe=pe,N=gray.Ny(N0,r,k,y,Ca)) est$AIC=fit$value*2+2*length(fit$par) est$fit=fit return(est) } gray.negllik.pe=function(pars,y,Nhat,sehat,Ca){ # Evaluate logistic model negative log-likelihood nll=sum(log(sehat+exp(pars[4]))+(Nhat-gray.Ny(exp(pars[1]),exp(pars[2]),exp(pars[3]),y,Ca))^2/ (2*(sehat^2+exp(pars[4])^2))) return(nll) } gray.Ny=function(N0,r,k,y,Ca){ # Logistic model forward projections of all Ns from year min(y)-1 (taken to be N0) to max(y) ys=(min(y)-1):max(y) nys=length(ys) Ns=rep(N0,nys) for(i in 2:nys){ Nleft=Ns[i-1]-Ca[i-1] # remove catch first Ns[i]=Nleft+r*Nleft*(1-Nleft/k) # then project forward a year } N=Ns[which(is.element(ys,y))] return(N) } # Matrix state model functions: # ---------------------------- gray.make.Leslie=function(phi.j,phi.a,f0,fmax,K,N,z=1,nj=5,na=1){ m=(nj+na) if(length(N)!=1) stop(paste("N must be a scalar, but is of length:",length(N))) M=matrix(rep(0,m^2),nrow=m) f=f0+(fmax-f0)*(1-(N/K)^z) # M[(2:m),(1:(m-1))]=M[(2:m),(1:(m-1))]+diag(rep(phi.j,nj)f) M[1,((m-1):m)]=f*c(phi.j,phi.a) M[(2:m),(1:(m-1))]=diag(rep(phi.j,nj)) M[m,m]=phi.a # f=f0+(fmax-f0)*(1-(N/K)^z) # M[1,(nj+1):m]=f return(M) } fzero=function(phi.j,phi.a,amax=200){ # Equation (3) of Brandon & Wade (2006) return((1-phi.a)/(phi.j^5*(1-phi.a^(amax-5-1)))) } # functions to do bounded transformation and un-tranfromation of survival parameters: bound.logit.phi.a=function(p,pmin=0.94,pmax=0.99999999){ return(logit((p-pmin)/(pmax-pmin))) } bound.logit.phi.j=function(p,pmin=0.50,pmax=0.99999999){ return(logit((p-pmin)/(pmax-pmin))) } bound.invlogit.phi.a=function(x,pmin=0.94,pmax=0.99999999){ return((inv.logit(x)*(pmax-pmin)+pmin)) } bound.invlogit.phi.j=function(x,pmin=0.50,pmax=0.99999999){ return((inv.logit(x)*(pmax-pmin)+pmin)) } gray.tfm.pars=function(phi.j,phi.a,fmax,K,sigma.pe,N0){ tfmpars=c(bound.logit.phi.j(phi.j),bound.logit.phi.a(phi.a),logit(fmax),log(K),log(sigma.pe),log(N0)) names(tfmpars)=c("bound.logit(phi.j)","bound.logit(phi.a)","logit(fmax)","log(K)","log(sigma.pe)","log(N0)") return(tfmpars) } gray.untfm.pars=function(tfmpars){ pars=c(bound.invlogit.phi.j(tfmpars[1]), bound.invlogit.phi.a(tfmpars[2]), inv.logit(tfmpars[3]), exp(tfmpars[4]), exp(tfmpars[5]), exp(tfmpars[6])) pars=list(phi.j=pars[1],phi.a=pars[2],fmax=pars[3],K=pars[4],sigma.pe=pars[5],N0=pars[6]) return(pars) } gray.fit.MM=function(phi.j,phi.a,fmax,K,sigma.pe,N0,y,Nhat,sehat,Ca,amax=200,control=list(trace=0,maxit=1000),hessian=FALSE) { if(phi.j<0.50) stop("phi.j must be between 0.50 and 1, but it is:",phi.j) if(phi.j>1) stop("phi.j must be between 0.50 and 1, but it is:",phi.j) if(phi.a<0.94) stop("phi.a must be between 0.94 and 1, but it is:",phi.a) if(phi.a>1) stop("phi.j must be between 0.94 and 1, but it is:",phi.a) tfmpars=gray.tfm.pars(phi.j,phi.a,fmax,K,sigma.pe,N0) fit=optim(par=tfmpars,fn=gray.MM.negllik,y=y,Nhat=Nhat,sehat=sehat,Ca=Ca,amax=amax,control=control,hessian=hessian) pars=gray.untfm.pars(fit$par) f0=fzero(pars$phi.j,pars$phi.a,amax) # asymptotic birth rate N=gray.MM.Ny(pars$phi.j,pars$phi.a,pars$f0,pars$fmax,pars$K,z,pars$N0,y,Ca,amax,nj=5,na=1) ys=min(y):max(y) all.N=gray.MM.Ny(pars$phi.j,pars$phi.a,pars$f0,pars$fmax,pars$K,z,pars$N0,ys,Ca,amax,nj=5,na=1) # get rid of some inapplicable names that were carried over: names(pars$phi.j)=names(pars$phi.a)=names(pars$fmax)=names(pars$K)=names(pars$sigma.pe)=names(pars$N0)=names(f0)="" est=list(phi.j=pars$phi.j,phi.a=pars$phi.a,fmax=pars$fmax,K=pars$K,z=z,sigma.pe=pars$sigma.pe,N0=pars$N0, f0=f0,y=y,N=N,ys=ys,all.N=all.N,amax=amax) est$AIC=fit$value*2+2*length(fit$par) est$fit=fit est$amax=amax return(est) } gray.MM.negllik=function(tfmpars,y,Nhat,sehat,Ca,amax=200,nj=5,na=1){ pars=gray.untfm.pars(tfmpars) f0=fzero(pars$phi.j,pars$phi.a,amax) Ns=gray.MM.Ny(pars$phi.j,pars$phi.a,pars$f0,pars$fmax,pars$K,z=1,pars$N0,y,Ca,amax,nj,na) # model abundances if(any(Ns<0)) { nll=Inf } else { nll=sum(0.5*log(sehat^2+pars$sigma.pe^2)+(Nhat-Ns)^2/(2*(sehat^2+pars$sigma.pe^2))) } return(nll) } gray.MM.Ny=function(phi.j,phi.a,f0,fmax,K,z,N0,y,Ca,amax,nj=5,na=1){ ys=(min(y)-1):max(y) nys=length(ys) m=nj+na if(length(N0)!=1) stop(paste("N0 must be scalar but it is:",N0)) f0=fzero(phi.j,phi.a,amax) M=gray.make.Leslie(phi.j,phi.a,f0,fmax,K,N0,z,nj,na) # set initial Leslie matrix ev1=Re(eigen(M)$vectors[,1]) #first eigenvector sdbn=ev1/sum(ev1) # stable age distribution N0vec=N0*sdbn # intial abundance by age Nsa=matrix(rep(N0vec,nys),nrow=m) # intialise age-structured abundances Ns=rep(0,nys) # intialise total abundances negabund=FALSE for(i in 2:nys){ Nsa[,i]=Nsa[,(i-1)] # number at start of year i is number at end of year (i-1) Nsa[m,i]=Nsa[m,i]-Ca[i] # first remove adult catch Ntot=sum(Nsa[,i]) # calculate resulting total before survival, ageing and births M=gray.make.Leslie(phi.j,phi.a,f0,fmax,K,Ntot,z,nj,na) # set Leslie matrix for survival, ageing and births Nsa[,i]=M%*%Nsa[,i] # project ahead one step, implementing survival, ageing and births if(Nsa[m,i]<0) negabund=TRUE # catch is of adults, so only adults can go negative Ns[i]=sum(Nsa[,i]) # record total abundanc at end of year i } N=Ns[which(is.element(ys,y))] if(negabund) N=N*0-1 # flag for impossible abundances return(N) } gray.MM.bs.trajectory=function(est,y,Nhat,sehat,Ca,B,seed){ mu=est$fit$par amax=est$amax phi.j=est$phi.j phi.a=est$phi.a f0=fzero(phi.j,phi.a,amax) sigma=solve(est$fit$hessian) ys=min(y):max(y) nys=length(ys) N=matrix(rep(-1,B*nys),nrow=B) parmat=matrix(rep(0,B*5),nrow=B) colnames(parmat)=c("phi.j","phi.a","fmax","K","N0") if(!is.null(seed)) set.seed(seed) # calculate trajectories: b=bad.bs=bad.f0=0 while(b < B){ tfmpars=rmvnorm(1,mu,sigma) # sample from multivariate normal parameters # transform parameters: pars=gray.untfm.pars(tfmpars) phi.j=pars$phi.j phi.a=pars$phi.a fmax=pars$fmax K=pars$K N0=pars$N0 if(phi.j==1 & phi.a==1) f0=0 else f0=fzero(phi.j,phi.a,amax) # cat("phi.j, phi.a, f0:",phi.j, phi.a, f0,"\n") # Impose bounds on fmax: if(fmax<f0) { fmax=1.001*f0 # can't have fmax<f0 bad.f0=bad.f0+1 } try.N=gray.MM.Ny(phi.j,phi.a,f0,fmax,K,z=1,N0,ys,Ca,amax,nj=5,na=1) # abundances for all years if(!any(try.N<0) & (try.N[1]<=try.N[nys-1])) { # remove stupid trajectories b=b+1 N[b,]=try.N parmat[b,]=c(phi.j,phi.a,fmax,K,N0) } else { bad.bs=bad.bs+1 } } return(list(N=N,pars=parmat,bad.bs=bad.bs,bad.f0=bad.f0)) } #sim.MM=function(est,nsim,phi.j,phi.a,fmax,K,sigma.pe,N0,y,Nhat,sehat,Ca,amax=200, # control=list(trace=0,maxit=1000),hessian=FALSE)){ #}
Graywhale analysis
library(boot) # need logit() and inv.logit() functions, and library boot has them library(mvtnorm) # needed for parameteric bootstrap # Fit logistic model with additional error (extension .pe for "process error") # ------------------------------------------------------------------------- # set starting values r=0.17 k=25 sigma.pe=1 N0=14 N=gray.Ny(N0,r,k,y,Ca) # get series of estimates with starting parameters, just to check they are not stupid control=list(trace=0,maxit=5000) est.logistic=gray.fit.logistic(N0,r,k,sigma.pe,y,Nhat,sehat,Ca,control=control) # fit model est.logistic$fit$convergence # check convergence=0 est.logistic$AIC est.logistic[1:4] # look at estimates # plot: par(mar=c(5,5,4,2)) # else y-label falls off LHS of plot plot(y+68,Nhat,ylim=range(lower,upper),xlab="Year",ylab=expression(paste(hat(italic(N))," (thousands)"))) segments(x0=y+68,y0=lower,y1=upper,lwd=2) lines(y+68,est.logistic$N,lty=2) title("Gray whale logistic model") # Fit state model with additional error # ------------------------------------- # set starting values phi.j=0.999 phi.a=0.99 amax=200 fmax=0.4 K=22 z=1 sigma.pe=2 N0=12 f0=fzero(phi.j,phi.a,amax) N=gray.MM.Ny(phi.j,phi.a,f0,fmax,K,z,N0,y,Ca,amax,nj=5,na=1) # get series of estimates with starting parameters, just to check they are not stupid plot(x=67:86, y=N, xlab="Year", ylab="N (thousands)",main="Abundance trajectory using starting parameters") control=list(trace=0,maxit=5000) est.MM=gray.fit.MM(phi.j,phi.a,fmax,K,sigma.pe,N0,y,Nhat,sehat,Ca,amax,control=control,hessian=TRUE) est.MM$fit$convergence # check convergence=0 est.MM$AIC est.MM[1:8] # look at estimates # plot plot(y+68,Nhat,ylim=range(lower,upper),xlab="Year",ylab=expression(paste(hat(italic(N))," (thousands)"))) segments(x0=y+68,y0=lower,y1=upper,lwd=2) lines(est.MM$ys+68,est.MM$all.N,lty=1) title("Gray whale matrix model") # ====================== Plots for Chapter 6 ============================= # par(mar=c(5,5,4,2)) # else y-label falls off LHS of plot bs=gray.MM.bs.trajectory(est.MM,y,Nhat,sehat,Ca,B=100,seed=123456789) # bootstrap trajectories plot(y+68,Nhat,ylim=range(lower,upper),xlab="Year",ylab=expression(paste(hat(italic(N))," (thousands)"))) for(i in 1:dim(bs$N)[1]){ lines(est.MM$ys+68,bs$N[i,],col="gray") } points(y+68,Nhat) segments(x0=y+68,y0=lower,y1=upper,lwd=2) lines(est.MM$ys+68,est.MM$all.N,lty=1) lines(y+68,est.logistic$N,lty=2) # add logistic model estimates title("Gray whale plot for Chapter 6")
Section 6.4.2.2 Kalman filter estimates for gray whale analysis
The following code produces Figure 6.5 for a normal dynamic state-space model fitted to gray whale data.
Code for Kalman filter
Modified to apply to the Gray Whale Example. This includes handling missing data.
#Chapter 6 R Code for # Kalman Filter, (Negative log) Likelihood, and MLE code # See Chapter 4 R Code for KF, etc. # State process matrix is time-specific and is a function of n_{t-1} # hence have conditional Gaussian NDLM # n[t] = A_t x n[t-1] + "Q_t" # y[t] = 1' x n[t] + "R_t" where y[t] is scalar, sum of abundances # 28 October 2013 - 20 November 2013 # 8 November 2013: allowing for different optimization approaches # 1. Use Nelder-Mead with transformation of params: w/ phi.j and phi.a # possibly constrained to [0.5, 0.999] and [0.94, 0.999] # 2. Use BFGS w/ same transformation options # 3. Tried L-BFGS-B w/ no transformations, just box constraints [DIDN'T WORK: # due to the beta0 calculations going "out of bounds"; e.g., get beta0 < 0 # 16 November 2013 # occasional non-positive definite Q[t] when lambda_t < 0 which in turn happens # when beta0 < beta1 and K < Ntotal # - Trying different bounds on beta1 (max fecundity), K, and N0 # - Forcing lambda_t to be non-negative: this does not allow estimation # of Hessian, however; main problem is Ntotal > K in lambda.t calculation # 17 November 2013 # Creating a separate file which calls these routines to ease de-bugging # Whale-call.r # 19-20 November 2013: matching Catch years with Obs'n Error only model work # Bring in the data files and R functions library(boot) # need logit() and inv.logit() functions, and library boot has them library("mvtnorm") # ---------- Part 0: Read in the data ----------------------- # ----- Data from Punt & Wade (2010) ------------------------------- # My year numbering is 1 higher than David's my.y <- c(68:80,85,86,88,93,94,96,98) # survey years my.y <- (my.y-min(my.y))+1 # start years at 1, i.e., 1968 = 1, 1998 = 31 # abundance estimates and CVs Nhat <- c(13426,14548,14553,12771,11079,17365,17375,15290,17564,18377,19538, 15384,19763,23499,22921,26916,15762, 20103,20944,21135) cvhat <- c(0.094,0.080,0.083,0.081,0.092,0.079,0.082,0.084,0.086,0.080, 0.088,0.080,0.083,0.089,0.081,0.058,0.067, 0.055,0.061,0.068) names(Nhat) <- names(cvhat) <- as.character(my.y) sehat <- cvhat*Nhat # convert to thousands Nhat <- Nhat/1000 sehat <- sehat/1000 # Fill in missing values full.y <- min(my.y):max(my.y) num.years <- length(full.y) Nhat.filled <- rep(NA,num.years) names(Nhat.filled) <- as.character(full.y) cvhat.filled <- sehat.filled <- Nhat.filled ok <- match(names(Nhat),names(Nhat.filled)) Nhat.filled[ok] <- Nhat cvhat.filled[ok] <- cvhat sehat.filled[ok] <- sehat obs <- matrix(data=Nhat.filled,nrow=1,ncol=num.years,dimnames=list(NULL,full.y)) # Harvest: my.ys=68:98 # harvest years # male catches: Cm=c(92,93,70,62,66,98,94,58,69,86,94,57,53,36,56,46,59,55,46,47,43, 61,67,69,0,0,21,48,18,48,64) # female catches: Cf=c(109,121,81,91,116,80,90,113,96,101,90,126,129,100,112,125,110,115, 125,112,108,119,95,100,0,0,23,43,25,31,61) # convert to thousands Cm=Cm/1000 Cf=Cf/1000 my.Ca=Cm+Cf # total catches names(my.Ca) <- as.character(full.y)# Harvest begins with year 1 # Obs'n model matrix B <- matrix(data=rep(1,6),nrow=1,ncol=6) # Obs'n variance Rt <- sehat.filled^2 # ---------- Part 1: Utility functions ----------------------- # functions to do bounded transformation and un-transformation # of survival parameters: bound.logit.phi.j=function(p,pmin=0.50,pmax=0.9999){ return(logit((p-pmin)/(pmax-pmin))) } bound.invlogit.phi.j=function(x,pmin=0.50,pmax=0.9999){ return((inv.logit(x)*(pmax-pmin)+pmin)) } bound.logit.phi.a=function(p,pmin=0.70,pmax=0.9999){ # bound.logit.phi.a=function(p,pmin=0.94,pmax=0.9999){ return(logit((p-pmin)/(pmax-pmin))) } bound.invlogit.phi.a=function(x,pmin=0.70,pmax=0.9999){ # bound.invlogit.phi.a=function(x,pmin=0.94,pmax=0.9999){ return((inv.logit(x)*(pmax-pmin)+pmin)) } bound.logit.beta1 <- function(p,pmin=0.2,pmax=0.9999){ return(logit((p-pmin)/(pmax-pmin))) } bound.invlogit.beta1 <- function(x,pmin=0.2,pmax=0.9999){ return((inv.logit(x)*(pmax-pmin)+pmin)) } bound.log.K <- function(x,Kmin=10) { return(log(x-Kmin)) } bound.invlog.K <- function(x,Kmin=10) { return(exp(x)+Kmin) } bound.log.N0 <- function(x,N0max=20) { return(logit(x/N0max)) } bound.invlog.N0 <- function(x,N0max=20) { return(N0max*inv.logit(x)) } # ---------- Part 2: Preliminary functions for NDLM ------------------ gray.tfm.pars=function(params,tfm.option=1) { # slight modification of obs'n error only model phi.j <- params$phi.j phi.a <- params$phi.a beta1 <- params$beta1 K <- params$K N0 <- params$N0 sigma.pe <- params$sigma.pe if(tfm.option==1) { #tfmpars=c(bound.logit.phi.j(phi.j),bound.logit.phi.a(phi.a), # bound.logit.beta1(beta1),bound.log.K(K),bound.log.N0(N0),log(sigma.pe)) tfmpars=c(bound.logit.phi.j(phi.j),bound.logit.phi.a(phi.a), logit(beta1),log(K),log(N0),log(sigma.pe)) } else { # not bounding logit tfmpars=c( logit(phi.j),logit(phi.a), logit(beta1),log(K),log(N0),log(sigma.pe)) } names(tfmpars) <- c("inv.phi.j","inv.phi.a","inv.beta1", "inv.K","inv.N0","inv.sigma.pe") return(tfmpars) } gray.untfm.pars=function(tfmpars,tfm.option=1){ # slight modification of obs'n error only model if(tfm.option==1) { params=c( bound.invlogit.phi.j(tfmpars["inv.phi.j"]), bound.invlogit.phi.a(tfmpars["inv.phi.a"]), #bound.invlogit.beta1(tfmpars["inv.beta1"]), #bound.invlog.K(tfmpars["inv.K"]), #bound.invlog.N0(tfmpars["inv.N0"]), inv.logit(tfmpars["inv.beta1"]), exp(tfmpars["inv.K"]), exp(tfmpars["inv.N0"]), exp(tfmpars["inv.sigma.pe"])) } else { #not bounding logit params=c( inv.logit(tfmpars["inv.phi.j"]), inv.logit(tfmpars["inv.phi.a"]), inv.logit(tfmpars["inv.beta1"]), exp(tfmpars["inv.K"]), exp(tfmpars["inv.N0"]), exp(tfmpars["inv.sigma.pe"])) } names(params) <- NULL params=list(phi.j=params[1],phi.a=params[2],beta1=params[3],K=params[4], N0=params[5],sigma.pe=params[6]) return(params) } cor.matrix.fun <- function(VCOV) { temp <- VCOV se <- sqrt(diag(VCOV)) temp1 <- VCOV/se out <- t(t(temp1)/se) return(out) } beta0.calc <- function(phi.j,phi.a,amax=200){ # Equation (3) of Brandon & Wade (2006) # identical to fzero() function out <- (1-phi.a)/(phi.j^5*(1-phi.a^(amax-5-1))) return(out) } lambdat.calc <- function(beta0,beta1,Ntotal,K,z=1,nonneg.lambda=FALSE, verbose.option=FALSE) { # the fecundity rate parameter # modified 17 Nov 2013 to ensure lambda_t is at least beta0 # => Hessian cannot be calculated, however out <- beta0 + (beta1-beta0)*(1-(Ntotal/K)^z) if(nonneg.lambda) { #if(verbose.option & out<0) if(out<0) cat("lambda=",out,"Nttl=",Ntotal,"K=",K,"beta0=",beta0, "beta1",beta1,"\n") out <- max(beta0, out) } return(out) } At.calculation <- function(phi.j,phi.a,lambdat,dim.state) { # this calculates the A_t projection matrix in the State eq'n At <- matrix(data=0,nrow=dim.state,ncol=dim.state) At[1,5] <- lambdat*phi.j At[1,6] <- lambdat*phi.a At[2,1] <- At[3,2] <- At[4,3] <- At[5,4] <- At[6,5] <- phi.j At[6,6] <- phi.a return(At) } Qt.calculation <- function(phi.j,phi.a,lambdat,n.past.vector,harvest,dim.state){ # this calculates the Q_t variance-covariance matrix in the State eq'n Qt <- matrix(data=0,nrow=dim.state,ncol=dim.state) # Ages 2 thru 5 are simply binomial survival with binomial variance Qt[2,2] <- phi.j*(1-phi.j)*n.past.vector[1] Qt[3,3] <- phi.j*(1-phi.j)*n.past.vector[2] Qt[4,4] <- phi.j*(1-phi.j)*n.past.vector[3] Qt[5,5] <- phi.j*(1-phi.j)*n.past.vector[4] # Age 6 variance is sum of 2 independent binomial survival processes Qt[6,6] <- phi.j*(1-phi.j)*n.past.vector[5] + phi.a*(1-phi.a)*(n.past.vector[6]-harvest) # Age 1 variance is more involved: dealing with sum of 2 binomial survivals # (ages 5 and 6) and binomial birth process E.n.6 <- phi.j*n.past.vector[5]+ phi.a*(n.past.vector[6]-harvest) V.E.part <- lambdat^2*Qt[6,6] E.V.part <- lambdat*(1-lambdat)*E.n.6 Qt[1,1] <- V.E.part + E.V.part # Covariance between Ages 1 and 6 E.n1.x.n6 <- lambdat*(Qt[6,6]+E.n.6^2) E.n1.E.n6 <- lambdat*E.n.6*E.n.6 Qt[1,6] <- Qt[6,1] <- E.n1.x.n6-E.n1.E.n6 return(Qt) } # Deterministic model predictions using estimated initial abundance, mles, # and Leslie matrix Leslie.projection <- function(params,harvest,nj=5,na=1,verbose=FALSE) { # This year's harvest is subtracted from last year's abundance to create # an "initial" abundance in the current year, to which the Leslie matrix # is applied phi.j <- params$"phi.j" phi.a <- params$"phi.a" beta1 <- params$"beta1" K <- params$"K" N0 <- params$"N0" dim.state <- nj+na num.times <- length(harvest) pred.n.matrix <- matrix(data=NA,nrow=dim.state,ncol=num.times+1, dimnames=list(paste("age",1:dim.state),0:num.times)) beta0 <- beta0.calc(phi.j=phi.j,phi.a=phi.a,amax=200) lambdat <- lambdat.calc(beta0=beta0,beta1=beta1,Ntotal=N0,K=K,z=1) A0 <- At.calculation(phi.j=phi.j,phi.a=phi.a,lambdat=lambdat, dim.state=dim.state) ev1 <- Re(eigen(A0)$vectors[,1]) # first eigenvector sdbn <- ev1/sum(ev1) # stable age distribution N0vec <- N0*sdbn # initial abundance by age pred.n.matrix[,"0"] <- N0vec # initial abundance vector (year = 0) # the iterations for(i in 1:num.times) { if(verbose) {cat("\n-------Time period",i,"\n") } last.years.abundance <- pred.n.matrix[,as.character(i-1)] #Subtract last year's harvest from last year's adults last.years.abundance[nj+na] <- last.years.abundance[nj+na] - harvest[as.character(i)] Ntotal <- sum(last.years.abundance) lambdat <- lambdat.calc(beta0=beta0,beta1=beta1,Ntotal=Ntotal,K=K,z=1) At <- At.calculation(phi.j=phi.j,phi.a=phi.a,lambdat=lambdat, dim.state=dim.state) pred.n.matrix[,as.character(i)] <- At %*% last.years.abundance } out <- pred.n.matrix return(out) } # ---------- Part 3: Kalman Filter Function ----------------------- kf.missing <- function(params,obs,B,Rt,harvest,nj=5,na=1,verbose=FALSE) { #-- Kalman filter estimates of state vectors # Tailored to gray whale example # Added handling of missing obs'ns # and dynamic state projection and covariance matrices # Output: # Kalman filter estimates of state vector for static SSM # Also returns matrix of 1 step ahead predicted states and corresponding # array of 1 step ahead predicted covariance matrices # Input: # params is a vector of parameters # obs is a p x T matrix # At is a dynamic state transition matrix, here the "Leslie" matrix # Qt is a dynamic state covariance matrix, n[t] = A x n[t-1] + "Q" # B is the observation "linkage" matrix, here a simple summation across # the age classes # Rt is the obs'n covariance "matrix", y[t] = B x n[t] + "Rt", here a # simple scalar, "std error" squared # N0 is the initial state "total" # 16 Nov 2013, var.stabilizer is a constant added to the diagonal of Qt # to try to avoid non-singularity; by phi.j <- params$"phi.j" phi.a <- params$"phi.a" beta1 <- params$"beta1" K <- params$"K" N0 <- params$"N0" sigma.pe <- params$"sigma.pe" dim.state <- nj+na num.times <- dim(obs)[2] pred.n <- update.n <- matrix(0,nrow=dim.state,ncol=1) cov.pred.n <- cov.update.n <- diag(x=rep(0,dim.state)) filter.n.matrix <- matrix(data=0,nrow=dim.state,ncol=num.times) #-- the following are used for likelihood evaluation pred.n.matrix <- matrix(data=0,nrow=dim.state,ncol=num.times) cov.pred.array <- array(data=0,dim=c(dim.state,dim.state,num.times)) #Initialization beta0 <- beta0.calc(phi.j=phi.j,phi.a=phi.a,amax=200) lambdat <- lambdat.calc(beta0=beta0,beta1=beta1,Ntotal=N0,K=K,z=1) if(verbose) cat("0",lambdat,"Nttl=",N0,"K=",K,"\n") A0 <- At.calculation(phi.j=phi.j,phi.a=phi.a,lambdat=lambdat, dim.state=dim.state) ev1 <- Re(eigen(A0)$vectors[,1]) # first eigenvector sdbn <- ev1/sum(ev1) # stable age distribution N0vec <- N0*sdbn # initial abundance by age # the iterations update.n <- N0vec identity.mat <- diag(x=1,nrow=dim.state) for(i in 1:num.times) { if(verbose) {cat("\n-------Time period",i,"\n") } last.years.abundance <- update.n # Subtract this year's harvest from last year's adult abundance last.years.abundance[nj+na] <- last.years.abundance[nj+na] - harvest[i] Ntotal <- sum(last.years.abundance) lambdat <- lambdat.calc(beta0=beta0,beta1=beta1,Ntotal=Ntotal,K=K,z=1) if(verbose) cat(i,"lambda=",lambdat,"Nttl=",Ntotal,"K=",K,"\n") At <- At.calculation(phi.j=phi.j,phi.a=phi.a,lambdat=lambdat, dim.state=dim.state) Qt <- Qt.calculation(phi.j=phi.j,phi.a=phi.a,lambdat=lambdat, n.past.vector=update.n,harvest=harvest[i],dim.state=dim.state) #Prediction step pred.n <- At %*% last.years.abundance cov.pred.n <- At %*% cov.update.n %*% t(At) + Qt x.p <- det(cov.pred.n) # the following problem "almost always" is due to negative lambda_t, from # time period 22, and this resulted from K < Ntotal # if(x.p <= 0) { # cat("\n ------------singular prediction covariance matrix at obs=",i,"det=",x.p,"\n") # } # Updating ("filtering") step if(all(!is.na(obs[,i]))) { Kalman.gain <- cov.pred.n %*% t(B) %*% solve(B %*% cov.pred.n %*% t(B) + Rt[i] + sigma.pe^2) update.n <- pred.n + Kalman.gain %*% (obs[,i,drop=FALSE]-B%*%pred.n) cov.update.n <- (identity.mat-Kalman.gain %*% B) %*% cov.pred.n } else { #Missing Obs'ns update.n <- pred.n cov.update.n <- cov.pred.n } filter.n.matrix[,i] <- update.n pred.n.matrix[,i] <- pred.n cov.pred.array[,,i] <- cov.pred.n } out <- list(N0vec=N0vec,filter.n.matrix=filter.n.matrix, pred.n.matrix=pred.n.matrix, cov.pred.array=cov.pred.array) return(out) } # ---------- Part 4: Negative Log-Likelihood Function ------ like.fn.gray.whale.SSM <- function(theta,obs,B,Rt,harvest,nj=5,na=1, tfm.option=1,verbose.opt=FALSE) { #theta is a vector of the parameters transformed to real number line params <- gray.untfm.pars(theta,tfm.option=tfm.option) #print("params");print(unlist(params)) sigma.pe.sq <- params$sigma.pe^2 if(verbose.opt) print(unlist(params)) num.times <- length(obs) out <- kf.missing(params=params,obs=obs,B=B,Rt=Rt,harvest=harvest,nj=nj,na=na, verbose=verbose.opt) pred.n.matrix <- out$pred.n.matrix cov.pred.array <- out$cov.pred.array N0vec <- out$N0vec log.l <- 0 for(i in 1:num.times) { if(!is.na(obs[,i])) { #skip all years with missing abundance estimates if(i > 1) { # THIS WAS WHERE THE BUG WAS, if(i > i), thus never satisfied cond.mean.vector <- B %*% pred.n.matrix[,i] cond.cov.matrix <- B %*% cov.pred.array[,,i] %*% t(B) + Rt[i] + sigma.pe.sq } else { cond.mean.vector <- B %*% cbind(N0vec) cond.cov.matrix <- Rt[i] + sigma.pe.sq } if(verbose.opt) cat(i,as.vector(cond.mean.vector), as.vector(cond.cov.matrix),obs[,i],"\n") # if multivariate obs'ns #log.l.comp <- dmvnorm(x=obs[,i+1], mean=cond.mean.vector, # sigma=cond.cov.matrix, log=TRUE) # The following is just an error check; in MVN case check for non-pos def if(is.na(sqrt(cond.cov.matrix))) { cat("problem with variance ",i,cond.cov.matrix,"\n") print("cov pred array"); print(cov.pred.array[,,i]) print(unlist(params)) } log.l.comp <- dnorm(x=obs[,i], mean=cond.mean.vector, sd=sqrt(cond.cov.matrix), log=TRUE) log.l <- log.l + log.l.comp } } out <- -log.l return(out) } ``` ### Fitting NDLM to gray whale data This code sources the functions, data, etc. and then "executes" the code for fitting the NDLM to gray whale data and displays results ```{r sect6.4.2.2, echo=TRUE, warnings=FALSE} # - Testing the KF function with values estimated in the obs'n error only model --# params.obs.err <- as.list(c("phi.j"=0.9998996,"phi.a"=0.9400000,"beta1"=0.3418253, "K"=21.6390683,"N0"=11.6035525,"sigma.pe"=2.4018017)) out <- kf.missing(params=params.obs.err,obs=obs,B=B,Rt=Rt,harvest=my.Ca,nj=5,na=1, verbose=FALSE) print(as.data.frame(out$filter.n.matrix)) obs.pred <- Leslie.projection(params=params.obs.err,harvest=my.Ca, nj=5,na=1,verbose=FALSE) #-- testing log likelihood function my.tfm.option <- 1 # option =1 is bounded parameters, option=2 is unbounded theta<- gray.tfm.pars(params.obs.err,my.tfm.option) like.fn.gray.whale.SSM(theta,obs,B,Rt, harvest=my.Ca,nj=5,na=1,tfm.option=my.tfm.option,verbose.opt=FALSE) # ------- Calculating MLEs: with different sets of initial values ------ initial.1 <- list(phi.j=0.88,phi.a=0.90,beta1=0.4,K=28,N0=12,sigma.pe=2.4) # check on legitimacy of initial values wrt beta0 = minimum fecundity cat("min fecundity=", beta0.calc(phi.j=initial.1$phi.j, phi.a=initial.1$phi.a,amax=200), "max fecundity=",initial.1$beta1,"\n") initial <- initial.1 my.tfm.option <- 1 # option =1 is bounded parameters, option=2 is unbounded my.optim.option <- "Nelder-Mead" #"Nelder-Mead" # "BFGS" #When use BFGS sometimes get singular Hessian my.control <- list(trace=0,maxit=5000) my.lower <- -Inf my.upper <- Inf theta <- gray.tfm.pars(initial,tfm.option=my.tfm.option) like.fn.gray.whale.SSM(theta,obs,B,Rt, harvest=my.Ca,nj=5,na=1,tfm.option=my.tfm.option,verbose.opt=FALSE) mle.1 <- optim(par=theta,fn=like.fn.gray.whale.SSM, hessian=TRUE, method=my.optim.option,lower=my.lower,upper=my.upper, tfm.option=my.tfm.option, control=my.control, obs=obs,B=B,Rt=Rt,harvest=my.Ca,nj=5,na=1,verbose.opt=FALSE) param.est.1 <- unlist(gray.untfm.pars(tfmpars=mle.1$par, tfm.option=my.tfm.option)) vcov.1 <- solve(mle.1$hessian) out.1 <- kf.missing(params=as.list(param.est.1), obs=obs,B=B,Rt=Rt,harvest=my.Ca,nj=5,na=1,verbose=FALSE) cat("min fecundity=", beta0.calc(phi.j=param.est.1["phi.j"], phi.a=param.est.1["phi.a"],amax=200), "max fecundity=",param.est.1["beta1"],"\n") print(mle.1$convergence) print(mle.1$value) print(param.est.1) print(vcov.1) round(sqrt(diag(vcov.1))/mle.1$par*100,2) round(cor.matrix.fun(VCOV=vcov.1),2) filtered.states.1 <- out.1$filter.n.matrix print(filtered.states.1) print(out.1$pred.n.matrix) temp <- ts(t(filtered.states.1),start=68) plot(temp) total.states.1 <- apply(filtered.states.1,2,sum) total.ts.1 <- ts(total.states.1,start=68) avg.pred.1 <- Leslie.projection(params=as.list(param.est.1),harvest=my.Ca, nj=5,na=1,verbose=FALSE) avg.pred.1 <- apply(avg.pred.1[,-1],2,sum) # add Obs'n error only predictions obs.err.only <- Leslie.projection(params=params.obs.err,harvest=my.Ca,nj=5,na=1,verbose=FALSE) obs.err.only <- apply(obs.err.only[,-1],2,sum) my.ylim <- range(c(total.states.1,obs,avg.pred.1,obs.err.only),na.rm=TRUE) plot(total.ts.1,xlab="Year",ylab="Abundance",ylim=my.ylim,type='b',pch=1) lines(68:98,obs,col="red",lty=2,type='b',pch=2) lines(68:98,avg.pred.1,col="blue",lty=3,type='b',pch=3) lines(68:98,obs.err.only,col="purple",lty=4,type='b', pch=4) legend("topleft",legend=c("Filtered","Obsns","Projected SSM", "Projected Obsn Error Only"), col=c("black","red","blue","purple"),lty=1:4,pch=1:4)