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

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.103567e+02 6.392210e+00 4.488389e-02 3.151592e-04 2.212939e-06
#> [6] 1.553849e-08
predict(fit, dpar = "sigma", type = "link")
#> [1]   6.813837   1.855080  -3.103676  -8.062433 -13.021189 -17.979946
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.0664096 0.2 1784.466410
#> [2,]   -11.9285020 0.6   13.128502
#> [3,]     0.9120292 1.0    1.087971
#> [4,]     1.3993823 1.4    1.400618
#> [5,]     1.7999957 1.8    1.800004
#> [6,]     2.2000000 2.2    2.200000
#> attr(,"calibrated")
#> [1] FALSE
#> attr(,"prob")
#> [1] 0.025 0.500 0.975
#> attr(,"label")
#> [1] "distributional (plug-in) interval"