Transformation models for ROC analysis

To fix ideas, we first consider the case of two samples without covariates. Let the shift term take the form $\mu_{d}(\text{\boldmath$x$})=\delta d$. The CDF of the test results in the nondiseased population is given by $F_{0}(y)=F_{Z}(h(y))$ and in the diseased population by $F_{1}(y)=F_{Z}(h(y)-\delta)$. Using Equation 3, the induced ROC curve can be expressed as

$\displaystyle\text{ROC}(p)$

$\displaystyle=1-F_{Z}\left(F_{Z}^{-1}(1-p)-\delta\right).$

(4)

The model assumption implies that a monotone function $h$ exists to transform both $Y_{1}$ and $Y_{0}$ into the same distribution, $Z\sim F_{Z}$ separated by a shift parameter, $\delta$. The induced ROC curve from this model doesn’t assume a particular distribution of the test result, rather, it quantifies the difference between the test result distributions on the scale of a user-defined $F_{Z}$. In this sense, the difference between the test result distributions is described by $\delta$. Each choice of $F_{Z}$ leads to a different interpretation of $\delta$. For example, when $F_{Z}$ is selected to be the standard normal distribution function, $\delta$ is the difference in means of the transformed test results comparing the diseased and nondiseased groups, $\E[h(Y_{1})-h(Y_{0})]$. Similarly, when $F_{Z}$ is the standard logistic distribution function, $\exp(\delta)$ is the ratio of odds of having a positive test result comparing diseased and nondiseased groups. Closed-form expressions can be derived for summary indices of the ROC curve by solving the appropriate integrals. The expressions of AUC, $J$, the optimal threshold $c^{*}$, sensitivity and specificity at $c^{*}$ are given for some well known choices of $F_{Z}$ in Table 1.

The accuracy of a diagnostic test may be influenced by a set of covariates $X$. To evaluate this effect on the ROC curve and its summary indices, we assume a linear transformation model with a shift term that takes the form

$\displaystyle\mu_{d}(\text{\boldmath$x$})=\delta d+\text{\boldmath$x$}^{\top}\text{\boldmath$\xi$}+d\text{\boldmath$x$}^{\top}\text{\boldmath$\gamma$},$

(5)

where $\text{\boldmath$\xi$},\text{\boldmath$\gamma$}\in\mathbb{R}^{P}$ are the coefficient vectors for the covariates and interaction term, respectively. Under this model, the resulting covariate-specific ROC curve is

$\displaystyle\text{ROC}(p\mid\text{\boldmath$x$})=1-F_{Z}\left(F_{Z}^{-1}(1-p)-(\delta+\text{\boldmath$x$}^{\top}\text{\boldmath$\gamma$})\right),$

where the covariate effect on the ROC curve is given by the difference in shift terms between diseased and nondiseased subjects, $\delta(\text{\boldmath$x$})=\delta+\text{\boldmath$x$}^{\top}\text{\boldmath$\gamma$}$. Similarly, the covariate-specific AUC is given by

$\displaystyle\text{AUC}(\text{\boldmath$x$})=\Prob(Y_{0}<Y_{1}\mid\text{\boldmath$X$}=\text{\boldmath$x$})=a\left(\delta(\text{\boldmath$x$})\right)=a\left(\delta+\text{\boldmath$x$}^{\top}\text{\boldmath$\gamma$}\right),$

(6)

where $a:\mathbb{R}\mapsto[0,1]$ is the AUC function from the first row of Table 1 for different choices of $F_{Z}$. The expressions for $J$, $c^{*}$, sensitivity and specificity can analogously be adjusted to account for covariates, with $\delta$ replaced by $\delta+\text{\boldmath$x$}^{\top}\text{\boldmath$\gamma$}$ in Table 1. Interpretation of the interaction coefficient is similar to Pepe (1998). For example, in the case of a continuous covariate $X=x\in\mathbb{R}$, for each possible specificity value $1-p\in(0,1)$, a unit increase in $x$ results in a $\gamma$-unit increase in the ROC curve (or an increase in the sensitivity) on the scale of $F_{Z}$. If $\gamma$ is positive, an increase in $x$ corresponds to an increase in the ROC curve, indicating that a test is better able to discriminate the two populations for larger values of $x$ and, vice versa. Note that the ROC curve varies with the covariate contingent upon the presence of an interaction between $d$ and $x$ (Pepe, 2003). For $\gamma=0$, the covariate affects the distribution of the test results from the diseased and nondiseased population, but not the ROC curve. That is, for all levels of $x$, the difference between the transformed distributions $h(Y_{1})$ and $h(Y_{0})$ is given by $\delta$, thus the ROC curve is unchanged. Analogous interpretations hold when we are interested in modeling a set of covariates $X$, which could possibly include categorical covariates.

Standard regression techniques have also been proposed as an alternative to assess the effect of covariates on summary indices rather than deriving the induced ROC curve. For example, Dodd and Pepe (2003) model the partial AUC as a regression function of covariates. Our model equivalently results in a regression model for the AUC when the shift term is given from Equation 5 and a particularly chosen scale $F_{Z}$. Specifically, note that Equation 6 is the model proposed by Dodd and Pepe (2003), where $\delta+\text{\boldmath$x$}^{\top}\text{\boldmath$\gamma$}$ is in the form of a usual linear predictor and $a$ is a monotonically increasing inverse link function which defines the scale for the regression coefficients. This relationship provides another interpretation of the interaction term. For example, when $F_{Z}=\operatorname{cloglog}^{-1}$ or $F_{Z}=\operatorname{loglog}^{-1}$, $\exp(\gamma)$ is a ratio of AUC odds comparing a subject with covariate value $x+1$ to $x$. The AUC odds conditional on $x$ are defined as ${\Prob(Y_{1}>Y_{0}\mid X=x)}/{\Prob(Y_{1}\leq Y_{0}\mid X=x)}$. As will be shown in Section 2.2, an advantage of our approach is that we estimate the regression parameters of the transformation model using maximum likelihood estimation and do not rely on less efficient binary regression techniques. Thus, we retain parameter interpretation while allowing for efficient estimation of a variety of covariate-specific summary indices of the ROC curve.

We show that our method is additionally related to the probabilistic index model (PIM) of Thas et al. (2012); De Neve et al. (2019), who propose to model the AUC as a function of $\text{\boldmath$X$}^{*}-\text{\boldmath$X$}$. Here, the superscript (*) indicates an alternate configuration of the modeled covariates. We discuss the PIM of the form

$\displaystyle\Prob(Y<Y^{*}\mid D,D^{*},X,X^{*})=a\left(\text{$\tilde{\delta}$}(D^{*}-D)+(\text{\boldmath$X$}^{*}-\text{\boldmath$X$})^{\top}\text{\boldmath$\tilde{\xi}$}+(D^{*}\text{\boldmath$X$}^{*}-D\text{\boldmath$X$})^{\top}\text{\boldmath$\tilde{\gamma}$}\right)$

where $a$ is the function as defined previously and the tilde ($\texttildelow$) is used to distinguish the parameters of the PIM from those of the transformation model. Contrasting a diseased and nondiseased subject with the same covariate value, we let $D^{*}=1$, $D=0$, $\text{\boldmath$X$}^{*}=\text{\boldmath$X$}$ and $a=\operatorname{logit}^{-1}$, then the probabilistic index is given by $\operatorname{expit}(\text{$\tilde{\delta}$}+\text{\boldmath$X$}^{\top}\text{\boldmath$\tilde{\gamma}$})$. Thus, this is equivalent to the AUC resulting from a linear transformation model with $F_{Z}=\operatorname{cloglog}^{-1}$ and a shift term as in Equation 5. Similarly, when $F_{Z}=\operatorname{probit}^{-1}$, we have the following relationship between the parameters of the linear transformation model and the PIM, $\text{$\tilde{\delta}$}=\delta/\sqrt{2}$ and $\text{\boldmath$\tilde{\gamma}$}=\text{\boldmath$\gamma$}/\sqrt{2}$ . As was the case for ROC curves modeled by transformation models, the interaction between $D$ and $X$ allows the AUC to vary with covariates. The main effect $\tilde{\xi}$, however, does not effect the AUC for configurations where the interest is to compare the distributions of $Y_{1}$ and $Y_{0}$ among subjects with a specific value of covariates.

We can also consider more general and potentially nonlinear formulations of the shift and scale terms in our framework. For the special case of $F_{Z}=\operatorname{probit}^{-1}$, the AUC from a shift-scale transformation model (Siegfried et al., 2022) is given by

$\displaystyle\Prob(Y_{0}<Y_{1}\mid\text{\boldmath$X$}_{0}=\text{\boldmath$x$}_{0},\text{\boldmath$X$}_{1}=\text{\boldmath$x$}_{1})=\Phi\left(\frac{\mu_{1}(\text{\boldmath$x$}_{1})-\ \mu_{0}(\text{\boldmath$x$}_{0})}{\displaystyle\sqrt{\sigma_{0}(\text{\boldmath$x$}_{0})^{2}+\sigma_{1}(\text{\boldmath$x$}_{1})^{2}}}\right)$

where $\text{\boldmath$X$}_{0}$ and $\text{\boldmath$X$}_{1}$ are the corresponding (potentially different) sets of covariates in the nondiseased and diseased populations, respectively. When the scale term depends only on a single set of covariates, $\sigma_{0}(\text{\boldmath$x$}_{0})=\sigma_{1}(\text{\boldmath$x$}_{1})=\sigma(\text{\boldmath$x$})$, despite varying sets of covariates in the shift terms, all the expressions hold from Table 1 with $\delta$ replaced by $\frac{\mu_{1}(\text{\boldmath$x$}_{1})-\mu_{0}(\text{\boldmath$x$}_{0})}{\sigma(\text{\boldmath$x$})}$. However, such closed-form expressions cannot be derived for other choices of $F_{Z}$ when the scale term depends on the disease indicator or on different sets of covariates. In such cases, AUCs and other summary indices can be derived using numeric techniques on the induced ROC curve.

If two (or more) diagnostic tests are measured on the same set of subjects (paired design), the test results are likely to be correlated leading to paired or correlated ROC curves (Swets and Pickett, 1982). To compare the accuracy of these tests, an approach that is frequently used in practice is to consider a summary index of the ROC curve and evaluate a hypothesis for the equivalence of the indices. For example, DeLong et al. (1988) developed an asymptotically normal nonparametric statistic to compare correlated AUCs based on the theory of U-statistics; see also Hanley and McNeil (1983); Wieand et al. (1989); Thompson and Zucchini (1989); McClish (1989); Jiang et al. (1996); Bandos et al. (2005, 2006) who provide similar methods. The problem with using the AUC as a summary index is that diagnostic tests with different ROC curves could potentially have the same value for the AUC. To overcome this limitation, we construct a distribution-free test for differences in the marginal ROC curve parameters using multivariate transformation models that inherently account for the correlations between outcomes in their dependency structure. The null hypothesis represents the equality of entire ROC curves instead of comparing correlated summary indices. Additionally, our test allows for covariate adjusted comparisons to be made.

Let the continuous random vector $\text{\boldmath$Y$}=(Y_{1}$, $Y_{2},\dots,Y_{J})$ denote the results of $J$ diagnostic tests, possibly dependent on a set of covariates $X$. This notation should not be confused with subscripts denoting quantities conditional on $D=1$. To account for the correlation between the test results, we consider a multivariate transformation model with the strictly monotonically increasing conditional transformation function $g:\mathbb{R}^{J}\mapsto\mathbb{R}^{J}$. This function maps the unknown conditional distribution of $Y$ depending on the disease indicator and covariates, to a random vector with a known distribution $\text{\boldmath$Z$}=g(\text{\boldmath$Y$}\mid d,\text{\boldmath$x$})$. Specifically, the vector $\text{\boldmath$Z$}=(Z_{1},\dots,Z_{J})$ follows a zero mean multivariate normal distribution $\text{\boldmath$Z$}\sim N_{J}(\mathbf{0},\mathbf{\Sigma})$ with covariance matrix $\mathbf{\Sigma}\in\mathbb{R}^{J\times J}$. Consequently, the marginal distribution of a multivariate normal vector is itself normally distributed with $Z_{j}\sim N(0,\varsigma^{2}_{j})$, where the diagonal elements of $\mathbf{\Sigma}$ give the variances $\varsigma_{j}^{2}$ for $j=1,\dots,J$. Thus, the conditional joint CDF of $Y$ is given by

	$\displaystyle\Prob(\text{\boldmath$Y$}\leq\text{\boldmath$y$}\mid D=d,\text{\boldmath$X$}=\text{\boldmath$x$})$	$\displaystyle=\Prob(g(\text{\boldmath$Y$}\mid d,\text{\boldmath$x$})\leq g(\text{\boldmath$y$}\mid d,\text{\boldmath$x$}))$
		$\displaystyle=\Prob(\text{\boldmath$Z$}\leq g(\text{\boldmath$y$}\mid d,\text{\boldmath$x$}))$
		$\displaystyle={\Phi}_{\mathbf{0},\mathbf{\Sigma}}\left(g_{1}(y_{1}\mid d,\text{\boldmath$x$}),\dots,g_{J}(y_{J}\mid d,\text{\boldmath$x$})\right)$

where $\Phi_{\mathbf{0},\mathbf{\Sigma}}$ is the joint CDF of a multivariate normal distribution with a zero mean vector and covariance matrix $\mathbf{\Sigma}$. Letting

$\displaystyle g_{j}(y_{j}\mid d,\text{\boldmath$x$})$

$\displaystyle=\Phi^{-1}_{0,\varsigma_{j}^{2}}\left(F_{Z}(h_{j}(y_{j})-\mu_{d}(\text{\boldmath$x$}))\right)$

where $\Phi_{0,\varsigma^{2}}$ is the CDF of a normal distribution with variance $\varsigma^{2}$, the CDF of $Y$ coincides to the Gaussian copula and the marginal distribution of $Y_{j}$ is given by

	$\displaystyle\Prob(Y_{j}\leq y_{j}\mid D=d,\text{\boldmath$X$}=\text{\boldmath$x$})$	$\displaystyle={\Phi}_{\mathbf{0},\mathbf{\Sigma}}\left(\infty,\dots,\infty,g_{j}(y_{j}\mid d,\text{\boldmath$x$}),\infty,\dots,\infty\right)$
		$\displaystyle=\Phi_{0,\varsigma_{j}^{2}}\left(g_{j}(y_{j}\mid d,\text{\boldmath$x$})\right)$
		$\displaystyle=F_{Z}\left(h_{j}(y_{j})-\mu_{d,j}(\text{\boldmath$x$})\right).$

By using the multivariate model specified we allow $Y$ to have a nonlinear dependence structure through the marginal transformation functions $h_{j}$. Analogous to our univariate transformation model, this multivariate model implies that the $j$th marginal transformed test results follow

$\displaystyle h_{j}(Y_{j})=\mu_{d,j}(\text{\boldmath$x$})+Z$

where the function $\mu_{d,j}$ is the shift term for the $j$th test result model and $Z\sim F_{Z}$, as described previously. The marginal ROC curve is then given by

$\displaystyle\text{ROC}_{j}(p\mid\text{\boldmath$x$})=1-F_{Z}\left(F_{Z}^{-1}(1-p)-\delta_{j}(\text{\boldmath$x$})\right)$

where $\delta_{j}(\text{\boldmath$x$})=\mu_{1,j}(\text{\boldmath$x$})-\mu_{0,j}(\text{\boldmath$x$})$ and values of marginal summary indices can be computed from the closed-form expressions given in Table 1 using the marginal ROC curve parameters. Specifically, $\delta$ can be replaced by the covariate effect on the summary indices given by $\delta_{j}(\text{\boldmath$x$})$. From this result, we see that the null hypothesis of equal ROC curves at a given set of covariate values $x$ for the $j$th and $k$th test is equivalent to testing if $\delta_{j}(\text{\boldmath$x$})-\delta_{k}(\text{\boldmath$x$})=0$. Moreover, this hypothesis test is distribution-free as we make no assumptions about $h_{j}$.

Klein et al. (2022) establish conditions which decompose general transformation functions into the marginal ones described above. Specifically, the covariance matrix takes the form $\mathbf{\Sigma}=\mathbf{\Lambda}^{-1}\mathbf{\Lambda}^{-\top}$ where $\mathbf{\Lambda}$ is a lower triangular $(J\times J)$ matrix whose coefficients characterize the dependence structure given by the Gaussian copula. The matrix $\mathbf{\Lambda}$ could also vary between diseased and nondiseased subjects and potentially depend on a set of covariates $x$. This allows the dependence structure of $Y$ to change as a function of $x$. For the special case of two diagnostic tests, the correlation between $Y_{1}$ and $Y_{2}$ for a given set of covariates $x$ is given by $\frac{-\lambda_{d}(\text{\boldmath$x$})}{\sqrt{\lambda_{d}(\text{\boldmath$x$})^{2}+1}}$. For more details and examples of multivariate transformation models, see Klein et al. (2022); Barbanti et al. (2022).

We parameterize the transformation function as

$\displaystyle h(y\mid\text{\boldmath$\vartheta$})=\text{\boldmath$b$}(y)^{\top}\text{\boldmath$\vartheta$}=\sum_{m=0}^{M}\vartheta_{m}b_{m}(y)\quad\text{for}\quad y\in\mathbb{R},$

(7)

where $\text{\boldmath$b$}(y)=(b_{0}(y),\dots,b_{M}(y))^{\top}$ is vector of $M+1$ basis functions with coefficients $\text{\boldmath$\vartheta$}\in\mathbb{R}^{M+1}$. Polynomials in Bernstein form offer a choice of basis that provides a flexible way of estimating the underlying function. The Bernstein basis polynomial of order $M$ is defined on the interval $[l,u]$ as

$\displaystyle b_{m}(y)={M\choose m}\tilde{y}^{m}(1-\tilde{y})^{M-m},\quad m=0,\dots,M,$

(8)

where $\tilde{y}=\frac{y-l}{u-y}\in[0,1]$. The restriction $\vartheta_{m}\leq\vartheta_{m+1}$ for $m=0,\dots,M-1$, guarantees the monotonicity of $h$. Observe that the transformation function is linear in the parameters that define it and, any nonlinearity of the test results is modeled by the basis functions. If the order $M$ is chosen to be sufficiently large, Bernstein polynomials can uniformly approximate any real-valued continuous function on an interval (Farouki, 2012).

Denote the complete parameter vector as $\text{\boldmath$\theta$}=(\text{\boldmath$\beta$}^{\top},\text{\boldmath$\vartheta$}^{\top})^{\top}$, where $\text{\boldmath$\beta$}=(\delta,\text{\boldmath$\xi$}^{\top},\text{\boldmath$\gamma$}^{\top})^{\top}\in\mathbb{R}^{2P+1}$ are the vector of regression coefficients parameterizing the function $\mu_{d}$ from Section 2.1 and $\text{\boldmath$\vartheta$}\in\mathbb{R}^{M+1}$ are the basis coefficients. We follow the maximum likelihood approach proposed by Hothorn et al. (2018) to jointly estimate $\beta$ and $\vartheta$. The advantages of embedding the model in the likelihood framework are as follows. (i) All forms of random censoring (right, left, interval) as well as truncation can directly be incorporated into likelihood contributions defined in terms of the distribution function (Lindsey et al., 1996). (ii) If the given model is correctly specified, under regularity conditions, the maximum likelihood estimator (MLE) will be asymptotically the most efficient estimator (minimum variance and unbiased). (iii) The MLE is asymptotically normally distributed and has a sample variance that can be computed from the inverse of the Fisher information matrix. This can be used to generate confidence intervals for the estimated parameters. (iv) The MLE is equivariant which implies invariance of the score test (or the Lagrange multiplier test) to reparameterizations (Rao, 1948; Aitchison and Silvey, 1958; Silvey, 1959; Dagenais and Dufour, 1991). Specifically, as will be shown in Section 2.4.1, by inverting the score test, our method produces confidence bands for the ROC curve and appropriate score intervals for its summary indices.

The likelihood contribution of a single observation $O=(Y,D,\text{\boldmath$X$})$ where $Y\in(\underline{y},\overline{y}]=\{y\in\mathbb{R}:\underline{y}<y\leq\overline{y}\}$ is given by

$$L(\text{\boldmath$\theta$}\mid O)=\begin{cases}\begin{aligned} &f_{Z}(h(y\mid\text{\boldmath$\vartheta$})-\mu_{d}(\text{\boldmath$x$}\mid\text{\boldmath$\beta$}))h^{\prime}(y\mid\text{\boldmath$\vartheta$})&&y\in\mathbb{R}&\quad&\text{`exact continuous'}\\ &1-F_{Z}(h(\underline{y}\mid\text{\boldmath$\vartheta$})-\mu_{d}(\text{\boldmath$x$}\mid\text{\boldmath$\beta$}))&&y\in(\underline{y},\infty)&\quad&\text{`right censored'}\\ &F_{Z}(h(\overline{y}\mid\text{\boldmath$\vartheta$})-\mu_{d}(\text{\boldmath$x$}\mid\text{\boldmath$\beta$}))&&y\in(-\infty,\overline{y})&\quad&\text{`left censored'}\\ &\begin{split}F_{Z}(h(\overline{y}\mid\text{\boldmath$\vartheta$})-\mu_{d}(\text{\boldmath$x$}\mid\text{\boldmath$\beta$}))-\\ F_{Z}(h(\underline{y}\mid\text{\boldmath$\vartheta$})-\mu_{d}(\text{\boldmath$x$}\mid\text{\boldmath$\beta$}))\end{split}&&y\in(\underline{y},\overline{y}]&\quad&\text{`interval censored'},\\ \end{aligned}\end{cases}$$

where $f_{Z}$ is the density function of $Z$ and $h^{\prime}(y\mid\text{\boldmath$\vartheta$})$ is the first derivative of the transformation function with respect to $y$. Given a sample of $N$ independent and identically distributed observations $O_{i}$ for $i=1,\dots,N$, the log-likelihood is given by $\ell(\text{\boldmath$\theta$})=\sum_{i=1}^{N}\log(L_{i}(\text{\boldmath$\theta$}))$, where $L_{i}$ is the likelihood contribution of observation $i$. The (unconditional) maximum likelihood estimate of $\theta$ is the solution to the optimization problem

$\displaystyle\hat{\text{\boldmath$\theta$}}=(\text{\boldmath$\hat{\beta}$},\text{\boldmath$\hat{\vartheta}$})=\operatorname*{{arg\,max}}_{\text{\boldmath$\beta$},\text{\boldmath$\vartheta$}}\ell(\text{\boldmath$\beta$},\text{\boldmath$\vartheta$}),$

subject to the monotonicity constraint $\vartheta_{m}\leq\vartheta_{m+1}$ for $m=0,\dots,M-1$. The resulting ROC curve only depends on $\beta$ which is decoupled from the parameters needed to model the transformation function $\vartheta$. The score function is defined as the first derivative of the log-likelihood function with respect to each of the parameters and is given by

$\displaystyle S(\text{\boldmath$\theta$})=\left(\begin{array}[]{l}\frac{\partial\ell(\text{\boldmath$\theta$})}{\partial\text{\boldmath$\beta$}}\\ \frac{\partial\ell(\text{\boldmath$\theta$})}{\partial\text{\boldmath$\vartheta$}}\end{array}\right)=\left(\begin{array}[]{l}S_{\text{\boldmath$\beta$}}(\text{\boldmath$\theta$})\\ S_{\text{\boldmath$\vartheta$}}(\text{\boldmath$\theta$})\end{array}\right).$

We perform constrained optimization using the likelihood and score contributions to determine the maximum likelihood estimates for $\beta$ and $\vartheta$. For computational details, see Hothorn (2020). The asymptotic variance of the MLE can further be estimated by the expected Fisher information matrix which is the variance-covariance matrix of the score function and is defined as

$\displaystyle I(\text{\boldmath$\theta$})=-\E\left(\begin{array}[]{ll}\frac{\partial^{2}\ell(\text{\boldmath$\theta$})}{\partial\text{\boldmath$\beta$}\partial\text{\boldmath$\beta$}^{\top}}&\frac{\partial^{2}\ell(\text{\boldmath$\theta$})}{\partial\text{\boldmath$\beta$}\partial\text{\boldmath$\vartheta$}^{\top}}\\ \frac{\partial^{2}\ell(\text{\boldmath$\theta$})}{\partial\text{\boldmath$\vartheta$}\partial\text{\boldmath$\beta$}^{\top}}&\frac{\partial^{2}\ell(\text{\boldmath$\theta$})}{\partial\text{\boldmath$\vartheta$}\partial\text{\boldmath$\vartheta$}^{\top}}\end{array}\right)=\left(\begin{array}[]{ll}I_{\text{\boldmath$\beta$},\text{\boldmath$\beta$}}(\text{\boldmath$\theta$})&I_{\text{\boldmath$\beta$},\text{\boldmath$\vartheta$}}(\text{\boldmath$\theta$})\\ I_{\text{\boldmath$\beta$},\text{\boldmath$\vartheta$}}(\text{\boldmath$\theta$})^{\top}&I_{\text{\boldmath$\vartheta$},\text{\boldmath$\vartheta$}}(\text{\boldmath$\theta$})\end{array}\right).$

The matrix is partitioned such that the submatrix $I_{\text{\boldmath$\beta$},\text{\boldmath$\beta$}}(\text{\boldmath$\theta$})$ corresponds to the parameter related to the disease indicator and covariates. The variation in $\vartheta$ and its relationship to $\beta$ is not of direct relevance here but is nonetheless estimated.

Instrument precision can affect the evaluation of diagnostic biomarkers. For example, when biomarker levels are at or below the limit of detection (LOD) $y_{\text{LOD}}$, the observed value lies in an interval $(-\infty,y_{\text{LOD}})$ and the resulting measurement is left censored. Often a replacement value is substituted for such measurements. Alternatively, only biomarker values above the LOD are used for the ROC analysis. It has been shown that these approaches lead to biased estimation (Lynn, 2001; Singh and Nocerino, 2002). Various adjustments to ROC curves and its summary indices have been proposed to handle such censored measurements (Mumford et al., 2006; Perkins et al., 2007, 2009; Ruopp et al., 2008; Bantis et al., 2017; Xiong et al., 2022). However, these methods typically do not account for covariates. Our framework naturally accounts for such observations in the likelihood function for left censored test results. Similarly, the right censored likelihood accounts for measurements which are affected by an upper limit of detection. Thus, our method provides a smooth covariate-specific ROC curve for all values of specificity with estimates and inference appropriately incorporating the observed information.

Many diagnostic tests of interest to clinical researchers are measured on ordinal scales. These are rank-ordered categorical variables for which quantitative differences between levels are unknown. Examples include the Memory Impairment Screen (MIS) for the measurement of cognition to distinguish Alzheimer’s disease and other dementias (Buschke et al., 1999) or the risk assessment of disease from radiological images such as with the Breast Imaging Reporting and Data System (BI-RADS) (American College of Radiology (2013), ACR). A simple reparameterization of the transformation function facilitates an extension of ROC curves to such ordinal test results.

Let the random variable $Y\in\{y_{1},\dots,y_{K}\}$ denote the results of an ordered categorical diagnostic test with $K$ categories, where $y_{k}<y_{k+1}$ for $k=1,...,K-1$ . Define the covariate-specific ROC curve as a discrete function of the subject categories

$\displaystyle\left\{\left(1-F_{0}(y\mid\text{\boldmath$x$}),1-F_{1}(y\mid\text{\boldmath$x$})\right),y\in\{y_{1},\dots,y_{K}\}\right\}.$

This set of points plots the covariate-specific $1-\text{specificity}$ and sensitivity coordinates for each of the $K$ categories. Suppose that the following transformation model defines the distribution of the ordinal test results in the diseased and nondiseased populations

$\displaystyle F_{d}(y_{k}\mid\text{\boldmath$x$})=F_{Z}\left(h(y_{k})-\mu_{d}(\text{\boldmath$x$})\right),$

where the transformation function assigns one parameter to each of the categories of the test result except the last one, $h(y_{k})=\vartheta_{k}\in\mathbb{R}$ for $k=1,\dots,K-1$. Note that since $\Prob(Y\leq y_{K}\mid\text{\boldmath$X$}=\text{\boldmath$x$})=1$, only $K-1$ parameters need to be specified. The monotonicity constraint $\vartheta_{k}<\vartheta_{k+1}$ is also necessary so that $F_{d}(y_{k}\mid\text{\boldmath$x$})<F_{d}(y_{k+1}\mid\text{\boldmath$x$})$ for $k=1,\dots,K-1$. Observe that this is a special case of the general shift-scale model in Equation 2. Under a linear shift term as given in Equation 5 and $F_{Z}=\operatorname{logit}^{-1}$, this model is known as the proportional odds logistic regression (POLR).

Setting $p_{k}=1-F_{0}(y_{k}\mid\text{\boldmath$x$})$, the ROC curve under the ordinal transformation model is given by

$\displaystyle\text{ROC}(p_{k})=1-F_{Z}\left(F_{Z}^{-1}(1-p_{k})-\delta(\text{\boldmath$x$})\right),\quad k=1,\dots,K.$

The parameters of the model can be estimated by using the maximum likelihood scheme described previously with the interval censored likelihood contribution

$\displaystyle F_{Z}(\vartheta_{k}-\mu_{d}(\text{\boldmath$x$}\mid\text{\boldmath$\beta$}))-F_{Z}(\vartheta_{k-1}-\mu_{d}(\text{\boldmath$x$}\mid\text{\boldmath$\beta$})),$

where the observed subject has a test result of $y_{k}$ and covariate value $x$. After estimating the parameters, it is possible to plot a continuous ROC curve as a function of all $p\in[0,1]$. This curve’s shape depends on the choice of $F_{Z}$ and represents the underlying parametric form of the model. However, it does not represent the ROC curve for a continuous test. Specifically, the function is only defined at the $K$ discrete specificity values. Summary measures can be calculated using the expressions detailed in Table 1. Note that these would be based on a smooth interpolation given by the model rather a discrete one and lead to the interpretation of an ordinal test with a latent continuous scale.

In the two-sample univariate case where $\delta$ defines the ROC curve, as in Equation 4, we can construct score intervals for $\delta$. Unlike the Wald and other commonly used intervals, score intervals are especially desirable as they are invariant to transformations of the parameters. That is, a score CI for $G(\delta)$ (e.g., the AUC $a(\delta)$), provides the same level of coverage as would a score CI for $\delta$. In turn, under a correctly specified model, a score CI for $\delta$ allows the construction of accurately covered uniform confidence bands for the ROC curve as well as intervals for its summary indices such as the AUC and the Youden index.

We first generate score confidence intervals for $\delta$ by inverting the score test. In this case, the null hypothesis is given by $H_{0}:\delta=\delta_{0}$ where $\delta_{0}$ is a specific value of the parameter of interest. Under $H_{0}$, the restricted (conditional) maximum likelihood estimator for $\vartheta$ can be obtained by

$\displaystyle\text{\boldmath$\hat{\vartheta}$}(\delta_{0})=\operatorname*{{arg\,max}}_{\text{\boldmath$\vartheta$}}\ell(\delta_{0},\text{\boldmath$\vartheta$}),$

or as a solution of the $M+1$ score equations $S_{\text{\boldmath$\vartheta$}}(\delta_{0},\text{\boldmath$\vartheta$})=\mathbf{0}$. Note that this estimate is a function of $\delta_{0}$. Letting $\text{\boldmath$\tilde{\theta}$}=(\delta_{0},\text{\boldmath$\hat{\vartheta}$}(\delta_{0}))$, the quadratic (Rao) score statistic simplifies to

	$\displaystyle R(\delta_{0})$	$\displaystyle=S(\text{\boldmath$\tilde{\theta}$})^{\top}I^{-1}(\text{\boldmath$\tilde{\theta}$})S(\text{\boldmath$\tilde{\theta}$})$
		$\displaystyle=(S_{\delta}(\text{\boldmath$\tilde{\theta}$})^{\top},\mathbf{0}^{\top})I^{-1}(\text{\boldmath$\tilde{\theta}$})(S_{\delta}(\text{\boldmath$\tilde{\theta}$})^{\top},\mathbf{0}^{\top})^{\top}$
		$\displaystyle=S_{\delta}(\text{\boldmath$\tilde{\theta}$})^{\top}A_{\delta,\delta}(\text{\boldmath$\tilde{\theta}$})S_{\delta}(\text{\boldmath$\tilde{\theta}$}),$

where $A_{\delta,\delta}(\text{\boldmath$\tilde{\theta}$})$ denotes the submatrix corresponding to $\delta$ of the inverse Fisher information matrix and is given by the Schur complement $I_{\delta,\delta}(\text{\boldmath$\theta$})-I_{\delta,\text{\boldmath$\vartheta$}}(\text{\boldmath$\theta$})I^{-1}_{\text{\boldmath$\vartheta$},\text{\boldmath$\vartheta$}}(\text{\boldmath$\theta$})I_{\delta,\text{\boldmath$\vartheta$}}(\text{\boldmath$\theta$})^{\top}$. Under $H_{0}$, $R(\delta_{0})$ converges asymptotically to a chi-square distribution with 1 degree of freedom, $R(\delta_{0})\stackrel{{\scriptstyle D}}{{\longrightarrow}}\chi^{2}_{1}$. This result is explained by Rao (2005). Thus, inverting the score statistic by enumerating values of $\delta_{0}$ allows for the construction of $(1-\alpha)$ score confidence intervals for $\delta$ defined as

$\displaystyle\{\delta_{0}\in\mathbb{R}\mid R(\delta_{0})<\chi^{2}_{1}(1-\alpha)\},$

where $\chi^{2}_{P}(1-\alpha)$ is the $(1-\alpha)$ quantile value of the $\chi^{2}_{1}$ distribution. Equivalently, we can use the square root of the score statistic to construct a $(1-\alpha)$ score interval using quantiles of the standard normal distribution, $\{\delta_{0}\in\mathbb{R}\mid\Phi^{-1}(\alpha/2)<\sqrt{R(\delta_{0})}\leq\Phi^{-1}(1-\alpha/2)\}$. Finally, we apply the function $G$ to both the lower and upper limits of the interval to construct score confidence bands for the ROC curve or score confidence intervals for its summary indices

The score statistic is given by $R(\delta_{0})=S_{\delta}(\text{\boldmath$\tilde{\theta}$})^{2}A_{\delta,\delta}(\text{\boldmath$\tilde{\theta}$})$. Testing if there is a significant difference between the nondiseased and diseased populations coincides to the hypothesis test, $H_{0}:\delta=0$. This is computationally efficient because only the distribution of $R(0)$ needs to be computed. However, computing score CIs requires updating the restricted MLEs $\text{\boldmath$\hat{\vartheta}$}(\delta_{0})$ for an enumeration of $\delta_{0}$ values. This becomes computationally intractable when enumerating a higher dimensional grid of parameters.

Since the MLE satisfies

$\displaystyle\sqrt{n}(\text{\boldmath$\hat{\beta}$}-\text{\boldmath$\beta$})\xrightarrow{D}N_{P+1}\left(\mathbf{0},A_{\text{\boldmath$\beta$},\text{\boldmath$\beta$}}(\text{\boldmath$\theta$})\right),$

then by the multivariate delta method, $G(\text{\boldmath$\hat{\beta}$})$ also follows a normal distribution with

$\displaystyle\mathbb{V}(G(\text{\boldmath$\hat{\beta}$}))=\frac{1}{n}\nabla G(\text{\boldmath$\beta$})^{\top}A_{\text{\boldmath$\beta$},\text{\boldmath$\beta$}}(\text{\boldmath$\theta$})\nabla G(\text{\boldmath$\beta$}),$

where $\nabla G(\text{\boldmath$\beta$})$ is the gradient of $G$ evaluated at $\beta$ and the inverse Fisher information matrix $A_{\text{\boldmath$\beta$},\text{\boldmath$\beta$}}(\text{\boldmath$\theta$})$ is given by the Schur complement $I_{\text{\boldmath$\beta$},\text{\boldmath$\beta$}}(\text{\boldmath$\theta$})-I_{\text{\boldmath$\beta$},\text{\boldmath$\vartheta$}}(\text{\boldmath$\theta$})I^{-1}_{\text{\boldmath$\vartheta$},\text{\boldmath$\vartheta$}}(\text{\boldmath$\theta$})I_{\text{\boldmath$\beta$},\text{\boldmath$\vartheta$}}(\text{\boldmath$\theta$})^{\top}$. For example, when the shift term takes the linear form as in Equation 5 and $G$ defines the AUC function for $F_{Z}=\operatorname{probit}^{-1}$, the entries of $\nabla G(\text{\boldmath$\beta$})$ are given by

$\displaystyle\frac{\partial G(\text{\boldmath$\beta$})}{\partial\delta}=\frac{1}{\sqrt{2}}C,\quad\frac{\partial G(\text{\boldmath$\beta$})}{\partial\xi_{i}}=0\quad\text{and}\quad\frac{\partial G(\text{\boldmath$\beta$})}{\partial\gamma_{i}}=\frac{x_{i}}{\sqrt{2}}C,$

where $C=\phi\left(\frac{\delta+\text{\boldmath$x$}^{\top}\text{\boldmath$\gamma$}}{\sqrt{2}}\right)$, $\phi$ is the density of the standard normal distribution and $i$ indexes the $P$ covariates. In general, the gradient can be estimated by calculating such derivatives and evaluating the resulting function at the MLE. Similarly, the variance-covariance matrix of the estimated parameters $A_{\text{\boldmath$\beta$},\text{\boldmath$\beta$}}(\text{\boldmath$\theta$})$ can be computed by inverting the numerically evaluated Hessian matrix. Thus, a $(1-\alpha)$ level confidence interval for $G(\text{\boldmath$\beta$})$ is given by

$\displaystyle G(\text{\boldmath$\beta$})\pm\Phi^{-1}(\alpha/2)\sqrt{\hat{\mathbb{V}}(G(\text{\boldmath$\hat{\beta}$}))}.$

When the function $G$ has complex derivatives, as would be the case for nonlinear shift terms $\mu_{d}(\text{\boldmath$x$})$ or when calculating optimal thresholds $c^{*}$ where $G$ includes the inverse of the transformation function, constructing confidence intervals using the delta method becomes infeasible. For these cases, we apply a simple simulation-based algorithm which utilizes the asymptotic normality of the MLE to calculate confidence intervals for the ROC curve and its summary indices, which are functions of the parameters of interest. The steps of the algorithm to construct $(1-\alpha)$ level confidence intervals for $G(\text{\boldmath$\hat{\beta}$})$ can be summarized as follows:

1. 1.

   Generate $B$ independent samples from the asymptotic multivariate normal distribution of the parameter estimates $N_{P+1}\left(\text{\boldmath$\hat{\beta}$},\frac{1}{n}\hat{A}_{\text{\boldmath$\beta$},\text{\boldmath$\beta$}}(\mathbf{\hat{\Theta}})\right)$ and denote as $\text{\boldmath$\hat{\beta}$}^{*}_{1},\dots,\text{\boldmath$\hat{\beta}$}^{*}_{B}$.

2. 2.

   For each sample $b=1,\dots,B$, calculate the function of interest $G(\text{\boldmath$\hat{\beta}$}^{*}_{b})$.

3. 3.

   Construct the confidence interval $\left(Q_{G(\text{\boldmath$\hat{\beta}$}^{*})}(\alpha/2),Q_{G(\text{\boldmath$\hat{\beta}$}^{*})}(1-\alpha/2)\right)$, where $Q_{G(\text{\boldmath$\hat{\beta}$}^{*})}$ is the empirical quantile function of the sample $G(\text{\boldmath$\hat{\beta}$}^{*}_{1}),\dots,G(\text{\boldmath$\hat{\beta}$}^{*}_{B})$.

A similar algorithm is presented in Mandel (2013), who discuss its asymptotic validity and present several examples that show its empirical coverage adheres to nominal levels with results similar to the delta method.