---
title: "Response families"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Response families}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
# The graphics device canvas, not the ggplot theme, is what makes a figure's
# background opaque, so the transparent device is what lets the page colour show
# through. On the package website that matters twice over: the light page is
# warm off-white rather than pure white, and the dark page inverts the figure,
# where an opaque matte would read as a slab in either mode.
knitr::opts_chunk$set(
  collapse = FALSE, comment = "", fig.width = 6, fig.height = 3.4,
  dev.args = list(bg = "transparent")
)
# Console colour carries no meaning on a rendered page. pkgdown turns it on for
# its own build, and the escape sequences then reach the reader as literal text,
# so colour is switched off here for a plain vignette render and a site build
# alike. The fixed width keeps printed output inside the documentation column.
options(cli.num_colors = 1, cli.hyperlink = FALSE, crayon.enabled = FALSE,
        width = 80)
```

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

Behavioural outcomes are seldom Gaussian. pilotr maps the linear predictor
$\eta = \beta_0 + \sum_k \beta_k x_k + (\text{random effects})$ to an outcome through one of
eight response families. An important consequence is that the fixed intercept and effect live on
the family's own scale. That scale is the identity for the Gaussian and ex-Gaussian families,
the log scale for lognormal reaction times, reading times and counts, and the logit scale for
accuracy, ordinal and proportion outcomes. This vignette presents each family with
scale-appropriate parameters.

The helper below builds a two-group between-subjects design and reports the group means.

```{r}
demo <- function(family, intercept, effect, n = 4000, ...) {
  spec <- build_spec(list(
    name = family, seed = 1, design_kind = "between", n_subject = n,
    factor_name = "group", lev1 = "control", lev2 = "treatment",
    intercept = intercept, effect = effect, family = family,
    resp_name = "", ...))
  d <- simulate_design(spec)
  y <- d[[spec$response$name]]
  list(spec = spec, data = d, y = y, by_group = tapply(y, d$group, mean))
}

library(ggplot2)
fam_hist <- function(y, fill, title, xlab) {
  ggplot(data.frame(y = y), aes(y)) +
    geom_histogram(bins = 40, fill = fill, colour = NA) +
    labs(title = title, x = xlab, y = "count") +
    theme_minimal(base_size = 12) +
    # theme_minimal still paints a white plot.background over the transparent
    # device canvas, so both surfaces have to be cleared for the page colour to
    # reach the figure. The ink stays at its default, because the website
    # inverts the figure in dark mode, which turns the dark axis text light of
    # its own accord.
    theme(plot.background  = element_rect(fill = NA, colour = NA),
          panel.background = element_rect(fill = NA, colour = NA),
          panel.grid       = element_line(colour = "grey80"))
}
```

## Gaussian

This is the default, with $y = \eta + \varepsilon$ and residual standard deviation `sigma`. It
suits continuous, roughly symmetric outcomes such as ratings averaged over many trials or
standardised scores.

```{r}
g <- demo("gaussian", intercept = 100, effect = 5, sigma = 10)
round(g$by_group, 2)
fam_hist(g$y, "#2C6FB0", "Gaussian", "score")
```

## Shifted lognormal (reaction times)

Reaction times are right-skewed and bounded below. pilotr models them as
$y = \text{shift} + \exp(\eta + \varepsilon)$, with the effect on the log scale. An
intercept of 6 implies a typical RT near `exp(6)` ms above the shift.

```{r}
rt <- demo(
  "shifted_lognormal", intercept = 6, effect = 0.1, sigma = 0.3, shift = 200
)
round(rt$by_group, 1)
fam_hist(rt$y, "#B0402C", "Shifted lognormal (RT)", "RT (ms)")
```

## Lognormal (positive continuous)

The plain `lognormal` family is the shifted lognormal without the shift,
$y = \exp(\eta + \varepsilon)$, suited to positive continuous outcomes such as per-word reading
times. As with reaction times the effect is on the log scale.

```{r}
ln <- demo("lognormal", intercept = 6, effect = 0.1, sigma = 0.3)
round(ln$by_group, 1)
fam_hist(ln$y, "#7A4FB0", "Lognormal", "reading time (ms)")
```

## Bernoulli (accuracy)

Binary accuracy is modelled through a logit link. The intercept is the log-odds of a correct
response, and the effect is a log-odds difference between conditions.

```{r}
acc <- demo("bernoulli", intercept = 0, effect = 0.5)
round(acc$by_group, 3)   # P(correct) by group
```

## Poisson (counts)

Counts via a log link (e.g. number of fixations, errors or events). An intercept of 1.5
implies a base rate near `exp(1.5)`.

```{r}
cts <- demo("poisson", intercept = 1.5, effect = 0.3)
round(cts$by_group, 2)   # mean count by group
table(cts$y)[1:8]
```

## Ordinal (Likert)

Ordered categorical responses via a cumulative-logit model with user thresholds. The effect
shifts the latent distribution across the thresholds.

```{r}
ord <- build_spec(list(
  name = "likert", seed = 1, design_kind = "between", n_subject = 4000,
  factor_name = "group", lev1 = "control", lev2 = "treatment",
  intercept = 0, effect = 0.8, family = "ordinal", resp_name = "rating",
  thresholds = "-2, -0.6, 0.6, 2"))
r <- simulate_design(ord)
# category proportions by group
round(prop.table(table(r$group, r$rating), 1), 2)
```

## Beta (proportions)

Bounded proportions in (0, 1) are modelled through a mean–precision parameterisation. The mean
is `logit⁻¹(η)` and `phi` is the precision, with larger values giving a tighter distribution.

```{r}
bt <- demo("beta", intercept = 0, effect = 0.8, phi = 8)
round(bt$by_group, 3)   # mean proportion by group
fam_hist(bt$y, "#2E8B57", "Beta", "proportion")
```

## Choosing a family

The families above, together with the ex-Gaussian, cover the outcome types a behavioural study
usually produces. The table below sets each one against the scale its intercept and effect are
written on.

| Outcome | Family | Scale of the effect |
|---|---|---|
| Continuous, symmetric | `gaussian` | identity |
| Reaction time | `shifted_lognormal` | log |
| Positive continuous (e.g. reading time) | `lognormal` | log |
| Reaction time, on the response scale | `exgaussian` | identity |
| Accuracy (0/1) | `bernoulli` | logit |
| Counts | `poisson` | log |
| Likert / ordered categories | `ordinal` | logit (cumulative) |
| Proportions in (0, 1) | `beta` | logit (mean) |

These match the families that researchers fit in `lme4`, `glmmTMB` and `brms`, so a design
simulated here corresponds to the model that will later be fit. The point-and-click application
exposes six of the eight families. The full engine, including the plain lognormal, the
ex-Gaussian, continuous predictors, interactions, nesting and partial crossing, is reached by
writing the specification directly, as the last section of this vignette shows.

The `exgaussian` family, new in 0.3, takes `sigma` and `beta` and draws a normal plus an
exponential, mean-centred by subtracting the exponential's own mean so that $\eta$ remains the
mean of the response. That is brms's `exgaussian(mu, sigma, beta)` parameterisation, so a
specification and the model fitted to it agree on what the intercept means. A shifted lognormal
will not stand in for it, because `model_data()` logs the response back and leaves a symmetric
residual on the analysis scale. `build_spec()` does not cover the family either, so an
ex-Gaussian design has to be written out by hand.

## Writing the specification directly

The specification is a plain list, so richer designs than the builder covers can be assembled
by hand. Starting from a `build_spec()` result, the example below adds a continuous item-level
predictor, an interaction with the categorical effect, a per-subject item subset (partial
crossing) and an extra grouping factor that nests subjects.

```{r}
spec <- build_spec(list(
  name = "reading", seed = 1, design_kind = "within", include_items = TRUE,
  n_subject = 12, n_item = 24,
  factor_name = "condition", lev1 = "related", lev2 = "unrelated",
  intercept = 6, effect = 0.05,
  subj_int_sd = 0.12, subj_slope_sd = 0.04, subj_corr = 0.2,
  item_int_sd = 0.08, item_slope_sd = 0.02, item_corr = -0.1,
  family = "lognormal", resp_name = "RT", sigma = 0.25))

spec$predictors <- list(
  list(name = "freq", varies_by = "item", mean = 0, sd = 1)
)
# an interaction with the effect
spec$fixed$coefficients[["effect:freq"]] <- 0.02
# each subject sees 10 of the 24 items
spec$units$item$per_subject <- 10
# subjects nested in classes
spec$random$class <- list(over = "subject", n = 6, intercept_sd = 0.05)

head(simulate_design(spec))
```

The derived columns follow from the specification: `freq` is the continuous predictor and
`class` is the nesting factor. The auto-derived formula picks up the interaction as
`effect_freq` and the extra grouping factor as `(1 | class)`.

```{r}
model_formula(spec)
```

Ready-to-run specifications of this kind, including a continuous-predictor reading-time design,
a nested-clusters design and a partial-crossing design, ship in the repository's
[`spec/examples/`](https://github.com/pablobernabeu/pilotr/tree/main/spec/examples) directory,
and `load_spec()` reads any of them back.