---
title: "Getting Started with bgms"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Getting Started with bgms}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Introduction

The **bgms** package implements Bayesian methods for analyzing graphical
models. It supports three variable types:

- **ordinal** (including binary) — Markov random field (MRF) models,
- **Blume–Capel** — ordinal MRF with a reference category,
- **continuous** — Gaussian graphical models (GGM).

The package estimates main effects and pairwise interactions, with optional
Bayesian edge selection via spike-and-slab priors. It provides two main
entry points:

- `bgm()` for one-sample designs (single network),
- `bgmCompare()` for independent-sample designs (group comparisons).

This vignette walks through the basic workflow for ordinal data. For
continuous data, set `variable_type = "continuous"` in `bgm()` to fit a
Gaussian graphical model.


# Wenchuan dataset

The dataset `Wenchuan` contains responses from survivors of the 2008 Wenchuan
earthquake on posttraumatic stress items. Here, we analyze a subset of the
first five items as a demonstration.


```{r}
library(bgms)

# Analyse a subset of the Wenchuan dataset
?Wenchuan
data = Wenchuan[, 1:5]
head(data)
```

# Fitting a model

The main entry point is `bgm()` for single-group models and 
`bgmCompare()` for multiple-group comparisons.

```{r, eval=FALSE}
fit = bgm(data, seed = 1234)
```

```{r, include=FALSE}
fit = bgm(data, seed = 1234, chains = 2, display_progress = "none", verbose = FALSE)
```

# Posterior summaries

```{r}
summary(fit)
```

You can also access posterior means or inclusion probabilities directly:

```{r}
coef(fit)
```

# Network plot

To visualize the network structure, we threshold the posterior inclusion probabilities at 0.5 and plot the resulting adjacency matrix.

```{r, fig.width= 7, fig.height= 7}
library(qgraph)

median_probability_network = coef(fit)$pairwise
median_probability_network[coef(fit)$indicator < 0.5] = 0.0

qgraph(median_probability_network,
  theme = "TeamFortress",
  maximum = 1,
  fade = FALSE,
  color = c("#f0ae0e"), vsize = 10, repulsion = .9,
  label.cex = 1, label.scale = "FALSE",
  labels = colnames(data)
)
```

# Continuous data (GGM)

For continuous variables, `bgm()` fits a Gaussian graphical model when
`variable_type = "continuous"`. The workflow is the same:

```{r, eval=FALSE}
fit_ggm = bgm(continuous_data, variable_type = "continuous", seed = 1234)
summary(fit_ggm)
```

The pairwise effects are unstandardized partial associations derived from
the off-diagonal entries of the precision matrix (stored as `-0.5 * K_ij`).
They are not standardized partial correlations. Use
`extract_partial_correlations()` to convert them to actual partial
correlations, or `extract_precision()` to obtain the precision matrix.
Missing values can be imputed during sampling with `na_action = "impute"`.

# Next steps

- For comparing groups, see `?bgmCompare` or the *Model Comparison* vignette.
- For diagnostics and convergence checks, see the *Diagnostics* vignette.