Hi everybody,
I have set up a code for estimating a multivariate state-space model
with correlated process errors among states. The model is devised for
jointly estimating density-dependence and across-species environmental
correlation in a suite of 5 bird species. I suspect that something is
wrong with my code; although the estimates for the across-species
environmental correlation seem correct according to my expectations, the
estimates for the AR1 parameters and state estimates are clearly not
correct. Please, could anyone with expertise briefly check my code below
to see if I did something wrong in the coding? Really need help on this!
Thank you very much in advance,
Regards
Pablo
#### Here's the problem...
Specifically, my model takes the form:
X.mu_t+1 = beta0 + b1*X_t
X_t+1 ~ MN(X.mu_t+1, Omega)
Omega ~ Wis(R,5)
n_t+1 ~ N(X_t+1, tau2)
where:
1) X is a column vector of population sizes (size 5).
2) beta0 and b1 are column vectors of intercepts and AR1 coefficients,
respectively.
3) Omega is the square matrix (5 x 5) with the inverse variances of the
state-processes in the diagonal and inverse covariances among states in
the off-diagonal.
4) n_t+1 is a column vector with log-transformed population sizes.
Here`s my code:
model; {
# Observation model
for (i in 1:(N-1)) {
for(j in 1:K) {
n[i+1 , j] ~ dnorm(state[i+1 , j], itau2[j])
}
}
# Initial hidden states
for(j in 1:K){
state[1, j] ~ dunif(2, 13)
}
# Predictions of the multivariate state-space model
for (i in 1:(N-1)) {
state[i+1, 1:K] ~ dmnorm(state.mu[i+1 , 1:K], omega[ , ])
}
# Multivariate state-space model
for (i in 1:(N-1)) {
for(j in 1:K) {
state.mu[i+1, j] <- beta0[j] + b1[j] * state[i, j]
}
}
# Covariance matrix of the process errors
omega[1 : K , 1 : K] ~ dwish(R[ , ], 5)
sigma[1 : K , 1 : K] <- inverse(omega[ , ])
# Parameters of the multivariate state-space model
for (j in 1:K) {
beta0[j] ~ dnorm(0.0, 0.001)
b1[j] ~ dnorm(0.0, 0.001)
beta1[j] <- b1[j] - 1
itau2[j] ~ dgamma(0.001,0.001)
tau2[j] <- 1/itau2[j]
tau[j] <- pow(tau2[j], -2)
init.isigma2[j] ~ dgamma(0.1,1)
init.sigma2[j] <- 1/init.isigma2[j]
init.sigma[j] <- pow(init.sigma2[j], -2)
}
}
# Data
list(K = 5, N = 28, n = structure(.Data =
c(8.146,6.96,8.982,7.285,7.52,7.931,8.721,5.788,6.877,6.59,7.932,8.222,
7.019,8.997,6.515,2.319,4.467,10.746,5.963,6.495,7.89,6.918,6.527,7.602,
8.404,5.7,9.689,9.507,6.237,5.714,6.51,6.291,9.926,9.234,7.035,6.47,
6.051,8.121,10.045,8.491,7.726,5.924,6.59,9.413,6.037,8.356,8.917,10.538,
10.241,8.878,6.047,4.283,2.254,7.183,6.6,7.702,9.177,7.491,2.773,9.451,
7.798,6.141,10.533,9.18,4.621,5.146,9.181,6.37,5.265,9.352,8.876,4.877,
6.841,8.14,6.736,6.024,9.779,6.989,5.566,8.342,3.24,8.149,7.016,6.44,
7.889,6.535,7.816,8.673,6.461,8.76,7.643,7.152,5.704,7.841,8.37,11.059,
4.93,6.819,6.546,7.497,6.798,6.651,7.218,6.094,6.424,6.949,6.845,6.608,
8.375,6.496,7.214,4.23,6.237,7.684,7.747,3.09,6.371,9.261,5.659,7.846,
7.29,8.84,5.383,4.175,7.549,11.776,8.667,9.818,5.92,3.684,6.371,6.372,
7.656,5.92,7.126,8.071,5.271,5.345,5.864,6.33),.Dim = c(28, 5)),
R = structure(.Data = c(
1, 0, 0, 0, 0,
0, 1, 0, 0, 0,
0, 0, 1, 0, 0,
0, 0, 0, 1, 0,
0, 0, 0, 0, 1), .Dim = c(5, 5)))
# Inits
list(
beta0=c(0,0,0,0,0),b1=c(0,0,0,0,0),init.isigma2=c(10,10,10,10,10),
itau2=c(10,10,10,10,10),
omega=structure(.Data = c(
1, 0, 0, 0, 0,
0, 1, 0, 0, 0,
0, 0, 1, 0, 0,
0, 0, 0, 1, 0,
0, 0, 0, 0, 1), .Dim = c(5, 5)))
-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list
