\begin{document}

	
	\title{Supplemental Material to: \\ A Bayesian accelerated failure time model for interval censored three-state screening outcomes}
	\author{Thomas Klausch (\href{t.klausch@amsterdamumc.nl}{t.klausch@amsterdamumc.nl})}
	\author[1] {\\Eddymurphy U Akwiwu}
	\author[1] {Mark A. van de Wiel}
	\author[1] {\mbox{Veerle M. H. Coup\'e}}
	\author[1] {\mbox{Johannes Berkhof}}
	\affil[1]{Amsterdam University Medical Centers, \mbox{Department of Epidemiology and Data Science}, Amsterdam, The Netherlands} 

\maketitle

\tableofcontents

\section{Supplement A: Proofs and Information Criteria} 

\subsection{Proof for Section 2 (Informative Interval Censoring)} \label{sec::infcens}
We first define informative interval censoring in the standard interval censored survival analysis setting. For this consider the random survival time variable $Y$ with density $f_Y(y)$ which is observed for $n$ units $i=1,...,n$, subject to interval censored observations at times $\vf_i = (0, v_{i,1},...,v_{i,m_i}, \infty )$ with first and last entries set to zero and infinity, respectively, and $m_i$ random. In particular, the censoring interval $\phi_i = (l_i,r_i]$ is defined such that
\begin{align*}
	\Pr( l_i = v_{i,j}, r_i = v_{i,j+1} | y_i, \vf_i ) = \mathds{1}_{\{ v_{i,j} < y_i \le v_{i,j+1} \}}.
\end{align*}

Non-informative interval censoring is defined as \citep[][p. 79]{kalbfleisch_statistical_2002}
\begin{align} \label{eq::noninf}
	\Pr( y_i \in (l_i, r_i] | l_i, r_i) = \Pr( y_i \in (l_i, r_i]).
\end{align}

This situation occurs when the screening times $\vf$ are generated independently of $y_i$ or are fixed in advance of data collection. Informative censoring occurs when (\ref{eq::noninf}) does not hold. Standard survival analysis methods for interval censored data usually assume non-informative censoring and yield biased estimates and wrong inference when (\ref{eq::noninf}) does not hold. \\

Now consider the data generating mechanism of the three-state model for screening outcomes defined in Section 3 of the main manuscript. We have that $y_i = x_i + t_i$ and 
\begin{align*}  
	\Pr( l_i = v_{i,j}, r_i = v_{i,j+1} | x_i, \vf_i ) = \mathds{1}_{\{ v_{i,j} < x_i \le v_{i,j+1} \}}
\end{align*}	
with censoring states 
\begin{align*}  
	\delta_i = \begin{cases}
		1 \ &\text{if} \ \ l_i < x_i < r_i = \infty  \\
		2 \ &\text{if} \ \ l_i < x_i \le r_i < y_i \\
		3 \ &\text{if} \ \ l_i < x_i < y_i \le r_i < \infty,
	\end{cases}
\end{align*}
where we assume \textit{non-informative interval censoring in the three state model}, i.e.,
\begin{align} \label{eq::noninf3s}
	q(\vf_i | x_i, t_i) = q(\vf_i).
\end{align}

In the main manuscript and Section \ref{sec::proof4.1} of the Supplemental Material we additionally show that the estimation methods discussed in the main manuscript only require (\ref{eq::noninf3s}) to hold conditionally on a vector of covariates $\z_i$, which are suppressed here for brevity. \\

A naive approach to using the three-state data for modeling the time distribution of $y_i$ first transforms the censoring times $\phi_i$ to standard interval censoring of $y_i$ by the rules
\begin{align*}
	\tilde{l}_i = \tilde{l}(l_i,r_i,\delta_i, y_i) &= \begin{cases}
		l_i\ &\text{if} \ \ \delta_i = 1 \\
		r_i\ &\text{if} \ \ \delta_i = 2 \\
		l_i\ &\text{if} \ \ \delta_i = 3 
	\end{cases} \\
   &= \begin{cases}
		l_i\ &\text{if} \ \ r_i = \infty \\
		r_i\ &\text{if} \ \ r_i < y_i \\
		l_i\ &\text{if} \ \ y_i \le r_i 
	\end{cases}
\end{align*}
and 
\begin{align*}
	\tilde{r}_i = \tilde{r}(r_i,\delta_i,y_i) &= \begin{cases}
		\infty\ &\text{if} \ \ \delta_i = 1 \\
		\infty\ &\text{if} \ \ \delta_i = 2 \\
		r_i\ &\text{if} \ \ \delta_i = 3 
	\end{cases} \\ 
    &= \begin{cases}
	\infty\ &\text{if} \ \ r_i = \infty \\
	\infty\ &\text{if} \ \ r_i < y_i \\
	r_i\ &\text{if} \ \ y_i \le r_i.
\end{cases}
\end{align*}

However, we then have that

\begin{align*}
	\Pr( y_i \in (\tilde{l}_i, \tilde{r}_i] | \tilde{l}_i, \tilde{r}_i) \neq \Pr( y_i \in ( \tilde{l}_i,  \tilde{r}_i])
\end{align*}
in violation of condition (\ref{eq::noninf}) due to the dependence of $\tilde{l}(l_i,r_i,\delta_i, y_i) $ and $\tilde{r}(r_i,\delta_i,y_i)$ on $y_i$. Thus, we expect standard survival analysis methods using the interval $(\tilde{l}(l_i,r_i,\delta_i, y_i), \tilde{r}(r_i,\delta_i,y_i)$] to yield biased estimates of parameters of models of $y_i$. \\

Consider now the problem to model $x_i$ instead of $y_i$ by standard survival analysis methods for interval censored data. Then we recognize that $x_i$ is bounded by the interval $\phi_i$ so that no pre-processing is needed and the interval can directly be used for standard survival analysis methods. Now we notice that $l_i, r_i$ are independent of $x_i$ by the non-informative censoring assumption of the three-state model in (\ref{eq::noninf3s}) and we have
\begin{align*} 
	\Pr( x_i \in (l_i, r_i] | l_i, r_i) = \Pr( x_i \in (l_i, r_i]).
\end{align*}

Thus, we expect standard survival analysis methods using the interval $(l_i, r_i]$ to yield unbiased estimates of parameters of standard survival analysis models of $x_i$.

\subsection{Proofs for Section 4.1 (Truncation Bounds)} \label{sec::proof4.1}  
We first proof the proposition
\begin{align*} 
	q( x_i |  \betaf, \sigmaf, t_i, \mathcal{D}_i,\vf_i; \tau, \lambdaf) =
	f_{x}(x_i | a(\phi_i,\delta_i,t_i)<x_i<b(\phi_i,\delta_i,t_i), \z_i, \betaf_x, \sigma_x )
\end{align*}
with truncation bounds
\begin{align} \label{eq::bounds_ab}
	(a(\phi_i,\delta_i,t_i), b(\phi_i,\delta_i,t_i))= 
	\begin{cases}
		(l_i , \infty) \quad &\text{if} \ \delta_i=1 \\
		(\max(r_i-t_i,l_i), r_i ] \quad &\text{if} \ \delta_i=2 \\
		(l_i,r_i -t_i ] \quad &\text{if} \ \delta_i=3. \\
	\end{cases}
\end{align}

We have
\begin{align}
	&q( x_i |  \betaf, \sigmaf, t_i, \mathcal{D}_i,\vf_i; \tauf, \lambdaf)  = \frac{ q(x_i, \betaf, \sigmaf, t_i, \z_i, \phi_i, \delta_i, \vf_i; \tauf, \lambdaf)
	}{ \int q(x_i, \betaf, \sigmaf, t_i, \z_i, \phi_i, \delta_i, \vf_i; \tauf, \lambdaf) d x_i} \nonumber \\
	\label{eq:app1} = &\frac{f_x(x_i | \z_i, \betaf_x, \sigma_x )
		f_t(t_i | \z_i, \betaf_t, \sigma_t )\
		\Pr( \delta_i | x_i, t_i, \phi_i)
		q(\phi_i | x_i, \vf_i )
		q(\vf_i)
	}{ \int f_x(x_i | \z_i, \betaf_x, \sigma_x )
		f_t(t_i | \z_i, \betaf_t, \sigma_t )
		\Pr( \delta_i | x_i, t_i, \phi_i)
		q(\phi_i | x_i, \vf_i )
		q(\vf_i) d x_i} \\
	\label{eq:app2} = &\frac{f_x(x_i | \z_i, \betaf_x, \sigma_x )
		\Pr( \delta_i | x_i, t_i, \phi_i)
		q(\phi_i | x_i, \vf_i )
	}{ \int f_x(x | \z_i, \betaf_x, \sigma_x )
		\Pr( \delta_i | x_i, t_i, \phi_i)
		q(\phi_i | x_i, \vf_i )
		d x_i} \\
	\label{eq:app3} = &\frac{f_x(x_i | \z_i, \betaf_x, \sigma_x )
		\mathds{1}_{\{  a(\phi_i,\delta_i,t_i) < x_i < b(\phi_i,\delta_i,t_i) \}}
	}{ \int f_x(x | \z_i, \betaf_x, \sigma_x )
		\mathds{1}_{\{a(\phi_i,\delta_i,t_i) < x_i < b(\phi_i,\delta_i,t_i) \}}
		d x_i} \\

	\label{eq:app4} = &f_{x}(x_i | a(\phi_i,\delta_i,t_i)<x_i<b(\phi_i,\delta_i,t_i), \z_i, \betaf_x, \sigma_x ).
\end{align}

Equation (\ref{eq:app1}) follows from the factorization implied by the data generating process that is visualized by the DAG in Figure 2 of the main manuscript. In equation (\ref{eq:app2}), terms independent of $x_i$ cancel, among which the density of the screening history $q(\vf_i)$. This result still holds if $q(\vf_i)$ instead depends on covariates $\z_i$, i.e., $q(\vf_i|\z_i)$, as noted in Section "Censoring mechanism" in the main manuscript. Equation (\ref{eq:app4}) follows from the definition of a truncated distribution. Therefore, it remains to be shown that $\Pr( \delta_i | x_i, t_i, \phi_i)q(\phi_i | x_i, \vf_i ) = \mathds{1}_{\{  a(\phi_i,\delta_i,t_i) < x_i < b(\phi_i,\delta_i,t_i) \}}$ holds in (\ref{eq:app3}), where the bounds $a$ and $b$ are defined in (\ref{eq::bounds_ab}). For this we note that $\Pr( \delta_i | x_i, t_i, \phi_i) q(\phi_i | x_i, \vf_i )$ (densities defined in Section "Censoring mechanism" in the main manuscript) is a product of the degenerate distributions
\begin{align} \label{eq:app5} 
	\Pr( \delta_i | x_i, t_i, \phi_i) q(\phi_i | x_i, \vf_i ) = \Pr( \delta_i | x_i, t_i, \phi_i) =
	\begin{cases}
		\mathds{1}_{\{ l_i < x_i <\infty \}} \quad &\text{if} \ \delta_i=1 \\
		\mathds{1}_{\{ l_i < x_i \le r_i <  x_i+t_i  \}} \quad &\text{if} \ \delta_i=2 \\
		\mathds{1}_{\{ l_i < x_i < x_i+t_i \le r_i < \infty \}} \quad &\text{if} \ \delta_i=3. \\
	\end{cases} 
\end{align}

The first equation holds, because $q(\phi_i | x_i, \vf_i ) $  evaluates to one whenever $\Pr( \delta_i | x_i, t_i, \phi_i)$ does. To find the bounds $(a(\phi_i,\delta_i,t_i),b(\phi_i,\delta_i,t_i))$, the inequalities in (\ref{eq:app5}) are solved for $x_i$, yielding (\ref{eq::bounds_ab}). \\

Next, we proof the proposition that
\begin{align*} 
	q( t_i |  \betaf, \sigmaf, x_i, \mathcal{D}_i,\vf_i; \tauf, \lambdaf) &= f_{t}(t_i | c(\phi_i,\delta_i,x_i)<t_i<d(\phi_i,\delta_i,x_i), \z_i, \betaf_t, \sigma_t )
\end{align*}
with truncation bounds 
\begin{align} \label{eq::bounds_cd}
	(c(\phi_i,\delta_i,x_i), d(\phi_i,\delta_i,x_i))= 
	\begin{cases}
		(0 , \infty) \quad &\text{if} \ \delta_i=1 \\
		(r_i-x_i, \infty ) \quad &\text{if} \ \delta_i=2 \\
		(0, r_i-x_i ] \quad &\text{if} \ \delta_i=3. \\
	\end{cases}
\end{align}

Here we have
\begin{align*} 
	q( t_i |  \betaf, \sigmaf, x_i, \mathcal{D}_i,\vf_i; \tauf, \lambdaf) =  &\frac{f_t(t_i | \z_i, \betaf_t, \sigma_t )
		\Pr( \delta_i | x_i, t_i, \phi_i) 
	}{ \int f_t(t | \z_i, \betaf_t, \sigma_t )
		\Pr( \delta_i | x_i, t_i, \phi_i) 
		d t_i} \\
	 = &\frac{f_t(t_i | \z_i, \betaf_t, \sigma_t )
		\mathds{1}_{\{  c(\phi_i,\delta_i,x_i) < t_i < d(\phi_i,\delta_i,x_i) \}}
	}{ \int f_t(t | \z_i, \betaf_t, \sigma_t )
		\mathds{1}_{\{c(\phi_i,\delta_i,x_i) < t_i < d(\phi_i,\delta_i,x_i) \}}
		d t_i},
\end{align*}
where we again solve the inequalities in (\ref{eq:app5}), but now for $t_i$, to find the bounds $c(\phi_i,\delta_i,x_i) $ and $d(\phi_i,\delta_i,x_i)$ given in (\ref{eq::bounds_cd}).

\subsection{Proofs for Section 4.2 (Predictive Densities)} 
We first proof the proposition that

\begin{align*} 
	q(x_i|\mathcal{D}_i,\betaf,\sigmaf) \propto g_x(x_i|\z_i,\betaf_{x},\sigma_x) [F_t(d(\phi_i,\delta_i,x_i) | \z_i, \betaf_{t},\sigma_t ) - F_t(c(\phi_i,\delta_i,x_i) | \z_i, \betaf_{t},\sigma_t )],
\end{align*}
where $g_x(x_i|\cdot) = f_x(x_i|\cdot)$ if $l_i<x_i\le r_i$ and $g_x(x_i|\cdot) = 0$ else. We have by the data generating mechanism implied by the DAG in Figure 2 of the main manuscript
\begin{align}
	q(x_i|\phi_i,\delta_i,\z_i,\betaf, \sigmaf) &= \frac{ \int q(x_i, t_i, \z_i, \phi_i, \delta_i,  \betaf, \sigmaf) dt_i
	}{ \iint q(x_i, t_i, \z_i, \phi_i, \delta_i, \betaf, \sigmaf) d x_i dt_i} \nonumber \\
    &= \frac{\int f_x(x_i | \z_i, \betaf_x, \sigma_x )
	f_t(t_i | \z_i, \betaf_t, \sigma_t )\
	\Pr( \delta_i | x_i, t_i, \phi_i)
	q(\phi_i | x_i, \vf_i )
	q(\vf_i) dt_i
}{  \iint  f_x(x_i | \z_i, \betaf_x, \sigma_x )
	f_t(t_i | \z_i, \betaf_t, \sigma_t )
	\Pr( \delta_i | x_i, t_i, \phi_i)
	q(\phi_i | x_i, \vf_i )
	q(\vf_i) d x_i dt_i} \nonumber \\
	& = \frac{ f_x(x_i | \z_i, \betaf_x, \sigma_x ) \int 
		f_t(t_i | \z_i, \betaf_t, \sigma_t )\
		\Pr( \delta_i | x_i, t_i, \phi_i) dt_i
	}{ \int  f_x(x_i | \z_i, \betaf_x, \sigma_x)  \int
		f_t(t_i | \z_i, \betaf_t, \sigma_t )
		\Pr( \delta_i | x_i, t_i, \phi_i)  d x_i dt_i}. \label{eq::condpredx}
\end{align}

Now, the numerator of (\ref{eq::condpredx}) depends on $\delta_i$; in particular
\begin{align*}
	&f_x(x_i | \z_i, \betaf_x, \sigma_x ) \int 
	f_t(t_i | \z_i, \betaf_t, \sigma_t )\
	\Pr( \delta_i | x_i, t_i, \phi_i) d t_i \\ 
	&= f_x(x_i|\z_i,\betaf_x, \sigma_x) \mathds{1}_{\{ l_i < x_i <r_i \}} \times 
	\begin{cases}
		1 \quad \ &\text{if} \ \delta_i = 1 \\
		\int \mathds{1}_{\{r_i-x_i \le t_i\}} f_t(t_i | \z_i, \betaf_t, \sigma_t ) dt_i \quad \ &\text{if} \ \delta_i = 2 \\
		\int \mathds{1}_{\{r_i-x_i>t_i\}} f_t(t_i | \z_i, \betaf_t, \sigma_t ) dt_i \quad \ &\text{if} \ \delta_i = 3 \\
	\end{cases} \\
	&= g_x(x_i|\z_i,\betaf_x, \sigma_x) \times \begin{cases}
		1 \quad \ &\text{if} \ \delta_i = 1 \\
		[ 1 - F_t(r_i-x_i| \z_i, \betaf_t, \sigma_t )] \quad \ &\text{if} \ \delta_i = 2 \\
		F_t(r_i-x_i| \z_i, \betaf_t, \sigma_t ) \quad \ &\text{if} \ \delta_i = 3. \\ 
	\end{cases}  \\
	&= g_x(x_i|\z_i,\betaf_x, \sigma_x) \times \begin{cases}
	F_t(\infty) - F_t(0) \quad \ &\text{if} \ \delta_i = 1 \\
	[ F_t(\infty) - F_t(r_i-x_i| \z_i, \betaf_t, \sigma_t )] \quad \ &\text{if} \ \delta_i = 2 \\
	F_t(r_i-x_i| \z_i, \betaf_t, \sigma_t ) - F_t(0) \quad \ &\text{if} \ \delta_i = 3. \\ 
	\end{cases} 
\end{align*}

We proceed similarly to derive the predictive density of $t_i$, i.e.
\begin{align*} 
	q(t_i|\mathcal{D}_i,\betaf,\sigmaf) \propto f_t(t_i | \z_i, \betaf_{t},\sigma_t ) [F_x(b(\phi_i,\delta_i,t_i) | \z_i, \betaf_{x},\sigma_x ) - F_x(a(\phi_i,\delta_i,t_i) | \z_i, \betaf_{x},\sigma_x )].
\end{align*}

We have
\begin{align}
	q(t_i|\phi_i,\delta_i,\z_i,\betaf, \sigmaf) = & \frac{ \int q(x_i, t_i, \z_i, \phi_i, \delta_i,  \betaf, \sigmaf) dx_i
	}{ \iint q(x_i, t_i, \z_i, \phi_i, \delta_i, \betaf, \sigmaf) d x_i dt_i} \nonumber \\
	= & \frac{ f_t(t_i | \z_i, \betaf_t, \sigma_t ) \int 
		f_x(x_i | \z_i, \betaf_x, \sigma_x )\
		\Pr( \delta_i | x_i, t_i, \phi_i) dx_i
	}{ \int  f_t(t_i | \z_i, \betaf_t, \sigma_t)  \int
		f_x(x_i | \z_i, \betaf_x, \sigma_x )
		\Pr( \delta_i | x_i, t_i, \phi_i)  d x_i dt_i}. \label{eq::condpredt}
\end{align}

The numerator in (\ref{eq::condpredt}) then simplifies as follows
\begin{align*}
	&f_t(t_i | \z_i, \betaf_t, \sigma_t ) \int 
	f_x(x_i | \z_i, \betaf_x, \sigma_x )\
	\Pr( \delta_i | x_i, t_i, \phi_i) dx_i \\ 
	&= f_t(t_i | \z_i, \betaf_t, \sigma_t ) \times 
	\begin{cases}
		\int \mathds{1}_{\{l_i < x_i \le r_i\}} f_x(x_i|\z_i,\betaf_x, \sigma_x) dx_i \quad \ &\text{if} \ \delta_i = 1 \\
		\int \mathds{1}_{\{l_i < x_i \le r_i\}} \mathds{1}_{\{ x_i > r_i-t_i\}} f_x(x_i|\z_i,\betaf_x, \sigma_x) dx_i \quad \ &\text{if} \ \delta_i = 2 \\
		\int \mathds{1}_{\{l_i < x_i \le r_i\}} \mathds{1}_{\{ x_i \le r_i-t_i\}} f_x(x_i|\z_i,\betaf_x, \sigma_x) dx_i \quad \ &\text{if} \ \delta_i = 3 \\
	\end{cases} \\
	&= f_t(t_i | \z_i, \betaf_t, \sigma_t ) \times
	\begin{cases}
		F_x(r_i|\z_i,\betaf_x, \sigma_x) - F_x(l_i|\z_i,\betaf_x, \sigma_x)  \quad \ &\text{if} \ \delta_i = 1 \\
		F_x(r_i|\z_i,\betaf_x, \sigma_x) - F_x(\max(l_i,r_i-t_i)|\z_i,\betaf_x, \sigma_x) \quad \ &\text{if} \ \delta_i = 2 \\
		F_x(r_i-t_i|\z_i,\betaf_x, \sigma_x) - F_x(l_i|\z_i,\betaf_x, \sigma_x)  \quad \ &\text{if} \ \delta_i = 3. \\
	\end{cases} 
\end{align*}

\subsection{Details on the information criteria and their implementation}
The following exposition of the information criteria WAIC and DIC follows closely the discussion in \cite{gelman_understanding_2014}. WAIC and DIC require evaluating the posterior expectation or the posterior variance of the observed data likelihood that is given by
\begin{align} \label{eq::obs_LL}
	L(\betaf, \sigmaf | \bm{\mathcal{D}}) \propto &\prod_{i: \delta_i=1} \big[ 1- F_x( l_i | \z_i, \betaf_x, \sigma_x) \big] \times \\  \nonumber
	&\prod_{i: \delta_i=2} \big[\int_{l_i}^{r_i} f_x(x_i | \z_i, \betaf_x, \sigma_x) [1 - F_t(r_i - x_i | \z_i, \betaf_t, \sigma_t)] dx_i  \big] \times \\ \nonumber
	&\prod_{i: \delta_i=3} \big[\int_{l_i}^{r_i} f_x(x_i | \z_i, \betaf_x, \sigma_x) F_t(r_i - x_i | \z_i, \betaf_t, \sigma_t) dx_i  \big]. \nonumber
\end{align}

The integration in (\ref{eq::obs_LL}) was implemented by adaptive quadrature using \texttt{R} function \texttt{integrate}. Then, 
\begin{align*}
	\text{DIC} = 2 [\log L( \hat{\betaf}, \hat{\sigmaf} | \bm{\mathcal{D}}) - 2 \E_{\betaf, \sigmaf| \bm{\mathcal{D}} }(\log L(\betaf, \sigmaf | \bm{\mathcal{D}}))],
\end{align*}
where $\hat{\betaf}$, $\hat{\sigmaf}$ are the posterior means and $\E_{\betaf, \sigmaf| \bm{\mathcal{D}} }(\log L(\betaf, \sigmaf | \bm{\mathcal{D}}))]$ denotes the posterior expectation of the likelihood, computed as
\begin{align*}
	\frac{1}{K} \sum_{k=1}^{K} \sum_{i=1}^{n} L(\betaf^{(k)}, \sigmaf^{(k)} | \mathcal{D}_i).
\end{align*}

Furthermore,
\begin{align*}
	\text{WAIC}_1 = -2 \sum_{i=1}^n \bigg[ - \log \E_{\betaf, \sigmaf| \bm{\mathcal{D}} }(L(\betaf, \sigmaf | \mathcal{D}_i)) + 2 \E_{\betaf, \sigmaf| \bm{\mathcal{D}} }(\log L(\betaf, \sigmaf | \mathcal{D}_i)) \bigg],
\end{align*}
where $\log \E_{\betaf, \sigmaf| \bm{\mathcal{D}} }(L(\betaf, \sigmaf | \mathcal{D}_i))$ is computed as
\begin{align*}
	\log \bigg[ \frac{1}{K} \sum_{k=1}^{K} L(\betaf^{(k)}, \sigmaf^{(k)} | \mathcal{D}_i)\bigg]
\end{align*}
and $\E_{\betaf, \sigmaf| \bm{\mathcal{D}} }(\log L(\betaf, \sigmaf | \mathcal{D}_i))$ is computed as
\begin{align*}
	\frac{1}{K} \sum_{k=1}^{K} \log L(\betaf^{(k)}, \sigmaf^{(k)} | \mathcal{D}_i).
\end{align*}

Finally, 
\begin{align*}
	\text{WAIC}_2 = -2 \sum_{i=1}^n \bigg[  \log \E_{\betaf, \sigmaf| \bm{\mathcal{D}} }(L(\betaf, \sigmaf | \mathcal{D}_i))  - \V_{\betaf, \sigmaf| \bm{\mathcal{D}} }( \log L(\betaf, \sigmaf |\mathcal{D}_i)) \bigg],
\end{align*}
where $\V_{\betaf, \sigmaf| \bm{\mathcal{D}} }$ is computed by the sample variance across the posterior draws
\begin{align*}
	\frac{1}{K-1} \sum_{k=1}^{K} (\log L(\betaf^{(k)}, \sigmaf^{(k)} | \mathcal{D}_i) - \bar{L})^2,
\end{align*}
where 
\begin{align*}
	\bar{L} = \frac{1}{K} \sum_{k=1}^{K} \log L(\betaf^{(k)}, \sigmaf^{(k)} | \mathcal{D}_i).
\end{align*}

\section{Supplement B: EM algorithm for the lognormal-lognormal model}\label{appn}

We now present the EM algorithm for one specific choice of error distribution, the lognormal-lognormal model, which sets $f_{\epsilon}$ and $f_{\xi}$ as standard normal. In doing so, it is useful to introduce notation for the log of transition and censoring times as follows:

\begin{align*}
\tilde{x}_i := \log x_i \quad  \quad \tilde{t}_i := \log t_i \\
\tilde{l}_i := \log l_i \quad  \quad \tilde{r}_i := \log r_i
\end{align*}

\subsection{EM algorithm}
\textbf{M-step.}
The M-step of the lognormal-lognormal model updates:

\begin{align*}
&\betaf_x^{(k+1)} = (Z_x'Z_x)^{-1}Z_x \E_{x | \mathcal{D}, \betaf_x^{(k)}, \sigma_x^{(k)} } [ \tilde{\xf} ]\\
&\betaf_t^{(k+1)} = (Z_t'Z_t)^{-1}Z_t \E_{t | \mathcal{D}, \betaf_t^{(k)}, \sigma_t^{(k)} } [ \tilde{\tf} ]
\end{align*}

and

\begin{align*}
&\sigma_x^{(k+1)} = \bigg[ \frac{1}{n} \sum_{i=1}^{n} \big( \E_{x | \mathcal{D}, \betaf_x^{(k)}, \sigma_x^{(k)} } [\tilde{x}_i^2] - 2 \E_{x | \mathcal{D}, \betaf_x^{(k)}, \sigma_x^{(k)} } [\tilde{x}_i] \z_{x,i}' \betaf_x^{(k)} + (\z_{x,i}' \betaf_x^{(k)})^2 \big) \bigg]^{\frac{1}{2}} \\
&\sigma_t^{(k+1)} = \bigg[ \frac{1}{n} \sum_{i=1}^{n} \big( \E_{t | \mathcal{D}, \betaf_t^{(k)}, \sigma_t^{(k)} } [\tilde{t}_i^2] - 2 \E_{t | \mathcal{D}, \betaf_t^{(k)}, \sigma_t^{(k)} } [\tilde{t}_i] \z_{t,i}' \betaf_t^{(k)} + (\z_{t,i}' \betaf_t^{(k)})^2 \big) \bigg]^ {\frac{1}{2}}.
\end{align*}

The E-step only involves evaluation of the normal density and the normal distribution function and is, therefore, fast.  \\

\textbf{E-step.}
If $\delta_i = 1$, the E-step consists of computing the following quantities:

\begin{align*}
&\E_{x | \mathcal{D}, \betaf_x^{(k)}, \sigma_x^{(k)} } [\tilde{x}_i] = \mu_{x,i}^{(k)}  + \sigma_x^{(k)}  \Psi \big(\frac{\tilde{l}_i-\mu_{x,i}^{(k)} }{\sigma_x^{(k)} }\big)  \\
&\E_{x | \mathcal{D}, \betaf_x^{(k)}, \sigma_x^{(k)} } [\tilde{x}_i^2] = (\mu_{x,i}^{(k)})^2  + (\sigma_x^{(k)})^2 +  \sigma_x^{(k)} (\tilde{l}_i + \mu_{x,i}^{(k)}) \Psi \big(\frac{\tilde{l}_i-\mu_{x,i}^{(k)}}{\sigma_x^{(k)}} \big)  \\
&\E_{t | \mathcal{D}, \betaf_t^{(k)}, \sigma_t^{(k)}} [\tilde{t}_i] = \mu_{t,i}^{(k)} \\
&\E_{t | \mathcal{D}, \betaf_t^{(k)}, \sigma_t^{(k)} } [\tilde{t}_i^2] = (\mu_{t,i}^{(k)})^2  + (\sigma_t^{(k)})^2, 
\end{align*}
where $\mu_{x,i}^{(k)} = \z_{x,i}' \beta_x^{(k)}$, $\mu_{t,i}^{(k)} = \z_{t,i}' \beta_t^{(k)}$, $\Psi(x) = \phi(x|0,1)/\Phi(x|0,1)$, with  $\phi(x|0,1)$ and $\Phi(x|0,1)$ the standard normal density and standard  normal cumulative distributions functions (i.e., $N(0,1)$). It can, therefore, be seen that the E-step has partly closed form when $\delta_i=1$ and is generally fast. \\

If $\delta_i = 2$, we compute:
\begin{align}
&\E_{x | \mathcal{D}, \theta^{(t)} } [\tilde{x}_i] = h_{x,i}^{-1} \int_{\tilde{l}_i}^{\tilde{r}_i} \tilde{x}_i \phi(\tilde{x}_i|\mu_{x,i}^{(k)},\sigma_x^{(k)}) [1-\Phi(\log( r_i- e^{\tilde{x}_i})|\mu_{t,i}^{(k)},\sigma_t^{(k)})] d\tilde{x}_i   \label{eq::int_xd2} \\
&\E_{x | \mathcal{D}, \theta^{(t)} } [\tilde{x}_i^2] = h_{x,i}^{-1} \int_{\tilde{l}_i}^{\tilde{r}_i} \tilde{x}_i^2 \phi(\tilde{x}_i|\mu_{x,i}^{(k)},\sigma_x^{(k)}) [1-\Phi(\log( r_i-e^{\tilde{x}_i})|\mu_{t,i}^{(k)},\sigma_t^{(k)})] d\tilde{x}_i  \label{eq::int_x2d2} \\
&\E_{t | \mathcal{D}, \theta^{(t)} } [\tilde{t}_i] =  h_{t,i}^{-1} \int_{-\infty}^{\infty} \tilde{t}_i  \phi(\tilde{t}_i|\mu_{t,i}^{(k)},\sigma_t^{(k)}) [\Phi(\tilde{r}_i|\mu_{x,i}^{(k)},\sigma_x^{(k)})-\Phi(\log(\max(l_i,r_i-e^{\tilde{t}_i}))|\mu_{x,i}^{(k)},\sigma_x^{(k)})] d\tilde{t}_i  \label{eq::int_td2} \\ 
&\E_{t | \mathcal{D}, \theta^{(t)} } [\tilde{t}^2_i] = h_{t,i}^{-1} \int_{-\infty}^{\infty} \tilde{t}_i^2  \phi(\tilde{t}_i|\mu_{t,i}^{(k)},\sigma_t^{(k)}) [\Phi(\tilde{r}_i|\mu_{x,i}^{(k)},\sigma_x^{(k)})-\Phi(\log(\max(l_i,r_i-e^{\tilde{t}_i}))|\mu_{x,i}^{(k)},\sigma_x^{(k)})] d\tilde{t}_i. \label{eq::int_t2d2} 
\end{align}

Quantities $h_{x,i}$ and $h_{t,i}$ are the normalization constants defined as

\begin{align*}
	&h_{x,i} = \int_{\tilde{l}_i}^{\tilde{r}_i} \phi(\tilde{x}_i|\mu_{x,i}^{(k)},\sigma_x^{(k)}) [1-\Phi(\log( r_i- e^{\tilde{x}_i})|\mu_{t,i}^{(k)},\sigma_t^{(k)})] d\tilde{x}_i   \\
	&h_{t,i} = \int_{-\infty}^{\infty} \tilde{t}_i  \phi(\tilde{t}_i|\mu_{t,i}^{(k)},\sigma_t^{(k)}) [\Phi(\tilde{r}_i|\mu_{x,i}^{(k)},\sigma_x^{(k)})-\Phi(\log(\max(l_i,r_i-e^{\tilde{t}_i}))|\mu_{x,i}^{(k)},\sigma_x^{(k)})] d\tilde{t}_i.
\end{align*}

The integration of all quantities above is done numerically, as described below. \\

Furthermore, if $\delta_i = 3$, we compute
\begin{align}
	&\E_{x | \mathcal{D}, \theta^{(t)} } [\tilde{x}_i] = h_{x,i}^{-1} \int_{\tilde{l}_i}^{\tilde{r}_i} \tilde{x}_i \phi(\tilde{x}_i|\mu_{x,i}^{(k)},\sigma_x^{(k)}) \Phi(\log( r_i-e^{\tilde{x}_i})|\mu_t,\sigma_t^{(k)}) d\tilde{x}_i  \label{eq::int_xd3} \\
	&\E_{x | \mathcal{D}, \theta^{(t)} } [\tilde{x}_i^2] = h_{x,i}^{-1} \int_{\tilde{l}_i}^{\tilde{r}_i} \tilde{x}_i^2 \phi(\tilde{x}_i|\mu_{x,i}^{(k)},\sigma_x^{(k)}) \Phi(\log( r_i-e^{\tilde{x}_i})|\mu_t,\sigma_t^{(k)}) d\tilde{x}_i \label{eq::int_x2d3} \\
	&\E_{t | \mathcal{D}, \theta^{(t)} } [\tilde{t}_i] =  h_{t,i}^{-1} \int_{-\infty}^{\log(r_i-l_i)} \tilde{t}_i \phi(\tilde{t}_i|\mu_t,\sigma_t^{(k)}) [\Phi( \log(r_i-e^{\tilde{t}_i})|\mu_{x,i}^{(k)},\sigma_x^{(k)})-\Phi(\tilde{l}_i|\mu_{x,i}^{(k)},\sigma_x^{(k)})] d\tilde{t}_i  \label{eq::int_td3} \\ 
	&\E_{t | \mathcal{D}, \theta^{(t)} } [\tilde{t}^2_i] = h_{t,i}^{-1} \int_{-\infty}^{\log(r_i-l_i)} \tilde{t}_i^2  \phi(\tilde{t}_i|\mu_t,\sigma_t^{(k)}) [\Phi(\log(r_i-e^{\tilde{t}_i})|\mu_{x,i}^{(k)},\sigma_x^{(k)})-\Phi(\tilde{l}_i|\mu_{x,i}^{(k)},\sigma_x^{(k)})] d\tilde{t}_i. \label{eq::int_t2d3}
\end{align}

Quantities $h_x$ and $h_t$ are again the normalization constants defined as

\begin{align*}
	&h_{x,i} = \int_{\tilde{l}_i}^{\tilde{r}_i}\phi(\tilde{x}_i|\mu_{x,i}^{(k)},\sigma_x^{(k)}) \Phi(\log( r_i-e^{\tilde{x}_i})|\mu_t,\sigma_t^{(k)}) d\tilde{x}_i   \\
	&h_{t,i} =  \int_{-\infty}^{\log(r_i-l_i)} \phi(\tilde{t}_i|\mu_t,\sigma_t^{(k)}) [\Phi(\log(r_i-e^{\tilde{t}_i})|\mu_{x,i}^{(k)},\sigma_x^{(k)})-\Phi(\tilde{l}_i|\mu_{x,i}^{(k)},\sigma_x^{(k)})] d\tilde{t}_i .
\end{align*}

\subsection{Implementation and numerical integration}
We implemented the lognormal-lognormal EM algorithm in \texttt{R} as function \texttt{em.lognorm.r} which is available as part of the \texttt{Bayes3S} package. An important technical issue is the integration in the E-step when $\delta_i \in (2, 3)$. General-purpose integration software in \texttt{R}, such as function \texttt{integrate}, is often prohibitive for integrating functions which are zero or close to zero across ranges of the integration space, such as densities (see, e.g., \texttt{?integrate}). We therefore use importance sampling for integration. For this we first note that the expectations in  (\ref{eq::int_xd2} - \ref{eq::int_t2d3}) are all of the form

\begin{align*}
	h^{-1} \int l(x) u(x) dx,
\end{align*}

where $h$ is the normalization constant, $l(x) = x$ or $x^2$ and $u(x) = \phi(x) \Phi(x)$ or $\phi(x) (1-\Phi(x))$ the remainder terms (i.e., the products of a normal density and a cumulative normal distribution function), where we suppress the dependence on the parameters in the notation for brevity. The same holds for $t$. The expectations are approximated as follows:

\begin{enumerate}
	\item Generate samples $s_i$, $i=1,...,m$, from a truncated normal proposal distribution with density $\phi(s|a \le s \le b) = v(s) / [\Phi(b)-\Phi(a)] $ where $v(s) = \phi(s)$ if $a \le s \le b$ and 0 else,  $\phi(s)$ has equal mean and standard deviation as $\phi(x)$ in $u$, and the truncation bounds $(a,b)$ are
	\begin{itemize}
		\item $(\tilde{l}, \tilde{r})$ in (\ref{eq::int_xd2} - \ref{eq::int_x2d2}) and (\ref{eq::int_xd3} - \ref{eq::int_x2d3}), 
		\item $(-\infty,\infty)$ in (\ref{eq::int_td2} - \ref{eq::int_t2d2}), and
		\item $(-\infty,\log(r_i-l_i))$ in (\ref{eq::int_td3} - \ref{eq::int_t2d3}). \\
	\end{itemize}
\item Compute importance weights:
\begin{align*}
	w(s_i) = u(s_i)/\phi(s_i|a \le s \le b).
\end{align*} \\
\item Then approximate:
\begin{align*}
	\int l(x) u(x) dx \ \hat{=} \ m^{-1} \sum_{i=1}^{m} w(s_i) l(s_i)
\end{align*}
and 
\begin{align*}
	h^{(-1)} \ \hat{=} \ m^{-1} \sum_{i=1}^{m} w(s_i).
\end{align*}
\\

\end{enumerate}

\subsection{Accuracy of integration and convergence} \label{sec::mlem_remarks}

The E- and M-steps are iterated until the change in observed data likelihood is below a threshold value (the convergence limit). However, since we have to approximate the integrals in  (\ref{eq::int_xd2} - \ref{eq::int_t2d3}), the success of the EM algorithm in finding a maximum of the observed data likelihood at a given convergence limit strongly depends on the accuracy of the integration, which in turn depends on the number of importance samples $m$. The choice of $m$ strongly impacts computation time, since it has to be carried out for all $\delta_i \in (2,3)$. Therefore, $m$ and the convergence limit have to be chosen appropriately in any application. In general, choosing smaller limits and larger $m$ seems preferable, but can be computationally expensive.

\section{Supplement C: Additional details on the simulation set-up and simulation results}

\subsection{Details on the simulation set-up}
In this Section, we provide further details on the transition time and censoring distributions of the simulation study as described in Section 5.1 of the main paper.  We generated data according to the following models: 

\begin{align*}
	\log(x_i) &= 3 + \beta_{1,x}  z_{i,1}+ \beta_{2,x} z_{i,2} + 0.2  \epsilon_{i}   \\
	\log(t_i) &= 1.2 + \beta_{1,t}  z_{i,1}+ \beta_{2,t} z_{i,2} + 0.3 \xi_{i}
\end{align*}

with $p=0$, we have $\beta_x = \beta_t = 0$, or with $p=2$, we set $\beta_x=\beta_t=0.5$. A plot of the resulting marginal transition time densities $f_x(x|\betaf_x, \sigma_x)$ and $f_t(t|\betaf_t, \sigma_t)$ is given in Figure \ref{fig::times} and the log-times are given in Figure \ref{fig::logtimes}.

\begin{figure}
\caption{Transition time distributions of $x_i$ and $t_i$ based on a sample of size $n=10^6$ .}
\label{fig::times}
\end{figure}
\begin{figure}
\caption{Log transition time distributions of $x_i$ and $t_i$ based on a sample of size $n=10^6$.}
\label{fig::logtimes}
\end{figure}

Time $x_i$ had greater median (20.1)  than $t_i$ (3.3) under $p=0$; similarly, for $p=2$ the medians were: 25.8 and 4.3, respectively. However, the introduction of covariates translates into stronger heterogeneity resulting in distributions exhibiting stronger left skewed forms and longer right tails. Note that with $p=0$ the two distributions are lognormal, whereas with $p=2$ the distributions are a mixture. The times were marginally independent under $p=0$ but strongly correlated under $p=2$ (Pearson $r=.799$). This (marginal) correlation emerged, because the same covariates $\z_i$ were included in both models. Conditionally on the covariates, the times were independent. \\

The simulation was designed to mimic the practical setting where the disease has long onset time $x_i$ but progresses more rapidly once state 2 (e.g., CIN-2) has occurred. A similar finding is reported for the POBASCAM data in the main manuscript, Section 6. We did not intend to replicate the same distributions from the application, but we simulate a similar disease process in terms of a strong difference in progression speeds.  \\

We introduced the censoring times $(l_i, r_i)$ by a visiting process that first simulated $v_{i,1} \sim \text{uniform}(c_{min},c_{max})$, where $c_{min}<c_{max}$, and then recursively:
\begin{align}
	v_{i,t+1} \sim \text{uniform}(v_t+c_{min},v_t+c_{max})
\end{align}
until $v_{i,t+1} > v_{i,rc}$ where $v_{i,rc} \sim \text{exp}(\theta^{-1})$ was the transition time of right censoring of $x_i$ (end of follow-up). In particular, the bounds $(c_{min},c_{max})$ denote, respectively, the minimum and maximum time that elapsed between recurring visits. For example, subjects may be required not to be screened more frequently than every year for a particular disease, $c_{min} = 1$, but are expected or known to attend screening at least every seven years, $c_{max} = 7$. However, in the simulation study we needed to obtain approximately equal probabilities of censoring states $\delta_i \in \{1,2,3\}$ while the distributions of the underlying transition times of $x_i$ and $t_i$ differed across the number of covariates $p$. We tuned the visiting process parameters $c_{min},c_{max},\theta$ to achieve this goal. \\

\begin{table*} [h]
	\centering
	\caption{Parameters of the screening times process and resulting distributions of $\delta_i$ \\($^*$estimated by Monte Carlo with $10^6$ replications). Parameters $c_{min}$ and $c_{max}$ indicate the minimum and maximum time between screening moments and $\theta$ indicates the mean time to right censoring of $x_i$}
	\label{tab::simcens}
	\begin{tabular}{@{}l c c l c c @{}}
		& \multicolumn{2}{c}{Medium censoring} && \multicolumn{2}{c}{Strong censoring} \\
		\cline{2-3}
		\cline{5-6}
		& $p=0$ & $p=2$ &&  $p=0$ & $p=2$  \\ 
		\hline
		$c_{min}$ & 1 & 1 && 1 & 1 \\  
		$c_{max}$ & 8 & 8.7 && 7 & 6.5 \\  
		$\theta$ & 40 & 56 && 20 & 26.1 \\
		$\Pr(\delta_i = 1 | c_{min},c_{max},\theta )^*$ & 0.368 & 0.370 && 0.600 & 0.601 \\
		$\Pr(\delta_i = 2 | c_{min},c_{max},\theta)^*$ & 0.426 & 0.426 && 0.298 & 0.296 \\
		$\Pr(\delta_i = 3 | c_{min},c_{max},\theta)^*$ & 0.206 & 0.204 && 0.103 & 0.103 \\	
		\hline
	\end{tabular}
\end{table*}

The final choice of parameters is reproduced \ref{tab::simcens} from the main manuscript. As can be seen, the parameters differ between $p=0$ and $p=2$ for the medium and strong censoring conditions, respectively. This was done to achieve approximate equivalence in the probability mass of $\delta_i$ between the levels of $p$. The main difference introduced between the medium and strong censoring conditions is due to parameter $\theta$ which determines the mean time to right censoring. This time is substantially shorter for the strong censoring condition (i.e., 20/26.1 as opposed to 40/56) so that the probability of $\delta_i = 1$ is strongly elevated. \\ 

It is useful to note that the probabilities given in Table \ref{tab::simcens} are large-sample statistics that were determined in a sample of size $10^6$. In samples of size $n=1000$ or $n=2000$, as used in the simulation study, the observed proportion of cases with $\delta_i=j$ had a sampling variance. This is illustrated by Figure \ref{fig::delta_variance}. 

\begin{figure}
\caption{Histograms of the proportion of cases in categories of $\delta_i$ over 1,000 replicated data sets illustrating sampling variance of the sample proportions of the censoring events.}
\label{fig::delta_variance}
\end{figure}

\subsection{Setting the jumping distribution proposal variance}
In the main manuscript, we study convergence of the MCMC algorithm, where the Gibbs algorithm employed a Metropolis step with an independent multivariate normal jumping distribution. It is well-known that the proposal variance of the jumping distribution impacts the acceptance rate of the Metropolis sampler and thus its speed of convergence. It was, therefore, essential to have approximately equal acceptance rates across experimental conditions to have comparable convergence rates. In general, both, too large and too small acceptance rates are problematic, where for multivariate posterior distributions the 'optimal' proposal variance invokes an acceptance rate of approximately 23\% \citep[][p. 296]{gelman_bayesian_2013}. \\

We tuned the proposal variance to every data set separately. For this we employed an automated heuristic search as part of the simulation code, described as follows:

\begin{enumerate}
	\item Set target interval $J = [J_1,J_2]$ for the acceptance rate $w$ to $J=[0.22,0.24]$.
	\item Run the Gibbs sampler with $m=3$ chains for an initial $D=5,000$ draws with proposal standard deviation $s=0.1$ in the Metropolis step.
	\item Determine $w$ as the proportion of accepted draws from the proposal distribution across $m$.
	\item If $w \in J$ go to step 5. Else set: $s_{new} : = s ( 1-(J_1 - w)/J_1 )$ if $s \ge J_1$ and  $s_{new} : = s ( 1+(J_2 - w)/J_2 )$ if $s < J_2$.
	\item Run the Gibbs sampler for another $D = 3,000$ draws and calculate $w$ as the proportion of accepted draws from the proposal distribution across $m$ and the last $D$ draws.
	\item If $w \in J$, go to step 7. Else set $s_{new}$ as described in step 4 and go to step 5.
	\item Repeat steps 5 and 6 for $D=6,000$ and $D=12,000$. Stop if $D=12,000$ and $w \in J$. 
\end{enumerate}

We emphasize that this is a heuristic search algorithm that mimics manual tuning of the proposal variance with the goal to achieve an optimal sampling rate. The reason this step was automated is that it is run as part of the simulation study, which required tuning the acceptance rate to $8 \times 100$ data sets under different experimental conditions.  Alternatives to manual tuning of the proposal variance include adaptive Metropolis Hastings algorithms which could be considered in future work. However, in initial tests with adaptive algorithms we found degenerate MCMC chains when used within the Gibbs sampler of the 3-state model. We therefore gave preference to manual tuning of the proposal variance which proved to reliably achieve convergence. \\

After the search is completed the final proposal variance $s^2$ is used in a new run of the MCMC sampler (that is, discarding any draws from the search algorithm). The proposal variance thus is kept constant at the determined value across the whole run of the MCMC sampler. We compare the acceptance rates obtained from the various experimental conditions in the simulations in Figure \ref{fig::acc_rate}. As can be seen, the procedure was successful in obtaining acceptance rates close to 23\% in all conditions. 

\begin{figure}
\caption{ Metropolis sampling acceptance probability over the simulation replications (proportion of accepted draws from the proposal distribution after discarding half of samples for burn-in). Dashed line indicates 'optimal' acceptance rate of 23\%.}
\label{fig::acc_rate}
\end{figure}

\subsection{Implementation details for EM, MSM, and NPMLE} \label{sec::simimpl}

Besides the Bayesian MCMC algorithm, we modeled the simulated data also using the lognormal-lognormal maximum likelihood EM algorithm (Supplemental Material, Section B), the  \texttt{msm} package \citep{jackson_multi-state_2011} and a non-parametric maximum likelihood estimator \citep[NPMLE;][]{turnbull_empirical_1976}. \\

We implemented the EM algorithm for the lognormal-lognormal model in \texttt{R} as function \texttt{em.lognorm.r}; see Supplemental Material, Section B. For an initial, approximate run, we first set $m=100$ (the number of samples used in importance sampling) and considered the EM algorithm as converged when the change in observed data likelihood came below a value of $10^{-4}$. The approximate solution can typically be obtained fast but may be inaccurate. Therefore, we refined the approximate solution by more precise ML-EM runs with higher integration accuracy. For this, we set $m=5,000$ and started a new EM run at the values estimated by the initial approximate run. The convergence limit for this run was set to $10^{-5}$. \\

For \texttt{msm}, data were first transformed into the format required by the package. We then defined the following transition matrix (\texttt{qmatrix}):
\lstinputlisting[language=R]{Q.txt}
which lets \texttt{msm} estimate a 3-state model with allowed transitions 1 to 2 and 2 to 3; all other transitions (zero entries) are impossible. The specific values in the transition matrix are starting values. We passed the \texttt{qmatrix} and data to \texttt{msm} and left all other values at their defaults. To estimate the marginal cumulative incidence functions (CIF), we proceeded as follows. First we estimated the conditional transition probabilities using function \texttt{pmatrix.msm} for all observed covariate data ($\z$). Let $P_i$ denote the estimated transition probability matrix for individual $i$ with $\z_i$ at time \texttt{t1} with entry $p_{ijk}$, then the CIF percentile of $x_i$ is obtained by the sum of matrix entries $p_{i\cdot}=p_{i12} + p_{i13}$, the percentile of $t_i$ is obtained by $p_{i13}$ and the percentile of $y_i$ is obtained by $p_{i23}$. The corresponding percentile of the marginal CIF is finally obtained by the average of all $p_{i\cdot}$ and $p_{i23}$ across $i=1,...,n$. \\

Furthermore, we used function \texttt{icfit} from package \texttt{interval} \citep{fay_exact_2010} to obtain a NPMLE of the marginal CIFs of $x$ and $y$, see Supplemental Material Section \ref{sec::infcens}. In particular, NPMLE was implemented for $x$ by treating $\delta=2$ and $3$ as the same event (since $\delta=3$ occurring in an interval suggests that state 2 occurred in the same interval). This was achieved by simply passing $l$ to argument \texttt{L} and $r$ to \texttt{R} using the same coding applied in the main manuscript. For $y$, we applied the recoding scheme described in  Section \ref{sec::infcens}, before passing $\tilde{l}$ and $\tilde{r}$ to the arguments \texttt{L} and \texttt{R}. 

\subsection{MCMC convergence rates by parameters}
We evaluated the rate of convergence of the MCMC chains for each of the parameters of $f_x$ and $f_t$ in the 3-state model, i.e., the intercepts $\beta_0$, the slope coefficients $\beta_1$ and $\beta_2$, and the scale coefficients $\sigmaf$. We did so for the first $10^5$ MCMC draws using the Gelman Rubin convergence statistic (Figure \ref{fig::sim_convergence}). We found that the chains of the parameters of $f_x$ quickly converged whereas the chains of the parameters of $f_t$ needed substantially more draws, in particular in the $p=2$ condition. The MCMC chain of parameter $\sigma_t$ had the slowest rate of convergence. We also found that stronger censoring decreased the rate of convergence for all parameters. We note again that eventually all MCMC chains converged, as described above and shown in the main manuscript.

\begin{figure}
\caption{Gelman-Rubin convergence diagnostic $R$ across the first $10^5$ MCMC draws (medians of $R$ over 500 replications per condition). The speed of descent indicates the rate of convergence for the parameters of $f_x$ and $f_t$ ($\beta_1$ and $\beta_2$ reported jointly as 'betas').}
\label{fig::sim_convergence}
\end{figure}

\subsection{Frequentist coverage probability of the Bayesian credible intervals}

We studied the frequentist coverage properties of Bayesian credible intervals as follows. In the simulation, we computed the $95\%$ credible intervals for each parameter by determining the 0.025 and 0.975 quantiles of the posterior distribution (i.e., equal tailed credible intervals). We then determined per parameter and simulation condition the proportion out of all 500 credible intervals that included the true parameter (Figure \ref{fig::coverage}). \\

\begin{figure}
\caption{ True parameter coverage rate of 95\% credible intervals over 500 replicated data sets. Dashed line indicates nominal coverage rate (0.95). Dotted line at 0.908 indicates the minimum value at which a normal approximation 95\% confidence interval for the estimated coverage rate still includes the nominal (true) coverage rate of 0.95 after adjusting for multiple testing using Bonferroni correction. }
\label{fig::coverage}
\end{figure}

Figure \ref{fig::coverage} displays a horizontal dotted line per sub-plot which indicates the minimum value at which a normal approximation 95\% confidence interval for the estimated coverage rate still includes the nominal (true) coverage rate of 0.95. This cut-off value was determined by finding the zero point of
\begin{align*}
	\bigg(s + \Phi^{-1} \bigg (1- \frac{0.05}{48 \times 2}\bigg ) \sqrt{\frac{s(1-s)}{K}} - 0.95 \bigg)^2 \quad ; \quad s \in (0,1)
\end{align*}
across $s$, with $s$ the lower bound of interest, $\Phi^{-1}$ the quantile function of cumulative normal distribution and $K=500$ the number of replications in the simulation. The ratio $\frac{0.05}{48 \times 2}$ carries out a Bonferroni correction for multiple testing (48 credible intervals are evaluated). We determined $s$ at approximately $0.908$ by \texttt{optimize} in \texttt{R}. It can be seen that in almost all of the conditions the estimated proportion was greater than this lower bound, indicating that the coverage probability of the credible intervals was at nominal level for most parameters and under all conditions. However, for $\sigma_t$ we found slight under-coverage in the condition with $n=2000$ and $p=2$. In addition, coverage of this parameter seemed to be slightly lower across all conditions in the larger sample size setting ($n=2000$).

\subsection{Simulation results EM compared to Bayes}
We compared the relative error of the parameter estimates of the Bayesian Gibbs sampler to those of ML-EM (Figures \ref{fig::mcmcvsem_x}-\ref{fig::mcmcvsem_t_out}). Relative error was defined as $[\hat{\kappa}-\kappa]/\kappa$  where $\kappa$ is the true value of the parameter and $\hat{\kappa}$ a posterior median estimate or a ML estimate. Figure \ref{fig::mcmcvsem_x} shows the relative errors for the parameters of the model of $x_i$. As can be seen, EM and Bayes had largely the same error sizes. In particular, the error for all parameters was zero in central tendency (approximately unbiased estimators) under all conditions. The same result was found for the parameters of the model of $t_i$ (Figure \ref{fig::mcmcvsem_t_out}). In addition, we estimated the root of the mean squared error (root MSE)  for all parameters, defined as $\sqrt{\frac{1}{K}\sum_{k=1}^K(\hat{\kappa}_k-\kappa)^2}$. As can be seen from Figures \ref{fig::mcmcvsem_x_rmse} and \ref{fig::mcmcvsem_t_rmse}, the estimators had very similar root MSE sizes.

\begin{figure}
\caption{ Relative error of maximum likelihood EM compared to Bayesian estimation over 500 replications for the parameters of the transition time $x_i$.  }
\label{fig::mcmcvsem_x}
\end{figure}

\begin{figure}
\caption{ Relative error of maximum likelihood EM compared to Bayesian estimation over 500 replications for the parameters of the transition time $t_i$.}
\label{fig::mcmcvsem_t_out}
\end{figure}

\begin{figure}
\caption{ Root mean squared error of maximum likelihood EM compared to Bayesian estimation over 500 replications for the parameters of the transition time $x_i$.}
\label{fig::mcmcvsem_x_rmse}
\end{figure}

\begin{figure}
\caption{ Root mean squared error of maximum likelihood EM compared to Bayesian estimation over 500 replications for the parameters of the transition time $t_i$.}
\label{fig::mcmcvsem_t_rmse}
\end{figure}

\subsection{Comparison of cumulative incidence functions  models}

In this section, we give the cumulative incidence functions compared across estimation methods for all simulation scenarios. For the simulation setting, see the labeling on top of the figures.

\begin{figure}
\caption{Marginal cumulative incidence functions (CIFs) compared for Bayesian estimation, EM, Markov model, and NPMLE. CIFs are point-wise averaged over 500 replications. For simulation setting, see information on top of the plot. }
\end{figure}
\begin{figure}
\caption{Marginal cumulative incidence functions (CIFs) compared for Bayesian estimation, EM, Markov model, and NPMLE. CIFs are point-wise averaged over 500 replications. For simulation setting, see information on top of the plot. }
\end{figure}
\begin{figure}
\caption{Marginal cumulative incidence functions (CIFs) compared for Bayesian estimation, EM, Markov model, and NPMLE. CIFs are point-wise averaged over 500 replications. For simulation setting, see information on top of the plot. }
\end{figure}
\begin{figure}
\caption{Marginal cumulative incidence functions (CIFs) compared for Bayesian estimation, EM, Markov model, and NPMLE. CIFs are point-wise averaged over 500 replications. For simulation setting, see information on top of the plot. }
\end{figure}
\begin{figure}
\caption{Marginal cumulative incidence functions (CIFs) compared for Bayesian estimation, EM, Markov model, and NPMLE. CIFs are point-wise averaged over 500 replications. For simulation setting, see information on top of the plot. }
\end{figure}
\begin{figure}
\caption{Marginal cumulative incidence functions (CIFs) compared for Bayesian estimation, EM, Markov model, and NPMLE. CIFs are point-wise averaged over 500 replications. For simulation setting, see information on top of the plot. }
\end{figure}
\begin{figure}
\caption{Marginal cumulative incidence functions (CIFs) compared for Bayesian estimation, EM, Markov model, and NPMLE. CIFs are point-wise averaged over 500 replications. For simulation setting, see information on top of the plot. }
\end{figure}
\begin{figure}
\caption{Marginal cumulative incidence functions (CIFs) compared for Bayesian estimation, EM, Markov model, and NPMLE. CIFs are point-wise averaged over 500 replications. For simulation setting, see information on top of the plot. }
\end{figure}

\subsection{Estimation performance under misspecification}

In this section, we study the performance of the information criteria WAIC-1, WAIC-2, and DIC in selecting the correct distributions from a number of candidate models, where only one model is correctly specified. Subsequently, we study the bias in the model parameters and CIFs arising from distributional misspecification. For this we also used the simulation setting described in the main manuscript with lognormally distributed $x_i$ and lognormally distributed $t_i$ (in the following called the 'lognormal-lognormal' model). We specified and estimated the following models:
\begin{enumerate}
	\item lognormal - lognormal
	\item lognormal - Weibull
	\item lognormal - exponential
	\item Weibull - lognormal
	\item Weibull - Weibull
\end{enumerate}
We fit these models to the true lognormal-lognormal data from the simulation study. Models 2 to 5 misspecified the true model either with respect to $t_i$ (models 2 and 3),  $x_i$ (model 4), or both transition times (model 5). Subsequently, we ran the model until convergence and selected the model yielding smallest WAIC-1, WAIC-2, or DIC . We replicated this experiment 100 times on randomly simulated data sets for the setting with $p = 2$ covariates (Table \ref{tab::ICcomp}). Our results indicated that the true distribution of $x_i$ was very reliably identified, regardless of whether $t_i$ was correctly (model 4) or incorrectly specified (model 5). The results applied to all information criteria. The true distribution of $t_i$ was correctly selected in 66\% to 69\% of replicated runs correctly (WAIC-1), unless the sample was small and censoring was strong (50\%). WAIC-2 and DIC performed similarly. For $t_i$ we considered two alternative misspecified models, Weibull (model 2) and exponential (model 3). Among the two, the information criteria reliably selected the more complex misspecified model in 31\% to 34\% of all runs (Weibull / model 2) over the simpler misspecified model in none of the runs (exponential / model 3). 

\begin{table*} [h]
	\centering
	\caption{Performance of the information criteria WAIC-1, WAIC-2, and DIC in identifying the correct model among various misspecified models (1. lognormal - lognormal, 2. lognormal - Weibull, 3. lognormal - exponential, 4. Weibull - lognormal, 5. Weibull - Weibull). Setting with $p=2$ covariates and either $n = 1000$ or 2000 and medium or strong censoring. Numbers give the proportion out of all runs that a model is selected (in percent \% of 100 replicated runs). Model 1 is the correctly specified model (lognormal-lognormal). \\}
		\label{tab::ICcomp}
		\begin{tabular}{@{}l c c c c c l c c c c c l c c c c c @{}}
			& \multicolumn{5}{c}{WAIC-1} && \multicolumn{5}{c}{WAIC-2} && \multicolumn{5}{c}{DIC} \\
			\cline{2-6}
			\cline{8-12}
			\cline{14-18}
			Model & 1 & 2 & 3 & 4 & 5 && 1 & 2 & 3 & 4 & 5 && 1 & 2 & 3 & 4 & 5 \\ 
			\hline
			n = 1000, medium & 66 & 34 & 0 & 0 & 0 && 65 & 35 & 0 & 0 & 0 && 64 & 36 & 0 & 0 & 0 \\
			n = 1000, strong & 50 & 50 & 0 & 0 & 0 && 49 & 51 & 0 & 0 & 0 && 51 & 49 & 0 & 0 & 0 \\
			n = 2000, medium & 69 & 31 & 0 & 0 & 0 && 71 & 29 & 0 & 0 & 0 && 69 & 31 & 0 & 0 & 0 \\
			n = 2000, strong & 69 & 31 & 0 & 0 & 0 && 68 & 32 & 0 & 0 & 0 && 68 & 32 & 0 & 0 & 0 \\
			\hline
		\end{tabular}
\end{table*}

As the information criteria occasionally selected the wrong distribution for $t_i$, we also assessed parameter bias and accuracy of the estimated CIFs under misspecification. We found that the coefficients $\beta_1$ and $\beta_2$ were estimated with negligible bias (Figures \ref{fig::modcomp_par1}-\ref{fig::modcomp_par2}), unless the $t_i$ model was misspecified as exponential. The intercept $\beta_0$ and scale parameter $\sigma$ had larger error. However, it is important to note that these parameters are hard to compare across distributions. It is more important to assess, whether the cumulative incidence functions were comparable across misspecified models (Figure \ref{fig::modcomp_cif}). We found that the bias emerging from misspecification in the CIFs was, indeed, negligible, unless the transition was misspecified as exponentially distributed.

\begin{figure}
\caption{Mean relative error by misspecified models over 100 replicated data sets (500 sets for the lognormal-lognormal model). Parameters of the $x$-model shown. }
\label{fig::modcomp_par1}
\end{figure}
\begin{figure}
\caption{Mean relative error by misspecified models over 100 replicated data sets (500 sets for the lognormal-lognormal model). Parameters of the $t$-model shown. Mean relative error of the lognormal-exponential model off the scale.}
\label{fig::modcomp_par2}
\end{figure}
\begin{figure}
\caption{Cumulative Incidence functions (CIFs) averaged point-wise over 100 replicated data sets. }
\label{fig::modcomp_cif}
\end{figure}

\section{Supplement D: Additional information on the POBASCAM case study}
\subsection{MCMC chains of lognormal-loglogistic model}
To assess convergence of each models, we inspected the MCMC chains of the parameters of the models of $\log x_i$ and $\log t_i$. Figures \ref{fig::mcmc_chain.x} and \ref{fig::mcmc_chain.t} depict the MCMC chain of the lognormal-loglogistic model including the warm-up period. The figures depict the first five of fifteen MCMC chains, for clarity. The other chains looked similar.
\begin{figure}
\caption{MCMC chains 1 through 5 (out of 15) of the lognormal-loglogistic model for the parameters of $\log x_i$. The chains were thinned using a step length of 200 replications. }
\label{fig::mcmc_chain.x}
\end{figure}

\begin{figure}
\caption{MCMC chains 1 through 5 (out of 15) of the lognormal-loglogistic model for the parameters of $\log t_i$. The chains were thinned using a step length of 200 replications.}
\label{fig::mcmc_chain.t}
\end{figure}
 Convergence of all chains looked appropriate, as judged by a clear trend towards a single region of the posterior parameter space and good mixing of the MCMC chains. The thinned chains exhibited auto-correlation which is a well-known characteristic of MCMC chains when using data augmentation; see discussion Section of the main manuscript.

\subsection{Implementation notes for NPMLE and the Markov model}
The NPMLE for the CIF of $x_i$ and $y_i$ was implemented for the POBASCAM data using \texttt{R} package \texttt{interval} \citep{fay_exact_2010}, as in the simulation study (Section \ref{sec::simimpl}). The interval bounds were recoded before estimating $y_i$, as explained in Section \ref{sec::simimpl}. \\

A three-state Markov model was also implemented for the POBASCAM data using \texttt{R} package \texttt{msm} \citep{jackson_multi-state_2011}. The implementation was equivalent to that of the simulation study, Section  \ref{sec::simimpl}. 

\subsection{Lognormal-lognormal model fit by ML-EM }

We estimated the lognormal-lognormal model by the maximum likelihood expectation-maximization (ML-EM) algorithm described in Section B of the Supplemental Material.

\subsubsection{Implementation}
 We ran the ML-EM algorithm 100 times with random starting values. Each run first executed an approximate estimation by running the ML-EM algorithm with a low number of importance samples in the numerical integration ($m=100$) and a convergence limit of $10^{-4}$. By convergence limit we mean that the ML-EM algorithm was considered as converged when the change in the observed data likelihood came below the limit. The approximate solution can typically be obtained fast but may be inaccurate. Therefore, we refined the approximate solution by more precise ML-EM runs with higher integration accuracy, increasing $m$ to $2 \times 10^4$ importance samples and decreasing the convergence limit to $10^{-5}$. This run used the estimates from the approximate solution as starting values.

\subsubsection{Results}
The repeated runs of the ML-EM algorithm converged at observed data log-likelihood values between -977.874 and -977.698. This result indicated a flat and/or multi-modal likelihood function, because we used a convergence limit of $10^{-5}$ and the differences in the obtained log-likelihood values after convergence had a range of values of 0.3. By this we mean that the algorithm either converged too early in a flat region before it reached the optimum or the algorithm found a local maximum. We, therefore, inspected the parameter estimates of the converged EM runs (Figure \ref{fig::MLEs}). The estimates of the parameters of the model of $x_i$ were stable across runs, but parameter estimates of the model of $t_i$ showed strong variation. In particular, $\sigma_t$ (minimum: 3.36, maximum: 5.76) and the regression coefficient of HPV-16 (-3.99, -2.43) varied strongly. However, visual inspection of the likelihood values against the parameter values indicated that the maximum log-likelihood value was approximately -977.698. We, therefore, used the estimates from this run as the maximum likelihood estimates (MLE). \\
\begin{table*}
	\centering
	\caption{Maximum likelihood estimates obtained by EM for the lognormal-lognormal model (log-likelihood = -977.698) compared to posterior median estimates from the Bayesian lognormal-lognormal model ($x_i$: time to CIN-2; $t_i$: time from CIN-2 to CIN-3+).}
	\label{tab::MLEs}
	\begin{tabular}{@{} l r r r r r @{}}
		& \multicolumn{2}{c}{MLE} && \multicolumn{2}{c}{Bayes} \\
		\cline{2-3}
		\cline{5-6}
		Parameter & $x_i$ & $t_i$ && $x_i$ & $t_i$ \\
		\hline		
		Intercept          & 8.89   &  0.64  &&  7.94 &  0.06  \\
		Age (standardized) & 0.95   & -0.45  &&  0.84 & -0.34  \\
		HPV-16*            & -3.76  & -3.03  && -3.04 & -2.01  \\
		HPV-18*            & -2.25  & -1.53  && -1.38 & -0.53  \\
		HPV-31*            & -2.05  & -1.08  && -1.35 & -0.33  \\
		sigma              & 5.45   &  4.29  &&  4.96 &  4.07  \\
		\hline
		\multicolumn{6}{l}{ * Reference: other HPV sub-type} \\
	\end{tabular}
\end{table*}

Table \ref{tab::MLEs} compares the MLE with the posterior median estimates from the Bayesian lognormal-lognormal model. It is important to stress that the Bayesian estimates in Table \ref{tab::MLEs} are for the lognormal-lognormal model and they differ slightly from those in the main manuscript, because the model presented there is the lognormal-loglogistic model which fitted the data best. We found that the MLE and Bayesian estimates were comparable in magnitude and signs. In general, the  Bayesian estimates were slightly smaller, in absolute value than the MLE, which could be a result of the applied regularization. Figure \ref{fig::res_ppdplots_modcomp} compares the marginal cumulative incidence functions (CIF) of the Bayesian lognormal-lognormal model and the ML-EM fit to the Bayesian lognormal-logistic model (shown in the main manuscript). We find that the CIFs of the three models were almost identical on $t_i$ and $y_i$. The ML-EM CIF of $x_i$ had slightly stronger curvature than those of the Bayesian models.

\begin{figure}
\caption{Parameter estimates by the EM algorithm (convergence limit $10^{-5}
\label{fig::MLEs}
\end{figure}

\begin{figure}
\caption{Posterior median of predictive cumulative incidence functions (CIF) for the Bayesian lognormal-loglogistic model (best performance acc. to WAIC-1), Bayesian lognormal-lognormal model, and lognormal-lognormal model fit by ML-EM. 95\% credible intervals of the lognormal-loglogistic model shown as shaded region.  }
\label{fig::res_ppdplots_modcomp}
\end{figure}

\subsection{Conditional independence sensitivity analysis} \label{sec::application_sensitivity}

We assessed the sensitivity of our findings to the conditional independence assumption for the lognormal-loglogistic model. In particular, we specified a grid of parameters 
\begin{align*}
	(\beta_{x,w},\beta_{t,w}) = \{ (-3,3),(-2,2),(-1,1),(0,0),(1,1),(2,2),(3,3) \}
\end{align*}
in model 
\begin{align} 
	\log x_i = \z_{x,i}' \betaf_x + \beta_{x,w} w_i + \sigma_x \epsilon_{i} \label{eq::model_eq1_sens1} \\
	\log t_i = \z_{t,i}' \betaf_t + \beta_{t,w} w_i + \sigma_t \xi_{i} \label{eq::model_eq2_sens1},
\end{align}
where $(0,0)$ denotes the conditional independence model. We included both positive and negative $\beta_{x,w}$ to allow for an induced positive or negative partial association of $x_i$ and $t_i$. Condition $(1,1)$ generates, for example, a moderately strong unobserved confounding variable $w_i \sim N(0,1)$ for all $i$: individuals with $w_i$ one and two standard deviations above the mean have $\exp(1) = 2.71$ and $\exp(2) = 7.38$ times greater transitions times. Condition $(3,3)$ denotes a very strong effect, similar to that of HPV-16 in the conditional independence model for $\log x_i$.  All considered deviations from conditional independence caused a stronger curvature in the predictive CIF, especially that of $t_i$ (Figure \ref{fig::sens_ppd}). However, the conclusion of a rapid transition from CIN-2 to CIN-3+ by over half of the population was not changed, and the overall findings on the transition from baseline to CIN-2 and CIN-3+ remained very similar. The overlaid 95\% credible bands of all sensitivity conditions in Figure \ref{fig::sens_ppd} quantify the uncertainty added by relaxing conditional independence. \\

The sign and significance of the regression coefficients did not change as compared to conditional independence (Figure \ref{fig::sens_cind_par}). The posterior distribution of the strongest effect (HPV-16) was sensitive to $(\beta_{x,w},\beta_{t,w})$ but it moved away from zero (i.e., the estimated effect was stronger) for larger absolute values of  $(\beta_{x,w},\beta_{t,w})$, while uncertainty increased. This held for both positive and negative $\beta_{x,w}$. In addition, scale parameters $\sigmaf$ and the intercept were more extreme for larger absolute values of  $(\beta_{x,w},\beta_{t,w})$, resulting in the stronger curvature observed in Figure \ref{fig::sens_ppd}. 

\begin{figure}
\caption{Posterior median of predictive cumulative incidence functions for different assumptions on $(\beta_{x,w}
\label{fig::sens_ppd}
\end{figure}

\begin{figure}
\caption{Posterior parameter sensitivity to the conditional independence assumption. Black lines equivalent to main analysis model (lognormal-loglogistic).}
\label{fig::sens_cind_par}
\end{figure}

 

 

    

\end{document}