---
title: "Seroincidence package methodology"
author: "European Centre for Disease Prevention and Control (ECDC)"
date: '`r Sys.Date()`'
output:
  pdf_document:
    toc: yes
header-includes:
- \usepackage{subfigure}
- \usepackage[utf8]{inputenc}
references:
- DOI: 10.1002/sim.5322
  author:
  - family: Teunis
    given: P. F.
  - family: van Eijkeren
    given: J. C.
  - family: Ang
    given: C. W.
  - family: van Duynhoven
    given: Y. T.
  - family: Simonsen
    given: J. B.
  - family: Strid
    given: M. A.
  - family: van Pelt
    given: W.
  container-title: Statistics in Medicine
  id: teunis2012
  issue: 20
  issued:
    month: 9
    year: 2012
  page: 2240-2248
  title: 'Biomarker dynamics: estimating infection rates from serological data'
  type: article-journal
  volume: 31
  
- author:
  - family: Versteegh
    given: F. G.
  - family: Mertens
    given: P. L.
  - family: de Melker
    given: H. E.
  - family: Roord
    given: J. J.
  - family: Schellekens
    given: J. F.
  - family: Teunis
    given: P. F.
  container-title: Epidemiology and Infection
  id: versteegh2005
  issue: 4
  issued:
    month: 8
    year: 2005
  page: 737-748
  title: Age-specific long-term course of IgG antibodies to pertussis toxin after
    symptomatic infection with Bordetella pertussis
  type: article-journal
  volume: 133
  
- DOI: 10.1002/sim.3592
  author:
  - family: Simonsen
    given: J.
  - family: Mølbak
    given: K.
  - family: Falkenhorst
    given: G.
  - family: Krogfelt
    given: K. A.
  - family: Linneberg
    given: A.
  - family: Teunis
    given: P. F.
  container-title: Statistics in Medicine
  id: simonsen2009
  issue: 14
  issued:
    month: 6
    year: 2009
  page: 1882-1895
  title: Estimation of incidences of infectious diseases based on antibody measurements
  type: article-journal
  volume: 28
  
- author:
  - family: Feller
    given: W.
  container-title: John Wiley & Sons
  id: feller1968
  issued:
    year: 1968
  title: An introduction to probability theory and its applications
  type: book
  volume: 2
  
- author:
  - family: Higgins
    given: P. M.
  container-title: 'London: Springer Verlag'
  id: higgins2008
  issued:
    year: 2008
  title: Number Story. From counting to cryptography
  type: book
  
- author:
  - family: Scalia Tomba
    given: G. P.
  - family: Svenson
    given: Å.
  - family: Asikainen
    given: T.
  - family: Giesecke
    given: J.
  container-title: Mathematical Biosciences
  id: scalia2010
  issue: 1
  issued:
    year: 2010
  page: 24-31
  title: Some model based considerations on observing generation times for communicable
    diseases
  type: article-journal
  volume: 223
fontsize: 11pt
---





# 1. Introduction

The simplest possible model for serum antibody decay is an exponential
decay function (see [Appendix](#method_longitudinal)). This is the default 
longitudinal model used in **seroincidence** package. 
Various modifications can be made, however, to allow for
different interpretations of low antibody concentrations in cross-sectional data. 
The following three approaches are implemented in the current version of the package:

1. Exponential decay: seroconversion is instantaneous, and is followed by 
exponential decrease to zero (see [Appendix](#method_longitudinal) for more information).
  
2. Exponential decay with left censoring: seroconversion is instantaneous, 
and it is followed by exponential decay, but antibody concentrations below a
given (fixed) level are treated as censored (see [Appendix](#censoring) for 
more information). This implies that observations below censoring level are 
considered to have been  generated by seroconversions from at least a 
given time in the past (corresponding to the censoring level). In simple 
terms: for very low antibody concentrations we can only conclude that 
seroconversion occurred long ago, but not precisely how long.

3. Exponential decay with nonzero baseline: seroconversion is again instantaneous, 
and is again followed by exponential decay, but decay now is not towards zero 
but instead, approaches some baseline antibody concentration that may be 
different for any individual (see [Appendix](#baseline) for more information). 
Low antibody concentrations in the cross-sectional sample now may correspond to
baseline levels, in which case seroconversion would have occurred an infinite period ago.

Because of their different interpretations of low cross-sectional
serum antibody concentrations, both alternatives (left censoring and
decay to baseline) lead to decreased estimates of seroincidence,
compared to the default model.

Finally, note that the methods described in the [Appendix](#appendices) allow for
arbitrary antibody decay patterns, so that adaptations can be made to
accommodate any longitudinal model for the shape of the decaying arm
of the seroresponse.

# 2. Influence of censoring on seroincidence

In case serum antibody data are known to be censored, a cutoff level
can be defined in the seroincidence calculator script as a list, with
separate values for each of the three antibody classes (IgG, IgM,
and/or IgA).

\begin{figure}[h!]
\centering
\mbox{
\subfigure[Same cutoff]{
\includegraphics[width=0.48\linewidth]{./fig/cutoff-campylobacter-1.png}}
\subfigure[Different cutoffs]{
\includegraphics[width=0.48\linewidth]{./fig/cutoff-campylobacter-2.png}}}
\caption{\label{censoring_campy}Influence of cutoff levels on
  seroincidence estimates for a cross-sectional sample of {\em Campylobacter} antibodies. 
  (a) Identical cutoffs for all three antibody types; 
  (b) Different cutoffs for IgG, IgM and IgA.}
\end{figure}

\begin{figure}[h!]
\centering
\mbox{
\subfigure[Same cutoff]{
\includegraphics[width=0.48\linewidth]{./fig/cutoff-salmonella-1.png}}
\subfigure[Different cutoffs]{
\includegraphics[width=0.48\linewidth]{./fig/cutoff-salmonella-2.png}}}
\caption{\label{censoring_salmo}Influence of cutoff levels on
  seroincidence estimates for a cross-sectional sample of {\em Salmonella} antibodies (mixed ELISA). 
  (a) Identical cutoffs for all three antibody types; 
  (b) Different cutoffs for IgG, IgM and IgA.}
\end{figure}

Figures \ref{censoring_campy} and \ref{censoring_salmo} show the
change in seroincidence estimate, for a given set of cross-sectional
antibody data (*Campylobacter* and *Salmonella*), with
increasing cutoff levels. Cutoffs need not be identical for IgG, IgM
and IgA antibodies, as illustrated in these graphs.

# 3. Checks using the longitudinal data

The longitudinal data may be used to obtain a crude estimate of bias
in the estimated incidences (Figure \ref{verif}a, \ref{verif}e).

\begin{figure}[h!]
\centering
\mbox{
\subfigure[Cutoff 0.0]{
\includegraphics[width=0.24\linewidth]{./fig/delft-cutoff0.png}}
\subfigure[Cutoff 0.5]{
\includegraphics[width=0.24\linewidth]{./fig/delft-cutoff05.png}}
\subfigure[Cutoff 1.0]{
\includegraphics[width=0.24\linewidth]{./fig/delft-cutoff1.png}}
\subfigure[Cutoff 2.0]{
\includegraphics[width=0.24\linewidth]{./fig/delft-cutoff2.png}}}
\mbox{
\subfigure[Cutoff 0.0]{
\includegraphics[width=0.24\linewidth]{./fig/ssi-cutoff0.png}}
\subfigure[Cutoff 0.5]{
\includegraphics[width=0.24\linewidth]{./fig/ssi-cutoff05.png}}
\subfigure[Cutoff 1.0]{
\includegraphics[width=0.24\linewidth]{./fig/ssi-cutoff1.png}}
\subfigure[Cutoff 2.0]{
\includegraphics[width=0.24\linewidth]{./fig/ssi-cutoff2.png}}}
\caption{\label{verif}Estimated incidences for longitudinal sera at
  different times since onset of illness symptoms. delta(incid) =
  difference between estimated incidence and 365/(time since
  onset). Top: Campylobacter, Delft ELISA, Bottom: Campylobacter, SSI
  ELISA.}
\end{figure}

It is seen that for samples taken 100--250 days post symptom onset, the
estimated incidences appear unbiased (the 95\% interval includes
0). Increasing the cutoff (here applied identical values for all three
antibody classes) causes underestimation of the incidence (Figures
\ref{verif}b,c,d and \ref{verif}f,g,h). The high degree of variation
in antibody levels in late samples (more than 100 days post symptom
onset) renders estimation problematic.

\newpage

\newpage

# Appendices {#appendices}

This is a slightly modified description from [@teunis2012], extended with adaptations for
censoring and elevated baseline levels.

## A.1. Characteristics of the serum antibody response {#method_longitudinal}

The longitudinal model used for describing serum antibody responses
has been published [@simonsen2009, @teunis2012]. As
longitudinal data may include early samples taken shortly after
infection, possibly during the initial increase in antibody
concentrations (during seroconversion) the longitudinal model includes
the rising phase of the seroresponse. This avoids misspecification of
the peak antibody level, resulting in estimates for the (joint
distribution of) peak levels, time to peak, and decay rates of the
serum antibody concentration [@teunis2012, @versteegh2005].

However, seroconversion usually is rapid: the increase in antibody
levels is orders of magnitude faster than the subsequent decrease. For
seroconversion rate calculations we ignore the time to reach peak
levels. For simplicity it will be assumed that the antibody response
decays exponentially but this is not an essential assumption and it
can be relaxed to any monotonically decreasing function of time. So,
the antibody response $y$ to incidence at $t = 0$ is an instantaneous
increase to peak level ($A$) followed by an exponential decrease
(with rate $k$):

$$
y = A \text{e}^{-kt}
$$

## A.2. Random sampling of randomly repeating responses {#A2}

In a cross-sectional sample of responses, each subject will be sampled
at a random time within its current inter-event interval $\Delta t$
(interval, for short). The within--interval sampling distribution is
uniform:

\begin{equation}
\label{fundamental}
u_{f}(\tau\vert\Delta t)=\frac{[0\le \tau < \Delta t]}{\Delta t}
\end{equation}

The subscript $f$ is introduced to show that this distribution is
fundamental to sampling from the period between events that are
separated by a constant time interval $\Delta t$. As the intervals
$\Delta t$ generally will vary, the within--interval sampling
distribution and the distribution of intervals combine to a population
distribution of backward recurrence times. The backward recurrence
time distribution can be translated into a corresponding distribution
of response values (antibody concentrations), due to the monotonicity
of the response.

The intervals $\Delta t$ are supposed to vary randomly with
distribution $p(\Delta t)$. This introduces a subtle question about
the probability of sampling the response within an infinitesimal range
of intervals $[\Delta t, \Delta t + \text{d}\Delta t]$. One would
be inclined to assume that the probability is given only by the
probability density $p(\Delta t)$, but one should keep in mind that
responses with longer intervals are more likely to be sampled than
those with shorter intervals. Obviously, the probability of the range
of sampling intervals $[\Delta t, \Delta t + \text{d}\Delta t]$
should be proportional to the interval $\Delta t$ itself (size biased
sampling, [@feller1968, @higgins2008, @scalia2010]). Then,
the probability of sampling a response within an infinitesimal range
of intervals $[\Delta t, \Delta t + \text{d}\Delta t]$ is

$$
\frac{p(\Delta t)\Delta t}{\overline{\Delta t}_{p}}\text{d}\Delta t
$$

Here $\overline{\Delta t}_{p}$ is the expected interval length from
the distribution $p$. We define the incidence as its reciprocal
$1/\overline{\Delta t}_{p}$.

Combining this probability with the uniform distribution of sampling
time within the interval $\Delta t$ leads to a population distribution
of time since last event (backward recurrence time)

\begin{eqnarray*}
u(\tau) &=& \int_{\Delta t=0}^{\infty}u_{f}(\tau\vert\Delta t) \frac{p(\Delta t)\Delta t}{\overline{\Delta t}_{p}}\text{d}\Delta t\\
        &=& \frac{1}{\overline{\Delta t}_{p}} \int_{\Delta t=\tau}^{\infty} p(\Delta t)\text{d}\Delta t = \frac{1-P(\tau)}{\overline{\Delta t}_{p}}
\end{eqnarray*}

where $P$ is the cumulative distribution of $p$. In going from the
first to the second equality we used the fact that $\Delta t$ should
exceed $\tau$ (equation (\ref{fundamental})). Note that the
distribution $u(\tau)$ is monotonically non-increasing.

## A.3. Delta distribution {#A3}

If $p(\Delta t)$ is the delta-distribution $p(\Delta t) = \delta(\Delta t-\Delta t^{\ast})$ 
representing identical intervals, then the uniform distribution $u_{f}(\tau\vert \Delta t^{\ast})$ 
is regained: for the delta distribution $p(\Delta t) = \delta(\Delta t -\Delta t^{\ast})$ 
the expected value is $\overline{\Delta t}_{p} = \Delta t^{\ast}$. It follows that

$$
u(\tau\vert\Delta t^{\ast}) 
= \frac{1}{\Delta t^{\ast}} \int_{\tau}^{\infty} \delta (\Delta t-\Delta t^{\ast}) \text{d}\Delta t 
= \frac{\left[0\le \tau < \Delta t^{\ast}\right]}{\Delta t^{\ast}}
$$

## A.4. Distribution of response values given the distribution of backward recurrence times {#A4}

Given incidence at $\tau = 0$, the response decays as $y(\tau) = A \text{e}^{-k\tau}$, 
with $\tau$ backward recurrence time. Consider

\begin{eqnarray}
P(y^{\prime}\le y) & = & P\left(A\text{e}^{-k\tau^{\prime}}\le A\text{e}^{-k\tau}\right) [y\le A] = P(\tau \le \tau^{\prime})[y\le A] \nonumber \\
                   & = & \left(1-U(\tau)\right)[y\le A] = \left(1-U(\log(A/y)/k)\right) [y\le A] \label{cumul_distr}
\end{eqnarray}

Differentiating this expression leads to

\begin{eqnarray}
\rho(y) & = & \frac{\text{d}}{\text{d}y} P(y^{\prime}\le y) = -u\left(\log(A/y)/k\right) \left(\frac{1}{k}\right) \left(\frac{-1}{y}\right)[y\le A] \nonumber\\
        & = & \frac{[y\le A]}{ky}u\left(\log(A/y)/k\right)\label{distrib_y}
\end{eqnarray}

Given the antibody response $y$ with respect to time of incidence,
$\rho(y)$ is the corresponding density for sampled response values.

For the $n$th moment of the backward recurrence time one derives

\begin{eqnarray*}
E(\tau^{n}) & = & \int_{\tau = 0}^{\infty} \tau^{n} u(\tau)\text{d}\tau 
                    = \int_{\tau = 0}^{\infty}\left(\tau^{n}\frac{1}{\overline{\Delta t}_{p}} 
                    \int_{\Delta t = \tau}^{\infty} p(\Delta t)\text{d}\Delta t\right)\text{d}\tau \\
            & = & \frac{1}{\overline{\Delta t}_{p}}\int_{\Delta t=0}^{\infty}p(\Delta t) 
                    \left(\int_{\tau = 0}^{\Delta t}\tau^{n}\text{d}\tau\right)\text{d}\Delta t \\
            & = & \frac{1}{\overline{\Delta t}_{p}}\int_{\Delta t=0}^{\infty}p(\Delta t) 
                    \frac{(\Delta t)^{{n+1}}}{n+1}\text{d}\Delta t 
                    = \frac{\overline{(\Delta t)^{n+1}}_{p}}{(n+1)\overline{\Delta t}_{p}}
\end{eqnarray*}

## A.5. Family of Gamma distributions {#A5}

The distribution of intervals $p(\Delta t)$ may have any
form. Let the family of $\Gamma(\lambda, m)$ distributions be given:

$$
p_{m}(\Delta t\vert \lambda) = \frac{\lambda^{m+1}\Delta t^{m}}{m!}\text{e}^{-\lambda \Delta t}
$$

From

$$
\int_{x=0}^{\infty}x^{m+1}\text{e}^{-x}\text{d}x = (m+1)!
$$

one easily derives that $\overline{\Delta t}_{p} = (m+1)/\lambda$.

We shall need the cumulative distribution

\begin{eqnarray}
P_{m}(\tau\vert \lambda) & = & \int_{\Delta t=0}^{\tau} p_{m}(\Delta t\vert \lambda) \text{d}\Delta t 
                                = 1 - \int_{\Delta t=\tau}^{\infty}p_{m}(\Delta t\vert \lambda)d\Delta t \nonumber \\ 
                         & = & 1 - \frac{\Gamma(m+1,\lambda \tau)}{m!} \label{p_dist_gamma}
\end{eqnarray}

where $\Gamma(m + 1, \lambda\tau)$ denotes the (upper) incomplete gamma
function (not to be confused with the gamma distribution).

Note: working out the integral

$$
P_{m}(\tau\vert\lambda) = \int_{\Delta t = 0}^{\tau} p_{m}(\Delta t\vert\lambda)\text{d}\Delta t,
$$

it is straightforward to show that

\begin{equation}
\label{gamma_pj}
\Gamma(m+1,\lambda\tau) = \frac{m!}{\lambda}\sum_{j=0}^{m}p_{j}(\tau\vert\lambda)
\end{equation}

For the distribution of backward recurrence times one finds

\begin{equation}
\label{u_dist_gamma}
u_{m}(\tau\vert \lambda) = \frac{\lambda}{m+1} \int_{\Delta t=\tau}^{\infty}p_{m}(\Delta t\vert \lambda)\text{d}\Delta t 
    = \frac{\lambda}{(m+1)!}\Gamma(m+1,\lambda\tau)
\end{equation}

when $m = 0$

$$
u_{0}(\tau\vert \lambda) = \lambda \Gamma(1,\lambda\tau) = p_{0}(\tau\vert \lambda),
$$

the famous property of the Poisson process.

## A.6. Heterogeneity in serum antibody responses {#A6}

Incorporation of responses that also vary in amplitude $A$ and decay
rate $k$ can be easily described with a joint density $g(A, k)$. Of
course, such heterogeneity changes the population based distribution
of $y$ to

$$
\rho(y) = \frac{1}{y}\int_{k=0}^{\infty}\int_{A=y}^{\infty} \frac{1}{k}u\left(\frac{\log(A/y)}{k}\right)g(A,k)\text{d}A\text{d}k
$$

instead of corresponding equation (\ref{distrib_y}). Note that the factor $[y\le A]$ 
in equation ($\ref{distrib_y}$) forces the lower limit $y$ in the
integration over $A$. In case the peak antibody levels $A$ and
antibody decay rates $k$ are available as a Monte Carlo sample
$\{A_{n},k_{n}\}$ of their joint distribution the integration can be
approximated by

\begin{equation}
\label{mc_heterogeneity}
\rho(y)=\frac{1}{y}\frac{1}{N}\sum_{n=1}^{N}\frac{1}{k_{n}} u\left(\frac{\log(A_{n}/y)}{k_{n}}\right)[y\le A_{n}]
\end{equation}

where $N$ is the size of the Monte-Carlo sample, maintaining any
correlation between peak levels $A_{n}$ and decay rates $k_{n}$.
Assuming e.g. that infections occur with Gamma distributed intervals
with rate $\lambda$ and shape factor $m$:

\begin{eqnarray*}
\rho(y\vert \lambda, m) & = & \frac{1}{(m+1) N}\sum_{n=1}^{N}\left( \frac{\lambda}{k_{n}A_{n}} \left(\frac{y}{A_{n}}\right)^{-1+\lambda/k_{n}}\right.\\ 
                        &\times & \left.\sum_{j=0}^{m} \frac{\left(\lambda/k_{n} \log(A_{n}/y)\right)^{j}}{j!}[y\le A_{n}]\right)
\end{eqnarray*}

For a cross-sectional sample $\{Y_{1}, Y_{2},\dots , Y_{N_{c}}\}$ a
likelihood

$$
\ell (\lambda, m) = \prod_{i=1}^{N_{c}} \rho(Y_{i}\vert \lambda, m)
$$

can be calculated to allow estimation of the parameters $(\lambda, m)$ of 
the process generating infections. The corresponding incidence is
the average rate of infections $\lambda/(m+1)$.

## A.7. Censored observations {#censoring}

In case observations are censored at $y_{c}$ such that an observed $Y
= \text{max}(Y, y_{c})$, then for $y_{c} < Y$ the density $\rho(y)$
as in eq. ($\ref{distrib_y}$) holds, but the likelihood of any $Y \le y_{c}$
can be calculated

$$
\ell(\lambda, m \vert y \le y_{c}) = \int_{z=0}^{y_{c}}\rho(z\vert \lambda,m)\text{d}z = R(y_{c}\vert \lambda,m)
$$

and the likelihood function for a cross-sectional sample $\{Y_{1}, Y_{2}, \dots, Y_{N_{c}}\}$ becomes

$$
\ell (\lambda, m) = \prod_{i=1}^{N_{c}} \left(\rho(Y_{i}\vert \lambda, m)[y_{c} < Y_{i}] \times R(y_{c}\vert \lambda,m)[Y_{i}\le y_{c}]\right)
$$

We need the cumulative distribution $R(y)$ as given by
eq. ($\ref{cumul_distr}$)

\begin{equation}
\label{r_cumul_def}
R(y) = 1 - U\left(\frac{\log\left(\frac{A}{\min(y,A)}\right)}{k}\right) = 1 - U\left(\frac{\text{lm}(y,A)}{k}\right)
\end{equation}

with obvious definition of $\text{lm}(y, A)$.

For the family of gamma distributions for the intervals $\Delta t$,
eq. ($\ref{r_cumul_def}$) leads to

$$
R_{m}(y) = 1 - U_{m}\left(\frac{\text{lm}(y,A)}{k}\vert \lambda\right)
$$

Substituting eq. ($\ref{gamma_pj}$) into eq. ($\ref{u_dist_gamma}$), and
using the cumulative distribution $P_{j}$

$$
R_{m}(y) = 1 - \frac{1}{m+1}\sum_{j=0}^{m}P_{j} \left(\frac{\text{lm}(y,A)}{k}\vert\lambda\right)
$$

which can also be expressed in terms of incomplete Gamma
functions, using eq. ($\ref{p_dist_gamma}$)

$$
R_{m}(y) = 1 - \frac{1}{m+1}\sum_{j=0}^{m} 1 - \frac{\Gamma\left(j+1,\lambda\frac{\text{lm}(y,A)}{k}\right)}{j!}
$$

or, expanding $\text{lm}$

$$
R(y) = \frac{1}{m+1}\sum_{j=0}^{m}\frac{1}{j!}\Gamma\left(j+1,  \frac{\lambda}{k}\log\left(\frac{A}{\min(y,A)}\right)\right)
$$

The latter form is useful for numeric calculations as there is a
reliable and efficient implementation available in R.

## A.8. Nonzero baseline {#baseline}

The procedure outlined in [@teunis2012] is valid
for any monotonically non-increasing serum antibody response. If for
instance, after seroconversion, antibody levels decrease to some
baseline level greater than zero (however long ago seroconversion has
occurred), the antibody response is

$$
y = y_{0} + A \text{e}^{-kt}
$$

Such responses have been assumed for back-calculations [@simonsen2009].

Now the distribution of response values given the distribution of the
backward recurrence times becomes

\begin{eqnarray*}
P(y^{\prime}\le y) & = & P\left(y_{0}+A\text{e}^{-k\tau^{\prime}}\le y_{0}+A\text{e}^{-k\tau}\right)[y_{0} < y \le A] \\
                   & = & P(\tau \le \tau^{\prime})[y_{0} < y \le A] = \left(1 - U(\tau)\right)[y_{0} < y \le A] \\
                   & = & \left(1-U\left(\frac{1}{k}\log\left(\frac{A}{y-y_{0}}\right)\right)\right) [y_{0} < y \le A]
\end{eqnarray*}

Differentiation produces the density

$$
\rho(y) = \frac{\text{d}}{\text{d}y} P(y^{\prime}\le y) = \frac{1}{k(y-y_{0})} u\left(\frac{1}{k}\log\left(\frac{A}{y-y_{0}}\right)\right)
$$

Including heterogeneity in $A$, $k$ and $y_{0}$ via their joint distribution

$$
\rho(y) = \int_{k=0}^{\infty}\int_{A=y}^{\infty}\int_{y_{0}=0}^y 
    \frac{1}{k(y-y_{0})} u\left(\frac{1}{k}\log\left(\frac{A}{y-y_{0}}\right)\right) 
    g(A,K,y_{0})\text{d}y_{0}\text{d}A\text{d}k
$$

or, when a Monte-Carlo sample of $\{A_{n}, k_{n}, y_{0,n}\}$ triplets
is available

$$
\rho(y) = \frac{1}{N}\sum_{n=1}^{N}\frac{u\left(\frac{1}{k_{n}}  \log\left(\frac{A_{n}}{y-y_{0,n}}\right)\right)}{k_{n}(y-y_{0,n})}  [y_{0,n} < y < A_{n}]
$$

For the gamma family of distributions for $\Delta t$,

\begin{eqnarray*}
\rho(y) & = & \frac{1}{N(m+1)}\sum_{n=1}^{N}\left(\frac{\lambda}{k_{n}A_{n}} \left(\frac{y-y_{0,n}}{A_{n}}\right)^{-1+\lambda/k_{n}}\right. \\
        & \times & \left.\sum_{j=0}^{m}\frac{1}{j!}\left(\frac{\lambda}{k_{n}} \log\left(\frac{A_{n}}{y-y_{0,n}}\right)\right)^{j} [y_{0,n} < y \le A_{n}]\right)
\end{eqnarray*}

For the Poisson process ($m = 0$) this reduces to

$$
\rho(y) = \frac{1}{N}\sum_{n=1}^{N}\frac{\lambda}{k_{n}A_{n}} \left(\frac{y-y_{0,n}}{A_{n}}\right)^{-1+\lambda/k_{n}}[y_{0,n} < y \le A_{n}]
$$

# References
