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