extract_descriptives_mbp_jags: Extraction of descriptive statistics from a JAGS object for...

In hunansona/ed4bhm: Empirical Determinacy for Bayesian Hierarchical Models

Extraction of descriptive statistics from a JAGS object for base and weighted models

Description

This function takes a JAGS object as input and extracts the descriptive statistics (mean and standard deviation) of all the model parameters (without random effects) for the base and weighted models.

Usage

extract_descriptives_mbp_jags(jags.object.base, delta, prec.name)

Arguments

jags.object.base	character string, the name of the JAGS object fitted in rjags with the base model (non-weighted likelihood, w = 1)
delta	numeric, numerical differentiation step, the weighting factor w = 1 \pm δ, the default value is 0.01.
prec.name	character string, the name of the precision parameter to be used for logarithmic transformation of the posterior sampples for the precision parameter.

Details

The numerical value represents δ, the matrix has six columns with names ("mean.minus", "sd.minus", "mean.base", "sd.base", "mean.plus", "sd.plus") for the mean and standard deviation of the marginal posterior distributions from the models with w = 1 - 0.01, w = 1 (base) and w = 1 + 0.01, and has as many rows as there are parameters in the model (without random effects).

Value

list composed of a numerical matix and a numerical value

See Also

extract_descriptives_jags

Examples

data(eight_schools)
#prior settings
mean_mu<-0
prec_mu<-1/(4^2) 
prec_tau<-1/(5^2)

##########################################################
library(rjags)
library(coda)
#sessionInfo()
list.modules()
# In order to sample from the likelihood in JAGS one has to load the "dic" module and 
#mention "deviance" in variable.names
load.module("dic")
list.modules()

set.seed(44566)

###############
# rjags interface with code in a model string
###############

# In JAGS the normal distribution is parametrized by its mean and precision
jags_data<-list(y=eight_schools$y, sigma_2=eight_schools$prec, 
                mean_mu=mean_mu, prec_mu=prec_mu, 
                prec_tau=prec_tau)

# define parameters
schools.params.jags<-c("mu", "tau_2", "deviance")

schools_centered_modelString <- "
# likelihood
model{
for(i in 1:length(y)){
y[i] ~ dnorm(theta[i], sigma_2[i])
theta[i] ~ dnorm(mu, tau_2);
}

# priors
mu ~ dnorm(mean_mu, prec_mu);
tau ~ dnorm(0,prec_tau)T(0,);
tau_2 <- 1/(tau^2);
}
"

writeLines(schools_centered_modelString, con="TempModel.txt") # write to a file

# # model initiation
rjags.8schools  <- jags.model(
  file = "TempModel.txt",
  data = jags_data,
  n.chains = 4,
  n.adapt = 4000,
  inits = list( list(mu = 3.1, tau = 3.117181, 
                .RNG.name = "base::Wichmann-Hill", .RNG.seed = 314159),
                list(mu = 3.3, tau = 2.517181, 
                .RNG.name = "base::Marsaglia-Multicarry", .RNG.seed = 159314),
                list(mu = 3.6, tau = 3.317181, 
                .RNG.name = "base::Super-Duper", .RNG.seed = 413159),
                list(mu = 3.8, tau = 2.317181,  
                .RNG.name = "base::Mersenne-Twister", .RNG.seed = 143915)))

# burn-in
update(rjags.8schools, n.iter = 20000)
N.iter <- 2*(10^5)
n.thin <- 100
# sampling/monitoring
fit.rjags.8schools <- coda.samples(
  model = rjags.8schools ,
  variable.names = schools.params.jags,
  n.iter = N.iter,
  thin = n.thin)

del = 0.01

descriptives_jags_8schools <- extract_descriptives_mbp_jags(jags.obj = fit.rjags.8schools, 
                                                          prec.name = "tau_2", delta = del)

hunansona/ed4bhm documentation built on June 15, 2022, 6:42 p.m.