Minimum MSE Weighted Estimator to Make Inferences for a Common Risk Ratio across Sparse Meta-Analysis Data

Chukiat Viwatwongkasem^1*, Sutthisak Srisawad², Pichitpong Soontornpipit¹, Jutatip Sillabutra¹, Pratana Satitvipawee¹, Prasong Kitidamrongsuk¹, Hathaikan Chootrakool¹

¹Department of Biostatistics, Faculty of Public Health, Mahidol University, Bangkok, Thailand.
²Nursing Research Center, Faculty of Nursing, Mahidol University, Nakhon Pathom, Thailand.

DOI: 10.4236/ojs.2022.121004                

Abstract

The paper aims to discuss three interesting issues of statistical inferences for a common risk ratio (RR) in sparse meta-analysis data. Firstly, the conventional log-risk ratio estimator encounters a number of problems when the number of events in the experimental or control group is zero in sparse data of a 2 × 2 table. The adjusted log-risk ratio estimator with the continuity correction points  based upon the minimum Bayes risk with respect to the uniform prior density over (0, 1) and the Euclidean loss function is proposed. Secondly, the interest is to find the optimal weights of the pooled estimate  that minimize the mean square error (MSE) of  subject to the constraint on  where , , . Finally, the performance of this minimum MSE weighted estimator adjusted with various values of points  is investigated to compare with other popular estimators, such as the Mantel-Haenszel (MH) estimator and the weighted least squares (WLS) estimator (also equivalently known as the inverse-variance weighted estimator) in senses of point estimation and hypothesis testing via simulation studies. The results of estimation illustrate that regardless of the true values of RR, the MH estimator achieves the best performance with the smallest MSE when the study size is rather large  and the sample sizes within each study are small. The MSE of WLS estimator and the proposed-weight estimator adjusted by , or , or are close together and they are the best when the sample sizes are moderate to large (and) while the study size is rather small.

Keywords

Minimum MSE Weights, Adjusted Log-Risk Ratio Estimator, Sparse Meta-Analysis Data, Continuity Correction

1. Introduction

Frequently, biostatisticians would like to evaluate the effects of treatments or risk factors in terms of risk difference, relative risk (risk ratio), and/or odds ratio between two independent sample groups (e.g., treatment or control, presence or absence of a risk factor) and binary outcomes (e.g., disease or non-disease, success or failure, dead or alive) in a 2 × 2 table. Let $p_1$ be the probability of outcome in the treatment or exposed group and $p_2$ be the probability of outcome in the control or unexposed group. On comparing two independent groups, the popular effect parameters are defined by the risk difference $RD=p_1-p_2$, the relative risk $RR=p_1/p_2$, or the odds ratio $OR=\frac{p_1/\left(1-p_1\right)}{p_2/\left(1-p_2\right)}$. Obviously, all three effect sizes are related to the estimation of proportion p on each arm. It has widely been known that the conventional proportion estimator $\hat{p}=X/n$ is a good choice for estimating p in general, but may not be in sparse data. For example, in sparse data coping with the small number of events X and the small sample size n, the variance of $\hat{p}$, estimated by $\frac{\hat{p}\left(1-\hat{p}\right)}{n}$, can cause a problem with a value of 0 when $X=0$ or $X=n$. To solve the problem, a continuity correction term c is often added to each cell of each group in the 2 × 2 table, yielding to ${\hat{p}}_c=\left(X+c\right)/\left(n+2c\right)$ in each arm. Yate [1] first used the continuity correction of 0.5 in the approximation of a discrete distribution to a continuous one in 1934. And it seems that the correction value of 0.5 has been used extensively until now, for examples: Lane [2], Stijnen et al. [3], White et al. [4], Lui and Lin [5], Sankey et al. [6], Gart and Zwefel [7], Walter [8], and Cox [9] used value 0.5 adjustment for zero observations in each cell of the 2 × 2 table. Another choice of c for this class such as 0.25, 0.5 and 1 had been suggested by Li and Wang [10]; 1/6 by Turkey [11] and Sánchez-Meca and Marin-Martinez [12]; Böhning and Viwatwongkasem [13] showed that the simple adjusted estimate ${\hat{p}}_c=\left(X+c\right)/\left(n+2c\right)$ with $c=1$ performed surprisingly well with the smallest average MSE. In addition, Agresti and Caffo [14] also suggested a value of 1 to solve the zero observations; however, a correction value of 2 was recommended by MaClave and Sincich [15].

Our focus of interest is not on the simple adjusted estimate ${\hat{p}}_c=\left(X+c\right)/\left(n+2c\right)$, but on the logarithm of ${\hat{p}}_c$ instead, leading to the final interest to the logarithm of the risk ratio estimate, ${\hat{\theta }}_c=\log {\widehat{RR}}_c=\log {\hat{p}}_{c1}-\log {\hat{p}}_{c2}=\log \left(\frac{X_1+c_1}{n_1+2c_1}\right)-\log \left(\frac{X_2+c_2}{n_2+2c_2}\right)$. The reason is that we are interested in the logarithm function, because it is widely known that the risk ratio has a non-symmetric and rather right-skewed distribution about the null value of 1; consequently, the natural logarithm of the risk ratio is needed to transform to be more reliable for the normal distribution. For estimating $\log \left(\frac{p+c/n}{1+2c/n}\right)$, Pettigrew, Gart and Thomas [16] proposed the moment estimator $\log \left(\frac{X+c}{n+2c}\right)$ with the bias and the first four cumulants by means of asymptotic Taylor’s series expansion. Likewise, we re-derive and expand on the details of the solution properties for the adjusted estimator of the log-risk ratio, ${\hat{\theta }}_c=\log {\widehat{RR}}_c=\log {\hat{p}}_{c1}-\log {\hat{p}}_{c2}=\log \left(\frac{X_1+c_1}{n_1+2c_1}\right)-\log \left(\frac{X_2+c_2}{n_2+2c_2}\right)$.

Therefore, the first objective of the study is to find the optimal points $\left(c_1,c_2\right)$ based upon continuity correction of this adjusted log-risk ratio estimator under various settings, such as the smallest bias, the smallest average bias, the smallest MSE, and the smallest average MSE.

Secondly, after obtaining the optimal value of $\left(c_1,c_2\right)$ which is the best choice when considering the smallest bias and/or the smallest MSE of ${\hat{\theta }}_c$. The next focus of interest is concerned with sparse data in meta-analysis studies that combine the various risk ratios from k studies to produce a single summary risk ratio. Under a common risk ratio overall studies or homogeneity of risk ratios across k studies, the second aim of the study is to find the optimal weights ${\hat{f}}_j$ of the pooled estimate ${\hat{\theta }}_w={\displaystyle {\sum }_{j=1}^k{\hat{\theta }}_{cj}{\hat{f}}_j}$ that minimize the MSE of ${\hat{\theta }}_w$ subject to the constraint on ${\displaystyle {\sum }_{j=1}^k{\hat{f}}_j}=1$ where ${\hat{\theta }}_{cj}=\log \left({\widehat{RR}}_{cj}\right)=\log \left({\hat{p}}_{c1j}\right)-\log \left({\hat{p}}_{c2j}\right)$, ${\hat{p}}_{c1j}=\left(X_{1j}+c_1\right)/\left(n_{1j}+2c_1\right)$, ${\hat{p}}_{c2j}=\left(X_{2j}+c_2\right)/\left(n_{2j}+2c_2\right)$. Indeed, these optimal weights ${\hat{f}}_j$ actually come from the roughly similar formula to the work of Viwatwongkasem et al. [19] that was used for the risk difference in multi-center studies.

Finally, the final objective of the study of interest is about comparing the performance of well-known summary effect estimators, such as the Mantel-Haenszel (MH) estimator stated with its variance estimate by Greenland and Robin [20] and the weighted least square (WLS) estimator or equivalently known as the inverse-variance weighted estimator to the new adjusted summary effect estimator ${\hat{\theta }}_w$ with various adjustment points $\left(c_1,c_2\right)$ in terms of point estimation and hypothesis testing via simulation studies.

The rest of the paper is organized as follows. Section 2 contains the derivative methods and the results of the adjusted log-risk ratio estimator. Section 3 discusses on the methodology and the outcomes of the minimum MSE weights of the adjusted summary relative risk estimator. Section 4 states on the other well-known estimators and tests as the comparative candidates to compare the performances of making inference for a common relative risk across meta-analysis studies. Section 5 consists of the simulation plans for studying the performances in senses of the estimation, the type I errors and the power of the tests. Section 6 contains the results of Section 5. Section 7 is about the discussion and the recommendation.

2. Adjusted Estimator of the Logarithm of the Risk Ratio

In a single study ( $k=1$ ), since the conventionally popular estimator of the logarithm of risk ratio, obtained by $\hat{\theta }=\log \widehat{RR}=\log {\hat{p}}_1-\log {\hat{p}}_2=\log \left(\frac{X_1}{n_1}\right)-\log \left(\frac{X_2}{n_2}\right)$, may have a problem when data are sparse. To solve the problem, a continuity correction term $\left(c_1,c_2\right)$ is often added to each cell of the 2 × 2 table, leading to the adjusted estimate as ${\hat{\theta }}_c=\log {\widehat{RR}}_c=\log {\hat{p}}_{c1}-\log {\hat{p}}_{c2}=\log \left(\frac{X_1+c_1}{n_1+2c_1}\right)-\log \left(\frac{X_2+c_2}{n_2+2c_2}\right)$. In addition, the adjusted risk ratio estimator can be further found as ${\widehat{RR}}_c=\frac{{\hat{p}}_{c1}}{{\hat{p}}_{c2}}=\frac{\left(X_1+c_1\right)/\left(n_1+2c_1\right)}{\left(X_2+c_2\right)/\left(n_2+2c_2\right)}$. Due to the work of Pettigrew, Gart and Thomas [16], we re-derive and extend the results to get the expectation, bias, variance, and MSE of ${\hat{\theta }}_c$ in the following:

$E\left({\hat{\theta }}_c\right)=\left(\log p_1-\log p_2\right)+B\left({\hat{\theta }}_c\right)$ (1)

$\begin{array}{c}B\left({\hat{\theta }}_c\right)=\frac{2c_1-1+p_1\left\{1-2\left(2c_1\right)\right\}}{2n_1{p}_1}\\ +\frac{-6c_1^2+6{\left(2c_1{p}_1\right)}^2+12c_1{q}_1+4q_1\left(q_1-p_1\right)-9q_1^2}{12{\left(n_1{p}_1\right)}^2}\\ -[\frac{2c_2-1+p_2\left\{1-2\left(2c_2\right)\right\}}{2n_2{p}_2}\\ +\frac{-6c_2^2+6{\left(2c_2{p}_2\right)}^2+12c_2{q}_2+4q_2\left(q_2-p_2\right)-9q_2^2}{12{\left(n_2{p}_2\right)}^2}]+O\left(n^{-3}\right)\end{array}$ (2)

$V\left({\hat{\theta }}_c\right)=\frac{q_1}{n_1{p}_1}+\frac{q_1\left(1-2c_1\right)+\frac{1}{2}{q}_1^2}{{\left(n_1{p}_1\right)}^2}+\frac{q_2}{n_2{p}_2}+\frac{q_2\left(1-2c_2\right)+\frac{1}{2}{q}_2^2}{{\left(n_2{p}_2\right)}^2}+O\left(n^{-3}\right)$ (3)

$\begin{array}{c}MSE\left({\hat{\theta }}_c\right)=\frac{q_1}{n_1{p}_1}+\frac{\left(c_1-2c_1{p}_1\right)\left(c_1-2c_1{p}_1-q_1\right)+q_1\left(1-2c_1\right)+\frac{3}{4}{q}_1^2}{{\left(n_1{p}_1\right)}^2}+\frac{q_2}{n_2{p}_2}\\ +\frac{\left(c_2-2c_2{p}_2\right)\left(c_2-2c_2{p}_2-q_2\right)+q_2\left(1-2c_2\right)+\frac{3}{4}{q}_2^2}{{\left(n_2{p}_2\right)}^2}+O\left(n^{-3}\right)\end{array}$ (4)

when $q_1=1-p_1$ and $q_2=1-p_2$.

The first setting to find the optimal point $\left(c_1,c_2\right)$ that minimizes the smallest bias of ${\hat{\theta }}_c$ is investigated. The first derivatives of the bias of ${\hat{\theta }}_c$ with respect to $c_1$ and $c_2$ are given by

$\frac{\partial B\left({\hat{\theta }}_c\right)}{\partial c_1}=\frac{1+\left(-1+n_1\right)p_1-2n_1{p}_1^2+c_1\left(-1+4p_1^2\right)}{n_1^2{p}_1^2}$ (5)

$\frac{\partial B\left({\hat{\theta }}_c\right)}{\partial c_2}=-\left[\frac{1+\left(-1+n_2\right)p_2-2n_2{p}_2^2+c_2\left(-1+4p_2^2\right)}{n_2^2{p}_2^2}\right]$ (6)

Setting $\frac{\partial B\left({\hat{\theta }}_c\right)}{\partial c_1}=0$ and $\frac{\partial B\left({\hat{\theta }}_c\right)}{\partial c_2}=0$, the critical roots $\left(c_1,c_2\right)$ are obtained as follows:

$c_1=\frac{-1+p_1-n_1{p}_1+2n_1{p}_1^2}{-1+4p_1^2}$, $p_1\ne 0.5$ (7)

$c_2=\frac{-1+p_2-n_2{p}_2+2n_2{p}_2^2}{-1+4p_2^2}$, $p_2\ne 0.5$ (8)

The sufficient conditions of the second-order derivatives of the bias of ${\hat{\theta }}_c$ to guarantee the solutions of (7) and (8) being minimum points are that the determinants $D_1$ and $D_2$, evaluated at critical points, are all positive where $D_1=\frac{{\partial }^{2}B\left({\hat{\theta }}_c\right)}{\partial c_1^2}$ and $D_2=\left|\begin{array}{cc}\frac{{\partial }^{2}B\left({\hat{\theta }}_c\right)}{\partial c_1^2}& \frac{{\partial }^{2}B\left({\hat{\theta }}_c\right)}{\partial c_1\partial c_2}\\ \frac{{\partial }^{2}B\left({\hat{\theta }}_c\right)}{\partial c_1\partial c_2}& \frac{{\partial }^{2}B\left({\hat{\theta }}_c\right)}{\partial c_2^2}\end{array}\right|$. Unfortunately, it is impossible to find the minimum point $\left(c_1,c_2\right)$ such that ${\hat{\theta }}_c=\log \left(\frac{X_1+c_1}{n_1+2c_1}\right)-\log \left(\frac{X_2+c_2}{n_2+2c_2}\right)$ has the smallest bias for all $\left(p_1,p_2\right)$. In addition, the solution to (7) and (8) is practically inflexible because it cannot be applied when $p_1$ and $p_2$ equal 0.5. Therefore, we further consider the other settings, such as the smallest average bias and the smallest MSE, to find the optimal points $\left(c_1,c_2\right)$. However, these settings still fail, they cannot provide the optimal points $\left(c_1,c_2\right)$.

Another successful setting for finding the optimal values $\left(c_1,c_2\right)$ is coming from an average MSE or equivalently known as Bayes risk. Suppose that the squared error loss function is given by $\text{Loss}={\left({\hat{\theta }}_c-\theta \right)}^2$. The risk as the expectation of loss or the MSE of ${\hat{\theta }}_c$ in such this case, is given by $\text{Risk}=MSE\left({\hat{\theta }}_c\right)=E{\left({\hat{\theta }}_c-\theta \right)}^2$. The prior distribution is usually constructed based on the previously scientific knowledge or the prior observed data; especially, in this paper we set the prior distribution based on Table 1. For given the prior uniform density $g\left(p_1,p_2\right)=1_{\left[0,1\right]}\times 1_{\left[0,1\right]}$ over [0, 1] × [0, 1], the Bayes risk of ${\hat{\theta }}_c$ (or the average MSE of ${\hat{\theta }}_c$ ) with respect to the Euclidean loss function is obtained as

$\text{Bayesrisk}=m\left(c_1,c_2\right)={\displaystyle {\int }_0^1{\displaystyle {\int }_0^{1}MSE\left({\hat{\theta }}_c\right)g\left(p_1,p_2\right)\text{d}{p}_1\text{d}{p}_2}}$

Table 1. High sparse data of a multi-center clinical trial in CALGB study (Cancer and Leukemia Group B) from Lipsitz et al. [17] and Cooper et al. [18].

*Either $X_{1j}$ or $X_{2j}$ is 0, leading to the infinity value problem of the natural logarithm.

$\begin{array}{c}m\left(c_1,c_2\right)={\displaystyle {\int }_0^1{\displaystyle {\int }_0^1\frac{q_1}{n_1{p}_1}}}+\frac{\left(c_1-2c_1{p}_1\right)\left(c_1-2c_1{p}_1-q_1\right)+q_1\left(1-2c_1\right)+\frac{3}{4}{q}_1^2}{{\left(n_1{p}_1\right)}^2}\\ +\frac{q_2}{n_2{p}_2}+\frac{\left(c_2-2c_2{p}_2\right)\left(c_2-2c_2{p}_2-q_2\right)+q_2\left(1-2c_2\right)+\frac{3}{4}{q}_2^2}{{\left(n_2{p}_2\right)}^2}\text{d}{p}_1\text{d}{p}_2\end{array}$

$m\left(c_1,c_2\right)=\frac{\left(-1+c_1+3c_1^2\right)n_2^2-n_1{n}_2^2+n_1^2\left(3.25-c_2-3c_2^2+n_2\right)}{n_1^2{n}_2^2}$ (9)

The first and second order of partial derivatives of $m\left(c_1,c_2\right)$ with respect to $c_1$ and $c_2$ are as follows: $\frac{\partial m\left(c_1,c_2\right)}{\partial c_1}=\frac{1+6c_1}{n_1^2}$, $\frac{\partial m\left(c_1,c_2\right)}{\partial c_2}=\frac{-1-6c_2}{n_2^2}$, $\frac{{\partial }^{2}m\left(c_1,c_2\right)}{\partial c_1^2}=\frac{6}{n_1^2}$, $\frac{{\partial }^{2}m\left(c_1,c_2\right)}{\partial c_2^2}=\frac{-6}{n_2^2}$, $\frac{{\partial }^{2}m\left(c_1,c_2\right)}{\partial c_1\partial c_2}=0$. Unfortunately, the result of Bayes risk shows that the critical point $\left(c_1,c_2\right)=\left(-\frac{1}{6},-\frac{1}{6}\right)$ is not a minimum point. With the conditions of $D_1=\frac{{\partial }^{2}m\left(c_1,c_2\right)}{\partial c_1^2}$ and $D_2=\left|\begin{array}{cc}\frac{{\partial }^{2}m\left(c_1,c_2\right)}{\partial c_1^2}& \frac{{\partial }^{2}m\left(c_1,c_2\right)}{\partial c_1\partial c_2}\\ \frac{{\partial }^{2}m\left(c_1,c_2\right)}{\partial c_1\partial c_2}& \frac{{\partial }^{2}m\left(c_1,c_2\right)}{\partial c_2^2}\end{array}\right|$, the critical point $\left(c_1,c_2\right)$ is a saddle point since $D_2=\frac{-36}{n_1^2{n}_2^2}<0$. However, in particular case, we let $c=c_1=c_2$ and try again to find the minimum point. Fortunately, with the condition of $n_2>n_1$, the optimal point c with the smallest average MSE is obtained by $c=-\frac{1}{6}$ ; equivalently, with the condition of $n_1>n_2$, the minimum point is $c=\frac{1}{6}$. The solution of $c=\frac{-1}{6}$ or $c=\frac{1}{6}$ looks well and can be considered as an appropriate one in practice.

3. Minimum MSE Weights of the Adjusted Summary Relative Risk Estimator

Under a common risk ratio overall k studies or homogeneity of risk ratios across k studies, the optimal weights $f_j$ are investigated to minimize the MSE of ${\hat{\theta }}_w$ of the form ${\hat{\theta }}_w={\displaystyle {\sum }_{j=1}^k{\hat{\theta }}_{cj}{f}_j}$, subject to the constraint on ${\displaystyle {\sum }_{j=1}^k{f}_j}=1$, where ${\hat{\theta }}_{cj}=\log \left({\widehat{RR}}_{cj}\right)=\log \left({\hat{p}}_{c1j}\right)-\log \left({\hat{p}}_{c2j}\right)$, ${\hat{p}}_{c1j}=\left(X_{1j}+c_1\right)/\left(n_{1j}+2c_1\right)$ and ${\hat{p}}_{c2j}=\left(X_{2j}+c_2\right)/\left(n_{2j}+2c_2\right)$. The MSE of ${\hat{\theta }}_w={\displaystyle {\sum }_{j=1}^k{\hat{\theta }}_{cj}{f}_j}$, used under a true common risk ratio $\theta$ across all k studies, is given by

$MSE\left({\hat{\theta }}_w\right)=E{\left({\hat{\theta }}_w-\theta \right)}^2=E{\left({\displaystyle {\sum }_{j=1}^k{\hat{\theta }}_{cj}{f}_j}-\theta \right)}^2$

To find the optimal weights $f_j$ with the constraint on ${\displaystyle {\sum }_{j=1}^k{f}_j}=1$, we form the auxiliary function $\varphi$ in extending $\text{MSE}\left({\hat{\theta }}_w\right)$ under the Lagrange’s method where $\lambda$ is a Lagrange multiplier as follows:

$\varphi =E{\left({\displaystyle {\sum }_{j=1}^k{\hat{\theta }}_{cj}{f}_j}-\theta \right)}^2+\lambda \left({\displaystyle {\sum }_{j=1}^k{f}_j}-1\right)$

$\varphi =E{\left({\displaystyle {\sum }_{j=1}^k{\hat{\theta }}_{cj}{f}_j}\right)}^2-2\theta E\left({\displaystyle {\sum }_{j=1}^k{\hat{\theta }}_{cj}{f}_j}\right)+{\theta }^2+\lambda \left({\displaystyle {\sum }_{j=1}^k{f}_j}-1\right)$

$\begin{array}{c}\varphi =Var\left({\displaystyle {\sum }_{j=1}^k{\hat{\theta }}_{cj}{f}_j}\right)+{\left(E\left({\displaystyle {\sum }_{j=1}^k{\hat{\theta }}_{cj}{f}_j}\right)\right)}^2-2\theta E\left({\displaystyle {\sum }_{j=1}^k{\hat{\theta }}_{cj}{f}_j}\right)\\ +{\theta }^2+\lambda \left({\displaystyle {\sum }_{j=1}^k{f}_j}-1\right)\end{array}$

$\begin{array}{c}\varphi ={\displaystyle {\sum }_{j=1}^k{f}_j^{2}Var\left({\hat{\theta }}_{cj}\right)}+{\left({\displaystyle {\sum }_{j=1}^k{f}_{j}E\left({\hat{\theta }}_{cj}\right)}\right)}^2-2\theta \left({\displaystyle {\sum }_{j=1}^k{f}_{j}E\left({\hat{\theta }}_{cj}\right)}\right)\\ +{\theta }^2+\lambda \left({\displaystyle {\sum }_{j=1}^k{f}_j}-1\right)\end{array}$

$\varphi ={\displaystyle {\sum }_{j=1}^k{f}_j^{2}Var\left({\hat{\theta }}_{cj}\right)}+{\left({\displaystyle {\sum }_{j=1}^k{f}_{j}E\left({\hat{\theta }}_{cj}\right)}-\theta \right)}^2+\lambda \left({\displaystyle {\sum }_{j=1}^k{f}_j}-1\right)$

From the previous section, we have the following results:

$\begin{array}{l}{V}_j=Var\left({\hat{\theta }}_{cj}\right)=\frac{q_{1j}}{n_{1j}{p}_{1j}}+\frac{q_{1j}\left(1-2c_1\right)+\frac{1}{2}{q}_{1j}^2}{{\left(n_{1j}{p}_{1j}\right)}^2}+\frac{q_{2j}}{n_{2j}{p}_{2j}}\\ +\frac{q_{2j}\left(1-2c_2\right)+\frac{1}{2}{q}_{2j}^2}{{\left(n_{2j}{p}_{2j}\right)}^2}+O\left(n^{-3}\right)\end{array}$

and

$\begin{array}{c}{E}_j=E\left({\hat{\theta }}_{cj}\right)=\left(\log p_{1j}-\log p_{2j}\right)+(\frac{2c_1-1+p_{1j}\left\{1-2\left(2c_1\right)\right\}}{2n_{1j}{p}_{1j}}\\ -\frac{2c_2-1+p_{2j}\left\{1-2\left(2c_2\right)\right\}}{2n_{2j}{p}_{2j}})+O\left(n^{-2}\right)\end{array}$

The partial derivatives with respect to $\lambda$ and $f_j$ yield

$\frac{\partial \varphi }{\partial \lambda }={\displaystyle {\sum }_{j=1}^k{f}_j}-1$

$\frac{\partial \varphi }{\partial f_j}=2f_j{V}_j+2\left({\displaystyle {\sum }_{j=1}^k{f}_j{E}_j}-\theta \right)E_j+\lambda$

Setting $\frac{\partial \varphi }{\partial \lambda }=0$, $\frac{\partial \varphi }{\partial f_j}=0$ and solving the equations, it yields

$\lambda =-\frac{2}{a}-\frac{2b\left({\displaystyle {\sum }_{j=1}^k{f}_j{E}_j}-\theta \right)}{a}$

$f_j=\left(\frac{V_j^{-1}\left(1+{\tau }_j\theta \right)}{a}\right)-\left(\frac{V_j^{-1}{\tau }_j}{a+{\displaystyle {\sum }_{j=1}^k{\tau }_j{E}_j{V}_j^{-1}}}\right)\left(\frac{{\displaystyle {\sum }_{j=1}^k{V}_j^{-1}{E}_j\left(1+{\tau }_j\theta \right)}}{a}\right)$

where $a={\displaystyle {\sum }_{j=1}^k\frac{1}{V_j}}={\displaystyle {\sum }_{j=1}^k{V}_j^{-1}}$, $b={\displaystyle {\sum }_{j=1}^k\frac{E_j}{V_j}}$, ${\tau }_j=a{E}_j-b$.

More details of derivation can be found from the work of Viwatwongkasem et al. [19] that used the minimum MSE weights for the risk difference in multi-center studies. After replacing the unknown parameter quantities with their sample estimates, the minimum MSE weighted estimate is obtained in the following:

${\hat{\theta }}_w={\displaystyle {\sum }_{j=1}^k{\hat{\theta }}_{cj}{\hat{f}}_j}$ (10)

where

${\hat{\theta }}_{cj}=\log \left({\widehat{RR}}_{cj}\right)=\log \left({\hat{p}}_{c1j}\right)-\log \left({\hat{p}}_{c2j}\right)$, ${\hat{p}}_{c1j}=\left(X_{1j}+c_1\right)/\left(n_{1j}+2c_1\right)$,

${\hat{p}}_{c2j}=\left(X_{2j}+c_2\right)/\left(n_{2j}+2c_2\right)$, ${\hat{q}}_{c1j}=1-{\hat{p}}_{c1j}$, ${\hat{q}}_{c2j}=1-{\hat{p}}_{c2j}$

${\hat{f}}_j=\left(\frac{{\hat{V}}_j^{-1}\left(1-{\hat{\tau }}_j{\hat{\theta }}_{pool}\right)}{\hat{a}}\right)-\left(\frac{{\hat{V}}_j^{-1}{\hat{\tau }}_j}{\hat{a}+{\displaystyle {\sum }_{j=1}^k{\hat{\tau }}_j{\hat{E}}_j{\hat{V}}_j^{-1}}}\right)+\left(\frac{{\displaystyle {\sum }_{j=1}^k{\hat{E}}_j{\hat{V}}_j^{-1}\left(1+{\hat{\tau }}_j{\hat{\theta }}_{pool}\right)}}{\hat{a}}\right)$

$\begin{array}{c}{\hat{E}}_j=\left(\log {\hat{p}}_{c1j}-\log {\hat{p}}_{c2j}\right)+(\frac{\left\{c_1-0.5+\left(0.5-2c_1\right){\hat{p}}_{c1j}\right\}}{n_{1j}{\hat{p}}_{c1j}}\\ -\frac{\left\{c_2-0.5+\left(0.5-2c_2\right){\hat{p}}_{c2j}\right\}}{n_{2j}{\hat{p}}_{c2j}})+O\left(n^{-2}\right)\end{array}$

$\begin{array}{c}{\hat{V}}_j=\frac{{\hat{q}}_{c1j}}{n_{1j}{\hat{p}}_{c1j}}+\frac{{\hat{q}}_{c1j}\left(1-2c_1\right)+0.5{\hat{q}}_{c1j}^2}{{\left(n_{1j}{\hat{p}}_{c1j}\right)}^2}+\frac{{\hat{q}}_{c2j}}{n_{2j}{\hat{p}}_{c2j}}\\ +\frac{{\hat{q}}_{c2j}\left(1-2c_2\right)+0.5{\hat{q}}_{c2j}^2}{{\left(n_{2j}{\hat{p}}_{c2j}\right)}^2}+O\left(n^{-3}\right)\end{array}$

${\hat{\theta }}_{pool}=\log \left({\hat{p}}_1\right)-\log \left({\hat{p}}_2\right)$

where

${\hat{p}}_1=\frac{{\displaystyle {\sum }_{j=1}^k{n}_{1j}{\hat{p}}_{1j}}}{{\displaystyle {\sum }_{j=1}^k{n}_{1j}}}=\frac{{\displaystyle {\sum }_{j=1}^k{X}_{1j}}}{{\displaystyle {\sum }_{j=1}^k{n}_{1j}}}$, ${\hat{p}}_2=\frac{{\displaystyle {\sum }_{j=1}^k{n}_{2j}{\hat{p}}_{2j}}}{{\displaystyle {\sum }_{j=1}^k{n}_{2j}}}=\frac{{\displaystyle {\sum }_{j=1}^k{X}_{2j}}}{{\displaystyle {\sum }_{j=1}^k{n}_{2j}}}$

$\hat{a}={\displaystyle {\sum }_{j=1}^k\frac{1}{{\hat{V}}_j}}={\displaystyle {\sum }_{j=1}^k{\hat{V}}_j^{-1}}$, $\hat{b}={\displaystyle {\sum }_{j=1}^k\frac{{\hat{E}}_j}{{\hat{V}}_j}}$, ${\hat{\tau }}_j=\hat{a}{\hat{E}}_j-\hat{b}$

Assuming that a normal approximation is reliable, the minimum MSE weights Z-test for testing $H_0:\theta ={\theta }_0$, ( ${\theta }_0=\log R{R}_0$ ) is

$Z_{cw}=\frac{{\displaystyle {\sum }_{j=1}^k{\hat{f}}_j{\hat{\theta }}_{cj}}-{\theta }_0}{\sqrt{{\displaystyle {\sum }_{j=1}^k{\hat{f}}_j^2\hat{V}ar\left({\hat{\theta }}_{cj}|H_0\right)}}}=\frac{{\displaystyle {\sum }_{j=1}^k{\hat{f}}_j{\hat{\theta }}_{cj}}-{\theta }_0}{\sqrt{{\displaystyle {\sum }_{j=1}^k{\hat{f}}_j^2{\hat{V}}_j}}}$

We will reject $H_0$ at $\alpha$ level for two-sided test if $\left|Z_{cw}\right|>Z_{\alpha /2}$ where $Z_{\alpha /2}$ is the upper $100\left({\alpha }^{th}\right)$ percentile of the standard normal distribution. Alternatively, reject $H_0$ when the p-value (p) is less than or equal to $\alpha$ where $p=2\left(1-\Phi \left(\left|Z_{cw}\right|\right)\right)$ and $\Phi \left(Z\right)$ is the cumulative standard normal distribution.

4. Other Well-Known Estimators and Tests for Making Inference for a Common Relative Risk

Under a common risk ratio or homogeneity of risk ratios across k studies, we wish to compare the performance of the minimum MSE weighted estimator adjusted by various points $c=1/6$, $c=1/3$, $c=1/2$ with the other well-known summary risk ratio estimators, such as the Mantel-Haenszel (MH) estimator and the weighted least square (WLS) estimator or equivalently known as the inverse-variance weighted estimator via a simulation study. According to these well-known estimators, we will present briefly both estimators.

Mantel-Haenszel Weights (MH)

For Mantel-Haenszel (MH) relative risk estimator overall centers/studies from binomial data, the estimator has been proposed by

${\widehat{RR}}_{MH}=\frac{{\displaystyle {\sum }_{j=1}^k{n}_{2j}{X}_{1j}/N_j}}{{\displaystyle {\sum }_{j=1}^k{n}_{1j}{X}_{2j}/N_j}}$ (11)

where $N_j=n_{1j}+n_{2j}$.

The variance estimator of the log-relative risk of Mantel-Haenszel was proposed by Green and Robin [20] based on unconditional binomial distribution is given by

$\hat{V}\left[\log {\widehat{RR}}_{MH}\right]=\frac{{\displaystyle {\sum }_{j=1}^k{D}_j}}{\left({\displaystyle {\sum }_{j=1}^k{R}_j}\right)\left({\displaystyle {\sum }_{j=1}^k{S}_j}\right)}$ (12)

where:

$D_j=\left(n_{2j}{n}_{1j}{t}_j-x_{1j}{x}_{2j}{N}_j\right)/N_j^2$, $R_j=x_{1j}{n}_{2j}/N_j$, $S_j=x_{2j}{n}_{1j}/N_j$, $t_j=x_{1j}+x_{2j}$ and $N_j=n_{1j}+n_{2j}$. Note that under a binomial sparse-data model, the Mantel-Haenszel relative risk is consistent in sparse stratification.

Assuming that a normal approximation is valid, the Mantel-Haenszel’s Z-test for testing $H_0:\theta ={\theta }_0$, ( ${\theta }_0=\log R{R}_0$ ) is

$Z_{MH}=\frac{\log {\widehat{RR}}_{MH}-{\theta }_0}{\sqrt{\hat{V}\left(\log {\widehat{RR}}_{MH}|H_0\right)}}$

We will reject $H_0$ at $\alpha$ level for two-sided test if $\left|Z_{MH}\right|>Z_{\alpha /2}$ where $Z_{\alpha /2}$ is the upper $100\left({\alpha }^{th}\right)$ percentile of the standard normal distribution. Alternatively, reject $H_0$ when the p-value (p) is less than or equal to $\alpha$ where $p=2\left(1-\Phi \left(\left|Z_{MH}\right|\right)\right)$ and $\Phi \left(Z\right)$ is the cumulative standard normal distribution.

Weighted least square (WLS) estimator

The weighted least square (WLS) estimator or equivalently known as the inverse-variance weighted estimator of the log-relative risk overall centers/studies is

${\hat{\theta }}_{WLS}=\log {\widehat{RR}}_{WLS}={\displaystyle {\sum }_{j=1}^k{w}_j\log {\widehat{RR}}_j}/{\displaystyle {\sum }_{j=1}^k{w}_j}$ (13)

where

$\log {\widehat{RR}}_j=\log {\hat{p}}_{1j}-\log {\hat{p}}_{2j}=\log \left(\frac{X_{1j}}{n_{1j}}\right)-\log \left(\frac{X_{2j}}{n_{2j}}\right)$

$w_j=\frac{1}{V\left[\log {\widehat{RR}}_j\right]}={\left(\frac{1-p_{1j}}{n_{1j}{p}_{1j}}+\frac{1-p_{2j}}{n_{2j}{p}_{2j}}\right)}^{-1}$ (14)

Practically, the weights $w_j$ are often replaced by their estimates.

${\hat{w}}_j=\frac{1}{\hat{V}\left[\log {\widehat{RR}}_j\right]}={\left(\frac{1-{\hat{p}}_{1j}}{n_{1j}{\hat{p}}_{1j}}+\frac{1-{\hat{p}}_{2j}}{n_{2j}{\hat{p}}_{2j}}\right)}^{-1}=\left(\frac{1}{X_{1j}}-\frac{1}{n_{1j}}+\frac{1}{X_{2j}}-\frac{1}{n_{2j}}\right)$

The variance of the summary estimator ${\hat{\theta }}_{WLS}$ is given by

$V\left({\hat{\theta }}_{WLS}\right)=\frac{{\displaystyle {\sum }_{j=1}^k{w}_j^{2}V}\left[\log {\widehat{RR}}_j\right]}{{\left({\displaystyle {\sum }_{j=1}^k{w}_j}\right)}^2}=\frac{{\displaystyle {\sum }_{j=1}^k{w}_j^2}\left(\frac{1}{w_j}\right)}{{\left({\displaystyle {\sum }_{j=1}^k{w}_j}\right)}^2}=\frac{1}{{\displaystyle {\sum }_{j=1}^k{w}_j}}$ (15)

Assuming that a normal approximation is valid, the weighted least square Z-test for testing $H_0:\theta ={\theta }_0$, ( ${\theta }_0=\log R{R}_0$ ) is

$Z_{WLS}=\frac{{\displaystyle {\sum }_{j=1}^k{\hat{w}}_j\log {\widehat{RR}}_j}/\left({\displaystyle {\sum }_{j=1}^k{\hat{w}}_j}\right)-{\theta }_0}{\sqrt{1/\left({\displaystyle {\sum }_{j=1}^k{\hat{w}}_j}\right)}}$

We will reject $H_0$ at $\alpha$ level for two-sided test if $\left|Z_{WLS}\right|>Z_{\alpha /2}$ where $Z_{\alpha /2}$ is the upper $100\left({\alpha }^{th}\right)$ percentile of the standard normal distribution. Alternatively, reject $H_0$ when the p-value (p) is less than or equal to $\alpha$ where $p=2\left(1-\Phi \left(\left|Z_{WLS}\right|\right)\right)$ and $\Phi \left(Z\right)$ is the cumulative standard normal distribution.

5. Simulation Plan for the Estimation, Studying the Type I Error and the Power of the Test

We present here a simulation study using the following designs:

Parameters: Let the common relative risk be some constants (RR = 1, 2 and 4) and generate the baseline proportion $p_{2j}$ in the control arm for the j^th study from a uniform distribution in which the range corresponds to the values of RR. If $RR=1$, then $p_{2j}~U\left(0,0.9\right)$ ; if $RR=2$, then $p_{2j}~U\left(0,0.45\right)$ ; and if $RR=4$, then $p_{2j}~U\left(0,0.23\right)$. The correspondent proportion risks in the exposure arm is $p_{1j}=p_{2j}\times RR$. The sample size $n_{1j}$ and $n_{2j}$ in each study are fixed and varied as 4, 8, 16, 32 and 100. The number of studies is 1, 2, 4, 8, 16, 32 and 100.

Statistic 1 (estimation of RR): Generate $x_{1j}$ and $x_{2j}$ from the binomial distribution with parameter $\left(p_{1j},n_{1j}\right)$ and $\left(p_{2j},n_{2j}\right)$ for each study j ( $j=1,2,\cdots ,k$ ). All summary estimates are calculated. The procedure is replicated 5000 times. From these replicates, we compute bias, variance, and mean square error (MSE) for the adjusted relative risk estimator with the proposed weight to compare the performance with Mantel-Haenszel (MH) estimator and the weighted least square estimator.

Statistic 2 (studying the type I error under $H_0:RR=R{R}_0$ ): Generate $x_{1j}$ and $x_{2j}$ from the binomial distribution with parameter $\left(p_{1j},n_{1j}\right)$ and $\left(p_{2j},n_{2j}\right)$ and replicate these generations 5000 times for every procedure. From these replicates, the number of null hypothesis rejections when $H_0$ is true under three Z-tests is counted for the actual type I error.

$\text{TheactualtypeIerror}=\frac{\text{Numberofrejectionsof}{H}_0\text{when}{H}_0\text{istrue}}{\text{Numberofreplicates}\left(5000\text{times}\right)}$

Statistic 3 (studying the power of the test under $H_1$ ): Before comparing the power of tests, all test statistics should be calibrated to handle the same type I error rate under the null hypothesis. Under the alternative hypothesis with the random effects model, the powers of three candidate tests are compared. We need to revise the parameter setting for studying the power of the test. Let $p_{2j}$ be a uniform distribution over [0, 0.25] and we assume that ${\theta }_j=\log R{R}_j$ follows as ${\theta }_j=0.1+U_m=0.1+mm\left(2U-1\right)$ where $U_m$ is a uniform $\left[-mm,mm\right]$ random variable for a given $mm=0.2,0.4,0.6$, or equivalently, U is a uniform over (0, 1). Note that these parameter settings provide: $E\left({\theta }_j\right)=0.1$ and $Var\left({\theta }_j\right)={\left(2\times mm\right)}^2/12$. Consequently, we still have $p_{1j}=p_{2j}\times RR$. Binomial variates $x_{1j}$ and $x_{2j}$ are also generated with parameter $\left(p_{1j},n_{1j}\right)$ and $\left(p_{2j},n_{2j}\right)$, respectively. All proposed test statistics under this alternative hypothesis are computed and replicated 5000 times. From these replicates, the number of null hypothesis rejection is counted for the power of the test.

$\text{Theactualpowerofthetest}=\frac{\text{Numberofrejectionsof}{H}_0\text{when}{H}_1\text{istrue}}{\text{Numberofreplicates}\left(5000\text{times}\right)}$

6. Results from Simulation Studies

6.1. Comparative Performance for Point Estimation

Under a constant of relative risk (RR = 1, 2 and 4), the performance in terms of bias, variance, and mean square error (MSE) of several summary relative risk estimators are compared. Results show that increasing k can decrease the variance and the MSE of all estimators and the increase of both $n_{ij}$ can also decrease the variance of all estimators while fixing k. The unbalance cases of $n_{ij}$ ( $i=1,2$ and $j=1,\cdots ,k$ ) have less affected on the order performance of MSE estimators. The summary adjusted relative risk estimator in meta-analysis study of size k has shrinkage estimator to be a simple adjusted relative risk estimator in one single study case. The optimal point ( $c=c_1=c_2$ ) providing the bias, variance, and MSE of ${\hat{\theta }}_w$, adjusted by $c=-1/6$ is identical to $c=1/6$. By these, the numerical evidence has confirmed the derivation process of finding the root c and it is very useful in practice.

For a single center study ( $k=1$ ), regardless of a true value of RR, the proposed estimator adjusted by $c=1/3$ performs the best with smallest MSE.

For a multi-center study of size k, when $RR=1$, the WLS estimator is the best in sense of the smallest MSE ignoring the sample size $n_{ij}$. In case $RR=2$, $k=4$ the proposed estimator adjusted by c performs the best with the smallest MSE. Another issue, when $RR=2,k=16$, the MH estimator achieves the smallest MSE when sample sizes are small. The MSE of the WLS estimator and the proposed estimator adjusted by c are close together when the sample sizes are moderate and large ( $n_{ij}\ge 16$ ). For the case $RR=4$ and $k=4$ or 16, the MH estimator well performs generally with the smallest MSE. Some comparisons of the performances of all summary estimators when $n_{1j}=4,8$ and $n_{2j}=4,8,16,32$ are depicted in Figures 1(a)-(f) for bias; Figures 2(a)-(f) for variance; and Figures 3(a)-(f) for MSE.

In summary for the performance of estimators, regardless of the true values of RR, the MH estimator achieves the best performance with the smallest MSE when the study size is rather large ( $k\ge 16$ ) and the sample sizes within each study are small. The MSE of WLS estimator and the proposed-weight estimator adjusted by $c=1/6$, $c=1/3$, $c=1/2$ are close together and they are the best when the sample sizes are moderate to large ( $n_{ij}\ge 16$ ) while the study size is rather small.

6.2. Studying the Type I Error

A type I error comparison between the tests is considered by comparing the actual (empirical) type I error ( $\hat{\alpha }$ ) with the nominal level of significance. In this study, the evaluation of the ability to control type I error probability for two-sided tests is based on Cochran limits as follows.

At $\alpha =0.01$ significant level, the actual $\hat{\alpha }$ value is between [0.005, 0.015].

At $\alpha =0.05$ significant level, the actual $\hat{\alpha }$ value is between [0.04, 0.06].

At $\alpha =0.10$ significant level, the actual $\hat{\alpha }$ value is between [0.08, 0.12].

If the actual type I error or the empirical alpha lies within those Cochran limits, then the statistical test can control type I error rate.

For a single study ( $k=1$ ), regardless of the true values of RR, it is unfortunate that almost all tests shown in Table 2 cannot control type I error rate. There are some type I error rates lying in the Cochran limits.

For a meta-analysis study of size k, also displayed in Table 2, regardless of the true values of RR, the Mantel-Haenszel’s Z-Test can control type I error rate when the sample size either in treatment or control group is moderate to large ( $n_{1j}\ge 16$ or $n_{2j}\ge 16$ ). In addition, the weighted least square Z-test and the proposed weight Z-test adjusted by $c=1/6$, $c=1/2$, $c=1$ can handle type I error rate when both $n_{ij}$ are large. But the proposed Z-test adjusted by $c=2$ cannot control type I error rate into Cochran’s limit at almost all situations.

6.3. Studying the Power of the Test

Usually, the empirical power of the tests will be compared under the same type I error value. In summary, the Mantel-Haenszel Z-test performs best when $n_{ij}$ is moderate to large with satisfying the type I error value within Cochran’s range limit, regardless of the study size k. The inverse variance weighted Z-test is good when $n_{ij}$ is large. For the proposed Z-test adjusted by $c=1$, $c=2$, the Z-tests perform well under the same type I error value when k is small ( $k\le 4$ ) and $n_{ij}$ is large ( $n_{ij}\ge 32$ ). The results in Table 3 illustrate the comparison of the power of tests under the same rate of the type I error.

Figure 1. Bias comparison of relative risk between well-known estimators and adjusted relative risk estimators at k = 1, k = 4 and k = 16 for n₁ = 4 and n₂ = 4, 8, 16 and 32 ((a)-(c)) and at k = 1, k = 4 and k = 16 for n₁ = 8 and n₂ = 8, 16 and 32 ((d)-(f)).

Figure 2. Variance comparison of relative risk between well-known estimators and adjusted relative risk estimators at k = 1, k = 4 and k = 16 for n₁ = 4 and n₂ = 4, 8, 16 and 32 ((a)-(c)) and at k = 1, k = 4 and k = 16 for n₁ = 8 and n₂ = 8, 16 and 32 ((d)-(f)).

Figure 3. MSE comparison of relative risk between well-known estimators and adjusted relative risk estimators at k = 1, k = 4 and k = 16 for n₁ = 4 and n₂ = 4, 8, 16 and 32 ((a)-(c)) and at k = 1, k = 4 and k = 16 for n₁ = 8 and n₂ = 8, 16 and 32 ((d)-(f)).

Table 2. Comparison of the empirical type I error for testing $H_0:RR=R{R}_0$ in a center study ( $k=1$ ) and multi-center study ( $k=4$, 16 and 32) at 5% significant level.

Bold Values denote that the statistical tests can control the type I error.

Table 3. Comparison of the empirical power of test (percent).

^#denotes that the statistical tests cannot control type I error rate under $H_0:RR=1$.

7. Discussion and Recommendation

The main question rises which continuity correction values are the best choice for the adjusted relative risk estimator in a center study and multi-center study with sparse data. Due to the conventional continuity correction, most of investigators such as Yate [1], Lane [3], Stijnen et al., [4], Lui and Lin [6], Sankey et al., [7], Gart and Zwefel [8], Walter [9] and Cox [10], suggest to use $c=0.5$. However, in this study with the smallest average MSE of ${\hat{\theta }}_c=\log {\widehat{RR}}_c$, the optimal point $\left(c_1,c_2\right)=\left(1/6,1/6\right)$ can perform the best for the point estimation. The minimum point $\left(c_1,c_2\right)=\left(1/6,1/6\right)$ of ${\hat{\theta }}_c$ agrees with the suggestion of Turkey [12], which is very useful and the most appropriate in a practice way.

For estimation of fixing $\theta$ ( $\theta =\log RR$ ) in a center ( $k=1$ ), regardless of a true value of RR, the proposed estimator adjusted by $c=c_1=c_2=1/3$ performs the best with the smallest MSE. For a meta-analysis study of size k, in general the MH estimator achieves the smallest MSE when the sample size $n_{ij}$ is small while the study size is rather large ( $k\ge 16$ ). The MSE of the WLS estimator and the proposed estimator adjusted by the various values of c are closed together and they are the best when the sample sizes are moderate to large ( $n_{1j}\ge 16$ and $n_{2j}\ge 16$ ) while the study size is rather small. This finding is consonant with the work of Viwatwongkasem et al. [19]. Since the true value of RRis usually not available in practice as mentioned earlier, we suggest to choose the proposed relative risk estimator adjusted by $c=1/6$ that can minimize the Bayes risk with respect to uniform prior (0, 1) and Euclidean loss function.

For the empirical power of the test under $H_1:{\theta }_j=0.1+U_m\left(-mm,mm\right)$, regardless of the study size k, the MH Z-test performs the best with the highest power when both $n_{ij}$ are moderate to large. The inverse-variance weighted Z-test is good when $n_{ij}$ is large. In accordance with Soulakova and Bright [21], the empirical power of the test will increase when the sample sizes increase. For the power of the proposed Z-test adjusted by $c=1$, $c=2$, the Z-test has the higher performance under the same type I error when k is small ( $k\le 4$ ) and $n_{ij}$ is large ( $n_{ij}\ge 32$ ).

If we don’t know information about parameter RR, we recommend to use the adjusted estimator $\hat{\theta }$ by using continuity correction defined by $c_1=c_2=1/6$, or 1/3, or 1/2 in a center study. For a multi-center study of size k, we recommend to use adjusted ${\hat{\theta }}_w$ defined by $c_1=c_2=1/6$, or 1/3, or 1/2, including optimal weights ${\hat{f}}_j$.

Generally, the effect of exposures or the effect of treatments with binary outcomes covers the risk difference, the relative risk, and the odds ratio. Obviously, all three conventional effect estimators have the same problem of the zero values in sparse data. The conventional proportion estimate of $\hat{p}=X/n$ is in need of replacement by ${\hat{p}}_c=\left(X+c\right)/\left(n+2c\right)$ to solve this problem. Therefore, the recommendation for a further study is to use these ideas, such as the smallest MSE, the smallest Bayes risk to fine the appropriate point $\left(c_1,c_2\right)$ in estimating the odds ratio parameter.

Acknowledgements

The study was partially supported for publication by the Faculty of Public Health, Mahidol University, Bangkok, Thailand.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Yates, F. (1934) Contingency Tables Involving Small Numbers and the Chi-Squared Test. Supplement to the Journal of the Royal Statistical Society, 1, 217-235. https://doi.org/10.2307/2983604
[2]	Lane, P.W. (2013) Meta-Analysis of Incidence of Rare Events. Statistical Methods in Medical Research, 22, 117-132. https://doi.org/10.1177/0962280211432218
[3]	Stijnen, T., Hamza, T.H. and Özdemir, P. (2010) Random Effects Meta-Analysis of Event Outcome in the Framework of the Generalized Linear Mixed Model with Applications in Sparse Data. Statistics in Medicine, 29, 3046-3067. https://doi.org/10.1002/sim.4040
[4]	White, I.R., Daniel, R. and Royston, P. (2010) Avoiding Bias Due to Perfect Prediction in Multiple Imputation of Incomplete Categorical Variables. Computational Statistics and Data Analysis, 54, 2267-2275. https://doi.org/10.1016/j.csda.2010.04.005
[5]	Lui, K.J. and Lin, C.D. (2003) A Revisit on Comparing the Asymptotic Interval Estimators of Odds Ratio in a Single 2 × 2 Table. Biometrical Journal, 45, 226-237. https://doi.org/10.1002/bimj.200390008
[6]	Sankey, S.S., Weissfeld, L.A., Fine, M.J. and Kapoor, W. (1996) An Assessment of the Use of the Continuity Correction for Sparse Data in Meta-Analysis. Communications in Statistics—Simulation and Computation, 25, 1031-1056. https://doi.org/10.1080/03610919608813357
[7]	Gart, J.J. and Zweifel, J.R. (1967) On the Bias of Various Estimators of the Logit and Its Variance with Application to Quantal Bioassay. Biometrika, 54, 181-187. https://doi.org/10.1093/biomet/54.1-2.181
[8]	Walter, S.D. (1975) The Distribution of Levin’s Measure of Attributable Risk. Biometrika, 62, 371-372. https://doi.org/10.1093/biomet/62.2.371
[9]	Cox, D.R. (1970) The Continuity Correction. Biometrika, 57, 217-219. https://doi.org/10.1093/biomet/57.1.217
[10]	Li, L. and Wang, X. (2017) Meta-Analysis of Rare Binary Events in Treatment Groups with Unequal Variability. Statistical Methods in Medical Research, 28, 263-274. https://doi.org/10.1177/0962280217721246
[11]	Tukey, J.W. (1977) Exploratory Data Analysis. Addison-Wesley Publishing Company, Boston.
[12]	Sánchez-Meca, J. and Marín-Martínez, F. (2000) Testing the Significance of a Common Risk Difference in Meta-Analysis. Computational Statistics and Data Analysis, 33, 299-313. https://doi.org/10.1016/S0167-9473(99)00055-9
[13]	Böhning, D. and Viwatwongkasem, C. (2005) Revisiting Proportion Estimators. Statistical Methods in Medical Research, 14, 147-169. https://doi.org/10.1191/0962280205sm393oa
[14]	Agresti, A. and Caffo, B. (2000) Simple and Effective Confidence Intervals for Proportions and Differences of Proportions Result from Adding Two Successes and Two Failures. The American Statistician, 54, 280-288. https://doi.org/10.1080/00031305.2000.10474560
[15]	McClave, J.T. and Sincich, T.T. (2012) Statistics. Pearson Education, London.
[16]	Pettigrew, H.M., Gart, J.J. and Thomas, D.G. (1986) The Bias and Higher Cumulants of the Logarithm of a Binomial Variate. Biometrika, 73, 425-435. https://doi.org/10.1093/biomet/73.2.425
[17]	Lipsitz, S.R., Dear, K.B.G., Laird, N.M. and Molenberghs, G. (1998) Tests for Homogeneity of the Risk Difference When Data Are Sparse. Biometrics, 54, 148-160.
[18]	Cooper, M.R., Dear, K.B.G., McIntyre, O.R., Ozer, H., Ellerton, J.Cannellos, G., Duggan, B. and Schiffer, C. (1993) A Randomized Clinical Trial Comparing Melphalan/Prednisone with and without α-2b Interferon in Newly-Diagnosed Patients with Multiple Myeloma: A Cancer and Leukemia Group B Study. Journal of Clinical Oncology, 11, 155-160. https://doi.org/10.1200/JCO.1993.11.1.155
[19]	Viwatwongkasem, C., Jitthavech, J., Böhning, D. and Lorchirachoonkul, V. (2012) Minimum MSE Weights of Adjusted Summary Estimator of Risk Difference in Multi-Center Studies. Open Journal of Statistics, 2, 48-59. https://doi.org/10.4236/ojs.2012.21006
[20]	Greenland, S. and Robin, J.M. (1985) Estimation of a Common Effect Parameter from Sparse Follow-up Data. Biometrics, 45, 55-68.
[21]	Soulakova, J.N. and Bright, B.C. (2013) Applications of Asymptotic Confidence Intervals with Continuity Corrections for Asymmetric Comparisons in Noninferiority Trials. Pharmaceutical Statistics, 12, 147-155. https://doi.org/10.1002/pst.1566