SEQUENTIAL TESTING OF MEASUREMENT ERRORS IN INTER-RATER RELIABILITY STUDIES

Inter-rater reliability studies are conducted to investigate the reproducibility and level of agreement on assessments among different raters. The analysis of reliability is a common feature in research and practice since most measurements involve measurement errors, particularly those made by humans. Measurement error can seriously affect statistical analysis and interpretation; it therefore becomes important to assess the amount of such error by calculating a reliability index ( Shrout and Fleiss (1979 )). Inter-rater reliability refers to the reproducibility of the raters when repeated at random on the same subject or specimen under the same condition by the same rater or different ones, either simultaneously or at several time points. A common measure of the reliability of measurements is the so-called intraclass correlation coefficient, with larger values indicating higher level of consistency. There are various versions for intraclass correlation coefficient derived from different statistical models. Shrout and Fleiss (1979) proposed a set of guidelines for choosing the appropriate model in reliability studies. In actual inter-rater reliability studies, multiple raters usually evaluate multiple subjects and the two-way analysis of variance (ANOVA) model is considered to be the appropriate analytical model. Moreover, when raters can be considered as random sample from the target population of raters, the two-way random effects model can be applied ( Fleiss (1999) ; Zou and McDermott (1999) ; Cappelleri and Ting (2003) ; Tian and Cappelleri (2003) ). For fixed sample design, there are two main approaches available to determine the sample size. One approach is to find the size such that the power of statistical test is met ( Donner and Eliasziw (1987) ; Walter, Eliasziw, and Donner (1998) ) while the other is to assure the precision of estimation ( Shoukri, Asyali, and Walter (2003) ; Bonett (2002) ; Saito et al. (2006) ). Sequential methods were introduced in response to demands for more efficient testing of anti-aircraft gunnery during World War II, culminating in Wald’s development of the sequential probability ratio test which had an immediate impact on weapons testing ( Wald (1947) ; Siegmund (1985) ; Lai (2004) ). Gradually, in studies involving human subjects, specialized statistical methods were called for to balance the ethical and financial advantages of stopping a study early against the risk of an incorrect conclusion ( Jennison and Turnbull (2000) ). Armitage (1954 , 1958 ) and Bross (1952 , 1958 ) introduced the use of sequential methods in the medical field. However, the early development for group sequential methods came from Pocock (1977) , O’Brien and Fleming (1979) , and Lan and DeMets (1983) . Motivated by the idea of sequential testing that is widely used in clinical trials, it is natural to adopt and extend these sequential testing methods in the design and analysis of reliability studies to reduce the sample size and study cost. In reliability studies evaluating the measurement error by applying the one-way ANOVA model, the multistage group sequential designs were proposed. Under one-way ANOVA, the sums of squares in the estimation of the intraclass correlation coefficient possess independent increments, thus simplifying the calculation of stopping boundaries ( Liu, Schisterman, and Wu (2006) ). In this paper, we develop multistage testing procedures using two-way ANOVA for hypotheses concerning the intraclass correlation coefficient in a inter-rater reliability study. In Section 2, we state the hypotheses of interest, introduce the structure and assumption of the two-way ANOVA models, and propose simulation designs for the one-stage problem for sample size and power calculation. In Section 3, we develop methods to determine critical values, sample size, and power using Lan and DeMets’s (1983) error spending approach. Realizing that the between-rater sum of squares violates the independent increments assumption, we develop simulation techniques to effectively calculate the critical values. The performance of the proposed methods is examined in Section 4 using simulation studies, and critical values are tabulated for a range of two-stage design parameters. In Section 5, we exemplify the methods using data from the Physician Reliability Study for diagnosis of endometriosis. Finally, summary and discussion are given in Section 6.

The intraclass correlation coefficient ρ is commonly used to assess inter-rater reliability. The designs developed in this paper are based on testing a null hypothesis H 0 : ρ ≤ ρ 0 that the true intraclass correlation coefficient is less than some uninteresting level ρ 0 . We require that under the null hypothesis the probability of concluding that the rating is sufficiently promising be less than α , the level of significance of the test. We also require that under a specified alternative hypothesis H 1 : ρ ≥ ρ 1 , the probability of rejecting the rating should be less than β . If the null hypothesis is not rejected, then we may conclude that the raters are not well trained or that the variable to be measured possesses unacceptable measurement errors. Consider a random sample of n subjects for which a continuous variable y is measured independently by d raters randomly selected from a population of raters. Three sources of variation usually occur from such a design: subjects, raters, and random errors. Higher inter-rater reliability is achieved if the variations from random errors and raters are relatively smaller than the variation from subjects. Denote by y ij the measurement made on the i th subject by the j th rater, i = 1, …, n and j = 1, …, d . Then the two-way ANOVA model is (2.1) y i j = μ + s i + r j + ε i j , where μ is the grand mean of all measurements. The subject effect s i , the rater effect r j , and the random error ε ij are assumed to be independent and normally distributed with mean 0 and variances σ s r , σ r 2 and σ ε 2 , measuring the magnitude of the variation from the three resources, respectively. Note that these assumptions lead to Var ( y i j ) = σ s 2 + σ r 2 + σ ε 2 , Cov ( y i j , y i ′ j ) = σ r 2 if i ≠ i ′, Cov ( y i j , y i j ′ ) = σ s 2 if j ≠ j ′, and Cov( y ij , y i ′ j ′ ) = 0 if i ≠ i ′ and j ≠ j ′. The vector including measurements from all subjects y 1 × n d T = ( y 11 , ⋯ , y 1 d , ⋯ , y n 1 , ⋯ , y n d ) is then distributed as N nd ( u , Σ), (2.2) u = μ 1 n d , ∑ = σ ε 2 I n ⊗ I d + σ s 2 I n ⊗ J d + σ r 2 J n ⊗ I d , where, throughout, the superscription T stands for the transpose of a matrix or vector, 1 m is the vector of order m with all elements 1, I n is the identity matrix or unit matrix of size n , and ⊗ is the Kronecker product of two matrices ( Rao (1973) ), and J n × m is the n × m matrix with all elements 1, with a simpler notation J n (= J n × n ) being used instead if m = n . The intraclass correlation coefficient in a two-way ANOVA model is given by It provides a measure for the overall consistency among measurements from the same subject. The index reduces to the intraclass correlation coefficient in a one-way ANOVA model, if there is no variation among raters. Furthermore, it can be shown that ρ is the correlation coefficient between two measurements y ij 1 and y ij 2 on the same subject i by two different raters j 1 and j 2 . Thus, larger values of ρ indicate higher coherence among measurements on the same subject by different raters. Perfect consistency among these measurements occurs with ρ = 1 when subjects are the only source of variation, in which case σ r = σ ε = 0. A low reliability may indicate that the raters are not well trained or that the variable to be measured possesses unacceptable measurement errors. Consider the sums of squares S S subject = ∑ i = 1 n ∑ j = 1 d ( Y ¯ i . - Y ¯ .. ) 2 , S S rater = ∑ i = 1 n ∑ j = 1 d ( Y ¯ . j - Y ¯ .. ) 2 , S S error = ∑ i = 1 n ∑ j = 1 d ( Y i j - Y ¯ i . - Y ¯ . j + Y ¯ .. ) 2 , where Here ( Rao (1973) ) SS subject , SS rater , and SS error are independent and can be written as (2.4) S S subject ≅ ( d σ s 2 + σ ε 2 ) V 1 , S S rater ≅ ( n σ r 2 + σ ε 2 ) V 2 , S S error ≅ σ ε 2 V 3 , where V 1 , V 2 , and V 3 are independent chi-square variables with degrees of freedom n − 1, d − 1, and ( n − 1)( d − 1), respectively, and ≅ means that its two sides have the same distribution. Take MS subject as SS subject divided by the degrees of freedom of its corresponding χ 2 -distribution, and similarly define MS rater and MS error . Then E ( M S subject ) = d σ s 2 + σ ε 2 , E ( M S rater ) = n σ r 2 + σ ε 2 , E ( M S error ) = σ ε 2 , yielding unbiased estimators of σ s 2 , σ r 2 and σ ε 2 as An estimator of the intraclass correlation coefficient ρ can be obtained by substituting the variance estimators, so (2.6) ρ ^ = σ s 2 ^ σ s 2 ^ + σ r 2 ^ + σ ε 2 ^ = n ( d - 1 ) S - n n ( d - 1 ) S + d ( n - 1 ) R + ( n d - n - d ) , where (2.7) R = S S rater S S error ≅ ( n ω + 1 ) V 2 V 3 ~ n ω + 1 n - 1 F d - 1 , ( n - 1 ) ( d - 1 ) , (2.8) S = S S subject S S error ≅ d ρ ( ω + 1 ) + ( 1 - ρ ) 1 - ρ V 1 V 3 ~ d ρ ( ω + 1 ) + ( 1 - ρ ) ( d - 1 ) ( 1 - ρ ) F n - 1 , ( n - 1 ) ( d - 1 ) , where ω = σ r 2 / σ ε 2 and F a , b is the central F-distribution with degrees of freedom a and b . We refer to ρ̂ as the sample intraclass correlation coefficient. The two-way ANOVA model reduces to the popular one-way ANOVA model if σ r 2 = 0 , that is, if there is no variations among raters. In this case, observations from the same subject are treated as interchangeable repeated measurements, and observations from different subjects are independent. For the one-way ANOVA model, the sample intraclass correlation coefficient is a function of an F-statistic, and methods are directly available for computing the powers and sample sizes of the test based on the F-statistic; see, among others, Donner and Eliasziw (1987) , Walter, Eliasziw, and Donner (1998) , and Liu, Schisterman, and Wu (2006) . For the two-way ANOVA model observations from different subjects but the same rater are not independent, and methods for power and sample size calculations based on one-way ANOVA are not applicable because the sample intraclass correlation coefficient involves correlated F-statistics, S and R in ( 2.6 ), arising from the rater effects. Kraemer (1976) considered reliability in a two-way ANOVA model; however, the intraclass correlation coefficient there differed from the traditional one in that the rater variability was excluded. For a pre-specified significance level α , the fixed sample size design rejects the null hypothesis H 0 : ρ ≤ ρ 0 if the sample intraclass correlation coefficient exceeds some critical value 0 < c c ) depends, in addition to ρ , on n , d , and the variance ratio ω . It can be shown (see Appendix ) that, for given n , d , ω , Ψ is an increasing function in ρ , implying that we can control the type I error rate at the nominal level α by setting Ψ = α for ρ = ρ 0 , and set the power to be 1 − β at ρ 1 for ρ ≥ ρ 1 . To deal with the nuisance parameter ω , we can take the conventional approach of finding 0 < c < 1 such that sup ω Ψ( ρ 0 ; ω ) = α . Similarly, given the power 1 − β at an alternative value ρ = ρ 1 , the sample size n can be determined from inf ω Ψ( ρ 1 ; ω ) = 1 − β . Allowing the nuisance parameter ω to change freely in the entire range [0, ∞] often results in a very conservative test for values of ω encountered in practice. Very often, data from other similar studies can provide an empirical estimate of ω , thus making it possible to have an upper bound ω 1 for ω . Taking these arguments into account, we propose to determine the critical value and sample size via the equations Exact calculation of the operating characteristics of the test statistic ρ̂ appears to be complicated because it involves the joint distribution of two correlated F-statistics S and R in ( 2.6 ). As ρ̂ is a function of three independent χ 2 random variables, we propose to calculate the operating characteristics based on independently generating large numbers of replicates of the χ 2 variables. Specifically, for a given set ( n , d , ω , ρ , c ), the power Ψ( ρ; ω ) is computed as the proportion of simulated values of ρ̂ that are larger than c . To generate a random value ρ̂ from a simulation, we generate a random observation of ( V 1 , V 2 , V 3 ), which then yields random values of R and S according to ( 2.7 ) and ( 2.8 ). Substituting these values in ( 2.6 ) then yields the simulated value of ρ̂ . Using this approach, we can find the critical value and sample size that satisfy the error requirements ( 2.10 ). Note that the sample size is taken as the smallest n that meets the error requirements. Setting n = 104, d = 4, α = 0.05, ρ 0 = 0.5 and ω 1 = 0.5, we found c = 0.6173, using the approach. For ω ≤ ω 1 , Figure 1 plots the type I error Ψ( ρ 0 ; ω ) and Figure 2 plots the power Ψ( ρ 1 ; ω ) at ρ 1 = 0.7 and 0.8, both versus ω . Although no analytic proof is available, it appears from the figures that, as ω increases, the type I error increases and the power decreases (for alternative values of ρ sufficiently away from the null value ρ 0 ). Thus for ω ≤ ω 1 , ( 2.10 ) implies

Sequential statistical methods were originally developed for economic benefits. In any experiment or survey in which data accumulate steadily over a period of time, it is natural to monitor results as they occur with a view to taking such action as early termination or a modification of the study design. The many reasons for conducting interim analyses of accumulating data can be loosely categorized in three classes: ethical, administrative, and economic ( Jennison and Turnbull (2000) ). For example ( Rothstein (1990) ), in job performance rating, such critically important decisions as pay, promotions, and firing are made on the basis that the ratings are reliable enough to produce useful and valid performance appraisals. An interim analysis is needed to ensure that the rating process is being executed as planned, and that the raters are well trained and satisfy eligibility criteria. An early examination of interim results can reveal the presence of problems that can be remedied before too much expense is incurred. Early stopping ensures that resources are not wasted. Sequential methods typically lead to savings in sample size, time, and cost when compared with standard fixed sample procedures. Suppose the subjects are taken into groups. Denote by n k the number of subjects in the k th group, and let N k = n 1 + n 2 + ··· + n k . Denote by ρ̂ k the estimate of the intraclass correlation coefficient ρ computed from data observed up to the k th group. Let K be the prespecified maximum number of interim analyses planned for the study. A K -stage testing procedure for H 0 is carried out as follows. At each stage k = 1, …, K − 1, with critical value c k , stop if (3.1) ρ ^ i ≤ c i , i ≤ k - 1 , ρ ^ k > c k , and reject H 0 . Otherwise, d measurements are taken from the raters for each of the n k +1 subjects in the ( k + 1)th group, and ρ̂ k +1 is computed based on all available data up to the ( k + 1)th stage. If the test has not stopped, stop at the K th stage and reject H 0 if ρ̂ K > c K . Here c k , k = 1, …, K , are the critical values that satisfy the error requirements. The power function of the group sequential test is (3.2) ∑ k = 1 K P ρ ( ρ ^ i ≤ c i , i ≤ k - 1 , ρ ^ k > c k ) , which depends on the variance ratio ω . The design approach considered here is to specify the parameters ρ 0 , ρ 1 , α , and β , along with the nuisance parameter ω , and then determine the critical values and sample sizes that have power α at ρ 0 and 1 − β at ρ 1 . To determine sample sizes and critical values for testing the null hypothesis H 0 , we employ the Lan and DeMets error spending approach ( Lan and DeMets (1983) ). We focus on two-stage designs that are more appealing in terms of complexity of computation and potential applications to a reliability study. However, the simulation approach to a two-stage test can be straightforwardly extended to more stages. In a two-stage study, the null hypothesis H 0 is rejected if (i) at the first stage ρ̂ 1 > c 1 or (ii) ρ̂ 1 ≤ c 1 at the first stage and ρ̂ 2 > c 2 at the second. Suppose that the type I error to spend at the first stage is α 1 . We have (The errors spent at each stage then add up to the total error of α .) The requirement that the power at ρ 1 be 1 − β yields These three equations determine the critical values and the sample size assuming the allocation ratio, n 1 / N 2 , of subjects to the first group is specified. The critical value c 1 for the first stage and the first summand of ( 3.5 ) can be determined by the approach proposed in Section 2 for a fixed sample size test. Computation of ( 3.4 ) and the second summand in ( 3.5 ) involves the joint distribution of ρ̂ 1 and ρ̂ 2 or, equivalently, the joint distribution of two correlated χ 2 -statistics (apart from proper constants), For the one-way ANOVA model, Liu, Schisterman, and Wu (2006) showed that the independence between subjects warrants that the sums of squares involved in estimation of the intraclass correlation coefficient can be decomposed into independent increments. This independent-increment property substantially simplifies the simulation for calculation of critical values and power of the test. For the two-way ANOVA model, the between-subject sum of squares SS subject and the within sum of squares SS error have independent increments: SS subject (1), SS subject (2) − SS subject (1), SS error (1), and SS error (2) − SS error (1) are mutually independent. The four increments have the same distributions as ( d σ s 2 + σ ε 2 ) U 1 , ( d σ s 2 + σ ε 2 ) U 2 , σ ε 2 U 3 , and σ ε 2 U 4 , respectively, where U 1 , U 2 , U 3 and U 4 are independent chi-square random variables with degrees of freedom n 1 − 1, N 2 − n 1 , ( n 1 − 1)( d − 1), and ( N 2 − n 1 )( d − 1), respectively. These sums of squares at the second stage can be expressed as functions of independent chi-square variables, similar to the one-way ANOVA model in Liu, Schisterman, and Wu (2006) . A major complication arises from calculation of the between-rater sum of squares SS rater (2). Because observations from different subjects made by the same rater are not independent, the between-rater sums of squares do not have the independent increments property; see the Appendix . Therefore, unlike the other two sums of squares, SS rater (1) and SS rater (2) cannot be generated using independent chi-square variables. We propose a more efficient method for simulation based on two independent multivariate distributions. With specified parameters the probabilities in ( 3.3 ) – ( 3.5 ) can be estimated by simulating a large number of realizations of the sums of squares. The probability P ( ρ̂ 1 > c 1 ) can be derived using simulation techniques described in Section 2.4. To compute the probability P ( ρ̂ 1 ≤ c 1 , ρ̂ 2 > c 2 ), we first simulate a large number of observations of U = ( U 1 , U 2 , U 3 , U 4 ), defined as before. For each simulation, the between-subject and within sum of squares are Empirical values of SS rater (1) and SS rate (2) can be generated from two independent multivariate normal variables. Details are given in the Appendix . With each simulated value of SS subject (1), SS sugject (2), SS error (1), SS error (2), SS rater (1), and SS rater (2), we compute ρ̂ 1 and ρ̂ 2 . The probability is then estimated as the proportion of {( ρ̂ 1 , ρ̂ 2 ) : ρ̂ 1 ≤ c 1 , ρ̂ 2 > c 2 }. Although the probabilities can be estimated by simulating N 2 -variate multivariate normal distributions under model ( 2.1 ), our experiences show that the proposed simulation approach based on independent chi-square and multivariate normal variables is much more time-efficient.

To set up the simulation, set n 1 = N 2 /2. The error spending function of Kim and DeMets (1987) with α 1 = α /2 is applied to allocate the error. The values of design parameters considered were α = (0.025, 0.05), ρ 0 = (0.5, 0.6), and 1 − β = (0.8, 0.9) at ρ 1 = (0.7, 0.8, 0.9). Here, according to Landis and Koch (1977) , the chosen null values indicate moderate consistency in measurements. We considered the numbers of raters (4, 6, 8), and the values of ω to be (0.1, 0.5, 1, 10). For each set of design parameters ( α , β , ρ 0 , ρ 1 , ω ), we searched a potential range of ( c 1 , c 2 , N 2 ) to find a combination that meets the error requirements. For each ( c 1 , c 2 , N 2 ), c 1 could be derived using the simulation techniques described in Section 2.4. A total of m = 10, 000 random vectors U = ( U 1 , U 2 , U 3 , U 4 ) were generated. Simultaneously, two samples each with m = 10, 000 observations were drawn from independent multivariate normal variables as defined in the Appendix . The empirical probabilities in ( 3.3 ), ( 3.4 ), and ( 3.5 ) were computed as the proportion that fell into the rejection regions. The combination that satisfied the error requirements was thus the desired solution to the design. To help narrow the search range, we first limited the values of c 2 and N 2 close to those required in a fixed sample design computed by methods mentioned in Section 2.4, and then gradually expanded the range until a solution is found. Table 1 presents the simulation results of critical values and sample size required, for a selected range of designs. For each design, the final critical value c 2 is smaller than the interim critical value c 1 , partly reflecting the increase of sample size. The Average Sample Number (ASN) under the alternative is smaller than the fixed sample size. Note that for a two-stage design the ASN at ρ is ASN = n 1 P ( ρ̂ 1 > c ) + N 2 P ( ρ̂ 1 ≤ c 1 ). All simulations were conducted using the R software. Codes are available from the authors upon request.

We illustrate the proposed group sequential testing procedure on data from a study on endometriosis diagnosis. Endometriosis is a gynecological medical condition occurring in roughly 5%–10% of women. Despite its relatively high prevalence, substantive and methodological challenges exist, including diagnostic proficiency. The Physician Reliability Study ( Schliep et al. (2012) ) addressed this issue by examining agreement in diagnosing endometriosis and its association with both rater- and subject-specific information. We used the review results of four experts who independently examined study subjects’ intraoperative photos for signs of endometriosis or other gynecologic pathology. The specific outcome we used was the Revised American Society for Reproductive Medicine’s classification that ranges from 1 to 150 with higher scores corresponding to severe symptoms of endometriosis. The interim analysis was conducted based on the first 24 (about half of the 47 samples that have complete ratings) subjects’ scores from the four experts. The interim error α 1 was set to α /2 based on an error spending function of Kim and DeMets (1987) . Thus, n 1 = 24, N 2 = 47, and α 1 = α /2 = 0.025. The null value is ρ 0 = 0.4, and the variance ratio ω is 0.013 based on the real data. With these design parameters, we applied the simulation approach of Section 3.3, and found c 1 = 0.5953 and c 2 = 0.5359. The power of the test at ρ 1 = 0.55, 0.60, 0.65 was 0.6088, 0.8150, and 0.9548, respectively. The average sample size (ASN) at these alternative values was 40, 36, and 30, respectively, reflecting the cost-effective nature of the design. We then applied the test based on the scores on the first 24 subjects in the endometriosis study. The sample intraclass correlation coefficient computed at the first stage was ρ̂ 1 = 0.6077. We then could terminate the test at the first stage and reject the null hypothesis because the test statistic ρ̂ 1 exceeded the critical value c 1 . Early evidence is that the four experts demonstrated quite satisfactory agreement in their scoring.

In this paper, we develop multistage testing procedures using two-way ANOVA for hypothesis concerning the intraclass correlation coefficient in a inter-rater reliability study. The designs are proposed to determine critical values, sample size, and power using Lan and DeMets’s error spending function. To our best knowledge there have been no exact tests available for the intraclass correlation coefficient based on two-way ANOVA framework due to the involvement of two correlated F statistics. Realizing that the between-rater sum of squares violates the independent increments assumption, we develop simulation techniques to effectively calculate the operating characteristics of the sequential test. The performance of the proposed methods is examined using simulation studies and critical values are tabulated for a range of two-stage design parameters. The operating characteristics of the test are affected substantially by the variance ratio ω . Our proposed approach requires an upper bound of ω be specified. The methods discussed in this paper could be applied to the field of Genetics in the sense of the heritability analysis that estimates the relative contributions of differences in genetic and non-genetic factors to the total phenotypic variance in a population. Future research is needed to develop and compare various group sequential designs, and to propose methods for inference following the sequential testing such as the point and interval estimation of the intraclass correlation coefficient.