Estimation of Some Epidemiological Measures of Association in Multiple Comparative Trials

Estimation of Some Epidemiological Measures of Association in Multiple Comparative Trials Soumik Banerjee*and Krishna Saha† ‡ Abstract Multiple comparative trials with binary outcomes are commonly used in biomedical research and other disciplines for estimating epidemiological measures by combining information from multiple comparative trials. The epidemiological measures are commonly estimated as a weighted average of summary statistics based on the 2 × 2 table data from each trial. Three of the most important epidemiological measures are frequently used: the risk ratio (RR), odds ratio (OR), and risk difference (RD). The RD/RR are preferable due to a more meaningful and interpretable treatment measure for binary outcome. The estimation procedures for estimating the overall RD/RR in multiple comparative trials with binary outcomes are very challenging and difficult, especially when the number of patients in a single trial is small and when the number of events is zero for some trials. Considering the above situations, we develop some efficient estimation procedures for estimating the overall RD/RR in multiple comparative trials with binary outcomes. We illustrate those estimation procedures by analyzing two real-life data sets obtained from multiple comparative trials in biomedical research. Key Words: Clinical Trials, Binary Data, Likelihood Ratio Test, Risk Difference,Test of Homogeneity, Estimation. 1. Introduction Multiple comparative trials with binary outcomes (see, for example, Beitler and Landis, 1985; Cooper et al., 1993; Chu et al., 2010) are commonly used in biomedical research and other disciplines for estimating the epidemiological measures by combining the information from multiple comparative trials. In multiple comparative trials with binary outcome, the epidemiological measures can be estimated as a weighted average of summary statistics based on the 2 × 2 table data from each trial. Three of the most important epidemiological measures are frequently used: the risk ratio, odds ratio, and risk difference. The risk ratio and odds ratio are both multiplicative measures of effect, whereas the risk difference is an absolute measure of effect. The issue of deciding which effect measure to use in a particular application basically depends on how the study was designed. One would prefer to use the risk difference/ratio over the odds ratio due to a more meaningful and interpretable treatment measure for binary outcome. However, the estimation procedures for estimating the overall risk difference/ratio in multiple comparative trials with binary outcome are very challenging and difficult, especially when the number of patients in a single trial is small (see, for example, Cooper et al., 1993) and when the number of events is zero for some trials (see, for example, Beitler and Landis, 1985). *Department of Quantitative Sciences, Canisius University, Science Hall, Buffalo, NY 14208, USA †Department of Mathematical Sciences, Central Connecticut State University, 1615 Stanley Street, New Britain, CT 06050, USA ‡Correspondence However, the estimation procedures for estimating the overall risk difference/ratio in multiple comparative trials with binary outcome are very challenging and difficult, especially when the number of patients in a single trial is small (see, for example, Cooper et al., 1993) and when the number of events is zero for some trials (see, for example, Beitler and Landis, 1985). For example, data presented by Cooper et al. (1993) from a Cancer and Leukemia Group B randomized trial clearly showed that the sample sizes for both treatment groups are very small ranging from 2 to 12 and data presented by Beitler and Landis (1985) also showed that for clinic 5 and clinic 6, the number of events for a control group is zero. Another important issue is that binary outcomes from the same clinic or hospital may not be independent, and estimating the overall risk difference/ratio without counting the correlation between binary outcomes within the hospital may often lead to misleading conclusions. The choice of estimation methods for the odds ratio has been thoroughly studied and compared, but relatively little attention has been paid to developing efficient methods for the estimation of the risk difference/ratio. Taking into account all situations discussed above, we focus on developing estimation procedures for estimating the overall risk difference (RD) in multiple comparative trials with binary outcome. In this proceeding, we have discussed various methods of estimating RD under null and alternative hypothesis (Guo, 2022) and also the procedures of testing the null hypothesis (equality of the RDs across multiple centers) developed in Banerjee et al. (2025). We provided the simulation results derived in Banerjee et al. (2025) and the application of the test procedures on two real life datasets. Subsequently, we have also discussed the estimation of RD using a couple of methods (Guo, 2022). 1.1 Basic notations and definition Suppose, we wish to compare two treatments using the data combined for k centers. For each of k centers, we randomly assign patients to one of the two treatment groups. Let πij, i = 1,2, ..., k;j= 1,2be the corresponding cell probabilities of the outcome of interest. Then, ˆπij =Xij/nij is the sample estimate of this probability, where Xij and nij are the number of successes and the number of subjects in the ith center and the jth treatment group, respectively. The risk difference δiin the ith center is simply πi1−πi2. We are interested in testing H0:δ1=δ2=..... =δk=δ, (1) where δis a true common risk difference between the new and the standard treatments. 2. Estimation Methods Under H0 One simple meta-analysis approach is to assume a common RD across all studies, that is, the Common effect (CE) model, which is also commonly referred to as the Fixed effect (FE) model or equal-effect model. Let dibe the observed risk difference for the ith center. The estimate of the risk difference for center iis di= ˆπi1−ˆπi2=Xi1 ni1 −Xi2 ni2 .(2) Under H0, one can model for dias di=δH0+σiϵi, where δH0is a common risk difference across all centers and ϵiis the random variation which can follow as ϵi∼N(0,1), and the σidescribes the variation. Note that ˆσ2 i=dvar(ϵi) = dvar(di) = 1 n3 1i X1i(n1i−X1i) + 1 n3 2i X2i(n2i−X2i). 2.1 Weighted average method The common risk difference δH0can be estimated as a weighted average of individual risk difference diof each center, where the weight given to each center is the inverse of the variance of dias ˆ δH0=Pk i=1 widi Pk i=1 wi with wi= 1/dvar(di). 2.2 Mantel–Haenszel (MH) method The weights based on the Mantel–Haenszel (MH) method which uses sample sizes to define weights as w∗ i=n1in2i n1i+n2i . Similarly, the common risk difference δcan be estimated based on the above weights as ˆ δ∗ H0=Pk i=1 w∗ idi Pk i=1 w∗ i with w∗ i=n1in2i n1i+n2i . 3. Estimation Methods Under Ha: Random Effect Model If the null hypothesis is rejected, then the effects are different across centers due to their heterogeneous populations. To account for the heterogeneity of center-level true effects, the RE (random effect) model uses random effects to represent the between-center heterogeneity. Once rejecting H0, one needs to estimate the overall risk difference under Ha. Under Ha, one can model for dias di=δHa+ui+ϵi, where δHais the overall risk difference and ϵiis the random variation which can follow as ϵi∼N(0, σ2 i)and uiis the random effect due to between-center heterogeneity which can follow as ui∼N(0, σ2 u).The variance of the observed risk difference diis given by var(di) = var(ϵi) + var(ui) = σ2 i+σ2 u, and the estimated variance of diis given by dvar(di) = ˆσ2 i+ ˆσ2 u, where ˆσ2 iis the estimate of σ2 idiscussed earlier and ˆσ2 uis estimate of σ2 uwhich will be discussed in section 3. Once we obtain the estimated variance of diunder Hausing the estimates of σ2 iand σ2 u, we can similarly obtain the overall risk difference δHa, as a weighted average of individual risk difference diof each center, given by ˆ δHa=Pk i=1 wHa idi Pk i=1 wHa i , where wHa i= 1/(ˆσ2 i+ ˆσ2 u)is the weight given to each center as the inverse of the variance of diunder Ha. 3.1 Methods Under Ha: : Estimate of σ2 u The estimate of the between-center variance σ2 ucan be obtained based on the following methods (Guo et al. 2022) discussed as follows: • DL (DerSimonian and Laird, 1986) estimator • HO (Hedges and Olkin, 1985) estimator • PM (Paule and Mandel, 1982) estimator • ML (Maximum Likelihood) estimator • Restricted ML estimator • HS (Hunter and Schmidt, 2015) estimator • SJ (Sidik-Jonkman (2005)) estimator • EB (Empirical Bayes) estimator DL Estimator: The DL estimator is the most widely used method to estimate σ2 uin the current literature. The estimator is given by: ˆσ2 u,DL =    [Q−(k−1)]/[Pk iwi−Pk iw2 i Pk iwi], Q > k −1 0otherwise, where Q=Pk i=1 wi(di−ˆ δH0)2with wi= 1/dvar(di)is the Cochran’s Q heterogeneity statistic, kis the number of studies included in the meta-analysis. HO Estimator: The HO estimator is an unbiased estimator of σ2 ugiven by; ˆσ2 u,HE =Pk i(di−ˆ δH0)2 k−1−1 k k X i ˆσi2. PM Estimator: The PM estimator applies iterations to the generalized Q statistic. Q(r+1) P M = k X i w(r) i,P M (di−ˆ δ(r) P M ),ˆ δ(r) P M =Pk iw(r) i,P M di Pk iw(r) i,P M , w(r) i,P M = (ˆσ2 i+ ˆσ2,(r) P M )−1 The generalized Q statistic is expected to equal k−1, and Q(r+1) P M is the generalized Q statistic at the (r + 1)th iteration. The ˆ δ(r) P M and ˆσ2,(r) P M are the estimates of δ(r) P M and σ2,(r) P M at the r’th iteration. ML Estimator: The RE model is essentially a linear mixed model, which could be implemented by the ML estimation. The logarithmic likelihood function given by: logL(δML, σ2 ML|d1, d2, ..., dk) = −1 2 k X i=1 log(σ2 ML +σ2 i)−1 2 k X i=1 (di−δML)2 σ2 ML +σ2 i Taking partial derivatives to the log-likelihood function, δ2 ML, could be solved via the following equation: ˆ δ2 ML =Pk i=1 w2 i,ML[(di−ˆ δML)2−ˆσ2 i)] Pk i=1 w2 i,ML where wi,ML = (ˆ δ2 ML + ˆσ2 i)−1and ˆ δML is the estimated δML. We estimate ˆσML and δML through iterations. Restricted ML estimator: The ML estimator fails to account for the heterogeneity of the data and also,is sometimes negatively biased; the Restricted ML estimator (REML) solves these problems by a transformation. Specifically, the REML estimation modifies the log-likelihood function as logL(σ2 REML|d1, d2, d3, ...dk) = −1 2[ k X i=1 log(σ2 REML+σ2 i)+log k X i=1 1 σ2 REML +σ2 i +1 wi,REML ] where ˆ δML is the likelihood-based estimate of δ, which is identical to the obtained as the ML estimator. The weight wi,REML is (ˆσ2 REML + ˆσ2 i)−1. The estimation also requires iterations. HS Estimator: The HS estimator is an unbiased estimate of σ2 u, which is given by, ˆσ2 u,HS =Pk iwi(di−ˆ δH0)2 Pk iwi −Pk iwiˆσ2 i Pk iwi . SJ Estimator: The SJ estimator is a general between-study variance estimator that is based on the weighted residual sum of squares as follows: ˆσ2 SJ =1 k−1 k X i=1 (di−ˆ δˆv)2 ˆvi,SJ Here δˆv=Pk i=1 ˆv−1 i,SJ di/Pk i=1 ˆv−1 i,SJ ,ˆvi,SJ = ˆri,SJ + 1, and ˆri,SJ = ˆσ2 i/ˆσ2 0, where ˆσ0 is a rough estimate of τsuch that ˆσ2 0=1 kPk i=1(di−ˆ d)2with ˆ d=1 kPk i=1 di. EB estimator: The EB estimator is a multi-parameter shrinkage estimator of σu. If δk is the true effect size for study k and ˆ δH0is the effect estimate from the weighted average method, then di|δk∼N(δk, σ2 i), δk|ˆ δH0, σ2 EB ∼N(ˆ δH0, σ2 EB), i = 1,2, .., k The posterior distribution of δiis: δi|di,ˆ δH0, σ2 EB ∼N(δ∗ i, σ2 i(1 −Bi)) where δ∗ i= (1 −Bi)di+Biˆ δH0and Bi=σ2 i/(σ2 i+σ2 EB). Similar to ˆσ2 REML,ˆσ2 EB can be estimated by iterations using the following equations if the within-study variances σ2 iare homogeneous; ˆσ2 EB =Pk i=1 wi,EB [k k−1(di−ˆ δ2 H0)−ˆσ2 i] Pk i=1 Wi,EB where wi,EB = (ˆτ2 EB + ˆσ2 i)−1. 4. Test Procedures for testing the homogeneity of risk difference The test procedures for testing H0:δ1=δ2=.... =δK=δproposed by Lipsitz et al. (1998) and also developed by Banerjee et al. (2025) are discussed briefly as follows: 4.1 Non-parametric test procedures Let di= ˆπi1−ˆπi2, where the cell probabilities ˆπij =Xij/nij and Xs ij are assumed to follow binomial distribution with Xij ∼Bin(nij, πij ),i= 1, ..., k and j= 1,2. Then, ˆ δ=P(di/ˆwi)/P(1/ˆwi)where ˆwi=d V ar(di) = (ˆπi1(1 −ˆπi1))/(ni1−1) + (ˆπi2(1 −ˆπi2))/(ni2−1). As all nij → ∞, the weighted least squares test statistic QW LS =X(di−ˆ δ)2/ˆwi∼χ2 k′−1, where k′=k−qwith qbeing the number of centers in which ˆwi= 0. As all k→ ∞, they proposed the test statistic Z2 W LS = (QW LS −(k′−1))2/[2(k′−1)] ∼F1,(k′−1). To improve the normal approximation of a chi-square distribution of QW LS, Lui and Kelly (2000) proposed the following test statistic as ZLW LS =ln[QW LS /(k′−1)]/2+1/[2(k′−1)] p1/[2(k′−1)] ∼N(0,1) When all nij’s are not large enough, the performances of QW LS and Z2 W LS are expected not to be satisfactory. To improve this, Lipsitz et al. (1998) proposed the following three weighted test statistics as follows: Z2 W LS,R = (QW LS −k′)2/P[(di−ˆ δ)2/ˆwi−1]2 Z2 V= [P((di−ˆ δ)2−ˆwi)/ai]2/P[(di−ˆ δ)2−ˆwi]2/a2 i Z2 K= [P((di−ˆ δ)2−ˆwi)/((1/ni1) + (1/ni2))2]2/P[(di−ˆ δ)2−ˆwi]2/((1/ni1) + (1/ni2))4 4.2 Parametric test procedures Due to the heavy dependence on the number of clusters and the average sample size, Banerjee et al. (2025) developed two model based (parametric) test procedures. Binomial Likelihood ratio test (LRT) Assume Xij ∼Bin(nij, πij), i = 1,2, .., K;j= 1,2. Then the log-likelihood can be written as ln L=Piln ni1 πi1+Piln ni2 πi2+Pixi1ln πi1+Pixi2ln πi2+Pi(ni1−xi1) ln(1 − πi1)+Pi(ni2−xi2) ln(1 −πi2) Recall that δi=πi2−πi1=⇒πi2=δi+πi1. We can then rewrite the log-likelihood as ln L=Piln ni1 πi1+Pilnni2 πi2+Pixi1lnπi1 +Pixi2ln(δi+πi1)+Pi(ni1−xi1)ln(1 −πi1) +Pi(ni2−xi2)ln(1 −δi−πi1) Let ln L0and ln Labe the maximum log-likelihoods under H0and Ha, respectively. Then, the LRT statistic for testing H0against Hacan be obtained as LRT = 2[ln La−ln L0]∼χ2 K−1. Binomial C(α)test For the convenience of the derivation of the C(α)test statistic, we re-parameterize δiunder Haby δi=δ+αi,i= 1, ..., k−1, with αk= 0. Then testing H0:δ1=δ2=... =δk=δ reduces to testing H0:α1=... =αk−1= 0. Let α′= (α1, α2, .., αk−1)and ω′= (δ, πi1, πi2, ...., πik), with ω′being treated as nuisance parameter. Let’s denote U=δl δα |α= 0,V=δl δω |α= 0,A=E(−δ2l δαδα′|α= 0),C=E(−δ2l δαδω |α= 0), and D=E(−δ2l δωδω′|α= 0) The C(α)test (Neyman, 1959 and Saha, 2013) based on the residual S(ω) = U(ω)− β′V(ω)I, for testing H0against Ha, is given by ˆ S(ˆ A−ˆ Cˆ D−1ˆ C′)−1ˆ S′, which asymptotically follows χ2 k−1as nij → ∞. 5. Simulation Study for Test Procedures Following Lui and Kelly (2000), we conducted the simulation study to check the performance in terms of the Type I error rate of the previous seven test methods, along with our two proposed likelihood-based test statistics, (i) Binomial LRT and (ii) Binomial C(α)test methods. We follow the instructions of Lui and Kelly (2000) to generate random samples. As mentioned, the probability πi2varies between centers due to center effects, we set πi2= .10 + .70p, where pis generated from a beta distribution with mean π=α/(α+β) and variance π(1 −π)/(T+ 1) with T=α+β. We then set πi1=πi2+τ, here τ is the underlying common risk difference. We have performed simulations for τ= 0.1 and 0.2. Note that pfollows the beta distribution with T=2, and π= 0.5is equivalent to assuming that pfollows a uniform distribution on (0,1). We generate the sample size nij following the probability mass function f(nij)=0.20 for nij =n−s, where nis a given fixed constant and s= 2,1,0,−1,−2. We have taken the intraclass correlation ρ(= 1/[2(T+ 1)]) from approximately 0.05, 0.1 and 0.2; the number of centers Kequals 8, 16, 20 and 30; and the mean treatment group size per center E(nij) = nhas been considered 4,8,10,15,20,25,30,40 and 60. All results are obtained using R codes. For each configuration of the above sets of parameters, we generated 10,000 repeated samples under the null hypothesis. We calculated Type I error rates for the nine test procedures, including seven previous and two new proposed methods. Results have been provided in Table 1. From the results given in Table 1, it can be seen that the parametric test procedures, LRT and C(α), perform much better than the other non-parametric test procedures developed by Lipsitz et al. (1998) in terms of Type-I error rate. 6. Data Analysis After simulation, we have provided the detailed analysis of two separate datasets. First one is the example when the null hypothesis is accepted, whereas, for the second one, the null hypothesis of testing the homogeneity of risk difference is rejected. Example 1: Data from Multi-Center Randomized Clinical Trial This data set contains information on the distribution of favorable responses to active drug and control treatment in a multi-center randomized clinical trial This dataset contains information on the distribution of favorable responses to active drug and control treatment in a multi-center, randomized clinical trial. Data were collected from 8 clinics with 273 patients or 17 patients per clinic.Full data is available in Table 2. We first need to test the following hypothesis: H0:δ1=δ2=.... =δk=δ LRT = 8.583 with p-value = 0.263 C(α) = 13.065 with p-value = 0.075 Conclusion: do not reject H0, indicating a non-significant between-clinic heterogeneity, that is, the two methods discussed under H0will be suitable to estimate the common risk difference. The details of estimation of the common risk difference using the weighted average method is given in Table 3. Then, the common risk difference δH0based on weighted average is ˆ δH0=Pk i=1 widi Pk i=1 wi = 0.124. The details of estimation of the common risk difference using the weighted average (MH) method is given below in Table 4. Then, the common risk difference δH0based on MH method is ˆ δ∗ H0=Pk i=1 w∗ idi Pk i=1 w∗ i = 0.128. Example 2: Risk of stunting in left-behind children Fellmeth et al. (2018) presented a meta-analysis investigating whether left-behind children and adolescents had different health conditions compared with children and adolescents with non-migrant parents. We focus on stunting as a binary outcome. This research identified 16 studies without zero events, including 38,779 children (17,066 left-behind children and adolescents and 21,713 children and adolescents with non-migrant parents). Because this dataset contains no zero events, we can use it to investigate whether each method gives a consistent estimate of RDs. Full dataset is available below in Table 5. We first need to test the following hypothesis: H0:δ1=δ2=.... =δK=δ LRT = 52.81 with p-value = 0.0001 C(α) = 55.21 with p-value = 0.0001. Conclusion: H0is rejected, indicating a significant between-study heterogeneity so the DL and HO methods under Hawill be suitable to estimate the overall risk difference.Details of DL method has been provided below in Table 6. Estimating Risk Difference: DL Method Based on Table 6, we obtain Q= K X i=1 wi(di−ˆ δH0)2= 52.71. The DL estimator of σ2 uis given by ˆσ2 u,DL = [Q−(k−1)]/[ k X i wi−Pk iw2 i Pk iwi ] = 0.001. Then, the overall risk difference δHabased on DL method can be obtained as ˆ δHa=Pk i=1 wDL,idi Pk i=1 wDL,i = 0.014. Here wDL,i = 1/(ˆσ2 i+ ˆσ2 u,DL).