In this vignette, we’ll derive the analytic solutions for the primary censored delay distributions for a range of commonly used distributions.
In epidemiological analysis, it is common for primary events to be arriving at exponentially increasing or decreasing rates, for example, incidence of new infections in a growing epidemic. In this case, the distribution of the primary event time within its censoring window is biased by the exponential growth or decay. If we assume a reference uniform distribution within a primary censoring window \([k, k + w_P)\) then the distribution of the primary event time within the censoring window is the exponential tilted uniform distribution:
\[ f_P(t) \propto \exp(r t) \mathbb{1}_{[k, k + w_P]}(t). \]
In this case, the distribution function for \(C_P\), that is the length of time left in the primary censor window after the primary event time, is given by:
\[ F_{C_P}(p; r) = { 1 - \exp(-r p) \over 1 - \exp(-r w_P)}, \qquad p \in [0, w_P]. \tag{2.1} \] Note that taking the limit \(\lim_r \rightarrow 0\) in equation (2.1) gives the uniform distribution function \(F_{C_P}(p, 0) = p / w_P\).
In the following, it is convenient to use a (linear) difference operator defined as:
\[ \Delta_{w}f(t) = f(t + w) - f(t). \]
Applying a uniform primary event time distribution to equation 3.1 from “Why it works” gives:
\[ Q_{S_+}(t) = Q_T(t + w_P) + { 1 \over w_P} \int_0^{w_P} f_T(t+p) p~ dp. \]
This is analytically solvable whenever the upper distribution function of \(T\) is known and the mean of \(T\) is analytically solvable from its integral definition.
In each case considered below it is easier to change the integration variable:
\[ \begin{aligned} Q_{S_+}(t) &= Q_T(t + w_P) + { 1 \over w_P} \int_t^{t+w_P} f_T(z) (z-t)~ dz \\ &= Q_T(t + w_P) + { 1 \over w_P} \Big[ \int_t^{t+w_P} f_T(z) z~ dz - t \Delta_{w_P}F_T(t) \Big]. \end{aligned} \tag{3.1} \]
Note that for any distribution with an analytically available distribution function \(F_T\) equation (3.1) can be solved so long as the partial expectation
\[ \int_t^{t+w_P} f_T(z) z~ dz \tag{3.2} \]
can be reduced to an analytic expression.
The insight here is that this will be possible for any distribution where the average of the distribution can be calculated analytically, which includes commonly used non-negative distributions such as the Gamma, Log-Normal and Weibull distributions.
First, we note that equation 3.3 from “Why it works” can be written using the difference operator: \(f_n = -\Delta_1 Q_{S_+}(n-1)\). We can insert this expression into equation (3.1) to give the discrete censored delay distribution for uniformly distributed primary event times:
\[ \begin{aligned} f_n &= \Delta_1\Big[(n-1) \Delta_1F_T(n-1)\Big] - \Delta_1Q_T(n) - \Delta_1\Big[ \int_{n-1}^n f_T(z) z ~dz \Big] \\ &= (n+1)F_T(n+1) + (n-1)F_T(n-1) - 2nF_T(n) - \Delta_1\Big[ \int_{n-1}^n f_T(z) z ~dz \Big]. \end{aligned} \tag{3.3} \]
We now consider some specific delay distributions.
The Gamma distribution has the density function:
\[ f_T(z;k, \theta) = {1 \over \Gamma(k) \theta^k} z^{k-1} \exp(-z/\theta). \] Where \(\Gamma\) is the Gamma function.
The Gamma distribution has the distribution function:
\[ \begin{aligned} F_T(z;k, \theta) &= {\gamma(k, z/\theta) \over \Gamma(k)}, \qquad z\geq 0,\\ F_T(z;k, \theta) &= 0, \qquad z < 0. \end{aligned} \] Where \(\gamma\) is the lower incomplete gamma function.
We know that the full expectation of the Gamma distribution is \(\mathbb{E}[T] = k\theta\), which can be calculated as a standard integral. Doing the same integral for the partial expectation gives:
\[ \begin{aligned} \int_t^{t+w_P} f_T(z) z~ dz &= {1 \over \Gamma(k) \theta^k} \int_t^{t+w_P} \mathbb{1}(z \geq 0) z z^{k-1} \exp(-z/\theta)~dz \\ &= {\Gamma(k+1) \theta^{k+1} \over \Gamma(k) \theta^k} {1 \over \Gamma(k+1) \theta^{k+1}} \int_t^{t+w_P} \mathbb{1}(z \geq 0) z^{k} \exp(-z/\theta)~dz\\ &= k\theta \Delta_{w_P} F_T(t; k + 1, \theta). \end{aligned} \tag{3.4} \]
By substituting equation (3.4) into equation (3.3) we can solve for the survival function of \(S_+\) in terms of analytically available functions:
\[ Q_{S_+}(t; k, \theta) = Q_T(t + w_P; k, \theta) + { 1 \over w_P} \big[ k \theta \Delta_{w_P}F_T(t; k+1, \theta) - t \Delta_{w_P}F_T(t; k, \theta) \big]. \tag{3.5} \]
By substituting (3.5) into (3.3) we get the discrete censored delay distribution in terms of analytically available functions: \[ \begin{aligned} f_n &= (n+1) F_T(n+1; k, \theta) + (n-1) F_T(n-1; k, \theta) - 2 n F_T(n; k, \theta) - k \theta \Delta_1^{(2)}F_T(n-1; k+1, \theta)\\ &= (n+1) F_T(n+1; k, \theta) + (n-1) F_T(n-1; k, \theta) - 2 n F_T(n; k, \theta) \\ &+ k \theta \Big( 2 F_T(n; k+1, \theta) - F_T(n-1; k+1, \theta) - F_T(n+1; k+1,\theta) \Big) \qquad n = 0, 1, \dots. \end{aligned} \]
Which was also found by Cori et al.[1]
The Log-Normal distribution has the density function:
\[ \begin{aligned} f_T(z;\mu, \sigma) &= {1 \over z \sigma \sqrt{2\pi}} \exp\left( - {(\log(z) - \mu)^2 \over 2 \sigma^2} \right)\\ F_T(z;\mu, \sigma) &= 0, \qquad z < 0. \end{aligned} \] And distribution function:
\[ F_T(z;\mu, \sigma) = \Phi\left( {\log(z) - \mu \over \sigma} \right). \] Where \(\Phi\) is the standard normal distribution function.
We know that the full expectation of the Log-Normal distribution is \(\mathbb{E}[T] = e^{\mu + \frac{1}{2} \sigma^2}\), which can be calculated by integration with the integration substitution \(y = (\ln z - \mu) / \sigma\). This has transformation Jacobian:
\[ \frac{dz}{dy} = \sigma e^{\sigma y + \mu}. \]
Doing the same integral for the partial expectation, and using the same integration substitution gives:
\[ \begin{aligned} \int_t^{t+w_P} z~ f_T(z; \mu, \sigma) dz &= {1 \over \sigma \sqrt{2\pi}} \int_t^{t+w_P} \mathbb{1}(z \geq 0) \exp\left( - {(\log(z) - \mu)^2 \over 2 \sigma^2} \right) dz \\ &= {1 \over \sqrt{2\pi}} \int_{(\ln t - \mu)/\sigma}^{(\ln(t+w_P) - \mu)/\sigma} e^{\sigma y + \mu} e^{-y^2/2} dy\\ &= {e^{\mu + \frac{1}{2} \sigma^2} \over \sqrt{2 \pi} } \int_{(\ln t - \mu)/\sigma}^{(\ln(t+w_P) - \mu)/\sigma} e^{-(y- \sigma)^2/2} dy \\ &= e^{\mu + \frac{1}{2} \sigma^2} \Big[\Phi\Big({\ln(t+w_P) - \mu \over \sigma} - \sigma\Big) - \Phi\Big({\ln(t) - \mu \over \sigma} - \sigma\Big) \Big]\\ &= e^{\mu + \frac{1}{2} \sigma^2} \Delta_{w_P}F_T(t; \mu + \sigma^2, \sigma). \end{aligned} \tag{3.6} \]
By substituting equation (3.6) into equation (3.3) we can solve for the survival function of \(S_+\) in terms of analytically available functions:
\[ Q_{S+}(t ;\mu, \sigma) = Q_T(t + w_P;\mu, \sigma) + { 1 \over w_P} \Big[ e^{\mu + \frac{1}{2} \sigma^2} \Delta_{w_P}F_T(t; \mu + \sigma^2, \sigma) - t\Delta_{w_P}F_T(t; \mu, \sigma) \Big] \]
By substituting (3.6) into (3.3) we get the discrete censored delay distribution in terms of analytically available functions:
\[ \begin{aligned} f_n &= (n+1) F_T(n+1; \mu, \sigma) + (n-1) F_T(n-1; \mu, \sigma) - 2 n F_T(n; \mu, \sigma) \\ &- e^{\mu + \frac{1}{2} \sigma^2} \Delta_1^{(2)}F_T(n-1;\mu + \sigma^2, \sigma) \\ &= (n+1) F_T(n+1; \mu, \sigma) + (n-1) F_T(n-1; \mu, \sigma) - 2 n F_T(n; \mu, \sigma) \\ &+ e^{\mu + \frac{1}{2} \sigma^2} \Big( 2 F_T(n; \mu + \sigma^2, \sigma) - F_T(n+1; \mu + \sigma^2, \sigma) - F_T(n-1; \mu + \sigma^2, \sigma) \Big)\qquad n = 0, 1, \dots. \end{aligned} \]
The Weibull distribution has the density function:
\[ f_T(z;\lambda,k) = \begin{cases} \frac{k}{\lambda}\left(\frac{z}{\lambda}\right)^{k-1}e^{-(z/\lambda)^{k}}, & z\geq0 ,\\ 0, & z<0, \end{cases} \] And distribution function:
\[ F_T(z;\lambda,k))=\begin{cases}1 - e^{-(z/\lambda)^k}, & z\geq0,\\ 0, & z<0.\end{cases} \] Where \(\Phi\) is the standard normal distribution function.
We know that the full expectation of the Weibull distribution is \(\mathbb{E}[T] = \lambda \Gamma(1 + 1/k)\), which can be calculated by integration using the integration substitution \(y = (z / \lambda)^k\), which has transformation Jacobian:
\[ \frac{dz}{dy} = \frac{\lambda}{k}y^{1/k - 1}. \]
Doing the same integral for the partial expectation, and using the same integration substitution gives:
\[ \begin{aligned} \int_{t}^{t+w_P} z~ f_T(z; \lambda,k) dz &= \int_t^{t+w_P} \mathbb{1}(z \geq 0) \frac{kz}{\lambda}\left(\frac{z}{\lambda}\right)^{k-1}e^{-(z/\lambda)^{k}} dz \\ &= k\int_t^{t+w_P} \mathbb{1}(z \geq 0) \left(\frac{z}{\lambda}\right)^{k}e^{-(z/\lambda)^{k}} dz \\ &= \lambda k \int_{(t / \lambda)^k}^{((t + w_P) / \lambda)^k} \mathbb{1}(y \geq 0) y y^{1/k - 1} e^{-y} dy \\ &= \lambda\int_{(t / \lambda)^k}^{((t + w_P) / \lambda)^k} \mathbb{1}(y \geq 0) y^{1/k} e^{-y} dy\\ &= \lambda \Delta_{w_P} g(t; \lambda,k) \end{aligned} \tag{3.7} \]
Where
\[ g(t; \lambda, k) = \gamma\left(1 + 1/k, \left({t\vee 0 \over \lambda}\right)^k\right) = \frac{1}{k}\gamma\left(1/k, \left({t\vee 0 \over \lambda}\right)^k\right) - \frac{t}{\lambda}\exp\left(-\left({t\vee 0 \over \lambda}\right)^k\right) \] is a reparametrisation of the lower incomplete gamma function. Note that the \(\vee\) operator \(t \vee 0 = \text{max}(0, t)\) comes into the expression due to \(\mathbb{1}(y \geq 0)\) term in the integrand.
By substituting equation (3.7) into equation (3.3) we can solve for the survival function of \(S_+\) in terms of analytically available functions:
\[ Q_{S+}(t ;\lambda,k) = Q_T(t + w_P; \lambda,k) + { 1 \over w_P} \Big[ \lambda \Delta_{w_P} g(t; \lambda,k) - t\Delta_{w_P}F_T(t; \lambda,k)\Big]. \]
By substituting (3.7) into (3.3) we get the discrete censored delay distribution in terms of analytically available functions:
\[ \begin{aligned} f_n &= (n+1)F_T(n+1) + (n-1)F_T(n-1) - 2nF_T(n) - \Delta_1\Big[ \int_{n-1}^n f_T(z) z ~dz \Big] \\ &= (n+1)F_T(n+1) + (n-1)F_T(n-1) - 2nF_T(n) - \lambda \Delta_1^{(2)} g(n-1; \lambda,k) \\ &= (n+1)F_T(n+1) + (n-1)F_T(n-1) - 2nF_T(n) \\ &+ \lambda [2 g(n; \lambda,k) - g(n+1; \lambda,k) - g(n-1; \lambda,k)] \qquad n = 0, 1, \dots. \end{aligned} \]
Which was also found by Cori et al.[1]
vignette("why-it-works").vignette("primarycensored").1. Cori, A., Ferguson, N. M., Fraser, C., & Cauchemez, S. (2013). A new framework and software to estimate time-varying reproduction numbers during epidemics. American Journal of Epidemiology, 178(9), 1505–1512.