Chapter 6: Modelling population dynamics using closed-population abundance estimates

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)