In this vignette, we’ll explain the underlying statistical model that
the primarycensored package. We’ll cover the following key
points:
primarycensored.If you are new to the package, we recommend that you start with the
vignette("primarycensored") vignette.
Time to event analysis, also known as survival analysis, concerns
estimating the distribution of delay times between events. A distinctive
feature of the field are the methodological techniques used to deal with
the missing data problems common in data sets of delay times[1]. In
primarycensored we focus on two particular missing data
problems:
In statistical epidemiology, these missing data problems occur frequently in both data analysis and theoretical modelling. For a more detailed description of these problems in an epidemiological context see[2,3].
In data analysis, events in epidemiology are commonly reported as occurring on a particular day or week (interval censoring). In an emerging outbreak, datasets can be incompletely observed (right truncation) and their can be a great deal of uncertainty around the precise timing of events (interval censoring).
In theoretical epidemiological modelling, it is often appropriate to
model the evolution of an infectious disease as occurring in discrete
time, for example in the EpiNow2 and EpiEstim
modelling packages. This means appropriately discretising continuous
distributions, such as the generation interval distribution. In
primarycensored we treat the discretisation of
intrinsically continuous distributions as an interval censoring
problem which allows us to simultaneously provide methods for both
applied and theoretical contexts.
primarycensoredAs described in Getting Started with
primarycensored, primarycensored focuses
on a subset of methods from time to event analysis that address data
missingness problems commonly found in epidemiological datasets. We
present the statistical problem as a double interval censoring problem,
where both the primary event time and the secondary event times are
interval censored. We can recover single interval censoring problems by
reducing one of the intervals to a point. In particular, a key
assumption we make is that the censoring window for events is known and
independent of the event time. This is known as non-informative
censoring.
The target for inference is the distribution of the delay time between the primary and secondary events. We assume that the delay time is a random variable T = S − Pu with distribution function FT(t) = Pr(T < t) and density function fT(t). In this treatment we assume that the delay time is shift-invariant, that is, the distribution of the delay time is the same regardless of the primary event time.
The (unconditional) primary event time is a random variable Pu with distribution function FPu(t) and density function fPu(t). The secondary event time is a random variable S, but in this treatment we construct the secondary event time from the primary time and delay, therefore the marginal distribution of S is not considered.
The censoring window for each event is the interval within which each event is known to have occurred, respectively, P ∈ [tP, tP + wP] and S ∈ [tS, tS + wS]. The lengths of the censoring windows are wP and wS respectively. The precise event times within their windows are unobserved. Note that since the primary event time is known up to the censoring window, we are predominantly interested in the conditional primary time P = Pu|{Pu ∈ [tP, tP + wP]} which has density function:
$$ \begin{aligned} f_{P}(p) &= {f_{P_{u}}(p) \over F_{P_{u}}(t_P + w_P) - F_{P_{u}}(t_P)}, \qquad &p \in [t_P, t_P + w_P],\\ f_{P}(p) &= 0, \qquad &\text{otherwise}. \end{aligned} $$
In this note, we measure the censored delay time Tc between the primary and secondary event windows from endpoint to endpoint: Tc = tS + wS − (tP + wP). Note that in the generative model for delays from “Getting started” the truncated delay tvalid is measured from startpoint to startpoint of event windows.
In our treatment below we focus on the survival function of the time after the end of the primary window to the secondary event time which we denote S+. We then use this to derive the distribution of the censored delay time Tc. This is equivalent to, but differs in mathematical approach from other treatments of the censoring problems in epidemiology, such as[2], see section Connections to other approaches for details.
In this section, we explain how to derive the distribution of the censored delay time Tc from the distribution of the delay time T and the condition distribution of the primary event time P.
In order to reason upon the distribution of the censored delay time Tc, it is useful to consider the time from the end (right) point of the primary censoring interval to the secondary time as a random variable,
S+ = S − (tP + wP) = T − ((tP + wP) − P) = T − CP.
Where T is the delay distribution of interest and CP = (tP + wP) − P is interval between the end (right) point of the primary censoring window and the primary event time; note that by definition CP is not observed but we can relate its distribution to the distribution of P: FCP(p) = Pr(CP < p) = Pr(P > tP + wP − p).
With non-informative censoring, it is possible to derive the upper distribution function of S+, or survival function of S+, from the distribution of T and the distribution of CP.
$$ \begin{equation} \begin{split} Q_{S_+}(t) &= Pr(S_+ > t) \\ &= Pr(T > C_P + t) \\ &= \mathbb{E}_{C_P} \Big[Q_T(t + C_P)\Big] \\ &= \int_0^{w_P} Q_T(t + p) f_{C_P}(p) dp. \end{split} \end{equation} $$
Using integration by parts gives: QS+(t) = QT(t + wP) + ∫0wPfT(t + p)FCP(p)dp.(#eq : survivalfunc)
Where we have used that QT′ = −fT, QT is the survival function of the actual delay distribution of interest and wP is the length of the primary censoring window.
Equation @ref(eq:survivalfunc) is the key equation in this note and is used to derive the distribution of the censored delay time Tc. It has the interpretation that the probability that the secondary event time is greater than t after the end of the primary censoring window is the sum of two disjoint event probabilities:
Note that in “Getting started” the target for numerical quadrature is the cumulative distribution function of the sum of the primary time within the primary censor window and the delay time.
Having constructed the survival function of S+ with equation @ref(eq:survivalfunc), using numerical quadrature or in some other way, we can calculate the probability mass of a secondary event time falling within a observed secondary censoring window of length wS that begins at time n − wS after the primary censoring window. This is the probability that the censored delay time Tc is n.
This gives the censored delay time probability by integrating over censored values:
Pr(S+ ∈ [n − wS, n)) = QS+(n − wS) − QS+(n).(#eq : seccensorprob)
Note that the censored secondary event time can also occur within the primary censoring window. This happens with probability, QS+(−wP) − QS+(0) = 1 − QT(wP) − ∫0wPfT(p)FCP(p)dp = Pr(T < C).
A common situation is when every primary and secondary event window is a single time unit, wP = wS = 1, and data arrives in non-overlapping intervals. For example, in the context of infectious disease modelling, we might have a data set of symptom onsets (primary event) and time of seeking health care (secondary event) which are both reported in days.
In this situation, we are interested in the probability of censored delays Tc in units of the event window. In this case equation @ref(eq:seccensorprob) is:
fn = Pr(Tc = n) = Pr(S+ ∈ [n − 1, n)) = QS+(n − 1) − QS+(n) n = 0, 1, ….(#eq : disccensprob)
Which is a discrete probability mass function of the secondary event time relative to the primary event time in units of the event window.
In this section, we discuss how the approach taken in
primarycensored relates to some other approaches to the
censored data problems in time to event analysis.
Using the notation from Park et al[2], we write the conditional probability of the secondary event time S ∈ (SL, SR) given the primary event time P ∈ (PL, PR) as:
$$ \begin{aligned} \mathrm{Pr}(S_L < S < S_R | P_L < P < P_R) &= \frac{\mathrm{Pr}(P_L < P < P_R, S_L < S < S_R)}{\mathrm{Pr}(P_L < P < P_R)} \\ &= \frac{\int_{P_L}^{P_R} \int_{S_L}^{S_R} g_P(x) f_x(y-x) dy dx}{\int_{P_L}^{P_R} g_P(x) dx}\\ &= \int_{P_L}^{P_R} \int_{S_L}^{S_R} g_P(x|P_L, P_R) f_x(y-x)dy dx \end{aligned} $$
In this note, we assume that the forward distribution doesn’t vary over time (such that fx = f), then
∫PLPR∫SLSRgP(x|PL, PR)fx(y − x)dydx = ∫PLPRgP(x|PL, PR)[F(SR − x) − F(SL − x)]dx
Then, by using integration by parts, we get:
$$ \begin{split} \int_{P_L}^{P_R} g_P(x|P_L, P_R) \big[F(S_R - x) - F(S_L - x)\big] dx &= F(S_R - P_R) - F(S_L - P_R) \\ & - \int_{P_L}^{P_R} G_P(x|P_L, P_R) \big[f(S_L - x) - f(S_R - x)\big] dx \end{split} (\#eq:park) $$ Where we have used that ∂xF(SR − x) = −f(SR − x) and ∂xF(SL − x) = −f(SL − x).
We can now compare this to equation @ref(eq:seccensorprob) by considering the following transformations:
Then equation @ref(eq:park) becomes:
$$ \begin{aligned} \mathrm{Pr}(S_L < S < S_R | P_L < P < P_R) &= F(n) - F(n-w_S) - \int_{-w_P}^{0} G_P(x|-w_P, 0) \big[f(n - w_S - x) - f(n - x)\big] dx \end{aligned} $$ Making the transformation x = −p, and rewriting in the notation of this note gives: $$ \begin{aligned} &= F(n) - F(n-w_S) + \int_{w_P}^{0} G_P(-p|-w_P, 0) \big[f_T(n - w_S + p) - f_T(n +p)\big] dp \\ &= F(n) - F(n-w_S) + \int_{0}^{w_P} G_P(-p|-w_P, 0) \big[f_T(n + p) - f_T(n - w_S +p)\big] dp\\ &= F(n) - F(n-w_S) + \int_{0}^{w_P} (1 - F_{C_P}(p)) \big[f_T(n + p) - f_T(n - w_S +p)\big] dp\\ &= F(n + w_P) - F(n-w_S + w_P) + \int_{0}^{w_P} [f_T(n + p - w_S) - f_T(n + p)] F_{C_P}(p) dp\\ &= Q_T(n-w_S + w_P) - Q_T(n + w_P) + \int_{0}^{w_P} [f_T(n + p - w_S) - f_T(n + p)] F_{C_P}(p) dp \\ &= Q_{S_+}(n-w_S) - Q_{S_+}(n ). \end{aligned} $$ which is same as equation @ref(eq:seccensorprob).
In this derivation, we have used that GP(x|−wP, 0) is the distribution function from the time from the start of the primary interval until primary event time, and FCP is the distribution function of the time until the end of the primary event window from the primary event time. Therefore, GP(−p|−wP, 0) = Pr(P < −p|P ∈ (−wP, 0)) = 1 − Pr(CP ≤ p) = 1 − FCP(p).
vignette("analytic-solutions").vignette("primarycensored").