Estimating epidemiological delay distributions

From R/Stan to Julia

Sam Abbott

London School of Hygiene & Tropical Medicine

14 August 2026

Talk plan

  • What an epidemiological delay is, and where the estimate ends up
  • Both events arrive as windows, and the pairs that have not finished are missing
  • primarycensored, and what shipping the same likelihood in Stan cost
  • Delays vary, so partial pooling, brms and epidist
  • No line list, so epinowcast
  • The same three adjustments in Julia, as composable wrappers

samabbott.co.uk/JuliaCon2026/delays

JuliaCon 2026, Muschel — N3, Friday 14 August 2026, 14:45.

Delays are not observed directly

A delay is the time between two events

  • Infection to symptom onset, onset to hospitalisation, onset to death
  • Nowcasts, forecasts and transmission models take a delay distribution as given
  • Neither event is seen at a point. Both arrive as a day, sometimes a week
  • In real time the longer delays have not happened yet

Diagram for this talk, scripts/delays-double-censoring.py

A timeline. A blue box marks the primary event window of width w_P, a red box marks the secondary event window of width w_S, an arrow between them is the delay T = S - P, and a dashed line further right marks the observation cutoff C, beyond which longer delays have not yet been seen.

Estimate one naively and it is biased

  • In a growing epidemic short delays are oversampled, so the delay comes out too short
  • Interval censoring of both events applies to almost any line list
  • Right truncation is the correction when the observation cutoff sits near the primary events
  • How cases were ascertained matters. On the primary event the delay is right censored instead

Fig 3 of Charniga et al. (2024), PLOS Comput Biol, CC-BY 4.0

A decision tree. Questions about time variation, whether the analysis is retrospective or real time, whether the cutoff was near the primary events, and which event cases were ascertained through, lead to actions: adjust for interval censoring and right censoring, adjust for interval censoring and right truncation, or interval censoring alone.

1. The primary event sits in a window

  • The primary event falls somewhere in a window \([0, w_P]\)
  • It need not be uniform, but a uniform primary is fair when the window is short
  • Mix the delay over where in the window the event sat, a convolution on the CDF

A timeline. The primary event sits at an unknown point tau inside a blue window of width w_P, the secondary event is a single red line, and the delay is f(t - tau) averaged over tau.

\[ F_{\mathrm{PEC}}(t) = \int_{0}^{w_P} g(\tau)\, F(t - \tau)\, d\tau \]

Diagram from How to serial interval?, labelled there for a serial interval · maths after primarycensored

Note

Gamma, lognormal and Weibull against a uniform primary have closed forms.

2. So does the secondary event

  • The secondary event is recorded to the day as well
  • The probability of a daily bin is a difference of the primary-censored CDF
  • It doubles as the probability mass function for daily discrete-time models

The primary-censored delay density, with the daily bin from n to n + w_S shaded red and labelled as the difference of the CDF at its two edges.

\[ \Pr\!\big(T \in [n, n + w_S]\big) = F_{\mathrm{PEC}}(n + w_S) - F_{\mathrm{PEC}}(n) \]

Diagram from How to serial interval?

Note

In the past this was often written as a double integral, which is typically more complicated to evaluate.

3. In real time you only see the pairs that finished

  • A pair enters the data only if the secondary event happened before the cutoff \(C\)
  • Condition on being observed, so divide by the CDF at the horizon
  • Short delays are over-represented early in an epidemic

Five delays from one primary event. The three that end before the cutoff C are drawn in red and observed, the two that end after it are dashed and never seen, and the likelihood is divided by the CDF at C minus P.

\[ L^{\mathrm{RT}} = \frac{F_{\mathrm{PEC}}(n + w_S) - F_{\mathrm{PEC}}(n)} {F_{\mathrm{PEC}}(C - P)} \]

Diagram from How to serial interval?

What the truncation costs, in days

  • The shortfall grows the nearer the cutoff is to now
  • A mean of 3.9 days rather than 5.9 is what the forecast and the transmission model are handed

LogNormal(1.6, 0.6) with daily windows, computed for this talk with CensoredDistributions.jl v0.2.22, scripts/delays-right-truncation.jl

Two panels. On the left, the delay distribution seen by day 7 and by day 14 sits well to the left of the true distribution. On the right, the mean estimated delay rises from under 2 days towards the true mean of 5.9 days as the cutoff moves away from now, passing 3.9 days at a cutoff of 7 and 5.3 days at a cutoff of 14.

Software for the three adjustments

I have a line list and I want a delay out of it

  • primarycensored wraps any delay family in the three adjustments
  • It picks the closed form where the pair has one and integrates numerically otherwise
  • It ships the same likelihood twice, in R and as Stan functions
delay_data <- data.frame(
  left = samples, right = samples + swindow,
  pwindow = rep(pwindow, n), D = rep(D, n)
)

fit_norm <- fitdistdoublecens(
  delay_data, distr = "norm",
  start = list(mean = 0, sd = 1)
)

R/fitdistdoublecens.R:93-103, the roxygen example, reflowed · primarycensored

Note

Stan has no import statement, so the 1,081 line function bundle is spliced into the downstream model as text.

The integral solver would not hold

// before
result = integrate_1d(
  primary_censored_integrand, lower_bound, d,
  theta, {d, pwindow}, ids, 1e-2
);

// after
result = ode_rk45(
  primary_censored_ode, y0, lower_bound, {d},
  theta, {d, pwindow}, ids
)[1, 1];
  • Stan’s quadrature errored inside sampling rather than rejecting, so the chain died instead of the proposal
  • 243 of 256 chains gone. Nine days of tuning got it to 12 per cent. Recast as an ODE, 0 of 600
  • The integral was fine. Stan’s autodiff and the missing reject were what made it unusable

epinowcast/primarycensored#34 and PR #64, commit 1c4e6f1, September 2024. The scenario was made harder in the same PR, so the 600 chain runs are the like for like pair

What Stan has no vocabulary for

real dist_lcdf(real delay, array[] real params,
               int dist_id) {
  if (dist_id == 1)
    return lognormal_lcdf(delay | params[1], params[2]);
  else if (dist_id == 2)
    return gamma_lcdf(delay | params[1], params[2]);
  // ... 15 more branches ...
  else if (dist_id == 25)
    return von_mises_lcdf(delay | params[1], params[2]);
  else reject("Invalid distribution identifier: ",
              dist_id);
}
  • No user types, so a distribution is an integer. 18 branches for the 25 families R knows, and the other seven have no Stan CDF to call
  • A second table of ten says which have positive support, kept in step by hand
  • dist_id = 3 meant Normal in Stan and Weibull in R. Ask for a Weibull on the numerical path and you fitted a normal

primarycensored_ode.stan:38-51 and :52-77 · epinowcast/primarycensored#277, commit 4aaf249, February 2026

Delays vary by district, by age and by wave

  • Fitting each stratum on its own runs out of data. Partial pooling shares what the small strata cannot support
  • brms metaprograms Stan and hands you the whole regression apparatus, so epidist puts the censored likelihood behind a brms formula

Sierra Leone Ebola line list, epidist/vignettes/ebola.Rmd:244-250, renamed and sampler arguments cut · epidist · brms

fit <- epidist(
  data = obs_prep,
  formula = bf(
    mu ~ 1 + sex + (1 | district),
    sigma ~ 1 + sex + (1 | district)
  ),
  family = lognormal()
)

Note

The delay family has to be one brms already knows.

Extending brms took more than a formula

  • family.R, formula.R, stancode.R and prior.R exist only to bend brms around a censored likelihood. 774 lines with marginal_model.R
  • Five S3 generics were defined to hang the extensions on, from epidist_family_model to epidist_family_prior
  • The Stan source is rewritten by string substitution before it compiles, five gsub() passes over a template

epidist/R/marginal_model.R:339-348, one of five passes at :332-368 · epidist

dist_id <- primarycensored::pcd_stan_dist_id(family_name)

# Replace the dist_id passed to primarycensored
stanvars_functions[[1]]$scode <- gsub(
  "dist_id",
  dist_id,
  stanvars_functions[[1]]$scode,
  fixed = TRUE
)

Does Julia need a brms

  • A Julia one would very likely be more extensible. The metaprogramming would not be fighting a second language with no user types
  • brms exists partly because in R you cannot drop the delay machinery into a model you wrote yourself
  • In Julia a censored distribution is an ordinary Distributions.jl object. It goes into any model somebody writes, with no formula interface and no code generation
  • A few people here are probably working on one

Note

Where is a formula interface still needed, once the components compose on their own?

Composing components rather than generating models is the 16:45 talk in this room, samabbott.co.uk/JuliaCon2026/composable

Sometimes there is no line list

  • Counts arrive by reference date and by report date, and today’s count is not finished. Right truncation on a table, not on a pair
  • epinowcast writes the delay as discrete-time reporting hazards, so effects and random walks can modify it by delay and by date
  • It estimates the reporting process and the epidemic curve together
  • It is Stan, as are primarycensored and epidist. One backend the whole way up

epinowcast

The count reported at delay \(d\)

\[ \mathbb{E}[n_{t,d}] = \lambda_t\, h_d \textstyle\prod_{d'<d}(1-h_{d'}) \]

Note

The README says the default lognormal reporting delay can fail on multimodal delays.

What about in Julia

using CensoredDistributions, Distributions, Turing

@model function double_censored_model(values, weights)
    α ~ truncated(Normal(1, 2), 0, Inf)
    θ ~ truncated(Normal(1, 2), 0, Inf)
    d = double_interval_censored(
        Gamma(α, θ); upper = 15, interval = 1
    )
    values ~ weight(d, weights)
end

model = double_censored_model(values, weights)
chain = sample(model, NUTS(), MCMCThreads(), 1000, 2)

Model code from README.md:71-85, run for this talk against v0.2.22 · CensoredDistributions.jl

  • Changing language here changes the backend. Turing, not Stan
  • Each adjustment is a wrapper. primary_censored is the \(w_P\) integral, upper truncates at \(C\), interval_censored is \(w_S\). double_interval_censored stacks all three in that order
  • The delay is built inside the model from the sampled parameters
  • cdf(primary_censored(Frechet(2, 1), Uniform(0, 1)), 3.0) returns 0.847. Nobody wrote that method

R wins this, by about a hundred to one

R Julia
Downloads, all time 11,244 116
Downloads, last month 528 49

CRAN and JuliaPkgStats, 10 August 2026. CRAN counts requests to one mirror, JuliaPkgStats counts unique IP addresses, so the ratio is an upper bound. Both include continuous integration, so neither is a count of people.

  • Nine repositories on GitHub depend on primarycensored, two of them on CRAN. Seven are from my own orgs. Two are CDC’s, both as suggests
  • Everything that depends on CensoredDistributions.jl is mine, Sam Brand’s or Sebastian Funk’s. I looked for one from outside and did not find it

Note

Twenty-one months on CRAN against twelve in General. I do not know how much of the gap is R having the users and how much is the package.

Where this goes

  • ConvolvedDistributions.jl adds sums and differences of delays
  • The generation time is the same shape of problem, a convolution and a difference
  • A composed Julia epidist, with Turing submodels doing what the brms formula does now

Important

The epidist layer on top is not written yet. Open issue CensoredDistributions.jl#749 says the joint fit is too slow.

A graph of a composed delay. A latent period and a transmission delay combine into a generation time, the generation time and a difference of incubation periods combine into a serial interval, and double interval censoring turns that into the observed delay.

using CensoredDistributions
using ConvolvedDistributions, Distributions

incubation = Gamma(2.0, 1.0)
reporting = LogNormal(1.0, 0.5)

d = convolved(incubation, reporting)
observed = double_interval_censored(
    d; upper = 15, interval = 1
)

cdf(observed, 5.0)   # 0.4561

Diagram from How to serial interval?, labelled there for a serial interval · Plain Distributions.jl dispatch, composed for this talk rather than lifted from a README. Run against CensoredDistributions.jl v0.2.22 and ConvolvedDistributions.jl. Neither package appears in the other’s [deps], and both list Distributions

Thank you

Important

Tell me about a delay you cannot write down as a composition.

samabbott.co.uk/JuliaCon2026/delays

Sam Abbott, LSHTM · epiaware.org · github.com/seabbs

Backup

Three closed forms, picked by type

primarycensored_cdf(dist::Gamma, primary_event::Uniform,
    x::Real, ::AnalyticalSolver)

primarycensored_cdf(dist::LogNormal, primary_event::Uniform,
    x::Real, ::AnalyticalSolver)

primarycensored_cdf(dist::Weibull, primary_event::Uniform,
    x::Real, ::AnalyticalSolver)

primarycensored_cdf(dist::D1, primary_event::D2,
    x::Real, method::NumericSolver
) where {D1 <: UnivariateDistribution,
         D2 <: UnivariateDistribution}
  • Three closed forms picked by the type you passed, and one fallback for any UnivariateDistribution
  • Stan has the same three, and in Stan Frechet needs a new branch, a new integer and a rebuild

primarycensored_cdf.jl:184-201, :323, :373, :415, v0.2.22 · primarycensored_analytical_cdf.stan

Everything I tried in between

Horizontal bar chart of ten runs reported in epinowcast/primarycensored issue 34. The first, integrate_1d at tolerance 1e-6, lost 243 of 256 chains. Eight tuning attempts sit between 12 and 27 per cent, one of them a model with the truncation adjustment removed. The last, recast as ode_rk45, lost 0 of 600.

Every run is one Sam Abbott reported in epinowcast/primarycensored#34 between 5 and 12 September 2024, read back from the issue, scripts/delays-stan-failures.py. The runs are not one dataset. The scenario was made harder in the same PR that landed the ODE, so the last two rows on 600 chains are the pair closest to like for like.

The observation process is the type

julia> typeof(double_interval_censored(
           Gamma(2.0, 3.0); upper = 15, interval = 1))

IntervalCensored{
  Truncated{
    PrimaryCensored{
      Gamma{Float64},
      Uniform{Float64},
      AnalyticalSolver{GaussLegendre{…}}}}}
  • The tree is the type. Nothing dispatches on a name or an integer, so there is nothing to keep in step
  • Swap Gamma for any UnivariateDistribution and the rest of the tree is unchanged

Real output, CensoredDistributions.jl v0.2.22. Reflowed onto several lines, and the type parameters of the solver and of Truncated are elided at …

The numerical problem moved, it did not go away

Autodiff backend Broken
ForwardDiff 0
ReverseDiff, tape 0
Mooncake forward 3
Mooncake reverse 3
Enzyme forward 4
Enzyme reverse 8

Scenarios registered broken out of 63, in test/ADFixtures/src/ADFixtures.jl

  • Six autodiff backends, each with its own CI workflow. Two of them pass every scenario, four carry scenarios registered broken
  • The default solver is GaussLegendre(; n = 64) because fixed nodes trace through all six backends and adaptive ones do not
  • Stan gave up adaptive quadrature for an ODE. Julia gave it up for a dot product

Six per-backend workflows in .github/workflows/, over a shared ad-backend.yaml · why an ecosystem needs that CI was the 14:30 talk in this room, samabbott.co.uk/JuliaCon2026/roadmap

The same component, in two languages

primarycensored, R and Stan CensoredDistributions.jl
Adding a distribution New branch, new integer, rebuild Pass the type
Distributions covered 18, listed in a switch Any univariate one
Choosing one An integer, set before sampling The object itself
Sharing the code Splice 1,081 Stan lines in at build Project.toml
The integral integrate_1d, then an ODE recast 64 fixed nodes
Worst failure 243 of 256 chains died 8 scenarios broken on one backend

Important

The Stan file has to know every distribution it will ever support. The Julia file does not have to know any of them.

primarycensored_ode.stan:52-77 and pcd_functions.stan · primarycensored_cdf.jl:184-201, v0.2.22 · epinowcast/primarycensored#34 and PR #64, September 2024

What the move fixed, and what it did not

Removed

0

tables left to keep in step

A class of bug. dist_id = 3 meant Normal in Stan and Weibull in R, so asking for a Weibull fitted a normal and returned numbers. The two integer tables agreed on lognormal, gamma and exponential, and on nothing else.

primary_censored(Weibull(2, 1), Uniform(0, 1)) carries the type. There is no second table for it to disagree with.

Moved

4 of 6

autodiff backends carry a broken scenario

The numerical problem, from the language to the autodiff backends. Six backends run in CI, each with its own workflow. Two pass every scenario, four carry scenarios registered broken, out of 63.

The default solver is GaussLegendre(; n = 64) because fixed nodes trace through every backend. Stan gave up adaptive quadrature for an ODE, Julia for a dot product.

Unchanged

100×

more downloads in R, all time

Who uses it. 11,244 downloads all time against 116, and 528 last month against 49.

Nine GitHub repositories depend on primarycensored, two of them on CRAN. Seven are mine, two are CDC’s. Everything depending on CensoredDistributions.jl is mine, Sam Brand’s or Sebastian Funk’s.

Twenty-one months on CRAN against twelve in General, so the gap is not all the package.

CRAN and JuliaPkgStats, 10 August 2026. CRAN counts requests to one mirror, JuliaPkgStats counts unique IP addresses, so the ratio is an upper bound · epidist/R/marginal_model.R:339-347