Analytic solutions for censored delay distributions
Source:vignettes/analytic-solutions.Rmd
analytic-solutions.Rmd1 Analytic solutions for exponentially tilted primary event times
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{1.1} \] Note that taking the limit \(\lim_r \rightarrow 0\) in equation (1.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). \]
2 Uniform primary event time (\(r=0\))
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{2.1} \]
General partial expectation
Note that for any distribution with an analytically available distribution function \(F_T\) equation (2.1) can be solved so long as the partial expectation
\[ \int_t^{t+w_P} f_T(z) z~ dz \tag{2.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.
General Discrete censored delay distribution
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 (2.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{2.3} \]
We now consider some specific delay distributions.
2.1 Gamma distributed delay times
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.
Gamma partial expectation
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}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}z^{k} \exp(-z/\theta)~dz\\ &= k\theta \Delta_{w_P} F_T(t; k + 1, \theta). \end{aligned} \tag{2.4} \]
Survival function of \(S_{+}\) for Gamma distribution
By substituting equation (2.4) into equation (2.3) we can solve for the survival function of \(S_+\) in terms of analytically available functions:
\[ \begin{aligned} Q_{S_+}(t; k, \theta) &= Q_T(t + w_P; k, \theta) + { 1 \over w_P} \left( {1 \over \Gamma(k) \theta^k} \int_t^{t+w_P} z^{k-1} \exp(-z/\theta) (z-t)~ dz \right), \\ &= 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]. \end{aligned} \tag{2.5} \]
Gamma discrete censored delay distribution
By substituting (2.5) into (2.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].
2.2 Log-Normal distribution
The Log-Normal distribution has the density function:
\[ f_T(z;\mu, \sigma) = {1 \over z \sigma \sqrt{2\pi}} \exp\left( - {(\log(z) - \mu)^2 \over 2 \sigma^2} \right). \] And distribution function:
\[ F_T(z;\mu, \sigma) = \Phi\left( {\log(z) - \mu \over \sigma} \right). \] Where \(\Phi\) is the standard normal distribution function.
Log-Normal partial expectation
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} \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{2.6} \]
Survival function of \(S_{+}\) for Log-Normal distribution
By substituting equation (2.6) into equation (2.3) we can solve for the survival function of \(S_+\) in terms of analytically available functions:
\[ \begin{aligned} Q_{S+}(t ;\mu, \sigma) &= Q_T(t + w_P;\mu, \sigma) + { 1 \over w_P} \int_t^{t+w_P} {1 \over z \sigma \sqrt{2\pi}} \exp\left( - {(\log(z) - \mu)^2 \over 2 \sigma^2} \right) (z-t)~ dz, \\ &= 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] \end{aligned} \]
Log-Normal discrete censored delay distribution
By substituting (2.6) into (2.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} \]
2.3 Weibull distribution
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.
Weibull partial expectation
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} \frac{kz}{\lambda}\left(\frac{z}{\lambda}\right)^{k-1}e^{-(z/\lambda)^{k}} dz \\ &= k\int_t^{t+w_P} \left(\frac{z}{\lambda}\right)^{k}e^{-(z/\lambda)^{k}} dz \\ &= \lambda k \int_{(t / \lambda)^k}^{((t + w_P) / \lambda)^k} y y^{1/k - 1} e^{-y} dy \\ &= \lambda\int_{(t / \lambda)^k}^{((t + w_P) / \lambda)^k} y^{1/k} e^{-y} dy\\ &= \lambda \Delta_{w_P} g(t; \lambda,k) \end{aligned} \tag{2.7} \]
Where
\[ g(t; \lambda, k) = \gamma(1 + 1/k, (t/\lambda)^k) = \frac{1}{k}\gamma(1/k, (t/\lambda)^k) - \frac{t}{\lambda}e^{-(t/\lambda)^k}. \] is a reparametrisation of the lower incomplete gamma function.
Survival function of \(S_{+}\) for Weibull distribution
By substituting equation (2.7) into equation (2.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]. \]
Weibull discrete censored delay distribution
By substituting (2.7) into (2.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]. # References