After fitting a distributional model, the next question is not only “what are the coefficients?” A useful workflow checks the fit, extracts the fitted distributional parameters, predicts to meaningful covariate values, examines residuals, and simulates from the fitted model.
This article uses a small growth example, but the same workflow
applies to other Gaussian location-scale models in
drmTMB.
Fit the model
Suppose growth is measured along a temperature gradient in two habitat types. The model asks whether temperature and habitat change expected growth, and whether temperature also changes residual variability.
Symbolically:
Matching R syntax:
library(drmTMB)
#>
#> Attaching package: 'drmTMB'
#> The following object is masked from 'package:base':
#>
#> beta
set.seed(12)
n <- 120
fish <- data.frame(
temperature = runif(n, -1.5, 1.5),
habitat = factor(sample(c("reef", "kelp"), n, replace = TRUE))
)
eta_mu <- 1 + 0.7 * fish$temperature + 0.35 * (fish$habitat == "kelp")
eta_sigma <- -0.5 + 0.25 * fish$temperature
fish$growth <- rnorm(n, mean = eta_mu, sd = exp(eta_sigma))
fit <- drmTMB(
bf(growth ~ temperature + habitat, sigma ~ temperature),
family = gaussian(),
data = fish
)Here growth ~ temperature + habitat defines the location
model, and sigma ~ temperature defines the residual
standard deviation model on the log scale.
Check before interpreting
Run check_drm() before interpreting estimates:
check_drm(fit)
#> <drm_check: 10 checks>
#> ok: 10; notes: 0; warnings: 0; errors: 0
#> check status
#> optimizer_convergence ok
#> optimizer_budget ok
#> finite_objective ok
#> fixed_gradient ok
#> sdreport_status ok
#> hessian_positive_definite ok
#> standard_errors_finite ok
#> dropped_rows ok
#> positive_scale ok
#> fixed_effect_design_size ok
#> value
#> 0
#> iterations=21; function=32; gradient=22
#> 110.6
#> max=0.0005009; component=beta_sigma[2]
#> ok
#> TRUE
#> range=[0.06459,0.1109]
#> nobs=120; dropped=0
#> min=0.4388
#> total_mb=0.02132; max_cols=3; largest=mu; largest_class=matrix; largest_density=0.8639
#> message
#> nlminb convergence code is 0.
#> Optimizer evaluation counts recorded; no eval.max or iter.max control was supplied.
#> Objective and log-likelihood are finite.
#> Maximum absolute fixed gradient is <= 0.001; largest component is beta_sigma[2].
#> TMB::sdreport() completed successfully.
#> sdreport reports a positive-definite Hessian.
#> All fixed-effect standard errors are finite.
#> No rows were dropped by model-frame or known-covariance filtering.
#> All fitted scale values are finite and positive.
#> Dense fixed-effect design matrices are modest for this fit.A clean diagnostic table does not prove that the model is
scientifically correct, but it catches common problems early: optimizer
convergence, gradients, optimizer evaluation counts, finite objective
values, Hessian status, dropped rows, finite fixed-effect standard
errors, positive scale values, residual rho12 values near
the boundary, Student-t nu values near the boundary or
close to Gaussian behaviour, known sampling covariance summaries,
random-effect standard deviations near zero, and random-effect design
issues when random effects are present.
Read the ordinary summary first
The ordinary R habit is still the right first habit in
drmTMB:
summary(fit)
#> <summary.drmTMB>
#> estimate std_error
#> mu:(Intercept) 1.3217192 0.08692854
#> mu:temperature 0.6576071 0.06883055
#> mu:habitatreef -0.3776491 0.11093835
#> sigma:(Intercept) -0.4914499 0.06458580
#> sigma:temperature 0.2251913 0.08324771
#> Distributional, random-effect, scale, and correlation parameters:
#> component dpar term estimate minimum
#> fitted:sigma distributional-scale sigma fitted range 0.6188036 0.4388434
#> maximum scale
#> fitted:sigma 0.8465488 response
#> logLik: -110.6
#> convergence: 0summary() puts the fixed-effect estimates,
response-scale scale, shape, random-effect SD, correlation summaries,
fitted covariance summaries, derived variance ratios, log likelihood,
and convergence code in one place when those quantities are present. The
fixed-effect table includes standard errors when
TMB::sdreport() was computed successfully.
Read the printed summary as a map, then use a more specific extractor only when you need more detail:
| Summary component | What it reports | Use it for |
|---|---|---|
coefficients |
Fixed-effect coefficients for each distributional parameter, usually on the model’s link scale. | Direction, sign, and approximate uncertainty of terms in
mu, sigma, rho12, shape, or
zero-inflation formulas. |
parameters |
Direct response-scale parameters, fitted parameter ranges, random-effect SDs, and fitted correlations. | Constant sigma, residual rho12,
sd:mu:(1 | group), and quick response-scale uncertainty
when std_error is available. |
covariance |
Fitted random-effect covariance and correlation rows for supported covariance blocks. | Latent correlation layers such as group-level, phylogenetic, or structured-effect covariance summaries. |
derived |
Simple variance-ratio summaries whose ingredients are unambiguous. | Repeatability or phylogenetic signal point estimates, including the random-effect and residual variances used in the ratio. |
confint |
Confidence-interval rows only when
conf.int = TRUE. |
Checking which requested fixed-effect or profile-likelihood intervals were actually computed. |
The printed summary is meant to answer “what did this fitted model
estimate?” For workflows that need the full object behind a row, use
fixef(), sigma(), rho12(),
ranef(), corpairs(), or
profile_targets() after reading the summary.
For models with random effects, summary() is also the
first place to look for variance components. This compact example adds a
site random intercept to the location model:
set.seed(13)
n_site <- 18
n_per_site <- 6
site <- factor(rep(seq_len(n_site), each = n_per_site))
site_effect <- rnorm(n_site, sd = 0.45)
fish_site <- data.frame(
site = site,
temperature = runif(n_site * n_per_site, -1.5, 1.5),
habitat = factor(
sample(c("reef", "kelp"), n_site * n_per_site, replace = TRUE)
)
)
site_mu <- site_effect[as.integer(fish_site$site)]
fish_site$growth <- rnorm(
nrow(fish_site),
mean = 1 + 0.7 * fish_site$temperature +
0.35 * (fish_site$habitat == "kelp") + site_mu,
sd = 0.35
)
fit_site <- drmTMB(
bf(growth ~ temperature + habitat + (1 | site), sigma ~ 1),
family = gaussian(),
data = fish_site
)
summary(fit_site)
#> <summary.drmTMB>
#> estimate std_error
#> mu:(Intercept) 1.3588788 0.10385484
#> mu:temperature 0.6814229 0.04764308
#> mu:habitatreef -0.2416812 0.07768774
#> sigma:(Intercept) -0.9603029 0.07455733
#> Distributional, random-effect, scale, and correlation parameters:
#> component dpar term estimate std_error
#> sigma distributional-scale sigma (constant) 0.3827769 0.02853883
#> sd:mu:(1 | site) random-effect-sd mu (1 | site) 0.3721619 0.07323934
#> scale
#> sigma response
#> sd:mu:(1 | site) response
#> Derived summaries:
#> quantity level group term estimate
#> derived:repeatability(site) repeatability group site (1 | site) 0.4859419
#> random_effect_variance residual_variance
#> derived:repeatability(site) 0.1385045 0.1465182
#> logLik: -66.61
#> convergence: 0In that output, the sd:mu:(1 | site) row is the fitted
random-intercept standard deviation on the response scale. Squaring it
gives the corresponding random-effect variance component. When the
ingredients are unambiguous, the derived component also
reports repeatability and shows both the random-effect variance and the
residual variance used in that ratio.
When TMB::sdreport() succeeds, direct response-scale
parameter rows such as constant sigma and
sd:mu:(1 | site) also include delta-method standard errors.
Treat those as quick local approximations. For fitted direct SD or
correlation targets, profile-likelihood confidence intervals remain the
safer interval route. drmTMB is a frequentist package, so
these are confidence intervals, not Bayesian credible intervals:
profile_targets(fit_site)
#> parm target_class dpar term
#> 1 fixef:mu:(Intercept) fixed-effect mu (Intercept)
#> 2 fixef:mu:temperature fixed-effect mu temperature
#> 3 fixef:mu:habitatreef fixed-effect mu habitatreef
#> 4 fixef:sigma:(Intercept) fixed-effect sigma (Intercept)
#> 5 sigma distributional-scale sigma (constant)
#> 6 sd:mu:(1 | site) random-effect-sd mu (1 | site)
#> 7 derived:repeatability(site) derived-summary mu (1 | site)
#> tmb_parameter index estimate link_estimate scale transformation
#> 1 beta_mu 1 1.3588788 1.3588788 link linear_predictor
#> 2 beta_mu 2 0.6814229 0.6814229 link linear_predictor
#> 3 beta_mu 3 -0.2416812 -0.2416812 link linear_predictor
#> 4 beta_sigma 1 -0.9603029 -0.9603029 link linear_predictor
#> 5 beta_sigma 1 0.3827769 -0.9603029 response exp
#> 6 log_sd_mu 1 0.3721619 -0.9884264 response exp
#> 7 <NA> NA 0.4859419 NA response variance_ratio
#> target_type profile_ready profile_note
#> 1 direct TRUE ready
#> 2 direct TRUE ready
#> 3 direct TRUE ready
#> 4 direct TRUE ready
#> 5 direct TRUE ready
#> 6 direct TRUE ready
#> 7 derived FALSE derived_target
summary(
fit_site,
conf.int = TRUE,
method = "profile",
ci_parm = "sd:mu:(1 | site)"
)Each distributional parameter has its own coefficient vector:
coef(fit, "mu")
#> (Intercept) temperature habitatreef
#> 1.3217192 0.6576071 -0.3776491
coef(fit, "sigma")
#> (Intercept) temperature
#> -0.4914499 0.2251913The sigma coefficients are on the log-standard-deviation
scale. A positive temperature coefficient in sigma means
that residual variability increases with temperature. Exponentiating the
coefficient gives a multiplicative change in residual standard deviation
for a one-unit temperature increase.
Use the scale of the fitted parameter when deciding what to report:
| Output | Fitted scale | Report when the question is… |
|---|---|---|
coef(fit, "mu") |
response units | expected growth, treatment effects, or mean trait differences |
coef(fit, "sigma") |
log residual SD | multiplicative residual-SD or residual-variance contrasts |
sigma(fit) |
residual SD | fitted predictability or extra heterogeneity on the response scale |
rho12(fit) |
residual correlation | residual coupling between two responses after means and residual SDs are modelled |
corpairs(fit) |
fitted correlation rows | group-level, phylogenetic, or residual correlation layers, depending
on level
|
profile_targets(fit) |
interval target inventory | whether a fitted SD or correlation target can receive a profile interval |
For confidence intervals, keep the fast default summary separate from the expensive profile step. Wald 95% confidence intervals are available for fixed effects:
summary(fit, conf.int = TRUE)
#> <summary.drmTMB>
#> confidence intervals: wald, level = 0.95
#> estimate std_error conf.low conf.high conf.level
#> mu:(Intercept) 1.3217192 0.08692854 1.15134239 1.4920960 0.95
#> mu:temperature 0.6576071 0.06883055 0.52270175 0.7925125 0.95
#> mu:habitatreef -0.3776491 0.11093835 -0.59508428 -0.1602140 0.95
#> sigma:(Intercept) -0.4914499 0.06458580 -0.61803570 -0.3648640 0.95
#> sigma:temperature 0.2251913 0.08324771 0.06202876 0.3883538 0.95
#> conf.method conf.status profile.boundary profile.message
#> mu:(Intercept) wald wald NA <NA>
#> mu:temperature wald wald NA <NA>
#> mu:habitatreef wald wald NA <NA>
#> sigma:(Intercept) wald wald NA <NA>
#> sigma:temperature wald wald NA <NA>
#> Distributional, random-effect, scale, and correlation parameters:
#> component dpar term estimate minimum
#> fitted:sigma distributional-scale sigma fitted range 0.6188036 0.4388434
#> maximum scale conf.status
#> fitted:sigma 0.8465488 response wald_unavailable
#> logLik: -110.6
#> convergence: 0Fixed effects include distributional-parameter coefficients such as
sigma, Student-t nu, zero-inflation, hurdle,
and residual rho12 formula terms. Those Wald intervals are
on the fitted link scale. For example, a Student-t
fixef:nu:x row is an interval for the
log(nu - 2) coefficient, not a response-scale
degrees-of-freedom interval at a particular covariate value.
Profile-likelihood 95% confidence intervals are available for
selected direct targets such as constant sigma,
random-effect SDs, random-effect correlations, and constant residual
rho12. Fixed effects keep Wald confidence intervals in the
same profile summary unless fixed-effect profile targets are selected.
For example, fit_site has sigma ~ 1, so its
residual standard deviation is a direct fitted-object target:
summary(fit_site, conf.int = TRUE, method = "profile", ci_parm = "sigma")For the first model in this article, sigma depends on
temperature. There is no single response-scale
sigma target until you name the covariate row you want to
profile:
sigma_grid <- data.frame(
temperature = c(-1, 1),
habitat = factor("reef", levels = levels(fish$habitat))
)
row.names(sigma_grid) <- c("cool_reef", "warm_reef")
confint(fit, parm = "sigma", method = "profile", newdata = sigma_grid)The same row-specific rule applies to bivariate Gaussian scale and
residual correlation models. In a two-response fit with
sigma1 ~ z1, sigma2 ~ z2, and
rho12 ~ w, profile the biologically meaningful rows
directly:
biv_grid <- data.frame(
x = 0,
z1 = c(-1, 1),
z2 = c(1, -1),
w = c(-0.5, 0.5)
)
row.names(biv_grid) <- c("cool_low_coupling", "warm_high_coupling")
confint(fit_biv, parm = "sigma1", method = "profile", newdata = biv_grid)
confint(fit_biv, parm = "sigma2", method = "profile", newdata = biv_grid)
confint(fit_biv, parm = "rho12", method = "profile", newdata = biv_grid)For random-effect SDs and correlations, copy the exact target name
from profile_targets(). Random-slope terms keep the
coefficient suffix in the target name, and residual rho12
remains separate from group-level or phylogenetic correlation rows:
profile_targets(fit_slope)
summary(
fit_slope,
conf.int = TRUE,
method = "profile",
ci_parm = "sd:mu:(1 + temperature | p | site):temperature"
)
summary(
fit_slope,
conf.int = TRUE,
method = "profile",
ci_parm = "cor:mu:cor((Intercept),temperature | p | site)"
)
summary(
fit_biv_group,
conf.int = TRUE,
method = "profile",
ci_parm = "cor:mu:cor(mu1:(Intercept),mu2:(Intercept) | p | id)"
)Use the profile call when a report needs fixed-effect confidence
intervals and profile-likelihood intervals for variance components, SDs,
or correlations in one summary. Successful profile rows also include
profile.boundary and profile.message; values
other than "ok" tell you to check boundary or
near-correlation-limit behaviour before interpreting the interval.
For correlations, use profile intervals for fitted-model inference
when the target is profile-ready. The Phase 18 simulation helpers also
include a Fisher-z back-transformed Wald route for coverage tables with
supplied correlation standard errors, but that helper is a simulation
producer contract, not a public confint() default for
fitted correlation rows.
Use profile_targets(fit) before asking for a profile
interval. A row is profile-ready only when it maps to a direct TMB
target and the fitted object kept its TMB object. Interval-aware output
uses conf.status to explain rows without intervals:
newdata_required means profile a supplied row with
confint(..., newdata = ...), while
derived_interval_unavailable means the point estimate
combines several quantities and needs a later derived-profile
method.
Read conf.status as an action column:
conf.status |
What it means | What to do next |
|---|---|---|
wald |
A Wald interval was returned for a fixed-effect row. | Use it for routine fixed-effect summaries, while remembering it is asymptotic. |
profile |
A profile-likelihood interval was returned. | Inspect profile.boundary and
profile.message before treating the interval as
stable. |
profile_ready |
The row is a direct profile target, but the current summary call did not request it. | Use the row’s parm value in
confint(..., method = "profile") or in
summary(..., method = "profile", ci_parm = ...). |
newdata_required |
The fitted surface needs a concrete row before it can be profiled. | Build a biologically meaningful newdata row and call
confint(..., newdata = ...). |
derived_interval_unavailable |
The point estimate combines multiple fitted quantities and has no validated interval method yet. | Report the point estimate as exploratory, or simplify to a direct target if interval support is essential. |
wald_unavailable |
Wald uncertainty is not available for this row, or
sdreport() was skipped or failed. |
Use a profile-ready target when one exists; otherwise report the missing interval explicitly. |
target_unavailable or
profile_unavailable
|
The row is descriptive or has no current direct interval target. | Treat the interval as unsupported and check
profile_targets(fit) before trying another route. |
not_requested |
The table is reporting a point estimate without an interval request. | Refit is not needed; request intervals only for the rows you plan to interpret. |
Parametric-bootstrap intervals are still a planned interval route,
not a public method value. Calls such as
confint(fit, method = "bootstrap"),
summary(fit, conf.int = TRUE, method = "bootstrap"), or
corpairs(fit, conf.int = TRUE, method = "bootstrap") stop
before interval work begins. For now, use Wald intervals for routine
fixed-effect summaries, profile intervals for direct profile-ready
targets, and explicit simulation studies when bootstrap coverage is the
scientific target.
Predict distributional parameters
predict() returns one distributional parameter at a
time. For fitted rows, fitted(fit) and
predict(fit, dpar = "mu") return the same location values
for this fixed-effect model:
head(fitted(fit))
#> [1] 0.4721452 1.5709838 2.1949328 0.8667508 0.2917530 0.4021785
head(predict(fit, dpar = "mu"))
#> [1] 0.4721452 1.5709838 2.1949328 0.8667508 0.2917530 0.4021785
head(sigma(fit))
#> [1] 0.4573171 0.7582294 0.8249552 0.5234837 0.4892755 0.4464902For new covariate values, supply newdata. These are
fixed-effect, population-level predictions. You can write the grid
directly:
new_fish <- data.frame(
temperature = c(-1, 0, 1),
habitat = factor("reef", levels = levels(fish$habitat))
)
predict(fit, newdata = new_fish, dpar = "mu")
#> [1] 0.2864629 0.9440701 1.6016772
predict(fit, newdata = new_fish, dpar = "sigma")
#> [1] 0.4883899 0.6117388 0.7662410
predict(fit, newdata = new_fish, dpar = "sigma", type = "link")
#> [1] -0.7166411 -0.4914499 -0.2662586Or build an explicit prediction grid. This keeps the focal term, conditioned nuisance predictors, and grid rule attached to the data:
temperature_grid <- prediction_grid(
fit,
focal = "temperature",
at = list(temperature = c(-1, 0, 1)),
condition = list(habitat = "reef")
)
attr(temperature_grid, "prediction_grid")
#> $focal_terms
#> [1] "temperature"
#>
#> $conditioned_terms
#> [1] "habitat"
#>
#> $margin
#> [1] "mean_reference"
#>
#> $weights
#> [1] "equal"
#>
#> $grid_source
#> [1] "mean_reference_prediction_grid"
#>
#> $reference_terms
#> character(0)
#>
#> $predictor_terms
#> [1] "temperature" "habitat"
#>
#> $n_source_rows
#> [1] 120
#>
#> $n_grid_rows
#> [1] 3Use predict_parameters() when the interpretation task
needs several distributional parameters on the same grid:
temperature_surface <- predict_parameters(
fit,
newdata = temperature_grid,
dpar = c("mu", "sigma"),
conf.int = TRUE
)
temperature_surface
#> row row_label dpar component type estimate std.error
#> 1 1 1 mu location response 0.2864629 0.08099454
#> 2 2 2 mu location response 0.9440701 0.07419006
#> 3 3 3 mu location response 1.6016772 0.11799790
#> 4 1 1 sigma distributional-scale response 0.4883899 0.05061759
#> 5 2 2 sigma distributional-scale response 0.6117388 0.03950964
#> 6 3 3 sigma distributional-scale response 0.7662410 0.08203217
#> conf.low conf.high conf.level conf.status interval_source temperature
#> 1 0.1277166 0.4452093 0.95 wald wald -1
#> 2 0.7986602 1.0894799 0.95 wald wald 0
#> 3 1.3704056 1.8329489 0.95 wald wald 1
#> 4 0.3986086 0.5983933 0.95 wald wald -1
#> 5 0.5390022 0.6942911 0.95 wald wald 0
#> 6 0.6212086 0.9451337 0.95 wald wald 1
#> habitat
#> 1 reef
#> 2 reef
#> 3 reef
#> 4 reef
#> 5 reef
#> 6 reef
unique(temperature_surface[, c(
"dpar",
"conf.status",
"interval_source",
"conf.level"
)])
#> dpar conf.status interval_source conf.level
#> 1 mu wald wald 0.95
#> 4 sigma wald wald 0.95Here both rows receive Wald 95% confidence intervals because the
table uses an explicit newdata grid and both
mu and sigma have ordinary fixed-effect bases
in this fitted model. The interval is a confidence band for the fitted
distributional-parameter surface, not a prediction interval for future
raw growth observations.
If you request intervals for fitted rows without
newdata, the table keeps the point estimates and says what
is missing:
fitted_interval_status <- predict_parameters(
fit,
dpar = "mu",
conf.int = TRUE
)
unique(fitted_interval_status[, c(
"dpar",
"conf.status",
"interval_source",
"conf.level"
)])
#> dpar conf.status interval_source conf.level
#> 1 mu newdata_required not_available 0.95Read newdata_required as an instruction to choose a row
on purpose. For a surface, build a compact grid such as
temperature_grid instead of asking the plot to decide which
fitted rows deserve intervals.
Use an empirical grid when the question is about the fitted response at focal temperatures while averaging over the fitted-row covariate distribution. Here each requested temperature is crossed with the fitted model rows, so the habitat distribution is the one in the fitted data:
temperature_empirical <- prediction_grid(
fit,
focal = "temperature",
at = list(temperature = c(-1, 0, 1)),
margin = "empirical"
)
attr(temperature_empirical, "prediction_grid")
#> $focal_terms
#> [1] "temperature"
#>
#> $conditioned_terms
#> NULL
#>
#> $margin
#> [1] "empirical"
#>
#> $weights
#> [1] "equal"
#>
#> $grid_source
#> [1] "empirical_prediction_grid"
#>
#> $reference_terms
#> [1] "habitat"
#>
#> $predictor_terms
#> [1] "temperature" "habitat"
#>
#> $n_source_rows
#> [1] 120
#>
#> $n_grid_rows
#> [1] 360Then reduce the empirical grid by the focal term:
marginal_parameters(
fit,
newdata = temperature_empirical,
dpar = c("mu", "sigma"),
by = "temperature"
)
#> dpar component type temperature estimate n conf.status
#> 1 mu location response -1 0.4406697 120 not_requested
#> 2 mu location response 0 1.0982768 120 not_requested
#> 3 mu location response 1 1.7558840 120 not_requested
#> 4 sigma distributional-scale response -1 0.4883899 120 not_requested
#> 5 sigma distributional-scale response 0 0.6117388 120 not_requested
#> 6 sigma distributional-scale response 1 0.7662410 120 not_requested
#> interval_source
#> 1 not_available
#> 2 not_available
#> 3 not_available
#> 4 not_available
#> 5 not_available
#> 6 not_availableThese helpers are tables, not plotting functions. They are the
current Phase 17 surface for future visualization work: each row names
the distributional parameter, fitted component, prediction scale,
estimate, interval status, and grid values. By default, prediction
tables are point-estimate summaries with
conf.status = "not_requested" and
interval_source = "not_available". When you request
conf.int = TRUE on an explicit newdata grid,
ordinary fixed-effect distributional parameters receive Wald
fixed-effect intervals with conf.status = "wald" and
interval_source = "wald". Rows that cannot receive that
interval, such as direct random-effect SD surfaces, keep an explicit
unavailable status instead of silently drawing a band.
Predict a random-effect SD surface
Random-effect scale formulas answer a different question from
residual sigma. In the next example,
sd(site) ~ reef_cover asks whether the among-site spread in
expected growth changes along a site-level reef-cover gradient. The
predictor is constant within each site:
set.seed(14)
n_site <- 14
n_per_site <- 5
site <- factor(rep(seq_len(n_site), each = n_per_site))
reef_cover_site <- seq(-1, 1, length.out = n_site)
reef_cover <- reef_cover_site[as.integer(site)]
temperature <- runif(n_site * n_per_site, -1.5, 1.5)
site_sd <- exp(-0.65 + 0.45 * reef_cover_site)
site_effect <- rnorm(n_site, sd = site_sd)
fish_site_scale <- data.frame(
site = site,
temperature = temperature,
reef_cover = reef_cover
)
fish_site_scale$growth <- rnorm(
nrow(fish_site_scale),
mean = 1 + 0.6 * fish_site_scale$temperature +
site_effect[as.integer(fish_site_scale$site)],
sd = 0.35
)
fit_site_scale <- drmTMB(
bf(
growth ~ temperature + (1 | site),
sigma ~ 1,
sd(site) ~ reef_cover
),
family = gaussian(),
data = fish_site_scale
)The fitted equation for the site random-effect SD is
log(sd_site) = alpha_0 + alpha_1 reef_cover. Use a grid
when you want the site-level SD surface on the response scale:
sd_site_grid <- prediction_grid(
fit_site_scale,
focal = "reef_cover",
at = list(reef_cover = c(-1, 0, 1)),
condition = list(temperature = 0)
)
sd_site_surface <- predict_parameters(
fit_site_scale,
newdata = sd_site_grid,
dpar = "sd(site)"
)
sd_site_surface
#> row row_label dpar component type estimate
#> 1 1 1 sd(site) random-effect-sd-model response 0.1626495
#> 2 2 2 sd(site) random-effect-sd-model response 0.3942447
#> 3 3 3 sd(site) random-effect-sd-model response 0.9556064
#> conf.status interval_source temperature reef_cover
#> 1 not_requested not_available 0 -1
#> 2 not_requested not_available 0 0
#> 3 not_requested not_available 0 1The table reports the random-effect-sd-model component.
Its estimates are random-intercept SDs, not residual SDs and not raw
responses. If you request a Wald interval for this direct random-effect
SD surface, the table keeps the point estimate and marks the interval
status explicitly:
sd_site_interval_status <- predict_parameters(
fit_site_scale,
newdata = sd_site_grid,
dpar = "sd(site)",
conf.int = TRUE
)
unique(sd_site_interval_status[, c(
"dpar",
"component",
"conf.status",
"interval_source",
"conf.level"
)])
#> dpar component conf.status interval_source conf.level
#> 1 sd(site) random-effect-sd-model wald_unavailable not_available 0.95Read that status as a boundary, not as a plotting failure. Direct random-effect SD surfaces need a profile-likelihood or bootstrap route before a confidence band is available, so the current plotting helper should show the line without a ribbon. If the report needs one row per reef-cover value, reduce the same explicit grid:
marginal_parameters(
fit_site_scale,
newdata = sd_site_grid,
dpar = "sd(site)",
by = "reef_cover"
)
#> dpar component type reef_cover estimate n conf.status
#> 1 sd(site) random-effect-sd-model response -1 0.1626495 1 not_requested
#> 2 sd(site) random-effect-sd-model response 0 0.3942447 1 not_requested
#> 3 sd(site) random-effect-sd-model response 1 0.9556064 1 not_requested
#> interval_source
#> 1 not_available
#> 2 not_available
#> 3 not_availableEstimate marginal means for mu
When the question is an adjusted mean of the location parameter
mu, the first emmeans bridge can build an
estimated marginal mean table directly from a fixed-effect univariate
fit. Here the reference grid holds temperature at zero and
estimates mu for each habitat:
if (requireNamespace("emmeans", quietly = TRUE)) {
emmeans::emmeans(
fit,
~habitat,
at = list(temperature = 0)
)
}
#> habitat emmean SE df asymp.LCL asymp.UCL
#> kelp 1.322 0.0869 Inf 1.151 1.49
#> reef 0.944 0.0742 Inf 0.799 1.09
#>
#> Confidence level used: 0.95Read this output as estimated marginal means of the native
distributional parameter mu, not as summaries of
sigma or as a random-effect, bivariate, zero-inflated,
hurdle, ordinal, or slope workflow. Generic emmeans
contrasts can be formed from this returned mu grid, but
broader drmTMB-specific contrast helpers remain a separate future
contract. For other distributional targets, use explicit
prediction_grid() and predict_parameters()
tables until the corresponding reference-grid algebra and tests
exist.
The first optional plotting helper consumes that table directly. For a figure, use a denser grid than the small printed table and request intervals before plotting:
temperature_plot_grid <- prediction_grid(
fit,
focal = c("temperature", "habitat"),
at = list(temperature = seq(-1.5, 1.5, length.out = 40))
)
pred_temperature <- predict_parameters(
fit,
newdata = temperature_plot_grid,
dpar = c("mu", "sigma"),
conf.int = TRUE
)
unique(pred_temperature[, c(
"dpar",
"conf.status",
"interval_source",
"conf.level"
)])
#> dpar conf.status interval_source conf.level
#> 1 mu wald wald 0.95
#> 81 sigma wald wald 0.95
if (requireNamespace("ggplot2", quietly = TRUE)) {
plot_parameter_surface(pred_temperature, x = "temperature", colour = "habitat")
}
For a discrete focal predictor, the same interval columns draw
interval bars rather than a ribbon. This example asks for 90% Wald
intervals, so the table records conf.level = 0.9 on rows
where conf.status = "wald":
habitat_plot_grid <- prediction_grid(
fit,
focal = "habitat",
condition = list(temperature = 0)
)
pred_habitat <- predict_parameters(
fit,
newdata = habitat_plot_grid,
dpar = "mu",
conf.int = TRUE,
conf.level = 0.90
)
unique(pred_habitat[, c(
"dpar",
"conf.status",
"interval_source",
"conf.level"
)])
#> dpar conf.status interval_source conf.level
#> 1 mu wald wald 0.9
if (requireNamespace("ggplot2", quietly = TRUE)) {
plot_parameter_surface(pred_habitat, x = "habitat", line = FALSE)
}
Treat conf.level as the requested level for rows where
an interval was computed or attempted. It is not a success flag by
itself; read it together with conf.status and
interval_source.
The figure is a convenience layer, not a new estimand. Keep reference grids, marginalization choices, and interval sources visible instead of hiding those decisions inside a plot.
For reports, pair model surfaces with the raw response pattern. The raw-data panel uses the observed response scale:
if (requireNamespace("ggplot2", quietly = TRUE)) {
ggplot2::ggplot(
fish,
ggplot2::aes(x = temperature, y = growth, colour = habitat)
) +
ggplot2::geom_point(alpha = 0.65) +
ggplot2::labs(
x = "Temperature",
y = "Observed growth",
colour = "Habitat"
)
}
The fitted panels then use the prediction table. Keep mu
and sigma on their own y-axes because they are different
biological quantities. When a plot keeps one distributional parameter,
plot_parameter_surface() labels the y-axis with the
parameter and prediction scale. The bands in these panels come from the
Wald fixed-effect intervals in pred_temperature; if those
columns were absent or marked
interval_source = "not_available", the helper would draw
only lines and points:
if (requireNamespace("ggplot2", quietly = TRUE)) {
plot_parameter_surface(
pred_temperature,
x = "temperature",
colour = "habitat",
dpar = "mu",
facet = NULL
)
}
if (requireNamespace("ggplot2", quietly = TRUE)) {
plot_parameter_surface(
pred_temperature,
x = "temperature",
colour = "habitat",
dpar = "sigma",
facet = NULL
)
}
Do not put raw growth points on the sigma
panel. Raw observations show the response; fitted sigma
shows the residual standard deviation implied by the model. The same
rule applies to uncertainty: draw a band only when the table names a
real interval source such as "wald" or
"profile".
For this Gaussian model, the response-scale sigma
predictions are residual standard deviations. The link-scale predictions
are the corresponding log-standard deviations. Report residual variance
as sigma(fit)^2 or
predict(fit, newdata = new_fish, dpar = "sigma")^2; do not
square the sigma coefficient, because coefficients are
fitted on the log-SD scale.
In meta-analytic Gaussian models with
meta_known_V(V = V), sigma(fit) still returns
the modelled residual heterogeneity scale. Simulation and Pearson
residuals combine the known sampling covariance with that residual scale
internally.
Inspect residuals
Response residuals are observed minus fitted location:
Pearson residuals divide by the fitted observation scale:
head(residuals(fit, type = "pearson"))
#> [1] -0.2250256 0.8277832 0.5261101 -1.3961680 -1.2289151 1.2745363For bivariate Gaussian fits,
residuals(fit, type = "pearson") returns a two-column
matrix whitened by the fitted sigma1, sigma2,
and residual correlation rho12.
Simulate from the fitted model
Simulation is a compact way to ask whether the fitted model generates data on the same scale as the observed response:
sims <- simulate(fit, nsim = 3, seed = 1)
head(sims)
#> sim_1 sim_2 sim_3
#> 1 0.18565716 0.2407622 0.7956104
#> 2 1.71022758 2.5893154 2.3550747
#> 3 1.50557669 2.0179144 2.3792942
#> 4 1.70185429 0.7727559 0.4067617
#> 5 0.45297304 0.2427321 0.8607630
#> 6 0.03584737 0.7203770 -0.4908756For this univariate model, the output has one column per simulated
data set. For bivariate Gaussian models, simulations are returned in
paired columns such as sim_1_y1 and
sim_1_y2.
A compact post-fit checklist
For routine fitted models, the post-fit loop is:
-
check_drm(fit)before interpreting estimates. -
summary(fit)for fixed effects, response-scale parameters, variance components, derived summaries, interval status, and the extractor to use next. -
profile_targets(fit)before asking for profile-likelihood intervals on fitted SD or correlation targets. -
predict(fit, dpar = ...)for fitted or new-data distributional parameters. -
residuals(fit, type = "pearson")to inspect standardized lack of fit. -
simulate(fit)to compare fitted-model behaviour with the observed data.
The same loop will extend as more families arrive. The parameter names change, but the discipline stays the same: check the fit, then interpret each distributional parameter on its own scale.