Bayesian inference for a discrete Markov random field (Ising or ordinal) by Markov chain Monte Carlo. By default, the posterior is approximated by the coordinate-rescaled (CoRe) pseudo-posterior of Arena and Marsman (2026): the pseudo-posterior rescaled around its mode so that its covariance matches the Godambe-Huber-White (GHW) covariance, which corrects the underestimated posterior variability of the pseudo-likelihood. All samplers use the Fisher adaptive Metropolis-adjusted Langevin algorithm (FisherMALA; Titsias, 2024).
dmrfit_bayes(
data,
parinit = NULL,
method = c("core", "adacore", "exact", "dmh"),
scale = c("ghw", "mch", "rm"),
nsim = 1000,
burnin = 1000,
ncores = 1,
thresholds_alpha = 0.5,
thresholds_beta = 0.5,
interactions_location = 0,
interactions_scale = 2.5,
sigma2 = 0.1,
seed = 123,
progress = TRUE,
control = list()
)data matrix, with rows as samples and columns as variables. Each variable should be rescaled to the range of 0 to m-1, where m is the number of categories for that variable. The baseline category is always the minimum value in the variable. The internal processing will check if the variables are rescaled and will rescale them if necessary. If there are any NAs in the data, they will be removed before optimization (listwise deletion).
initial parameter values for the optimization that finds the pseudo-posterior mode, a vector of length equal to the number of parameters in the model, sum(n_categories - 1) + P * (P - 1) / 2, where P is the number of variables. If NULL (default), a vector of zeros.
the sampler: "core" (default) samples the CoRe pseudo-posterior with a fixed rescaling
matrix; "adacore" adapts the rescaling matrix during burn-in at the running mean of the draws;
"exact" samples the posterior based on the full likelihood, whose normalizing constant is computed by
enumerating all response patterns (only for small networks, see control); "dmh" samples the
posterior based on the full likelihood with the double Metropolis-Hastings algorithm (Liang, 2010). DMH simulates auxiliary data at every
iteration and is much slower than the other methods (minutes to hours instead of seconds, depending on the
network size and control$dmh_aux); "core" runs at the cost of the pseudo-posterior sampler.
the covariance to which the pseudo-posterior is rescaled, only used when method = "core":
"ghw" (default) the Godambe-Huber-White (GHW) covariance; "mch" the inverse of the negative
Hessian of the full log-posterior at the pseudo-posterior mode, estimated by Monte Carlo simulation; "rm"
the same Hessian at the full-posterior mode, estimated by a Newton-type Robbins-Monro algorithm (the sampler
remains centered at the pseudo-posterior mode). "mch" and "rm" simulate data from the model and are
slower.
number of posterior draws kept after burn-in. Default is 1000.
number of burn-in iterations, after an initial adaptive stage of control$adaptive_stage iterations. Default is 1000.
number of cores used to compute the gradient, Hessian and pseudo-likelihood in parallel when finding the pseudo-posterior mode. It is capped at the number of available cores minus one, and has no effect when the package was built without OpenMP support. Default is 1.
alpha parameter for the Beta-Prime prior on thresholds (default is 0.5).
beta parameter for the Beta-Prime prior on thresholds (default is 0.5).
location parameter for the Cauchy prior on pairwise interactions (default is 0.0).
scale parameter for the Cauchy prior on pairwise interactions (default is 2.5).
initial value for the adaptive variance parameter in the FisherMALA sampler (default is 0.1).
random seed for reproducibility of the MCMC sampler (default is 123). The caller's random number stream is restored on exit.
logical, whether to show a progress bar while sampling (default is TRUE).
a list of tuning settings, each with a default: adaptive_stage (500) iterations of the
initial adaptive stage; max_states (1e6) largest number of response patterns for method = "exact";
dmh_aux (25000) and dmh_gibbs_iter (5) auxiliary draws and Gibbs sweeps per iteration of
method = "dmh"; mc_draws (1e5) and mc_gibbs_iter (5) simulated observations and Gibbs sweeps
for scale = "mch" and "rm"; rm_iter (50) Robbins-Monro iterations for scale = "rm".
an object of class dmrfit_bayes (also of class dmrfit), including the pseudo-posterior mode (argument), the posterior draws (draws, a matrix with parameters by rows and iterations by columns), the acceptance rate (acceptance), the sampling method and scale (method, scale), the Savage-Dickey Bayes factors with the effective sample size of each parameter (savage_dickey), and the multivariate effective sample size of the draws (mess; Vats, Flegal and Jones, 2019), which is NA when there are fewer than P + 1 batches of floor(sqrt(nsim)) draws per parameter.
Arena, G. and Marsman, M. (2026). Bayesian inference for discrete Markov random fields through coordinate rescaling. arXiv preprint. doi:10.48550/arXiv.2601.17205
Liang, F. (2010). A double Metropolis-Hastings sampler for spatial models with intractable normalizing constants. Journal of Statistical Computation and Simulation, 80(9), 1007-1022.
Titsias, M. K. (2024). Optimal preconditioning and Fisher adaptive Langevin sampling. In Proceedings of the 37th International Conference on Neural Information Processing Systems (NIPS '23). Red Hook, NY, USA: Curran Associates.
Vats, D., Flegal, J. M., and Jones, G. L. (2019). Multivariate output analysis for Markov chain Monte Carlo. Biometrika, 106(2), 321-337.
dmrfit for point estimation and Bayes factors by sampling importance resampling.
# binary responses of 1000 observations on 3 nodes (independent, for illustration)
set.seed(123)
n <- 1000
P <- 3
prob <- c(0.2, 0.5, 0.3)[rep(1:P, each = n)]
data <- matrix(rbinom(n * P, size = 1, prob = prob), nrow = n, ncol = P)
fit <- dmrfit_bayes(data = data, nsim = 500, burnin = 500, progress = FALSE)
fit
#>
#> Call:
#> dmrfit_bayes(data = data, nsim = 500, burnin = 500, progress = FALSE)
#>
#> Thresholds (posterior mode):
#> mu[1,1] mu[2,1] mu[3,1]
#> -1.4437 0.0164 -0.8830
#>
#> Pairwise interactions (posterior mode):
#> theta[2,1] theta[3,1] theta[3,2]
#> -0.1373 0.3452 -0.0040
#>
#> Posterior: coordinate-rescaled pseudo-posterior (Godambe-Huber-White covariance)
#> MCMC samples: 500 Acceptance rate: 0.6 Elapsed: 0 sec
#> Multivariate effective sample size: 236
#>
summary(fit)
#>
#> Call:
#> dmrfit_bayes(data = data, nsim = 500, burnin = 500, progress = FALSE)
#>
#> Discrete Markov Random Field (Ising)
#> Estimation method: Bayesian (FisherMALA), coordinate-rescaled pseudo-posterior (Godambe-Huber-White covariance)
#> Nodes: 3 Observations: 1000 Parameters: 6
#> MCMC samples: 500 Acceptance rate: 0.6 Elapsed: 0 sec
#> Multivariate effective sample size: 236
#> ------------------------------------------------------------
#>
#> Thresholds:
#> Post.Mode Post.Mean Post.SD 2.5% 97.5%
#> mu[1,1] -1.4437 -1.4451 0.1262 -1.6943 -1.2143
#> mu[2,1] 0.0164 0.0099 0.0798 -0.1300 0.1659
#> mu[3,1] -0.8830 -0.8824 0.1072 -1.0950 -0.6627
#>
#> Pairwise interactions:
#> Post.Mode Post.Mean Post.SD 2.5% 97.5%
#> theta[2,1] -0.1373 -0.1309 0.1577 -0.4166 0.1759
#> theta[3,1] 0.3452 0.3467 0.1651 0.0266 0.7016
#> theta[3,2] -0.0040 -0.0033 0.1405 -0.2926 0.2595
#>
#> Negative log pseudo-likelihood at posterior mode: 1814.438
#>
#> Savage-Dickey density ratio [prior: Cauchy( 0 , 2.5 )]
#> H0: theta = 0 for each pairwise interaction
#>
#> Post.Mode Post.SD BF_01 Pr(=0|x)
#> theta[2,1] -0.1373 0.1577 14.13 0.9339
#> theta[3,1] 0.3452 0.1651 2.469 0.7118
#> theta[3,2] -0.0040 0.1405 22.69 0.9578
#>
#> MCMC samples: 500 Effective sample size of the interactions: min 181.9 , median 201
# small network: the exact posterior is also available
fit_exact <- dmrfit_bayes(data = data, method = "exact", nsim = 500, burnin = 500, progress = FALSE)
summary(fit_exact)
#>
#> Call:
#> dmrfit_bayes(data = data, method = "exact", nsim = 500, burnin = 500,
#> progress = FALSE)
#>
#> Discrete Markov Random Field (Ising)
#> Estimation method: Bayesian (FisherMALA), full-likelihood posterior (exact normalizing constant)
#> Nodes: 3 Observations: 1000 Parameters: 6
#> MCMC samples: 500 Acceptance rate: 0.592 Elapsed: 0 sec
#> Multivariate effective sample size: 225.9
#> ------------------------------------------------------------
#>
#> Thresholds:
#> Post.Mode Post.Mean Post.SD 2.5% 97.5%
#> mu[1,1] -1.4437 -1.4425 0.1273 -1.7090 -1.1983
#> mu[2,1] 0.0164 0.0105 0.0792 -0.1427 0.1588
#> mu[3,1] -0.8830 -0.8810 0.1072 -1.0900 -0.6708
#>
#> Pairwise interactions:
#> Post.Mode Post.Mean Post.SD 2.5% 97.5%
#> theta[2,1] -0.1373 -0.1353 0.1608 -0.4456 0.1609
#> theta[3,1] 0.3452 0.3366 0.1589 0.0230 0.6419
#> theta[3,2] -0.0040 -0.0015 0.1370 -0.2828 0.2825
#>
#> Negative log pseudo-likelihood at posterior mode: 1814.438
#>
#> Savage-Dickey density ratio [prior: Cauchy( 0 , 2.5 )]
#> H0: theta = 0 for each pairwise interaction
#>
#> Post.Mode Post.SD BF_01 Pr(=0|x)
#> theta[2,1] -0.1373 0.1608 15.11 0.9379
#> theta[3,1] 0.3452 0.1589 2.506 0.7148
#> theta[3,2] -0.0040 0.1370 23.54 0.9592
#>
#> MCMC samples: 500 Effective sample size of the interactions: min 154.2 , median 178.8