# Homogeneity Test of Many-to-One Risk Differences for Correlated Binary Data under Optimal Algorithms

Authors: Keyi Mou; Zhiming Li
Journal: Complexity
 (2021)

---

## Abstract

In clinical studies, it is important to investigate the effectiveness of different therapeutic designs, especially, multiple treatment groups to one control group. The paper mainly studies homogeneity test of many-to-one risk differences from correlated binary data under optimal algorithms. Under Donner’s model, several algorithms are compared in order to obtain global and constrained MLEs in terms of accuracy and efficiency. Further, likelihood ratio, score, and Wald-type statistics are proposed to test whether many-to-one risk differences are equal based on optimal algorithms. Monte Carlo simulations show the performance of these algorithms through the total averaged estimation error, SD, MSE, and convergence rate. Score statistic is more robust and has satisfactory power. Two real examples are given to illustrate our proposed methods.


---

## Body

## 1. Introduction

Binary data are often encountered for paired organs (e.g., eyes and ears) in medical clinical studies. The responses of each patient are collected and recorded as paired data at the end of the study. The outcome can be none, unilateral, or bilateral cured. Data from all patients can be summarized in a contingency table. The correlation between responses from paired parts should be taken into account to avoid biased or misleading results.Some probability models have been proposed for analyzing correlated paired data. Rosner introduced a constantR model under the assumption that the conditional probability of a response at one side of the paired body parts given response at the other side was R times the unconditional probability [1]. Under Rosner’s model, asymptotic and exact tests were discussed [2–5]. However, Dallal pointed out that Rosner’s model could give a poor fit if the characteristic was almost certain to occur bilaterally with widely varying group-specific prevalence [6]. He assumed that each group had a constant conditional probability γ and derived the likelihood ratio statistic to test prevalence equality. M’lan and Chen presented three objective Bayesian methods for bilateral data under Dallal’s model [7]. Donner proposed an alternative model by assuming that the correlation coefficient was a fixed constant ρ in each of the groups [8]. Thompson proved that Donner’s model could make full use of single and two-organ data to optimize the power of study [9]. Pei et al. applied the model into stratified paired data and assumed that the correlation coefficients of responses were the same to all subjects in two groups of each stratum [10].Testing the homogeneity has received considerable attention for bilateral binary data. In ophthalmologic studies, Rosner [1] proposed two statistics to test the equality of rate difference in (in)dependence models. Tang et al. [2] developed exact and approximate unconditional procedures for the aforementioned statistics in small sample designs or sparse data structures. Further, Tang et al. [11] developed several statistics for testing the equality of cure rates, including likelihood ratio, score, and two Wald-type statistics in the (in)dependence models. Ma et al. [3] extended these tests to multigroup cases and investigated whether the response rates of the g groups g≥2 were identical under Rosner’s model. From the above results, we note that it is crucial to derive the global and constrained MLEs under the hypotheses. However, there are usually no closed-form solutions for maximum likelihood estimates (MLEs). Under Donner’s model, Ma and Liu [12] used two-step algorithm to obtain MLEs and developed several tests for the proportion equality among g groups g≥2. Liu et al. [13] also used the method for constrained MLEs. Peng et al. [14] constructed confidence intervals (CIs) of proportion ratio under Rosner’s model. They introduced Fisher scoring algorithm for constrained MLEs. Many algorithms were proposed to obtain MLEs for correlated binary data. However, there are few research studies on comparison of different algorithms for MLEs in multigroup binary design.Under Donner’s model, this paper aims to provide several algorithms for calculating global and constrained MLEs and extends the homogeneity tests of Tang et al. [11] to many-to-one case under optimal algorithms. Fisher scoring algorithm, two-step method, and generalized expectation-maximization (GEM) algorithm are taken into account, since they are widely used in calculating MLEs. Optimal algorithms for MLEs required by the objective test can be found through comparing these algorithms. The rest of this article is organized as follows. In Section 2, we review data structure and establish Donner’s model for multigroup correlated binary data. Global and constrained MLEs are derived by various algorithms in Section 3. Based on the optimal algorithms, the likelihood ratio, score, and Wald-type statistics are constructed for testing the equality of many-to-one risk differences. The performance of algorithms is compared by the total averaged estimation error, SD of the averaged estimation error, MSE, and convergence rate in Section 4. Monte Carlo simulations show the empirical type I error rate and power of these tests. Two real examples are provided to illustrate the proposed methods in Section 5. Conclusions and further work are given in Section 6.


## 2. Preliminaries

Suppose there areg groups involving M individuals in the clinical trial, where the first group is control group and other g−1 groups are treatment groups. Let mli be the number of patients with l responses l=0,1,2 in the i-th group i=1,2,…,g and Ni=∑l=02mli be the total number of patients in the i-th group, which is assumed to be fixed. The data structure is shown in Table 1.Table 1
Data structure of multigroup binary data.
Number of responseslGroupiTotal12…g0m01p01m02p02…m0gp0gS01m11p11m12p12…m1gp1gS12m21p21m22p22…m2gp2gS2TotalN1N2…NgMLetplil=0,1,2,i=1,2,…,g be the probabilities of none, unilateral, and bilateral response(s) in the i-th group, where p0i+p1i+p2i=1 for any fixed i. Denote mi=m0i,m1i,m2iT and pi=p0i,p1i,p2iT. For the i-th group, mi follows a trinomial distribution. Thus, the probability density of mi is expressed as follows:(1)fpi|mi=Ni!m0i!m1i!m2i!p0im0ip1im1ip2im2i,i=1,2,…,g.LetZhik be an indicator of the k-th organ’s response k=1,2 for the h-th patient h=1,…,Ni in the i-th group. If there is a response, then Zhik=1, and 0 otherwise. Suppose that PrZhik=1=πi0≤πi≤1, and CorrZhik,Zhi3−k=ρ0≤ρ≤1 under Donner’s model. Thus, the probabilities pli can be obtained by(2)p0i=ρ1−πi+1−ρ1−πi2,p1i=2πi1−ρ1−πi,p2i=ρπi+1−ρπi2,for i=1,2,…,g. Based on the observed data m=m1,m2,…,mg, the log-likelihood function can be given by(3)ℓπ,ρ|m=∑i=1gliπi,ρ|mi+logC=∑i=1gm0ilogρ1−πi+1−ρ1−πi2+m1ilog2πi1−ρ1−πi+m2ilogρπi+1−ρπi2+logC,where C=∏i=1gNi!/m0i!m1i!m2i! is a constant.LetδT=δ2,…,δg, where δi=πi−π1i=2,…,g is the risk difference between the first group and the i-th group. We are interested to test the hypotheses below.(4)H0:δ2=⋯=δg≜δversusHa:δiis not all the same.UnderH0, that is πi=π1+δi=2,…,g, the log-likelihood function can be rewritten as(5)ℓ0π1,δ,ρ|m=l01π1,ρ|m1+∑i=2gl0iπ1,δ,ρ|mi+logC,where(6)l01π1,ρ|m1=m01logρ1−π1+1−ρ1−π12+m11log2π11−ρ1−π1+m21logρπ1+1−ρπ12,l0iπ1,δ,ρ|mi=m0ilogρ1−π1−δ+1−ρ1−π1−δ2+m1ilog2π1+δ1−ρ1−π1−δ+m2ilogρπ1+δ+1−ρπ1+δ2,i=2,…,g,and C=∏i=1gNi!/m0i!m1i!m2i! is a constant.


## 3. Test Methods

In this section, the global and constrained MLEs are first derived by various algorithms. Then, likelihood ratio, score, and Wald-type tests are constructed based on the optimal algorithms.


### 3.1. Global MLEs

Letπ^i and ρ^ be the global MLEs of πi and ρ. For the unknown parameters πii=1,2,…,g and ρ, their global MLEs are the solutions of the following equations:(7)∂ℓ∂πi=0,∂ℓ∂ρ=0,i=1,2,…,gwhere(8)∂ℓ∂πi=2πi−1m1iπiπi−1+m2iρ+2πi−2ρπiπiρ+πi−ρπi−m0iρ+2πi−2ρπi−2πi−1ρπi−πi+1,∂ℓ∂ρ=∑i=1gm1iρ−1−πi−1m2iρ+πi−ρπi+πim0iρπi−πi+1.However, there are no closed-form solutions for the above equations. Thus, we need to obtain the global MLEsπ^ii=1,2,…,g and ρ^ by different algorithms.


#### 3.1.1. Global MLEs Based on Fisher Scoring Algorithm

The initial values ofπi and ρ can be given by(9)πi0=m1i+2m2i2Ni,i=1,2,…,g,ρ0=4S0S2−S12S1+2S0S1+2S2.Thet+1-th approximates πt+1=π1t+1,π2t+1,…,πgt+1T and ρt+1 can be obtained by Fisher scoring algorithm:(10)πt+1ρt+1=πtρt+I1−1π,ρ×∂ℓ∂π∂ℓ∂ρ|π=πt,ρ=ρt,where ∂ℓ/∂π=∂l1/∂π1,∂l2/∂π2,…,∂lg/∂πgT and I1 is a g+1×g+1 Fisher information matrix (see Appendix A). Repeat the process until the result converges.


#### 3.1.2. Global MLEs Based on Two-Step Method

The two-step method is described by a third-order polynomial and Newton–Raphson algorithm. The detailed procedure is provided below.(i)
Take the initial valueρ0=4S0S2−S12/S1+2S0S1+2S2. Moreover, equation (8) can be simplified as a third-order polynomial:(11)4ρ−2ρ2−2Niπi3+3ρ2Ni−ρ5m0i+6m1i+7m2i+2m0i+3m1i+4m2iπi2+4ρ−ρ2Ni−2ρm0i−m1i−2m2iπi−ρm1i+m2i=0,i=1,2,…,g.Putρtt=0,1,2,… into the polynomial and solve its real root to obtain the t-th approximates of πii=1,2,…,g, denoted by πit.(ii)
Update thet+1-th approximate of ρ by Newton–Raphson algorithm:(12)ρt+1=ρt−∂2ℓ∂ρ2πt,ρt−1∂ℓ∂ρπt,ρtwhere(13)∂2ℓ∂ρ2=−∑i=1gπi2m0iρπi−πi+12+m1iρ−12+πi−12m2iρ+πi−ρπi2,andπt=π1t,π2t,…,πgtT. Repeat (i)-(ii) until convergence.


#### 3.1.3. Global MLEs Based on GEM Algorithm

According to equation (1), we have(14)fπi,ρ|mi∝ρ1−πi+1−ρ1−πi2m0i2πi1−ρ1−πim1iρπi+1−ρπi2m2i.Suppose patients with 0 response can be divided into two parts, whose responsibilities areρ1−πi and 1−ρ1−πi2, respectively. Let latent variables z1i and m0i−z1i be their total numbers. Observable variable m2i can also be split into two latent variables z2i and m2i−z2i. Suppose the probability of result happening in z2i (or m2i−z2i) individuals is ρπi (or 1−ρπi2). When mi,ρ, and πi are given, z1i and z2i follow binomial distributions:(15)z1i∼Bm0i,ρ1−πiρ1−πi+1−ρ1−πi2,z2i∼Bm2i,ρπiρπi+1−ρπi2.Denotez=z1,z2,…,zg and zi=z1i,z2iT. Then, m,z are complete data and observable data are incomplete data. Based on complete data, we have(16)fπi,ρ|mi,zi∝ρ1−πiz1i1−ρ1−πi2m0i−z1i2πi1−ρ1−πim1i×ρπiz2i1−ρπi2m2i−z2i.Thus, the log-likelihood function about complete data is(17)ℓπ,ρ|m,z∝∑i=1gliπi,ρ|mi,zi,where(18)liπi,ρ|mi,zi=z1ilogρ1−πi+m0i−z1ilog1−ρ1−πi2+z2ilogρπi+m1ilog2πi1−ρ1−πi+m2i−z2ilog1−ρπi2.The initial values ofρ,πi are defined in equation (9). The process of GEM algorithm is described by the expectation (E) and maximization (M) steps as follows.(i)
E Step. Given m and the current approximates πt,ρt. Denote(19)Ez1iz1i|m0i,πit,ρt=m0iρt1−πitρt1−πit+1−ρt1−πit2,Ez2iz2i|m2i,πit,ρt=m2iρtπitρtπit+1−ρtπit2.LetQπ,ρ|∣m,πt,ρt be the expected value of the log-likelihood function of π,ρ, with respect to the current conditional distribution of z as follows:(20)Qπ,ρ|m,πt,ρt=Ezℓπ,ρ|m,z∝∑i=1gEziliπi,ρ|mi,ziwhere(21)Eziliπi,ρ|mi,zi=Ez1iz1i|m0i,πit,ρtlogρ1−πi+m1ilog2πi1−ρ1−πi+m0i−Ez1iz1i|m0i,πit,ρtlog1−ρ1−πi2+Ez2iz2i|m2i,πit,ρtlogρπi+m2i−Ez2iz2i|m2i,πit,ρtlog1−ρπi2.(ii)
M Step. Update π1,π2,…,πg and ρ in order to successively maximize expected value of the log-likelihood function in E step. The new approximate of parameters can be obtained by maximizing Qπ,ρ|m,πt,ρt when other parameters are given as their latest approximates. Repeat E and M steps until the result converges.


### 3.2. Constrained MLEs

Letπ˜1,δ˜ and ρ˜ be the constrained MLEs of π1,δ, and ρ under H0. Under H0: δ2=⋯=δg≜δ, it is obvious that πi=π1+δi=2,…,g. Thus, the constrained MLEs satisfy the following equations:(22)∂ℓ0∂π1=2π1−1m11π1π1−1+m21ρ+2π1−2ρπ1π1ρ+π1−ρπ1−m01ρ+2π1−2ρπ1−2π1−1ρπ1−π1+1+∑i=2gm2i21−ρδ+π1+ρδ+π12−ρδ+π1δ+π1−1−m0i2ρ−1δ+π1+2−ρδ+π1−1δ+π11−ρ−1+m1i2δ+2π1−1δ+π1δ+π1−1=0,∂ℓ0∂δ=∑i=2gm2i21−ρδ+π1+ρδ+π12−ρδ+π1δ+π1−1−m0iρ−12δ+2π1−1+1δ+π1−1δ+π11−ρ−1+m1i2δ+2π1−1δ+π1δ+π1−1=0,∂ℓ0∂ρ=m11ρ−1−π1−1m21ρ+π1−ρπ1+π1m01ρπ1−π1+1+∑i=2gm11ρ−1−π1+δi−1m21ρ+1−ρπ1+δi+π1+δim01ρ−1π1+δi+1=0.However, their closed-form solutions are not given by these equations. Thus, we introduce Fisher scoring algorithm, two-stage procedure, and GEM algorithm to obtain the constrained MLEsδ˜, π˜1, and ρ˜.


#### 3.2.1. Constrained MLEs Based on Fisher Scoring Algorithm

The initial values ofρ,π1 are defined in equation (9), and δ0=1. The t+1-th updates of δ, π1, and ρ can be calculated by Fisher scoring algorithm as follows:(23)δt+1π1t+1ρt+1=δtπ1tρt+I2−1δ,π1,ρ×∂ℓ0∂δ∂ℓ0∂π1∂ℓ0∂ρ|δ=δt,π1=π1t,ρ=ρt,where I2 is a 3×3 Fisher information matrix (see Appendix B).


#### 3.2.2. Constrained MLEs Based on Two-Stage Procedure

The two-stage procedure is different from the two-step method in Section3.1.2. Firstly, the MLE δ˜ of δ is given by Newton–Raphson algorithm. Then, π˜1 and ρ˜ are obtained by Fisher scoring algorithm under given MLE δ˜. The detailed process is described as follows.(i)
The initial valuesπ10 and ρ0 are defined in the equation (9), and δ0=1. The t+1-th approximate δt+1 of δ is obtained by Newton–Raphson algorithm:(24)δt+1=δt−∂2ℓ0∂δ2−1∂ℓ0∂δ|δ=δt,π1=π1t,ρ=ρt,where(25)∂2ℓ0∂δ2=∑i=2gm0iρ−12δ+2π1−1+1δ+π1−12δ+π11−ρ−1−2ρ−1δ+π1−1δ+π11−ρ−1−ρ−1ρ−12δ+2π1−1+1δ+π1−1δ+π11−ρ−12+m2i−2δ+π11−ρ+ρ2δ+π121−ρ+ρδ+π12−2ρ−1δ+π121−ρ+ρδ+π1+m1i2δ+π1δ+π1−1−2δ+2π1−1δ+π1δ+π1−12−2δ+2π1−1δ+π12δ+π1−1.(ii)
Givenδ=δ˜, the t+1-th approximates of π1 and ρ can be calculated as(26)π1t+1ρt+1=π1tρt+I3−1π1,ρ×∂ℓ0∂π1∂ℓ0∂ρ|δ=δ˜,π1=π1t,ρ=ρt,whereI3 is a 2×2 Fisher information matrix of π1 and ρ (see Appendix B).


#### 3.2.3. Constrained MLEs Based on GEM Algorithm

Similar to global MLEs, the GEM algorithm is used to calculate constrained MLEs underH0. The initial values of ρ,π1 are defined in equation (9), and δ0=∑i=2gπi0−π10/g−1. The detailed process is shown as follows.(i)
E Step. Let(27)Qπ1,δ,ρ|m,π1t,δt,ρt=Ezℓ0π1,δ,ρ|m,z=Ez1l01π1,ρ|m1,z1+∑i=2gEzil0iπ1,δ,ρ|mi,zi,where(28)Ez1l01π1,ρ|m1,z1=Ez11z11|m01,π1t,ρtlogρ1−π1+m01−Ez11z11|m01,π1t,ρtlog1−ρ1−π12+m11log2π11−ρ1−π1+Ez21z21|m21,π1t,ρtlogρπ1+m21−Ez21z21|m21,π1t,ρtlog1−ρπ12,Ezil0iπ1,δ,ρ|mi,zi=Ez1iz1i|m0i,π1t,δt,ρtlogρ1−π1−δ+m0i−Ez1iz1i|m0i,π1t,δt,ρtlog1−ρ1−π1−δ2+m1ilog2π1+δ1−ρ1−π1−δ+m2i−Ez2iz2i|m2i,π1t,δt,ρtlog1−ρπ1+δ2+Ez2iz2i|m2i,π1t,δt,ρtlogρπ1+δ,i=2,…,g,Ez11z11|m01,π1t,ρt=m01ρt1−π1tρt1−π1t+1−ρt1−π1t2,Ez21z21|m21,π1t,ρt=m21ρtπ1tρtπ1t+1−ρtπ1t2,Ez1iz1i|m0i,π1t,δt,ρt=m0iρt1−π1t−δtρt1−π1t−δt+1−ρt1−π1t−δt2,Ez2iz2i|m2i,π1t,δt,ρt=m2iρtπ1t+δtρtπ1t+δt+1−ρtπ1t+δt2,i=2,…,g.(ii)
M Step. Similar to Section 3.1.3, we choose parameters π1,δ,ρ so that Qπ1,δ,ρ|m,π1t,δt,ρt increases. Repeat E and M steps until the result converges.


### 3.3. Likelihood Ratio Test

Likelihood ratio test statistic can be constructed through the global and constrained MLEs as follows:(29)TL=2ℓπ^,ρ^−ℓ0π˜1,δ˜,ρ˜=2∑i=1gliπ^i,ρ^−l01π˜1,ρ˜−∑i=2gl0iπ˜1,δ˜,ρ˜where π^i,ρ^i=1,2,…,g are the global MLEs and π˜1,δ˜,ρ˜ are the constrained MLEs.


### 3.4. Score Test

Note thatH0:δ2=⋯=δg≜δ is equivalent to π2=⋯=πg. The homogeneity test of many-to-one risk differences can be achieved by testing the equality of risks in treatment groups. Denote U=∂l2/∂π2,…,∂lg/∂πg,0,0. Score test statistic is derived as(30)TSC=UI4π2,…,πg,π1,ρ−1UT|π=π˜,ρ=ρ˜,where π=π1,π2,…,πgT and π˜=π˜1,π˜2,…,π˜gT. π˜i and ρ˜ are constrained MLEs under H0. I4 is a g+1×g+1 information matrix (see Appendix C for more information). Thus, TSC can be simplified as(31)TSC=∑i=2gui2π˜i,ρ˜aiπ˜i,ρ˜+∑i=2guiπ˜i,ρ˜biπ˜i,ρ˜/aiπ˜i,ρ˜2Dπ˜,ρ˜−∑i=1gbi2π˜i,ρ˜/aiπ˜i,ρ˜,where(32)uiπi,ρ=∂li∂πi=2πi−1m1iπiπi−1+m2iρ+2πi−2ρπiπiρ+πi−ρπi−m0iρ+2πi−2ρπi−2πi−1ρπi−πi+1,aiπi,ρ=E−∂2li∂πi2=Ni−4ρ2πi2+4ρ2πi−ρ2+6ρπi2−6ρπi+2ρ−2πi2+2πiπi1−πiρ+πi−ρπiρπi−πi+1,biπi,ρ=E−∂2li∂πi∂ρ=Niρ2πi−1ρ+πi−ρπiρπi−πi+1,Dπ,ρ=E−∂2ℓ∂ρ2=∑i=1gNiπiρ+11−πi1−ρρ+πi−ρπiρπi−πi+1.


### 3.5. Wald-Type Test

Letβ=π2,…,πg,π1,ρ and(33)C=1−10……00001−1000⋮⋱⋱⋱⋮⋮⋮⋮⋱⋱⋱⋮⋮⋮0……01−100g−2×g+1.The null hypothesisH0:δ2=⋯=δg≜δ is equivalent to CβT=0. Thus, Wald-type statistic is(34)TW=βCTCI4−1CT−1CβT|π=π^,ρ=ρ^,where the Fisher information matrix I4 is the same as that of score test (see Appendix D). It can be simplified as(35)TW=∑i=1g−2∑j=1g−2π^i+1−π^i+2π^j+1−π^j+2∑k=1g−2rkiπ^,ρ^rkjπ^,ρ^,where(36)rijπ,ρ=−1liiπ,ρ∑k=ji−1likπ,ρrkjπ,ρ,i>j,1lijπ,ρ,i=j,0,i<j,lijπ,ρ=fijπ,ρ−∑k=1j−1likπ,ρljkπ,ρljjπ,ρ,i>j,fijπ,ρ−∑k=1j−1ljkπ,ρ21/2,i=j,0,i<j,fijπ,ρ=Iijπ,ρ−Ii+1jπ,ρ−Iij+1π,ρ−Ii+1j+1π,ρ,Iijπ,ρ=1ai+1πi+1,ρ+bi+1πi+1,ρbj+1πj+1,ρ/ai+1πi+1,ρaj+1πj+1,ρDπ,ρ−∑i=1gbi2πi,ρ/aiπi,ρ,i=j,bi+1πi+1,ρbj+1πj+1,ρ/ai+1πi+1,ρaj+1πj+1,ρDπ,ρ−∑i=1gbi2πi,ρ/aiπi,ρ,i≠j.UnderH0, test statistics TL,TSC, and TW are asymptotically distributed as chi-square distribution with g−2 degrees of freedom. Thus, H0 should be rejected if the value of test statistic is larger than χg−2,1−α2 at the significant level α, where χg−2,1−α2 is the 1001−α percentile of the chi-square distribution with g−2 degrees of freedom.


## 3.1. Global MLEs

Letπ^i and ρ^ be the global MLEs of πi and ρ. For the unknown parameters πii=1,2,…,g and ρ, their global MLEs are the solutions of the following equations:(7)∂ℓ∂πi=0,∂ℓ∂ρ=0,i=1,2,…,gwhere(8)∂ℓ∂πi=2πi−1m1iπiπi−1+m2iρ+2πi−2ρπiπiρ+πi−ρπi−m0iρ+2πi−2ρπi−2πi−1ρπi−πi+1,∂ℓ∂ρ=∑i=1gm1iρ−1−πi−1m2iρ+πi−ρπi+πim0iρπi−πi+1.However, there are no closed-form solutions for the above equations. Thus, we need to obtain the global MLEsπ^ii=1,2,…,g and ρ^ by different algorithms.


### 3.1.1. Global MLEs Based on Fisher Scoring Algorithm

The initial values ofπi and ρ can be given by(9)πi0=m1i+2m2i2Ni,i=1,2,…,g,ρ0=4S0S2−S12S1+2S0S1+2S2.Thet+1-th approximates πt+1=π1t+1,π2t+1,…,πgt+1T and ρt+1 can be obtained by Fisher scoring algorithm:(10)πt+1ρt+1=πtρt+I1−1π,ρ×∂ℓ∂π∂ℓ∂ρ|π=πt,ρ=ρt,where ∂ℓ/∂π=∂l1/∂π1,∂l2/∂π2,…,∂lg/∂πgT and I1 is a g+1×g+1 Fisher information matrix (see Appendix A). Repeat the process until the result converges.


### 3.1.2. Global MLEs Based on Two-Step Method

The two-step method is described by a third-order polynomial and Newton–Raphson algorithm. The detailed procedure is provided below.(i)
Take the initial valueρ0=4S0S2−S12/S1+2S0S1+2S2. Moreover, equation (8) can be simplified as a third-order polynomial:(11)4ρ−2ρ2−2Niπi3+3ρ2Ni−ρ5m0i+6m1i+7m2i+2m0i+3m1i+4m2iπi2+4ρ−ρ2Ni−2ρm0i−m1i−2m2iπi−ρm1i+m2i=0,i=1,2,…,g.Putρtt=0,1,2,… into the polynomial and solve its real root to obtain the t-th approximates of πii=1,2,…,g, denoted by πit.(ii)
Update thet+1-th approximate of ρ by Newton–Raphson algorithm:(12)ρt+1=ρt−∂2ℓ∂ρ2πt,ρt−1∂ℓ∂ρπt,ρtwhere(13)∂2ℓ∂ρ2=−∑i=1gπi2m0iρπi−πi+12+m1iρ−12+πi−12m2iρ+πi−ρπi2,andπt=π1t,π2t,…,πgtT. Repeat (i)-(ii) until convergence.


### 3.1.3. Global MLEs Based on GEM Algorithm

According to equation (1), we have(14)fπi,ρ|mi∝ρ1−πi+1−ρ1−πi2m0i2πi1−ρ1−πim1iρπi+1−ρπi2m2i.Suppose patients with 0 response can be divided into two parts, whose responsibilities areρ1−πi and 1−ρ1−πi2, respectively. Let latent variables z1i and m0i−z1i be their total numbers. Observable variable m2i can also be split into two latent variables z2i and m2i−z2i. Suppose the probability of result happening in z2i (or m2i−z2i) individuals is ρπi (or 1−ρπi2). When mi,ρ, and πi are given, z1i and z2i follow binomial distributions:(15)z1i∼Bm0i,ρ1−πiρ1−πi+1−ρ1−πi2,z2i∼Bm2i,ρπiρπi+1−ρπi2.Denotez=z1,z2,…,zg and zi=z1i,z2iT. Then, m,z are complete data and observable data are incomplete data. Based on complete data, we have(16)fπi,ρ|mi,zi∝ρ1−πiz1i1−ρ1−πi2m0i−z1i2πi1−ρ1−πim1i×ρπiz2i1−ρπi2m2i−z2i.Thus, the log-likelihood function about complete data is(17)ℓπ,ρ|m,z∝∑i=1gliπi,ρ|mi,zi,where(18)liπi,ρ|mi,zi=z1ilogρ1−πi+m0i−z1ilog1−ρ1−πi2+z2ilogρπi+m1ilog2πi1−ρ1−πi+m2i−z2ilog1−ρπi2.The initial values ofρ,πi are defined in equation (9). The process of GEM algorithm is described by the expectation (E) and maximization (M) steps as follows.(i)
E Step. Given m and the current approximates πt,ρt. Denote(19)Ez1iz1i|m0i,πit,ρt=m0iρt1−πitρt1−πit+1−ρt1−πit2,Ez2iz2i|m2i,πit,ρt=m2iρtπitρtπit+1−ρtπit2.LetQπ,ρ|∣m,πt,ρt be the expected value of the log-likelihood function of π,ρ, with respect to the current conditional distribution of z as follows:(20)Qπ,ρ|m,πt,ρt=Ezℓπ,ρ|m,z∝∑i=1gEziliπi,ρ|mi,ziwhere(21)Eziliπi,ρ|mi,zi=Ez1iz1i|m0i,πit,ρtlogρ1−πi+m1ilog2πi1−ρ1−πi+m0i−Ez1iz1i|m0i,πit,ρtlog1−ρ1−πi2+Ez2iz2i|m2i,πit,ρtlogρπi+m2i−Ez2iz2i|m2i,πit,ρtlog1−ρπi2.(ii)
M Step. Update π1,π2,…,πg and ρ in order to successively maximize expected value of the log-likelihood function in E step. The new approximate of parameters can be obtained by maximizing Qπ,ρ|m,πt,ρt when other parameters are given as their latest approximates. Repeat E and M steps until the result converges.


## 3.1.1. Global MLEs Based on Fisher Scoring Algorithm

The initial values ofπi and ρ can be given by(9)πi0=m1i+2m2i2Ni,i=1,2,…,g,ρ0=4S0S2−S12S1+2S0S1+2S2.Thet+1-th approximates πt+1=π1t+1,π2t+1,…,πgt+1T and ρt+1 can be obtained by Fisher scoring algorithm:(10)πt+1ρt+1=πtρt+I1−1π,ρ×∂ℓ∂π∂ℓ∂ρ|π=πt,ρ=ρt,where ∂ℓ/∂π=∂l1/∂π1,∂l2/∂π2,…,∂lg/∂πgT and I1 is a g+1×g+1 Fisher information matrix (see Appendix A). Repeat the process until the result converges.


## 3.1.2. Global MLEs Based on Two-Step Method

The two-step method is described by a third-order polynomial and Newton–Raphson algorithm. The detailed procedure is provided below.(i)
Take the initial valueρ0=4S0S2−S12/S1+2S0S1+2S2. Moreover, equation (8) can be simplified as a third-order polynomial:(11)4ρ−2ρ2−2Niπi3+3ρ2Ni−ρ5m0i+6m1i+7m2i+2m0i+3m1i+4m2iπi2+4ρ−ρ2Ni−2ρm0i−m1i−2m2iπi−ρm1i+m2i=0,i=1,2,…,g.Putρtt=0,1,2,… into the polynomial and solve its real root to obtain the t-th approximates of πii=1,2,…,g, denoted by πit.(ii)
Update thet+1-th approximate of ρ by Newton–Raphson algorithm:(12)ρt+1=ρt−∂2ℓ∂ρ2πt,ρt−1∂ℓ∂ρπt,ρtwhere(13)∂2ℓ∂ρ2=−∑i=1gπi2m0iρπi−πi+12+m1iρ−12+πi−12m2iρ+πi−ρπi2,andπt=π1t,π2t,…,πgtT. Repeat (i)-(ii) until convergence.


## 3.1.3. Global MLEs Based on GEM Algorithm

According to equation (1), we have(14)fπi,ρ|mi∝ρ1−πi+1−ρ1−πi2m0i2πi1−ρ1−πim1iρπi+1−ρπi2m2i.Suppose patients with 0 response can be divided into two parts, whose responsibilities areρ1−πi and 1−ρ1−πi2, respectively. Let latent variables z1i and m0i−z1i be their total numbers. Observable variable m2i can also be split into two latent variables z2i and m2i−z2i. Suppose the probability of result happening in z2i (or m2i−z2i) individuals is ρπi (or 1−ρπi2). When mi,ρ, and πi are given, z1i and z2i follow binomial distributions:(15)z1i∼Bm0i,ρ1−πiρ1−πi+1−ρ1−πi2,z2i∼Bm2i,ρπiρπi+1−ρπi2.Denotez=z1,z2,…,zg and zi=z1i,z2iT. Then, m,z are complete data and observable data are incomplete data. Based on complete data, we have(16)fπi,ρ|mi,zi∝ρ1−πiz1i1−ρ1−πi2m0i−z1i2πi1−ρ1−πim1i×ρπiz2i1−ρπi2m2i−z2i.Thus, the log-likelihood function about complete data is(17)ℓπ,ρ|m,z∝∑i=1gliπi,ρ|mi,zi,where(18)liπi,ρ|mi,zi=z1ilogρ1−πi+m0i−z1ilog1−ρ1−πi2+z2ilogρπi+m1ilog2πi1−ρ1−πi+m2i−z2ilog1−ρπi2.The initial values ofρ,πi are defined in equation (9). The process of GEM algorithm is described by the expectation (E) and maximization (M) steps as follows.(i)
E Step. Given m and the current approximates πt,ρt. Denote(19)Ez1iz1i|m0i,πit,ρt=m0iρt1−πitρt1−πit+1−ρt1−πit2,Ez2iz2i|m2i,πit,ρt=m2iρtπitρtπit+1−ρtπit2.LetQπ,ρ|∣m,πt,ρt be the expected value of the log-likelihood function of π,ρ, with respect to the current conditional distribution of z as follows:(20)Qπ,ρ|m,πt,ρt=Ezℓπ,ρ|m,z∝∑i=1gEziliπi,ρ|mi,ziwhere(21)Eziliπi,ρ|mi,zi=Ez1iz1i|m0i,πit,ρtlogρ1−πi+m1ilog2πi1−ρ1−πi+m0i−Ez1iz1i|m0i,πit,ρtlog1−ρ1−πi2+Ez2iz2i|m2i,πit,ρtlogρπi+m2i−Ez2iz2i|m2i,πit,ρtlog1−ρπi2.(ii)
M Step. Update π1,π2,…,πg and ρ in order to successively maximize expected value of the log-likelihood function in E step. The new approximate of parameters can be obtained by maximizing Qπ,ρ|m,πt,ρt when other parameters are given as their latest approximates. Repeat E and M steps until the result converges.


## 3.2. Constrained MLEs

Letπ˜1,δ˜ and ρ˜ be the constrained MLEs of π1,δ, and ρ under H0. Under H0: δ2=⋯=δg≜δ, it is obvious that πi=π1+δi=2,…,g. Thus, the constrained MLEs satisfy the following equations:(22)∂ℓ0∂π1=2π1−1m11π1π1−1+m21ρ+2π1−2ρπ1π1ρ+π1−ρπ1−m01ρ+2π1−2ρπ1−2π1−1ρπ1−π1+1+∑i=2gm2i21−ρδ+π1+ρδ+π12−ρδ+π1δ+π1−1−m0i2ρ−1δ+π1+2−ρδ+π1−1δ+π11−ρ−1+m1i2δ+2π1−1δ+π1δ+π1−1=0,∂ℓ0∂δ=∑i=2gm2i21−ρδ+π1+ρδ+π12−ρδ+π1δ+π1−1−m0iρ−12δ+2π1−1+1δ+π1−1δ+π11−ρ−1+m1i2δ+2π1−1δ+π1δ+π1−1=0,∂ℓ0∂ρ=m11ρ−1−π1−1m21ρ+π1−ρπ1+π1m01ρπ1−π1+1+∑i=2gm11ρ−1−π1+δi−1m21ρ+1−ρπ1+δi+π1+δim01ρ−1π1+δi+1=0.However, their closed-form solutions are not given by these equations. Thus, we introduce Fisher scoring algorithm, two-stage procedure, and GEM algorithm to obtain the constrained MLEsδ˜, π˜1, and ρ˜.


### 3.2.1. Constrained MLEs Based on Fisher Scoring Algorithm

The initial values ofρ,π1 are defined in equation (9), and δ0=1. The t+1-th updates of δ, π1, and ρ can be calculated by Fisher scoring algorithm as follows:(23)δt+1π1t+1ρt+1=δtπ1tρt+I2−1δ,π1,ρ×∂ℓ0∂δ∂ℓ0∂π1∂ℓ0∂ρ|δ=δt,π1=π1t,ρ=ρt,where I2 is a 3×3 Fisher information matrix (see Appendix B).


### 3.2.2. Constrained MLEs Based on Two-Stage Procedure

The two-stage procedure is different from the two-step method in Section3.1.2. Firstly, the MLE δ˜ of δ is given by Newton–Raphson algorithm. Then, π˜1 and ρ˜ are obtained by Fisher scoring algorithm under given MLE δ˜. The detailed process is described as follows.(i)
The initial valuesπ10 and ρ0 are defined in the equation (9), and δ0=1. The t+1-th approximate δt+1 of δ is obtained by Newton–Raphson algorithm:(24)δt+1=δt−∂2ℓ0∂δ2−1∂ℓ0∂δ|δ=δt,π1=π1t,ρ=ρt,where(25)∂2ℓ0∂δ2=∑i=2gm0iρ−12δ+2π1−1+1δ+π1−12δ+π11−ρ−1−2ρ−1δ+π1−1δ+π11−ρ−1−ρ−1ρ−12δ+2π1−1+1δ+π1−1δ+π11−ρ−12+m2i−2δ+π11−ρ+ρ2δ+π121−ρ+ρδ+π12−2ρ−1δ+π121−ρ+ρδ+π1+m1i2δ+π1δ+π1−1−2δ+2π1−1δ+π1δ+π1−12−2δ+2π1−1δ+π12δ+π1−1.(ii)
Givenδ=δ˜, the t+1-th approximates of π1 and ρ can be calculated as(26)π1t+1ρt+1=π1tρt+I3−1π1,ρ×∂ℓ0∂π1∂ℓ0∂ρ|δ=δ˜,π1=π1t,ρ=ρt,whereI3 is a 2×2 Fisher information matrix of π1 and ρ (see Appendix B).


### 3.2.3. Constrained MLEs Based on GEM Algorithm

Similar to global MLEs, the GEM algorithm is used to calculate constrained MLEs underH0. The initial values of ρ,π1 are defined in equation (9), and δ0=∑i=2gπi0−π10/g−1. The detailed process is shown as follows.(i)
E Step. Let(27)Qπ1,δ,ρ|m,π1t,δt,ρt=Ezℓ0π1,δ,ρ|m,z=Ez1l01π1,ρ|m1,z1+∑i=2gEzil0iπ1,δ,ρ|mi,zi,where(28)Ez1l01π1,ρ|m1,z1=Ez11z11|m01,π1t,ρtlogρ1−π1+m01−Ez11z11|m01,π1t,ρtlog1−ρ1−π12+m11log2π11−ρ1−π1+Ez21z21|m21,π1t,ρtlogρπ1+m21−Ez21z21|m21,π1t,ρtlog1−ρπ12,Ezil0iπ1,δ,ρ|mi,zi=Ez1iz1i|m0i,π1t,δt,ρtlogρ1−π1−δ+m0i−Ez1iz1i|m0i,π1t,δt,ρtlog1−ρ1−π1−δ2+m1ilog2π1+δ1−ρ1−π1−δ+m2i−Ez2iz2i|m2i,π1t,δt,ρtlog1−ρπ1+δ2+Ez2iz2i|m2i,π1t,δt,ρtlogρπ1+δ,i=2,…,g,Ez11z11|m01,π1t,ρt=m01ρt1−π1tρt1−π1t+1−ρt1−π1t2,Ez21z21|m21,π1t,ρt=m21ρtπ1tρtπ1t+1−ρtπ1t2,Ez1iz1i|m0i,π1t,δt,ρt=m0iρt1−π1t−δtρt1−π1t−δt+1−ρt1−π1t−δt2,Ez2iz2i|m2i,π1t,δt,ρt=m2iρtπ1t+δtρtπ1t+δt+1−ρtπ1t+δt2,i=2,…,g.(ii)
M Step. Similar to Section 3.1.3, we choose parameters π1,δ,ρ so that Qπ1,δ,ρ|m,π1t,δt,ρt increases. Repeat E and M steps until the result converges.


## 3.2.1. Constrained MLEs Based on Fisher Scoring Algorithm

The initial values ofρ,π1 are defined in equation (9), and δ0=1. The t+1-th updates of δ, π1, and ρ can be calculated by Fisher scoring algorithm as follows:(23)δt+1π1t+1ρt+1=δtπ1tρt+I2−1δ,π1,ρ×∂ℓ0∂δ∂ℓ0∂π1∂ℓ0∂ρ|δ=δt,π1=π1t,ρ=ρt,where I2 is a 3×3 Fisher information matrix (see Appendix B).


## 3.2.2. Constrained MLEs Based on Two-Stage Procedure

The two-stage procedure is different from the two-step method in Section3.1.2. Firstly, the MLE δ˜ of δ is given by Newton–Raphson algorithm. Then, π˜1 and ρ˜ are obtained by Fisher scoring algorithm under given MLE δ˜. The detailed process is described as follows.(i)
The initial valuesπ10 and ρ0 are defined in the equation (9), and δ0=1. The t+1-th approximate δt+1 of δ is obtained by Newton–Raphson algorithm:(24)δt+1=δt−∂2ℓ0∂δ2−1∂ℓ0∂δ|δ=δt,π1=π1t,ρ=ρt,where(25)∂2ℓ0∂δ2=∑i=2gm0iρ−12δ+2π1−1+1δ+π1−12δ+π11−ρ−1−2ρ−1δ+π1−1δ+π11−ρ−1−ρ−1ρ−12δ+2π1−1+1δ+π1−1δ+π11−ρ−12+m2i−2δ+π11−ρ+ρ2δ+π121−ρ+ρδ+π12−2ρ−1δ+π121−ρ+ρδ+π1+m1i2δ+π1δ+π1−1−2δ+2π1−1δ+π1δ+π1−12−2δ+2π1−1δ+π12δ+π1−1.(ii)
Givenδ=δ˜, the t+1-th approximates of π1 and ρ can be calculated as(26)π1t+1ρt+1=π1tρt+I3−1π1,ρ×∂ℓ0∂π1∂ℓ0∂ρ|δ=δ˜,π1=π1t,ρ=ρt,whereI3 is a 2×2 Fisher information matrix of π1 and ρ (see Appendix B).


## 3.2.3. Constrained MLEs Based on GEM Algorithm

Similar to global MLEs, the GEM algorithm is used to calculate constrained MLEs underH0. The initial values of ρ,π1 are defined in equation (9), and δ0=∑i=2gπi0−π10/g−1. The detailed process is shown as follows.(i)
E Step. Let(27)Qπ1,δ,ρ|m,π1t,δt,ρt=Ezℓ0π1,δ,ρ|m,z=Ez1l01π1,ρ|m1,z1+∑i=2gEzil0iπ1,δ,ρ|mi,zi,where(28)Ez1l01π1,ρ|m1,z1=Ez11z11|m01,π1t,ρtlogρ1−π1+m01−Ez11z11|m01,π1t,ρtlog1−ρ1−π12+m11log2π11−ρ1−π1+Ez21z21|m21,π1t,ρtlogρπ1+m21−Ez21z21|m21,π1t,ρtlog1−ρπ12,Ezil0iπ1,δ,ρ|mi,zi=Ez1iz1i|m0i,π1t,δt,ρtlogρ1−π1−δ+m0i−Ez1iz1i|m0i,π1t,δt,ρtlog1−ρ1−π1−δ2+m1ilog2π1+δ1−ρ1−π1−δ+m2i−Ez2iz2i|m2i,π1t,δt,ρtlog1−ρπ1+δ2+Ez2iz2i|m2i,π1t,δt,ρtlogρπ1+δ,i=2,…,g,Ez11z11|m01,π1t,ρt=m01ρt1−π1tρt1−π1t+1−ρt1−π1t2,Ez21z21|m21,π1t,ρt=m21ρtπ1tρtπ1t+1−ρtπ1t2,Ez1iz1i|m0i,π1t,δt,ρt=m0iρt1−π1t−δtρt1−π1t−δt+1−ρt1−π1t−δt2,Ez2iz2i|m2i,π1t,δt,ρt=m2iρtπ1t+δtρtπ1t+δt+1−ρtπ1t+δt2,i=2,…,g.(ii)
M Step. Similar to Section 3.1.3, we choose parameters π1,δ,ρ so that Qπ1,δ,ρ|m,π1t,δt,ρt increases. Repeat E and M steps until the result converges.


## 3.3. Likelihood Ratio Test

Likelihood ratio test statistic can be constructed through the global and constrained MLEs as follows:(29)TL=2ℓπ^,ρ^−ℓ0π˜1,δ˜,ρ˜=2∑i=1gliπ^i,ρ^−l01π˜1,ρ˜−∑i=2gl0iπ˜1,δ˜,ρ˜where π^i,ρ^i=1,2,…,g are the global MLEs and π˜1,δ˜,ρ˜ are the constrained MLEs.


## 3.4. Score Test

Note thatH0:δ2=⋯=δg≜δ is equivalent to π2=⋯=πg. The homogeneity test of many-to-one risk differences can be achieved by testing the equality of risks in treatment groups. Denote U=∂l2/∂π2,…,∂lg/∂πg,0,0. Score test statistic is derived as(30)TSC=UI4π2,…,πg,π1,ρ−1UT|π=π˜,ρ=ρ˜,where π=π1,π2,…,πgT and π˜=π˜1,π˜2,…,π˜gT. π˜i and ρ˜ are constrained MLEs under H0. I4 is a g+1×g+1 information matrix (see Appendix C for more information). Thus, TSC can be simplified as(31)TSC=∑i=2gui2π˜i,ρ˜aiπ˜i,ρ˜+∑i=2guiπ˜i,ρ˜biπ˜i,ρ˜/aiπ˜i,ρ˜2Dπ˜,ρ˜−∑i=1gbi2π˜i,ρ˜/aiπ˜i,ρ˜,where(32)uiπi,ρ=∂li∂πi=2πi−1m1iπiπi−1+m2iρ+2πi−2ρπiπiρ+πi−ρπi−m0iρ+2πi−2ρπi−2πi−1ρπi−πi+1,aiπi,ρ=E−∂2li∂πi2=Ni−4ρ2πi2+4ρ2πi−ρ2+6ρπi2−6ρπi+2ρ−2πi2+2πiπi1−πiρ+πi−ρπiρπi−πi+1,biπi,ρ=E−∂2li∂πi∂ρ=Niρ2πi−1ρ+πi−ρπiρπi−πi+1,Dπ,ρ=E−∂2ℓ∂ρ2=∑i=1gNiπiρ+11−πi1−ρρ+πi−ρπiρπi−πi+1.


## 3.5. Wald-Type Test

Letβ=π2,…,πg,π1,ρ and(33)C=1−10……00001−1000⋮⋱⋱⋱⋮⋮⋮⋮⋱⋱⋱⋮⋮⋮0……01−100g−2×g+1.The null hypothesisH0:δ2=⋯=δg≜δ is equivalent to CβT=0. Thus, Wald-type statistic is(34)TW=βCTCI4−1CT−1CβT|π=π^,ρ=ρ^,where the Fisher information matrix I4 is the same as that of score test (see Appendix D). It can be simplified as(35)TW=∑i=1g−2∑j=1g−2π^i+1−π^i+2π^j+1−π^j+2∑k=1g−2rkiπ^,ρ^rkjπ^,ρ^,where(36)rijπ,ρ=−1liiπ,ρ∑k=ji−1likπ,ρrkjπ,ρ,i>j,1lijπ,ρ,i=j,0,i<j,lijπ,ρ=fijπ,ρ−∑k=1j−1likπ,ρljkπ,ρljjπ,ρ,i>j,fijπ,ρ−∑k=1j−1ljkπ,ρ21/2,i=j,0,i<j,fijπ,ρ=Iijπ,ρ−Ii+1jπ,ρ−Iij+1π,ρ−Ii+1j+1π,ρ,Iijπ,ρ=1ai+1πi+1,ρ+bi+1πi+1,ρbj+1πj+1,ρ/ai+1πi+1,ρaj+1πj+1,ρDπ,ρ−∑i=1gbi2πi,ρ/aiπi,ρ,i=j,bi+1πi+1,ρbj+1πj+1,ρ/ai+1πi+1,ρaj+1πj+1,ρDπ,ρ−∑i=1gbi2πi,ρ/aiπi,ρ,i≠j.UnderH0, test statistics TL,TSC, and TW are asymptotically distributed as chi-square distribution with g−2 degrees of freedom. Thus, H0 should be rejected if the value of test statistic is larger than χg−2,1−α2 at the significant level α, where χg−2,1−α2 is the 1001−α percentile of the chi-square distribution with g−2 degrees of freedom.


## 4. Monte Carlo Simulations

In this section, the performance of several algorithms are compared with respect to average errors of MLEs, the number of iteration, and time cost. For convenience, we denote Fisher scoring algorithm, two-step method, and GEM algorithm for global MLEs as FSA, TSM, and GEM and Fisher scoring algorithm, two-stage procedure, and GEM algorithm for constrained MLEs as FSA, TSP, and GEM for tables and figures. Then, we investigate the type I error rates (TIEs) and power of the likelihood ratio, score, and Wald-type tests. In simulations,g and N=N1,N2,…,NgT are arranged as shown in Table 2, where the scenarios 4, 8, and 12 are unbalanced designs.Table 2
Scenarios of simulation.
gScenarioN1,N2,…,Ng3130,30,30270,70,703100,100,100430,70,1006530,30,30,30,30,30670,70,70,70,70,707100,100,100,100,100,100830,40,50,60,70,809930,30,30,30,30,30,30,30,301070,70,70,70,70,70,70,70,7011100,100,100,100,100,100,100,100,1001230,40,50,60,70,80,90,100,110


### 4.1. Selection of Algorithms

UnderHa or H0, we randomly select 1000 sets of πii=1,2,…,g and ρ for each scenario in Table 2. Further, 10,000 samples are randomly produced for each parameter setting.


#### 4.1.1. Evaluation of MLEs

Lete¯i be averaged error among estimators and true values of parameters from the 10000 random samples under the i-th parameter setting. Denote e¯¯=1/1000∑i=11000e¯i in each scenario, which is the total averaged error for 1000 parameter settings. The dispersion of e¯i in each scenario can be reflected by standard deviation (SD) value of e¯i. The MSEs of the global MLEs can be evaluated by the differences from the true parameters to the corresponding estimated values under various parameter settings. The global MLEs can be calculated based on Fisher scoring algorithm, two-step method, and GEM algorithm. The constrained MLEs are obtained by Fisher scoring algorithm, two-stage procedure, and GEM algorithm. The random samples for the former are generated under Ha and the latter are generated under H0. The convergence accuracy ϵ is defined by the differences from two close iterations and fixed as 1×10−5. The MSEs of the three algorithms for global MLEs have no significant difference as shown in Tables 3–5 . That is to say, the global MLEs are identical by these algorithms. In Tables 6–8 , the values of e¯¯, SD, and MSEs from Fisher scoring algorithm are usually smaller than other two algorithms for constrained MLEs. So, Fisher scoring algorithm has higher accuracy for constrained MLEs. All MSEs become smaller and close to each other when sample size increases. Algorithms for MLEs have better MSEs in balanced designs than unbalanced designs.Table 3
Errors of the global MLEs forg=3.
IndexMLEScenario 1Scenario 2Scenario 3Scenario 4FSATSMGEMFSATSMGEMFSATSMGEMFSATSMGEMe¯¯π^1−0.000071−0.000109−0.000076−0.000037−0.000125−0.000036−0.000058−0.000078−0.000064−0.000544−0.001051−0.000616π^2−0.000030−0.000043−0.000026−0.000147−0.000243−0.000147−0.000030−0.000049−0.000030−0.000081−0.000584−0.000116π^30.0000250.0000110.0000280.000036−0.0000690.000036−0.000047−0.000061−0.000054−0.000373−0.000824−0.000385ρ^−0.007755−0.006677−0.007174−0.003747−0.003360−0.003342−0.002615−0.002376−0.002213−0.046909−0.040828−0.040980SD ofe¯iπ^10.0024030.0024790.0024540.0014920.0016150.0015230.0012540.0012770.0012810.0024010.0024660.002459π^20.0023470.0024130.0024020.0015420.0016800.0015700.0012510.0012850.0012740.0017760.0018420.001825π^30.0023620.0024810.0024260.0014790.0017010.0015090.0012550.0012740.0012840.0014800.0015600.001515ρ^0.0119820.0158990.0138700.0062040.0079190.0078240.0047850.0059110.0064770.0086150.0109310.010636MSEπ^10.0054280.0054380.0054280.0022900.0023190.0022890.0015990.0016050.0015990.0053270.0053450.005327π^20.0053320.0053390.0053320.0022860.0023210.0022860.0016090.0016160.0016090.0031920.0032050.003192π^30.0053320.0053380.0053310.0022780.0023140.0022780.0016020.0016080.0016020.0022890.0023030.002289ρ^0.0093170.0093170.0093140.0040990.0040980.0040970.0028360.0028360.0028350.0057160.0057150.005714Table 4
Errors of the global MLEs forg=6.
IndexMLEScenario 5Scenario 6Scenario 7Scenario 8FSATSMGEMFSATSMGEMFSATSMGEMFSATSMGEMe¯¯π^10.000010−0.0000040.0000110.0000120.0000080.000016−0.000064−0.000064−0.000064−0.000022−0.000026−0.000029π^20.000009−0.0000070.000009−0.000013−0.000021−0.0000130.0000670.0000670.000066−0.000027−0.000029−0.000027π^30.000006−0.0000090.000007−0.000030−0.000036−0.000030−0.000038−0.000039−0.000039−0.000096−0.000097−0.000095π^40.0000890.0000800.000092−0.000068−0.000075−0.0000700.0000770.0000770.0000770.0000470.0000450.000047π^5−0.000121−0.000134−0.0001200.0000560.0000480.0000560.0000510.0000510.000051−0.000010−0.000012−0.000011π^60.0000620.0000480.000063−0.000011−0.000019−0.0000150.0000550.0000550.000056−0.000071−0.000075−0.000073ρ^−0.011639−0.011557−0.011491−0.005095−0.005052−0.004959−0.003580−0.003578−0.003538−0.006436−0.006387−0.006288SD ofe¯iπ^10.0023780.0023950.0023950.0015390.0015430.0015540.0012560.0012560.0012580.0022980.0023020.002313π^20.0023600.0023780.0023750.0014900.0014990.0015000.0013170.0013180.0013210.0019840.0019910.002004π^30.0023530.0023670.0023760.0015220.0015290.0015380.0012650.0012660.0012700.0017570.0017610.001777π^40.0024580.0024700.0024850.0014730.0014780.0014850.0012790.0012790.0012820.0015710.0015750.001583π^50.0023160.0023330.0023430.0014610.0014700.0014700.0012490.0012500.0012530.0014990.0015010.001512π^60.0022940.0023080.0023110.0014950.0015060.0015030.0012460.0012460.0012480.0013800.0013880.001395ρ^0.0049620.0055200.0058990.0026080.0029460.0035620.0016870.0017110.0018730.0033620.0037200.004300MSEπ^10.0053530.0053590.0053530.0023000.0023030.0023000.0016040.0016040.0016040.0053020.0053020.005302π^20.0053870.0053930.0053870.0022940.0022970.0022940.0015970.0015970.0015970.0040130.0040140.004013π^30.0053310.0053360.0053300.0022930.0022950.0022930.0016140.0016140.0016140.0032410.0032420.003241π^40.0053790.0053830.0053790.0022800.0022830.0022800.0015970.0015970.0015970.0026800.0026810.002680π^50.0053730.0053780.0053720.0022940.0022960.0022940.0016080.0016080.0016080.0023000.0023010.002300π^60.0053020.0053070.0053020.0023020.0023040.0023020.0016170.0016170.0016170.0020050.0020060.002005ρ^0.0049070.0049070.0049060.0020630.0020630.0020620.0014300.0014300.0014300.0026400.0026400.002639Table 5
Errors of the global MLEs forg=9.
IndexMLEScenario 9Scenario 10Scenario 11Scenario 12FSATSMGEMFSATSMGEMFSATSMGEMFSATSMGEMe¯¯π^1−0.000013−0.000018−0.0000140.0000080.0000080.000010−0.000043−0.000043−0.000045−0.000006−0.000006−0.000003π^2−0.000106−0.000106−0.0001030.0000470.0000470.000048−0.000009−0.000009−0.0000090.0000080.0000080.000010π^3−0.000043−0.000050−0.000048−0.000065−0.000066−0.000065−0.000029−0.000029−0.0000270.0000960.0000960.000093π^40.0000180.0000140.0000190.0000040.0000040.000003−0.000010−0.000010−0.000010−0.000052−0.000053−0.000054π^50.0000300.0000260.000032−0.000021−0.000021−0.000021−0.000024−0.000024−0.000024−0.000024−0.000024−0.000022π^6−0.000012−0.000016−0.000013−0.000021−0.000021−0.0000210.0000710.0000710.000072−0.000118−0.000118−0.000118π^7−0.000106−0.000108−0.000102−0.000036−0.000036−0.0000350.0000770.0000770.000077−0.000019−0.000020−0.000019π^8−0.000151−0.000155−0.000151−0.000005−0.000006−0.0000050.0000120.0000130.000014−0.000022−0.000022−0.000024π^9−0.000043−0.000044−0.0000430.0000580.0000570.000058−0.000008−0.000008−0.000008−0.000018−0.000018−0.000017ρ^−0.011963−0.011923−0.011885−0.005189−0.005188−0.005158−0.003624−0.003622−0.003582−0.005150−0.005138−0.005082SD ofe¯iπ^10.0023440.0023540.0023570.0015110.0015110.0015170.0013550.0013550.0013580.0024190.0024200.002428π^20.0023440.0023450.0023530.0015800.0015800.0015830.0012930.0012930.0012950.0020620.0020630.002073π^30.0022720.0022800.0022790.0015780.0015770.0015810.0012810.0012810.0012860.0017370.0017380.001745π^40.0023540.0023620.0023620.0015180.0015180.0015200.0012700.0012700.0012740.0016120.0016120.001616π^50.0023440.0023500.0023500.0014620.0014630.0014640.0012590.0012590.0012600.0014880.0014880.001491π^60.0022880.0022900.0022980.0015120.0015130.0015140.0012840.0012840.0012860.0014570.0014570.001461π^70.0023080.0023140.0023200.0015130.0015130.0015170.0012630.0012630.0012650.0013400.0013400.001343π^80.0022470.0022500.0022570.0014970.0014970.0014980.0012320.0012320.0012340.0012370.0012380.001241π^90.0024150.0024160.0024250.0015680.0015680.0015690.0012680.0012680.0012750.0012060.0012060.001208ρ^0.0032640.0036380.0038960.0015590.0015750.0017100.0013130.0013360.0015720.0016550.0017660.002223MSEπ^10.0053270.0053290.0053270.0023270.0023270.0023270.0016180.0016180.0016180.0053700.0053700.005370π^20.0053650.0053660.0053650.0023580.0023580.0023580.0016030.0016030.0016030.0040340.0040340.004034π^30.0054050.0054070.0054050.0023220.0023230.0023220.0016130.0016130.0016130.0032050.0032050.003205π^40.0053320.0053330.0053320.0023120.0023120.0023120.0016220.0016220.0016210.0027050.0027050.002705π^50.0053930.0053940.0053930.0023070.0023080.0023070.0016180.0016180.0016180.0023250.0023250.002325π^60.0053930.0053940.0053930.0023390.0023390.0023390.0016090.0016090.0016090.0020150.0020150.002015π^70.0053440.0053450.0053440.0023140.0023150.0023140.0016180.0016180.0016180.0018110.0018110.001811π^80.0053700.0053710.0053700.0023250.0023260.0023250.0015990.0015990.0015990.0016270.0016270.001627π^90.0053750.0053760.0053750.0023240.0023240.0023240.0016230.0016230.0016230.0014710.0014710.001471ρ^0.0033220.0033220.0033210.0013550.0013550.0013550.0009490.0009490.0009490.0013670.0013670.001367Table 6
Errors of the constrained MLEs forg=3.
IndexMLEScenario 1Scenario 2Scenario 3Scenario 4FSATSPGEMFSATSPGEMFSATSPGEMFSATSPGEMe¯¯π˜10.0000220.0000100.0000100.000002−0.000002−0.000002−0.000019−0.000017−0.000017−0.000007−0.000005−0.000005δ˜−0.000028−0.000089−0.000643−0.000019−0.000001−0.0002020.0000190.000018−0.001885−0.0000060.0000120.000919ρ˜−0.0089810.0241310.024078−0.0039950.0258560.025855−0.0028050.0284250.028424−0.0042020.0158390.015835SD ofe¯iπ˜10.0010570.0010570.0010400.0005360.0005440.0005400.0004300.0004350.0004330.0009660.0008040.000804δ˜0.0014920.0031360.0805220.0006690.0024520.0788930.0005540.0027880.0821190.0010030.0019080.081664ρ˜0.0016280.0352850.0350900.0010060.0338350.0338190.0007700.0367030.0366950.0010870.0220450.022035MSEπ˜10.0052270.0052640.0052640.0022180.0022340.0022340.0015640.0015750.0015750.0051680.0052180.005218δ˜0.0078360.0088370.0100140.0033340.0033910.0077150.0023490.0023490.0077940.0060980.0061430.009410ρ˜0.0096060.0104200.0104200.0042170.0056540.0056540.0029230.0047980.0047980.0044240.0048760.004876Table 7
Errors of the constrained MLEs forg=6.
IndexMLEScenario 5Scenario 6Scenario 7Scenario 8FSATSPGEMFSATSPGEMFSATSPGEMFSATSPGEMe¯¯π˜1−0.000041−0.000029−0.000029−0.000018−0.000013−0.000013−0.000011−0.000010−0.000010−0.000037−0.000030−0.000030δ˜0.0000400.0001250.0022680.0000200.0001050.0009820.000014−0.000047−0.0002350.0000400.0000520.001381ρ˜−0.0045730.0169950.016998−0.0019520.0165970.016597−0.0013990.0188780.018878−0.0025440.0109350.010935SD ofe¯iπ˜10.0009600.0008390.0008350.0005060.0004920.0004900.0004260.0004200.0004190.0009660.0007810.000781δ˜0.0010810.0023510.0408530.0005630.0017040.0392400.0004720.0019020.0411400.0009890.0014530.041265ρ˜0.0009850.0254330.0254290.0006190.0225970.0225960.0004750.0256470.0256470.0007400.0156970.015697MSEπ˜10.0051820.0052310.0052310.0022430.0022640.0022640.0015760.0015890.0015890.0051460.0052020.005202δ˜0.0062350.0062230.0060460.0026950.0026910.0034290.0018930.0018890.0030180.0056720.0056930.005703ρ˜0.0048620.0054650.0054650.0020780.0027410.0027410.0014590.0023790.0023790.0026810.0029270.002927Table 8
Errors of the constrained MLEs forg=9.
IndexMLEScenario 9Scenario 10Scenario 11Scenario 12FSATSPGEMFSATSPGEMFSATSPGEMFSATSPGEMe¯¯π˜10.000003−0.000005−0.0000050.0000090.0000070.000007−0.000003−0.000007−0.000007−0.000032−0.000010−0.000010δ˜0.000006−0.000017−0.000751−0.000015−0.000020−0.000252−0.000003−0.000036−0.0004990.0000350.0000190.001520ρ˜−0.0030090.0119710.011972−0.0013540.0135960.013596−0.0009230.0136030.013604−0.0013330.0062050.006205SD ofe¯iπ˜10.0009690.0008270.0008250.0005110.0004910.0004900.0004240.0004120.0004120.0010140.0008070.000807δ˜0.0010220.0017150.0268040.0005460.0013680.0271700.0004450.0012380.0268980.0010250.0011140.027407ρ˜0.0008100.0177770.0177750.0004620.0180450.0180450.0003840.0172870.0172870.0004720.0093520.009352MSEπ˜10.0051900.0052430.0052430.0022120.0022350.0022350.0015470.0015620.0015620.0051250.0051840.005184δ˜0.0058470.0058600.0053830.0024910.0024960.0027240.0017440.0017480.0021140.0053870.0054280.005078ρ˜0.0032050.0035120.0035120.0014090.0018520.0018520.0009770.0014180.0014180.0013990.0014890.001489Tables9 and 10 show the number of failures for these algorithms to converge when ɛ=1×10−5,1×10−4 within k=100,200 iterations. Fisher scoring algorithm has lower failure rate for convergence when calculating global MLEs. GEM algorithm for global MLEs is sensitive to ϵ and k. Both reducing ϵ and adding k can markedly improve the convergence possibility of GEM algorithm for global MLEs. Fisher scoring and GEM algorithms for constrained MLEs hardly fail to converge within 100 iterations. As g increases, the number of failures becomes small.Table 9
Failures of algorithms for global MLEs converging withink iterations.
Algorithmϵkg=3g=6g=9Scenario 1Scenario 2Scenario 3Scenario 4Scenario 5Scenario 6Scenario 7Scenario 8Scenario 9Scenario 10Scenario 11Scenario 12FSA1×10−51001000000000001×10−52001000000000001×10−4100100000000000TSM1×10−510024331713508310611612595453483256112014336410389607298683160321×10−520024331713508310611612595453483256112014336410389607298683160321×10−41002363361332701040181240325266225070199153580937998713967425896GEM1×10−510023356419944218253318256794508850986570386451498113777131859224371×10−5200422083503833405350592128619705164261915412860125751142299411×10−4100616118469410010000Table 10
Failures of algorithms for constrained MLEs converging withink iterations.
Algorithmϵkg=3g=6g=9Scenario 1Scenario 2Scenario 3Scenario 4Scenario 5Scenario 6Scenario 7Scenario 8Scenario 9Scenario 10Scenario 11Scenario 12FSA1×10−51003000000000001×10−52001000000000001×10−4100100000000000TSP1×10−5100133968540653993922555061923128152476953143415364105504041×10−5200133968540653993922555061923128152476953143415364105504041×10−410013396854065399392255506192312815247695314341536410550404GEM1×10−51000000000000001×10−52000000000000001×10−4100000000000000


#### 4.1.2. Number of Iteration

Since the algorithms for global MLEs have identical accuracy in terms of MSEs, we further compare their efficiency by calculating the number of iteration and time cost when the convergence accuracy isϵ=1×10−5. The average number of iteration is recorded for every parameter setting in Figure 1.Figure 1
Number of iteration of algorithms for global MLEs from 1000 parameter settings.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)As shown in Figure1, GEM algorithm takes more iterations to converge than Fisher scoring algorithm and two-step method. The number of iterations is most intensive between 20 and 30 and the upper limit increases when Nii=1,2,…,g is bigger. For Fisher scoring algorithm, the number of iterations has slight superiority over two-step method.


#### 4.1.3. Time Cost

The time required to achieve convergence can also reflect the convergence rate of an algorithm. Figure2 shows that GEM algorithm has the worst performance in terms of time cost. Fisher scoring algorithm always takes less time, especially when g=6,9. As Nii=1,2,…,g increases, the time distribution is more clustered for Fisher scoring algorithm and two-step method. When g is larger, difference of time cost is bigger between the two algorithms.Figure 2
Time cost of algorithms for global MLEs from 1000 parameter settings.
(a)(b)(c)Based on the above results, Fisher scoring algorithm, two-step method, and GEM algorithm can produce the same global MLEs when the number of iterations is large enough. However, Fisher scoring algorithm for global MLEs has a better convergence rate in terms of number of iteration and time cost. Constrained MLEs from Fisher scoring algorithm has smaller MSEs. Therefore, it is advisable to choose Fisher scoring algorithms to obtain the global and constrained MLEs required by the homogeneity test in this paper.


### 4.2. Evaluation of Test Statistics

Since Fisher scoring algorithm has higher accuracy for global MLEs and is more efficient for constrained MLEs, we construct test statisticsTL,TSC, and TW based on global and constrained MLEs from the algorithm. The empirical TIEs and power of three test statistics can be demonstrated by Monte Carlo simulations. All tests are conducted at the significance level of α=0.05.


#### 4.2.1. Type I Error Rate

TIE is the probability of rejecting a null hypothesis that is true. According to Tang et al. [2], a test is robust if its empirical TIE is between 0.04 and 0.06; liberal if it is greater than 0.06; and conservative if it is less than 0.04 when the significance level is α=0.05. Power is the probability of correctly rejecting a null hypothesis when it is false in a statistical test. A good test should not only be robust but also make power as high as possible.1,000 random parameter settings involvingπi and ρi=1,2,…,g are generated under H0, where g and N are shown in Table 2. 10,000 samples are randomly produced for every parameter setting. The empirical TIEs can be computed by dividing the number of times of rejecting H0 with 10,000. Figure 3 further reflects the comprehensive performance of TL,TSC, and TW. The empirical TIEs of TSC are smaller and close to 0.05, which means score test is more robust. Wald-type test statistic tends to be liberal when Nii=1,2,…,g is small (i.e., Ni=30,70) and unbalanced. Wald-type test statistic is more liberal when g is bigger. Empirical TIEs of three tests are closer to 0.05 as Ni increases. Balanced designs make test statistics have better performance than unbalanced designs.Figure 3
Empirical TIEs of test statistics.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)


#### 4.2.2. Power

LetN1=N2=⋯=Ng=m, g=3,6,9, π1=0.2, ρ=0.2,0.5,0.8, δp=0.1,0.2,0.3, πi=π1+δp when i=2k, and πi=π1 when i=2k+1k=1,2,…. 10,000 samples will be randomly generated under every parameter setting when m changes from 10 to 160 at intervals of 10. Empirical power can be computed by dividing the number of times of rejecting H0 with 10,000. Figure 4 reflects how empirical power changes as m changes. TW has higher power and TL and TSC have close power. Empirical power of three test statistics is close to each other as m increases. Increasing the number of groups can produce higher power. It means that all test statistics work well under multigroup cases. The moderately and highly relevant data (i.e., ρ=0.5 and ρ = 0.8) have lower power than the mildly relevant data (i.e., ρ = 0.2).Figure 4
Empirical power of test statistics.
(a)(b)(c)According to the above results, score test is more robust and has satisfactory power. Thus, it should be recommended for the homogeneity test about many-to-one comparison of risk differences.


## 4.1. Selection of Algorithms

UnderHa or H0, we randomly select 1000 sets of πii=1,2,…,g and ρ for each scenario in Table 2. Further, 10,000 samples are randomly produced for each parameter setting.


### 4.1.1. Evaluation of MLEs

Lete¯i be averaged error among estimators and true values of parameters from the 10000 random samples under the i-th parameter setting. Denote e¯¯=1/1000∑i=11000e¯i in each scenario, which is the total averaged error for 1000 parameter settings. The dispersion of e¯i in each scenario can be reflected by standard deviation (SD) value of e¯i. The MSEs of the global MLEs can be evaluated by the differences from the true parameters to the corresponding estimated values under various parameter settings. The global MLEs can be calculated based on Fisher scoring algorithm, two-step method, and GEM algorithm. The constrained MLEs are obtained by Fisher scoring algorithm, two-stage procedure, and GEM algorithm. The random samples for the former are generated under Ha and the latter are generated under H0. The convergence accuracy ϵ is defined by the differences from two close iterations and fixed as 1×10−5. The MSEs of the three algorithms for global MLEs have no significant difference as shown in Tables 3–5 . That is to say, the global MLEs are identical by these algorithms. In Tables 6–8 , the values of e¯¯, SD, and MSEs from Fisher scoring algorithm are usually smaller than other two algorithms for constrained MLEs. So, Fisher scoring algorithm has higher accuracy for constrained MLEs. All MSEs become smaller and close to each other when sample size increases. Algorithms for MLEs have better MSEs in balanced designs than unbalanced designs.Table 3
Errors of the global MLEs forg=3.
IndexMLEScenario 1Scenario 2Scenario 3Scenario 4FSATSMGEMFSATSMGEMFSATSMGEMFSATSMGEMe¯¯π^1−0.000071−0.000109−0.000076−0.000037−0.000125−0.000036−0.000058−0.000078−0.000064−0.000544−0.001051−0.000616π^2−0.000030−0.000043−0.000026−0.000147−0.000243−0.000147−0.000030−0.000049−0.000030−0.000081−0.000584−0.000116π^30.0000250.0000110.0000280.000036−0.0000690.000036−0.000047−0.000061−0.000054−0.000373−0.000824−0.000385ρ^−0.007755−0.006677−0.007174−0.003747−0.003360−0.003342−0.002615−0.002376−0.002213−0.046909−0.040828−0.040980SD ofe¯iπ^10.0024030.0024790.0024540.0014920.0016150.0015230.0012540.0012770.0012810.0024010.0024660.002459π^20.0023470.0024130.0024020.0015420.0016800.0015700.0012510.0012850.0012740.0017760.0018420.001825π^30.0023620.0024810.0024260.0014790.0017010.0015090.0012550.0012740.0012840.0014800.0015600.001515ρ^0.0119820.0158990.0138700.0062040.0079190.0078240.0047850.0059110.0064770.0086150.0109310.010636MSEπ^10.0054280.0054380.0054280.0022900.0023190.0022890.0015990.0016050.0015990.0053270.0053450.005327π^20.0053320.0053390.0053320.0022860.0023210.0022860.0016090.0016160.0016090.0031920.0032050.003192π^30.0053320.0053380.0053310.0022780.0023140.0022780.0016020.0016080.0016020.0022890.0023030.002289ρ^0.0093170.0093170.0093140.0040990.0040980.0040970.0028360.0028360.0028350.0057160.0057150.005714Table 4
Errors of the global MLEs forg=6.
IndexMLEScenario 5Scenario 6Scenario 7Scenario 8FSATSMGEMFSATSMGEMFSATSMGEMFSATSMGEMe¯¯π^10.000010−0.0000040.0000110.0000120.0000080.000016−0.000064−0.000064−0.000064−0.000022−0.000026−0.000029π^20.000009−0.0000070.000009−0.000013−0.000021−0.0000130.0000670.0000670.000066−0.000027−0.000029−0.000027π^30.000006−0.0000090.000007−0.000030−0.000036−0.000030−0.000038−0.000039−0.000039−0.000096−0.000097−0.000095π^40.0000890.0000800.000092−0.000068−0.000075−0.0000700.0000770.0000770.0000770.0000470.0000450.000047π^5−0.000121−0.000134−0.0001200.0000560.0000480.0000560.0000510.0000510.000051−0.000010−0.000012−0.000011π^60.0000620.0000480.000063−0.000011−0.000019−0.0000150.0000550.0000550.000056−0.000071−0.000075−0.000073ρ^−0.011639−0.011557−0.011491−0.005095−0.005052−0.004959−0.003580−0.003578−0.003538−0.006436−0.006387−0.006288SD ofe¯iπ^10.0023780.0023950.0023950.0015390.0015430.0015540.0012560.0012560.0012580.0022980.0023020.002313π^20.0023600.0023780.0023750.0014900.0014990.0015000.0013170.0013180.0013210.0019840.0019910.002004π^30.0023530.0023670.0023760.0015220.0015290.0015380.0012650.0012660.0012700.0017570.0017610.001777π^40.0024580.0024700.0024850.0014730.0014780.0014850.0012790.0012790.0012820.0015710.0015750.001583π^50.0023160.0023330.0023430.0014610.0014700.0014700.0012490.0012500.0012530.0014990.0015010.001512π^60.0022940.0023080.0023110.0014950.0015060.0015030.0012460.0012460.0012480.0013800.0013880.001395ρ^0.0049620.0055200.0058990.0026080.0029460.0035620.0016870.0017110.0018730.0033620.0037200.004300MSEπ^10.0053530.0053590.0053530.0023000.0023030.0023000.0016040.0016040.0016040.0053020.0053020.005302π^20.0053870.0053930.0053870.0022940.0022970.0022940.0015970.0015970.0015970.0040130.0040140.004013π^30.0053310.0053360.0053300.0022930.0022950.0022930.0016140.0016140.0016140.0032410.0032420.003241π^40.0053790.0053830.0053790.0022800.0022830.0022800.0015970.0015970.0015970.0026800.0026810.002680π^50.0053730.0053780.0053720.0022940.0022960.0022940.0016080.0016080.0016080.0023000.0023010.002300π^60.0053020.0053070.0053020.0023020.0023040.0023020.0016170.0016170.0016170.0020050.0020060.002005ρ^0.0049070.0049070.0049060.0020630.0020630.0020620.0014300.0014300.0014300.0026400.0026400.002639Table 5
Errors of the global MLEs forg=9.
IndexMLEScenario 9Scenario 10Scenario 11Scenario 12FSATSMGEMFSATSMGEMFSATSMGEMFSATSMGEMe¯¯π^1−0.000013−0.000018−0.0000140.0000080.0000080.000010−0.000043−0.000043−0.000045−0.000006−0.000006−0.000003π^2−0.000106−0.000106−0.0001030.0000470.0000470.000048−0.000009−0.000009−0.0000090.0000080.0000080.000010π^3−0.000043−0.000050−0.000048−0.000065−0.000066−0.000065−0.000029−0.000029−0.0000270.0000960.0000960.000093π^40.0000180.0000140.0000190.0000040.0000040.000003−0.000010−0.000010−0.000010−0.000052−0.000053−0.000054π^50.0000300.0000260.000032−0.000021−0.000021−0.000021−0.000024−0.000024−0.000024−0.000024−0.000024−0.000022π^6−0.000012−0.000016−0.000013−0.000021−0.000021−0.0000210.0000710.0000710.000072−0.000118−0.000118−0.000118π^7−0.000106−0.000108−0.000102−0.000036−0.000036−0.0000350.0000770.0000770.000077−0.000019−0.000020−0.000019π^8−0.000151−0.000155−0.000151−0.000005−0.000006−0.0000050.0000120.0000130.000014−0.000022−0.000022−0.000024π^9−0.000043−0.000044−0.0000430.0000580.0000570.000058−0.000008−0.000008−0.000008−0.000018−0.000018−0.000017ρ^−0.011963−0.011923−0.011885−0.005189−0.005188−0.005158−0.003624−0.003622−0.003582−0.005150−0.005138−0.005082SD ofe¯iπ^10.0023440.0023540.0023570.0015110.0015110.0015170.0013550.0013550.0013580.0024190.0024200.002428π^20.0023440.0023450.0023530.0015800.0015800.0015830.0012930.0012930.0012950.0020620.0020630.002073π^30.0022720.0022800.0022790.0015780.0015770.0015810.0012810.0012810.0012860.0017370.0017380.001745π^40.0023540.0023620.0023620.0015180.0015180.0015200.0012700.0012700.0012740.0016120.0016120.001616π^50.0023440.0023500.0023500.0014620.0014630.0014640.0012590.0012590.0012600.0014880.0014880.001491π^60.0022880.0022900.0022980.0015120.0015130.0015140.0012840.0012840.0012860.0014570.0014570.001461π^70.0023080.0023140.0023200.0015130.0015130.0015170.0012630.0012630.0012650.0013400.0013400.001343π^80.0022470.0022500.0022570.0014970.0014970.0014980.0012320.0012320.0012340.0012370.0012380.001241π^90.0024150.0024160.0024250.0015680.0015680.0015690.0012680.0012680.0012750.0012060.0012060.001208ρ^0.0032640.0036380.0038960.0015590.0015750.0017100.0013130.0013360.0015720.0016550.0017660.002223MSEπ^10.0053270.0053290.0053270.0023270.0023270.0023270.0016180.0016180.0016180.0053700.0053700.005370π^20.0053650.0053660.0053650.0023580.0023580.0023580.0016030.0016030.0016030.0040340.0040340.004034π^30.0054050.0054070.0054050.0023220.0023230.0023220.0016130.0016130.0016130.0032050.0032050.003205π^40.0053320.0053330.0053320.0023120.0023120.0023120.0016220.0016220.0016210.0027050.0027050.002705π^50.0053930.0053940.0053930.0023070.0023080.0023070.0016180.0016180.0016180.0023250.0023250.002325π^60.0053930.0053940.0053930.0023390.0023390.0023390.0016090.0016090.0016090.0020150.0020150.002015π^70.0053440.0053450.0053440.0023140.0023150.0023140.0016180.0016180.0016180.0018110.0018110.001811π^80.0053700.0053710.0053700.0023250.0023260.0023250.0015990.0015990.0015990.0016270.0016270.001627π^90.0053750.0053760.0053750.0023240.0023240.0023240.0016230.0016230.0016230.0014710.0014710.001471ρ^0.0033220.0033220.0033210.0013550.0013550.0013550.0009490.0009490.0009490.0013670.0013670.001367Table 6
Errors of the constrained MLEs forg=3.
IndexMLEScenario 1Scenario 2Scenario 3Scenario 4FSATSPGEMFSATSPGEMFSATSPGEMFSATSPGEMe¯¯π˜10.0000220.0000100.0000100.000002−0.000002−0.000002−0.000019−0.000017−0.000017−0.000007−0.000005−0.000005δ˜−0.000028−0.000089−0.000643−0.000019−0.000001−0.0002020.0000190.000018−0.001885−0.0000060.0000120.000919ρ˜−0.0089810.0241310.024078−0.0039950.0258560.025855−0.0028050.0284250.028424−0.0042020.0158390.015835SD ofe¯iπ˜10.0010570.0010570.0010400.0005360.0005440.0005400.0004300.0004350.0004330.0009660.0008040.000804δ˜0.0014920.0031360.0805220.0006690.0024520.0788930.0005540.0027880.0821190.0010030.0019080.081664ρ˜0.0016280.0352850.0350900.0010060.0338350.0338190.0007700.0367030.0366950.0010870.0220450.022035MSEπ˜10.0052270.0052640.0052640.0022180.0022340.0022340.0015640.0015750.0015750.0051680.0052180.005218δ˜0.0078360.0088370.0100140.0033340.0033910.0077150.0023490.0023490.0077940.0060980.0061430.009410ρ˜0.0096060.0104200.0104200.0042170.0056540.0056540.0029230.0047980.0047980.0044240.0048760.004876Table 7
Errors of the constrained MLEs forg=6.
IndexMLEScenario 5Scenario 6Scenario 7Scenario 8FSATSPGEMFSATSPGEMFSATSPGEMFSATSPGEMe¯¯π˜1−0.000041−0.000029−0.000029−0.000018−0.000013−0.000013−0.000011−0.000010−0.000010−0.000037−0.000030−0.000030δ˜0.0000400.0001250.0022680.0000200.0001050.0009820.000014−0.000047−0.0002350.0000400.0000520.001381ρ˜−0.0045730.0169950.016998−0.0019520.0165970.016597−0.0013990.0188780.018878−0.0025440.0109350.010935SD ofe¯iπ˜10.0009600.0008390.0008350.0005060.0004920.0004900.0004260.0004200.0004190.0009660.0007810.000781δ˜0.0010810.0023510.0408530.0005630.0017040.0392400.0004720.0019020.0411400.0009890.0014530.041265ρ˜0.0009850.0254330.0254290.0006190.0225970.0225960.0004750.0256470.0256470.0007400.0156970.015697MSEπ˜10.0051820.0052310.0052310.0022430.0022640.0022640.0015760.0015890.0015890.0051460.0052020.005202δ˜0.0062350.0062230.0060460.0026950.0026910.0034290.0018930.0018890.0030180.0056720.0056930.005703ρ˜0.0048620.0054650.0054650.0020780.0027410.0027410.0014590.0023790.0023790.0026810.0029270.002927Table 8
Errors of the constrained MLEs forg=9.
IndexMLEScenario 9Scenario 10Scenario 11Scenario 12FSATSPGEMFSATSPGEMFSATSPGEMFSATSPGEMe¯¯π˜10.000003−0.000005−0.0000050.0000090.0000070.000007−0.000003−0.000007−0.000007−0.000032−0.000010−0.000010δ˜0.000006−0.000017−0.000751−0.000015−0.000020−0.000252−0.000003−0.000036−0.0004990.0000350.0000190.001520ρ˜−0.0030090.0119710.011972−0.0013540.0135960.013596−0.0009230.0136030.013604−0.0013330.0062050.006205SD ofe¯iπ˜10.0009690.0008270.0008250.0005110.0004910.0004900.0004240.0004120.0004120.0010140.0008070.000807δ˜0.0010220.0017150.0268040.0005460.0013680.0271700.0004450.0012380.0268980.0010250.0011140.027407ρ˜0.0008100.0177770.0177750.0004620.0180450.0180450.0003840.0172870.0172870.0004720.0093520.009352MSEπ˜10.0051900.0052430.0052430.0022120.0022350.0022350.0015470.0015620.0015620.0051250.0051840.005184δ˜0.0058470.0058600.0053830.0024910.0024960.0027240.0017440.0017480.0021140.0053870.0054280.005078ρ˜0.0032050.0035120.0035120.0014090.0018520.0018520.0009770.0014180.0014180.0013990.0014890.001489Tables9 and 10 show the number of failures for these algorithms to converge when ɛ=1×10−5,1×10−4 within k=100,200 iterations. Fisher scoring algorithm has lower failure rate for convergence when calculating global MLEs. GEM algorithm for global MLEs is sensitive to ϵ and k. Both reducing ϵ and adding k can markedly improve the convergence possibility of GEM algorithm for global MLEs. Fisher scoring and GEM algorithms for constrained MLEs hardly fail to converge within 100 iterations. As g increases, the number of failures becomes small.Table 9
Failures of algorithms for global MLEs converging withink iterations.
Algorithmϵkg=3g=6g=9Scenario 1Scenario 2Scenario 3Scenario 4Scenario 5Scenario 6Scenario 7Scenario 8Scenario 9Scenario 10Scenario 11Scenario 12FSA1×10−51001000000000001×10−52001000000000001×10−4100100000000000TSM1×10−510024331713508310611612595453483256112014336410389607298683160321×10−520024331713508310611612595453483256112014336410389607298683160321×10−41002363361332701040181240325266225070199153580937998713967425896GEM1×10−510023356419944218253318256794508850986570386451498113777131859224371×10−5200422083503833405350592128619705164261915412860125751142299411×10−4100616118469410010000Table 10
Failures of algorithms for constrained MLEs converging withink iterations.
Algorithmϵkg=3g=6g=9Scenario 1Scenario 2Scenario 3Scenario 4Scenario 5Scenario 6Scenario 7Scenario 8Scenario 9Scenario 10Scenario 11Scenario 12FSA1×10−51003000000000001×10−52001000000000001×10−4100100000000000TSP1×10−5100133968540653993922555061923128152476953143415364105504041×10−5200133968540653993922555061923128152476953143415364105504041×10−410013396854065399392255506192312815247695314341536410550404GEM1×10−51000000000000001×10−52000000000000001×10−4100000000000000


### 4.1.2. Number of Iteration

Since the algorithms for global MLEs have identical accuracy in terms of MSEs, we further compare their efficiency by calculating the number of iteration and time cost when the convergence accuracy isϵ=1×10−5. The average number of iteration is recorded for every parameter setting in Figure 1.Figure 1
Number of iteration of algorithms for global MLEs from 1000 parameter settings.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)As shown in Figure1, GEM algorithm takes more iterations to converge than Fisher scoring algorithm and two-step method. The number of iterations is most intensive between 20 and 30 and the upper limit increases when Nii=1,2,…,g is bigger. For Fisher scoring algorithm, the number of iterations has slight superiority over two-step method.


### 4.1.3. Time Cost

The time required to achieve convergence can also reflect the convergence rate of an algorithm. Figure2 shows that GEM algorithm has the worst performance in terms of time cost. Fisher scoring algorithm always takes less time, especially when g=6,9. As Nii=1,2,…,g increases, the time distribution is more clustered for Fisher scoring algorithm and two-step method. When g is larger, difference of time cost is bigger between the two algorithms.Figure 2
Time cost of algorithms for global MLEs from 1000 parameter settings.
(a)(b)(c)Based on the above results, Fisher scoring algorithm, two-step method, and GEM algorithm can produce the same global MLEs when the number of iterations is large enough. However, Fisher scoring algorithm for global MLEs has a better convergence rate in terms of number of iteration and time cost. Constrained MLEs from Fisher scoring algorithm has smaller MSEs. Therefore, it is advisable to choose Fisher scoring algorithms to obtain the global and constrained MLEs required by the homogeneity test in this paper.


## 4.1.1. Evaluation of MLEs

Lete¯i be averaged error among estimators and true values of parameters from the 10000 random samples under the i-th parameter setting. Denote e¯¯=1/1000∑i=11000e¯i in each scenario, which is the total averaged error for 1000 parameter settings. The dispersion of e¯i in each scenario can be reflected by standard deviation (SD) value of e¯i. The MSEs of the global MLEs can be evaluated by the differences from the true parameters to the corresponding estimated values under various parameter settings. The global MLEs can be calculated based on Fisher scoring algorithm, two-step method, and GEM algorithm. The constrained MLEs are obtained by Fisher scoring algorithm, two-stage procedure, and GEM algorithm. The random samples for the former are generated under Ha and the latter are generated under H0. The convergence accuracy ϵ is defined by the differences from two close iterations and fixed as 1×10−5. The MSEs of the three algorithms for global MLEs have no significant difference as shown in Tables 3–5 . That is to say, the global MLEs are identical by these algorithms. In Tables 6–8 , the values of e¯¯, SD, and MSEs from Fisher scoring algorithm are usually smaller than other two algorithms for constrained MLEs. So, Fisher scoring algorithm has higher accuracy for constrained MLEs. All MSEs become smaller and close to each other when sample size increases. Algorithms for MLEs have better MSEs in balanced designs than unbalanced designs.Table 3
Errors of the global MLEs forg=3.
IndexMLEScenario 1Scenario 2Scenario 3Scenario 4FSATSMGEMFSATSMGEMFSATSMGEMFSATSMGEMe¯¯π^1−0.000071−0.000109−0.000076−0.000037−0.000125−0.000036−0.000058−0.000078−0.000064−0.000544−0.001051−0.000616π^2−0.000030−0.000043−0.000026−0.000147−0.000243−0.000147−0.000030−0.000049−0.000030−0.000081−0.000584−0.000116π^30.0000250.0000110.0000280.000036−0.0000690.000036−0.000047−0.000061−0.000054−0.000373−0.000824−0.000385ρ^−0.007755−0.006677−0.007174−0.003747−0.003360−0.003342−0.002615−0.002376−0.002213−0.046909−0.040828−0.040980SD ofe¯iπ^10.0024030.0024790.0024540.0014920.0016150.0015230.0012540.0012770.0012810.0024010.0024660.002459π^20.0023470.0024130.0024020.0015420.0016800.0015700.0012510.0012850.0012740.0017760.0018420.001825π^30.0023620.0024810.0024260.0014790.0017010.0015090.0012550.0012740.0012840.0014800.0015600.001515ρ^0.0119820.0158990.0138700.0062040.0079190.0078240.0047850.0059110.0064770.0086150.0109310.010636MSEπ^10.0054280.0054380.0054280.0022900.0023190.0022890.0015990.0016050.0015990.0053270.0053450.005327π^20.0053320.0053390.0053320.0022860.0023210.0022860.0016090.0016160.0016090.0031920.0032050.003192π^30.0053320.0053380.0053310.0022780.0023140.0022780.0016020.0016080.0016020.0022890.0023030.002289ρ^0.0093170.0093170.0093140.0040990.0040980.0040970.0028360.0028360.0028350.0057160.0057150.005714Table 4
Errors of the global MLEs forg=6.
IndexMLEScenario 5Scenario 6Scenario 7Scenario 8FSATSMGEMFSATSMGEMFSATSMGEMFSATSMGEMe¯¯π^10.000010−0.0000040.0000110.0000120.0000080.000016−0.000064−0.000064−0.000064−0.000022−0.000026−0.000029π^20.000009−0.0000070.000009−0.000013−0.000021−0.0000130.0000670.0000670.000066−0.000027−0.000029−0.000027π^30.000006−0.0000090.000007−0.000030−0.000036−0.000030−0.000038−0.000039−0.000039−0.000096−0.000097−0.000095π^40.0000890.0000800.000092−0.000068−0.000075−0.0000700.0000770.0000770.0000770.0000470.0000450.000047π^5−0.000121−0.000134−0.0001200.0000560.0000480.0000560.0000510.0000510.000051−0.000010−0.000012−0.000011π^60.0000620.0000480.000063−0.000011−0.000019−0.0000150.0000550.0000550.000056−0.000071−0.000075−0.000073ρ^−0.011639−0.011557−0.011491−0.005095−0.005052−0.004959−0.003580−0.003578−0.003538−0.006436−0.006387−0.006288SD ofe¯iπ^10.0023780.0023950.0023950.0015390.0015430.0015540.0012560.0012560.0012580.0022980.0023020.002313π^20.0023600.0023780.0023750.0014900.0014990.0015000.0013170.0013180.0013210.0019840.0019910.002004π^30.0023530.0023670.0023760.0015220.0015290.0015380.0012650.0012660.0012700.0017570.0017610.001777π^40.0024580.0024700.0024850.0014730.0014780.0014850.0012790.0012790.0012820.0015710.0015750.001583π^50.0023160.0023330.0023430.0014610.0014700.0014700.0012490.0012500.0012530.0014990.0015010.001512π^60.0022940.0023080.0023110.0014950.0015060.0015030.0012460.0012460.0012480.0013800.0013880.001395ρ^0.0049620.0055200.0058990.0026080.0029460.0035620.0016870.0017110.0018730.0033620.0037200.004300MSEπ^10.0053530.0053590.0053530.0023000.0023030.0023000.0016040.0016040.0016040.0053020.0053020.005302π^20.0053870.0053930.0053870.0022940.0022970.0022940.0015970.0015970.0015970.0040130.0040140.004013π^30.0053310.0053360.0053300.0022930.0022950.0022930.0016140.0016140.0016140.0032410.0032420.003241π^40.0053790.0053830.0053790.0022800.0022830.0022800.0015970.0015970.0015970.0026800.0026810.002680π^50.0053730.0053780.0053720.0022940.0022960.0022940.0016080.0016080.0016080.0023000.0023010.002300π^60.0053020.0053070.0053020.0023020.0023040.0023020.0016170.0016170.0016170.0020050.0020060.002005ρ^0.0049070.0049070.0049060.0020630.0020630.0020620.0014300.0014300.0014300.0026400.0026400.002639Table 5
Errors of the global MLEs forg=9.
IndexMLEScenario 9Scenario 10Scenario 11Scenario 12FSATSMGEMFSATSMGEMFSATSMGEMFSATSMGEMe¯¯π^1−0.000013−0.000018−0.0000140.0000080.0000080.000010−0.000043−0.000043−0.000045−0.000006−0.000006−0.000003π^2−0.000106−0.000106−0.0001030.0000470.0000470.000048−0.000009−0.000009−0.0000090.0000080.0000080.000010π^3−0.000043−0.000050−0.000048−0.000065−0.000066−0.000065−0.000029−0.000029−0.0000270.0000960.0000960.000093π^40.0000180.0000140.0000190.0000040.0000040.000003−0.000010−0.000010−0.000010−0.000052−0.000053−0.000054π^50.0000300.0000260.000032−0.000021−0.000021−0.000021−0.000024−0.000024−0.000024−0.000024−0.000024−0.000022π^6−0.000012−0.000016−0.000013−0.000021−0.000021−0.0000210.0000710.0000710.000072−0.000118−0.000118−0.000118π^7−0.000106−0.000108−0.000102−0.000036−0.000036−0.0000350.0000770.0000770.000077−0.000019−0.000020−0.000019π^8−0.000151−0.000155−0.000151−0.000005−0.000006−0.0000050.0000120.0000130.000014−0.000022−0.000022−0.000024π^9−0.000043−0.000044−0.0000430.0000580.0000570.000058−0.000008−0.000008−0.000008−0.000018−0.000018−0.000017ρ^−0.011963−0.011923−0.011885−0.005189−0.005188−0.005158−0.003624−0.003622−0.003582−0.005150−0.005138−0.005082SD ofe¯iπ^10.0023440.0023540.0023570.0015110.0015110.0015170.0013550.0013550.0013580.0024190.0024200.002428π^20.0023440.0023450.0023530.0015800.0015800.0015830.0012930.0012930.0012950.0020620.0020630.002073π^30.0022720.0022800.0022790.0015780.0015770.0015810.0012810.0012810.0012860.0017370.0017380.001745π^40.0023540.0023620.0023620.0015180.0015180.0015200.0012700.0012700.0012740.0016120.0016120.001616π^50.0023440.0023500.0023500.0014620.0014630.0014640.0012590.0012590.0012600.0014880.0014880.001491π^60.0022880.0022900.0022980.0015120.0015130.0015140.0012840.0012840.0012860.0014570.0014570.001461π^70.0023080.0023140.0023200.0015130.0015130.0015170.0012630.0012630.0012650.0013400.0013400.001343π^80.0022470.0022500.0022570.0014970.0014970.0014980.0012320.0012320.0012340.0012370.0012380.001241π^90.0024150.0024160.0024250.0015680.0015680.0015690.0012680.0012680.0012750.0012060.0012060.001208ρ^0.0032640.0036380.0038960.0015590.0015750.0017100.0013130.0013360.0015720.0016550.0017660.002223MSEπ^10.0053270.0053290.0053270.0023270.0023270.0023270.0016180.0016180.0016180.0053700.0053700.005370π^20.0053650.0053660.0053650.0023580.0023580.0023580.0016030.0016030.0016030.0040340.0040340.004034π^30.0054050.0054070.0054050.0023220.0023230.0023220.0016130.0016130.0016130.0032050.0032050.003205π^40.0053320.0053330.0053320.0023120.0023120.0023120.0016220.0016220.0016210.0027050.0027050.002705π^50.0053930.0053940.0053930.0023070.0023080.0023070.0016180.0016180.0016180.0023250.0023250.002325π^60.0053930.0053940.0053930.0023390.0023390.0023390.0016090.0016090.0016090.0020150.0020150.002015π^70.0053440.0053450.0053440.0023140.0023150.0023140.0016180.0016180.0016180.0018110.0018110.001811π^80.0053700.0053710.0053700.0023250.0023260.0023250.0015990.0015990.0015990.0016270.0016270.001627π^90.0053750.0053760.0053750.0023240.0023240.0023240.0016230.0016230.0016230.0014710.0014710.001471ρ^0.0033220.0033220.0033210.0013550.0013550.0013550.0009490.0009490.0009490.0013670.0013670.001367Table 6
Errors of the constrained MLEs forg=3.
IndexMLEScenario 1Scenario 2Scenario 3Scenario 4FSATSPGEMFSATSPGEMFSATSPGEMFSATSPGEMe¯¯π˜10.0000220.0000100.0000100.000002−0.000002−0.000002−0.000019−0.000017−0.000017−0.000007−0.000005−0.000005δ˜−0.000028−0.000089−0.000643−0.000019−0.000001−0.0002020.0000190.000018−0.001885−0.0000060.0000120.000919ρ˜−0.0089810.0241310.024078−0.0039950.0258560.025855−0.0028050.0284250.028424−0.0042020.0158390.015835SD ofe¯iπ˜10.0010570.0010570.0010400.0005360.0005440.0005400.0004300.0004350.0004330.0009660.0008040.000804δ˜0.0014920.0031360.0805220.0006690.0024520.0788930.0005540.0027880.0821190.0010030.0019080.081664ρ˜0.0016280.0352850.0350900.0010060.0338350.0338190.0007700.0367030.0366950.0010870.0220450.022035MSEπ˜10.0052270.0052640.0052640.0022180.0022340.0022340.0015640.0015750.0015750.0051680.0052180.005218δ˜0.0078360.0088370.0100140.0033340.0033910.0077150.0023490.0023490.0077940.0060980.0061430.009410ρ˜0.0096060.0104200.0104200.0042170.0056540.0056540.0029230.0047980.0047980.0044240.0048760.004876Table 7
Errors of the constrained MLEs forg=6.
IndexMLEScenario 5Scenario 6Scenario 7Scenario 8FSATSPGEMFSATSPGEMFSATSPGEMFSATSPGEMe¯¯π˜1−0.000041−0.000029−0.000029−0.000018−0.000013−0.000013−0.000011−0.000010−0.000010−0.000037−0.000030−0.000030δ˜0.0000400.0001250.0022680.0000200.0001050.0009820.000014−0.000047−0.0002350.0000400.0000520.001381ρ˜−0.0045730.0169950.016998−0.0019520.0165970.016597−0.0013990.0188780.018878−0.0025440.0109350.010935SD ofe¯iπ˜10.0009600.0008390.0008350.0005060.0004920.0004900.0004260.0004200.0004190.0009660.0007810.000781δ˜0.0010810.0023510.0408530.0005630.0017040.0392400.0004720.0019020.0411400.0009890.0014530.041265ρ˜0.0009850.0254330.0254290.0006190.0225970.0225960.0004750.0256470.0256470.0007400.0156970.015697MSEπ˜10.0051820.0052310.0052310.0022430.0022640.0022640.0015760.0015890.0015890.0051460.0052020.005202δ˜0.0062350.0062230.0060460.0026950.0026910.0034290.0018930.0018890.0030180.0056720.0056930.005703ρ˜0.0048620.0054650.0054650.0020780.0027410.0027410.0014590.0023790.0023790.0026810.0029270.002927Table 8
Errors of the constrained MLEs forg=9.
IndexMLEScenario 9Scenario 10Scenario 11Scenario 12FSATSPGEMFSATSPGEMFSATSPGEMFSATSPGEMe¯¯π˜10.000003−0.000005−0.0000050.0000090.0000070.000007−0.000003−0.000007−0.000007−0.000032−0.000010−0.000010δ˜0.000006−0.000017−0.000751−0.000015−0.000020−0.000252−0.000003−0.000036−0.0004990.0000350.0000190.001520ρ˜−0.0030090.0119710.011972−0.0013540.0135960.013596−0.0009230.0136030.013604−0.0013330.0062050.006205SD ofe¯iπ˜10.0009690.0008270.0008250.0005110.0004910.0004900.0004240.0004120.0004120.0010140.0008070.000807δ˜0.0010220.0017150.0268040.0005460.0013680.0271700.0004450.0012380.0268980.0010250.0011140.027407ρ˜0.0008100.0177770.0177750.0004620.0180450.0180450.0003840.0172870.0172870.0004720.0093520.009352MSEπ˜10.0051900.0052430.0052430.0022120.0022350.0022350.0015470.0015620.0015620.0051250.0051840.005184δ˜0.0058470.0058600.0053830.0024910.0024960.0027240.0017440.0017480.0021140.0053870.0054280.005078ρ˜0.0032050.0035120.0035120.0014090.0018520.0018520.0009770.0014180.0014180.0013990.0014890.001489Tables9 and 10 show the number of failures for these algorithms to converge when ɛ=1×10−5,1×10−4 within k=100,200 iterations. Fisher scoring algorithm has lower failure rate for convergence when calculating global MLEs. GEM algorithm for global MLEs is sensitive to ϵ and k. Both reducing ϵ and adding k can markedly improve the convergence possibility of GEM algorithm for global MLEs. Fisher scoring and GEM algorithms for constrained MLEs hardly fail to converge within 100 iterations. As g increases, the number of failures becomes small.Table 9
Failures of algorithms for global MLEs converging withink iterations.
Algorithmϵkg=3g=6g=9Scenario 1Scenario 2Scenario 3Scenario 4Scenario 5Scenario 6Scenario 7Scenario 8Scenario 9Scenario 10Scenario 11Scenario 12FSA1×10−51001000000000001×10−52001000000000001×10−4100100000000000TSM1×10−510024331713508310611612595453483256112014336410389607298683160321×10−520024331713508310611612595453483256112014336410389607298683160321×10−41002363361332701040181240325266225070199153580937998713967425896GEM1×10−510023356419944218253318256794508850986570386451498113777131859224371×10−5200422083503833405350592128619705164261915412860125751142299411×10−4100616118469410010000Table 10
Failures of algorithms for constrained MLEs converging withink iterations.
Algorithmϵkg=3g=6g=9Scenario 1Scenario 2Scenario 3Scenario 4Scenario 5Scenario 6Scenario 7Scenario 8Scenario 9Scenario 10Scenario 11Scenario 12FSA1×10−51003000000000001×10−52001000000000001×10−4100100000000000TSP1×10−5100133968540653993922555061923128152476953143415364105504041×10−5200133968540653993922555061923128152476953143415364105504041×10−410013396854065399392255506192312815247695314341536410550404GEM1×10−51000000000000001×10−52000000000000001×10−4100000000000000


## 4.1.2. Number of Iteration

Since the algorithms for global MLEs have identical accuracy in terms of MSEs, we further compare their efficiency by calculating the number of iteration and time cost when the convergence accuracy isϵ=1×10−5. The average number of iteration is recorded for every parameter setting in Figure 1.Figure 1
Number of iteration of algorithms for global MLEs from 1000 parameter settings.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)As shown in Figure1, GEM algorithm takes more iterations to converge than Fisher scoring algorithm and two-step method. The number of iterations is most intensive between 20 and 30 and the upper limit increases when Nii=1,2,…,g is bigger. For Fisher scoring algorithm, the number of iterations has slight superiority over two-step method.


## 4.1.3. Time Cost

The time required to achieve convergence can also reflect the convergence rate of an algorithm. Figure2 shows that GEM algorithm has the worst performance in terms of time cost. Fisher scoring algorithm always takes less time, especially when g=6,9. As Nii=1,2,…,g increases, the time distribution is more clustered for Fisher scoring algorithm and two-step method. When g is larger, difference of time cost is bigger between the two algorithms.Figure 2
Time cost of algorithms for global MLEs from 1000 parameter settings.
(a)(b)(c)Based on the above results, Fisher scoring algorithm, two-step method, and GEM algorithm can produce the same global MLEs when the number of iterations is large enough. However, Fisher scoring algorithm for global MLEs has a better convergence rate in terms of number of iteration and time cost. Constrained MLEs from Fisher scoring algorithm has smaller MSEs. Therefore, it is advisable to choose Fisher scoring algorithms to obtain the global and constrained MLEs required by the homogeneity test in this paper.


## 4.2. Evaluation of Test Statistics

Since Fisher scoring algorithm has higher accuracy for global MLEs and is more efficient for constrained MLEs, we construct test statisticsTL,TSC, and TW based on global and constrained MLEs from the algorithm. The empirical TIEs and power of three test statistics can be demonstrated by Monte Carlo simulations. All tests are conducted at the significance level of α=0.05.


### 4.2.1. Type I Error Rate

TIE is the probability of rejecting a null hypothesis that is true. According to Tang et al. [2], a test is robust if its empirical TIE is between 0.04 and 0.06; liberal if it is greater than 0.06; and conservative if it is less than 0.04 when the significance level is α=0.05. Power is the probability of correctly rejecting a null hypothesis when it is false in a statistical test. A good test should not only be robust but also make power as high as possible.1,000 random parameter settings involvingπi and ρi=1,2,…,g are generated under H0, where g and N are shown in Table 2. 10,000 samples are randomly produced for every parameter setting. The empirical TIEs can be computed by dividing the number of times of rejecting H0 with 10,000. Figure 3 further reflects the comprehensive performance of TL,TSC, and TW. The empirical TIEs of TSC are smaller and close to 0.05, which means score test is more robust. Wald-type test statistic tends to be liberal when Nii=1,2,…,g is small (i.e., Ni=30,70) and unbalanced. Wald-type test statistic is more liberal when g is bigger. Empirical TIEs of three tests are closer to 0.05 as Ni increases. Balanced designs make test statistics have better performance than unbalanced designs.Figure 3
Empirical TIEs of test statistics.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)


### 4.2.2. Power

LetN1=N2=⋯=Ng=m, g=3,6,9, π1=0.2, ρ=0.2,0.5,0.8, δp=0.1,0.2,0.3, πi=π1+δp when i=2k, and πi=π1 when i=2k+1k=1,2,…. 10,000 samples will be randomly generated under every parameter setting when m changes from 10 to 160 at intervals of 10. Empirical power can be computed by dividing the number of times of rejecting H0 with 10,000. Figure 4 reflects how empirical power changes as m changes. TW has higher power and TL and TSC have close power. Empirical power of three test statistics is close to each other as m increases. Increasing the number of groups can produce higher power. It means that all test statistics work well under multigroup cases. The moderately and highly relevant data (i.e., ρ=0.5 and ρ = 0.8) have lower power than the mildly relevant data (i.e., ρ = 0.2).Figure 4
Empirical power of test statistics.
(a)(b)(c)According to the above results, score test is more robust and has satisfactory power. Thus, it should be recommended for the homogeneity test about many-to-one comparison of risk differences.


## 4.2.1. Type I Error Rate

TIE is the probability of rejecting a null hypothesis that is true. According to Tang et al. [2], a test is robust if its empirical TIE is between 0.04 and 0.06; liberal if it is greater than 0.06; and conservative if it is less than 0.04 when the significance level is α=0.05. Power is the probability of correctly rejecting a null hypothesis when it is false in a statistical test. A good test should not only be robust but also make power as high as possible.1,000 random parameter settings involvingπi and ρi=1,2,…,g are generated under H0, where g and N are shown in Table 2. 10,000 samples are randomly produced for every parameter setting. The empirical TIEs can be computed by dividing the number of times of rejecting H0 with 10,000. Figure 3 further reflects the comprehensive performance of TL,TSC, and TW. The empirical TIEs of TSC are smaller and close to 0.05, which means score test is more robust. Wald-type test statistic tends to be liberal when Nii=1,2,…,g is small (i.e., Ni=30,70) and unbalanced. Wald-type test statistic is more liberal when g is bigger. Empirical TIEs of three tests are closer to 0.05 as Ni increases. Balanced designs make test statistics have better performance than unbalanced designs.Figure 3
Empirical TIEs of test statistics.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)


## 4.2.2. Power

LetN1=N2=⋯=Ng=m, g=3,6,9, π1=0.2, ρ=0.2,0.5,0.8, δp=0.1,0.2,0.3, πi=π1+δp when i=2k, and πi=π1 when i=2k+1k=1,2,…. 10,000 samples will be randomly generated under every parameter setting when m changes from 10 to 160 at intervals of 10. Empirical power can be computed by dividing the number of times of rejecting H0 with 10,000. Figure 4 reflects how empirical power changes as m changes. TW has higher power and TL and TSC have close power. Empirical power of three test statistics is close to each other as m increases. Increasing the number of groups can produce higher power. It means that all test statistics work well under multigroup cases. The moderately and highly relevant data (i.e., ρ=0.5 and ρ = 0.8) have lower power than the mildly relevant data (i.e., ρ = 0.2).Figure 4
Empirical power of test statistics.
(a)(b)(c)According to the above results, score test is more robust and has satisfactory power. Thus, it should be recommended for the homogeneity test about many-to-one comparison of risk differences.


## 5. Real Examples

In this section, we introduce two real examples to illustrate the aforementioned methods. The first example was presented by Rosner [1] to illustrate the new proposed methods. As shown in Table 11, 216 persons aged 20–39 with retinitis pigmentosa (RP) were classified into four genetic types including autosomal dominant RP (DOM), autosomal recessive RP (AR), sex-linked RP (SL), and isolate RP (ISO). The results from the four groups were assessed by the Snellen visual acuity (VA). An eye was considered affected if VA was 20/50 or worse and normal if VA was 20/40 or better.Table 11
The number of affected eyes in each genetic type.
Number of affected eyeslGenetic typeDOMARSLISO01573671652242791457Through Fisher scoring algorithm, two-step method, and GEM algorithm, global MLEsπ^1,π^2, π^3,π^4, and ρ^ are listed in Table 12. We observe that three algorithms can produce the same results. According to global MLEs of proportions in AR, SL, and ISO groups, estimated risk differences can be calculated as δ^2=0.5455−0.3625=0.1830, δ^3=0.7926−0.3625=0.4301, and δ^4=0.4658−0.3625=0.1033. Constrained MLEs π˜1,δ˜,ρ˜ vary slightly in Fisher scoring algorithm, two-stage procedure, and GEM algorithm as shown in Table 12.Table 12
Global and constrained MLEs.
Algorithmπ^1π^2π^3π^4ρ^π˜1δ˜ρ˜FSA0.36250.54550.79260.46580.64160.36310.14490.6537TSM/TSP0.36250.54550.79260.46580.64160.35710.15090.6572GEM0.36250.54550.79260.46580.64160.33340.17790.6558Through Fisher scoring algorithms, the values ofTL,TSC, and TW are 9.4332, 8.7076, and 13.0135, respectively, bigger than 95% percentile of the chi-square distribution with 2 degrees of freedom. Three tests all have p values smaller than 0.05 (see Table 13). Thus, there is strong evidence to reject the null hypothesis H0: δ2=δ3=δ4≜δ at a significance level of 0.05.Table 13
Test statistics’ values andp values.
StatisticTLTSCTWValue9.43328.707613.0135p value0.00890.01290.0015The second example is a cross-sectional population-based study about avoidable blindness in Iran by Rajavi et al. [15]. Nearly 3000 persons were examined, where blindness is assessed for seven age groups presented in Table 14.Table 14
Prevalence of avoidable blindness from a sample population in Iran.
BlindnessAge group (unit: yrs)50−5455−5960−6465−6970−7475−7980+0964541469257242127104123171816323029228453910Table15 provides global MLEs and constrained MLEs. The common risk differences δ from three algorithms are estimated to be 0.0419, 0.0433, and 0.0462. Table 16 shows that p values of test statistics TL,TSC, and TW are 1.4474×10−21, 8.9818×10−26, and 3.4369×10−14. So we have enough evidence to reject H0: δ2=δ3=δ4=δ5=δ6=δ7≜δ at a significance level of 0.05.Table 15
Global and constrained MLEs.
Algorithmπ^1π^2π^3π^4π^5π^6π^7ρ^π˜1δ˜ρ˜FSA0.01460.02620.02620.04490.07360.14460.17060.27010.01480.04190.2942TSM/TSP0.01460.02620.02620.04490.07360.14460.17060.27010.01370.04330.3024GEM0.01470.02620.02630.04500.07370.14480.17080.27390.01110.04620.3036Table 16
Test statistics’ values andp values.
StatisticTLTSCTWValue107.4057127.281472.2784p value1.4474×10−218.9818×10−263.4369×10−14


## 6. Conclusion

This paper mainly studies the many-to-one comparison of risk differences for correlated binary data. Three algorithms, including Fisher scoring algorithm, two-step method, and GEM algorithm, are introduced for global MLEs. Constrained MLEs are calculated by Fisher scoring algorithm, two-stage procedure, and GEM algorithm. Simulations show that Fisher scoring algorithm, two-step method, and GEM algorithm can produce global MLEs with close MSEs. However, Fisher scoring algorithm for global MLEs behaves better in terms of convergence rate. Fisher scoring algorithm can obtain constrained MLEs with lower MSEs. Thus, it is recommended for getting constrained MLEs. Likelihood ratio, score, and Wald-type tests are constructed based on global and constrained MLEs from Fisher scoring algorithms. The results of simulations show that score test statistic is robust and has high empirical power.In particular, Fisher scoring algorithm may not work well with high dimensional Fisher information matrix when the number of groups is massive. Then, the two-step method or two-stage procedure can be adopted to solve the high dimensional problems. As the modified Fisher scoring algorithm or Newton–Raphson algorithm, they only need to calculate a small part of parameters at every step to reduce dimension. GEM algorithm is also an alternative approach regardless of time cost. The biggest obstacle for GEM algorithm may occur inM step when the local optimal solution does not exist or cannot be found. Furthermore, the convergent performance of GEM algorithm is sensitive for the convergence accuracy and the iteration number.In this article, we investigate the homogeneity test of many-to-one risk differences in multigroup design. If samples can be stratified by some control variables (e.g., age, gender, etc.), the treatment-by-stratum interaction should be considered. Thus, it is of great significance to develop many-to-one comparison in stratified correlated binary data. Other indexes such as risk ratio can be adopted for evaluating differences of proportions instead of risk difference.



---

*Source: 6685951-2021-04-08.xml*