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

BACKGROUND: 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.
METHODS: 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.
RESULTS: 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.
CONCLUSION: 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 intra-individual 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 intra-individual 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 inter-individual 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.

Methods

Statistical model

Let denote the kth observation of the jth subject by the ith measure, . Let and be and . The repeated measurements for the subject consist of times examinations by the first measure and times examinations by the second measure. If there are missing values, may be different from . Thus, missing values are allowed in the model. In our mixed effect model, is expressed bywhere is a design vector for fixed effects, is a parameter vector of the fixed effects regression coefficients, is the random effect of the subject for the measure and assumed to be normally distributed and are correlated by a covariance and and are assumed to be independent if for , is the residual error and is assumed that . The fixed effect part includes a factor of the measure, where and denote means for the first measure and for the second measure, respectively.

Using notations of vector, the observations of the subject and all observations can be expressed by and , where and , respectively. Let denote a diagonal matrix for the residuals of the subject. It can be seen thatDefine as . The covariance matrix for the random effect of the subject can be expressed by , whereNote that is the same structure for all . Define the matrix as . In the matrix notation, the model (1) can be represented bywhere and are and design matrices for the fixed effects and the random effect of subjects, respectively, is a parameter vector for the random effect and , is a vector of the random residual error and it follows a , and are independent. is the parameter vector mentioned above. We can make full rank. consists of column vectors of and , , where , has 1 for the corresponding elements to the subject by the measure and 0 otherwise.

The expectation and the variance of are expressed byThe components of 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-11].

According to [11], the log-likelihood can be partitioned into two parts, and . The is the part including and not including . The variance and covariance parameters are estimated using an iterative algorithm based on , and fixed effect parameters are estimated based on the asFor the detail of the parameter estimation, see “Appendix 2”.

Evaluation of the difference of intra-individual CVs

As mentioned in statistical model section, and represent means of the first and second measure, respectively. Let denote the difference of CVs between measures, where , . Our objective is to test H0: . Using the delta method [12] and , the expectation and variance of are obtained asUnder H0, has the unit normal distribution by applying the large sample theory, where is the estimator of 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 .

Another method to test H0 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 of the original data set. The resampling will be straightforward if the numbers of observations by measure are the same across the subjects, namely for all . 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 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 , and () 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 with replacement from the subgroup .

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 (see Additional file 1). Summary statistics of age and sex and those of PAS and PV were shown in Tables 1 and 2, respectively.

Table 1

Summary statistics of age and sex

	Mean (SD)/N (%)a	Min–Max
Age	65.4 (8.4)	43–81
Sex
Female	11 (25.6)
Male	32 (74.4)

Sample size: 43 patients

aMean (SD) for age and N (%) for sex

Table 2

Summary statistics of the pain intensity data for 43 patients

	VAS	PV
Meana	21.2	26.1
Variance
Intrab	165.7	1725.1
Interc	477.3	274.9
Individual SDd
Mean (SD)	9.2 (8.7)	21.7 (30.7)
Min–Max	0–33	0–172

VAS and PV were repeatedly measured for 43 patients. The number of data records were 158

aCalculate mean for each of 43 patients (m subj), and calculate mean of 43 m subj values

bUnbiased variance, i.e. , SS is a sum of squares for patient j

cVariance of m subj

dSD for each of 43 patients was calculated

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 and (PV) as least square means (LSMs, [15]), let , , and 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 and denote the proportion of male and mean age, respectively. The MLEs of and are and , where and and is obtained by (4). SEs of LSMs are estimated as . Estimated parameters were shown in Table 2.

To obtain variance of , and should be estimated as shown in (6) in addition to the results in Table 3. Using Fisher’s information matrix we see that and . The difference of CVs and its 95 % confidence interval was estimated as −0.96 (−1.59, −0.33) by (5) and (6). Thus the intra-individual variation was significantly smaller in the VAS than PV.

Table 3

Result of parameter estimation for the pain intensity data

	β_1	β_2	σe12	σe22	σb12	σb22	σb122
MLE	21.0	26.3	166	1707	474	305	423
SE	3.50	4.25	21.8	220.1	117.6	176.2	120.9

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.

Fig. 1

Q–Q plots by measure

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.

Fig. 2

Residual plots by measure

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 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 inter-individual 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 inter-individual 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 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 () 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 () by analyzing if lognormal distributions can be assumed for and the geometric CV [17] of the original data can be approximated by , where is the residual error variance of 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 and and by deriving derivatives of .