tests/regression_tests/hojsgaard_model_tests/triangle_model/triangle-model_bayes-loglin_Notes-example_JAGS.R

In npetraco/CRFutil: Handy Utility Functions for Use with CRF and gRbase packages

library(CRFutil)
library(CRFutilRcppComponents)
library(R2jags)

# Triangle model
grphf <- ~1:2+1:3+2:3
adj <- ug(grphf, result="matrix")
gp <- ug(grphf, result = "graph")
dev.off()
iplot(gp)

f0 <- function(y){ as.numeric(c((y==1),(y==2)))}
n.states <- 2

# Instantiate model and pick a "true" theta:
knm <- make.empty.field(adj.mat = adj, parameterization.typ = "standard", plotQ = T)
set.seed(6)
knm$par <- runif(6,-1.5,1.1)
knm$par
# "true" theta:
#0.07629757  0.93786913 -0.81268462 -0.51175581  0.59945681  1.04299689
out.pot <- make.pots(parms = knm$par,  crf = knm,  rescaleQ = T, replaceQ = T)
# So now sample from the model as if we obtained an experimental sample:
num.samps <- 25
set.seed(1)
samps <- sample.exact(knm, num.samps)
mrf.sample.plot(samps)

# Make state count freq table from samples and compute frequencies of all possible state configs.
# State configs not observed will have 0 freq
ftab  <- data.frame(ftable(data.frame(samps)))
X.all <- ftab[,1:ncol(samps)]
X.all <- sapply(1:ncol(X.all), function(xx){as.numeric(X.all[,xx])}) # Needed for C routines
freqs <- ftab[,ncol(ftab)]

# Model Matrix with respect to graph
theta.pars     <- fix_node_and_edge_par(node_par = knm$node.par, edge_par = knm$edge.par) # Needed for C routines
knm$node.par.c <- theta.pars$node_par
knm$edge.par.c <- theta.pars$edge_par
M              <- compute_model_matrix(
  configs       = X.all,
  edge_mat      = knm$edges,
  node_par      = knm$node.par.c,
  edge_par      = knm$edge.par.c,
  num_params_in = knm$n.par)

# Check C vs R model matrix routines:
X.all2 <- ftab[,1:ncol(samps)]
M2     <- compute.model.matrix(configs=X.all2, edges.mat=knm$edges, node.par = knm$node.par, edge.par = knm$edge.par, ff = f0)
M-M2

dat <- list(
  p = ncol(M),
  N = nrow(M),
  y = freqs,
  Mmodl = M
)
dat

fpth <- "inst/poisson_model.bug"
model_parameters <- c("alpha","theta")

# Run the model
model_run = jags(data = dat, #inits = init,
                 parameters.to.save = model_parameters,
                 model.file = fpth,
                 n.chains = 4,
                 n.iter = 50000,
                 n.burnin = 5000,
                 n.thin = 50)

# Check the output - are the true values inside the 95% CI?
# Also look at the R-hat values - they need to be close to 1 if convergence has been achieved
plot(model_run)
print(model_run)
traceplot(model_run)

# Prepare a fit to obtain joint dist. This time we will use JAGS instead of Stan
fit <- make.empty.field(adj.mat = adj, parameterization.typ = "standard", plotQ = T)

# Here we remove the extra index put in by CRF using the C function. We need to do this for the the RcppComponent functions
#theta.pars   <- fix_node_and_edge_par(node_par = fit$node.par, edge_par = fit$edge.par)
#fit$edge.par # Should we replace with what's in theta.pars??

post    <- print(model_run)
fit$par <- post$mean$theta

out.pot2 <- make.pots(parms = fit$par,  crf = fit,  rescaleQ = T, replaceQ = T)
fit$node.pot
fit$edge.pot

# Configuration probabilities:
pot.info <- make.gRbase.potentials(fit, node.names = gp@nodes, state.nmes = c("1","2"))
pot.info

gR.dist.info    <- distribution.from.potentials(pot.info$node.potentials, pot.info$edge.potentials)
logZ            <- gR.dist.info$logZ
joint.dist.info <- as.data.frame(as.table(gR.dist.info$state.probs))
joint.dist.info

# True model:
pot.info.true.model <- make.gRbase.potentials(knm, node.names = gp@nodes, state.nmes = c("1","2"))
gR.dist.info.true.model    <- distribution.from.potentials(pot.info.true.model$node.potentials, pot.info.true.model$edge.potentials)
logZ.true.model            <- gR.dist.info.true.model$logZ
joint.dist.info.true.model <- as.data.frame(as.table(gR.dist.info.true.model$state.probs))
joint.dist.info.true.model

# Compare:
fit.cp  <- round(joint.dist.info[,4]*100, 1)     # Bayes logistic from JAGS
true.model.cp  <- round(joint.dist.info.true.model[,4]*100, 1) # True
cbind(joint.dist.info[,c(2,3,1)], fit.cp, true.model.cp)

# X1 X2 X3 fit.cp true.model.cp
# 1  1  1  1   10.2          18.4
# 2  1  1  2    1.1           8.0
# 3  2  1  1   22.2          15.6
# 4  2  1  2   23.4          22.6
# 5  1  2  1    8.9           4.2
# 6  1  2  2   10.1          14.9
# 7  2  2  1    2.0           1.3
# 8  2  2  2   22.1          15.0

# Double check
cbind(
  joint.dist.info.true.model,
  joint.dist.info
)