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")
    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
    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