Fitting with phylo(1 | species)
The phylo() term goes inside the mean (mu)
formula. It takes a random-effect specification,
1 | species, and the tree as tree = tree. Wrap
the formulas in bf() (the drmTMB formula
builder) and pass a family, exactly as for any other drmTMB
model.
fit <- drmTMB(
bf(trait ~ x + phylo(1 | species, tree = tree),
sigma ~ 1),
family = gaussian(),
data = dat
)
check_drm(fit)
#> <drm_check: 16 checks>
#> ok: 15; notes: 1; warnings: 0; errors: 0
#> check status
#> optimizer_convergence ok
#> optimizer_budget ok
#> finite_objective ok
#> logsigma_clamp_active ok
#> fixed_gradient ok
#> sdreport_status ok
#> hessian_positive_definite ok
#> standard_errors_finite ok
#> standard_errors_inflated ok
#> dropped_rows ok
#> positive_scale ok
#> random_effect_sd_boundary ok
#> interval_reliability_scope note
#> fixed_effect_design_size ok
#> phylo_mu_replication ok
#> phylo_mu_diagnostics ok
#> value
#> 0
#> iterations=19; function=29; gradient=20
#> 26.37
#> <NA>
#> max=0.0000009013; component=beta_mu[2]
#> ok
#> TRUE
#> range=[0.03125,0.4704]
#> n_inflated=0; max_se=0.4704; median_se=0.07706
#> nobs=96; dropped=0
#> min=0.2693
#> min=0.7795; boundary=0.0001000; term=mu.phylo(1 | species)
#> sd_targets=1; assessed_here=0
#> total_mb=0.01537; max_cols=2; largest=mu; largest_class=matrix; largest_density=1.000
#> min_species_n=6
#> group=species; n_species=16; min_species_n=6; phylo_sd=0.7795; sd_ratio=2.894
#> 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.
#> The log(sigma) clamp is not active at the optimum.
#> Maximum absolute fixed gradient is <= 0.001; largest component is beta_mu[2].
#> TMB::sdreport() completed successfully.
#> sdreport reports a positive-definite Hessian.
#> All fixed-effect standard errors are finite.
#> No fixed-effect standard error is inflated relative to the others.
#> No rows were dropped by model-frame or known-covariance filtering.
#> All fitted scale values are finite and positive.
#> All fitted random-effect standard deviations are finite, positive, and above the requested lower-boundary warning threshold.
#> This fit has 1 random-effect standard-deviation target. `check_drm()` assesses the fit, not interval reliability: a target can pass every check above and still return an interval that `confint()` warns about at a variance boundary. Before reporting an interval, call `confint()` and read `conf.status` and any boundary warning.
#> Dense fixed-effect design matrices are modest for this fit.
#> Every observed species has at least two fitted observations.
#> The phylogenetic random effect has replicated species and a non-negligible fitted SD relative to residual scale.Start with check_drm(): it reports convergence, Hessian,
scale, and phylogenetic-replication diagnostics in one public table.
Then use the public coefficient and target extractors to see the fixed
slope and the two SDs:
coef(fit, "mu")
#> (Intercept) x
#> 1.2837197 -0.4055049
sd_targets <- profile_targets(fit)
sd_targets[
sd_targets$parm %in% c("sigma", "sd:mu:phylo(1 | species)"),
c("parm", "estimate", "scale", "profile_ready", "profile_note")
]
#> parm estimate scale profile_ready profile_note
#> 4 sigma 0.2693226 response TRUE ready
#> 5 sd:mu:phylo(1 | species) 0.7795200 response TRUE readyThe recovered slope (x) is close to its true value of
-0.4, and the phylogenetic and residual SDs are in the
neighbourhood of 0.8 and 0.3. The fixed
intercept need not match the simulated 0.5: with only 16
species the species-level deviations carry a non-zero sample mean that
the random effect absorbs, so the intercept and the deviations trade
off. This is expected behaviour, not a fitting error – the
scientifically interpretable quantities are the slope and the two
SDs.
Reading the phylogenetic signal
summary() adds a derived row for the phylogenetic signal
\(\lambda\), alongside the
random-effect and residual variances it is built from.
summary(fit)$derived[, c(
"quantity", "estimate",
"random_effect_variance", "residual_variance"
)]
#> quantity estimate
#> derived:phylogenetic_signal(species) phylogenetic_signal 0.8933605
#> random_effect_variance residual_variance
#> derived:phylogenetic_signal(species) 0.6076515 0.07253469estimate is \(\lambda\). A value well above zero says
that, after accounting for the covariate, related species really do
resemble each other more than unrelated ones. The fitted object also
contains the augmented phylogenetic state used by the sparse-precision
representation. To inspect only the species-tip deviations, match that
state by the tree’s tip labels:
Uncertainty for the two SDs
Variance components have their own confint() target. The
interval is on the response (SD) scale, obtained by transforming the
Wald interval for the log-SD parameter. The table contains both the
residual SD and the phylogenetic location SD.
confint(fit, parm = "variance_components")[, c(
"parm", "lower", "upper", "scale"
)]
#> parm lower upper scale
#> 1 sigma 0.2315703 0.3132297 response
#> 2 sd:mu:phylo(1 | species) 0.5046785 1.3699263 responseA transformed log-SD interval is necessarily positive, so its lower
bound should not be read as a test of a zero variance
component. Use it to describe uncertainty in the magnitude of each SD.
The earlier check_drm(fit) table reports diagnostics for
the recognised phylogenetic layer.
The figure uses Confidence Eyes rather than flat interval bars. Each pale eye is the finite 95% Wald confidence region shaped on the log-SD scale: compatibility is greatest near its centre and tapers towards the endpoints. The eye is a frequentist compatibility display, not a posterior density.
vc <- confint(fit, parm = "variance_components")
sd_targets <- profile_targets(fit)
target <- c("sd:mu:phylo(1 | species)", "sigma")
interval_row <- match(target, vc$parm)
target_row <- match(target, sd_targets$parm)
stopifnot(!anyNA(interval_row), !anyNA(target_row))
sd_tab <- data.frame(
label = factor(
c("Phylogenetic SD", "Residual SD"),
levels = c("Residual SD", "Phylogenetic SD")
),
estimate = c(
sd_targets$estimate[target_row]
),
lower = vc$lower[interval_row],
upper = vc$upper[interval_row]
)
stopifnot(
all(is.finite(unlist(sd_tab[c("estimate", "lower", "upper")]))),
all(sd_tab$lower > 0),
all(sd_tab$lower <= sd_tab$estimate),
all(sd_tab$estimate <= sd_tab$upper)
)
sd_eye <- do.call(rbind, lapply(seq_len(nrow(sd_tab)), function(i) {
log_lower <- log(sd_tab$lower[i])
log_upper <- log(sd_tab$upper[i])
log_centre <- 0.5 * (log_lower + log_upper)
log_value <- seq(log_lower, log_upper, length.out = 401L)
half_width <- 0.5 * (log_upper - log_lower)
height <- pmax(1 - ((log_value - log_centre) / half_width)^2, 0)
data.frame(
label = as.character(sd_tab$label[i]),
value = exp(log_value),
height = height
)
}))
sd_eye$label <- factor(sd_eye$label, levels = levels(sd_tab$label))
sd_eye$y <- as.numeric(sd_eye$label)
sd_tab$y <- as.numeric(sd_tab$label)
ggplot2::ggplot() +
ggplot2::geom_ribbon(
data = sd_eye,
ggplot2::aes(
x = value,
ymin = y - 0.20 * height,
ymax = y + 0.20 * height,
group = label
),
fill = "#0072B2",
alpha = 0.24,
colour = NA
) +
ggplot2::geom_point(
data = sd_tab,
ggplot2::aes(x = estimate, y = y),
shape = 21,
fill = "white",
colour = "#0072B2",
size = 3.0,
stroke = 1.0
) +
ggplot2::scale_y_continuous(
breaks = seq_along(levels(sd_tab$label)),
labels = levels(sd_tab$label),
expand = ggplot2::expansion(add = 0.38)
) +
ggplot2::scale_x_continuous(
limits = c(0, NA),
expand = ggplot2::expansion(mult = c(0, 0.04))
) +
ggplot2::labs(
x = "Standard deviation (response scale)",
y = NULL
) +
ggplot2::theme_minimal(base_size = 12.5) +
ggplot2::theme(
axis.line.x = ggplot2::element_line(colour = "grey40", linewidth = 0.35),
axis.ticks.x = ggplot2::element_line(colour = "grey40", linewidth = 0.35),
panel.grid.major.y = ggplot2::element_blank(),
panel.grid.minor = ggplot2::element_blank(),
axis.text.y = ggplot2::element_text(colour = "grey15")
)
Confidence Eyes for the two response-scale SDs. Pale shapes are the default finite 95% Wald confidence regions, constructed on the log-SD scale; hollow circles are the raw fitted SDs. The default small-sample correction shifts the phylogenetic eye slightly relative to its raw estimate. The data-generating values were 0.8 and 0.3.