Modelling Population Dynamics: Chapter 4

Fitting state space models

Kalman Filter, (Negative log) Likelihood, and MLE code

for NDLM with 2 component state and observation vectors (See Table 4.1 and Section 4.4.1)

# ---------- Part 1:  Simulating the NDLM -----------------------
#   n[t] = A x n[t-1] + "Q"
#   y[t] = B x n[t] + "R"
set.seed(22)
library(mvtnorm)

# State process matrix
A <- matrix(data=c(0.0, 1.5, 0.3, 0.55),nrow=2,ncol=2,byrow=TRUE)
# State variance
sigma.p <- 0.5

# Obs'n model matrix
B <- diag(x=c(0.8,1.1),nrow=2)
# Obs'n variance
sigma.o <- 2.0

# Initial state values
n0 <- cbind(c(30,50))
num.times <- 10
n.vec <- y.vec <- matrix(data=0,nrow=2,ncol=num.times+1,
                         dimnames=list(c("Juveniles","Adults"),NULL))
n.vec[,1] <- n0

#-- Simulate states and observations
for(i in 1:num.times) {
  n.vec[,i+1] <- A %*% n.vec[,i] + rnorm(n=2,mean=0,sd=sigma.p)
  y.vec[,i+1] <- B %*% n.vec[,i+1] + rnorm(n=2,mean=0,sd=sigma.o)
}
print(n.vec)
print(y.vec)

# ---------- Part 2:  Kalman Filter Function -----------------------
kf <- function(obs,A,B,Q,R,n0) {
  #-- Kalman filter estimates of state vectors
  
  # 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:
  # obs is a p x [T+1] matrix
  # A is the state transition matrix
  # Q is the state covariance matrix, n[t] = A x n[t-1] + "Q"
  # B is the observation "linkage" matrix
  # R is the obs'n covariance matrix, y[t] = B x n[t] + "R"
  # n0 is the initial state vector
  
  num.times <- dim(obs)[2]-1  # start with time = 0
  dim.state <- dim(Q)[1]
  filter.n <- matrix(data=0,nrow=dim.state,ncol=num.times+1)
  pred.n <- update.n <- matrix(0,nrow=dim.state,ncol=1)
  cov.pred.n <- cov.update.n <- diag(x=rep(0,dim.state),nrow=dim.state)
  
#-- the following are used for likelihood evaluation
  pred.n.matrix <- matrix(data=0,nrow=dim.state,ncol=num.times+1)
  cov.pred.array <- array(data=0,dim=c(dim.state,dim.state,num.times+1))
  pred.n.matrix[,1] <- n0
  
# the iterations
  update.n <- n0
  filter.n[,1] <- update.n
  identity.mat <- diag(x=1,nrow=dim.state)
  for(i in 1:num.times) {
    pred.n <- A %*% update.n
    cov.pred.n <- A %*% cov.update.n %*% t(A) + Q
    
    Kalman.gain <- cov.pred.n %*% t(B) %*% solve(B %*% cov.pred.n %*% t(B) + R)
    
    update.n <- pred.n + Kalman.gain %*% (y.vec[,i+1,drop=FALSE]-B%*%pred.n)
    cov.update.n <- (diag(x=1,nrow=2)-Kalman.gain %*% B) %*% cov.pred.n
    
    filter.n[,i+1] <- update.n
    
    pred.n.matrix[,i+1] <- pred.n
    cov.pred.array[,,i+1] <- cov.pred.n
  }
  out <- list(filter.n=filter.n,pred.n.matrix=pred.n.matrix,
              cov.pred.array=cov.pred.array)
  return(out)
}

# ---------- Part 3: Applying KF with true parameter values  ------# 
Q <- diag(sigma.p^2,nrow=2)
R <- diag(sigma.o^2,nrow=2)
out <- kf(obs=y.vec,A=A,B=B,Q=Q,R=R,n0=n0)
print(out$filter.n, digits=3)

# Juvenile Obs'ns, True Values, and Filtered Values
temp <- round(rbind(y.vec[1,],n.vec[1,],out$filter.n[1,]))
dimnames(temp) <- list(c("y-juv","n-juv","n.kf"),0:(dim(y.vec)[2]-1))
print(temp)

# adult Obs'ns, True Values, and Filtered Values
temp <- round(rbind(y.vec[2,],n.vec[2,],out$filter.n[2,]))
dimnames(temp) <- list(c("y-ad","n-ad","n.kf"),0:(dim(y.vec)[2]-1))
print(temp)

# ---------- Part 3: Negative Log-Likelihood Function  ------ 
#  For state matrix A, with independent state components.
#  Need to have mvtnorm package installed
library(package=mvtnorm)
like.fn.2state.SSM <- function(theta,obs,B,Q,R,n0) {
  num.times  <- dim(obs)[2]-1  # time starts at 0
  num.states <- dim(Q)[1] 
  A <- matrix(data=theta,nrow=num.states,ncol=num.states,byrow=FALSE)
  out <- kf(obs=obs,A=A,B=B,Q=Q,R=R,n0=n0)
  pred.n.matrix  <- out$pred.n.matrix
  cov.pred.array <- out$cov.pred.array
  
  log.l <- 0
  for(i in 1:num.times) { 
    cond.mean.vector <- B %*% pred.n.matrix[,i]
    cond.cov.matrix  <- B %*% cov.pred.array[,,i] %*% t(B) + R
    log.l.comp <- dmvnorm(x=obs[,i+1], mean=cond.mean.vector, 
                          sigma=cond.cov.matrix, log=TRUE)
    log.l      <- log.l + log.l.comp
  }
  out <- -log.l 
  return(out)
}

#-- testing log likelihood function
Q <- diag(sigma.p^2,nrow=2)
R <- diag(sigma.o^2,nrow=2)
like.fn.2state.SSM(theta=as.vector(A),obs=y.vec,B=B,Q=Q,R=R,n0=n0)  

A.alt <- jitter(A)
like.fn.2state.SSM(theta=as.vector(A.alt),obs=y.vec,B=B,Q=Q,R=R,n0=n0)  

# ---------- Part 4: Calculating MLE of State Transition Matrix  ------ 
A.alt <- jitter(A)
mle <- optim(par=as.vector(A.alt),fn=like.fn.2state.SSM,
             obs=y.vec,B=B,Q=Q,R=R,n0=n0)

A.est <- matrix(data=mle$par,nrow=dim(Q)[1],ncol=dim(Q)[1],byrow=FALSE)
cat("MLE of State Transition Matrix \n"); print(A.est)
cat("True value of S.T.M. \n"); print(A)

# checks on optimization
like.fn.2state.SSM(theta=as.vector(A.est),obs=y.vec,B=B,Q=Q,R=R,n0=n0)
like.fn.2state.SSM(theta=as.vector(A),obs=y.vec,B=B,Q=Q,R=R,n0=n0)

# Applying Kalman Filter with mle for A   
out <- kf(obs=y.vec,A=A.est,B=B,Q=Q,R=R,n0=n0)

# Juvenile Obs'ns, True Values, and Filtered Values
temp <- round(rbind(y.vec[1,],n.vec[1,],out$filter.n[1,]))
dimnames(temp) <- list(c("y-juv","n-juv","n.kf"),0:(dim(y.vec)[2]-1))
print(temp)

# adult Obs'ns, True Values, and Filtered Values
temp <- round(rbind(y.vec[2,],n.vec[2,],out$filter.n[2,]))
dimnames(temp) <- list(c("y-ad","n-ad","n.kf"),0:(dim(y.vec)[2]-1))
print(temp)

Section 4.5.2.1 (WinBUGS code)

#(i) Model Statement
model {
#Prior distribution for parameters
beta1  ~ dnorm(0,0.01)
beta0  ~ dnorm(0,0.01)
sigma  ~ dunif(0.01,10)
tau <- 1/(sigma*sigma) #the precision
#Likelihood defined
for(i in 1:n) {
mu[i] <- beta0+beta1*x[i]
y[i]  ~ dnorm(mu[i],tau)
}
}

#(ii) Data Statement
list(n=15,
y=c( 29.4, 33.8, 24.8, 26.1, 31.6, 29.4, 14.6, 12.6, 29.3, 27.5,
28.9, 28.3 ,24.7, 28.5, 14.9),
x=c(8, 10, 7, 7, 9, 9, 3, 3, 8, 7, 8, 8, 7, 8, 3))		
#(iii) Initial Values Statement
list(beta0=0.5, beta1=1.0, sigma=3.0)	

Section 4.5.3 salmon SSM in WinBUGS

Save this block of WinBUGS code locally as C:/PopDynBook/BUGSCode/Ricker-Poisson-logN.txt)

model {
 #WinBUGS model code for Ricker model Poisson-LogN
 # assuming known obs'n noise CV

 # Priors
 alpha    ~ dunif(1,2.5)
 beta     ~ dunif(0.00001,0.001)
 n.init   ~ dunif(50,500)
 sigma.sq <- log(CV.obs*CV.obs+1)
 tau.obs  <- 1/sigma.sq

 mu[1]    <- alpha*n.init*exp(-beta*n.init)
 state[1] ~ dpois(mu[1])
 temp[1]  <- log(state[1])-sigma.sq/2
 obs[1]   ~ dlnorm(temp[1],tau.obs)

 for(i in 2:T) {
   mu[i]    <- alpha*state[i-1]*exp(-beta*state[i-1])
   state[i] ~  dpois(mu[i])
   temp[i]  <- log(state[i])-sigma.sq/2
   obs[i]   ~  dlnorm(temp[i],tau.obs)
 }
}
}

Load the following code into R

library(R2WinBUGS)
input.data <- list(T=num.times, obs=estimates,CV.obs=CV.obs)
init.value.generator <- function(num.starts=1,n,alpha.params=c(1,2),
      beta.params=c(0.00001,0.001),n.init.params=c(10,2000)) {
 init.values <- NULL
 for(i in 1:num.starts) {
     alpha<- runif(1,alpha.params[1],alpha.params[2])
     beta <- runif(1,beta.params[1],beta.params[2])
     n.init <- runif(1,n.init.params[1],n.init.params[2])
     state <- numeric(n)
     state[1] <- rpois(1,alpha*n.init*exp(-beta*n.init))
     for(j in 2:n)
      state[j] <-rpois(1,alpha*state[j-1]*exp(-beta*state[j-1]))
     init.values[[i]] <-
          list(alpha=alpha, beta =beta,n.init=n.init,state=state)
  }
 return(init.values)
}

init.values <- init.value.generator(num.starts=3,n=num.times,
      n.init.params=c(0.5*estimates[1],1.2*estimates[1]))

params <- c("alpha","beta","n.init","state")

out.ssm <- bugs(data=input.data, inits=init.values,
  parameters.to.save=params, model.file=
  "C:/PopDynBook/BUGSCode/Ricker-Poisson-logN.txt",
 n.chains=3, n.iter=50000, n.burnin=10000,n.thin=10,debug=TRUE)

Section 4.5.4 WinBUGS model statement for the BRS SSM

model {
 #WinBUGS model code for BRS model (Bin-Bin-Bin)-LogN assuming known obs'n noise CV

 # Priors
 #Survival parameters
 phi1     ~ dunif(0.05,0.95)
 phi2     ~ dunif(0.05,0.95)

 #Growth
 Pi       ~ dunif(0.05,0.95)
 not.Pi   <- 1-Pi

 #Birth
 lambda   ~ dunif(0.05,0.95)

 n1.init  ~ dunif(10,500)
 n2.init  ~ dunif(10,500)

 sigma.sq <- log(CV.obs*CV.obs+1)
 tau.obs  <- 1/sigma.sq

 #Year 1 calculations
 u1s[1]   ~ dbin(phi1,n1.init)
 u2s[1]   ~ dbin(phi2,n2.init)

 u1r[1]   ~ dbin(not.Pi,u1s[1])
 n2[1]    <- u2s[1] + u1s[1]-u1r[1]

 births[1] ~ dbin(lambda,n2[1])
 n1[1]    <- u1r[1] + births[1]

 temp1[1]    <- log(n1[1]) - sigma.sq/2
 obs.imm[1]   ~ dlnorm(temp1[1],tau.obs)
 temp2[1]    <- log(n2[1]) - sigma.sq/2
 obs.mat[1]  ~ dlnorm(temp2[1],tau.obs)

 #Iterate over remaining years
 for(i in 2:T) {
   u1s[i] ~ dbin(phi1,n1[i-1])
   u2s[i] ~ dbin(phi2,n2[i-1])

   u1r[i] ~ dbin(not.Pi,u1s[i])
   n2[i] <- u2s[i] + u1s[i]-u1r[i]

   births[i] ~ dbin(lambda,n2[i])
   n1[i] <- u1r[i] + births[i]

  temp1[i]    <- log(n1[i]) - sigma.sq/2
  obs.imm[i]   ~ dlnorm(temp1[i],tau.obs)
  temp2[i]    <- log(n2[i]) - sigma.sq/2
  obs.mat[i]  ~ dlnorm(temp2[i],tau.obs)
 }
}

R code for generating initial values to fit BRS SSM with WinBUGS.

init.value.generator <- function(num.starts=1,n,
    phi1.params=c(0.1,0.9),  phi2.params=c(0.1,0.9),
    Pi.params = c(0.1,0.9), lambda.params=c(0.1,0.9),
    n1.init.params=c(10,500), n2.init.params=c(10,500)) {
 init.values <- NULL
 for(i in 1:num.starts) {
     phi1   <- runif(1,phi1.params[1],phi1.params[2])
     phi2   <- runif(1,phi2.params[1],phi1.params[2])
     Pi     <- runif(1,Pi.params[1],Pi.params[2])
     lambda <- runif(1,lambda.params[1],lambda.params[2])

     n1.init <- round(runif(1,n1.init.params[1],n1.init.params[2]))
     n2.init <- round(runif(1,n2.init.params[1],n2.init.params[2]))
     n1 <- n2 <- u1s <- u2s <- u1r <- births <- numeric(n)

     u1s[1]    <- rbinom(1,size= n1.init,prob=phi1)
     u2s[1]    <- rbinom(1,size= n2.init,prob=phi2)

     u1r[1]    <-rbinom(1,size=u1s[1],prob=1-Pi)
     n2[1]     <- u2s[1] + u1s[1]-u1r[1]
     births[1] <- rbinom(1,size=n2[1],prob=lambda)
     n1[1]     <- u1r[1]+ births[1]

     for(j in 2:n) {
      u1s[j] <- rbinom(1,n1[j-1],phi1)
      u2s[j] <- rbinom(1,n2[j-1],phi2)

      u1r[j] <- rbinom(1,u1s[j],1-Pi)
      n2[j]  <- u2s[j] + (u1s[j]-u1r[j])

      births[j] <- rbinom(1,n2[j],lambda)
      n1[j] <- u1r[j] + births[j]
    }

     init.values[[i]] <-
          list(phi1=phi1, phi2 =phi2,Pi=Pi,lambda=lambda,
             n1.init=n1.init,n2.init=n2.init,
            u1s=u1s,u2s=u2s,u1r=u1r,births=births)
  }
 return(init.values)
}

init.values <- init.value.generator(num.starts=3,n=num.yrs,
    phi1.params=c(0.45,0.9),  phi2.params=c(0.55,0.9),
    Pi.params = c(0.5,0.9), lambda.params=c(0.7,0.85),
    n1.init.params=c(40,100), n2.init.params=c(40,100))

Section 4.5.5.2 Sequential importance sampling

R code for univariate coho salmon SSM (revised January 2016)

# Chapter SIS example complete
# Section 4.5.5.2 Sequential importance sampling
# R code for using SIS to fit univariate coho salmon SSM

ricker.sis <- function(N=1000,params,obs,CV.obs=0.3,smth.param=0.95) {
 r.a <- params[,1]; r.b <- params[,2]; n0 <- params[,3]
 obs.sd <- sqrt(log(CV.obs^2+1))
 num.times <- length(obs)
 states <- matrix(NA,nrow=N,ncol=num.times)
 uni.tally <- 1:N

 states[,1] <- rpois(n=N,lambda=r.a*n0*exp(-r.b*n0))
 wts <- dlnorm(x=obs,meanlog=log(states[,1])-obs.sd^2/2,obs.sd)
 for(i in 2:num.times) {
   states[,i] <- rpois(n=N,lambda=r.a*states[,i-1]*exp(-r.b*states[,i-1]))
   wts <- wts*dlnorm(x=obs[i],meanlog=log(states[,i])-obs.sd^2/2,obs.sd)
   new.set <- sample(1:N,N,replace=TRUE,prob=wts)
   uni.tally <- uni.tally[new.set]
   wts <- 1 #re-initialize wts
   states <- states[new.set,]
   params <- params[new.set,]

   #kernel smoothing of params
   mean.param <- apply(params,2,mean)
   cov.param  <- cov(params)
   smth.params <- smth.param*params + (1-smth.param)*mvrnorm(n=N,mu=cbind(mean.param),Sigma=cov.param)
   params <- smth.params
   r.a <- params[,1]; r.b <- params[,2]; n0 <- params[,3]
 }
 out <- list(states=states,params=params,tally=uni.tally)
 return(out)
}
 
# simulate coho data
set.seed(2345)
num.times <- 20
Ricker.a <- 1.5
Ricker.b <- 0.0003
equilibrium <-  log(Ricker.a)/Ricker.b
cat("equilibrium=", round(equilibrium),"\n")

initial.n <- round(0.1*equilibrium)
juveniles <- estimates <- numeric(num.times)
juveniles[1] <- rpois(1,lambda = Ricker.a*initial.n*exp(-Ricker.b*initial.n))
CV.obs <- .30
obs.sd <- sqrt(log(CV.obs^2+1))
estimates[1] <- round(rlnorm(1,meanlog=log(juveniles[1]),sdlog=obs.sd))
for(i in 2:num.times) {
 juveniles[i] <- rpois(1,lambda=Ricker.a*juveniles[i-1]*exp(-Ricker.b*juveniles[i-1]))
 estimates[i] <-  round(rlnorm(1,meanlog=log(juveniles[i])-obs.sd^2/2,sdlog=obs.sd))
 }

# Look at state and obs'n time series
my.ylim <- range(c(juveniles,estimates) )
plot(as.ts(juveniles),xlab="Time",ylab="",pch="s",col="blue",type="b",
     main="Simulated Coho Salmon Population and Estimates",ylim=my.ylim)
lines(as.ts(estimates),pch="o",col="red",type="b")

# set-up for SIS
num.particles <- 200000
set.seed(814)
init.params <-  matrix(NA,nrow=num.particles,ncol=3,dimnames=list(NULL,c("alpha","beta","n0")))
init.params[,1] <- runif(num.particles,min=1,max=2.5)
init.params[,2] <- runif(num.particles,min= .00001,max=0.001 )
init.params[,3] <- runif(num.particles,min=50,max=500)
dimnames(init.params) <- list(NULL,c("alpha","beta","n0"))

par(mfrow=c(2,2))
for(i in 1:3) hist(init.params[,i],main=dimnames(init.params)[[2]][i])
par(mfrow=c(1,1))

# run SIS
# load mvrnorm package
library(MASS)

temp <- ricker.sis(N=num.particles,params=init.params,obs=estimates,CV.obs=CV.obs,smth.param=0.95)

# Examine results
cat("particle depletion=",1-length(unique(temp$tally))/num.particles,"\n")

post.states <- apply(temp$states,2,mean)

par(mfrow=c(2,2),oma=c(0,0,3,0))
true.par <- c(Ricker.a,Ricker.b,initial.n)
for(i in 1:3) {
 my.xlim <- range(c(init.params[,i],temp$params[,i]))
 my.ylim <- range(c(density(init.params[,i])$y, density(temp$params[,i])$y))
 plot(density(init.params[,i]),xlim=my.xlim,ylim=my.ylim,
      main=dimnames(init.params)[[2]][i],xlab="")
 lines(density(temp$params[,i]),col=2,lty=2)
 text(c(true.par[i],mean(temp$params[,i])),c(0,0),"X",col=1:2)
 }
my.ylim <- range(c(juveniles,estimates,post.states))
plot(1:num.times,juveniles,xlab="Year",ylab="",ylim=my.ylim,type="l")
lines(1:num.times,estimates,lty=2,col=2)
lines(1:num.times,post.states,lty=3,col=3,lwd=2)
legend(7,0.45*max(my.ylim),legend=c("State","Observation","Smooth"),lty=1:3,col=1:3)
par(mfrow=c(1,1))
print(rbind(true.par,apply(temp$params,2,mean)))