## Clarifications regarding state-space 1st order autoregression time-series analysis

2013-01-15
2014-10-01
• I'm new to both bayesian and state-space modeling. My goal is to develop a state-space first order autoregressive population model:

*The State-Space approach combines two models: 1) ecological process model and 2) observation model

1)The ecological process model is: X[t] = B0 + (1+B1)X[t-1] + epsilon
X[t] = loge(abundance) in year t
X[t-1] = loge(abundance) in year t-1
B0 = Intinsic growth rate
B1 = 1st order density dependent autoregressive term
epsilon = process variance = a time-dependent stochastic term with mean=0

2)The Observation process model is: Y[t] = X[t] + observation error
Y[t] ~ normal(X[t], error)*

The model I have developed so far is:

```#model    {
Parameters and Priors ###
alpha.mu ~ dnorm(0.0, 1.0E-6);        #Prior on intrinsic rate of increase
tau.alpha ~ dgamma(1.0E-3, 1.0E-3);
sigma.alpha <- 1/sqrt(tau.alpha)      #Prior on precision of intrinsic rate of increase
beta ~ dnorm(0.0, 1.0E-6);            #Prior on 1st order density dependence
epsilon ~ dnorm(0.0, 1.0E-6);         #Prior on process variation
tau.pop ~ dgamma(1.0E-3, 1.0E-3);
sigma.pop <- 1/sqrt(tau.pop);         # Prior on precision of intrinsic rate of increase
tau.obs ~ dgamma(1.0E-3, 1.0E-3);
sigma.obs <- 1/sqrt(tau.obs)          # Prior on normal precision of observation error

### State-Space model ###

# System process model
for (t in 2:T){
alpha[t] ~ dnorm(alpha.mu, sigma.alpha);
log.pop.mu[t] <- alpha[t] + beta*log(pop[t-1]) + epsilon[t];
log.pop[t] ~ dnorm(log.pop.mu[t], sigma.pop);
}

# Observation process model
for (t in 2:T){
Y[t] ~ dnorm(log.pop[t], sigma.obs);
}
}
```

Sample input data:
"T" <- 45
"pop" <- c(8150, 5725, NA, NA, 4476, NA, 3842, 3172, 4305, 5543, 7281,
8215, 9981, 10529, 12607, 12014, 12828, 12680, 10838, NA, NA,
16019, NA, NA, NA, 16286, 17007, 18913, 16536, 14829, 12027,
12859, 17585, 19045, 16791, NA, NA, 11692, 11742, 14538, 13400,
11969, 9215, 8335, 9545)

Sample initial values:
"alpha.mu" <- 0.5
"tau.alpha" <- 0.5
"tau.obs" <- 1
"tau.pop" <- 1
"beta" <- 1
"epsilon" <- 1

```#Using the following rjags code:
library("rjags", "coda")
parameters = c("pop") # The parameter(s) to be monitored.
adaptSteps = 500 # Number of steps to "tune" the samplers.
burnInSteps = 1000 # Number of steps to "burn-in" the samplers.
nChains = 2 # Number of chains to run.
numSavedSteps=50000 # Total number of steps in chains to save.
thinSteps=1 # Number of steps to "thin" (1=keep every step).
nIter = ceiling( ( numSavedSteps * thinSteps ) / nChains ) # Steps per chain.
m <- jags.model("state-space.bug", data, inits, n.chains=nChains, n.adapt=adaptSteps)
```

I'm getting the following error: Compilation error on line (log.pop.mu[t] <- alpha[t] + beta*log(pop[t-1])). Unable to resolve node pop[3]

My clarifications are:
1. What is wrong with my code? Is the error message because of missing(NA) data?
2. Have I defined the precisions for the distributions of the parameters correctly?
3. Do I need the epsilon term in the code? (Other examples I have seen seem to ignore using it?)
4. Don't I need to reconvert my logarithmic estimates to population counts to be able to compare with the estimates provided in the input?
5. To compare the original time series with the bayesian estimated time-series, is the following code an appropriate approach:

```#coda code
codaSamples <- coda.samples(m, variable.names=parameters , n.iter=nIter , thin=thinSteps )
mcmcChain = as.matrix( codaSamples )
mcmc.means <- colMeans(mcmcChain, na.rm=TRUE, dims=1) #Do I need to take the means?
years <- seq(1961,2005, by=1)
plot(years,mcmc.means, type='l')
lines(years, data\$pop, col='blue')
```

Thank you.

• Francy Lisboa
2014-10-01

Hi. I'm also new to JAGS, but I realised you might have forgotten of specifying the node log.pop.mu. Where is it came from?

# System process model

```for (t in 2:T){
alpha[t] ~ dnorm(alpha.mu, sigma.alpha);
log.pop.mu[t] <- alpha[t] + beta*log(pop[t-1]) + epsilon[t];
log.pop[t] ~ dnorm(log.pop.mu[t], sigma.pop);
}
```

Maybe you should look at the this part of the model and see that log.pop.mu[t] is not defined anywhere.

Last edit: Francy Lisboa 2014-10-01