# Predict values from complex Rjags model

It's the first time I'm working with R2Jags, MCM chains and Bayesian models and I'm having trouble to compute the predicted values for my model. The model is based on research by Hallmann et al. 2017, specifically the basic model (R code can be found in the appendix of the paper). Just as a little background, I'm doing this for a class that's focused on reproducibility in science.

Here's my jags code (slightly edited/renamed)

cat("model{
## Likelihood function for the latent expected daily biomass
for (i in 1:n) {
m_bio[i] ~ dnorm(sum(y[tau1[i]:tau2[i]]), sig_sq[i])
sig_sq[i] <- 1/Var[i]
Var[i] <- sum(vr[tau1[i]:tau2[i]])
}

## Likelihood function for muHat, it's dependent function and variance
for (i in 1:ndaily) {
z[i] <- exp(y[i])
y[i] <- g_intcp + log.lambda * year[i] + c[1] * daynr[i] + c[2] * daynr2[i] +
c[3] * daynr[i] * year[i] + c[4] * daynr2[i] * year[i] + b[loctype[i]] +
eps[plot[i]]
vr[i] <- exp(2 * y[i] + lvar) * (exp(lvar) - 1)
}

## Priors
g_intcp ~ dnorm(0, .01)
log.lambda ~ dnorm(0, .01)
b[1] <- 0
for( i in 2:3) {b[i] ~ dnorm(0, .01)}
for( i in 1:4) {c[i] ~ dnorm(0, .01)}
sdhat ~ dunif(0, 5)
lvar <- pow(sdhat, 2)
for (i in 1:nrandom) {
eps[i] ~ dnorm(0, tau.re)
}
tau.re <- pow(sd.re, -2)
sd.re ~ dunif(0, 1)
}


The mathematical equations can be found in the paper (eqs. 2, 3, 4, 6)

I'm trying to get values for m_bio. I've looked into other posts like this one: How to predict values using estimates from rjags / JAGS but couldn't 'extrapolate' the solution from the simple regression formula (a + b * x) to mine.

I would start by extracting the mean of every parameter, with something along the lines of this:

jags$$BUGSoutput$$mean


The part where I'm stuck is the matrix multiplication. What changes when I have more than one predictor ?

I'm also unsure about what I have to do after that. If I understood the model correctly, I need to do that with all three equations, i.e. z[i], y[i], and vr[i] to then be able to calculate the values of m_bio (and the dependent eqs. sig_sq and Var), is that correct ?

Hope I adequately conveyed my problem. Appreciate any and all replies :)

I believe the standard way to do predictions when using MCMC is to use each simulated set of parameters to form predictions, then study the distribution of the predictions across all of the simulations. If the model is linear in the parameters then this will give the same results as using the average parameters to form predictions, but in many interesting models (e.g. all GLMs) you'll get different results.

With a large dataset this might not be practical: if you have 10000 MCMC steps and 100000 data points, you'd need to store a billion values. In that case you might be able to do predictions only at a subset of "typical" data points, and you can certainly thin the output without much loss.