I am trying to convert a model from the BUGS language (implemented in JAGS) but am having trouble, maybe in part because I’m not clear exactly what’s happening in the JAGS model. The model is from Schofield et al. 2016: https://www.tandfonline.com/doi/abs/10.1080/01621459.2015.1110524
The model is:
The current (with a bit more hierarchical complexity than described above) JAGS model is:
model{
for(i in 1:M){
for(j in f[i]:l[i]) {
y[i,j] ~ dnorm(alpha0[i] + alpha1[i]*age[i,j] + eta[j], 1/sd_y[i]/sd_y[i])
}
alpha0[i] ~ dnorm(mu_a0, 1/sd_a0/sd_a0)
alpha1[i] ~ dnorm(mu_a1, 1/sd_a1/sd_a1)T( , 0) # truncated norm to force negative
sd_y[i] ~ dt(0, 0.04, 3)T(0, ) # variance differs by tree
}
mu_a0 ~ dnorm(0, 0.0001)
mu_a1 ~ dnorm(0, 0.0001)T( , 0)
sd_a0 ~ dt(0, 0.04, 3)T(0, )
sd_a1 ~ dt(0, 0.04, 3)T(0, )
for(t in 1:Tea) {
x[t] ~ dnorm(mu_x, 1 / sd_x / sd_x)
eta[t] ~ dnorm(beta0 * x[t], 1 / sd_eta / sd_eta)
}
mu_x ~ dnorm(0, 0.0001)
sd_x ~ dt(0, 0.04, 3)T(0, )
beta0 ~ dnorm(0, 0.0001)
sd_eta ~ dt(0, 0.04, 3)T(0, )
}
y_{it} is a matrix with observations across all 500 values of t. The primary interest of the model is to backcast X_t (reconstruct historic climate). In the JAGS model X_t come in as data with values t 1-400 as NA
and then 401-500 as observed values. With this formulation the missing values are estimated as the historic climate. I am not sure if sd_x is being estimated from just the observed values or if it’s somehow being estimated across all the potential values X_t given the observed values of y_{it} and age_{it}.
However, in Stan with the separate data
, model
, and generated quantities
block I am not sure how to handle estimating the historic (unobserved) X_t values. I assume since X_t are brought in as data then I can’t then have X_t ~ normal(mu_x, sd_x)
. But if I wanted to put it in the generated quantities
, I’m not sure where sd_x would come from. Maybe the model would need to be reparameterized. Any suggestions would be appreciated.