Skip to contents

Predict distributional parameters

Source: R/methods.R
predict.drmTMB.Rd

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"),
  prob = c(0.025, 0.5, 0.975),
  ...
)

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", or "quantile".

prob

Numeric vector of probabilities in (0, 1), used only when type = "quantile".

...

Reserved for future prediction options.

Value

A numeric vector for type "response"/"link". 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 fitted.drmTMB() for fitted-row response means or expected scores; use explicit distributional-parameter predictions when a new-data mean must be assembled outside the fitted rows.

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).

See also

fitted.drmTMB(), rho12(), stats::sigma(), fitted_distribution(), exceedance()

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)
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"