Analyzing a psychological network with dmrfit

This vignette illustrates a complete network analysis with dmrfit. It covers point estimation with confidence intervals, the evidence for each edge, Bayesian estimation, centrality, and the choice between the sampling methods. Each section first describes how to run the analysis and interpret its output. The paragraphs titled Details explain the methods and can be skipped on a first reading. The notation follows Arena and Marsman (2026).

Discrete Markov random fields and the pseudo-likelihood

A discrete Markov random field models the conditional dependencies between discrete random variables, such as the items of a questionnaire, as a network: two variables are connected when they are dependent given all the other variables. Let $\mathbf{X} = (\mathbf{X}_1, \ldots, \mathbf{X}_n)$ denote the observed data for $n$ observations on $p$ ordinal variables, where each variable takes one of the categories $0, 1, \ldots, m$, and let $\mathbf{X}_{\nu} = (X_{\nu 1}, \ldots, X_{\nu p})$ denote the $\nu$-th observation. The ordinal Markov random field (OMRF) is defined by the parameter vector $\boldsymbol{\eta} = (\boldsymbol{\mu}, \boldsymbol{\theta})$, which comprises the threshold parameters $\mu_{ih}$, the tendency of variable $i$ to take the non-baseline category $h = 1, \ldots, m$ that is not explained by its interactions with the other variables, and the interaction parameters $\theta_{ij}$, the pairwise conditional dependency between variables $i$ and $j$. An interaction $\theta_{ij} = 0$ means that variables $i$ and $j$ are conditionally independent, that is, there is no edge between them. When the variables are binary, the OMRF reduces to the Ising model.

The full likelihood of the OMRF contains a normalizing constant that sums over all $(m + 1)^p$ possible response patterns, which is infeasible beyond a few variables. dmrfit uses the pseudo-likelihood instead, the product of the conditional distributions of each variable given the remaining variables, \[ \tilde{f}(\mathbf{X}; \boldsymbol{\eta}) = \prod_{\nu=1}^{n} \prod_{i=1}^{p} f\left(X_{\nu i} \mid \mathbf{X}_{\nu(-i)}, \boldsymbol{\eta}\right), \qquad f\left(X_{\nu i} = h \mid \mathbf{X}_{\nu(-i)}, \boldsymbol{\eta}\right) = \frac{\exp\left(\mu_{ih} + h \sum_{j \neq i} \theta_{ij} X_{\nu j}\right)} {1 + \sum_{k=1}^{m} \exp\left(\mu_{ik} + k \sum_{j \neq i} \theta_{ij} X_{\nu j}\right)}, \] with $\mu_{i0} = 0$, whose normalizing constants are computationally tractable. Combined with a prior $\pi(\boldsymbol{\eta})$, it defines the pseudo-posterior $\tilde{\pi}(\boldsymbol{\eta} \mid \mathbf{X}) \propto \tilde{f}(\mathbf{X}; \boldsymbol{\eta})\, \pi(\boldsymbol{\eta})$. The pseudo-likelihood estimates are consistent, but the curvature of the pseudo-likelihood overstates the information in the data: used naively, it underestimates the variability of the estimates and the posterior variance. The package corrects this underestimation in two ways.

dmrfit() estimates the parameters by maximizing the pseudo-likelihood (or the pseudo-posterior), with robust standard errors, confidence intervals and Bayes factors for the edges.
dmrfit_bayes() samples from the pseudo-posterior rescaled to the scale of the full posterior (coordinate rescaling, CoRe, Arena and Marsman, 2026).

library(dmrfit)

The data

rads2 contains the responses of 917 adolescents from Lima (Peru) to the 25 items of the four-factor version of the Reynolds Adolescent Depression Scale (RADS-2), each answered on a four-point scale (Ramos-Vera et al., 2023). The positively worded items were reverse-coded, so that higher values mean more symptoms for every item. The items belong to four clusters, stored with the data. Here we use the 10 items of two clusters, dysphoria and anhedonia/negative affect, with $p = 10$ variables, $m + 1 = 4$ categories and $n = 917$ observations.

data(rads2)
clusters <- attr(rads2, "clusters")
items <- names(clusters)[clusters %in% c("Dysphoria", "Anhedonia/Negative Affect")]
X <- rads2[, items]
clusters[items]
#>                          D1                          D3 
#> "Anhedonia/Negative Affect"                 "Dysphoria" 
#>                          D5                          D6 
#> "Anhedonia/Negative Affect"                 "Dysphoria" 
#>                          D7                          D8 
#>                 "Dysphoria"                 "Dysphoria" 
#>                         D12                         D16 
#> "Anhedonia/Negative Affect"                 "Dysphoria" 
#>                         D21                         D26 
#>                 "Dysphoria"                 "Dysphoria"

The response frequencies show that the highest category is rarely chosen for some items. This matters below, because a threshold of a rare category is estimated from few observations.

sapply(X, function(item) table(factor(item, levels = 1:4)))
#>    D1  D3  D5  D6  D7  D8 D12 D16 D21 D26
#> 1 500 275 257 417 282 323 384 325 495 177
#> 2 325 231 348 243 265 275 331 324 208 339
#> 3  66 295 200 175 281 219 137 209 151 275
#> 4  26 116 112  82  89 100  65  59  63 126

Point estimation and confidence intervals

dmrfit() estimates the thresholds and interactions. With with_prior = TRUE it maximizes the pseudo-posterior, which gives the maximum a posteriori (MAP) estimate, with beta-prime priors on the exponentiated thresholds, $\exp(\mu_{ih}) \sim \text{Beta-Prime}(0.5, 0.5)$, and Cauchy priors on the interactions, $\theta_{ij} \sim \text{Cauchy}(0, 2.5)$. Without the prior it gives the maximum pseudo-likelihood estimate (MPLE). lrt_intervals = TRUE adds likelihood-ratio intervals for all parameters, and savage_dickey = TRUE the Bayes factors of the next section.

fit <- dmrfit(X, with_prior = TRUE, savage_dickey = TRUE, lrt_intervals = TRUE)

In the output, mu[i,h] is the threshold $\mu_{ih}$ and theta[i,j] the interaction $\theta_{ij}$, with the variables numbered in the column order of the data. summary() prints the estimates with their standard errors, and the Wald and likelihood-ratio intervals next to each other.

summary(fit)
#> 
#> Call:
#> dmrfit(data = X, with_prior = TRUE, savage_dickey = TRUE, lrt_intervals = TRUE)
#> 
#> Discrete Markov Random Field
#> Estimation method: maximum pseudo-likelihood (with prior)
#> Nodes: 10  Observations: 917  Free parameters: 75 
#> Standard errors: Godambe-Huber-White (GHW) estimator
#> ------------------------------------------------------------ 
#> 
#> Thresholds:
#>          Post.Mode   Post.SD  z value  Pr(>|z|)
#> mu[1,1]  -1.820877  0.144418 -12.6084 < 2.2e-16
#> mu[1,2]  -5.576649  0.370759 -15.0412 < 2.2e-16
#> mu[1,3]  -9.452950  0.705495 -13.3990 < 2.2e-16
#> mu[2,1]  -1.169270  0.115919 -10.0870 < 2.2e-16
#> mu[2,2]  -2.417688  0.205029 -11.7920 < 2.2e-16
#> mu[2,3]  -5.524127  0.387496 -14.2560 < 2.2e-16
#> mu[3,1]  -0.209862  0.105177  -1.9953   0.04601
#> mu[3,2]  -1.582957  0.180676  -8.7613 < 2.2e-16
#> mu[3,3]  -3.399964  0.323605 -10.5065 < 2.2e-16
#> mu[4,1]  -1.451465  0.114853 -12.6376 < 2.2e-16
#> mu[4,2]  -3.087075  0.231522 -13.3338 < 2.2e-16
#> mu[4,3]  -5.610251  0.404348 -13.8748 < 2.2e-16
#> mu[5,1]  -1.468367  0.140364 -10.4611 < 2.2e-16
#> mu[5,2]  -3.827827  0.327542 -11.6865 < 2.2e-16
#> mu[5,3]  -8.693257  0.650305 -13.3680 < 2.2e-16
#> mu[6,1]  -1.445533  0.129159 -11.1918 < 2.2e-16
#> mu[6,2]  -3.850868  0.282119 -13.6498 < 2.2e-16
#> mu[6,3]  -7.703633  0.506745 -15.2022 < 2.2e-16
#> mu[7,1]  -1.112222  0.128882  -8.6298 < 2.2e-16
#> mu[7,2]  -3.438122  0.283897 -12.1105 < 2.2e-16
#> mu[7,3]  -6.222909  0.568813 -10.9402 < 2.2e-16
#> mu[8,1]  -0.850789  0.108832  -7.8175 5.389e-15
#> mu[8,2]  -2.513389  0.216437 -11.6125 < 2.2e-16
#> mu[8,3]  -5.457662  0.402758 -13.5507 < 2.2e-16
#> mu[9,1]  -2.282995  0.131021 -17.4247 < 2.2e-16
#> mu[9,2]  -4.633978  0.279444 -16.5829 < 2.2e-16
#> mu[9,3]  -8.242312  0.510441 -16.1474 < 2.2e-16
#> mu[10,1]  0.045836  0.107724   0.4255   0.67048
#> mu[10,2] -1.056172  0.173051  -6.1032 1.039e-09
#> mu[10,3] -3.082789  0.305565 -10.0888 < 2.2e-16
#> 
#> Pairwise interactions:
#>              Post.Mode    Post.SD z value  Pr(>|z|)    
#> theta[2,1]   0.1121620  0.0703940  1.5933 0.1110827    
#> theta[3,1]   0.4062825  0.0694625  5.8490 4.947e-09 ***
#> theta[4,1]   0.0689947  0.0617306  1.1177 0.2637062    
#> theta[5,1]   0.2416338  0.0884372  2.7323 0.0062900 ** 
#> theta[6,1]  -0.0561661  0.0794795 -0.7067 0.4797695    
#> theta[7,1]   0.4131106  0.0758050  5.4496 5.047e-08 ***
#> theta[8,1]   0.1234728  0.0690276  1.7887 0.0736560 .  
#> theta[9,1]   0.2035943  0.0701258  2.9033 0.0036929 ** 
#> theta[10,1] -0.0476660  0.0672855 -0.7084 0.4786879    
#> theta[3,2]   0.1099157  0.0517141  2.1255 0.0335490 *  
#> theta[4,2]   0.2402308  0.0504967  4.7574 1.961e-06 ***
#> theta[5,2]   0.2744462  0.0614688  4.4648 8.014e-06 ***
#> theta[6,2]   0.2313118  0.0564412  4.0983 4.162e-05 ***
#> theta[7,2]   0.1433946  0.0575102  2.4934 0.0126535 *  
#> theta[8,2]   0.0939737  0.0524526  1.7916 0.0731982 .  
#> theta[9,2]   0.1510427  0.0572081  2.6402 0.0082850 ** 
#> theta[10,2]  0.0902188  0.0507052  1.7793 0.0751938 .  
#> theta[4,3]   0.0343960  0.0504221  0.6822 0.4951364    
#> theta[5,3]   0.0113068  0.0603851  0.1872 0.8514689    
#> theta[6,3]  -0.0468243  0.0572661 -0.8177 0.4135508    
#> theta[7,3]   0.3783185  0.0600302  6.3021 2.936e-10 ***
#> theta[8,3]  -0.0029214  0.0531458 -0.0550 0.9561629    
#> theta[9,3]   0.1144986  0.0540580  2.1181 0.0341691 *  
#> theta[10,3] -0.0168471  0.0506416 -0.3327 0.7393810    
#> theta[5,4]   0.1831356  0.0584058  3.1356 0.0017152 ** 
#> theta[6,4]   0.0861031  0.0541438  1.5903 0.1117747    
#> theta[7,4]   0.0657230  0.0555872  1.1823 0.2370703    
#> theta[8,4]   0.0894205  0.0495450  1.8048 0.0711011 .  
#> theta[9,4]   0.2424225  0.0527873  4.5924 4.381e-06 ***
#> theta[10,4]  0.0534099  0.0480940  1.1105 0.2667697    
#> theta[6,5]   0.8163969  0.0866513  9.4216 < 2.2e-16 ***
#> theta[7,5]   0.1360569  0.0740861  1.8365 0.0662881 .  
#> theta[8,5]   0.1658596  0.0679567  2.4407 0.0146602 *  
#> theta[9,5]   0.1407857  0.0779097  1.8070 0.0707566 .  
#> theta[10,5]  0.1813855  0.0623835  2.9076 0.0036423 ** 
#> theta[7,6]  -0.0278154  0.0694366 -0.4006 0.6887246    
#> theta[8,6]   0.1779921  0.0614563  2.8962 0.0037766 ** 
#> theta[9,6]   0.1920124  0.0705658  2.7210 0.0065076 ** 
#> theta[10,6]  0.2100408  0.0593717  3.5377 0.0004036 ***
#> theta[8,7]  -0.0286764  0.0594765 -0.4821 0.6297012    
#> theta[9,7]   0.1635554  0.0672080  2.4336 0.0149508 *  
#> theta[10,7]  0.0152829  0.0590018  0.2590 0.7956162    
#> theta[9,8]   0.2149666  0.0587718  3.6576 0.0002545 ***
#> theta[10,8]  0.1885638  0.0529913  3.5584 0.0003731 ***
#> theta[10,9]  0.1542764  0.0574723  2.6844 0.0072668 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Intervals (95%): Wald (GHW standard errors) and likelihood-ratio
#>             Wald lower Wald upper LRT lower LRT upper
#> mu[1,1]        -2.1039    -1.5378   -2.1084   -1.5420
#> mu[1,2]        -6.3033    -4.8500   -6.3232   -4.8697
#> mu[1,3]       -10.8357    -8.0702  -10.8810   -8.1158
#> mu[2,1]        -1.3965    -0.9421   -1.3976   -0.9430
#> mu[2,2]        -2.8195    -2.0158   -2.8260   -2.0219
#> mu[2,3]        -6.2836    -4.7646   -6.3009   -4.7818
#> mu[3,1]        -0.4160    -0.0037   -0.4156   -0.0031
#> mu[3,2]        -1.9371    -1.2288   -1.9406   -1.2319
#> mu[3,3]        -4.0342    -2.7657   -4.0485   -2.7794
#> mu[4,1]        -1.6766    -1.2264   -1.6787   -1.2282
#> mu[4,2]        -3.5409    -2.6333   -3.5508   -2.6428
#> mu[4,3]        -6.4028    -4.8177   -6.4247   -4.8395
#> mu[5,1]        -1.7435    -1.1933   -1.7456   -1.1951
#> mu[5,2]        -4.4698    -3.1859   -4.4842   -3.1999
#> mu[5,3]        -9.9678    -7.4187  -10.0023   -7.4533
#> mu[6,1]        -1.6987    -1.1924   -1.7008   -1.1943
#> mu[6,2]        -4.4038    -3.2979   -4.4153   -3.3092
#> mu[6,3]        -8.6968    -6.7104   -8.7206   -6.7343
#> mu[7,1]        -1.3648    -0.8596   -1.3668   -0.8614
#> mu[7,2]        -3.9945    -2.8817   -4.0077   -2.8943
#> mu[7,3]        -7.3378    -5.1081   -7.3762   -5.1463
#> mu[8,1]        -1.0641    -0.6375   -1.0649   -0.6381
#> mu[8,2]        -2.9376    -2.0892   -2.9453   -2.0964
#> mu[8,3]        -6.2471    -4.6683   -6.2679   -4.6889
#> mu[9,1]        -2.5398    -2.0262   -2.5436   -2.0298
#> mu[9,2]        -5.1817    -4.0863   -5.1954   -4.0999
#> mu[9,3]        -9.2428    -7.2419   -9.2697   -7.2691
#> mu[10,1]       -0.1653     0.2570   -0.1640    0.2585
#> mu[10,2]       -1.3953    -0.7170   -1.3974   -0.7186
#> mu[10,3]       -3.6817    -2.4839   -3.6938   -2.4955
#> theta[2,1]     -0.0258     0.2501   -0.0250    0.2512
#> theta[3,1]      0.2701     0.5424    0.2724    0.5450
#> theta[4,1]     -0.0520     0.1900   -0.0518    0.1904
#> theta[5,1]      0.0683     0.4150    0.0699    0.4170
#> theta[6,1]     -0.2119     0.0996   -0.2123    0.0996
#> theta[7,1]      0.2645     0.5617    0.2668    0.5644
#> theta[8,1]     -0.0118     0.2588   -0.0116    0.2594
#> theta[9,1]      0.0662     0.3410    0.0671    0.3424
#> theta[10,1]    -0.1795     0.0842   -0.1798    0.0842
#> theta[3,2]      0.0086     0.2113    0.0089    0.2119
#> theta[4,2]      0.1413     0.3392    0.1423    0.3404
#> theta[5,2]      0.1540     0.3949    0.1553    0.3964
#> theta[6,2]      0.1207     0.3419    0.1218    0.3432
#> theta[7,2]      0.0307     0.2561    0.0313    0.2570
#> theta[8,2]     -0.0088     0.1968   -0.0085    0.1973
#> theta[9,2]      0.0389     0.2632    0.0399    0.2644
#> theta[10,2]    -0.0092     0.1896   -0.0088    0.1901
#> theta[4,3]     -0.0644     0.1332   -0.0645    0.1333
#> theta[5,3]     -0.1070     0.1297   -0.1071    0.1298
#> theta[6,3]     -0.1591     0.0654   -0.1595    0.0652
#> theta[7,3]      0.2607     0.4960    0.2627    0.4983
#> theta[8,3]     -0.1071     0.1012   -0.1073    0.1012
#> theta[9,3]      0.0085     0.2205    0.0089    0.2211
#> theta[10,3]    -0.1161     0.0824   -0.1164    0.0824
#> theta[5,4]      0.0687     0.2976    0.0695    0.2987
#> theta[6,4]     -0.0200     0.1922   -0.0196    0.1928
#> theta[7,4]     -0.0432     0.1747   -0.0433    0.1749
#> theta[8,4]     -0.0077     0.1865   -0.0075    0.1870
#> theta[9,4]      0.1390     0.3459    0.1400    0.3472
#> theta[10,4]    -0.0409     0.1477   -0.0407    0.1480
#> theta[6,5]      0.6466     0.9862    0.6517    0.9917
#> theta[7,5]     -0.0091     0.2813   -0.0086    0.2822
#> theta[8,5]      0.0327     0.2991    0.0336    0.3003
#> theta[9,5]     -0.0119     0.2935   -0.0106    0.2952
#> theta[10,5]     0.0591     0.3037    0.0600    0.3048
#> theta[7,6]     -0.1639     0.1083   -0.1643    0.1082
#> theta[8,6]      0.0575     0.2984    0.0584    0.2996
#> theta[9,6]      0.0537     0.3303    0.0549    0.3319
#> theta[10,6]     0.0937     0.3264    0.0947    0.3277
#> theta[8,7]     -0.1452     0.0879   -0.1457    0.0877
#> theta[9,7]      0.0318     0.2953    0.0324    0.2963
#> theta[10,7]    -0.1004     0.1309   -0.1005    0.1310
#> theta[9,8]      0.0998     0.3302    0.1007    0.3314
#> theta[10,8]     0.0847     0.2924    0.0855    0.2934
#> theta[10,9]     0.0416     0.2669    0.0424    0.2680
#> 
#> Negative log pseudo-likelihood: 9654.417 
#> 
#> Savage-Dickey density ratio  [prior: Cauchy( 0 , 2.5 )]
#> H0: theta = 0 for each pairwise interaction
#> 
#>             Estimate     SE     BF_01  Pr(=0|x)
#> theta[2,1]    0.1122 0.0704     12.35    0.9251
#> theta[3,1]    0.4063 0.0695 1.394e-06 1.394e-06
#> theta[4,1]    0.0690 0.0617     24.28    0.9604
#> theta[5,1]    0.2416 0.0884    0.7111    0.4156
#> theta[6,1]   -0.0562 0.0795     29.92    0.9677
#> theta[7,1]    0.4131 0.0758 2.232e-05 2.232e-05
#> theta[8,1]    0.1235 0.0690     10.83    0.9155
#> theta[9,1]    0.2036 0.0701    0.6426    0.3912
#> theta[10,1]  -0.0477 0.0673     30.68    0.9684
#> theta[3,2]    0.1099 0.0517     9.208     0.902
#> theta[4,2]    0.2402 0.0505 0.0002968 0.0002967
#> theta[5,2]    0.2744 0.0615  0.004975  0.004951
#> theta[6,2]    0.2313 0.0564   0.02635   0.02568
#> theta[7,2]    0.1434 0.0575     2.121    0.6796
#> theta[8,2]    0.0940 0.0525     12.54    0.9261
#> theta[9,2]    0.1510 0.0572     1.681     0.627
#> theta[10,2]   0.0902 0.0507     15.88    0.9408
#> theta[4,3]    0.0344 0.0504        51    0.9808
#> theta[5,3]    0.0113 0.0604     58.36    0.9832
#> theta[6,3]   -0.0468 0.0573     38.94     0.975
#> theta[7,3]    0.3783 0.0600 1.316e-07 1.316e-07
#> theta[8,3]   -0.0029 0.0531     61.41     0.984
#> theta[9,3]    0.1145 0.0541     6.342    0.8638
#> theta[10,3]  -0.0168 0.0506     60.07    0.9836
#> theta[5,4]    0.1831 0.0584     0.311    0.2372
#> theta[6,4]    0.0861 0.0541     21.56    0.9557
#> theta[7,4]    0.0657 0.0556     26.24    0.9633
#> theta[8,4]    0.0894 0.0495     8.474    0.8944
#> theta[9,4]    0.2424 0.0528 0.0006504   0.00065
#> theta[10,4]   0.0534 0.0481     41.02    0.9762
#> theta[6,5]    0.8164 0.0867 6.701e-17 6.701e-17
#> theta[7,5]    0.1361 0.0741     10.09    0.9098
#> theta[8,5]    0.1659 0.0680       2.1    0.6774
#> theta[9,5]    0.1408 0.0779      9.86    0.9079
#> theta[10,5]   0.1814 0.0624     1.076    0.5183
#> theta[7,6]   -0.0278 0.0694     46.33    0.9789
#> theta[8,6]    0.1780 0.0615    0.9622    0.4904
#> theta[9,6]    0.1920 0.0706     1.051    0.5125
#> theta[10,6]   0.2100 0.0594   0.03755   0.03619
#> theta[8,7]   -0.0287 0.0595     42.37    0.9769
#> theta[9,7]    0.1636 0.0672     1.809     0.644
#> theta[10,7]   0.0153 0.0590     45.56    0.9785
#> theta[9,8]    0.2150 0.0588   0.07817    0.0725
#> theta[10,8]   0.1886 0.0530   0.08349   0.07705
#> theta[10,9]   0.1543 0.0575     1.058    0.5141
#> 
#> BF_01 of 22 interaction(s) from a normal approximation of the draws (few draws near zero).
#> 
#> SIR samples: 1000   Effective sample size of the importance sample: 2648.7 (of 10000 proposal draws)

The two kinds of intervals agree for most parameters, but not for the thresholds of the rarely chosen highest categories. For item D1, whose highest category was chosen 26 times, the likelihood-ratio interval of $\mu_{13}$ extends further below the estimate than above it, which the symmetric Wald interval cannot do:

confint(fit, parm = c("mu[1,3]", "theta[2,1]"))
#>                   2.5 %     97.5 %
#> mu[1,3]    -10.88102304 -8.1158110
#> theta[2,1]  -0.02504397  0.2512157
confint(fit, parm = c("mu[1,3]", "theta[2,1]"), method = "wald")
#>                   2.5 %     97.5 %
#> mu[1,3]    -10.83569524 -8.0702043
#> theta[2,1]  -0.02580772  0.2501316

plot(fit, type = "intervals") shows the estimates and intervals of all parameters, the thresholds and the interactions in separate panels, each sorted by estimate.

plot(fit, type = "intervals")

Estimates and likelihood-ratio intervals of the thresholds and interactions of the 10 items, sorted by estimate

plot of chunk intervals

When the sample is small relative to the number of parameters and the likelihood-ratio intervals are not computed, dmrfit() warns that the Wald intervals may be inaccurate:

fit_small <- dmrfit(X[1:200, ])
#> Warning in dmrfit(X[1:200, ]): The sample size (200) is small relative to the
#> number of free parameters (75): the Wald intervals assume a symmetric sampling
#> distribution and may be inaccurate. Likelihood-ratio intervals are suggested:
#> refit with lrt_intervals = TRUE (or a vector of parameter names).

Details. The standard errors come from the Godambe-Huber-White (GHW, or sandwich) covariance \[ \boldsymbol{\Sigma}_{\text{GHW}} = \left(-\mathbf{H}^{-1}\right) \mathbf{U} \left(-\mathbf{H}^{-1}\right), \] where $\mathbf{H}$ is the Hessian of the log pseudo-likelihood and $\mathbf{U} = \sum_{\nu=1}^{n} \mathbf{u}_{\nu} \mathbf{u}_{\nu}^{\intercal}$ the variance of the score, both evaluated at the estimate $\hat{\boldsymbol{\eta}}$. The naive covariance $-\mathbf{H}^{-1}$ would underestimate the variability. The likelihood-ratio interval of a parameter $\eta_k$ is based on its profile: for every value of $\eta_k$, the other parameters are set to the values that maximize the pseudo-likelihood. Because it is computed from a pseudo-likelihood, the profile likelihood-ratio statistic is not $\chi^2_1$-distributed but $\lambda_k \chi^2_1$-distributed, with $\lambda_k = [\boldsymbol{\Sigma}_{\text{GHW}}]_{kk} / [-\mathbf{H}^{-1}]_{kk}$ (Geys, Molenberghs and Ryan, 1999; Pace, Salvan and Sartori, 2011). The interval collects the values of $\eta_k$ for which \[ \frac{2}{\lambda_k} \left[\log \tilde{f}(\mathbf{X}; \hat{\boldsymbol{\eta}}) - \max_{\boldsymbol{\eta}_{-k}} \log \tilde{f}(\mathbf{X}; \eta_k, \boldsymbol{\eta}_{-k})\right] \leq \chi^2_{1, 1 - \alpha}, \] with $\chi^2_{1, 1 - \alpha}$ the quantile of the $\chi^2_1$ distribution at the level of the interval. When a prior is used, the pseudo-posterior takes the place of the pseudo-likelihood, and $\lambda_k$ is computed from its curvature $-(\mathbf{H} + \mathbf{H}_{\boldsymbol{\eta}})$, with $\mathbf{H}_{\boldsymbol{\eta}}$ the Hessian of the log prior, and from the posterior GHW covariance $\boldsymbol{\Gamma}\boldsymbol{\Gamma}^{\intercal}_{\text{GHW}} = (\boldsymbol{\Sigma}_{\text{GHW}}^{-1} - \mathbf{H}_{\boldsymbol{\eta}})^{-1}$. The interval follows the shape of the pseudo-likelihood, so it can be asymmetric, and in large samples it agrees with the Wald interval.

A network structure can also be imposed with structure, a matrix with 1 for the edges to estimate and 0 for the edges fixed at zero. For instance, without edges between the two clusters:

S <- outer(clusters[items], clusters[items], "==") * 1
diag(S) <- 0
fit_within <- dmrfit(X, structure = S, with_prior = TRUE)
plot(fit_within, all_edges = TRUE, groups = clusters[items])

Network of the constrained fit without edges between the two clusters, with the nodes grouped by cluster

plot of chunk network-within

Is the edge there? Bayes factors

An interaction estimated close to zero does not tell whether the edge is absent or the data are not informative enough to tell. The Savage-Dickey Bayes factor $\text{BF}_{01}$ (denoted $\text{BF}_{\text{cu}}$ in Arena and Marsman, 2026) compares the hypothesis of conditional independence, $\theta_{ij} = 0$, with the hypothesis $\theta_{ij} \neq 0$: $\text{BF}_{01} > 1$ is evidence for conditional independence (no edge), and $\text{BF}_{01} < 1$ evidence for an interaction. The plots classify the evidence for each edge as

$\text{BF}_{01}$	Evidence
$< 1/10$	included
$1/10$ to $1/3$	weak included
$1/3$ to $3$	inconclusive
$3$ to $10$	weak excluded
$\geq 10$	excluded

evidence <- table(cut(fit$savage_dickey$bf_01, c(0, 1/10, 1/3, 3, 10, Inf), right = FALSE,
                      labels = c("included", "weak included", "inconclusive", "weak excluded", "excluded")))
evidence
#> 
#>      included weak included  inconclusive weak excluded      excluded 
#>            11             1            10             4            19

Of the 45 possible edges, 11 are included and 19 excluded, and the data leave 10 inconclusive. plot(fit, type = "bf") shows the evidence for every edge, with the size of the circle given by the size of the interaction, and the network plot draws the included edges.

plot(fit, type = "bf")

Evidence for every edge from its Savage-Dickey Bayes factor, with the size of the interaction

plot of chunk bf-network

plot(fit, groups = clusters[items])

Network of the included edges, with the nodes grouped by cluster

plot of chunk bf-network

Details. The Savage-Dickey Bayes factor is the ratio of the marginal posterior and the prior density of the interaction at zero, $\text{BF}_{01} = \pi(\theta_{ij} = 0 \mid \mathbf{X}) / \pi(\theta_{ij} = 0)$. dmrfit() estimates the posterior density from draws of the coordinate-rescaled pseudo-posterior (see the next section), obtained by sampling importance resampling (Skare et al., 2003), with a Gaussian kernel estimate. When few draws lie near zero (22 interactions here), a kernel estimate rests on the handful of draws that happen to fall there, or on none, and varies by orders of magnitude between seeds. The posterior density at zero is then that of a normal distribution with the mean and standard deviation of the draws, and summary() reports how many Bayes factors use this approximation.

Bayesian estimation with coordinate rescaling

dmrfit_bayes() samples from the posterior of the parameters. By default, it samples the coordinate-rescaled pseudo-posterior, which has the location of the pseudo-posterior and the scale of the GHW covariance.

fit_bayes <- dmrfit_bayes(X, nsim = 10000, burnin = 1000, progress = FALSE)

We first assess the convergence and mixing of the chain. The trace plot shows the draws of the interactions with the largest absolute posterior modes, which should fluctuate around a stable level, without trends. The effective sample sizes (ESS) summarize the mixing.

plot(fit_bayes, type = "trace")

Trace plots of the posterior draws of the four interactions with the largest absolute posterior modes

plot of chunk trace

min(fit_bayes$savage_dickey$ess)  # smallest effective sample size of a parameter
#> [1] 766.0416
fit_bayes$mess                    # multivariate effective sample size of all draws
#> [1] 1933.92

The posterior densities show the posterior mode (solid line) and the 95% highest posterior density interval (dashed lines), the shortest interval that contains 95% of the posterior.

plot(fit_bayes, type = "density", pars = c("mu[1,3]", "mu[1,1]", "theta[3,1]", "theta[6,5]"))

Posterior densities of two thresholds and two interactions, with their posterior modes and highest posterior density intervals

plot of chunk density

summary() reports the posterior summaries and the Bayes factors computed from the draws, and confint() the highest posterior density intervals. The Bayes factors from the draws and from dmrfit() agree closely: they classify 84% of the edges in the same evidence class.

summary(fit_bayes)
#> 
#> Call:
#> dmrfit_bayes(data = X, nsim = 10000, burnin = 1000, progress = FALSE)
#> 
#> Discrete Markov Random Field
#> Estimation method: Bayesian (FisherMALA), coordinate-rescaled pseudo-posterior (Godambe-Huber-White covariance) 
#> Nodes: 10  Observations: 917  Parameters: 75 
#> MCMC samples: 10000  Acceptance rate: 0.58  Elapsed: 10.1 sec
#> Multivariate effective sample size: 1933.9 
#> ------------------------------------------------------------ 
#> 
#> Thresholds:
#>          Post.Mode Post.Mean Post.SD     2.5%   97.5%
#> mu[1,1]    -1.8209   -1.8188  0.1395  -2.0979 -1.5464
#> mu[1,2]    -5.5766   -5.5751  0.3616  -6.2710 -4.8439
#> mu[1,3]    -9.4529   -9.4878  0.6919 -10.8343 -8.0737
#> mu[2,1]    -1.1693   -1.1668  0.1153  -1.3927 -0.9417
#> mu[2,2]    -2.4177   -2.4079  0.2053  -2.8119 -2.0039
#> mu[2,3]    -5.5241   -5.5080  0.3852  -6.2659 -4.7490
#> mu[3,1]    -0.2099   -0.2047  0.1082  -0.4187  0.0051
#> mu[3,2]    -1.5830   -1.5685  0.1866  -1.9406 -1.1932
#> mu[3,3]    -3.4000   -3.3860  0.3279  -4.0245 -2.7330
#> mu[4,1]    -1.4515   -1.4440  0.1178  -1.6751 -1.2134
#> mu[4,2]    -3.0871   -3.0716  0.2364  -3.5449 -2.6303
#> mu[4,3]    -5.6103   -5.6014  0.4100  -6.3932 -4.8180
#> mu[5,1]    -1.4684   -1.4693  0.1432  -1.7518 -1.1898
#> mu[5,2]    -3.8278   -3.8309  0.3338  -4.4747 -3.1541
#> mu[5,3]    -8.6933   -8.6997  0.6575  -9.9772 -7.3785
#> mu[6,1]    -1.4455   -1.4359  0.1329  -1.6923 -1.1761
#> mu[6,2]    -3.8509   -3.8286  0.2882  -4.3750 -3.2610
#> mu[6,3]    -7.7036   -7.6699  0.5119  -8.6602 -6.6777
#> mu[7,1]    -1.1122   -1.1064  0.1281  -1.3507 -0.8618
#> mu[7,2]    -3.4381   -3.4339  0.2844  -3.9912 -2.8898
#> mu[7,3]    -6.2229   -6.2410  0.5738  -7.3498 -5.1137
#> mu[8,1]    -0.8508   -0.8462  0.1072  -1.0557 -0.6330
#> mu[8,2]    -2.5134   -2.5002  0.2152  -2.9206 -2.0884
#> mu[8,3]    -5.4577   -5.4548  0.4087  -6.2590 -4.6650
#> mu[9,1]    -2.2830   -2.2905  0.1310  -2.5469 -2.0352
#> mu[9,2]    -4.6340   -4.6486  0.2780  -5.1889 -4.1002
#> mu[9,3]    -8.2423   -8.2962  0.5092  -9.3031 -7.2752
#> mu[10,1]    0.0458    0.0581  0.1088  -0.1577  0.2726
#> mu[10,2]   -1.0562   -1.0324  0.1743  -1.3760 -0.6863
#> mu[10,3]   -3.0828   -3.0432  0.3065  -3.6423 -2.4444
#> 
#> Pairwise interactions:
#>             Post.Mode Post.Mean Post.SD    2.5%  97.5%
#> theta[2,1]     0.1122    0.1194  0.0696 -0.0167 0.2629
#> theta[3,1]     0.4063    0.4076  0.0704  0.2682 0.5482
#> theta[4,1]     0.0690    0.0723  0.0618 -0.0505 0.1925
#> theta[5,1]     0.2416    0.2439  0.0897  0.0616 0.4203
#> theta[6,1]    -0.0562   -0.0600  0.0795 -0.2136 0.0949
#> theta[7,1]     0.4131    0.4228  0.0765  0.2742 0.5764
#> theta[8,1]     0.1235    0.1232  0.0681 -0.0086 0.2553
#> theta[9,1]     0.2036    0.1976  0.0701  0.0604 0.3357
#> theta[10,1]   -0.0477   -0.0567  0.0671 -0.1881 0.0739
#> theta[3,2]     0.1099    0.1112  0.0516  0.0102 0.2137
#> theta[4,2]     0.2402    0.2414  0.0508  0.1391 0.3382
#> theta[5,2]     0.2744    0.2779  0.0623  0.1554 0.3986
#> theta[6,2]     0.2313    0.2298  0.0577  0.1183 0.3453
#> theta[7,2]     0.1434    0.1386  0.0582  0.0245 0.2515
#> theta[8,2]     0.0940    0.0943  0.0526 -0.0052 0.1980
#> theta[9,2]     0.1510    0.1516  0.0578  0.0376 0.2624
#> theta[10,2]    0.0902    0.0861  0.0504 -0.0140 0.1848
#> theta[4,3]     0.0344    0.0327  0.0511 -0.0666 0.1294
#> theta[5,3]     0.0113    0.0116  0.0623 -0.1108 0.1299
#> theta[6,3]    -0.0468   -0.0473  0.0589 -0.1611 0.0698
#> theta[7,3]     0.3783    0.3797  0.0598  0.2625 0.4985
#> theta[8,3]    -0.0029   -0.0040  0.0551 -0.1145 0.1017
#> theta[9,3]     0.1145    0.1135  0.0541  0.0085 0.2210
#> theta[10,3]   -0.0168   -0.0193  0.0512 -0.1218 0.0786
#> theta[5,4]     0.1831    0.1835  0.0597  0.0668 0.2999
#> theta[6,4]     0.0861    0.0844  0.0557 -0.0248 0.1951
#> theta[7,4]     0.0657    0.0637  0.0546 -0.0460 0.1672
#> theta[8,4]     0.0894    0.0921  0.0500 -0.0058 0.1890
#> theta[9,4]     0.2424    0.2442  0.0530  0.1392 0.3488
#> theta[10,4]    0.0534    0.0504  0.0494 -0.0460 0.1457
#> theta[6,5]     0.8164    0.8184  0.0863  0.6444 0.9837
#> theta[7,5]     0.1361    0.1326  0.0729 -0.0089 0.2772
#> theta[8,5]     0.1659    0.1669  0.0666  0.0382 0.2952
#> theta[9,5]     0.1408    0.1391  0.0797 -0.0214 0.2917
#> theta[10,5]    0.1814    0.1795  0.0633  0.0521 0.3006
#> theta[7,6]    -0.0278   -0.0305  0.0704 -0.1677 0.1046
#> theta[8,6]     0.1780    0.1764  0.0611  0.0531 0.2947
#> theta[9,6]     0.1920    0.1944  0.0725  0.0549 0.3348
#> theta[10,6]    0.2100    0.2076  0.0610  0.0905 0.3291
#> theta[8,7]    -0.0287   -0.0306  0.0595 -0.1505 0.0873
#> theta[9,7]     0.1636    0.1695  0.0646  0.0414 0.2971
#> theta[10,7]    0.0153    0.0188  0.0597 -0.1008 0.1371
#> theta[9,8]     0.2150    0.2140  0.0582  0.0956 0.3253
#> theta[10,8]    0.1886    0.1873  0.0532  0.0849 0.2924
#> theta[10,9]    0.1543    0.1595  0.0564  0.0516 0.2716
#> 
#> Negative log pseudo-likelihood at posterior mode: 9654.417 
#> 
#> 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.1122  0.0696     10.16    0.9104
#> theta[3,1]     0.4063  0.0704 2.322e-06 2.322e-06
#> theta[4,1]     0.0690  0.0618     26.44    0.9636
#> theta[5,1]     0.2416  0.0897    0.8679    0.4646
#> theta[6,1]    -0.0562  0.0795     31.81    0.9695
#> theta[7,1]     0.4131  0.0765 9.397e-06 9.397e-06
#> theta[8,1]     0.1235  0.0681     8.536    0.8951
#> theta[9,1]     0.2036  0.0701    0.8434    0.4575
#> theta[10,1]   -0.0477  0.0671     31.96    0.9697
#> theta[3,2]     0.1099  0.0516     5.954    0.8562
#> theta[4,2]     0.2402  0.0508 0.0007707 0.0007701
#> theta[5,2]     0.2744  0.0623  0.002429  0.002423
#> theta[6,2]     0.2313  0.0577   0.01941   0.01904
#> theta[7,2]     0.1434  0.0582     3.146    0.7588
#> theta[8,2]     0.0940  0.0526     11.76    0.9216
#> theta[9,2]     0.1510  0.0578     1.748    0.6361
#> theta[10,2]    0.0902  0.0504     13.65    0.9317
#> theta[4,3]     0.0344  0.0511     49.41    0.9802
#> theta[5,3]     0.0113  0.0623     47.42    0.9793
#> theta[6,3]    -0.0468  0.0589     36.62    0.9734
#> theta[7,3]     0.3783  0.0598 9.326e-08 9.326e-08
#> theta[8,3]    -0.0029  0.0551     56.11    0.9825
#> theta[9,3]     0.1145  0.0541     6.286    0.8628
#> theta[10,3]   -0.0168  0.0512     55.62    0.9823
#> theta[5,4]     0.1831  0.0597    0.4648    0.3173
#> theta[6,4]     0.0861  0.0557     18.56    0.9489
#> theta[7,4]     0.0657  0.0546     28.29    0.9659
#> theta[8,4]     0.0894  0.0500      11.1    0.9174
#> theta[9,4]     0.2424  0.0530  0.001444  0.001442
#> theta[10,4]    0.0534  0.0494     38.16    0.9745
#> theta[6,5]     0.8164  0.0863 1.086e-18 1.086e-18
#> theta[7,5]     0.1361  0.0729     7.823    0.8867
#> theta[8,5]     0.1659  0.0666     2.038    0.6709
#> theta[9,5]     0.1408  0.0797     8.712     0.897
#> theta[10,5]    0.1814  0.0633    0.8821    0.4687
#> theta[7,6]    -0.0278  0.0704      41.2    0.9763
#> theta[8,6]     0.1780  0.0611    0.7921     0.442
#> theta[9,6]     0.1920  0.0725     1.184    0.5422
#> theta[10,6]    0.2100  0.0610    0.1559    0.1348
#> theta[8,7]    -0.0287  0.0595     47.88    0.9795
#> theta[9,7]     0.1636  0.0646     1.554    0.6084
#> theta[10,7]    0.0153  0.0597      49.3    0.9801
#> theta[9,8]     0.2150  0.0582     0.063   0.05926
#> theta[10,8]    0.1886  0.0532    0.1206    0.1076
#> theta[10,9]    0.1543  0.0564     1.021    0.5053
#> 
#> BF_01 of 22 interaction(s) from a normal approximation of the draws (few draws near zero).
#> 
#> MCMC samples: 10000   Effective sample size of the interactions: min 861.1 , median 1066.2
confint(fit_bayes, parm = c("mu[1,3]", "theta[2,1]"))
#>                   2.5 %     97.5 %
#> mu[1,3]    -10.73963450 -7.9904352
#> theta[2,1]  -0.02117394  0.2553808

Details. The pseudo-posterior $\tilde{\pi}(\boldsymbol{\eta} \mid \mathbf{X})$ is approximately centered at the right location, but its scale is too small. CoRe defines the transformed parameters \[ \boldsymbol{\beta} = \mathbf{A}\left(\boldsymbol{\eta} - \boldsymbol{\eta}^{\star}\right) + \boldsymbol{\eta}^{\star}, \qquad \mathbf{A} = \boldsymbol{\Gamma} \mathbf{L}^{\intercal}, \] where $\boldsymbol{\eta}^{\star}$ is the pseudo-posterior mode, $\mathbf{L}^{\intercal}$ the upper-triangular Cholesky factor of the negative pseudo-posterior Hessian, $\mathbf{L}\mathbf{L}^{\intercal} = -(\mathbf{H} + \mathbf{H}_{\boldsymbol{\eta}})$, and $\boldsymbol{\Gamma}$ the lower-triangular Cholesky factor of the GHW posterior covariance $\boldsymbol{\Gamma}\boldsymbol{\Gamma}^{\intercal}_{\text{GHW}}$. The transformed pseudo-posterior has the covariance $\boldsymbol{\Gamma}\boldsymbol{\Gamma}^{\intercal}_{\text{GHW}}$ to leading order, and it is sampled with the Fisher adaptive Metropolis-adjusted Langevin algorithm (FisherMALA, Titsias, 2024). The plots use the marginal posterior mode as the point summary by default (estimate = "mode"): in small samples the posterior of a threshold of a rarely chosen category is asymmetric, which moves the mean and the median away from the mode (Arena and Marsman, 2026).

Centrality

The expected influence of a variable is the sum of its interactions, $\sum_{j} \theta_{ij}$ (Robinaugh, Millner and McNally, 2016), which measures how strongly it is connected to the rest of the network, with the sign of each connection. For a dmrfit_bayes fit, it is computed on every draw, which also gives the posterior probability that each variable is the most central.

plot(fit_bayes, type = "centrality")

Expected influence of the 10 items with highest posterior density intervals and the posterior probability of being the most central item

plot of chunk centrality

Item D7 has the largest expected influence, and is the most central item in 99% of the posterior draws. For a dmrfit fit, plot(fit, type = "centrality") shows the estimates with Wald intervals.

Which method when

CoRe approximates the posterior based on the full likelihood. For small networks, this posterior can be sampled exactly, because its normalizing constant can be computed by enumerating all $(m + 1)^p$ response patterns. This gives a direct check of the approximation, here on the first five items ($4^5 = 1{,}024$ response patterns):

five <- items[1:5]
fit_core <- dmrfit_bayes(rads2[, five], nsim = 5000, burnin = 1000, progress = FALSE)
fit_exact <- dmrfit_bayes(rads2[, five], method = "exact", nsim = 5000, burnin = 1000, progress = FALSE)
sd_ratio <- apply(fit_core$draws, 1, sd) / apply(fit_exact$draws, 1, sd)
summary(sd_ratio)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.9697  1.0100  1.0556  1.0512  1.0821  1.1636

The posterior standard deviations of CoRe are close to those of the exact posterior, and slightly larger on average (median ratio 1.06, largest 1.16). This agrees with the theory, because the GHW covariance is at least as large as the inverse Fisher information, with equality when the pseudo-likelihood is efficient, and the residual gap grows for small $n$ relative to the number of parameters (Arena and Marsman, 2026). The other options of dmrfit_bayes() are the following.

method = "core" (default): the coordinate-rescaled pseudo-posterior, as fast as sampling the pseudo-posterior.
method = "adacore": adapts the rescaling during burn-in, at the running mean of the draws.
method = "exact": the posterior based on the full likelihood, by enumeration, only for small networks (up to control$max_states response patterns).
method = "dmh": the posterior based on the full likelihood with the double Metropolis-Hastings algorithm (Liang, 2010), for any network size, but much slower: a reference rather than a routine choice.
scale, for method = "core": the covariance the pseudo-posterior is rescaled to, the GHW covariance ("ghw", default), or the inverse of a Monte Carlo estimate of the negative Hessian of the full log-posterior, at the pseudo-posterior mode ("mch") or at the Robbins-Monro estimate of the full-posterior mode ("rm"). The last two simulate data from the model and are slower.

References

Arena, G. and Marsman, M. (2026). Bayesian inference for discrete Markov random fields through coordinate rescaling. arXiv preprint. https://doi.org/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. https://doi.org/10.1080/01621459.1999.10474176
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.
Pace, L., Salvan, A., and Sartori, N. (2011). Adjusting composite likelihood ratio statistics. Statistica Sinica, 21(1), 129-148.
Ramos-Vera, C., Quispe Callo, G., Basauri Delgado, M., Vallejos Saldarriaga, J., and Saintila, J. (2023). Factorial and network structure of the Reynolds Adolescent Depression Scale (RADS-2) in Peruvian adolescents. PLOS ONE, 18(5), e0286081.
Robinaugh, D. J., Millner, A. J., and McNally, R. J. (2016). Identifying highly influential nodes in the complicated grief network. Journal of Abnormal Psychology, 125(6), 747-757.
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.
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.

\(\text{BF}_{01}\)	Evidence
\(< 1/10\)	included
\(1/10\) to \(1/3\)	weak included
\(1/3\) to \(3\)	inconclusive
\(3\) to \(10\)	weak excluded
\(\geq 10\)	excluded