--- title: "Analytic solutions for censored delay distributions" format: html bibliography: library.bib output: bookdown::html_document2: fig_caption: yes code_folding: show pkgdown: as_is: true csl: https://raw.githubusercontent.com/citation-style-language/styles/master/apa-numeric-superscript-brackets.csl link-citations: true vignette: > %\VignetteIndexEntry{Analytic solutions} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- # What are we going to do in this Vignette In this vignette, we'll derive the analytic solutions for the primary censored delay distributions for a range of commonly used distributions. 1. Examine the case of exponentially tilted primary event times 2. Analyse the uniform primary event time case ($r=0$) 3. Explore general partial expectation and its role in solving equation \@ref(eq:unifprim) 4. Provide specific analytic solutions for distributions where the partial expectation can be reduced to an analytic expression # 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]. (\#eq:fcp) $$ Note that taking the limit $\lim_r \rightarrow 0$ in equation \@ref(eq:fcp) 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). $$ # Uniform primary event time ($r=0$) Applying a uniform primary event time distribution to equation 3.1 from ["Why it works"](why-it-works.html) 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} (\#eq:unifprim) $$ ## General partial expectation Note that for any distribution with an analytically available distribution function $F_T$ equation \@ref(eq:unifprim) can be solved so long as the _partial expectation_ $$ \int_t^{t+w_P} f_T(z) z~ dz (\#eq:partexp) $$ 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"](why-it-works.html) can be written using the difference operator: $f_n = -\Delta_1 Q_{S_+}(n-1)$. We can insert this expression into equation \@ref(eq:unifprim) 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} (\#eq:disccensunifprim) $$ We now consider some specific delay distributions. ## 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} \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} (\#eq:gammapartexp) $$ ### Survival function of $S_{+}$ for Gamma distribution By substituting equation \@ref(eq:gammapartexp) into equation \@ref(eq:disccensunifprim) 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]. (\#eq:survgammaunifprim) $$ ### Gamma discrete censored delay distribution By substituting \@ref(eq:survgammaunifprim) into \@ref(eq:disccensunifprim) 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_ [@cori2013new]. ## Log-Normal distribution 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. ### 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} \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} (\#eq:lognormpartexp) $$ ### Survival function of $S_{+}$ for Log-Normal distribution By substituting equation \@ref(eq:lognormpartexp) into equation \@ref(eq:disccensunifprim) 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] $$ ### Log-Normal discrete censored delay distribution By substituting \@ref(eq:lognormpartexp) into \@ref(eq:disccensunifprim) 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} $$ ## 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} \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} (\#eq:weibullpartexp) $$ 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. ## Survival function of $S_{+}$ for Weibull distribution By substituting equation \@ref(eq:weibullpartexp) into equation \@ref(eq:disccensunifprim) 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 \@ref(eq:weibullpartexp) into \@ref(eq:disccensunifprim) 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_ [@cori2013new]. # Learning more - For more mathematical background on the analytic solutions see the `vignette("why-it-works")`. - For a more introductory explanation of the primary event censored distribution see the `vignette("primarycensored")`. # References