model {
#likelihood
for(t in 1:T) {
y[t] ~ dnorm(mu[t], tau.err)
mu[t] <- beta + theta[t]
}
# prior for temporal effects
# RW prior for theta[t] - specified using car.normal with
#neighbours (t-1) and (t+1)
# for theta[2],....,theta[T-1], and neighbours (t+1) for
#theta[1] and (t-1) for theta[T]
theta[1:T] ~ car.normal(adj[], weights[], num[], tau)
beta ~ dflat()
# Specify weight matrix and adjacency matrix corresponding to RW(1) prior
# (Note - this could be given in the data file instead)
for(t in 1:1) {
weights[t] <- 1;
adj[t] <- t+1;
num[t] <- 1
}
for(t in 2:(T-1)) {
weights[2+(t-2)*2] <- 1;
adj[2+(t-2)*2] <- t-1
weights[3+(t-2)*2] <- 1;
adj[3+(t-2)*2] <- t+1;
num[t] <- 2
}
for(t in T:T) {
weights[(T-2)*2 + 2] <- 1;
adj[(T-2)*2 + 2] <- t-1;
num[t] <- 1
}
# Alternatively, a weight matrix and adjacency matrix
#corresponding to RW(2) prior can
# be specified or given in the data file
#(note, no need to change the prior distribution
# on theta, just the weights/adjacencies)
# for(t in 1:1) {
# weights[t] <- 2; adj[t] <- t+1
# weights[t+1] <- -1; adj[t+1] <- t+2; num[t] <- 2
# }
# for(t in 2:2) {
# weights[t+1] <- 2; adj[t+1] <- t-1
# weights[t+2] <- 4; adj[t+2] <- t+1
# weights[t+3] <- -1; adj[t+3] <- t+2; num[t] <- 3
# }
#for(t in 3:(T-2)) {
# weights[6+(t-3)*4] <- -1; adj[6+(t-3)*4] <- t-2
# weights[7+(t-3)*4] <- 4; adj[7+(t-3)*4] <- t-1
# weights[8+(t-3)*4] <- 4; adj[8+(t-3)*4] <- t+1
# weights[9+(t-3)*4] <- -1; adj[9+(t-3)*4] <- t+2; num[t] <- 4
# }
# for(t in (T-1):(T-1)) {
# weights[(T-4)*4 + 6] <- 2; adj[(T-4)*4 + 6] <- t+1
# weights[(T-4)*4 + 7] <- 4; adj[(T-4)*4 + 7] <- t-1
# weights[(T-4)*4 + 8] <- -1; adj[(T-4)*4 + 8] <- t-2; num[t] <- 3
# }
# for(t in T:T) {
# weights[(T-4)*4 + 9] <- 2; adj[(T-4)*4 + 9] <- t-1
# weights[(T-4)*4 + 10] <- -1; adj[(T-4)*4 + 10] <- t-2; num[t] <- 2
# }
# other priors
tau.err ~ dgamma(0.01, 0.01) # measurement error precision
sigma.err <- 1 / sqrt(tau.err)
sigma2.err <- 1/tau.err
tau ~ dgamma(0.01, 0.01) # random walk precision
sigma <- 1 / sqrt(tau)
sigma2 <- 1/tau
# include this variable to use in time series (model fit) plot
for(t in 1:T) { day[t] <- t }
}