Predict distributional parameters

predict() returns fitted or predicted values for one distributional parameter of a drmTMB fit.

Usage

# S3 method for class 'drmTMB'
predict(
  object,
  newdata = NULL,
  dpar = NULL,
  type = c("response", "link", "quantile", "mean"),
  prob = c(0.025, 0.5, 0.975),
  re.form = NULL,
  ...
)

Arguments

object: A drmTMB fit.
newdata: Optional data frame for prediction. If omitted, fitted rows are used. When supplied, newdata must include the predictors used by the requested dpar; required predictor values must be complete, required numeric predictors must be finite, and factor predictors must use fitted levels. Transformed predictor terms, such as log(size), must also evaluate to finite design-matrix values.
dpar: Distributional parameter to predict. If NULL, the first fitted distributional parameter is used. For type = "quantile" on a bivariate biv_gaussian fit, dpar instead selects which response's marginal quantile to compute (see Details).
type: Prediction scale: "response", "link", "quantile", or "mean". "mean" returns the response mean E[Y] and takes no dpar (see Details).
prob: Numeric vector of probabilities in (0, 1), used only when type = "quantile".
re.form: Must be NULL (the default), which is the only behaviour: training rows are predicted conditionally on the fitted random effects and newdata rows at the population level. Any other value is an error.
...: Reserved for future prediction options; any value is an error (for example se.fit, which predict() does not compute; use predict_parameters() with conf.int = TRUE for standard errors).

Value

A numeric vector for type "response"/"link"/"mean". For type = "quantile", a numeric matrix with one row per observation and one column per prob (columns named as percentages, e.g. "2.5%"), with attr(., "calibrated") == FALSE.

Details

By default, predictions are returned on the distributional parameter's response scale. For positive scale parameters such as sigma, this means the exponentiated value. For bivariate residual correlation rho12 or a fitted corpair() model, this means the correlation scale. Use type = "link" to return the linear predictor instead.

When newdata = NULL, predictions are for the fitted rows and include currently implemented conditional random-effect contributions for mu, including registry-backed q > 2 ordinary covariance blocks, bivariate mu1/mu2, phylogenetic mu, and residual-scale sigma including bivariate sigma1/sigma2 blocks. Fitted-row predictions also include ordinary zero-one-beta zoi and coi random-effect contributions. When newdata is supplied, predictions are fixed-effect, population-level predictions for the supplied rows.

type = "quantile" returns per-row conditional quantiles of the fitted RESPONSE distribution (not of a linear predictor): qnorm-style inverse CDF evaluation via fitted_distribution()$q() at predict_parameters()'s fixed-effect, population-level parameter estimates – see fitted_distribution() for the fixed-effect-only scope and the "spike"/"unimplemented" status rules this inherits (a "spike"-status family emits a one-time warning; an unregistered model_type raises a clear error). For bivariate biv_gaussian fits, dpar selects which response's MARGINAL quantile to return ("mu1"/"sigma1" for response 1, "mu2"/"sigma2" for response 2); the joint rho12 correlation, any joint tail structure, and any known bivariate sampling covariance are ignored – this is a marginal-only computation, not a joint bivariate quantile. The result carries attr(., "calibrated") <- FALSE: this is a distributional (plug-in) interval at the point estimate theta_hat, not a calibrated-coverage interval, and it does not propagate theta_hat uncertainty.

For retained training data with a modelled missing predictor, prediction uses that predictor's finalized conditional plug-in in the fixed-effect design. This is a distributional-parameter prediction, not an integrated response mean, and fixed-effect covariance does not propagate imputation-parameter uncertainty through the design basis.

More generally, type = "response" means "on the response scale of the requested distributional parameter", not always E[Y]. For zero-inflated and hurdle count families it returns the component count mean when dpar = "mu", not the unconditional mean after structural-zero or hurdle mixing. For lognormal fits it returns mu on the log-response scale, not exp(mu + sigma^2 / 2). For cumulative-logit fits it returns the location linear predictor on its response scale, not an expected ordinal score. Use type = "mean" for E[Y].

type = "mean" returns the response mean E[Y] for each row, so it differs from type = "response" exactly where the two answer different questions: exp(mu + sigma^2 / 2) for lognormal fits; (1 - zi) * mu for zero-inflated count fits; the hurdle mean for hurdle fits; and the expected score sum(k * Pr(y = k)) on the category indices 1:K for cumulative-logit fits. Truncated and zero-one-beta fits use their own mean formulas, and families whose mean is mu return mu. With newdata = NULL it equals fitted.drmTMB(); with newdata it uses the population-level parameter predictions for the supplied rows. dpar must be left NULL, and bivariate fits are an error.

Binary predictor labels follow the fitted encoding, regardless of the levels present in the new batch. Numeric new values match raw numeric fitted labels; for nonnumeric fitted labels, numeric 0/1 values specify the encoded states. If distinct labels become identical as numbers (for example, "01" and "1"), supply character or factor values to preserve their identity.

For families in drm_clamped_scale_families(), the log(sigma) linear predictor is soft-clamped inside the TMB likelihood to keep the objective finite. predict(dpar = "sigma") (and the bivariate sigma1/sigma2) reports that same clamped scale, on both type = "link" and type = "response", so the value matches what the likelihood actually evaluated. When check_drm() reports the clamp active, this makes the fit's own diagnostics honest, not correct: the estimate itself is still unreliable near the clamp (see the rescaling advice there).

For predictor-dependent random-effect SD parameters (for example dpar = "sd(id)" or "sd_animal(id)"), the value is the SD multiplier exp(W beta) on the latent effect, not a marginal SD. When the effect is structured (animal(), relmat(), phylo(), spatial()), the multiplier applies to the supplied kernel as given, so the marginal SD of group g is multiplier * sqrt(C_gg); predict() never returns that product. For a unit-diagonal kernel the two agree. type = "link" is the log SD after the same soft clamp the likelihood applies (the identity inside [-12, 12], saturating towards [-15, 15]), not the raw W beta; type = "response" is its exponential. With newdata = NULL the result has one value per fitted group, named by group level. With newdata, each row is one group-level row (not an observation row): the result has one value per row, named by rownames(newdata), and a structured SD for such a row is a population multiplier with no kernel node and no sqrt(C_gg) factor.

Examples

dat <- data.frame(
  y = c(0.2, 0.5, 1.1, 1.4, 1.8, 2.2),
  x = c(-1, -0.5, 0, 0.5, 1, 1.5)
)
fit <- drmTMB(bf(y ~ x, sigma ~ x), data = dat)
#> Warning: `drmTMB()`: the fitted `log(sigma)` spans 21.7 across observations, more than
#> `log(1e4)` = 9.2.
#> ℹ A scale range wider than 10^4 usually means the scale is running off (or the
#>   design is extreme), so the reported convergence, estimates and standard
#>   errors may be unreliable even when the clamp is not hit. Rescale the
#>   response, simplify the `sigma` formula, or use a penalized/MAP fit;
#>   `check_drm()` reports the fit state.
predict(fit, dpar = "mu")
#> [1] 0.2 0.6 1.0 1.4 1.8 2.2
predict(fit, dpar = "sigma")
#> [1] 9.103562e+02 6.392210e+00 4.488391e-02 3.151595e-04 2.297952e-06
#> [6] 3.412479e-07
predict(fit, dpar = "sigma", type = "link")
#> [1]   6.813836   1.855080  -3.103676  -8.062432 -12.983492 -14.890657
predict(fit, newdata = data.frame(x = c(0, 1)), dpar = "mu")
#> [1] 1.0 1.8
predict(fit, type = "quantile", prob = c(0.025, 0.5, 0.975))
#>               2.5% 50%       97.5%
#> [1,] -1784.0653817 0.2 1784.465382
#> [2,]   -11.9285012 0.6   13.128501
#> [3,]     0.9120292 1.0    1.087971
#> [4,]     1.3993823 1.4    1.400618
#> [5,]     1.7999955 1.8    1.800005
#> [6,]     2.1999993 2.2    2.200001
#> attr(,"calibrated")
#> [1] FALSE
#> attr(,"prob")
#> [1] 0.025 0.500 0.975
#> attr(,"label")
#> [1] "distributional (plug-in) interval"
predict(fit, type = "mean")
#> [1] 0.2 0.6 1.0 1.4 1.8 2.2

Usage

Arguments

Value

Details

See also

Examples