Handling delays better?

[zsusswein]

@seabbs and I have been thinking about ways to better handle the delay from infection to case observation. We have a current thing that's only meh and two additional things that probably don't work and we're stuck on.

What we currently do

We have a fitted model that takes a user-requested vector of dates $\hat f(\tilde t)$. We ask the user for the dates they want to see $R_t$.

The model expects to be provided dates that correspond to cases, but the user-specified dates corresponds to an infection based measure. So we map the user-requested infection-corresponding dates to the model-expected case-corresponding dates by adding a user-specified mean delay $\delta$.

Then we simulate draws of the expected posterior incident cases out of $\hat f(\tilde t + \delta)$ and convolve them with the user-supplied GI PMF to get $R_t$.

This has the obvious problems of not correcting for the uncertainty induced by the delay convolution. It's been a hard problem to solve because we don't have a latent process; we're just fitting to cases.

Option 1: signal regression

mgcv has a "linear functional terms" feature that enables signal regression approaches of the form

$$g(\mu_i) = ... + \sum_j L_{ij}f(x_ij) + ...$$

where the weighted sum of $f(x)$ and $L$ is part of the model.

Hallelujah! Except the convolution happens on the scale of the linear predictor, not the mean response scale. And we're stuck.

One bad workaround would be to switch to an identity link, but that feels like an unpleasant road to walk down. Out of the "poor delay handling" frying pan and into the "negative cases" fire. This also would have a static variance over time, rather than having the variance be a function of the mean.

Option 2: Variance propagation a la density surface models

Inspired by the work here, we could implement their approach from section 3. In that work they have spatial data and some kind of "detection function" indicating likelihood of observation.

Which they smuggle in through a random effect via the Hessian in a way I don't get

(from these slides)

As far as I can tell, this would involve trying to do something clever with a Taylor series and the gamma distribution of the delay? I'm a little worried (1) we would lose something in the tail with the linearization of the convolution from the Taylor series approximation and (2) I'm not totally understanding the implementation. I think this is where the magic happens but it's a little scary.

Option 3: Put a tensor product smooth on the lag

Skip knowing the delay and instead throw more splines at the wall until you get something out. This feels dangerous and like a bad idea in our situation but see here because it's cool: https://ecogambler.netlify.app/blog/distributed-lags-mgcv/

More spaghetti to throw at the wall

Write a custom mgcv family that implements the convolution. Looks like we'd need a bunch of higher order derivatives of a convolution.

[SamuelBrand1]

Re Option 2 (only one I've really looked at).

The linked paper doesn't really work (at a couple of reads) because the data in the paper does have a distance to detection (when observed). E.g. the canonical modelling problem is moving along a track and then getting distance to target data observed from the track, then back-inferring the likely spatial density given an estimate of how likely you are to spot a target at a range.

I think the analogy to disease modelling a latent time series of infections would be if we had case data coming in at some cadence, and when we get it we also know the infection time of the cases but the chance of being an asserted case depends on time since infection. Thats not a typical situation.

Laplace Method

A vaguely analogous approx to the paper would be to try and approximate the integral

$$\mathbb{E}[n(t) | \cdot ] = \int_0^\infty D(s) \exp\left(\beta_0 + \sum_k f_k(t - s)\right) ds.$$

using Laplace method. Here I've mainly used the linked paper notation: $n(t)$ is observed counts, $D$ is a delay kernel distribution, the rest is a standard GAM formulation (detection probability it absorbed into offset).

Approx 1

Assuming that we can rewrite $D(s) = \exp\left(\tau g(s)\right)$ where $\tau$ is a precision parameter for the delay kernel then at high precision you basically get back the mean shift

$$\mathbb{E}[n(t) | \cdot ] \approx \exp\left(\beta_0 + \sum_k f_k(t - \mu)\right).$$

where $\mu$ is the argmax of $g(s)$. When this is a reasonable approx this will be the mean of the distribution. The reason is that at moderately high precision the delay kernel will be approx Gaussian and the Laplace Method just returns the above. There might be some non-Gaussian correction terms of interest though.

Approx 2

We try and exploit the GAM structure to do a joint optimisation. This time we write (nb that Laplace method approx does depend on how you choose your expansion):

$$\begin{aligned} \mathbb{E}[n(t) | \cdot ] &= \int_0^\infty \exp\left( \beta_0 + h(s) \right) ds \\\ h(s) &= \tau \left[g(s) + \frac{1}{\tau} \sum_k f_k(t-s) \right]. \end{aligned}$$

Again, at high precision this limits to mean shift however at moderate precision you get a shift correction. Define $\mu$ as the maximising solution to

$$g'(\mu) = \frac{1}{\tau} \sum_k f'_k(t-\mu)$$

then the Laplace approximation is

$$\mathbb{E}[n(t) | \cdot ] \approx \sqrt{2\pi} \exp\left(\beta_0 + \sum_k f_k(t-\mu) + \tau g(\mu) - \frac{1}{2} [\ln \Big| \tau g''(\mu) + \sum_k f''_k(t-\mu)\Big| ]\right).$$

In this version of the approximation the mean shift approach is corrected a bit by terms that involve:

• The curvature of the GAM smooths.
• The delay kernel parameters (though the function $g$).

For example, note that if we have a growing incidence rate (in the GAM representation this is $\sum_k f'_k(t-\mu) > 0$ then the "optimal" look back $\mu$ will occur when $g'(\mu) > 0$. For a unimodal delay distribution this will occur before the peak of the delay distribution e.g. this approx seems to get the phase effect on the backwards delays.

EDIT: PS this will already exist somewhere in the literature so I'll keep my eyes out.

[bbolker]

See #116 for one approach.

A completely different approach is to try to deconvolve the case reporting time series (either up-front with raw data, or based on the smooth estimate?)

Goldstein, Edward, Jonathan Dushoff, Junling Ma, Joshua B. Plotkin, David J. D. Earn, and Marc Lipsitch. 2009. “Reconstructing Influenza Incidence by Deconvolution of Daily Mortality Time Series.” Proceedings of the National Academy of Sciences 106 (51): 21825–29. https://doi.org/10.1073/pnas.0902958106.