Meta-Analysis of Multi-Arm Trials Using Binomial Approach

Hathaikan Chootrakool, Pichet Treewai

Department of Biostatistics, Faculty of Public Health, Mahidol University, Bangkok, Thailand.

DOI: 10.4236/ojs.2022.121002                

Abstract

Most meta-analysis has concentrated on combining of treatment effect measures based on comparisons of two treatments. Meta-analysis of multi-arm trials is a key component of submission to summarize evidence from all possible studies. In this paper, an exact binomial model is proposed by using logistic regression model to compare different treatment in multi-arm trials. Two approaches such as unconditional maximum likelihood and conditional maximum likelihood have been determined and compared for the logistic regression model. The proposed models are performed using the data from 27 randomized clinical trials (RCTs) which determine the efficacy of antiplatelet therapy in reduction venous thrombosis and pulmonary embolism.

Keywords

Meta-Analysis, Multi-Arm Trials, Binomial Approach

1. Introduction

There has been a massive growth in the number of randomized clinical trials (RCTs) since the first RCT was introduced in the well-known streptomycin trial in 1946 [1] . In making some of this information more readily available, an attempt is made to pull together the existing evidence in a form that can be used by researchers or statisticians; this is called systematic review. The aim of systematic reviews is to identify and summarize the findings of all possible studies addressing the clinical question of the review. Systematic review reduces the large quantities of information to a manageable size. An important aspect of most reviews is the quantitative synthesis of results; thus, meta-analysis is the statistical part of systematic review. The main purpose of meta-analysis is to increase the precision of the conclusions of a review. With statistical perspective, it is able to detect treatment effects with greater power and estimate these effects with greater precision than any single study. A meta-analysis from systematic reviews of Antiplatelet Trialists’ Collaboration [2] was applied in the paper.

Most meta-analysis has focused on summarizing treatment effect measures based on the comparison of two treatments or called arms, sometimes also called interventions or exposures. In this comparison, two groups of individual studies are exposed to two different treatments. Standard two-arm RCTs are frequently used in clinical research due in part to its relative simplicity of design and interpretation. At its most basic, one power, one significance level and one magnitude of difference are analyzed for two-arm comparisons. Conclusions are straightforward: either the two arms are shown to be different or they are not. The implementation for the model is not complicated. When more than two arms are included in meta-analysis, complexity ensues. These types of dataset are called multi-arm trials [3] [4] [5] although some authors call it mixed treatment comparison (MTC) [6] [7] , or multiple treatment comparison [4] [8] . Additionally, some authors call network meta-analysis for indirect head-to-head evidence for multiple evidences from RCTs [9] - [14] .

For multi-arm trials data, Chootrakool and Shi [3] introduced the normal approximation model using an empirical logistic transform requiring a large number of individual observations and the probability of being the case to be not too near zero or one. If the number of individual observations is small, the normal approximation model is not suitable. The number of samples in a single study should usually be larger than 20 [4] . In this paper, an exact binomial model is introduced to fit the binary multi-arm trials data. There are two alternative maximum likelihood approaches that can be used to make inferences for the unknown parameters in the logistic regression model. These are unconditional method and conditional method. The logistic regression model has become increasingly popular with the easy availability of appropriate computer routines. Many authors have described maximum likelihood estimation procedures which turn out to be iterative [15] . There have been a large number of studies about unconditional and conditional methods, for example in Tritchler [16] , Sartori [17] , Lee et al. [18] and Caterina et al. [19] .

Most existing methods for meta-analysis of multi arm trials use the logistic regression model with the unconditional approach. Lu and Ades [20] introduced the Bayesian hierarchical model for multi-arm trials using the unconditional method to estimate unknown parameters. Weber et al. [21] used random effect model using the unconditional maximum likelihood method for the meta-analysis. They conducted a simulation study comparing two zero-cell corrections under the ordinary random-effects model. Logistic regression model is introduced using the conditional approach in the paper.

There has the direct and indirect comparison for RCTs in a meta-analysis [9] [20] , this has led to network meta-analysis. Seide et al. [12] perform simulation study for sparse network of trial including multi-arm trials. They compared Bayesian and Frequentist methods in random effects network meta-analysis. Jenkins et al. [13] introduced methods for inclusion of evidence in network meta-analysis. More applications of network meta-analysis can be found in Greco et al. [22] , Wang [23] and Zhao et al. [14] .

The structure of the paper is arranged as follows. The data structure of multi-arm trials is firstly introduced in Section 2. Fitting the logistic regression model is proposed in the next section. Unconditional maximum likelihood approach for the model including the standard error of maximum likelihood estimators (MLEs) are described in Section 4. Similarly, conditional maximum likelihood approach for the model is presented in Section 5. In Section 6, the logistic regression model is illustrated with the unconditional and conditional approaches with the data. The final section is discussion and conclusion including the advantages and the limitations of the two approaches.

2. The Data Structure of Multi-Arm Trials

Let the indices $i=1,\cdots ,M$ and $j=0,1,\cdots ,K$ stand for the studies and the treatments respectively, where the index $j=0$ stands for the control group. To make multi-arm comparisons, suppose that there is M RCTs comparing K + 1 treatments. For the ith study, let $r_{ij}$ represent the number of cases on treatment j and let $n_{ij}$ be the total number in the corresponding group. Let ${\pi }_{ij}$ be the probability of being the case for a patient received the treatment j in the ith study. The $r_{ij}$ has a binomial distribution

$r_{ij}~Bin\left({\pi }_{ij},n_{ij}\right);i=1,\cdots ,M;\text{and}j=0,1,\cdots ,K.$

A data structure of multi-arm trials shall be defined by introducing an index set $J_i$ comprising the treatments involved in the ith study. The data structure is shown as

$\mathcal{D}=\left\{\left(r_{ij},n_{ij}\right):i=1,\cdots ,M;j\in J_i\right\}.$

Some studies of meta-analysis might not have all the treatments available. For example, some studies might compare fewer than K + 1 treatments, or some baseline treatments may be different, or both cases could occur simultaneously. The data structure is analogous to an incomplete-blocks design. For example, from Antiplatelet data, the 8^th - 17th studies are composed of two treatment comparison unlike the 1^st - 7th studies comparing three treatments.

Let b(i) denote baseline treatment or called study-specific reference treatment [4] corresponding to the ith study, which can be the control group or any other treatments. As mentioned earlier about indirect comparison, in a situation that the treatments in some studies cannot be compared directly to the control group, we need to use evidence from the external studies. To make it clear, if b(i) = 0 then the direct comparison is involved in this study. Conversely, if $b\left(i\right)\ne 0$ then the study makes indirect comparison. Let $J_{\left(i\right)}=J_i\backslash \left\{b\left(i\right)\right\}$ represent the set of treatments involved in the ith study but excluding the baseline treatment b(i). Let $k_i$ and $k_{i+1}$ denote the number of treatments in the sets $J_i$ and $J_i$ respectively. The $r_{ib\left(i\right)}$ and $r_{ij}$ are binomially distributed, respectively as $Bin\left(n_{ib\left(i\right)},{\pi }_{ib\left(i\right)}\right)$ and $Bin\left(n_{ij},{\pi }_{ij}\right)$ for $i=1,\cdots ,M$ and $j\in J_{\left(i\right)}$ . Let D be set of studies that make the direct comparison of treatment thus let I be set of indirect comparison of treatment.

The multi-arm trials data is consisted of 27 RCTs from systematic reviews of Antiplatelet Trialists’ Collaboration [2] in total as shown in Table 1. The studies compare three treatments: aspirin plus dipyridamole (A), aspirin alone (B) and control group (C), where seven trials compare aspirin plus dipyridamole, aspirin alone and control group (i.e. comparing all A, B and C), ten trials compare aspirin plus dipyridamole and control group (i.e. comparing A and C) and ten trials compare aspirin alone and control group (i.e. comparing B and C) The “event” in Table 1 represents the number of patients in whom deep venous thrombosis was detected by systematic fibrinogen scans or venography, or both, after general and orthopedic surgery and in high risk medical patients. The “total” represents the number of patients controlled in each group.

3. Fitting the Logistic Regression Model

This section illustrates how to fit the logistic regression model to the binary data related to multi-arm trials including the direct and indirect comparisons. Logistic regression is a regression model for a binomially distributed response/dependent variable. It is useful for modelling the probability of an event occurring as a function of other factors. Logistic regression is part of a category of statistical models called generalized linear models and uses the logit as its link function. Logistic regression can be used only with two types of dependent variables: one is a categorical dependent variable that has exactly two categories (i.e. a binary or dichotomous variable). The other is a continuous dependent variable that has values in the range 0 to 1 representing the probability values or the proportions. The names for logistic regression used in various other application areas are logistic model or logit model. Logistic regression is similar to linear regression. In linear regression, ultimate objective for the study may be either estimation of the coefficient values, or prediction of the response value. One significant difference between logistic and linear models is that the linear model has a continuous response variable and the logistic model uses a binary or dichotomous response. Logistic regression models for the ith study can be defined by

$\log \left(\frac{{\pi }_{ib\left(i\right)}}{1-{\pi }_{ib\left(i\right)}}\right)={\alpha }_i$ (1)

$\log \left(\frac{{\pi }_{ij}}{1-{\pi }_{ij}}\right)={\alpha }_i+{\delta }_{i,b\left(i\right)j},j\in J\left(i\right)$ (2)

where ${\alpha }_i$ and ${\delta }_{i,b\left(i\right)j}$ are trial effect and treatment effect, they will be detailed next.

Table 1. RCTs of Antiplatelet data.

3.1. Trial Effect

Two assumptions are usually made about the trial effect ${\alpha }_i$ . The first one is that the trial effects are assumed to be study-level effects, which means that the ${\alpha }_i$ ’s are different parameters and are treated as nuisance parameters in the model. The M different unknown parameters are needed to include in the model. The second one is that the model may be assumed for the ${\alpha }_i$ ’s. A special case is to assume that the trial effect is a fixed effect, defined by ${\alpha }_i={\alpha }_0$ . Conversely, it may be assumed to be a random effect, given by ${\alpha }_i~N\left({\mu }_{\alpha 0},{\tau }_{\alpha 0}^2\right)$ , where ${\mu }_{\alpha 0}$ is the overall mean of the trial effect and ${\tau }_{\alpha 0}^2$ measures the magnitude of the variation between the studies. Most of the existing methods therefore used the first assumption. However, the number of unknown parameters is the same as the number of studies if the first assumption of the trial effect is made. It will result in some theoretical and computational problems. The accuracy of the estimation depends on the sample size of each study not the overall sample size of the pool in the meta-analysis. Therefore, the ${\alpha }_i$ ’s are assumed to be different in the paper.

3.2. Treatment Effect

The treatment effect ${\delta }_{i,b\left(i\right)j}$ can be direct treatment effect if $i\in D$ or indirect treatment effect if $i\in I$ : they are defined as follows.

${\delta }_{i,b\left(i\right)j}=\{\begin{array}{l}{\delta }_{i,0j}\sim N\left({\mu }_{0j},{\tau }_{0j}^2\right)\text{if}i\in D,\\ {\delta }_{i,0j}-{\delta }_{i,0b\left(i\right)}~N\left({\mu }_{0j}-{\mu }_{0b\left(i\right)},{\tau }_{0j}^2+{\tau }_{0b\left(i\right)}^2-2{\rho }_{jb\left(i\right)}{\tau }_{0j}{\tau }_{0b\left(i\right)}\right)\text{if}i\in I.\end{array}$

where ${\rho }_{jb\left(i\right)}$ is the correlation coefficient between ${\delta }_{i,0j}$ and ${\delta }_{i,0b\left(i\right)}$ . For example, if the data, includes 3 treatments such as treatment A, B, C and C stands for control group. The baseline treatment for some studies exclude the control group C (such as A and B) thus the indirect treatment effect can be written as

${\delta }_{i,AB}={\delta }_{i,AC}-{\delta }_{i,BC}\sim N\left({\mu }_{AC}-{\mu }_{BC},{\tau }_{AC}^2+{\tau }_{BC}^2-2{\rho }_{AB}{\tau }_{AC}{\tau }_{BC}\right),i=I.$

Next, we shall consider the treatment effect in a matrix form. Let ${\delta }_{i,0}$ and ${\mu }_0$ represent the vectors of ${\left({\delta }_{i,0j},j=1,\cdots ,K\right)}^t$ and ${\left({\mu }_{0j},j=1,\cdots ,K\right)}^t$ respectively where the superscript t stands for matrix transposition and let ${\Omega }_0$ represent the $K\times K$ covariance matrix. The model ${\delta }_{i,0}$ can be written as

${\delta }_{i,0}~MVN\left({\mu }_0,{\Omega }_0\right).$

This is called the basic model of random treatment effect. Let $F_{ij}$ be the index vector of length K consisting of elements 0 and 1 corresponding to ${\delta }_{i,b\left(i\right)j}$ , given by

$F_{ij}=\{\begin{array}{l}\left(0,\cdots ,\underset{b\left(i\right)\text{th}}{\underset{︸}{0}},\cdots \underset{j\text{th}}{\underset{︸}{1}},\cdots ,0\right)\text{if}i\in D,\\ \left(0,\cdots ,\underset{b\left(i\right)\text{th}}{\underset{︸}{-1}},\cdots ,\underset{j\text{th}}{\underset{︸}{1}},\cdots ,0\right)\text{if}i\in I.\end{array}$

Now, the random effect ${\delta }_{i,b\left(i\right)j}$ can be written in the form of ${\delta }_{i,0}$ and $F_{ij}$ :

${\delta }_{i,b\left(i\right)j}=F_{ij}{\delta }_{i,0}~N\left(F_{ij}{\mu }_0,F_{ij}{\Omega }_0{F}_{ij}^t\right).$

As before, the covariance between the treatment effects ${\delta }_{i,b\left(i\right)j}$ and ${\delta }_{i,b\left(i\right)k}$ for $j\ne k$ and $j,k\in J_{\left(i\right)}$ may be dependent. For the ith study, let $F_i$ be the following $k_i\times K$ matrix

$F_i={\left(F_{ij}\right)}_{k_i\times K}\text{for}j\in J_{\left(i\right)},$

Let ${\delta }_i$ denote the vector ${\left({\delta }_{i,b\left(i\right)j},j\in J_{\left(i\right)}\right)}^t$ then we have

${\delta }_i=F_i{\delta }_{i,0}~MVN\left({\mu }_i,{\Omega }_i\right),$

where ${\mu }_i=F_i{\mu }_0$ and ${\Omega }_i=F_i{\Omega }_0{F}_i^t.$

The ${\alpha }_i$ are assumed to be different and the ${\delta }_{i,b\left(i\right)j}$ are assumed to be a random effect as ${\delta }_i$ . The models in the Equation (1) and Equation (2) can be used for both treatment comparisons. From model in Equation (2), log ${\pi }_{ij}/1-{\pi }_{ij}$ is called the logistic transform of probability ${\pi }_{ij}$ or alternatively log odds ${\pi }_{ij}$ or logit ${\pi }_{ij}$ . Having considered the properties of logit ${\pi }_{ij}$ , the term ${\pi }_{ij}/1-{\pi }_{ij}$ is the odds of an unsuccessful outcome from a patient treated with treatment j and so logit ${\pi }_{ij}$ is the log odds of being the case. It is easily seen that a value of ${\pi }_{ij}$ in the range (0, 1) corresponds to a value of logit ${\pi }_{ij}$ in (−∞, ∞). As ${\pi }_{ij}$ approaches to 0, logit ${\pi }_{ij}$ approaches to −∞; as ${\pi }_{ij}$ approaches to 1, logit ${\pi }_{ij}$ approaches to ∞ and for ${\pi }_{ij}=0.5$ , logit ${\pi }_{ij}=0$ . After some rearrangement, the logistic regression models in Equation (1) and Equation (2) have respectively equivalent formulations as

${\pi }_{ib\left(i\right)}=\left(\frac{{\text{e}}^{{\alpha }_i}}{1+{\text{e}}^{{\alpha }_i}}\right)$ and ${\pi }_{ij}=\left(\frac{{\text{e}}^{{\alpha }_i+{\delta }_{i,b\left(i\right)j}}}{1+{\text{e}}^{{\alpha }_i+{\delta }_{i,b\left(i\right)j}}}\right)$ . (3)

There are two alternative maximum likelihood (ML) approaches, the unconditional and conditional approaches that can be used to estimate the unknown parameters in a logistic regression model.

4. Unconditional Maximum Likelihood Approach

Generally, unconditional ML estimation is preferred if the number of parameters in the model is small relative to the number of studies in a meta-analysis [25] .

4.1. Probability Functions

To demonstrate the unconditional ML estimation, let $p\left(r_{ib\left(i\right)}|{\alpha }_i\right)$ and $p\left(r_{ij}|{\alpha }_i,{\delta }_{i,b\left(i\right)j}\right)$ denote the probability functions associated with the distributions of $r_{ib\left(i\right)}|{\alpha }_i$ and $r_{ij}|{\alpha }_i,{\delta }_{i,b\left(i\right)j}$ respectively for $i=\text{1},\cdots ,M$ and $j\in J_{\left(i\right)}$ , defined as follows.

For the baseline treatment,

$p\left(r_{ib\left(i\right)}|{\alpha }_i\right)=\left(\begin{array}{c}{n}_{ib\left(i\right)}\\ r_{ib\left(i\right)}\end{array}\right){\pi }_{ib\left(i\right)}^{r_{ib\left(i\right)}}{\left(1-{\pi }_{ib\left(i\right)}\right)}^{n_{ib\left(i\right)}-r_{ib\left(i\right)}}=\left(\begin{array}{c}{n}_{ib\left(i\right)}\\ r_{ib\left(i\right)}\end{array}\right)\frac{{\text{e}}^{{\alpha }_i{r}_{ib\left(i\right)}}}{{\left(1+{\text{e}}^{{\alpha }_i}\right)}^{n_{ib\left(i\right)}}}$ (4)

For the treatments j, $j\in J_{(\; i\; )}$

$p\left(r_{ij}|{\alpha }_i,{\delta }_{i,b\left(i\right)j}\right)=\left(\begin{array}{c}{n}_{ij}\\ r_{ij}\end{array}\right){\pi }_{ij}^{r_{ij}}{\left(1-{\pi }_{ij}\right)}^{n_{ij}-r_{ij}}=\left(\begin{array}{c}{n}_{ij}\\ r_{ij}\end{array}\right)\frac{{\text{e}}^{\left({\alpha }_i+{\delta }_{i,b\left(i\right)j}\right)r_{ij}}}{{\left(1+{\text{e}}^{{\alpha }_i+{\delta }_{i,b\left(i\right)j}}\right)}^{n_{ij}}}.$ (5)

The combination in Equation (4) represents the number of possible combinations of $n_{ib\left(i\right)}$ taken $r_{ib\left(i\right)}$ at a time. The ${\pi }_{ib\left(i\right)}$ in the middle term of Equation (4) is substituted from Equation (3) and $\left(1-{\pi }_{ib\left(i\right)}\right)$ becomes $1/1+{\text{e}}^{{\alpha }_i}$ . The combination in Equation (5) can be analyzed in the same way.

4.2. The Unconditional Likelihood

From the probability functions $p\left(r_{ib\left(i\right)}|{\alpha }_i\right)$ and $p\left(r_{ij}|{\alpha }_i,{\delta }_{i,b\left(i\right)j}\right)$ , the trial effects ${\alpha }_i$ ’s are study-level effects. They are assumed to be different and also included in both probability functions. While the ${\delta }_{i,b\left(i\right)j}$ is a random effect, thus the $p\left(r_{ij}|{\alpha }_i,{\delta }_{i,b\left(i\right)j}\right)$ involves the vector of random effects, ${\delta }_i$ . The standard method of handling a probability function which involves random variables that have a fully specified probability is to integrate the probability function with respect to the distribution of those variables. To deal with the random effects ${\delta }_i$ , let $r_{\left(i\right)}$ be the vector ${\left(r_{ij},j\in J_{\left(i\right)}\right)}^t$ . The probability function $p\left(r_{ij}|{\delta }_i\right)$ will be integrated with respect to ${\delta }_i$ . The $p\left(r_{\left(i\right)}\right)$ contains $k_i$ integrals, which is the number of treatments in the set $J_{\left(i\right)}$ , and is given by

$p\left(r_{\left(i\right)}\right)={\displaystyle {\int }_{{\delta }_i}{\displaystyle {\prod }_{j\in J_{\left(i\right)}}p}}\left(r_{ij}|{\delta }_i\right)\varphi \left({\delta }_i;{\mu }_i,{\Omega }_i\right)\text{d}{\delta }_i,$ (6)

where $\varphi \left({\delta }_i;{\mu }_i,{\Omega }_i\right)$ is the probability density function of the normal distribution with mean ${\mu }_i$ and covariance ${\Omega }_i$ given by

$\varphi \left({\delta }_i;{\mu }_i,{\Omega }_i\right)=\frac{1}{{\left(2\pi \right)}^{k_i/2}{\left|{\Omega }_i\right|}^{1/2}}{\text{e}}^{-{\left({\delta }_i-{\mu }_i\right)}^{\prime }{\Omega }_i^{-1}\left({\delta }_i-{\mu }_i\right)/2}$ (7)

The integral in Equation (6) can be calculated numerically; one way to do it is to use the Gauss-Hermite method. To apply Gauss-Hermite approximation, the probability function $p\left(r_{\left(i\right)}\right)$ for the ith study can be estimated by

$p\left(r_{\left(i\right)}\right)\approx {\pi }^{k_i/2}{\displaystyle {\sum }_{n_1=1}^{l_1}{w}_{n_1}^{\left(1\right)}}\cdots {\displaystyle {\sum }_{n_{k_i}=1}^{l_{k_i}}{w}_{n_{k_i}}^{\left(k_i\right)}}\left\{{\displaystyle {\prod }_{j\in J_{\left(i\right)}}\left(\begin{array}{c}{n}_{ij}\\ r_{ij}\end{array}\right)}\frac{{\text{e}}^{\left({\alpha }_i+\left({\mu }_i+\sqrt{2}{\Omega }_i^{1/2}{d}_{i,n}\right)\right)r_{ij}}}{{\left(1+{\text{e}}^{{\alpha }_i+\left({\mu }_i+\sqrt{2}{\Omega }_i^{1/2}{d}_{i,n}\right)}\right)}^{n_{ij}}}\right\}$ , (8)

where the sampling nodes are at ${\mu }_i+\sqrt{2}{\Omega }_i^{1/2}{d}_{i,n}$ and $d_{i,n}=\left(x_{n_1}^{\left(1\right)},\cdots ,x_{n_{k_i}}^{\left(k_i\right)}\right)$ . The vector $d_{i,n}$ depends on the number $k_i$ , which is the number of treatments comprising in the ith study. The resulting function $p\left(r_{\left(i\right)}\right)$ does not depend on the ${\delta }_i$ . For most practical purposes, $l_{k_i}$ need not be greater than 20, although some authors suggest using even smaller values [26] .

As before, let $\theta$ be the collection of all unknown parameters for the meta-analysis including all trial effects $\left({\alpha }_i,\cdots ,{\alpha }_M\right)$ , $\mu$ and $\Omega$ and let $r_i$ be the vector $\left(r_{ij},j\in J_i\right)$ . The likelihood function for the ith study can be written as

$L\left(\theta |r_i\right)={\displaystyle {\prod }_{j\in J_i}p\left(r_{ij}\right)}=p\left(r_{ib\left(i\right)}|{\alpha }_i\right)p\left(r_{\left(i\right)}\right).$ (9)

where $p\left(r_{ib\left(i\right)}|{\alpha }_i\right)$ and $p\left(r_{\left(i\right)}\right)$ are given in Equation (4) and Equation (8) respectively. Let $l_{u,i}=\log L\left(\theta |r_i\right)$ standing for the unconditional log-likelihood function of the logistic regression model for the ith study. The log-likelihood function of $\theta$ for the logistic regression models is given by

$l_u\left(\theta \right)={\displaystyle {\sum }_{i=1}^M{l}_{u,i}}\text{.}$ (10)

The number of ${\alpha }_i$ ’s is the same as the number of studies. The computation of MLEs may be quite unstable if the number of studies is large while the sample size of each study is small. As discussed earlier, this may also result in a biased or misleading estimate. A conditional approach is suggested to eliminate all nuisance parameters.

4.3. Asymptotic Variance-Covariance Matrix

This section shows how to calculate the standard errors for the MLEs of the logistic regression model using the unconditional approach. Since there are random effects in the model, some integrals are involved in the likelihood function. The unconditional log-likelihood function $l_u\left(\theta \right)$ in Equation (10) can be written as

$\begin{array}{c}{l}_u\left(\theta \right)={\displaystyle {\sum }_{i=1}^M\log p\left(r_{ib\left(i\right)}\right)}+{\displaystyle {\sum }_{i=1}^M\log p\left(r_{\left(i\right)}\right)}\\ ={\displaystyle {\sum }_{i=1}^M\log p\left(r_{ib\left(i\right)}\right)}+{\displaystyle {\sum }_{i=1}^M\log }{\displaystyle {\int }_{{\delta }_i}{\displaystyle {\prod }_{j\in J_{\left(i\right)}}p}}\left(r_{ij}|{\delta }_i\right)\varphi \left({\delta }_i;{\mu }_i,{\Omega }_i\right)\text{d}{\delta }_i.\end{array}$ (11)

Let $l_1$ and $l_2$ stand for the first and second terms of the above log-likelihood function, given by $l_u\left(\theta \right)=l_1+l_2$ . Three types of unknown parameters are involved in $\theta$ ; the trial effects, ${\alpha }_i$ ’s, the overall mean effects μ’s (for $\mu$ ), and the variances τ’s and the correlation coefficients ρ’s in the covariance matrix $\Omega$ . Let τ represent a parameter (either τ or ρ) involved in $\Omega$ . There is no random effect involved in $l_1$ .

First, the second-order partial derivative ${\partial }^2{l}_1/\partial {\alpha }_i^2$ can be calculated in the usual way; while the other terms are

$\frac{{\partial }^2{l}_1}{\partial {\alpha }_i{\alpha }_j}=\frac{{\partial }^2{l}_1}{\partial {\alpha }_i\mu }=\frac{{\partial }^2{l}_1}{\partial {\alpha }_i\tau }=0$ , $i\ne j$ and $i,j\in \left\{1,\cdots ,M\right\}$ .

Next, consider the second term of Equation (11), for notational convenience, let $P_i\left({\delta }_i\right)$ represent the function ${\displaystyle {\prod }_{j\in J_{\left(i\right)}}p\left(r_{ij}|{\delta }_i\right)}$ in $l_2$ . Now the term $l_2$ takes the form

$l_2\left({\alpha }_1,\cdots ,{\alpha }_M,\mu ,\Omega \right)={\displaystyle \underset{i=1}{\overset{M}{\sum }}\log }{\displaystyle {\int }_{{\delta }_i}p\left({\delta }_i\right)\varphi \left({\delta }_i\right)\text{d}{\delta }_i}={\displaystyle \underset{i=1}{\overset{M}{\sum }}{l}_{2i}},$

where $\varphi \left(\delta \right)$ is the density of the multivariate normal distribution with mean ${\mu }_i$ and variance matrix ${\Omega }_i$ , and $l_{2i}$ is a summand of the log-likelihood involving the integrals ( $\log {\displaystyle {\int }_{{\delta }_i}{P}_i\left({\delta }_i\right)\varphi \left({\delta }_i\right)\text{d}{\delta }_i}$ ). The first-order partial derivatives relating to $l_2$ are shown as follows

$\frac{\partial l_2}{\partial {\alpha }_i}={\displaystyle \underset{i=1}{\overset{M}{\sum }}{\text{e}}^{-l_{2i}}}{\displaystyle {\int }_{{\delta }_i}\frac{\partial P_i\left({\delta }_i\right)}{\partial {\alpha }_i}\varphi \left({\delta }_i\right)\text{d}{\delta }_i},$

$\frac{\partial l_2}{\partial \mu }={\displaystyle \underset{i=1}{\overset{M}{\sum }}{\text{e}}^{-l_{2i}}}{\displaystyle {\int }_{{\delta }_i}{P}_i\left({\delta }_i\right)\frac{\partial \varphi \left({\delta }_i\right)}{\partial \mu }\text{d}{\delta }_i},$

$\frac{\partial l_2}{\partial \tau }={\displaystyle \underset{i=1}{\overset{M}{\sum }}{\text{e}}^{-l_{2i}}}{\displaystyle {\int }_{{\delta }_i}{P}_i\left({\delta }_i\right)\frac{\partial \varphi \left({\delta }_i\right)}{\partial \tau }\text{d}{\delta }_i}.$

Similarly, the second-order partial derivatives are

$\frac{{\partial }^2{l}_2}{\partial {\alpha }_i^2}=\underset{i=1}{\overset{M}{{\displaystyle \sum }}}\left({\text{e}}^{-l_{2i}}{\displaystyle {\int }_{{\delta }_i}\frac{{\partial }^2{P}_i\left({\delta }_i\right)}{\partial {\alpha }_i^2}\varphi \left({\delta }_i\right)\text{d}{\delta }_i}-{\left({\text{e}}^{-l_{2i}}{\displaystyle {\int }_{{\delta }_i}\frac{\partial P\left({\delta }_i\right)}{\partial {\alpha }_i}\varphi \left({\delta }_i\right)\text{d}{\delta }_i}\right)}^2\right),$ (12)

$\frac{{\partial }^2{l}_2}{\partial {\mu }^2}=\underset{i=1}{\overset{M}{{\displaystyle \sum }}}\left({\text{e}}^{-l_{2i}}{\displaystyle {\int }_{{\delta }_i}{P}_i}\left({\delta }_i\right)\frac{{\partial }^2\varphi \left({\delta }_i\right)}{\partial {\mu }^2}\text{d}{\delta }_i-{\left({\text{e}}^{-l_{2i}}{\displaystyle {\int }_{{\delta }_i}{P}_i}\left({\delta }_i\right)\frac{\partial \varphi \left({\delta }_i\right)}{\partial \mu }\text{d}{\delta }_i\right)}^2\right),$ (13)

$\frac{{\partial }^2{l}_2}{\partial {\tau }^2}=\underset{i=1}{\overset{M}{{\displaystyle \sum }}}\left({\text{e}}^{-l_{2i}}{\displaystyle {\int }_{{\delta }_i}{P}_i}\left({\delta }_i\right)\frac{{\partial }^2\varphi \left({\delta }_i\right)}{\partial {\tau }^2}\text{d}{\delta }_i-{\left({\text{e}}^{-l_{2i}}{\displaystyle {\int }_{{\delta }_i}{P}_i}\left({\delta }_i\right)\frac{\partial \varphi \left({\delta }_i\right)}{\partial \tau }\text{d}{\delta }_i\right)}^2\right),$ (14)

$\frac{{\partial }^2{l}_2}{\partial \mu \partial \tau }=\underset{i=1}{\overset{M}{{\displaystyle \sum }}}\left({\text{e}}^{-l_{2i}}{\displaystyle {\int }_{{\delta }_i}{P}_i}\left({\delta }_i\right)\frac{{\partial }^2\varphi \left({\delta }_i\right)}{\partial \mu \partial \tau }\text{d}{\delta }_i-{\left({\text{e}}^{-l_{2i}}{\displaystyle {\int }_{{\delta }_i}{P}_i}\left({\delta }_i\right)\frac{\partial \varphi \left({\delta }_i\right)}{\partial \mu }\text{d}{\delta }_i\right)}^2\right).$ (15)

Note that the second-order partial derivative ${\partial }^2{l}_2/\partial {\alpha }_i\partial {\alpha }_j$ is equal to zero. The second-order partial derivatives of ${\partial }^2{l}_2/\partial {\mu }_i\partial {\mu }_j$ and ${\partial }^2{l}_2/\partial {\tau }_i\partial {\tau }_j$ can be expressed in similar equations to Equation (13) and Equation (14) respectively. The integrals in the first-order and second-order partial derivatives can be approximated by Gaussian quadrature.

From the log-likelihood $l_u\left(\theta \right)$ in Equation (11), the second-order partial derivatives for the observed Fisher information matrix can be calculated as

$\frac{{\partial }^2{l}_u}{\partial {\alpha }_i^2}=\frac{{\partial }^2{l}_1}{\partial {\alpha }_i^2}+\frac{{\partial }^2{l}_2}{\partial {\alpha }_i^2}$ , $\frac{{\partial }^2{l}_u}{\partial {\mu }^2}=\frac{{\partial }^2{l}_2}{\partial {\mu }^2}$ , $\frac{{\partial }^2{l}_u}{\partial {\mu }_i\partial {\mu }_j}=\frac{{\partial }^2{l}_2}{\partial {\mu }_i\partial {\mu }_j}$ ,

$\frac{{\partial }^2{l}_u}{\partial {\tau }^2}=\frac{{\partial }^2{l}_2}{\partial {\tau }^2}$ , $\frac{{\partial }^2{l}_u}{\partial {\tau }_i\partial {\tau }_j}=\frac{{\partial }^2{l}_2}{\partial {\tau }_i\partial {\tau }_j}$ , $\frac{{\partial }^2{l}_u}{\partial \mu \partial \tau }=\frac{{\partial }^2{l}_2}{\partial \mu \partial \tau }$ .

As set earlier, the second partial derivatives of ${\partial }^2{l}_u/\partial {\rho }^2$ and ${\partial }^2{l}_u/\partial \mu \partial \rho$ (and ${\partial }^2{l}_u/\partial \tau \partial \rho$ ) can be calculated in similar equations to ${\partial }^2{l}_u/\partial {\tau }^2$ and ${\partial }^2{l}_u/\partial \mu \partial \tau$ respectively. Notice that the second-order derivative of $l_1$ is only related in ${\partial }^2{l}_u/\partial {\alpha }_i^2$ . The matrix of second partial derivatives can be partitioned into a block matrix with null matrices in the off diagonals:

$\mathcal{H}\left(\theta \right)=\left(\begin{array}{cc}{\mathcal{H}}_{\alpha }\left(\theta \right)& 0\\ 0& {\mathcal{H}}_{\mu \tau \rho }\left(\theta \right)\end{array}\right)$ ,

where ${\mathcal{H}}_{\alpha }\left(\theta \right)$ and ${\mathcal{H}}_{\mu \tau \rho }\left(\theta \right)$ are the second-order partial derivatives about ${\alpha }_j$ , and μ, τ and ρ respectively. By multiplying $\mathcal{H}\left(\theta \right)$ by −1, the observed Fisher information matrix $\mathcal{J}\left(\theta \right)$ is obtained. The inverse of $\mathcal{J}\left(\theta \right)$ is the asymptotic variance-covariance matrix of MLEs and their standard errors are the square roots of the diagonal of $\mathcal{J}{\left(\theta \right)}^{-1}$ .

5. Conditional Maximum Likelihood Approach

Conditional likelihood is widely used in logistic regression models with binary data. In particular, this leads to accurate inferences for the parameters of interest and eliminates all nuisance parameters [25] . The conditional likelihood will be defined and the maximum likelihood estimation will be described in this section.

5.1. Conditional Likelihood

From the logistic regression models in Equation (1) and Equation (2), the conditional likelihood $r_i$ given that $C_i={\displaystyle {\sum }_{j\in J_i}{r}_{ij}}=c_i$ for the ith study, is given by

$f\left(r_i|C_i=c_i;{\delta }_i\right)=f\left(r_i|{\displaystyle {\sum }_{j\in J_i}{r}_{ij}}=c_i;{\delta }_i\right)=\frac{f\left(r_i|{\delta }_i\right)}{f\left({\displaystyle {\sum }_{j\in J_i}{r}_{ij}}=c_i|{\delta }_i\right)}.$ (16)

The conditional likelihood reflects the probability of the observed data configuration relative to the probability of all possible configurations of the given data. The numerator $f\left(r_i|{\delta }_i\right)$ is exactly the same as the unconditional likelihood obtained from Equation (4) and Equation (5). The denominator is what makes the conditional likelihood different from the unconditional likelihood; it sums the joint probability for all possible configurations. To derive the Equation (16), the conditional likelihood $r_i$ given $C_i$ can be simplified as

$f\left(r_i|C_i=c_i;{\delta }_i\right)=\frac{{{\displaystyle \prod }}_{j\in J_i}\left(\begin{array}{c}{n}_{ij}\\ r_{ij}\end{array}\right){\text{e}}^{{\delta }_{i{r}_{ij}}}}{{{\displaystyle \sum }}_{u_i\in U_i}\left(\begin{array}{c}{n}_{i0}\\ c_i-{{\displaystyle \sum }}_{j\in J_{\left(i\right)}}{u}_{ij}\end{array}\right){{\displaystyle \prod }}_{j\in J_{\left(i\right)}}\left(\begin{array}{c}{n}_{ij}\\ u_{ij}\end{array}\right){\text{e}}^{{\delta }_{i{u}_{ij}}}},$ (17)

where $u_i={\left(u_{ij},j\in J_{\left(i\right)}\right)}^t$ and

$U_i=\left\{u_i:0\le u_{ij}\le n_{ij},j\in J_{\left(i\right)}\text{and}{c}_i-n_{i0}\le \underset{j\in J_{\left(i\right)}}{{\displaystyle \sum }}{u}_{ij}\le c_i\right\}.$

Notice that this likelihood function does not involve any nuisance parameters ${\alpha }_i$ ’s and is a function of ${\delta }_i$ alone. The removal of the trial effects from the conditional likelihood is important because when the conditional likelihood is used, estimates are obtained only for the parameters of interest in the model and not for the ${\alpha }_i$ ’s.

5.2. Estimation

The conditional likelihood in Equation (17) has $k_i$ random effects so the likelihood $f\left(r_i|{\displaystyle {\sum }_{j\in J_i}{r}_{ij}}=c_i\right)$ involves $k_i$ integrations:

$f\left(r_i|{\displaystyle {\sum }_{j\in J_i}{r}_{ij}}=c_i\right)={\displaystyle {\int }_{{\delta }_i}f}\left(r_i|{\displaystyle {\sum }_{j\in J_i}{r}_{ij}}=c_i;{\delta }_i\right)\varphi \left({\delta }_i;{\mu }_i,{\Omega }_i\right)\text{d}{\delta }_i,$ (18)

where $\varphi \left({\delta }_i;{\mu }_i,{\Omega }_i\right)$ is the probability density function of multivariate normal distribution with mean ${\mu }_i$ and covariance ${\Omega }_i$ given in Equation (7). Gauss-Hermite approximation is applied to Equation (18) and obtain:

$f\left(r_i|{\displaystyle \underset{j\in J_i}{\sum }{r}_{ij}}=c_i\right)\approx {\pi }^{-k_i/2}\underset{n_1=1}{\overset{l_1}{{\displaystyle \sum }}}{w}_{n_1}^{\left(1\right)}\cdots \underset{n_{k_i}=1}{\overset{l_{k_i}}{{\displaystyle \sum }}}{w}_{n_{k_i}}^{\left(k_i\right)}f\left(r_i|{\displaystyle \underset{j\in J_i}{\sum }{r}_{ij}}=c_i;{\delta }_{i,n}\right),$ (19)

where $f\left(r_i|{\displaystyle {\sum }_{j\in J_i}{r}_{ij}}=c_i;{\delta }_{i,n}\right)$ is obtained from Equation (17) where the

sampling nodes is ${\delta }_{i,n}={\mu }_i+\sqrt{2}{\Omega }_i^{1/2}{d}_{i,n}$ and $d_{i,n}=\left(d_{n_1}^{\left(1\right)},\cdots ,d_{n_{k_i}}^{\left(k_i\right)}\right)$ . The likelihood for the ith study $L\left(\theta |r_i\right)$ can be written as

$L\left(\theta |r_i\right)=f\left(r_i|{\displaystyle \underset{j\in J_i}{\sum }{r}_{ij}}=c_i;{\delta }_i\right).$

The log-likelihood function of the logistic regression models using the conditional approach is

$l_c\left(\theta \right)=\log L\left(\theta |r\right)={\displaystyle {\sum }_{i=1}^M\log L}\left(\theta |r_i\right)$ (20)

By maximizing the conditional likelihood function over $\theta$ , an exact parameter estimate is obtained for $\theta$ , called the conditional maximum likelihood estimate. To calculate the standard error of their MLEs, the log-likelihood function in Equation (20) can be written as

$\begin{array}{c}{l}_c\left(\theta \right)={\displaystyle {\sum }_{i=1}^M\log f}\left(r_i|{\displaystyle {\sum }_{j\in J_i}{r}_{ij}}=c_i\right)\\ ={\displaystyle {\sum }_{i=1}^M\log }{\displaystyle {\int }_{{\delta }_i}f}\left(r_i|{\displaystyle {\sum }_{j\in J_i}{r}_{ij}}=c_i,{\delta }_i\right)\varphi \left({\delta }_i;{\mu }_i,{\Omega }_i\right)\text{d}{\delta }_i.\end{array}$ (21)

Let $P_i\left({\delta }_i\right)$ represent $f\left(r_i|{\displaystyle {\sum }_{j\in J_i}{r}_{ij}}=c_i,{\delta }_i\right)$ in the above equation. The second-order partial derivatives of ${\partial }^2{l}_2/\partial {\mu }^2$ , ${\partial }^2{l}_2/\partial {\tau }^2$ and ${\partial }^2{l}_2/\partial \mu \partial \tau$ are similar to Equation (13) - Equation (15) respectively. In a similar way to the previous section, the standard errors for the MLEs are obtained.

6. Application to Antiplatelet Therapy Data

From the data given in Table 1, the total number of individual observations is small thus the empirical log-odds model is not appropriate. The logistic regression model will be applied using the unconditional and conditional approaches with the data.

6.1. Unconditional Inference

From the data, there are 27 studies investigating the use of aspirin plus dipyridamole or aspirin alone in comparison with the control group. The studies compare three treatments: aspirin plus dipyridamole (A), aspirin alone (B) and the control treatment (C). Seven studies compare A, B and C, ten studies compare A and C and ten studies compare B and C. There is no indirect comparison for this dataset, so the set D is $\left\{1,\cdots ,27\right\}$ . The baseline treatment for all studies is the control group (b(i) = 0).

The indices $i=1,\cdots ,27$ and $j=0,1,2$ stand for the studies and the treatments C, A and B, respectively. The data is partitioned into three groups: $G_1=\left\{1,\cdots ,7\right\}$ , $G_2=\left\{8,\cdots ,17\right\}$ and $G_3=\left\{18,\cdots ,27\right\}$ . The sets $J_i$ and $J_{\left(i\right)}$ are given by

$\{\begin{array}{l}{J}_i=\left\{0,1,2\right\},J_{\left(i\right)}=\left\{1,2\right\}\text{for}i\in G_1,\\ J_i=\left\{0,1\right\},J_{\left(i\right)}=\left\{1\right\}\text{for}i\in G_2,\\ J_i=\left\{0,2\right\},J_{\left(i\right)}=\left\{2\right\}\text{for}i\in G_3.\end{array}$ (22)

Let $r_{i0},r_{i1}$ and $r_{i2}$ be the numbers of patients who suffered reocclusions on treatments C, A and B respectively, where the ith study is in $G_1\cup G_2\cup G_3$ , $G_1\cup G_2$ and $G_1\cup G_3$ , respectively. The total numbers of patients are $n_{i0},n_{i1}$ and $n_{i2}$ Let ${\pi }_{i0},{\pi }_{i1}$ and ${\pi }_{i2}$ be the probabilities that patients have reocclusions on treatments C, A and B respectively in the ith study. The $r_{i0},r_{i1}$ and $r_{i2}$ are binomially distributed as

$r_{i0}~Bin\left({\pi }_{i0},n_{i0}\right),i\in G_1\cup G_2\cup G_3,$

$r_{i1}~Bin\left({\pi }_{i1},n_{i1}\right),i\in G_1\cup G_2,$

$r_{i2}~Bin\left({\pi }_{i2},n_{i2}\right),i\in G_1\cup G_3.$

The treatment effect ${\delta }_i$ for $G_1$ is defined as

$\left(\begin{array}{c}{\delta }_{i,01}\\ {\delta }_{i,02}\end{array}\right)~MVN\left(\left(\begin{array}{c}{\mu }_{01}\\ {\mu }_{02}\end{array}\right),\left(\begin{array}{cc}{\tau }_{01}^2& \rho {\tau }_{01}{\tau }_{02}\\ \rho {\tau }_{01}{\tau }_{02}& {\tau }_{02}^2\end{array}\right)\right)$ (23)

Logistic regression models for the data can be fitted using the Equation (1) and Equation (2) where b(i) = 0 and $J_{\left(i\right)}$ is given in Equation (22). Note that the trial effects are assumed to be different in each study. To define the unconditional likelihood function, let $r_{\left(i\right)}$ represent the vector $\left(r_{i1},r_{i2}\right)$ . The probability functions $p\left(r_{i0}\right)$ and $p\left(r_{\left(i\right)}\right)$ are formulated from Equation (4) and Equation (8) respectively.

From ${\delta }_i$ for $G_1$ , the correlation coefficient ρ between ${\delta }_{i,01}$ and ${\delta }_{i,02}$ is in the form ${\tau }_0^2/{\tau }_{01}{\tau }_{02}$ , where is obtained from ${\delta }_{i0}~N\left({\mu }_0,{\tau }_0^2\right)$ . Note that ${\mu }_0$ and ${\tau }_0^2$ are not estimable unless some other information is used. The assumption of homogeneity variance will be considered here. Suppose that all heterogeneity parameters are the same: ${\tau }_{01}={\tau }_{02}=\tau$ and the correlation coefficient takes the value 1/2. The unknown parameter $\theta$ for the models is $\left\{{\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_{27},{\mu }_{01},{\mu }_{02},{\tau }^2\right\}$ . The log-likelihood function $l_u\left(\theta \right)$ is obtained from Equation (10). By maximizing the log-likelihood function, the MLEs can be estimated. Also their standard errors are given by the observed Fisher information matrix.

The results for the treatment effects ${\delta }_{01}$ and ${\delta }_{02}$ are given in Table 2. SD and CI stand for standard deviation and confidence interval respectively. The trial effects are presented in Table 3. The overall means on the log odds ratio

Table 2. The results of the treatment effects for the model using the unconditional method.

Table 3. The trial effect of the model using the unconditional method.

(LOR) scale for ${\delta }_{01}$ and ${\delta }_{02}$ are −1.17849 (SD 0.08499) and −0.63700 (SD 0.03728), and the heterogeneity parameter is 0.0372 (SD 0.04752). On the odds ration (OR) scale, the means are 0.30774 and 0.52800 respectively. Their confidence interval (CI) can be calculated from the related CI on the LOR scale. To be concluded, treatments aspirin plus dipyridamole and aspirin only in antiplatelet therapy reduce deep venous thrombosis by over 70% and 45% respectively. The average of both treatments reduces deep venous thrombosis by over 55%.

6.2. Conditional Inference

The models and other parameters are similar to those defined in the unconditional method. The function $C_i$ for the data can be defined by

$\{\begin{array}{l}{C}_i=r_{i0}+r_{i1}+r_{i2}\text{for}i\in G_1,\\ C_i=r_{i0}+r_{i1}\text{for}i\in G_2,\\ C_i=r_{i0}+r_{i2}\text{for}i\in G_3.\end{array}$ (24)

Let $r_i$ denote the vector $\left(r_{i0},r_{i1},r_{i2}\right)$ . The conditional likelihood $f\left(r_i|C_i\right)$ for the ith study is given in Equation (18). To handle the random treatment effect ${\delta }_i$ , the likelihood function is approximated by Gaussian-Hermite approximation as defined in Equation (19). The unknown parameter $\theta$ for the models is $\left\{{\mu }_{01},{\mu }_{02},\tau \right\}$ . By using the log-likelihood function in Equation (20), the results of the models are given in Table 4. On the LOR scale, the overall mean effects for both treatment effects are −0.87516 (SD 0.04340) and −0.39000 (SD 0.31160) while their variation between studies is 0.37000 (SD 0.03900). Those means on the OR scale are 0.41679 and 0.67434. As before their confidence intervals are obtained from the related CI on the LOR scale. The results indicate that treatments aspirin plus dipyridamole and aspirin only produce a reduction in deep venous thrombosis by over 55% and 30% respectively. The average of both treatments in antiplatelet therapy reduces deep venous thrombosis by over 40%.

As seen from Table 2 and Table 4, the results from using the unconditional likelihood (on the LOR scale) are smaller than from using conditional likelihood. Note that those results are negative. That is to say that estimation with unconditional likelihood may cause underestimation or bias. Collaboration [2] summarized that antiplatelet therapy produced a highly significant (2p ≤ 0.00001) reduction in deep venous thrombosis of about 40%. The results from the model using the conditional likelihood support this.

7. Conclusions and Discussion

The logistic regression model was introduced for the exact binomial distribution. Two alternative approaches for making inferences were presented. The unconditional likelihood involves nuisance parameters (from the trial effects). If the number of studies (M) is large, it may lead to inconsistent estimators. Cox and Snell [15] concluded for the unconditional likelihood that if the number of studies (M) is large and the number of individual observations ( $n_{ij}$ ) is small then it makes estimation inaccurate and inconsistent. The conditional maximum likelihood approach was introduced for the model to eliminate all nuisance parameters. In making a choice between the two approaches, it is needed to consider the number of studies and the number of individual observations. However, the use

Table 4. The results of the treatment effects for the model using the conditional method.

of this method can be expensive in term of the cost of computer running time, especially if the number of individual observations is large. Using the unconditional maximum likelihood approach, note that if the number of studies is large and the number of individual observations is small then the estimate may be biased or misleading.

Some other methods can be used in the logistic regression model, for example, using a pseudo-loglikelihood, (see Severini [27] ); or the modified profile likelihood (see Bellio and Sartori [28] , Hamza, Houwelingen and Stijnen [29] and Caterina et al. [19] ). Gaussian-Hermite quadrature was used to calculate the integral forms of the probabilities including random effects in the likelihood functions for both approaches [30] . The approximation is reasonably effective for low-order integrations [31] . Implementing Gaussian-Hermite approximation, we used the function “gauss.quad” in the software R to estimate MLEs for the model. The number of integrands depends on the number of treatments involved in those studies. If this number is large then it makes the dimensionality of the integral large and so it cannot be approximated accurately. Other approximations such as Laplace approximation or Monte Carlo method can be used, see Shi and Copas [32] and Caterina [19] . Laplace approximation could make the calculation of second-order derivatives for the observed Fisher information matrix easier than using Gaussian approximation since there is no weight term in the approximation [33] [34] . Additionally, simulation studies could be conducted in further work to compare these two approaches in different scenarios including improving the approximation if multi-arm trials are more than three treatment comparisons.

The established binomial tree method could be applied in the further work [24] .

8. Conclusion

Simulation studies could be conducted in further work to compare these two approaches in different scenarios including improving the approximation if multi-arm trials are more than three treatment comparisons.

Acknowledgements

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

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Hill, S.A.B. (1990) Memories of the British Streptomycin Trial in Tuberculosis: The First Randomized Clinical Trial. Controlled Clinical Trials, 11, 77-79. https://doi.org/10.1016/0197-2456(90)90001-I
[2]	Antiplatelet Trialists’ Collaboration (1994) Collaborative Overview of Randomised Trials of Antiplatelet Therapy—III: Reduction in Venous Thrombosis and Pulmonary Embolism by Antiplatelet Prophylaxis among Surgical and Medical Patients. British Medical Journal, 308, 235-246. https://doi.org/10.1136/bmj.308.6923.235
[3]	Chootrakool, H. and Shi, J.Q. (2008) Meta-Analysis of Multi-Arm Trials Using Empirical Logistic Transform. The Open Medical Informatics Journal, 2, 112-116. https://doi.org/10.2174/1874431100802010112
[4]	Chootrakool, H., Shi, J.Q. and Yue, R. (2011) Meta-Analysis and Sensitivity Analysis for Multi-Arm Trials with Selection Bias. Statistics in Medicine, 30, 1183-1198. https://doi.org/10.1002/sim.4143
[5]	Rücker, G. (2017) Methods for Including Information from Multi-Arm Trials in Pairwise Meta-Analysis Multi-Arm Trials in Pairwise Meta-Analysis. Research Synthesis Methods, 8, 392-403. https://doi.org/10.1002/jrsm.1259
[6]	Lu, G. and Ades, A.E. (2004) Combination of Direct and Indirect Evidence in Mixed Treatment Comparisons. Statistics in Medicine, 23, 3105-3124. https://doi.org/10.1002/sim.1875
[7]	Lee, A.W. (2014) Review of Mixed Treatment Comparisons in Published Systematic Reviews Shows Marked Increase since 2009. Journal of Clinical Epidemiology, 67, 138-143.
[8]	Salanti, G. (2012) Indirect and Mixed-Treatment Comparison, Network, or Multiple-Treatments Meta-Analysis: Many Names, Many Benefits, Many Concerns for the Next Generation Evidence Synthesis Tool. Research Synthesis Methods, 3, 80-97. https://doi.org/10.1002/jrsm.1037
[9]	Lumley, T. (2002) Network Meta-Analysis for Indirect Treatment Comparisons. Statistics in Medicine, 21, 2313-2324. https://doi.org/10.1002/sim.1201
[10]	Mbuagbaw, L., et al. (2017) Approaches to Interpreting and Choosing the Best Treatments in Network Meta-Analyses. Systematic Reviews, 6, Article No. 79. https://doi.org/10.1186/s13643-017-0473-z
[11]	Veroniki, A.A., et al. (2018) Is Providing Uncertainty Intervals in Treatment Ranking Helpful in a Network Meta-Analysis? Journal of Clinical Epidemiology, 100, 122-129. https://doi.org/10.1016/j.jclinepi.2018.02.009
[12]	Seide, S.E., Jensen, K. and Kieser, M. (2019) Simulation and Data-Generation for Random-Effects Network Meta-Analysis of Binary Outcome. Statistics in Medicine, 38, 3288-3303. https://doi.org/10.1002/sim.8193
[13]	Jenkins, D.A., et al. (2021) Methods for the Inclusion of Real-World Evidence in Network Meta-Analysis. BMC Medical Research Methodology, 21, Article No. 207. https://doi.org/10.1186/s12874-021-01399-3
[14]	Zhao, Y.P., et al. (2021) A Network Meta-Analysis for Neoadjuvant and Adjuvant Treatments for Resectable Squamous Cell Carcinoma of Esophagus. Scientific Reports, 11, Article No. 6800. https://doi.org/10.1038/s41598-021-86102-8
[15]	Cox, D.R. and Snell, E.J. (1989) Analysis of Binary Data. 2nd Edition, Chapman & Hall/CRC, Boca Raton.
[16]	Tritchler, D. (1984) An Algorithm for Exact Logistic Regression. Journal of the American Statistical Association, 79, 709-711. https://doi.org/10.2307/2288420
[17]	Sartori, N. (2003) Modified Profile Likelihoods in Models with Stratum Nuisance Parameters. Biometrika, 90, 533-549. http://www.jstor.org/stable/30042064 https://doi.org/10.1093/biomet/90.3.533
[18]	Lee, W.J., Shi, J.Q. and Lee, Y.J. (2010) Approximate Conditional Inference in Mixed-Effects Models with Binary Data. Computational Statistics & Data Analysis, 54, 173-184. https://doi.org/10.1016/j.csda.2009.07.027
[19]	Di Caterina, C., Cortese, G. and Sartori, N. (2018) Monte Carlo Modified Profile Likelihood for Panel Data Models. arXiv:1801.02597v1.
[20]	Lu, G.B. and Ades, A.E. (2006) Assessing Evidence Inconsistency in Mixed Treatment Comparisons. Journal of the American Statistical Association, 101, 447-459. https://doi.org/10.1198/016214505000001302
[21]	Weber, F., et al. (2020) Zero-Cell Corrections in Random-Effects Meta-Analyses. Research Synthesis Methods, 11, 913-919. https://doi.org/10.1002/jrsm.1460
[22]	Greco, T., et al. (2015) The Attractiveness of Network Meta-Analysis: A Comprehensive Systematic and Narrative Review. Heart Lung and Vessels, 7, 133-142. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4476767/
[23]	Wang, F.F., et al. (2018) Pharmacologic Therapy to Induce Weight Loss in Women Who Have Obesity/Overweight with Polycystic Ovary Syndrome: A Systematic Review and Net-Work Meta-Analysis. Obesity Reviews, 19, 1424-1445. https://doi.org/10.1111/obr.12720
[24]	Zhu, L.K., et al. (2019) A New Binomial Tree Method for European Options under the Jump Diffusion Model. Journal of Applied Mathematics and Physics, 7, 3012-3021. https://doi.org/10.4236/jamp.2019.712211
[25]	Kleinbaum, D.G. (1994) Logistic Regression. Springer, New York. https://doi.org/10.1007/978-1-4757-4108-7
[26]	Collett, D. (1991) Modelling Binary Data (Hardback). Taylor & Francis, Oxfordshire.
[27]	Severini, T.A. (1998) Likelihood Functions for Inference in the Presence of a Nuisance Parameter. Biometrika, 85, 507-522. https://doi.org/10.1093/biomet/85.3.507
[28]	Bellio, R. and Sartori, N. (2003) Extending Conditional Likelihood in Models for Stratified Binary Data. Statistical Methods and Applications, 12, 121-132. https://doi.org/10.1007/s10260-003-0055-1
[29]	Hamza, T.H., van Houwelingen, H.C. and Stijnen, T. (2008) The Binomial Distribution of Meta-Analysis Was Preferred to Model Within-Study Variability. Journal of Clinical Epidemiology, 61, 41-51. https://doi.org/10.1016/j.jclinepi.2007.03.016
[30]	Adeniran, A.T., et al. (2020) Derivation of Gaussian Probability Distribution: A New Approach. Applied Mathematics, 11, 436-446. https://doi.org/10.4236/am.2020.116031
[31]	Crouch, E.A.C. and Spiegelman, D. (1990) The Evaluation of Integrals of the Form ∫+∞ -∞ f(t)exp(-t2) dt: Application to Logistic-Normal Models. Journal of the American Statistical Association, 85, 464-469. https://doi.org/10.1080/01621459.1990.10476222
[32]	Shi, J.Q. and Copas, J. (2002) Publication Bias and Meta-Analysis for 2 × 2 Tables: An Average Markov Chain Monte Carlo EM Algorithm. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 64, 221-236. https://doi.org/10.1111/1467-9868.00334
[33]	Liu, Q. and Pierce, D.A. (1994) A Note on Gauss-Hermite Quadrature. Biometrika, 81, 624-629. https://doi.org/10.2307/2337136
[34]	Shaobo, J. and Björn, A. (2020) A Note on the Accuracy of Adaptive Gauss-Hermite Quadrature. Biometrika, 107, 737-744. https://doi.org/10.1093/biomet/asz080