Analysis of proportions for clustered binary data in the presence of unequal dispersions

Journal of Statistical Planning and Inference 141 (2011) 2933–2945

Contents lists available at ScienceDirect

Journal of Statistical Planning and Inference journal homepage: www.elsevier.com/locate/jspi

Analysis of proportions for clustered binary data in the presence of unequal dispersions Krishna K. Saha Department of Mathematical Sciences, Central Connecticut State University, New Britain, CT 06050, United States

a r t i c l e i n f o

abstract

Article history: Received 20 April 2010 Received in revised form 8 March 2011 Accepted 10 March 2011 Available online 16 March 2011

This paper investigates the test procedures for testing the homogeneity of the proportions in the analysis of clustered binary data in the context of unequal dispersions across the treatment groups. We introduce a simple test procedure based on adjusted proportions using a sandwich estimator of the variance of the proportion estimators obtained by the generalized estimating equations approach of Zeger and Liang (1986) [Biometrics 42, 121–130]. We also extend the exiting test procedures of testing the hypothesis of proportions in this context. These test procedures are then compared, by simulations, in terms of size and power. Moreover, we derive the score test for testing the homogeneity of the dispersion parameters among several groups of clustered binary data. An illustrative application of the recommended test procedures is also presented. & 2011 Elsevier B.V. All rights reserved.

Keywords: Dispersion parameter Homogeneity of the proportions Generalized estimating equations Sandwich estimator

1. Introduction There has been considerable interest in analyzing clustered binary data that often exhibit extra-dispersion in many applications (Weil, 1970; Kleinman, 1973; Williams, 1975; Prentice, 1986; Donovan et al., 1994; Gibson and Austin, 1996; Keen, 1996; Paul and Islam, 1998; Ridout et al., 1999; Paul et al., 2003; Zhu et al., 2003; Kuk, 2004; Zou and Donner, 2004; Yu and Zelterman, 2008; Xiao et al., 2010). Standard approaches of analyzing such data are based on the assumption that the binary responses of the individuals within cluster are independent. However, in many real-life applications in biological investigations (Crowder, 1978; Paul, 1982; Otake and Prentice, 1984; Donner and Banting, 1989; Bowman and George, 1995; Lindsey and Altham, 1998), it often shows that these responses within a cluster are likely correlated. For example, in studies where the experimental unit is a litter (cluster), there is a tendency of littermates to respond more alike than animals from different litters (Weil, 1970). Therefore, analyzing such data based on the standard approaches, which ignore the litter effect, would give misleading inferences by producing biased estimates of the parameters and inflating the Type I error rates of the hypothesis testing about the parameters. A number of approaches of analyzing clustered binary data have been developed by taking into account the cluster effect or the within cluster correlations on the proposed parametric models, for example, the beta-binomial model, correlated-binomial model, and multiplicativebinomial models. Due to its simplicity, the beta-binomial model has been widely used. In some situations, the full parametric assumption may be too restrictive in which cases a more flexible model can be used that only specifies the mean and variance of the data distribution. However, in many real-life applications, the specific form of the mean– variance relationship of the data may not be justified (Moore and Tsiatis, 1991). In addition, these parametric and

E-mail address: [email protected] 0378-3758/$ - see front matter & 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.jspi.2011.03.018

2934

K.K. Saha / Journal of Statistical Planning and Inference 141 (2011) 2933–2945

semiparametric models involve a dispersion parameter, which requires numerical methods for its estimation (Paul and Islam, 1998; Lee, 2004). In this paper we introduce a simple test procedure for testing the homogeneity of proportions of several treatment groups of clustered binary data in the presence of unequal dispersions based on adjusted proportions using a sandwich estimator of the variance of the proportion estimators obtained from the generalized estimating equations approach of Zeger and Liang (1986). Paul and Islam (1994, 1995) developed Cða) test procedures based on beta-binomial likelihood, quasi-likelihood and extended quasi-likelihood for testing the equality of proportions with the assumption of common dispersion among several treatment groups. However, this may be too restrictive, which is true for the application discussed in Section 4, where the assumption of common dispersion is rejected. So we extend the test procedures just mentioned with the assumption of unequal dispersion among treatment groups and include them in our investigation. Rao and Scott (1992) analyzed the clustered binary data based on the adjusted proportions using an estimator of the variance of a ratio estimator obtained by applying the concepts of design and sample size. We also include this approach for the unequal dispersion case and compare this with our proposed test through simulations. This paper is organized as follows. The test procedures based on adjusted proportions using the generalized estimating equations approach as well as a ratio estimator approach are developed in Section 2. Section 3 extends the Cða) tests based on the beta-binomial likelihood, the quasi-likelihood, and the extended quasi-likelihood in the context of unequal dispersions. A simulation study is conducted in Section 4. Section 5 provides the score test for testing the equality of the dispersion parameters among several treatment groups of clustered binary data. A real-life example is given in Section 5, and a discussion follows in the last section.

2. The test statistics based on adjusted data 2.1. The test based on generalized estimating equations Let yij ðj ¼ 1, . . . ,mi Þ be a random sample of the number of affected individuals among the nij ðj ¼ 1, . . . ,mi Þ individuals drawn from the ith population, i ¼ 1, . . . ,k, with the expected value of the number of affected individuals EðYij Þ ¼ nij pi , i ¼ 1, . . . ,k, j ¼ 1, . . . ,mi , where pi is the expected proportion of affected individuals in the ith population. In this study we are interested in making inferences about the proportions pi , i ¼ 1, . . . ,k, using the estimators of the proportions obtained based on the generalized estimating equation approach of Zeger and Liang (1986). That is, we wish to test H0 : p1 ¼    ¼ pk ¼ p, where p is unspecified, against H1: at least two p’s are unequal. A usual estimator of pi is the overall sample proportion P i Pmi p^ i ¼ yi: =ni: , where yi: ¼ m variance of p^ i is given by j yij and ni: ¼ j nij . Now, if yi. follows binomial(ni: , pi ), the P varðp^ i Þ ¼ pi ð1pi Þ=ni: . Then, for testing H0 against H1, the usual statistic is given by W ¼ ki ¼ 1 ðyi: ni: p^ Þ2 =½ni: p^ ð1p^ Þ with P P p^ ¼ ki¼ 1 yi: = ki¼ 1 ni: , which is asymptotically distributed as chi-square with k 1 degrees of freedom under H0. Ignoring ^ p^ i Þ may cause under/over-estimating the variance of the proportions, the cluster effects, the estimated binomial variance varð which results too liberal or too conservative behavior of the test statistic W. In order to obtain a consistent estimator of varðp^ i Þ, we use an sandwich estimator based on the generalized estimating equation approach proposed by Zeger and Liang (1986) as follows. Let Yijl be the lth binary response of the jth group in the ith population. Now, let EðYijl Þ ¼ mijl ¼ gðXi bi Þ, where g is the link function, Xi is an nij  p design matrix, and bi is a p  1 vector of regression parameters. Further, let yij ¼ ðyij1 , . . . ,yijnij Þ0 and Rij ðaÞ be the nij  nij working correlation matrix for yij. Then, following Zeger and Liang (1986) bi is estimated by solving the estimating equation: mi X

C0ij U1 ij Lij ¼ 0,

ð2:1Þ

j¼1

where Lij ¼ yij mij , mij ¼ ðmij1 , . . . , mijnij Þ0 , Cij ¼ @mij =@bi , and Uij is the covariance matrix for yij given by

Uij ¼ Gij1=2 Rij ðaÞG1=2 =fi , ij where Gij is the nij  nij diagonal matrix with the variance function uðmill Þ as the lth diagonal element. Then a sandwich estimator of the variance of b^ i is given by 1 ð2:2Þ varðb^ i Þ ¼ D1 i0 Di1 Di0 , Pmi P mi 1 0 0 C0ij U1 where Di0 ¼ j ¼ 1 C0ij U1 ij Cij and Di1 ¼ ij covðyij ÞUij Cij with covðyij Þ ¼ Lij Lij . j¼1 In this case, there is no covariate, EðYijl Þ ¼ mijl ¼ pi ¼ expðb0i Þ=½1þ expðb0i Þ, Cij ¼ pi ð1pi Þ, and Rij ðaÞ is an exchangeable working correlation matrix with diagonal elements as 1 and off-diagonal elements as a. Then Eq. (2.1) becomes mi X

pi ð1pi Þfyij nij pi g ¼ 0,

j¼1

where yij ¼

Pnij

l¼1

yijl , which gives p^ i ¼

Pmi

j¼1

yij =

Pmi

j¼1

nij ¼ yi: =ni: .

K.K. Saha / Journal of Statistical Planning and Inference 141 (2011) 2933–2945

2935

P i P i Also, we obtain Di0 ¼ ½pi ð1pi Þ2 m n and Di1 ¼ ½pi ð1pi Þ2 m fy nij pi g2 . Now using Di0 and Di1 , a sandwich j ¼ 1 ij j ¼ 1 ij ^ estimator of the variance of b 0i is given by Eq. (2.2) as Pmi fy nij pi g2 j ¼ 1 ij ^ : varðb 0i Þ ¼ ½pi ð1pi Þni: 2 Using the delta method, we obtain a sandwich estimator of the variance of p^ i as Pmi   ^ 2 @pi 2 j ¼ 1 fyij nij p i g VGEE ¼ varðb^ 0i Þjpi ¼ p^ i ¼ : 2 @b0i ni: As mi -1,VGEE is a consistent estimator of varðp^ i Þ in the sense that mi ½VGEE varðp^ i Þ converges in probability to zero, for ^ p^ i Þ ¼ p^ i ð1p^ i Þ=ni: is the variance inflation or each i. Notice that the ratio of VGEE to the estimated binomial variance varð deflation factor due to clustering and is given by di ¼ ni: VGEE =½p^ i ð1p^ i Þ. Following Rao and Scott (1992), we now adjust the ith population data of yi: and ni: by yi: ¼ yi: =di and ni: ¼ ni: =di , respectively, and treat yi: as a binomial random variable with parameters ni: and pi . Then, the statistic for testing H0 against H1 is given by TGEE ¼

k X ðy n p^ Þ2 i:

i:

ni: p^ ð1p^ Þ

i¼1

,

which is asymptotically distributed as chi-square with k 1 degrees of freedom under H0. 2.2. The test based on a ratio estimator As we have seen in Section 2.1, the estimator of pi is the overall sample proportion, which can be written as p^ i ¼ yi: =ni: ¼ y i: =n i: , where y i: ¼ yi: =mi and n i: ¼ ni: =mi . Since p^ i is the ratio of the two sample means, using a result provided by Cochran (1954) we obtain an estimator of the variance of p^ i as

WR ¼

mi

Pmi

j¼1

ðyij nij p^ i Þ2

ðmi 1Þn2i:

:

Note that WR is also a consistent estimator of varðp^ i Þ in the sense that mi ½WR varðp^ i Þ converges in probability to zero, as mi -1 for each i. Similar to the previous case, the inflation or deflation factor in this case is di ¼ ni: Wi =p^ i ð1p^ i Þ. The adjusted data of yi. and ni. for population i are, respectively, given by y~ i: ¼ yi: =di and n~ i: ¼ ni: =di , and treat y~ i: as binomial(n~ i: , pi ). Then, the statistic for testing H0 against H1 is given by TR ¼

k X ðy~ i: n~ i: p^ Þ2 , n~ i: p^ ð1p^ Þ

i¼1

which is asymptotically distributed as chi-square with k 1 degrees of freedom under H0 (for a detailed derivation, see Rao and Scott, 1992). 3. The CðaÞ tests based on BB, QL, and EQL 3.1. The CðaÞ test based on the beta-binomial likelihood Let yij ðj ¼ 1, . . . ,mi Þ be a sample from the beta-binomial distribution BBðpi , fi Þ with probability mass function ! Qy 1 Qnij yij 1 ij nij ½ð1pi Þð1fi Þ þr fi  r ¼ 0 ½ð1fi Þpi þr fi  r¼0 Pðyij jpi , fi Þ ¼ Qnij 1 yij ½ð1 f Þ þ r fi  i r¼0 for yij ¼ 0,1,2, . . ., 0 r pi r1 and maxð1=ðni 1ÞÞ o f o 1. The mean and variance of Yij are EðYij Þ ¼ pi and var(Yij Þ ¼ nij pi ð1pi Þf1 þðni 1Þfi g. The parameter fi is called the over-dispersion parameter, which can also be interpreted as an intraclass correlation parameter. With the above range of fi ,Pðyij jpi , fi Þ is a valid probability function (see Prentice, 1986). Obviously, when fi -0, the BBðpi , fi Þ becomes binomial ðpi Þ. In order to develop procedures for testing the equality of proportions of the BBðpi , fi Þ populations, that is, testing H0 : p1 ¼    ¼ pk against H1: at least two p’s are unequal, with the assumption fi af for all i ¼ 1, . . . ,k based on the CðaÞ test, we reparameterize pi , under H1, by pi ¼ p þ ai with ak ¼ 0. Let a0 ¼ ða1 ; :::; ak1 Þ and n0 ¼ ðp, f1 , . . . , fk Þ. Then testing H0 is equivalent to testing H00 : a ¼ 0 against H01 : at least two a’s are unequal, with n0 being treated as nuisance parameters. Under H01 , the log-likelihood, apart from a constant, can be written as "y 1 # nij yij 1 nij 1 mi ij k X X X X X l¼ lnfð1fi Þðp þ ai Þ þ r fi g þ lnfð1ðp þ ai ÞÞð1fi Þ þ r fi g lnfð1fi Þ þ r fi g : i¼1j¼1

r¼0

r¼0

r¼0

2936

K.K. Saha / Journal of Statistical Planning and Inference 141 (2011) 2933–2945

Write UðnÞ ¼

 @l  , @aa ¼ 0

VðnÞ ¼

 2  @ l  CðnÞ ¼ E @a@na ¼ 0

 @l  , @na ¼ 0

and

   @2 l  , AðnÞ ¼ E  0 @a@a a ¼ 0

 2   @ l  DðnÞ ¼ E : @n@n0 a ¼ 0

The dimensions of the matrices U, V, A, C, and D are, respectively, (k  1)  1, (kþ1)  1, (k 1)  (k  1), (k 1)  (kþ1), and (kþ1)  (kþ1). After replacing the nuisance parameter n by its On-consistent estimate n^ under H00 in U, V, A, C, and D, the 0 0 Cða) test (Neyman, 1959), based on the residual LðnÞ ¼ UðnÞb VðnÞI with b ¼ ðb1 , . . . , bk þ 1 Þ and I ¼ ð1, . . . ,1Þ0 being the 0 0 0 column vector of (k 1) order, for testing H0 against H1 is given by Lðn^ Þ ½Aðn^ ÞCðn^ ÞD1 ðn^ ÞC 0 ðn^ ÞLðn^ Þ, which follows a chisquare distribution with k 1 degrees of freedom under H00 . Note that b is the partial regression coefficient of U on V. After some straightforward algebra and by using the maximum likelihood estimate of the nuisance parameter n0 under H00 , the Cða) test becomes TML ¼

k X c2 ðn^ Þ i

i¼1

oi ðn^ Þ

,

where

ci11 ðn^ Þ ¼ ð1f^ i Þ

" y mi ij X X j¼1

r¼1

1 R1iðr1Þ ðn^ Þ

nij yij

X

1

r¼1

R2iðr1Þ ðn^ Þ



# and

oi ðn^ Þ ¼ si11 ðn^ Þs2i12 ðn^ Þ=si22 ðn^ Þ, where ^ Þ2 si11 ðn^ Þ ¼ ð1f i

si12 ðn^ Þ ¼

" n # nij mi ij X X X PðYij ZrÞ PðYij r nij rÞ þ , R2 ðn^ Þ r ¼ 1 R22iðr1Þ ðn^ Þ j ¼ 1 r ¼ 1 1iðr1Þ

" # nij nij mi ^ ÞX X X PðYij Z rÞ PðYij r nij rÞ ð1f i ^ ^ þ ð1p Þ p R2 R22iðr1Þ ðn^ Þ ðn^ Þ f^ i j ¼ 1 r ¼ 1 1iðr1Þ r¼1

and si22 ðn^ Þ ¼

mi 1 X

"

2

f^ i

j¼1

p^ 2

# nij nij nij X X PðYij ZrÞ PðYij r nij rÞ X 1 2 ^ þð1  p Þ , R2 R2 R22iðr1Þ ðn^ Þ ðn^ Þ ðn^ Þ r ¼ 1 1iðr1Þ r¼1 r ¼ 1 3iðr1Þ

where R1ir ðnÞ ¼ pð1fi Þ þ r fi , R2ir ðnÞ ¼ ð1pÞð1fi Þ þ r fi , and R3ir ðnÞ ¼ 1fi þ r fi . The maximum likelihood estimates p^ and f^ i ði ¼ 1 . . . ,kÞ are obtained by solving the maximum likelihood estimating equations: " y # nij yij mi ij k X X X X 1 1  ð1fi Þ ¼0 R ðnÞ r ¼ 1 R2iðr1Þ ðnÞ j ¼ 1 r ¼ 1 1iðr1Þ ii ¼ 1 and " y # nij yij nij mi ij X X ðr1Þð1pÞ X X ðr1Þp r2 þ  ¼ 0, R2iðr1Þ ðnÞ R R ðnÞ ðnÞ r¼1 r ¼ 1 3iðr1Þ j ¼ 1 r ¼ 1 1iðr1Þ

i ¼ 1, . . . ,k,

simultaneously. 3.2. The CðaÞ test based on quasi-likelihood The quasi-log-likelihood using the mean–variance relationship of the counts Yij given above under H01 becomes     p þ ai 1pai þ ðnij yij Þln mi yij ln k X X zij 1zij , Q¼ 1 þ ðnij 1Þfi i¼1j¼1 where zij ¼ yij =nij . Now, given fi ði ¼ 1, . . . ,kÞ, we can easily obtain the unbiased estimating functions for the parameters a1 , . . . , ak1 , p as Ui ða, nÞ ¼ @Q =@ai ,i ¼ 1, . . . ,k1 and V1 ða, nÞ ¼ @Q =@p. However, the estimating function for fi ði ¼ 1; :::; kÞ based on Q does not exist. In such a case, following Breslow (1984, 1990) and Moore and Tsiatis (1991), given a1 , . . . , ak1 , p, an unbiased estimating function for fi ði ¼ 1, . . . ,kÞ can be obtained as V2,t þ 1 ða, nÞ ¼

mt X

ntj ðztj pat Þ2 ðmt 2Þ, ðp þ at Þð1 þ p þ at Þ½1 þ ðnit 1Þft  j¼1

t ¼ 1, . . . ,k:

K.K. Saha / Journal of Statistical Planning and Inference 141 (2011) 2933–2945

2937

We then treat Ui, V1, and V2,t þ 1 as the likelihood score analogs and following the procedure above, we obtain the CðaÞ statistic based on the quasi-likelihood as TQL ¼

k X Z2i ðn^ Þ 1 , p^ ð1p^ Þ i ¼ 1 di ðn^ Þ

where mi X

Zi ðn^ Þ ¼

^ 1 nij ðzij p^ Þ½1þ ðnij 1Þf i

j¼1

and mi X

di ðn^ Þ ¼

^ 1 : nij ½1 þðnij 1Þf i

j¼1

^ ði ¼ 1, . . . ,kÞ are obtained by solving The method of moment estimates p^ and f i mi k X X

nij ðzij pÞ½1þ ðnij 1Þfi 1 ¼ 0

i¼1j¼1

and mi X

nij ðzij pÞ½pð1pÞf1 þ ðnij 1Þfi g1 ¼ ðmi 2Þ,

i ¼ 1, . . . ,k,

j¼1

simultaneously. 3.3. The CðaÞ statistic based on extended quasi-likelihood The extended quasi-likelihood under H01 with the parameters a and n, apart from a constant term, is     2 p þ ai 1pai 3 y ln y Þln þ ðn ij ij ij m k i 6 XX 7 zij 1zij 6 1 lnf1þ ðnij 1Þfi g þ 7: Qþ ¼ 4 2 5 1Þ f 1 þðn ij i i¼1j¼1 Now, based on Q þ one can easily obtain the estimating a1 . . . , ak1 , p, f1 , . . . , fk1 , and fk using Ui ða, nÞ ¼ @Q þ =@ai , i ¼ 1, . . . ,k1; V1 ða, nÞ ¼ @Q þ =@p; and V2,t þ 1 ða, nÞ ¼ @Q þ =@fi , t ¼ 1, . . . ,k. Then, similar to the procedure in Section 3.1 and treating Ui, V1, and V2,t þ 1 as the likelihood score analogs, it can be shown that the CðaÞ statistic, TEQL, based on the extended quasi-likelihood, becomes TQEL ¼

k X Z2i ðn~ Þ 1 , p~ ð1p~ Þ i ¼ 1 di ðn~ Þ

where

Zi ðn~ Þ ¼

mi X

~ 1 , nij ðzij p~ Þ½1þ ðnij 1Þf i

j¼1

di ðn~ Þ ¼

mi X

~ 1 , nij ½1 þðnij 1Þf i

j¼1

~ ; :::; f ~ Þ is the maximum extended quasi-likelihood estimate of n under H . The maximum extended quasiand n~ ¼ ðp~ , f 0 1 k ~ ði ¼ 1, . . . ,kÞ under H can be obtained by solving likelihood estimates p~ and f 0 i mi k X X

nij ðzij pÞ½1þ ðnij 1Þfi 1 ¼ 0

i¼1j¼1

and       zij 1zij 1 þ ðnij yij Þln  f1 þ ðnij 1Þfi g mi ðnij 1Þ yij ln X 2 p 1p j¼1

simultaneously.

½1þ ðnij 1Þfi 2

¼ 0,

i ¼ 1, . . . ,k,

2938

K.K. Saha / Journal of Statistical Planning and Inference 141 (2011) 2933–2945

4. The score test for testing the equality of the dispersion parameters The score test is a special case of the more general CðaÞ test discussed in Section 3.1 in which the nuisance parameters are replaced by their maximum likelihood estimates. We follow the same procedure and notation as in Section 3.1 to derive the score test for testing the equality of the dispersion parameters among the treatment groups. Let H0: f1 ¼    ¼ fk ¼ f against H1: at least two fi ’s are unequal. We rewrite the alternative hypothesis using fi ¼ f þci with ck ¼0. Then testing the null hypothesis of common f is equivalent to testing H0 : ci ¼ 0, for all i, with pi ði ¼ 1, . . . ,kÞ and f treated as nuisance parameters. Now, define a0 ¼ ða1 , . . . , ak1 Þ ¼ ðc1 , . . . ,ck1 Þ, n0 ¼ ðn1 , . . . , nk þ 1 Þ ¼ ðp1 , . . . , pk , fÞ, and ci ðnÞ ¼ ½@l=@ai c ¼ 0 , where l is the beta-binomial likelihood using fi ¼ f þci under H1. The (i, j)th element of A, C, and D are obtained based on the expected mixed partial derivatives similar to those of Section 3.1 and evaluated under H0. We then replace the nuisance parameter n by the maximum likelihood estimate of n and hence we get Li ðn^ Þ ¼ ci ðn^ Þ. After some derivations and simple algebra, we obtain the score test as ST f ¼

k X c2 ðn^ Þ i

i¼1

Zi ðn^ Þ

,

where

ci ðn^ Þ ¼

Zi ðn^ Þ ¼

" y # nij yij nij mi ij X ðr1Þð1p^ i Þ X X X ðr1Þp^ i r2 þ  , R2iðr1Þ ðn^ Þ R R ðn^ Þ ðn^ Þ r¼1 r ¼ 1 3iðr1Þ j ¼ 1 r ¼ 1 1iðr1Þ mi 1 X

f^

2 j¼1

"

p^ 2i

# nij nij nij X X PðYij Z rÞ PðYij r nij rÞ X 1 2 ^ þ ð1  p Þ , i R2 R2 R22iðr1Þ ðn^ Þ ðn^ Þ ðn^ Þ r ¼ 1 1iðr1Þ r¼1 r ¼ 1 3iðr1Þ

^ are the maximum with S1ir ðnÞ ¼ pi ð1fÞ þr f, S2ir ðnÞ ¼ ð1pi Þð1fÞ þr f, and S3ir ðnÞ ¼ 1f þ r f, and p^ i ði ¼ 1, . . . ,kÞ and f likelihood estimates of pi ði ¼ 1, . . . ,kÞ and f under H0, obtained by solving the maximum likelihood estimating equations: " y # nij yij mi ij X X X 1f 1f  ¼ 0, i ¼ 1, . . . ,k S S r ¼ 1 2iðr1Þ j ¼ 1 r ¼ 1 1iðr1Þ and " y mi ij k X X X ðr1Þpi i¼1j¼1

r¼1

R1iðr1Þ

# nij X ðr1Þð1pi Þ X r2  ¼ 0, R2iðr1Þ R r¼1 r ¼ 1 3iðr1Þ

nij yij

þ

simultaneously. 5. Simulations This section reports on a simulation study conducted to investigate the small and moderate sample behavior of the test statistics TML, TQL, TEQL, TR, and TGEE in terms of Type I error rate and power. For simplicity, we consider k¼2 groups and the number of clusters for k¼2 groups are chosen as (i) m1 ¼m2 ¼16; (ii) m1 ¼19 and m2 ¼27; and (iii) m1 ¼ 73 and m2 ¼87. In this study, we consider the cluster sizes for case (i) from the control group (m1 ¼16) and the treatment group (m2 ¼ 16) for the data in Section 3 of Williams (1975). For case (ii), the cluster sizes are used as those of the low dose treatment group (m1 ¼19) and the control group (m2 ¼ 27) of the data in Table 1 of Paul (1982) and for case (iii), the cluster sizes are the same as those of the Dose: 0 group (m1 ¼73) and Dose: 30 group (m2 ¼87) for the data in Table 4. For size comparison, we consider p1 ¼ p2 ¼0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50 and (f1 , f2 )¼ (0.1, 0.2), (0.1, 0.4), (0.2, 0.5), (0.3, 0.6). For power comparison, we consider p1 ¼ p and p2 ¼ p þ d, with p ¼ 0:1 and d ¼0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50. All results are obtained using FORTRAN 90 code. We generate data from the beta-binomial distribution using the IMSL random number generators RNBET and RNBIN. We then compute the empirical Type I error rates and the empirical power as the proportion of the test statistics based on 1000 samples generated from the beta-binomial distribution for the above sets of parameters under the null and alternative hypotheses, respectively. Tables 1–3 present the empirical Type I error rates of all five test statistics corresponding to nominal level a ¼ 0:05. As seen from Tables 1 and 2, Type I error rates for the TR and TGEE statistics are close to the nominal level, but show some conservative behavior for smaller proportions ( o 0:20) and larger dispersions. The statistics TML, TQL, and TEQL in general show some conservative behavior, however, the statistics TML and TQL become extremely conservative for small proportions. For larger number of clusters, we can see from Table 3, the performances of all five test statistics are in general similar and hold the nominal level well. However, the statistics TEQL show extremely liberal behavior for small proportions ( o 0:15). The performances of the statistics TR and TGEE are reasonably good in that both maintain nominal level quite well in almost all data situations studied here. Empirical powers for the number of clusters (i) m1 ¼ m2 ¼16 as well as (ii) m1 ¼19 and m2 ¼27 at a ¼ 0:05 are presented, respectively, in Figs. 1 and 2 and those for the number of clusters (iii) m1 ¼73 and m2 ¼87 at a ¼ 0:05 are presented in

K.K. Saha / Journal of Statistical Planning and Inference 141 (2011) 2933–2945

2939

Table 1 Type I error rates (%) for all five tests TML, TQL, TEQL, TR, and TGEE based on 1000 replications for data generated from beta-binomial with m1 ¼m2 ¼16 for k¼ 2 groups at a ¼ 5% level of significance.

f1

f2

Method

p1 ¼ p2 0.10

0.15

0.20

0.25

0.30

0.35

0.40

0.45

0.50

0.1

0.2

TML TQL TEQL TRS TZL

1.6 1.5 3.3 3.2 3.8

2.4 1.9 3.7 5.1 6.5

3.7 3.7 3.9 5.2 6.0

5.3 4.2 4.2 5.9 6.5

4.7 3.9 3.4 5.4 6.2

5.2 4.2 3.4 5.6 6.7

5.3 4.2 3.1 5.0 5.9

4.6 4.1 3.7 5.4 6.1

5.2 4.2 3.4 5.8 6.8

0.1

0.4

TML TQL TEQL TRS TZL

1.1 3.2 4.9 4.0 4.3

1.7 3.2 3.5 3.6 4.2

3.2 5.1 4.5 5.4 6.3

3.9 4.3 3.6 6.1 7.0

4.1 4.6 3.0 6.2 6.9

5.3 4.6 3.5 5.7 6.3

5.7 5.4 3.9 6.5 7.4

7.0 5.3 3.7 6.6 7.4

5.3 4.7 2.7 5.6 6.1

0.2

0.5

TML TQL TEQL TRS TZL

1.0 0.8 3.4 1.6 2.1

1.6 2.2 3.6 2.8 3.5

1.1 2.5 2.8 3.6 4.2

2.6 3.0 2.6 4.1 4.3

3.7 3.1 2.4 4.5 5.2

4.9 4.3 2.8 5.4 6.4

4.3 3.8 2.2 4.8 5.9

5.4 4.1 2.9 5.4 6.1

5.4 4.4 3.2 6.2 7.5

0.3

0.6

TML TQL TEQL TRS TZL

1.0 1.5 1.6 1.5 1.9

1.4 1.8 2.4 2.2 2.4

1.4 1.7 2.1 2.2 2.6

1.0 1.6 1.7 1.7 2.0

2.7 2.4 1.8 3.2 4.0

4.5 4.1 3.0 5.0 5.7

5.1 4.4 2.9 5.6 6.2

4.2 4.0 2.8 4.5 4.8

4.2 4.1 2.6 4.9 6.0

Table 2 Type I error rates (%) for all five tests TML, TQL, TEQL, TR, and TGEE based on 1000 replications for data generated from beta-binomial with m1 ¼19 and m2 ¼27 for k ¼2 groups at a ¼ 5% level of significance.

f1

f2

Method

p1 ¼ p2 0.10

0.15

0.20

0.25

0.30

0.35

0.40

0.45

0.50

0.1

0.2

TML TQL TEQL TR TGEE

1.6 1.3 4.0 3.9 4.7

3.7 3.2 4.4 5.5 5.9

4.2 4.7 4.1 6.5 7.0

5.5 5.3 4.6 6.9 7.4

4.7 4.1 3.0 5.8 6.0

4.1 3.5 2.7 5.3 6.0

5.4 4.8 4.0 6.8 6.8

5.3 4.6 3.6 5.7 6.4

4.5 4.0 3.1 4.7 5.0

0.1

0.4

TML TQL TEQL TR TGEE

1.2 2.0 4.6 3.4 3.5

1.8 2.4 4.0 5.0 5.8

3.3 3.3 3.5 5.2 5.7

4.4 3.9 3.5 5.5 5.9

4.8 4.4 3.4 5.0 5.3

4.7 4.0 3.3 5.5 5.6

3.6 3.7 2.0 4.3 4.8

4.9 4.5 2.8 5.1 5.6

4.5 4.1 2.2 4.6 5.0

0.2

0.5

TML TQL TEQL TR TGEE

1.8 1.6 3.5 2.6 3.1

1.1 1.8 2.7 3.2 3.6

3.6 3.5 4.2 5.2 5.8

4.2 3.6 3.4 4.5 5.0

4.8 4.1 3.1 5.9 6.2

3.1 2.8 2.0 3.8 4.2

4.3 4.1 2.6 5.0 5.6

4.7 3.8 2.6 5.7 6.1

4.7 3.5 1.8 4.9 5.0

0.3

0.6

TML TQL TEQL TR TGEE

1.5 1.2 2.9 1.8 2.1

1.6 1.9 3.4 2.8 2.9

1.1 1.4 1.9 2.1 2.4

2.6 2.0 2.2 3.1 3.2

3.9 3.0 2.3 4.3 5.2

4.9 3.8 2.5 5.2 5.5

4.9 5.0 3.3 6.1 6.6

5.4 4.4 2.9 5.4 6.0

4.2 4.3 2.1 5.0 5.5

Fig. 3. From Figs. 1 and 2, we see that the power performance of the statistic TGEE is uniformly best. As expected, power of the statistic TR is slightly smaller than that of the statistic TGEE. Power of the statistic TML is somewhat smaller than that of the other four. However, for the larger proportions ( 40:4) the power performances of all five statistics are similar. From Fig. 3, we see that powers of all five statistics are very similar irrespective of the values of the proportions and dispersions. From all three Figs. 1–3, we see that powers of all five statistics become higher when the number of clusters are larger. Note that the results in this limited simulation study are based on data generated from the beta-binomial model and the cluster sizes taken from the real-life data sets that are available in the existing literatures in this area of studies. However, in order to generalize the results above, one could further investigate the performances of these methods more details by generating data from other competing over-dispersed proportion models such as probit-normal binomial and

2940

K.K. Saha / Journal of Statistical Planning and Inference 141 (2011) 2933–2945

Table 3 Type I error rates (%) for all five tests TML, TQL, TEQL, TR, and TGEE based on 1000 replications for data generated from beta-binomial with m1 ¼73 and m2 ¼87 for k¼2 groups at a ¼ 5% level of significance.

f1

p1 ¼ p2

Method

f2

0.10

0.15

0.20

0.25

0.30

0.35

0.40

0.45

0.50

0.1

0.2

TML TQL TEQL TR TGEE

4.3 3.7 5.6 5.0 5.2

5.7 5.5 5.6 5.4 5.5

5.1 5.3 4.4 5.4 5.6

5.2 5.3 4.0 6.0 6.2

5.6 5.5 4.4 5.8 6.0

5.3 4.7 3.4 5.7 5.8

6.0 6.1 4.7 5.8 5.9

4.3 4.2 3.4 4.9 5.0

4.1 4.1 2.7 4.3 4.5

0.1

0.4

TML TQL TEQL TR TGEE

4.6 5.2 8.9 6.6 6.8

5.6 5.3 6.1 5.8 6.0

4.9 4.4 4.2 4.9 5.1

4.1 3.7 3.2 4.2 4.3

4.1 4.2 3.4 4.9 5.0

5.6 5.6 4.0 6.4 6.6

5.5 5.4 3.4 5.8 6.1

4.7 4.9 3.3 5.0 5.1

4.7 4.5 2.5 5.6 5.6

0.2

0.5

TML TQL TEQL TR TGEE

4.0 3.5 8.1 4.5 4.5

5.3 5.2 6.3 5.8 6.1

4.6 4.3 4.5 5.2 5.3

4.8 4.8 3.7 5.4 5.8

5.0 5.0 3.1 5.2 5.4

6.1 5.5 4.2 6.5 6.8

6.1 5.2 3.7 5.3 5.4

6.7 6.0 4.0 5.7 5.8

5.3 5.4 3.2 5.6 5.6

0.3

0.6

TML TQL TEQL TR TGEE

3.0 3.3 6.1 2.7 2.7

4.4 4.0 5.5 4.7 4.7

3.6 3.0 3.1 3.8 3.9

7.4 6.3 5.2 7.0 7.1

4.4 4.4 3.6 5.1 5.2

6.7 4.9 3.4 4.9 4.9

5.1 4.2 2.7 4.3 4.4

6.1 5.0 3.1 5.3 5.3

6.0 6.0 4.1 5.9 5.9

phi1 = 0.1, phi2 = 0.2, pi = 0.1

1.0

0.8

0.8

0.6

0.6

Power

Power

1.0

0.4 0.2

phi1 = 0.1, phi2 = 0.4, pi = 0.1

ML QL EQL R GEE

0.4 0.2

0.0

0.0 0.1

0.2

0.4

0.3

0.5

0.1

0.6

0.2

0.5

0.6

d

1.0 phi1 = 0.2, phi2 = 0.5, pi = 0.1

1.0

0.8

0.8

0.6

0.6

Power

Power

0.4

0.3

d

0.4 0.2

phi1 = 0.3, phi2 = 0.6, pi = 0.1

0.4 0.2

0.0

0.0 0.1

0.2

0.4

0.3 d

0.5

0.6

0.1

0.2

0.4

0.3

0.5

0.6

d

Fig. 1. Empirical power curves corresponding to nominal level a ¼ 5% for the number of clusters m1 ¼16 and m2 ¼16.

K.K. Saha / Journal of Statistical Planning and Inference 141 (2011) 2933–2945

1.0

0.8

0.8

0.6

0.6

Power

Power

phi1 = 0.1, phi2 = 0.2, pi = 0.1 1.0

0.4

0.2

0.0

0.0 0.1

0.2

0.3

0.4

0.5

phi1 = 0.1, phi2 = 0.4, pi = 0.1

ML QL EQL R GEE

0.4

0.2

0.1

0.6

0.2

0.3

phi1 = 0.2, phi2 = 0.5, pi = 0.1

1.0

0.8

0.8

0.6

0.6

Power

Power

0.4

0.5

0.6

d

d

1.0

2941

0.4 0.2

phi1 = 0.3, phi2 = 0.6, pi = 0.1

0.4 0.2

0.0

0.0 0.1

0.2

0.3

0.4 d

0.5

0.6

0.1

0.2

0.3

0.4

0.5

0.6

d

Fig. 2. Empirical power curves corresponding to nominal level a ¼ 5% for the number of clusters m1 ¼ 19 and m2 ¼ 27.

log-normal binomial models following the three steps procedure given by Ochi and Prentice (1984) and using IMSL subroutines CHFAC and RNMVN. Furthermore, the cluster sizes can be generated from the empirical distribution of 523 litter sizes quoted by Kupper et al. (1986) and from the truncated negative binomial distribution with mean 3.12 and variance 4.52, which correspond to the US sibship size distribution in 1950 (Brass, 1958).

6. Application We consider the developmental toxicology study data, originally given and analyzed by Bowman and George (1995). The data in Table 4 consist of response as the number of combined end points of death, resorptions or malformation (y) and the number of implants (n) for each of k¼6 dose groups in a large-scale developmental toxicology study of the herbicide 2,4,5-trichlorophenoxiacetic acid (2,4,5-T) conducted at the US National Center for Toxicological Research. From days 6 to 14 of gestation, pregnant mice from several strains were given daily does of 2,4,5-T and recorded here only the number of implantation sites, fetal deaths, resorptions and cleft palate malformations from the outbreed strain CD-1 for each female mouse. For this outbreed strain CD-1, there were six dose groups corresponding to exposure levels of 0, 30, 45, 60, 75 and 90 mg/kg of 2,4,5-T. Bowman et al. (1995) and Paul and Saha (2007) analyzed these data by fitting a separate beta-binomial distribution for each dose group. In the analysis of the data, we first check the assumption of unequal dispersions among treatment groups, using the score test statistic given in Section 5. The maximum likelihood estimates of the parameters and the values of the score test statistic ST f are given in Table 5, which indicates that the dispersion parameters among the all six dose groups are not homogeneous. The results in Table 5 also indicate that the treatment combinations 12, 13, 14, 15, 16, 23, 24, 25, 26, 34, and 35 have significantly different dispersions, whereas the treatment combinations 36, 45, 46, and 56 have significantly similar dispersions. For testing homogeneity of proportions across all groups or any pair of groups, we now apply the recommended test statistics TR and TGEE. The results are presented in Table 6, which show strong evidence by both statistics that the proportions for all six groups are different. The results in Table 6 also show strong evidence that the proportions for any pair of groups are also different.

2942

K.K. Saha / Journal of Statistical Planning and Inference 141 (2011) 2933–2945

phi1 = 0.1, phi2 = 0.4, pi = 0.1 1.0

0.8

0.8

0.6

0.6

Power

Power

phi1 = 0.1, phi2 = 0.2, pi = 0.1 1.0

0.4

ML QL EQL R GEE

0.4 0.2

0.2 0.0

0.0 0.1

0.2

0.4

0.3

0.5

0.1

0.6

0.2

0.8

0.8

0.6

0.6

Power

Power

1.0

0.4

0.4

0.2

0.2

0.0

0.0 0.4

0.3

0.6

phi1 = 0.3, phi2 = 0.6, pi = 0.1

phi1 = 0.2, phi2 = 0.5, pi = 0.1

0.2

0.5

d

1.0

0.1

0.4

0.3

d

0.5

0.1

0.6

0.2

0.4

0.3

0.5

0.6

d

d

Fig. 3. Empirical power curves corresponding to nominal level a ¼ 5% for the number of clusters m1 ¼73 and m2 ¼87.

Table 4 Distribution of the number of implants (n) and number of combined endpoints (y) following exposure to 2, 4, 5-T. Dose: 0 y1/n1 0/1 0/3 0/5 1/5 0/6 0/7 0/7 0/8 0/9 0/9 0/9 0/9 0/9 1/9 1/9 2/9 0/10 0/10 0/10 0/10 1/10 2/10 0/11 0/11

Dose: 30 y2/n2 1/12 1/12 1/12 2/12 2/12 2/12 0/13 0/13 1/13 1/13 1/13 1/13 1/13 3/13 0/14 0/14 1/14 1/14 2/14 4/14 3/16 1/17 2/17

0/2 0/3 0/4 0/4 0/6 0/6 0/7 0/7 1/7 2/7 2/7 0/8 1/8 2/8 2/8 0/9 0/9 1/9 1/9 6/9 0/10 0/10 0/10 0/10

Dose: 45 y3/n3 1/12 1/12 1/12 1/12 2/12 2/12 3/12 3/12 3/12 4/12 0/13 0/13 0/13 1/13 1/13 1/13 1/13 1/13 2/13 3/13 3/13 4/13 7/13 0/14

0/1 0/3 1/3 0/4 1/5 2/5 2/6 0/7 1/7 0/8 1/8 1/9 1/9 1/9 2/9 4/9 8/9 0/10 1/10 1/10 5/10 6/10 9/10 0/11

Dose: 60 y4/n4 2/12 2/12 2/12 3/12 3/12 4/12 5/12 6/12 6/12 7/12 11/12 12/12 12/12 1/13 2/13 2/13 2/13 3/13 3/13 3/13 5/13 5/13 5/13 8/13

0/1 3/3 3/3 1/4 4/5 5/6 5/6 6/6 2/7 7/7 0/8 2/8 3/8 8/8 8/8 1/9 1/9 2/9 4/9 6/9 8/9 9/9 0/10 1/10

12/12 12/12 12/12 1/13 2/13 2/13 2/13 2/13 3/13 3/13 3/13 5/13 7/13 8/13 8/13 11/13 13/13 13/13 1/14 3/14 14/14 14/14 0/15 2/15

Dose: 75 y5/n5

Dose: 90 y6/n6

5/5 5/5 1/8 5/8 8/8 0/9 3/9 5/9 7/9 9/9 9/9 9/9 4/10 8/10 9/10 10/10 0/11 2/11 2/11 8/11 11/11 11/11 11/11 11/11

0/3 3/3 4/4 6/6 7/7 4/8 10/10 10/10 10/10 10/10 10/11 10/11 11/11 11/11 12/12 12/12 12/12 12/12 12/12 12/12 10/13 13/13 13/13 14/14

K.K. Saha / Journal of Statistical Planning and Inference 141 (2011) 2933–2945

2943

Table 4 (continued ) Dose: 0 y1/n1

Dose: 30 y2/n2

0/11 0/11 0/11 0/11 1/11 1/11 1/11 1/11 1/11 1/11 1/11 2/11 2/11 2/11 2/11 2/11 0/12 0/12 0/12 0/12 0/12 0/12 0/12 1/12 1/12 1/12

1/10 1/10 1/10 1/10 1/10 1/10 1/10 1/10 2/10 0/11 0/11 0/11 0/11 0/11 0/11 0/11 0/11 1/11 1/11 1/11 1/11 2/11 3/11 0/12 0/12 1/12

Dose: 45 y3/n3 0/14 1/14 1/14 2/14 2/14 8/14 0/15 1/15 1/15 2/15 3/15 5/15 15/18

0/11 0/11 0/11 0/11 0/11 0/11 0/11 1/11 1/11 1/11 1/11 2/11 2/11 2/11 2/11 2/11 3/11 5/11 10/11 0/12 0/12 0/12 1/12 1/12 1/12 2/12

Dose: 60 y4/n4 13/13 13/13 1/14 2/14 2/14 3/14 4/14 4/14 5/14 6/14 6/14 6/14 13/14 14/14 0/15 1/15 2/15 3/15 4/15 5/15 6/15 3/16 6/16 18/21

1/10 2/10 10/10 10/10 10/10 10/10 10/10 1/11 1/11 1/11 1/11 4/11 4/11 9/11 9/11 10/11 11/11 11/11 0/12 0/12 0/12 1/12 1/12 5/12 5/12 9/12

7/15 5/16

Dose: 75 y5/n5

Dose: 90 y6/n6

2/12 2/12 10/12 11/12 11/12 12/12 12/12 12/12 10/13 12/13 12/13 12/13 12/13 13/13 13/13 13/13 7/14 14/14 14/14 15/15

14/14

Table 5 The maximum likelihood estimates of the parameters under the null hypothesis of the equality of the dispersion parameters among treatment groups and the score test statistic in Section 4 for the data in Table 4. Group

p^ 1

p^ 2

p^ 3

p^ 4

p^ 5

p^ 6

f^

ST f (p-value)

All 12 13 14 15 16 23 24 25 26 34 35 36 45 46 56

0.1378 0.0813 0.1132 0.1522 0.1219 0.0773

0.1780 0.1211

0.3159

0.5311

0.7951

0.9439

0.3335 0.0666 0.2249 0.3922 0.2644 0.0417 0.2358 0.3750 0.2711 0.1345 0.4189 0.3697 0.3218 0.5206 0.5244 0.5074

58.2702 10.4566 18.8850 25.9264 31.0256 20.5972 11.2677 28.2822 25.1480 4.7666 10.9418 5.6547 0.4843 0.0374 0.0074 0.0008

0.2935 0.5366 0.8066 0.9642 0.1539 0.1889 0.1623 0.1319

0.2956 0.5351 0.8056 0.9609 0.3356 0.3241 0.3133

0.5388 0.7884 0.9451 0.5456 0.5458

0.7564 0.7594

0.9228 0.9249

(0.0000) (0.0012) (0.0000) (0.0000) (0.0000) (0.0000) (0.0008) (0.0000) (0.0000) (0.0290) (0.0009) (0.0174) (0.4865) (0.8466) (0.9315) (0.9771)

7. Conclusion In this article, we have investigated five test methods for testing the homogeneity of the proportions of several treatment groups of clustered binary data when dispersions among treatment groups are different. In this case, we have derived the test statistic based on adjusted proportions using a sandwich estimator of the variance of the proportion estimators employed by the generalized estimating equations of Zeger and Liang (1986). In addition, we have extended the test procedures based on the beta-binomial likelihood, the quasi-likelihood, and the extended quasi-likelihood studied by Paul and Islam (1994, 1995) for the unequal dispersions case. Moreover, we have included the test procedure studied by Rao and Scott (1992) based on the adjusted proportions using estimated variance of a ratio estimator. We have then compared all five test procedures, through simulations, in terms of Type I error rate and power. The test statistics based on the generalized estimating equations and a ratio estimator hold level quite well in all data situations investigated here. The test statistics based on the beta-binomial likelihood, the quasi-likelihood, and the extended quasi-likelihood generally hold level well for the larger number of cluster sizes, whereas for the small number of cluster sizes, they show somewhat conservative behavior and become extremely conservative particularly for small

2944

K.K. Saha / Journal of Statistical Planning and Inference 141 (2011) 2933–2945

Table 6 The test statistics TGEE and TR for testing the equality of the proportions in the presence of unequal dispersions for the data in Table 4. Group

p^ 0

TGEE (p-value)

TR (p-value)

All 12 13 14 15 16 23 24 25 26 34 35 36 45 46 56

0.3449 0.1043 0.2059 0.2792 0.3357 0.2850 0.2223 0.2924 0.3459 0.3035 0.3768 0.4415 0.4203 0.5916 0.5962 0.8342

552.5047 10.6165 102.4337 202.4789 323.8141 379.0254 22.2307 58.8862 120.6432 218.5859 11.6907 51.8886 116.8276 16.1119 53.3921 11.1644

540.7050 10.4756 101.0735 199.7191 318.7947 368.7199 21.9854 58.1729 118.5739 211.0215 11.5536 50.9393 112.5543 15.8024 51.4301 10.7547

(3.68E  117) (1.12E  03) (4.46E  24) (6.01E  46) (2.14E  72) (2.03E  84) (2.42E  06) (1.67E  14) (4.57E  28) (1.84E  49) (6.28E  04) (5.87E  13) (3.13E  27) (5.97E  05) (2.73E  13) (8.34E  04)

(1.30E  114) (1.21E 03) (8.86E  24) (2.41E 45) (2.65E 71) (3.56E 82) (2.75E 06) (2.40E 14) (1.30E 27) (8.22E  48) (6.76E 04) (9.53E 13) (2.70E 26) (7.03E  05) (7.42E  13) (1.04E 03)

proportions. For small number of cluster sizes, these three test statistics are less powerful than those of the test procedures based on the adjusted proportions. However, for large number of cluster sizes, all five test procedures have the same power property for all data situations studied here. Overall, the test statistics based on adjusted proportions perform quite well, holds the nominal level well, are powerful especially for small number of cluster sizes, and they do not require iterative estimates of the dispersion parameters. Although these two tests have very simple forms and hold nominal level quite well, the test based on the generalized estimating equations has some edge in power compared to the test based on a ratio estimate, particularly for small number of cluster sizes. Therefore, the test statistic based on generalized estimating equations might be the method of choice.

Acknowledgments This research was partially supported by the AAUP University Research Grants and the Mathematics Department Kenneth G. Fuller Endowment Funds. This paper was presented to the 2010 International Indian Statistical Association Conference, held in Visakhapatnam, India and the 2010 Statistical Society of Canada Annual Meeting, held in Quebec City, Canada. The referees made thoughtful comments and suggestions which improved the presentation of this manuscript. References Bowman, D., Chen, J.J., George, E.O., 1995. Estimating variance functions in developmental toxicity studies. Biometrics 51, 1523–1528. Bowman, D., George, E.O., 1995. A saturated model for analyzing exchangeable binary data: application to clinical and developmental toxicity studies. Journal of the American Statistical Association 90, 871–879. Brass, W., 1958. Models of birth distribution in human populations. Bulletin of the International Statistical Institute 36, 165–167. Breslow, N.E., 1984. Extra-Poisson variation in log-linear models. Applied Statistics 33, 38–44. Breslow, N.E., 1990. Tests of hypotheses in overdispersed Poisson regression and other quasi-likelihood models. Journal of the American Statistical Association 85, 565–571. Cochran, W.G., 1954. Some methods for strengthening the common w2 tests. Biometrics 10, 417–451. Crowder, M.J., 1978. Beta-binomial ANOVA for proportions. Applied Statistics 27, 34–37. Donner, A., Banting, D., 1989. Adjustment of frequently used chi-square procedures for the effect of site-to-site dependencies in the analysis of dental data. Journal of Dental Research 68, 1350–1354. Donovan, A.M., Ridout, M.S., James, D.J., 1994. Assessment of somaclonal variation in apple. II. Rooting ability and shoot proliferation. Journal of Horticultural Science 69, 115–122. Gibson, G.J., Austin, E.J., 1996. Fitting and testing spatio-temporal stochastic models with applications in plant pathology. Plant Pathology 45, 172–184. Kleinman, J.C., 1973. Proportions with extraneous variance: single and independent samples. Journal of the American Statistical Association 8, 46–54. Keen, K.J., 1996. Limiting the effects of single-member families in the estimation of the intraclass correlation. Biometrics 52, 823–832. Kuk, A.Y.C., 2004. A litter-based approach to risk assessment in developmental toxicity studies via a power family of completely monotone functions. Applied Statistics 53, 369–386. Kupper, L.L., Portier, C., Hogan, M.D., Yamamoto, E., 1986. The impact of litter effects on dose-response modeling in teratology. Biometrics 42, 85–98. Lee, Y., 2004. Estimating intraclass correlation for binary data using extended quasi-likelihood. Statistical Modelling 4, 113–126. Lindsey, J.K., Altham, P.M.E., 1998. Analysis of the human sex ratio by using overdispersion models. Applied Statistics 47, 149–157. Moore, D.F., Tsiatis, A., 1991. Robust estimation of the standard error in moment methods for extra-binomial and extra-Poisson variation. Biometrics 47, 383–401. Neyman, J., 1959. Optimal asymptotic tests for composite hypotheses. In: Grenander, U. (Ed.), Probability and Statistics. Wiley, New York, pp. 213–234. Ochi, Y., Prentice, R.L., 1984. Likelihood inference in a correlated probit regression model. Biometrika 71, 531–543. Otake, M., Prentice, R.L., 1984. The analysis of chromosomally aberrant cells based on beta-binomial distribution. Radiation Research 98, 456–470. Paul, S.R., 1982. Analysis of proportions of affected foetuses in teratological experiments. Biometrics 38, 361–370. Paul, S.R., Islam, A.S., 1994. CðaÞ tests for homogeneity of proportions in toxicology in presence of beta-binomial over-dispersion. In: Williams, D. (Ed.), Statistics in Toxicology. Oxford University Press, New York, pp. 66–74. Paul, S.R., Islam, A.S., 1995. Analysis of proportions in the presence of over-/under-dispersion. Biometrics 51, 1400–1411.

K.K. Saha / Journal of Statistical Planning and Inference 141 (2011) 2933–2945

2945

Paul, S.R., Islam, A.S., 1998. Joint estimation of the mean and dispersion parameters in the analysis of proportions: a comparison of efficiency and bias. The Canadian Journal of Statistics 26, 83–94. Paul, S.R., Saha, K.K., Balasooriya, U., 2003. An empirical investigation of different operating characteristics of several estimators of the intraclass correlation in the analysis of binary data. Journal of Statistical Computation and Simulation 73, 507–523. Paul, S.R., Saha, K.K., 2007. The generalized linear model and extensions: a review and some biological and environmental applications. Environmetrics 18, 421–443. Prentice, R.L., 1986. Binary regression using an extended beta-binomial distribution, with discussion of correlation induced by covariate measurement errors. Journal of the American Statistical Association 81, 321–327. Rao, J.N.K., Scott, A.J., 1992. A simple method for the analysis of clustered binary data. Biometrics 48, 577–585. Ridout, M.S., Demetrio, C.G.B., Firth, D., 1999. Estimating intraclass correlation for binary data. Biometrics 55, 137–148. Weil, C.S., 1970. Selection of valid number of sampling units and a consideration of their combination in toxicological studies involving reproduction, teratogenesis or carcinogenesis reproduction, teratogenesis. Food and Cosmetic Toxicology 8, 177–182. Williams, D.A., 1975. Analysis of binary responses from toxicological experiments involving reproduction and teratogenicity. Biometrics 31, 949–952. Xiao, Y., Liu, J., Bhandary, M., 2010. Profile likelihood based confidence intervals for common intraclass correlation coefficient. Communications in Statistics—Simulation and Computation 39, 111–118. Yu, C., Zelterman, D., 2008. Sums of exchangeable Bernoulli random variables for family and litter frequency data. Computational Statistics & Data Analysis 52, 1636–1649. Zeger, S.L., Liang, K.Y., 1986. Longitudinal data analysis for discrete and continuous outcomes. Biometrics 42, 121–130. Zhu, J., Eickhoff, J.C., Kaiser, M.S., 2003. Modeling the dependence between number of trials and success probability in beta-binomial-Poisson mixture distributions. Biometrics 59, 955–961. Zou, G., Donner, A., 2004. Confidence interval estimation of the intraclass correlation coefficient for binary outcome data. Biometrics 60, 807–811.