- Email: [email protected]
Confidence regions for comparison of two normal samples
Statistics and Probability Letters 119 (2016) 273–280
Contents lists available at ScienceDirect
Statistics and Probability Letters journal homepage: www.elsevier.com/locate/stapro
Confidence regions for comparison of two normal samples Andrew L. Rukhin Statistical Engineering Division, National Institute of Standards and Technology, 100 Bureau Drive, Gaithersburg, MD 20899, USA
article
info
Article history: Received 12 April 2016 Received in revised form 18 August 2016 Accepted 19 August 2016 Available online 30 August 2016
abstract The comparison of two normal populations is considered when one is interested not only in the difference between the means but also in the ratio of the variances. Several joint confidence regions for these two parameters are obtained. Published by Elsevier B.V.
Keywords: Bartlett correction Behrens–Fisher problem Higher order inference Minimum area region
1. Introduction The problem of two normal populations raises important statistical issues and is often encountered in practice when two different methods (instruments, treatments, materials, etc.) are to be compared. The classical monograph by Lehmann and Romano (2003) dedicates a whole section to this subject whose theoretical importance is due to the absence of a single pivotal quantity on which to base statistical inference in the classical Behrens–Fisher problem. Indeed there seems to be no unique universally accepted method either to test equality of the means or to derive a confidence interval covering their difference. However two samples even with the same mean can be quite dissimilar unless their variances are not too different. Therefore the variances must enter in a meaningful comparison calling for a hypothesis test about the value of their ratio or for a joint confidence region for the difference in the two means and this ratio. We review several such confidence regions in Section 3. These include the exact regions produced by the likelihood ratio and its Bartlett and two-moments based corrections, the minimum area set, as well as the Mood-type confidence set. One of the suggested procedures is obtained from the modern higher order asymptotic statistical inference (Brazzale et al., 2007). The situation when the variances ratio is known is briefly discussed in Section 2. Section 4 concludes the paper with some recommendations. 2. Comparing means and variances: known variances ratio Here we assume given two normal observables x¯ ∼ N (µ + ∆, θ ) and y¯ ∼ N (µ, η), as well as two statistics s2 and u2 which are independent of x¯ and y¯ . Their probability distributions are of the form θ χ 2 (n − 1)/(n − 1) and ηχ 2 (m − 1)/(m − 1). For example, s2 is an unbiased estimator of θ which is the variance of the first sample mean. Thus the reduction from two random samples to sufficient statistics (¯x, y¯ , s2 , u2 ) is made, θ = σ12 /n, η = σ22 /m, where σ12 , σ22 are the samples’ variances, and n, m the respective sample sizes.
E-mail address: [email protected]. http://dx.doi.org/10.1016/j.spl.2016.08.011 0167-7152/Published by Elsevier B.V.
274
A.L. Rukhin / Statistics and Probability Letters 119 (2016) 273–280
The log-likelihood function corresponding to x¯ , y¯ , s2 and u2 involves a four-dimensional parameter ψ = (∆, µ, θ , η). It has the form
(¯y − µ)2 + (m − 1)u2 + n log θ + m log η . (1) 2 θ η In this section we look at the simplest case when the ratio of two variances is given, ρ = η/θ . Then λ = (µ, θ ) is the nuisance parameter and the likelihood for ψ = (∆, λ) can be obtained from (1) by replacing η by ρθ . The maximum ˆ = x¯ − y¯ , and that of µ is µ likelihood estimator θˆ of θ is [(n − 1)s2 + (m − 1)u2 /ρ]/(m + n), while that of ∆ is ∆ ˆ = y¯ . The constrained maximum likelihood estimator of µ for fixed ∆, L(ψ) = −
µ ˜ =
1 (¯x − µ − ∆)2 + (n − 1)s2
x¯ − ∆ + y¯ /ρ 1 + 1/ρ
+
,
(2)
is the weighted mean of x¯ − ∆ and y¯ . The constrained maximum likelihood estimator θ˜ satisfies the relation, θ˜ = ˆ − ∆)2 , and the signed likelihood root statistic has the form θˆ + (n + m)−1 (ρ + 1)−1 (∆
ˆ − ∆) (n + m) log θ˜ /θˆ r = r (∆) = sign(∆
1/2
.
The statistic r is known to be approximately normally distributed although this approximation typically can be poor for 2 ˆ − ∆|2 ≤ [exp(zα/ small sample sizes. For example, the confidence interval {∆ : |r | ≤ zα/2 } = {∆ : |∆ 2 /(n + m)) − 1](n +
m)(ρ + 1)θˆ } is too short and does not maintain the nominal confidence coefficient well. A much more accurate normal approximation can be obtained via the Bartlett correction discussed in more detail in Section 3. Another possibility is to employ a third order accurate statistic (Section 8.5 of Brazzale et al., 2007). One of these, the
ˆ − ∆)θˆ / (ρ + 1)θ˜ 3 , provides third order improvement on r, with the confidence interval, Skovgaard (1996) statistic, v = (∆ {∆ : |r + log(v/r )| ≤ zα/2 }.
(3)
This interval turns out to be very close to
ˆ ± tα/2 (n + m − 2) (n + m)(n + m − 2)−1 (ρ + 1)θˆ , ∆ which is the usual t-interval. Here zα/2 and tα/2 (n + m − 2) denote the familiar critical points. The tests of the hypothesis: ∆ = 0, ρ = 1 were studied by Perng and Littell (1976), and more recently by Zhang et al. (2012). The goal of this work is to construct a joint confidence region for (∆, ρ). 3. Comparing means and variances: unknown variances ratio With ψ = (∆, ρ; θ , µ) = (∆, ρ; λ), the score vector ℓ = ℓ(ψ; x¯ , y¯ , s2 , u2 ), has the form, x¯ − µ − ∆ y¯ − µ , + , θ θ η (¯x − µ − ∆)2 + (n − 1)s2 n (¯y − µ)2 + (m − 1)u2 m T − , − . 2θ 2 2θ 2η 2 2η ˆ = x¯ − y¯ , the maximum likelihood estimator of µ is y¯ , ρˆ = n(m − The maximum likelihood estimator of ∆ remains ∆ 1)u2 /[m(n − 1)s2 ], θˆ = [(n − 1)s2 + (m − 1)u2 /ρ]/( ˆ m + n) = (n − 1)s2 /n. The Fisher information matrix has the form 1 1
x¯ − µ − ∆
θ 0 I (ψ) = 1 θ 0
0
m 2ρ 2 0 m 2ρθ
0
0
ρ+1 ρθ 0
θ
2ρθ . 0 n + m m
2θ 2
Now the parameter of interest has two components (∆, ρ); the constrained maximum likelihood estimators of µ and θ for fixed ∆, ρ determine the likelihood ratio test of the hypothesis postulating these null values. The first of these estimators µ ˜ has already been found in (2). The estimator θ˜ of θ can be written as
(¯x − y¯ − ∆)2 /(ρ + 1) + (n − 1)s2 + (m − 1)u2 /ρ n+m ˆ − ∆)2 (∆ (mρˆ + nρ)θˆ = + , (n + m)(ρ + 1) (n + m)ρ
θ˜ =
(4)
A.L. Rukhin / Statistics and Probability Letters 119 (2016) 273–280
275
leading to the likelihood ratio (Neyman–Pearson) statistic,
ˆ , ρ, R = R(∆, ρ) = 2[L(∆ ˆ y¯ , θˆ ) − L(∆, ρ, µ, ˜ θ˜ )] = (n + m) log(θ˜ /θˆ ) + m log(ρ/ρ). ˆ Thus the joint confidence region for ∆ and ρ based on R is
(∆, ρ) : (n + m) log(θ˜ /θˆ ) + m log(ρ/ρ) ˆ ≤ χα2 (2) ˆ − ∆)2 (∆ m(ρˆ − ρ) ρ 2 = (∆, ρ) : (n + m) log 1 + + ≤ χα (2) . + m log (n + m)ρ ρˆ (n + m)(ρ + 1)θˆ
(5)
Here χα2 (2) is the α -critical point of a χ 2 -distribution with 2 degrees of freedom, use of which is recommended by the first order asymptotic theory. However unless both n and m are large, the cut-off constant χα2 (2) of the region (5) is too small leading to insufficient coverage probability. ˆ and Y = mρ/( ˆ − ∆)2 /[n(ρ + 1)θ] Clearly R depends only on Z = (∆ ˆ nρ), which is a two-dimensional pivot. The joint ˆ − ∆)2 /[(ρ + 1)θ], nθˆ /θ , and distribution of (Z , Y ) can be obtained from that of three independent χ 2 -random variables, (∆ ˆ mρˆ θ/[ρθ ] which have degrees of freedom 1, n − 1 and m − 1 respectively. This distribution provides the joint density of Z and Y , p(z , y) =
Γ ((n + m − 1)/2)y(m−3)/2 . Γ (1/2)Γ ((n − 1)/2)Γ ((m − 1)/2)z 1/2 (1 + z + y)(n+m−1)/2
One can considerably improve accuracy of the χ 2 -approximation in (5) by using the Bartlett correction (Cordeiro and Cribari-Neto, 2014). In our situation this correction is simplified by the fact that the distribution of the likelihood ratio statistic does not depend on unknown parameters,
n−1 m−1 − nΨ − mΨ − (n + m) log(n + m) + n log(n) + m log(m) 2 2 2 11(n2 + nm + m2 ) 1 =2 1+ +O . 12nm(n + m) min(n, m)2
ER = (m + n)Ψ
n+m−1
Here Ψ denotes the well-known log-derivative of the gamma-function. Thus the Bartlett corrected likelihood ratio statistic,
−1 11(n2 + nm + m2 ) 11(n2 + nm + m2 ) ≈R 1+ ≈R 1− , ER 12nm(n + m) 12nm(n + m)
2R
admits much better χ 2 -approximation than R, so that the nominal confidence coefficient of the confidence region
(∆, ρ) : R ≤ χα2 (2)ER/2
(6)
is considerably closer to 1 − α than that of (5). For example, the variance of the Bartlett corrected likelihood ratio statistic, which can be obtained from the formula Var(R) = n Ψ 2
′
n−1
2
+m Ψ 2
′
m−1
2
− ( m + n) Ψ 2
′
n+m−1
2
,
coincides with Var(χ 2 (2)) = 4 up to the order O min(n, m)−2 .
More general confidence sets (∆, ρ) : Q ≤ K , where K is a threshold constant and Q is a linear combination of log(1 + Z + Y ) and log(Y ), Q = a log(1 + Z + Y ) − b log(Y ) − c ;
a.b > 0,
are worthy of consideration. For example, Q = R, when a = m + n, b = m, c = (n + m) log(n + m)− n log(n)− m log(m). The mean EQ and the variance Var(Q ) admit explicit parameter-free formulas. A possible choice of constants is a = a0 , b = b0 which equate the first two moments of Q to those of the χ 2 (2)-distribution,
2 = EQ = a Ψ
n+m−1 2
− log
n+m
2
n m m−1 n−1 − ( a − b) Ψ − log −b Ψ − log 2
2
2
(7)
2
and 4 = Var(Q ) = (a − b)2 Ψ ′
n−1 2
+ b2 Ψ ′
m−1 2
− a2 Ψ ′
n+m−1 2
.
(8)
276
A.L. Rukhin / Statistics and Probability Letters 119 (2016) 273–280
This choice leads to the region
n ρ (x − ∆)2 mρˆ n+m + b0 log ≤ χα2 (2) + a0 log + b0 log . (∆, ρ) : Q0 = a0 log 1 + + nρ ρˆ n m n(ρ + 1)θˆ
(9)
When n = m, one gets b0 = a0 /2, which is determined by the first moment identity so that the second equation is not needed. In this situation Q0 coincides with the Bartlett corrected likelihood ratio statistic above. When min(n, m) ≥ 5, there are two positive solutions a0 to the simultaneous equations above; the largest root is recommended if m > n. √ The linear transform of Q , 2[(Q − EQ )/ Var(Q ) + 1] admits a much better χ 2 (2)-approximation than Q itself. This fact gives an approximate formula for the threshold constant
K = K (Q ) = EQ + (χα2 (2)/2 − 1) Var(Q ),
(10)
where EQ and Var(Q ) are found from (7), (8). Under this choice the confidence coefficient of the region
(∆, ρ) : a log(1 + Z + Y ) − b log(Y ) ≤ K
(11)
is close to 1 − α . The confidence sets in the estimation problem for a normal mean and variance were studied by Arnold and Shavelle (1998). Frey et al. (2009) discuss the minimum area regions for these parameters, whose counterparts can be derived in our problem. Indeed if C is (1−α)-confidence set for (∆, ρ), let C˜ denote the set of (z , y)-values such that (∆, ρ) ∈ C if and only if (z , y) ∈ C˜ . The formula for the area of C can be obtained from the Jacobian J = J (z , y) of the transformation (∆, ρ) → (z , y),
θˆ (y + mρ/ ˆ n)/(nzy5 ). For all confidence sets C considered here, C˜ = (z , y) : z ≤ g (y, ρ), ˆ y− ≤ y ≤ y+ with a
J = mρˆ
positive function g and positive y− , y+ . The area of such a region depends on ρˆ (and on θˆ ),
d∆ dρ =
C
C˜
J (z , y) dz dy =
2mρˆ
θˆ
n
y+
g (y, ρ)( ˆ ny + mρ) ˆ
1/2
dy.
y5
y−
(12)
For (11), g (y, ρ) ˆ = g (y) = Ayκ − 1 − y, A = exp(K /a), κ = b/a and y− , y+ are the roots of the equation g (y) = 0. Since C˜ p(z , y) dz dy = 1 − α , the fundamental Neyman–Pearson lemma gives the form of the minimum area region,
C˜ = {(z , y) : J (z , y) ≤ K ′ p(z , y)}, i.e.
ˆ − ∆)2 (∆ mρˆ C = (∆, ρ) : (n + m − 1) log 1 + + ˆ nρ n(ρ + 1)θ
+ (m + 1) log
nρ mρˆ
+ log(1 + ρ) ≤ K .
(13)
When ρˆ is small, C has approximately the form (11) with a = m + n − 1, b = m + 1, for large values of ρˆ it has the similar form with a = m + n − 1, b = m + 2. In these two cases the threshold constant K can be found from (10) by using the identities (7) and (8). If ρˆ = n/m, the explicit form of these two moments follows from the fact that Z /(1 + Z + Y ) and 1/(1 + Y ) are independent and beta-distributed with parameters 1/2, (m + n)/2 − 1 and (n − 1)/2, (m − 1)/2 respectively, but in the general case to determine the constant K after the confidence coefficient 1 − α , one has to use numerical integration or Monte Carlo evaluations. When max(n, m) ≥ 8, min(n, m) ≥ 3, numerical results favor the region (11) with a = m + n − 1, b = m + 2, say, Q1 = (m + n − 1) log(1 + Z + Y ) − (m + 2) log(Y ). The minimization of the area is directly related to the given parameterization so that these regions depend on the labeling of the samples. However the regions based on the likelihood ratio (5) and (6) are invariant. The formulas in this section suggest to declare the sample with fewer observations as the first one. Indeed if ρ 4 > n/m, one gets larger Fisher information for ρ defined that way. For t , 0 ≤ t ≤ 1, let β(t ; k, ℓ) denote the density of a beta-distribution with parameters k and ℓ, and let B(t ; k, ℓ) be the corresponding cumulative distribution function (incomplete beta-function). To calculate the confidence coefficient corresponding to C˜ = (z , y) : z ≤ g (y, ρ), ˆ y− ≤ y ≤ y+ ,
y+ y−
g (y,ρ) ˆ
p(z , y) dz dy =
0
t+
B 1 − 1 + tg t−
1−t t
, ρˆ
−1
1 m+n
; , 2
2
n−1 m−1 − 1 β t; , dt , (14) 2
2
one needs only to evaluate one-dimensional integrals with t− = 1/(1 + y+ ) and t+ = 1/(1 + y− ). The exact Mood region in our situation, C˜ = {(z , y) : z ≤ a1 , a2 ≤ y ≤ a3 }, i.e.,
ˆ − ∆| ≤ (∆, ρ) : |∆
na1 (ρ + 1)θˆ , mρ/( ˆ na3 ) ≤ ρ ≤ mρ/( ˆ na2 ) .
(15)
A.L. Rukhin / Statistics and Probability Letters 119 (2016) 273–280
277
Fig. 1. The shape of five exact confidence regions: (6) (dotted line), (9) (dash-dotted line), (13) (line marked by points), Q1 (dashed line), (15) (solid line), ˆ = 0, ρˆ = θˆ = 1, α = 0.05. when n = 8, m = 13, ∆
The constants a1 , a2 , and a3 are chosen so that the coverage probability is 1 − α . The area of the Mood region is
4 na1 θˆ 3
1+
mρˆ
3/2
− 1+
na2
mρˆ
3/2
na3
,
which agrees with (12). To minimize this area for a fixed a1 with g (y, ρ) ˆ = a1 , y− = a2 , y+ = a3 , one has to find t+ = 1/(1 + a2 ), t− = 1/(1 + a3 ) which satisfy the equation,
a1 t+ 1 m+n n−1 m−1 (1 − t+ )5/2 B ; , − 1 β t+ ; , [n + (mρˆ − n)t+ ]1/2 1 + a1 t+ 2 2 2 2 5/2 (1 − t− ) a1 t− 1 m+n n−1 m−1 = B ; , − 1 β t− ; , . [n + (mρˆ − n)t− ]1/2 1 + a1 t− 2 2 2 2 For example, if n = 8, m = 13, α = 0.05, the minimum area Mood region is defined by a1 = 1.328, a2 = 0.508, a3 = 29.937; when ρˆ = 1, θˆ = 1 this area is 32.73. In contrast, the area of (6) is 30.88, that of (9) is 29.67 the minimum area is 23.29, and that of Q1 is 23.46. Thus the Mood region (15) is more voluminous than the confidence set (9) and the Bartlett corrected set (6). This fact ˆ = 0, ρˆ = 1, θˆ = 1, n = 8, m = is confirmed by Fig. 1, depicting four exact confidence regions discussed above when ∆ 13, α = 0.05. Clearly Q1 is the smallest region, with (6) and (9) quite close one to another, while the noticeably larger set (15) of a very different shape. The classical likelihood ratio based region (5) with the cut-off constant χα2 (2) is not shown here since it does not maintain the nominal confidence coefficient 1 − α unless n and m are quite large. The area of these regions when normalized by the area of the Mood set, is shown in Fig. 2 as a function of ρˆ , 0.1 < ρˆ < 5. The relative area of Q1 is quite close to that of the minimal area region (13); they both decrease as ρˆ increases while these of (6) and (9) are approaching one. The confidence ellipsoid for (∆, ρ) based on the Wald procedure,
ˆ − ∆)2 nm(ρˆ − ρ)2 (∆ + ≤ χα2 (2), 2(n + m)ρˆ 2 (ρˆ + 1)θˆ can be interpreted as a Taylor series approximation to (5). This set is an ellipse truncated so that ρ is positive. As the result, this region does not maintain the nominal confidence coefficient 1 − α . The score statistic based confidence set which uses the form of the inverse of the Fisher information matrix I (ψ),
ˆ − ∆)2 (∆ (n + m) + 2 (ρ + 1) θ˜ 2nmθ˜ 2
ˆ − ∆) 2 ρ 2 (∆ − mθ˜ (ρ + 1)2
2 ≤ χα2 (2),
is unbounded. The confidence regions based on a linear combination of 1 + Z + Y and Y are typically unbounded as well as the set obtained from the collection of all confidence intervals (3).
278
A.L. Rukhin / Statistics and Probability Letters 119 (2016) 273–280
Fig. 2. The normalized areas of four exact confidence regions (6) (dotted line), (9) (dash-dotted line), (13) (line marked by points), Q1 (dashed line).
These regions are not considered further. Instead we compare the exact confidence regions with the third order accurate modification of the likelihood ratio statistic, namely a two-parameter version of the Skovgaard statistic (2001) which extends the confidence interval for ∆ in Section 2. Straightforward calculations show that the observed Fisher information matrix J (ψ) when evaluated at the maximum ˆ , is equal to the classical Fisher information matrix I (ψ). Their common determinant is then likelihood estimator ψ ˆ = det J (ψ) ˆ = nm/[4ρˆ 3 θˆ 4 ]. The determinant of the lower 2 × 2 sub-matrix of the observed Fisher information det I (ψ) ˜ = det Iλλ (ψ) ˜ = (n + m)(ρ + 1)/[2ρ θ˜ 3 ]. It coincides corresponding to the nuisance parameter λ = (µ, θ ), is det Jλλ (ψ) 3 ˜4 ˜ ˜ with det Iλλ (ψ), while det I (ψ) = nm/(4ρ θ ). The needed matrix S is
1 0 1 S= θ˜ 1
0
0
0
mθˆ 2ρ
1+
mρˆ 2ρ 2
0
mθˆ
0
2
0
1
1
ρ
2ρ θ˜ , 0 nρ + mρˆ ρ θ˜
with
q=
ˆ −∆ m ∆ , (ρ + 1)θ˜ 2ρˆ
T ρˆ θˆ n 1 1 m ρˆ θˆ − 1 , 0, − −1 . + 2 θ˜ 2ρˆ ρ θ˜ ρ θ˜ θˆ
˜ is The score vector ℓ evaluated at the restricted maximum likelihood estimator ψ ℓ=
ˆ − ∆)2 ˆ − ∆ (∆ m ∆ , + 2 2ρ (ρ + 1)θ˜ 2(ρ + 1) θ˜
T ρˆ θˆ − 1 , 0, 0 . ρ θ˜
Therefore with adj(S ) denoting the adjoint matrix, nm ℓT adj(S )q = ρ(2ρˆ + ρ ρˆ + ρ 2 )(x − ∆)2 + m(ρ + 1)2 (ρˆ − ρ)(ρˆ θˆ − ρ θ˜ ) , 8ρ 4 (ρ + 1)2 ρˆ θ˜ 5 and
˜ = ℓT I −1 (ψ)ℓ
ˆ − ∆)2 (∆ ( n + m) ˆ − ∆)2 + m(ρ + 1)2 (ρˆ θˆ − ρ θ˜ )]2 . + [ρ 2 (∆ ˜ (ρ + 1)θ 2nmρ 2 (ρ + 1)4 θ˜ 2
The Skovgaard (2001) (Eq. (12)) statistic for a full two-parameter exponential model has the form V =
ℓT adj(S )q 4(ρ ρ) ˆ 3/2 (θˆ θ˜ )2 ℓT adj(S )q = . ˜ mnℓT I −1 (ψ)ℓ ˜ ˆ det I (ψ) ˜ ℓT I −1 (ψ)ℓ det I (ψ)
A.L. Rukhin / Statistics and Probability Letters 119 (2016) 273–280
279
ˆ = 0, ρˆ = 1, θˆ = 1. Fig. 3. The surface plot and contours for V − χ02.05 (2) when n = 8, m = 13, ∆
Here we have used the fact that the determinants of sub-blocks of observed information matrix corresponding to the ˜ = det Jλλ (ψ) ˆ . nuisance parameter λ = (µ, θ ) are equal, det Jλλ (ψ) All quantities entering in the formula for v have been evaluated so that one can get the final form of v , as well as that of the corresponding modified likelihood ratio test statistic, V = V (∆, ρ) = R2 + 2 log v , or its asymptotic equivalent R2 (1 + R−2 log v)2 . Using the statistic V and taking into account (5), one gets the confidence region for ∆ and ρ based on the third order accurate likelihood ratio test,
(∆, ρ) : (n + m − 2) log
θ˜ θˆ
ρ 2 + 2 log w ≤ χα (2) . + (m − 1) log ρˆ
(16)
Here
ˆ − ∆)2 m(ρ + 1)2 (ρˆ − ρ)(ρˆ θˆ − ρ θ˜ ) ρ(2ρˆ + ρ ρˆ + ρ 2 )(∆ + w= θˆ θˆ −1 ˆ − ∆)2 + m(ρ + 1)2 (ρˆ θˆ − ρ θ˜ )]2 ˆ − ∆)2 (n + m)[ρ 2 (∆ 2ρ 2 (ρ + 1)θ˜ (∆ + × θˆ 2 nm(ρ + 1)2 θˆ [(2n + m)nm−1 Y 2 + (nρˆ + m)Y + mn−1 (n + 2m)ρˆ 2 ]Z (nY − m)2 = + mρˆ + nY nm [n(2n + m)Y 2 + 2nmρˆ Y + m(n + 2m)ρˆ 2 ]Z 2 2(n + m)(ρˆ + 1)YZ (nY − m)2 −1 × + + , 2 (mρˆ + nY ) mρˆ + nY nm
with pivots Z and Y . The surface and contours of the function V (∆, ρ) are shown in Fig. 3. Although the coverage probability of (16) depends on the parameter ρ , this region maintains the nominal 95% confidence coefficient almost perfectly. However the area of (16) when ρˆ = 1 is 33.66 which exceeds that of (9) or that of Q1 . The shape of that region is shown in Fig. 4 for three values of ρˆ when n = 8, m = 13, θˆ = 1, α = 0.05. The larger is ρˆ , the larger is the volume of the region. The derived procedures can be used to test the null hypothesis H0 : ∆ = ∆0 , ρ = ρ0 . Depending on the test’s outcome one may decide to employ the statistic based on the given ρ0 , cf. Sprott and Farewell (1993). See also She et al. (2011) for further testing procedures. Bédart et al. (2007), and Fraser et al. (2009) give some results derived by using third order accurate tests for ∆ in the Behrens–Fisher problem One can entertain testing a more difficult hypothesis: ∆ = ∆0 , ρ ≥ ρ0 , where ρ0 represents the smallest permissible ratio of variances. For example, in applications discussed in Rukhin (2015), one of the variances can be expected to be (much) larger than another. A test can be derived from the intersection-union principle. Indeed, the likelihood ratio test of the null hypothesis H0 for a known ρ0 in our notation rejects this hypothesis when R(∆0 , ρ0 ) ≥ χα2 (2), and {minρ:ρ≥ρ0 R(∆0 , ρ) ≥ χα2 (2)} is the rejection region of the test corresponding to (5). The level of the resulting test is bounded above by α . Similar regions can be derived for other considered statistics like Q0 or Q1 . For instance that for (16) has the form {minρ:ρ≥ρ0 V (∆0 , ρ) ≥ χα2 (2)}. However none of these tests is practical as the functions R(∆0 , ρ) and V (∆0 , ρ) commonly are not unimodal with the rejection region which becomes a union of disjoint intervals.
280
A.L. Rukhin / Statistics and Probability Letters 119 (2016) 273–280
Fig. 4. Confidence region (16) for three values of ρˆ = 1/2: ρˆ = 1/2 (dotted line), ρˆ = 1 (dashed line), ρˆ = m/n = 1.625 (solid line).
4. Conclusions We summarize here four recommended confidence regions for ∆ and ρ .
• • • •
The Bartlett corrected likelihood ratio test statistic given in (6), The set in (9) based on statistic (11) whose two first moments coincide with that of χ 2 (2)-distribution, The set Q1 with the cut-off constant K found from (10) for a = m + n − 1, b = m + 2, The region based on the third-order optimal modified likelihood ratio test statistic W given by (16).
One can get their confidence coefficients from the formula (14) and the areas by using (12). The determination of the threshold constant K of the minimum area region (13) is not automatic. The exact Mood region which is given in (15) may be used if large values of ρ are not anticipated, but the Wald statistic or the score statistic do not present viable alternatives. Acknowledgment The author is grateful to the referee for helpful comments. References Arnold, B.C., Shavelle, R.M., 1998. Joint confidence sets for the mean and variance of a normal distribution. Amer. Statist. 52, 133–140. Bédart, M., Fraser, D.A.S., Wong, A.C.M., 2007. Higher accuracy for Bayesian and frequentist inference: large sample theory for small sample likelihood. Statist. Sci. 22, 301–121. Brazzale, A.R., Davison, A.C., Reid, N., 2007. Applied Asymptotics: Case Studies in Small Sample Statistics. Cambridge University Press, Cambridge, UK. Cordeiro, G., Cribari-Neto, F., 2014. An Introduction to Bartlett Correction and Bias Reduction. Springer, New York. Fraser, D.A.S., Wong, A., Sun, J., 2009. Three enigmatic examples and inference from likelihood. Canad. J. Statist. 37, 161–181. Frey, J., Marrero, O., Norton, D., 2009. Minimum area confidence sets for a normal distribution. J. Statist. Plann. Inference 139, 1023–1032. Lehmann, E.L., Romano, J.P., 2003. Testing Statistical Hypotheses, third ed. Springer, New York. Perng, S.K., Littell, R.C., 1976. A test of equality of two normal population means and variances. J. Amer. Statist. Assoc. 71, 968–971. Rukhin, A.L., 2015. Random effects model for bias estimation: higher order asymptotic inference. Stat 4, 130–139. She, Y., Wong, A., Zhou, X., 2011. Revisit the two sample t-test with a known ratio of variances. Open J. Statist. 1, 151–156. Skovgaard, I., 1996. An explicit large-deviation approximation to one-parameter tests. Bernoulli 2, 145–165. Skovgaard, I., 2001. Likelihood asymptotics. Scand. J. Stat. 28, 3–32. Sprott, D.A., Farewell, V.T., 1993. The difference between two normal means. Amer. Statist. 47, 126–128. Zhang, L., Xu, X., Chen, G., 2012. The exact likelihood ratio test for equality of two normal populations. Amer. Statist. 66, 180–184.
Contents lists available at ScienceDirect
Statistics and Probability Letters journal homepage: www.elsevier.com/locate/stapro
Confidence regions for comparison of two normal samples Andrew L. Rukhin Statistical Engineering Division, National Institute of Standards and Technology, 100 Bureau Drive, Gaithersburg, MD 20899, USA
article
info
Article history: Received 12 April 2016 Received in revised form 18 August 2016 Accepted 19 August 2016 Available online 30 August 2016
abstract The comparison of two normal populations is considered when one is interested not only in the difference between the means but also in the ratio of the variances. Several joint confidence regions for these two parameters are obtained. Published by Elsevier B.V.
Keywords: Bartlett correction Behrens–Fisher problem Higher order inference Minimum area region
1. Introduction The problem of two normal populations raises important statistical issues and is often encountered in practice when two different methods (instruments, treatments, materials, etc.) are to be compared. The classical monograph by Lehmann and Romano (2003) dedicates a whole section to this subject whose theoretical importance is due to the absence of a single pivotal quantity on which to base statistical inference in the classical Behrens–Fisher problem. Indeed there seems to be no unique universally accepted method either to test equality of the means or to derive a confidence interval covering their difference. However two samples even with the same mean can be quite dissimilar unless their variances are not too different. Therefore the variances must enter in a meaningful comparison calling for a hypothesis test about the value of their ratio or for a joint confidence region for the difference in the two means and this ratio. We review several such confidence regions in Section 3. These include the exact regions produced by the likelihood ratio and its Bartlett and two-moments based corrections, the minimum area set, as well as the Mood-type confidence set. One of the suggested procedures is obtained from the modern higher order asymptotic statistical inference (Brazzale et al., 2007). The situation when the variances ratio is known is briefly discussed in Section 2. Section 4 concludes the paper with some recommendations. 2. Comparing means and variances: known variances ratio Here we assume given two normal observables x¯ ∼ N (µ + ∆, θ ) and y¯ ∼ N (µ, η), as well as two statistics s2 and u2 which are independent of x¯ and y¯ . Their probability distributions are of the form θ χ 2 (n − 1)/(n − 1) and ηχ 2 (m − 1)/(m − 1). For example, s2 is an unbiased estimator of θ which is the variance of the first sample mean. Thus the reduction from two random samples to sufficient statistics (¯x, y¯ , s2 , u2 ) is made, θ = σ12 /n, η = σ22 /m, where σ12 , σ22 are the samples’ variances, and n, m the respective sample sizes.
E-mail address: [email protected]. http://dx.doi.org/10.1016/j.spl.2016.08.011 0167-7152/Published by Elsevier B.V.
274
A.L. Rukhin / Statistics and Probability Letters 119 (2016) 273–280
The log-likelihood function corresponding to x¯ , y¯ , s2 and u2 involves a four-dimensional parameter ψ = (∆, µ, θ , η). It has the form
(¯y − µ)2 + (m − 1)u2 + n log θ + m log η . (1) 2 θ η In this section we look at the simplest case when the ratio of two variances is given, ρ = η/θ . Then λ = (µ, θ ) is the nuisance parameter and the likelihood for ψ = (∆, λ) can be obtained from (1) by replacing η by ρθ . The maximum ˆ = x¯ − y¯ , and that of µ is µ likelihood estimator θˆ of θ is [(n − 1)s2 + (m − 1)u2 /ρ]/(m + n), while that of ∆ is ∆ ˆ = y¯ . The constrained maximum likelihood estimator of µ for fixed ∆, L(ψ) = −
µ ˜ =
1 (¯x − µ − ∆)2 + (n − 1)s2
x¯ − ∆ + y¯ /ρ 1 + 1/ρ
+
,
(2)
is the weighted mean of x¯ − ∆ and y¯ . The constrained maximum likelihood estimator θ˜ satisfies the relation, θ˜ = ˆ − ∆)2 , and the signed likelihood root statistic has the form θˆ + (n + m)−1 (ρ + 1)−1 (∆
ˆ − ∆) (n + m) log θ˜ /θˆ r = r (∆) = sign(∆
1/2
.
The statistic r is known to be approximately normally distributed although this approximation typically can be poor for 2 ˆ − ∆|2 ≤ [exp(zα/ small sample sizes. For example, the confidence interval {∆ : |r | ≤ zα/2 } = {∆ : |∆ 2 /(n + m)) − 1](n +
m)(ρ + 1)θˆ } is too short and does not maintain the nominal confidence coefficient well. A much more accurate normal approximation can be obtained via the Bartlett correction discussed in more detail in Section 3. Another possibility is to employ a third order accurate statistic (Section 8.5 of Brazzale et al., 2007). One of these, the
ˆ − ∆)θˆ / (ρ + 1)θ˜ 3 , provides third order improvement on r, with the confidence interval, Skovgaard (1996) statistic, v = (∆ {∆ : |r + log(v/r )| ≤ zα/2 }.
(3)
This interval turns out to be very close to
ˆ ± tα/2 (n + m − 2) (n + m)(n + m − 2)−1 (ρ + 1)θˆ , ∆ which is the usual t-interval. Here zα/2 and tα/2 (n + m − 2) denote the familiar critical points. The tests of the hypothesis: ∆ = 0, ρ = 1 were studied by Perng and Littell (1976), and more recently by Zhang et al. (2012). The goal of this work is to construct a joint confidence region for (∆, ρ). 3. Comparing means and variances: unknown variances ratio With ψ = (∆, ρ; θ , µ) = (∆, ρ; λ), the score vector ℓ = ℓ(ψ; x¯ , y¯ , s2 , u2 ), has the form, x¯ − µ − ∆ y¯ − µ , + , θ θ η (¯x − µ − ∆)2 + (n − 1)s2 n (¯y − µ)2 + (m − 1)u2 m T − , − . 2θ 2 2θ 2η 2 2η ˆ = x¯ − y¯ , the maximum likelihood estimator of µ is y¯ , ρˆ = n(m − The maximum likelihood estimator of ∆ remains ∆ 1)u2 /[m(n − 1)s2 ], θˆ = [(n − 1)s2 + (m − 1)u2 /ρ]/( ˆ m + n) = (n − 1)s2 /n. The Fisher information matrix has the form 1 1
x¯ − µ − ∆
θ 0 I (ψ) = 1 θ 0
0
m 2ρ 2 0 m 2ρθ
0
0
ρ+1 ρθ 0
θ
2ρθ . 0 n + m m
2θ 2
Now the parameter of interest has two components (∆, ρ); the constrained maximum likelihood estimators of µ and θ for fixed ∆, ρ determine the likelihood ratio test of the hypothesis postulating these null values. The first of these estimators µ ˜ has already been found in (2). The estimator θ˜ of θ can be written as
(¯x − y¯ − ∆)2 /(ρ + 1) + (n − 1)s2 + (m − 1)u2 /ρ n+m ˆ − ∆)2 (∆ (mρˆ + nρ)θˆ = + , (n + m)(ρ + 1) (n + m)ρ
θ˜ =
(4)
A.L. Rukhin / Statistics and Probability Letters 119 (2016) 273–280
275
leading to the likelihood ratio (Neyman–Pearson) statistic,
ˆ , ρ, R = R(∆, ρ) = 2[L(∆ ˆ y¯ , θˆ ) − L(∆, ρ, µ, ˜ θ˜ )] = (n + m) log(θ˜ /θˆ ) + m log(ρ/ρ). ˆ Thus the joint confidence region for ∆ and ρ based on R is
(∆, ρ) : (n + m) log(θ˜ /θˆ ) + m log(ρ/ρ) ˆ ≤ χα2 (2) ˆ − ∆)2 (∆ m(ρˆ − ρ) ρ 2 = (∆, ρ) : (n + m) log 1 + + ≤ χα (2) . + m log (n + m)ρ ρˆ (n + m)(ρ + 1)θˆ
(5)
Here χα2 (2) is the α -critical point of a χ 2 -distribution with 2 degrees of freedom, use of which is recommended by the first order asymptotic theory. However unless both n and m are large, the cut-off constant χα2 (2) of the region (5) is too small leading to insufficient coverage probability. ˆ and Y = mρ/( ˆ − ∆)2 /[n(ρ + 1)θ] Clearly R depends only on Z = (∆ ˆ nρ), which is a two-dimensional pivot. The joint ˆ − ∆)2 /[(ρ + 1)θ], nθˆ /θ , and distribution of (Z , Y ) can be obtained from that of three independent χ 2 -random variables, (∆ ˆ mρˆ θ/[ρθ ] which have degrees of freedom 1, n − 1 and m − 1 respectively. This distribution provides the joint density of Z and Y , p(z , y) =
Γ ((n + m − 1)/2)y(m−3)/2 . Γ (1/2)Γ ((n − 1)/2)Γ ((m − 1)/2)z 1/2 (1 + z + y)(n+m−1)/2
One can considerably improve accuracy of the χ 2 -approximation in (5) by using the Bartlett correction (Cordeiro and Cribari-Neto, 2014). In our situation this correction is simplified by the fact that the distribution of the likelihood ratio statistic does not depend on unknown parameters,
n−1 m−1 − nΨ − mΨ − (n + m) log(n + m) + n log(n) + m log(m) 2 2 2 11(n2 + nm + m2 ) 1 =2 1+ +O . 12nm(n + m) min(n, m)2
ER = (m + n)Ψ
n+m−1
Here Ψ denotes the well-known log-derivative of the gamma-function. Thus the Bartlett corrected likelihood ratio statistic,
−1 11(n2 + nm + m2 ) 11(n2 + nm + m2 ) ≈R 1+ ≈R 1− , ER 12nm(n + m) 12nm(n + m)
2R
admits much better χ 2 -approximation than R, so that the nominal confidence coefficient of the confidence region
(∆, ρ) : R ≤ χα2 (2)ER/2
(6)
is considerably closer to 1 − α than that of (5). For example, the variance of the Bartlett corrected likelihood ratio statistic, which can be obtained from the formula Var(R) = n Ψ 2
′
n−1
2
+m Ψ 2
′
m−1
2
− ( m + n) Ψ 2
′
n+m−1
2
,
coincides with Var(χ 2 (2)) = 4 up to the order O min(n, m)−2 .
More general confidence sets (∆, ρ) : Q ≤ K , where K is a threshold constant and Q is a linear combination of log(1 + Z + Y ) and log(Y ), Q = a log(1 + Z + Y ) − b log(Y ) − c ;
a.b > 0,
are worthy of consideration. For example, Q = R, when a = m + n, b = m, c = (n + m) log(n + m)− n log(n)− m log(m). The mean EQ and the variance Var(Q ) admit explicit parameter-free formulas. A possible choice of constants is a = a0 , b = b0 which equate the first two moments of Q to those of the χ 2 (2)-distribution,
2 = EQ = a Ψ
n+m−1 2
− log
n+m
2
n m m−1 n−1 − ( a − b) Ψ − log −b Ψ − log 2
2
2
(7)
2
and 4 = Var(Q ) = (a − b)2 Ψ ′
n−1 2
+ b2 Ψ ′
m−1 2
− a2 Ψ ′
n+m−1 2
.
(8)
276
A.L. Rukhin / Statistics and Probability Letters 119 (2016) 273–280
This choice leads to the region
n ρ (x − ∆)2 mρˆ n+m + b0 log ≤ χα2 (2) + a0 log + b0 log . (∆, ρ) : Q0 = a0 log 1 + + nρ ρˆ n m n(ρ + 1)θˆ
(9)
When n = m, one gets b0 = a0 /2, which is determined by the first moment identity so that the second equation is not needed. In this situation Q0 coincides with the Bartlett corrected likelihood ratio statistic above. When min(n, m) ≥ 5, there are two positive solutions a0 to the simultaneous equations above; the largest root is recommended if m > n. √ The linear transform of Q , 2[(Q − EQ )/ Var(Q ) + 1] admits a much better χ 2 (2)-approximation than Q itself. This fact gives an approximate formula for the threshold constant
K = K (Q ) = EQ + (χα2 (2)/2 − 1) Var(Q ),
(10)
where EQ and Var(Q ) are found from (7), (8). Under this choice the confidence coefficient of the region
(∆, ρ) : a log(1 + Z + Y ) − b log(Y ) ≤ K
(11)
is close to 1 − α . The confidence sets in the estimation problem for a normal mean and variance were studied by Arnold and Shavelle (1998). Frey et al. (2009) discuss the minimum area regions for these parameters, whose counterparts can be derived in our problem. Indeed if C is (1−α)-confidence set for (∆, ρ), let C˜ denote the set of (z , y)-values such that (∆, ρ) ∈ C if and only if (z , y) ∈ C˜ . The formula for the area of C can be obtained from the Jacobian J = J (z , y) of the transformation (∆, ρ) → (z , y),
θˆ (y + mρ/ ˆ n)/(nzy5 ). For all confidence sets C considered here, C˜ = (z , y) : z ≤ g (y, ρ), ˆ y− ≤ y ≤ y+ with a
J = mρˆ
positive function g and positive y− , y+ . The area of such a region depends on ρˆ (and on θˆ ),
d∆ dρ =
C
C˜
J (z , y) dz dy =
2mρˆ
θˆ
n
y+
g (y, ρ)( ˆ ny + mρ) ˆ
1/2
dy.
y5
y−
(12)
For (11), g (y, ρ) ˆ = g (y) = Ayκ − 1 − y, A = exp(K /a), κ = b/a and y− , y+ are the roots of the equation g (y) = 0. Since C˜ p(z , y) dz dy = 1 − α , the fundamental Neyman–Pearson lemma gives the form of the minimum area region,
C˜ = {(z , y) : J (z , y) ≤ K ′ p(z , y)}, i.e.
ˆ − ∆)2 (∆ mρˆ C = (∆, ρ) : (n + m − 1) log 1 + + ˆ nρ n(ρ + 1)θ
+ (m + 1) log
nρ mρˆ
+ log(1 + ρ) ≤ K .
(13)
When ρˆ is small, C has approximately the form (11) with a = m + n − 1, b = m + 1, for large values of ρˆ it has the similar form with a = m + n − 1, b = m + 2. In these two cases the threshold constant K can be found from (10) by using the identities (7) and (8). If ρˆ = n/m, the explicit form of these two moments follows from the fact that Z /(1 + Z + Y ) and 1/(1 + Y ) are independent and beta-distributed with parameters 1/2, (m + n)/2 − 1 and (n − 1)/2, (m − 1)/2 respectively, but in the general case to determine the constant K after the confidence coefficient 1 − α , one has to use numerical integration or Monte Carlo evaluations. When max(n, m) ≥ 8, min(n, m) ≥ 3, numerical results favor the region (11) with a = m + n − 1, b = m + 2, say, Q1 = (m + n − 1) log(1 + Z + Y ) − (m + 2) log(Y ). The minimization of the area is directly related to the given parameterization so that these regions depend on the labeling of the samples. However the regions based on the likelihood ratio (5) and (6) are invariant. The formulas in this section suggest to declare the sample with fewer observations as the first one. Indeed if ρ 4 > n/m, one gets larger Fisher information for ρ defined that way. For t , 0 ≤ t ≤ 1, let β(t ; k, ℓ) denote the density of a beta-distribution with parameters k and ℓ, and let B(t ; k, ℓ) be the corresponding cumulative distribution function (incomplete beta-function). To calculate the confidence coefficient corresponding to C˜ = (z , y) : z ≤ g (y, ρ), ˆ y− ≤ y ≤ y+ ,
y+ y−
g (y,ρ) ˆ
p(z , y) dz dy =
0
t+
B 1 − 1 + tg t−
1−t t
, ρˆ
−1
1 m+n
; , 2
2
n−1 m−1 − 1 β t; , dt , (14) 2
2
one needs only to evaluate one-dimensional integrals with t− = 1/(1 + y+ ) and t+ = 1/(1 + y− ). The exact Mood region in our situation, C˜ = {(z , y) : z ≤ a1 , a2 ≤ y ≤ a3 }, i.e.,
ˆ − ∆| ≤ (∆, ρ) : |∆
na1 (ρ + 1)θˆ , mρ/( ˆ na3 ) ≤ ρ ≤ mρ/( ˆ na2 ) .
(15)
A.L. Rukhin / Statistics and Probability Letters 119 (2016) 273–280
277
Fig. 1. The shape of five exact confidence regions: (6) (dotted line), (9) (dash-dotted line), (13) (line marked by points), Q1 (dashed line), (15) (solid line), ˆ = 0, ρˆ = θˆ = 1, α = 0.05. when n = 8, m = 13, ∆
The constants a1 , a2 , and a3 are chosen so that the coverage probability is 1 − α . The area of the Mood region is
4 na1 θˆ 3
1+
mρˆ
3/2
− 1+
na2
mρˆ
3/2
na3
,
which agrees with (12). To minimize this area for a fixed a1 with g (y, ρ) ˆ = a1 , y− = a2 , y+ = a3 , one has to find t+ = 1/(1 + a2 ), t− = 1/(1 + a3 ) which satisfy the equation,
a1 t+ 1 m+n n−1 m−1 (1 − t+ )5/2 B ; , − 1 β t+ ; , [n + (mρˆ − n)t+ ]1/2 1 + a1 t+ 2 2 2 2 5/2 (1 − t− ) a1 t− 1 m+n n−1 m−1 = B ; , − 1 β t− ; , . [n + (mρˆ − n)t− ]1/2 1 + a1 t− 2 2 2 2 For example, if n = 8, m = 13, α = 0.05, the minimum area Mood region is defined by a1 = 1.328, a2 = 0.508, a3 = 29.937; when ρˆ = 1, θˆ = 1 this area is 32.73. In contrast, the area of (6) is 30.88, that of (9) is 29.67 the minimum area is 23.29, and that of Q1 is 23.46. Thus the Mood region (15) is more voluminous than the confidence set (9) and the Bartlett corrected set (6). This fact ˆ = 0, ρˆ = 1, θˆ = 1, n = 8, m = is confirmed by Fig. 1, depicting four exact confidence regions discussed above when ∆ 13, α = 0.05. Clearly Q1 is the smallest region, with (6) and (9) quite close one to another, while the noticeably larger set (15) of a very different shape. The classical likelihood ratio based region (5) with the cut-off constant χα2 (2) is not shown here since it does not maintain the nominal confidence coefficient 1 − α unless n and m are quite large. The area of these regions when normalized by the area of the Mood set, is shown in Fig. 2 as a function of ρˆ , 0.1 < ρˆ < 5. The relative area of Q1 is quite close to that of the minimal area region (13); they both decrease as ρˆ increases while these of (6) and (9) are approaching one. The confidence ellipsoid for (∆, ρ) based on the Wald procedure,
ˆ − ∆)2 nm(ρˆ − ρ)2 (∆ + ≤ χα2 (2), 2(n + m)ρˆ 2 (ρˆ + 1)θˆ can be interpreted as a Taylor series approximation to (5). This set is an ellipse truncated so that ρ is positive. As the result, this region does not maintain the nominal confidence coefficient 1 − α . The score statistic based confidence set which uses the form of the inverse of the Fisher information matrix I (ψ),
ˆ − ∆)2 (∆ (n + m) + 2 (ρ + 1) θ˜ 2nmθ˜ 2
ˆ − ∆) 2 ρ 2 (∆ − mθ˜ (ρ + 1)2
2 ≤ χα2 (2),
is unbounded. The confidence regions based on a linear combination of 1 + Z + Y and Y are typically unbounded as well as the set obtained from the collection of all confidence intervals (3).
278
A.L. Rukhin / Statistics and Probability Letters 119 (2016) 273–280
Fig. 2. The normalized areas of four exact confidence regions (6) (dotted line), (9) (dash-dotted line), (13) (line marked by points), Q1 (dashed line).
These regions are not considered further. Instead we compare the exact confidence regions with the third order accurate modification of the likelihood ratio statistic, namely a two-parameter version of the Skovgaard statistic (2001) which extends the confidence interval for ∆ in Section 2. Straightforward calculations show that the observed Fisher information matrix J (ψ) when evaluated at the maximum ˆ , is equal to the classical Fisher information matrix I (ψ). Their common determinant is then likelihood estimator ψ ˆ = det J (ψ) ˆ = nm/[4ρˆ 3 θˆ 4 ]. The determinant of the lower 2 × 2 sub-matrix of the observed Fisher information det I (ψ) ˜ = det Iλλ (ψ) ˜ = (n + m)(ρ + 1)/[2ρ θ˜ 3 ]. It coincides corresponding to the nuisance parameter λ = (µ, θ ), is det Jλλ (ψ) 3 ˜4 ˜ ˜ with det Iλλ (ψ), while det I (ψ) = nm/(4ρ θ ). The needed matrix S is
1 0 1 S= θ˜ 1
0
0
0
mθˆ 2ρ
1+
mρˆ 2ρ 2
0
mθˆ
0
2
0
1
1
ρ
2ρ θ˜ , 0 nρ + mρˆ ρ θ˜
with
q=
ˆ −∆ m ∆ , (ρ + 1)θ˜ 2ρˆ
T ρˆ θˆ n 1 1 m ρˆ θˆ − 1 , 0, − −1 . + 2 θ˜ 2ρˆ ρ θ˜ ρ θ˜ θˆ
˜ is The score vector ℓ evaluated at the restricted maximum likelihood estimator ψ ℓ=
ˆ − ∆)2 ˆ − ∆ (∆ m ∆ , + 2 2ρ (ρ + 1)θ˜ 2(ρ + 1) θ˜
T ρˆ θˆ − 1 , 0, 0 . ρ θ˜
Therefore with adj(S ) denoting the adjoint matrix, nm ℓT adj(S )q = ρ(2ρˆ + ρ ρˆ + ρ 2 )(x − ∆)2 + m(ρ + 1)2 (ρˆ − ρ)(ρˆ θˆ − ρ θ˜ ) , 8ρ 4 (ρ + 1)2 ρˆ θ˜ 5 and
˜ = ℓT I −1 (ψ)ℓ
ˆ − ∆)2 (∆ ( n + m) ˆ − ∆)2 + m(ρ + 1)2 (ρˆ θˆ − ρ θ˜ )]2 . + [ρ 2 (∆ ˜ (ρ + 1)θ 2nmρ 2 (ρ + 1)4 θ˜ 2
The Skovgaard (2001) (Eq. (12)) statistic for a full two-parameter exponential model has the form V =
ℓT adj(S )q 4(ρ ρ) ˆ 3/2 (θˆ θ˜ )2 ℓT adj(S )q = . ˜ mnℓT I −1 (ψ)ℓ ˜ ˆ det I (ψ) ˜ ℓT I −1 (ψ)ℓ det I (ψ)
A.L. Rukhin / Statistics and Probability Letters 119 (2016) 273–280
279
ˆ = 0, ρˆ = 1, θˆ = 1. Fig. 3. The surface plot and contours for V − χ02.05 (2) when n = 8, m = 13, ∆
Here we have used the fact that the determinants of sub-blocks of observed information matrix corresponding to the ˜ = det Jλλ (ψ) ˆ . nuisance parameter λ = (µ, θ ) are equal, det Jλλ (ψ) All quantities entering in the formula for v have been evaluated so that one can get the final form of v , as well as that of the corresponding modified likelihood ratio test statistic, V = V (∆, ρ) = R2 + 2 log v , or its asymptotic equivalent R2 (1 + R−2 log v)2 . Using the statistic V and taking into account (5), one gets the confidence region for ∆ and ρ based on the third order accurate likelihood ratio test,
(∆, ρ) : (n + m − 2) log
θ˜ θˆ
ρ 2 + 2 log w ≤ χα (2) . + (m − 1) log ρˆ
(16)
Here
ˆ − ∆)2 m(ρ + 1)2 (ρˆ − ρ)(ρˆ θˆ − ρ θ˜ ) ρ(2ρˆ + ρ ρˆ + ρ 2 )(∆ + w= θˆ θˆ −1 ˆ − ∆)2 + m(ρ + 1)2 (ρˆ θˆ − ρ θ˜ )]2 ˆ − ∆)2 (n + m)[ρ 2 (∆ 2ρ 2 (ρ + 1)θ˜ (∆ + × θˆ 2 nm(ρ + 1)2 θˆ [(2n + m)nm−1 Y 2 + (nρˆ + m)Y + mn−1 (n + 2m)ρˆ 2 ]Z (nY − m)2 = + mρˆ + nY nm [n(2n + m)Y 2 + 2nmρˆ Y + m(n + 2m)ρˆ 2 ]Z 2 2(n + m)(ρˆ + 1)YZ (nY − m)2 −1 × + + , 2 (mρˆ + nY ) mρˆ + nY nm
with pivots Z and Y . The surface and contours of the function V (∆, ρ) are shown in Fig. 3. Although the coverage probability of (16) depends on the parameter ρ , this region maintains the nominal 95% confidence coefficient almost perfectly. However the area of (16) when ρˆ = 1 is 33.66 which exceeds that of (9) or that of Q1 . The shape of that region is shown in Fig. 4 for three values of ρˆ when n = 8, m = 13, θˆ = 1, α = 0.05. The larger is ρˆ , the larger is the volume of the region. The derived procedures can be used to test the null hypothesis H0 : ∆ = ∆0 , ρ = ρ0 . Depending on the test’s outcome one may decide to employ the statistic based on the given ρ0 , cf. Sprott and Farewell (1993). See also She et al. (2011) for further testing procedures. Bédart et al. (2007), and Fraser et al. (2009) give some results derived by using third order accurate tests for ∆ in the Behrens–Fisher problem One can entertain testing a more difficult hypothesis: ∆ = ∆0 , ρ ≥ ρ0 , where ρ0 represents the smallest permissible ratio of variances. For example, in applications discussed in Rukhin (2015), one of the variances can be expected to be (much) larger than another. A test can be derived from the intersection-union principle. Indeed, the likelihood ratio test of the null hypothesis H0 for a known ρ0 in our notation rejects this hypothesis when R(∆0 , ρ0 ) ≥ χα2 (2), and {minρ:ρ≥ρ0 R(∆0 , ρ) ≥ χα2 (2)} is the rejection region of the test corresponding to (5). The level of the resulting test is bounded above by α . Similar regions can be derived for other considered statistics like Q0 or Q1 . For instance that for (16) has the form {minρ:ρ≥ρ0 V (∆0 , ρ) ≥ χα2 (2)}. However none of these tests is practical as the functions R(∆0 , ρ) and V (∆0 , ρ) commonly are not unimodal with the rejection region which becomes a union of disjoint intervals.
280
A.L. Rukhin / Statistics and Probability Letters 119 (2016) 273–280
Fig. 4. Confidence region (16) for three values of ρˆ = 1/2: ρˆ = 1/2 (dotted line), ρˆ = 1 (dashed line), ρˆ = m/n = 1.625 (solid line).
4. Conclusions We summarize here four recommended confidence regions for ∆ and ρ .
• • • •
The Bartlett corrected likelihood ratio test statistic given in (6), The set in (9) based on statistic (11) whose two first moments coincide with that of χ 2 (2)-distribution, The set Q1 with the cut-off constant K found from (10) for a = m + n − 1, b = m + 2, The region based on the third-order optimal modified likelihood ratio test statistic W given by (16).
One can get their confidence coefficients from the formula (14) and the areas by using (12). The determination of the threshold constant K of the minimum area region (13) is not automatic. The exact Mood region which is given in (15) may be used if large values of ρ are not anticipated, but the Wald statistic or the score statistic do not present viable alternatives. Acknowledgment The author is grateful to the referee for helpful comments. References Arnold, B.C., Shavelle, R.M., 1998. Joint confidence sets for the mean and variance of a normal distribution. Amer. Statist. 52, 133–140. Bédart, M., Fraser, D.A.S., Wong, A.C.M., 2007. Higher accuracy for Bayesian and frequentist inference: large sample theory for small sample likelihood. Statist. Sci. 22, 301–121. Brazzale, A.R., Davison, A.C., Reid, N., 2007. Applied Asymptotics: Case Studies in Small Sample Statistics. Cambridge University Press, Cambridge, UK. Cordeiro, G., Cribari-Neto, F., 2014. An Introduction to Bartlett Correction and Bias Reduction. Springer, New York. Fraser, D.A.S., Wong, A., Sun, J., 2009. Three enigmatic examples and inference from likelihood. Canad. J. Statist. 37, 161–181. Frey, J., Marrero, O., Norton, D., 2009. Minimum area confidence sets for a normal distribution. J. Statist. Plann. Inference 139, 1023–1032. Lehmann, E.L., Romano, J.P., 2003. Testing Statistical Hypotheses, third ed. Springer, New York. Perng, S.K., Littell, R.C., 1976. A test of equality of two normal population means and variances. J. Amer. Statist. Assoc. 71, 968–971. Rukhin, A.L., 2015. Random effects model for bias estimation: higher order asymptotic inference. Stat 4, 130–139. She, Y., Wong, A., Zhou, X., 2011. Revisit the two sample t-test with a known ratio of variances. Open J. Statist. 1, 151–156. Skovgaard, I., 1996. An explicit large-deviation approximation to one-parameter tests. Bernoulli 2, 145–165. Skovgaard, I., 2001. Likelihood asymptotics. Scand. J. Stat. 28, 3–32. Sprott, D.A., Farewell, V.T., 1993. The difference between two normal means. Amer. Statist. 47, 126–128. Zhang, L., Xu, X., Chen, G., 2012. The exact likelihood ratio test for equality of two normal populations. Amer. Statist. 66, 180–184.