Confidence Intervals for Multinomial Proportions

In a sample of size N, let nj for $j=1,2,\ldots, m$be the observed frequencies in the j-th cell of a m-cell multinomial distribution with corresponding parameter ${\bf p}=(p_1, p_2, \ldots, p_m)'$. The maximum likelihood estimator of ${\bf p}$ is $\hat {\bf p}=({\hat p}_1, {\hat p}_2, \ldots, {\hat p}_m)'$ where $ {\hat p}_j=n_j/N$ for $j=1,2,\ldots, m$. Asymptotically, $\sqrt N (\hat {\bf p} - {\bf p})$has m-variate normal distribution with zero mean vector and variance-covariance matrix $\Sigma$ with elements

\begin{eqnarray*}\sigma_{jj} = & p_j(1-p_j)& (j=1,2,\ldots,m) \\ \sigma_{jk} = &-p_jp_k. & (j\neq k) \end{eqnarray*}


Therefore, the corresponding correlation matrix has singular negative product structure with $\{ \rho^{(m)}_{jk}=-\rho_j\rho_k\}$, where $\rho_j=\sqrt{p_j/(1-p_j)}$, for $j=1,2,\ldots, m$, and rank m-1. Gold (1963) and Quesenberry and Hurst (1964) provided two alternative asymptotic simultaneous confidence intervals for ${\bf p}$ as follows:

\begin{displaymath}p_j \in {\hat p_j} \pm \chi\biggl[\frac{\hat p_j(1-\hat p_j)}{N}\biggr]^{\frac{1}{2}} \end{displaymath} (4)

and

\begin{displaymath}p_j \in \frac {\chi^2+2n_j \pm \biggl\{ \chi^2 \biggl[ \chi^2 + 4 n_j(N-n_j)/N \biggr]\biggr\}^{\frac {1}{2}}}{2(N+\chi^2)} \end{displaymath} (5)

for $j=1,2,\ldots, m$, respectively, while Bailey (1980) suggested the other two confidence intervals for ${\bf p}$ as

\begin{displaymath}p_j \in \biggl\{ \sin \biggl[Y_{1j} \pm \frac{\chi }{\sqrt {4N+2}} \biggr] \biggr\}^2 \end{displaymath} (6)

and

\begin{displaymath}p_j \in \frac{\biggl\{ Y_{2j} \pm \biggl[ C(C+1-Y_{2j}^2 )\biggr]^{\frac{1}{2}} \biggr\}^2 }{(C+1)^2} \end{displaymath} (7)

for $j=1,2,\ldots, m$, where $C=\chi^2/(4N)$, $Y_{1j}=\sin^{-1} \sqrt{(n_j+3/8)/(N+3/4)}$, and $Y_{2j}=\sqrt{(n_j+3/8)/(N+1/8)}$. Based on the Bonferroni inequality, without consideration of the singular correlation structure of $\hat {\bf p}$, Goodman's (1965) criterion (hereafter called GC) is usually used to determine $\chi^2 =\chi_1^2(\alpha /m)$ for the above four techniques in (4)-(7), where $\chi_1^2(\alpha /m)$ is defined as the $100(1-\alpha/m)$percentage point of the chi-square distribution with 1 degree of freedom.

After taking into account the singular correlation structure of ${\bf\hat p}$, Kwong and Iglewicz (1996) proposed another criterion (hereafter called KIC) of setting $\chi^2=z_m^2(\alpha)$ in (4)-(7),where $z_m(\alpha)$ is defined as the two-sided $100(1-\alpha)$ percentage point of the standardized m-variate normal distribution with the singular negative equi-correlated structure. Kwong and Iglewicz (1996) conducted a simulation study to compare the performances of these two criteria. The study concluded that KIC provides less conservative and shorter confidence intervals than GC does when the sample size is sufficiently large.

However, based on the conjecture that the least favorable parameter vector having all the elements equal to one another, KIC uses only part of singular correlation structure. Intuitively, it is possible to improve KIC if the singular negative product correlation structure, instead of singular equi-correlated structure, is incorporated into the evaluation of critical values.

For a given confidence level $1-\alpha$, let

\begin{eqnarray*}H_m[t_{\alpha}({\bf p})] &=& \Pr \biggl[ \bigcap_{j=1}^m \vert ... ...iggr\{ \rho^{(m)}_{jk}=-\rho_i\rho_j \biggr\} \biggr] =1-\alpha. \end{eqnarray*}


By the central limit theorem,

\begin{eqnarray*}\lim _{N \rightarrow \infty} \Pr\biggl[\bigcap _{j=1}^m \frac{\... ..._{jk}=-\rho_i\rho_j \biggr\} \biggr] = H_m[t_{\alpha}({\bf p})]. \end{eqnarray*}


Since the correlation structure of $\{ \rho^{(m)}_{jk} \}$ is dependent on the unknown parameter ${\bf p}$, the sample proportions are used to replace the parameters in order to estimate the correlation structure. Therefore, based on the result in Section 2.2, we numerically solve

\begin{eqnarray*}H_m[t_{\alpha}(\hat {\bf p})] =1-\alpha. \end{eqnarray*}


for $t_{\alpha}(\hat {\bf p})$ and then use it to construct the simultaneous confidence intervals for multinomial proportions when the sample size is relatively large. The numerical solution of this equation requires the use of a nonlinear equation solving method. We considered the use of variations on the secant method for this problem because the computation of the derivative of Hm requires the evaluation of m m-1-dimensional integrals, so the use of more rapidly convergent method, like Newton's method, is infeasible. When $\alpha$ is small, the objective function $h(t) = H_m(t)-1+\alpha$ is very flat near the values of t where $h(t)\approx 0$, and we found that the ordinary secant method frequently diverged unless the starting values for t were very close to the solution. We therefore tried some modified secant methods, which provide at each iteration a bounding interval for the final solution. One of these methods, called the ``Pegasus'' method (see Ralston and Rabinowitz, 1978, p. 341), which we found to be robust, reliable and efficient, was used for the results presented in the remainder of this paper.

The Pegasus method requires an initial bounding interval for the solution point t* ( $t^* = t^*_\alpha(\hat{\bf p})$, with $H(t^*_\alpha(\hat{\bf p}))=1-\alpha$). A crude initial bounding interval can be determined from the simple Bonferroni bound. Using notation similar to that in Section 2.1, we let $A_j(t) = \{X_j: \vert X_j\vert \leq t\}$, and then $P(\cup_{j=1}^m A_j^c(t)) = 1 - H_m(t)$. If we let $S_1(t) = \sum_{j=1}^m P(A_j^c(t))$, then the simple Bonferroni bound is $P(\cup_{j=1}^m A_j^c(t)) \leq S_1(t) $. A simple lower bound is $\min_jP(A_j^c(t)) \leq P(\cup_{j=1}^m A_j^c(t))$. Combining these bounds, and using $P(A_j^c(t)) = 2\Phi(-t)$, a crude bounding interval for t* is given by (tCL, tGC), where $t_{GC} = \Phi^{-1}(1-\alpha/(2m))$ (equivalent to GC previously discussed) and $t_{CL} = \Phi^{-1}(1-\alpha/2)$. When m is large, the repeated evaluation of h(t) requires expensive numerical integrations, so in order to reduce the number of evaluations of h(t), we looked for better bounds that could be used to reduce the size of the initial bounding interval, thereby reducing the number of steps required for the convergence of the Pegasus method. Better bounds can be inexpensively computed from combinations of univariate and bivariate probabilities. The best bounds using these probabilities are the Hunter-Worsley bound for an upper t* limit (see Hsu, 1996, p. 229), and a modified Bonferroni bound for a lower t* limit (see Kwerel, 1975). If we let $S_2(t) = \sum_{j<i} P(A_j^c(t)\cap A_i^c(t))$ , the Hunter-Worsley and modified Bonferroni bounds guarantee that

\begin{displaymath}2(S_1(t)-S_2(t)/k)/(k+1) \leq P(\cup_{j=1}^m A_j^c(t)) \leq S_1(t) - \sum_{(i,j)\in T^*} P(A_j^c(t)\cap A_i^c(t)), \end{displaymath}

where $k = 1 + 2\lfloor S_2/S_1\rfloor$, and T* is a maximal spanning tree for the complete graph of order m with edge weights given by $P(A_j^c(t)\cap A_i^c(t))$. If we use tMB to denote the t value where $2(S_1(t)-S_2(t)/k)/(k+1) = \alpha$ and tHW to denote the t value where $S_1(t) - \sum_{(i,j)\in T^*} P(A_j^c(t)\cap A_i^c(t)) = \alpha$, then $t_{CL} \leq t_{MB} \leq t^* \leq t_{HW} \leq t_{GC}$. A combined procedure for determining t* is to first use the Pegasus method, starting with (tCL, tGC), to determine tHW, then use the Pegasus method, starting with (tCL, tHW), to determine tMB, and finally to use the Pegasus method starting with (tMB, tHW) and combined with numerical integration for computing h(t), to determine t*. For the tests reported in Table 2, we found that the interval (tMB, tHW) is usually very small, particularly for the smaller $\alpha$ values, and in some cases (depending on the accuracy required for t*) it was not necessary to complete the final step in this procedure.

Another problem that arises when an iterative numerical method like the Pegasus method is combined with a numerical integration method, is the problem of balancing the errors from the two methods. We let $\hat{t} = t^*+\epsilon_t$ be an approximation to t*, and $\hat{h(t)}= h(t) + \epsilon_I$ be the numerical integration estimate of h(t). When carrying out the iterations for the Pegasus method, if we attempt to compute h(t) at a point $\hat{t}$ close to t*, we actually compute

\begin{displaymath}\hat{h(\hat{t})}= h(\hat{t}) + \epsilon_I \approx \epsilon_th'(t^*) + \epsilon_I. \end{displaymath}

If we want to be able to stop the hybrid numerical procedure when $\epsilon_t < \delta$ for some small $\delta$, then we need an estimate of h'(t*) so that we can set the absolute error tolerance for our numerical integration method to be less than $\delta/h'(t^*)$. As the Pegasus iterations are carried out we can compute a simple estimate $\hat{h'}$ for h'(t*) with a difference quotient. We found $\hat{h'}$ as small as .01 for some problems for the smallest $\alpha$ values considered, but in these cases we usually have very small starting intervals, and often no numerical integrations are necessary. The numerical integration error tolerance was set at $\delta/(10\hat{h'})$ for the tests that we carried out.

Selected test results are summarized in Table 2. We used the same sets of pj's that were used for the tests that were done to produce Table 1. We also added two other pj sets for m=5. The error tolerance $\delta$ was set at 0.001 for all of the tests. In the cases where $t_{HW}-t_{MB} < 2\delta$, no numerical integration was necessary, and we set t* = (tHW+tMB)/2. For most of these problems, the computation time required to compute the required t* values to within the $\delta = 0.001$ requested accuracy was not significant. In most cases the modified Bonferroni and Hunter-Worseley bonds could be used to quickly provide t* values to within the $\delta = 0.01$, and the KIC value was very accurate. The two additional pj sets were added for m=5 to provide examples where the KIC value is not as accurate. This occurs when there is significant variation in the pj values.


 
Table: Bounds and Estimated Values for $t_\alpha$
pj's $\alpha$ tCL tMB t* tKIC tHW tGC f-Values
(.2, .1, .4, .3) .100 1.645 2.189 2.190 2.193 2.208 2.241 30016
  .050 1.960 2.465 2.466 2.468 2.476 2.498 38864
  .010 2.576 3.011 3.013 3.013 3.014 3.023 23520
  .005 2.807 3.219 3.220 3.221 3.221 3.227 0

(.1, .2, .2, .2, .3)

 .100 1.645 2.288 2.288 2.289 2.308 2.326 44512
  .050 1.960 2.554 2.554 2.555 2.565 2.576 95200
  .010 2.576 3.084 3.084 3.084 3.086 3.090 46768
  .005 2.807 3.286 3.287 3.287 3.288 3.291 0

(.3, .1, .05, .5, .05)

 .100 1.645 2.288 2.288 2.289 2.308 2.326 44512
  .050 1.960 2.554 2.554 2.555 2.565 2.576 95200
  .010 2.576 3.084 3.084 3.084 3.086 3.090 46768
  .005 2.807 3.286 3.287 3.287 3.288 3.291 0
(.1, .1, .2, .2, .2, .2) .100 1.645 2.362 2.363 2.363 2.383 2.394 38256
  .050 1.960 2.621 2.621 2.621 2.632 2.638 50960
  .010 2.576 3.139 3.139 3.139 3.142 3.144 70336
  .005 2.807 3.339 3.340 3.339 3.340 3.341 0

(.10, .30, .05, .50, .05)

 .100 1.645 2.280 2.282 2.289 2.296 2.326 15728
  .050 1.960 2.547 2.549 2.555 2.555 2.576 13712
  .010 2.576 3.079 3.079 3.084 3.080 3.090 0
  .005 2.807 3.282 3.283 3.287 3.283 3.291 0

(.05, .05, .05, .05, .80)

 .100 1.645 2.277 2.280 2.289 2.290 2.326 97088
  .050 1.960 2.546 2.547 2.555 2.551 2.576 137872
  .010 2.576 3.079 3.079 3.084 3.080 3.090 0
  .005 2.807 3.283 3.283 3.287 3.283 3.291 0

(.1, .1, .2, .2, .2, .1,

 .100 1.645 2.421 2.422 2.422 2.441 2.450 72272
.1) .050 1.960 2.675 2.675 2.676 2.685 2.690 112208
  .010 2.576 3.185 3.185 3.185 3.187 3.189 18544
  .005 2.807 3.382 3.383 3.382 3.383 3.384 0

(.1, .1, .1, .1, .15, .05,

 .100 1.645 2.471 2.472 2.473 2.490 2.498 160144
.2, .2) .050 1.960 2.721 2.722 2.721 2.730 2.734 154880
  .010 2.576 3.224 3.225 3.224 3.226 3.227 20048
  .005 2.807 3.419 3.419 3.419 3.420 3.421 0

(.01, .02, .07, .1, .15, .05,

 .100 1.645 2.514 2.515 2.516 2.531 2.539 33680
.3, .2, .1) .050 1.960 2.760 2.761 2.761 2.768 2.773 50000
  .010 2.576 3.258 3.259 3.258 3.259 3.261 0
  .005 2.807 3.451 3.451 3.451 3.452 3.452 0

(.1, .05, .05, .04, .06, .1,

 .100 1.645 2.553 2.554 2.554 2.570 2.576 55536
.15, .15, .1, .2) .050 1.960 2.796 2.796 2.796 2.804 2.807 484416
  .010 2.576 3.288 3.289 3.288 3.290 3.291 0
  .005 2.807 3.479 3.480 3.479 3.480 3.481 0

(.02, .08, .04, .06, .1, .1,

 .100 1.645 2.586 2.588 2.587 2.604 2.609 691008
.16, .14, .15, .1, .05) .050 1.960 2.827 2.827 2.827 2.835 2.838 489952
  .010 2.576 3.315 3.316 3.315 3.317 3.317 0
  .005 2.807 3.505 3.505 3.505 3.506 3.506 0

(.01, .03, .06, .05, .05, .1,

 .100 1.645 2.617 2.618 2.618 2.634 2.638 501248
.15, .05, .1, .14, .16, .1) .050 1.960 2.855 2.857 2.855 2.863 2.865 74544
  .010 2.576 3.339 3.340 3.339 3.341 3.341 0
  .005 2.807 3.528 3.529 3.528 3.529 3.529 0


Alan C Genz
2001-02-09