### Likelihood based methodology

Suppose that we are interested in comparing the reproducibility of two instruments. Let *x*
_{
ijl
}be the *jth* measurement of the *ith* subject by the *lth* instrument, *j* = 1,2,... *m*
_{
l
}, *i* = 1,2,... *n*, and *l* = 1, 2. To evaluate the WSCV we consider the one-way random effects model

*x*
_{
ijl
}= *μ*
_{
l
}+ *b*
_{
i
}+ *e*
_{
ijl
}

where *μ*
_{
l
}is the mean value of measurements made by the *lth* instrument, *b*
_{
i
}are independent random subject effects with *b*
_{
i
}~ *N*(0, {\sigma}_{b}^{2}), and *e*
_{
ijl
}are independent *N*(0, {\sigma}_{l}^{2}). Many authors have used the intra-class correlation coefficient (ICC), *ρ*
_{
l
}defined by the ratio {\rho}_{l}={\sigma}_{b}^{2}/\left({\sigma}_{b}^{2}+{\sigma}_{l}^{2}\right) as measure of reproducibility/reliability [18, 23]. Quan and Shih [8] argued that *ρ*
_{
l
}is study-population based since it involves between-subject variation. Meaning that the more heterogeneity in the population, the larger the *ρ*
_{
l
}. Alternatively, they proposed the within-subject coefficient of variation (WSCV) *θ*
_{
l
}= *σ*
_{
l
}/*μ*
_{
l
}as a measure of reproducibility. It determines the degree of closeness of repeated measurements taken on the same subject either by the same instruments or on different occasions under the same conditions. It is clear that, the smaller the WSCV, the better the reproducibility. We distinguish the WSCV from the coefficient of variation C{V}_{l}={\left({\sigma}_{b}^{2}+{\sigma}_{l}^{2}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}/{\mu}_{l} since *CV*
_{
l
}involves {\sigma}_{b}^{2} in the numerator and similar to *ρ*
_{
l
}is population based. Therefore, more heterogeneity in the population would result in a large value of *CV*
_{
l
}. For that reason we shall focus our work on the WSCV rather than the CV. We also note that there is an inverse relationship between the ICC (*ρ*
_{
l
}) and the corresponding within subject variance {\sigma}_{l}^{2}. Clearly, larger values of ICC (higher reliability) would be associated with smaller WSCV (better reproducibility). The focus of this paper is on aspects of statistical inference on the difference between two correlated WSCV. The inferential procedure depends on the multivariate normality of the measurements and is mainly likelihood based. The following set-up is to facilitate the construction of the likelihood function.

Let

{X}_{i}={({X}_{i1},{X}_{i2},\mathrm{.......}{X}_{i{m}_{1}},{X}_{i,{m}_{1}+1},{X}_{i,{m}_{1}+2},\mathrm{....},{X}_{i,{m}_{1}+{m}_{2}})}^{\text{'}}

denote the measurements on the *i*
^{th}subject, *i* = 1,2,....,*n* where {X}_{i1},{X}_{i2},\mathrm{....},{X}_{i{m}_{1}} are the *m*
_{1} measurements obtained by the first method (platform), {X}_{i{m}_{1}+1},{X}_{i{m}_{1}+2},\mathrm{....},{X}_{i{m}_{1}+{m}_{2}} are the *m*
_{2} measurements obtained the second method (platform). We assume that *X*
_{
i
}~ *N*(*μ*, Σ), where {\mu}^{T}=({\mu}_{1}{1}_{{m}_{1}}^{T},{\mu}_{2}{1}_{{m}_{2}}^{T}) and,

\Sigma =\left[\begin{array}{cc}{\sigma}_{1}^{2}{I}_{{m}_{1}}+\frac{{\rho}_{1}}{1-{\rho}_{1}}{\sigma}_{1}^{2}{J}_{{m}_{1}}& {\rho}_{12}{\sigma}_{1}{\sigma}_{2}{J}_{{m}_{1}x{m}_{2}}\\ {\rho}_{12}{\sigma}_{1}{\sigma}_{2}{J}_{{m}_{1}x{m}_{2}}& {\sigma}_{2}^{2}{I}_{{m}_{2}}+\frac{{\rho}_{2}}{1-{\rho}_{2}}{\sigma}_{2}^{2}{J}_{{m}_{2}}\end{array}\right]

(2)

In these expressions 1_{
k
}is a column vector with all *k* elements equal to 1, *I*
_{
k
}is a *k* × *k* identity matrix and *J*
_{
k
}and *J*
_{
kxt
}are *k* × *k* and *k* × *t* matrices with all the elements equal to 1. Thus the model assumes that the *m*
_{1} observations taken by the first platform have common mean *μ*
_{1}, common variance {\sigma}_{1}^{2}, and common intra-class correlation *ρ*
_{1}, whereas the *m*
_{2} measurements taken by the second platform have common mean *μ*
_{2}, common variance {\sigma}_{2}^{2}, and common intra-class correlation *ρ*
_{2}. Moreover, *ρ*
_{12} denotes the interclass correlation between any pair of measurements *x*
_{
ij
}(*j* = 1,2,... *m*
_{1}) and {x}_{im1+t}\left(t=1,2,\dots {m}_{2}\right), and also assumed constant across all subjects in the population.

For the *l*
^{th}method, the WSCV, which will be denoted as *θ*
_{
l
}in the remainder of the paper is defined as

*θ*
_{
l
}= *σ*
_{
l
}/*μ*
_{
l
}, *l* = 1, 2.

Our primary aim is to develop and evaluate methods of testing *H*
_{0}:*θ*
_{1} = *θ*
_{2} taking into account dependencies induced by a positive value of *ρ*
_{12}. We restrict our evaluation to reproducibility studies having *m*
_{1} = *m*
_{2} = *m*.

### Methods for testing the null hypothesis

#### Wald test (WT)

If *X*
_{1}, *X*
_{2},.... *X*
_{
n
}is a sample from the above multivariate normal distribution, then the *log*-likelihood function *l*, as a function of *ψ* = (*μ*
_{1}, *μ*
_{2}, {\sigma}_{1}^{2}, {\sigma}_{2}^{2}, *ρ*
_{1}, *ρ*
_{2}, *ρ*
_{12}) is given by:

-2L=Q+nm\mathrm{log}\left({\sigma}_{1}^{2}{\sigma}_{2}^{2}\right)-n\mathrm{log}\left(\left(1-{\rho}_{1}\right)\left(1-{\rho}_{2}\right)\right)+n\mathrm{log}w

(3)

where,

*w* = *u*
_{1}
*u*
_{2} - *m*
^{2} {\rho}_{12}^{2},

*u*
_{
l
}= 1 + (*m* - 1)*ρ*
_{
l
}, *l* = 1, 2 and,

\begin{array}{c}Q=\frac{{S}_{1}^{2}}{{\sigma}_{1}^{2}}+\frac{m\left(1-{\rho}_{1}\right){u}_{2}}{w{\sigma}_{1}^{2}}{\displaystyle \sum _{i=1}^{n}{\left({\overline{x}}_{i1}-{\mu}_{1}\right)}^{2}}+\frac{{S}_{2}^{2}}{{\sigma}_{2}^{2}}+\frac{m\left(1-{\rho}_{2}\right){u}_{1}}{w{\sigma}_{2}^{2}}{\displaystyle \sum _{i=1}^{n}{\left({\overline{x}}_{i2}-{\mu}_{2}\right)}^{2}}\\ -\frac{2{m}^{2}{\rho}_{12}}{w{\sigma}_{1}{\sigma}_{2}}{\left(\left(1-{\rho}_{1}\right)\left(1-{\rho}_{2}\right)\right)}^{1/2}{\displaystyle \sum _{i=1}^{n}\left({\overline{x}}_{i1}-{\mu}_{1}\right)\left({\overline{x}}_{i2}-{\mu}_{2}\right)}\end{array}

From [24] the conditions {1 + (*m* - 1)*ρ*
_{1}}{1 + (*m* - 1)*ρ*
_{2}} > *m*
^{2} {\rho}_{12}^{2} and -1/(*m* - 1) <*ρ*
_{
l
}< 1 must be satisfied for the likelihood function to be a sample from a non-singular multivariate normal distribution.

The summary statistics given in (3) are defined as:

\begin{array}{ll}{\overline{x}}_{ij}={\displaystyle \sum _{k=1}^{m}{x}_{ijk}/m}\hfill & i=1,2,\dots n;j=1,2\hfill \\ {S}_{j}^{2}={\displaystyle \sum _{i=1}^{n}{\displaystyle \sum _{k=1}^{m}{\left({x}_{ijk}-{\overline{x}}_{ij}\right)}^{2}}}\hfill \end{array}

The maximum likelihood estimates (MLE) for *μ*
_{
l
}and {\sigma}_{l}^{2} are given respectively by {\stackrel{\u2322}{\mu}}_{l}={\overline{x}}_{l},{\stackrel{\u2322}{\sigma}}_{l}^{2}={S}_{l}^{2}/n\left(m-1\right), where {\overline{x}}_{l}=\frac{1}{n}{\displaystyle \sum _{i=1}^{n}{\overline{x}}_{ij}} and *l* = 1, 2. Clearly, {\widehat{\sigma}}_{l}^{2} exists for values of *m* > 1. Therefore we shall assume that *m* > 1 throughout this paper. From [24], we obtain {\widehat{\rho}}_{1} and {\widehat{\rho}}_{2} by computing Pearson's product-moment correlation over all possible pairs of measurements that can be constructed within platforms 1 and 2 respectively, with {\widehat{\rho}}_{12} similarly obtained by computing this correlation over the *nm*
^{2} pairs (*x*
_{
ij
}, *x*
_{
i
},_{
m+l
}).

The WT of *H*
_{0}:*θ*
_{1} = *θ*
_{2} requires the evaluation of variance of {\stackrel{\u2322}{\theta}}_{l}, *l* = 1, 2, and \mathrm{cov}\left({\stackrel{\u2322}{\theta}}_{1},{\stackrel{\u2322}{\theta}}_{2}\right). To obtain these values we use elements of Fisher's information matrix, along with the delta method [26, 27]. On writing:

*ψ* = (*ψ*
_{1}, *ψ*
_{2})',*ψ*
_{1} = (*μ*
_{1}, *μ*
_{2})', *and* {\psi}_{2}={\left({\sigma}_{1}^{2},{\sigma}_{2}^{2},{\rho}_{1},{\rho}_{2},{\rho}_{12}\right)}^{\prime}, the Fisher's information matrix *I* = -*E*⌊∂^{2}
*l*/∂*ψ*∂*ψ*'⌋ has the following structure:

I=\left[\begin{array}{cc}{I}_{11}& O\\ O& {I}_{22}\end{array}\right].

(4)

This is based on a result from [26] (page 239) indicating that, *I*
_{12} = {I}_{21}^{\text{'}} = -*E*(∂^{2}
*l*/∂*ψ*
_{1}∂*ψ*'_{2}) = 0. Therefore, from the asymptotic theory of maximum likelihood estimation we have:

{I}_{11}^{-1}=\left[\begin{array}{cc}\mathrm{var}\left({\stackrel{\u2322}{\mu}}_{1}\right)& \mathrm{cov}\left({\stackrel{\u2322}{\mu}}_{1},{\stackrel{\u2322}{\mu}}_{2}\right)\\ \mathrm{cov}\left({\stackrel{\u2322}{\mu}}_{1},{\stackrel{\u2322}{\mu}}_{2}\right)& \mathrm{var}\left({\stackrel{\u2322}{\mu}}_{2}\right)\end{array}\right]

And the elements of *I*
_{22} are given in the Appendix.

The elements of {I}_{22}^{-1} are the asymptotic variance- covariance matrix of the maximum likelihood estimators of the covariance parameters. Inverting Fisher's information matrices we get:

\mathrm{var}\left({\stackrel{\u2322}{\mu}}_{l}\right)=\frac{{\sigma}_{l}^{2}}{nm\left(1-{\rho}_{l}\right)}\left[1+\left(m-1\right){\rho}_{l}\right].

(5)

Applying the delta method [27], we can show, to the first order of approximation that:

\begin{array}{cc}\mathrm{var}\left({\stackrel{\u2322}{\sigma}}_{l}\right)\approx {\sigma}_{l}^{2}/2n\left(m-1\right).& l=1,2\end{array}

(6)

The maximum likelihood estimator of *θ*
_{
l
}is {\widehat{\theta}}_{l}=\frac{{\widehat{\mu}}_{l}}{{\widehat{\sigma}}_{l}}. Again, by application of the delta method, we can show to the first order of approximation that:

\mathrm{var}\left({\stackrel{\u2322}{\theta}}_{l}\right)\approx \frac{{\theta}_{l}^{4}\left[1+\left(m-1\right){\rho}_{l}\right]}{nm\left(1-{\rho}_{l}\right)}+\frac{{\theta}_{l}^{2}}{2n\left(m-1\right)},

(7)

as was shown by Quan and Shih [8].

Again using the delta method we show approximately that:

\mathrm{cov}\left({\stackrel{\u2322}{\theta}}_{1},{\stackrel{\u2322}{\theta}}_{2}\right)\approx \frac{2{\theta}_{1}^{2}{\theta}_{2}^{2}{\rho}_{12}}{n\sqrt{\left(1-{\rho}_{1}\right)\left(1-{\rho}_{2}\right)}}.

(8)

From [28] we apply the large sample theory of maximum likelihood to establish that:

Z=\frac{{\widehat{\theta}}_{1}-{\widehat{\theta}}_{2}}{\sqrt{\mathrm{var}({\widehat{\theta}}_{1})+\mathrm{var}({\widehat{\theta}}_{2})-2\mathrm{cov}({\widehat{\theta}}_{1},{\widehat{\theta}}_{2})}}

(9)

is approximately distributed under *H*
_{0} as a standard normal deviate. The denominator of Z is the standard error of {\widehat{\theta}}_{1}-{\widehat{\theta}}_{2} and is denoted by *SE* {\widehat{\theta}}_{1}-{\widehat{\theta}}_{2}. Since the standard error of {\widehat{\theta}}_{1}-{\widehat{\theta}}_{2} contains unknown parameters, its maximum likelihood estimate \widehat{S}E({\widehat{\theta}}_{1}-{\widehat{\theta}}_{2}) is obtained by substituting {\widehat{\theta}}_{l} for *θ*
_{
l
}, {\tilde{\rho}}_{l} for *ρ*
_{
l
}and {\tilde{\rho}}_{12} for *ρ*
_{12}. Moreover, we may construct an approximate (1-*α*)100% confidence interval on (*θ*
_{1} - *θ*
_{2}) given as:

{\widehat{\theta}}_{1}-{\widehat{\theta}}_{2}\pm {z}_{\alpha /2}\widehat{S}E({\widehat{\theta}}_{1}-{\widehat{\theta}}_{2}), where *z*
_{
α/2 }is the (1-*α*/2)100% cut-off point of the standard normal distribution.

#### Likelihood ratio test (LRT)

An LRT of *H*
_{0} : *θ*
_{1} = *θ*
_{2} was developed numerically, and computed by first setting *μ*
_{
l
}= *σ*
_{
l
}/*θ*
_{
l
}, *l* = 1,2 in Equation (3), and then adopting the following algorithm:

1- Set *μ*
_{
l
}= *σ*
_{
l
}/*θ*
_{
l
}, *l* = 1,2 in Equation (3), thereafter;

2- Set *θ*
_{1} = *θ*
_{2} = *θ* in (3)

3- Minimize the resulting expression with respect to all six parameters (*σ*
_{1}, *σ*
_{2}, *ρ*
_{1}, *ρ*
_{2}, *ρ*
_{12}, *θ*) and;

4- Subtract the minimum from the minimum of -2*L* as computed over all seven parameters (*σ*
_{1}, *σ*
_{2}, *ρ*
_{1}, *ρ*
_{2}, *ρ*
_{12}, *θ*
_{1}, *θ*
_{2}) in the model.

It then follows from standard likelihood theory that the resulting test statistic is approximately chi-square distributed with 1 degree of freedom under H_{0}.

#### Score test

One of the advantages of likelihood based inference procedure is that in addition to the WT and the LRT "Rao's score test" can also be readily developed. The motivation for it is that it can sometimes be easier to maximize the likelihood function under the null hypothesis than under the alternative hypothesis. A standard procedure for performing the score test of *H*
_{0} : *θ*
_{1} = *θ*
_{2} is to set *θ*
_{2} = *θ*
_{1} + Δ, so that the null hypothesis is equivalent to *H*
_{0} : Δ = 0, where Δ is unrestricted. Replacing *μ*
_{
l
}by *σ*
_{
l
}/*θ*
_{
l
}, the log-likelihood function *L* is then independent of *μ*
_{
l
}.

Let *L* = *L*(Δ; *ψ*
^{•}) = *L*(Δ; *θ*
_{1}, *σ*
_{1}, *σ*
_{2}, *ρ*
_{1}, *ρ*
_{2}, *ρ*
_{12}) and {l}_{1}=\frac{\partial L}{\partial \Delta},{l}_{2}=\frac{\partial L}{\partial {\psi}^{\ast}}.

From [28] the score statistic is given by:

S={\dot{l}}_{1}^{T}{A}_{1\u20222}^{-1}{\dot{l}}_{1},

where

{l}_{1}^{\u2022}=\frac{\partial L}{\partial \Delta}\underset{\Delta =0}{|}=\frac{nm}{\widehat{w}{\widehat{\theta}}_{1}^{2}}\left[{\widehat{\mu}}_{1}\left(1-{\widehat{\rho}}_{2}\right)\left(\frac{{\widehat{\theta}}_{1}-{\widehat{\theta}}_{2}}{{\widehat{\theta}}_{1}{\widehat{\theta}}_{2}}\right)\right]

(10)

and {A}_{1\u20222}={A}_{11}-{A}_{12}{A}_{22}^{-1}{A}_{21}. The matrices on the right hand side of *A*
_{1•2} are obtained from partitioning the Fisher's information matrix A so that A=\left(\begin{array}{cc}{A}_{11}& {A}_{12}\\ {A}_{21}& {A}_{22}\end{array}\right) where {A}_{11}=E\left(-\frac{{\partial}^{2}l}{\partial {\Delta}^{2}}\right),{A}_{12}={A}_{21}^{T}=E\left(-\frac{{\partial}^{2}l}{\partial \Delta \partial {\psi}^{\ast}}\right), and {A}_{22}=E\left(-\frac{{\partial}^{2}l}{\partial {\psi}^{\ast}\partial {\psi}^{\ast T}}\right) with all the matrices on the right hand side of *A*
_{1•2} evaluated at Δ = 0. When an estimator other than the MLE is used for the nuisance parameters *ψ**, provided that the estimator {\widehat{\psi}}^{\ast} is \sqrt{n} consistent, it was shown that the asymptotic distribution of *S* is that of a chi-square with 1 degree of freedom [29, 30].

The score test has been applied in many situations and has been proven to be locally powerful. Unfortunately, the inversion of *A*
_{1•2} is quite complicated and we cannot obtain a simple expression for *S* that can be easily used. Moreover, we have also found through extensive simulations that while the score test holds its levels of significance, it is less powerful than LRT and WT across all parameter configurations. We therefore focus our subsequent discussion of power to LRT and WT.

#### Regression test

Pitman [1] and Morgan [2] introduced a technique to test the equality of variances of two correlated normally distributed random variables. It is constructed to simply test for zero correlation between the sums and differences of the paired data. Bradley and Blackwood [31] extended Pitman and Morgan's idea to a regression context that affords a simultaneous test for both the means and the variances. The test is applicable to many paired data settings, for example, in evaluating the reproducibility of lab test results obtained from two different sources. The test could also be used in repeated measures experiments, such as in comparing the structural effects of two drugs applied to the same set of subjects. Here we generalize the results of Bradley and Blackwood to establish the simultaneous equality of means and variances of two correlated variables, implying the equality of their coefficients of variations.

Let {\overline{X}}_{ij}={\displaystyle \sum _{k=1}^{m}{X}_{ijk}}/m, and define {d}_{i}={\overline{X}}_{i1}-{\overline{X}}_{i2}, and {s}_{i}={\overline{X}}_{i1}+{\overline{X}}_{i2}.

Direct application of the multivariate normal theory shows that the conditional expectation of *d*
_{
i
}on *s*
_{
i
}is linear [32]. That is

*E*(*d*
_{
i
}| *s*
_{
i
}) = *α* + *βs*
_{
i
},

where

\alpha =\left({\mu}_{1}-{\mu}_{2}\right)-\left({\mu}_{1}+{\mu}_{2}\right)\left({\sigma}_{1}^{2}-{\sigma}_{2}^{2}\right)k

(11.a)

\beta =\left({\sigma}_{1}^{2}-{\sigma}_{2}^{2}\right)k

(11.b)

where

{k}^{-1}={\sigma}_{1}^{2}{\left(1-{\rho}_{1}\right)}^{-1}\left(1+\left(2m-1\right){\rho}_{1}\right)+{\sigma}_{2}^{2}{\left(1-{\rho}_{2}\right)}^{-1}\left(1+\left(2m-1\right){\rho}_{2}\right) is strictly positive.

The proof is straightforward and is therefore omitted.

It can be shown then from direct application of the multivariate normal theory that the conditional expectation (11) is linear, and does not depend on the parameter *ρ*
_{12}.

From (11.a) and (11.b), it is clear that *α* = *β* = 0 if and only if *μ*
_{1} = *μ*
_{2} and *σ*
_{1} = *σ*
_{2} simultaneously. Therefore, testing the equality of two correlated coefficients of variations is equivalent to testing the significance of the regression equation (11). From the theory of least squares, if we define:

\text{TSS}={\displaystyle \sum _{i=1}^{n}{\left({d}_{i}-\overline{d}\right)}^{2}},RSS={\widehat{\beta}}_{1}^{2}{\displaystyle \sum _{i=1}^{n}{\left({s}_{i}-\overline{s}\right)}^{2}} and EMS = (*TSS* - *RSS*)/(*n* - 2),

the hypothesis *H*
_{0} : *α* = *β* = 0 is rejected when *RSS/EMS* exceeds *F*
_{
v,1,(n-2)}, the (1 - *v*) 100% percentile value of the *F*-distribution with 1 and *(n-2)* degrees of freedom [32].