Generalized Roe and Metz receiver operating characteristic model: analytic link between simulated decision scores and empirical AUC variances and covariances

Abstract. Modeling and simulation are often used to understand and investigate random quantities and estimators. In 1997, Roe and Metz introduced a simulation model to validate analysis methods for the popular endpoint in reader studies to evaluate medical imaging devices, the reader-averaged area under the receiver operating characteristic (ROC) curve. Here, we generalize the notation of the model to allow more flexibility in recognition that variances of ROC ratings depend on modality and truth state. We also derive and validate equations for computing population variances and covariances for reader-averaged empirical AUC estimates under the generalized model. The equations are one-dimensional integrals that can be calculated using standard numerical integration techniques. This work provides the theoretical foundation and validation for a Java application called iRoeMetz that can simulate multireader multicase ROC studies and numerically calculate the corresponding variances and covariances of the empirical AUC. The iRoeMetz application and source code can be found at the “iMRMC” project on the google code project hosting site. These results and the application can be used by investigators to investigate ROC endpoints, validate analysis methods, and plan future studies.


Introduction
The area under the receiver operating characteristic (ROC) curve, denoted AUC, is a common endpoint in reader studies to evaluate medical imaging devices. 1 ROC data from reader studies are confidence-of-disease ratings from clinicians (readers) evaluating images (cases). Therefore, the endpoint in a multireader multicase (MRMC) ROC study is affected by two important sources of variability-the readers and the casesand in studies past the exploratory stage, we often want to account for both sources in our analyses.
Modeling and simulation are tools that help us to understand and investigate the distribution and statistical behavior of random quantities and estimators. Roe and Metz 2 (R&M) proposed a simulation model that launched the study of different analysis methods and endpoints related to MRMC ROC studies. Their model was developed to validate the Dorfman, Berbaum, and Metz (DBM) method that compares AUCs from two modalities. 3 The R&M model simulates ROC ratings according to a binormal model for each reader and generates data from a "fully crossed" study design where each patient is imaged by two or more modalities, with the resulting cases evaluated once by each reader.
With few modifications, the R&M model has been used to validate and characterize many other MRMC ROC analysis methods for almost two decades. [4][5][6][7][8][9][10] The R&M model has also been used to investigate power and sizing methods for ROC studies, 11,12 and adapted to yield discrete ROC ratings 13,14 and explore alternative study designs. 15 Finally, while AUC has been the primary reader performance measure analyzed with the R&M model to date, the R&M model has also been used to analyze binary performance measures 16,17 and utility. 18 The R&M model assumes a four-factor modality × reader × case × truth mixed-effects analysis of variance model for the ROC ratings. In the original work, reader and case effects and their interactions were assumed to be random and the variance components were chosen to be the same across truth states and across modalities. As we show below, the assumptions applied to the R&M model in the original work can be relaxed. In the current work, we generalize the notation to clarify that the variance components can depend on truth state and modality.
As discussed by Hillis, 19 when the variance components are assumed to be the same across truth states and modalities, the R&M model has the following interpretations: (1) ROC ratings for each reader are generated from an equal-variance binormal model (i.e, a binormal model such that variances of the nondiseased and diseased ROC ratings are equal); and (2) the expected differences (or separations) between the nondiseased and diseased ROC ratings vary across readers, with the separations having the same variance for each modality. This last result implies that for a simulation study that assumes equal AUCs across modalities (i.e., a null-hypothesis study), the resulting AUC estimates will have the same variance for each modality. It is natural to question this assumption, especially when comparing an imaging modality with and without a computer aid. Beiden et al. 20 found that the reader variability of readers' AUCs was much smaller when using a computer aid in classifying microcalcifications in mammograms compared to those without the aid.
When fitting an ROC curve with a binormal model, it is generally recognized that for real data the distributions of the latent decision variables of the diseased and nondiseased ROC ratings will often have different variances, with the diseased distribution typically wider. This causes the ROC curve to be unsymmetric about the negative diagonal. Such unsymmetric ROC curves have been seen as far back as the early psychophysical experiments of the 1960s 21,22 and in recent studies evaluating medical imaging modalities. 19,[23][24][25][26] Unsymmetric ROC curves have motivated other models of ROC ratings [27][28][29] and are sometimes characterized by a mean-to-sigma ratio defined as the difference of the binormal means divided by the difference of the binormal standard deviations across truth states in an unequal-variance binormal model. 19,21,26 For this reason, Hillis 19 introduced an unequal-variance binormal model by allowing some of the variance components to depend on truth but with some additional constraints.
In this paper, our main purpose is to present the exact nature of the relationship between the R&M model inputs (these include the model parameters and numbers of readers and cases) and the means, variances and covariances of the resulting reader-averaged empirical AUC estimates. R&M note that they were not able to determine such a relationship. Knowledge of this relationship builds on earlier work by Gallas et al. 16 and allows an investigator to, among other things, (1) verify that the simulation model has been correctly programmed by comparing parameter estimates based on the simulations to the true values of the parameters; and (2) quantify the bias of AUC variance estimates by similarly comparing simulation results to the true values. We consider this paper to be an important first step toward our ultimate aim of being able to calibrate a simulation model that will produce data that matches a real data set with respect to the estimated parameters from an analysis method, such as that proposed by Obuchowski and Rockette. 30 In this paper, we begin by generalizing the notation of the R&M model by allowing all of its variance components to depend on both modality and truth state. Then, we present and validate equations for computing the population variances and covariances for empirical AUC estimates computed from data simulated from the generalized R&M model. The generalized model includes the original R&M model (with the equal variance-components assumption) and the unequal-variance model proposed by Hillis 19 as special cases. Although these two special cases are sufficient for many situations, we anticipate that researchers may want to use other special cases of the generalized R&M model, or the generalized model itself for simulating data. Presenting equations for the generalized model eliminates the need to derive equations for each special case in the future.

Generalized Roe and Metz Model
The R&M simulation model is for simulating rating data that emulate an ROC reader-performance study that has N 0 nondiseased cases and N 1 diseased cases that are interpreted and rated by N R readers. Here, we focus on a fully crossed study design to compare two modalities, denoted by A and B, where "fully crossed" refers to the data collection: all the readers rate all the cases in both modalities with respect to confidence of disease. The result of such a study is a dataset with 2 × ðN 0 þ N 1 Þ × N R ROC ratings ½ðmodalitiesÞ × ðcasesÞ × ðreadersÞ. R&M denote the ROC ratings by X ijkt , where i denotes modality (i ¼ "A" or "B"), j denotes reader, k denotes case, and t denotes truth ðt ¼ 0;1 ≡ nondiseased; diseasedÞ.
Using the notation of R&M, the model is given by where X ijkt denotes the value of the ROC rating for modality i, reader j, case k, and truth state t. Modality and truth are fixed factors and reader and case are random factors, i.e., effects involving reader or case are random effects and all other effects are fixed. Consequently, the Greek terms μ t ; τ it are fixed effects and the remaining seven terms are random. The random terms in the R&M model are all independent zero-mean Gaussian random variables: R is the reader effect, C is the case effect, ½RC is the reader × case effect, ½τR is the modality × reader effect, ½τC is the modality × case effect, ½τRC is the modality × reader × case effect, and E is an independent random error term. The corresponding variance components are denoted σ 2 R ; σ 2 C ; σ 2 RC ; σ 2 τR ; σ 2 τC ; σ 2 τRC , and σ 2 E . The independent error term E ijkt can be attributed to a reader's inability to exactly reproduce their ROC rating for a case; E is sometimes referred to as internal noise or reader jitter. R&M pointed out that the pure error and the three-way interaction variance components cannot be separately estimated "without multiple readings of each case in each modality by all readers." As such, they defined an aggregate term ε ijkt ¼ ½τRC ijkt þ E ijkt .

Updated notation
Here, we generalize the R&M model, clarifying that the model is, in fact, a four-factor model. In addition to modality, reader, and case, truth is a factor. In particular, truth has a hierarchical relationship with cases: 31 each case can have only one truth state and cases are "nested" within truth. Mathematically, we will use α to express the truth factor and rewrite the R&M model as (2) where we have omitted the pure error term since it cannot be distinguished from the four-way interaction term without replications. Following statistical conventions, we could surround α and its subscript t with parentheses (for terms that include case) to indicate the nesting of cases within truth. However, we will not follow this convention for simplicity. Now, the four-factor structure of the model is clearly shown in Eq. (2). We point out that this model includes only interactions with truth or effects nested within truth; in particular, there are no effects for modality alone or reader alone. The rationale for omitting the modality and reader effects is that these terms would have no effect on the ROC curve for a given reader and test, since the ROC curve is invariant to location shifts of the decision variable. We note that the fixed interaction ½τα it allows different modalities to have different ROC curves.

Allow variances to depend on modality and truth
Given the new notation for the R&M model, we now generalize it to allow the variance components to depend on modality and truth.
There are three random-effect terms that do not include modality: ½Rα jt ; ½Cα kt ; ½RCα jkt . They correspond to six variance components: σ 2 R0 ; σ 2 C0 ; σ 2 RC0 are for nondiseased cases ðt ¼ 0Þ and σ 2 R1 ; σ 2 C1 ; σ 2 RC1 are for diseased cases ðt ¼ 1Þ. The sum of these six variance components that do not depend on modality is There are three random-effect terms that include modality: ½τRα ijt , ½τCα ikt , and ½τRCα ijkt . They correspond to 12 variance components: The sums of the variance components that are specific to modality A and B are The original R&M model is a special case of the generalized model. All we need to do is assume, as R&M did, that the variance components of the ROC ratings do not depend on modality or truth. We can replicate the original R&M model by setting generalized R&M model variance components equal to the original ones as given in Table 1. Other simplifications can be similarly handled.

Expected AUCs
Here, we examine the expected value of the empirical estimate of AUC, also known as the trapezoidal estimate, 32 given the R&M model. Without loss of generality, we focus on modality A. A similar discussion can be derived for modality B.
The estimated reader-averaged AUC for modality A is where we shall refer to sðxÞ as the "success function;" sðxÞ equals 1.0 when reader j successfully rates diseased case k 0 higher than nondiseased case k, sðxÞ equals 0.0 if the ratings are in the wrong order, and sðxÞ equals 0.5 if the ratings are tied.
As we consider the expected reader-averaged AUC for modality A, we are averaging over readers and cases as we average over the ROC ratings. We express the expectation of d where we have introduced the random variable s Ajkk 0 ¼ sðX Ajk 0 1 − X Ajk0 Þ that we refer to as a "success observation." The equivalences above are true because we can pull the summations out of the expected value and all the summands yield the same result. To be clear, E½s Ajkk 0 is the expected value of s Ajkk 0 for a randomly selected reader reading randomly selected diseased and nondiseased cases. In a sense, the expected value averages over the subscripts, and is not dependent on them. In Eq. (6), X Ajk 0 1 − X Ajk0 is nothing more than the difference of a few fixed effects and many zero-mean random effects. Specifically, it is a normal random variable with a mean and variance given by where σ 2 Ω is the sum of the six variance components that do not depend on modality [Eq. (3)] and σ 2 A is the sum of the six variance components that are specific to modality A [Eq. (4)]. Therefore, we can express the expected value of d AUC A analytically as where Φð·Þ is the cumulative distribution function of the standard normal distribution.

Variances and Covariances
Here, we turn our attention to variances and covariances. Both appear in the variance of the difference of estimated reader-averaged AUCs: Without loss of generality, we first discuss the variance of d AUC A . A similar discussion can be given for the variance d AUC B . We will then discuss the covariance of d AUC A and d AUC B .

Variance
The variance of a single modality can be decomposed into different representations. We shall use the success-moment representation that can be derived using U-statistics; 33 where c A is a vector of coefficients (see Table 2 for the coefficients of the fully-crossed study design considered in this paper) and M A is a vector of eight product moments (see Table 3, every other row of column 1). Each element, ½M A l ðl ¼ 1;2; : : : ; 8Þ, is the expected value of the product of two success observations from that modality. There are eight moments because the two success observations may come from the same or different readers (when reader subscripts match or not), the same or different nondiseased cases (when the nondiseased case subscripts match or not), and the same or different diseased cases (when the diseased case subscripts match or not). To see this clearly or to understand any of the discussion below, it may be useful to write a success observation in terms of the success function acting on a difference in ratings; recall s Ajkk 0 ¼ sðX Ajk 0 1 − X Ajk0 Þ. Furthermore, it may even be necessary to write out the difference in ratings in terms of the constituent random effects.
The particular moment will be clear from the subscripts as, from this point forward, we will not allow different subscripts to take on the same value. The interpretation of each moment is driven by the unique subscripts that appear. For example, in the expression ½M A 1 ¼ Eðs Ajkk 0 s Ajkk 0 Þ we see that the subscripts are identical on both success observations; therefore, the expression is the expected value for a randomly selected reader reading a randomly selected diseased and nondiseased case. In contrast, in the expression ½M A 8 ¼ Eðs Ajkk 0 s Aj 0 k 0 0 k 0 0 0 Þ we see that every subscript on the first success observation is different from its Table 2 Coefficients of the moments that are found in the variance, Eq. (11) and in the covariance (Eq. 13). For the fully-crossed study design considered in this paper, Table 3 Each constituent moment of varðAUC A Þ and covðAUC A ; AUC B Þ calculated using the generalized R&M variance components. Each row is an equation for the variances σ 2 Að·Þ and σ 2 Ωð·Þ that appear in Eq. (12) and (15). For example, in the row corresponding to ½M AB 5 , . Note that the generalized R&M variance components are organized in columns, leaving spaces when variance components are not included.

Moments
Variance components Special case, refer to the text ½M A 8 ¼ E ðs Ajkk 0 s Aj 0 k 0 0 k 0 0 0 Þ∶ Special case, refer to the text ½M AB 8 ¼ E ðs Ajkk 0 s Bj 0 k 0 0 k 0 0 0 Þ∶ Special case, refer to the text counterpart on the second success observation; therefore, the expression is an expected value over a pair of randomly selected readers (that are unique), reading randomly selected diseased and nondiseased cases (that are all unique).
To derive the eight moments, we first treat two special cases, then we discuss the rest. The two special cases are the examples above: ½M A 1 and ½M A 8 . In the first special case, the success observations are identical and both are equal to one or zero (ROC ratings from the generalized R&M model are continuous); consequently, their product equals one or zero. Therefore, In the second special case, the success observations are independent because they come from different readers reading different nondiseased and diseased cases. Therefore, the expected value can be factored and ½M A 8 ¼ Eðs Ajkk 0 ÞEðs Aj 0 k 0 0 k 0 0 0 Þ ¼ AUC 2 A . The remaining moments of M A (moments 2 to 7) can be written as where expressions for the variances σ 2 ΩðlÞ and σ 2 AðlÞ ðl ¼ 2; : : : ; 7Þ are listed in alternate rows of Table 3, and ϕð·Þ is the probability density function of the standard normal distribution. We point out that AðlÞ is the sum of the variances of the random effects that are common to the two success observations within each moment.
Given a fully-specified R&M model, we can compute the moments above (and the moments corresponding to the variance of d AUC B ) using basic numerical integration methods. The derivation of the expression above follows a common outline. We exemplify the derivation for one of the moments in Appendix A. Here is the outline: 1. Identify the random effects that are common to the two success observations. These will be the effects corresponding to shared subscripts.
2. Write the moment as the expected value of the product of two conditional expected values, one for each success observation. Each conditional expected value assumes that the random effects found in Step 1 are fixed.
3. Recognize that each conditional expected value equals a probability that can be expressed in the form ΦðZÞ, where Z is a normal random variable.
4. Integrate over the random variables that were initially fixed in Step 1. This is a one-dimensional (1-D) integral of the product of two normal cumulative distribution functions and one normal probability distribution function.

Covariance
The covariance of d AUC A and d AUC B can be decomposed into a success-moment representation analogous to the variance; namely, where c AB ¼ c A for the fully-crossed study design considered in this paper (see Table 2) and M AB is a vector of eight moments (see Table 3, every other row of column 1). Here, each moment is the expected value of the product of one success observation from modality A and another from B. For example, where the two success observations come from different random readers reading the same cases in different modalities. There is only one special case this time: ½M AB 8 . The success observations for this moment are again independent because they again come from different readers reading different cases. Therefore, the expected value can be factored and ½M AB 8 ¼ The remaining moments of M AB (moments 1 to 7) can be written as where the variances σ 2 ΩðlÞ are listed in every other row of Table 3. We point out that σ 2 Ω − σ 2 ΩðlÞ is the sum of the variances of the random effects that are common to the two success observations within each moment.
Given a fully specified R&M model, we can compute the moments above using basic numerical integration methods. The derivation of the expression above follows the same common outline as above. We illustrate the derivation for ½M AB 5 in Appendix A.

Design of Simulation Studies
In this paper, we mimic the Monte Carlo (MC) simulation experiments run by R&M; however, we modify their experiments by perturbing the input variance components so that the variance components of the ROC ratings depend on modality and truth state. The purpose of the experiments here is to validate the numerical calculations given above against simulated results. Like R&M, we simulate data sets that are fully crossed: all readers read all images in both modalities with no re-reading. In what follows, we describe the original R&M simulations and detail how the simulations in this paper are built on the original. R&M set to zero the effect for nondiseased cases and all four modality-truth interaction effects: Consequently, the expected differences in ROC ratings for both modalities are equal to the effect for diseased cases ðΔ A ¼ Δ B ¼ α 1 Þ. R&M considered three levels of this effect; namely, α 1 ¼ 0.75, 1.50, and 2.50. These three levels ultimately control the expected AUCs via Eqs. (8) and (9). Now, since the variances of the random effects did not depend on modality, they effectively simulated a null-hypothesis experiment: The experiments in this paper will use the same values for the fixed effects. However, because we perturb some variance components in order to generate ROC ratings that depend on modality and truth state (described below), AUC A will not equal AUC B and they will both be different from the AUCs of the original R&M simulation.
R&M investigated experiments with three or five readers and case-sets that ranged from 50 to 400. They also explored ratios of case mixes (nondiseased to diseased) of 1:1 and 9:1. We shall only investigate experiments with five readers, 50 nondiseased cases, and 50 diseased cases. We feel that limiting the investigations related to the size of the experiment is appropriate since the purpose here is more modest than the original purpose (validating a hypothesis test).
R&M based their original simulation experiments on correlation estimates found in actual ROC analyses using the CORROC algorithm, 34 which assumes an underlying bivariate binormal distribution. They investigated high versus low data correlation between and within readers (due to reading the same cases) at both high and low values of total reader variability. We explore a similar set of variance structures, but we perturb them to get variance structures that are different across modality and truth state. Specifically, we halve the variances of effects that involve modality A and truth ¼ diseased, and we double the variances of effects that involve modality B and truth ¼ diseased: This perturbation causes σ 2 A to be different from σ 2 B , which then causes AUC A to be different from AUC B [Eq. (9)].
In total, we explore two levels of between-and within-reader data correlation, four levels of reader variability, and three levels of performance in a factorial fashion for a total of 24 simulation configurations. Details are provided in Appendix B that will allow the interested reader to replicate our experiments. Figure 1 shows the variances and covariances of d AUC A and d AUC B (numerical results versus estimates) of the 24 simulation experiments described above: 8 configurations × 3 performance levels. The x-axis shows the variances by numerical integration. The y-axis shows the MC estimates of variance given 1000 MC trials. There is some variability along the line of equality due to the relatively small number of MC trials. The absolute value of the relative differences between the numerical results and the estimates averaged over all 24 simulation configurations is 5.2%. For 100,000 MC trials, the averaged absolute value of the relative differences is 0.5% and the variability along the line of equality is not visible (not shown).

Conclusions
In this paper, we generalized the R&M model by allowing all of its variance components to depend on both modality and truth state. This will allow investigators to model and simulate MRMC ROC studies that better fit their data. Additionally, we presented and validated equations for computing the population variances and covariances for empirical AUC estimates computed from data simulated from the generalized R&M model. These equations show the core relationships between the ROC data and the reader-averaged AUCs. These equations and relationships should help investigators to validate new MRMC variance and covariance estimation methods, explore novel study designs, size future trials, and model MRMC ROC data.
For the interested investigator, the first author has made available a stand-alone Java application called iRoeMetz. iRoeMetz and its source code can be found in Ref. 35. iRoeMetz can simulate ROC studies according to the generalized R&M model and numerically calculate the expected moments and variances found in this paper. Additionally, all the original R&M simulation configurations and those used in this paper can be downloaded and tested.
Future work is needed in this area, investigating real data sets and calibrating the different possible R&M models of decision scores to be consistent with real data estimates of AUC variance components. To cover the full range, datasets that study different anatomical locations and diseases, and different imaging modalities (with and without computer aids as appropriate) are needed. The results of such future work would provide investigators with more information and examples to help them calibrate the R&M model (or any other model) to their situation so that they can appropriately size future reader studies. We are collecting datasets for this work for public sharing toward this goal; we welcome any and all contributions.
There are other interesting future efforts that could support the need mentioned above. One direction of study could link the success moments at the core of the work here to the parameters of the Obuchowski and Rockette 30 model for AUC observations. Another direction of study could be to uncover the inverse mapping: given the parameters of the Obuchowski and Rockette model for AUC observations, solve for the variance components of the ROC ratings that are consistent with them. It is not currently known whether the solution is unique or if it can be treated efficiently and effectively.

Appendix A
In this appendix, we derive one of the product moments, ½M AB 5 , that enters into the computation of the covariance between other moments that are required for computing the variance of d AUC A and the variance of d AUC B is straightforward. Although the generalized R&M model has many parameters and the variances and the covariance of the reader-averaged AUCs are fairly complex, all of the success moments are derived the same way and the results have a common form. The derivations follow a common outline. The results are straightforward 1-D integrals involving the normal distribution that can be calculated using basic numerical integration methods. We illustrate the derivation for ½M AB 5 .
Here is the derivation outline: 1. Identify the random effects that are common to the two success observations. These will be the effects corresponding to shared subscripts.
2. Write the moment as the expected value of the product of two conditional expected values, one for each success observation. Each conditional expected value assumes that the random effects found in Step 1 are fixed.
3. Recognize that each conditional expected value equals a probability that can be expressed in the form ΦðZÞ, where Z is a normal random variable equal to a difference in ROC ratings with the random effects found in Step 1 temporarily fixed.
4. Integrate over the random variables that were initially fixed in Step 1. This is a 1-D integral of the product of two normal cumulative distribution functions and one normal probability distribution function.
Step 1: The fifth product moment is The random effects that are common to both success functions are the ones that do not depend on reader or modality: ½Cα k 0 1 ; ½Cα k0 . Note that these random effects are the effects corresponding to shared subscripts. This can be seen if we write out the differences in ROC ratings in terms of the constituent random effects.
Step 2: Here, we first write the moment as the expected value of a conditional expected value, where the random effects found in Step 1 are fixed; namely, ½̱ M AB 5 ¼ E½Eðs Ajkk 0 s Bj 0 kk 0 j½Cα k 0 1 ; ½Cα k0 Þ With the common random effects fixed, s Ajkk 0 is independent of s Bj 0 kk 0 . Therefore, we can factor the conditional expected value to obtain ½̱ M AB 5 ¼ E½Eðs Ajkk 0 j½Cα k 0 1 ; ½Cα k0 Þ × Eðs Bj 0 kk 0 j × ½Cα k 0 1 ; ½Cα k0 Þ: Step 3: In the expected value concerning modality A, the argument of the success function is a normal random variable Z with mean and conditional variance given by Using the definitions given in Eq. (4) and Table 3, we see that σ 2 Z ¼ σ 2 A þ σ 2 Ωð5Þ . We have a similar result for the argument of the success function in the expected value concerning modality B. Since the argument to each conditional expected value is a normal random variable, we can write Step 4: Last, we see that the only randomness that remains above is the difference ½Cα k 0 1 − ½Cα k0 . This difference is a normal random variable Y ¼ ½Cα k 0 1 − ½Cα k0 with mean zero and variance Ωð5Þ q we can effectively integrate over the random variables that were initially fixed, ½Cα k 0 1 ; ½Cα k0 , with the following: Ωð5Þ q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi σ 2 B þ σ 2 Ωð5Þ q 1 C AϕðxÞdx: This last expression can be calculated using basic numerical integration methods. In our software, we simply sample the 1-D integral at 256 points on the interval (−10; 10) and use the midpoint rule (rectangle rule).