---
title: "Plotting lag-sequential models"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Plotting lag-sequential models}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", message = FALSE,
                      warning = FALSE, dpi = 130, fig.width = 6.8,
                      fig.height = 5.4, out.width = "100%",
                      fig.align = "center")
set.seed(2026)
library(lagseq)
options(digits = 3)
has <- function(p) requireNamespace(p, quietly = TRUE)
fit <- lsa(engagement)
```

`lagseq` has one plotting verb, `plot(fit, type = )`, plus dedicated
plots for the resampling and comparison results. Every view of a fit
displays the same quantity -- the adjusted residual, the departure of a
transition from chance -- and differs only in geometry. The single fit
used throughout is the bundled `engagement` data (138 students, three
engagement states).

| view | call | backend |
|---|---|---|
| residual heatmap | `plot(fit)` | ggplot2 |
| residual network | `plot(fit, type = "network")` | cograph |
| transition (TNA) network | `plot(fit, type = "network", weights = "prob")` | tna |
| chord diagram | `plot(fit, type = "chord")` | cograph |
| polar sunburst | `plot(fit, type = "sunburst")` | ggplot2 |
| uncertainty forest | `plot(bootstrap_lsa(fit))`, `plot(certainty_lsa(fit))` | ggplot2 |
| group barrel | `plot(compare_lsa(g))`, `plot(bayes_compare_lsa(g))` | ggplot2 |

**Colour.** Two conventions are used, by purpose. The residual *heatmap*,
*chord*, and *sunburst* use the residual diverging scale (warm = over-
represented, cool = avoided). The *network* and *group comparison* plots
follow the wider TNA / Nestimate convention -- **blue = more than chance**,
**red = less** (avoided edges dashed) -- so an `lsa` network reads like any
other transition network in that ecosystem.

# Residual heatmap

The default. Rows are the current state, columns the next; colour is the
adjusted residual. `which` selects the matrix.

```{r heatmap, eval = has("ggplot2")}
plot(fit)                       # adjusted residuals (default)
plot(fit, which = "prob")       # transition probabilities
```

# Residual network

The same residuals as a directed graph: **blue** edges (solid) are
over-represented, **red** edges (dashed, with a soft halo) are avoided.

```{r net-residual, eval = has("cograph")}
plot(fit, type = "network")
```

# Transition network (a TNA model)

Weight the edges by probability and the view becomes the familiar
transition network: `lagseq` builds a `tna` object on the fly and renders
it with tna's styling. This shows what happens, where the residual network
shows what is surprising.

```{r net-transition, eval = has("tna")}
plot(fit, type = "network", weights = "prob")
```

# Chord and sunburst

A chord diagram of the transition flow, and a polar sunburst of each
state's outgoing distribution.

```{r chord, eval = has("cograph")}
plot(fit, type = "chord")
```

```{r sunburst, eval = has("ggplot2")}
plot(fit, type = "sunburst")                 # rose (default)
plot(fit, type = "sunburst", style = "wedge")  # frequency wedges
```

# Uncertainty forests

A fitted edge is one estimate; the forest shows its interval. Both the
resampling bootstrap and the analytic certainty plot as a circular forest
of per-edge intervals.

```{r forest, eval = has("ggplot2"), fig.height = 6.5}
plot(bootstrap_lsa(fit, R = 200))            # resampling CIs
plot(certainty_lsa(fit))                     # analytic credible intervals
```

# Group comparison

For a real grouping we use a long event log, `tna::group_regulation_long`,
and fit one model per achievement group.

```{r gfit, eval = has("tna")}
gfit <- lsa(tna::group_regulation_long, actor = "Actor",
            action = "Action", time = "Time", group = "Achiever")
```

`plot()` on a comparison draws a back-to-back **barrel**: each group's bar
runs to one side, the higher group's bar is bordered (darker for a larger
difference), and the centre chip carries the difference p-value.

```{r barrel, eval = has("tna") && has("ggplot2"), fig.height = 7}
cmp <- compare_lsa(gfit, R = 500)
plot(cmp)                                    # frequency-ranked rows
plot(cmp, rank = "effect")                   # surface avoided (red) edges
plot(cmp, style = "heatmap")                 # full-grid difference
```

The Bayesian comparison reuses the same barrel, with the credible
posterior in place of the permutation p-value.

```{r bayes-barrel, eval = has("tna") && has("ggplot2"), fig.height = 7}
plot(bayes_compare_lsa(gfit, seed = 1))
```

# Grouped fits

Every `type` works on a grouped fit, drawing one panel per group
(`combined = FALSE`, the default) or a single tiled figure
(`combined = TRUE`).

```{r grouped, eval = has("tna") && has("ggplot2")}
plot(gfit)                                   # one heatmap per group
```

# The worker functions

`plot(fit, type = )` is a front door over exported workers, each with its
own full argument list:

```{r workers, eval = FALSE}
plot_transitions(fit, weights = "residuals", significant = TRUE)  # network
plot_chords(fit, width = "prob")                                  # chord
plot_polar(fit, style = "wedge")                                  # sunburst
plot_forest(bootstrap_lsa(fit), n_top = 12)                       # forest
```

# In short

```r
plot(fit)                                   # residual heatmap (default)
plot(fit, type = "network")                 # residual network (blue = more)
plot(fit, type = "network", weights = "prob")  # transition network (a TNA model)
plot(fit, type = "chord"); plot(fit, type = "sunburst")
plot(bootstrap_lsa(fit)); plot(certainty_lsa(fit))   # uncertainty forests
plot(compare_lsa(gfit)); plot(bayes_compare_lsa(gfit))  # group barrels
plot(gfit)                                  # grouped: one panel per group
```