Gaussian Graphical Model
The Gaussian graphical model (GGM) is a Markov random field for continuous variables. It is the most widely studied undirected graphical model (Lauritzen, 1996) and forms the continuous component of the Mixed MRF in bgms. To fit a GGM, set variable_type = "continuous".
The model as a regression
For continuous variables, the model answers the same question as the Ordinal MRF: given all the other variables, what does this one do? Each variable is an ordinary linear regression on the rest. The threshold’s job from the ordinal model is handled before the regression starts: bgms first subtracts each variable’s mean, so that every score becomes a deviation from that variable’s own average. A variable’s overall level is thereby removed from the data and cannot masquerade as association; the model only describes how variables move together around their averages. What remains is the pull of the neighbors:
\[ E\bigl(x_i \mid \text{rest}\bigr) \;=\; \frac{1}{\theta_{ii}}\, \underbrace{\textstyle\sum_{j \neq i} 2\,\omega_{ij}\, x_j}_{\text{the pull of the neighbors}} . \]
The pull is the same expression as in the ordinal regression: each neighboring value \(x_j\) enters with twice its edge weight \(\omega_{ij}\), conditional on all remaining variables. Only the noise around the regression changes, Gaussian rather than categorical. The factor \(1/\theta_{ii}\) is the residual variance of \(x_i\), the variance that remains after its neighbors have done their work, and it converts the pull into the units of \(x_i\).
Together, the variable-level regressions define one joint model, a multivariate normal distribution whose precision matrix collects all the parameters. That joint form is the model below, and it is also the model behind the familiar EBICglasso analysis: readers arriving from that workflow have been fitting it all along, and what changes here is not the model but the question asked of it.
The model
A GGM assumes that \(p\) continuous variables follow a multivariate normal distribution:
\[ \mathbf{x} \sim \mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Sigma}) \]
where \(\boldsymbol{\mu}\) is the mean vector and \(\boldsymbol{\Sigma}\) is the \(p \times p\) covariance matrix. The covariance matrix describes marginal associations between pairs of variables, but marginal associations confound direct and indirect effects. Two variables may be strongly correlated because both depend on a third variable, not because they influence each other directly.
The GGM is defined not through the covariance matrix but through its inverse, the precision matrix \(\boldsymbol{\Theta} = \boldsymbol{\Sigma}^{-1}\). Dempster (1972) introduced this idea under the name covariance selection: the precision matrix separates direct from indirect associations, and a zero entry \(\theta_{ij} = 0\) means that variables \(i\) and \(j\) are conditionally independent given all other variables. In the graph, this corresponds to the absence of the edge between nodes \(i\) and \(j\).
Connection to the MRF framework
The Graphical Models page introduced the general pairwise MRF:
\[ p(\mathbf{x}) \propto \exp\!\left(\sum_{i} \mu_i(x_i) + \mathbf{x}^{\sf T}\boldsymbol{\Omega}\, \mathbf{x}\right) \]
The GGM is a special case of this form. In bgms, continuous variables are centered before estimation, so \(\boldsymbol{\mu} = \mathbf{0}\) and the node potentials vanish. The multivariate normal density then reduces to:
\[ p(\mathbf{x}) \propto \exp\!\left( - \tfrac{1}{2}\,\mathbf{x}^{\sf T} \boldsymbol{\Theta}\, \mathbf{x}\right) \]
Comparing this with the general MRF identifies the interaction matrix as \(\boldsymbol{\Omega} = -\tfrac{1}{2}\,\boldsymbol{\Theta}\). The off-diagonal entry \(\omega_{ij} = -\tfrac{1}{2}\,\theta_{ij}\) is the partial association between variables \(i\) and \(j\), and when \(\theta_{ij} = 0\) the partial association is zero and the edge is absent.
Partial associations
As described in Graphical Models, bgms reports all pairwise effects on a unified partial association scale. For the GGM, the partial association between variables \(i\) and \(j\) is:
\[ \omega_{ij} = -\tfrac{1}{2}\, \theta_{ij} \]
A positive \(\omega_{ij}\) means that, holding all other variables constant, higher values of \(x_i\) are associated with higher values of \(x_j\).
The partial association and the precision element carry the same information — one is a linear rescaling of the other. The rescaling ensures that the GGM partial association occupies the same structural role as the partial association in the ordinal MRF, where \(\omega_{ij}\) is half the log adjacent-category odds ratio (see Ordinal MRF). This unified scale means that coef(fit)$pairwise returns partial associations regardless of model type, and the spike-and-slab prior (Prior Basics) operates on the same quantity across models. To obtain the precision matrix directly, use extract_precision(fit).
Partial correlations
The partial correlation is a standardized measure of conditional association. It removes the influence of the variables’ scales and is bounded between \(-1\) and \(1\):
\[ \rho_{ij \cdot \text{rest}} = -\frac{\theta_{ij}}{\sqrt{\theta_{ii}\, \theta_{jj}}} \]
Because the partial correlation standardizes the off-diagonal precision element by the conditional standard deviations of both variables, it is invariant to rescaling. Two datasets that differ only in their units of measurement produce the same partial correlations.
The partial correlation is related to the partial association by:
\[ \rho_{ij \cdot \text{rest}} = \frac{2\,\omega_{ij}}{\sqrt{\theta_{ii}\, \theta_{jj}}} \]
When comparing the strength of associations across variable pairs — identifying which edges are strongest relative to each variable’s conditional variability — partial correlations are the natural summary. When the goal is edge selection or comparison with discrete graphical models, the partial associations are more appropriate. To obtain the partial correlation matrix directly, use extract_partial_correlations(fit).
The precision matrix
The diagonal entries of the precision matrix, \(\theta_{ii}\), are the inverse of the residual variance of variable \(i\) — the variance that remains after regressing \(x_i\) on all other variables. A small residual variance \(1 / \theta_{ii}\) means the variable is well predicted by the rest of the network.
In bgms, the diagonal entries are nuisance parameters. They receive an exponential prior by default (exponential_prior(eta = 1), with the rate given in the standardized frame described in Prior Basics) and are estimated during MCMC sampling, but they are not subject to edge selection. The off-diagonal entries — or equivalently the partial associations — are the primary inferential targets.
Likelihood
Unlike the ordinal MRF, whose normalizing constant is intractable and requires a pseudolikelihood approximation (see Ordinal MRF), the GGM has a closed-form likelihood — estimates for continuous data involve no approximation. The exact expression and how it is evaluated are given in GGM Internals.
Edge selection
When edge_selection = TRUE (the default), each partial association \(\omega_{ij}\) — and therefore each precision entry \(\theta_{ij}\) — receives a spike-and-slab prior, so an edge is either exactly absent or estimated from the data. The prior specification is described in Prior Basics, and how to read the resulting inclusion probabilities and Bayes factors in Edge Selection.
For continuous data, the edge prior interacts with the precision-matrix prior, because excluding an edge constrains a matrix that must remain positive definite. The precision_graph_prior argument of bgm() controls how the two priors compose:
"hierarchical"(the default) normalizes the precision prior per graph, so the prior over structures is exactly the edge prior and nominal prior inclusion probabilities apply. Each edge move then evaluates a normalizer ratio through a fast approximation, and a trust gauge audits that approximation while the chains run (see MCMC Diagnostics). That approximation is built once per analysis before the chains start, which is the pause you see on a first hierarchical fit; the build announces itself with the size it covers and the number of cores it uses, then reports its elapsed time. It is cached in the session and on disk, so a repeat fit of the same configuration reuses it silently. Setoptions(bgms.verbose = FALSE)to mute the announcement."joint"combines them in un-normalized form instead. The prior probability that an edge is included then differs from the nominal edge-prior value; useextract_prior_inclusion_probabilities()for the effective values when computing Bayes factors by hand (see Edge Selection). A fit reports this once when it starts. Withbeta_bernoulli_prior()orsbm_prior(), the hyperparameter updates additionally carry a normalizing-constant correction. The correction table is built from the prior at the first fit of a model configuration — a one-time cost of the order of minutes, announced when it actually builds — and cached on disk, so later fits of the same configuration skip the build. Withbeta_bernoulli_prior(), the sampled inclusion probability is returned per chain infit$inclusion_parameter_samples.
The two specifications differ only where the sampler moves between graphs. With edge_selection = FALSE they coincide exactly, and on data without a continuous block the argument simply does not apply. bgm() accepts it in both cases rather than stopping with an error, and mentions this when you set the argument yourself and verbose = TRUE.
Sampling with the Gibbs sampler
Besides the default NUTS sampler and adaptive Metropolis, all-continuous data admit a third option: update_method = "gibbs". It updates the precision matrix row by row from exact conjugate conditional distributions, and — under edge selection — updates the graph with a full-conditional birth/death step between parameter updates. It requires a Normal or Cauchy interaction prior and a Gamma-family prior on the precision diagonal (the defaults qualify).
Because every draw comes from an exact conditional, there are no step sizes or acceptance targets to tune — target_accept is ignored — and the warmup is spent settling rather than adapting: with edge selection, the first 15% of warmup runs the full model so the precision matrix settles, and selection is active for the remaining 85%. See GGM Internals for the update formulas.
Paired with precision_graph_prior = "joint", this sampler is the fast route for large models — two hundred variables in minutes rather than hours. Large Gaussian Graphical Models shows what the pair costs and buys, including a measured comparison with the neighboring CRAN packages.
Missing data
Missing values can be handled via listwise deletion (the default) or imputed within the Gibbs sampler by setting na_action = "impute". For continuous data, missing entries are drawn from their conditional normal distribution at each MCMC iteration. See Missing Data for details on both approaches and recommendations.
Fitting a GGM
To fit a GGM, set variable_type = "continuous":
fit = bgm(x = data, variable_type = "continuous", seed = 1234)The summary() and coef() methods work the same as for ordinal and mixed models. coef(fit)$pairwise returns the posterior means of the partial associations \(\omega_{ij}\), and coef(fit)$indicator returns the posterior inclusion probabilities.
To obtain model-specific parameterizations, use the dedicated extractors:
extract_precision(fit)— returns the full precision matrix \(\boldsymbol{\Theta}\).extract_partial_correlations(fit)— returns the partial correlation matrix.