Point estimation and Bayes factors for discrete Markov random fields

Fits a discrete Markov random field (Ising or ordinal) through the pseudo-likelihood. The estimates maximize the pseudo-likelihood, or the pseudo-posterior with with_prior = TRUE, and are obtained with a trust region algorithm; the standard errors come from the Godambe-Huber-White (GHW, or sandwich) estimator. Optionally, Savage-Dickey Bayes factors are computed for each pairwise interaction from the coordinate-rescaled pseudo-posterior (Arena and Marsman, 2026).

dmrfit(
  data,
  parinit = NULL,
  structure = NULL,
  with_prior = FALSE,
  savage_dickey = FALSE,
  M = 1000,
  oversampling = 10,
  ncores = 1,
  thresholds_alpha = 0.5,
  thresholds_beta = 0.5,
  interactions_location = 0,
  interactions_scale = 2.5,
  proposal_df = 5,
  lrt_intervals = FALSE,
  level = 0.95,
  seed = 123
)

Arguments

data: 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).
parinit: initial parameter values for the optimization, 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.
structure: network structure, a P x P symmetric matrix with 1 for an edge and 0 for no edge. If NULL (default), the network is fully connected.
with_prior: logical, whether to include the prior in the optimization: a Beta-Prime prior on the thresholds and a Cauchy prior on the pairwise interactions. Default is FALSE.
savage_dickey: logical, whether to compute the Savage-Dickey density ratio Bayes factor for each pairwise interaction via Bayesian Sampling Importance Resampling (BSIR) from the coordinate-rescaled pseudo-posterior: the pseudo-posterior rescaled around its mode to the Godambe-Huber-White (GHW) covariance, which corrects the underestimated posterior variability of the pseudo-likelihood. The proposal is a multivariate t with the GHW covariance. When zero lies beyond all resampled draws of an interaction, its density at zero cannot be estimated from the draws: the Bayes factor is then reported at the floor 1/(10 * M) and flagged in savage_dickey$zero_beyond_draws. Only available when with_prior = TRUE. Default is FALSE.
M: number of draws resampled in the BSIR step (default is 1000). Only used when savage_dickey = TRUE.
oversampling: multiplier for the number of proposal draws in the BSIR step (default is 10). The total number of proposal draws is M * oversampling. Only used when savage_dickey = TRUE.
ncores: number of cores used to compute the gradient, Hessian and pseudo-likelihood in parallel. 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.
thresholds_alpha: alpha parameter for the Beta-Prime prior on thresholds (default is 0.5). Only used when with_prior = TRUE.
thresholds_beta: beta parameter for the Beta-Prime prior on thresholds (default is 0.5). Only used when with_prior = TRUE.
interactions_location: location parameter for the Cauchy prior on pairwise interactions (default is 0.0). Only used when with_prior = TRUE.
interactions_scale: scale parameter for the Cauchy prior on pairwise interactions (default is 2.5). Only used when with_prior = TRUE.
proposal_df: degrees of freedom for the multivariate t proposal distribution in the BSIR step (default is 5). Only used when savage_dickey = TRUE.
lrt_intervals: whether to compute likelihood-ratio intervals (default is FALSE; see Details): TRUE for all free parameters, or a vector of parameter names (as in names(fit$argument)) for some of them. Each interval needs a series of constrained optimizations, which take a while for large networks. They are returned by confint(fit, method = "lrt"). When they are not computed and there are fewer than 10 observations per free parameter, a warning suggests them.
level: confidence level of the intervals (default is 0.95): the Wald intervals printed by summary() and the likelihood-ratio intervals, so that both are at the same level. It is also the default level of confint().
seed: random seed for reproducibility of the BSIR step (default is 123). The caller's random number stream is restored on exit. Only used when savage_dickey = TRUE.

Value

an object of class dmrfit, a list including the estimates (argument), the value, gradient, Hessian and GHW covariance of the objective at the estimates (utils), the matched call (call), the data dimensions (P, N, n_categories), and, with savage_dickey = TRUE, the Savage-Dickey Bayes factors with the importance-sampling effective sample size (savage_dickey).

Details

The likelihood-ratio interval of a parameter $\theta_k$ is based on its profile pseudo-likelihood: for each value of $\theta_k$, the other parameters are set to the values that maximize the pseudo-likelihood (or the pseudo-posterior with with_prior = TRUE). The interval collects the values where $$C_k \cdot 2 [\log PL(\hat\theta) - \max_{\theta_{-k}} \log PL(\theta_k, \theta_{-k})] \le \chi^2_{1, level}.$$ For a full likelihood the statistic would be $\chi^2_1$-distributed ($C_k = 1$); for the pseudo-likelihood it is $\lambda_k \chi^2_1$-distributed, with $\lambda_k = \Sigma_{kk} / (H^{-1})_{kk}$, because the curvature $H$ of the pseudo-likelihood overstates the information in the data ($\Sigma$ is the GHW covariance). The scaling $C_k = 1 / \lambda_k$ calibrates the statistic (Geys, Molenberghs and Ryan, 1999; Pace, Salvan and Sartori, 2011; Varin, Reid and Firth, 2011). Unlike the Wald intervals, the likelihood-ratio intervals can be asymmetric around the estimate. The bounds are found by uniroot on each side of the estimate. The asymmetry matters when the sample size is small relative to the number of parameters (Arena and Marsman, 2026), so with fewer than 10 observations per free parameter dmrfit() warns that the Wald intervals may be inaccurate.

References

Arena, G. and Marsman, M. (2026). Bayesian inference for discrete Markov random fields through coordinate rescaling. arXiv preprint. doi:10.48550/arXiv.2601.17205

Geys, H., Molenberghs, G., and Ryan, L. M. (1999). Pseudolikelihood modeling of multivariate outcomes in developmental toxicology. Journal of the American Statistical Association, 94(447), 734-745. doi:10.1080/01621459.1999.10474176

Pace, L., Salvan, A., and Sartori, N. (2011). Adjusting composite likelihood ratio statistics. Statistica Sinica, 21(1), 129-148.

Skare, Ø., Bølviken, E., and Holden, L. (2003). Improved sampling-importance resampling and reduced bias importance sampling. Scandinavian Journal of Statistics, 30(4), 719-737.

Varin, C., Reid, N., and Firth, D. (2011). An overview of composite likelihood methods. Statistica Sinica, 21(1), 5-42.

Examples