\section{Introduction}

Medicine is in a period of transition. An ever-increasing amount of information is available on patients ranging from genetic and epigenetic profiles enabled by next-generation sequencing to moment-to-moment data collected by physical activity monitors. With this wealth of information comes the opportunity to provide more targeted healthcare including, for example, prediction of pre-clinical atherosclerosis \citep{McGeachie2009}, individualized cancer screening (Saini, van Hees, and Vijan, 2014)\nocite{Saini2014}, sub-typing of scleroderma (Schulam, Wigley, and Saria, 2015)\nocite{Schulam2015},  and personalized cancer treatment \citep{Hayden2009}. In order to fully realize the promise of patient-focused medicine, principled statistical methods are needed that integrate data from a variety of sources in order to provide physicians and patients with relevant syntheses to inform their decision-making. These methods must also accommodate limitations common to data generated in an observational setting including measurement error and informative missing data patterns.

An excellent example of this challenge is low-risk prostate cancer diagnosis. Tumor lethality is an aspect of an individual's health state that is not directly observable but is manifest in multiple types of measurements including biomarkers, histology of biopsied tissue, genetic markers, and family history of the disease. Individualized predictions of the latent disease state are critical to guide treatment decisions. If the tumor is potentially lethal, immediate treatment (including surgery or radiation) can be life-saving. Yet, some tumors are indolent and not life-threatening. In this case, treatment is not recommended due to the risk of lasting side effects including urinary incontinence and erectile dysfunction \citep{Chou2011b}.


Active surveillance (AS) offers an alternative to early treatment for individuals with lower risk disease  \citep{DallEra2012}. Though AS regimes vary, the approach generally entails regular biopsies (e.g., annually) with intervention recommended upon detection of higher risk histological features, as determined by the Gleason grading system \citep{Gleason1992}. Biopsies with a Gleason score of 6 (the minimum for prostate cancer diagnosis) indicate low risk disease while a subsequent Gleason score of 7 or above is considered ``grade reclassification'' \citep{Tosoian2015}; treatment is recommended once grade reclassification is observed. Prostate-specific antigen (PSA), a blood serum biomarker of inflammation in the prostate, is also routinely measured and may be used as the basis for a biopsy recommendation. 

The success of AS programs depends on clinicians' ability to identify tumors with metastatic potential with sufficient time for curative intervention to be effective. Yet, biopsies used to characterize tumors typically sample less than one percent of the prostate tissue and so have imperfect sensitivity and specificity \citep{Epstein2012}. Existing decision support tools that predict biopsy outcomes for AS patients (including, most recently, \citet{Ankerst2015}) provide patients and physicians with valuable information to guide decisions about biopsy timing and frequency but are insufficient to directly address patients' primary concerns about their tumors' lethality. Patients and clinicians need predictions of the pathological make-up of the entire prostate to guide their decision-making.

With this application in mind, we have developed a Bayesian hierarchical model that enables prediction of an individual's underlying disease state via joint modeling of repeated PSA measurements and biopsies. Specifically, we predict a binary cancer state-- \textit{indolent} or \textit{aggressive}-- with the latter defined as a ``true" Gleason score of 7 or higher. Predictions are informed by a subset of patients for whom the true state is observed-- patients who, either before or after biopsy grade reclassification, chose to undergo prostatectomy and have post-surgery, entire-prostate Gleason score determinations. In this sense, cancer state operates as a partially-latent class in the proposed model \citep{Wu2015}. 

An individual's cancer state is assumed to be manifest in both the level and trajectory of PSA measurements as well as in the outcomes from repeated biopsies. These relationships are illustrated by the directed acyclic graph (DAG) in Figure \ref{fig:dag1}. In the model we are proposing, PSA measurements follow a multilevel model with mean intercept and age effects varying across latent classes. Then, repeated annual biopsies constitute a time-to-event outcome since patients exit AS after grade reclassification on biopsy. So, time until reclassification on biopsy is modeled using pooled logistic regression under the assumption that biopsy results are independent conditional on cancer state and covariates \citep{Cupples1988}. Pooled logistic regression provides survival estimates equivalent to those of a time-varying Cox model for discrete event times and conditionally independent intervals \citep{D'Agostino1990}. As indicated in Figure \ref{fig:dag1}, PSA and biopsy results are also assumed to be conditionally independent given latent class. 

\begin{figure}
\begin{center}
\begin{subfigure}[b]{0.45\textwidth}
\includegraphics[width=\textwidth]{dag2}
\caption{Time-varying scenario}
\label{fig:dag1}
\end{subfigure}
\begin{subfigure}[b]{0.45\textwidth}
\includegraphics[width=\textwidth]{dag3}
\caption{Informative observation process scenarios}
\label{fig:dag2}
\end{subfigure}
\label{fig:dags}
\caption{DAGs describing the relationships between latent class (circled) and clinical outcomes (squared).}
\end{center}
\end{figure}

The model depicted in Figure \ref{fig:dag1} is related to previous work by \citet{Lin2002}, who proposed a joint latent class model (JLCM) to analyze longitudinal PSA and time-to-diagnosis of prostate cancer, extending earlier joint models by \citet{Schluchter1992}, \citet{DeGruttola1994}, and Henderson, Diggle, and Dobson (2000)\nocite{Henderson2000}. \citet{Inoue2008} used a Bayesian approach to jointly model PSA and time-to-diagnosis at various stages of disease in order to estimate the underlying natural history process for prostate cancer initiation and progression. \citet{Proust2009} developed a dynamic extension of the JLCM to predict prostate cancer recurrence after radiation therapy.

A credible statistical solution to the active surveillance problem requires three extensions of existing latent variable models for multivariate outcomes (such as the JLCM). First, the model must accommodate measurement error inherent in monitoring disease state. In our approach, we focus prediction on a partially-observed true Gleason score, instead of relying on biopsy Gleason scores for accurate characterization of the latent health state, and model a stochastic, rather than deterministic, relationship between the two.

Second, the model must allow disease monitoring to reflect patterns of clinical practice, including discrete, possibly informative, observation times. Specifically, prostate biopsies are scheduled to occur annually, but patients may opt to forgo the procedure. Our method replaces the JLCM's usual survival model for right-censored outcomes in favor of a pooled logistic regression model for biopsy grade reclassification where the possibility of reclassification in any year is conditional on a biopsy being performed. Furthermore, it is possible that the choice to receive a biopsy depends on the true cancer state or, more generally, that unobserved confounding exists, as shown by the dotted arrow from true cancer state to the ``Biopsy Performed" node in Figure \ref{fig:dag2}. If so, biopsy results are missing not at random (MNAR), and predictions of the true state that ignore the MNAR mechanism will be biased \citep{Little2014}. In response, our approach also includes a regression model for the probability of receiving a biopsy in each interval; the occurrence of a biopsy is allowed to depend on the latent health state, as well as previous biopsy and PSA observations.

Third, the active surveillance model must allow surgical removal of the prostate and subsequent observation of the underlying cancer state to be informative of that latent state. Consider the dotted arrow from true cancer state to the ``Surgical Removal" node in Figure \ref{fig:dag2}. If, after conditioning on clinical observations, an individual's true cancer state is associated with his choice to undergo surgery, whether through direct causation or unmeasured confounding, then informative missingness is present and failure to accommodate this in the model will result in biased predictions of the cancer state. While the association between the true cancer state and a binary indicator of its observation is not identifiable, we propose to model the time until surgery (and true state observation) conditional on the latent state. Evidence of this relationship among the subset of patients with surgery and mild assumptions on the structure of the hazard function (such as an additive or smooth effect) provide identifiability. This approach shares a similar intuition with missing data models for repeated attempt designs in which the estimated association between the number of attempts needed to elicit a response and its value is used to account for outcomes suspected to be MNAR \citep{Jackson2012}. In this application, patients have the opportunity to elect surgery throughout their participation. For simplicity, we also model the time until surgery with a pooled logistic regression model.

This paper is organized as follows. In Section 2, a hierarchical model for latent class prediction is described and estimation procedures are outlined. In Section 3, we specify our model to predict latent cancer states for patients in the Johns Hopkins Active Surveillance cohort and outline a simulation study based on this application. Results are presented in Section 4. We close with a discussion.

\vspace{-14pt}

\section{Hierarchical Latent Class Model}

We propose a Bayesian hierarchical model of the underlying cancer state, measurement process, and clinical outcomes of patients enrolled in active surveillance (AS). Predictions are made by incorporating information from repeated PSA and biopsy measurements for all patients and true cancer state observations in a potentially non-random subset of the cohort. Predictions are also informed by the presence of some observations, which we refer to as an \textit{informative observation process} (IOP). In this section, we introduce notation and conditional distributions for the observed data given the latent variables and parameters, then give the likelihood function. The model is completed by specifying appropriate priors and defining the joint posterior distribution. Overall model structure is summarized in Figure \ref{fig:model-summary}.

\subsection{Latent cancer state $\eta_i$ for patient $i$, $i=1,\dots,n$}

Define individual $i$'s true cancer state, $\eta_i$, as either indolent, $\eta_i=0$, or aggressive, $\eta_i=1$, $i=1,\dots,n.$ We use the Gleason score that would be assigned if his entire prostate were to be surgically removed and analyzed to define $\eta_i=0$ if Gleason $=$ 6 and $\eta_i=1$ if Gleason $\geq$ 7. Note that this definition assumes that cancer state is constant during the time under consideration. This assumption is discussed in more detail in Section 5.

True cancer state is then modeled as a Bernoulli random variable, $\eta_i \sim Bern(\rho_i)$. We assume a shared underlying probability of aggressive cancer, $\rho_i=\rho$, for simplicity in initial presentation. We observe this true cancer state on a possibly non-random subset of patients who choose surgical removal of the prostate and, hence, $\eta_i$ is a partially-latent variable.

\subsection{Longitudinal data $Y_{im}$ given latent class $\eta_i$, $m=1,\dots,M_i$}

Next, we consider PSA, which is influenced by the true cancer state $\eta_i$ as well as covariates including age and prostate volume. Unlike biopsies, PSA measurements are a routine part of each clinic visit so the times of observation are assumed to be independent of $\eta_i$. We use a multilevel model to estimate the linear trend (on a log scale) of an individual's PSA as he ages \citep{Gelman2006}. Patient-level coefficients, $\mathbf{b}_i$, vary about an $\eta_i$-specific mean intercept and slope ($\boldsymbol{\mu}_{\eta_i}$). Specification follows that of a hierarchically-centered multilevel model to speed convergence of the posterior sampling algorithm (Gelfand, Sahu, and Carlin, 1995)\nocite{Gelfand1995}. Specifically, given $\mathbf{b}_i$, the log-transformed PSA for patient $i$'s $m$th visit, $Y_{im}$, is assumed equal to $\mathbf{X}_{im}\boldsymbol{\beta} + \mathbf{Z}_{im}\mathbf{b}_i + \epsilon_{im}$ where $\mathbf{X}_{im}$ and $\mathbf{Z}_{im}$ are covariate vectors for individual $i$ at visit $m$, $\boldsymbol{\beta}$ is a parameter vector of population-level coefficients, and residual $\epsilon_{im}$ is assumed to follow a Gaussian distribution with mean zero and variance $\sigma^2$. In comparison to the commonly used mixed effects model of \citet{Laird1982}, covariates in $\mathbf{Z}_{im}$ are not a subset of covariates in $X_{im}$; covariates corresponding to patient-level effects $\mathbf{b}_i$ are only included in $\mathbf{Z}_{im}$, and the $\mathbf{b}_i$ are not centered at zero. In our application, $\mathbf{Z}_{im}$ includes an intercept and age so that PSA intercepts and slopes vary across individuals. $\mathbf{X}_{im}$ includes prostate volume, and $\boldsymbol{\beta}$ is the population-level association between volume and log-PSA.  

Modeling of patient-level coefficients follows the recommendation of \citet{Gelman2006} who advocate the use of a scaled inverse Wishart prior on the covariance matrix. The inverse Wishart prior, which is commonly used for Bayesian estimation of multilevel models \citep{Gelfand1995}, imposes dependence between variance and correlation components of the covariance matrix. To reduce prior dependence and allow for a flat prior on the correlation between individual-level intercepts and slopes, \citet{OMalley2008} introduce a scale parameter, $\boldsymbol{\xi}$, for the patient-level random effects:  $\mathbf{b}_i = diag(\check{\mathbf{b}_i}\boldsymbol{\xi}^T)$. Unscaled random effects, $\check{\mathbf{b}_i}$, are assumed to follow a  latent-class specific multivariate Gaussian distribution with mean vector $\boldsymbol{\mu}_{\eta_i}$ and covariance matrix $\Sigma_{\eta_i}$. 



\subsection{Biopsy Occurrence $B_{ij}$ and Result $R_{ij}$ for patient $i$ in time interval $j$, $j=1,\dots,J_i$}
We then consider information about the true cancer state contained in the occurrence and results of prostate biopsies. Biopsy data are categorized into discrete time intervals with $(B_{ij}, R_{ij})$ denoting binary outcomes for individual $i$ in time interval $j$. $B_{ij}$ indicates whether a biopsy was performed ($B_{ij}=1$) or not ($B_{ij}=0$) and, when it was performed, $R_{ij}$ indicates if grade reclassification occurred ($R_{ij}=1$) or not ($R_{ij}=0$). $B_{ij}$ and $R_{ij}$ are defined for $j=1,\dots, J_i$,  where $J_i$ is the time interval of reclassification or censoring for patient $i$. For each time interval, we use logistic regression to model the occurrence of a biopsy and, when a biopsy was performed, its result; both outcomes are conditional on true cancer state:
\bea
\label{eq:p_bx}
\text{logit} \{ \, P(B_{ij}=1 | \eta_i, \mathbf{U}_{ij}, \boldsymbol{\nu}) \, \}= \mathbf{U}_{ij}\boldsymbol{\nu}_1 + \eta_i\boldsymbol{\nu}_2 + \mathbf{U}_{ij}\eta_i\boldsymbol{\nu}_3 \\
\label{eq:p_rc}
\text{logit}\{\, P(R_{ij}=1 | \eta_i, \mathbf{V}_{ij}, B_{ij}=1, \boldsymbol{\gamma}) \,\}= \mathbf{V}_{ij}\boldsymbol{\gamma}_1 + \eta_i\boldsymbol{\gamma}_2 + \mathbf{V}_{ij}\eta_i\boldsymbol{\gamma}_3
\eea
where $\mathbf{U}_{ij}$ and $\mathbf{V}_{ij}$ are covariate vectors including time-varying predictors and $\boldsymbol{\nu}=(\boldsymbol{\nu}_1, \boldsymbol{\nu}_2, \boldsymbol{\nu}_3)$ and $\boldsymbol{\gamma}=(\boldsymbol{\gamma}_1, \boldsymbol{\gamma}_2, \boldsymbol{\gamma}_3)$ are parameter vectors to be estimated that include the main effects of covariates $\mathbf{U}_{ij}$ or $\mathbf{V}_{ij}$, $\eta_i$, and the possible interactions $\mathbf{U}_{ij}\eta_i$ and $\mathbf{V}_{ij}\eta_i$, respectively. Since reclassification occurs at most once, Equation (\ref{eq:p_rc}) corresponds to a modified pooled logistic regression model for time-to-reclassification in which only intervals with biopsies contribute.

This model specification represents three important aspects of data generated in active surveillance:  whether a biopsy is performed may be informative of true cancer state, time-to-reclassification depends on a patient's decision to receive a biopsy, and biopsy outcomes are prone to measurement error. In this application, $U_{ij}$ and $V_{ij}$ may include age, time since diagnosis, and calendar date. Previous PSA and biopsy results may also influence the decision to get a biopsy, but they do not influence biopsy findings. 

\vspace{-8pt}

\subsection{Surgical Removal of Prostate $S_{ij}$ and its Cancer Lethality $\eta_i$}
Lastly, to allow for the possibility that surgical removal of the prostate (and subsequent observation of the true cancer state) is informative, we define $S_{ij}$ to be a binary indicator of surgery ($S_{ij}=1$) or not ($S_{ij}=0$) for individual $i$ during time interval $j$ for $j=1,\dots,J_{S_i}$, where $J_{S_i}$ is the time of surgery or other censoring for patient $i$ and $J_{S_i} \geq J_i$ for all $i$. The probability of surgery in each time interval is modeled with logistic regression and conditional on the true cancer state: $\text{logit} \{\,P(S_{ij}=1 | \eta_i, \mathbf{W}_{ij}, \boldsymbol{\omega})\,\} =   \mathbf{W}_{ij}\boldsymbol{\omega}_1 + \eta_i\boldsymbol{\omega}_2 + \mathbf{W}_{ij}\eta_i\boldsymbol{\omega}_3 \nonumber$ where $\mathbf{W}_{ij}$ is a vector of time-varying predictors and $\boldsymbol{\omega}=(\boldsymbol{\omega}_1, \boldsymbol{\omega}_2, \boldsymbol{\omega}_3)$ is a parameter vector to be estimated. Age, time since diagnosis, calendar date, and previous PSA and biopsy results may all be considered as possible predictors of surgery. 


\vspace{-8pt}

\subsection{Posterior Distribution Estimation}
Having specified models for each information source, we define the likelihood of the latent states, patient-level coefficients, and population-level parameters given the observed data as the product of the contribution of each component described above:
\bea
&&L\big( \, \textbf{parameters, patient-level latent class and coefficients}\, | \,\textbf{data}\,\big) \nonumber \\ 
&& \quad =  \prod_{i=1}^{n}  \,  \rho^{\eta_i}\,(1-\rho)^{1-\eta_i}\, f(\mathbf{Y}_i | \underline{\mathbf{X}_i}, \underline{\mathbf{Z}_i}, \boldsymbol{\beta}, \boldsymbol{\xi}, \check{\mathbf{b}_i}, \sigma^2) \, g(\check{\mathbf{b}_i} |  \boldsymbol{\mu}_{\eta_i}, \Sigma_{\eta_i})  \nonumber \\
&& \qquad \qquad \prod_{j=1}^{J_i}  P(B_{ij}=1 | \eta_i, \mathbf{U}_{ij}, \boldsymbol{\nu}) ^{B_{ij}} P(B_{ij}=0 | \eta_i, \mathbf{U}_{ij}, \boldsymbol{\nu})^{1-B_{ij}}  \nonumber \\
&& \qquad \qquad \qquad  \qquad \{P(R_{ij}=1 | \eta_i, \mathbf{V}_{ij}, \boldsymbol{\gamma}) ^{R_{ij}} P(R_{ij}=0 |  \eta_i, \mathbf{V}_{ij}, \boldsymbol{\gamma})^{1-R_{ij}} \}^{B_{ij}}\nonumber \\
\label{eq:lik-inf}
&& \qquad \qquad  \prod_{j=1}^{J_{S_i}}  P(S_{ij}=1 | \eta_i, \mathbf{W}_{ij}, \boldsymbol{\omega}) ^{S_{ij}} P(S_{ij}=0 |  \eta_i, \mathbf{W}_{ij}, \boldsymbol{\omega})^{1-S_{ij}} 
\eea
where $f$ and $g$ are multivariate normal densities for the vector of log-transformed PSAs $\mathbf{Y}_i$ and unscaled patient-level effects $\check{\mathbf{b}_i}$, respectively, each with mean and covariance as defined in Section 2.2. $\underline{\mathbf{X}_i}$ denotes the matrix of covariate vectors $[\mathbf{X}_{i1},\dots, \mathbf{X}_{iM_i}]$; $\underline{\mathbf{Z}_i}$ is similarly defined.


We use a Bayesian approach for model estimation. The prior distributions for the latent states and patient-level coefficients have been described (Sections 2.1 and 2.2, respectively). Standard prior distributions are used for model parameters, including a beta prior on the probability of having aggressive cancer ($\rho$) and minimally informative Gaussian priors on logistic regression model coefficients, as shown in the model summary given in Figure \ref{fig:model-summary}. The joint posterior distribution of the parameters, latent states, and patient-level effects is proportional to the product of the likelihood and joint prior density of model parameters and is given explicitly in the online supplement.

\begin{figure}
\begin{center}
\includegraphics[width=\textwidth]{model-summary-table.pdf}
\caption{Model summary with priors used for application to Johns Hopkins Active Surveillance data. $D_X$ is the length of vector $\mathbf{X}$ and $\mathbf{I}_{D_X}$ is the identity matrix with dimension $D_X \times D_X$. $D_Z, \, D_U, \, D_V,$ and $D_W$ and the associated identity matrices are similarly defined for covariate vectors $\mathbf{Z},\,\mathbf{U},\,\mathbf{V},$ and $\mathbf{W}$.}
\label{fig:model-summary}
\end{center}
\end{figure}


For those patients without $\eta_i$ observed, a data augmentation approach is used to sample the true unobserved cancer state from its full conditional posterior at each iteration of the MCMC sampling algorithm \citep{Tanner1987}. Averaging the resulting posterior sample produces the posterior probability that a patient has a true Gleason 7 or higher prostate cancer, $P(\eta_i=1)$.



\section{Johns Hopkins Active Surveillance Cohort}

\subsection{The Data}
From January 1995 to June 2014, the Johns Hopkins Active Surveillance (JHAS) cohort enrolled 1,298 prostate cancer patients \citep{Tosoian2015}. This study prospectively follows patients with very-low-risk or low-risk prostate cancer diagnoses (according to criteria outlined in \citet{Epstein1994}) who elect to delay curative intervention in favor of active surveillance (AS). Results of all prior PSA tests and diagnostic biopsies are collected at enrollment. As part of the surveillance regimen, PSA tests are performed every six months and biopsies are performed annually, though biopsy intervals may vary based upon patient preferences and clinician recommendations. Treatment is recommended upon biopsy grade reclassification, that is, when the Gleason score assigned on a biopsy first exceeds six. Some patients also choose to undergo treatment prior to reclassification. For patients who elect surgical removal of the prostate, the true Gleason score assigned to the entire prostate after pathologic assessment is recorded when available.

A total of 874 patients who met study criteria and had at least two PSA measurements and at least one post-diagnosis biopsy as of October 1, 2014 were included in the analysis. Patient outcomes are given in Figure \ref{fig:consort}. Grade reclassification was observed in 160 patients (18$\%$ of the analysis cohort). Notably, over a quarter of patients with grade reclassification who underwent prostatectomy were downgraded after surgery (17/65) while nearly a third of patients who underwent prostatectomy in the absence of grade reclassification were upgraded (30/96). Further details on the analysis dataset are given in the online supplement. 

\begin{figure}
\begin{center}
\includegraphics[width=\textwidth]{consort-diagram-biometrics.pdf}
\caption{CONSORT diagram for Johns Hopkins Active Surveillance prospective cohort patients included in this analysis. Post-surgery full prostate Gleason score (GS) observations are also given (circled). (Six patients who underwent prostatectomy did not have true GS observations available.)}
\label{fig:consort}
\end{center}
\end{figure}


\subsection{Model Specification}
We applied models with and without the biopsy and surgery informative observation process (IOP) components to data from the JHAS cohort. 

PSA observations were modeled with a hierarchically-centered multilevel model, as described in Section 2.2. Patient-level coefficients for intercept and age were estimated for each patient. A shared covariance matrix was assumed for the unscaled patient-level effects, that is,  $\Sigma_0=\text{Var}(\check{\mathbf{b}_i}|\eta_i=0) = \text{Var}(\check{\mathbf{b}_i}|\eta_i=1) = \Sigma_1,$ in order to reduce model complexity. (The plausibility of this assumption was checked by examining estimated covariance matrices in the subset of patients with known cancer state.) The PSA model also included a population-level coefficient for prostate volume.

Biopsy, reclassification, and surgery observations were categorized into annual intervals, and exploratory data analysis was performed to identify predictors of each. Covariates were selected that lowered the Akaike Information Criterion (AIC) of a multivariable logistic regression model for each outcome \citep{Akaike1998} and are listed in Figure \ref{fig:model-summary}. Natural splines with up to four degrees of freedom and knots at percentiles of the predictor variable were used when doing so lowered the AIC. Additional details are given in the online supplement. 

Model parameters and their minimally informative priors are presented in the model summary given in Figure \ref{fig:model-summary}. Posterior sampling was performed with \texttt{JAGS} \citep{Plummer2015} via the \texttt{R} package \texttt{R2JAGS} \citep{R2jags2015}. Parallel chains were run to confirm model estimates converged to similar values. Cumulative quantile and trace plots were also used to monitor convergence. Analysis code, posterior sampler settings, and diagnostic plots are included in the supplementary material. 



\subsection{Model Assessment}

Predictive accuracy was assessed among patients with post-surgery true Gleason score observations. Out-of-sample posterior predictions of $\eta$ were obtained for each patient by removing his true state observation from the analysis dataset and re-running the posterior sampler with an additional data augmentation step for the patient of interest. Out-of-sample predictions of $\eta_i$ were then compared to known values with receiver operating characteristic (ROC) curves \citep{Hanley1982} and calibration plots \citep{Steyerberg2010}. For the former, the area under the curve (AUC) and associated 95$\%$ bootstrapped intervals were calculated. For the latter, a plot comparing posterior predictions to observed rates of class membership was constructed by performing logistic regression of the observed true state on a natural spline representation of out-of-sample posterior predictions (degrees of freedom = 2). The mean squared error (MSE) between observed and predicted cancer state was also calculated. For comparison, posterior predictions were obtained from a logistic regression model fit with data from patients with post-surgery observations of $\eta$; covariates included age, time since diagnosis, and PSA and biopsy results. We also compared specificity of model predictions to the specificity of using final biopsy results to predict the true cancer state by fixing sensitivity at the observed true positive rate of biopsy Gleason score (dichotomized $<$ 7 or $\geq 7$).

Calibration plots were also drawn to assess model fit for outcomes observed on all patients: the occurrence of a biopsy, grade reclassification on biopsy, and the occurrence of surgery.
Code for reproducing all plots is available in the supplementary material





\subsection{Simulations}
We performed a simulation study to examine model performance in this application. 200 simulated datasets were generated using posterior estimates of model parameters obtained from the biopsy and surgery IOP analysis of JHAS data. For each dataset, the proposed model was estimated under four settings: unadjusted (no IOP components), biopsy IOP only, surgery IOP only, and both biopsy and surgery IOP components. Posterior predictions of the latent state were obtained for all simulated patients and compared to known (data-generating) values. For patients without surgery, posterior samples of $\eta$ were generated with a data augmentation step as a matter of course in model estimation. For patients with surgery, the posterior probability of $\eta=1$ was estimated via an importance sampling algorithm performed on the joint posterior \citep{Bishop2006}, which is less computationally intensive than the out-of-sample methods used in Section 3.3. (See technical report of \citet{Fisher2015} for further details.) Posterior predictions of $\eta$ were also compared to fitted probabilities from a logistic regression model. Code for generating data, estimating the joint posterior, and obtaining predictions is included in the supplementary material.

\vspace{-14pt}

\section{Results}



The estimated marginal probability of harboring a prostate cancer with Gleason score above 7 was 0.23 (95$\%$ CI: 0.16, 0.33) for the proposed model with biopsy and surgery IOP components, 0.20 (0.14, 0.28) with surgery IOP only, 0.31 (0.24, 0.39) with biopsy IOP only and 0.30 (0.23, 0.38) with no IOP components. Patients with $\eta=1$ were less likely to receive biopsies--leading to underestimation of $\rho$ in models without the biopsy IOP component--and more likely to elect surgery, such that $\rho$ was overestimated when not accounting for informative observation. Parameter estimates and credible intervals from all models are given in the online supplement (Appendix Tables A3-A7). 

A histogram of predictions of $\eta$ from the model with biopsy and surgery IOP components is given in Figure \ref{fig:hist-eta-pred}. Patients with posterior predictions above 60$\%$ are primarily those who both experienced grade reclassification (solid bars) and elected prostatectomy (red and green). 95$\%$ of AS patients who neither reclassified (diagonal shading) nor underwent surgery (black) have posterior predictions that are lower than 50$\%$; a majority have predictions below 20$\%$. Figure \ref{fig:scatter-eta-pred} shows a scatterplot comparing posterior probabilities of aggressive cancer, $P(\eta=1)$, between models with and without IOP components. The models produce similar posterior predictions for most patients, particularly those patients for which the non-IOP model assigns very low risk. Inclusion of biopsy and surgery IOP components decreases posterior predictions of $\eta$ most markedly for patients with frequent biopsies and no surgery (larger black circles below the x=y axis) and tends to increase posterior predictions for patients who elect surgery or have infrequent biopsies (colored or smaller black circles, respectively, above the x=y axis). These trends are further illustrated by density plots of predictions of $\eta$ stratified by reclassification and surgery given in the online supplement.



\begin{figure}
\begin{center}
\begin{subfigure}[b]{0.49\textwidth}
\includegraphics[width=\textwidth]{hist-pred-eta-stratified}
\caption{Histogram of $B,\,S$ IOP predictions}
\label{fig:hist-eta-pred}
\end{subfigure}
\begin{subfigure}[b]{0.49\textwidth}
\includegraphics[width=\textwidth]{pred-eta-iop-vs-unadj}
\caption{$B,\,S$ IOP vs Non-IOP predictions}
\label{fig:scatter-eta-pred}
\end{subfigure}
\caption{Posterior predictions of true prostate cancer state. On both plots, coloring indicates whether $\eta$ was observed and, if so, its value. In the histogram (a), diagonal shading represents patients whose final biopsy was assigned a Gleason score of 6 while solid bars represent patients whose final biopsy was assigned a Gleason score of 7 or higher (i.e., grade reclassification). In the scatterplot (b), circle size indicates the frequency with which a patient received prostate biopsies; larger circles represent more frequent biopsies while smaller circles represent less frequent biopsies.}
\label{fig:post-preds-eta}
\end{center}
\end{figure}



Posterior predictions of $\eta$ from the proposed model with biopsy and surgery IOP components were more accurate than those from the proposed model with a single or no IOP components or from logistic regression (Figures \ref{fig:jhas-auc} and \ref{fig:roc}). The out-of-sample AUC among patients with observed true cancer state is highest for the biopsy and surgery IOP model (0.75, 95$\%$ bootstrapped interval: 0.67, 0.83), and the MSE from this model was also the lowest (0.201, 95$\%$ Int: 0.17, 0.24; Table 9 in online supplement). While this improvement is slight, it is widely recognized that drastic increases in classification accuracy are rare to achieve \citep{Pepe2004}. The false positive rate (FRP) of predictions from the biopsy and surgery IOP model (0.14, 95$\%$ Int: 0.07, 0.23) is also lower than that of the binary classifier based on final biopsy results (FPR=0.20, 95$\%$ Int: 0.12, 0.29) at a fixed true positive rate of 0.62 (the sensitivity of final biopsies in patients with eventual surgery). The improvement in specificity offered by the biopsy and surgery IOP model corresponds to avoiding, on average, 30$\%$ of unnecessary diagnoses of more aggressive cancer in comparison to a diagnosis based solely on a patient's most recent biopsy (95$\%$ bootstrapped interval of FPR($R$) - FPR($\eta$): -7.5$\%$, 13$\%$). These comparisons are limited because accuracy of posterior predictions can only be assessed among patients with $\eta$ observed. Yet, we expect the predictive accuracy gained by incorporating IOP components to be seen more definitively in patients without true state observations if biopsy and surgery results are indeed MNAR. This is explored further with simulations.

Posterior predictions of $\eta$ from the IOP model also appear to accurately estimate a patient's risk of having more aggressive cancer. The calibration plot in Figure \ref{fig:calibration-eta} shows that, for patients with known values of $\eta$, the average posterior predicted probability of $\eta=1$ is close to the average observed value of $\eta$, indicating that the model reasonably reproduces the mean of observations. The risks of clinical outcomes (biopsy results) and choices (occurrence of biopsy and surgery) for all patients appear to be accurately estimated by the IOP model as well, as demonstrated by calibration plots in the online supplement (Appendix Figure A3).


\begin{figure}
\begin{center}
\begin{subfigure}[b]{0.9\textwidth}
\includegraphics[width=\textwidth]{JHAS-pred-accuracy-eta}
\caption{}
\label{fig:jhas-auc}
\end{subfigure}

\begin{subfigure}[b]{0.9\textwidth}
\vspace{20pt}
\includegraphics[width=\textwidth]{SIM-pred-accuracy-eta}
\caption{}
\label{fig:auc-sim}
\end{subfigure}
\caption{Estimated AUC for predictions of $\eta$ from the JHAS cohort (a) and simulation studies (b). Intervals in (a) are quantile-based 95$\%$ intervals from 10,000 bootstrap samples of patients with post-surgery observations of $\eta$. Intervals in (b) are quantile-based 95$\%$ intervals from the estimated AUC in 200 simulation studies.}
\end{center}
\end{figure}



\begin{figure}
\begin{center}
\begin{subfigure}[b]{0.475\textwidth}
\includegraphics[width=\textwidth]{pred-vs-obs-eta-roc-annotated}
\caption{ROC curves: $B,S$ IOP $\&$ non-IOP predictions}
\label{fig:roc}
\end{subfigure}
\begin{subfigure}[b]{0.475\textwidth}
\includegraphics[width=\textwidth]{pred-vs-obs-eta-calibration-annotated}
\caption{Calibration plot: $B,\,S$ IOP predictions}
\label{fig:calibration-eta}
\end{subfigure}
\caption{Predictive accuracy of out-of-sample predictions of $\eta$ among patients with true state observed in the JHAS cohort. In (a), the specificity of predictions from each model is highlighted at the sensitivity of a binary classifier defined by final biopsy result (*). In (b), the dark line shows the empirical rate of observing a true Gleason score of 7 or above (y-axis) given an out-of-sample posterior probability of true state (x-axis) under the model with biopsy and surgery IOP components; shading gives the 95$\%$ point-wise confidence interval. Perfect agreement lies on the x=y axis (dotted line). Hashmarks at y=0 and y=1 correspond to observed cancer states ($\eta=0,1$, respectively) for patients with post-surgery true state observations. Hashmarks are located along the x-axis at each patient's out-of-sample posterior probability of the true state.}
\end{center}
\end{figure}



Estimation of the time-to-surgery model depends on sufficient evidence among patients with surgery of a relationship between $\eta$ and surgery time. To assess the strength of this evidence, we re-ran the biopsy and surgery IOP model with more informative priors on the coefficients capturing this relationship and found posterior estimates to be robust. Details are given in the online supplement (Appendix Section A3.1.6 and Figure A18).

Through simulation studies, we found that our sampling procedure produced unbiased estimates with nominal coverage (Tables 10-14 in the online supplement) and that the proposed model with biopsy and surgery IOP components outperforms other model variations and logistic regression when it correctly reflects the data-generating mechanism. The AUC for each model among patients with and without post-surgery true state observations are given in Figure \ref{fig:auc-sim}. Differences in predictive accuracy across models are similar in magnitude to those observed in the application (Figure \ref{fig:jhas-auc}).  As expected, we also see that the logistic regression model, which is only estimated with data from patients with true state observations, has the poorest predictive accuracy among patients without surgery. Predictions from models incorporating the surgery IOP component show appropriate calibration while those without overestimate the risk of aggressive prostate cancer (Figure 20, online supplement).


\section{Discussion}
In this paper, we have presented a hierarchical Bayesian model for predicting latent cancer state among low risk prostate cancer patients. Multiple models have been developed to predict biopsy results in this population \citep{Ankerst2015, Truong2013}. However, our model predicts the outcome of chief interest--the true underlying state of an individual's prostate cancer. Focusing on the actual health state, even when latent, is equivalent to subsetting patients into subgroups for which optimal treatments differ. Subsetting is the goal of precision medicine \citep{Saria2015}.


The proposed model integrates four sources of information about whether a tumor is aggressive or indolent: repeated measures of the biomarker PSA; repeated results from tissue biopsies; repeated decisions to have a biopsy; and the time to surgical removal of the prostate. In the subset of patients that have their prostates removed, the true tumor pathology state is observed. This data-integrating method is an example of semi-supervised learning because patients both with and without true state observations are included in model estimation (Chapelle, Sch{\"o}lkopf, and Zien, 2006)\nocite{Chapelle2006}. We adjust for possible informative missingness by modeling the time until surgery depending on the true state. While it is ideal to assess model sensitivity to parametric assumptions embedded in selection models for missing not at random data mechanisms \citep{Daniels2008}, existing methods for re-parameterizing selection models as pattern mixture models do not accommodate event time outcomes with the possibility of censoring. Further research is needed to develop these methods. 

The methods proposed here are tailored to available measurements that address the clinical questions arising from active surveillance of prostate cancer: should I have a biopsy this year; what is the chance my tumor is indolent; should I undertake removal or irradiation of my prostate despite the known serious side effects? The model extends naturally to provide improved answers to these questions as additional data become available. For example, when genetic markers for prostate cancer risk are identified, the probability distribution for latent state ($\rho_i$) could easily be informed by subgroups defined by their expression. Or, when MRI or ultrasound images are commonly used before biopsy, these data will be included in the model as well. In the case that some measurements are not available for all patients, the proposed framework is also able to adjust for informative observation of predictors and outcomes. 

The proposed model can also be modified in response to advancement in scientific understanding about the relationship between clinical measurements and the underlying cancer state. In particular, in the event of new research findings on the rate of progression in this population, the model could be extended to allow an indolent cancer to transition to a lethal one, for example, as a Markov process. Because an individual's true cancer state can only be observed once, the current data contain insufficient information to simultaneously support identifiability of both the rate of biopsy misclassification and the rate of pathological progression in the underlying state. The model currently assumes that an individual's cancer categorization (Gleason score) does not change over the time period under surveillance while allowing for imperfect sensitivity and specificity of biopsies. This assumption reflects the current clinical understanding that biopsy upgrading in AS is more frequently due to misdiagnosis rather than true grade progression \citep{Porten2011}. A more recent analysis by \citet{Inoue2014} suggests a rate of disease progression in the JH AS cohort of 12-14$\%$ within a decade of enrollment, but this estimate is sensitive to prior specification. A dynamic state extension would require strong prior knowledge about the progression rate parameter in order to be identifiable from the current data. The effect of allowing for a state transition would be to give greater weight to more recent PSA and biopsy outcomes when predicting the underlying state rather than giving equal weight to all observations. 


The proposed prediction model exemplifies the statistical underpinnings of a learning health care system (Goolsby, Olsen, and McGinni, 2012; Smith et al., 2013)\nocite{IOM2012,Smith2013}, a system with the ability to continuously integrate patient data and medical knowledge to optimize patient care. As more patients enroll in the Johns Hopkins active surveillance cohort, and as more information is collected on existing patients, our ability to predict underlying health states and the likely trajectory of clinical outcomes will improve. Furthermore, importance sampling methods can be used to obtain real-time prediction updates based on the most current information in order to support decision-making in a clinical setting. An example interactive decision-support tool that provides fast predictions of a patient's latent prostate cancer state is demonstrated at \texttt{https://rycoley.shinyapps.io/} \texttt{dynamic-prostate-surveillance}.





\section*{Supplementary Materials}
Simulated data, \texttt{JAGS} scripts, and \texttt{R} code to reproduce the analysis and figures are provided at \texttt{http://github.com/rycoley/prediction-prostate-surveillance}.

\bibliographystyle{biom}
