Original Research

Assessing inter-rater agreement between multiple medical instruments with heteroscedastic measurements

© 2017 Nawarathna et al; licensee Herbert Publications Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

---

Introduction

In the field of clinical medicine, measurements of the human body take a major part of diagnostic, prognostic and therapeutic evaluations. Due to the rapid advancement technology, new methods and instruments are introduced into this field. These introduced instruments or methods might be more advance, cheaper and easier to use than the old standard instruments. Before these newly introduced methods or instruments put into use, the accuracy and precision of the measurements need to be verified. If these instruments agree sufficiently well with the already existing methods, they can be used interchangeably. To get an understanding of these measurements, they must be compared against with a well-established technique. Then we need to assess the degree of agreement between these methods [1-5].

Experts in this field have already introduced many techniques to test the degree of agreement between two methods. But most of these techniques can be applied under certain assumptions.

In order to assess the agreement between the assays, the first step is to model the data. The Linear mixed effects model is most commonly used in modeling the method comparison data [6-15]. Because of the flexibility in modeling of within subject dependence, linear mixed models are popular. In this model, normality is assumed for the error terms, but in practice sometimes the normality assumption is violated. In such cases, linear mixed models cannot be used to model the data [16].

To overcome this problem [17] proposed a robust mixed model called "General skew-t mixed model (GSTMM)" that assumes a multivariate skew-t distribution for random effects and an independent multivariate t-distribution for the errors. But this general skew-t mixed model (GSTMM) and most of the other existing models available for method comparison studies are based on the assumption that the variability of the continuous measurement remains a constant throughout the range of the measurement. Though, this is not the case in some practical situations and the variability of the measurements might change with the magnitude, i.e., the 'Heteroscedasticity' of error terms [18].

A novel model to method comparison data with heteroscedastic error variances is proposed in [15], to evaluate the agreement between two methods of measurements measuring continuous data. However, this model cannot accommodate the comparison of multiple methods of heteroscedastic measurements. Therefore, the main objective of this study is to propose a heteroscedastic mixed effects model to analyze heteroscedastic method comparison data with more than two methods of measurements and to adapt the agreement evaluation methodology for multiple methods, taking heteroscedasticity into account.

The lung tumor size measurements data from [19] which motivated this work and are analyzed later in this article, provide a specific example of this phenomenon. The information was gathered from August 2000 to May 2001. In the dataset, out of 33 patients there are 40 lung tumors. These 40 lung tumors belong to 16 men and 17 women from the ages 43 to 78 years. All the lung tumors are larger than 1.5 cm in maximal diameter. Computed Tomography (CT) images of these 40 lung tumors are distributed among five radiologists. All five radiologists are with Thoracic Fellowship training and have more than four years of post-training experience. The five radiologists measured the lung tumors on the CT images by using a ruler or with calipers. Each of these measurements is performed independently and also each of these images was measured twice by each inspector. Here we will consider the measurements taken by the five radiologists as different methods of measurements. Therefore, the dataset contains 40 lung tumors (subjects), 5 readers (methods) and 2 replications for each measurement. In total 40*5*2=400 records.

The rest of this article is organized as follows. Section 2 presents the proposed methodology to deal with heteroscedastic method comparison data measured by multiple methods or raters. Section 3 discusses agreement evaluation under the proposed model using two techniques of model fitting and the last section discusses the results and conclusions of the study. The model fitting and the analysis was carried out by using the R statistical software.

Methodology

In this study, the methodology was divided into two main parts. The first step is to model the data set using an appropriate model and the second step is to measure the agreement among the methods of measurements. The most popular statistical modeling technique in method comparison data is the "Liner mixed effects model". This model is an extension of linear regression models [12] and contains both fixed effects and random effects. In practice the mixed effect models are important when there are repeated measurements. The standard mixed model in the matrix form can be represented as follows.

$$Y_i=X_i\beta +Z_i{b}_i+e_i;i=1,\dots ,n,$$ (1) (1)

where, Y_i is the vector of observed responses on the i^th subject, X_i, Z_i are design matrices for fixed and random effects, ß is the vector of fixed effects, b_i is the vector of random effects, e_i is the error vector and e is the vector of all unknown parameters in the model.

The main assumptions of this model are the random effects, random errors have normal distributions and they are mutually independent and the model has a constant error variance which depends only on the method.

$$\begin{array}{l}{b}_i~independentN\left(0,\text{Ψ}\right);e_i~independentN\left(0,R_i\right);\\ R_i=diag\left({\sigma }_1^2,\dots ,{\sigma }_1^2,{\sigma }_2^2,\dots ,{\sigma }_2^2,\dots ,{\sigma }_n^2,\dots ,{\sigma }_n^2\right)\end{array}$$

where Ψ is a lxl positive definite matrix with diagonals Ψ₁²,........,Ψ₁², σ_j² is the error variances of method j. After fitting the dataset with the mixed effects model, the residuals were analyzed to identify whether there are any model assumption violations.

Heteroscedastic Mixed Effects Model for multiple methods
This model is used when the assumptions (i.e., the random effects and the random errors are normally distributed, independent and has a constant variance) of the standard mixed effect model (homoscedastic model) are violated where error variance changes with the magnitude of the measurements. The heteroscedastic model is as follows.

Let Y_ijk is the k^th replicate measurement by the j^th method on the ith subject, where k=1,...,n_ij, j=1,...., and i=1,.....,m. Here n_ij is the number of measurement from the j^th method on the ith subject, ß_j is the fixed mean of the j^th method, b_ijis the random effect on the i^th subject on the j^th method, e_ijk is the within subject random error, ν is the variance covariate and the ν_i is values for the ith subject, z(ρ_ij) is the Fisher's z-transformation of the correlation coefficient between methods i and lastly δ is the heteroscedastic parameter, when δ=0 this assumes homoscedasticity, µ_ij=β_j+b_ij denotes the conditional mean of E (Y_ijk|b_i) and ν_i=h (μ_i) denotes that the νi is a function of μi and Σ_ij(ν_i) is n_ij×n_ij diagonal matrix and Σ_ij(ν_i) is n_i×n_i diagonal matrix where, b_i =(b_i1, b_i2,....b_il)

$${\text{Σ}}_{ij}\left(v_i\right)=diag\left\{{\sigma }_j^2{g}_j^2\left(v_i,{\delta }_j\right),\dots .,{\sigma }_j^2{g}_j^2\left(v_i,{\delta }_j\right)\right\}$$

and

$${\text{Σ}}_i\left(v_i\right)=diag\left\{{\text{Σ}}_{i1}\left(v_i\right),{\text{Σ}}_{i2}\left(v_i\right),\dots .,{\text{Σ}}_{ij}\left(v_i\right)\right\}$$

Variance covariate is a function of magnitude of measurement μ_i that will be used to model the error variances. A model for conditional error variance is as follows,

$$\mathrm{var}\left(e_{ijk}\text{|}{b}_i\right)={\sigma }_j^2{g}_j^2\left(v_i,{\delta }_j\right),{\mu }_i=\beta +b_i,j=1,2,\dots ,l,$$

where g_jis the variance function and δ_j is the heteroscedasticity parameter of the method j. According to [12] there exists some common variance functions such as, power model: (ν,δ)=| ν |^δ, constant plus power model:

$$\begin{array}{l}\left(v,\delta \right)=|v{|}^{\delta }\text{constant plus power model}:\\ g\left(v,\delta \right)={\delta }_0+|v{|}^{{\delta }_1}\text{and exponential model}:\end{array}$$

g(ν,δ)= exp (δ_ν). The heteroscedastic model in the matrix form is as follows.

$$\begin{array}{l}{Y}_i=X_i\beta +Z_i{b}_i+e_i;i=1,\dots ,m\\ {\text{e}}_{\text{i}}|{\text{b}}_{\text{i}}\text{~ independent}{N}_{n_i}\left(0,{\text{Σ}}_i\left(v_i\right)\right),b_i~independent{N}_l\left(0,\text{Ψ}\right).\\ \end{array}$$

Due to the scarcity of closed form for the likelihood functions, the exact modeling approaches will be troublesome. Therefore, the model fitting is carried out by two model approximations. First using the mean vector of the reference method as the covariate and the second is using the Best Linear Unbiased Predictor (BLUP) as the covariate. In this case the likelihood functions will be possible in a closed form. Score test and likelihood ratio tests have been carried out to assess the validity of these models. After confirming the models validity, the agreement evaluation is done by using Concordance correlation coefficient (CCC) and Total deviation index (TDI).

Model fitting
The model fitting for this heteroscedastic mixed effects models is carried out by selecting an appropriate value for the variance covariate. In this study, there are two main options for the variance covariate, namely mean of the reference method and the BLUP as the variance covariate. The power model function g(ν,δ)=| ν |^δ was selected to fit the model [12].

Using observable mean measurement as variance covariate
As the first approach, we select the mean vector of the reference method as the variance covariate. In this case μ_i^*, an observable quantity is fixed and can fit the model by maximum likelihood method. Here the μ_i^* is expected to be close to μ_i and μ_i=Y_i.

Using BLUP as the covariate
As the second approach, we use the best linear unbiased predictors (BLUP) as the variance covariate and fit the data using the heteroscedastic model. Random effects and error terms are independent in this heteroscedastic model and according to the [12] the BLUP can be written as,

$${\mu }_{i,blup}=E\left(b_i\text{|}{Y}_i\right)=\beta +b_{i,blup};b_{i,blup}=E\left(b_i\text{|}{y}_i\right);v_i^{*}=h\left({\mu }_i{}^{*}\right).$$

In order to calculate the μi,blup any statistical software can be used. ν_i^* depends on the unknown parameter θ. The method of calculating the μ_i,blup is an iteratively reweighted method. This method has two steps and it is continued until it converges [12]. More details of this can be found on [15].

When considering the two approaches, the method of using the true mean as the covariate is much simpler than the method of using BLUP as the covariate, because of the complexity of calculating the BLUP than calculating the mean. However, method of using BLUP is more accurate than the other [14].

Agreement evaluation under the proposed
Methodology
In health science researches, agreement evaluation is a topic which has considerable interest. This is assessing the agreement between two or more methods measuring the same response [20,21]. In this Section we discuss two measurements that are used to assess the agreement in this study [22,23].

Concordance correlation coefficient (CCC)
Concordance correlation coefficient (CCC) is one of the most common measurements used in order to assess the agreement between methods of measurements. This was introduced by [22]. CCC value ranges between -1 to 1. The higher the values it gives a better agreement. Concordance correlation coefficient under the proposed heteroscedastic mixed model is follows for the lung tumor size measurements,

$$CCC\left(v_0\right)=\frac{2{\psi }_{1j}}{{({\beta }_1-{\beta }_j)}^2+{\psi }_1^2+{\sigma }_1^2{g}_1^2\left(v_0,{\delta }_1\right)+{\psi }_j^2+{\sigma }_j^2{g}_j^2\left(v_0,{\delta }_j\right)};j=2,3,4,5.$$

This represents the concordance correlation coefficient between the reader1 and the reader j. For greater accuracy of the measurement CCC was first calculated with Fisher's ztransformation and then converted into the CCC [24].

Total Deviation Index (TDI)
Total deviation index (TDI) is another common measurement of evaluating agreement between two methods of measurements. TDI is the π₀^th percentile of absolute value of the differences between the methods (μ_d), for a given large probability π₀. TDI always takes a positive value and the smaller value for TDI indicates good agreement. Under the proposed heteroscedastic model, Total deviation index (TDI) between the reader1 and the reader j is defined as follows,

$${\text{τ}}^2\left(v_0\right)={\psi }_1^2+{\psi }_j^2-2{\psi }_{1j}+{\sigma }_1^2{g}_1^2\left(v_0,{\delta }_1\right)+{\sigma }_j^2{g}_j^2\left(v_0,{\delta }_j\right);$$

$$TDI\left(v_0,p_0\right)=\tau \left(v_0\right){\left\{{\chi }_1^2\left(p_0{\left\{\frac{{\beta }_1-{\beta }_j}{\tau \left(v_0\right)}\right\}}^2\right)\right\}}^{1/2};j=2,3,4,5;p_0=0.95$$

Results

The initial model fitting was carried out by using the standard mixed effects model. To check the normality of the residuals for each method, the quantile-quintile plot and the Shapiro- Wilk normality test was used. Table 1 represents the results of the Shapiro-Wilk normality test, which was applied to the residuals of each reader of the mixed effect model. According to the results, it is obvious that the residuals do not follow a normal distribution as all p-values are small. Therefore, the main assumption for the standard mixed effects model (i.e., the residuals have a constant variance) is violated. Moreover, Figure 1 shows separate quantile-quantile plot for each reader when standard mixed effects is fit to the data. The circles cross the line three times indicates that the hump is not the right shape for these data to be normal. These data are therefore not exactly normal. Hence, we model the error variability using two approximations for μ_i, namely mean of the reference method and the BLUP of μ_i.

We consider reader1 as the reference method. As the first approach, let's take the mean vector of the reader1 as the variance covariate and fit the data into a heteroscedastic model.

$${\mu }_{i1}=\stackrel{-}{Y_1}$$

In the dataset there are 40 subjects and each subject has 2 replicates. Since there are 5 readers in the dataset, we need to compare 5 different methods of measurements. Therefore, as for the model stated in the methodology, i=1,2,.....,40; k=1,2; j=1,2,3,4,5.

The Table 2 represents the estimates and standard errors of the fitted heteroscedastic model parameters, which was fitted by using mean of the reference method and BLUP as the variance covariate. Considering the results from the Table 2, both approaches take approximately equal values when compared with the other methods. Estimator of the reader3 has the largest men value, and the smallest estimator is from reader2.

The Likelihood ratio test has been carried out to identify the most appropriate model among the standard mixed effects model and the proposed heteroscedastic model. The test rejected the null hypothesis, which implies that the reduced model is more appropriate, confirming the heteroscedastic model is the best among these two models. Furthermore, the score test is also performed on these models. The result was the same as before confirming that the heteroscedastic model is more suitable for the given dataset. Since the fitted heteroscedastic model is appropriate the next step is to assess the agreement between the methods of measurements.

Agreement evaluation under the Heteroscedastic model
The Figure 2 represents the estimates and their confidence intervals of CCC and TDI using the mean method. The solid line indicates the estimate and the shaded region indicates the estimate±standard error. The dotted line represents the lower bound of the concordance correlation coefficient and the upper bound of the Total deviation index. The upper two graphs represent the CCC with Fisher's z transformation and estimated CCC and lower two graphs represent the TDI with log transformation and estimated TDI values. The CCC gradually decreases from 0.82 to 0.79 through the diameter of 1.45 to 8.5 cm. The TDI value increases from 1.505 to 1.69 with the increase of the magnitude of the measurement. From both CCC and TDI values we can observe that for small diameter values the both methods (readers) have satisfactory agreement, while with the increase of the magnitude of the lung tumor diameter, the agreement decreases. Therefore, it is safe to conclude that the reader2 sufficiently agrees with the reference method (reader1). Likewise, the readings from the readers 3, 4 and 5 derived almost similar results as the reader2 by confirming that all the methods (readers) agree well with the reference method. i.e., reader1.

The Figure 3 represents the CCC and TDI values for the reader1 and reader2, using the BLUP based method. The Figure 2 and the Figure 3 are almost identical. Therefore, reader2 is well agreed with the reader1.

Discussion

The proposed heteroscedastic mixed effects model was fitted by using two methods. One was by using true mean as the variance covariate and the other method is by using the Best Linear Unbiased Predictor (BLUP) as the variance covariate. When fitting the model by using both methods, it was confirmed that the heteroscedastic model is better than the standard mixed effect model. The results from the likelihood ratio test and the score test confirms that conclusion. When considering the means from the Table 2, we can observe that the means from both the approaches are identical if round up to 2 decimal places. This means that the mean method and the BLUP method generate almost identical estimates.

CCC and TDI were used in evaluating the agreement between the reference method (reader1) and the other 4 methods (readers). Both mean based method and the BLUP based method produced almost identical (with slight differences) results for CCC and for TDI. Figures 2 and 3 both show that reader2 agrees well with reader1. Not only the reader2, but also reader3, reader4 and reader5 are also agreed well with the reader1. Moreover, Figures 2 and 3 show that the CCC and TDI values varies with the magnitude. If the data set was modeled using a homoscedastic model, constant values for CCC (0.8069) and TDI (1.6140) can be obtained. Therefore, if a dataset with heteroscedasticity is modeled with standard mixed effects model, the resulting outcome will be misleading and inaccurate. To overcome this issue, the data must be modeled with a model that accounts the heteroscedasticity.

Conclusions

We propose a method successfully in order to fit 5 methods of heteroscedastic clinical measurements and this proposed model can be easily extended to deal with any number of multiple methods which has replicated data. Our approach can accommodate balanced or unbalanced data designs and it works well with any scalar measure of agreement. Moreover, the two model fitting methods discussed give almost identical estimates to the model. The limitation of this study is that the proposed method can be applied only with replicated measurements.

List of abbreviations

BLUP: Best Linear Unbiased Predictor
CCC: Concordance Correlation Coefficient
TDI: Total Deviation Index

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

Authors' contributions	IDD	LSN
Research concept and design	√	√
Collection and/or assembly of data	√	√
Data analysis and interpretation	√	√
Writing the article	√	√
Critical revision of the article	√	√
Final approval of article	√	√
Statistical analysis	√	√

Acknowledgement

The authors are grateful to Professor L. Broemeling for providing the lung tumor size measurements dataset.

Publication history

Editor: Jimmy Efird, East Carolina University, USA.
Received: 18-Oct-2016 Final Revised: 26-Dec-2016
Accepted: 09-Jan-2017 Published: 27-Jan-2017