Use this guide when you have a dated, ultrametric phylogenetic tree and want to ask whether closely related species depart from a regression in similar ways. Closely related species are not independent observations: they share ancestry. A phylogenetic mixed model represents that history as a species-level random effect whose covariance comes from the tree, so sister taxa can tend to depart from a regression in the same direction.
drmTMB exposes this through the phylo()
term. The runnable example starts with a Gaussian trait measured
repeatedly within species. Those repeats help separate phylogenetic
variation from residual variation. Count examples appear later, but
their uncertainty support is narrower. If you only have one measurement
per species, treat separation of these two sources of variation with
particular caution.
If you are choosing between animal(),
phylo(), spatial(), and relmat(),
read the structural-dependence
overview first. When observations belong to a pair of
species drawn from two different trees, read two-tree phylogenetic
interactions.
The phylogenetic covariance
A phylo() term adds a vector of species-level location
deviations
with a mean of zero and a covariance proportional to the phylogenetic
covariance matrix
,
where is the species of observation in group . Under a Brownian motion model of trait evolution, is the height of the most recent common ancestor of species and : the longer two species have shared an evolutionary path, the more strongly their deviations covary. For an ultrametric tree (all tips equidistant from the root), the diagonal of is constant and the off-diagonal entries are the shared root-to-ancestor path lengths.
Two standard deviations appear, and they answer different questions:
- is the phylogenetic SD: how much species deviate from the fixed-effect prediction in a way that tracks the tree.
- is the residual SD: variation among observations within a species, independent of ancestry.
Their ratio is the phylogenetic signal,
the proportion of the random-plus-residual variance attributable to
phylogeny. drmTMB reports this as a derived quantity (see
below). Internally drmTMB does not invert the dense
.
It uses the Hadfield and Nakagawa (2010) sparse phylogenetic precision
(an
built from the tree’s branch lengths), which keeps the augmented-state
Laplace approximation fast as the number of species grows.
A small tree with ape::rcoal()
We simulate a tree small enough to fit in a fraction of a second.
ape::rcoal() draws a random ultrametric
coalescent tree, which is exactly the shape phylo() expects
(branch lengths present, all tips contemporaneous).
library(ape)
set.seed(2026)
n_species <- 16L
tree <- rcoal(n_species, tip.label = paste0("sp", seq_len(n_species)))
c(ultrametric = is.ultrametric(tree), n_tip = length(tree$tip.label))
#> ultrametric n_tip
#> 1 16The covariance
implied by this tree is used here only to simulate a trait with known
phylogenetic structure. You do not compute it by hand for a real fit:
provide the tree to phylo() instead.
Simulating a Gaussian trait with phylogenetic signal
We draw one species-level deviation per tip from using a Cholesky factor of , add a fixed covariate effect, and add independent residual noise. Several observations per species let the model separate the phylogenetic SD from the residual SD.
sd_phylo_true <- 0.8 # phylogenetic SD
sigma_true <- 0.3 # residual SD
n_per_species <- 6L
# One deviation per species, correlated along the tree.
u <- as.vector(t(chol(A)) %*% rnorm(n_species, sd = sd_phylo_true))
names(u) <- tree$tip.label
species <- rep(tree$tip.label, each = n_per_species)
x <- rnorm(length(species))
# y = intercept + slope * x + phylogenetic deviation + residual noise
trait <- 0.5 - 0.4 * x + u[species] + rnorm(length(species), sd = sigma_true)
dat <- data.frame(
trait = unname(trait),
x = x,
species = species
)
head(dat)
#> trait x species
#> 1 0.4962920 0.2284787 sp7
#> 2 0.0187960 1.1096678 sp7
#> 3 0.3970166 -0.4624204 sp7
#> 4 1.0612423 -1.6297240 sp7
#> 5 1.2921510 -1.0717278 sp7
#> 6 0.4733090 0.7681771 sp7Fitting 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: 19 checks>
#> ok: 18; notes: 1; warnings: 0; errors: 0
#> check status
#> optimizer_convergence ok
#> convergence_status ok
#> optimizer_budget ok
#> finite_objective ok
#> logsigma_clamp_active ok
#> fixed_gradient ok
#> sdreport_status ok
#> hessian_positive_definite ok
#> hessian_conditioning ok
#> standard_errors_finite ok
#> standard_errors_inflated ok
#> observations_per_parameter 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
#> converged
#> iterations=19; function=29; gradient=20
#> 26.37
#> <NA>
#> max=0.0000000000001891; component=beta_sigma
#> ok
#> TRUE
#> min_eig=4.519; cond=226.6
#> range=[0.03125,0.4704]
#> n_inflated=0; max_se=0.4704; reference_median=0.07706
#> n_obs=96; n_par=4; ratio=24.00
#> 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 convergence and uncertainty diagnostics are consistent with a proper interior optimum.
#> 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_sigma.
#> TMB::sdreport() completed successfully.
#> sdreport reports a positive-definite Hessian.
#> Minimum eigenvalue and condition number of TMB's sdreport() fixed-effect covariance (sdr$cov.fixed), inverted. These are a genuinely different read of the fit's conditioning than TMB's internal pdHess flag -- comparable across fits, not claimed to be numerically identical to any raw TMB gradient or Hessian quantity. This fit's Hessian conditioning is within the requested threshold.
#> All fixed-effect standard errors are finite.
#> No fixed-effect standard error is inflated relative to the others.
#> Observations per estimated parameter are at or above the small-samplenote threshold (10).
#> 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
,
alongside the random-effect and residual variances it is built from.
summary(fit)$derived[, c(
"quantity", "estimate",
"random_effect_variance", "residual_variance"
)]
#> quantity
#> derived:phylo_total_variance_share(species) phylo_total_variance_share
#> estimate random_effect_variance
#> derived:phylo_total_variance_share(species) 0.8933605 0.6076515
#> residual_variance
#> derived:phylo_total_variance_share(species) 0.07253469estimate is
.
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.
Non-Gaussian responses
The same phylo(1 | species, tree = tree) term works for
ordinary Poisson and negative-binomial (NB2) location models. The
phylogenetic deviation now acts on the linear predictor of the count
mean through the log link, so related species share a baseline
abundance. Two differences from the Gaussian case:
- Poisson has no residual
sigmaformula. NB2 instead has a modelled overdispersionsigma; its exact q1 phylogeneticsigmamodel accepts an intercept plus one independent slope, separately from themufield. Parameter recovery has been checked for this form; this does not establish the reliability of its intervals; - the simulated deviation is built from the correlation form of (unit diagonal) scaled by the phylogenetic SD, so the SD is interpretable on the log-mean scale.
set.seed(11)
# Standardise A to a correlation matrix, then scale by the phylogenetic SD.
A_cor <- A / outer(sqrt(diag(A)), sqrt(diag(A)))
sd_phylo_count <- 0.5
u_count <- as.vector(t(chol(A_cor)) %*% rnorm(n_species)) * sd_phylo_count
names(u_count) <- tree$tip.label
species_c <- rep(tree$tip.label, each = n_per_species)
x_c <- rep(seq(-1, 1, length.out = n_per_species), times = n_species)
eta <- log(3) - 0.3 * x_c + u_count[species_c] # log mean
count <- rpois(length(eta), lambda = exp(eta))
dat_count <- data.frame(count = count, x = x_c, species = species_c)
range(dat_count$count)
#> [1] 0 8A Poisson fit uses the same call with
family = poisson():
fit_pois <- drmTMB(
bf(count ~ x + phylo(1 | species, tree = tree)),
family = poisson(link = "log"),
data = dat_count
)
check_drm(fit_pois)
#> <drm_check: 19 checks>
#> ok: 18; notes: 1; warnings: 0; errors: 0
#> check status
#> optimizer_convergence ok
#> convergence_status ok
#> optimizer_budget ok
#> finite_objective ok
#> logsigma_clamp_active ok
#> fixed_gradient ok
#> sdreport_status ok
#> hessian_positive_definite ok
#> hessian_conditioning ok
#> standard_errors_finite ok
#> standard_errors_inflated ok
#> observations_per_parameter 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
#> converged
#> iterations=11; function=13; gradient=12
#> 158.9
#> <NA>
#> max=0.0000000000006847; component=log_sd_phylo
#> ok
#> TRUE
#> min_eig=4.243; cond=21.65
#> range=[0.1053,0.1801]
#> n_inflated=0; max_se=0.1801; reference_median=0.1427
#> n_obs=96; n_par=3; ratio=32.00
#> nobs=96; dropped=0
#> min=1.000
#> min=0.2530; boundary=0.0001000; term=mu.phylo(1 | species)
#> sd_targets=1; assessed_here=0
#> total_mb=0.008080; 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.2530
#> message
#> nlminb convergence code is 0.
#> Optimizer convergence and uncertainty diagnostics are consistent with a proper interior optimum.
#> Optimizer evaluation counts recorded; no eval.max or iter.max control was supplied.
#> Objective and log-likelihood are finite.
#> The log(sigma) clamp does not apply to this family.
#> Maximum absolute fixed gradient is <= 0.001; largest component is log_sd_phylo.
#> TMB::sdreport() completed successfully.
#> sdreport reports a positive-definite Hessian.
#> Minimum eigenvalue and condition number of TMB's sdreport() fixed-effect covariance (sdr$cov.fixed), inverted. These are a genuinely different read of the fit's conditioning than TMB's internal pdHess flag -- comparable across fits, not claimed to be numerically identical to any raw TMB gradient or Hessian quantity. This fit's Hessian conditioning is within the requested threshold.
#> All fixed-effect standard errors are finite.
#> No fixed-effect standard error is inflated relative to the others.
#> Observations per estimated parameter are at or above the small-samplenote threshold (10).
#> 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 finite positive fitted SD.
coef(fit_pois, "mu") # log-mean intercept near log(3) ~ 1.10, slope near -0.3
#> (Intercept) x
#> 0.5749885 -0.3875214
summary(fit_pois)$parameters # phylogenetic SD on the log-mean scale
#> component dpar term estimate
#> sd:mu:phylo(1 | species) random-effect-sd mu phylo(1 | species) 0.2529665
#> std_error minimum maximum scale
#> sd:mu:phylo(1 | species) 0.1223278 NA NA response
#> parm profile_ready profile_note
#> sd:mu:phylo(1 | species) sd:mu:phylo(1 | species) TRUE readyIf the counts are overdispersed relative to a Poisson, swap in
nbinom2(). The NB2 family adds an overdispersion (scale)
parameter while keeping the identical phylogenetic location term.
fit_nb <- drmTMB(
bf(count ~ x + phylo(1 | species, tree = tree)),
family = nbinom2(),
data = dat_count
)
check_drm(fit_nb)
#> <drm_check: 20 checks>
#> ok: 15; notes: 4; warnings: 1; errors: 0
#> check status
#> optimizer_convergence ok
#> convergence_status note
#> optimizer_budget ok
#> finite_objective ok
#> logsigma_clamp_active ok
#> fixed_gradient ok
#> sdreport_status ok
#> hessian_positive_definite ok
#> hessian_conditioning note
#> standard_errors_finite ok
#> standard_errors_inflated note
#> weakly_identified_fixed_effects warning
#> observations_per_parameter 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
#> degenerate
#> iterations=34; function=39; gradient=35
#> 158.9
#> <NA>
#> max=0.000000005557; component=beta_sigma
#> ok
#> TRUE
#> min_eig=0.00000001111; cond=8266765167.
#> range=[0.1053,9487.]
#> n_inflated=1; max_se=9487.; reference_median=0.1427; example=sigma:(Intercept)
#> n_inflated=1; max_se=9487.; reference_median=0.1427; example=sigma:(Intercept)
#> n_obs=96; n_par=4; ratio=24.00
#> nobs=96; dropped=0
#> min=0.00001544
#> min=0.2530; 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.2530
#> message
#> nlminb convergence code is 0.
#> The optimizer reported success but the likelihood geometry or uncertainty is degenerate (no reliable maximum or non-finite standard errors). multi_start cannot clear this state; is_converged() is FALSE.
#> 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_sigma.
#> TMB::sdreport() completed successfully.
#> sdreport reports a positive-definite Hessian.
#> Minimum eigenvalue and condition number of TMB's sdreport() fixed-effect covariance (sdr$cov.fixed), inverted. These are a genuinely different read of the fit's conditioning than TMB's internal pdHess flag -- comparable across fits, not claimed to be numerically identical to any raw TMB gradient or Hessian quantity. The condition number exceeds 1e+08, signalling a near-flat, weakly identified direction; a clean pdHess is necessary, not sufficient. Treat standard errors and correlations among the affected parameters with caution.
#> All fixed-effect standard errors are finite.
#> At least one standard error is finite but extremely large despite a positive-definite Hessian (pdHess = TRUE), which signals a near-flat, weakly identified direction the Hessian did not resolve -- often a correlation or SD running to a boundary. Confirm with a likelihood profile (profile()) and consider a simpler model before interpreting the affected parameter; a clean Hessian is necessary, not sufficient.
#> At least one fixed-effect standard error is pathologically large relative to the reference scale; the fit is not identified enough for plain is_converged() to read TRUE and multi_start cannot clear this.
#> Observations per estimated parameter are at or above the small-samplenote threshold (10).
#> 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 finite positive fitted SD.
summary(fit_nb)$parameters # phylogenetic SD, NB2 mean model
#> component dpar term
#> sigma distributional-scale sigma (constant)
#> sd:mu:phylo(1 | species) random-effect-sd mu phylo(1 | species)
#> estimate std_error minimum maximum scale
#> sigma 1.544417e-05 0.1465222 NA NA response
#> sd:mu:phylo(1 | species) 2.529665e-01 0.1223278 NA NA response
#> parm profile_ready profile_note
#> sigma sigma TRUE ready
#> sd:mu:phylo(1 | species) sd:mu:phylo(1 | species) TRUE readyThe phylogenetic SD is reported on the log-mean (link) scale for count families, so it is not directly comparable to the Gaussian response-scale SD; compare it instead to other log-scale effects in the same model.
A larger fit, for reference
The fits above are tiny by design. For intuition about runtime on a
more realistic tree, the chunk below sketches a 200-species fit. It is
marked eval = FALSE so the vignette never blocks on it; the
numbers in the comments are illustrative of the shape of the output, not
a benchmarked claim.
set.seed(99)
big_tree <- ape::rcoal(200, tip.label = paste0("t", 1:200))
A_big <- drmTMB:::drm_phylo_tip_covariance(big_tree)
u_big <- as.vector(t(chol(A_big)) %*% rnorm(200, sd = 0.7))
names(u_big) <- big_tree$tip.label
sp <- rep(big_tree$tip.label, each = 4L)
xb <- rnorm(length(sp))
yb <- 0.2 + 0.5 * xb + u_big[sp] + rnorm(length(sp), sd = 0.3)
big <- data.frame(y = yb, x = xb, species = sp)
fit_big <- drmTMB(
bf(y ~ x + phylo(1 | species, tree = big_tree), sigma ~ 1),
family = gaussian(),
data = big
)
summary(fit_big)$parameters
# The sparse-precision path keeps this on the order of a second on a laptop;
# cost grows roughly linearly in the number of species rather than cubically.Choose a phylogenetic route
Start with the Gaussian example in this article if you have repeated continuous measurements within species. It is the clearest route for learning what the tree term, phylogenetic SD, and residual SD each mean.
| If your question is… | Eligible route | Read or do next | Key limitation |
|---|---|---|---|
| Do related species differ around a continuous-response regression? |
phylo(1 | species, tree = tree) in the Gaussian mean
formula |
Work through the example above, then run
check_drm(). |
With little within-species replication, phylogenetic and residual SDs can be hard to distinguish. |
| Do predictors change phylogenetic variation or residual variation? | Gaussian phylo() terms in the mean and/or scale
formula |
When variance carries signal, Part II | Interpret each scale separately; a fitted term is not by itself evidence that every interval is reliable. |
| Are count means structured by ancestry? | One ordinary Poisson or NB2
phylo(1 + x | species, tree = tree) mean term |
Use the count-family guide and inspect diagnostics. | This route has a narrower evidence base than the Gaussian example; do not generalize its uncertainty to richer count structures. |
| Do two continuous responses share phylogenetic deviations? | Matching phylo() terms in two Gaussian mean
formulas |
Use Structural dependence overview. | Phylogenetic covariance is distinct from bivariate residual
correlation rho12. |
For the complete syntax boundary, including unsupported combinations,
see What can I fit today?. Do not treat a
successful fit as a general interval or coverage guarantee. The current
public routes do not cover simultaneous phylogenetic and spatial layers,
structured effects in rho12, or arbitrary additional slope
and label combinations.