Introduction
A fitted graphical model raises questions of two kinds. Which
edges does the evidence settle? is a question about the posterior,
and verdicts() answers it. Does the model describe the
data? is a question about fit, and it splits in two: whether the
conditional prediction for each variable is calibrated, which
calibration_check() answers, and whether the joint
distribution is reproduced, which you check by building a display on
simulate().
This vignette works through all of them on the Wenchuan data. The
convergence diagnostics – R-hat, effective sample size, the trust gauge
– are a separate matter and live in the diagnostics vignette.
library(bgms)
data = Wenchuan[, 1:9]
fit = bgm(data, seed = 1234)
What the evidence settles
Each edge carries an inclusion Bayes factor, and reading it at a
threshold gives one of three verdicts: evidence of presence, evidence of
absence, or undecided. The third is not a failure of the analysis; it is
the honest answer when the data do not separate the two hypotheses.
edges = verdicts(fit)
edges
#> Edge verdicts at an inclusion Bayes factor of 10 (and 0.1 for absence):
#> presence: log BF > 2.30; absence: log BF < -2.30
#>
#> presence 11 | undecided 15 | absence 10 (36 indicators)
#>
#> parameter pip log_bf verdict fragile
#> intrusion-dreams 1.000 291.671 presence FALSE
#> intrusion-flash 1.000 16.221 presence FALSE
#> intrusion-upset 0.947 2.887 presence TRUE
#> intrusion-physior 0.632 0.539 undecided FALSE
#> intrusion-avoidth 0.026 -3.615 absence FALSE
#> intrusion-avoidact 0.038 -3.244 absence FALSE
#> intrusion-amnesia 0.036 -3.275 absence FALSE
#> intrusion-lossint 0.165 -1.620 undecided FALSE
#> dreams-flash 1.000 132.236 presence FALSE
#> dreams-upset 0.667 0.695 undecided FALSE
#> ... (26 more rows)
#>
#> 7 verdicts are Monte-Carlo fragile: a verdict boundary lies within two standard
#> errors of the evidence, so the verdict could change on a rerun. Consider a
#> longer run.
The fragile column is the part that cannot be read off
the Bayes factor alone. A verdict is fragile when a threshold sits
within two standard errors of the estimated evidence, which is exactly
the regime in which Monte Carlo noise, not the data, decides the answer.
In a known-truth calibration study of 37,010 edge-fits across ordinal,
binary, and Gaussian graphical models, every verdict error sat within
0.58 of a threshold on the log Bayes factor scale, and no edge further
out was ever misclassified.
A fragile verdict is not a wrong verdict; it is a verdict the run is
too short to settle. The remedy is more iterations, not a different
model.
edges[edges$fragile, c("parameter", "log_bf", "verdict")]
#> parameter log_bf verdict
#> 3 intrusion-upset 2.887300 presence
#> 12 dreams-avoidth -2.208473 undecided
#> 13 dreams-avoidact -2.524920 absence
#> 16 flash-upset -2.464686 absence
#> 20 flash-amnesia 2.837967 presence
#> 24 upset-avoidact -2.064197 undecided
#> 30 physior-lossint 2.232483 undecided
The default picture encodes the same three verdicts, so what a reader
sees and what the table says cannot drift apart.
There is one panel per verdict – the pairs the data support, the
pairs the data rule out, and the pairs the data cannot decide – drawn on
one shared layout, so a variable sits in the same place in all three.
Only the first panel is weighted, with line width the posterior mean
association and colour its sign; the other two are drawn at uniform
width, dashed for a ruled-out pair and dotted for an undecided one,
because there the classification is the result. A single drawing would
have to leave the last two classes as the same blank, which is exactly
what a reader cannot interpret. For one edge in detail, the
spike-and-slab posterior shows how the mass is actually divided:
plot_edge_posterior(fit, "intrusion", "upset")
Centrality is a function of the graph, so it inherits the same
uncertainty. Evaluating it on every posterior draw gives its posterior
directly, which is what extract_centrality() returns.
strength = extract_centrality(fit)
head(summary(strength))
#> node mean lower upper p_most_central
#> 6 avoidth 0.8754598 0.7461277 1.0248335 0.5605
#> 7 avoidact 0.8570654 0.7355831 0.9941833 0.3965
#> 4 upset 0.7399619 0.6218007 0.8629425 0.0310
#> 2 dreams 0.7051229 0.5926109 0.8192193 0.0085
#> 5 physior 0.6776066 0.5664244 0.8046862 0.0025
#> 1 intrusion 0.6626375 0.5596425 0.7838838 0.0010
plot(strength)
Are the conditional predictions calibrated?
The model’s person-level predictions are the conditional
distributions P(x_ij | the rest of case i's responses),
which is the model’s own regression unit.
calibration_check() compares those predictions with what
was observed, one variable at a time: an isotonic fit of observed
frequency on predicted probability for a discrete variable, and the
probability integral transform for a continuous one.
calibration = calibration_check(fit, nrep = 200, seed = 1)
calibration
#> Calibration of the conditional predictions, 95% consistency band from 200 resamples:
#>
#> variable kind mean_dev max_dev share_outside_band
#> avoidth pav 0.075 0.269 0.248
#> dreams pav 0.052 0.192 0.168
#> intrusion pav 0.044 0.148 0.109
#> avoidact pav 0.060 0.142 0.267
#> flash pav 0.047 0.141 0.129
#> lossint pav 0.036 0.121 0.149
#> amnesia pav 0.039 0.108 0.109
#> physior pav 0.030 0.088 0.059
#> upset pav 0.036 0.087 0.040
#>
#> A calibrated variable tracks the diagonal and stays inside its band;
#> 'share_outside_band' is the proportion of the curve that does not, on a
#> 0 to 1 scale (0.3 is 30% of the curve, not 0.3%).
One panel per variable is not a display preference. Pooling across
variables lets a variable predicted too high and one predicted too low
cancel, and the pooled curve then tracks the diagonal while neither
variable does.
The grey band is the wander of a model that is calibrated by
construction. It is built by resampling each case’s category
from its own predicted distribution, rather than resampling the
threshold events independently: the cumulative events of one case are
nested, and treating them as independent would make the band too narrow
and an ordinary wander look like miscalibration.
A continuous variable has no category to have fallen in, so it gets
the same question through the probability integral transform:
u = F(y | rest), which is uniform exactly when the
conditional predictive distribution is right. Its panel is the
distribution of those u values against the uniform
diagonal, and F is the mixture over posterior draws, so the
parameter uncertainty is inside the distribution the observation is
transformed by. That band is not resampled from the model at all – under
the transform the null is uniform whatever the density was – so it
depends only on the number of cases. A mixed fit draws both panel kinds
on one figure and records which is which in the kind
column.
Checking the joint distribution with simulate()
Calibrated conditional predictions do not imply that the model
reproduces the joint distribution: predicting each variable from its
neighbours is an easier task than reproducing how all of them covary,
and a model can do the first while falling short on the second.
bgms ships no joint-level check, and the omission is deliberate. What
does have diagnostic value is reading a display for systematic pattern,
and that is an analyst’s judgement rather than a pass/fail
statistic.
Which statistic you display matters more than it looks. Pairwise
dependence statistics – correlations, covariances, cross-tabulations –
sit close to the model’s sufficient set, so how they read depends on how
the estimator is anchored rather than on whether the model fits. The
sum score lies outside that set: no anchoring makes it
match by construction, so a departure is about the model. It also has a
one-sentence reading, which a pairwise statistic does not.
simulate() gives you the material. One replicated
dataset per posterior draw propagates both parameter and structure
uncertainty into the replicates:
observed = as.matrix(data[complete.cases(data), ])
replicates = simulate(fit,
nsim = nrow(observed), method = "posterior-sample",
ndraws = 200, seed = 1, display_progress = "none"
)
scores = 0:(ncol(observed) * max(observed))
distribution = function(x) {
as.numeric(table(factor(rowSums(x), levels = scores))) / nrow(x)
}
observed_share = distribution(observed)
replicated_share = vapply(replicates, distribution, observed_share)
band = apply(replicated_share, 1, quantile, probs = c(0.025, 0.975))
Draw the observed distribution against its pointwise band rather than
counting how often the two agree. The count is the weak reading; the
informative one is whether the observed curve leaves the band in a
direction:
plot(scores, observed_share,
type = "n", las = 1, xlab = "Sum score", ylab = "Share of cases",
ylim = range(c(band, observed_share))
)
polygon(c(scores, rev(scores)), c(band[1, ], rev(band[2, ])),
col = adjustcolor("grey55", 0.25), border = NA
)
lines(scores, observed_share, col = "#0072B2", lwd = 1.8)
mean(observed_share < band[1, ] | observed_share > band[2, ])
#> [1] 0.1086957
The grey band is where a dataset generated by the fitted model puts
its sum scores; the blue curve is where the real one puts them. Read the
shape, not the count: a curve that wanders in and out of the band is
ordinary sampling variation, whereas one that sits above the band across
a whole stretch of the scale and below it across another is the model
getting the spread of total severity wrong, which is what a graphical
model can miss while predicting each symptom from its neighbours
perfectly well.
Here the departures are scattered – a few isolated scores at the
sparse low end of the scale and one at the ceiling, with no run in
either direction. That is the ordinary-variation reading, and the
printed share is what makes the point about counts: about a tenth of the
scale sits outside a 95% pointwise band, which sounds alarming and is
not, because the band is pointwise and the scale has 46 points. The
shape is what carries the information.
Any statistic outside the model’s sufficient set works the same way –
a tail probability, the number of endorsed symptoms, whatever the
substantive question turns on.