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)