Comparison of intra-individual coefficients of variation on the paired sampling data when inter-individual variations are different between measures

Pain intensities of patients are repeatedly measured by Visual Analog Scale (VAS) and Pain Vision (PV) in a clinical research. Two measurements by VAS and PV are performed at the same time. In order to evaluate within patient consistency, intra-individual coefficient of variations (CVs) are compared between measures assuming that the pain status of each patient is stable during the research period. The correlated samples and different inter-individual variation due to different scales of the measures should be taken into account in statistical analysis. The adjustment of covariates will improve the estimation of population mean values of the measures. In this paper, statistical approach to compare the intra-individual CVs is proposed. The approach consists of two steps: (1) estimating population mean values and intra-individual variances of the pain intensities by measure in a mixed effect model framework, (2) computing intra-individual CVs and comparing them between measures. The mixed effect model includes measure and some variables as fixed effects and subject by measure as a random effect. The different inter-individual variations between measures and their covariance reflect the paired sampling in the variance component. The confidence interval of the difference of intra-individual CVs is constructed using the asymptotic normality and the delta method. Bootstrap method is available if sample size is small. The proposed approach is illustrated using pain research data. Measure (VAS and PV), age and sex are included in the model as fixed effects. The confidence intervals of the difference of intra-individual CVs between measures are estimated by the asymptotic theory and by bootstrap using a subgroup resampling, respectively. Both confidence intervals are similar. The proposed approach is useful to compare two intra-individual CVs taking it into account to reflect the paired sampling, different inter-individual variations between measures and some covariates. Although the inclusion of covariates did not improve the goodness-of-fit in the illustration, the proposed model with covariates will improve the accuracy and/or precision if covariates truly influence response variable. This approach can be applicable with small modification to various situations.


Background
A clinical research of oxaliplatin chemotherapy was conducted in colorectal cancer patients. Pain intensities of oxaliplatin-induced peripheral neuropathy were assessed by two measures of Visual Analog Scale (VAS) and Pain Vision PS-2100 (PV). The two pain assessments were performed at the same time on the same set of patients and the assessments were repeated during some time interval. VAS is a widely used subjective scale to measure the pain intensity, and patients assess their pain intensities along a continuous line from 0 to 100. PV is an analytical instrument that is designed to assess sense perception and nociception quantitatively and objectively. The pain intensity by PV is defined as a percentage change of two electrical stimulation values with and without pain. Thus measurements by PV are non-negative values (range of the actual data was 0 to about 400). The objective of the research was to compare these subjective and objective measures, and the intra-individual variation was compared in order to evaluate the within patient consistency as the reliability index assuming that the pain status of each patient has been stable during the research period. Since the scales of two measures are different, the intraindividual coefficient of variation (CV) is used as a variation index.
In this situation, there are two conditions to be considered in the statistical analysis: • correlated samples by paired sampling (assessments by VAS and PV were performed at the same time) • different inter-individual variation between measures (different scales: 0-100 for VAS and non-negative values of 0 to about 400 for PV) Furthermore, there exists an additional requirement to improve the estimation of population mean values of measures: • taking covariates into account, for example sex and/or age In many clinical studies, the standard deviation, the coefficient of variation or a composite index such as the mean amplitude of glycemic excursion (MAGE) in diabetic area are used as intra-individual variation indexes, and those indexes are calculated for each subject and sometimes compared between groups or treatments by two sample t test, linear regression analysis or analysis of variance [e.g. 1-3]. Onishi et al. [4] compared the intraindividual CVs in pre-breakfast self-measured plasma glucose between treatments including some variables as covariates and subject as a random effect. There are many statistical discussions (e.g. 5, 6) relating to comparing intra-individual CVs. Forkman [5] gave a method to compare intra-individual CVs (although "intra-individual" is not explicitly written). Iron [6] compared several methods to calculate intra-class correlation coefficient (ICC). For the paired sampling, Shoukri et al. [7] gave several procedures for testing the equality of two dependent intra-individual CVs. They used a mixed model including one fixed effect and subject as a random effect and they gave a covariance structure reflecting the paired sampling to the model. Their procedures, however, did not clearly model the different inter-individual variation and covariates were not included either. Furthermore, they assumed that the number of repetition of both measurements was the same for all subjects.
In this article, we propose an approach to compare two intra-individual CVs considering the two conditions (correlated samples by paired sampling and different interindividual variation between measures) and adjustment of the covariates. A mixed effect model will be used to estimate population mean values and intra-individual variances for measures. The mixed effect model includes measure and some variables as fixed effects and subject by measure as a random effect. Different inter-individual variation between measures and the paired sampling are taken into account in this model. The number of repetitions may be varied among subjects due to missing. By using iterative algorithms in the mixed effect models, we obtain the estimates without restriction of the same number of repetitions. The intra-individual CVs will be calculated and compared using these estimates.

Statistical model
Let y ijk ≥ 0 denote the kth observation of the jth subject by the ith measure, i = 1, 2, j = 1, 2, . . . , n, k = 1, 2, . . . , m ij . Let m j and m be m j = m 1j + m 2j and m = i,j m ij . The repeated measurements for the jth subject consist of m 1j times examinations by the first measure and m 2j times examinations by the second measure. If there are missing values, m 1j may be different from m 2j . Thus, missing values are allowed in the model. In our mixed effect model, y ijk is expressed by where x ijk is a p × 1 design vector for fixed effects, β = β 1 , . . . , β p ′ is a p × 1 parameter vector of the fixed effects regression coefficients, γ ij is the random effect of the jth subject for the ith measure and assumed to be normally distributed γ ij ∼ N 0, σ 2 bi , γ 1j and γ 2j are correlated by a covariance σ 2 b12 and γ ij and γ i ′ j ′ are assumed to be independent if j � = j ′ for i, i ′ = 1, 2, e ijk is the residual error and is assumed that e ijk ∼ i.i.d.N 0, σ 2 ei . The fixed effect part x ′ ijk β includes a factor of the measure, where β 1 and β 2 denote means for the first measure and for the second measure, respectively.
Using notations of vector, the observations of the jth subject and all observations can be expressed by y j and y, where y j = y 1j1 , . . . , y 1jm 1j , y 2jm 1j +1 , . . . , y 2jm 1j +m 2j ′ and y = y ′ 1 , y ′ 2 , . . . , y ′ n ′ , respectively. Let R j denote a m j × m j diagonal matrix for the residuals of the jth subject. It can be seen that Define R as block diag(R 1 , . . . , R n ). The covariance matrix for the random effect of the jth subject can be expressed by g, where Note that g is the same structure for all j. Define the 2n × 2n matrix G as G = block diag (g, . . . , g). In the matrix notation, the model (1) can be represented by where X and Z are m × p and m × 2n design matrices for the fixed effects and the random effect of subjects, respectively, γ is a 2n × 1 parameter vector for the random effect and γ ∼ N (0, G), e is a m × 1 vector of the random residual error and it follows a N (0, R), e and γ are independent. β is the p × 1 parameter vector mentioned above. We can make X full rank. Z consists of m × 1 column vectors of Z 1j and Z 2j , j = 1, 2, . . . , n, where Z ij , i = 1 2, has 1 for the corresponding elements to the jth subject by the ith measure and 0 otherwise.
The expectation and the variance of y are expressed by The components of V are explained in "Appendix 1".

Parameter estimation based on likelihood
The parameters in the model (2) are estimated based on the usual maximum likelihood (ML) approach for linear mixed models. The covariance matrix is estimated by the method of restricted maximum likelihood (REML) where degrees of freedom is adjusted in estimating fixed effects [8][9][10][11]. According to [11], the log-likelihood can be partitioned into two parts, L 1 and L 2 . The L 1 is the part including β and L 2 not including β. The variance and covariance parameters are estimated using an iterative algorithm based on L 2 , and fixed effect parameters β are estimated based on the L 1 as For the detail of the parameter estimation, see "Appendix 2".

Evaluation of the difference of intra-individual CVs
As mentioned in statistical model section, β 1 and β 2 represent means of the first and second measure, respectively. Let g â denote the difference of CVs between measures, where â = β 1 ,β 2 ,σ 2 e1 ,σ 2 Our objective is to test H 0 : g(a) = 0. Using the delta method [12] and E ∂ 2 L 1 /∂β∂θ l = 0, the expectation and variance of g â are obtained as Under H 0 , g â / Var g â has the unit normal distribution by applying the large sample theory, where Var g â is the estimator of Var g â obtained by substituting estimators for all parameters. Alternatively the equality of CVs between measures can be evaluated by an interval estimation of the difference of CVs using the asymptotic normality of g â / Var g â .
Another method to test H 0 is bootstrap method [13] which is particularly useful if sample size is not sufficiently large to apply the asymptotic theory. A simple bootstrapping is the subject resampling method with replacement. Since the data is obtained based on the paired sampling, the resampling is taken place by subject (kind of block or cluster bootstrap, [14]) in order to capture the dependence structure. The resampled number of subjects is equal to the size n of the original data set. The resampling will be straightforward if the numbers of observations by measure are the same across the subjects, namely m ij = r i , i = 1, 2 for all j = 1, . . . , n . Otherwise, a subgroup resampling classified by number of observations by measure will be one way to make the resample data set with the same records m as in the original data set if two measures (e.g. VAS and PV) are always observed paired although the number of observations may be different among subjects. For example, let n 1 , n 2 and n 3 (n 1 + n 2 + n 3 = n) be defined as the number of subjects repeatedly measured twice, thrice and four times, respectively. Thus there are three subgroups. The resampling is performed in the same size n i with replacement from the ith subgroup (i = 1, 2, 3).
There exist another bootstrap methods such as s parametric bootstrap. Since the objective of this article, however, is not to discuss the bootstrapping, only the blocked and subgrouped way will be used in the later section.

Application result
Our proposed model was applied to the clinical research data mentioned in the background section. Forty-three patients were treated with oxaliplatin chemotherapy. The pain intensity of oxaliplatin-induced peripheral neuropathy was assessed by two measures of VAS and PV. All patients with repeated measurements were included in the analysis, and the number of records of the data set was m = 158 (see Additional file 1). Summary statistics of age and sex and those of PAS and PV were shown in Tables 1 and 2, respectively.
A stopping criterion of numerical iterations by the Fisher's scoring method to estimate parameters was set to 10 −6 . Sex, age, interaction between measure and sex and interaction between measure and age were included in the mixed model as covariates (fixed effects in full model). To obtain the MLEs of β 1 (VAS) and β 2 (PV) as least square means (LSMs, [15]), let β 3 , β 4 , β 5 and β 6 be regression parameters for sex (0: female, 1: male), interaction between measure and sex (1: VAS and male, 0: otherwise), age and interaction between measure and age (age in year for VAS and 0 for PV). Let  Table 2.
To obtain variance of g â , Cov σ 2 e1 ,σ 2 e2 and Cov β 1 ,β 2 should be estimated as shown in (6) in addition to the results in Table 3. Using Fisher's information matrix we see that Cov σ 2 e1 ,σ 2 e2 = −13.2 and  (5) and (6). Thus the intra-individual variation was significantly smaller in the VAS than PV. Bootstrapping was also conducted as sensitivity analysis of the asymptotic inference. 1000 data sets with the same size as the original sample were generated with replacement by the subgroup method, and the model parameters and the difference of CVs between measures were estimated for each data set. A Bootstrap estimate of the difference of intra-individual CVs was obtained as the mean of the 1000 estimates, and a 95 % confidence interval was constructed by mean of the 25th and 26th values and mean of the 975th and 976th values out of the 1000 estimates. The bootstrap estimate was −0.90 (−1.50, −0.31), and this estimate was similar as in the asymptotic theory.

Modeling process and diagnostics
The assumption of normality was checked by Q-Q plots (Fig. 1) of residuals. There were several patients whose pain intensities were close to the maximum pain intensity, and the tails on the right side of distributions of VAS and PV were longer, particularly in PV. Since patients  with high pain intensities cannot be excluded from the analysis population for our research objective and a logarithmically transformation is not possible for zero values, all patients were included in the statistical analyses. We performed a sensitivity analysis excluding two patients with large values and confirmed similar result as shown in the former analysis. The residual plots for fitted value (Fig. 2) also show several large values in PV and the residual errors do not seem to be totally random. The Q-Q plots and residual plots may suggest to use appropriate data transformations, other distributional assumption and/or robust methods. These methods generally improve model fittings and make comparison of means precise. It is, however, not certain if these methods are relevant to compare intra-individual variations. Therefore, we still assume normality in the analysis. Residual plot for each measure, P-P plot by measure and plot of Cook's distance against leverage/(1-leverage) are presented for further diagnostic information in Additional file 2: Figures S1, S2 and S3, respectively.
The goodness-of-fit of the full model, where fixed effects were measure, sex, age, interactions between measure and sex and measure and age, was checked in comparison with the null model including only measure as a fixed effect using a likelihood ratio test. The test statistics was 3.048 and P-value according to χ 2 distribution with four degrees of freedom was 0.550. Thus the addition of sex and age did not contribute to improvement of the model in the pain research data, and the null  . 2 Residual plots by measure model will be applied in practice. Since the objective of this application section is to illustrate the availability of the proposed model, we presented the results of the full model assuming normality. Summary statistics of individual standard deviations (SDs) in VAS and PV were shown in Table 2. The range of individual SDs were vary from 0 to 33 in VAS and from 0 to 172 in PV. The comparison of the intra-individual CV between VAS and PV corresponds to compare the average individual SDs by mean.

Conclusion
Statistical approach to compare two intra-individual CVs considering the paired sampling, different interindividual variations between measures and some covariates was proposed. The approach is (1) estimating intra-individual variances and mean values of measures by the mixed effect model including measure and some covariates as fixed effects and subject by measure as a random effect to reflect different interindividual variation between measures and the paired sampling, (2) computing intra-individual CVs and comparing these CVs between measures. The utility of the approach was shown by applying to the real pain research data. Although the inclusion of covariates did not improve the goodness-of-fit in the illustration, the proposed model with covariates will improve the accuracy and/or precision if covariates truly influence response variable.
The proposed approach is applicable for many situations. One example is to give σ 2 b1 = σ 2 b2 = σ 2 b if the same inter individual variations can be assumed between measures. In a study [16], the frequently-sampled intravenous glucose tolerance test was carried out for subjects by two protocols in order to estimate insulin sensitivity, and intra-individual CVs of the insulin sensitivity were compared between protocols. If our proposed model is applied to this data, the same inter individual variation (σ 2 b1 = σ 2 b2 ) for two protocols will be observed since the same scaled variable was measured by the two different protocols.
Another example is to compare intra-individual variances σ 2 ei(z) (i = 1, 2) by analyzing z ijk = log y ijk if lognormal distributions can be assumed for y ijk and the geometric CV [17] of the original data can be approxi- is the residual error variance of z ijk and it is assumed to be small enough to hold the approximation. The approximation of the geometric CV means that CVs in original scale data can be approximated by variances in the logarithmically transformed data. This procedure suggests to omit the CV calculation step. For example, area under the pharmacokinetic curve and triglyceride have such characteristics. These variables have skewed distributions with longer tail upper side, and are frequently logarithmically transformed before statistical analyses. However, note that possible biases should be discussed before the transformation since the logarithmic transformation can give a larger benefit to a group with large values than the other group with small values when estimating population intra-individual variances.
An extension to multi-group comparisons is straightforward by giving corresponding R j and G j and by deriving derivatives of V .
where I u is a u × u unit matrix and J u×v is a u × v matrix with all the elements equal to 1.

Appendix 2: Parameter estimation
According to [11], the log-likelihood can be partitioned into two parts, L 1 and L 2 , where L 1 is the part including β and L 2 not including β. L 1 and L 2 are expressed as below.
where C is a certain constant and The variance and covariance parameters θ = (θ 1 , θ 2 , θ 3 , θ 4 , θ 5 ) are estimated based on L 2 , while the fixed effect parameters β are estimated based on the L 1 . The first order derivatives of L 2 is ∂V /∂θ l are obtained by simple manipulations as ML estimates θ are the solution of the estimating equations ∂L 2 /∂θ l = 0. Usually ML estimates θ are not expressed explicitly and therefore an iterative algorithm such as Newton-Raphson method is used to obtain θ . Fisher's scoring method, in which Fisher's information matrix I (θ ) is used instead of Hessian in Newton-Raphson method, is now employed as the iterative algorithm, since I(θ ) is a simpler form than Hessian in mixed models and the asymptotic variance of θ can be estimated by I −1 (θ ) (not observed information [18]). Thus second order derivatives are needed to apply the algorithm where θ k is the kth estimates, θ k+1 is an improved one, 0 < c < 1 is a constant value to adjust an improvement quantity and S(θ k ) = (∂L 2 /∂θ l )| θ =θ k . Simple manipulations give second order derivatives of L 2 and I(θ ), where the (θ l , θ l ′ ) element of the former and the latter as can be expressed as respectively. With respect to β, the ML estimator is derived as the solution of the first order derivatives of (7) equal to 0 as θ k+1 = θ k + cI −1 (θ k )S(θ k ), E − ∂ 2 L 2 ∂θ l ∂θ l ′ = 1 2 tr −P ∂V ∂θ l ′ P ∂V ∂θ l , l, l ′ = 1, 2, . . . , 5,