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. Student-t
Pearson residuals use the fitted Student-t scale; this is not the marginal
residual standard deviation when nu > 2. Skew-normal Pearson residuals use
the public response-standard-deviation parameter. 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 beta models, response residuals are y - mu
and Pearson residuals divide by
sqrt(mu * (1 - mu) * sigma^2 / (1 + sigma^2)). 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. All 18 currently fitted
model_type values have status = "reference" 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()). Families whose distributional
diagnostics have not been evaluated by simulation emit a one-time warning
that the residual is exploratory; 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. The evidence used 400 seeds across
18 fitted families. Tweedie evidence is incomplete: 99 of 400 intended seeds
were run, and 66 of the 99 dispersion-arm fits did not converge. For the
completed family scenarios, 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 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 result for a family
does not establish interval coverage or other inferential performance for
that family.
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.03 -0.09 0.09 -0.03
residuals(fit, type = "pearson")
#> [1] 0.4472136 -1.3416408 1.3416408 -0.4472136
residuals(fit, type = "quantile")
#> [1] 0.4472136 -1.3416408 1.3416408 -0.4472136
if (requireNamespace("DHARMa", quietly = TRUE)) {
sims <- simulate(fit, nsim = 20)
dh <- DHARMa::createDHARMa(
simulatedResponse = as.matrix(sims),
observedResponse = fit$model$y,
fittedPredictedResponse = fitted(fit),
integerResponse = FALSE
)
DHARMa::testUniformity(dh)
}
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: simulationOutput$scaledResiduals
#> D = 0.3, p-value = 0.7708
#> alternative hypothesis: two-sided
#>