---
title: "Random Effects"
author: "Gilles Colling"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Random Effects}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  eval = identical(Sys.getenv("NOT_CRAN"), "true") || identical(Sys.getenv("IN_PKGDOWN"), "true")
)
```

## Installation

```{r, eval = FALSE}
install.packages("pak")                   # development version
pak::pak("gcol33/tulpaRatio")
```

## Overview

This vignette covers random effects in tulpaRatio, including random intercepts, random slopes, crossed effects, and nested effects. tulpaRatio uses [lme4](https://cran.r-project.org/package=lme4)-style formula syntax with extensions for correlated vs. uncorrelated slopes.

The examples use one simulated data set: counts and sampling effort from three plots at each of 16 sites in four regions, surveyed in six years. Region, site, plot and observation-level effects enter both the count and the effort; year and observer effects and the site-specific slopes of the covariates enter only the count.

```{r data}
library(tulpaRatio)

set.seed(42)
df <- expand.grid(plot = 1:3, site = 1:4, region = 1:4, year = 1:6)
df$site <- factor(paste(df$region, df$site, sep = "_"))
df$plot <- factor(df$plot)
df$region <- factor(df$region)
df$year <- factor(df$year)
n <- nrow(df)
df$depth <- rnorm(n)
df$temp <- rnorm(n)
df$x <- rnorm(n)
df$z <- rnorm(n)
df$observer <- factor(sample(1:6, n, replace = TRUE))

b_region <- rnorm(4, 0, 0.4)[df$region]
b_site <- rnorm(16, 0, 0.5)[df$site]
b_plot <- rnorm(48, 0, 0.3)[interaction(df$site, df$plot)]
b_year <- rnorm(6, 0, 0.4)[df$year]
b_observer <- rnorm(6, 0, 0.3)[df$observer]
slopes <- matrix(rnorm(16 * 5, 0, 0.2), nrow = 16)[df$site, ]
u <- rnorm(n, 0, 0.4)
df$effort <- rgamma(n, shape = 5, rate = 5 / exp(0.5 + b_region + b_site + b_plot + u))
df$count <- rpois(n, exp(
  1 + b_region + b_site + b_plot + b_year + b_observer + u +
    (0.4 + slopes[, 1]) * df$depth + (-0.2 + slopes[, 2]) * df$temp +
    (0.3 + slopes[, 3]) * df$x + (-0.2 + slopes[, 4]) * df$z +
    (0.1 + slopes[, 5]) * df$x * df$z
))

ctl <- list(seed = 1, adapt_delta = 0.95)
```

# Random Intercepts

## Single Grouping Factor

```{r}
# Random intercept for site
fit <- tratio(
  count | effort ~ depth + (1 | site),
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)
summary(fit)
```

This estimates a site-specific intercept $b_{\text{site}[i]}$ for each observation $i$:

$$
\log(\mu_i) = \beta_0 + \beta_1 \cdot \text{depth}_i + b_{\text{site}[i]}
$$

where $b_j \sim \text{Normal}(0, \sigma_{\text{site}}^2)$.

## Multiple (Crossed) Grouping Factors

When observations are classified by multiple independent grouping factors:

```{r}
# Crossed random effects: site and year
fit <- tratio(
  count | effort ~ depth + (1 | site) + (1 | year),
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)
```

This estimates:
- $b^{\text{site}}_j \sim \text{Normal}(0, \sigma_{\text{site}}^2)$
- $b^{\text{year}}_k \sim \text{Normal}(0, \sigma_{\text{year}}^2)$

The linear predictor becomes:

$$
\log(\mu_i) = \beta_0 + \beta_1 \cdot \text{depth}_i + b^{\text{site}}_{\text{site}[i]} + b^{\text{year}}_{\text{year}[i]}
$$

# Nested Random Effects

When groups are hierarchically nested (plots within sites):

```{r}
# Nested: plots within sites
fit_nested <- tratio(
  count | effort ~ depth + (1 | site/plot),
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)
```

This expands to:

```{r}
# Equivalent explicit form
fit <- tratio(
  count | effort ~ depth + (1 | site) + (1 | site:plot),
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)
```

The `site:plot` term groups by each plot-within-site combination.

### Deep Nesting

For deeper hierarchies:

```{r}
# region/site/plot expands to:
# (1 | region) + (1 | region:site) + (1 | region:site:plot)
# Every region, site and plot holds several observations, so each effect
# is well informed by the data and the centred parameterisation suits it
fit <- tratio(
  count | effort ~ x + (1 | region/site/plot),
  data = df,
  family = ratiod_poisson_gamma(),
  control = modifyList(ctl, list(re_param = "centered"))
)
```

# Random Slopes

Random slopes allow the effect of a predictor to vary by group.

## Correlated Intercept and Slope

Use single `|` for correlated random effects:

```{r}
# Intercept and depth effect vary by site, with correlation
fit_slope <- tratio(
  count | effort ~ depth + (1 + depth | site),
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)
```

This estimates:

$$
\begin{pmatrix} b_{0j} \\ b_{1j} \end{pmatrix} \sim \text{MVN}\left(\begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} \sigma_0^2 & \rho\sigma_0\sigma_1 \\ \rho\sigma_0\sigma_1 & \sigma_1^2 \end{pmatrix}\right)
$$

The correlation $\rho$ captures whether sites with higher intercepts also tend to have steeper slopes. The summary reports it as `cor_re[1,intercept:depth]` beside the two standard deviations:

```{r}
coefs <- coef(summary(fit_slope))
coefs[grepl("^(sigma_re|cor_re)", coefs$parameter), c("parameter", "mean", "sd")]
```

## Uncorrelated Intercept and Slope

Use double `||` for independent random effects:

```{r}
# Intercept and depth effect vary by site, independently
fit <- tratio(
  count | effort ~ depth + (1 + depth || site),
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)
```

This is computationally simpler and appropriate when you don't expect correlation between intercepts and slopes.

## Multiple Slopes

```{r}
# Multiple correlated slopes
fit <- tratio(
  count | effort ~ depth + temp + (1 + depth + temp | site),
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)

# Multiple uncorrelated slopes
fit <- tratio(
  count | effort ~ depth + temp + (1 + depth + temp || site),
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)
```

# Formula Expansion

tulpaRatio supports lme4-style formula expansion for random slopes.

## Interactions

```{r}
# x*z expands to x + z + x:z
fit <- tratio(
  count | effort ~ x*z + (1 + x*z | site),
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)
# Random slopes for x, z, and their interaction
```

## Polynomials

```{r}
# Quadratic random slope
fit <- tratio(
  count | effort ~ x + I(x^2) + (1 + x + I(x^2) | site),
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)
```

## Parameter Naming

tulpaRatio cleans parameter names for readability:

| Formula Term | Parameter Name |
|--------------|----------------|
| `x:z` | `x_z` |

# Shared Random Effects

By default, tulpaRatio shares random effects between numerator and denominator when the same grouping factor appears in both. This is the core of the tulpaRatio philosophy (see `vignette("philosophy")`).

## Default: Shared Structure

```{r}
# Site effect shared between count and effort processes
fit <- tratio(
  count | effort ~ (1 | site),
  formula_num = ~ depth,
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)
```

The site effect enters both linear predictors identically:

$$
\log(\mu_{\text{count}}) = \beta_0 + \beta_1 \cdot \text{depth} + \mathbf{b_{\text{site}}}
$$
$$
\log(\mu_{\text{effort}}) = \alpha_0 + \mathbf{b_{\text{site}}}
$$

Shared effects cancel in the ratio, removing spurious variation from common unmeasured drivers.

## Explicit Shared Structure

```{r}
# Only site is shared, observer is process-specific
fit <- tratio(
  count | effort ~ (1 | site) + (1 | observer),
  shared = ~ (1 | site),
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)
```

## Independent Structure (Use with Caution)

```{r}
set.seed(7)
trials_df <- data.frame(
  block = factor(rep(1:8, each = 6)),
  treatment = factor(rep(c("control", "treated"), 24)),
  trials = 20
)
b_block <- rnorm(8, 0, 0.5)[trials_df$block]
trials_df$successes <- rbinom(
  48, size = 20,
  prob = plogis(-0.5 + 0.8 * (trials_df$treatment == "treated") + b_block)
)

# Request independence (triggers warning)
fit <- tratio(
  successes | trials ~ treatment + (1 | block),
  shared = ~ 0,
  data = trials_df,
  family = ratiod_binomial(),
  control = ctl
)
```

The warning is issued every time `shared = ~ 0` is requested, because independence between the two processes is an assumption the data rarely support.

# Extracting Random Effects

Random-effect draws are named `re[term,group]`, with terms numbered in formula order (for `site/plot`, term 1 is `site` and term 2 is `site:plot`) and a third index naming the component of a random slope. `as_draws()` returns them in the [posterior](https://cran.r-project.org/package=posterior) format.

```{r}
# Posterior summaries for random effects
re <- posterior::subset_draws(as_draws(fit_nested), variable = "re")
head(posterior::summarise_draws(re))

# For random slopes: intercept and depth components of the first site
re_slope <- posterior::subset_draws(
  as_draws(fit_slope),
  variable = c("re[1,1,intercept]", "re[1,1,depth]")
)
posterior::summarise_draws(re_slope)
```

# Prior Specification

## Default PC Priors

tulpaRatio uses penalized complexity (PC) priors that favor simpler models:

```{r}
# Default: P(sigma > 1) = 0.01
fit <- tratio(
  count | effort ~ depth + (1 | site),
  data = df,
  family = ratiod_poisson_gamma(),
  priors = ratiod_priors(),
  control = ctl
)
```

## Custom Priors

```{r}
# Expect smaller random effects
fit <- tratio(
  count | effort ~ depth + (1 | site),
  data = df,
  family = ratiod_poisson_gamma(),
  priors = ratiod_priors(
    sigma = prior_half_cauchy(0.5)  # Half-Cauchy with scale 0.5
  ),
  control = ctl
)
```

# Latent Factors for Unmeasured Confounders

When you suspect unmeasured confounders are driving both numerator and denominator but don't have a natural grouping structure, latent factors provide an alternative to random effects.

## Single Latent Factor

```{r}
# Capture shared unmeasured variation
# Observation-level factors mix slowly, so the chains run longer
fit_latent <- tratio(
  count | effort ~ (1 | site),
  formula_num = ~ depth,
  data = df,
  family = ratiod_poisson_gamma(),
  latent = latent_factor(n_factors = 1),
  control = list(seed = 1)
)
head(latent_factors(fit_latent))
```

Latent factors are observation-level scores on which each process loads with its own loading:

$$
\log(\mu_{\text{count}}) = \beta_0 + \beta_1 \cdot \text{depth} + b_{\text{site}} + \lambda_1 f_i
$$
$$
\log(\mu_{\text{effort}}) = \alpha_0 + b_{\text{site}} + \lambda_2 f_i
$$

where $f_i \sim \text{Normal}(0, 1)$ with a sum-to-zero constraint and $\lambda_1 > 0$. The factor adds the covariance $\lambda_1 \lambda_2$ between an observation's two linear predictors and moves the ratio by $(\lambda_1 - \lambda_2) f_i$.

## Two Latent Factors

Two factors with a lower-triangular loading matrix give the two linear predictors an unrestricted residual covariance. A third factor would have nothing left to identify, so `latent_factor()` stops at two (one when `shared = FALSE`).

```{r, eval = FALSE}
fit <- tratio(
  count | effort ~ depth,
  data = df,
  family = ratiod_poisson_gamma(),
  latent = latent_factor(n_factors = 2),
  control = list(seed = 1)
)
```

## Latent Factors vs Random Effects

| Feature | Random Effects | Latent Factors |
|---------|---------------|----------------|
| Grouping | Required (site, year, etc.) | None needed |
| Level | Group-level | Observation-level |
| Use case | Known grouping structure | Residual confounding |
| Dimensionality | 1 per group | K factors × N observations |

Use latent factors when:
- You suspect unmeasured confounders but don't know the grouping structure
- Residual variation after random effects is still correlated
- You want to capture observation-level shared structure

## Prior Specification

```{r}
# Stronger regularization on latent factors
fit <- tratio(
  count | effort ~ depth,
  data = df,
  family = ratiod_poisson_gamma(),
  latent = latent_factor(
    n_factors = 1,
    prior = prior_pc(U = 0.5, alpha = 0.01)  # P(loading > 0.5) = 0.01
  ),
  control = list(seed = 1)
)
```

## See Also

- `vignette("getting-started")` - Basic usage
- `vignette("philosophy")` - Why shared structure matters
- `vignette("spatial-temporal")` - Spatial and temporal extensions