Confidence interval estimating procedures for standardized incidence rates

Computational Statistics and Data Analysis 52 (2008) 3501–3516 www.elsevier.com/locate/csda

Confidence interval estimating procedures for standardized incidence rates Hon Keung Tony Ng a,∗ , Giovanni Filardo a,b , Gang Zheng c a Department of Statistical Science, Southern Methodist University, Dallas, TX 75275, USA b Institute for Health Care Research and Improvement, Baylor Research Institute, Dallas, TX 75206, USA c Office of Biostatistics Research, National Heart, Lung and Blood Institute, Bethesda, MD 20892, USA

Received 3 May 2007; received in revised form 27 August 2007; accepted 13 November 2007 Available online 19 November 2007

Abstract Performances of some common confidence interval construction procedures for direct standardized incidence or mortality rate were compared. The asymptotic method, the method based on matching moments, the method based on gamma distribution and the method based on beta distribution were investigated through Monte Carlo simulations based on data from an epidemiologic study of cardiovascular disease and randomly generated scenarios. These interval estimation procedures, as well as some mixed procedures developed for the present study, were evaluated in terms of their coverage probability and expected length. c 2007 Elsevier B.V. All rights reserved.

1. Introduction In epidemiology, age-standardized incidence rates of occurrence of disease or mortality are measures of primary interest. Standardization is an elemental tool used to minimize distortion during the comparison of rates. Often standardization is achieved with age adjustment of the disease incidence (or mortality whichever is the measure of interest) rate through a weighted sum of age-specific rates. Since the occurrences of specific diseases are often modeled using a Poisson distribution (Selvin, 1996; Rothman and Greenland, 1998; Gail and Benichou, 2000), direct standardizations of disease incidence rates are examples of weighted sums of Poisson rate parameters. Several methods can be used to estimate confidence intervals for standardized rates. The most commons are the asymptotic method, the moment matching method based on a single Poisson parameter, and the methods based on gamma and beta distributions. The asymptotic method is an approach frequently used to construct confidence intervals for rates. This method is based on the normal approximation of the Poisson distribution. Although the method is relatively straightforward, the performance might be inconsistent. The moment matching method was proposed by Dobson et al. (1991) for obtaining approximate confidence limits for the weighted sums of Poisson parameters as linear functions of confidence limits for a single Poisson parameter. ∗ Corresponding author. Tel.: +1 214 7682465; fax: +1 214 7684035.

E-mail address: [email protected] (H.K.T. Ng). c 2007 Elsevier B.V. All rights reserved. 0167-9473/$ - see front matter doi:10.1016/j.csda.2007.11.004

3502

H.K.T. Ng et al. / Computational Statistics and Data Analysis 52 (2008) 3501–3516

Constructions of interval estimation of a single Poisson rate have long been used (Garwood, 1936; Stevens, 1957; Crow and Gardner, 1959). In recent years, the development of many novel interval methods of estimation of a single Poisson parameter (Cohen and Yang, 1994; Kabaila and Byrne, 2001; Johnson et al., 2005) permitted the construction of new and possibly more efficient interval estimation methods for weighted sums of Poisson parameter based on the Dobson et al. (1991) approach. An approach alternative to the one proposed by Dobson et al. (1991) for obtaining confidence intervals for weighted sums of Poisson parameters was proposed by Fay and Feuer (1997) and Tiwari et al. (2006). Both methods designed to obtain confidence intervals for direct standardized rates were based on gamma distribution. Tiwari et al. (2006) also developed an alternative approach based on beta distribution. Although the aforementioned approaches are commonly used by researchers in comparing standardized rates and ultimately drawing conclusions, a systematic comparison of confidence interval construction procedures based on the asymptotic method, the moment matching method of Dobson et al. (1991), the method based on gamma distribution of Fay and Feuer (1997) and Tiwari et al. (2006), and the method based on beta distribution of Tiwari et al. (2006) has never been conducted. In this article we are presenting a number of confidence intervals procedures based on the combinations of the four aforementioned methods. Performances of these procedures were compared by means of the Monte Carlo method. Notation as regards the interval estimation of direct standardized incidence rates is defined in Section 2. In this section there are also reported details on the four different confidence interval construction methods considered: the asymptotic method, the method based on matching moments (Dobson et al., 1991), the method based on gamma distribution (Fay and Feuer, 1997), and the method based on beta distribution (Tiwari et al., 2006). The various confidence interval construction procedures based on the aforesaid methods and their computational formulae are also presented in the same section. In Section 3, a large scale Monte Carlo simulation study carried out to investigate the performance of the confidence interval construction procedures in terms of their coverage probabilities and the expected lengths of confidence intervals is presented. Findings and discussion are also reported in Section 3. Data from a cardiovascular disease study are provided in Section 4 to show procedure performance. Lastly, conclusions are reported in Section 5. 2. Confidence intervals for standardized incidence rates Incidence (or mortality) rates are commonly studied using direct standardized rates. The numbers of events observed in each age group are independent and are assumed to be Poisson distributed. Use of standardized rates raises the question of how best to construct their confidence intervals and thus a weighted sum of Poisson parameters. Let the random variables X i , i = 1, 2, . . . , k, be the number of occurrences of an event (disease or death) in the i-th subgroups (e.g. age groups). Suppose each X i is independent and has a Poisson distribution with mean n i λi , where λi is the incidence rate and n i is the sampling frame (number of person-years) for subgroup i. We further define θi = n i λi , that is, X i ∼ Poisson(θi ) for i = 1, 2, . . . , k. A standardized incidence (or mortality) rate is defined as θ=

k X

ci λi =

i=1

k X

wi θi

i=1

where wi = ci /n i and the ci ’s are the corresponding numbers of person-years of the standard population (for example, this can be obtained from the Segi World Standard Population) which are known constants for i = 1, 2, . . . , k. The maximum likelihood estimator (MLE) of θ is Y =

k X

wi X i ,

i=1

and we have the expected value E(Y ) = Var(Y ) is d )= V = Var(Y

k X i=1

wi2 X i .

Pk

i=1 wi θi

and variance Var(Y ) =

Pk

2 i=1 wi θi .

An estimate of the variance

H.K.T. Ng et al. / Computational Statistics and Data Analysis 52 (2008) 3501–3516

3503

We are interested in the interval estimation of the unknown parameter θ. 2.1. Asymptotic method A common approach for constructing confidence intervals for the parameter θ is using the normal approximation for the Poisson distribution for each X i while the linear combination of the independent normal random variables is considered normal. Different normal approximations of the MLE or transformed MLE distribution can also be used as these are good standards for comparison. (A1) Normal approximation of the MLE Y −E(Y ) On the basis of the normal approximation of the MLE, we assumed that √ is asymptotically normally Var(Y ) distributed with zero mean and unit variance. Therefore, a 100(1 − α)% confidence interval for θ is v ! ! u k k u X X √ t Y ±z V = w X ±z w2 X , (1) 1−α/2

i

1−α/2

i

i=1

i

i

i=1

where z q is the q-th upper percentile of a standard normal distribution. (A2) Normal approximation of the log-transformed MLE The problem with applying A1 is that when θi ’s are small, the normal approximation may be poor. However a different transformation of the MLE can be used to correct the inadequate performance of the normal approximation. Based on the normal approximation of the log-transformed MLE, a 100(1 − α)% confidence interval for θ is h i p exp ln(Y ) ± z 1−α/2 Var(ln Y ) (2) where Var(ln Y ) can be approximated by the delta method as k P

d Y ) = V = i=1 Var(ln  k Y2 P

wi2 X i 2 .

wi X i

i=1

(A3) Normal approximation of the cubic-root-transformed MLE Based on the normal approximation of the cubic-root-transformed MLE, a 100(1 − α)% confidence interval for θ is i3 h p (3) Y 1/3 ± z 1−α/2 Var(Y 1/3 ) where Var(Y 1/3 ) can be approximated by the delta method as k P

d 1/3 ) = V = Var(Y 9Y 4/3

i=1

 9

k P

wi2 X i

wi X i

4/3 .

i=1

(A4) Edgeworth correction for skewness Y −E(Y ) If we assume that √ is asymptotically normally distributed with zero mean and unit variance, the Var(Y ) Y −E(Y ) distribution function of √ can be written in a series form as Var(Y )   Y − E(Y ) Pr √ ≤ t ≈ Φ(t) + φ(t)g(t) Var(Y ) where √ √     β2 − 3 15β1 β2 − 3 5β1 3 β1 5 β1 β1 2 + + t− t − + t − t , g(t) = 6 8 72 6 24 36 72

(4)

3504

H.K.T. Ng et al. / Computational Statistics and Data Analysis 52 (2008) 3501–3516

√ β1 and β2 are the coefficients of skewness and kurtosis of Y , respectively, and Φ(·) is the c.d.f. of the standard normal distribution with corresponding p.d.f. φ(·). Formula (4) is known as the Edgeworth expansion (Hall, 1992; Johnson et al., 1994). The corrected confidence interval can be constructed as h  √  √ i Y − z 1−α/2 − g(z ˆ 1−α/2 ) V , Y + z 1−α/2 − g(z ˆ α/2 ) V √ where g(t) ˆ is identical to g(t) except that β1 and β2 are replaced by their estimates based on the sample moments k P

q

βˆ1 = 

i=1 k P

i=1

k P

wi3 X i

wi2 X i

3/2

and

βˆ2 =

i=1

wi4 X i (1 + 3X i )



k P

i=1

wi2 X i

2

,

(5)

respectively. 2.2. Moment matching method Dobson et al. (1991) proposed a method for obtaining approximate confidence limits for the weighted sum of Poisson Pkparameters as linear functions of confidence limits for a single Poisson parameter (the unweighted sum X = i=1 X i ) by matching the first and second moments. Suppose the confidence interval of the unweighted sum Pk of Poisson parameters i=1 θi is (X L , X U ); then an approximate confidence interval of the weighted sum of Poisson parameters θ is (TL , TU ) with  1/2 V TL = Y + (X L − X ) (6) X  1/2 V TU = Y + (X U − X ). (7) X Pk Since the construction of confidence intervals for i=1 θi is the problem inherent to constructing Poisson confidence intervals, a multitude of methods have been proposed to address this issue. The expectation of a Poisson random variable for each confidence interval method can have the corresponding approximate confidence interval for θ based on the moment matching method. Here, weP considered the following methods for constructing confidence intervals for the unweighted sum of Poisson k parameters i=1 θi : (M1) Exact confidence interval (Johnson et al., 2005, Ch. 4) On the basis of the relationship between the Poisson and chi-squared distributions, we can construct the exact confidence interval for the unweighted sum of Poisson parameters as 1 2 1 2 X L = χα/2 X U = χ1−α/2 (2X ) , (2(X + 1)) , 2 2 where χq2 (d) is the q-th quantile of the chi-squared distribution with d degrees of freedom. (M2) Boice–Monson method (Rothman and Greenland, 1998)     1 1 X L = exp ln X − z 1−α/2 √ , X U = exp ln X + z 1−α/2 √ . X X (M3) Normal approximation (Johnson et al., 2005, Ch. 4) Solving the inequalities k P

−z 1−α/2

k P Xi − θi i=1 i=1 ≤ < z 1−α/2 s k P θi i=1

H.K.T. Ng et al. / Computational Statistics and Data Analysis 52 (2008) 3501–3516

yields 100(1 − α)% confidence limits for

Pk

i=1 θi .

2 z 1−α/2 1 2 X L = X + z 1−α/2 − z 1−α/2 X + 2 4

XU

3505

2 z 1−α/2 1 2 + z 1−α/2 X + = X + z 1−α/2 2 4

We solve the quadratic equations and obtain !1/2 , !1/2 .

(M4) More complicated normal approximation (Rothman and Boice, 1979) Based on the Wilson–Hilferty approximation to the chi-squared distribution (Johnson et al., 1994, Ch. 18), a 100(1 − α)% confidence limits for θ is  !1/2  2 2 z 1−α/2 +2 1 (2z 1−α/2 + 1) 1 , −  + z 1−α/2 X + − XL = X + 6 2 18 2  !1/2  2 2 (2z 1−α/2 + 1) z 1−α/2 +2 1 1 . +  + z 1−α/2 X + + XU = X + 6 2 18 2 (M5) Byar’s method (Rothman and Boice, 1979)     z 1−α/2 3 z 1−α/2 3 1 1 XL = X 1 − − √ , X U = (X + 1) 1 − + √ . 9X 9(X + 1) 3 X + 1 3 X (M6) Exact short confidence interval of Crow and Gardner (1959) The exact short confidence interval of Crow and Gardner (1959) has provided a method for computing exact confidence intervals. They also tabulated the two-sided confidence intervals for the Poisson mean for different confidence level and observed values from 0 to 300. (M7) Exact mid- p confidence interval (Cohen and Yang, 1994) The mid- p 100(1 − α)% confidence interval is defined as (X L , X U ) where X L is the smallest non-negative µ such that 1 α Pr(X = x|µ) ≥ 2 2 and X U is the largest µ such that Pr(X > x|µ) +

1 α Pr(X = x|µ) ≥ , 2 2 where X is Poisson distributed with mean µ and x is the observed value of random variable X . (M8) Approximation to mid- p confidence interval  3   z 1−α/2 1 1   , 1− q − q XL = X + 2 1 1 9 X+ 3 X+ Pr(X < x|µ) +

2

 XU =





1  1 X+ 1− q 2 9 X+

1 2

2

3 z 1−α/2  . + q 1 3 X+2

(M9) Simple approximation to mid- p confidence interval (Vandenbroucke, 1982) !2 !2 r r 1 z 1−α/2 1 z 1−α/2 , XU = . XL = X+ − X+ + 2 2 2 2 (M10) Exact short confidence intervals of Kabaila and Byrne (2001)

3506

H.K.T. Ng et al. / Computational Statistics and Data Analysis 52 (2008) 3501–3516

Kabaila and Byrne (2001) proposed a confidence interval construction method for the expectation of a Poisson random variable. Suppose s is the smallest positive integer such that maxµ Pr(x −s ≤ X ≤ x −1|µ) > 1 − α and q is the smallest positive integer such that maxµ Pr(x + 1 ≤ X ≤ x + q|µ) > 1 − α.  0, for x = 0, XL = max{µ : Pr(x − s ≤ X ≤ x − 1|µ) = 1 − α}, otherwise, and X U = max{µ : Pr(x ≤ X ≤ x + q − 1|µ) = 1 − α}. 2.3. Methods based on gamma distribution Fay and Feuer (1997) proposed an approximation of confidence intervals for standardized mortality rates based on gamma distribution similar to the method proposed by Dobson et al. (1991). The 100(1 − α)% confidence interval (L , U ) is  2 2Y V 2 L= χα/2 , (8) 2Y V   2(Y + w ∗ )2 V + w ∗2 2 χ , (9) U= 2(Y + w ∗ ) 1−α/2 V + w ∗2 where w ∗ = w(k) = max(w1 , w2 , . . . , wk ) (say, method G1). Since Fay and Feuer (1997) showed that these confidence intervals are conservative, we also considered w ∗ = w(1) = min(w1 , w2 , . . . , wk ) (say, method G2) w +w and w ∗ = (1) 2 (k) (say, method G3) to account for the conservativeness. Tiwari et al. (2006) proposed a modification on the Fay and Feuer (1997) method based on continuity correction which assumed distributing a count of 1 uniformly to all subgroups. The proposed 100(1 − α)% confidence interval is (L , U ∗ ) where L is given in (8) and  ∗2  2Y V∗ 2 U∗ = χ , 1−α/2 ∗ 2Y V∗ with Y∗ = Y +

k 1X wi k i=1

and

V∗ = V +

k 1X w2 . k i=1 i

We named this method as method G4. Pk It is worth mentioning that Tiwari et al. (2006) have also considered the upper confidence limit U with w ∗ = k1 i=1 wi and found that these two intervals performed very similarly. Note that although there are an infinite number of choices for w ∗ which will result in different upper confidence limits, method G4 proposed by Tiwari et al. (2006) is proposed based on the idea of continuity correction. 2.4. Methods based on beta distribution Tiwari et al. (2006) proposed a method of constructing confidence intervals for age-adjusted incidence rates which can be adopted here to construct confidence intervals for the overall standardized incidence rates. An approximate 100(1 − α)% confidence interval (L b , Ub ) can be obtained by solving the following beta integrals: Z Lb α h(x|a, b)dx = , 2 0 Z Ub α h(x|a, b)dx = 1 − , 2 0 where h(·|a, b) is the probability density function of beta distribution with parameters a and b. Following the method proposed by Tiwari et al. (2006), the values of a and b are given by     Y (1 − Y ) Y (1 − Y ) −1 , b = (1 − Y ) −1 , a=Y V V

H.K.T. Ng et al. / Computational Statistics and Data Analysis 52 (2008) 3501–3516

3507

and we named this method as method B1. Tiwari et al. (2006) considered the case when the beta parameters a and b are defined as  ∗   ∗  Y (1 − Y ∗ ) Y (1 − Y ∗ ) ∗ , b = (1 − Y ) , a = Y∗ − 1 − 1 V∗ V∗ namely method B2. Like for the methods based on gamma distribution, there are an infinite number of choices for Y ∗ and V ∗ which will result different confidence intervals; however, the method B2 is well motivated by continuity correction. 3. Monte Carlo simulation studies 3.1. Specific settings In this section, Monte Carlo simulation was employed to investigate the performance of the interval estimation procedures. Criteria appropriate to the evaluation of the various methods under scrutiny include: closeness of the coverage probability to its nominal value; and expected interval width. The International Mathematical and Statistical Libraries (IMSL) subroutine RNPOI was employed for generating the desired Poisson samples. Simulation study parameters were based on Dobson et al. (1991) corresponding to the Segi World Standard Population for six age groups (35-39, 40-44, 45-49, 50-54, 55-59, 60-64). Note that Dobson et al. (1991) compared the performance of procedures M1 and M6 only. Here, we considered three cases involving small values of θi which may apply for death rates per 10,000 people. Suppose there are six groups each of n i = 10, 000 and the values of ci are {6/31, 6/31, 6/31, 5/31, 4/31, 4/31}; we considered the following settings: P6 P6 (1) The event rates λi are {1,1,2,4,5,7} per 10,000, i.e. θ = i=1 wi θi = 2.968 and i=1 θi = 20. P6 P6 (2) The event rates λi are {0.5,0.5,1,2,2.5,3.5} per 10,000, i.e. θ = i=1 wi θi = 1.484 and i=1 θ = 10. P6 P6 i (3) The event rates λi are {2.5,2.5,5,10,12.5,17.5} per 10,000, i.e. θ = i=1 wi θi = 7.419 and i=1 θi = 50. We also considered a special case: when all the wi ’s are equal, i.e. ci are {1/6, 1/6, 1/6, 1/6, 1/6, 1/6} and P6 P6 (4) The event rates λi are {1,1,2,4,5,7} per 10,000, i.e. θ = i=1 wi θi = 3.333 and i=1 θi = 20. P6 P6 (5) The event rates λi are {0.5,0.5,1,2,2.5,3.5} per 10,000, i.e. θ = i=1 wi θi = 1.667 and i=1 θ = 10. P6 P6 i (6) The event rates λi are {2.5,2.5,5,10,12.5,17.5} per 10,000, i.e. θ = i=1 wi θi = 8.333 and i=1 θi = 50. When all weights are equal or when some the weights are equal and the rest are zero, methods M1, G1, G2, G3 and G4 are equivalent, as are methods A2 and M2. In this situation, the standardized incidence rate can be written as a scaled Poisson variate and we can obtain an exact confidence interval (Hirji, 2006) by methods M1 (which is equivalent to G1–G4), M6 and M10. Note that method M1 gives an exact central confidence interval while methods M6 and M10 produce exact non-central confidence intervals. These exact confidence intervals are conservative, and maintain the coverage probability at or above the desired confidence level (see, for example, Casella and Robert (1989); Fleiss et al. (2003) and Hirji (2006)). For the purpose of comparison, we considered the following settings with a wider range of wi values and all the θi ’s are equal: P6 P6 (7) wi are {0.75,0.05,0.05,0.05,0.05,0.05}, θi = 2, i = 1, 2, . . . , 6, i.e. θ = i=1 wi θi = 2 and i=1 θi = 10. P6 P6 (8) wi are {0.75,0.05,0.05,0.05,0.05,0.05}, θi = 4, i = 1, 2, . . . , 6, i.e. θ = i=1 wi θi = 4 and i=1 θ = 20. P6 P6 i (9) wi are {0.75,0.05,0.05,0.05,0.05,0.05}, θi = 10, i = 1, 2, . . . , 6, i.e. θ = i=1 wi θi = 10 and i=1 θi = 50. For each simulated sample under a particular setting, we computed confidence intervals and checked whether the true value lay within the interval and recorded the length of the confidence interval. This procedure was repeated 10,000 times. The estimated probability coverage was computed as the number of confidence intervals that covered the true values divided by 10,000 while the estimated conditional expected width of the confidence interval was computed as the sum of the lengths for all intervals that covered the true values divided by the number of confidence intervals that covered the true values. The simulation study considered 90%, 95% and 99% confidence intervals and the corresponding simulation standard errors (SEs) were 0.0030, 0.0022 and 0.0010, respectively. The simulation results are summarized and presented in Tables 1–5 for settings (1) to (9) with 90%, 95% and 99% confidence intervals.

3508

H.K.T. Ng et al. / Computational Statistics and Data Analysis 52 (2008) 3501–3516

Table 1 Empirical coverage probabilities for weighted sums of Poisson parameters, based on 10,000 simulations for settings (1) and (2) Setting (1) Method

90% C.I. (SE = 0.0030) Left Right Coverage (in %) (in %) (in %)

Exp. width

95% C.I. (SE = 0.0022) Left Right Coverage (in %) (in %) (in %)

Exp. width

99% C.I. (SE = 0.0010) Left Right Coverage (in %) (in %) (in %)

Exp. width

A1 A2 A3 A4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 G1 G2 G3 G4 B1 B2

2.81 5.42 4.40 2.11 3.81 5.27 5.45 3.81 3.81 3.89 4.86 4.86 5.08 5.88 3.81 3.81 3.81 3.81 3.81 6.55

2.229 2.251 2.221 2.321 2.370 2.251 2.238 2.370 2.370 2.257 2.235 2.217 2.230 2.363 2.429 2.343 2.385 2.361 2.225 2.275

0.97 2.84 2.20 1.29 1.82 2.80 2.87 1.82 1.82 1.82 2.37 2.37 2.37 2.88 1.82 1.82 1.82 1.82 1.82 3.39

2.654 2.709 2.650 2.789 2.800 2.708 2.685 2.800 2.801 2.679 2.665 2.646 2.661 2.761 2.862 2.769 2.814 2.784 2.625 2.732

0.05 0.72 0.40 1.47 0.36 0.70 0.75 0.36 0.36 0.36 0.43 0.45 0.40 0.52 0.36 0.36 0.36 0.36 0.37 0.75

3.471 3.649 3.496 3.734 3.643 3.643 3.588 3.642 3.648 3.516 3.514 3.495 3.496 3.631 3.717 3.609 3.660 3.613 3.441 3.660

Setting (2)

90% C.I. (SE = 0.0030)

A1 A2 A3 A4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 G1 G2 G3 G4 B1 B2

2.08 5.41 3.99 1.50 3.33 5.31 5.41 3.33 3.33 3.12 4.32 4.36 4.80 4.62 3.33 3.33 3.33 3.33 3.33 7.42

7.85 4.03 5.24 7.51 3.87 4.26 4.51 3.87 3.87 5.41 5.23 4.26 4.52 3.35 3.28 4.39 3.82 3.87 6.68 3.87

9.62 3.40 5.86 8.12 3.28 3.40 4.35 3.28 3.40 5.32 4.72 3.40 4.64 3.38 2.37 4.35 3.40 4.34 7.81 4.34

89.34 90.55 90.36 90.38 92.32 90.47 90.04 92.32 92.32 90.70 89.91 90.88 90.40 90.77 92.91 91.80 92.37 92.32 89.51 89.58

88.30 91.19 90.15 90.38 93.39 91.29 90.24 93.39 93.27 91.56 90.96 92.24 90.56 92.00 94.30 92.32 93.27 92.33 88.86 88.24

5.24 1.80 2.86 3.35 2.04 1.90 2.01 2.04 2.04 2.58 2.47 1.90 2.53 1.65 1.40 2.05 1.76 2.05 3.35 2.05

93.79 95.36 94.94 95.36 96.14 95.30 95.12 96.14 96.14 95.60 95.16 95.73 95.10 95.47 96.78 96.13 96.42 96.13 94.83 94.56

95% C.I. (SE = 0.0022) 1.585 1.616 1.575 1.671 1.720 1.614 1.598 1.720 1.720 1.597 1.592 1.565 1.589 1.668 1.785 1.693 1.737 1.704 1.575 1.670

0.73 2.86 2.05 0.89 1.68 2.85 2.88 1.68 1.66 1.72 2.12 2.27 2.24 2.70 1.68 1.68 1.68 1.68 1.68 3.66

6.56 1.48 3.28 4.42 1.56 1.48 1.56 1.56 1.56 2.80 2.35 1.48 2.35 1.56 1.12 2.03 1.56 2.05 4.35 2.05

92.71 95.66 94.67 94.69 96.76 95.67 95.56 96.76 96.78 95.48 95.53 96.25 95.41 95.74 97.20 96.29 96.76 96.27 93.97 94.29

2.05 0.13 0.47 0.31 0.35 0.13 0.22 0.31 0.31 0.47 0.43 0.13 0.57 0.31 0.24 0.41 0.28 0.41 0.88 0.41

97.90 99.15 99.13 98.22 99.29 99.17 99.03 99.33 99.33 99.17 99.14 99.42 99.03 99.17 99.40 99.23 99.36 99.23 98.75 98.84

99% C.I. (SE = 0.0010) 1.882 1.966 1.881 2.017 2.025 1.962 1.927 2.024 2.026 1.906 1.900 1.872 1.893 2.005 2.099 1.993 2.043 1.996 1.856 1.994

0.03 0.83 0.44 0.85 0.35 0.82 0.89 0.37 0.30 0.46 0.46 0.49 0.42 0.54 0.37 0.37 0.37 0.37 0.37 0.89

3.28 0.00 0.61 1.12 0.41 0.00 0.12 0.41 0.41 0.60 0.60 0.00 0.61 0.41 0.12 0.46 0.19 0.61 1.12 0.61

96.69 99.17 98.95 98.03 99.24 99.18 98.99 99.22 99.29 98.94 98.94 99.51 98.97 99.05 99.51 99.17 99.44 99.02 98.51 98.50

2.455 2.715 2.488 2.720 2.634 2.707 2.616 2.632 2.642 2.517 2.514 2.485 2.487 2.637 2.723 2.594 2.653 2.580 2.409 2.733

We defined a confidence interval construction procedure as being conservative when the value of (1 − simulated coverage probability) was 20% lower than the nominal level (i.e., > 99.2% for 99% confidence intervals, > 96% for 95% confidence intervals and > 92% for 90% confidence intervals), liberal when the value of (1−simulated significance level) was 20% higher than the nominal level (i.e., < 98.8% for 99% confidence intervals, < 94% for 95% confidence intervals and < 88% for 90% confidence intervals), and robust otherwise. A liberal procedure should be used with caution as it provides coverage probability below the pre-chosen nominal level. Conversely, conservative procedures are of less concern. Confidence interval construction procedures within each group, asymptotic method (say, A-group), moment matching method (say, M-group), based on gamma (say, G-group) and based on beta distribution (say, B-group) were compared.

3509

H.K.T. Ng et al. / Computational Statistics and Data Analysis 52 (2008) 3501–3516 Table 2 Empirical coverage probabilities for weighted sums of Poisson parameters, based on 10,000 simulations for settings (3) and (4) Setting (3) Method

90% C.I. (SE = 0.0030) Left Right Coverage (in %) (in %) (in %)

Exp. width

95% C.I. (SE = 0.0022) Left Right Coverage (in %) (in %) (in %)

Exp. width

99% C.I. (SE = 0.0010) Left Right Coverage (in %) (in %) (in %)

Exp. width

A1 A2 A3 A4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 G1 G2 G3 G4 B1 B2

3.33 5.05 4.42 2.66 4.03 5.04 5.05 4.03 4.03 3.50 4.55 4.55 4.73 5.32 4.03 4.03 4.03 4.03 4.03 5.58

3.511 3.525 3.507 3.604 3.656 3.525 3.517 3.656 3.656 3.549 3.515 3.504 3.513 3.614 3.707 3.630 3.668 3.656 3.508 3.536

1.45 2.60 2.17 1.90 1.97 2.55 2.60 1.97 1.97 1.95 2.26 2.27 2.27 2.76 1.97 1.97 1.97 1.97 1.97 2.84

4.179 4.216 4.179 4.317 4.329 4.214 4.201 4.328 4.329 4.611 4.189 4.176 4.186 4.309 4.385 4.302 4.342 4.326 4.163 4.229

0.11 0.58 0.39 2.34 0.34 0.58 0.59 0.34 0.34 0.33 0.39 0.40 0.36 0.49 0.34 0.34 0.34 0.34 0.34 0.55

5.482 5.595 5.501 5.746 5.653 5.592 5.559 5.649 5.653 5.530 5.513 5.501 5.501 5.637 5.714 5.623 5.668 5.642 5.463 5.599

Setting (4) Method

90% C.I. (SE = 0.0030) Left Right Coverage (in %) (in %) (in %)

Exp. width

95% C.I. (SE = 0.0022) Left Right Coverage (in %) (in %) (in %)

Exp. width

99% C.I. (SE = 0.0010) Left Right Coverage (in %) (in %) (in %)

Exp. width

A1 A2, M2 A3 A4 M1, G1–G4 M3 M4 M5 M6 M7 M8 M9 M10 B1 B2

3.35 5.07 5.07 2.02 3.35 5.07 3.35 3.35 5.07 5.07 5.07 5.07 7.65 3.35 7.65

2.457 2.492 2.461 2.559 2.627 2.478 2.627 2.627 2.487 2.468 2.452 2.467 2.600 2.457 2.565

1.22 3.35 2.02 1.22 2.02 3.35 2.02 2.02 2.02 2.02 2.02 2.02 3.35 2.02 3.35

2.927 2.997 2.930 3.085 3.100 2.971 3.099 3.101 2.961 2.951 2.934 2.946 3.049 2.952 2.997

0.08 0.73 0.44 1.25 0.44 0.73 0.44 0.44 0.44 0.44 0.44 0.44 0.44 0.44 0.73

3.843 4.034 3.868 4.124 4.034 3.973 4.030 4.037 3.891 3.891 3.870 3.869 4.023 3.811 4.044

6.72 4.56 5.41 6.64 4.35 4.59 4.65 4.35 4.35 5.70 5.12 4.59 4.95 3.62 3.74 4.57 4.20 4.31 5.94 4.31

6.40 3.97 6.40 6.40 3.97 3.97 3.97 3.97 3.97 3.97 3.97 3.97 2.17 6.40 3.97

89.95 90.39 90.17 90.70 91.62 90.37 90.30 91.62 91.62 90.80 90.33 90.86 90.32 91.06 92.23 91.40 91.77 91.66 90.03 90.11

90.25 90.96 88.53 91.58 92.68 90.96 92.68 92.68 90.96 90.96 90.96 90.96 90.18 90.25 88.38

3.94 2.03 2.68 2.57 2.10 2.05 2.10 2.10 2.10 2.81 2.55 2.05 2.55 1.83 1.86 2.23 2.02 2.10 3.05 2.10

3.97 2.17 2.17 3.97 2.17 2.17 2.17 2.17 2.17 2.17 2.17 2.17 1.00 3.97 2.17

94.61 95.37 95.15 95.53 95.93 95.40 95.30 95.93 95.93 95.24 95.19 95.68 95.18 95.41 96.17 95.80 96.01 95.93 94.98 95.06

94.81 94.48 95.81 94.81 95.81 94.48 95.81 95.81 95.81 95.81 95.81 95.81 95.65 94.01 94.48

1.16 0.31 0.55 0.10 0.42 0.33 0.33 0.42 0.42 0.54 0.53 0.33 0.55 0.34 0.33 0.43 0.40 0.42 0.59 0.42

2.17 0.19 0.43 0.31 0.43 0.19 0.19 0.19 0.43 0.43 0.19 0.43 0.43 1.00 0.43

98.73 99.11 99.06 97.56 99.24 99.09 99.08 99.24 99.24 99.13 99.08 99.27 99.09 99.17 99.33 99.23 99.26 99.24 99.07 99.03

97.75 99.08 99.13 98.44 99.13 99.08 99.37 99.37 99.13 99.13 99.37 99.13 99.13 98.56 98.84

A-group: those procedures based on log-transformed MLE (A2) and cubic-root-transformed MLE (A3) were found to be robust and/or conservative while A1 and A4 were found to be liberal in the majority of the cases. Therefore, both A2 and A3 are recommended in the A-group but if one also considers the expected widths of the confidence intervals, A3 produced shorter confidence intervals on average compared to A2. M-group: the procedures based on the mid- p confidence interval (M7, M8 and M9) had similar performances and provided better coverage probabilities than the other methods in the M-group. Moreover, we observed that M8 was conservative most of the time and provided a smaller expected width. Furthermore, M8 is easier to compute than M7 since it uses an approximation to the mid- p confidence interval. Therefore, M8 might be considered the best procedure in the M-group in terms of the coverage probabilities and expected width. When the symmetric of the confidence interval (i.e. left-side error rate equal to right-tail error rate) is taken into account, M7 is a better procedure in the M-group. G-group: we observed that the choice of the value w ∗ affected the conservativeness of the procedures: the larger the value of w ∗ the more conservative the procedure. The expected width of the confidence interval was found to be

3510

H.K.T. Ng et al. / Computational Statistics and Data Analysis 52 (2008) 3501–3516

Table 3 Empirical coverage probabilities for weighted sums of Poisson parameters, based on 10,000 simulations for settings (5) and (6) Setting (5) Method

90% C.I. (SE = 0.0030) Left Right Coverage (in %) (in %) (in %)

Exp. width

95% C.I. (SE = 0.0022) Left Right Coverage (in %) (in %) (in %)

Exp. width

99% C.I. (SE = 0.0010) Left Right Coverage (in %) (in %) (in %)

Exp. width

A1 A2, M2 A3 A4 M1, G1–G4 M3 M4 M5 M6 M7 M8 M9 M10 B1 B2

2.47 4.38 4.38 1.32 4.38 4.38 4.38 4.38 2.47 4.38 4.38 4.38 4.38 4.38 8.22

1.739 1.792 1.747 1.840 1.899 1.769 1.899 1.899 1.781 1.774 1.733 1.754 1.848 1.776 1.869

0.63 2.47 2.47 0.63 1.32 2.47 1.32 1.32 1.32 2.47 2.47 2.47 2.47 1.32 4.38

2.085 2.174 2.079 2.219 2.243 2.134 2.242 2.244 2.116 2.108 2.070 2.101 2.220 2.021 2.264

0.02 0.63 0.33 0.70 0.33 0.63 0.33 0.33 0.63 0.33 0.33 0.33 0.63 0.33 0.63

2.716 3.003 2.753 3.002 2.918 2.902 2.916 2.926 2.784 2.784 2.756 2.762 2.921 2.672 3.011

Setting (6) Method

90% C.I. (SE = 0.0030) Left Right Coverage (in %) (in %) (in %)

Exp. width

95% C.I. (SE = 0.0022) Left Right Coverage (in %) (in %) (in %)

Exp. width

99% C.I. (SE = 0.0010) Left Right Coverage (in %) (in %) (in %)

Exp. width

A1 A2, M2 A3 A4 M1, G1–G4 M3 M4 M5 M6 M7 M8 M9 M10 B1 B2

3.01 5.25 4.06 3.01 4.06 5.25 4.06 4.06 3.01 5.25 5.25 5.25 5.25 4.06 5.25

3.887 3.902 3.880 3.982 4.049 3.893 4.049 4.049 3.936 3.887 3.877 3.886 4.002 3.932 3.902

1.59 3.01 2.21 1.85 2.21 3.01 2.21 2.21 2.21 2.21 2.21 2.21 3.01 2.21 3.01

4.622 4.665 4.624 4.771 4.793 4.650 4.792 4.793 5.096 4.637 4.626 4.634 4.768 4.652 4.665

0.11 0.55 0.34 2.14 0.34 0.55 0.34 0.34 0.25 0.34 0.34 0.34 0.55 0.34 0.55

6.067 6.191 6.088 6.350 6.258 6.155 6.254 6.258 6.124 6.103 6.090 6.089 6.240 6.059 6.200

6.59 2.95 6.59 6.59 2.95 2.95 2.95 2.95 6.59 6.59 2.95 2.95 2.95 6.59 2.95

6.49 4.76 4.76 6.49 4.76 4.76 4.76 4.76 6.49 4.76 4.76 4.76 3.36 6.49 4.76

90.94 92.67 89.03 92.09 92.67 92.67 92.67 92.67 90.94 89.03 92.67 92.67 92.67 89.03 88.83

90.50 89.99 91.18 90.50 91.18 89.99 91.18 91.18 90.50 89.99 89.99 89.99 91.39 89.45 89.99

6.59 1.02 2.95 2.97 1.02 1.02 1.02 1.02 2.95 2.95 1.02 2.95 1.02 2.95 2.94

3.36 2.32 2.32 2.32 2.32 2.32 2.32 2.32 2.32 2.32 2.32 2.32 1.65 3.36 2.32

92.78 96.51 94.58 96.40 97.66 96.51 97.66 97.66 95.73 94.58 96.51 94.58 96.51 95.73 92.68

95.05 94.67 95.47 95.83 95.47 94.67 95.47 95.47 95.47 95.47 95.47 95.47 95.34 94.43 94.67

2.95 0.00 0.25 1.02 0.25 0.05 0.25 0.25 0.25 0.25 0.00 1.02 0.25 1.02 0.25

0.98 0.31 0.46 0.14 0.46 0.31 0.46 0.46 0.63 0.46 0.31 0.46 0.31 0.63 0.46

97.03 99.37 99.42 98.28 99.42 99.32 99.42 99.42 99.12 99.42 99.67 98.65 99.12 98.65 99.12

98.91 99.14 99.20 97.72 99.20 99.14 99.20 99.20 99.12 99.20 99.35 99.20 99.14 99.03 98.99

increasing with the value of w ∗ . In sum, if w ∗ = max(w1 , . . . , wk ) (G1) is used, the resulting confidence intervals are always conservative (this was the most conservative procedure among the ones considered) but with wider confidence interval on average. The procedure with w ∗ = min(w1 , . . . , wk ) (G2) is a less conservative one (it can be liberal sometimes), but its expected width was the smallest among G1–G3. Perhaps taking w ∗ as the average of the maximum and minimum weights (G3) or the modification of Tiwari et al. (2006) (G4), which will give a compromise between the coverage probability and the expected width of the interval, is the better option. When we compared methods G3 and G4, we found that they performed similarly in Settings (1), (2) and (3) while G3 has slightly better probability coverages in Settings (7), (8) and (9) in which the variation of the weights wi was large. However, we should keep in mind that G4 has a nice reasoning and it is motivated by continuity correction which used all the values of the weight. B-group: there is not much difference on using the adjusted version (B2) compared to the unadjusted one (B1). Therefore, both B1 and B2 will be compared with other groups in the subsequent discussions. Comparison across groups: For settings (1)–(3), when coverage probabilities for procedures A3, M7, M8, G3, G4, B1 and B2 were investigated, M8, G3 and G4 were shown to be generally robust and in some cases conservative, while A3, M7, B1 and B2 were found to be robust most of the time and liberal in a few instances (see, for example, in Table 3, 90% confidence intervals under setting (5); the coverage probabilities for procedure A3, M7 and B1 are 89.03%, 89.03%, 89.03% and 88.83%, respectively). Expected widths of the confidence intervals were quite similar for A3 and M8 while G3 and G4 produced wider confidence intervals in comparison.

3511

H.K.T. Ng et al. / Computational Statistics and Data Analysis 52 (2008) 3501–3516 Table 4 Empirical coverage probabilities for weighted sums of Poisson parameters, based on 10,000 simulations for settings (7) and (8) Setting (7)

90% C.I. (SE = 0.0030) Right Coverage (in %) (in %)

Method

Left (in %)

A1 A2 A3 A4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 G1 G2 G3 G4 B1 B2

0.84 5.39 3.47 0.56 1.50 2.98 3.45 1.50 1.50 1.90 2.07 2.07 2.46 3.94 2.07 2.07 2.07 2.07 2.07 3.99

Setting (8)

90% C.I. (SE = 0.0030)

Method

Left (in %)

A1 A2 A3 A4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 G1 G2 G3 G4 B1 B2

1.37 5.27 3.76 0.95 2.15 3.29 3.29 2.15 2.15 2.11 2.77 2.77 3.01 3.76 2.75 2.75 2.75 2.75 2.75 4.06

13.70 13.70 13.70 13.70 13.70 13.70 13.70 13.70 13.70 13.70 13.70 13.70 13.70 13.70 0.00 13.70 13.16 13.70 13.70 13.70

Right. (in %). 10.95 9.20 9.28 9.34 9.28 9.29 9.29 9.28 9.28 9.38 9.29 9.29 9.29 9.26 1.94 9.31 8.22 9.26 9.35 9.26

85.46 80.91 82.83 85.74 84.80 83.32 82.85 84.80 84.80 84.40 84.23 84.23 83.84 82.36 97.93 84.23 84.77 84.23 84.23 82.31

Exp. width

95% C.I. (SE = 0.0022) Left Right Coverage (in %) (in %) (in %)

Exp. width

99% C.I. (SE = 0.0010) Left Right Coverage (in %) (in %) (in %)

Exp. width

3.645 4.004 3.680 4.069 3.987 3.735 3.685 3.986 3.986 3.700 3.692 3.654 3.682 3.839 4.365 3.546 4.001 3.666 3.698 3.936

0.15 2.46 1.43 0.12 0.56 1.43 1.43 0.56 0.56 0.82 0.84 1.00 1.00 1.30 1.13 1.13 1.13 1.13 1.13 1.48

4.368 5.096 4.491 5.064 4.734 4.566 4.498 4.734 4.737 4.456 4.449 4.401 4.430 4.656 5.149 4.280 4.626 4.393 4.419 5.037

0.00 0.75 0.29 0.00 0.09 0.15 0.26 0.09 0.09 0.13 0.14 0.14 0.14 0.14 0.20 0.20 0.20 0.20 0.20 0.33

5.751 7.602 6.113 7.127 6.185 6.303 6.129 6.183 6.200 5.916 5.917 5.875 5.863 6.179 6.756 5.774 5.730 5.866 5.764 7.564

13.70 13.70 13.70 13.70 13.70 13.70 13.70 13.70 13.70 13.70 13.70 13.70 13.70 13.70 0.00 13.70 9.51 13.70 13.70 13.70

86.15 83.84 84.87 86.18 85.74 84.87 84.87 85.74 85.74 85.48 85.46 85.30 85.30 85.00 98.87 85.17 89.36 85.17 85.17 84.82

Exp. width

Left (in %)

87.68 85.53 86.96 89.71 88.57 87.42 87.42 88.57 88.57 88.51 87.94 87.94 87.70 86.98 95.31 87.94 89.03 87.99 87.90 86.68

5.078 5.282 5.064 5.478 5.371 5.100 5.075 5.371 5.370 5.124 5.070 5.043 5.063 5.359 5.828 4.982 5.389 5.100 5.066 5.212

0.40 2.51 1.63 0.23 0.82 1.35 1.48 0.82 0.82 0.82 1.04 1.04 1.04 1.33 1.26 1.26 1.26 1.26 1.26 1.64

Right (in %) 9.30 1.91 8.59 8.59 9.00 9.00 9.00 9.00 9.00 9.20 9.20 9.00 9.20 8.39 1.82 9.00 3.49 8.95 9.15 8.59

86.30 85.55 86.01 86.30 86.21 86.15 86.04 86.21 86.21 86.17 86.16 86.16 86.16 86.16 99.80 86.10 99.80 86.11 86.10 85.98

99% C.I. (SE = 0.0010)

95% C.I. (SE = 0.0022)

Coverage (in %)

13.70 13.70 13.70 13.70 13.70 13.70 13.70 13.70 13.70 13.70 13.70 13.70 13.70 13.70 0.00 13.70 0.00 13.69 13.70 13.69

Coverage (in %)

Exp. width

Left (in %)

Right (in %)

Coverage (in %)

Exp. width

90.30 95.58 89.78 91.18 90.18 89.65 89.52 90.18 90.18 89.98 89.76 89.96 89.76 90.28 96.92 89.74 95.25 90.15 89.59 89.77

6.050 6.303 6.105 6.728 6.381 6.173 6.125 6.380 6.383 6.125 6.097 6.061 6.088 6.285 6.876 5.986 6.272 6.089 6.098 6.710

0.04 0.71 0.33 0.01 0.12 0.28 0.30 0.12 0.12 0.13 0.14 0.14 0.13 0.16 0.23 0.23 0.23 0.23 0.23 0.33

8.22 1.91 1.91 2.16 2.04 1.91 1.91 2.04 2.04 3.49 2.64 1.91 3.49 1.94 0.00 1.91 1.91 1.91 1.91 1.91

91.74 97.38 97.76 97.83 97.84 97.81 97.79 97.84 97.84 96.38 97.22 97.95 96.38 97.90 99.77 97.86 97.86 97.86 97.86 97.76

7.925 8.899 7.924 9.057 8.057 8.026 7.924 8.056 8.067 7.846 7.813 7.756 7.810 8.005 8.865 7.700 8.182 7.806 7.429 8.864

For setting (4)–(6) when all weights are equal, methods M1 (which is equivalent to G1–G4), M6 and M10 produce exact confidence intervals which are conservative and have relatively larger expected width. Since most of the methods considered here produce similar and satisfactory results in terms of their coverage probabilities, we compare these methods based on their expected widths. Among all the methods which have simulated coverage probability always above the nominal level, method M8 gives the shortest expected width in most cases. For setting (7)–(9), the variation between the weights wi was large (Var(wi ) = 0.0817) and the values of θi was small (θi = 2). G1 was the only procedure that maintained the coverage probabilities above the nominal level for all 90%, 95% and 99% confidence intervals under these conditions. This finding is in agreement with the conclusion

3512

H.K.T. Ng et al. / Computational Statistics and Data Analysis 52 (2008) 3501–3516

Table 5 Empirical coverage probabilities for weighted sums of Poisson parameters, based on 10,000 simulations for setting (9) Setting (9) Method

90% C.I. (SE = 0.0030) Left Right Coverage (in %) (in %) (in %)

Exp. width

95% C.I. (SE = 0.0022) Left Right Coverage (in %) (in %) (in %)

Exp. width

99% C.I. (SE = 0.0010) Left Right Coverage (in %) (in %) (in %)

Exp. width

A1 A2 A3 A4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 G1 G2 G3 G4 B1 B2

2.30 5.05 4.16 1.75 2.90 3.90 3.90 2.90 2.90 2.78 3.55 3.55 3.58 4.42 3.52 3.52 3.52 3.52 3.52 4.16

7.978 8.038 7.922 8.372 8.240 7.953 7.938 8.240 8.240 8.023 7.941 7.918 7.937 8.137 8.711 7.903 8.267 8.009 7.926 8.065

0.74 2.56 1.75 0.58 1.09 1.55 1.55 1.09 1.09 1.00 1.19 1.22 1.22 1.56 1.44 1.44 1.44 1.44 1.44 1.75

9.474 9.663 9.443 10.129 9.748 9.491 9.472 9.747 9.748 11.111 9.463 9.427 9.456 9.688 10.263 9.403 9.771 9.501 9.380 9.691

0.08 0.64 0.38 0.05 0.17 0.30 0.31 0.17 0.17 0.19 0.21 0.21 0.20 0.26 0.30 0.30 0.30 0.30 0.30 0.39

12.353 13.008 12.415 13.657 12.653 12.494 12.434 12.646 12.653 12.413 12.376 12.326 12.359 12.607 13.291 12.303 12.713 12.407 12.167 13.020

9.87 6.05 7.09 8.08 7.38 7.58 7.58 7.38 7.38 8.81 7.89 7.58 7.76 6.83 3.83 7.85 6.07 7.24 8.08 7.24

87.83 88.90 88.75 90.17 89.72 88.52 88.52 89.72 89.72 88.41 88.56 88.87 88.66 88.75 92.65 88.63 90.41 89.24 88.40 88.60

6.69 3.12 4.04 5.06 4.89 4.65 4.89 4.89 4.89 3.32 5.23 4.65 5.23 4.06 2.11 4.89 3.34 4.20 5.06 4.20

92.57 94.32 94.21 94.36 94.02 93.80 93.56 94.02 94.02 95.68 93.58 94.13 93.55 94.38 96.45 93.67 95.22 94.36 93.50 94.05

3.20 0.84 1.07 1.50 1.86 1.50 1.50 1.86 1.86 2.33 2.20 1.50 2.32 1.50 0.46 1.41 0.97 1.18 1.50 1.18

96.72 98.52 98.55 98.45 97.97 98.20 98.19 97.97 97.97 97.48 97.59 98.29 97.48 98.24 99.24 98.29 98.73 98.52 98.20 98.43

regarding the conservativeness of the procedures presented in the G-group and the observations of Fay and Feuer (1997). Conceivably, when the variation between the weights is large (say, > 0.05) and the values of θi are small (say ≤ 2), G1 can be considered as the best procedure to use in terms of coverage probability. A further investigation of the performance of the confidence intervals under different values of Var(wi ) will be done in the next subsection. 3.2. Randomly generated scenarios In order to investigate the performance of the confidence intervals Pk in different situations, we randomly generate the values of θi from a uniform distribution and rescale them to get i=1 θi = 20 and randomly generate the weights wi , Pk i = 1, . . . , k from a uniform distribution and standardize them such that i=1 wi = 1. (Note that the simulation study P18 in Fay and Feuer (1997) generated the weights and standardized them such that i=1 wi = 18 and the Var(weights) 2 ). We have also considered in their figures are corresponding to the Var(w ) that we are considering multiplied by 18 i Pk i=1 θi = 40 to examine how the larger incidence rates in each subgroup affect the performance of the interval estimation and we found that the results are similar. For the sake of saving space, we only presented the results for Pk i=1 θi = 20 here. We consider k = 6 and 95% confidence level. This process is repeated to obtain 500 different sets of randomly generated θi ’s and wi ’s and 1,000 simulations are used to estimate the coverage probabilities and the conditional expected widths for different methods. Since the performances of the confidence intervals depend on the variance of the weights and in most real applications the weights are known (and so is the variance of the weights), we study the performance of the confidence intervals for large (say > 0.05), moderate (say 0.01 to 0.05) and small (say < 0.01) values of Var(wi ) and make our recommendations based on these results. It is noteworthy that the variations between weights (wi = ci /n i ) are usually very small ( 0.01) in epidemiologic studies since epidemiologic studies usually deal with large values of n i . Following Tang and Ng (2004), we present the mean coverage probability (MCP), minimum coverage probability (MinCP), mean left-tail error rate (MLTER), mean righttail error rate (MRTER), proportion Pkof conservativeness (PC), proportion of liberality (PL) over of the 500 sets of randomly generated scenarios for i=1 θi = 20 in Tables 6–8, where PC and PL are defined as PC =

500 1 X I (estimate coverage probability for the j-th generated scenarios ≥ 0.96) 500 j=1

3513

H.K.T. Ng et al. / Computational Statistics and Data Analysis 52 (2008) 3501–3516 Table 6 Simulated coverage properties (in %) of various 95% confidence intervals with

P6

i=1 θi = 20 for Var(wi ) < 0.01

Method

MCP

MinCP

MLTER

MRTER

PC

PL

MEW

A1 A2 A3 A4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 G1 G2 G3 G4 B1 B2

93.5 95.1 94.7 95.1 95.9 95.1 94.9 95.9 95.9 95.1 94.9 95.6 94.8 95.2 96.7 94.8 95.8 95.5 94.3 94.4

89.4 92.8 91.6 91.6 92.9 92.1 92.1 92.9 92.9 92.2 91.9 92.6 91.9 92.6 94.6 92.2 93.9 93.4 90.8 92.2

1.0 2.9 2.2 1.0 1.8 2.8 2.9 1.8 1.8 1.9 2.3 2.3 2.3 2.8 1.9 1.9 1.9 1.9 1.9 3.1

5.5 2.0 3.1 3.9 2.3 2.1 2.3 2.3 2.3 3.0 2.9 2.1 2.9 1.9 1.4 3.3 2.3 2.6 3.8 2.5

0.0 13.4 5.4 12.8 49.0 12.4 8.2 49.2 49.6 12.8 8.0 30.4 7.2 18.4 86.0 9.2 43.6 29.2 1.4 1.6

70.0 9.0 20.0 6.8 1.6 7.8 14.0 1.4 1.4 8.6 14.2 2.6 16.6 6.2 0.0 15.4 0.6 2.0 36.4 31.4

3.492 3.570 3.488 3.689 3.683 3.561 3.531 3.683 3.685 3.524 3.508 3.480 3.502 3.631 3.804 3.508 3.639 3.609 3.453 3.447

Table 7 Simulated coverage properties (in %) of various 95% confidence intervals with

P6

i=1 θi = 20 for 0.01 ≤ Var(wi ) ≤ 0.05

Method

MCP

MinCP

MLTER

MRTER

PC

PL

MEW

A1 A2 A3 A4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 G1 G2 G3 G4 B1 B2

92.9 94.8 94.3 94.6 95.4 94.6 94.4 95.4 95.4 94.5 94.3 95.1 94.3 94.8 96.9 93.9 95.7 95.0 93.8 94.2

75.3 75.0 75.7 77.0 77.0 76.2 76.2 77.0 77.0 76.5 76.4 76.8 76.4 76.4 95.0 75.6 90.9 79.0 75.6 78.0

0.9 2.9 2.1 0.8 1.6 2.6 2.6 1.6 1.6 1.7 2.1 2.1 2.1 2.6 1.8 1.8 1.8 1.8 1.8 2.5

6.3 2.3 3.6 4.6 3.0 2.8 2.9 2.9 2.9 3.8 3.6 2.8 3.6 2.6 1.3 4.3 2.6 3.3 4.4 3.3

0.0 12.8 6.4 9.0 34.6 9.6 5.8 34.8 35.0 10.0 5.6 21.0 4.4 12.4 86.2 1.2 31.4 16.4 1.0 2.2

79.8 15.2 28.2 18.4 8.4 20.8 26.4 8.4 8.4 23.2 29.8 11.4 30.0 15.6 0.0 45.6 1.8 14.4 50.0 36.6

5.236 5.373 5.228 5.566 5.516 5.331 5.287 5.516 5.518 5.279 5.255 5.212 5.245 5.436 5.746 5.179 5.416 5.317 5.166 5.158

PL =

500 1 X I (estimate coverage probability for the j-th generated scenarios ≤ 0.94) 500 j=1

where I (A) is the indicator function of event A. Most statisticians and practitioners want guaranteed coverage with the shortest width confidence interval and some want central confidence intervals (Fay and Feuer, 1997) while others will allow non-central intervals in order to get a shorter width (Casella and Robert, 1989); thus we will compare the performance of the confidence interval

3514

H.K.T. Ng et al. / Computational Statistics and Data Analysis 52 (2008) 3501–3516

Table 8 Simulated coverage properties (in %) of various 95% confidence intervals with

P6

i=1 θi = 20 for Var(wi ) > 0.05

Method

MCP

MinCP

MLTER

MRTER

PC

PL

MEW

A1 A2 A3 A4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 G1 G2 G3 G4 B1 B2

90.9 93.3 92.7 92.9 93.5 92.8 92.6 93.5 93.5 92.6 92.4 93.2 92.4 93.0 97.0 92.1 95.4 94.2 92.0 93.0

42.4 39.7 39.7 42.4 42.4 42.2 42.2 42.4 42.4 42.2 42.2 42.2 42.2 42.2 94.0 39.7 84.6 60.7 39.7 69.2

0.8 2.9 2.1 0.6 1.5 2.4 2.5 1.5 1.5 1.6 1.9 2.0 2.0 2.4 1.8 1.8 1.8 1.8 1.8 2.3

7.3 2.8 4.2 5.5 3.9 3.7 3.9 3.9 3.9 4.8 4.6 3.7 4.6 3.6 1.3 5.1 2.7 3.0 5.2 4.0

0.2 16.0 6.0 9.2 26.2 8.6 4.6 26.0 26.4 5.6 4.6 17.0 4.4 9.0 82.6 1.4 28.6 7.0 1.4 3.0

82.6 25.4 40.8 37.2 27.2 41.8 46.0 27.2 27.2 43.8 49.2 32.8 50.0 36.0 0.2 58.0 11.2 25.6 58.4 52.0

8.654 8.983 8.643 9.282 9.106 8.799 8.725 9.106 9.110 8.718 8.675 8.602 8.659 8.972 9.453 8.506 8.870 9.459 8.502 8.481

construction procedures in terms of (i) their coverage probabilities (MCP, MinCP, PC and PL); (ii) their expected width (MEW); and (iii) their symmetric of left-tail and right-tail error rates (MLTER and MRTER). For large values of the variance of weights (say ≥ 0.05), even though most of the procedures have MCP close to or above the desired nominal level 95%, their MinCP can be as low as 39.7%. Therefore, one should use any of these procedures to construct confidence intervals for the standardized incidence rate with great caution when the variance of the weights is large. The procedure which gives a reasonable result (MinCP > 90%) is G1; as a result, we would recommend the use of G1 for the situation when the variance of weights is large, > 0.05. For moderate values of the variance of weights (say 0.01–0.05), similar to the case for large values of Var(wi ), although most of the procedures have MCP close to or above the desired nominal level 95%, their MinCP can be lower than 90% and the procedures which give reasonable results (MinCP > 90%) are G1 and G3. G1 has better MCP and MinCP, but larger MEW compared to G3. Since the differences in MCP and MinCP between G1 and G3 are not much but G3 gives smaller MEW (about 6% shorter), we would recommend the use of G3 for moderate values of the variance of weights. For small values of the variance of weights (say, < 0.01), on the basis of the simulation results in Table 6, we observe a pattern in the performance of the procedures similar to the one we obtained in Section 3.1. First, when we compared the coverage probabilities of the procedures, we found that most of the procedures have MCP close to or above 95%. As we mentioned before, a conservative procedure is of less concern compared to a liberal procedure; therefore, we would prefer those procedures with high value of PC and small values of PL. This gives the top seven candidates to be procedures M1, M4, M5, M8, G1, G3 and G4. All of these seven procedures have MCP > 95% but the procedures in the G-group have larger MinCP compared to those in the M-group. Second, when we compare the procedures in terms of expected width, among M1, M4, M5 and M8, M8 gives the smallest MEW (MEW = 3.480); among G1, G3 and G4, G4 gives the smallest MEW (MEW = 3.609) while G3 gives a very close result (MEW = 3.639). Third, when we compare the symmetric of left-tail and right-tail error rates of the three procedures M8, G3 and G4, M8 is the most symmetric one and G4 is the least symmetric one. Although the performances of G3 and G4 are similar, as we mentioned before, G4 is well motivated by continuity correction and used all the weights in the computation. Therefore, in conclusion, G4 is the desirable candidate when the coverage probabilities of the confidence intervals are preferred to be above the nominal level. Otherwise, M8 provides a reliable alternative with its coverage probability being only occasionally slightly lower than the nominal level. M8 provides shorter expected length than G4 but does not require computer software to evaluate the probability distribution of a chi-squared distribution with non-integer degree of freedom.

3515

H.K.T. Ng et al. / Computational Statistics and Data Analysis 52 (2008) 3501–3516 Table 9 Incidence rates of myocardial infarction in women in Augsburg in 1986 Age group

Person-year, n i

Event, X i

Weight, ci

35–39 40–44 45–49 50–54 55–59 60–64

7971 7084 9291 7743 7798 8809

0 0 1 2 4 10

6/31 6/31 6/31 5/31 4/31 4/31

Table 10 95% confidence intervals for age-standardized incidence rates per 10,000 based on different procedures and the width of the confidence intervals Method

Confidence interval

Width

A1 A2 A3 A4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 G1 G2 G3 G4 B1 B2

1.4306, 4.0726 1.7025, 4.4471 1.6307, 4.2952 1.4512, 4.2530 1.5915, 4.4220 1.7002, 4.4428 1.7078, 4.4233 1.5923, 4.4225 1.5906, 4.4222 1.5416, 4.2728 1.6456, 4.3317 1.6530, 4.3214 1.6500, 4.3306 1.7106, 4.5309 1.5931, 4.6080 1.5931, 4.3981 1.5931, 4.4968 1.5931, 4.4320 1.5931, 4.2213 1.7805, 4.4319

2.6420 2.7446 2.6645 2.8018 2.8305 2.7426 2.7156 2.8303 2.8316 2.7312 2.6861 2.6684 2.6805 2.8204 3.0148 2.8049 2.9037 2.8388 2.6282 2.6514

In general, regardless of the variation of the weights, if one wants a symmetric confidence interval which maintains nominal coverage, we would recommend the use of method G1. However, once again, the trade-off is that method G1 has a larger expected width. 4. Illustrative example In this section we illustrate the confidence interval construction procedures discussed in this paper using data from the World Health Organization Monitoring Trends and Determinants in Cardiovascular Diseases (WHO MONICA, (Dobson et al., 1991)) project. The incidence rates for non-fatal myocardial infarction in women aged 35–64 years stratified by 5-year age group in 1986 are presented in Table 9. The weights ci correspond to the Segi World Standard Population (Waterhouse et al., 1976). P6 ci For this data set, the MLE of the age-standardized incidence rates per 10,000 is (10,000 × i=1 n i X i ) = 2.75158 and the resulting 95% confidence intervals for age-standardized incidence rates per 10,000 people based on different procedures and their widths are listed in Table 10. Note that confidence intervals A1, M1 and M6 correspond to the results of ‘crude normal approximation’, ‘chi-squared method’ and ‘Crow and Gardner method’ presented in Table III of Dobson et al. (1991). Findings presented here are consistent with those obtained through the Monte Carlo simulation. Specifically, M8 provided short confidence intervals.

3516

H.K.T. Ng et al. / Computational Statistics and Data Analysis 52 (2008) 3501–3516

5. Conclusions In this article we presented a comparative study of some common non-iterative confidence interval construction methods (asymptotic, the method based on matching moments, the method based on gamma distribution and the method based on beta distribution) as regards the estimation of direct standardized incidence (or mortality) rates. The Monte Carlo simulation showed that the procedures based on gamma distribution maintained the coverage probabilities above the desired nominal level. Nonetheless, this desirable feature was somewhat mitigated by the larger than average width of the confidence intervals. When all weights are equal or when some the weights are equal and the rest are zero, methods M1 (which is equivalent to G1–G4), M6 and M10 produce exact confidence intervals which are conservative but have larger expected width. In this case, M8 is recommended since it gives the shortest expected width among all the methods which can control the coverage probability above the nominal level. w +w In summary, the method based on gamma distribution with w ∗ = (1) 2 (k) (G3) and the version with w ∗ = w(k) (G1) are recommended for moderate and large values of Var(wi ), respectively, while the other procedures should be used with caution as they can produce extremely low coverage probabilities. For small values of Var(wi ) (say < 0.01), we would recommend the method proposed by Tiwari et al. (2006) based on gamma distribution (G4) when the coverage probabilities of the confidence intervals are preferred to be above the nominal level and the matching moments method with mid- p approximation confidence interval (M8) as the best alternative when both coverage probability and average length are of interest. In general, regardless of the variation of the weights, if one is not concerned about the length of interval and wants a symmetric confidence interval which maintains nominal coverage, we would recommend the use of method G1. Acknowledgments We express our sincere thanks to the editor, the associate editor and the two anonymous reviewers for their comments. The authors would also like to thank Dr. Ram C. Tiwari for his useful suggestions which greatly improved this article. References Casella, G., Robert, C., 1989. Refining Poisson confidence intervals. The Canadian Journal of Statistics 17, 45–57. Cohen, G.R., Yang, S.-Y., 1994. Mid- p confidence intervals for the Poisson expectation. Statistics in Medicine 13, 2189–2203. Crow, E.L., Gardner, R.S., 1959. Confidence intervals for the expectation of a Poisson variable. Biometrika 46, 441–453. Dobson, A.J., Kuulasmaa, K., Eberle, E., Scherer, J., 1991. Confidence intervals for weighted sums of Poisson parameters. Statistics in Medicine 10, 457–462. Fay, M.P., Feuer, E.J., 1997. Confidence intervals for directly standardized rates: A method based on the gamma distribution. Statistics in Medicine 16, 791–801. Fleiss, J.L., Levin, B., Paik, M.C., 2003. Statistical Methods for Rates and Proportions, third ed. John Wiley & Sons, Hoboken, New Jersey. Gail, M.H., Benichou, J., 2000. Encyclopedia of Epidemiologic Methods. John Wiley & Sons, Chichester. Garwood, F., 1936. Fiducial limits for the Poisson distribution. Biometrika 28, 437–442. Hall, P., 1992. The Bootstrap and Edgeworth Expansion. Springer-Verlag, New York. Hirji, K.F., 2006. Exact Analysis of Discrete Data. Chapman & Hall/CRC, Boca Raton. Johnson, N.L., Kemp, A., Kotz, S., 2005. Univariate Discrete Distributions, third ed. John Wiley & Sons, New York. Johnson, N.L., Kotz, S., Balakrishnan, N., 1994. Continuous Univariate Distributions Vol. 1, Second ed. John Wiley & Sons, New York. Kabaila, P., Byrne, J., 2001. Exact short Poisson confidence intervals. The Canadian Journal of Statistics 29, 99–106. Rothman, K.J., Boice, J.D., 1979. Epidemiologic Analysis with a Programmable Calculator. National Institutes of Health, Washington. Rothman, K.J., Greenland, S., 1998. Modern Epidemiology, second ed. Lippincott-Raven, Philadelphia. Selvin, S., 1996. Statistical Analysis of Epidemiologic Data, second ed. Oxford University Press, Oxford. Stevens, W.L., 1957. Shorter intervals for the parameter of the binomial and Poisson distributions. Biometrika 44, 436–440. Tang, M.L., Ng, H.K.T., 2003. Letter to the Editor, comments on “Confidence limits for the ratio of two rates based on likelihood scores: Noniterative method” by Graham P.L., Mengersen K., Morton A.P. Statistics in Medicine 22, 2071–2083. Statistics in Medicine 23, 685–693. Tiwari, R.C., Clegg, L.X., Zou, Z., 2006. Efficient interval estimation for age-adjusted cancer rates. Statistical Methods in Medical Research 15, 547–569. Vandenbroucke, J.P., 1982. A shortcut method for calculating the 95 percent confidence interval of the standardized mortality ratio. (Letter). American Journal of Epidemiology 115, 303–304. Waterhouse, J., Muir, C., Correa, P., Powell, J. (Eds.), 1976. Cancer Incidence in Five Continents, vol. III. International Agency for Research on Cancer, Lyon.