Search code examples
rbayesianjagsrjags

Exponentiation of vector in JAGS model, Bayesian analysis in R


I have the following JAGS model for use in a Bayesian model in R. I am trying to estimate the posterior distribution for my variable "R". All variables but R are supposed to be deterministic nodes. Each variable, s_A, z_A, z_W, and d are vectors. While tau_s is a data.frame. TTD_aquifer and O2s_all are therefore expected to be a vector for each i.

model {
for (i in 1:N){
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- sum(O2s_all)/2


tau_s_bar[i] = (s_A[i]*z_A[i])/R[i]*log(z_A[i]/(z_A[i]-z_W[i]))


TTD_aquifer <- t((d[i]*sqrt(tau_s_bar[i]))/sqrt(4*3.14*d[i]*t(tau_s[,i]^3))*exp(-1*((d[i]*tau_s_bar[i])/(4*t(tau_s[,i])))*
                                                                      (1-t(tau_s[,i])/tau_s_bar[i])^2))

O2s_all <- t(O2_o[i]-k_o[i]*t(tau_s[,i]))*TTD_aquifer

# prior on R
 R[i] ~ dlnorm(-2, 1/(0.6)^2)
 }
# prior on tau and sigma
 tau <- pow(sigma, -2)
 sigma ~ dunif(0, 100)
}

When I run this in jags.model() I get the following error: RUNTIME ERROR: Invalid vector argument to exp. So it looks like I cannot input a vector into exp() like you can in R. The equations for TTD_aquifer and O2s_all run fine in R for a deterministic example. How should I write my equation for TTD_aquifer in JAGS to avoid the exp issue?


Solution

  • In JAGS, inverse link functions like exp only take scalar arguments. You could change your model to this in order to use exp. Note that you will need to include an object in your data list that denotes how many rows are in the tau_s data frame. Since I do not know what your model is doing I have not checked to determine if your parentheses are in the correct location across all of your divisions and multiplications.

    model {
    for (i in 1:N){
        y[i] ~ dnorm(mu[i], tau)
        mu[i] <- sum(O2s_all[,i])/2
    
    tau_s_bar[i] <- (s_A[i]*z_A[i])/R[i]*log(z_A[i]/(z_A[i]-z_W[i]))
    
    for(j in 1:K){ # K = nrow of tau_s
    
    TTD_aquifer[j,i] <- t((d[i]*sqrt(tau_s_bar[i]))/
    sqrt(4*3.14*d[i]*t(tau_s[j,i]^3))*
    exp(-1*((d[i]*tau_s_bar[i])/(4*t(tau_s[j,i])))*
    (1-t(tau_s[j,i])/tau_s_bar[i])^2))
    
    O2s_all[j,i] <- t(O2_o[i]-k_o[i]*t(tau_s[j,i]))*TTD_aquifer[j,i]
    
    } # close K loop
    # prior on R
     R[i] ~ dlnorm(-2, 1/(0.6)^2)
     }
    # prior on tau and sigma
     tau <- pow(sigma, -2)
     sigma ~ dunif(0, 100)
    }
    

    As TTD_aquifer and 02s_all should be a vector for each i, they should then be a two dimensional matrix of the same size as tau_s for each step in an MCMC chain. If you have a big dataset (i.e., big N and K) and are running this model for many iterations, tracking those derived parameters will take up considerable memory. Just something to keep in mind if you are running this on a computer without sufficient RAM. Thinning the chain is one way to help deal with the computational intensity of tracking said parameters.