---
title: "Left truncation with delay_min"
description: |
  "How to use the delay_min argument to account for delays that are only
  observed above a minimum value."
output:
  bookdown::html_vignette2:
    fig_caption: yes
    code_folding: show
    number_sections: true
pkgdown:
  as_is: true
link-citations: true
vignette: >
  %\VignetteIndexEntry{Left truncation with delay_min}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: references.bib
---



Some delays are only ever recorded above a minimum value.
Generation intervals are often defined to exclude same day transmission.
Case definitions sometimes require a minimum time between two events before a pair is recorded.
In both settings the observed delays are a left truncated sample from the delay distribution of interest.

Write $L$ for that minimum and $F$ for the delay distribution.
Only delays above $L$ can enter the sample, so the observed delays follow the conditional distribution
$$
F(t \mid T > L) = \frac{F(t) - F(L)}{1 - F(L)}, \qquad t > L .
$$
Fitting $F$ directly to such a sample biases the estimate.
The model has to explain the missing short delays by shifting the whole distribution.

With right truncation at $D$ as well, the likelihood renormalises over $[L, D]$
$$
\frac{F(t) - F(L)}{F(D) - F(L)} .
$$

The `delay_min` argument of `as_epidist_marginal_model()` sets $L$.
It is passed to the [`primarycensored`](https://primarycensored.epinowcast.org/) likelihood as its `L` argument.
`delay_min` is a duration in the same units as the delay, not a date.

`delay_min` truncates on the same scale as the response, that is the observed delay between the primary window lower bound and the secondary event.
Data should therefore be excluded on the observed delay rather than on the underlying continuous delay.

# Setup


``` r
library(epidist)
library(ggplot2)
library(dplyr)
```

# Simulate left truncated data

We simulate an outbreak and a lognormal delay from each primary event.
The parameters are chosen so that a large share of the distribution falls below the truncation point.


``` r
set.seed(101)

true_meanlog <- 1.0
true_sdlog <- 0.8
delay_min <- 2

obs <- simulate_gillespie(r = 0.2, seed = 101) |>
  simulate_secondary(
    dist = rlnorm,
    meanlog = true_meanlog,
    sdlog = true_sdlog
  )
```

Roughly this share of the delay distribution lies below the truncation point.


``` r
plnorm(delay_min, true_meanlog, true_sdlog)
#> [1] 0.3506501
```

`simulate_dates()` adds daily reporting windows and returns date.
We then drop the pairs whose observed delay is below `delay_min`.


``` r
obs_trunc <- obs |>
  simulate_dates() |>
  filter(as.numeric(sdate_lwr - pdate_lwr) >= delay_min) |>
  slice_sample(n = 500)

linelist <- as_epidist_linelist_data(obs_trunc)
```

The observed delays are truncated below `delay_min`.


``` r
plot_delays(linelist, delay_min = delay_min)
```

<div class="figure">
<img src="figures/left-truncation-hist-1.png" alt="Observed delays after left truncation. The dashed line marks the minimum delay."  />
<p class="caption">(\#fig:hist)Observed delays after left truncation. The dashed line marks the minimum delay.</p>
</div>

# Fit models with and without the adjustment


``` r
marginal_no_trunc <- as_epidist_marginal_model(linelist)
marginal_trunc <- as_epidist_marginal_model(linelist, delay_min = delay_min)
```

A `delay_min` column already in the data is picked up without the argument.
Pass a column name instead of a number when the minimum varies between observations.


``` r
linelist$delay_min <- delay_min

identical(
  as_epidist_marginal_model(linelist)$delay_min,
  marginal_trunc$delay_min
)
#> [1] TRUE
```

`delay_min` is stored as a column and passed to the likelihood.


``` r
marginal_trunc
#> # A tibble: 500 × 19
#>    ptime_lwr ptime_upr stime_lwr stime_upr obs_time  case pdate_lwr  pdate_upr
#>        <dbl>     <dbl>     <dbl>     <dbl>    <dbl> <int> <date>     <date>
#>  1        33        34        38        39       84  7022 2024-02-03 2024-02-04
#>  2        40        41        53        54       84  8003 2024-02-10 2024-02-11
#>  3        13        14        17        18       84  1000 2024-01-14 2024-01-15
#>  4        20        21        28        29       84  2739 2024-01-21 2024-01-22
#>  5        18        19        22        23       84  2301 2024-01-19 2024-01-20
#>  6        26        27        30        31       84  5070 2024-01-27 2024-01-28
#>  7        22        23        26        27       84  3556 2024-01-23 2024-01-24
#>  8        32        33        35        36       84  6810 2024-02-02 2024-02-03
#>  9        11        12        15        16       84   720 2024-01-12 2024-01-13
#> 10        29        30        31        32       84  6148 2024-01-30 2024-01-31
#> # ℹ 490 more rows
#> # ℹ 11 more variables: sdate_lwr <date>, sdate_upr <date>, obs_date <date>,
#> #   pwindow <dbl>, swindow <dbl>, relative_obs_time <dbl>,
#> #   orig_relative_obs_time <dbl>, delay_lwr <dbl>, delay_upr <dbl>, n <dbl>,
#> #   delay_min <dbl>
```


``` r
fit_no_trunc <- epidist(
  marginal_no_trunc,
  chains = 2, cores = 2, refresh = ifelse(interactive(), 250, 0)
)

fit_trunc <- epidist(
  marginal_trunc,
  chains = 2, cores = 2, refresh = ifelse(interactive(), 250, 0)
)
```

# Compare parameter estimates

`add_delay_parameter_draws()` returns posterior draws of the distributional parameters.
For the lognormal family `mu` is the log scale mean and `sigma` is the log scale standard deviation.
We use `epidist_strata()` to take a single row of the transformed data because there are no covariates here.
`add_summaries()` adds the natural scale mean and standard deviation, and records the family so that the draws can be plotted.


``` r
param_draws <- list(
  "No adjustment" = fit_no_trunc,
  "With delay_min" = fit_trunc
) |>
  lapply(\(fit) delay_parameter_draws(fit, newdata = epidist_strata(fit))) |>
  bind_rows(.id = "model") |>
  add_summaries(family = fit_trunc)

true_values <- data.frame(mu = true_meanlog, sigma = true_sdlog) |>
  add_summaries(family = "lognormal") |>
  unlist()
```


``` r
plot(param_draws, by = "model", true_values = true_values)
```

<div class="figure">
<img src="figures/left-truncation-params-plot-1.png" alt="Posterior draws of the lognormal parameters, and of the natural scale mean and standard deviation. Dashed lines are the simulation values. The unadjusted model is biased away from them."  />
<p class="caption">(\#fig:params-plot)Posterior draws of the lognormal parameters, and of the natural scale mean and standard deviation. Dashed lines are the simulation values. The unadjusted model is biased away from them.</p>
</div>

The unadjusted model shifts the distribution to the right and understates its spread.
The adjusted model recovers the simulation values.

# Compare the estimated delay distributions

Plotting the same draws with `type = "delay"` draws the delay distribution each model implies, as the posterior median density with a ribbon between the 5% and 95% quantiles.


``` r
plot(param_draws, type = "delay", by = "model") +
  geom_function(
    fun = dlnorm,
    args = list(meanlog = true_meanlog, sdlog = true_sdlog),
    linewidth = 1
  )
```

<div class="figure">
<img src="figures/left-truncation-pdf-plot-1.png" alt="Estimated delay distributions against the true lognormal density (black line)."  />
<p class="caption">(\#fig:pdf-plot)Estimated delay distributions against the true lognormal density (black line).</p>
</div>

These are the underlying delay distributions rather than the truncated ones, because the distributional parameters do not describe the observation process.
Predictions from `tidybayes::add_predicted_draws()` do the same only when `newdata` sets `delay_min` to 0.

# Summary

Set `delay_min` whenever delays below a threshold could not have been observed.
The default of 0 leaves the likelihood unchanged.
`delay_min` must not exceed the smallest observed delay.
Above it the data are impossible under the model and `as_epidist_marginal_model()` errors.