---
title: "Bayesian SSM Analysis"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Bayesian SSM Analysis}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
set.seed(12345)
library(ggplot2)
# Output that follows the console width wraps at 77 columns, so with its "#> "
# prefix it fits in 80. Text printed with cat() does not wrap.
options(width = 77)
```

```{r setup}
library(circumplex)
```

**Level:** Advanced. Read "Intermediate SSM Analysis" first.

::: {.alert .alert-warning role="alert"}
**Not yet peer reviewed.** The Bayesian SSM recipe this page teaches is the
package's own proposal. Its authors have not yet published it in a
peer-reviewed venue. Read it as a research tool, and state that status when
you report results from it.
:::

## 1. Overview

This vignette estimates the SSM parameters with a Bayesian model and
summarizes its posterior draws. Section 2, "Why a Bayesian SSM?", motivates
the approach. Section 3, "The cosine model as a linear regression", rewrites
the cosine model so that a regression can fit it. Section 4, "Fitting the
model with brms", fits it to the `jz2017` octant scores. Section 5, "From
posterior draws to SSM summaries", turns the draws into SSM summaries with
`ssm_draws()`. Section 6, "The induced prior on amplitude", shows what the
priors on $x$ and $y$ imply for the amplitude. The Wrap-up names related
functions and the next page to read.

## 2. Why a Bayesian SSM?

The Structural Summary Method (SSM) describes a circumplex profile with an
elevation $e$, an amplitude $a$, and a displacement $d$.
`vignette("introduction-to-ssm-analysis")` defines these three parameters.
The package's `ssm_analyze()` estimates them with bootstrap or Monte Carlo
confidence intervals.

A Bayesian alternative is attractive when you want prior information or
hierarchical structure. Prior information is what you know about the
parameters before you see the data. An example of hierarchical structure is partial
pooling across persons or groups, where each person's or group's estimate is
pulled toward the overall estimate. A Bayesian alternative is also attractive
when you want full posterior distributions for derived quantities, such as
amplitude and displacement. A posterior distribution combines the prior and
the data into one distribution over the parameters.

The division of labor is deliberate. A general-purpose Bayesian package, such
as **brms**, fits the model. Then `ssm_draws()` converts the resulting
posterior draws into SSM parameter draws. Posterior draws are samples
from the posterior distribution. `ssm_draws()` summarizes them with the same
circular statistics (statistics for angles) that the rest of the package uses. Displacement is an
angle, so its posterior needs circular treatment. In particular, an interval
that straddles the 0°/360° boundary must wrap across 0° rather than invert.
An inverted interval runs the long way around the circle. `ssm_draws()`
handles that by construction.

## 3. The cosine model as a linear regression

The SSM's cosine model for a profile of scores $S_j$ observed at scale
angles $\theta_j$ is

$$S_j = e + a \cos(\theta_j - d).$$

Expanding the cosine of a difference gives

$$S_j = e + \underbrace{a \cos d}_{x} \cos \theta_j +
\underbrace{a \sin d}_{y} \sin \theta_j,$$

a linear regression of the scores on $\cos \theta_j$ and $\sin \theta_j$.
The intercept is the elevation $e$, the *cosine* coefficient is $x$, and the
*sine* coefficient is $y$. The amplitude $a$ and displacement $d$ are
recovered by

$$a = \sqrt{x^2 + y^2}, \qquad d = \operatorname{atan2}(y,\, x),$$

with the displacement wrapped into $[0°, 360°)$.

Note the argument order: `atan2(y, x)` takes the **sine coefficient
first**. Swapping the arguments is a classic silent error. It returns a
valid-looking angle that is wrong for almost every profile. The example below
pins the convention with a profile whose displacement is known to be 90°.
The swapped call would return 0° instead:

```{r known-direction}
# Truth: e = 1, a = 2, d = 90 degrees
theta <- as.numeric(octants()) * pi / 180
scores <- 1 + 2 * cos(theta - pi / 2)
fit <- lm(scores ~ cos(theta) + sin(theta))
x_hat <- coef(fit)[["cos(theta)"]]
y_hat <- coef(fit)[["sin(theta)"]]
d_hat <- atan2(y_hat, x_hat) * 180 / pi
round(c(x = x_hat, y = y_hat, d = d_hat), 6)
```

The printed `d` is 90. A hidden check stops the vignette build unless
`atan2(y_hat, x_hat)` gives 90° and the swapped `atan2(x_hat, y_hat)` does not.

```{r known-direction-check, echo = FALSE}
stopifnot(isTRUE(all.equal(d_hat, 90)))          # atan2(y, x): correct
stopifnot(!isTRUE(all.equal(atan2(x_hat, y_hat) * 180 / pi, 90)))  # swapped
```

## 4. Fitting the model with brms

We model raw octant scores from the `jz2017` data. Octant scores are scores
on eight scales placed 45° apart around the circle. The data are in long
format, with one row per person-scale observation. A random intercept per
person gives each person their own baseline level. It is included to absorb
the dependence among a person's eight scores. The fixed effects are the
coefficients shared by all persons (intercept, cosine coefficient, sine
coefficient). They are then the group-level $(e, x, y)$. A random subsample
of 200 persons, drawn with a fixed seed, keeps the example light.

The long data frame is `dat`, with 1600 rows: one per person and octant scale.
Its columns are the person `id`, `cos_theta` and `sin_theta` (the cosine and
sine of the scale's angle from `octants()`), and the `score`. (The code that
draws the subsample and builds `dat` is omitted. It builds the same data as
the seeded script named below.) Its first rows are:

```{r data-prep, echo = FALSE}
data("jz2017")
set.seed(12345)
sub <- jz2017[sample(nrow(jz2017), 200), ]
scales <- c("PA", "BC", "DE", "FG", "HI", "JK", "LM", "NO")
theta <- as.numeric(octants()) * pi / 180
dat <- data.frame(
  id = rep(seq_len(200), times = length(scales)),
  cos_theta = rep(cos(theta), each = 200),
  sin_theta = rep(sin(theta), each = 200),
  score = unlist(sub[scales], use.names = FALSE)
)
```

```{r data-head}
head(dat)
```

**brms** fits its models with Stan, and Bayesian sampling with Stan requires a
working Stan toolchain. So the model below is
not re-fitted when this vignette is rebuilt. Its posterior draws were
generated once by the seeded script `data-raw/bayesian_ssm_draws.R`, and
they ship with the package. The `normal(0, 1)` prior on the regression
coefficients is a deliberate modeling choice. Section 6 examines its
consequences for the amplitude.

```{r brms-fit, eval = FALSE}
library(brms)
bfit <- brm(
  score ~ cos_theta + sin_theta + (1 | id),
  data = dat,
  prior = set_prior("normal(0, 1)", class = "b"),
  chains = 4, iter = 2000, cores = 4, seed = 12345
)
draws <- as.matrix(bfit,
                   variable = c("b_Intercept", "b_cos_theta", "b_sin_theta"))
```

## 5. From posterior draws to SSM summaries

The draws form a matrix with one row per posterior draw and three columns
interpreted **in column order** as $(e, x, y)$. The vignette reads the shipped
draws into the matrix `draws`. (The code that reads the file is omitted.)

```{r load-draws, echo = FALSE}
draws <- readRDS("bayesian_ssm_draws.rds")
```

```{r draws-shape}
dim(draws)
head(round(draws, 3))
```

`ssm_draws()` accepts two draw shapes. The first is three-column
*parameter* draws like these. The second is *profile* draws: posterior draws
of a whole profile of scale scores, with one column per scale and `angles`
supplied. A three-column matrix without angles is
ambiguous, because it could also be profile draws from a three-scale
instrument. So the shape must be stated explicitly with
`type = "parameters"`:

```{r adapter}
res <- ssm_draws(draws, type = "parameters")
summary(res)
```

Each posterior draw of $(x, y)$ was converted to a draw of $(a, d)$. Because
$d$ is an angle, the displacement's credible interval comes from circular
quantiles. A credible interval is the Bayesian counterpart of a confidence
interval. The circular mean is the direction of the average of the angles
taken as unit vectors. The draws are centered on their circular mean and
unwrapped (laid on a line with no jump near the mean). Then the quantiles are
taken, and the bounds are re-wrapped. The interval never comes from naive
linear quantiles that would misbehave near 0°/360°.

Point estimates are posterior medians for the linear parameters: elevation,
$x$, $y$, amplitude, and fit where it is available. The amplitude posterior is right-skewed, so a mean
would overstate it. For displacement, the point estimate is the circular
mean.

One caveat to keep in mind: these marginal summaries are not jointly
coherent. A marginal summary describes one parameter's draws on its own, so
the summaries need not agree with each other. The reported amplitude is the median of the amplitude draws. It is
not $\sqrt{x^2 + y^2}$ evaluated at the reported $x$ and $y$. Likewise, the
reported displacement is not the direction of the reported $(x, y)$. Each is
the honest marginal summary of its own posterior.

Model fit ($R^2$) is not reported here, because parameter draws carry no
profile to measure fit against. Profile draws (the second shape) do yield
fit draws from `ssm_draws()`.

## 6. The induced prior on amplitude

Independent priors on $x$ and $y$ do not induce a flat prior on the
structural parameters. The induced prior is the prior on $a$ and $d$ that
the priors on $x$ and $y$ imply. With $x, y \sim \mathrm{Normal}(0, 1)$, the implied
prior on $a = \sqrt{x^2 + y^2}$ is Rayleigh-shaped: zero density at $a = 0$,
then a single peak, with a right skew. Its mass is pushed *away* from $a = 0$.
But the implied prior on $d$ is uniform. A prior-predictive simulation makes
the amplitude prior visible. It draws 10000 values each of $x$ and $y$ from
$\mathrm{Normal}(0, 1)$ and computes $a = \sqrt{x^2 + y^2}$ for each pair. The
histogram shows those draws of $a$. The two printed numbers are their median
and the share of them below 0.1. (The simulation and plotting code is omitted.)

```{r induced-prior, echo = FALSE, fig.width = 5, fig.height = 3}
x_prior <- rnorm(10000, 0, 1)
y_prior <- rnorm(10000, 0, 1)
a_prior <- sqrt(x_prior^2 + y_prior^2)
ggplot(data.frame(a = a_prior), aes(x = a)) +
  geom_histogram(bins = 60, fill = "grey35") +
  labs(
    x = "Amplitude implied by the priors on x and y",
    y = "Prior draws",
    title = "Rayleigh-shaped induced prior on amplitude"
  ) +
  theme_minimal()
round(c(prior_median = median(a_prior), prior_mass_below_0.1 =
          mean(a_prior < 0.1)), 3)
```

The Rayleigh shape is not a defect. It is a modeling choice to be aware of:
the prior mildly disfavors exactly-flat profiles ($a = 0$). If your application needs prior
mass concentrated near $a = 0$, place priors on $(a, d)$ directly in a
custom Stan model instead. Its posterior draws can still be summarized
here. Convert them to $(x, y) = (a \cos d, a \sin d)$, or pass profile
draws.

## Wrap-up

* Per-person descriptive SSM parameters (no pooling):
  `ssm_parameters_id()` and its `summary()` method.
* Hierarchical pooling across persons or groups is exactly what the
  brms/Stan route is for. Fit the model you need, and feed the posterior
  draws (parameter or profile shape) back through `ssm_draws()`.
* Frequentist inference on group profiles and contrasts: `ssm_analyze()`.

The next page to read is "Growth Models on SSM Parameters". It fits a growth
model to SSM parameters measured at several waves and turns its fixed
effects into amplitude and displacement curves with `ssm_draws()`.