José Moral de la Rubia
School of Psychology, Universidad Autónoma de Nuevo León, Monterrey, Mexico.
DOI: 10.4236/ojs.2026.162004   PDF    HTML   XML   64 Downloads   312 Views   Citations

Abstract

Measures of interindividual variation are rarely used or studied in descriptive/inference statistics, despite constituting a third approach to describing data variability, alternative to ranges (the difference between two extreme measures of position) and to the mean or median of the distances of the data from a measure of central tendency. The aim of this article is to present point and interval estimation (asymptotic and bootstrap) for four absolute measures of variation based on interindividual differences and one relative measure, to develop an R script for their computation, and to demonstrate the performance of the script using two examples. One example is based on a large sample (n = 1000) drawn from a bounded continuous distribution with negative skewness (PERT (0, 8, 10)), which may simulate data from a knowledge–ability test; the other is based on a medium-sized sample (n = 100) drawn from a bounded discrete distribution with positive skewness (BN (n = 10, p = 0.25)), which may simulate the recording (present/absent) of a behavior in independent observations. Both distributions are slightly platykurtic and non-normal. With the measure based on the median of interindividual differences, difficulties arise when generating its sampling distribution for discrete data, but not for continuous data. The absolute measures converge to a normal distribution, whereas the relative measure does not. It is concluded that the script enables the appropriate computation of interindividual variation statistics. Its use is recommended, and further investigation of these measures through simulation studies is encouraged.

Keywords

Descriptive Measures, Variability, Interindividual Differences, Confidence Intervals, Sampling Distribution

Share and Cite:

1. Introduction

The absolute measures of variation describe the dispersion of data and are expressed in the units of the variable. They can be grouped into three types. The first set is based on the distance between two extreme points, defined by order statistics: minimum (min), maximum (max), and percentiles (P10, P25, P75, and P90) [1]. These are known as ranges (R), which include the total range (R = max − min), the percentile range (PR = P90 − P10), the interquartile range (IQR = P75 − P25), and the semi-interquartile range (SIQR = (P75 − P25)/2).

A second set is based on the mean or median of the absolute or squared distances of the scores from a measure of central tendency such as the mean, median, or mode [2]. This group includes the variance (S²), the standard deviation (S), the mean absolute deviation about the mean (MAD), the median absolute deviation (MEAD), and the standard deviation from the mode (SDM). The latter was proposed by Kvalseth [3].

A third set of absolute variation measures is based on the mean or median of the differences between each sample data point and all other data points in the sample, also known as interindividual differences [4]. Given a sample of n data points, the n × (n − 1)/2 positive pairwise differences are computed. If their mean is calculated, the resulting statistic is the Mean Interindividual Difference (MID). If the median is computed, the result is the Median Interindividual Difference (MEID). If the differences are squared and their mean is taken, the resulting measure is the Interindividual Variance (IVAR). Taking the square root of the interindividual variance yields the Interindividual Standard Deviation (ISD). These measures typically produce values that are greater than or equal to those obtained from variation measures based on the distances from a central tendency statistic. A correspondence can also be established between the two groups: MID and MAD, MEID and MEAD, IVAR and S², as well as ISD and S.

Regarding notation, these measures are written in uppercase when used in descriptive statistics or when referring to the population parameter or value, and in lowercase or uppercase with a circumflex accent when used in inferential statistics to denote the point estimate obtained from a sample of size n [5].

2. Theoretical Foundations

2.1. Definition of the Measures of Interindividual Variation

The procedure to calculate the variation measures based on interindividual differences begins by computing these differences. To do this, the data are first sorted in ascending order, and the absolute difference between each value and all subsequent values in the sequence yields the n × (n − 1)/2 interindividual differences (Equation (1)). The cardinality, or number of elements, in this set corresponds to the number of combinations without repetition of the n sample values taken two at a time.

$x={\left\{x_i\right\}}_{i=1}^n=x_1,x_2,\dots ,x_n$

$x_{\left(1\right)}\le x_{\left(2\right)}{x}_{\left(1\right)}\le \cdots \le x_{\left(n\right)}$

$\begin{array}{c}di\left(x\right)={\left\{\left|x_{\left(i\right)}-x_{\left(i^{\prime }\right)}\right|\right\}}_{i=1i^{\prime }>1}^{n-1n}\\ =\left\{\frac{1\left(n-1\right)/2}{\left|x_{\left(1\right)}-x_{\left(2\right)}\right|,\left|x_{\left(1\right)}-x_{\left(3\right)}\right|,\left|x_{\left(1\right)}-x_{\left(4\right)}\right|,\dots ,\left|x_{\left(n-1\right)}-x_{\left(n\right)}\right|}\right\}\end{array}$ (1)

Second, the arithmetic mean of the interindividual differences is calculated to obtain the Mean Interindividual Difference (MID). See Equation (2).

$MID=\frac{{\displaystyle \sum {}_{i=1}^{n-1}{\displaystyle \sum {}_{i^{\prime }>i}^n\left|x_{\left(i\right)}-x_{\left(i^{\prime }\right)}\right|}}}{{}_n{C}_2}=\frac{{\displaystyle \sum {}_{i=1}^{n-1}{\displaystyle \sum {}_{i^{\prime }>i}^n\left|x_{\left(i\right)}-x_{\left(i^{\prime }\right)}\right|}}}{n\left(n-1\right)/2}$ (2)

Third, the median of the interindividual differences is computed to obtain the Median Interindividual Difference (MEID). See Equation (3).

$\begin{array}{c}MEID=mdn\left({\left\{\left|x_{\left(i\right)}-x_{\left(i^{\prime }\right)}\right|\right\}}_{\begin{array}{cc}i= 1& i^{\prime }>i\end{array}}^{\begin{array}{cc}n-1& n\end{array}}\right)\\ =mdn\left(\stackrel{n\left(n-1\right)/2}{\overbrace{\left|x_{\left(1\right)}-x_{\left(2\right)}\right|,\left|x_{\left(1\right)}-x_{\left(3\right)}\right|,\left|x_{\left(1\right)}-x_{\left(4\right)}\right|,\dots \left|x_{\left(n-1\right)}-x_n\right|}}\right)\end{array}$ (3)

Fourth, the arithmetic mean of the squared interindividual differences is calculated to obtain the Interindividual Variance (IVAR). See Equation (4).

$IVAR=\frac{{{\displaystyle \sum }}_{i=1}^{n-1}{{\displaystyle \sum }}_{i^{\prime }>i}^n{\left(x_{\left(i\right)}-x_{\left(i^{\prime }\right)}\right)}^2}{n\left(n-1\right)/2}$ (4)

It can be shown algebraically that IVAR is equivalent to twice the sample variance of the n scores (Equation (5)).

$\begin{array}{c}IVAR=\frac{2{{\displaystyle \sum }}_{i=1}^{n-1}{{\displaystyle \sum }}_{i^{\prime }>i}^n\left(x_{\left(i\right)}^2+x_{\left(i^{\prime }\right)}^2-2x_{\left(i\right)}{x}_{\left(i^{\prime }\right)}\right)}{n\left(n-1\right)}\\ =\frac{{{\displaystyle \sum }}_{i=1}^n{{\displaystyle \sum }}_{i^{\prime }=1}^n\left(x_{\left(i\right)}^2+x_{\left(i^{\prime }\right)}^2-2x_{\left(i\right)}{x}_{\left(i^{\prime }\right)}\right)}{n\left(n-1\right)}\\ =\frac{{{\displaystyle \sum }}_{i=1}^n{{\displaystyle \sum }}_{i^{\prime }=1}^n{x}_{\left(i\right)}^2+{{\displaystyle \sum }}_{i=1}^{n-1}{{\displaystyle \sum }}_{i^{\prime }=1}^n{x}_{\left(i^{\prime }\right)}^2-2{{\displaystyle \sum }}_{i=1}^n{{\displaystyle \sum }}_{i^{\prime }=1}^n{x}_{\left(i\right)}{x}_{\left(i^{\prime }\right)}}{n\left(n-1\right)}\\ =\frac{n{{\displaystyle \sum }}_{i=1}^n{x}_{\left(i\right)}^2+n{{\displaystyle \sum }}_{i^i=1}^n{x}_{\left(i^{\prime }\right)}^2-2{{\displaystyle \sum }}_{i=1}^n{x}_{\left(i\right)}{{\displaystyle \sum }}_{i^i=1}^n{x}_{\left(i^{\prime }\right)}}{n\left(n-1\right)}\\ =\frac{2\left(n{{\displaystyle \sum }}_{i=1}^n{x}_{\left(i\right)}^2-{\left({{\displaystyle \sum }}_{i=1}^n{x}_{\left(i\right)}\right)}^2\right)}{n\left(n-1\right)}=\frac{2\left(n{{\displaystyle \sum }}_{i=1}^n{x}_{\left(i\right)}^2-{\left(n\overline{x}\right)}^2\right)}{n\left(n-1\right)}\\ =\frac{2n\left({{\displaystyle \sum }}_{i=1}^n{x}_{\left(i\right)}^2-n{\overline{x}}^2\right)}{n\left(n-1\right)}=2\frac{{{\displaystyle \sum }}_{i=1}^n{x}_{\left(i\right)}^2-n{\overline{x}}^2}{n-1}\\ =2\times s_{n-1}^2\left(x\right)\end{array}$ (5)

Fifth, the square root of IVAR is taken to obtain the Interindividual Standard Deviation (ISD), and it can be shown algebraically that ISD is equivalent to sample standard deviation (of the n scores) multiplied by squared root of two (Equation (6)).

$ISD=\sqrt{IVAR}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^{n-1}{{\displaystyle \sum }}_{i^{\prime }>i}^n{\left(x_{\left(i\right)}-x_{\left(i^{\prime }\right)}\right)}^2}{n\left(n-1\right)/2}}=\sqrt{2s_{n-1}^2}=\sqrt{2}\times s_{n-1}\left(x\right)$ (6)

The absolute variation measures based on interindividual differences share several properties: they are expressed in the unit of the variable, they are always positive, and larger values indicate greater variability. A value of 0 corresponds to a constant random variable (i.e., a sample in which the same value was obtained in all n observations of variable X). Like other absolute measures of variation, they are invariant to changes in origin but are affected by changes in scale in the values of X. Under a linear transformation, interindividual variation (IVAR) is multiplied by the square of the scale (slope) constant, whereas the other three statistics (MID, MEID, and ISD) are multiplied by the absolute value of that constant.

A measure of interindividual relative variation can be defined by adapting the modification that Merce [6] proposed for the Pearson coefficient of variation [7]. The ratio of the interindividual standard deviation to the total range, multiplied by 100, yields the Interindividual Coefficient of Variation (ICV). Precisely, it corresponds to the Merce coefficient of variation multiplied by the square root of two (Equation (7)).

$\begin{array}{c}ICV=100\times \frac{ISD}{Range\left(x\right)}=100\times \frac{\sqrt{2}\times s_{n-1}}{\max \left(x\right)-\min \left(x\right)}\\ =100\times \sqrt{2}\frac{s_{n-1}}{x_{\left(n\right)}-x_{\left(1\right)}}=100\times \sqrt{2}\times C{V}_{Merce}\end{array}$ (7)

For a constant variable, the ICV is indeterminate because the sample range, which appears in the denominator, equals zero. The infimum of ICV is 0, which are approaching degenerate variables. As an illustration, a sample consisting of 4999 identical values and a single different observation yields an ICV close to 2, reflecting very low relative interindividual variability. At the sample level, ICV values depend on both the distributional form and sampling variability. For instance, a sample of 5000 observations drawn from a standard Cauchy distribution (seed = 123) yields an ICV of 2.09. This value does not characterize the underlying distribution because it lacks fine moments; however, it illustrates a very low value.

Higher ICV values are typically observed in samples where the data points cluster near the boundaries of the support, particularly in highly bimodal or discrete distributions with many repeated values. For bounded variables, the ratio of the standard deviation to the range is bounded and has a possible maximum. Extremely bimodal discrete distributions with all probability mass concentrated at the boundaries, such as a symmetric Bernoulli distribution, provide an extreme benchmark for which the ICV reaches $100\sqrt{2}$ , or approximately 70.7.

Among continuous distributions supported on a finite interval, the maximum value is 50, reached by the arcsine distribution, which corresponds to the Beta (1/2, 1/2) distribution. For unbounded distributions, no finite upper bound exists for the ICV at the sample level, as the sample range may become arbitrarily small while the sample standard deviation remains positive.

For samples drawn from a standard normal distribution, the ICV is usually less than one-third and converges toward approximately 23 for large sample sizes. These reference values, together with the theoretical benchmarks described above for bounded distributions, provide a practical basis for interpreting the magnitude of the ICV. Accordingly, the following heuristic interpretation bands are proposed: very low variability for ICV values in the interval (0, 10), low for values in [10, 20), medium for values in [20, 40), high for values in [40, 50), and very high for values equal to or greater than 50. These categories should be regarded as descriptive guidelines rather than strict statistical thresholds, as is common in the interpretation of coefficients of variation and other relative dispersion measures in applied statistical practice [8].

2.2. Asymptotic Standard Error and Wald-Type Confidence Intervals

2.2.1. Mean Interindividual Difference (MID)

The distribution of the interindividual differences is non-negative and typically right-skewed, with support determined by the range of the original variable. If the underlying variable has finite moments, the absolute differences also have finite moments. Although the distribution of a single absolute difference does not converge to normality, statistics based on averages of such differences—such as the MID, a second-order non-degenerate U-statistic—are asymptotically normal under very general conditions, holding for a wide class of underlying distributions, e.g., logistic, gamma, beta, Laplace, and binomial [9].

As noted in the previous paragraph, MID is a symmetric, second-order U-statistic and is non-degenerate; that is, it is an unbiased estimator of the functional parameter θ = E [h (X₁, X₂)], defined from a symmetric second-order kernel (h (x_i, x_j) = |x_i – x_j|), whose firstorder projection has non-zero variance [10]. Consequently, it satisfies Hoeffding’s Central Limit Theorem [9] [11], which allows for the definition of an asymptotic standard error (Equation (8)), from which a Waldtype confidence interval can be obtained to estimate MID (Equation (9)).

$\begin{array}{l}x=\left\{x_1,x_2,\dots ,x_n\right\}\\ \widehat{MID}=\frac{{{\displaystyle \sum }}_{i=1}^{n-1}{{\displaystyle \sum }}_{i^{\prime }>i}^n\left|x_{\left(i\right)}-x_{\left(i^{\prime }\right)}\right|}{n\left(n-1\right)/2}\\ {\widehat{h}}_1={\left\{{\widehat{h}}_1\left(x_i\right)\right\}}_{i=1}^n={\left\{\frac{1}{n-1}{\displaystyle \sum }_{i\ne i^{\prime }}\left|x_i-x_{i^{\prime }}\right|-\widehat{MID}\right\}}_{i=1}^n\\ ase\left(\widehat{MID}\right)=2\frac{s_{n-1}\left({\widehat{h}}_1\right)}{\sqrt{n}}=2\times \sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{\left({\widehat{h}}_1\left(x_i\right)-\overline{{\widehat{h}}_1}\right)}^2}{n\left(n-1\right)}}\\ \overline{{\widehat{h}}_1}=\frac{{{\displaystyle \sum }}_{i=1}^n{\widehat{h}}_1\left(x_i\right)}{n}\end{array}$ (8)

$P\left(\widehat{MID}-z_{1-\frac{\alpha }{2}}ase\left(\widehat{MID}\right)\le MID\le \widehat{MID}+z_{1-\frac{\alpha }{2}}ase\left(\widehat{MID}\right)\right)=1-\alpha$ (9)

2.2.2. Median Interindividual Difference (MEID)

The MEID is a symmetric second-order U-quantile, defined as the sample quantile of the kernel $h\left(x_i,x_{i^{\prime }}\right)=\left|x_i-x_{i^{\prime }}\right|$ computed over all unordered pairs of observations [11]. Whereas U-statistics provide unbiased estimators of population parameters by averaging a kernel function across all possible subsets of the sample, U-quantiles extend this framework to robust estimation by replacing the average with a sample quantile of the kernel distribution, thereby reducing sensitivity to extreme values. U-quantiles are particularly suitable for data exhibiting weak forms of dependence, under which their asymptotic properties remain valid [12].

Once the estimator class to which the MEID belongs has been established, its asymptotic standard error can be derived under standard regularity conditions. First, the density at the MEID is estimated by kernel density estimation, ${\widehat{f}}_h\left(\widehat{MEID}\right)$ , using a Gaussian kernel (Equation (10)). The bandwidth h is selected according to Silverman’s rule of thumb (implemented as bw = “nrd0” in R, the default option) [13]. Specifically, the bandwidth equals 0.9 times the minimum of the sample standard deviation and the ratio—provided it is strictly positive—of the sample interquartile range to the interquartile range of the standard normal distribution (1.349), multiplied by $n^{-1/5}$ , as shown in Equation (10).

Second, the influence function of the MEID is specified [14], as presented in Equation (11), and its variance is computed. This allows application of the standard asymptotic variance formula for quantile estimators [15], yielding the expression given in Equation (12). The resulting asymptotic standard error can then be used to construct a Wald-type confidence interval, as illustrated in Equation (9).

The influence function of a quantile is proportional to the difference between the quantile level and the indicator function of the observation being less than or equal to the quantile, scaled by the reciprocal of the density evaluated at the quantile. It quantifies the effect of an infinitesimal contamination at a given point on the estimator and corresponds to the first-order directional derivative of the associated statistical functional [14].

$\begin{array}{l}x=\left\{x_1,x_2,\dots ,x_n\right\}\\ \widehat{MEID}=mdn\left({\left\{\left|x_{\left(i\right)}-x_{\left(i^{\prime }\right)}\right|\right\}}_{\begin{array}{cc}i=1& i^{\prime }>i\end{array}}^{\begin{array}{cc}n-1& n\end{array}}\right)\\ {\widehat{f}}_h\left(x=\widehat{MEID}\right)=\frac{1}{nh\sqrt{2\pi }}{\displaystyle \sum }_{i=1}^n{e}^{-\frac{1}{2}{\left(\frac{x-x_i}{h}\right)}^2}\\ h=0.9\times \min \left(s_{n-1}\left(x\right),\frac{q_x\left(0.75\right)-q_x\left(0.25\right)}{{\Phi }^{-1}\left(0.75\right)-{\Phi }^{-1}\left(0.25\right)}>0\right)\times n^{-\frac{1}{5}}\\ =\frac{0.9}{\sqrt[5]{n}}\times \min \left(s_{n-1}\left(x\right),\frac{R_{IQ}\left(x\right)}{1.349}>0\right)\end{array}$ (10)

$\widehat{\psi }={\left\{{\widehat{\psi }}_i\right\}}_{i=1}^n={\left\{\frac{1}{n-1}{\displaystyle \sum }_{i\ne i^{\prime }}\mathbb{I}\left(\left|x_i-x_{i^{\prime }}\right|\le \widehat{MEID}\right)-0.5\right\}}_{i=1}^n$ (11)

$ase\left(\widehat{MEID}\right)=\sqrt{\frac{s^2\left(\widehat{\psi }\right)}{n{\left({\widehat{f}}_h\left(\widehat{MEID}\right)\right)}^2}}$ (12)

It should be noted, however, that the validity of this asymptotic approximation relies on regularity conditions, including the existence of a positive and continuous density at the target quantile. For discrete distributions, or in situations with substantial mass points, these conditions may not hold. In such cases, the density at the quantile may be ill-defined or unstable, potentially leading to inflated or undefined asymptotic standard errors and poor coverage of Wald-type confidence intervals [16]. Therefore, caution is warranted when applying this approach to discrete or highly discretized data, and resampling-based alternatives may provide more reliable inference in those settings [17].

2.2.3. Interindividual Variance (IVAR)

Because IVAR corresponds to twice the sample variance, its asymptotic standard error is equal to twice the asymptotic standard error of the sample variance, as shown in Equation (13), where n denotes the sample size, m₄ the fourth central moment, and s² the sample variance. This error term can then be used to construct a Wald-type confidence interval, as presented in Equation (9).

$\begin{array}{c}ase\left(\widehat{IVAR}\right)=2\times ase\left(s_{n-1}^2\left(x\right)\right)=2\times \sqrt{\frac{m_4\left(x\right)-\frac{n-3}{n-1}{s}_{n-1}^2\left(x\right)}{n}}\\ =2\times \sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{\left(x_i-\overline{x}\right)}^4/n-\left(n-3\right)/\left(n-1\right){{\displaystyle \sum }}_{i=1}^n{\left(x_i-\overline{x}\right)}^2/\left(n-1\right)}{n}}\end{array}$ (13)

2.2.4. Interindividual Standard Deviation (ISD)

Similarly, because ISD corresponds to the sample standard deviation multiplied by the square root of two, its asymptotic standard error is equal to the asymptotic standard error of the sample standard deviation multiplied by the square root of two, as shown in Equation (14), where n denotes the sample size, s² the sample variance and m₄ the fourth central moment. This error term can then be used to construct a Waldtype confidence interval, as illustrated in Equation (9).

$\begin{array}{c}ase\left(\widehat{ISD}\right)=\sqrt{2}\times ase\left(s_{n-1}\left(x\right)\right)=\sqrt{2}\times \frac{ase\left(s_{n-1}^2\left(x\right)\right)}{4{\widehat{\sigma }}^2\left(x\right)}\\ =\sqrt{2}\times \sqrt{\frac{m_4\left(x\right)-\frac{n-3}{n-1}{s}_{n-1}^4\left(x\right)}{4n{s}_{n-1}^2\left(x\right)}}=\sqrt{\frac{m_4\left(x\right)-s_{n-1}^4\left(x\right)}{2n{s}_{n-1}^2\left(x\right)}}\\ =\frac{ase\left(\widehat{IVAR}\right)}{2\times \widehat{ISD}}\end{array}$ (14)

2.2.5. Interindividual Coefficient of Variation (ICV)

The dependence of the ICV on the sample range—an extreme-value functional and therefore a non-regular quantity—precludes the derivation of a closed-form asymptotic standard error via the delta method [18]. Although the ICV is a smooth function of the sample variance and can be expressed as a normalized scale statistic obtained by dividing the interindividual standard deviation by the sample range. A simplification may be possible when the sample is drawn from a bounded distribution with known support [a, b]. Examples of such distributions include the beta, Program Evaluation and Review Technique (PERT), uniform, triangular, Bernoulli, and binomial distributions in this case, the population range (b − a) can be treated as constant.

Under this assumption, a formal delta-method expression can be derived by scaling the asymptotic standard error of the ISD by 1/(b − a) [18] [19]. For distributions with unit support, such as the standard uniform U [0, 1], beta or Bernoulli, this scaling leaves the asymptotic standard error of the ISD unchanged. However, this approximation is unreliable in practice and yields confidence intervals with very poor coverage when the population range is large.

This issue is illustrated by a simple example. When 100 data points are sampled from a PERT (0, 8, 10) distribution using a fixed seed (123), the asymptotic standard error of the ICV is 0.0147 (ICV = 31.8197; 95% Wald-type Confidence Interval (CI): 31.7909, 31.8485; width (w) = 0.0576), whereas the bootstrap standard error is 2.1651 (95% Bias-Corrected and accelerated (BCa) CI: 27.0781, 35.2336; w = 8.1555).

This limitation is intrinsic to the definition of the ICV: dividing by a fixed constant does not stabilize the sampling variability of the statistic, even when the range is known. In particular, the resulting standard error depends on the absolute magnitude of the population range rather than on the relative scale of the variable. Therefore, inference for the ICV should rely on resampling procedures, such as the bootstrap, regardless of whether the population support is known [20].

3. Objectives

It is worth noting that variation statistics based on interindividual differences are seldom used as descriptive measures of variability [4]. Their limited use is mainly due to the fact that they are not widely known, are not included in standard statistical software, and are labor-intensive to compute. Nonetheless, these statistics can be calculated using Excel [21] or programmed in the R language [22]. As indicators of intraindividual variability [23], they have been applied in fields such as personality research [24] and population genetics [25].

The present methodological study does not conceptualize these statistics as field-specific indicators tied to particular areas of applied research. Rather, it considers them as general descriptive statistical measures applicable to any quantitative variable, whether continuous (e.g., reaction time, age, height, weight) or discrete (e.g., frequency of a behavior, number of correct answers). This perspective underscores the importance of having practical tools for their computation and a more detailed understanding of their statistical properties. Moreover, MEID and ICV are proposed as two new statistics of this type.

Addressing this gap in descriptive and inferential statistics, the objective of this study is to develop an R script that computes point and interval estimates for five measures of interindividual variation: MID, MEID, IVAR, ISD, and ICV. For interval estimation, both asymptotic and bootstrap approaches are implemented. For MEID, the jackknife method is included to estimate bias, standard error, and confidence intervals. Because MEID is based on the median of interindividual differences, the generation of bootstrap distributions can be problematic when applied to data drawn from a discrete distribution, such as the binomial distribution. Consequently, an alternative to the bootstrap method is considered. In addition, the bootstrap-based sampling distributions are examined in terms of asymmetry, mesokurtosis, and normality.

The R script is illustrated using two random samples. One sample consists of 1000 observations drawn from a left-skewed PERT distribution, which may simulate scores on a knowledge-ability test. The other sample consists of 100 observations drawn from a right-skewed binomial distribution, which may simulate the recording of a behavior. The R environment was selected because it is a comprehensive statistical platform widely endorsed by the statistical and mathematical communities and is freely available [26].

To provide greater assurance in the interpretation of the presented examples, two well-known distributions with finite moments that are widely used in simulation studies were selected [27]. One distribution is continuous and based on a large sample size, allowing the examination of convergence between asymptotic and bootstrap standard errors and confidence intervals. The other distribution is discrete and based on a medium-sized sample, serving to demonstrate the difficulty of generating the sampling distribution for the MEID [28] [29] and to observe how the distributional behavior of other statistics manifests, with particular attention to the relative measure. Both distributions were non-normal.

A simulation study, presented in a second article, evaluates whether the bootstrap and asymptotic approaches exhibit equivalent performance by comparing their standard errors and confidence interval widths and by testing whether the sampling distributions of interindividual variation statistics converge to normality.

4. Method

Confidence intervals for five variation statistics based on interindividual comparisons (MID, MEID, IVAR, ISD, and ICV) were computed using a bootstrap approach, specifically the biascorrected and accelerated (BCa) percentile method. This method draws 2000 bootstrap samples of n data points with replacement from the original random sample of size n [30]. In addition, the asymptotic approach is included for the first four statistics because it is applicable, which is not the case for ICV.

The bootstrap approach was selected because it does not rely on parametric assumptions about the population distribution and, under regularity conditions, provides consistent estimates of standard errors and confidence intervals by approximating the sampling distribution of the statistic through resampling from the empirical distribution [31]. In addition, as sample size increases, the bootstrap distribution converges to the true sampling distribution for a broad class of smooth statistics. The BCa method was employed because it adjusts for both median bias and skewness in the sampling distribution via the bias-correction and acceleration parameters, thereby improving interval accuracy under asymmetric conditions [32]. The use of 2000 bootstrap resamples follows the recommendations of Efron and Narasimhan [33] and represents a practical compromise between computational cost and the stability of confidence interval estimation.

For the MEID statistic, bootstrap percentile (perc) method was included, which computes confidence interval more easily than BCa method. The jackknife method was also incorporated to estimate bias, standard error, and TypeT confidence intervals as an alternative approach in cases where bootstrap estimates are difficult to obtain. This issue may arise with large samples drawn from discrete distributions (e.g., binomial, hypergeometric, or Poisson), because MEID is based on the median of interindividual differences. As a result, the variability of the estimates can be greatly reduced, in some cases even yielding a constant value (degenerated distribution). The jackknife is a more stable and straightforward resampling method; however, it does not guarantee greater variability in median-based estimates [34].

The operation of the R script was illustrated using two random samples, which were displayed graphically by means of box plots. One sample consisted of 100 observations drawn from a binomial distribution (n = 10, p = 0.25), which could represent 10 independent records of bullying behavior (present or absent) among adolescents in compulsory secondary education. The second sample consisted of 1000 observations drawn from a PERT (0, 8, 10) distribution, which could represent scores on a knowledge-ability test ranging from 0 to 10 collected from university students. A single seed (123) was set to generate the random sample and ensure reproducibility.

For the five variation measures based on interindividual differences, point and bootstrap estimates, bootstrap bias, bootstrap standard errors, and the lower limits, upper limits, and widths of the 95% bootstrap confidence intervals were calculated. The symmetry of the corresponding bootstrap distribution was examined using the D’Agostino test [35]; mesokurtosis was evaluated using the Anscombe-Glynn test [36]; and goodness-of-fit to the normal distribution was assessed with the Shapiro-Francia W test [37] [38], the K² test [39], and a normal Q-Q plot. The bootstrap distribution was graphically represented using a histogram. The asymptotic standard error and the 95% Wald-type confidence interval were provided for the MID, MEID, IVAR, and ISD statistics.

5. Results

The R script is provided in the Appendix and is also available as a Word file in the following GitHub repository: https://github.com/josemoraldelarubia579/R_script_for_IV_Measures.git.

By default, the script renders the graphics in the console. To save the graphics as JPEG files, the hash symbols (#) must be removed from the commands: #jpeg (…), #par (…), and #dev.off (). The script can be executed online (with graphics displayed in the console) at https://rdrr.io/snippets/. When run with a sample of 1000 observations drawn from a PERT distribution (Example 1), the user execution time is 55.130 s, the system time is 0.340 s, and the elapsed time is 55.789 s. When the same script and sample are executed on a desktop computer (with graphics saved as JPEG files) equipped with an Intel® Celeron® J4025 CPU @ 2.00 GHz, 12.0 GB of RAM (11.8 GB usable), and a 64-bit operating system (x64-based processor), the user execution time is 65.92 s, the system time is 5.91 s, and the elapsed time is 91.91 s, based on a single execution.

The script can be run directly in R (version 3.0 or later; preferably the most recent version, 4.5.3) or using the RStudio environment (R version 3.6 or later). Three packages are required: boot [40], moments [41], and nortest [42]. R can be downloaded from https://cran.r-project.org/bin/windows/base/, and RStudio from https://posit.co/download/rstudio-server/. The mc2d package [43] was used to generate the sample drawn from the PERT distribution.

To measure the execution time of the script applied to a sample of 100 observations drawn from a binomial distribution (Example 2), the BCa confidence interval codes for MEID were disabled: #boot_ci_MEID_bca < -…, #boot_width_MEID_bca < -..., and #cat ("95% BCa confidence interval…). Otherwise, an error message was generated because the BCa method could not be applied due to a degenerate bootstrap distribution. Under these conditions, online execution yields a user time of 1.960 s, a system time of 0.140 s, and an elapsed time of 2.173 s. When executed locally on the desktop computer, the user execution time is 1.22 s, the system time is 0.58 s, and the elapsed time is 1.84 s, based on a single execution.

5.1. Applying the Script to a Random Sample Drawn from a Continuous Distribution

A sample is composed of 1000 observations and was generated from a PERT distribution with minimum 0, mode 8, and maximum 10. The corresponding threshold parameters are 0 and 10, and the shape parameters are α = 4.2 and β = 1.8, as shown in Equation (15). This bounded, non-normal distribution is defined on the interval [0, 10] and exhibits negative or left-tailed skewness (β₁ = −0.5774 < 0) and slight platykurtosis (β₂ = 2.778 < 3). Its mathematical expectation is 7, median is 7.25, mode is 8, and its standard deviation is $\sqrt{3}=1.732$ , which is approximately one-sixth of the range: 10/6 = 1.667 [44].

$\begin{array}{l}X~PERT\left(a=0,b=8,c=10\right)\\ \alpha =1+4\frac{b-a}{c-a}=1+4\frac{8-0}{10-0}=4.2\\ \beta =1+4\frac{c-b}{c-a}=1+4\frac{10-8}{10-0}=1.8\end{array}$ (15)

For illustrative purposes, the values of this first sample may plausibly represent the grades of psychology students in a statistics course. The sample distribution is represented using a box plot, as shown in Figure 1, which is generated by the R script.

Sample size: n = 1000.

Sample minimum: min = 1.4106.

Sample maximum: max = 9.9199.

Figure 1. Box plot for sample distribution.

After executing the R script, the following results are obtained. The asymptotic and bootstrap approaches for estimating the standard errors of MID, MEID, IVAR, and ISD yield remarkably similar results, indicating that asymptotic approximations perform well even though the underlying population distribution is bounded and skewed. Likewise, the widths of the Wald-type and BCa confidence intervals are very similar. This suggests that sampling variability is captured comparably by both methods and that corrections (of bias and skewness) incorporated in the BCa procedure do not meaningfully inflate uncertainty for these statistics. Asymptotic standard errors were only marginally wider than bootstrap standard errors in most cases, as well as the width of confidence intervals, except to MEID.

Even the jackknife standard error of MEID, which is added in the R script to be computed when the bootstrap procedure is unable to generate the 1000 resamples for a discrete distribution, is very similar to the asymptotic and bootstrap standard errors, being slightly smaller than both (jackknife SE = 0.0342 < asymptotic SE = 0.0373 < bootstrap SE = 0.0384). Consequently, the jackknife confidence interval for MEID is slightly narrower (w = 0.134) than the asymptotic and bootstrap confidence intervals (w = 0.1462 and w = 0.1460, respectively).

The value of ICV was 28.8329, for which only bootstrap estimation was available (95% BCa confidence interval: 27.2483 - 30.0138), yielded within the interval [20, 40), showing a medium variability for this dataset. The standard error associated with ICV and width of BCa confidence interval were comparatively larger than that of the absolute measures due to the leptokurtosis of its distribution.

The histograms show good agreement with the overlaid normal density, and the Normal Q-Q plots display points closely aligned with the reference line (Figures 2-11). Specifically, the histograms are presented in the even-numbered figures (2, 4, 6, 8, and 10), whereas the corresponding Q-Q plots appear in the odd-numbered figures (3, 5, 7, 9, and 11). The results of normality K² and W’ tests, consistently fail to reject the null hypothesis of normality for MID, MEID, IVAR, and ISD. This indicates that their empirical bootstrap distributions are approximately Gaussian despite the non-normality of the original data, supporting the adequacy of Wald-type confidence intervals in these cases. In contrast, ICV departs from normality: its bootstrap distribution exhibits significant positive skewness and slight leptokurtosis. This right-skewed asymmetry accounts for the larger bootstrap bias and the wider BCa interval observed for ICV, confirming that bootstrap-based inference is essential for this measure given its non-linear dependence on the sample range. The Q-Q plot reveals this departure from normality through a convex curvature at the upper tail, whereas the histogram with superimposed density curves does not, indicating that the deviation from normality is present but not particularly pronounced.

Mean Interindividual Difference: MID = mean (|x_i − x_j|):

MID = 1.9587.

Asymptotic standard error: 0.0418.

95% Wald-type confidence interval: (1.8768, 2.0407); width = 0.1639.

Bootstrap standard error: 0.0432.

Bootstrap bias: −0.0042.

95% BCa confidence interval (1.878, 2.0434); width = 0.1655.

D’Agostino skewness test for the bootstrap distribution of MID:

$\sqrt{b_1}=-0.0379,z=-0.4922,\text{ p-value}=\mathrm{0.6226.}$

Anscombe-Glynn kurtosis test for the bootstrap distribution of MID:

b₂ = 2.7786, z = −1.5084, p-value = 0.1315.

Shapiro-Francia normality test for the bootstrap distribution of MID:

W’ = 0.9989, z = −0.7169, p-value = 0.7633.

D’Agostino-Belanger-D’Agostino normality test for the bootstrap distribution of MID: K² = 2.5175, p-value = 0.2840.

Figure 2. Histogram of the bootstrap distribution for MID.

Figure 3. Normal Q-Q plot of the bootstrap distribution for MID.

Median Interindividual Difference: MEID = median (|x_i − x_j|):

MEID = 1.6528.

Asymptotic standard error: 0.0373.

95% Wald-type confidence interval: (1.5797, 1.7259); width = 0.1462.

Bootstrap standard error: 0.0384.

Bootstrap bias: −0.003.

95% BCa confidence interval (1.5855, 1.7315); width = 0.1460.

Jackknife standard error: 0.0342.

Jackknife bias: 0.0323.

95% jackknife T-type confidence interval: (1.5535, 1.6875); width = 0.1340.

Figure 4. Histogram of the bootstrap distribution for MEID.

Figure 5. Normal Q-Q plot of the bootstrap distribution for MEID.

D’Agostino skewness test for the bootstrap distribution of MEID:

$\sqrt{b_1}=0.0053,z=0.0696,\text{ p-value}=\mathrm{0.9445.}$

Anscombe-Glynn kurtosis test for the bootstrap distribution of MEID:

b₂ = 2.8078, z = −1.2711, p-value = 0.2037.

Shapiro-Francia normality test for the bootstrap distribution of MEID:

W’ = 0.9986, z = −0.121, p-value = 0.5482.

D’Agostino-Belanger-D’Agostino normality test for the bootstrap distribution of MEID:

K² = 1.6205, p-value = 0.4448.

Interindividual Variance: IVAR = mean ((x_i − x_j)²) = 2 × var (x):

IVAR = 6.0196.

Asymptotic standard error: 0.2559.

95% Wald-type confidence interval: (5.5181, 6.5211); width = 1.003.

Bootstrap standard error: 0.2633.

Bootstrap bias: −0.0197.

95% BCa confidence interval (5.5425, 6.5708); width = 1.0283.

D’Agostino skewness test for the bootstrap distribution of IVAR:

$\sqrt{b_1}=-0.002,z=-0.0254,\text{ p-value}=\mathrm{0.9797.}$

Anscombe-Glynn kurtosis test for the bootstrap distribution of IVAR:

b₂ = 2.7174, z = −2.0323, p-value = 0.0421.

Shapiro-Francia normality test for the bootstrap distribution of IVAR:

W’ = 0.9983, z = 0.2838, p-value = 0.3883.

D’Agostino-Belanger-D’Agostino normality test for the bootstrap distribution of IVAR:

K² = 4.131, p-value = 0.1268.

Figure 6. Histogram of the bootstrap distribution for IVAR.

Figure 7. Normal Q-Q plot of the bootstrap distribution for IVAR.

Interindividual Standard Deviation: $ISD=\sqrt{IVAR}=\sqrt{\left[\text{mean}\left({\left(x_i-x_j\right)}^2\right)\right]}=\sqrt{\left[2\times sd\left(x\right)\right]}$ :

ISD = 2.4535.

Asymptotic standard error: 0.0521.

95% Wald-type confidence interval: (2.3513, 2.5557); width = 0.2044.

Bootstrap standard error: 0.0538.

Boot strap bias: −0.0046.

95% BCa confidence interval (2.3542, 2.5628); width = 0.2086.

D’Agostino skewness test for the bootstrap distribution of ISD:

$\sqrt{b_1}=-0.0587,z=-0.7632,\text{ p-value}=\mathrm{0.4453.}$

Anscombe-Glynn kurtosis test for the Bootstrap distribution of ISD:

b₂ = 2.7297, z = −1.9242, p-value = 0.0543.

Shapiro-Francia normality test for the Bootstrap distribution of ISD:

W’ = 0.9981, z = 0.5151, p-value = 0.3032.

D’Agostino-Belanger-D’Agostino normality test for the Bootstrap distribution of ISD:

K² = 4.2848, p-value = 0.1174.

Figure 8. Histogram of the bootstrap distribution for ISD.

Figure 9. Normal Q-Q plot of the bootstrap distribution for ISD.

Interindividual Coefficient of Variation: ICV = 100 × ISD/(max (x) – min (x)).

ICV = 28.8329.

Bootstrap standard error: 0.7582.

Bootstrap bias: 0.2898.

95% BCa confidence interval (27.2483, 30.0138); width = 2.7655.

D’Agostino skewness test for the Bootstrap distribution of ICV:

$\sqrt{b_1}=0.3625,z=4.5744,\text{ p-value}=0.$

Anscombe-Glynn kurtosis test for the Bootstrap distribution of ICV:

b₂ = 3.5441, z = 2.9395, p-value = 0.0033.

Shapiro-Francia normality test for the Bootstrap distribution of ICV:

W’ = 0.9916, z = 3.9734, p-value = 0.

D’Agostino-Belanger-D’Agostino normality test for the Bootstrap distribution of ICV:

K² = 29.5655, p-value = 0.

Figure 10. Histogram of the bootstrap distribution for ISD.

Figure 11. Normal Q-Q plot of the bootstrap distribution for ISD.

5.2. Applying the Script to a Random Sample Drawn from a Discrete Distribution

The second random sample is composed of 100 data points and was drawn from a binomial distribution with parameters n (number of independent trials) equal to 10 and p (probability of success) equal to one quarter. This bounded discrete distribution is positively skewed ($\sqrt{{\beta }_1\left(X\right)}=\left(q-p\right)/\sqrt{\left(n\times p\times q\right)}=\left(0.75-0.25\right)/\sqrt{\left(10\times 0.25\times 0.75\right)}=0.3651>0$ , where q = 1 − p), and is practically mesokurtic (β₂ (X) = 3 + (1 − (6 × p × q))/(n × p × q) = 3 + (1 − (6 × 0.25 × 0.75))/(10 × 0.25 × 0.75) = 2.9333 ≈ 3). Its standard deviation $\left(\sigma \left(X\right)=\sqrt{\left(n\times p\times q\right)}=\sqrt{\left(10\times 0.25\times 0.75\right)}=1.3693\right)$ represents approximately one seventh of its range (R (X)/7 = (10 − 0)/7 = 1.4286).

The binomial distribution (n = 10, p = 0.25) exhibits lower skewness than the PERT (0, 8, 10) distribution. Skewness is positive, or toward the right tail, rather than negative, or toward the left tail (β₁ (X) = 0.3651 versus β₁ (X) = −0.5774). The binomial distribution also exhibits lower platykurtosis (β₂ (X) = 2.9333 versus β₂ (X) = 2.7778) and lower variability (σ (X) = 1.3693 versus $\sqrt{3}=1.7321$ ), even though both distributions share the same bounded range (from 0 to 10) and the PERT distribution has lighter tails.

For illustrative purposes, the values of this second sample could plausibly represent records from ten randomly selected days of bullying behaviors during recess among one hundred male adolescents aged 12 to 16 years who are enrolled in compulsory secondary education, have been diagnosed with Conduct Disorder [45], and are receiving cognitive behavioral therapy [46] from school psychologists. The adolescents and their parents provided informed consent for video recording for research purposes, and the data were anonymized to ensure confidentiality.

In Figure 12, an outlier can be observed in the right tail according to Tukey’s criterion [47]—values lying one and a half times the interquartile range above the third quartile—and it can also be seen that the right box and whisker are longer than the left box and whisker, revealing the positive skewness, or right tail, exhibited by the binomial (n = 10, p = 0.25) distribution.

Figure 12. Box plot of a random sample of 100 data points drawn from a binomial (n = 10, p = 0.25) distribution.

In the random sample of 100 observations, the asymptotic and bootstrap standard errors and interval widths for MID, IVAR, and ISD were similar, with the bootstrap values being slightly larger than the asymptotic ones, whereas for the PERT (0, 8, 10) distribution the bootstrap estimates were slightly larger.

For MEID, the bootstrap confidence interval could not be computed using the BCa method, whereas it could be computed for the PERT distribution. This limitation arose because the sampling distribution of MEID was dichotomous, nearly degenerate at a constant, with most values equal to 1 and a few zeros, which resulted in an infinite upper limit. The jackknife standard error and the width of its confidence interval were zero. The bootstrap percentile method also yielded both lower and upper bounds equal to 1, producing a confidence interval of zero width. The bootstrap bias was small (0.011); however, the bootstrap standard error (0.1044) was approximately ten times larger than the asymptotic one (0.0282), rendering it unreliable. Consequently, the software does not compute the Gaussian confidence interval when requested. As a result, the point estimate of MEID is equal to 1 with no associated error, and no confidence interval is required. Indeed, the asymptotic standard error was very small (0.0282).

The relative measure indicated medium variability, slightly greater than that of PERT (0, 8, 10), with a difference of 2.86. This finding is consistent with the fact that the distribution generating the sample is less platykurtic (β₂ (X) =2.9333 for X~BN (n = 10, p = 0.25) versus 2.7778 for X~PERT (0, 9, 10)) and with the presence of an outlier in the sample. The ICV value was 31.6901 (95% BCa confidence interval: 27.4722, 36.0655), with a bootstrap bias of 0.855 and a bootstrap standard error of 2.3437. The error associated with ICV was comparatively larger than that of the absolute measures due to the leptokurtosis of its distribution. It was also larger than that observed for the PERT distribution because of the discrete nature of the binomial distribution.

The sampling distributions of MID exhibited symmetry according to D’Agostino’s test and mesokurtosis according to the Anscombe-Glynn test, conforming to normality under both the Shapiro-Francia W’ test (1972) and D’Agostino et al.’s K² test, as also illustrated in Figures 13-14.

Figure 13. Histogram of the bootstrap distribution for MID.

Figure 14. Normal Q-Q plot of the bootstrap distribution for MID.

The distribution of MEID was dichotomous, strongly positively skewed, and leptokurtic, as most values were equal to 1 and a few equal to 2; consequently, it departed completely from the normal distribution (Figures 15-16).

Figure 15. Histogram of the bootstrap distribution for MEID.

Figure 16. Normal Q-Q plot of the bootstrap distribution for MEID.

The IVAR distribution showed mesokurtosis but slight positive skewness; the null hypothesis of normality was retained under the K² test, though not under the W’ test. The plots (Figures 17-18) reveal a slight departure from normality due to a modest elongation of the right tail.

Figure 17. Histogram of the bootstrap distribution for IVAR.

Figure 18. Normal Q-Q plot of the bootstrap distribution for IVAR.

The sampling distribution of ISD was very similar to that of MID and showed a close fit to a normal distribution (Figures 19-20).

Figure 19. Histogram of the bootstrap distribution for ISD.

Figure 20. Normal Q-Q plot of the bootstrap distribution for ISD.

The ICV distribution also failed to meet normality assumptions due to positive skewness and leptokurtosis; however, being continuous and unimodal, it was not, as shown in Figures 21-22, as far from normality as the MEID distribution (Figures 15-16).

Figure 21. Histogram of the bootstrap distribution for ICV.

Figure 22. Normal Q-Q plot of the bootstrap distribution for ICV.

Mean Interindividual Difference: MID = mean (|x_i − x_j|):

MID = 1.4986.

Asymptotic standard error: 0.1031.

95% Wald-type confidence interval: (1.2965, 1.7007); width = 0.4043.

Bootstrap standard error: 0.1.

Bootstrap bias: −0.002.

95% BCa confidence interval (1.3224, 1.7104); width = 0.3880.

D’Agostino skewness test for the bootstrap distribution of MID:

$\sqrt{b_1}=-0.0273,z=-0.355,\text{ p-value}=\mathrm{0.7226.}$

Anscombe-Glynn kurtosis test for the bootstrap distribution of MID:

b₂ = 3.0965, z = 0.7151, p-value = 0.4745.

Shapiro-Francia normality test for the bootstrap distribution of MID:

W’ = 0.9985, z = 0.0163, p-value = 0.4935.

D’Agostino-Belanger-D’Agostino normality test for the bootstrap distribution of MID:

K² = 0.6374, p-value = 0.7271.

Median Interindividual Difference: MEID = median (|x_i − x_j|):

MEID = 1.

Asymptotic standard error: 0.0282.

95% Wald-type confidence interval: (0.9447, 1.0553); width = 0.1106.

Bootstrap standard error: 0.1044.

Bootstrap bias: 0.011.

95% BCa confidence interval: boot_ci_MEID < - boot.ci (bd, type = “bca”, index = 2). Error in bca.ci (boot.out, conf, index [1L], L = L, t = t.o, t0 = t0.o: estimated adjustment “w” is infinite. Error: object “boot_ci_MEID” did not find

95% bootstrap percentile confidence interval (1, 1); width = 0.

Jackknife standard error: 0.

Jackknife bias: 0.

95% jackknife T-type confidence interval: (1, 1); width = 0.

D’Agostino skewness test for the bootstrap distribution of MEID:

$\sqrt{b_1}=9.3766,z=33.5541,\text{ p-value}=0.$

Anscombe-Glynn kurtosis test for the bootstrap distribution of MEID:

b₂ = 88.9202, z = 20.151, p-value = 0.

Shapiro-Francia normality test for the bootstrap distribution of MEID:

W’ = 0.0761, z = 14.7359, p-value = 0.

D’Agostino-Belanger-D’Agostino normality test for the bootstrap distribution of MEID:

K² = 1531.94, p-value = 0.

Interindividual Variance: IVAR = mean ((x_i − x_j)²) = 2 * var (x):

IVAR = 3.6154.

Asymptotic standard error: 0.4615.

95% Wald-type confidence interval: (2.7108, 4.5199); width = 1.8091.

Bootstrap standard error: 0.4591.

Bootstrap bias: 0.0201.

95% BCa confidence interval (2.8369, 4.619); width = 1.7821.

D’Agostino skewness test for the bootstrap distribution of IVAR:

$\sqrt{b_1}=0.179,z=2.3103,\text{ p-value}=\mathrm{0.0209.}$

Anscombe-Glynn kurtosis test for the bootstrap distribution of IVAR:

b₂ = 3.1034, z = 0.7563, p-value = 0.4495.

Shapiro-Francia normality test for the bootstrap distribution of IVAR:

W’ = 0.9967, z = 1.8186, p-value = 0.0345.

D’Agostino-Belanger-D’Agostino normality test for the bootstrap distribution of IVAR:

K² = 5.9095, p-value = 0.0521.

Interindividual Standard Deviation: $ISD=\sqrt{IVAR}=\sqrt{\left[\text{mean}\left({\left(x_i-x_j\right)}^2\right)\right]}=\sqrt{\left[2\times sd\left(x\right)\right]}$ :

ISD = 1.9014.

Asymptotic standard error: 0.1214.

95% Wald-type confidence interval: (1.6635, 2.1393); width = 0.4757.

Bootstrap standard error: 0.1206.

Boot strap bias: 0.0015.

95% BCa confidence interval (1.6841, 2.1492); width = 0.4651.

D’Agostino skewness test for the bootstrap distribution of ISD:

$\sqrt{b_1}=-0.0171,z=-0.2219,\text{ p-value}=\mathrm{0.8244.}$

Anscombe-Glynn kurtosis test for the Bootstrap distribution of ISD:

b₂ = 3.0707, z = 0.56, p-value = 0.5755.

Shapiro-Francia normality test for the Bootstrap distribution of ISD:

W’ = 0.9985, z = 0.0744, p-value = 0.4704.

D’Agostino-Belanger-D’Agostino normality test for the Bootstrap distribution of ISD:

K² = 0.3628, p-value = 0.8341.

Interindividual Coefficient of Variation: ICV = 100 × ISD/(max (x) − min (x)):

ICV = 31.6901.

Bootstrap standard error: 2.3437.

Bootstrap bias: 0.855.

95% BCa confidence interval (27.4722, 36.0655); width = 8.5934.

D’Agostino skewness test for the Bootstrap distribution of ICV:

$\sqrt{b_1}=0.5258,z=6.4444,\text{ p-value}=0.$

Anscombe-Glynn kurtosis test for the Bootstrap distribution of ICV:

b₂ = 3.5675, z = 3.0362, p-value = 0.0024.

Shapiro-Francia normality test for the Bootstrap distribution of ICV:

W’ = 0.9838, z = 5.4658, p-value = 0.

D’Agostino-Belanger-D’Agostino normality test for the Bootstrap distribution of ICV:

K² = 50.7486, p-value = 0.

6. Discussion

The four absolute statistics based on interindividual differences are primarily intended for quantitative variables and are most naturally interpreted when the underlying measurement scale is continuous. When the variables are ordinal, particularly with a limited number of ordered categories—such as Likert-type items with five to seven response options—the median interindividual difference can be computed without strong metric assumptions, and the mean interindividual difference is often used as an approximation under an assumed interval structure. In contrast, the interindividual variance and the interindividual standard deviation, which rely on squared differences and thus on stronger assumptions about equal intervals and metric properties, are less straightforward to interpret for such ordinal variables and should therefore be used with caution [48] [49].

The most novel measure proposed is the interindividual coefficient of variation (ICV), defined as the ratio of the interindividual standard deviation to the sample range. As with other coefficients of variation, it is multiplied by 100. This coefficient belongs to the class of measures proposed by Merce [6], with which a direct correspondence can be established. Like other coefficients of variation, its meaningful application presupposes quantitative variables with a well-defined metric structure. Consequently, it is not appropriate for ordinal variables, particularly those with a very limited number of categories—such as Likert-type items—where both the standard deviation and the range are strongly constrained by the response scale.

Unlike classical coefficients of variation, the ICV does not, under general conditions, attain a theoretical upper bound of 100. For bounded variables, the ratio between the standard deviation and the sample range is itself bounded. Extremely bimodal discrete distributions, with all probability mass concentrated at the boundaries—such as a symmetric Bernoulli distribution (p = 0.5)—reach this bound and therefore provide a natural extreme benchmark for the interindividual coefficient of variation: $ICV=100\times ISD/Range=100\times \sqrt{2}\times s_n{}_{-1}/Range=100/\sqrt{2}$ . For bounded distributions, the statistic could be adjusted by multiplying it by $\sqrt{2}$ so that its potential range extends from 0 to 100. This rescaling $\left(ICVr=\surd 2\times 100\times ISD/Range=100\times 2\times s_n{}_{-1}/Range\right)$ does not alter the relative ordering of variability across samples but facilitates interpretability and comparability with other normalized dispersion indices, such as the quartile coefficient of variation.

In contrast, for unbounded distributions no such bound exists at the sample level. In finite samples drawn from distributions with unbounded support, the sample range may become arbitrarily small due to sampling variability, while the sample standard deviation remains positive, allowing the ratio s_{n − 1}/Range to exceed the Bernoulli benchmark. Consequently, it is not possible to define an absolute upper limit, and the rescaling obtained by multiplying by $\sqrt{2}$ lacks universal validity, particularly in samples from unknow distributions.

With large samples drawn from continuous distributions with finite moments, at least up to the fourth order, asymptotic and bootstrap-based approaches for the four absolute measures of interindividual variation are expected to yield similar results, as illustrated by the first example with 1000 observations, where the conditions for asymptotic approximations are well satisfied [50]. For relatively small random samples (approximately 30 ≤ n ≤ 50), bootstrap methods may offer a useful alternative for estimating standard errors and confidence intervals, although their performance depends on the smoothness of the statistic and the underlying distribution. When sample sizes fall below this range, bootstrap estimates may become unstable and their inferential accuracy may deteriorate [51]. Similar considerations apply to samples of any size drawn from distributions lacking finite fourth-order or higher moments, such as Student’s t distributions with fewer than four degrees of freedom, for which both asymptotic and bootstrap-based inference should be interpreted with caution.

With medium and large samples drawn from discrete distributions, the computation of the standard error and confidence interval of the MEID using resampling methods is problematic or impossible [11] [28]. The resulting sampling distribution is degenerate because it converges to a constant with a null standard error [29]. This suggests that the central limit theorem does not apply, as demonstrated by the sample of 100 observations in the second example, and it is therefore preferable to report only the point estimate.

In the discrete-data MEID case, it should be clarified that the nearly degenerate bootstrap distribution results from the limited set of possible pairwise differences and the resampling process. With discrete data containing few distinct values, interindividual differences cluster around a small number of outcomes, restricting the population median to a narrow range. Bootstrap resamples repeatedly reproduce these same configurations, yielding identical or nearly identical MEID estimates and thus a nearly degenerate distribution.

Practically, this means the bootstrap offers too little variability for meaningful interval estimation. A useful diagnostic is the proportion of unique MEID bootstrap values: if only a very small fraction (e.g., < 1 - 5%) of resamples produce distinct estimates, the distribution can be considered degenerate. In such cases, interval estimates should be suppressed or reported as undefined because the bootstrap cannot approximate the statistic’s sampling variability.

The sampling distribution of the relative measure, for which the delta method cannot define the asymptotic standard error [52], can be accessed via the bootstrap procedure [20]. This procedure is valid for estimating both the standard error and the bias, as well as for computing confidence intervals using either the BCa or percentile method [33]. These intervals can be computed for both continuous and discrete distributions, as demonstrated by the two examples presented.

The present article is limited to introducing an R script for the complete computation of interindividual variation measures and illustrating its performance with two samples: a large sample drawn from a continuous distribution and a medium-sized sample drawn from a discrete distribution. A more comprehensive evaluation of these measures requires dedicated simulation studies, such as those reported in three complementary articles accompanying the present proposal.

The derivation of the asymptotic standard errors was based on general statistical principles. However, further methodological work grounded in formal statistical theory would facilitate more rigorous derivations under restrictive assumptions, particularly for the MEID and, more specifically, for the Merce-type relative measure (ICV).

7. Conclusions

It is concluded that the R script is appropriate for its intended purposes and functions correctly with samples drawn from both continuous and discrete distributions. The large sample drawn from a non-normal continuous distribution with finite moments allows verification of the convergence of asymptotic and bootstrap standard errors and confidence intervals. The medium-sized sample drawn from a discrete distribution illustrates the difficulty of generating the sampling distribution for the MEID, which in this case becomes nearly degenerate (i.e., almost constant), a problem that did not arise with the continuous data sample. This discrete-data example also shows approximate convergence to normality for several absolute measures of interindividual variation.

In the two empirical illustrations presented here, the relative measure based on the range (ICV) does not exhibit convergence to normality. This behavior is consistent with the fact that the statistic depends on the sample range, whose sampling distribution is governed by extreme-order statistics and may remain highly skewed, especially in moderate samples or in the presence of discrete data. Nevertheless, approximate normality may be plausible under more regular conditions, such as large samples drawn from smooth continuous distributions with stable extremes. Therefore, the present observation should be interpreted as illustrative of the examples analyzed rather than as a general property of the ICV.

8. Suggestions

It is recommended to apply these measures of variation in the descriptive analysis of quantitative variables using the R script developed for this purpose and to explore their extension to ordinal variables. Likewise, it is advisable to deepen their study by employing simulated samples drawn from known continuous and discrete probability distributions, as well as empirical samples derived from hypothesized or unknown distributions. The type of distribution, the type of bounds (bounded, semi-bounded, and unbounded), the nature of the data (continuous quantitative, discrete quantitative, or ordinal), the degree of skewness, and the level of kurtosis may be considered design factors for examining the convergence between the asymptotic and bootstrap approaches in estimating the standard error and constructing the confidence interval, as well as for assessing the approximation to normality of the sampling distribution of the statistics. Moreover, this type of study would allow for the establishment of more precise interpretive intervals for the relative measure (ICV).

Acknowledgements

The author thanks the reviewers and the editor for their valuable comments and suggestions, which substantially improved the quality of the manuscript.

Appendix: Script to Compute Interindividual Variation Statistics

# This R script is also available in the following GitHub repository: https://github.com/josemoraldelarubia579/R_script_for_IV_Measures.git

# Generate the random sample of the first example presented in the article.

library(mc2d) # Required package to generate a sample from the PERT distribution.

set.seed(123) # A fixed seed is used to ensure reproducibility of the random sample.

x <- rpert(n = 1000, min = 0, mode = 8, max = 10) # PERT (0, 8, 10) distribution.

# Generate the random sample of the second example presented in the article.

# set.seed(123) # A fixed seed is used to ensure reproducibility of the random sample.

# x <- rbinom(n = 100, size = 10, prob = 0.25) # Binomial (10, 0.25) distribution.

# Replacing the example data with your own sample

# Option 1: Define your sample manually. For example: x = {8.097046, 7.485967, ..., 6.630036} ⊂ X = "Test scores".

# x <- c(8.097046, 7.485967, ..., 6.630036)

# Option 2: Load a plain text file (e.g., college_data.csv).

# data <- read.csv("C:/Users/YourName/Documents/college_data.csv")

# x <- data$test_score # Extract the random variable test_score and rename it x.

# Option 3: Load an Excel file (e.g., mydata.xlsx) using the ‘readxl’ library.

# library(readxl)

# data <- read_excel("C:/Users/YourName/Documents/college_data.xlsx")

# x <- data$test_score # Extract the random variable test_score and rename it x.

# Option 4: Load an SPSS file (e.g., college_data) using the ‘haven’ library.

# library(haven)

# data <- read_sav("C:/Users/YourName/Documents/college_data.sav")

# x <- data$test_score # Extract the random variable test_score and rename it x.

# Load the required R packages.

# If you don’t have one, install it from a nearby Comprehensive R Archive Network (CRAN) mirror. Select a mirror from the R toolbar (e.g. Mexico (Mexico City)).

# Remove the hash symbol (#) to this line to enable the code.

library(boot)

# install.packages("boot")

# install.packages("moments")

library(moments)

# install.packages("nortest")

library(nortest)

# Round results to 4 decimal places.

# Change 4 to the desired decimals, such as 2, 3, or 6.

rnd <- function(z, k = 4) round(z, k)

# Represent sample distribution using a box plot

# The graph is displayed in the console by default.

# Remove the next three hash symbols to save it as a JPG file.

# jpeg("Boxplot_sample_distribution.jpg", width = 3, height = 2, units = "in", res = 300)

# par(mar = c(4, 1, 1, 1))

boxplot(x, horizontal = TRUE, notch = FALSE, col = "darkolivegreen2", xlab = "X values", cex.axis = 0.7, cex.lab = 0.7)

# dev.off()

# ------------------------------------------------------------------------------------------------------------

# Point estimates of the variation statistics based on interindividual differences

# ------------------------------------------------------------------------------------------------------------

# Compute the interindividual differences

diffs <- as.vector(dist(x))

# Compute the variation statistics based on interindividual differences

# Mean Interindividual Difference (MID)

MID <- mean(diffs)

# Median Interindividual Difference (MEID)

MEID <- median(diffs)

# Interindividual Variance (IVAR)

IVAR <- mean(diffs^2)

# Interindividual Standard Deviation (ISD)

ISD <- sqrt(IVAR)

# Interindividual Coefficient of Variation (ICV)

ICV = 100 * ISD /(max(x) - min(x))

# ------------------------------------------------------------------------------------------------------------

# Asymptotic Standard Errors and Wald-type 95% confidence intervals

# ------------------------------------------------------------------------------------------------------------

alpha <- 0.05 # Significance level. This value can be changed to 0.1 for small samples.

z <- qnorm(1 - alpha/2) # Critical value

n <- length(x) # Sample size

## MID = mean(diffs)

# Kernel h(x_i, x_j) = |x_i - x_j|, a symmetric second-order U-statistic

# Estimation of the influence function:

# phi_i = E[h(X_i, X_j) | X_i = x_i] - theta ≈ (1 / (n - 1)) * sum_{j != i} |x_i - x_j| - MID

abs_row_means <- vapply(seq_len(n), function(i) {

mean(abs(x[i] - x[-i]))

}, numeric(1))

phi <- abs_row_means - MID

# Asymptotic variance: Var(MID) ≈ 4 * Var(phi) / n

asympt_se_MID <- sqrt(4 * var(phi) / n)

asympt_ci_MID <- c(MID - z * asympt_se_MID, MID + z * asympt_se_MID)

asympt_width_MID = asympt_ci_MID[2] - asympt_ci_MID[1]

## MEID = median(diffs)

den <- density(diffs, bw = "nrd0", kernel = "gaussian")

fhat_MEID <- approx(den$x, den$y, xout = MEID, rule = 2)$y

# Symmetric second-order U-quantile.

# Influence function estimator: psi_i = P(|x_i - X’| <= MEID | X_i = x_i) - 1/2

psi <- vapply(seq_len(n), function(i) {

mean(abs(x[i] - x[-i]) <= MEID) - 0.5

}, numeric(1))

asympt_se_MEID <- sqrt(4 * var(psi) / (n * fhat_MEID^2))

asympt_ci_MEID <- c(MEID - z * asympt_se_MEID, MEID + z * asympt_se_MEID)

asympt_width_MEID = asympt_ci_MEID[2] - asympt_ci_MEID[1]

## IVAR = mean(diffs^2) = 2 * var(x)

# SE(var) ≈ sqrt((m4 - (n - 3)/(n - 1) * s^4) / n)

m4 <- mean((x - mean(x))^4)

se_var <- sqrt(max(m4 - (n - 3)/(n - 1) * var(x)^2, 0) / n)

asympt_se_IVAR <- 2 * se_var

asympt_ci_IVAR <- c(IVAR - z * asympt_se_IVAR, IVAR + z * asympt_se_IVAR)

asympt_width_IVAR = asympt_ci_IVAR[2] - asympt_ci_IVAR[1]

## ISD = √IVAR = √[mean(diffs^2)] = sqrt[2*sd(x)]

asympt_se_ISD <- asympt_se_IVAR / (2 * ISD)

asympt_ci_ISD <- c(ISD - z * asympt_se_ISD, ISD + z * asympt_se_ISD)

asympt_width_ISD = asympt_ci_ISD[2] - asympt_ci_ISD[1]

# ------------------------------------------------------------------------------------------------------------

# Bootstrap standard error and 95% BCa confidence intervals

# ------------------------------------------------------------------------------------------------------------

# Compute all pairwise absolute differences

boot_stats <- function(data, index) {

xb <- data[index]

d <- as.vector(dist(xb))

c(

MID = mean(d),

MEID = median(d),

IVAR = 2 * var(xb),

ISD = sqrt(2) * sd(xb),

ICV = 100 * (sqrt(2) * sd(xb)) / (max(xb) - min(xb))

)}

# Generate the bootstrap distribution of each statistic

set.seed(123) # # A fixed seed is used to ensure reproducibility of the results.

bd <- boot(data = x, statistic = boot_stats, R = 1000)

## ----------------------------------------------------------------------------------------------------------

# Jackknife standard error and jackknife confidence interval (for discrete distributions)

## ----------------------------------------------------------------------------------------------------------

# Leave-one-out MEID

T_loo <- numeric(n)

for (i in seq_len(n)) {

xi <- x[-i]

T_loo[i] <- median(as.vector(dist(xi)))

}

# Pseudovalues and t-type CI

pseudo <- n * MEID - (n - 1) * T_loo

jack_MEID <- mean(pseudo)

jack_se_MEID <- sd(pseudo) / sqrt(n)

t_crit <- qt(1 - alpha/2, df = n - 1) # Hendi (2017)

jack_bias_MEID <- (n-1) * (mean(T_loo) - MEID)

jack_ci_MEID <- c(jack_MEID - t_crit * jack_se_MEID,

jack_MEID + t_crit * jack_se_MEID)

jack_width_MEID <- diff(jack_ci_MEID)

# -------------------------------------------------------------------------------------------------------

# Mean Interindividual Difference (MID)

# -------------------------------------------------------------------------------------------------------

# Tests for symmetry, kurtosis, and normality for the bootstrap distribution of MID

# D’Agostino skewness test

sq_b1_MID <- agostino.test(bd$t[,1])

# Anscombe–Gynn kurtosis test

b2_MID <- anscombe.test(bd$t[,1])

# Shapiro–Francia normality test

w_MID <- sf.test(bd$t[,1])

# D’Agostino–Belanger–D’Agostino normality test

k2_MID <- sq_b1_MID$statistic[2]^2 + b2_MID$statistic[2]^2

p_value_MID <- pchisq(k2_MID, df = 2, lower.tail = FALSE)

# Display or save the normal QQ plot to evaluate normality graphically

# QQ-plot normal. The graph is displayed in the console by default.

# Remove the next three hash symbols to save it as a JPG file.

# jpeg("QQ_plot_MID.jpg", width = 3, height = 2, units = "in", res = 300)

# par(mar = c(4, 4, 1, 1), cex.main = 0.8, cex.axis = 0.75, cex.lab = 0.7)

qqnorm(bd$t[,1], main = "Normal QQ plot for MID", xlab = "Theoretical quantiles: normal", ylab = "Sample quantiles: boot", cex.main = 0.8, cex.axis = 0.7, cex.lab = 0.7)

qqline(bd$t[,1], col = "darkblue")

# dev.off()

# Display or save density histogram with overlaid density curve.

# Remove the next three hash symbols to save the graph as a JPG file.

# jpeg("Histogram_MID.jpg", width = 3, height = 2, units = "in", res = 300)

# par(mar = c(4, 4, 1, 1))

hist(bd$t[,1], breaks = 50, probability = TRUE, col = "lightyellow", main = "", xlab = "MID", cex.axis = 0.8, cex.lab = 0.7)

lines(density(bd$t[,1]), col = "orange", lwd = 2)

x_values <- seq(mean(bd$t[,1]) - 4 * sd(bd$t[,1]), mean(bd$t[,1]) + 4 * sd(bd$t[,1]), length = 1000)

y_values <- dnorm(x_values, mean = mean(bd$t[,1]), sd = sd(bd$t[,1]))

lines(x_values, y_values, col = "darkred", lwd = 2)

legend("topright", title = "Density curves", legend = c("Empiric", "Normal"), fill = c("orange", "darkred"), cex = 0.6, bty = "n")

# dev.off()

# 95% BCa confidence interval for MID

boot_ci_MID <- boot.ci(bd, type = "bca", index = 1)

boot_width_MID <- boot_ci_MID$bca[5] - boot_ci_MID$bca[4]

# Print the results

cat("\nSample size: n =", n, ".\n")

cat("Sample minimum: min =", rnd(min(x)), ".\n")

cat("Sample maximum: max =", rnd(max(x)), ".\n")

# Print the results about the MID on the console.

cat("\nMean Interindividual Difference: MID = mean(|xi - xj|):\n")

cat(" MID =", rnd(MID), ".\n")

cat(" Asymptotic standard error:", rnd(asympt_se_MID), ".\n")

cat(" 95% Wald-type confidence interval: (", rnd(asympt_ci_MID[1]), "," , rnd(asympt_ci_MID[2]), "); width =", rnd(asympt_width_MID), ".\n")

cat(" Bootstrap standard error:", rnd(sd(bd$t[,1])), ".\n")

cat(" Bootstrap bias:", rnd(mean(bd$t[,1]) - MID), ".\n")

cat(" 95% BCa confidence interval (", rnd(boot_ci_MID$bca[4]), ",", rnd(boot_ci_MID$bca[5]), "); width =", rnd(boot_width_MID), ".\n")

cat(" D’Agostino skewness test for the bootstrap distribution of MID:\n")

cat(" √b1 =", rnd(sq_b1_MID$statistic[1]), ", z =", rnd(sq_b1_MID$statistic[2]), ", p-value =", rnd(sq_b1_MID$p.value), ".\n")

cat(" Anscombe-Glynn kurtosis test for the bootstrap distribution of MID:\n")

cat(" b2 =", rnd(b2_MID$statistic[1]), ", z =", rnd(b2_MID$statistic[2]), ", p-value =", rnd(b2_MID$p.value), ".\n")

cat(" Shapiro-Francia normality test for the bootstrap distribution of MID:\n")

cat(" W’ =", rnd(w_MID$statistic[1]), ", z =", rnd(qnorm(w_MID$p.value, mean = 0, sd = 1, lower.tail = FALSE)), ", p-value =", rnd(w_MID$p.value), ".\n")

cat(" D’Agostino-Belanger-D’Agostino (DBD) normality test for the bootstrap distribution of MID:\n")

cat(" K^2 =", rnd(k2_MID), ", p-value =", rnd(p_value_MID), ".\n")

# ------------------------------------------------------------------------------------------------------------

# Median Interindividual Difference (MEID)

# ------------------------------------------------------------------------------------------------------------

# Tests for symmetry, kurtosis, and normality for the bootstrap distribution of MEID

# D’Agostino skewness test

sq_b1_MEID <- agostino.test(bd$t[,2])

# Anscombe–Gynn kurtosis test

b2_MEID <- anscombe.test(bd$t[,2])

# Shapiro–Francia normality test

w_MEID <- sf.test(bd$t[,2])

# D’Agostino–Belanger–D’Agostino normality test

k2_MEID <- sq_b1_MEID$statistic[2]^2 + b2_MEID$statistic[2]^2

p_value_MEID <- pchisq(k2_MEID, df = 2, lower.tail = FALSE)

# Display or save the normal QQ plot to evaluate normality graphically

# Remove the next three hash symbols to save the graph as a JPG file.

# jpeg("QQ_plot_MEID.jpg", width = 3, height = 2, units = "in", res = 300)

# par(mar = c(4, 4, 1, 1))

qqnorm(bd$t[,2], main = "Normal QQ plot for MEID", xlab = "Theoretical quantiles: normal", ylab = "Sample quantiles: boot", cex.main = 0.8, cex.axis = 0.7, cex.lab = 0.7)

qqline(bd$t[,2], col = "darkblue")

# dev.off()

# Display or save density histogram with overlaid density curve.

# Remove the next three hash symbols to save the graph as a JPG file.

# jpeg("Histogram_MEID.jpg", width = 3, height = 2, units = "in", res = 300)

# par(mar = c(4, 4, 1, 1))

hist(bd$t[,2], breaks = 50, probability = TRUE, col = "lightyellow", main = "", xlab = "MEID", cex.axis = 0.8, cex.lab = 0.7)

lines(density(bd$t[,2]), col = "orange", lwd = 2)

x_values <- seq(mean(bd$t[,2]) - 4 * sd(bd$t[,2]), mean(bd$t[,2]) + 4 * sd(bd$t[,2]), length = 1000)

y_values <- dnorm(x_values, mean = mean(bd$t[,2]), sd = sd(bd$t[,2]))

lines(x_values, y_values, col = "darkred", lwd = 2)

legend("topright", title = "Density curves", legend = c("Empiric", "Normal"), fill = c("orange", "darkred"), cex = 0.6, bty = "n")

# dev.off()

# 95% bootstrap and jackknife confidence interval for MEID

# Add the hash symbol (#) to the beginning of the following two lines to disable the code in case of an error caused by a degenerate bootstrap distribution. This error typically occurs when samples are drawn from discrete distributions or ordinal variables.

boot_ci_MEID_bca <- boot.ci(bd, type = "bca", index = 2)

boot_width_MEID_bca <- boot_ci_MEID_bca$bca[5] - boot_ci_MEID_bca$bca[4]

# 95% bootstrap PERC confidence interval for MEID

boot_ci_MEID_perc <- boot.ci(bd, type = "perc", index = 2)

boot_width_MEID_perc <- boot_ci_MEID_perc$perc[5] - boot_ci_MEID_perc$perc[4]

# Print the results about the MEID on the console.

cat("\nMedian Interindividual Difference: MEID = median(|xi - xj|):\n")

cat(" MEID =", rnd(MEID), ".\n")

cat(" Asymptotic standard error:", rnd(asympt_se_MEID), ".\n")

cat(" 95% Wald-type confidence interval: (", rnd(asympt_ci_MEID[1]), "," , rnd(asympt_ci_MEID[2]), "); width =", rnd(asympt_width_MEID), ".\n")

cat(" Bootstrap standard error:", rnd(sd(bd$t[,2])), ".\n")

cat(" Bootstrap bias:", rnd(mean(bd$t[,2]) - MEID), ".\n")

cat(" 95% BCa confidence interval (", rnd(boot_ci_MEID_bca$bca[4]), ",", rnd(boot_ci_MEID_bca$bca[5]), "); width =", rnd(boot_width_MEID_bca), ".\n")

# Prepend the hash symbol (#) to the previous line to disable the code in case of an error.

cat(" 95% bootstrap percentile confidence interval (", rnd(boot_ci_MEID_perc$perc[4]), ",", rnd(boot_ci_MEID_perc$perc[5]), "); width =", rnd(boot_width_MEID_perc), ".\n")

# Prepend the hash symbol (#) to the previous to disable the code in case of an error.

cat(" Jackknife standard error:", rnd(jack_se_MEID), ".\n")

cat(" Jackknife bias:", rnd(jack_bias_MEID), ".\n")

cat(" 95% jackknife T-type confidence interval: (", rnd(jack_ci_MEID[1]), "," , rnd(jack_ci_MEID[2]), "); width =", rnd(jack_width_MEID), ".\n")

cat(" D’Agostino skewness test for the bootstrap distribution of MEID:\n")

cat(" √b1 =", rnd(sq_b1_MEID$statistic[1]), ", z =", rnd(sq_b1_MEID$statistic[2]), ", p-value =", rnd(sq_b1_MEID$p.value), ".\n")

cat(" Anscombe-Glynn kurtosis test for the bootstrap distribution of MEID:\n")

cat(" b2 =", rnd(b2_MEID$statistic[1]), ", z =", rnd(b2_MEID$statistic[2]), ", p-value =", rnd(b2_MEID$p.value), ".\n")

cat(" Shapiro-Francia normality test for the bootstrap distribution of MEID:\n")

cat(" W’ =", rnd(w_MEID$statistic[1]), ", z =", rnd(qnorm(w_MEID$p.value, mean = 0, sd = 1, lower.tail = FALSE)), ", p-value =", rnd(w_MEID$p.value), ".\n")

cat(" DBD normality test for the bootstrap distribution of MEID:\n")

cat(" K^2 =", rnd(k2_MEID), ", p-value =", rnd(p_value_MEID), ".\n")

# ------------------------------------------------------------------------------------------------------------

# Interindividual Variance (IVAR)

# ------------------------------------------------------------------------------------------------------------

# Tests for symmetry, kurtosis, and normality for the bootstrap distribution of IVAR

# D’Agostino skewness test

sq_b1_IVAR <- agostino.test(bd$t[,3])

# Anscombe–Gynn kurtosis test

b2_IVAR <- anscombe.test(bd$t[,3])

# Shapiro–Francia normality test

w_IVAR <- sf.test(bd$t[,3])

# D’Agostino–Belanger–D’Agostino normality test

k2_IVAR <- sq_b1_IVAR$statistic[2]^2 + b2_IVAR$statistic[2]^2

p_value_IVAR <- pchisq(k2_IVAR, df = 2, lower.tail = FALSE)

# Display or save the normal QQ plot to evaluate normality graphically

# Remove the next three hash symbols to save the graph as a JPG file.

# jpeg("QQ_plot_IVAR.jpg", width = 3, height = 2, units = "in", res = 300)

# par(mar = c(4, 4, 1, 1))

qqnorm(bd$t[,3], main = "Normal QQ plot for IVAR", xlab = "Theoretical quantiles: normal", ylab = "Sample quantiles: boot", cex.main = 0.8, cex.axis = 0.7, cex.lab = 0.7)

qqline(bd$t[,3], col = "darkblue")

# dev.off()

# Display or save density histogram with overlaid density curve.

# Remove the next three hash symbols to save the graph as a JPG file.

# jpeg("Histogram_IVAR.jpg", width = 3, height = 2, units = "in", res = 300)

# par(mar = c(4, 4, 1, 1))

hist(bd$t[,3], breaks = 50, probability = TRUE, col = "lightyellow", main = "", xlab = "IVAR", cex.axis = 0.8, cex.lab = 0.7)

lines(density(bd$t[,3]), col = "orange", lwd = 2)

x_values <- seq(mean(bd$t[,3]) - 4 * sd(bd$t[,3]), mean(bd$t[,3]) + 4 * sd(bd$t[,3]), length = 1000)

y_values <- dnorm(x_values, mean = mean(bd$t[,3]), sd = sd(bd$t[,3]))

lines(x_values, y_values, col = "darkred", lwd = 2)

legend("topright", title = "Density curves", legend = c("Empiric", "Normal"), fill = c("orange", "darkred"), cex = 0.6, bty = "n")

# dev.off()

# 95% BCa confidence interval for IVAR

boot_ci_IVAR <- boot.ci(bd, type = "bca", index = 3)

boot_width_IVAR <- boot_ci_IVAR$bca[5] - boot_ci_IVAR$bca[4]

# Print the results about the IVAR on the console.

cat("\nInterindividual Variance: IVAR = mean((xi - xj)^2) = 2 * var(x):\n")

cat(" IVAR =", rnd(IVAR), ".\n")

cat(" Asymptotic standard error:", rnd(asympt_se_IVAR), ".\n")

cat(" 95% Wald-type confidence interval: (", rnd(asympt_ci_IVAR[1]), ", ", rnd(asympt_ci_IVAR[2]), "); width =", rnd(asympt_width_IVAR), ".\n")

cat(" Bootstrap standard error:", rnd(sd(bd$t[,3])), ".\n")

cat(" Bootstrap bias:", rnd(mean(bd$t[,3]) - IVAR), ".\n")

cat(" 95% BCa confidence interval (", rnd(boot_ci_IVAR$bca[4]), ",", rnd(boot_ci_IVAR$bca[5]), "); width =", rnd(boot_width_IVAR), ".\n")

cat(" D’Agostino skewness test for the bootstrap distribution of IVAR:\n")

cat(" √b1 =", rnd(sq_b1_IVAR$statistic[1]), ", z =", rnd(sq_b1_IVAR$statistic[2]), ", p-value =", rnd(sq_b1_IVAR$p.value), ".\n")

cat(" Anscombe-Glynn kurtosis test for the bootstrap distribution of IVAR:\n")

cat(" b2 =", rnd(b2_IVAR$statistic[1]), ", z =", rnd(b2_IVAR$statistic[2]), ", p-value =", rnd(b2_IVAR$p.value), ".\n")

cat(" Shapiro-Francia normality test for the bootstrap distribution of IVAR:\n")

cat(" W’ =", rnd(w_IVAR$statistic[1]), ", z =", rnd(qnorm(w_IVAR$p.value, mean = 0, sd = 1, lower.tail = FALSE)), ", p-value =", rnd(w_IVAR$p.value), ".\n")

cat(" DBD normality test for the bootstrap distribution of IVAR:\n")

cat(" K^2 =", rnd(k2_IVAR), ", p-value =", rnd(p_value_IVAR), ".\n")

# ------------------------------------------------------------------------------------------------------------

# Interindividual Standard Deviation (ISD)

# ------------------------------------------------------------------------------------------------------------

# Tests for symmetry, kurtosis, and normality for the bootstrap distribution of ISD

# D’Agostino skewness test

sq_b1_ISD <- agostino.test(bd$t[,4])

# Anscombe–Gynn kurtosis test

b2_ISD <- anscombe.test(bd$t[,4])

# Shapiro–Francia normality test

w_ISD <- sf.test(bd$t[,4])

# D’Agostino–Belanger–D’Agostino normality test

k2_ISD <- sq_b1_ISD$statistic[2]^2 + b2_ISD$statistic[2]^2

p_value_ISD <- pchisq(k2_ISD, df = 2, lower.tail = FALSE)

# Display or save the normal QQ plot to evaluate normality graphically

# Remove the next three hash symbols to save the graph as a JPG file.

# jpeg("QQ_plot_ISD.jpg", width = 3, height = 2, units = "in", res = 300)

# par(mar = c(4, 4, 1, 1))

qqnorm(bd$t[,4], main = "Normal QQ plot for ISD", xlab = "Theoretical quantiles: normal", ylab = "Sample quantiles: boot", cex.main = 0.8, cex.axis = 0.7, cex.lab = 0.7)

qqline(bd$t[,4], col = "darkblue")

# dev.off()

# Display or save density histogram with overlaid density curve.

# Remove the next three hash symbols to save the graph as a JPG file.

# jpeg("Histogram_ISD.jpg", width = 3, height = 2, units = "in", res = 300)

# par(mar = c(4, 4, 1, 1))

hist(bd$t[,4], breaks = 50, probability = TRUE, col = "lightyellow", main = "", xlab = "ISD", cex.axis = 0.8, cex.lab = 0.7)

lines(density(bd$t[,4]), col = "orange", lwd = 2)

x_values <- seq(mean(bd$t[,4]) - 4 * sd(bd$t[,4]), mean(bd$t[,4]) + 4 * sd(bd$t[,4]), length = 1000)

y_values <- dnorm(x_values, mean = mean(bd$t[,4]), sd = sd(bd$t[,4]))

lines(x_values, y_values, col = "darkred", lwd = 2)

legend("topright", title = "Density curves", legend = c("Empiric", "Normal"), fill = c("orange", "darkred"), cex = 0.6, bty = "n")

# dev.off()

# 95% BCa confidence interval for ISD

boot_ci_ISD <- boot.ci(bd, type = "bca", index = 4)

boot_width_ISD <- boot_ci_ISD$bca[5] - boot_ci_ISD$bca[4]

# Print the results about the ISD on the console.

cat("\nInterindividual Standard Deviation: ISD = √(IVAR) = √[mean((xi - xj)^2)] = √(2) * sd(x):\n")

cat(" ISD =", rnd(ISD), ".\n")

cat(" Asymptotic standard error:", rnd(asympt_se_ISD), ".\n")

cat(" 95% Wald-type confidence interval: (", rnd(asympt_ci_ISD[1]), "," , rnd(asympt_ci_ISD[2]), "); width =", rnd(asympt_width_ISD), ".\n")

cat(" Bootstrap standard error:", rnd(sd(bd$t[,4])), ".\n")

cat(" Bootstrap bias:", rnd(mean(bd$t[,4]) - ISD), ".\n")

cat(" 95% BCa confidence interval (", rnd(boot_ci_ISD$bca[4]), ",", rnd(boot_ci_ISD$bca[5]), "); width =", rnd(boot_width_ISD), ".\n")

cat(" D’Agostino skewness test for the bootstrap distribution of ISD:\n")

cat(" √b1 =", rnd(sq_b1_ISD$statistic[1]), ", z =", rnd(sq_b1_ISD$statistic[2]), ", p-value =", rnd(sq_b1_ISD$p.value), ".\n")

cat(" Anscombe-Glynn kurtosis test for the Bootstrap distribution of ISD:\n")

cat(" b2 =", rnd(b2_ISD$statistic[1]), ", z =", rnd(b2_ISD$statistic[2]), ", p-value =", rnd(b2_ISD$p.value), ".\n")

cat(" Shapiro-Francia normality test for the Bootstrap distribution of ISD:\n")

cat(" W’ =", rnd(w_ISD$statistic[1]), ", z =", rnd(qnorm(w_ISD$p.value, mean = 0, sd = 1, lower.tail = FALSE)), ", p-value =", rnd(w_ISD$p.value), ".\n")

cat(" DBD normality test for the Bootstrap distribution of ISD:\n")

cat(" K^2 =", rnd(k2_ISD), ", p-value =", rnd(p_value_ISD), ".\n")

# ------------------------------------------------------------------------------------------------------------

# Interindividual Coefficient of Variation (ICV)

# ------------------------------------------------------------------------------------------------------------

# Tests for symmetry, kurtosis, and normality for the bootstrap distribution of ICV

# D’Agostino skewness test

sq_b1_ICV <- agostino.test(bd$t[,5])

# Anscombe–Gynn kurtosis test

b2_ICV <- anscombe.test(bd$t[,5])

# Shapiro–Francia normality test

w_ICV <- sf.test(bd$t[,5])

# D’Agostino–Belanger–D’Agostino normality test

k2_ICV <- sq_b1_ICV$statistic[2]^2 + b2_ICV$statistic[2]^2

p_value_ICV <- pchisq(k2_ICV, df = 2, lower.tail = FALSE)

# Display or save the normal QQ plot to evaluate normality graphically

# Remove the next three hash symbols to save the graph as a JPG file.

# jpeg("QQ_plot_ICV.jpg", width = 3, height = 2, units = "in", res = 300)

# par(mar = c(4, 4, 1, 1))

qqnorm(bd$t[,5], main = "Normal QQ plot for ICV", xlab = "Theoretical quantiles: normal", ylab = "Sample quantiles: boot", cex.main = 0.8, cex.axis = 0.7, cex.lab = 0.7)

qqline(bd$t[,5], col = "darkblue")

# dev.off()

# Display or save density histogram with overlaid density curve.

# Remove the next three hash symbols to save the graph as a JPG file.

# jpeg("Histogram_ICV.jpg", width = 3, height = 2, units = "in", res = 300)

# par(mar = c(4, 4, 1, 1))

hist(bd$t[,5], breaks = 50, probability = TRUE, col = "lightyellow", main = "", xlab = "ICV", cex.axis = 0.8, cex.lab = 0.7)

lines(density(bd$t[,5]), col = "orange", lwd = 2)

x_values <- seq(mean(bd$t[,5])-4*sd(bd$t[,5]), mean(bd$t[,5])+4*sd(bd$t[,5]), length=1000)

y_values <- dnorm(x_values, mean = mean(bd$t[,5]), sd = sd(bd$t[,5]))

lines(x_values, y_values, col = "darkred", lwd = 2)

legend("topright", title = "Density curves", legend = c("Empiric", "Normal"), fill = c("orange", "darkred"), cex = 0.6, bty = "n")

# dev.off()

# 95% BCa confidence interval for ICV

boot_ci_ICV <- boot.ci(bd, type = "bca", index = 5)

boot_width_ICV <- boot_ci_ICV$bca[5] - boot_ci_ICV$bca[4]

# Print the results about the ICV on the console.

cat("\nInterindividual Coefficient of Variation: ICV = 100 × ISD / (max(x) - min(x)):\n")

cat(" ICV =", rnd(ICV), ".\n")

cat(" Bootstrap standard error:", rnd(sd(bd$t[,5])), ".\n")

cat(" Bootstrap bias:", rnd(mean(bd$t[,5]) - ICV), ".\n")

cat(" 95% BCa confidence interval (", rnd(boot_ci_ICV$bca[4]), ",", rnd(boot_ci_ICV$bca[5]), "); width =", rnd(boot_width_ICV), ".\n")

cat(" D’Agostino skewness test for the Bootstrap distribution of ICV:\n")

cat(" √b1 =", rnd(sq_b1_ICV$statistic[1]), ", z =", rnd(sq_b1_ICV$statistic[2]), ", p-value =", rnd(sq_b1_ICV$p.value), ".\n")

cat(" Anscombe-Glynn kurtosis test for the Bootstrap distribution of ICV:\n")

cat(" b2 =", rnd(b2_ICV$statistic[1]), ", z =", rnd(b2_ICV$statistic[2]), ", p-value =", rnd(b2_ICV$p.value), ".\n")

cat(" Shapiro-Francia normality test for the Bootstrap distribution of ICV:\n")

cat(" W’ =", rnd(w_ICV$statistic[1]), ", z =", rnd(qnorm(w_ICV$p.value, mean = 0, sd = 1, lower.tail = FALSE)), ", p-value =", rnd(w_ICV$p.value), ".\n")

cat(" DBD normality test for the Bootstrap distribution of ICV:\n")

cat(" K^2 =", rnd(k2_ICV), ", p-value =", rnd(p_value_ICV), ".\n")

Conflicts of Interest

The author declares that there are no conflicts of interest.

References

[1]	Hinton, P.R. (2024) Statistics Explained. 4th Edition, Routledge.[CrossRef]
[2]	Alabi, O. and Bukola, T. (2023) Introduction to Descriptive Statistics. In: Kumar, S., Ed., Recent Advances in Biostatistics, IntechOpen, 1-12.[CrossRef]
[3]	Kvålseth, T.O. (1988) Measuring Variation for Nominal Data. Bulletin of the Psychonomic Society, 26, 433-436.[CrossRef]
[4]	San Martin-Castellanos, R., Espinosa-Gil, L. and Fernández-Pedreira, L. (1990) Psicoestadística Descriptive. 2nd Edition, Pirámide.
[5]	Casella, G. and Berger, R. (2024) Statistical Inference. 2nd Edition, Chapman and Hall/CRC.[CrossRef]
[6]	Merce, E. (2009) Pearson’s Coefficient of Variation, an Erroneous History in Assessing the Degree of Significance of the Mean Value. Bulletin UASVM Horticulture, 66, 295 300.
[7]	Pearson, K. and Filon, L.N.G. (1898) Mathematical Contributions to the Theory of Evolution. IV. On the Probable Errors of Frequency Constants and on the Influence of Random Selection on Variation and Correlation. Proceedings of the Royal Society of London, 62, 173-176.[CrossRef]
[8]	Sokal, R.R. and Rohlf, F.J. (1995) Biometry: The Principles and Practice of Statistics in Biological Research. 3rd Edition, Freeman and Company.
[9]	Henze, N. (2024) Limit Theorems for U-Statistics. In: Henze, N., Ed., Asymptotic Stochastics, Springer, 103-141.[CrossRef]
[10]	Korolyuk, V.S. and Borovskich, Y.V. (2013) Theory of U-Statistics (Mathematics and Its Applications, Vol. 273). Springer Science & Business Media.
[11]	Hoeffding, W. (1948) A Class of Statistics with Asymptotically Normal Distribution. The Annals of Mathematical Statistics, 19, 293-325.[CrossRef]
[12]	Wendler, M. (2011) Bahadur Representation for U-Quantiles of Dependent Data. Journal of Multivariate Analysis, 102, 1064-1079.[CrossRef]
[13]	Silverman, B.W. (1998) Density Estimation for Statistics and Data Analysis. Routledge.[CrossRef]
[14]	Ichimura, H. and Newey, W.K. (2022) The Influence Function of Semiparametric Estimators. Quantitative Economics, 13, 29-61.[CrossRef]
[15]	Stuart, A. and Ord, K. (1994) Kendall’s Advanced Theory of Statistics. Volume 1. Distribution Theory. 6th Edition, Edward Arnold.
[16]	Andrianov, I.V. and Awrejcewicz, J. (2024) Asymptotic Methods for Engineers. CRC Press.[CrossRef]
[17]	Clark, J.M. and Warr, R.L. (2023) Bootstrapping through Discrete Convolutional Methods. Applied Stochastic Models in Business and Industry, 40, 144-160.[CrossRef]
[18]	Van der Vaart, A.W. (1998) Asymptotic Statistics. Cambridge University Press.[CrossRef]
[19]	Leadbetter, R.M., Lindgren, G. and Rootzén, H. (1983) Extremes and Related Properties of Random Sequences and Processes. Springer.[CrossRef]
[20]	Efron, B. and Tibshirani, R.J. (1993) An Introduction to the Bootstrap. Chapman & Hall.
[21]	Microsoft Corporation (2024) Microsoft Excel (Versión 2409) [Software]. Microsoft. https://www.microsoft.com/es-mx/microsoft-365/get-started-with-office-2024
[22]	R Core Team (2025) R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. https://www.R-project.org/
[23]	Estabrook, R., Grimm, K.J. and Bowles, R.P. (2012) A Monte Carlo Simulation Study of the Reliability of Intraindividual Variability. Psychology and Aging, 27, 560-576.[CrossRef] [PubMed]
[24]	Kohoutová, L., Atlas, L.Y., Büchel, C., Buhle, J.T., Geuter, S., Jepma, M., et al. (2022) Individual Variability in Brain Representations of Pain. Nature Neuroscience, 25, 749-759.[CrossRef] [PubMed]
[25]	Chen, L., Zhernakova, D.V., Kurilshikov, A., Andreu-Sánchez, S., Wang, D., Augustijn, H.E., et al. (2022) Influence of the Microbiome, Diet and Genetics on Inter-Individual Variation in the Human Plasma Metabolome. Nature Medicine, 28, 2333-2343.[CrossRef] [PubMed]
[26]	Giorgi, F.M., Ceraolo, C. and Mercatelli, D. (2022) The R Language: An Engine for Bioinformatics and Data Science. Life, 12, Article 648.[CrossRef] [PubMed]
[27]	Wehrspohn, U. and Ernst, D. (2022) When Do I Take Which Distribution? A Statis-tical Basis for Entrepreneurial Applications. Springer.[CrossRef]
[28]	Schultzberg, M. and Ankargren, S. (2022) Resampling-free Bootstrap Inference for Quantiles. In: Arai, K., Ed., Proceedings of the Future Technologies Conference (FTC) 2022, Volume 1. FTC 2022 2022, Springer, 548-562.[CrossRef]
[29]	Wegner, L. and Wendler, M. (2024) Robust Change-Point Detection for Functional Time Series Based on U-Statistics and Dependent Wild Bootstrap. Statistical Papers, 65, 4767-4810.[CrossRef]
[30]	Efron, B. and Narasimhan, B. (2025) Packages ‘Bcaboot’. Bias Corrected Bootstrap Confidence Intervals. https://cran.r-project.org/web/packages/bcaboot/bcaboot.pdf
[31]	Mokhtar, S.F., Md Yusof, Z. and Sapiri, H. (2023) Confidence Intervals by Bootstrapping Approach: A Significance Review. Malaysian Journal of Fundamental and Applied Sciences, 19, 30-42.[CrossRef]
[32]	Grün, B. and Miljkovic, T. (2023) The Automated Bias-Corrected and Accelerated Bootstrap Confidence Intervals for Risk Measures. North American Actuarial Journal, 27, 731-750.[CrossRef]
[33]	Efron, B. and Narasimhan, B. (2020) The Automatic Construction of Bootstrap Confidence Intervals. Journal of Computational and Graphical Statistics, 29, 608-619.[CrossRef] [PubMed]
[34]	Wu, S., Zhu, X. and Wang, H. (2023) Subsampling and Jackknifing: A Practically Convenient Solution for Large Data Analysis with Limited Computational Resources. Statistica Sinica, 33, 2041-2064.[CrossRef]
[35]	D’Agostino, R.B. (1970) Transformation to Normality of the Null Distribution of G1. Biometrika, 57, 679-681.[CrossRef]
[36]	Anscombe, F.J. and Glynn, W.J. (1983) Distribution of the Kurtosis Statistic b2 for Normal Samples. Biometrika, 70, 227-234.[CrossRef]
[37]	Royston, P. (1993) A Toolkit for Testing for Non-Normality in Complete and Censored Samples. The Statistician, 42, 37-43.[CrossRef]
[38]	Shapiro, S.S. and Francia, R.S. (1972) An Approximate Analysis of Variance Test for Normality. Journal of the American Statistical Association, 67, 215-216.[CrossRef]
[39]	D’Agostino, R.B., Belanger, A. and D’agostino, R.B. (1990) A Suggestion for Using Powerful and Informative Tests of Normality. The American Statistician, 44, 316-321.[CrossRef]
[40]	Canty, A., Ripley, B. and Brazzale, A.R. (2025) ‘Boot’: Bootstrap Functions.[CrossRef]
[41]	Komsta, L. and Novomestky, F. (2022) ‘Moments’: Moments, Cumulants, Skewness, Kurtosis and Related Tests.[CrossRef]
[42]	Gross, J. and Ligges, U. (2015) ‘Nortest’: Tests for Normality.[CrossRef]
[43]	Pouillot, R. Delignette-Muller, M.L., Denis, J.B., Chen, Y and Havelaar, A. (2024) ‘MC2D’: Tools for Two-Dimensional Monte-Carlo Simulations.[CrossRef]
[44]	Clark, C.E. (1962) Letter to the Editor—The PERT Model for the Distribution of an Activity Time. Operations Research, 10, 405-406.[CrossRef]
[45]	American Psychiatric Association (2022) Diagnostic and Statistical Manual of Mental Disorders. 5th Edition, Text Revised, American Psychiatric Association Publishing.[CrossRef]
[46]	Matthys, W. and Schutter, D.J.L.G. (2022) Improving Our Understanding of Impaired Social Problem-Solving in Children and Adolescents with Conduct Problems: Implications for Cognitive Behavioral Therapy. Clinical Child and Family Psychology Review, 25, 552-572.[CrossRef] [PubMed]
[47]	Tukey, J.W. (1977) Exploratory Data Analysis. Addison-Wesley.
[48]	Chattamvelli, R. and Shanmugam, R. (2023) Descriptive Statistics for Scientists and Engineers (Synthesis Lectures on Mathematics and Statistics). 2th Edition, Springer. [Google Scholar] [CrossRef]
[49]	Nikolic, B. and Popovic, T. (2024) Variables, Data, and Descriptive Statistical Methods—Part I. Medicinski pregled, 77, 325-330.[CrossRef]
[50]	Zhang, X., Astivia, O.L.O., Kroc, E. and Zumbo, B.D. (2023) How to Think Clearly about the Central Limit Theorem. Psychological Methods, 28, 1427-1445.[CrossRef] [PubMed]
[51]	Chernozhuokov, V., Chetverikov, D., Kato, K. and Koike, Y. (2022) Improved Central Limit Theorem and Bootstrap Approximations in High Dimensions. The Annals of Statistics, 50, 2562-2586.[CrossRef]
[52]	Bera, A.K. and Koley, M. (2023) A History of the Delta Method and Some New Results. Sankhya B, 85, 272-306.[CrossRef]