Testing on the common mean of several normal distributions

Abstract

Point estimation of the common mean of several normal distributions with unknown and possibly unequal variances has attracted the attention of many researchers over the last five decades. Relatively less attention has been paid to the hypothesis testing problem, presumably due to the complicated sampling distribution(s) of the test statistics(s) involved. Taking advantage of the computational resources available nowadays there has been a renewed interest in this problem, and a few test procedures have been proposed lately including those based on the generalized $p$-value approach. In this paper we propose three new tests based on the famous Graybill–Deal estimator (GDE) as well as the maximum likelihood estimator (MLE) of the common mean, and these test procedures appear to work as good as (if not better than) the existing test methods. The two tests based on the GDE use respectively a first order unbiased variance estimate proposed by Sinha [Sinha, B.K., 1985. Unbiased estimation of the variance of the Graybill–Deal estimator of the common mean of several normal populations. The Canadian Journal of Statistics 13 (3), 243–247], as well as the little known exact unbiased variance estimator proposed by Nikulin and Voinov [Nikulin, M.S., Voinov, V.G., 1995. On the problem of the means of weighted normal populations. Qüestiió (Quaderns d’Estadistica, Sistemes, Informatica i Investigació Operativa) 19 (1–3), 93–106] (after we’ve fixed a small mistake in the final expression). On the other hand, the MLE, which doesn’t have a closed expression, uses a parametric bootstrap method proposed by Pal, Lim and Ling [Pal, N., Lim, W.K., Ling, C.H., 2007b. A computational approach to statistical inferences. Journal of Applied Probability & Statistics 2 (1), 13–35]. The extensive simulation results presented in this paper complement the recent studies undertaken by Krishnamoorthy and Lu [Krishnamoorthy, K., Lu, Y., 2003. Inferences on the common mean of several normal populations based on the generalized variable method. Biometrics 59, 237–247], and Lin and Lee [Lin, S.H., Lee, J.C., 2005. Generalized inferences on the common mean of several normal populations. Journal of Statistical Planning and Inference 134, 568–582].

Introduction

One of the oldest and most interesting problems in statistical sciences is the inference on a common mean of several normal distributions with unknown and possibly unequal variances. In literature, this is also known as Meta-Analysis where data from multiple sources are combined or integrated with a common objective. This problem arises in situations where different instruments or methods are used repeatedly to measure, say, blood alcohol level, or lead content in gasoline. One wishes to use multiple datasets, assuming normality, for an improved inference of the common mean, rather than relying on individual samples. As an application of the Meta-Analysis in clinical trials see Vazquez et al. (2007).

To formulate the present problem, assume that we have iid observations $X_{i1},\dots ,X_{i{n}_i}$ from $N(\mu ,{\sigma }_i^2),1\le i\le k$, where all parameters are assumed to be unknown. Define ${\bar{X}}_i$ and $S_i$ as

$${\bar{X}}_i=\sum _{j=1}^{n_i}{X}_{ij}/n_i\text{,}{S}_i=\sum _{j=1}^{n_i}{(X_{ij}-{\bar{X}}_i)}^2\text{;}$$

where ${\bar{X}}_i\sim N(\mu ,{\sigma }_i^2/n_i),S_i\sim {\sigma }_i^2{\chi }_{(n_i-1)}^2(1\le i\le k)$, and these statistics are all mutually independent. Throughout this paper it is assumed that $n_i\ge 2(1\le i\le k)$ unless mentioned otherwise.

Inference on the common mean $\mu$ has its genesis in a balanced incomplete block design (BIBD) with uncorrelated random block effects. For the $l$th treatment effect (say, ${\tau }_{\ell }$) one has two point estimates — namely, the intra-block estimate and the interblock estimate (say, ${\hat{\tau }}_{\ell }$ and ${\hat{\tau }}_{\ell }^{\ast }$ respectively). Under the usual design assumptions, ${\hat{\tau }}_{\ell }$ and ${\hat{\tau }}_{\ell }^{\ast }$ are independent, have normal distributions with common mean ${\tau }_{\ell }$, but with unknown and possibly unequal variances. Hence we have a common mean problem with $k=2$ (see also Montgomery (1997), page-216–218). Apart from estimating the $l$th treatment effect, one may wish to verify whether it is 0 or not, and thereby getting into a hypothesis testing problem. Zacks, 1966, Zacks, 1970, who contributed significantly to the theoretical foundation of the above mentioned common mean problem, described how he got interested in this problem by an application in soil science (see Kempthorne et al. (1991)).

Point estimation of $\mu$ has drawn the attention of many researchers over the past five decades. When ${\sigma }_i$’s are all known, the optimal estimator (MLE, BLUE as well as UMVUE) is

$$\hat{\mu }=\sum _{i=1}^k(n_i/{\sigma }_i^2){\bar{X}}_i/\sum _{i=1}^k(n_i/{\sigma }_i^2)\text{.}$$

But when ${\sigma }_i$’s are all unknown, one encounters the real challenge to propose an efficient estimator of $\mu$ combining the individual sample means ${\bar{X}}_i$’s.

Note that the minimal sufficient statistic $({\bar{X}}_i,S_i,i=1,2,\dots ,k)$ is not complete. As a result, one can not get the UMVUE (if it exists) using the standard Rao–Blackwell theorem on an unbiased estimator for estimating $\mu$.

In our present problem, where ${\sigma }_i$’s are all unknown, the most appealing unbiased estimator of $\mu$ has been the Graybill–Deal estimator (GDE) given as

$${\hat{\mu }}_{\mathrm{GDE}}=\sum _{i=1}^k(n_i/s_i^2){\bar{X}}_i/\sum _{i=1}^k(n_i/s_i^2)\text{,}$$

where $s_i^2=S_i/(n_i-1),1\le i\le k$. Graybill and Deal (1959) obtained conditions on the $n_i$’s for $k=2$ only for which ${\hat{\mu }}_{\mathrm{GDE}}$ has a smaller variance than each ${\bar{X}}_i(i=1,2)$.

It is important to remember that for estimating $\mu$ the above ${\hat{\mu }}_{\mathrm{GDE}}$ is not the MLE. The MLE, which does not have any closed expression, is

$${\hat{\mu }}_{\mathrm{MLE}}=\sum _{i=1}^k(n_i/{\hat{\sigma }}_{i\left(\mathrm{MLE}\right)}^2){\bar{X}}_i/\sum _{i=1}^k(n_i/{\hat{\sigma }}_{i\left(\mathrm{MLE}\right)}^2)\text{,}$$

where ${\hat{\sigma }}_{i\left(\mathrm{MLE}\right)}^2$ is the MLE of ${\sigma }_i^2(1\le i\le k)$ obtained by solving the following system of equations:

$${\hat{\sigma }}_{i\left(\mathrm{MLE}\right)}^2=(S_i/n_i)+{[\{\sum _{j=1}^k{n}_j({\bar{X}}_i-{\bar{X}}_j)/{\hat{\sigma }}_{j\left(MLE\right)}^2\}/\{\sum _{j=1}^k{n}_j/{\hat{\sigma }}_{j\left(MLE\right)}^2\}]}^2\text{,}i=1,2,\dots ,k\text{.}$$

It is because of this complicated expression that the MLE has dampened the interest of many researchers. In a recent paper Pal et al. (2007a) compared the MLE with the GDE for $k=2$, and observed that the MLE can be advantageous in a heavily unbalanced case. For a comprehensive review on the point estimation of $\mu$ till 1990 one can refer to Kubokawa (1991).

The main objective of this paper is to consider testing

$$H_0:\mu ={\mu }_0\text{ vs. }{H}_a:\mu \ne {\mu }_0$$

based on the statistics in (1.1). For a better understanding, we may occasionally focus on the special case of $k=2$. This special case is easy to understand in the context of BIBD mentioned earlier, but can be extended suitably for general $k$. Cohen and Sackrowitz (1977) showed that each individual $t$-test based only on one sample is admissible (which has the highest power at some boundary region of the parameter space). This is somewhat surprising in the light of inadmissibility results obtained in point estimation where, as expected, the combined data is supposed to help in inferences pertaining to the common mean. For the special case $k=2$ and $n_1=n_2=n$, Cohen and Sackrowitz (1984) considered several tests, including a normal approximate test based on ${\hat{\mu }}_{\mathrm{GDE}}$ in (1.3). Their study of power comparison showed that there is no outright winner among the tests considered. It is worth mentioning that the variance of ${\hat{\mu }}_{\mathrm{GDE}}$, which is

$$\mathrm{Var}\left({\hat{\mu }}_{\mathrm{GDE}}\right)=E[\sum _{i=1}^k{n}_i{\sigma }_i^2/s_i^4/{(\sum _{i=1}^k{n}_i/s_i^2)}^2]\text{,}$$

does not have a closed expression. As a result, it is not easy to get an estimate of the variance of ${\hat{\mu }}_{\mathrm{GDE}}$. Sinha (1985), using a special case of Haff’s (1979) Wishart identity, provided a way to obtain an approximate unbiased estimate of $\mathrm{Var}\left({\hat{\mu }}_{\mathrm{GDE}}\right)$ up to any desired order. The first order unbiased estimate of $\mathrm{Var}\left({\hat{\mu }}_{\mathrm{GDE}}\right)$ thus obtained, and hence denoted by ${\widehat{\mathrm{Var}}}_{\left(1\right)}\left({\hat{\mu }}_{\mathrm{GDE}}\right)$, is given as

$${\widehat{\mathrm{Var}}}_{\left(1\right)}\left({\hat{\mu }}_{\mathrm{GDE}}\right)={(\sum _{i=1}^k{n}_i/s_i^2)}^{-1}[1+4\sum _{i=1}^k{(n_i+1)}^{-1}(n_i/s_i^2){\{\sum _{i=1}^k(n_i/s_i^2)-{(n_i/s_i^2)}^2/{(\sum _{i=1}^k{n}_i/s_i^2)}^2\}}^{-1}]\text{,}$$

which is comparable to the approximation due to Meier (1953). The above result is helpful because the studentized version $(({\hat{\mu }}_{\mathrm{GDE}}-\mu )/{\left\{\widehat{\mathrm{Var}}\left({\hat{\mu }}_{\mathrm{GDE}}\right)\right\}}^{1/2})$, which follows a $N(0,1)$ asymptotically, can be used for testing as well as for interval estimation of $\mu$. Cohen and Sackrowitz (1984) used Meier’s (1953) approximation along with a suitable critical value in their studentized test statistic.

In the context of the BIBD mentioned earlier, Cohen and Sackrowitz (1989) proposed a test combining individual tests by weighing with respect to their sample variances. This idea was extended by Zhou and Mathew (1993) who proposed two tests and compared their power functions with that of Fisher’s (1932) test. Though no clear-cut winner could be found, one of the tests proposed by Zhou and Mathew (1993) can work well provided one has some prior knowledge about the ratio of two variances.

The test methods of Fisher (1932) and Mathew et al. (1993) are worth mentioning here. Let $t_{0i}^2$ be an observed value of $n_i{({\bar{X}}_i-{\mu }_0)}^2/s_i^2(1\le i\le k)$. Let $P_i=-lnP(t_{n_i-1}^2>t_{0i}^2)$. Fisher’s test is based on the fact that $2{\sum }_{i=1}^k{P}_i\sim {\chi }_{2k}^2$, and it rejects $H_0$ in (1.6) if $P({\chi }_{2k}^2>c)<\alpha$, where $c$ is an observed value of $2{\sum }_{i=1}^k{P}_i$. On the other hand, the approximate test due to Mathew et al. (1993) for $k=2$ suggests using $T^{\ast }=\{n_1{B}_1+n_2{B}_2{\hat{\theta }}^2+2\eta \hat{\theta }\sqrt{n_1{n}_2{B}_1{B}_2}\}/{(\sqrt{n_1}+\hat{\theta }\sqrt{n_2})}^2$, where $B_i=n_i{({\bar{X}}_i-{\mu }_0)}^2/S_i^{\ast },S_i^{\ast }=n_i{({\bar{X}}_i-{\mu }_0)}^2+S_i(i=1,2)$; and $\hat{\theta }=(n_2/n_1){(S_1^{\ast }/S_2^{\ast })}^{1/2}$. For an observed value $t^{\ast }$ of $T^{\ast }$ this test rejects $H_0$ in (1.6) whenever $P(Y>t^{\ast })<\alpha$, where $Y$ is a $\mathrm{Beta}(a,b)$ random variable with specified $a=a(n_1,n_2,\hat{\theta })$ and $b=b(n_1,n_2,\hat{\theta })$ provided by the above authors.

Recently Krishnamoorthy and Lu (2003) proposed a new test based on the generalized $p$-value approach suggested by Tsui and Weerahandi (1989). Though the generalized $p$-value approach doesn’t guarantee that the size would attain the nominal level $\alpha$, the simulation results provided by Krishnamoorthy and Lu (2003) show that the nominal level is attained most of the time.

The generalized $p$-value approach, as initiated by Tsui and Weerahandi (1989), depends heavily on the selection of a suitable pivotal element which is not necessarily unique. The pivot (or, the ‘generalized variable’) used by Krishnamoorthy and Lu (2003) is

$$T_{KL}=\sum _{i=1}^k{W}_i\{{\bar{x}}_i-T_i\sqrt{s_i/\left(n_i(n_i-1)\right)}\}\text{,}$$

with $W_i=(n_i{U}_i/s_i)/({\sum }_{j=1}^k{n}_j{U}_j/s_j)$, where $T_i$’s ($\sim t_{(n_i-1)}$) and $U_i$’s ($\sim {\chi }_{(n_i-1)}^2$) are independent $(1\le i\le k)$, and $({\bar{x}}_i,s_i)$ are the observed values of $({\bar{X}}_i,S_i),(1\le i\le k)$. For the observed data (i.e., $({\bar{x}}_i,s_i)$ of (1.1)), one can generate $T_{KL}$ a large number of times by replicating $T=(T_i,1\le i\le k)$ and $U=(U_i,1\le i\le k)$ from the known distributions. Let $T_{KL,p}$ be the $\left(100p\right)$th percentile point of the above mentioned generated values of $T_{KL}$. Then, an approximate $(1-\alpha )$ level two sided confidence interval for $\mu$ is found as $I_{KL}=(T_{KL,(\alpha /2)},T_{KL,1-(\alpha /2)})$. If ${\mu }_0\in I_{KL}$, then accept $H_0$; and reject $H_0$ otherwise. Alternatively, find $p_1^{\ast }$ ($p_2^{\ast }$) as the proportion of times $T_{KL}$ is less (more) than ${\mu }_0$. Then the $p$-value of the test is: $p$-$\text{value }=2min\{p_1^{\ast },p_2^{\ast }\}$.

Krishnamoorthy and Lu (2003) presented some simulation results comparing five test procedures discussed earlier (Zhou and Mathew’s (1993) two tests, Fisher’s test, Mathew et al.’s (1993) test, and their test based on the pivot $T_{KL}$ given above) for small sample sizes. Again, the results showed no clear-cut winner.

Lin and Lee (2005) used the same generalized $p$-value approach of Tsui and Weerahandi (1989), but with a different pivot given as

$$T_{LL}=\{(\sum _{i=1}^k{n}_i{\bar{x}}_i{U}_i/s_i)-Z\sqrt{\sum _{i=1}^k(n_i{U}_i/s_i)}\}/\sum _{i=1}^k(n_i{U}_i/s_i)$$

where $Z\sim N(0,1)$ and is independent of $U_i$’s $(\sim {\chi }_{(n_i-1)}^2)$, as given after (1.9). For the observed data $({\bar{x}}_i,s_i,1\le i\le k)$, the confidence interval and the subsequent testing for $\mu$ can be carried out by generating a large number of values of $T_{LL}$, similar to what mentioned earlier with $T_{KL}$.

Lin and Lee (2005) compared their generalized $p$-value test with a few others, including the one suggested by Krishnamoorthy and Lu (2003), in terms of power (without mentioning the exact size of their own test). However, it is not clear how closely the size of Lin and Lee’s (2005) test follow the nominal level. We will point this out in our comprehensive simulation study.

In this paper we propose three simple tests which, in our opinion, are natural contenders to those mentioned above. We also include the asymptotic likelihood ratio test (LRT) as a benchmark. Thus, the tests considered here are:

(a) The likelihood ratio test (LRT) which works fairly well for sample sizes greater or equal to 25;

(b) A studentized version of the GDE using the first order unbiased variance estimator proposed by Sinha (1985) (or Meier (1953), since there is hardly much difference between these two first order expressions);

(c) A studentized version of the GDE using the exact unbiased variance estimator proposed by Nikulin and Voinov (1995) (after correcting a mistake in their final expression); and

(d) Using the MLE as a test statistic with the help of a parametric bootstrap method.

The rest of the paper is organized as follows. In Section 2 we discuss the above three test methods in detail–especially the newly developed method using the MLE. In Section 3 we carry out a comprehensive simulation study to compare these four tests with the others discussed earlier. It is seen that the above three tests, the two based on the GDE and the other based on the MLE, exhibit good size and power behavior.