---
title: "Distributional regression with drmTMB"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Distributional regression with drmTMB}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
if (!"package:drmTMB" %in% search()) {
  library(drmTMB)
}
```

<style>
@media (max-width: 575.98px) {
  main h1 {
    hyphens: none;
    overflow-wrap: normal;
    word-break: normal;
  }

  main table {
    display: block;
    max-width: 100%;
    overflow-x: auto;
  }
}
</style>

`drmTMB` fits distributional regression models for one or two responses. This
introduction explains the model components, fits a first model, and shows how
to interpret it. If you only need to find a function or tutorial, use the
shorter [Get started function map](https://itchyshin.github.io/drmTMB/articles/function-map-cheatsheet.html).

The central workflow is simple: write one formula for each parameter you want
to estimate, then check that the fitted model matches the symbolic model you had
in mind. The broader implemented map lives in
[What can I fit today?](model-map.html), and family choice lives in
[Choosing response families](distribution-families.html).

The parameter names stay consistent across the site. Location is a
family-specific centre or location parameter; for Gaussian models it is the
expected response, but it need not be the unconditional response mean. Scale
describes residual variability; shape describes distribution features beyond
location and scale; and coscale means residual correlation such as `rho12`
between two responses. A family may expose only a subset of these components.

For example, an applied user might ask: do mean trait values and residual
variability change with an environmental predictor, after accounting for
repeated measures from the same site or species? In `drmTMB`, that question is
written as one formula for the mean and one formula for the residual scale.

## Install the pre-CRAN release

`drmTMB` 0.7.0 is the first CRAN-targeted release. At the time this
documentation was built, it had not yet been accepted by CRAN. Install the
current development source from GitHub with `pak`:

```r
install.packages("pak")
pak::pak("itchyshin/drmTMB")
```

After CRAN accepts the package, install the released version with
`install.packages("drmTMB")`.

You need R 4.1.0 or newer and a working compiler toolchain because TMB models
are compiled during installation. If installation fails while compiling C++,
install the usual R build tools for your platform: Rtools on Windows, Xcode
Command Line Tools on macOS, or the R development toolchain on Linux.

The core runtime dependencies are installed automatically by `pak`: `cli`,
`Matrix`, `TMB`, and the compiled headers from `RcppEigen` and `TMB`. The
articles and development checks also use optional packages such as `glmmTMB`,
`lme4`, `MASS`, `metafor`, `knitr`, `rmarkdown`, `testthat`, `withr`, and
`pkgdown`.

## Fit your first model

Start with a Gaussian location-scale model when the response is continuous and
the scientific question is about both the expected value and predictability.
In the small example below, habitat and temperature affect mean growth, while
habitat also changes residual variation:

```{r first-model-data}
set.seed(13)
n <- 120
dat <- data.frame(
  habitat = factor(rep(c("forest", "grassland"), each = n / 2)),
  temperature = rnorm(n)
)
mu <- 1 + 0.6 * (dat$habitat == "grassland") + 0.4 * dat$temperature
sigma <- exp(-0.5 + 0.45 * (dat$habitat == "grassland"))
dat$growth <- rnorm(n, mean = mu, sd = sigma)
```

The fitted model uses one formula for `mu` and one formula for `sigma`:

```{r first-model-fit}
fit <- drmTMB(
  drm_formula(growth ~ habitat + temperature, sigma ~ habitat),
  family = gaussian(),
  data = dat
)

check_drm(fit)
```

Read the `sigma` coefficient as a log residual-SD contrast. Exponentiating it
gives an SD ratio; exponentiating twice the coefficient gives a residual
variance ratio:

```{r first-model-interpret}
sigma_habitat <- coef(fit, "sigma")["habitatgrassland"]
data.frame(
  residual_sd_ratio = exp(sigma_habitat),
  residual_variance_ratio = exp(2 * sigma_habitat)
)
```

For a fuller walkthrough of fitted means, residual SDs, and residual
variances, read [When variance carries signal, Part
1](location-scale.html). Continue to [Part
2](location-scale-scale.html) when predictors model a grouped or
phylogenetic random-effect SD through `sd()`.

## Learning path

Start with the question you want to answer, not with the full list of features.
The tutorials are arranged so that each article pairs symbolic equations, R
syntax, fitted output, and interpretation.

The opening rows below follow the Get started menu: check the reporting
boundary, choose a family, find the smallest useful function sequence, and
then check the fitted model. The remaining rows route common modelling and
diagnostic questions to the matching guide.

| If your question is... | Read this first | Main parameter |
| --- | --- | --- |
| Can I fit this model, report its point estimate or interval, and what should I try if I cannot? | [Can I fit and report this model?](capability-and-limits.html) | reporting scope and named fallback |
| Which family fits my response — continuous, count, proportion, or robust? | [Choosing response families](distribution-families.html) | family-specific `mu`, `sigma`, `nu`, `zi`, `hu` |
| I know my task but not which function or shortest workflow to use. | [Function map and cheat sheet](https://itchyshin.github.io/drmTMB/articles/function-map-cheatsheet.html) | task-to-function route and minimal workflow |
| After fitting, which checks and interval-status columns should I read? | [Checking and using fitted models](https://itchyshin.github.io/drmTMB/articles/model-workflow.html) | `check_drm()`, `profile_targets()`, `conf.status` |
| How do I profile a random-effect SD and know when not to trust the interval? | [First-week intervals: fit, profile, and boundary](first-week-intervals.html) | `profile_targets()`, `confint(method = "profile")`, `profile.boundary` |
| Am I modelling residual variation, group-level variation, or likelihood weights? | [Which scale are you modelling?](which-scale.html) | `sigma`, `sd(group)`, `weights` |
| Do group or phylogenetic covariance rows answer a bivariate question beyond residual coupling? | [What can I fit today?](model-map.html) | `corpairs()`, `corpair()`, `sd1(group, level = "phylogenetic")`, `sd2(group, level = "phylogenetic")` |
| Do predictors change the mean and residual variability of one response? | [When variance carries signal, Part 1](location-scale.html) | `mu`, `sigma` |
| Do predictors change residual variability and the SD of a grouped or phylogenetic random effect? | [When variance carries signal, Part 2](location-scale-scale.html) | `mu`, `sigma`, `sd(group)` |
| Do counts show extra-Poisson variation or structural zeros? | [Count abundance and extra zeros](count-nbinom2.html) | NB2 `mu`, `sigma`, `zi` |
| Are proportions successes out of trials, continuous rates inside `(0, 1)`, or continuous rates with structural exact boundaries? | [Proportions and success rates](proportion-beta-binomial.html) | beta-binomial, beta, and zero-one beta `mu`, `sigma`, `zoi`, `coi` |
| Does a predictor change residual coupling between two responses? | [Changing residual coupling with rho12](bivariate-coscale.html) | `rho12` |
| Do repeated individuals have correlated average responses after accounting for residual coupling? | [Changing residual coupling with rho12](bivariate-coscale.html) | `corpairs(level = "group")` |
| Do effect sizes have known sampling variances or covariance? (a specialist route, not a family choice for raw response data) | [Mean effects and residual heterogeneity in meta-analysis](meta-analysis.html) | `meta_V(V = V)`, `sigma` (implemented/tested; tier unregistered, with no interval or coverage claim) |
| Do species trait means remain similar after accounting for shared ancestry, do sites share coordinate-structured deviations, or do individuals have known relatedness? | [Structural dependence overview](structural-dependence.html) | `animal()`, `phylo()`, `spatial()`, `relmat()` |
| I have more than one candidate fit — which one should I report? | [Model selection with AIC and BIC](model-selection.html) | `AIC()`, `BIC()` |
| Is the fitted distribution itself adequate, not just the mean structure? | [Distributional outputs and adequacy](https://itchyshin.github.io/drmTMB/articles/distributional-outputs-and-adequacy.html) | `worm_plot()`, `qq_plot()`, `centile_chart()` |
| The optimizer did not converge cleanly — what do I check? | [Errors, warnings, and convergence](convergence.html) | `optimizer_convergence`, `fixed_gradient` |
| My fit is slow or memory-heavy — how do I keep fitting at this scale? | [Working with large data](large-data.html) | `keep_data`, `keep_model_frame`, `se = FALSE`, `se_group_sd` |

Here `nu` is a shape parameter, `zi` and `hu` are zero-process parameters, and
`sd(group)` refers to a group-level random-effect standard deviation rather than
the residual scale `sigma`. Use the model guide "What can I fit today?" when you
need the longer status map before choosing syntax.

For a first applied analysis, fit the simplest model that answers the question,
run `check_drm()`, and then read the coefficient table on the parameter scale
used by the model. For example, `sigma` coefficients are on a log scale in
Gaussian location-scale models, while `rho12(fit)` returns residual
correlations on the response scale.

For slope and variance-component questions, name the estimand before reporting
the number. A `mu` slope is an expected-response effect, a `sigma` slope is a
log residual-SD effect, a random-slope SD is among-group variation in a
reaction norm, and `sd(group) ~ x_group` is a model for the SD of a group-level
mean effect. Those four quantities can all involve a predictor, but they answer
different biological questions.

## Check before interpreting

After fitting any model, run `check_drm()` before interpreting the estimates:

```r
check_drm(fit)
```

The diagnostic table checks convergence, gradients, Hessian status, standard
errors, dropped rows, scale values, random-effect replication, and relevant
parameter boundaries. Inspect a `note`; resolve a `warning` or `error` before
treating estimates as stable. The [errors, warnings, and convergence
guide](convergence.html) explains what to try next.

The key rule is to keep correlation layers separate. A bivariate residual
`rho12` is within-observation coupling after the two means and residual SDs are
modelled. A random-effect correlation is a group-level quantity. Phylogenetic
rows from `corpairs(..., level = "phylogenetic")` are structured-effect
quantities, and the first coordinate-spatial mean-mean row is reported by
`corpairs(..., level = "spatial")` when matching labelled spatial terms are
fitted in `mu1` and `mu2`. The constant coordinate-spatial q=4 location-scale
block is also fitted when matching labelled spatial terms appear in all four
`mu1`/`mu2`/`sigma1`/`sigma2` endpoints. Predictor-dependent spatial
`corpair()` rows remain planned. Those layers should not be reported as if they
were the same estimand.

Use [Can I fit and report this model?](capability-and-limits.html) for the
current reporting boundary and named fallback; use the [model map](model-map.html)
when you need syntax detail. Once this first fit is checked, continue with
[Checking and using fitted models](https://itchyshin.github.io/drmTMB/articles/model-workflow.html) or choose a scientific
question from the Learning path table above.