residuals() returns response residuals or Pearson-style residuals from a
drmTMB fit.
Arguments
- object
A
drmTMBfit.- type
Residual type:
"response","pearson", or"quantile".- ...
Reserved for future residual options; for
type = "quantile", forwarded todrm_quantile_residuals()(seed,nsim,response).
Value
A numeric vector for univariate models, or a two-column matrix for
bivariate Gaussian models. For type = "quantile" with nsim > 1, an
n-by-nsim matrix (see drm_quantile_residuals()).
Details
For univariate Gaussian models, response residuals are y - mu. Pearson
residuals divide by the fitted observation standard deviation. If a dense
known sampling covariance was used, Pearson residuals are whitened by the
fitted total observation covariance.
For lognormal models, response residuals are y - fitted_mean. Pearson
residuals are computed on the log-response scale as (log(y) - mu) / sigma.
For Gamma models, response residuals are y - mu and Pearson residuals
divide by the fitted Gamma standard deviation mu * sigma, where sigma is
the coefficient of variation. For Tweedie models, response residuals are
y - mu and Pearson residuals divide by
sqrt(sigma^2 * mu^nu). For zero-one beta models, response residuals
are observed proportions minus the unconditional fitted mean, including exact
zero-one boundary mass. For beta-binomial models, response residuals are
observed success proportions minus fitted mu, and Pearson residuals
divide by the fitted beta-binomial proportion standard deviation. For
binomial models, response residuals are observed success proportions minus
fitted event probability, and Pearson residuals divide by
sqrt(mu * (1 - mu) / trials). For
cumulative-logit ordinal models, response residuals are the observed
ordered-category score minus the fitted expected score, and Pearson
residuals divide by the fitted category-score standard deviation. For
Poisson models, response residuals are y - mu and Pearson residuals divide
by sqrt(mu). For zero-inflated
Poisson models, response residuals are y - (1 - zi) * mu, and Pearson
residuals divide by sqrt((1 - zi) * mu * (1 + zi * mu)). For
negative-binomial 2 models, Pearson residuals divide by
sqrt(mu + sigma^2 * mu^2). For zero-truncated NB2 models, response
residuals are y - mu / (1 - Pr_NB2(0)), and Pearson residuals divide by
the conditional positive-count standard deviation. For hurdle NB2 models,
response residuals are y - (1 - hu) * mu / (1 - Pr_NB2(0)), and Pearson
residuals divide by the unconditional standard deviation implied by the
hurdle-zero mixture. For zero-inflated NB2 models, response residuals are
y - (1 - zi) * mu, and Pearson residuals divide by the unconditional
standard deviation implied by the structural-zero mixture.
For bivariate Gaussian models, response residuals are returned as a
two-column matrix. Pearson residuals are standardized and whitened using the
fitted residual sigma1, sigma2, and rho12, or using the full row-paired
observation covariance when a dense bivariate known V was supplied.
type = "quantile" returns Dunn-Smyth randomized quantile residuals from
drm_quantile_residuals(): qnorm(F(y; theta_hat)) at the fitted,
fixed-effect distributional parameters. As of DO-T3 batch D, all 18 fitted
model_type values have a promoted entry in drm_family_dpq(), so
type = "quantile" is available for every family, including bivariate
biv_gaussian (where response – 1 or 2 – is REQUIRED and selects
which response's MARGINAL residuals to compute; omitting it errors
clearly, as does supplying it for a univariate fit; see
fitted_distribution()). Spike-status families (feasibility-only, not yet
DG2/DG3-verified) emit a one-time warning that the residual is
exploratory, not DG-verified; none currently have this status. See
drm_quantile_residuals() for the fixed-effect-only adequacy caveat: for
random-effect or structured fits, these residuals are conditional on the
fixed-effect prediction, not marginal, so a departure (or its absence) is
evidence about fixed-effect adequacy only – never a general validity
claim. Pass seed, nsim, and/or response through ... to
drm_quantile_residuals().
What this detects – and does not. Evidence: a 400-seed gated
simulation campaign across all 18 fitted families (tweedie: 99 of 400
seeds locally, 66/99 dispersion-arm non-convergence, full run deferred to
Totoro)
(docs/dev-log/simulation-artifacts/2026-07-12-dg3-power-arm-gated/).
Under a correctly specified fixed-effect model, type-I error stays at or
below the nominal rate (0.0025-0.025 at alpha = 0.05 across families); the
underlying KS+PIT statistic is conservative, so power is understated, not
overstated. type = "quantile" detects distributional shape and atom
mis-specification that a family cannot reabsorb through its own free
parameters – heavy tails fit as Gaussian, overdispersion or
zero-inflation ignored by a family with no free dispersion parameter,
truncation ignored, a missing zero/one atom – with gated power >= 0.8 at
n = 300-400 per arm (commonly 0.9-1.0). It has a genuine structural
blind spot, not a bug: a mis-specification that a fitted family's own
free nuisance/dispersion/inflation parameter absorbs leaves the
fitted-model residual marginally N(0,1) and is not detectable here –
for example heteroscedasticity absorbed by Student-t nu (power 0.035 at
n = 300, versus 1.0 for the same heteroscedasticity under Gaussian, which
has no absorbing parameter), missing zero-inflation absorbed by nbinom2
sigma (power 0.035, versus 0.9625 under Poisson), and zero-inflation/
hurdle/zero-one-inflation mechanism mis-specification (a constant
inflation probability fit when it truly varies with a covariate), which
splits into two patterns under the n-ladder (tested to n = 3000): for
hurdle_nbinom2/zero_one_beta, power stays flat at or below about 0.01
at every n – a genuine structural blind spot; for
zi_nbinom2/zi_poisson, power rises with n (to about 0.11/0.06 at
n = 3000), so the marginal is not identical under the mechanism mis-spec,
but power stays far below the >= 0.8 detectable benchmark even at
n = 3000, so it remains impractical to detect at realistic sample sizes.
A mean-structure diagnostic, not this one, is what catches
an absorbed mis-specification. Separately, gamma-vs-lognormal
wrong-family detection is sample-size limited rather than structurally
blind: power rises from about 0.19 at n = 300 to 0.79 at n = 1000 and 1.0
at n = 3000, so that specific mis-specification needs n well above 1000 to
be reliably caught. A distributional-output/adequacy tick on a family
never changes or implies anything about that family's own inference-tier
status; see tests/testthat/test-dg-firewall.R.
Examples
dat <- data.frame(y = c(0.2, 0.5, 1.1, 1.4), x = c(-1, 0, 1, 2))
fit <- drmTMB(bf(y ~ x, sigma ~ 1), data = dat)
residuals(fit)
#> [1] 0.03000001 -0.09000002 0.08999995 -0.03000008
residuals(fit, type = "pearson")
#> [1] 0.4472135 -1.3416404 1.3416394 -0.4472146
residuals(fit, type = "quantile")
#> [1] 0.4472135 -1.3416404 1.3416394 -0.4472146