Phylogenetic stable effects and temporal OU deviations
Source:vignettes/phylogenetic-temporal-effects.Rmd
phylogenetic-temporal-effects.RmdDevelopment status. This paired phylogenetic–OU model has passed its parser, dense-likelihood, method, and fixed-mean profile checks. Its first 24-fixture point-recovery gate missed the predeclared fixed-effect threshold; a replacement point-recovery study (24 contrast fixtures plus 300 independently generated trees) then met all of its predeclared point criteria. No interval-calibration study has been run. Treat the syntax below as a development workflow for checking a model and its diagnostics, not as an inference-qualified routine analysis. In particular, do not report confidence intervals from this paired model until its calibration gates are run and passed.
A comparative time-series data set can contain three different kinds of variation. Closely related species can have similar stable baselines because of shared evolutionary history. Each species can also make a short-lived departure from its baseline that persists across nearby observation times. Finally, observations have independent residual noise. The first paired phylogenetic-temporal route keeps these quantities separate.
For observation time in species , it fits
Here
is the covariance derived from the phylogenetic tree.
sd_phylo_stable describes stable between-species variation
that follows this tree. sd_temporal describes the
stationary SD of an OU departure within the same species, and
decay_temporal is its positive decay rate in the units of
the elapsed-time column. sigma is independent
observation-level noise.
The model is additive. It does not say that a temporal departure is shared between related species. In particular, for two different species, their covariance is , even when their observations occur close together in time. A separable phylogeny-by-time field is a different future model with a different scientific question.
A small comparative time-series example
The example uses genuinely irregular elapsed times. The gap from 0.5
to 2.5 is two time units; drmTMB retains that gap and does
not replace these values with ranks.
library(ape)
set.seed(20260909)
n_species <- 60L
elapsed_values <- c(0, 0.5, 2.5, 5)
tree <- rcoal(n_species, tip.label = paste0("sp", seq_len(n_species)))
# Simulate the stable phylogenetic component once per species.
A <- drmTMB:::drm_phylo_tip_covariance(tree)
stable <- drop(t(chol(A)) %*% rnorm(n_species, sd = 0.45))
names(stable) <- tree$tip.label
# Simulate independent OU departures for each species.
R_ou <- exp(-0.45 * abs(outer(elapsed_values, elapsed_values, "-")))
ou_by_species <- lapply(tree$tip.label, function(species) {
drop(t(chol(R_ou)) %*% rnorm(length(elapsed_values), sd = 0.55))
})
names(ou_by_species) <- tree$tip.label
dat <- do.call(rbind, lapply(tree$tip.label, function(species) {
treatment <- rep(c(-0.5, 0.5), length.out = length(elapsed_values))
data.frame(
species = species,
elapsed_days = elapsed_values,
treatment = treatment,
y = 1 + 0.4 * treatment + stable[[species]] +
ou_by_species[[species]] + rnorm(length(elapsed_values), sd = 0.35)
)
}))
dat$species <- factor(dat$species, levels = tree$tip.label)
head(dat)
#> species elapsed_days treatment y
#> 1 sp33 0.0 -0.5 0.1017484
#> 2 sp33 0.5 0.5 -0.8825828
#> 3 sp33 2.5 -0.5 0.2274549
#> 4 sp33 5.0 0.5 -0.5410026
#> 5 sp19 0.0 -0.5 -1.7771153
#> 6 sp19 0.5 0.5 0.3865659Fit the phylogenetic stable intercept and OU process with the same
species column. The tree object is supplied by
name, so its tips can be checked against the observed species before the
response is analysed.
fit <- drmTMB(
bf(
y ~ treatment +
phylo(1 | species, tree = tree) +
temporal(1 | species, time = elapsed_days, structure = "ou"),
sigma ~ 1
),
data = dat,
family = gaussian(),
REML = FALSE
)
summary(fit)
#> <summary.drmTMB>
#> estimator: ML
#> standard errors: unavailable; OU coefficient table is point estimates only (OU
#> Wald inference is deferred pending the inherited AR1 calibration prerequisite)
#> estimate std_error std_error.status
#> mu:(Intercept) 0.7636514 NA temporal_wald_unqualified
#> mu:treatment 0.3870943 NA temporal_wald_unqualified
#> sigma:(Intercept) -1.1208776 NA temporal_wald_unqualified
#> Distributional, random-effect, scale, and correlation parameters:
#> component dpar term
#> sigma distributional-scale sigma (constant)
#> sd:mu:sd_temporal random-effect-sd mu sd_temporal
#> sd:mu:sd_phylo_stable random-effect-sd mu sd_phylo_stable
#> decay:temporal:decay_temporal temporal-decay temporal decay_temporal
#> estimate scale
#> sigma 0.3259936 response
#> sd:mu:sd_temporal 0.5429868 response
#> sd:mu:sd_phylo_stable 0.4610330 response
#> decay:temporal:decay_temporal 0.4424672 response
#> Temporal fit: sigma, SD, persistence, decay, and lag-correlation rows are point
#> estimates only; their standard errors are withheld.
#> logLik: -220.6
#> convergence_status: converged
#> convergence: 0
check_drm(fit)
#> <drm_check: 23 checks>
#> ok: 20; notes: 3; 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
#> mle_existence ok
#> observations_per_parameter ok
#> dropped_rows ok
#> positive_scale ok
#> random_effect_sd_boundary ok
#> interval_reliability_scope note
#> temporal_boundary ok
#> temporal_mean_wald note
#> temporal_mean_profile note
#> fixed_effect_design_size ok
#> phylo_mu_replication ok
#> phylo_mu_diagnostics ok
#> value
#> 0
#> converged
#> iterations=24; function=29; gradient=25
#> 220.6
#> <NA>
#> max=0.000000000005493; component=log_sd_phylo
#> ok
#> TRUE
#> min_eig=4.209; cond=70.21
#> range=[0.05818,0.4383]
#> n_inflated=0; max_se=0.4383; reference_median=0.2823
#> finite
#> n_obs=240; n_par=6; ratio=40.00
#> nobs=240; dropped=0
#> min=0.3260
#> min=0.4610; boundary=0.0001000; term=mu.sd_phylo_stable
#> sd_targets=2; assessed_here=0
#> <NA>
#> unavailable; reason=calibration_deferred
#> available_for_this_fit; calibration=unqualified
#> total_mb=0.03625; max_cols=2; largest=mu; largest_class=matrix; largest_density=1.000
#> min_species_n=4
#> group=species; n_species=60; min_species_n=4; phylo_sd=0.4610; sd_ratio=1.414
#> 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 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.
#> The screened fixed-effect geometry is consistent with a finite maximum-likelihood estimate.
#> 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 2 random-effect standard-deviation targets. `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.
#> Residual sigma and the temporal persistence or decay are away from their boundaries.
#> OU mean-coefficient Wald intervals are intentionally unavailable while the inherited AR1 calibration prerequisite remains unresolved.
#> Temporal mean-coefficient profile intervals are available for this fit. Their general coverage calibration remains unresolved; do not treat this fit-level status as a coverage claim.
#> 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.The paired route is intentionally narrow: it needs a univariate
Gaussian ML model, sigma ~ 1, finite numeric elapsed times,
one response for every species–time key, at least three observed
species, at least two distinct times per species, and at least three
positive distinct lags across the data. It rejects an ordinary
(1 | species) intercept, different IDs for
phylo() and temporal(), temporal slopes, a
second structured term, non-Gaussian families, REML, non-unit weights,
forecasts, and newdata prediction.
Read the three sources of variation separately
The fit keeps the phylogenetic and temporal components distinct in
its structured metadata. The exact labels are intended to prevent a
temporal SD from being read as evolutionary variation, or vice versa.
The parameters table returned by summary()
reports each component under that label in its term
column.
parameter_table <- summary(fit)$parameters
component_estimates <- setNames(parameter_table$estimate, parameter_table$term)
component_estimates[c("sd_temporal", "sd_phylo_stable")]
#> sd_temporal sd_phylo_stable
#> 0.5429868 0.4610330
component_estimates["decay_temporal"]
#> decay_temporal
#> 0.4424672
sigma(fit)[1]
#> [1] 0.3259936The data were simulated with sd_phylo_stable = 0.45,
sd_temporal = 0.55, decay_temporal = 0.45,
sigma = 0.35, intercept 1, and treatment effect 0.4. One
compact vector collects the point estimates for comparison:
key_estimates <- c(
logLik = as.numeric(logLik(fit)),
coef(fit, "mu"),
component_estimates[c("sd_temporal", "sd_phylo_stable", "decay_temporal")],
sigma = unname(sigma(fit)[1])
)
round(key_estimates, 4)
#> logLik (Intercept) treatment sd_temporal sd_phylo_stable
#> -220.6328 0.7637 0.3871 0.5430 0.4610
#> decay_temporal sigma
#> 0.4425 0.3260With a smaller panel the residual SD can run to zero while the OU
term absorbs it. With this same generator and 16 species,
sigma falls to about 3e-5 although the optimizer, Hessian
and standard errors look regular; check_drm() reports this
in its temporal_boundary row and
convergence_status(fit) returns "boundary". At
8 species the optimizer also stops with false convergence. Such fits
should not be read component by component.
fitted(fit) and residuals(fit) use the
estimated stable and temporal conditional deviations for these observed
rows. simulate(fit) draws a new phylogenetic stable vector
and a new independent OU path for each species;
simulate(fit, re.form = NULL) instead retains the fitted
conditional latent components. These are in-sample tools only.
head(fitted(fit))
#> [1] -0.26288067 -0.33686029 -0.02288582 -0.31882141 -1.08574984 -0.05657510
head(residuals(fit))
#> [1] 0.3646291 -0.5457225 0.2503407 -0.2221812 -0.6913655 0.4431410
head(simulate(fit, nsim = 2))
#> sim_1 sim_2
#> 1 -0.48593847 1.7462285
#> 2 0.07711415 1.2323020
#> 3 -1.33625009 1.1059569
#> 4 0.18325630 1.0324259
#> 5 0.56111522 0.5540526
#> 6 1.30999371 1.5748198
head(simulate(fit, nsim = 2, re.form = NULL))
#> `simulate()` now follows the lme4 `re.form` convention.
#> • `re.form = NA` (the default) or `~0` draws fresh random effects.
#> • `re.form = NULL` holds the fitted random effects fixed.
#> ℹ Earlier drmTMB versions gave the explicit values the reverse meaning; check
#> code written for them.
#> This message is displayed once per session.
#> sim_1 sim_2
#> 1 -0.41434102 -0.1713432
#> 2 -0.09406648 -0.5407837
#> 3 0.26898707 -0.4104471
#> 4 -0.13101827 -0.2223853
#> 5 -0.96448534 -0.6511343
#> 6 -0.43419289 -0.1236500What can be reported now
This development slice has an independently checked fixed-mean
likelihood profile calculation, but its recovery evidence is not yet
sufficient to qualify an interval workflow. For now, use the example to
inspect the fitted covariance components and the diagnostics; do not use
confint(), Wald covariance, bootstrap intervals, decay or
variance-component intervals, forecasts, or newdata
prediction as a reported paired-model result.
check_drm(fit) is a necessary diagnostic but cannot turn an
unqualified recovery result into calibrated inference.
The existing temporal AR1 and OU guide remains the supported guide for independent-series temporal models. A later phylogeny-by-OU field will require a separate formula, likelihood, simulation study, and reader guide because it represents a different covariance model.