The R Package metaLik for Likelihood Inference in Meta-Analysis

Meta-analysis is a statistical method for combining information from different studies about the same issue of interest. Meta-analysis is widely diffuse in medical investigation and more recently it received a growing interest also in social disciplines. Typical applications involve a small number of studies, thus making ordinary inferential methods based on first-order asymptotics unreliable. More accurate results can be obtained by exploiting the theory of higher-order asymptotics. This paper describes the metaLik package which provides an R implementation of higher-order likelihood methods in meta-analysis. The extension to meta-regression is included. Two real data examples are used to illustrate the capabilities of the package.


Introduction
Meta-analysis is a statistical method for pooling the results from multiple separate studies about the same issue of interest.It has a considerable impact on medical and epidemiological investigation (Sutton, Jones, Abrams, Sheldon, and Song 2000), as testified by thousands of papers in specialized journals.More recently, the scope of application of meta-analysis reached other fields, such as sociology and economics (Roberts 2005;Sutton and Higgins 2008).
The traditional approach to meta-analysis exploits a random-effect formulation, with parameters estimated according to the method-of-moments procedure due to DerSimonian and Laird (1986).Despite its feasibility, the approach is prone to substantial disadvantages, giving rise to unreliable inferential results.These can be experienced mainly in case of small sample size.Several alternatives have been proposed to improve the results, see Van Houwelingen, Arends, and Stijnen (2002) for a review.Some authors suggest to resort to a likelihood approach, e.g., Hardy and Thompson (1996).However, standard first-order likelihood approximations can be inaccurate because of the typical small sample sizes of meta-analysis.To overcome such difficulties, Guolo (2012) shows that accurate inferential conclusions can be restored by exploiting the theory of higher-order asymptotics.Some general references about higher-order procedures are Severini (2000), Skovgaard (2001), and Reid (2003).
There exist already several packages for meta-analysis in the R system for statistical computing (R Development Core Team 2012).Packages meta (Schwarzer 2012) and rmeta (Lumley 2009) allow fitting fixed-and random-effects models for meta-analysis by relying on DerSimonian and Laird's procedure.An extension of the traditional analysis is provided by metafor (Viechtbauer 2010), which includes meta-regression problems and, in the meanwhile, considers approaches different from those by DerSimonian and Laird for estimating the betweenstudy heterogeneity.Maximum likelihood estimation in multivariate meta-analysis and metaregression is implemented in mvmeta (Gasparrini 2012).The package metatest (Huizenga, Visser, and Dolan 2011) is addressed to perform hypothesis testing via a variety of procedures including likelihood ratio test with Bartlett correction and permutations.
The package metaLik described in this paper extends the likelihood approach to meta-analysis and meta-regression to guarantee higher accuracy of the asymptotic results, which is better appreciable in case of small sample sizes.In particular, the package implements the secondorder likelihood method based on Skovgaard's statistic (Skovgaard 1996).The package met-aLik is available from the Comprehensive R Archive Network at http://CRAN.R-project.org/package=metaLik.After a brief summary of the meta-analysis problem in Section 2, the paper describes the theory underlying metaLik in Section 3. Section 4 contains a detailed description of the functionalities of metaLik, while its application to real data sets is illustrated in Section 5.The paper ends with final remarks and discussion about possible future developments in Section 6.

Meta-analysis
The main aim of meta-analysis is to perform inference on a true effect β, on the basis of summary information obtained from K studies.Let Y i be the summary measure of β provided by the i-th study, i = 1, . . ., K, such as the log odds ratio.Let σ 2 i denote the associated measure of the uncertainty.The common meta-analysis model is a linear random-effects model (DerSimonian and Laird 1986) with e i and ε i assumed to be independent.The variance component τ 2 accounts for the possibility of heterogeneity among the studies.The usual assumption is that the sample size of each study is sufficiently large to judge the within-study variance σ 2 i as known and equal to the variance reported in each study, σ2 i .The most common inference method in meta-analysis is due to DerSimonian and Laird (1986).This is a two-steps method-of-moments estimation procedure.First, τ 2 is estimated by τ 2 = t I (0,∞) (t), where I A (z) is the indicator of the event {z ∈ A} and In the expression above, q w is the observed value of i=1 w i being the estimator of β in the fixed-effects model obtained as a weighted average of Y i with optimal weights w i .Then, the interest parameter β is estimated by where w i (τ 2 ) = 1/(w −1 i + τ 2 ) and w i = 1/σ 2 i .The variance of the estimator is VAR( βDL ) = 1/ K i=1 w i (τ 2 ).This method suffers two inconveniences.First, it provides a biased estimate of τ 2 .Second, inference on β does not account for the uncertainty in estimating τ 2 .As a consequence, confidence intervals for β are usually narrower on average than they should be.
In order to take account of possible sources of heterogeneity, the meta-analysis model can be extended to include covariates at the study level.Let x i denote the vector of p covariates available at the aggregated meta-analysis level for each study and let β be the associated p-dimensional vector of effects.The meta-regression model is (Thompson and Higgins 2002;Knapp and Hartung 2003) Y i ∼ Normal x i β, σ2 i + τ 2 , which reduces to the meta-analysis model with scalar β when x i is one.DerSimonian and Laird's estimator (1) has a straightforward extension to the meta-regression case, e.g., Knapp and Hartung (2003).
Despite the feasibility of DerSimonian and Laird's approach, several authors emphasize that its application gives rise to unreliable inferential results, mainly in case of small sample sizes.This partly explains the proposal of alternatives in recent years, both from a frequentist and a Bayesian perspective, see Van Houwelingen et al. (2002) for a review.

Likelihood methods in meta-analysis
The log-likelihood function for the whole parameter vector ψ = (β , τ 2 ) based on a sample of size K and conditionally on the observed within-study variances is The maximum likelihood estimate ψ = ( β , τ 2 ) can be obtained iteratively (Brockwell and Gordon 2001).Suppose that ψ is partitioned into a component of interest θ and a nuisance component λ.In meta-analysis and meta-regression θ is usually a component of or the whole fixed-effects vector β, while the between-study variance τ 2 is a nuisance component.We denote by ψ = ( θ , λ ) the unconstrained maximum likelihood estimate of ψ and by ψ = (θ , λ θ ) the constrained maximum likelihood estimate of ψ for a fixed value of θ.Thus, inference on θ can rely on the profile log-likelihood P ( ψ) = (θ; λθ ), e.g., Hardy and Thompson (1996).For scalar θ, inference can be performed through the signed profile log-likelihood ratio Under mild regularity conditions, r P (θ) is asymptotically standard normally distributed up to an error of order O(n −1/2 ), see Severini (2000, Section 4.4).First-order likelihood results can be inaccurate in many applications, mainly in the case of small sample sizes, as it typically occurs in meta-analysis.A valid solution is resorting to the theory of higher-order asymptotics, see Severini (2000, Chapter 7).We refer in this paper to the modified version of r P (θ) suggested by Skovgaard (Skovgaard 1996), This term involves the expected information matrix Ĵ and the observed information matrix Î, both evaluated at the maximum likelihood estimate ψ, and the subblock J λλ of J corresponding to the parameter λ, evaluated at the constrained maximum likelihood estimate ψ.
The term [S −1 q] θ represents the component of the vector S −1 q corresponding to θ, with S and q being covariances of likelihood quantities, where ∇ z f (z) denotes the gradient of f (z) with respect to z.Despite its apparent complexity, Skovgaard's statistic r P (θ) can be rather easily derived, since it involves components whose evaluation is similar to that required by the expected information matrix.Skovgaard's statistic r P (θ) has still a standard normal approximate distribution as r P (θ), but with a higher-order accuracy.In fact, the standard normal approximation is up to an error of order O(n −1 ) for moderate deviations and O(n −1/2 ) for large deviations.Moreover, Skovgaard's statistic achieves a third-order accuracy in a full exponential family.This particular situation corresponds to the meta-analysis with equal within-study variances σ2 i = σ2 , i = 1, . . ., K. See Guolo (2012) for further details about how to derive Skovgaard's statistic components in meta-analysis and meta-regression.

Score test for heterogeneity
If the estimate of the heterogeneity parameter τ 2 is not significantly different from zero, then the random-effects model reduces to the fixed-effects model, where first-order results are exact under model assumptions, making higher-order adjustments unnecessary.The assessment of the significance of τ 2 estimate entails the one-sided test Since the null model is embedded in the boundary of the parameter space, classical inference results do not hold.For example, the likelihood-based test statistics do not follow the usual limiting χ 2 1 distribution.A practical solution is to evaluate the test statistic under the null hypothesis via parametric bootstrap.Sinha (2009) applies this idea to the score statistic which is particularly convenient since it does not require maximum likelihood estimation within each bootstrap sample.The algorithm is summarized below.Let ψ0 = ( β 0 , 0) be the maximum likelihood estimate of ψ under the null hypothesis and let .
be the score statistic.Hence, 1. Compute the observed value of the score statistic W obs .

The metaLik package
The main function in metaLik is metaLik(), which implements first-order and higher-order likelihood methods for inference in meta-analysis and meta-regression models described in Guolo (2012).The arguments of metaLik() are metaLik(formula, data, subset, contrasts = NULL, offset, sigma2, weights = 1/sigma2) This function uses standard model-frame specification via formula, data, subset, contrasts, and offset, following Chambers and Hastie (1992).A typical model has the form y ~x1 + x2 + ... + xJ, where y is a continuous response term and xj is the j-th covariate available at the aggregated meta-analysis level for each study, j = 1, . . ., p.If a model with no covariates (only intercept) is specified, then the classical meta-analysis model is fitted, by typing y~1.
An unusual argument of metaLik() with respect to the linear regression model specification of lm() is sigma2, which is needed to specify the vector of estimated within-study variances σ2 i .Alternatively, the inverse quantities weights can be specified.
The maximum likelihood estimation of the model parameters needs numerical optimization.
To this aim, we chose the "BFGS" algorithm implemented in optim(), with the analytic gradient of the log-likelihood provided.As pointed out by Cribari-Neto and Zeileis (2010), "BFGS" is generally thought to be the best performing quasi-Newton method.Starting values for the optimization procedure are the DerSimonian and Laird's estimates.The score test for the significance of the maximum likelihood estimate of the heterogeneity parameter τ 2 described in Section 3 is implemented with B = 1, 000 bootstrap samples.
Function metaLik() returns a fitted-object of class metaLik.This is a list structured similarly to other R regression functions, as for example glm().There are standard methods for objects of class metaLik that can be used to extract quantities of interest.They are coef() and vcov() to extract maximum likelihood estimates of fixed-effects parameters and their estimated variance-covariance matrix, confint() to obtain confidence intervals, logLik() to extract the value of the maximized log-likelihood and residuals() to compute Pearson residuals.
Method summary() provides summary information about the maximum likelihood parameter estimates, either for the fixed-effects components and the heterogeneity parameter, their associated standard errors, and first-and higher-order log-likelihood ratio statistics about the significance of the fixed-effects components.In particular, the values of r P and r P are reported, together with the associated p value.Method summary() internally calls function test.metaLik(),which performs hypothesis testing on a scalar component of the fixed-effects vector in meta-analysis and meta-regression models, using the signed profile log-likelihood ratio test and its higher-order Skovgaard's adjustment.In case the bootstrap score test indicates that the maximum likelihood estimate of τ 2 is not significant at 5% level, then τ 2 is set to zero and Skovgaard's adjustment vanishes.Function test.metaLik() has a syntax similar to that of other R functions for significance tests, test.metaLik(object,param = 1, value = 0, alternative = c("two.sided","less", "greater"), print = TRUE) The arguments of test.metaLik() are a fitted metaLik object (object), the index or the name of the parameter to test (param), the value of the parameter under the null hypothesis (value) and whether the user is interested in a two-sided or an one-sided alternative hypothesis (alternative).The output is a list of information about the value of the first-order statistic and Skovgaard's statistic, together with the associated p value.
Function test.metaLik() is also used by the profile() method, which computes and plots confidence intervals for a scalar component of the fixed-effects vector profile(fitted, param = 1, level = 0.95, display = TRUE, ...) This function returns a matrix with columns the endpoints of the confidence intervals at level for the specified parameter according to r P and r P .If display = TRUE, the signed square root of the profile likelihood is plotted together with the horizontal dashed lines relative to the specified level and the corresponding confidence intervals.

Illustrations
In this section, we illustrate the application of the metaLik package on two data sets available through the package.

Diuretics data
We start by analyzing data about prevention of pre-eclampsia with diuretics previously analyzed by Biggerstaff and Tweedie (1997).The data consists of the logarithm of the risk ratio of eclampsia in nine randomized trials, see Table 1.
A linear random-effects meta-analysis model can be fitted through function metaLik() R> library("metaLik") R> data("diuretics") The first-order likelihood approach (signed logLRT) suggests a significant effect of diuretics at 5% level, with an associated p value equal to 0.035.The Skovgaard's statistic (Skovgaard) instead yields a slightly larger p value 0.065, thus making declaration of significance doubtful.
Confidence interval for the model parameter β based on the Wald statistic is obtained by the standard method confint() 2.5 % 97.5 % (Intercept) -0.922 -0.113 while profile likelihood confidence intervals are provided by function profile()

R> profile(m)
Confidence interval for parameter (Intercept) 2.5 % 97.5 % signed logLRT -0.983 -0.0484 Skovgaard -1.064 0.0433 Signed square root profile likelihoods are visualized in Figure 1.The plot shows the small but significant difference since the first-order 95% confidence interval does not include the zero, while the higher-order version does.
Since the inferential interest is evaluating the reduction of the risk of pre-eclampsia, then it is more appropriate to test the null hypothesis H 0 : β = 0 against the one-sided alternative H 1 : β < 0. Hypothesis testing can be performed using r P or r P , through function test.metaLik()R> test.metaLik(m,value = 0, alternative = "less") Signed profile log-likelihood ratio test for parameter (Intercept) First-order statistic r:-2.11,p-value:0.0174Skovgaard's statistic rSkov:-1.85,p-value:0.0325alternative hypothesis: parameter is less than 0 Both the approaches indicate an effect of diuretics in reducing the risk of pre-eclampsia, although the level of support is sensibly different using first-(p value 0.0174) or higher-order (p value 0.0325) asymptotics.

BCG vaccine data
The second illustration regards a data set on tuberculosis prevention originally examined by Berkey, Hoaglin, Mosteller, and Colditz (1995) and Knapp and Hartung (2003) (Berkey et al. 1995).
others.The data refer to thirteen clinical studies evaluating the efficacy of the Bacillus Calmette-Guérin (BCG) vaccine for the prevention of tuberculosis.Data contain information about the log odds ratio in each study and the corresponding distance of each study from the equator, i.e., the latitude, which is considered as a surrogate for the presence of environmental mycobacteria providing a level of natural immunity against tuberculosis.Moreover, the year of the study and the estimated within-study variance are available, see Table 2.
We fit a linear mixed-effects model with covariate latitude, R> library("metaLik") R> data("vaccine") R> m <-metaLik(y ~latitude, data = vaccine, sigma2 = sigma2) R> summary(m) Likelihood inference in random-effects meta-analysis models Call: metaLik(formula = y ~latitude, data = vaccine, sigma2 = sigma2) Estimated heterogeneity parameter tau^2: 0. The first-order likelihood approach (signed logLRT) suggests a significant effect of the latitude on the disease risk at 5% level, with an associated p value equal to 0.034.The Skov-gaard's statistic (Skovgaard) does not fully support the same conclusion given the p value equal to 0.069.
Confidence intervals for the model parameters based on the Wald statistic are R> confint(m) 2.5 % 97.5 % (Intercept) -0.7442 0.13414 latitude -0.0279 -0.00294 The profile likelihood 95% confidence intervals for the parameter of interest latitude with first-and higher-order approximations are R> profile(m, param = "latitude") Confidence interval for parameter latitude 2.5 % 97.5 % signed logLRT -0.0290 -0.00141 Skovgaard -0.0315 0.00146 The profile likelihood confidence intervals with the two approximations are also displayed in the left panel of Figure 2, further illustrating how the improved accuracy of higher-order asymptotics leads to different conclusions about the parameter significance at the 95% level.
As one of the referees pointed out, an alternative analysis of these data substitutes the covariate latitude with its absolute value, so to highlight the distance of each study from the equator R> m2 <-metaLik(y ~abs(latitude), data = vaccine, sigma2 = sigma2) R> summary(m2)  model.In this case, results from metaLik refer only to the ordinary signed log-likelihood ratio statistic, which provides a strong indication of negative association between the distance from equator and the effectiveness of the vaccine.

Conclusions
Higher-order asymptotics literature has had an important impact on methodological journals in the last years, as the detailed review by Reid (2003) pointed out.The advantages of higher-order solutions methods are shown to be substantial with respect to their first-order counterparts, mainly in case of small sample sizes.Nevertheless, the application of the methods in practical contexts is still limited, the reason being the complexity that the expression of the higher-order solutions usually achieves.The book by Brazzale, Davison, and Reid (2007) illustrates some examples and case studies where the application of higher-order solutions is successful.Guolo (2012) shows that the meta-analysis problem is well-suited to deal with higher-order asymptotics, since the resulting statistics have a simple and compact form, while providing substantial improvements over first-order counterparts, especially for small sample sizes.
In this paper, we addressed the R implementation of the methods as illustrated in Guolo (2012).The resulting package metaLik supplies the first implementation of higher-order solutions in meta-analysis, and in particular of the second-order accurate inference on a scalar component of the fixed effects through the Skovgaard's statistic.The attention has been paid to meta-analysis models where summary information is provided for each study.Future research and implementation will focus on the extension of the meta-analysis model to individual patients data, where data for each subject are used in place of the summary information provided at the study level.The analysis of multiple outcome data is a further issue of investigation.

3.
Compute the score statistic for each bootstrap data set, W (b) , b = 1, . . ., B. 4. Estimate the p value by the proportion of W (b) exceeding W obs

Figure 2 :
Figure 2: BCG vaccine data: Profile likelihoods and 95% confidence intervals for latitude covariate (left panel) and absolute value of latitude (right panel).

Table 2 :
Thirteen clinical studies for evaluating the efficacy of the Bacillus Calmette-Guérin vaccine for the prevention of tuberculosis The maximum likelihood estimate of the heterogeneity parameter τ 2 is 0.004.The p value of the bootstrap score test for significance, 0.968, is a strong indication in favor of a fixed-effects