Framework for simulation study to validate bayesian models

Question

I am looking for a framework that would allow to take JAGS/bugs model and on many sets of simulated data test if there is a bias (or not) in the parameter estimates (the real parameters would be known because the data are simulated). Is there any such framework? I've seen only run.jags.study() in R package runjags, but this is focused on crossvalidation. I need to check the bias in parameter estimates.

Do you know such framework? If not, I am especially interested how to evaluate the bias across the generated data sets. How would you do that?

Thanks in advance!

Answer 1

I have been developing the framework you are referring to. See

https://github.com/audrey-b/simanalyse

The R package has separate functions to simulate the datasets (using JAGS or R code), analyze (using JAGS) and evaluate the performance of the model. It can easily save all the results to files (thus saving on RAM) and parallelize on a local or cluster computer. Setting the seed is also streamlined, simply using set.seed(). The MCMC chains are automatically ran until convergence, while being thinned and burned-in to a specified size, so you do not need to guess a thinning factor.

There are a lot of evaluation measures to choose from, including bias. The currently implemented measures are: bias, relative bias, bias ratio, coverage probability, expected interval length, expected posterior variance, expected posterior standard deviation, variance, standard error, mean square error, root mean square error, relative root mean square error, coefficient of variation and power.

Here is a simple example (code subject to change since it's not on CRAN yet):

remotes::install_github("audrey-b/simanalyse")
remotes::install_github("poissonconsulting/sims")

library(simanalyse)
library(sims)

set.seed(123)

# Set up the model

code <- "for(i in 1:10){  
          y[i] ~ dnorm(mu, 1/sigma^2)}"

params <- list(sigma = 2)

constants <- list(mu = 0)

# Simulate 100 datasets

data <- sims::sims_simulate(code, 
                            parameters = params, 
                            constants = constants,
                            nsims = 100)

# Analyse 100 datasets

prior <- "sigma ~ dunif(0, 6)"

analyses <- sma_analyse(data,
                        code = code,
                        code.add = prior,
                        mode = sma_set_mode("report"))

# Evaluate the bias

sma_evaluate(analyses, 
             measures = "bias",
             parameters = params)

 term      bias
sigma 0.2695341

Answer 2

Cross-validation is only one potential use of the run.jags.study function - it can also be used for estimating model performance in terms of bias and/or coverage of 95% confidence intervals. For example:

library('runjags')

# Target values:
mu_target <- 0.3
tau_target <- 0.2

N <- 10

# Function to simulate data:
datafun <- function(simulation_number){
    y <- rnorm(N, mean=mu_target, sd=1/tau_target^0.5)
    return(list(y=y))
}

# Simple model:
model <- '
model{

    for(i in 1:N){
        y[i] ~ dnorm(mu, tau)
    }
    #data# N

    mu ~ dnorm(0, 10^-6)
    tau ~ dgamma(0.01, 0.01)
}'

# Get bias and coverage:
results <- run.jags.study(240, model, datafun, targets=list(mu=mu_target, tau=tau_target), n.cores=6)
results

Average values obtained from a JAGS study with a total of 240 simulations:

    Target Av.Median Av.Mean Av.Lower95%CI Av.Upper95%CI Av.Range95%CI Prop.Within95%CI Av.AutoCorr(Lag10) Simulations
mu     0.3   0.30215 0.30247       -1.2433        1.8309        3.0742           0.9375          -0.000946         240
tau    0.2   0.23197 0.25005      0.053889       0.48251       0.42862          0.95833       -0.000085487         240

Average time taken:  0.2 seconds (range: 0.2 seconds - 0.3 seconds)
Average adapt+burnin required:  5000 (range: 5000 - 5000)
Average samples required:  10000 (range: 10000 - 10000)

So there is a small positive bias in mean and a slightly larger positive bias for tau, but both sets of 95% confidence intervals contain the true value approximately 95% of the time.

See also section 3 of the following for a more complete/useful example:

Denwood, M.J. 2016. runjags: An R Package Providing Interface Utilities, Model Templates, Parallel Computing Methods and Additional Distributions for MCMC Models in JAGS. J. Stat. Softw. 71. doi:10.18637/jss.v071.i09.

Matt