\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}


\section{Introduction} \label{Introduction}
Spatial resolution is always an important parameter for position-sensitive particle detectors, such as Micro-Pattern Gaseous Detectors~(MPGD). It is well known that a Gas Electron Multiplier~(GEM)~\cite{Sauli97} detector with straight readout strips of 400~$\mu$m pitch can reach 50-70~$\mu$m resolution \cite{COMPASSGEM, KondoNIMA}. In applications with large-area GEM detectors, the use of readout strips with such small pitch will quickly lead to a large number of required electronic channels and consequently to significant system cost incurred by the readout electronics. Readout structures employing pads or strips with a chevron or zigzag geometry were proposed and studied in the past to address this issue for MPGDs and other gaseous detectors~\cite{zz85, zz89, zz97, zzBYu03, zzBYu05}. A zigzag strip covers a larger area than a standard straight strip, so that zigzag strips can reduce the number of electronic channels needed to read out a given detector area. The challenge then is to make sure that the spatial resolution is still adequate and that the crosstalk between strips is well controlled.

We have previously demonstrated for small Triple-GEM detectors that parallel zigzag strips achieve good spatial resolution as the sensitivity of the charge sharing among strips to the hit position is enhanced between the interleaved ``zigs'' and ``zags'' of adjacent strips~\cite{IEEEzigzag}. Here we expand this approach to large-area GEM detectors with a radial readout strip geometry that is appropriate for detector systems with a disk or ring geometry as is commonly used in collider experiments. We report results from a beam test of a 1-meter-long GEM detector read out with wide radial zigzag strips.

This paper is organized as follows: In section \ref{theDetector} the design of the zigzag readout board is described in detail; in section \ref{setup} the beam test configuration and data acquisition system are introduced; in section \ref{basic_results} the data analysis methods and basic performances of the zigzag GEM detector are discussed; in section \ref{resolution} the angular resolution, i.e.\ the spatial resolution in azimuthal $\phi$-coordinate, of the zigzag GEM detector is studied under different operating conditions.

\section{The GEM detector and the zigzag readout strips} \label{theDetector}
The 1-meter-long Triple-GEM detector has the shape of a trapezoid with 22~cm width at the narrow end and 45~cm width at the wide end (Fig.~\ref{zigzagGEM}). It is based on a prototype for the CMS muon GEM upgrade at LHC~\cite{CMSGEM, CMSGEMTDR}, but has a modified readout board. The gas gaps (drift, transfer 1, transfer 2, induction) in this GEM detector are 3/1/2/1~mm.

On the printed circuit board (pcb) used for the readout, the zigzag strips are radially divided into eight sectors with lengths 10~cm (15~cm) at the narrow (wide) end of the trapezoid and with 128 strips per sector (Fig.~\ref{zigzagGEM}). In a collider geometry these sectors would correspond to different ranges of pseudorapidity $\eta$. From the narrow end to the wide end, the sectors are labeled from one to eight. The reason for this partitioning is to provide some coarse information on the radial position of hits and to have readout strips with reasonably low capacitance. In Fig.~\ref{zigzagDetail}, the details of the zigzag structure are shown. The zigzag strips run in radial direction and measure the azimuthal $\phi$ coordinate. The opening angle between the first strip and the last strip in a sector is 10$^\circ$ and the angular pitch between two neighboring strips is 1.37~mrad. The spacing between zigzag tips on adjacent strips is 0.1~mm and the distance between two tips in one strip is 0.5~mm. These geometric parameters are the same for all strips in all eight sectors on the readout board; consequently the strip width increases from $\sim$2.5~mm at the narrow end to $\sim$4.5~mm at the wide end. The total number of zigzag strips on the readout board is 1,024, to be compared with, e.g., 3,072 straight strips in the standard CMS GEM detector that the zigzag detector is derived from~\cite{CMSGEMTDR}. This means that 2/3 of the electronics can be eliminated by using zigzag strips.

A miniaturized ceramic high voltage divider is used for powering the GEM detector. The HV divider has a resistor chain with a total resistance of 4.4 M$\Omega$ and connects to each electrode of the GEM detector through eight pads on the drift board. One HV input channel is connected to the drift electrode and potentials between any two neighboring detector electrodes are provided through resistors in the HV divider. More details of the HV divider can be found in reference~\cite{CMSGEMTDR}.

\section{Beam test configuration and data acquisition system}\label{setup}
The zigzag GEM detector was tested as a tracking detector at the Fermilab Test Beam Facility~(FTBF) in October 2013. The eRD6-FLYSUB consortium~\footnote{eRD6-FLYSUB is a consortium of researchers from \textbf{FL}orida Institute of Technology, \textbf{Y}ale University, \textbf{S}tony Brook University, \textbf{U}niversity of Virginia, \textbf{B}rookhaven National Lab, formed to carry out R\&D of tracking and particle identification detectors for a future electron-ion collider.} conducted this beam test with a variety of GEM detectors. Fig.~\ref{detectorConfig}~(top) shows a diagram of the detector configurations in the tracking system, where a set of four GEM detectors with 2-D readout acts as a reference tracker (noted as ``Ref'' in the figure). The zigzag GEM detector was installed on a movable table in the center between the tracker detectors (Fig.~\ref{detectorConfig}, bottom). All detectors were operated with an Ar/CO$_2$ 70:30 gas mixture during the entire test period.

The RD51 Scalable Readout System (SRS)~\cite{SRS} was used to read out all ten detectors in the beam test. With 64 frontend hybrids carrying APV25 chips~\cite{APVUserGuide} and four SRS Frontend-Concentrator/ADC combinations, 8,192 channels were read out simultaneously~(Fig.~\ref{fig4}). Strip signals were digitized at 40 MHz. Readout was triggered by coincidence signals from plastic scintillators placed in the beam line. Data were transferred via Gigabit ethernet to a PC that acquired them using the DATE software and monitored them online in the \textbf{A}utomatic \textbf{MO}nito\textbf{R}ing \textbf{E}nvironment~(AMORE) (both software suites were originally developed for the ALICE experiment and later adapted by the RD51 collaboration for the SRS)~\cite{ALICEGuide}.

\section{Data analysis and offline results} \label{basic_results}
The data are analyzed with the AMORE software package as well; raw events are decoded and basic information on the events are retrieved. Events with multiple hits in any detector are excluded in the analysis for simplicity. Beam profiles are checked first to characterize the impinging beams. Then basic performance characteristics of the zigzag GEM detector are analyzed, such as strip multiplicity of strip clusters, i.e.\ the number of strips with induced charge on the strips above a certain threshold, strip-cluster charge distribution and detection efficiency. For simplicity, we will refer to a strip cluster as a ``cluster" throughout the remainder of this paper. Crosstalk between zigzag strips is also investigated in this section.

\subsection{Beam profiles}
During the beam test, we tested detectors with the primary 120~GeV/c proton beam and with secondary beams containing mixed hadrons~(mainly pions with an admixture of kaons and protons) at 20, 25, and 32~GeV/c momenta. Most of the data for the zigzag GEM detector were taken with mixed hadron beams. The 2D hit maps of the tracker detectors~(Fig.~\ref{fig5}) show that the beam spot of the pure proton beam is $\sim$2~cm in diameter, while in the case of mixed hadron beams it is $\sim$7~cm in X direction (horizontal) and $\sim$3~cm in Y direction (vertical).

\subsection{Basic performance characteristics of the zigzag GEM detector}

\subsubsection{Cluster charge vs.\ HV}
Before taking data for any scenario, we took a pedestal run with 5k events at standby detector voltages so that there was no gas gain. To find a cluster in data runs with full voltage, we set a 5$\sigma$ threshold cut where $\sigma$ is the width of the pedestal distribution. If N contiguous strips are fired in a detector with individual strip charges higher than the 5$\sigma$ cut, then we will consider this as a cluster with strip multiplicity N in the detector and accept this event.

The zizgzag GEM detector was tested at different high voltages $\mathrm{V_{drift}}$ applied to the drift electrode when a 25~GeV/c mixed hadron beam was impinging on central readout sector number five (Fig.~\ref{zigzagGEM}). Fig.~\ref{fig6}~(top) shows the total cluster charge distribution  at $\mathrm{V_{drift}}$ = 3200~V in terms of ADC counts (1 ADC $\simeq$ 0.0371~fC). The measured cluster charge distribution fits well to a Landau function. Fig.~\ref{fig6}~(bottom) shows the most probable values~(MPV) of the Landau fits as a function of $\mathrm{V_{drift}}$; this yields an approximately exponential curve as expected for an MPGD. The gas gain during the beam test is estimated to range from 660~(at $\mathrm{V_{drift}=3000~V}$) to $\mathrm{10^4}$~(at $\mathrm{V_{drift}=3400~V}$) based on the measured mean cluster charge. After the beam test, we measured the gas gain of this detector operated with Ar/CO$_{2}$ 70:30 in our lab in Florida (sea level) with X rays (AMPTEK Mini-X X-ray generator with gold anode, set at 10~kV and 5~$\mu$A) impinging on  central sector five. A good rate plateau is obtained and the gain reaches close to $\mathrm{3\times10^4}$ at $\mathrm{V_{drift}=3600~V}$ (Fig.~\ref{figGain}) under those conditions.

\subsubsection{Strip multiplicity vs.\ HV}
The strip multiplicity distribution of clusters measured at $\mathrm{V_{drift}}$ = 3200~V in central sector five is shown in Fig.~\ref{fig8}~(top); the mean strip multiplicity is 1.63. This mean value is found to increase quadratically with $\mathrm{V_{drift}}$ (Fig.~\ref{fig8}, bottom), and at the highest tested voltage~(3400~V) its maximum is just below three strips. This is expected since the strip width in that radial region is $\sim$3.6 mm.

\subsubsection{Detection efficiency vs.\ HV}
The detection efficiency of the zigzag GEM detector for mixed hadrons, which are close to minimally ionizing in this test beam, is obtained from the ratio of the number of observed hits to the total number of triggered events. The measured efficiency is (98.4~$\pm$~0.2)\% on the plateau when a 5$\sigma$ threshold cut is applied to the measured strip charge.  Fig.~\ref{fig9} shows efficiencies with different threshold cuts demonstrating that the plateau efficiency is not affected significantly by the applied threshold.

\subsubsection{Response uniformity from position scan data}
The zigzag GEM detector was tested (at $\mathrm{V_{drift}}$ = 3200~V) with 20~GeV/c mixed hadron beams impinging on two different spots in each readout sector, approximately 60~mm apart in the vertical direction, so that the response uniformity of the chamber could be examined.  Fig.~\ref{fig10} shows the MPVs of Landau fits to cluster charge distributions measured in sectors one to seven at these different positions. The sector eight was not tested because it could not be reached by moving the stage that the detector was placed on. The response varies by $\sim$25\% and is lower in sectors six and  seven, which are on the wide side of the trapezoidal chamber. This non-uniformity is likely caused by a known slight bending of the drift board that was occurred when the GEM foils were stretched during assembly.

\subsubsection{Crosstalk among zigzag strips}
We were able to investigate the crosstalk between zigzag strips with the beam test data due to a minor mistake that was made in the production procedure of the readout pcb, which accidentally connected strips number 63 and 127 in each sector to the ground plane. The large input capacitance presented by these grounded strips to the amplifier in the APV chip causes large noise as can be seen, for example, in the distribution of the pedestal widths for one of the sectors (Fig.~\ref{fig11}). These plots also show that the adjacent strips 62, 64 and 126 are victims of crosstalk as their pedestal widths are slightly higher than the average across all other strips.

We estimate the crosstalk as $\mathrm{\sqrt{rms_{victim}^{2}-rms_{avg}^{2}}/rms_{aggressor}}$, where $\mathrm{rms_{avg}}$ is the average pedestal width observed for strips not affected by crosstalk, $\mathrm{rms_{aggressor}}$ is the pedestal width on the noise source strips (63 or 127), and $\mathrm{rms_{victim}}$ is the pedestal width on a strip adjacent to the noise source strip. We find an average crosstalk among the zigzag strips of ($5.5~\pm~0.2$)\% with an rms width of 1.3\% (Fig.~\ref{figOverallCrossTalk}).

\section{Angular resolution studies}\label{resolution}
\subsection{Multiple Coulomb scattering in detector materials during beam test}
In the beam test, the total radiation length of the detector materials in the setup composed of ten GEM detectors was about 14\% (Table~\ref{mytable}), which impacts the spatial resolution analysis due to multiple Coulomb scattering (MCS) of the tracks in the material. For example, the rms width of the scattering angle distribution due to MCS in that amount of materal is estimated to be 147~$\mu$rad for 32~GeV/c hadrons using the standard MCS formula~\cite{PDG_multipleScattering}.

\begin{table}[h]
\caption[Table caption text]{Material estimate for the tracking detectors in the FNAL beam test.}
\begin{tabular}{|r||c|c|c|c|}
  \hline
  Detector & Gas gaps & Window mat./     & Readout mat./       & Rad. Len.  \\
           &     [mm] &   thickness [mm] & thickness [mm]                & [\%$X_0$] \\ \hline \hline
  Tracker 1 & 3/2/2/2 & Mylar/$\sim$~0.1 & G10/kapton/honeycomb & 0.32 \\ \hline
  Tracker 2 & 3/2/2/2 & Mylar/$\sim$~0.1  & G10/kapton/honeycomb & 0.32 \\ \hline
  SBS 1     & 3/2/2/2 & Al+kapton         & G10/kapton/honeycomb & 0.345 \\ \hline
  UVA-1m-GEM & 3/2/2/2       & Mylar/$\sim$~0.1 & G10/kapton/Rohacell foam & 0.42 \\ \hline
  FIT-1m-zigzag-GEM & 3/1/2/1 & PCB/3.175        & G10/3.175                       & 3.88 \\ \hline
 FIT-30cm & 3/2/2/2 & PCB/3.175      & G10/2.362                        &3.42 \\ \hline
 FIT-10cm-1&3/2/2/2& Mylar/$\sim$~0.1 & G10/2.362                   & 1.5 \\ \hline
 FIT-10cm-2&3/2/2/2& Honycomb/3.175&G10/2.362                   & 1.48 \\ \hline
 Tracker 3 & 3/2/2/2 & Al+kapton     & G10/kapton/honeycomb   & 0.345 \\ \hline
 Tracker 4 & 3/2/2/2 & Mylar/$\sim$~0.1 & G10/kapton/honeycomb& 0.32 \\ \hline
 Ar/CO$_{2}$ & 88 mm &                 &                                      & $\sim$~0.66 \\ \hline
 Air       & $\sim$~3 m &                  &                                        & $\sim$~1 \\ \hline \hline
 Total Material        &                  &       &                                 &  \textbf{14\% $X_{0}$}    \\
  \hline
\end{tabular}
\label{mytable}
\end{table}

A stand-alone Geant4~\cite{G4Reference} simulation was created to study the impact of MCS and to extract the intrinsic detector resolution precisely. The beam test setup as simulated in Geant4 is shown in Fig.~\ref{fig13}. The material distribution from Table~\ref{mytable} is fully implemented in the simulation, while details such as GEM holes and readout strip geometries are not. For the Geant4 physics list, FTFP\_BERT is used, which includes the MCS model based on Lewis theory. Perpendicular point-like beams start 20 mm in front of the first tracker detector (REF1) at position (x,y) = (0 mm, 0 mm) and different beam momenta and particles are simulated.

We estimate the magnitude of the MCS effect in each detector by running the simulation with perfect intrinsic detector resolutions, but with MCS turned on. The simulated hit positions are histogrammed for each reference tracker detector and for the zigzag GEM detector. In the simulation, the zigzag GEM can be treated as a 2D detector in a Cartesian system. In Fig.~\ref{fig14}, the top plot shows the hit position distributions in the horizontal X direction for 25 GeV/c pions, while the bottom plot is for 120 GeV/c protons. Beams are widened significantly at lower momenta due to MCS and the hit distribution widens as the beam particles move downstream. The scattering angles between detectors can be calculated (Fig.~\ref{ScatterAngle}). The scattering is small at the beginning of the beam and it becomes large at the end of the 3-m long tracking system; the overall mean scattering angle between tracker detectors one and four is about 97~$\mu$rad with a 66~$\mu$rad rms. The exclusive residual width for the zigzag GEM detector due to the MCS effect alone is found to be about 80~$\mu$m in Cartesian coordinates~(Fig.~\ref{G4ResFITEICX_Ex}) when simulating 25~GeV/c pions with all tracker detectors set to perfect intrinsic resolution. Here `exclusive' means that the hit of the zigzag detector is excluded from the track fit.

\subsection{Alignment of reference trackers}

The first step in the resolution measurement is an alignment of the four small tracking detectors that have a Cartesian X-Y strip readout. The trackers are first aligned to each other in Cartesian coordinates. The data sample with the highest statistics is used for this alignment, which was taken with 32~GeV/c mixed hadron beams. Only events with a single cluster observed in each tracker are used for the alignment. Any strip multiplicity is allowed in the clusters including single-strip clusters.

The first alignment step is to shift each of the four tracking detectors iteratively in the XY-plane to make their origins match each other in that plane. The initial shift parameters are mean values from position distributions in X and Y coordinates. In each iteration, straight lines are fitted to the hits in X and Y. Residuals are histogrammed for each detector and the residual distributions are fitted with a double-Gaussian function. Ten percent of the residual mean value of each detector is taken as the shift parameter in the next iteration to avoid overcorrections. The resulting residual mean values converge quickly towards zero after 40 iterations, as can be seen in Fig.~\ref{fig17}. This provides a first coarse alignment. In a second alignment step, we correct also for relative rotations of the tracking detectors around the beam in the XY-plane. We again fit straight lines to the hits in X and Y and iterate through a succession of offsets and rotations around the beam axis relative to the first tracking detector until the residual means from the track fits are very close to zero and the $\chi^2$ of the track fits are minimized. In each iteration, the detectors are first shifted and then rotated; then new residuals and rotation angles are calculated. After about 20 iterations the residual mean values become flat within $\pm$0.4~$\mu$m around zero (Fig.~\ref{fig18}). From this step, we get eight shift parameters and three relative rotation angles as refined alignment parameters for the tracker system. The final step is to optimize the three rotation angles one-by-one. With the shift parameters and two rotation angles kept fixed, the third rotation angle is changed in a small range around the value from the second step with a 1~mrad step. We find parabolic curves of mean $\chi^{2}$ of tracks vs.\ angle and calculate the optimized rotation angle from the minimum of parabolic fits to those curves (Fig.~\ref{fig19}). The final rotation angle is taken as the average of these angles obtained for X and Y directions.

\subsection{Spatial resolution of reference trackers}\label{TrackErrMethod}

The Geant4 simulation of pions traversing the setup allows us to estimate track errors and to measure intrinsic detector resolutions when MCS effects are included. The resolutions of the tracker detectors are studied to ensure that their resolutions are good enough to serve as precision reference detectors. In the first simulation step, the intrinsic detector resolutions are first set to zero (perfect resolution) but with MCS turned on, and then the obtained hit positions are additionally smeared by hand using a Gaussian to simulate the intrinsic detector resolution.
We initially assume that the four tracker detectors have the same resolution and smear the simulated hit data with resolutions from 50-80 $\mu$m in 5 $\mu$m steps. For each smeared resolution, we calculate the exclusive track-hit residuals for each detector and get the residual widths, so we can compare exclusive residual widths from simulation with those observed in experimental data (Fig.~\ref{TrackersRes}). When the residual width matches for a tracker detector, the corresponding input resolution used for the Gaussian smearing is taken as the intrinsic resolution of that detector. We do this in both X and Y coordinates and the average of the two is taken as the final resolution for a tracker detector. The resolutions of the tracking detectors 1-4 are found to be 73, 70, 59, and 68 $\mu$m, respectively. The statistical uncertainties of these resolutions are smaller than 2.5~$\mu$m (half of the step size). These measured resolutions fall into a typical range for standard GEM detectors with 400 $\mu$m strip pitch.

\subsection{Alignment of GEM detector with zigzag strips}

Since the zigzag strips in the large GEM detector run in radial direction and measure only the azimuthal $\phi$-coordinate, we need to perform the resolution study for this detector in polar coordinates. The tracker hit positions are transferred to the polar coordinate system that is naturally given by the geometry of the large zigzag detector with the vertex of the trapezoid used as the origin of this polar system (Fig.~\ref{fig21}), and the tracker tracks are then refit in these polar coordinates.

In the polar coordinate system, the $\phi$-coordinate of the center of each zigzag strip in the large trapezoidal detector is simply given by $\phi_{\mathrm{n}} = -0.5\cdot\alpha~+~n \cdot \alpha/(128-1)$, where $\alpha$~=~10$^o$ is the opening angle of the trapezoidal shape and n~=~0,1,$\dots$,127 is the strip number. The $\phi$-position of a hit is initially determined from the barycenter, or centroid, of the cluster using the strip charges as weights.

For transforming the Cartesian tracker coordinates into the polar coordinate system, we need to find the correct X$\mathrm{_{offset}}$ and Y$\mathrm{_{offset}}$ from the new origin to the center of the tracker (Fig.~\ref{fig21}) using the tracks in the $\phi$-coordinate. This has to be done for each test scenario. Note that since the tracker is already internally aligned, the centers of the individual trackers already match, so there is only a single pair (X$\mathrm{_{offset}}$, Y$\mathrm{_{offset}}$) of offsets to be found for each scenario. The idea behind our alignment procedure is that any misalignment will shift the residual means away from zero and will increase the residual widths when fitting tracks in the $\phi$-coordinate with hits in the trapezoidal GEM.

We use a two-fold iteration loop to find the X and Y offsets for each data-taking scenario. First, we keep a fixed X$\mathrm{_{offset}}$ value and step through Y$\mathrm{_{offset}}$ values with a fixed step of 0.2~mm within a reasonable physical range, then change to another X$\mathrm{_{offset}}$ value with a step of 1~mm and vary the Y$\mathrm{_{offset}}$ again until all X$\mathrm{_{offset}}$ values have been covered. For each (X$\mathrm{_{offset}}$, Y$\mathrm{_{offset}}$) pair, tracks are fitted in the $\phi$-coordinate including and excluding the $\phi$-hit in the zigzag detector and the corresponding inclusive and exclusive residuals for the zigzag GEM detector are histogrammed. We record residual means and widths, as well as the mean $\chi^2$ of the track fits in the $\phi$-coordinate.

We use a three-step procedure to find the best estimate of the offsets. Fig.~\ref{fig22} shows an example of this procedure for the data recorded in $\eta$-sector five of the zigzag GEM at 3300~V. For a given X$\mathrm{_{offset}}$ value we record that Y$\mathrm{_{offset}}$ value which produces a well-centered $\phi$-residual for the zigzag detector with a mean of zero (Fig.~\ref{fig22}, left). Then, for a given Y$\mathrm{_{offset}}$, we plot the residual widths as function of X$\mathrm{_{offset}}$ and find the X$\mathrm{_{offset}}$ value that minimizes the residual width using a parabolic fit (Fig.~\ref{fig22}, center). This yields two sets of (X$\mathrm{_{offset}}$, Y$\mathrm{_{offset}}$) pairs that are plotted as two curves. The (X$\mathrm{_{offset}}$, Y$\mathrm{_{offset}}$) values for which those two curves intersect is taken as the best estimate of the alignment offsets (Fig.~\ref{fig22}, right). Best (X$\mathrm{_{offset}}$, Y$\mathrm{_{offset}}$) pairs based on track-$\chi^2$ are also plotted and confirm the results.

Finally, it is not guaranteed that tracker and zigzag detector are installed without a relative rotation in the X-Y plane (Fig.~\ref{fig23}, top). Since aligning the zigzag GEM detector also requires matching the X direction to the trackers so that all calculations of $\phi$ are correct, the rotation angle of the zigzag GEM detector needs to be checked. The mean $\chi^2$ of tracks in $\phi$ vs.\ rotation angles is a parabola and the angle can be calculated from the minimum; it turns out to be very close to zero (Fig.~\ref{fig23}, bottom).

\subsection{Corrections for non-linear strip response of the GEM detector with zigzag strips}

Before proceeding to the measurement of the angular resolution for the zigzag GEM detector, we correct its hit data. The reason to do so is that there is a non-linear response of the zigzag strips when the hit positions are calculated from the centroids of clusters. Specifically, we use track information to apply the non-linear corrections to the hit positions. For each  cluster with a strip multiplicity N~$>$~1, we define a quantity $\eta = \mathrm{s_c-s_{max}}$, where $\mathrm{s_c} = \Sigma_{i=1}^{n}\mathrm{q_i\cdot s_i} / \Sigma_{i=1}^{n} \mathrm{q_i}$ is the centroid position of the cluster in terms of strip number; $\mathrm{s_i}$ and $\mathrm{q_i}$ are the number and the induced charge (in ADC counts) for the i$^\mathrm{th}$ strip in the cluster and $\mathrm{s_{max}}$ is the number of the strip in the cluster on which the maximum charge is induced. The quantity $\eta$ is then a measure of the difference between centroid position and the maximum-charge strip independent of the absolute strip number. We treat different strip multiplicities separately because $\eta$ has different characteristics for odd and even strip multiplicities. For even multiplicities, $s_c$ tends to be between two strips, whereas for odd multiplicities it is close to the center of the strip with maximum charge. In practice, we only need to correct clusters with strip multiplicities N~=~2 and N~=~3 using corresponding quantities $\eta_2$ and $\eta_3$. Corrections for N~$>$~3 are not necessary since there are very few of such clusters (see also Fig.~\ref{fig8}).

Using reference tracks given by the tracker, exclusive residuals of hits in the zigzag GEM are then plotted vs.\ $\eta$. We observe that the resulting distribution is not flat, which indicates a non-linear response of the zigzag strips to the true hit position. As shown in Fig.~\ref{fig24}, exclusive residuals are plotted vs.\ $\eta_2$ and $\eta_3$ and their profiles are fitted separately to appropriate fit functions. A 10$^{\mathrm{th}}$-degree polynomial produces a good fit for the $\eta_2$ distribution (Fig.~\ref{fig24}, top), while a serpentine function is used for the $\eta_3$ distribution (Fig.~\ref{fig24}, bottom). We then globally correct the orginal hit positions by subtracting the $\eta$-dependent offsets obtained from the fit functions from the original hit positions. Fig.~\ref{fig25} shows that the exclusive residual vs.\ $\eta$ plot is much flatter after this correction. Consequently, the overall residual distribution becomes significantly narrower and the spatial resolution is improved when the corrected hit positions are used. Note that the tracker detectors do not need to be corrected since their resolutions are already good enough as they have straight strips with much smaller pitch and consequently much more linear responses. For the position scan data, we apply individual correction functions for each sector since the strip width changes along the radius. For the HV scan data, which were taken in a single sector, the correction functions are obtained from all HV scan data combined together and applied equally to different voltage points. We have checked that the difference is very small if we find correction functions using only data taken at one voltage point.

\subsection{Angular resolution measurement for GEM detector with zigzag strips}
In polar coordinates, we smear the simulated tracker hits with the realistic intrinsic tracker resolutions obtained above and with MCS fully taken into account. We fit a linear track to the azimuthal coordinates of the four smeared hits from the tracker detectors and simulate an unsmeared hit for the zigzag GEM detector. The resulting width of the exclusive residual $\phi_{\mathrm{interpol. track}} - \phi_{\mathrm{unsmear.zigzag}}$ distribution for the zigzag GEM is then a measure of the interpolated-track error (IE) at the position of the zigzag GEM with MCS taken into account. For the experimental data, we also fit linear tracks to the $\phi$-coordinate and calculate the exclusive residual width (ER) for the zigzag GEM detector. The intrinsic resolution $\sigma$ for the zigzag GEM can then be calculated by subtracting the two quantities in quadrature: $\sigma = \mathrm{\sqrt{ER^2-IE^2}}$. We check that this method gives accurate results using the simulation by also smearing the resolution of the probed zigzag GEM detector from 10 to 390 $\mu$rad in 10 $\mu$rad steps and by then applying the same method. We find almost perfect agreement between input resolution from Gaussian smearing and calculated intrinsic resolution $\sigma$ over the full range (Fig.~\ref{figG4SimFITEIC}).

Fig.~\ref{fig27} shows the resulting measured intrinsic angular resolutions for the zigzag GEM detector in central sector five at different $\mathrm{V_{drift}}$ values. The resolution is around 180~$\mu$rad on the efficiency plateau for 2-strip and 3-strip clusters when the non-linear responses of zigzag strips are corrected. In the bottom right plot, we see that the overall angular resolution for the zigzag GEM including single-strip clusters, which cannot be corrected, is $\sigma =$~(193~$\pm$~3)~$\mu$rad at highest tested voltage. This value corresponds to 14\% of the angular pitch of the radial strips. The resolutions are also measured for different positions on the zigzag GEM operated at 3200~V (Fig.~\ref{fig28}).

\section{Summary and Conclusion}
The large-area zigzag GEM detector performed very well in the Fermilab beam test. The $\sim$5.5\% average crosstalk between the zigzag strips is small and does not show significant impact on signal collection. The detection efficiency is above 98\% for charged particles as expected for a GEM detector. The zigzag GEM detector achieves an angular spatial resolution of 193~$\mu$rad or $\sim$14\% of the angular strip pitch after the centroid positions of clusters are corrected for non-linear strip response.

In conclusion, the zigzag readout design is considerably more cost-effective for large-area GEM detectors than conventional straight strip readout structures while preserving good performance. Consequently, GEM detectors with zigzag readout strips are a viable option for application as affordable tracking detectors at future colliders such as the Electron-Ion Collider~(EIC).

\section*{Acknowledgments}
This work is supported by Brookhaven National Laboratory under the EIC eRD-6 consortium. We acknowledge the Fermilab Test Beam Facility staff for their great assistance during the beam test. We are grateful to Alexander Kiselev (Brookhaven National Lab) and Kondo Gnanvo (University of Virginia) for helpful discussions and suggestions on estimating track errors for the resolution studies. Finally, we also would like to thank Bob Azmoun (Brookhaven National Lab) for providing the double-Gaussian fit script and Thomas Hemmick (Stony Brook University) for the discussion on implementing the non-linear hit corrections.

\section*{References}