Closely related species are not independent observations. They share ancestry, so a trait value carried by one species is partly inherited from the same nodes that shaped its relatives. A phylogenetic mixed model encodes that shared history as a structured random effect: a species-level deviation whose covariance is read off the tree, so that sister taxa are expected to deviate from the regression line in the same direction.
drmTMB exposes this through the phylo()
term. This article explains where the covariance comes from, how to
write a phylo() model for a Gaussian response, and how the
same syntax carries over to non-Gaussian responses such as counts.
Everything below runs at build time on a small simulated tree, so you
can read the fitted numbers next to the values that generated them.
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 available from the same internal helper the fitting code uses. We only need it here to simulate a trait with known phylogenetic structure; you never compute it by hand for a real fit.
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: 17 checks>
#> ok: 16; 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
#> hessian_conditioning 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.0000000000001391; component=beta_sigma
#> ok
#> TRUE
#> min_eig=4.519; cond=226.6
#> 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_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.
#> 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 phylogeneticsigmagate accepts an intercept plus one independent slope at recovery grade, separately from themufield; - 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: 17 checks>
#> ok: 16; 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
#> hessian_conditioning 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=11; function=13; gradient=12
#> 158.9
#> <NA>
#> max=0.0000000000001308; 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; median_se=0.1427
#> 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 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.
#> 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: 17 checks>
#> ok: 14; notes: 3; 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
#> hessian_conditioning note
#> standard_errors_finite ok
#> standard_errors_inflated note
#> 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=34; function=39; gradient=35
#> 158.9
#> <NA>
#> max=0.000000005557; component=beta_sigma
#> ok
#> TRUE
#> min_eig=0.00000001112; cond=8260119122.
#> range=[0.1053,9483.]
#> n_inflated=1; max_se=9483.; median_se=0.1801; example=sigma:(Intercept)
#> 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.
#> 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.
#> 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.544381e-05 0.1464599 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.What phylo() fits today
This is a reader-level summary of the current capability ledger. The tier words matter: diagnostic-only means a fit/extractor smoke, point-fit recovery means simulation or oracle evidence for fitted values, and inference-ready with caveats means interval evidence exists only inside the stated design. The model map gives the wider package-level boundary.
| Question | Syntax | Status |
|---|---|---|
| Gaussian location SD |
phylo(1 | species, tree = tree) or one independent
slope in mu
|
ML intercept: inference-ready with caveats. ML intercept-plus-slope and REML intercept: point-fit recovery. |
| Gaussian residual-scale phylogeny |
sigma ~ phylo(1 | species, tree = tree) or one
independent slope |
ML intercept: point-fit recovery. ML intercept-plus-slope: inference-ready with caveats. REML pure scale route: point-fit recovery. |
| Gaussian location and residual scale together | matching labelled phylo() terms in mu and
sigma
|
Constant and first one-slope q=2 blocks have ML point-fit recovery; the constant REML q=2 block also has point-fit recovery. |
| Gaussian predictor-dependent phylogenetic SD |
phylo(1 | species, tree = tree) in mu plus
sd(species, level = "phylogenetic") ~ x
|
Fitted direct-SD route; REML admission is included in the q=1 point-fit-recovery cell. This is not a blanket interval/coverage claim. |
Poisson/NB2 q1 phylogenetic mu
intercept-plus-one-slope |
phylo(1 + x | species, tree = tree) in
mu
|
The exact unlabelled intercept plus one independent slope at recovery grade has point-fit recovery. |
| NB2 overdispersion deviations | NB2 q1 phylogenetic sigma via
phylo(1 + x | species, tree = tree)
|
The exact unlabelled intercept-plus-one-slope route has point-fit recovery; no interval or coverage promotion. |
| Gamma or lognormal location SD |
phylo(1 | species, tree = tree) in mu
|
Exact ML intercept routes have point-fit recovery; slopes,
structured sigma, REML, and intervals are outside those
cells. |
| Beta phylogenetic location-SD regression |
phylo(1 | spp_id, tree = tree) in mu plus
sd(spp_id, level = "phylogenetic") ~ x
|
Exact unlabelled q=1 ML route is inference-ready with caveats only
for the validated interior-response, 1,024-species, four-replicate
designs. Family sigma remains fixed-effect and
. |
| Student-t tail-weight deviations | Student-t q1 phylogenetic nu via
phylo(1 | species, tree = tree)
|
Diagnostic-only fit/extractor evidence; not point-fit recovery. |
| cumulative-logit location deviations | cumulative-logit q1 phylogenetic mu via
phylo(1 | species, tree = tree)
|
Diagnostic-only fit/extractor evidence; not point-fit recovery. |
| Two Gaussian means sharing phylogeny | matching phylo(1 | p | species, tree = tree) in
mu1 and mu2
|
Intercept q=2 has ML and REML point-fit recovery; the exact ML
slope-only q=2 cell is inference-ready with caveats. Read correlations
with corpairs(level = "phylogenetic"). |
| Bivariate phylogenetic location-scale blocks | matching labelled phylo() terms across
mu1, mu2, sigma1, and
sigma2
|
Exact q=2+q=2 and broader all-four cells are fitted, but their evidence ranges from diagnostic-only to point-fit recovery; they are not blanket interval/coverage claims. |
When a predictor changes the phylogenetic SD itself, or when the model contains location plus two scale submodels, continue with When variance carries signal, Part II: location-scale-scale models. For the bivariate covariance layouts and the broader parity ladder, use the structural-dependence tutorial.
Anything outside these exact rows remains unpromoted, including other
non-Gaussian families, additional slope/label combinations, public
matrix input such as phylo(1 | species, A = A),
simultaneous phylo() and spatial() layers, and
structured random effects inside residual rho12.