Chapter 3: Tests for Comparing Several Normal Quantiles and Pairwise Confidence Intervals

1. Introduction

Comparing several groups or populations is a fundamental problem in statistics and is often addressed by testing equality of means. However, the mean or median does not fully characterize a distribution, and important differences may occur in other parts of the distribution. In many applications – particularly in medical and reliability studies—the behavior of a large proportion of the population, reflected in specific percentiles, may be more relevant than the mean. For example, a smaller upper quartile for treatment time under one therapy indicates that most patients experience faster recovery, even if mean times are similar. As noted by [cox1985testing] and [li2012comparison], distributions may share similar means while differing substantially in their tails, motivating inference procedures based on quantiles. The problem of estimating/testing the difference between percentiles of two independent normal distributions has received some attention in the literature. However, the paper by [li2012comparison] seems to be the first one considered the problem of testing equality of quantiles of several normal populations. [malekzadeh2023simultaneous] have proposed simultaneous CIs for quantile differences of several normal populations. They presumed that the pivotal quantity found in [chakraborti2007confidence] is different or better than the classical NCT CI, and proposed methods of computing critical values to find the Chakraborti-Li CI. We see that the pairwise fiducial CIs and parametric bootstrap CIs based on the classical NCT pivotal quantity for a quantile and those proposed in [malekzadeh2023simultaneous] are essentially the same.

The rest of this chapter is organized as follows. In the following section, we describe the generalized variable test (GVT) by [li2012comparison]. We have enhanced the GVT by deriving theoretical expressions of some quantities. This closed-form expressions are easy to compute and thereby avoid additional simulation used in [li2012comparison]. Then we describe the modified LRT proposed in [abdollahnezhad2018testing] and a new modified MLRT. We also outline a parametric bootstrap (PB) approach for testing the equality of the normal quantiles.

In Section 3, we evaluate and compare the tests in terms of type I error rates and powers. Available pairwise CIs and simplified version of them are described in Section 4. These pairwise CIs are evaluated and compared using Monte Carlo simulation. In Section 5, we illustrate the tests and simultaneous CIs using two examples with real data.

2. Tests for Equality of Quantiles

Let $X_{i1},\mathrm{\dots },X_{i{n}_i}$ be a sample from a normal distribution with parameters ${\mu }_i$ and ${\sigma }_i^2$, say, $N({\mu }_i,{\sigma }_i^2)$, $i=1,\mathrm{\dots },k$. The $p$th quantile of the $i$th normal distribution is given by ${\xi }_i={\mu }_i+z_p{\sigma }_i$, where $z_p$ is the standard normal quantile. The problem of interest is to test the equality of the quantiles. Specifically, we like to test

(1) $$H_0:{\xi }_1=\mathrm{\dots }={\xi }_k\mathrm{vs}.H_a:{\xi }_i\ne {\xi }_j\text{for some }i\ne j.$$

Let $({\overline{X}}_i,S_i^2)$ denote the (mean, variance) based on a sample of size $n_i$, and let $({\overline{x}}_i,s_i^2)$ be an observed value of $({\overline{X}}_i,S_i^2)$ and $m_i=n_i-1$, $i=1,\mathrm{\dots },k$.

2.1. Fiducial Quantities for Quantiles

Let $(\overline{X},S^2)$ denote the (mean, variance) based on a sample of size $n$ from a $N(\mu ,{\sigma }^2)$ distribution. Let $(\overline{x},s^2)$ be an observed value of $(\overline{X},S^2)$, and let $\xi =\mu +z_p\sigma$. Using the stochastic representations that $\overline{X}=\mu +Z\sqrt{\sigma /n}$ and $S^2={\sigma }^2{U}^2$, where $Z\sim N(0,1)$ independently of $U^2={\chi }_m^2/m$, we see that

$$T=\frac{\xi -\overline{X}}{S/\sqrt{n}}\stackrel{d}{=}\frac{z_p\sigma +Z\frac{\sigma }{\sqrt{n}}}{\sigma U}=\frac{Z+z_p\sqrt{n}}{U}\stackrel{d}{=}{t}_m(z_p\sqrt{n}),$$

where $t_m(\delta )$ denotes the noncentral $t$ random variable with df = $m$ and the noncentrality parameter $\delta$. To get the 2nd equation, we used the results that $\mu -\overline{X}\sim \frac{\sigma }{\sqrt{n}}Z$.

Applying the [dawid1982functional] approach, we find a fiducial quantity (FQ) for $\xi$ by solving the “equation” (2.1) for $\xi$ and then replacing $(\overline{X},S)$ by the observed value $(\overline{x},s)$, and is given by

(2) $$Q_{\xi }=\overline{x}+t_m(z_p\sqrt{n})\frac{s}{\sqrt{n}}.$$

The FQ for $\xi$ given in [li2012comparison] is given by

(3) $$Q_{\xi }=\overline{x}-\frac{Z}{U}\frac{s}{\sqrt{n}}+z_p\frac{s}{U}=\overline{x}+\left(\frac{Z+z_p\sqrt{n}}{U}\right)\frac{s}{\sqrt{n}},$$

where $Z\sim N(0,1)$ independently of $U^2\sim {\chi }_m^2/m$. To get the 2nd equation in the above, we used the fact that $Z$ and $-Z$ are identically distributed. Note that the term within the brackets has a $t_m(z_p\sqrt{n})$ distribution. [li2012comparison] have used the expression (3) for $Q_{\xi }$ and unable to find the mean and the variance of $Q_{\xi }$ which are required to develop a test for equality of quantiles. The expression for $Q_{\xi }$ in (2) is simple and its moments can be found using the moments of a $t_m(\delta )$ distribution.

2.2. Generalized Variable Test

We shall now describe the generalized variable test of [li2012comparison] using the FQ (2). Let

(4) $$𝝃={({\xi }_1,\mathrm{\dots },{\xi }_k)}^{\prime },{𝑸}_{\xi }={(Q_{{\xi }_1},\mathrm{\dots },Q_{{\xi }_k})}^{\prime }\mathrm{and}𝐇={({𝐈}_{k-1},-\mathrm{𝟏})}_{(k-1)\times k},$$

where

(5) $$Q_{{\xi }_i}={\overline{x}}_i+t_{m_i}(z_p\sqrt{n_i})\frac{s_i}{\sqrt{n_i}},i=1,\mathrm{\dots },k,$$

and ${𝐈}_l$ denotes the identity matrix of order $l$ and $\mathrm{𝟏}$ is a $(k-1)\times 1$ vector of ones. In terms of these notations, the null hypothesis in (1) can be written as

(6) $$H_0:𝐇𝝃=\mathrm{𝟎}\mathrm{vs}.𝐇𝝃\ne \mathrm{𝟎}.$$

For a given $({\overline{x}}_1,\mathrm{\dots },{\overline{x}}_k)$ and $(s_1^2,\mathrm{\dots },s_k^2)$, let ${𝝁}_{\xi }=E\left({𝑸}_{\xi }\right)\mathrm{and}{𝚺}_{\xi }=\mathrm{Cov}({𝑸}_{\xi })$. The mean vector and covariance matrix can be found using the moments of noncentral $t$ random variables. Recall that

(7) $$E(t_m(\delta ))=\sqrt{m/2}\frac{\mathrm{\Gamma }((m-1)/2)}{\mathrm{\Gamma }(m/2)}\delta \mathrm{and}\mathrm{Var}(t_m(\delta ))=\frac{m}{m-2}(1+{\delta }^2)-{[E(t_m(\delta ))]}^2.$$

Let $a_i=\sqrt{m_i/2}\frac{\mathrm{\Gamma }((m_i-1)/2)}{\mathrm{\Gamma }(m_i/2)},i=1,\mathrm{\dots },k.$ Then

(8) $${\mu }_{{\xi }_i}=E\left({𝑸}_{{\xi }_i}\right)={\overline{x}}_i+E\left(t_{m_i}(z_p\sqrt{n_i})\right)\frac{s_i}{\sqrt{n_i}}={\overline{x}}_i+z_p{a}_i{s}_i,i=1,\mathrm{\dots },k,$$

and

(9) $${\sigma }_{{\xi }_i}^2=\mathrm{Var}\left({𝑸}_{{\xi }_i}\right)=\mathrm{Var}\left(t_{m_i}(z_p\sqrt{n_i})\right)\frac{s_i^2}{n_i}=\left[\frac{m_i}{m_i-2}(1+z_p^2{n}_i)-a_i^2{z}_p^2{n}_i\right]\frac{s_i^2}{n_i},i=1,\mathrm{\dots },k,$$

In terms of these quantities, the generalized test variable $T$ is given by

(10) $$T={\left({𝑸}_{\xi }-{𝝁}_{\xi }\right)}^{\prime }{𝐇}^{\prime }{\left(𝐇{𝚺}_{\xi }{𝐇}^{\prime }\right)}^{-1}𝐇\left({𝑸}_{\xi }-{𝝁}_{\xi }\right).$$

Note that, for a given $(\overline{𝐱},𝐬)$, ${𝝁}_{\xi }={({\mu }_{{\xi }_1},\mathrm{\dots },{\mu }_{{\xi }_k})}^{\prime }$ and ${𝚺}_{\xi }=\mathrm{diag}({\sigma }_{{\xi }_1}^2,\mathrm{\dots },{\sigma }_{{\xi }_k}^2)$ can be computed numerically. Let $T_0={𝝁}_{\xi }^{\prime }{𝐇}^{\prime }{\left(𝐇{𝚺}_{\xi }{𝐇}^{\prime }\right)}^{-1}𝐇{𝝁}_{\xi }$. The generalized p-value for testing $H_0$ in (6) is given by

$$\text{p-value}=P\left(T\ge T_0|H_0\right).$$

The above generalized p-value can be estimated using the following Algorithm 3.1.

Algorithm 3.1

For a given set of $k$ independent samples from normal populations,

1. compute the means ${\overline{x}}_i$’s and variances $s_i^2$’s.

2. Generate noncentral $t$ random variables $t_{m_1}(z_p\sqrt{n_1}),\mathrm{\dots },t_{m_k}(z_p\sqrt{n_k})$, and compute ${𝑸}_{\xi }={(Q_{{\xi }_1},\mathrm{\dots },Q_{{\xi }_k})}^{\prime }$ using (5).

3. Compute ${𝝁}_{\xi }$ using (8) and ${𝚺}_{\xi }$ using (9).

4. Compute $T={\left({𝑸}_{\xi }-{𝝁}_{\xi }\right)}^{\prime }{𝐇}^{\prime }{\left(𝐇{𝚺}_{\xi }{𝐇}^{\prime }\right)}^{-1}𝐇\left({𝑸}_{\xi }-{𝝁}_{\xi }\right).$

5. Repeat steps 2, 3, and 4 for large number of times, say, $10,000$.

6. Compute $T_0={𝝁}_{\xi }^{\prime }{𝐇}^{\prime }{\left(𝐇{𝚺}_{\xi }{𝐇}^{\prime }\right)}^{-1}𝐇{𝝁}_{\xi }$.

7. Find the proportion of times $T\ge T_0$. This proportion is an estimate of the generalized p-value.

The null hypothesis (1) is rejected at the level $\alpha$, if the generalized p-value is less than $\alpha$.

2.3. Likelihood Ratio Test

Define

(11) $${\overline{X}}_i=\frac{1}{n_i}\sum _{j=1}^{n_i}{X}_{ij}\mathrm{and}{\widehat{\sigma }}_i^2=\frac{1}{n_i}\sum _{j=1}^{n_i}{(X_{ij}-{\overline{X}}_i)}^2,i=1,\mathrm{\dots },k.$$

The log-likelihood function is given by

(12) $$\ln L({\mu }_1,\mathrm{\dots },{\mu }_k;{\sigma }_1^2,\mathrm{\dots },{\sigma }_k^2)=-\frac{1}{2}\sum _{i=1}^k{n}_i\ln ({\sigma }_i^2)-\frac{1}{2}\sum _{i=1}^k\frac{n_i{\widehat{\sigma }}_i^2+n_i{({\overline{X}}_i-{\mu }_i)}^2}{{\sigma }_i^2}.$$

The MLEs that maximize the above $\ln L$ are given by ${\widehat{\mu }}_i={\overline{X}}_i$ and ${\widehat{\sigma }}_i^2$, $i=1,\mathrm{\dots },k$.

The log-likelihood function under $H_0:{\xi }_1=\mathrm{\dots }={\xi }_k$ is given by

(13) $$\ln L(\xi -z_p{\sigma }_1,\mathrm{\dots },\xi -z_p{\sigma }_k;{\sigma }_1^2,\mathrm{\dots },{\sigma }_k^2)=-\frac{1}{2}\sum _{i=1}^k{n}_i\ln ({\sigma }_i^2)-\frac{1}{2}\sum _{i=1}^k\frac{n_i{\widehat{\sigma }}_i^2+n_i{({\overline{X}}_i-\xi +z_p{\sigma }_i)}^2}{{\sigma }_i^2},$$

where $\xi$ is the common unknown quantile under $H_0$ in (1). The values of $(\xi ,{\sigma }_1^2,\mathrm{\dots },{\sigma }_k^2)$ that maximize (13) are the constrained MLEs, and let us denote the constrained MLEs by $({\widehat{\xi }}_c,{\widehat{\sigma }}_{1c}^2,\mathrm{\dots },{\widehat{\sigma }}_{kc}^2)$. Details on calculation of the constrained MLEs and an algorithm are given in the appendix. The LRT statistic is expressed as

(14) $\mathrm{\Lambda }$ $=$ $2\left[\ln L({\widehat{\mu }}_1,\mathrm{\dots },{\widehat{\mu }}_k;{\widehat{\sigma }}_1^2,\mathrm{\dots },{\widehat{\sigma }}_k^2)-\ln L({\widehat{\xi }}_c-z_p{\widehat{\sigma }}_{1c},\mathrm{\dots },{\widehat{\xi }}_c-z_p{\widehat{\sigma }}_{kc};{\widehat{\sigma }}_{1c}^2,\mathrm{\dots },{\widehat{\sigma }}_{kc}^2)\right]$

$=$ ${\displaystyle \sum _{i=1}^k}{\displaystyle \frac{n_i}{{\widehat{\sigma }}_{ic}^2}}[{\widehat{\sigma }}_i^2+{({\overline{X}}_i+z_p{\widehat{\sigma }}_{ic}-{\widehat{\xi }}_c)}^2]+{\displaystyle \sum _{i=1}^k}{n}_i\ln \left({\displaystyle \frac{{\widehat{\sigma }}_{ic}^2}{{\widehat{\sigma }}_i^2}}\right)-{\displaystyle \sum _{i=1}^k}{n}_i.$

For a given level of significance $\alpha$, the LRT rejects the null hypothesis when $\mathrm{\Lambda }>{\chi }_{k-1;1-\alpha }^2$, where ${\chi }_{m;q}^2$ denotes the 100$q$ percentile of the chi-square distribution with df = $m$.

AJ Modified LRT

In general, the LRT is not accurate for small samples. To improve the LRT, [abdollahnezhad2018testing] have proposed a modification using the general theory of [skovgaard2001likelihood]. To describe this modification, let

$$𝜷={({\mu }_1/{\sigma }_1^2,\mathrm{\dots },{\mu }_k/{\sigma }_k^2,1/{\sigma }_1^2,\mathrm{\dots },1/{\sigma }_k^2)}^{\prime },$$

and

$$𝝉={(n_1{\mu }_1,\mathrm{\dots },n_k{\mu }_k,-0.5n_1({\mu }_1^2+{\sigma }_1^2),\mathrm{\dots },-0.5n_k({\mu }_k^2+{\sigma }_k^2))}^{\prime }$$

Let

$$\widehat{}𝜷={({\overline{x}}_1/{\widehat{\sigma }}_1^2,\mathrm{\dots },{\overline{x}}_k/{\widehat{\sigma }}_k^2,1/{\widehat{\sigma }}_1^2,\mathrm{\dots },1/{\widehat{\sigma }}_k^2)}^{\prime }$$

be the MLE of $𝜷$ and

$$\widehat{}𝝉={(n_1{\overline{x}}_1,\mathrm{\dots },n_k{\overline{x}}_k,-0.5n_1({\overline{x}}_1^2+{\widehat{\sigma }}_1^2),\mathrm{\dots },-0.5n_k({\overline{x}}_k^2+{\widehat{\sigma }}_k^2))}^{\prime }$$

be the MLE of $𝝉$. The constrained MLE of $𝜷$ is given by $\widehat{}{𝜷}_c$, this can be found by replacing the parameters with the constrained MLEs ${\widehat{\sigma }}_{ic}$ and ${\widehat{\mu }}_{ic}={\widehat{\xi }}_c-z_p{\widehat{\sigma }}_{ic}$, $i=1,\mathrm{\dots },k$. Similarly, we can find the constrained MLE $\widehat{}{𝝉}_c$. Furthermore, define

Let $\widehat{}𝚺$ and $\widehat{}{𝚺}_c$ denote the MLE and the constrained MLE of $𝚺$, respectively. Let

$$\gamma =\frac{{\left\{{\left(\widehat{}𝝉-\widehat{}{𝝉}_c\right)}^{\prime }\widehat{}{𝚺}_c^{-1}\left(\widehat{}𝝉-\widehat{}{𝝉}_c\right)\right\}}^{k/2}}{{\mathrm{\Lambda }}^{k/2-1}{(\widehat{}𝜷-\widehat{}{𝜷}_c)}^{\prime }\left(\widehat{}𝝉-\widehat{}{𝝉}_c\right)}{\left(\frac{|\widehat{}{𝚺}_c|}{|\widehat{}𝚺|}\right)}^{1/2}.$$

Then ${\widehat{\mathrm{\Lambda }}}^{\ast }=\mathrm{\Lambda }{\left(1-\frac{1}{\mathrm{\Lambda }}\ln (\gamma )\right)}^2\sim {\chi }_{k-1}^2,\mathrm{approximately}.$ The AJ-MLRT rejects the null hypothesis if ${\widehat{\mathrm{\Lambda }}}^{\ast }>{\chi }_{k-1;1-\alpha }^2.$

DK Modified LRT

We can also find another improved version of the likelihood approach, referred to as the DK modified LRT (DK-MLRT), which can be obtained by approximating the distribution of $\mathrm{\Lambda }/{\mu }_{\mathrm{\Lambda }}$ by the moment matching method, where ${\mu }_{\mathrm{\Lambda }}=E(\mathrm{\Lambda })$. Let ${\sigma }_{\mathrm{\Lambda }}^2=\mathrm{var}(\mathrm{\Lambda })$. This approximation is based on the results by [diciccio2001simple] on improving the usual LRT and Welch’s approximate degrees of freedom solution for the Behrens-Fisher problem. We approximate the distribution of $\mathrm{\Lambda }/{\mu }_{\mathrm{\Lambda }}$ by the distribution of ${\chi }_{\nu }^2/\nu$, where the degrees of freedom $\nu$ is found so that $\mathrm{var}\left(\mathrm{\Lambda }/{\mu }_{\mathrm{\Lambda }}\right)={\sigma }_{\mathrm{\Lambda }}^2/{\mu }_{\mathrm{\Lambda }}^2=\mathrm{var}({\chi }_{\nu }^2/\nu ).$ This moment matching method yields $\nu =2{\mu }_{\mathrm{\Lambda }}^2/{\sigma }_{\mathrm{\Lambda }}^2$. In general, it is difficult to find ${\mu }_{\mathrm{\Lambda }}$ and ${\sigma }_{\mathrm{\Lambda }}^2$ theoretically, but they can be replaced by Monte Carlo estimate $({\widehat{\mu }}_{\mathrm{\Lambda }},{\widehat{\sigma }}_{\mathrm{\Lambda }}^2)$ as shown in Algorithm 3.2. For an observed value ${\mathrm{\Lambda }}_0^{\ast }$ of $\mathrm{\Lambda }/{\widehat{\mu }}_{\mathrm{\Lambda }}$, the p-value of the DK-MLRT is given by

(15) $$P\left(\frac{1}{\widehat{\nu }}{\chi }_{\widehat{\nu }}^2>\frac{\mathrm{\Lambda }}{{\widehat{\mu }}_{\mathrm{\Lambda }}}\right),$$

which can be estimated as follows.

Algorithm 3.2

For a given set of $k$ samples, calculate $({\overline{X}}_i,{\widehat{\sigma }}_i^2)$, $i=1,\mathrm{\dots },k$.

1. Calculate the constrained MLEs ${\widehat{\xi }}_c$ and $({\widehat{\sigma }}_{1c}^2,\mathrm{\dots },{\widehat{\sigma }}_{kc}^2)$, and the LRT statistic $\mathrm{\Lambda }$ in (14).

2. Generate a sample of size $n_i$ from the $N({\widehat{\xi }}_c-z_p{\widehat{\sigma }}_{ic},{\widehat{\sigma }}_{ic}^2)$, $i=1,\mathrm{\dots },k$.

3. Calculate the LRT statistic for the samples generated in the previous step.

4. Repeat steps 2 and 3 for a large number of times, say, $10,000$.

5. Find the mean ${\widehat{\mu }}_{\mathrm{\Lambda }}$ and the variance ${\widehat{\sigma }}_{\mathrm{\Lambda }}^2$ of these 10,000 simulated LRT statistics, and compute $\widehat{\nu }=2{\widehat{\mu }}_{\mathrm{\Lambda }}^2/{\widehat{\sigma }}_{\mathrm{\Lambda }}^2$ and let ${\mathrm{\Lambda }}_0^{\ast }=\mathrm{\Lambda }/{\widehat{\mu }}_{\mathrm{\Lambda }}$.

6. The p-value of the MLRT is estimated by $P\left(\frac{1}{\widehat{\nu }}{\chi }_{\widehat{\nu }}^2>\frac{\mathrm{\Lambda }}{{\widehat{\mu }}_{\mathrm{\Lambda }}}\right).$

2.4. The Parametric Bootstrap Test

Let ${\widehat{\xi }}_i={\overline{X}}_i+z_p{S}_i$, $i=1,\mathrm{\dots },k$. It is easy to see that an estimate of the variance of ${\widehat{\xi }}_i$ is given by

$$\widehat{\mathrm{V}}({\widehat{\xi }}_i)=S_i^2\left(\frac{1}{n_i}+z_p^2(1-c_{m_i}^2)\right)\mathrm{with}{c}_{m_i}=\sqrt{\frac{2}{m_i}}\frac{\mathrm{\Gamma }((m_i+1)/2)}{\mathrm{\Gamma }(m_i/2)},i=1,\mathrm{\dots },k.$$

Following the lines of [krishnamoorthy2007parametric], who have developed a PB test for equality of normal means, we can develop a test statistic for testing $H_0$ in (1) as

(16) $T_{\xi }({\overline{X}}_1,\mathrm{\dots },{\overline{X}}_k;S_1^2,\mathrm{\dots },S_k^2)$ $=$ ${\displaystyle \sum _{i=1}^k}{\displaystyle \frac{{\widehat{\xi }}_i^2}{\widehat{\mathrm{V}}({\widehat{\xi }}_i)}}-{\displaystyle \frac{{\left({\sum }_{i=1}^k\frac{1}{\widehat{\mathrm{V}}({\widehat{\xi }}_i)}{\widehat{\xi }}_i\right)}^2}{{\sum }_{i=1}^k\frac{1}{\widehat{\mathrm{V}}({\widehat{\xi }}_i)}}}$

$=$ ${\displaystyle \sum _{i=1}^k}{\displaystyle \frac{1}{\widehat{\mathrm{V}}({\widehat{\xi }}_i)}}\left({\widehat{\xi }}_i^2-{\overline{\widehat{\xi }}}^2\right)$

where $\overline{\widehat{\xi }}={\sum }_{i=1}^k{W}_i{\widehat{\xi }}_i$, $W_i=\frac{{[v({\widehat{\xi }}_i)]}^{-1}}{{\sum }_{j=1}^k{[v({\widehat{\xi }}_i)]}^{-1}}$ and

The parametric bootstrap involves sampling from the estimated models where the model parameters are replaced with the sample estimates. The PB test statistic is the same as the statistic $T_{\xi }$ in (16) with $({\overline{X}}_i,S_i)$ replaced by $({\overline{X}}_i^{\ast },S_i^{\ast })$, $i=1,\mathrm{\dots },k$, which are the statistics based on bootstrap samples generated from $N({\widehat{\mu }}_{ic},{\widehat{\sigma }}_{ic}^2)$, $i=1,\mathrm{\dots },k$, distributions. Here, ${\widehat{\mu }}_{ic}={\widehat{\xi }}_c-z_p{\widehat{\sigma }}_{ic}$, $i=1,\mathrm{\dots },k$, where ${\widehat{\xi }}_c$ and ${\widehat{\sigma }}_{ic}^2$ are the constrained MLEs derived in Section 2.3. To express the PB statistic, let

$${\overline{X}}_i^{\ast }\sim N({\widehat{\mu }}_{ic},{\widehat{\sigma }}_{ic}^2)\text{independently of}{S}_i^{\ast 2}\sim {\widehat{\sigma }}_{ic}^2\frac{{\chi }_{m_i}^2}{m_i},i=1,\mathrm{\dots },k.$$

Also, let ${\widehat{\xi }}_i^{\ast }={\overline{X}}_i^{\ast }+z_p{S}_i^{\ast }$, $i=1,\mathrm{\dots },k$. Then the PB statistic can be expressed as

(17) $$T_{\xi }^{\ast }({\overline{X}}_1^{\ast },\mathrm{\dots },{\overline{X}}_k^{\ast };S_1^{\ast 2},\mathrm{\dots },S_k^{\ast 2})=\sum _{i=1}^k\frac{1}{\widehat{\mathrm{V}}({\widehat{\xi }}_i^{\ast })}\left({\widehat{\xi }}_i^{\ast 2}-{\overline{\widehat{\xi }}}^{\ast 2}\right),$$

where ${\overline{\widehat{\xi }}}^{\ast }={\sum }_{i=1}^k{W}_i^{\ast }{\widehat{\xi }}_i^{\ast }$, $W_i^{\ast }=\frac{{[\widehat{\mathrm{V}}({\widehat{\xi }}_i^{\ast })]}^{-1}}{{\sum }_{j=1}^k{[v({\widehat{\xi }}_i^{\ast })]}^{-1}}$ and $\widehat{\mathrm{V}}({\widehat{\xi }}_i^{\ast })=S_i^{\ast 2}\left(\frac{1}{n}+z_p^2(1-c_{m_i}^2)\right)$ with $c_{m_i}=\sqrt{\frac{2}{m_i}}\frac{\mathrm{\Gamma }((m_i+1)/2)}{\mathrm{\Gamma }(m_i/2)}$, $i=1,\mathrm{\dots },k$.

The PB test for the equality of the quantiles rejects the null hypothesis (1) if

where $({\overline{x}}_i,s_i^2)$ is an observed value of $({\overline{X}}_i,S_i^2)$, $i=1,\mathrm{\dots },k$.

The p-value of the above PB test can be computed using the following Algorithm 3.3.

Algorithm 3.3

For a given set of $({\overline{x}}_1,\mathrm{\dots },{\overline{x}}_k,s_1^2,\mathrm{\dots },s_k^2)$ and the sample sizes,

1. compute the test statistic $T_{\xi }({\overline{x}}_1,\mathrm{\dots },{\overline{x}}_k;s_1^2,\mathrm{\dots },s_k^2)$ using (16).

2. Compute the constrained MLEs, ${\widehat{\xi }}_c$ and ${\sigma }_{ic}^2$, $i=1,\mathrm{\dots },k$. Set ${\widehat{\mu }}_{ic}={\widehat{\xi }}_c-z_p{\widehat{\sigma }}_{ic}$, $i=1,\mathrm{\dots },k$.

3. Generate ${\overline{X}}_i^{\ast }\sim N({\widehat{\mu }}_{ic},{\widehat{\sigma }}_{ic}^2)$ and $S_i^{\ast 2}\sim {\widehat{\sigma }}_{ic}^2\frac{{\chi }_{m_i}^2}{m_i},$ $i=1,\mathrm{\dots },k.$

4. Compute the PB statistic $T_{\xi }^{\ast }({\overline{X}}_1^{\ast },\mathrm{\dots },{\overline{X}}_k^{\ast };S_1^{\ast 2},\mathrm{\dots },S_k^{\ast 2})$ in (17).

5. Repeat steps 3 and 4 for a large number of times, say, 10,000.

6. The proportion of $T_{\xi }^{\ast }$’s that are greater than $T_{\xi }$ is an estimate of the PB p-value for testing $H_0$ in (1)

3. Type I Error Rates and Power Studies

To assess the statistical properties of the tests, we estimated the type I error rates and powers of the tests using Monte Carlo simulation. To estimate the error rates of the GVT, we first generated 10,000 sets of $k$ samples from independent normal populations with assumed means and variances. For each set of samples, we used another 10,000 simulation runs to estimate the p-values using Algorithms 3.1 and 3.3. The PB test was evaluated similarly. To estimate the error rates of the DK MLRT, we used the estimates of ${\widehat{\mu }}_{\mathrm{\Lambda }}$ and ${\widehat{\sigma }}_{\mathrm{\Lambda }}^2$ based on simulation with 1000 runs. The estimates of the type I error rates of the AJ MLRT are based on 10,000 simulation runs. The error rates and powers of (1) generalized variable test (GVT), (2) Abdollahnezhada and Jafari’s LRT (AJ MLRT), (3) Our new MLRT (DK MLRT) and (4) the parametric bootstrap test (PB test).

In Table 1, we reported the estimated sizes of all the tests for $p=0.10,0.25,0.75,0.90$, and $k=3,5,10$. We first observe from this table that the GVT test is conservative having type I error rates less than the nominal level 0.05 for almost all cases. In particular, the GVT is too conservative for $k=10$. The AJ MLRT and DK MLRT perform similar in controlling type I error rates around 0.05, except that the former test is little conservative for $k=10$. In general, we see that the DK MLRT controls the error rates around the nominal level for all sample size and parameter configurations considered in Table 1.

The powers of tests were also estimated for $p=0.25,0.50,0.75$, $k=3,5,10$ and for some values of sample sizes ranging from 5 to 20. The powers are presented in Table 2. Since the GVT is in general conservative, the powers of the GVT are smaller than those of the other tests in most cases. For example, see the powers of the GVT for the case of $k=3$. There is no clear-cut winner between the AJ-MLRT (2) and the DK-MLRT (3) for $k\le 5$. In particular, the AJ-MLRT seems to have larger powers than the DK-MLRT when $p=0.25$ but smaller powers than DK-MLRT when $p=0.75$. These tests perform similar for the case $p=0.5$. However, for $k=10$, the DK-MLRT appears to have better power property than the AJ-MLRT for $p=0.25,0.5$ and 0.75. The DK-MLRT is preferable to the AJ-MLRT when the number of groups being compared is ten or more. The PB test (4) appears to be better than the AJ-MLRT and the DK-MLRT for $k\le 10$ and $p\le 0.5$. However, the PB test is much less powerful than the DK-MLRT when $p=0.75$. For the cases where $p\le 0.5$, the PB test can be recommended for applications.

We now review the performances of the tests when $p=0.5$, that is, for testing equality of the normal means. The PB test simplifies to the one for testing the equality of the means given in [krishnamoorthy2007parametric]. The Monte Carlo simulation studies in their paper indicated that the PB test controls the type I error rates always very close to the nominal level. The PB test is better than the popular Welch test ([welch1951comparison]) for $k\ge 10$. All the tests control the type I error rates within the nominal level 0.05 for $k\le 10$. The PB test is more powerful than the AJ MLRT for almost all the cases. The PB test is also somewhat more powerful than DK MLRT. To compare several normal means, the PB test is referrable to all other tests.

4. Simultaneous CIs for Pairwise Differences

Once the null hypothesis in (1) is rejected, then it is desired to find the population quantiles that are significantly different. This can be found by examining simultaneous confidence intervals for all possible pairs of differences ${\xi }_i-{\xi }_j$, , where ${\xi }_i={\mu }_i+z_p{\sigma }_i$, $i=1,\mathrm{\dots },k$.

4.1. One-Sample Confidence Intervals

Before proceeding to develop simultaneous CIs, we need to choose an appropriate CI for a normal quantile. The classical CI, referred to as the noncentral $t$ interval, is based on the pivotal quantity $\sqrt{n}(\overline{X}-\xi )/S$ which follows a noncentral $t$ distribution (e.g., see Owen, 1968 and Section 5.3.1.1 of [lawless2011statistical]). Since the CI based on $\sqrt{n}(\overline{X}-\xi )/S$ and the one based on $\sqrt{n}(\xi -\overline{X})/S$ are the same, we consider the latter version of the pivotal quantity. Letting $m=n-1$ and using the stochastic representations that $\overline{X}\stackrel{d}{=}\mu +Z\frac{\sigma }{\sqrt{n}}\mathrm{and}{S}^2\stackrel{d}{=}{\sigma }^2{U}^2,$ where $Z\sim N(0,1)$ independently of $U^2\sim {\chi }_m^2/m$ distribution, we see that

(18) ${\displaystyle \frac{\xi -\overline{X}}{S/\sqrt{n}}}\stackrel{d}{=}{\displaystyle \frac{z_p\sqrt{n}-Z}{U}}\stackrel{d}{=}{t}_m(z_p\sqrt{n}),$

where $t_m(\delta )$ denotes the noncentral $t$ random variable with degrees of freedom (df) $m$ and the noncentrality parameter $\delta$. To get the 2nd equation in (18), we used the result that $Z$ and $-Z$ are identically distributed. On the basis of the above distributional result, the $1-2\alpha$ CI for ${\eta }_p$ is given by

(19) $$(\overline{X}+t_{m;\alpha }(z_p\sqrt{n})\frac{S}{\sqrt{n}},\overline{X}+t_{m;1-\alpha }(z_p\sqrt{n})\frac{S}{\sqrt{n}}),$$

where $t_{m;\alpha }(\delta )$ denotes the 100$\alpha$ percentile of $t_m(\delta )$.

[chakraborti2007confidence] have proposed a pivotal quantity $T$ based on the minimum variance unbiased estimator of $\xi$ given by ${\widehat{\xi }}_u=\overline{X}+z_p{c}_{n}S,\mathrm{with}{c}_n=\sqrt{m/2}\mathrm{\Gamma }\left(m/2\right)/\mathrm{\Gamma }(n/2).$ In Chapter LABEL:ch2, we have shown that the CI based on $T$ is the same as the one in (19). To show this, we first note that the variance estimate of the estimator ${\widehat{\xi }}_u$ is given by $\widehat{\mathrm{V}}({\widehat{\xi }}_u)=\frac{S^2}{n}\left(1+n{z}_p^2(c_n^2-1)\right),$ and the pivotal quantity $T$ can be expressed as

(20) $$T=\frac{\xi -(\overline{X}+z_p{c}_{n}S)}{\frac{S}{\sqrt{n}}\sqrt{1+n{z}_p^2(c_n^2-1)}}=\frac{\sqrt{n}(\xi -\overline{X})/S-z_p\sqrt{n}{c}_n}{\sqrt{1+n{z}_p^2(c_n^2-1)}}.$$

Since $T$ is a one-to-one function of the usual pivotal quantity $\sqrt{n}(\xi -\overline{X})/S$, the CI based on $T$ should be the same as the classical CI in (19). For more details, see Chapter LABEL:ch2.

[malekzadeh2023simultaneous] have proposed several methods of constructing pairwise CIs for ${\xi }_i-{\xi }_j$, . Most of these CIs are based on the pivotal quantities of the type in (20). In view of the above discussion, we see that these simultaneous CIs based on the noncentral pivotal quantity (18) and those based on $T$ in (20) are essentially the same.

4.2. FG Method for Simultaneous Pairwise CIs

This procedure is based on the fiducial approach given in [malekzadeh2023simultaneous], and these authors call these resulting CIs as fiducial generalized (FG) CIs. To describe this method, let $a_i=\sqrt{m_i/2}\mathrm{\Gamma }(m_i/2)/\mathrm{\Gamma }(n_i/2)$ and $b_i=1/n_i+z_p^2(a_i^2-1)$, $i=1,\mathrm{\dots },k$. Let ${\widehat{\xi }}_{iu}={\overline{X}}_i+z_p{a}_i{S}_i$ so that $E({\widehat{\xi }}_{iu})={\xi }_i={\mu }_i+z_p{\sigma }_i$ and $\mathrm{var}({\widehat{\xi }}_{iu})=b_i^2{\sigma }^2$. Let ${\widehat{\xi }}_{iju}={\widehat{\xi }}_{iu}-{\widehat{\xi }}_{ju}$, . Using the relation that $t_m(\delta )\sim -t_m(-\delta )$, it is not difficult to check that the FQ $R_i$ for ${\xi }_i$ given in their paper can be expressed as $F_i={\overline{x}}_i+t_{m_i}(z_p\sqrt{n_i})\frac{s_i}{\sqrt{n_i}}$. Let ${\widehat{\xi }}_{iju}={\widehat{\xi }}_{iu}-{\widehat{\xi }}_{ju}$ and $F_{ij}=F_i-F_j$, . The difference ${\widehat{\xi }}_{iju}-F_{ij}$ can be simplified as

(21)

Let , where $v_{ij}=b_i{s}_i^2+b_j{s}_j^2.$ Let $t_{1-\alpha }^F$ denote the $100(1-\alpha )$ percentile of the conditional distribution of $T^F$ given $({\overline{x}}_i,s_i)$’s. Then a $1-\alpha$ simultaneous CIs for ${\xi }_{ij}$’s are given by

(22)

We refer to the above CIs as FGU CIs, because they are based on unbiased estimates of ${\xi }_i$’s.

In view of our discussion in Section 4.1, we can develop pairwise CIs on the basis of the noncentral $t$ pivotal quantity in (18) instead of the one in (20) based on ${\widehat{\xi }}_u$. Such pairwise simultaneous CIs are obtained from (21) and (22) by substituting $a_i=1$ and $b_i=1/n_i$, $i=1,\mathrm{\dots },k$, and the simultaneous CIs can be expressed as

(23)

where $t_{1-\alpha }^{\ast F}$ is the 100$(1-\alpha )$ percentile of and ${\widehat{\xi }}_{ij}-F_{ij}$ is the expression (21) with $a_i=1$, $i=1,\mathrm{\dots },k$. The above CIs, referred at as the FG CIs, are somewhat simpler than the ones in (22), because computation of them does not involve repeated calculation of the gamma function. Our simulation studies in Section 4.5 indicate that these pairwise CIs in (22) and those in (23) are very similar. Also, see Examples in Section 5. However, it seems to be prove theoretically.

Table 1. Sizes of the tests at $\alpha =0.05$

	$p=0.1$				$p=0.25$				$p=0.75$				$p=0.90$
	$({\mu }_1,{\mu }_2,{\mu }_3)=(1,1,1)$
$(n_1,\mathrm{\dots },n_k)$	(1)	(2)	(3)	(4)	(1)	(2)	(3)	(4)	(1)	(2)	(3)	(4)	(1)	(2)	(3)	(4)
	$\xi =-4$				$\xi =-2$				$\xi =2$				$\xi =2$
(5, 10, 11)	.025	.055	.050	.049	.026	.050	.049	.053	.032	.051	.049	.051	.026	.056	.052	.048
(8, 9, 13)	.030	.053	.049	.051	.028	.051	.052	.050	.033	.048	.050	.049	.024	.052	.050	.052
(9,20,5)	.037	.056	.049	.051	.032	.053	.051	.051	.030	.048	.049	.054	.032	.051	.047	.050
(10,10,10)	.023	.053	.049	.049	.033	.045	.049	.047	.038	.048	.048	.047	.022	.046	.048	.051
(10,50,100)	.045	.049	.049	.052	.052	.049	.047	.047	.049	.051	.049	.049	.042	.053	.050	.048
	$({\mu }_1,{\mu }_2,{\mu }_3)=(0,1,1)$
	$\xi =-3$				$\xi =-4$				$\xi =3$				$\xi =2$
(5,10,11)	.024	.052	.050	.051	.028	.051	.051	.050	.037	.051	.049	.048	.029	.052	.047	.048
(8,9,13)	.023	.051	.054	.051	.032	.044	.050	.048	.039	.045	.051	.053	.032	.054	.049	.053
(9,20,5)	.030	.054	.050	.053	.034	.053	.046	.052	.043	.056	.052	.051	.037	.053	.049	.047
(10,10,10)	.025	.053	.048	.049	.028	.049	.051	.046	.046	.049	.049	.048	.025	.049	.051	.050
(10,50,100)	.043	.049	.051	.046	.042	.051	.048	.051	.055	.054	.050	.050	.044	.054	.047	.047
	$({\mu }_1,{\mu }_2,{\mu }_3)=(0,1,2)$
	$\xi =-1$				$\xi =-3$				$\xi =4$				$\xi =3$
(5,10,11)	.021	.053	.050	.047	.025	.045	.048	.048	.027	.053	.047	.047	.032	.053	.051	.049
(8,9,13)	.032	.053	.045	.047	.034	.046	.049	.048	.033	.050	.052	.050	.029	.054	.053	.045
(9,20,5)	.029	.059	.050	.051	.033	.055	.047	.054	.039	.045	.054	.048	.027	.050	.053	.049
(10,10,10)	.033	.052	.048	.050	.029	.045	.051	.049	.032	.044	.047	.047	.024	.052	.047	.052
(10,50,100)	.044	.049	.049	.051	.042	.051	.051	.052	.041	.050	.048	.051	.049	.050	.048	.053
	$({\mu }_1,\mathrm{\dots },{\mu }_5)=(1,1,1,1,1)$
	$\xi =-3$				$\xi =-2$				$\xi =4$				$\xi =5$
(10, 30, 50, 80, 100)	.047	.052	.049	.049	.025	.047	.048	.048	.041	.047	.051	.048	.049	.054	.050	.053
(10, 10, 10, 10, 10)	.030	.050	.047	.049	.034	.052	.049	.051	.025	.047	.047	.049	.032	.047	.050	.049
(5, 10, 15, 20, 30)	.035	.056	.049	.048	.033	.053	.047	.047	.042	.056	.051	.052	.038	.050	.049	.050
	$({\mu }_1,\mathrm{\dots },{\mu }_5)=(1,0,1,0,1)$
(10, 30, 50, 80, 100)	.047	.051	.051	.056	.028	.047	.052	.051	.048	.047	.048	.053	.046	.048	.049	.047
(10, 10, 10, 10, 10)	.025	.048	.049	.051	.030	.049	.051	.051	.031	.046	.049	.050	.031	.048	.053	.051
(5, 10, 15, 20, 30)	.033	.052	.050	.053	.042	.053	.051	.052	.037	.053	.051	.051	.031	.051	.048	.050
	$({\mu }_1,\mathrm{\dots },{\mu }_5)=(2,1,3,2,3)$
(10, 30, 50, 80, 100)	.039	.049	.051	.051	.028	.049	.052	.045	.045	.048	.047	.053	.051	.049	.049	.047
(10, 10, 10, 10, 10)	.024	.054	.051	.052	.030	.046	.053	.051	.050	.046	.046	.049	.032	.046	.051	.054
(5, 10, 15, 20, 30)	.031	.049	.050	.053	.042	.053	.051	.047	.045	.053	.049	.047	.043	.057	.051	.049

1 – GVT; 2 – AJ MLRT; 3 – DK MLRT; 4 – PB test

Table 1 continued.
$p=0.1$ $p=0.25$ $p=0.75$ $p=0.90$ $({\mu }_1,\mathrm{\dots },{\mu }_{10})=(1,\mathrm{\dots },1)$ $(n_1,\mathrm{\dots },n_k)$ (1) (2) (3) (4) (1) (2) (3) (4) (1) (2) (3) (4) (1) (2) (3) (4) $\xi =-3$ $\xi =-2$ $\xi =7$ $\xi =5$ ${𝐧}_{\mathrm{𝟏}}$ .008 .034 .052 .039 .015 .036 .051 .040 .018 .034 .051 .043 .009 .035 .041 .052 ${𝐧}_{\mathrm{𝟐}}$ .030 .043 .049 .049 .032 .044 .050 .051 .032 .048 .047 .049 .030 .046 .048 .052 ${𝐧}_{\mathrm{𝟑}}$ .016 .038 .050 .050 .025 .043 .036 .050 .031 .040 .051 .049 .021 .035 .048 .049 ${𝐧}_{\mathrm{𝟒}}$ .033 .044 .047 .049 .038 .043 .049 .051 .041 .043 .051 .051 .031 .040 .051 .049 $({\mu }_1,\mathrm{\dots },{\mu }_{10})=(1,0,1,0,1,0,1,0,1,0)$ ${𝐧}_{\mathrm{𝟏}}$ .010 .033 .046 .037 .016 .039 .047 .044 .017 .038 .051 .041 .007 .035 .043 .052 ${𝐧}_{\mathrm{𝟐}}$ .030 .043 .048 .045 .041 .050 .052 .046 .037 .044 .047 .046 .032 .043 .047 .052 ${𝐧}_{\mathrm{𝟑}}$ .016 .037 .051 .053 .026 .042 .038 .051 .023 .042 .051 .053 .022 .034 .048 .048 ${𝐧}_{\mathrm{𝟒}}$ .031 .038 .047 .050 .045 .046 .049 .048 .038 .042 .051 .051 .032 .041 .051 .051 $({\mu }_1,\mathrm{\dots },{\mu }_{10})=(2,1,4,2,3,1,1,3,1,1)$ ${𝐧}_{\mathrm{𝟏}}$ .008 .039 .047 .041 .014 .037 .052 .042 .017 .040 .051 .039 .022 .033 .043 .049 ${𝐧}_{\mathrm{𝟐}}$ .032 .048 .050 .048 .047 .046 .049 .048 .044 .051 .047 .047 .041 .042 .050 .048 ${𝐧}_{\mathrm{𝟑}}$ .016 .035 .051 .050 .024 .040 .051 .051 .022 .041 .051 .051 .032 .038 .052 .048 ${𝐧}_{\mathrm{𝟒}}$ .031 .042 .051 .048 .037 .044 .051 .050 .031 .044 .051 .051 .033 .040 .051 .050

NOTE: ${𝐧}_{\mathrm{𝟏}}=(5,\mathrm{\dots },5)$; ${𝐧}_{\mathrm{𝟐}}=(15,\mathrm{\dots },15)$; ${𝐧}_{\mathrm{𝟑}}=(4,4,4,5,5,5,10,10,10,10)$; ${𝐧}_{\mathrm{𝟒}}=(4,4,4,12,12,12,15,15,15,15)$
1 – GVT; 2 – AJ MLRT; 3 – DK MLRT; 4 – PB test

Table 2. Powers of the tests for equality of quantiles at $\alpha =0.05$

$𝝈=(1,1,1)$
		$p=0.25$					$p=0.5$					$p=0.75$
$(n_1,\mathrm{\dots },n_k)$	$({\mu }_1,\mathrm{\dots },{\mu }_k)$	(1)	(2)	(3)	(4)		(1)	(2)	(3)	(4)		(1)	(2)	(3)	(4)
(5, 10, 11)	(1,1,1)	.031	.048	.051	.051		.037	.048	.049	.053		.029	.048	.048	.052
	(1.5,1,1)	.075	.089	.092	.131		.085	.099	.108	.105		.040	.095	.110	.074
	(2.5,1,1)	.359	.383	.357	.501		.516	.563	.557	.608		.328	.559	.662	.543
	(1,1,2.5)	.680	.757	.702	.744		.855	.865	.870	.859		.716	.754	.791	.752
	(1,3,1)	.899	.923	.893	.929		.977	.986	.986	.981		.932	.953	.973	.949
	(1,1,3)	.915	.944	.917	.949		.982	.989	.990	.987		.944	.950	.964	.950
(10,10,10)	(1,1,1)	.032	.046	.052	.048		.042	.046	.048	.052		.032	.049	.048	.047
	(1.5,1,1)	.101	.135	.142	.142		.134	.152	.174	.174		.103	.134	.156	.146
	(2.5,1,1)	.699	.744	.710	.791		.856	.888	.890	.891		.775	.838	.880	.853
	(1,1,2.5)	.720	.735	.717	.795		.875	.883	.896	.897		.800	.843	.879	.854
	(1,3,1)	.930	.934	.911	.956		.989	.990	.989	.992		.969	.979	.986	.984
	(1,1,3)	.938	.934	.907	.956		.989	.989	.992	.993		.967	.978	.989	.989
(10,15,20)	(1,1,1)	.034	.049	.053	.052		.040	.049	.052	.050		.038	.051	.052	.052
	(1.5,1,1)	.150	.148	.149	.194		.160	.181	.203	.203		.124	.161	.179	.150
	(2.5,1,1)	.750	.781	.749	.842		.923	.930	.924	.941		.890	.935	.950	.928
	(1,1,2.5)	.963	.972	.962	.975		.991	.993	.994	.991		.938	.971	.976	.977
	(1,3,1)	.994	.994	.991	.999		.999	1	1	1		.949	.999	1	1
	(1,1,3)	.999	.999	.999	1		1	1	1	1		.999	.999	1	1
$𝝈=(1,1,1,1,1)$
(5,5,5,10,10)	(1,1,1,1,1)	.029	.052	.050	.048		.033	.047	.052	.051		.025	.049	.047	.047
	(1.5,1,1.5,1,1)	.069	.094	.111	.116		.078	.111	.132	.112		.037	.101	.118	.076
	(2.5,1,1,1,1)	.242	.260	.244	.357		.379	.413	.386	.483		.208	.460	.540	.360
	(1,1,1,1,2.5)	.489	.591	.504	.548		.714	.789	.762	.730		.590	.692	.808	.629
	(1,3,1,1,1)	.392	.412	.368	.560		.604	.639	.573	.708		.395	.752	.808	.627
	(1,1,1,1,3)	.799	.836	.747	.835		.940	.969	.954	.954		.883	.938	.975	.910
(5,10,11,13,15)	(1,1,1,1,1)	.033	.048	.051	.051		.034	.047	.050	.050		.032	.051	.053	.048
	(1.5,1,1.5,1,1)	.111	.142	.152	.164		.163	.177	.180	.167		.083	.159	.174	.128
	(2.5,1,1,1,1)	.277	.286	.264	.463		.434	.432	.442	.575		.274	.532	.605	.440
	(1,1,1,1,2.5)	.851	.853	.780	.841		.954	.960	.959	.958		.906	.933	.962	.912
	(1,3,1,1,1)	.895	.875	.797	.915		.983	.986	.974	.983		.975	.991	.996	.979
	(1,1,1,1,3)	.982	.983	.961	.983		.999	.999	.999	.999		.995	.998	.999	.995
(5,5,7,15,15)	(1,1,1,1,1)	.030	.051	.048	.049		.037	.048	.047	.054		.034	.052	.052	.054
	(1.5,1,1.5,1,1)	.094	.107	.122	.146		.109	.136	.160	.142		.056	.127	.141	.084
	(2.5,1,1,1,1)	.259	.278	.258	.399		.402	.437	.395	.513		.255	.499	.566	.390
	(1,1,1,1,2.5)	.772	.746	.760	.772		.919	.951	.936	.961		.828	.884	.938	.845
	(1,3,1,1,1)	.425	.428	.432	.592		.631	.667	.610	.737		.503	.798	.837	.643
	(1,1,1,1,3)	.967	.980	.942	.966		.999	.999	.998	.997		.985	.994	.998	.986

NOTE: 1 – GVT; 2 – AJ MLRT; 3 – DK MLRT; 4 – PB test

Table 2 continued.
$𝝈=(1,\mathrm{\dots },1)$; $k=10$ $({\mu }_4,\mathrm{\dots },{\mu }_{10})=(1,\mathrm{\dots },1)$ $p=0.25$ $p=0.5$ $p=0.75$ $(n_1,\mathrm{\dots },n_k)$ $({\mu }_1,\mathrm{\dots },{\mu }_k)$ (1) (2) (3) (4) (1) (2) (3) (4) (1) (2) (3) (4) ${𝐧}_{\mathrm{𝟏}}$ $({\mu }_1,{\mu }_2,{\mu }_3)$ (1,1,1) .014 .034 .049 .047 .023 .040 .050 .046 .020 .040 .050 .044 (1.5,2,1) .047 .085 .119 .122 .107 .139 .175 .148 .047 .117 .188 .109 (2.5,2,1) .145 .191 .249 .261 .280 .337 .371 .375 .112 .345 .506 .282 (2.5,2,3) .327 .420 .521 .536 .600 .709 .761 .726 .321 .685 .858 .602 (3.5,2,3) .519 .591 .681 .743 .823 .899 .914 .910 .566 .920 .982 .841 (5,1,4) .927 .902 .901 .982 .998 .994 .996 .999 .989 .999 1 .998 ${𝐧}_{\mathrm{𝟐}}$ $({\mu }_1,{\mu }_2,{\mu }_3)$ (1,1,1) .028 .042 .051 .048 .028 .047 .050 .050 .032 .041 .050 .051 (1.5,2,1) .073 .082 .114 .160 .102 .128 .147 .171 .046 .120 .173 .090 (2.5,2,1) .157 .147 .233 .318 .239 .278 .320 .400 .097 .322 .459 .225 (2.5,2,3) .323 .289 .461 .573 .542 .578 .689 .719 .236 .729 .942 .538 (3.5,2,3) .472 .420 .594 .729 .735 .776 .853 .897 .461 .928 .968 .760 (5,1,4) .772 .719 .782 .989 .963 .969 .974 .998 .906 .999 1 .993

NOTE: ${𝐧}_{\mathrm{𝟏}}=(5,\mathrm{\dots },5)$; ${𝐧}_{\mathrm{𝟐}}=(4,4,4,7,7,7,8,9,10,11)$; NOTE: 1 – GVT; 2 – AJ MLRT; 3 – DK MLRT; 4 – PB test

4.3. Parametric Bootstrap Method

We shall now see the PB pairwise CIs developed in [malekzadeh2023simultaneous]. Let $T_i=\frac{{\widehat{\xi }}_{iu}-{\xi }_i}{\sqrt{b_i{S}_i}}$, where ${\widehat{\xi }}_{iu}$’s, $a_i$’s and $b_i$’s are as defined in the preceding Section 4.2. Let $Z_i\sim N(0,1)$ independently of $U_i^2={\chi }_{m_i}^2/m_i$, $i=1,..,k.$ Then $T_i$ is distributed as $T_i\stackrel{d}{=}\frac{Z_i/\sqrt{n_i}+z_p(a_i{U}_i-1)}{\sqrt{b_i}{U}_i}.$ Define

(24)

Let and let $t_{1-\alpha }^B$ denote the 100$(1-\alpha )$ percentile of $T$. Then

(25)

are $1-\alpha$ simultaneous CIs for ${\xi }_i-{\xi }_j$,

As in the preceding section, we can take $a_i=1$ and $b_i=1/n_i$ for $i=1,\mathrm{\dots },k.$ In this case ${\widehat{\xi }}_{iu}$ becomes ${\widehat{\xi }}_i={\overline{X}}_i+z_p{S}_i$ and the $T_i$ and $T_{ij}$ simplify to

respectively. Let $t_{1-\alpha }^{\ast }$ denote the 100$(1-\alpha )$ percentile of . Then

(26)

are $1-\alpha$ simultaneous CIs for ${\xi }_i-{\xi }_j$, .

4.4. Bonferroni Simultaneous Fiducial CIs

On the basis of a normal approximation to the noncentral $t$ distribution, in Chapter LABEL:ch2, we have developed a fiducial CI for the difference between two quantiles. This CI is not only simple, but also very satisfactory even for small samples. To describe their CI for ${\xi }_i-{\xi }_j$, let

(27) $$h(x;n_i,p_i)=\frac{c_{m_i}{z}_p+x\sqrt{\frac{c_{m_i}^2}{n_i}+\frac{1}{2m_i}\left(z_p^2-\frac{x^2}{n}\right)}}{c_{m_i}^2-\frac{x^2}{2m_i}}\mathrm{with}{c}_{m_i}=1+\frac{1}{4m_i}.$$

For , let

(28)

where $z_{\alpha }$ denotes the $\alpha$ quantile of the standard normal distribution. Furthermore, let

(29)

The 100$(1-2\alpha )$% approximate fiducial CI for ${\xi }_i-{\xi }_j$ is given by

(30) $$({\overline{x}}_i-{\overline{x}}_j+D_{ij;\alpha },{\overline{x}}_i-{\overline{x}}_j+D_{ij;1-\alpha }).$$

The simulation study in Chapter LABEL:ch2 indicated that the above CI is very satisfactory even for small sample sizes.

On the basis of the CI in (30), the 100$(1-\alpha )$% simultaneous Bonferroni CIs for all ${\xi }_i-{\xi }_j$, , can be expressed as

(31)

where ${\alpha }^{\ast }=\frac{\alpha }{k(k-1)}.$

4.5. Coverage and Precision Studies

To judge the properties of the pairwise CIs and to compare them, we carried out simulation studies as follows. To estimate the fiducial generalized CIs FGU and FG, we first generated 10,000 samples from different normal distributions, then we used simulation consisting of 10,000 runs to find simultaneous CIs based on each set of samples. The percentage of the 10,000 simultaneous CIs that include ${\xi }_i-{\xi }_j$ for all , is a Monte Carlo estimate of the coverage probability. To understand and to compare the precisions of the different methods, we also estimates expected volume^1/K, where $K$ is the number of simultaneous CIs. The coverage probabilities and precisions of the PB simultaneous CIs are estimated similarly. In our simulation study, we take the parameter values (without loss of generality) as $𝝁=(1,\mathrm{\dots },1)$ and $𝝈=({\sigma }_1,{\sigma }_2,\mathrm{\dots },{\sigma }_k)$ with ${\sigma }_1=1$ and ${\sigma }_j\le 1$, $j=2,\mathrm{\dots },k$. Furthermore, as argued in Chapter LABEL:ch2, we can take $p\in [0.5,1)$. This is because, if $(L_i,U_i)$ is a CI based on $({\overline{X}}_i,S_i^2,{\overline{X}}_j,S_j^2)$ for ${\xi }_{i,p}-{\xi }_{j,p}$, then $(-U_i,-L_i)$ is the CI based on $(-{\overline{X}}_i,S_i^2,-{\overline{X}}_j,S_j^2)$ for ${\xi }_{i,1-p}-{\xi }_{j,1-p}$.

In Table 3, we present the coverage probabilities and volumes^1/3 of FGU CIs based on unbiased estimators of ${\xi }_i$ and FG CIs based on ${\widehat{\xi }}_i={\overline{X}}_i+z_p{S}_i$. As expected, both estimated coverage probabilities and volumes^1/3 in Table 3 clearly indicate that the FGU and FG CIs are very similar. Minor differences between them could be due to simulation errors. In Table 4, we compare the Bonferroni (BF), fiducial CIs based on $\widehat{\xi }$ (FG) and the PB CIs. An examination of table values shows that the PB CIs are better than other CIs in terms of coverage probabilities and precisions. Their coverage probabilities are very close to the nominal level 0.90, and their averages of volumes^1/3 are smaller than those of the other CIs. The fiducial (FG) CIs and the PB CIs are practically the same for larger sample sizes. Bonferroni CIs are conservative having coverage probabilities larger than the nominal level for all the cases considered. The expected volumes of the BF CIs are somewhat larger than those of the PB abd FG CIs, but not much larger. Considering the simplicity of the BF CIs, they can be used in applications when sample sizes are 30 or more.

Coverage probabilities and averages of the volumes^1/6 of BF, FG and PB CIs for the case of $k=4$ are reported in Table 5. Note that, when $k=4$, there are six simultaneous CIs for the differences ${\xi }_i-{\xi }_j$, Performances and comparisons of these CIs are very similar as in the case of $k=3$. We once again see that the PB CIs are narrower than the other CIs for all cases.

Overall, we see that calculation of the PB CIs and the FG CIs involve simulation and they share similar computational difficulties. Nevertheless, the PB CIs are better than the FG CIs in terms coverage probability and precision. In fact, the PB CIs are preferable to other CIs for all cases. For large sample sizes or when the study involves a large number of groups, the BF CIs are straightforward to compute.

Table 3. Coverage probabilities and average of volume^1/3 of 90% simultaneous fiducial CIs

	$k=3$
	$𝐧=(10,10,10)$					$𝐧=(10,15,15)$
$𝝈$	$p=0.75$		$p=0.90$			$p=0.75$		$p=0.90$
	FGU CI	FG CI	FGU CI	FG CI		FGU CI	FG CI	FGU CI	FG CI
(1,1,1)	.931(2.31)	.930(2.31)	.946(2.95)	.947(2.95)		.926(1.99)	.928(1.99)	.934(2.46)	.937(2.51)
(1,.9,.8)	.931(2.11)	.933(2.11)	.947(2.64)	.950(2.64)		.921(1.79)	.922(1.79)	.937(2.29)	.940(2.26)
(1,.9,.6)	.931(1.96)	.931(1.96)	.939(2.46)	.941(2.47)		.926(1.71)	.926(1.71)	.926(2.13)	.929(2.14)
(1,.9,.1)	.899(1.73)	.900(1.74)	.898(2.17)	.899(2.22)		.897(1.51)	.898(1.52)	.899(1.92)	.900(1.95)
(1,.5,.1)	.904(1.34)	.904(1.34)	.900(1.68)	.900(1.71)		.902(1.18)	.903(1.19)	.905(1.49)	.904(1.53)
(1,.1,.1)	.918(0.83)	.929(0.83)	.924(1.05)	.927(1.06)		.900(0.73)	.899(0.73)	.918(0.94)	.918(0.94)
	$𝐧=(20,15,25)$					$𝐧=(30,35,30)$
(1,1,1)	.919(1.58)	.920(1.58)	.920(1.94)	.922(1.94)		.909(1.18)	.909(1.18)	.919(1.48)	.920(1.48)
(1,.9,.8)	.918(1.41)	.918(1.41)	.925(1.78)	.927(1.78)		.910(1.07)	.911(1.07)	.914(1.32)	.915(1.32)
(1,.9,.6)	.911(1.35)	.913(1.35)	.927(1.69)	.929(1.70)		.913(1.00)	.914(1.00)	.917(1.22)	.918(1.22)
(1,.9,.1)	.904(1.23)	.905(1.24)	.899(1.53)	.899(1.54)		.904(0.87)	.904(0.87)	.900(1.07)	.901(1.08)
(1,.5,.1)	.900(0.94)	.900(0.95)	.895(1.15)	.895(1.16)		.903(0.68)	.903(0.68)	.905(0.83)	.905(0.84)
(1,.1,.1)	.908(0.55)	.908(0.55)	.912(0.68)	.912(0.69)		.906(0.42)	.907(0.42)	.913(0.53)	.914(0.53)

Table 4. Coverage probabilities and average of volume^1/3 of 90% simultaneous CIs

	$k=3$
	$𝐧=(10,10,10)$
$𝝈$	$p=0.50$				$p=0.75$				$p=0.90$
	BF-CIs	FG CIs	PB CIs		BF-CIs	FG CIs	PB CIs		BF-CIs	FG CIs	PB CIs
(1,1,1)	.926(2.08)	.924(2.06)	.897(1.93)		.928(2.44)	.930(2.32)	.899(2.13)		.929(2.44)	.953(2.92)	.906(2.57)
(1,.9,.8)	.926(1.88)	.915(1.83)	.899(1.74)		.928(2.20)	.935(2.12)	.903(1.92)		.928(2.20)	.958(2.70)	.911(2.31)
(1,.9,.6)	.924(1.75)	.919(1.72)	.901(1.62)		.927(2.05)	.925(1.93)	.901(1.80)		.929(2.05)	.945(2.48)	.897(2.18)
(1,.9,.1)	.910(1.55)	.901(1.50)	.894(1.47)		.918(1.80)	.904(1.74)	.895(1.71)		.925(1.80)	.905(2.27)	.893(2.17)
(1,.5,.1)	.911(1.20)	.900(1.17)	.885(1.13)		.920(1.39)	.908(1.35)	.896(1.32)		.926(1.39)	.903(1.71)	.891(1.67)
(1,.1,.1)	.927(0.76)	.914(0.73)	.896(0.69)		.930(0.87)	.921(0.83)	.897(0.78)		.933(0.87)	.927(1.06)	.905(0.98)
	$𝐧=(10,15,15)$
(1,1,1)	.922(1.79)	.918(1.76)	.904(1.69)		.923(2.07)	.919(1.97)	.905(1.87)		.926(2.07)	.944(2.51)	.912(2.27)
(1,.9,.8)	.922(1.63)	.912(1.59)	.906(1.54)		.923(1.89)	.924(1.81)	.905(1.71)		.926(1.89)	.936(2.26)	.905(2.08)
(1,.9,.6)	.920(1.54)	.911(1.50)	.897(1.45)		.924(1.77)	.920(1.70)	.903(1.62)		.926(1.77)	.932(2.14)	.905(1.98)
(1,.9,.1)	.911(1.39)	.901(1.35)	.903(1.34)		.917(1.59)	.896(1.54)	.899(1.52)		.925(1.59)	.905(1.95)	.905(1.95)
(1,.5,.1)	.915(1.09)	.901(1.04)	.896(1.05)		.921(1.25)	.906(1.19)	.897(1.19)		.927(1.25)	.905(1.53)	.900(1.53)
(1,.1,.1)	.926(0.70)	.898(0.65)	.904(0.65)		.930(0.80)	.919(0.75)	.901(0.73)		.934(0.80)	.923(0.97)	.904(0.91)
	$𝐧=(20,15,25)$
(1,1,1)	.920(1.43)	.908(1.39)	.896(1.36)		.922(1.63)	.915(1.56)	.899(1.50)		.923(1.63)	.924(1.94)	.898(1.82)
(1,.9,.8)	.920(1.30)	.917(1.26)	.899(1.23)		.921(1.47)	.914(1.42)	.903(1.37)		.921(1.47)	.925(1.78)	.909(1.66)
(1,.9,.6)	.920(1.23)	.904(1.19)	.904(1.17)		.921(1.40)	.908(1.34)	.900(1.31)		.923(1.40)	.920(1.66)	.904(1.60)
(1,.9,.1)	.913(1.13)	.905(1.10)	.905(1.10)		.919(1.28)	.905(1.24)	.902(1.23)		.923(1.28)	.900(1.55)	.901(1.55)
(1,.5,.1)	.915(0.86)	.896(0.83)	.903(0.83)		.919(0.98)	.904(0.93)	.901(0.93)		.923(0.98)	.899(1.18)	.904(1.18)
(1,.1,.1)	.928(0.52)	.900(0.48)	.895(0.48)		.930(0.58)	.911(0.55)	.902(0.54)		.931(0.58)	.905(0.68)	.899(0.66)
	$𝐧=(30,35,30)$
(1,1,1)	.919(1.10)	.907(1.06)	.902(1.04)		.918(1.23)	.913(1.18)	.899(1.16)		.920(1.16)	.918(1.47)	.904(1.42)
(1,.9,.8)	.919(0.99)	.898(0.94)	.899(0.93)		.919(1.12)	.910(1.07)	.895(1.05)		.918(1.05)	.920(1.32)	.904(1.27)
(1,.9,.6)	.918(0.93)	.895(0.88)	.901(0.87)		.920(1.05)	.907(1.00)	.900(0.98)		.922(0.98)	.912(1.23)	.903(1.19)
(1,.9,.1)	.918(0.84)	.898(0.77)	.898(0.77)		.919(0.95)	.891(0.86)	.905(0.88)		.923(0.88)	.898(1.07)	.906(1.08)
(1,.5,.1)	.919(0.65)	.905(0.61)	.897(0.60)		.922(0.73)	.898(0.68)	.905(0.68)		.924(0.68)	.909(0.84)	.901(0.83)
(1,.1,.1)	.928(0.41)	.907(0.38)	.895(0.37)		.929(0.45)	.910(0.43)	.903(0.42)		.931(0.42)	.913(0.53)	.902(0.52)

Table 5. Coverage probabilities and average of volume^1/6 of 90% simultaneous CIs

	$k=4$
	$𝐧=(10,10,10,10)$
$𝝈$	$p=0.50$				$p=0.75$				$p=0.90$
	BF-CIs	FG CIs	PB CIs		BF-CIs	FG CIs	PB CIs		BF-CIs	FG CIs	PB CIs
(1,1,1,1)	.941(2.44)	.917(2.22)	.896(2.18)		.943(2.92)	.935(2.58)	.910(2.42)		.947(3.91)	.952(3.27)	.917(2.92)
(1,.9,.8,.9)	.941(2.20)	.912(2.03)	.890(1.97)		.944(2.63)	.934(2.32)	.909(2.19)		.947(3.52)	.949(2.95)	.915(2.63)
(1,.9,.6,.4)	.938(1.79)	.910(1.70)	.895(1.60)		.942(2.13)	.926(1.94)	.896(1.81)		.947(2.83)	.943(2.52)	.902(2.22)
(1,.9,.4,.1)	.930(1.53)	.903(1.48)	.898(1.39)		.938(1.80)	.911(1.73)	.905(1.62)		.945(2.36)	.910(2.25)	.888(2.05)
(1,.5,.1,.1)	.938(1.25)	.906(0.93)	.898(0.88)		.945(1.46)	.919(1.07)	.906(1.02)		.950(1.92)	.919(1.40)	.904(1.28)
(1,.1,.1,.1)	.943(0.64)	.914(0.59)	.899(0.56)		.947(0.75)	.928(0.68)	.906(0.64)		.950(0.99)	.944(0.87)	.915(0.79)
	$𝐧=(10,15,10,15)$
(1,1,1,1)	.936(2.16)	.905(2.01)	.899(1.99)		.940(2.56)	.923(2.28)	.905(2.20)		.943(3.37)	.944(2.89)	.902(2.61)
(1,.9,.8,.9)	.936(1.95)	.904(1.79)	.903(1.80)		.939(2.30)	.921(2.04)	.909(1.98)		.944(3.03)	.939(2.60)	.906(2.36)
(1,.9,.6,.4)	.934(1.61)	.915(1.53)	.902(1.50)		.937(1.89)	.928(1.75)	.906(1.67)		.943(2.49)	.927(2.21)	.892(2.02)
(1,.9,.4,.1)	.928(1.39)	.905(1.33)	.903(1.32)		.935(1.62)	.906(1.51)	.903(1.50)		.943(2.11)	.912(1.97)	.894(1.86)
(1,.5,.1,.1)	.937(1.11)	.905(0.83)	.909(0.82)		.941(1.28)	.907(0.94)	.907(0.92)		.947(1.66)	.917(1.22)	.901(1.15)
(1,.1,.1,.1)	.941(0.59)	.904(0.53)	.911(0.53)		.945(0.69)	.926(0.62)	.907(0.59)		.949(0.90)	.932(0.79)	.904(0.73)
	$𝐧=(20,15,20,25)$
(1,1,1,1)	.930(1.62)	.902(1.54)	.907(1.54)		.933(1.86)	.915(1.75)	.907(2.60)		.934(2.38)	.922(2.16)	.895(2.00)
(1,.9,.8,.9)	.930(1.47)	.901(1.38)	.905(1.38)		.931(1.68)	.916(1.56)	.901(2.34)		.935(2.14)	.926(1.93)	.898(1.81)
(1,.9,.6,.4)	.930(1.23)	.907(1.20)	.908(1.16)		.932(1.41)	.916(1.37)	.892(2.03)		.935(1.79)	.918(1.69)	.903(1.55)
(1,.9,.4,.1)	.929(1.08)	.903(1.05)	.902(1.01)		.933(1.22)	.904(1.17)	.901(1.89)		.938(1.55)	.905(1.51)	.898(1.40)
(1,.5,.1,.1)	.939(0.87)	.907(0.65)	.906(0.64)		.941(0.99)	.903(0.73)	.902(1.16)		.944(1.25)	.911(0.92)	.898(0.87)
(1,.1,.1,.1)	.940(0.43)	.915(0.40)	.904(0.39)		.942(0.49)	.915(0.44)	.903(0.73)		.944(0.61)	.918(0.55)	.898(0.52)
	$𝐧=(30,35,25,30)$
(1,1,1,1)	.928(1.29)	.895(1.22)	.892(1.21)		.930(1.46)	.891(1.33)	.898( 1.33)		.929(1.83)	.907(1.64)	.897(1.62)
(1,.9,.8,.9)	.927(1.15)	.897(1.07)	.892(1.08)		.928(1.31)	.891(1.19)	.901( 1.19)		.928(1.64)	.907(1.48)	.904(1.46)
(1,.9,.6,.4)	.928(0.94)	.905(0.90)	.898(0.87)		.929(1.06)	.909(1.01)	.902( 0.98)		.932(1.32)	.918(1.24)	.905(1.20)
(1,.9,.4,.1)	.928(0.80)	.902(0.78)	.892(0.74)		.933(0.90)	.904(0.87)	.900( 0.84)		.936(1.12)	.904(1.07)	.904(1.04)
(1,.5,.1,.1)	.937(0.64)	.901(0.48)	.896(0.47)		.939(0.72)	.912(0.54)	.895( 0.52)		.942(0.90)	.912(0.67)	.906(0.65)
(1,.1,.1,.1)	.938(0.34)	.908(0.31)	.896(0.31)		.939(0.38)	.909(0.35)	.900( 0.34)		.941(0.48)	.907(0.43)	.895(0.42)

5. Examples

Example 3.5.1. This example and the associated data were taken from [li2012comparison]. The hypothesis that Vitamin D protects against colon cancer emerged from a study by [garland1980sunlight]. The result of the study is appeared to support the hypothesis that there is a relationship between Vitamin D and Colorectal Cancer (CRC). However, the effects of Vitamin D supplementation on incidence and mortality of CRC remain inconclusive. To investigate further, a Vitamin D study was conducted in Roswell Park Cancer Center where CRC patients were given a 6-month treatment with vitamin D supplements. The purpose of the study was to find if the vitamin D supplement treatment could sufficiently increase the serum 1, 25-D3 and 24, 25-D3 levels at the end of the study period. Subjects were divided into three groups according to the baseline serum 25-D3 level of each subject, namely, (i) vitamin D3 deficient if serum 25-D3 level less than 20 ng/ml, (ii) insufficient if it was $\ge 20$ and , and (iii) sufficient if $\ge$ 32. Tests for normality of the data by [li2012comparison] indicated that the data fit normal distributions. The summary statistics of serum 1, 25-D3 and 24, 25-D3 vitamin D3 metabolites Li et al. (2012) are reproduced here in Table 6.

Table 6. Descriptive statistics for the data from a vitamin D study conducted in Roswell Park Cancer Institute

Group	Variable	Size ($n_i$)	Mean (${\overline{x}}_i$)	SD ($s_i$)
Vitamin D3 sufficient	1, 25-D3	16	62.39	17.99
	24, 25-D3	17	4.65	1.98
Vitamin D3 insufficient	1, 25-D3	22	72.60	23.52
	24, 25-D3	22	3.62	1.17
Vitamin D3 deficient	1, 25-D3	9	70.13	19.67
	24, 25-D3	9	2.66	1.35

A hypothesis of interest here is that if all of these three vitamin D groups would reach the same serum 1, 25-D3 and 24, 25-D3 vitamin levels. [li2012comparison] have noted that a thorough understanding about how the distributions differ across three groups is to compare the quantiles, especially the common quantiles such as 1st quartile, median, and third quartile.

Table 7. The p-values for testing the equality of percentiles for the vitamin D study

	serum level 1,25-D3					serum level 24,25-D3
$p$	GVT	PB	AJ-MLRT	DK-MLRT		GVT	PB	AJ-MLRT	DK-MLRT
.05	.973	.942	.851	.898		.297	.300	.276	.260
.10	.914	.886	.881	.844		.251	.235	.211	.174
.25	.653(.643)	.642	.689	.615		.123(.123)	.103	.092	.067
.50	.317(.306)	.316	.320	.316		.032(.034)	.027	.028	.025
.75	.197(.193)	.198	.196	.212		.019(.024)	.015	.008	.005
.90	.193	.194	.206	.208		.028	.019	.005	.003
.95	.211	.205	.227	.218		.031	.022	.005	.003

Note: The values in parentheses are the p-values of the GVT given in Table 7 of Li et al. (2012)

We estimated the p-values of the generalized variable test (GVT), the PB test and the MLRTs and reported them in Table 7. In serum level 1, 25-D3, no significant differences among the groups in terms of quartiles and medians. All the tests produced p-values larger than commonly used practical nominal levels. However, all the tests indicate that significant differences exist in medians (or means) and 3rd quartiles among groups for 24,25-D3. The DK-MLRT produced p-values that are smaller than the corresponding AJ-MLRT for testing all percentiles considered in the table. We also see that the GVT produced a larger p-value of 0.123 testing the equality of the first quartiles, because this test is conservative for most cases. All the tests, as indicated by our earlier simulation studies, produced similar p-values for testing the equality of group means ($p=0.5$) in both serum levels. In particular, all tests indicate that there is no significant difference among group means in 1,25-D3 and they provide evidence to conclude that the group means are significantly different in 24,25-D3 level.

We shall now compute various 95% simultaneous CIs for the differences ${\xi }_{1,.7}-{\xi }_{2,.7}$, ${\xi }_{1,.7}-{\xi }_{3,.7}$ and ${\xi }_{2,.7}-{\xi }_{3,.7}$. We chose 70th percentiles so that we can compare our results with those given in Table 9 of [malekzadeh2023simultaneous]. These simultaneous CIs are given in Table 8. We first observe that three simultaneous CIs by all methods include zero, and so percentiles from these three groups are not significantly different. Furthermore, we notice that the FG and FGU simultaneous CIs are very similar with same volume^1/3, and the PB and PBU CIs are in good agreement with practically the same volume^1/3.

Table 8. 95% simultaneous pairwise CIs

	Serum Vitamin Level 1,25-D3
Difference	FGU	FG	PBU	PB	BF
${\xi }_{1,.7}-{\xi }_{2,.7}$	(-32.38, 6.19)	(-32.44, 6.23)	(-30.89, 4.69 )	(-30.88, 4.66 )	(-31.99, 5.86 )
${\xi }_{1,.7}-{\xi }_{3,.7}$	(-31.64, 14.06)	(-31.44, 14.20)	(-29.86, 12.29)	(-29.83, 12.31)	(-36.28, 13.02)
${\xi }_{2,.7}-{\xi }_{3,.7}$	(-19.39, 28.01)	(-19.20, 28.18)	(-17.55, 26.17)	(-17.48, 26.16)	(-23.67, 26.70)
critical value	2.67	2.87	2.47	2.64	—–
volume1/3	43.72	43.72	40.32	40.18	45.47
	Serum Vitamin Level 24,25-D3
${\xi }_{1,.7}-{\xi }_{2,.7}$	(-0.09, 3.02)	(-0.11, 3.02)	(0.02, 2.91)	(0.01, 2.90)	(0.07, 3.12)
${\xi }_{1,.7}-{\xi }_{3,.7}$	(0.42, 4.21)	(0.42, 4.22)	(0.55, 4.08)	(0.57, 4.08 )	(0.23, 4.22)
${\xi }_{2,.7}-{\xi }_{3,.7}$	(-0.63, 2.33)	(-0.62, 2.35)	(-0.53, 2.23)	(-0.51, 2.24)	(-0.98, 2.20)
critical value	2.68	2.88	2.49	2.67	—–
volume1/3	3.27	3.28	3.04	3.03	3.38

Example 3.5.2. Operating room anesthesia involves clinical and managerial decision making that relies on communication over periods of less than 5 minutes. [ledolter2011analysis] have noted that the latency data including times for anesthesia providers to respond to messages well described by lognormal models. These authors have used the generalized variable approach to compare several lognormal means based on data consisting of 472 messages from four groups.

Group 1:	no prior message and the message was not anchored,
Group 2:	no prior message and the message was anchored,
Group 3:	with prior message and the message was not anchored, and
Group 4:	with prior message and the message was anchored.

The summary statistics reported in Table 8 of [li2012comparison] are reproduced here in Table 9. In Table 10, we reported the p-values of the tests for the equality of 100$p$ percentiles for several values of $p$ ranging from 0.05 to 0.95. All the tests indicate that differences exist among percentiles for $p\le 0.30$. As the GVT is conservative, it produced a little larger p-value than those of the other tests for testing the equality of medians.

Table 9. Descriptive statistics for the log time to acknowledge in [ledolter2011analysis]

Group	Size	Mean	$Q_1$	Median	$Q_3$	SD	Min	Max
1. No prior message, not anchored	245	0.117	-0.288	0.122	0.525	0.580	-1.715	1.564
2. No prior message, anchored	125	0.111	-0.248	0.068	0.470	0.607	-1.833	1.575
3. Prior message, not anchored	65	-0.019	-0.446	0.039	0.419	0.627	-1.427	1.515
4. Prior message, anchored	37	-0.112	-0.511	0.058	0.270	0.788	-2.040	1.358

Table 10. The estimated p-values for testing the equality of $p$th quantiles of the log time to acknowledge

$p$	GVT	PB	AJ-MLRT	DK-MLRT	Welch
.05	.030	.031	.008	.008
.10	.028	.030	.010	.010
.25	.042(.045)	.044	.027	.026	—
.50	.182(.186)	.177	.180	.177	.177
.75	.683(.680)	.659	.653	.648	—
.90	.876	.875	.834	.872
.95	.849	.860	.846	.873

Note: The values in parentheses are the p-values of the GVT given in Table 7 of Li et al. (2012)

We also constructed 95% pairwise CIs for all six differences of the 10th percentiles and reported them in Table 11. Since all the tests for equality indicated that the 10th percentiles are significantly different (see Table 10), we see that pairwise CIs due to all the methods indicate that the 10th percentiles of the groups 1 and 4 are significantly different at the level 0.05. The BF CIs also indicate that the groups 2 and 4 are somewhat different while other CIs indicate that they are not significantly different.

Table 11. 95% simultaneous pairwise CIs for the difference ${\xi }_{i,.10}-{\xi }_{j,.10}$

Difference	FGU	FG	PBU	PB	BF
${\xi }_{1,.10}-{\xi }_{2,.10}$	(-0.19, 0.27)	(-0.19, 0.27)	(-0.19, 0.27)	(-0.19, 0.27)	(-0.19, 0.29)
${\xi }_{1,.10}-{\xi }_{3,.10}$	(-0.11, 0.50	(-0.11, 0.50)	(-0.11, 0.50)	(-0.11, 0.50)	(-0.09, 0.55)
${\xi }_{1,.10}-{\xi }_{4,.10}$	( 0.02, 0.98)	( 0.01, 0.98)	( 0.02, 0.97)	( 0.02, 0.98)	( 0.07, 1.11)
${\xi }_{2,.10}-{\xi }_{3,.10}$	(-0.18, 0.49)	(-0.18, 0.49)	(-0.18, 0.49)	(-0.18, 0.49)	(-0.17, 0.53)
${\xi }_{2,.10}-{\xi }_{4,.10}$	(-0.04, 0.95)	(-0.04, 0.95)	(-0.04, 0.95)	(-0.04, 0.96)	( 0.00, 1.07)
${\xi }_{3,.10}-{\xi }_{4,.10}$	(-0.23, 0.84)	(-0.23, 0.83)	(-0.24, 0.84)	(-0.23, 0.84)	(-0.22, 0.94)
critical value	2.60	3.53	2.61	3.54
volume1/6	0.67	0.67	0.76	0.76	0.80