JAGS: Hyper-prior distributions not converging to true values despite priors working

Modeling
1

I have a question regarding JAGS and I couIdn’t find the answer in any other forums, but was recommended to post modeling questions in this forum. I am trying to use R2jags to estimate the parameters from simulated data of an exponentially decaying variable for 10 subjects.

While the posterior distributions for the priors seem to centered around accurate values, the distributions for the hyper-priors are very off (almost 0 for both parameters). I am adding reproducible code for what I am doing here

# Libraries
library(R2jags)
library(dplyr)
library(matrixStats)
library(MCMCvis)

#  R2JAGS modeling with simulated data-  using Exponential decay fits
# Creating simulated data

set.seed(1)

N = 10 # No of Subjects
M = 300 # No of Sessions

session_list= c(1:M)

Y_sim_e =matrix(data = NA, nrow = N, ncol = M)
parameters_sim_e = matrix(data =NA, nrow =N,ncol =2)

real_D_mu= 25
real_D_inv_var= 0.3
real_A_mu= 0.8
real_A_inv_var= 30

# For a time-step of one session
for(sub in 1:N){
  D_sub = rnorm(1,real_D_mu,1/real_D_inv_var) # initial starting point w norm distr
  A_sub = rnorm(1,real_A_mu,1/real_A_inv_var) # retention rate w normal dist
  Y_sim_e[sub,1]=  D_sub  
  
  for(session in 2:M){
    Y_sim_e[sub,session]=  A_sub*Y_sim_e[sub,session-1]  
    # + 0.01*rnorm(1)
  }
  parameters_sim_e[sub,]= c(D_sub,A_sub)
}

df_ysim_e= data.frame(t(Y_sim_e)); names(df_ysim_e)= paste("Subject",c(1:N)); df_ysim_e$Session = c(1:M)

#----------------------------------------------------------------------------------
# Creating model
model_random_e <- function(){
  
  # Likelihood
  for(i in 1:N){ # iterating through subjects
    
    for(j in 1:M) {
      Y_sim_e[i,j] ~dnorm(mu[i,j],inv.var[i])
    }
    for(j in 2:M) {
      mu[i,j] <- A[i]* mu[i,j-1]
    }
    
    A[i] ~ dnorm(mu.A[i],inv.var.A[i])
    mu[i,1] ~ dnorm(mu.D[i],inv.var.D[i])    
    
    # Priors
    log(mu.A[i])      <- mu.A_prior[i]  # Ensuring A is positive
    mu.A_prior[i]      ~ dnorm(mu.A_mean_hyperprior, 1)
    log(mu.D[i])       <-  mu.D_prior[i]  # Ensuring D is positive
    mu.D_prior[i]      ~ dnorm(mu.D_mean_hyperprior, 0.1)   
    
    inv.var[i]        ~ dgamma(0.1, 0.1) 
    inv.var.A[i]       ~ dgamma(0.1, 0.1) 
    inv.var.D[i]       ~ dgamma(0.1, 0.1)
  }
  
  # Hyper Priors
  log(mu.A_mean_hyperprior) <- mu.A_mean_hyperprior_pos
  mu.A_mean_hyperprior_pos ~ dnorm(1,1)
  log(mu.D_mean_hyperprior) <- mu.D_mean_hyperprior_pos
  mu.D_mean_hyperprior_pos ~ dnorm(2,1)
  
  # mu.A_mean_hyperprior ~ dnorm(0,1)
  # mu.D_mean_hyperprior ~ dnorm(2,1)
  
}

# specifies initial parameter values for the MCMC sampler
init_values <- function(){
  list()
}
# chooses the parameters whose posterior distributions will be reported  

params <- c("mu",
            "inv.var.A","inv.var.D", 
            "mu.A_prior", "mu.D_prior",
            "mu.A_mean_hyperprior", "mu.D_mean_hyperprior",
            "mu.A", "mu.D")

model_sim_e <- list(Y_sim_e=Y_sim_e,N=N,M=M ,session_list=session_list)

fit_exp <- jags(data = model_sim_e, inits = init_values, parameters.to.save = params, 
                model.file = model_random_e, n.chains = 3, n.iter = 8000, n.burnin = 4000, n.thin = 10, 
                DIC = F) #logical; if TRUE (default), compute deviance, pD, and DIC.

Any help will be greatly appreciated. Also, if there are more appropriate places to ask this question other than https://sourceforge.net or https://stats.stackexchange.com/ that would also be helpful

2

Welcome to discourse! I don’t have the energy to remember enough about JAGS to help you directly, and I think that help with JAGS code is probably off-topic here anyway. However, if you are open to using Stan and can write down the model mathematically (expanding your simulation code to simulate all the way down from your desired hyperpriors would be sufficient), this forum can certainly help you get it implemented.

1 Like