---
title: "Checking your fitted model"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Checking your fitted model}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: refs.bib
csl: apa.csl
link-citations: TRUE
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 3.4
)
```

# 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.

```{r, eval=FALSE}
library(bgms)
data = Wenchuan[, 1:9]
fit = bgm(data, seed = 1234)
```

```{r, include=FALSE}
library(bgms)
data = Wenchuan[, 1:9]
fit = bgm(data,
  seed = 1234, chains = 2, iter = 1000, warmup = 1000,
  display_progress = "none", verbose = FALSE
)
```

# 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.

```{r}
edges = verdicts(fit)
edges
```

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.

```{r}
edges[edges$fragile, c("parameter", "log_bf", "verdict")]
```

The default picture encodes the same three verdicts, so what a reader sees and
what the table says cannot drift apart.

```{r, fig.width = 12, fig.height = 4.5, out.width = "100%"}
plot(fit)
```

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:

```{r, fig.height = 3.6, fig.width = 5}
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.

```{r, fig.height = 3.6, fig.width = 6}
strength = extract_centrality(fit)
head(summary(strength))
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.

```{r}
calibration = calibration_check(fit, nrep = 200, seed = 1)
calibration
```

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.

```{r, fig.height = 5.5, fig.width = 7}
plot(calibration)
```

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:

```{r}
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:

```{r, fig.height = 3.8, fig.width = 6}
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, ])
```

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.

# What the three checks answer

They are three different questions, and answering one says nothing about the
other two.

- `verdicts()`: which edges the data settle, and which verdicts the run is too
  short to settle.
- `calibration_check()`: whether to trust a prediction for one variable.
- A `simulate()`-based display: whether to trust the model as a description
  of the joint distribution.

An edge can be firmly settled in a model that reproduces the joint
distribution poorly, and a model can predict each variable well while
understating the dependence among them.

# Next steps

- The diagnostics vignette covers convergence: R-hat, effective sample size,
  the NUTS diagnostics, and the hierarchical prior trust gauge.
- `prior_sensitivity_check()` reports whether a verdict depends on the prior
  rather than on the run length, which is the other way a verdict can be
  unstable.
