Received 10 February 2019, Accepted 9 December 2019, Available Online 23 June 2020.
1. INTRODUCTION
The problem of testing equality of several distributions within a family of distributions arises in many practical situations. If independent samples are available from different sources or collected from different life testing experiments, a homogeneity test may be used to check if the samples can be modeled by a single distribution. For example, in a placebo controlled randomized clinical trial of gamma interferon in chronic granulomatous disease (CGD), the data were collected on elapsed time (in days) from randomization to diagnosis of a serious infection [1]. The samples were taken from three hospitals, and the Kolmogorov–Smirnov test indicates that all three samples satisfy the model assumption of two-parameter exponential distribution. If a homogeneity test indicates that all three samples can be modeled by a single exponential distribution, then better inferences can be made on the basis of the pooled sample. For two-parameter exponential distributions, the problem of testing the equality of parameters of several distributions was considered by [2]. These authors have proposed an exact iterated test procedures for testing the equality of the location parameters under some conditions and for testing the equality of the scale parameters. Abebe [3] have addressed the problem of comparing several location parameters with more than one control. Li et al. [4] have proposed parametric bootstrap simultaneous confidence intervals for pairs of means. Recently, Krishna and Goel [5] have proposed some Bayesian solutions for estimating location parameters when the data are randomly censored.
To describe the problems that will be addressed in this paper and other related problems, we first note that the two-parameter exponential distribution has the probability density function (pdf) given by
f(x|a,b)=1bexp−x−ab,x>a>0,b>0,(1)
where
a is the location parameter (also called threshold parameter) and
b is the scale parameter. If it is assumed that
a=0, then the resulting distribution is referred to as the one-parameter exponential distribution with mean
b. Let us denote the two-parameter exponential distribution by
exp(a,b).
In life testing, the experiments are often discontinued as soon as a fixed number of items fail, yielding a type II censored sample. Let Xi1,…,Xiri be a type II censored sample from a life test of ni items whose lifetimes follow an exp (ai,bi) distribution, i=1,…,k. The problems that we consider here are (i) testing the equality of the scale parameters bi's, (ii) testing the equality of the location parameters ai's and (iii) the homogeneity test, that is, testing (a1,b1)=…=(ak,bk). The LRTs that are proposed in the sequel are applicable to the case where ri=ni for some or all i=1,…,k. In the following section we provide some preliminary results that are needed to derive the likelihood ratio tests (LRTs) for all three problems. In Section 3, we address the problem of testing the equality of scale parameters, derive the LRT and describe another two tests available in the literature. The problem of testing equality of the location parameters is addressed in Section 4, and the homogeneity test is described in Section 5. We also show that the null distributions of all three LRT statistics do not depend on any unknown parameters, and the percentiles or p-values of the LRTs can be estimated using Monte Carlo simulation. So the proposed LRTs are exact except for the simulation errors. For all three problems, we also propose convenient closed-form approximations to the null distributions, and they are very satisfactory even for small samples. Power studies and comparison of tests are given in Section 6. The proposed LRTs and other methods are illustrated using an example with real data in Section 7, and some concluding remarks are given in Section 8.
2. PRELIMINARIES
Let X(1)<…<X(r) be a type II censored sample obtained from a life test on n items whose lifetime have an exp (a,b) distribution, where a is the threshold parameter and b is the scale parameter. Let X¯r denote the mean of the uncensored lifetimes in the sample.
The MLEs (see [6]) based on the censored sample are given by
a^=X(1)andb^=rX¯r−X(1)+(n−r)X(r)−X(1)r. (2)
Furthermore, the MLEs are independent with
a^∼b2nχ22+a and b^∼b2rχ2r−22.(3)
For the uncensored case, that is, when r=n, the MLEs are given by
a^=X(1)andb^=X¯n−X(1),(4)
and
a^∼b2nχ22+aandb^∼b2nχ2n−22.
Log-likelihood function
Let a=(a1,…,ak) and b=(b1,…,bk). The log-likelihood function based on the k independent censored samples can be expressed as
l(a,b)=∑i=1k−rilnbi−∑j=1riXi(j)−aibi−(ni−ri)Xi(ri)−aibi+C,
where
C is a constant term. Following
(2), the MLEs of
ai and
bi can be expressed as
a^i=Xi1andb^i=riX¯ri−Xi1+ni−riXiri−Xi1∕ri,i=1,…,k,(5)
and
l(a,b) at
(a,b)=(a^,b^) can be simplified as
l(a^,b^)=−∑i=1kriln(b^i)−R+C,(6)
where
R=∑i=1kri.
3. TESTS FOR THE EQUALITY OF SCALE PARAMETERS
To test the equality of scale parameters, consider the hypotheses
H0:b1=…=bkvs.Ha:bi≠bjfor some i≠j.(7)
In the following, we first describe the LRT followed by an iterative test by [2], and a test based on union-intersection principle proposed in [7].
3.1. The Likelihood Ratio Test
Let b denote the common unknown parameter under H0 in (7). Then the log-likelihood function under this null hypothesis can be expressed as
l(a,b)=∑i=1k−rilnb−∑j=1riXi(j)−aib−(ni−ri)Xi(ri)−aib+C,
and is maximized at
a^ic=Xi(1)andb^c=1R∑i=1kriX¯ri−Xi(1)+(ni−ri)Xi(ri)−Xi(1)=∑i=1krib^iR,(8)
where
R=∑i=1kri and
b^i's are given in
(5). Thus, letting
Λb to denote the
−2ln(LRT statistic), it can be shown that
Λb=−2∑i=1krilnb^i∑i=1krib^iR, (9)
where
b^i is given in
(5). We note that the distribution of
Λb under
H0:b1=…=bk does not depend on any parameter, and
Λb∼−2∑i=1krilnWi∑i=1kriWiR,(10)
where
Wi's are independent
χ2ri−22∕(2ri) random variables. Notice also that the null distribution does not depend on the sample sizes and they depend only on
(r1,…,rk). For a given
(r1,…,rk), the exact percentiles of
Λb can be estimated using Monte Carlo simulation as shown in the following algorithm.
Algorithm 1
For a given (r1,…,rk) and (b^1,…,b^k),
Compute Λb=−2∑i=1krilnb^i∑i=1krib^iR,
Generate Wi=.5χ2ri−22∕ri, i=1,…,k,
Compute Λb∗=−2∑i=1krilnWi∑i=1kriWiR,
Repeat steps 2 and 3, for a large number of times, say, 100,000,
The proportion of Λb∗ that are greater than Λb is a Monte Carlo estimate of the p-value,
The 100 (1−α) percentile of these 100,000 Λb∗'s is a Monte Carlo estimate of the 100 (1−α) percentile of Λb.
The test statistic Λb is very similar to the one for testing b1=…=bk of several one-parameter exponential distributions, except that b^i in our present problem is distributed as χ2ri−22∕(2ri) whereas in the one-parameter case it is distributed like χ2ri2∕(2ri). This fact suggest that, for large samples, the null distribution of Λb can be approximated by χk−12 distribution. For large samples, the null hypothesis of equal scale parameters is rejected if Λb>χk−1;1−α2. A better approximation to the null distribution of Λb can be obtained by the moment matching method. Specifically, we determine the value of Qb so that E(Λb)=E(Qbχk−12) so that
ΛbQb∼χk−12,approximately.
Using the result that E(lnχm2)=ψ(m∕2)+ln(2), where ψ is the digamma function, and the expression (10), we see that
EΛb=2E∑i=1kriln∑i=1kriWiR−∑i=1kriln(Wi) =2∑i=1kriψ(R−k)−ψ(ri−1)+∑i=1kriln(ri∕R).(11)
Thus, Qb=EΛb∕(k−1), and the null hypothesis of equal shape parameters is rejected if Λb>Qbχk−1;1−α2.
To judge the accuracy of the above approximation, we estimated the upper percentiles of Λb in (10) using simulation with 1,000,000 runs and the percentiles of Qbχk−12 for k=3 and 5, and some small to large values of (r1,…,rk). Recall that the distribution of Λb does not depend on sample sizes, and so percentiles were estimated only for some values of (r1,…,rk). These percentiles are reported in Table 1. Comparison of Monte Carlo estimates of the percentiles and those of Qbχk−12 shows that the latter approximate percentiles are very close to those based on simulation even when (r1,…,r5)=(3,3,3,2,2). Thus, the improved approximation can be safely used for any sample sizes with two or more uncensored data. The usual χk−12 approximation may be used if ri≥100 for all i; see the last row of Table 1.
k=3
|
Percentiles
|
k=5
|
Percentiles
|
| (r1,r2,r3) |
99 |
95 |
90 |
(r1,…,r5) |
99 |
95 |
90 |
| (2, 2, 2) |
21.53 (22.18) |
14.23 (14.43) |
11.03 (11.09) |
(2, …, 2) |
30.73 (31.46) |
22.21 (22.48) |
18.29 (18.43) |
| (3, 3, 3) |
15.12 (15.31) |
9.92 (9.96) |
7.63 (7.66) |
(3, 3, 3, 2, 2) |
26.69 (26.71) |
19.07 (19.08) |
15.64 (15.65) |
| (4, 3, 4) |
14.08 (14.10) |
9.17 (9.17) |
7.05 (7.05) |
(4, 3, 4, 3, 4) |
20.30 (20.30) |
14.51 (14.51) |
11.91 (11.90) |
| (3, 2, 6) |
18.90 (19.13) |
12.40 (12.45) |
9.55 (9.57) |
(3, 2, 4, 3, 2) |
26.63 (26.64) |
19.06 (19.04) |
15.65 (15.61) |
| (5, 4, 3) |
13.97 (13.97) |
9.08 (9.09) |
6.99 (6.99) |
(5, 4, 3, 6, 2) |
23.45 (23.07) |
16.61 (16.48) |
13.57 (13.52) |
| (8, 7, 10) |
10.90 (10.90) |
7.10 (7.09) |
5.45 (5.45) |
(8, 7, 6, 8, 20) |
16.06 (16.08) |
11.48 (11.49) |
9.41 (9.42) |
| (10, 10, 10) |
10.46 (10.49) |
6.82 (6.82) |
5.24 (5.24) |
(4, 3, 2, 5, 6) |
23.41 (23.07) |
16.62 (16.48) |
13.59 (13.52) |
| (20, 20, 20) |
9.82 (9.81) |
6.37 (6.38) |
4.90 (4.90) |
(20, …, 20) |
14.14 (14.12) |
10.09 (10.09) |
8.27 (8.28) |
| (30, 30, 30) |
9.62 (9.60) |
6.25 (6.25) |
4.80 (4.80) |
(30, …, 30) |
13.85 (13.83) |
9.86 (9.88) |
8.10 (8.10) |
| (50, 50, 50) |
9.49 (9.44) |
6.14 (6.14) |
4.72 (4.72) |
(50, …, 50) |
13.62 (13.60) |
9.74 (9.72) |
7.97 (7.97) |
| (100, 100, 100) |
9.33 (9.32) |
6.07 (6.07) |
4.66 (4.66) |
(100, …, 100) |
13.43 (13.44) |
9.61 (9.60) |
7.88 (7.87) |
| χk−12 percentiles |
9.21 |
5.99 |
4.61 |
χk−12 percentiles |
13.28 |
9.49 |
7.78 |
Table 1Exact percentiles of Λb based on simulation and the percentiles of Qbχk−12.
3.2. An Iterative Test
We shall now describe a test for equality of scale parameters proposed by [2], which is based on an iterative procedure. Consider a sequence of nested hypotheses
H2:b1=b2,H3:b1=b2=b3,…,Hk:b1=…=bk.
A test for Hi is made if and only if Hi−1 is accepted for i=3,…,k. The null hypothesis (7) of equality scale parameters is accepted if and only if H2,…,Hk are accepted.
The test statistic for Hi:b1=…=bi is given by
Fi=rib^i∕(2ri−2)∑ji−1rjb^j∕S(i−1),(12)
where
S(i−1)=∑j=1i−1(2rj−2),
i=2,…,k. For a given level
α, choose
0<ηi<α so that
α=1−∏i=2k(1−ηi). The iterative test rejects
H0:b1=…=bk if
Fi<F2ri−2,S(i−1);ηi∕2orFi>F2ri−2,S(i−1);1−ηi∕2for any i∈{2,…,k},
where
Fm,n;q denotes the 100
q percentile of an
F distribution with the numerator df
m and the denominator df
n. If we choose
η2=…=ηk, then
ηi=1−(1−α)1∕(k−1),
i=2,…,k. Notice that the individual test for
Hi has size
ηi, but the overall test for
H0:b1=…=bk has size
α.
3.3. Union-Intersection Test
Shah and Rathod [7] have proposed the following test based on the union-intersection principle. The test rejects the null hypothesis of equality of scale parameters if
U=max{b^1,…,b^k}min{b^1,…,b^k}>u1−α,(13)
where
b^i is the MLE defined in
(5) and
u1−α is the 100
(1−α) percentile of
U.
Shah and Rathod have derived analytical expressions to the distribution function of U for k≤3 and r1=r2=r3. Even for these special cases, calculation of the critical value u1−α is numerically involved. As the null distribution of the test statistic U does not depend on any parameters, the critical value u1−α can be easily estimated using Monte Carlo simulation for any k≥2 and (r1,…,rk).
4. THE LRT FOR LOCATION PARAMETERS
Consider testing
H0:a1=…=akvs.Ha:ai≠ajfor some i≠j.(14)
Let a denote the common unknown parameter under the H0. Then the log-likelihood function under H0 in (14) can be expressed as
l(a,b)=∑i=1k−rilnbi−∑j=1riXi(j)−abi−(ni−ri)Xi(ri)−abi+C,
where
C is a constant term. The log-likelihood function is maximized at
a^=X(1)=min{X1(1),…,Xk(1)} (15)
and
b^ic=1ri∑j=1ri(Xi(j)−X(1))+(ni−ri)(Xi(ri)−X(1)) =1ririb^i+ni(Xi(1)−X(1)),i=1,…,k,(16)
where
b^i's are as defined in
(5). Letting
b^c=(b^1c,…,b^kc), we have
l(a^,b^c)=−∑i=1kriln(b^ic)−R+C.(17)
Let Λa denote −2ln(LRT statistic). Then, using l(a^,b^) defined in (6), it can be easily verified that
Λa=−2l(a^,b^c)−l(a^,b^) =−2∑i=1krilnb^ib^i+ni(Xi(1)−X(1))∕ri.(18)
As Λa is invariant under the transformation Xij→ciXij+d, where ci's and d are positive constants, its null distribution does not depend on any parameters. Specifically, the null distribution of Λa can be evaluated empirically assuming that a1=…=ak=0 and b1=…=bk=1. Under this assumption, using the distributional results in (3), we see that
Λa∼−2∑i=1krilnWiWi+ni(Ui−U(1)),(19)
where
Wi=rib^i∼χ2ri−22∕2,
Ui=χ22∕(2ni) and
U(1)=min{U1,…,Uk} and
Wi's and
Ui's are mutually independent. For an observed value of the LRT statistic
Λa, the p-value and the percentiles of
Λa can be estimated using the following algorithm.
Algorithm 2
For a given set of k samples, compute the LRT statistic Λa in (18),
Generate Ui∗=χ22∕ni and Wi∗=χ2ri−22, i=1,…,k,
Set U(1)∗=min{U1∗,…,Uk∗},
Compute Λa∗=−2∑i=1krilnWi∗Wi∗+ni(Ui∗−U(1)∗),
Repeat steps 2 – 4 for a large number of times, say, 100,000,
Proportion of Λa∗'s that are greater than Λa is an estimate of the p-value,
The 100 (1−α) percentile of the 100,000 Λa∗'s is a Monte Carlo estimate of the 100(1−α) percentile of Λa.
In Table 2, we present the estimated percentiles of Λa (based on 1,000,000 runs) for k=3 and 5, and for some values for (r1,…,rk) ranging from very small to large. Notice that these are exact percentiles, except for the simulation errors.
|
(n1,n2,n3)=(9,10,8)
|
|
(n1,…,n5)=(20,…,20)
|
| (r1,r2,r3) |
99 |
95 |
90 |
(r1,…,r5) |
99 |
95 |
90 |
| (3, 2, 3) |
22.43 (22.74) |
15.89 (16.25) |
12.96 (13.33) |
(3, 4, 3, 2, 4) |
31.25 (31.93) |
23.93 (24.64) |
20.58 (21.23) |
| (3, 3, 4) |
19.21 (19.38) |
13.69 (13.85) |
11.21 (11.36) |
(5, 4, 6, 5, 8) |
24.87 (25.19) |
19.18 (19.44) |
16.55 (16.75) |
| (4, 5, 4) |
17.38 (17.44) |
12.40 (12.46) |
10.18 (10.22) |
(6, 3, 7, 7, 10) |
24.93 (25.07) |
19.16 (19.35) |
16.46 (16.67) |
| (5, 5, 5) |
16.68 (16.26) |
11.88 (11.88) |
9.75 (9.74) |
(9, 8, 7, 10, 15) |
22.57 (22.74) |
17.41 (17.55) |
14.99 (15.12) |
| (6, 7, 8) |
15.59 (15.51) |
11.12 (11.08) |
9.11 (9.09) |
(10, …, 10) |
22.32 (22.47) |
17.25 (17.34) |
14.85 (14.94) |
| (9, 10, 8) |
14.99 (14.94) |
10.69 (10.68) |
8.77 (8.75) |
(5, 10, 15, 10, 15) |
22.60 (22.73) |
17.43 (17.55) |
15.02 (15.12) |
|
|
(n1,n2,n3) = (25, 30, 22) |
|
(n1,…,n5)=(60,50,50,70,60) |
|
| (10, 10, 10) |
14.77 (14.72) |
10.53 (10.52) |
8.64 (8.63) |
(20, …, 20) |
21.19 (21.21) |
16.35 (16.37) |
14.08 (14.11) |
| (15, 10, 20) |
14.33 (14.26) |
10.23 (10.19) |
8.39 (8.35) |
(4, 5, 6, 7, 7) |
24.60 (24.93) |
18.98 (19.24) |
16.34 (16.58) |
| (20, 20, 20) |
14.00 (13.95) |
9.98 (9.97) |
8.17 (8.18) |
(30, …, 30) |
20.80 (20.82) |
16.06 (16.07) |
13.84 (13.85) |
|
|
(n1,n2,n3) = (300, 200, 100) |
|
(n1,…,n5)=(200,…,200) |
|
| (30, 40, 50) |
13.61 (13.58) |
9.73 (9.71) |
7.97 (7.96) |
(50, 40, 50, 40, 50) |
20.61 (20.57) |
15.87 (15.88) |
13.68 (13.68) |
| (50, 50, 50) |
13.56 (13.53) |
9.68 (9.67) |
7.93 (7.93) |
(50, …, 50) |
20.45 (20.52) |
15.82 (15.84) |
13.64 (13.65) |
| (100, 100, 100) |
13.39 (13.40) |
9.59 (9.58) |
7.86 (7.85) |
(100, …, 100) |
20.25 (20.30) |
15.68 (15.67) |
13.50 (13.50) |
| (120, 150, 100) |
13.37 (13.38) |
9.54 (9.56) |
7.81 (7.84) |
(150, 130, 140, 130, 150) |
20.22 (20.24) |
15.61 (15.63) |
13.46 (13.46) |
| χ2k−22 percentiles |
13.28 |
9.49 |
7.78 |
χ2k−22 percentiles |
20.09 |
15.51 |
13.36 |
Table 2Monte Carlo estimates of the percentiles of Λa based on (19) and those of Qaχ2k−22.
To find a closed-form approximate distribution for Λa, we note that Wilks' theorem is not applicable to find an asymptotic χ2 null distribution, because the sample space of a two-parameter exponential distribution depends on an unknown parameter. However, our extensive simulation studies indicated that, for large samples, the null distribution of Λa is close to χ2k−22 distribution; see the last row in Table 2. This χ2k−22 approximation is satisfactory only when all ri's are 100 or more. This approximation can be improved along the lines of the preceding section. Specifically, we determine Qa so that E(Λa)=QaEχ2k−22, and then approximate the distribution of Λa by Qaχ2k−22. Toward this, we note that exact calculation of E(Λa) seems to be difficult, and so we resort to use an approximation of E(Λa). It is shown in Appendix that
E(Λa)≈−2∑i=1kriψ(ri−1)−ln(θi)+12θi2ri+ni2∑j=1knj2,
where
θi=ri−ni∕∑j=1knj. Finally, letting
Qa=E(Λa)∕(2k−2), we propose the distribution of
Qaχ2k−22 as an approximate null distribution of
Λa.
Monte Carlo estimates (based on one million runs) of the percentiles of Λa based on Algorithm 2 and those based on the above approximation are given in Table 2 for k=3 and 5 and for some sample sizes and (r1,…,rk). Comparison of the percentiles based on Monte Carlo method and tho approximate ones are very close even for small sample such as (r1,r2,r3)=(3,2,3) and (r1,…,r5)=(3,4,3,2,4). On an overall basis, we see that the approximate percentiles are very satisfactory and are safe to use in practical applications.
5. HOMOGENEITY TEST
For the homogeneity test, the hypotheses of interest are
H0:(a1,b1)=…=(ak,bk)vs.Ha:(ai,bi)≠(aj,bj)for some i≠j.(20)
Let (a,b) denote the common unknown parameter under H0 in (20). Then the log-likelihood function under the null hypothesis (20) can be expressed as
l(a,b)=∑i=1k−rilnb−∑j=1riXi(j)−ab−(ni−ri)Xi(ri)−ab+C,
where
C is a constant term. The log-likelihood function is maximized at
a^=X(1)=min{X1(1),…,Xk(1)} (21)
and
b^=1R∑i=1k∑j=1ri(Xi(j)−X(1))+(ni−ri)(Xi(ri)−X(1)) =1R∑i=1krib^i+ni(Xi(1)−X(1)),(22)
where
b^i is defined in
(5) and
R=∑i=1kri. Let
Λh denote the
−2ln(LRT statistic). It can be easily verified that
Λh=−2∑i=1krilnb^ib^ =−2∑i=1krilnb^i1R∑i=1krib^i+1R∑i=1kni(Xi(1)−X(1)).(23)
Since the above statistic is location-scale invariant, its null distribution does not depend on any parameters, and so the null distribution can be evaluated empirically assuming, without loss of generality, that a1=…=ak=0 and b1=…=bk=1. Under this assumption and using the distributional results in (3), we see that
Λh∼−2∑i=1krilnWi1R∑i=1kriWi+∑i=1kni(Ui−U(1))(24)
where
Wi=χ2ri−22∕(2ri),
Ui=χ22∕(2ni) and
U(1)=min{U1,…,Uk} and the random variables
Wi's and
Ui's are mutually independent. Thus, the percentiles of
Λh can be estimated using Monte Carlo simulation as in
Algorithms 1 and
2.
Alternatively, a closed-form approximation to the null distribution of Λh can be obtained along the lines of Section 3. As noted earlier, Wilks' theorem is not applicable to find a large sample approximate χ2 distribution. On the basis of simulation studies we found that Λh∼χ3(k−1)2 for large samples. This approximation can be improved along the lines of Section 3. That is, we determine the value of Qh so that E(Λh)=E(Qhχ3(k−1)2) and Λh∼Qhχ3(k−1)2, approximately. To find E(Λh), we note the b^i follows a χ2ri−22∕(2ri) distribution and b^ in (22) follows a χ2R−22∕(2R). Using this distributional results in (23), it can be verified that
EΛh=2∑i=1kriψ(R−1)−ψ(ri−1)+lnri∕R.(25)
Thus, Qh=EΛh∕(3(k−1)) and
Λh∼Qhχ3(k−1)2,approximately.
To judge the accuracy of the above approximation, we estimated the upper percentiles of Λh in (24) using simulation with 1,000,000 runs and the percentiles of Qhχ3(k−1)2 for k=3 and 4, and some small to large values of sample sizes. These percentiles are reported in Table 3. Comparison of Monte Carlo estimates of the percentiles and those of Qhχ2k2 shows that the latter approximate percentiles are very close to those based on simulation even when (r1,r2,r3)=(3,2,2). Thus, the improved approximation can be safely used for all values of ri's greater than or equal to three.
|
Percentiles
|
|
Percentiles
|
| (r1,r2,r3) |
99 |
95 |
90 |
(r1,…,r4) |
99 |
95 |
90 |
| (3, 2, 3) |
28.10 (27.70) |
20.92 (20.74) |
17.62 (17.54) |
(3, 2, 3, 2) |
36.73 (36.19) |
28.52 (28.26) |
24.69 (24.53) |
| (3, 3, 4) |
24.10 (24.04) |
18.04 (18.00) |
15.24 (15.22) |
(4, 3, 4, 5) |
28.73 (28.55) |
22.38 (22.30) |
19.40 (19.35) |
| (5, 5, 5) |
20.84 (20.86) |
15.62 (15.63) |
13.19 (13.21) |
(5, 6, 7, 9) |
25.37 (25.29) |
19.78 (19.75) |
17.16 (17.14) |
| (4, 5, 4) |
21.80 (21.77) |
16.29 (16.31) |
13.77 (13.78) |
(9, 8, 10, 9) |
24.12 (24.12) |
18.83 (18.84) |
16.35 (16.35) |
| (5, 4, 5) |
21.35 (21.32) |
15.98 (15.97) |
13.50 (13.50) |
(10, …, 10) |
23.79 (23.83) |
18.60 (18.61) |
16.14 (16.15) |
| (7, 7, 10) |
19.26 (19.22) |
14.42 (14.39) |
12.17 (12.17) |
(15, 10, 15, 10) |
23.43 (23.45) |
18.30 (18.31) |
15.90 (15.89) |
| (10, 10, 10) |
18.62 (18.61) |
13.94 (13.94) |
11.77 (11.78) |
(15, 15, 15, 15) |
23.06 (23.06) |
18.03 (18.01) |
15.65 (15.63) |
| (15, 15, 15) |
17.94 (17.97) |
13.46 (13.46) |
11.37 (11.38) |
(20, 15, 20, 15) |
22.89 (22.88) |
17.84 (17.87) |
15.49 (15.51) |
| (10, 15, 20) |
18.11 (18.09) |
13.55 (13.55) |
11.46 (11.46) |
(20, 20, 20, 20) |
22.67 (22.69) |
17.73 (17.72) |
15.38 (15.38) |
| (50, 50, 50) |
17.11 (17.14) |
12.83 (12.84) |
10.85 (10.85) |
(50, …, 50) |
22.11 (22.06) |
17.24 (17.23) |
14.96 (14.95) |
| (100, 100, 100) |
16.93 (16.98) |
12.70 (12.71) |
10.72 (10.75) |
(100, …, 100) |
21.86 (21.86) |
17.09 (17.07) |
14.82 (14.82) |
|
| χ3(k−1)2 percentiles |
16.81 |
12.59 |
10.64 |
χ3(k−1)2 percentiles |
21.67 |
16.92 |
14.68 |
Table 3Monte Carlo estimates of the percentiles of Λh in (9) and those of Qhχ2k2.
6. POWER STUDIES
6.1. Power Comparison of the Tests for Scale Parameters
To compare the LRT, iterative test (ITER) and the union-intersection test (UIT), we estimated the powers of these tests for k=3 and 5, and some moderate sample sizes using simulation. As all the tests are scale invariant, the powers were estimated assuming, without loss of generality, that b1=1 and b1>…>bk. The estimated powers at the level .05 are reported in Table 4. Examination of reported powers indicates that the UIT performs better than others only when the sample sizes are equal, one of the parameters is maximum and only one of the parameters is minimum; see the cases of (r1,r2,r3)=(10,10,10), (b1,b2,b3)=(1,1,.4) and (1,1,.3), (r1,…,r5)=(10,…,10), (b1,…,b5)=(1,.4,1,1,1), (1,.2,1,1,1). In these cases, the UIT has larger power than those of the other tests. However, for (r1,…,r5)=(5,15,10,5,15) and (b1,…,b5)=(1,.4,1,1,1), the power of the UIT is much smaller than those of other two tests. Specifically, we see under the column of (r1,…,r5)=(5,15,10,5,15) the powers of UIT are much smaller than the powers of other two tests. Between the LRT and ITER, the ITER is more powerful than the LRT in some cases; see the powers for the case (r1,…,r5)=(5,15,10,5,15). However, the LRT has appreciably larger powers than the ITER over larger parameter space, and so the LRT maybe preferred to other two tests.
Sample Sizes
|
|
(10, 10, 10)
|
(10, 10, 20)
|
|
(10, 10, 10, 10, 10)
|
(5, 15, 10, 5, 15)
|
| (b1,b2,b3) |
ITER |
UIT |
LRT |
ITER |
UIT |
LRT |
(b1,…,b5) |
ITER |
UIT |
LRT |
ITER |
UIT |
LRT |
| (1, 1, 1) |
.050 |
.050 |
.050 |
.050 |
.050 |
.050 |
(1, 1, 1, 1, 1) |
.050 |
.050 |
.050 |
.050 |
.050 |
.050 |
| (1, .8, .6) |
.121 |
.134 |
.138 |
.182 |
.108 |
.311 |
(1, .9, 1, .7, 1) |
.082 |
.089 |
.091 |
.064 |
.084 |
.079 |
| (1, .8, .5) |
.203 |
.228 |
.228 |
.342 |
.196 |
.377 |
(1, .6, 1, .9, .6) |
.153 |
.168 |
.193 |
.189 |
.047 |
.132 |
| (1, .5, .5) |
.265 |
.287 |
.310 |
.311 |
.293 |
.375 |
(1, .6, 1, .6,1) |
.173 |
.190 |
.210 |
.197 |
.111 |
.193 |
| (1, .4, .6) |
.373 |
.364 |
.377 |
.377 |
.411 |
.392 |
(1, .6, 1, .5, .6) |
.207 |
.253 |
.292 |
.185 |
.141 |
.190 |
| (1, .4, .7) |
.374 |
.365 |
.367 |
.375 |
.419 |
.409 |
(1, .4, 1, 1, 1) |
.367 |
.419 |
.384 |
.505 |
.100 |
.438 |
| (1, .4, .8) |
.385 |
.388 |
.387 |
.392 |
.457 |
.703 |
(1, .4, 1, 1, .6) |
.369 |
.410 |
.431 |
.462 |
.090 |
.347 |
| (1, .4, .5) |
.391 |
.395 |
.422 |
.409 |
.424 |
.504 |
(1, .5, 1, .3, .5) |
.547 |
.620 |
.675 |
.339 |
.360 |
.432 |
| (1, 1, .4) |
.443 |
.463 |
.453 |
.703 |
.473 |
.614 |
(1, .6, 1, .2, .6) |
.863 |
.879 |
.862 |
.399 |
.624 |
.524 |
| (1, .4, .4) |
.437 |
.473 |
.508 |
.504 |
.484 |
.709 |
(1, .2, 1, 1, 1) |
.874 |
.945 |
.917 |
.954 |
.633 |
.964 |
| (1, .3, .7) |
.611 |
.615 |
.601 |
.614 |
.677 |
.917 |
(1, .2, 1, .3, 1) |
.910 |
.949 |
.972 |
.959 |
.706 |
.978 |
| (1, .4, .3) |
.556 |
.627 |
.655 |
.709 |
.652 |
.937 |
(1, .2, .5, .2, 1) |
.934 |
.959 |
.979 |
.952 |
.745 |
.973 |
| (1, 1, .3) |
.706 |
.721 |
.710 |
.917 |
.775 |
.922 |
(1, .2, .4, .2, 1) |
.939 |
.957 |
.980 |
.962 |
.748 |
.973 |
| (1, .4, .2) |
.796 |
.851 |
.859 |
.937 |
.887 |
.995 |
(1, .2, .7, .2, 1) |
.937 |
.966 |
.985 |
.948 |
.768 |
.979 |
| (1, .3, .2) |
.819 |
.870 |
.890 |
.922 |
.932 |
.896 |
(1, .2, 1, .2, 1) |
.957 |
.981 |
.994 |
.964 |
.820 |
.990 |
| (1, 1, .2) |
.940 |
.941 |
.936 |
.995 |
.989 |
.969 |
(1, .2, 1, .2, .2) |
.958 |
.978 |
.993 |
.948 |
.970 |
.990 |
Table 4Powers of the iterative test (ITER), UIT and the LRT for testing the equality of the scale parameters.
6.2. Power of the LRT for Location Parameters
To judge the power properties of the LRT for the equality of location parameters, we estimated the powers by Monte Carlo simulation and reported them in Table 5. The powers were estimated for (b1,b2,b3)=(3,4,5) and (2,2,2) and sample sizes (n1,n2,n3)=(15,15,15) and (20,20,20). Given a set of sample sizes, we evaluated the powers at (r1,r2,r3)=(5,4,5) and (12,11,12) and for values of (a1,a2,a3) so that a1=4>a2>a3. We first observe from Table 5 that the powers are increasing with increasing sample sizes. For example, see the powers in the columns of (n1,n2,n3)=(15,15,15) and (r1,r2,r3)=(5,4,5), and (n1,n2,n3)=(20,20,20) and (r1,r2,r3)=(5,4,5). We also observe that the powers for larger values of (r1,r2,r3) are larger than those of smaller values of (r1,r2,r3) while the sample sizes and other parameters are fixed. For example, see the powers in the columns of (n1,n2,n3)=(15,15,15), (r1,r2,r3)=(5,4,5) and (12,11,12); (n1,n2,n3)=(20,20,20), (r1,r2,r3)=(5,4,5) and (12,11,12). Finally, we also notice that the powers for the case of (b1,b2,b3)=(3,4,5) are smaller than corresponding powers for (b1,b2,b3)=(2,2,2). This is expected because the variance of an exp(a,b) distribution is b2, and so the powers for the set of k populations with larger variances are expected to be smaller than those for populations with smaller variances while all other parameters and sample sizes are fixed. Thus, the LRT for the equality of location parameters possess all natural properties of an efficient test.
| (b1,b2,b3)=(3,4,5) |
(b1,b2,b3)=(2,2,2) |
|
| (n1,n2,n3) |
(15, 15, 15) |
(20, 20, 20) |
(15, 15, 15) |
(20, 20, 20) |
|
|
(r1,r2,r3)
|
(r1,r2,r3)
|
| (a1,a2,a3) |
(5, 4, 5) |
(12, 11, 12) |
(5, 4, 5) |
(12, 11, 12) |
(5, 4, 5) |
(12, 11, 12) |
(5, 4, 5) |
(12, 11, 12) |
| (4, 4, 4) |
.049 |
.051 |
.050 |
.049 |
.051 |
.050 |
.050 |
.050 |
| (4, 3.8, 4) |
.072 |
.080 |
.095 |
.104 |
.164 |
.195 |
.270 |
.341 |
| (4, 3.5, 3.5) |
.140 |
.151 |
.188 |
.233 |
.308 |
.453 |
.444 |
.673 |
| (4, 3.6, 4) |
.166 |
.195 |
.268 |
.340 |
.515 |
.683 |
.741 |
.877 |
| (4, 4, 3.6) |
.176 |
.198 |
.401 |
.518 |
.532 |
.677 |
.731 |
.876 |
| (4, 3, 3) |
.315 |
.450 |
.448 |
.681 |
.674 |
.942 |
.823 |
.994 |
| (4, 4, 3.5) |
.249 |
.301 |
.398 |
.513 |
.689 |
.838 |
.879 |
.952 |
| (4, 3, 3.5) |
.459 |
.618 |
.679 |
.839 |
.915 |
.974 |
.983 |
.996 |
| (4, 4, 3) |
.699 |
.836 |
.875 |
.953 |
.985 |
.996 |
.998 |
.999 |
Table 5Powers of the LRT for equality of location parameters.
6.3. Powers of the Homogeneity Test
Powers of the LRT for testing H0:(a1,b1)=…=(ak,bk) were estimated for sample sizes (n1,n2,n3)=(15,15,15) and (20,20,20), and for each set of sample sizes, (r1,r2,r3)=(7,8,11) and (14,13,10). The estimated powers are reported in Table 6. On the basis powers in Table 6, we see that the power properties of the homogeneity test is very similar to the LRT for location parameters discussed in the preceding section. In particular, the power is increasing with increasing sample sizes while other parameters are fixed and also increasing with increasing values of ri's while other values are fixed. For example, see the columns under (n1,n2,n3)=(15,15,15)(r1,r2,r3)=(7,8,11) and (14,13,10) in Table 6. We also notice that the power is increasing with increasing disparities among location parameters and/or scale parameters. Thus, the LRT has all natural properties of an efficient test.
| (n1,n2,n3) |
|
(15, 15, 15) |
(20, 20, 20) |
|
|
|
(r1,r2,r3)
|
| (a1,a2,a3) |
(b1,b2,b3) |
(7, 8, 11) |
(14, 13, 10) |
(7, 8, 11) |
(14, 13, 10) |
| (2, 2, 2) |
(3, 3, 3) |
.049 |
.048 |
.049 |
.048 |
| (2, 2, 2) |
(3, 3, 2) |
.083 |
.094 |
.084 |
.097 |
| (2, 1.8, 2) |
(3, 3, 3) |
.094 |
.090 |
.129 |
.133 |
| (2, 1.8, 2) |
(3, 3, 2) |
.141 |
.149 |
.195 |
.209 |
| (2, 1.8, 1.7) |
(3, 3, 2) |
.202 |
.218 |
.285 |
.301 |
| (2, 1.5, 1.7) |
(3, 3, 2) |
.314 |
.336 |
.506 |
.512 |
| (2, 1.5, 1.7) |
(3, 3, 1.5) |
.498 |
.520 |
.680 |
.699 |
| (2, 1.5, 1.7) |
(3, 3, 1) |
.792 |
.821 |
.902 |
.909 |
| (2, 1.3, 1.7) |
(3, 3, 1) |
.920 |
.928 |
.974 |
.975 |
| (2, 1.3, 1.7) |
(3, 1.5, 1) |
.982 |
.989 |
.998 |
.999 |
Table 6Powers of the homogeneity test.
7. AN EXAMPLE
The data were collected from a placebo controlled randomized clinical trial of gamma interferon in CGD. The data were subsets of large data given in Appendix D of [1]. They represent elapsed time (in days) from randomization to diagnosis of a serious infection. The samples were taken from three hospitals with ID codes 238, 204 and 332. We applied Kolmogorov–Smirnov test for the exponential distribution by [8], and the test indicates that all three samples satisfy the model assumption of two-parameter exponential distribution.
Notice that the samples are uncensored with (n1,n2,n3)=(10,14,9). The MLEs are a^i=Xi(1) and b^i=X¯i−Xi(1), i=1,2,3, and the MLEs for the three samples are reported in Table 7.
| Hospitals |
Time to Infection in Days |
MLEs (a^,b^) |
| 1 |
253 |
294 |
19 |
373 |
334 |
238 |
118 |
240 |
99 |
167 |
(19, 194.5) |
| 2 |
373 |
26 |
152 |
241 |
322 |
350 |
211 |
307 |
82 |
|
|
|
114 |
337 |
18 |
267 |
104 |
|
|
|
|
|
(18, 189.4) |
| 3 |
146 |
188 |
304 |
91 |
121 |
203 |
264 |
236 |
207 |
|
(91, 104.6) |
For testing H0:a1=a2=a3vs.Ha:ai≠aj for some i≠j, the LRT statistic Λa=9.63 and the 95th percentile of the null distribution is 10.51. To find the approximate 95th percentile, we evaluated Qa=E(Λa)∕(2k−2)=1.105 and 1.105×χ4;.952=1.105×9.488=10.48. Notice that the exact and the approximate percentiles are practically the same. Furthermore, we estimated the p-value using Algorithm 2 with 1000000 runs as .069 and the p-value based on the Qaχ42 distribution is also .069. Thus, at 5% level, the location parameters are not significantly different.
For testing H0:b1=b2=b3, the statistic Λb is 2.17 and the p-value based on Algorithm 1 is .385. If we use the level of significance .05, then the exact critical value is 6.83 and Qbχ2;.952=1.138×5.99 is 6.82, which is very close to the exact one. To apply the UIT, we estimated the 5% critical value as 3.056, and the test statistic max{b^1,b^2,b^3}∕min{b^1,b^2,b^3}=194.5∕104.6=1.859. To apply the iterative test procedure, we calculated the individual test statistics F2=.944 and F3=.563. We chose the levels for the individual test as η2=η3=η4=.01266 so that α=1−(1−2×.01266)2=.05. For these levels, the (lower, upper) critical values are (.383,2.84) for F2 and (.349,2.36) for F3. Note that both test statistics F2 and F3 fall in their corresponding acceptance intervals, and so the null hypothesis of equality of scale parameters is accepted. Thus, all three tests indicate that H0:b1=b2=b3 is tenable.
For testing H0:(a1,b1)=(a2,b2)=(a3,b3), the statistic Λh in (23) is 9.68 and the p-value based on (24) is .187. If we use the level of significance .05, then the exact critical value is 13.93; the approximate critical value Qhχ6;.952=1.101×12.59=13.86. Thus, homogeneity of exponential distributions is tenable. This means that all three samples may be pooled and the pooled sample can be modeled by a single two-parameter exponential distribution.
8. CONCLUDING REMARKS
Although several tests for comparing parameters of several two-parameter exponential distributions were proposed in the literature, none of them is based on the likelihood approach. As the sample space of a two-parameter exponential distribution depends on an unknown parameter, it does not satisfy all the regularity conditions. So an LRT statistic for comparing parameters does not have the asymptotic chi-square distribution with the degrees of freedom determined by the difference between the dimensions of the parameter spaces under the alternative and null hypotheses. We have shown in this article that the null distributions of the LRT statistics for all three problems do not depend on any unknown parameters, and so the LRTs are exact. However, calculation of the percentiles of the LRT statistics involves simulation. Even though calculation of percentiles based on Monte Carlo simulation is not difficult, we provided closed-form approximate chi-square distributions for all three problems. These approximate chi-square null distributions are not only simple, but also very accurate even for small samples. These approximate null distributions maybe warranted in situations where the number of distributions to be compared is large.
ACKNOWLEDGMENTS
The authors are grateful to a reviewer for providing useful comments and suggestions.
APPENDIX
Let Zi=Wi+ni(Ui−U(1)), θi=E(Zi), i=1,‥,k and g(x)=ln(x). Recall that Wi's and Ui's are independent with
Wi∼12χ2ri−22andUi∼12niχ22,i=1,…,k.
To find an approximation to Eln(Zi), we shall use the result (see Section 2.9, [9]) that
Eln(Zi)≈ln(θi)+Var(Zi)g″(θi)2.
Recall that niUi∼χ22∕2, or niUi follows a standard exponential distribution. So
P(U(1)≤u)=1−∏ikP(niUi≥niu)=1−exp−∑i=1kniu,
which is the distribution function of an exponential distribution with mean
1∕∑i=1kni. Thus,
E(U(1))=1∑i=1kniandVar(U(1))=1∑i=1kni2(A.1)
Using the above expectation, we have
E(Zi)=E(Wi)+niE(Ui−U(1))=ri−ni∑j=1knj,i=1,…,k.
To find the variance,
Var(Zi)=Var(Wi)+ni2Var(Ui)+Var(U(1))−2 Cov (Ui,U(1)).(A.2)
To find the Cov(Ui,U(1)), we find the joint distribution of Ui and U(1) as
FUi,U1u,v=PUi≤u,U1≤v =PUi≤u−PUi≤u,U1≥v =PUi≤u−PUi≤u,U1≥v,…,Uk≥v =PUi≤u−Pv≤Ui≤u∏j≠iPnjUj≥njv =1−e−niu−e−niv−e−niue−∑j≠injv.
By taking the derivative with respect to (u,v), we find the joint density as
fUi,U(1)(u,v)=nie−niu∑j≠inje−∑j≠injv,
which shows that
Ui and
U(1) are independent. Using this fact in
(A.2), we find
Var(Zi)=ri+ni2∑j=1knj2.
Thus,
Eln(Zi)≈ln(θi)−Var(Zi)2θi2 =ln(θi)−12θi2ri+ni2∑j=1knj2,
where
θi=E(Zi)=ri−ni∑j=1knj,
i=1,…,k. Using this approximation, we see that
E(Λa)=−E2∑i=1krilnχ2ri−22∕2Zi =−2∑i=1kriψ(ri−1)−Eln(Zi) ≈−2∑i=1kriψ(ri−1)−ln(θi)+12θi2ri+ni2∑j=1knj2.(A.3)
REFERENCES
1.T.R. Fleming and D.P. Harrington, Counting Processes and Survival Analysis, John Wiley & Sons, 2011. 
