Farshid Ahrestani
2013-01-15
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") data <- read.jagsdata("data.R") inits <- read.jagsdata("inits.R") 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?
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.