Temporal AR1, OU and Toeplitz effects
Status — Tested (Gaussian mean, ML)
Mirrors drmTMB's Temporal AR1, OU, and Toeplitz effects article (vignettes/temporal-random-effects.Rmd). In DRModels.jl today: one intercept-only temporal(1 | id, time, ar1) or temporal(1 | id, time, ou) term on the Gaussian mean, with sigma ~ 1, fitted by ML, optionally alongside an ordinary (1 | id) on the same id, or (OU only) alongside a phylogenetic stable intercept phylo(1 | species) on the same grouping; or a homogeneous Toeplitz temporal(1 | id, time, homtoep) covariance on a complete, equally spaced panel. Other families, temporal terms on sigma, slopes, REML and other structured terms are refused with an error. This page fits drmTMB's own two datasets and reproduces its numbers (table below).
Attribution
The explanatory prose on this page is adapted from the drmTMB article Temporal AR1 and OU random effects (vignettes/temporal-random-effects.Rmd; its Toeplitz section from the development version of that article, Temporal AR1, OU, and Toeplitz effects) by its copyright holder, Shinichi Nakagawa, and is reused here under the MIT licence. The code, the Julia-specific text and the comparisons are new.
Repeated measurements can vary for three different reasons. Individuals or sites may have stable baseline differences, their deviations may persist through time, and each observation has independent residual noise. Gaussian temporal AR1 and OU models separate those components while estimating the effects of treatments or other predictors.
For observation
where
Writing the term in Julia
drmTMB writes the term with keyword arguments, temporal(1 | site, time = occasion, structure = "ar1"). Julia's @formula cannot carry keyword arguments or strings inside a formula, so DRModels.jl uses the same three pieces by position, with a bare structure name:
| drmTMB (R) | DRModels.jl (Julia) |
|---|---|
temporal(1 | site, time = occasion, structure = "ar1") | temporal(1 | site, occasion, ar1) |
temporal(1 | site, time = elapsed_days, structure = "ou") | temporal(1 | site, elapsed_days, ou) |
If you call DRModels.jl from R through drm_bridge (see Coming from R?), you can keep drmTMB's keyword spelling: the bridge translates it to the positional form (shown at the end of this page).
A repeated-site example
The data are the ones drmTMB's article simulates: 48 sites visited on the irregular integer occasions 0, 1, 3, 4, 6 and 7, with a treatment that alternates between visits. The gap between occasions 1 and 3 remains a two-occasion gap when the model is fitted. The simulation is R code (set.seed(20260908)), so rather than re-simulate it in Julia we read the exact data frame it produced, which ships with the package's test fixtures.
using DRModels
# A tiny CSV reader, so this page needs nothing beyond DRModels.
function read_fixture(file)
lines = readlines(pkgdir(DRModels, "test", "fixtures", "temporal", file))
header = split(lines[1], ',')
rows = split.(lines[2:end], ',')
return name -> [r[findfirst(==(name), header)] for r in rows]
end
col = read_fixture("vignette_ar1_sites.csv")
dat = (y = parse.(Float64, col("y")),
treatment = parse.(Float64, col("treatment")),
occasion = parse.(Int, col("occasion")), # AR1 needs integer occasions
site = String.(col("site")))
(n = length(dat.y), sites = length(unique(dat.site)), occasions = sort(unique(dat.occasion)))(n = 288, sites = 48, occasions = [0, 1, 3, 4, 6, 7])Fit an ordinary random intercept and the temporal process together. Both terms use site; the first captures stable site-to-site differences and the second captures deviations through time within each site. Here treatment compares treated and untreated visits within sites after accounting for a linear occasion trend.
fit = drm(bf(@formula(y ~ treatment + occasion + (1 | site) +
temporal(1 | site, occasion, ar1)),
@formula(sigma ~ 1)),
Gaussian(); data = dat)
coeftable(fit)| Estimate | Std.Error | z | Pr(>|z|) | Lower 95% | Upper 95% | |
|---|---|---|---|---|---|---|
| mu: (Intercept) | 1.08377 | 0.127417 | 8.51 | <1e-16 | 0.834036 | 1.3335 |
| mu: treatment | 0.495548 | 0.072682 | 6.82 | <1e-11 | 0.353093 | 0.638002 |
| mu: occasion | 0.0638532 | 0.021044 | 3.03 | 0.0024 | 0.0226076 | 0.105099 |
| sigma: (Intercept) | -0.656984 | 0.201665 | -3.26 | 0.0011 | -1.05224 | -0.261728 |
| resd: site_iid | -0.575523 | 0.227628 | -2.53 | 0.0115 | -1.02167 | -0.12938 |
| resd: site | -0.542066 | 0.176301 | -3.07 | 0.0021 | -0.88761 | -0.196521 |
| temporal_phi: temporal(1 | site, time = occasion, structure = "ar1") | 0.815976 | 0.454056 | 1.80 | 0.0723 | -0.0739579 | 1.70591 |
The temporal process and the variance components are collected by temporal_parameters: the process SD (sd), the persistence (phi, or decay for OU), the ordinary-intercept SD (sd_iid) and the residual SD (sigma), all on their natural scales.
tp = temporal_parameters(fit)(label = "temporal(1 | site, time = occasion, structure = \"ar1\")", structure = :ar1, sd = 0.5815455889308447, phi = 0.6728736494884742, decay = nothing, sd_iid = 0.5624106363186523, sd_phylo = nothing, sigma = 0.5184125952024274, cor = nothing, pac = nothing)The mean-model rows report the intercept, treatment and occasion effects. Read a positive persistence estimate as deviations tending to continue in the same direction over adjacent occasions; read a negative value as alternating deviations.
Same numbers as drmTMB
The drmTMB article prints its key estimates rounded to four decimals. The reference values below are drmTMB's full-precision output on the same data, read from the parity cell test/parity/temporal/vignette-ar1-ri/ (generated by drmTMB 0.7.1 under R 4.6.1; the drmTMB commit is recorded in that cell's expected.meta.toml).
using TOML
parity_cell(name) = TOML.parsefile(pkgdir(DRModels, "test", "parity", "temporal",
name, "expected.toml"))
# drmTMB's numbers for one cell, as (quantity => value) pairs.
function drmtmb_estimates(ref)
t = ref["temporal"]
persistence = haskey(t, "phi") ? "phi" => t["phi"] : "decay" => t["decay"]
return ["logLik" => ref["fit"]["loglik"],
[replace(k, "mu_" => "") => v for (k, v) in ref["coef"]]...,
"sd_site" => t["sd_iid"], "sd_temporal" => t["sd"],
persistence, "sigma" => t["sigma"]]
end
# The same quantities from a DRModels.jl fit, matched by name.
function julia_estimates(fit)
tp = temporal_parameters(fit)
names = Dict(fit.coefnames)[:mu]
β = coef(fit, :mu)
return Dict("logLik" => loglik(fit), (names .=> β)...,
"sd_site" => tp.sd_iid, "sd_temporal" => tp.sd,
"phi" => tp.phi, "decay" => tp.decay, "sigma" => tp.sigma)
end
compare(ref, fit) = [(quantity = k, drmTMB = round(r; digits = 4),
DRModels = round(julia_estimates(fit)[k]; digits = 4),
abs_diff = abs(julia_estimates(fit)[k] - r))
for (k, r) in drmtmb_estimates(ref)]
ar1_ref = parity_cell("vignette-ar1-ri")
compare(ar1_ref, fit)8-element Vector{@NamedTuple{quantity::String, drmTMB::Float64, DRModels::Float64, abs_diff::Float64}}:
(quantity = "logLik", drmTMB = -345.3448, DRModels = -345.3448, abs_diff = 9.379164112033322e-12)
(quantity = "occasion", drmTMB = 0.0639, DRModels = 0.0639, abs_diff = 5.98965321785272e-14)
(quantity = "(Intercept)", drmTMB = 1.0838, DRModels = 1.0838, abs_diff = 1.6513457268274578e-12)
(quantity = "treatment", drmTMB = 0.4955, DRModels = 0.4955, abs_diff = 5.302702721365904e-12)
(quantity = "sd_site", drmTMB = 0.5624, DRModels = 0.5624, abs_diff = 1.1077472272802424e-11)
(quantity = "sd_temporal", drmTMB = 0.5815, DRModels = 0.5815, abs_diff = 1.9129475781198835e-11)
(quantity = "phi", drmTMB = 0.6729, DRModels = 0.6729, abs_diff = 5.0291992792494966e-11)
(quantity = "sigma", drmTMB = 0.5184, DRModels = 0.5184, abs_diff = 5.639155808978558e-12)Every difference is far below the fourth decimal: the two packages reach the same maximum-likelihood optimum.
Fitted values and conditional effects
One convention differs from drmTMB. In DRModels.jl, fitted(fit) (and predict) is the population-level mean fitted() adds the estimated site and temporal deviations. Those conditional deviations are available from ranef: ranef(fit)[:site] holds the temporal deviation for each row (in data order) and ranef(fit)[:site_iid] the stable intercept for each site.
conditional = fitted(fit) .+ ranef(fit)[:site] .+
ranef(fit)[:site_iid][indexin(dat.site, unique(dat.site))]
drmtmb_fitted = Float64.(ar1_ref["conditional"]["fitted"]) # drmTMB's fitted(fit)
(rows = length(conditional), first_rows = first(conditional, 3),
drmTMB_first_rows = first(drmtmb_fitted, 3),
max_abs_diff = maximum(abs.(conditional .- drmtmb_fitted)))(rows = 288, first_rows = [0.3428133449902613, 1.2082311574622335, 1.0921930259451367], drmTMB_first_rows = [0.34281340073803396, 1.2082311348353578, 1.0921930170833316], max_abs_diff = 1.490183669794476e-7)All 288 conditional values agree with drmTMB's fitted(fit) to the max_abs_diff shown. That agreement is looser than the agreement of the estimates because of drmTMB's side: for this fit, the modes drmTMB stores are not exactly the modes at its own reported optimum, while DRModels.jl's modes are exact. The residuals below inherit the same small difference.
Standardised residuals use the same conditional convention as drmTMB's Pearson and quantile residuals. residuals(fit; type = :quantile) is check_drm(fit).temporal_boundary first: at the residual-SD boundary (sigma_ratio, described below) the temporal path absorbs the data, so as residuals(fit) stays fitted.)
z = residuals(fit; type = :quantile)
drmtmb_pearson = Float64.(ar1_ref["residuals"]["pearson"]) # drmTMB's residuals(fit, type = "pearson")
(rms = sqrt(sum(abs2, z) / length(z)), first_rows = first(z, 3),
drmTMB_first_rows = first(drmtmb_pearson, 3),
max_abs_diff = maximum(abs.(z .- drmtmb_pearson)))(rms = 0.7141638575236557, first_rows = [-1.007915739730797, 0.2848137210449365, -0.3475001159900366], drmTMB_first_rows = [-1.007915847255362, 0.28481376468830066, -0.34750009889214034], max_abs_diff = 2.8743102120643016e-7)simulate(fit) draws a new realization from the fitted model: fresh stable and temporal effects plus residual noise, as drmTMB's default simulate().
using Random
first(simulate(fit; rng = Xoshiro(1)), 6)6-element Vector{Float64}:
0.8166877248213305
0.8899369928470595
1.9758120542419246
3.2502215723111467
2.449978066268154
3.454507480553618What the intervals cover
drmTMB deliberately exposes Wald intervals only for the mean coefficients of an AR1 fit, using the full observed marginal-likelihood Hessian, and treats the process SD, ordinary-intercept SD, persistence or decay and residual SD as point estimates. Read DRModels.jl the same way:
wald = confint(fit; parm = :mu)3-element Vector{@NamedTuple{param::Symbol, coef::String, estimate::Float64, lower::Float64, upper::Float64}}:
(param = :mu, coef = "(Intercept)", estimate = 1.0837697434692317, lower = 0.8340362186970499, upper = 1.3335032682414134)
(param = :mu, coef = "treatment", estimate = 0.49554761767492894, lower = 0.35309349755073327, upper = 0.6380017377991246)
(param = :mu, coef = "occasion", estimate = 0.06385316501431108, lower = 0.02260759824588969, upper = 0.10509873178273248)The same Wald limits from drmTMB (confint(fit, method = "wald"), stored in the parity cell), side by side:
# Pair DRModels.jl interval rows with drmTMB rows ("fixef:mu:<coef>").
function compare_ci(julia_rows, drmtmb_rows)
map(drmtmb_rows) do r
coefname = replace(r["parm"], "fixef:mu:" => "")
j = only(filter(x -> x.coef == coefname, julia_rows))
(coef = coefname, lower = j.lower, upper = j.upper,
drmTMB_lower = r["lower"], drmTMB_upper = r["upper"],
max_abs_diff = max(abs(j.lower - r["lower"]), abs(j.upper - r["upper"])))
end
end
compare_ci(wald, ar1_ref["wald"])3-element Vector{@NamedTuple{coef::String, lower::Float64, upper::Float64, drmTMB_lower::Float64, drmTMB_upper::Float64, max_abs_diff::Float64}}:
(coef = "(Intercept)", lower = 0.8340362186970499, upper = 1.3335032682414134, drmTMB_lower = 0.8340362181302124, drmTMB_upper = 1.3335032688049484, max_abs_diff = 5.668374658540642e-10)
(coef = "treatment", lower = 0.35309349755073327, upper = 0.6380017377991246, drmTMB_lower = 0.3530934975001707, drmTMB_upper = 0.6380017378390818, max_abs_diff = 5.056255414359612e-11)
(coef = "occasion", lower = 0.02260759824588969, upper = 0.10509873178273248, drmTMB_lower = 0.02260759823040296, drmTMB_upper = 0.105098731798339, max_abs_diff = 1.5606529957246096e-11)DRModels.jl's vcov, stderror and coeftable also print Wald standard errors for the variance and persistence coordinates (on their working scales), because every DRModels.jl route does. drmTMB does not report them, and their calibration has not been established in either package; do not use them as intervals for the temporal parameters.
Likelihood profiles are available for the mean coefficients of both AR1 and OU fits. Use them when reporting a treatment or other regression effect:
compare_ci(confint(fit; method = :profile, parm = :mu => "treatment"), ar1_ref["profile"])1-element Vector{@NamedTuple{coef::String, lower::Float64, upper::Float64, drmTMB_lower::Float64, drmTMB_upper::Float64, max_abs_diff::Float64}}:
(coef = "treatment", lower = 0.35213232259395677, upper = 0.6391770678043932, drmTMB_lower = 0.3521366880603612, drmTMB_upper = 0.6391744494852939, max_abs_diff = 4.365466404421259e-6)Each package locates profile endpoints with its own numerical root search, so expect profile limits to agree less tightly than the Wald limits; the differences on this page are measured, not asserted.
No coverage is claimed for AR1 intervals of either kind. drmTMB's first AR1 calibration pilot found a residual-SD boundary in a primary short-series case, so its AR1 profile and Wald intervals carry no coverage claim, and DRModels.jl makes none either.
drmTMB refuses every interval for sigma, the variance components and the persistence or decay: Wald, profile and bootstrap. DRModels.jl returns them as an extension, without a calibration claim, and does not warn that they are uncalibrated. Wald confint(fit) and bootstrap_ci cover every coordinate, and confint(fit; method = :profile, parm = :sigma) (or :resd, :temporal_phi, :temporal_decay) profiles one block and errors when it cannot locate an endpoint. Near a variance boundary a profile limit can be infinite and bootstrap refits warn that the Hessian is singular, so treat these intervals as diagnostics. Forecasting and newdata prediction of the temporal process are not available.
Supply a finite integer time column for AR1, keep its real gaps, and make sure every site–occasion pair is unique; if a site was measured more than once at an occasion, aggregate those records before fitting. Rows may be in any order. An AR1-only fit omits (1 | site) but keeps the same temporal() term.
Free correlation by discrete lag with homogeneous Toeplitz
Use homogeneous Toeplitz covariance when every site is measured at the same complete, equally spaced set of discrete occasions and the scientific question is whether correlation departs from AR1's exponential pattern. It estimates a separate correlation for each lag: visits one occasion apart share cor_lag1, visits two occasions apart share cor_lag2, and so on. This is more flexible than AR1, so it needs a common panel rather than the irregular or incomplete schedules accepted by AR1 and OU. In Julia the structure name is homtoep:
| drmTMB (R) | DRModels.jl (Julia) |
|---|---|
temporal(1 | site, time = occasion, structure = "homtoep") | temporal(1 | site, occasion, homtoep) |
drmTMB's article simulates 80 sites, each visited at the same six equally spaced occasions (set.seed(20261002)), with a treatment that alternates between visits. Its lag correlations (0.60, 0.45, 0.40, 0.30, 0.20) decay more slowly than an AR1 pattern with the same first-lag value would (0.60, 0.36, 0.22, 0.13, 0.08), and the total SD is 0.80. We read the exact data it produced.
col = read_fixture("vignette_homtoep_sites.csv")
panel = (y = parse.(Float64, col("y")),
treatment = parse.(Float64, col("treatment")),
site = String.(col("site")),
occasion = parse.(Int, col("occasion"))) # integer occasions
toeplitz_fit = drm(bf(@formula(y ~ treatment + temporal(1 | site, occasion, homtoep)),
@formula(sigma ~ 1)),
Gaussian(); data = panel)
toeplitz_tp = temporal_parameters(toeplitz_fit)
(sigma = toeplitz_tp.sigma, cor_lag = round.(toeplitz_tp.cor; digits = 3))(sigma = 0.8005278490582741, cor_lag = [0.642, 0.473, 0.371, 0.238, 0.178])Here sigma is the total within-site SD of the Toeplitz covariance. The model does not separately estimate a temporal-process SD, an independent residual SD or an ordinary (1 | site) intercept, because those components are not separately identifiable when the correlation at every lag is free; both packages refuse (1 | site) beside homtoep. The lag correlations are estimated through their partial autocorrelations (coef(toeplitz_fit, :temporal_pac), on the atanh scale), which keeps every fitted correlation matrix positive definite. There are no latent temporal states: fitted is simulate draws each site's six values jointly from the fitted
The same model fitted by drmTMB to the same data (parity cell test/parity/temporal/vignette-homtoep/; drmTMB's article prints these estimates rounded to four decimals):
toep_ref = parity_cell("vignette-homtoep")
toep_names = Dict(toeplitz_fit.coefnames)[:mu]
[(quantity = q, drmTMB = round(r; digits = 4), DRModels = round(j; digits = 4),
abs_diff = abs(j - r))
for (q, r, j) in [("logLik", toep_ref["fit"]["loglik"], loglik(toeplitz_fit)),
[(k, v, coef(toeplitz_fit, :mu)[findfirst(==(replace(k, "mu_" => "")), toep_names)])
for (k, v) in toep_ref["coef"]]...,
("sigma", toep_ref["temporal"]["sigma"], toeplitz_tp.sigma),
[("cor_lag$m", r, toeplitz_tp.cor[m])
for (m, r) in enumerate(toep_ref["temporal"]["cor"])]...]]9-element Vector{@NamedTuple{quantity::String, drmTMB::Float64, DRModels::Float64, abs_diff::Float64}}:
(quantity = "logLik", drmTMB = -465.3119, DRModels = -465.3119, abs_diff = 9.606537787476555e-12)
(quantity = "mu_(Intercept)", drmTMB = 1.0614, DRModels = 1.0614, abs_diff = 1.177502539917441e-12)
(quantity = "mu_treatment", drmTMB = 0.362, DRModels = 0.362, abs_diff = 1.2068124277675452e-13)
(quantity = "sigma", drmTMB = 0.8005, DRModels = 0.8005, abs_diff = 1.684763439868675e-12)
(quantity = "cor_lag1", drmTMB = 0.6424, DRModels = 0.6424, abs_diff = 5.613287612504791e-13)
(quantity = "cor_lag2", drmTMB = 0.4732, DRModels = 0.4732, abs_diff = 1.1403655797437295e-12)
(quantity = "cor_lag3", drmTMB = 0.3712, DRModels = 0.3712, abs_diff = 1.5349943538467414e-12)
(quantity = "cor_lag4", drmTMB = 0.238, DRModels = 0.238, abs_diff = 2.7781110745195292e-12)
(quantity = "cor_lag5", drmTMB = 0.1779, DRModels = 0.1779, abs_diff = 3.82233134033072e-12)drmTMB qualifies only likelihood-profile intervals for mean regression effects of this model (its article computes this one with profile_engine = "tmbprofile"):
compare_ci(confint(toeplitz_fit; method = :profile, parm = :mu => "treatment"), toep_ref["profile"])1-element Vector{@NamedTuple{coef::String, lower::Float64, upper::Float64, drmTMB_lower::Float64, drmTMB_upper::Float64, max_abs_diff::Float64}}:
(coef = "treatment", lower = 0.2886344333990577, upper = 0.4353838574633656, drmTMB_lower = 0.2886354358001363, drmTMB_upper = 0.43538255873387255, max_abs_diff = 1.2987294930599802e-6)A retained 4,000-fit drmTMB campaign qualified these mean-effect profiles in three predeclared 80-site, six-occasion panels: AR1-shaped, non-exponential, and negative first-lag correlation patterns. The 20-site stress panel had lower intercept coverage, so this is evidence for those primary panel designs, not a coverage claim for every Toeplitz analysis, and it was obtained with drmTMB: DRModels.jl reproduces the same likelihood, but no separate calibration study has been run on its intervals. Both packages refuse Wald covariance and intervals (vcov, stderror, confint(fit), predict(...; se = true)) and profile intervals for sigma or the lag correlations. confint(fit; method = :profile) without parm profiles the mean coefficients only, and coeftable prints the standard-error columns as NaN. profile_curve refuses the sigma and partial-autocorrelation coordinates, as drmTMB's profile() does, and so does parameter_surface, which has no drmTMB counterpart. drmTMB also refuses the bootstrap for temporal models. As an extension, DRModels.jl's bootstrap_ci, bootstrap_summary and bootstrap_result return percentile intervals for the mean coefficients only. They carry no calibration claim and do not warn that they are uncalibrated, so the caution above applies to them too. If elapsed gaps are genuinely irregular, use OU instead.
Standardised residuals account for the correlation within a site: residuals(toeplitz_fit; type = :quantile) returns the whitened type = "quantile" residuals for this model are not whitened: they are :quantile residuals with drmTMB's type = "pearson", not its type = "quantile".
The panel rules are drmTMB's: integer occasions, at least 3 and at most 12 common occasions, equally spaced, every site observing all of them, and at least as many sites as occasions (the K − 1 free lag correlations need that many independent series to be estimable at all; this is a floor, not a design recommendation). Missing responses are handled as drmTMB handles them. The rows with a missing response are dropped first, and then the panel rules apply to the rows that remain. A site that loses one occasion is therefore refused by name as incomplete, and is never silently dropped. An occasion that is missing for every site leaves a smaller common panel, which still fits if it stays equally spaced.
Irregular elapsed time with OU
Use an OU process when elapsed intervals carry meaning, such as visits at 0, 0.5, 2.5 and 5 days. It keeps the same three-way variance separation but estimates a positive decay rate, so the correlation over a gap set.seed(20260909), 24 sites).
col = read_fixture("vignette_ou_sites.csv")
dat_irregular = (y = parse.(Float64, col("y")),
treatment = parse.(Float64, col("treatment")),
site = String.(col("site")),
elapsed_days = parse.(Float64, col("elapsed_days")))
ou_fit = drm(bf(@formula(y ~ treatment + (1 | site) +
temporal(1 | site, elapsed_days, ou)),
@formula(sigma ~ 1)),
Gaussian(); data = dat_irregular)
ou_tp = temporal_parameters(ou_fit)(label = "temporal(1 | site, time = elapsed_days, structure = \"ou\")", structure = :ou, sd = 0.6784555417098643, phi = nothing, decay = 0.49279374272145376, sd_iid = 0.5976810254929501, sd_phylo = nothing, sigma = 0.1901118149997702, cor = nothing, pac = nothing)ou_ref = parity_cell("vignette-ou-ri")
compare(ou_ref, ou_fit)7-element Vector{@NamedTuple{quantity::String, drmTMB::Float64, DRModels::Float64, abs_diff::Float64}}:
(quantity = "logLik", drmTMB = -102.5806, DRModels = -102.5806, abs_diff = 1.4210854715202004e-14)
(quantity = "(Intercept)", drmTMB = 0.8494, DRModels = 0.8494, abs_diff = 1.2611023336717153e-12)
(quantity = "treatment", drmTMB = 0.3713, DRModels = 0.3713, abs_diff = 2.6015856136041293e-12)
(quantity = "sd_site", drmTMB = 0.5977, DRModels = 0.5977, abs_diff = 1.86204385244082e-11)
(quantity = "sd_temporal", drmTMB = 0.6785, DRModels = 0.6785, abs_diff = 1.7819856701351e-11)
(quantity = "decay", drmTMB = 0.4928, DRModels = 0.4928, abs_diff = 5.123956814401254e-11)
(quantity = "sigma", drmTMB = 0.1901, DRModels = 0.1901, abs_diff = 1.0602962952077633e-11)For OU, elapsed_days must be finite and numeric, genuine gaps are retained, and duplicate site–time records must be aggregated before fitting. decay is the positive rate in the units of elapsed_days: if you change days to hours, the numerical rate changes by the inverse factor while the fitted correlation over a physical gap does not.
drmTMB withholds OU Wald inference altogether and recommends a profile interval for a mean effect:
compare_ci(confint(ou_fit; method = :profile, parm = :mu => "treatment"), ou_ref["profile"])1-element Vector{@NamedTuple{coef::String, lower::Float64, upper::Float64, drmTMB_lower::Float64, drmTMB_upper::Float64, max_abs_diff::Float64}}:
(coef = "treatment", lower = 0.18352954458005696, upper = 0.5452422370344282, drmTMB_lower = 0.18353339959853773, drmTMB_upper = 0.5452361650682405, max_abs_diff = 6.07196618773731e-6)drmTMB attaches a condition to reporting this interval: report it only when check_drm(ou_fit) raises no temporal_mean_profile warning and the interval's conf.status shows no problem; finite endpoints alone are not calibration evidence. On these data drmTMB's check gives:
ou_ref["check_drm"]Dict{String, Any} with 2 entries:
"temporal_mean_profile_status" => "note"
"temporal_mean_profile_value" => "available_for_this_fit; calibration=unqual…The status is a note (profile intervals are available for this fit, with calibration unqualified), not a warning, so by drmTMB's rule the interval may be reported, without a coverage claim.
DRModels.jl has no equivalent check yet: its check_drm does not assess temporal profile intervals. Before reporting a DRModels.jl OU profile interval, run drmTMB's check on the same model, or treat the interval as uncalibrated.
drmTMB qualified fixed-mean profile intervals for OU in a retained 3,000-fit campaign across three predeclared scenarios (80 sites measured at 6 or 12 irregular occasions, with and without stable site intercepts). That result does not establish general temporal coverage, it was obtained with drmTMB, and this page does not extend it: DRModels.jl reproduces the same likelihood, but no separate calibration study has been run on its intervals. drmTMB refuses all OU Wald intervals, and every interval for sigma, the decay and the variance components. DRModels.jl returns them as described for AR1 above: as an extension, without a calibration claim, and without a warning that they are uncalibrated.
Phylogenetic stable effects with OU deviations
Attribution
This section adapts the drmTMB article Phylogenetic stable effects and temporal OU deviations (vignettes/phylogenetic-temporal-effects.Rmd, drmTMB development version) by its copyright holder, Shinichi Nakagawa, under the MIT licence.
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 paired phylogenetic–temporal model keeps these quantities separate. For observation time
with
The model is additive. It does not say that a temporal departure is shared between related species: for two different species the covariance is
The data are those simulated by drmTMB's article (set.seed(20260909)): 60 species on a random coalescent tree, each observed at the genuinely irregular elapsed times 0, 0.5, 2.5 and 5 with an alternating treatment, 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. As above, we read the exact data and tree it produced.
col = read_fixture("vignette_phylo_ou_species.csv")
dat_phylo = (y = parse.(Float64, col("y")),
treatment = parse.(Float64, col("treatment")),
species = String.(col("species")),
elapsed_days = parse.(Float64, col("elapsed_days")))
tree = read(pkgdir(DRModels, "test", "fixtures", "temporal",
"vignette_phylo_ou_species.newick"), String)
phylo_fit = drm(bf(@formula(y ~ treatment + phylo(1 | species) +
temporal(1 | species, elapsed_days, ou)),
@formula(sigma ~ 1)),
Gaussian(); data = dat_phylo, tree = tree)
phylo_tp = temporal_parameters(phylo_fit)(label = "temporal(1 | species, time = elapsed_days, structure = \"ou\")", structure = :ou, sd = 0.5429867728838157, phi = nothing, decay = 0.44246723280642586, sd_iid = nothing, sd_phylo = 0.4610330204884002, sigma = 0.3259935649984235, cor = nothing, pac = nothing)The two random terms must use the same species column, and the tree is passed to drm (drmTMB writes it inside the term, phylo(1 | species, tree = tree)). The three variance components are reported separately: sd_phylo is drmTMB's sd_phylo_stable, sd its sd_temporal, and decay its decay_temporal. Like drmTMB, DRModels.jl uses the tree's tip correlation matrix, so rescaling every branch length leaves the fit unchanged.
phylo_ref = parity_cell("vignette-phylo-ou")
pt = phylo_ref["temporal"]
phylo_names = Dict(phylo_fit.coefnames)[:mu]
[(quantity = q, drmTMB = round(r; digits = 4), DRModels = round(j; digits = 4),
abs_diff = abs(j - r))
for (q, r, j) in [("logLik", phylo_ref["fit"]["loglik"], loglik(phylo_fit)),
[(k, v, coef(phylo_fit, :mu)[findfirst(==(replace(k, "mu_" => "")), phylo_names)])
for (k, v) in phylo_ref["coef"]]...,
("sd_phylo_stable", pt["sd_phylo"], phylo_tp.sd_phylo),
("sd_temporal", pt["sd"], phylo_tp.sd),
("decay_temporal", pt["decay"], phylo_tp.decay),
("sigma", pt["sigma"], phylo_tp.sigma)]]7-element Vector{@NamedTuple{quantity::String, drmTMB::Float64, DRModels::Float64, abs_diff::Float64}}:
(quantity = "logLik", drmTMB = -220.6328, DRModels = -220.6328, abs_diff = 5.4569682106375694e-12)
(quantity = "mu_(Intercept)", drmTMB = 0.7637, DRModels = 0.7637, abs_diff = 6.45039577307216e-13)
(quantity = "mu_treatment", drmTMB = 0.3871, DRModels = 0.3871, abs_diff = 2.0650148258027912e-13)
(quantity = "sd_phylo_stable", drmTMB = 0.461, DRModels = 0.461, abs_diff = 1.1747269823558781e-12)
(quantity = "sd_temporal", drmTMB = 0.543, DRModels = 0.543, abs_diff = 2.0872192862952943e-13)
(quantity = "decay_temporal", drmTMB = 0.4425, DRModels = 0.4425, abs_diff = 1.0098588631990424e-12)
(quantity = "sigma", drmTMB = 0.326, DRModels = 0.326, abs_diff = 3.405054016525355e-13)The two packages reach the same maximum-likelihood optimum. With a smaller panel (drmTMB's article mentions 8 or 16 species at these four times), the residual SD can run to zero while the OU term absorbs it; a fit at that boundary should not be read component by component.
drmTMB's fitted() adds both conditional components. In DRModels.jl they are in ranef: ranef(phylo_fit)[:species_phylo] holds the stable effect of each species (first-seen order) and ranef(phylo_fit)[:species] the OU departure of each row.
re = ranef(phylo_fit)
phylo_conditional = fitted(phylo_fit) .+ re[:species] .+
re[:species_phylo][indexin(dat_phylo.species, unique(dat_phylo.species))]
(max_abs_diff = maximum(abs.(phylo_conditional .- Float64.(phylo_ref["conditional"]["fitted"]))),)(max_abs_diff = 1.6553425297161084e-12,)residuals(phylo_fit; type = :quantile) subtracts both conditional components and divides by
simulate(phylo_fit) draws a new phylogenetic stable vector, a new independent OU path for each species and new noise.
The pairing is intentionally narrow, and both packages refuse the same departures from it: it needs ou (not ar1), an unlabelled phylo(1 | species) intercept on the same grouping as temporal(), no ordinary (1 | species) (the stable between-species component is already the phylogenetic term), at least three species with at least two distinct times each, at least three distinct positive lags across the data, tree tips that are exactly the observed species, and an ultrametric tree (all tips at the same depth, to drmTMB's relative sqrt(eps) tolerance). Zero-length branches are fine.
Every temporal fit is also checked against drmTMB's boundary rules: a random-effect SD (here sd_phylo or sd) below 1e-4, a residual SD below 1e-3 of the response SD, an OU decay that leaves correlation at about 1 or 0 at every observed lag, or AR1 persistence beyond ±0.999. Such a fit warns when it is fitted and reports check_drm(fit).temporal_boundary.at_boundary = true; do not read its variance components separately. As in drmTMB, the residual-SD rule compares with the marginal SD of the response, so it can also fire when a covariate explains most of that SD.
What can be reported now. drmTMB's article marks this model as a development workflow: its parser, dense-likelihood, method and fixed-mean profile checks pass, and a replacement point-recovery study met its predeclared point criteria, but no interval-calibration study has been run, so it asks readers not to report confidence intervals from the paired model yet. DRModels.jl reproduces the same likelihood and makes no stronger claim: use the fit to inspect the variance components, not for interval inference. For the record only, the fixed-mean profile for treatment agrees between the two packages. Profile confint on the paired fit warns that the interval is not calibrated, as drmTMB's does. drmTMB refuses Wald and bootstrap intervals for this fit. DRModels.jl returns them (confint(fit) and bootstrap_ci) as an extension, without a calibration claim, and does not warn that they are uncalibrated, so the same caution applies to them without a reminder:
compare_ci(confint(phylo_fit; method = :profile, parm = :mu => "treatment"), phylo_ref["profile"])1-element Vector{@NamedTuple{coef::String, lower::Float64, upper::Float64, drmTMB_lower::Float64, drmTMB_upper::Float64, max_abs_diff::Float64}}:
(coef = "treatment", lower = 0.27177867502731057, upper = 0.5020291188445459, drmTMB_lower = 0.27177986075433835, drmTMB_upper = 0.5020281537416368, max_abs_diff = 1.1857270277859655e-6)The drmTMB spelling through the bridge
An R user calling DRModels.jl through drm_bridge (see Coming from R?) can send drmTMB's own keyword spelling as a formula string. The bridge translates it and reaches the same optimum:
out = drm_bridge(; formula = "y ~ treatment + (1 | site) + " *
"temporal(1 | site, time = elapsed_days, structure = \"ou\"); sigma ~ 1",
family = "gaussian", data = dat_irregular)
(bridge_loglik = out["loglik"], native_loglik = loglik(ou_fit))(bridge_loglik = -102.58057815043881, native_loglik = -102.58057815043881)Evidence
Eight drmTMB parity cells — AR1 and OU, each with and without (1 | id), on these two datasets and on two further simulated fixtures — are checked on every test run (test/test_parity_temporal.jl, cells in test/parity/temporal/). In the measured run (Julia 1.12.6 on Linux, Totoro), the log-likelihood agreed with drmTMB to better than test/test_temporal_ar1.jl, test/test_temporal_ou.jl).
Three cells cover homogeneous Toeplitz (drmTMB's article data above, a 40-site six-occasion panel with a non-exponential lag pattern, and a 30-site four-occasion panel with a negative first-lag correlation), generated from the drmTMB development version: logLik within test/test_temporal_homtoep.jl).
Two more cells cover the paired phylogenetic + OU model: drmTMB's article data above and a 24-species simulated fixture. On Julia 1.10.12 and 1.13 (Linux, Totoro) they agree with drmTMB to about residuals(fit; type = :quantile) against drmTMB's Pearson residuals) agree to 012258e9f (recorded in each cell's expected.meta.toml). The paired likelihood is checked against a dense oracle (test/test_temporal_phylo_ou.jl.