meta.pen {altmeta} | R Documentation |
A penalization approach to random-effects meta-analysis
Description
Performs the penalization methods introduced in Wang et al. (2022) to achieve a compromise between the common-effect and random-effects model.
Usage
meta.pen(y, s2, data, tuning.para = "tau", upp = 1, n.cand = 100, tol = 1e-10)
Arguments
y |
a numeric vector or the corresponding column name in the argument data, specifying the observed effect sizes in the collected studies. |
s2 |
a numeric vector or the corresponding column name in the argument data, specifying the within-study variances. |
data |
an optional data frame containing the meta-analysis dataset. If |
tuning.para |
a character string specifying the type of tuning parameter used in the penalization method. It should be one of "lambda" (use lambda as a tuning parameter) or "tau" (use the standard deviation as a tuning parameter). The default is "tau". |
upp |
a positive scalar used to control the upper bound of the range for the tuning parameter. Specifically, [0, T* |
n.cand |
the total number of candidate values of the tuning parameter within the specified range. The default is 100. |
tol |
the desired accuracy (convergence tolerance). The default is 1e-10. |
Details
Suppose a meta-analysis collects n
independent studies. Let \mu_{i}
be the true effect size in study i
(i
= 1, ..., n
). Each study reports an estimate of the effect size and its sample variance, denoted by y_{i}
and s_{i}^{2}
, respectively. These data are commonly modeled as y_{i}
from N(\mu_{i},s_{i}^{2})
. If study-specific true effect sizes \mu_{i}
are assumed i.i.d. from N(\mu, \tau^{2})
, this is the random-effects (RE) model, where \mu
is the overall effect size and \tau^{2}
is the between-study variance. If \tau^{2}=0
and thus \mu_{i}=\mu
for all studies, this implies that studies are homogeneous and the RE model is reduced to the common-effect (CE) model.
Marginally, the RE model yields y_{i} \sim N(\mu,s_{i}^2+\tau^{2})
, and its log-likelihood is
l(\mu, \tau^{2}) = -\frac{1}{2} \sum_{i=1}^{n} \left[\log(s_{i}^{2} + \tau^{2}) + \frac{(y_{i} - \mu)^2}{s_{i}^2 + \tau^2}\right] + C,
where C
is a constant. In the past two decades, penalization methods have been rapidly developed for variable selection in high-dimensional data analysis to control model complexity and reduce the variance of parameter estimates. Borrowing the idea from the penalization methods, we employ a penalty term on the between-study variance \tau^2
when the heterogeneity is overestimated. The penalty term increases with \tau^2
. Specifically, we consider the following optimization problem:
(\hat{\mu}(\lambda), \hat{\tau}^2(\lambda)) = \min_{\mu, \tau^2 \geq 0} \left\{\sum_{i=1}^{n} \left[ \log(s_{i}^{2} + \tau^{2}) + \frac{(y_{i} - \mu)^2}{s_{i}^2 + \tau^2}\right] + \lambda \tau^2\right\},
where \lambda \geq 0
is a tuning parameter that controls the penalty strength. Using the technique of profile likelihood by taking the target function's derivative in the above equation with respect to \mu
for a given \tau^2
. The bivariate optimization problem is reduced to a univariate minimization problem. When \lambda=0
, the minimization problem is equivalent to minimizing the log-likelihood without penalty, so the penalized-likelihood method is identical to the conventional RE model. By contrast, it can be shown that a sufficiently large \lambda
produces the estimated between-study variance as 0, leading to the conventional CE model.
As different tuning parameters lead to different estimates of \hat{\mu}(\lambda)
and \hat{\tau}^2(\lambda)
, it is important to select the optimal \lambda
among a set of candidate values. We perform the cross-validation process and construct a loss function of \lambda
to measure the performance of specific \lambda
values. The \lambda
corresponding to the smallest loss is considered optimal. The threshold, denoted by \lambda_{\max}
, based on the penalty function p(\tau^2)=\tau^2
can be calculated. For all \lambda > \lambda_{\max}
, the estimated between-study variance is 0. Consequently, we select a certain number of candidate values (e.g., 100) from the range [0,\lambda_{\max}]
for the tuning parameter. For a set of tuning parameters, the leave-one-study-out (i.e., n
-fold) cross-validation is used to construct the loss function. Specifically, we use the following loss function for the penalization method by tuning \lambda
:
\hat{L}(\lambda) = \left[\frac{1}{n} \sum_{i=1}^{n} \frac{(y_{i} - \hat{\mu}_{(-i)}(\lambda))^2}{s_{i}^2 + \hat{\tau}_{RE(-i)}^2 + \frac{\sum_{j \ne i}(s_{j}^2 + \hat{\tau}_{RE(-i)}^2) / (s_{j}^2 + \hat{\tau}^2_{(-i)}(\lambda))^2}{(\sum_{j \ne i} 1 / (s_{j}^2 + \hat{\tau}^2_{(-i)}(\lambda)))^2}}\right]^{1/2},
where the subscript (-i)
indicates that study i
is removed, \hat{\tau}^2_{(-i)}(\lambda)
is the estimated between-study variance for a given \lambda
, \hat{\tau}_{RE(-i)}^2
is the corresponding between-study variance estimate of the RE model, and
\hat{\mu}_{(-i)}(\lambda) = \frac{\sum_{j \ne i} y_{j} / (s_{j}^2 + \hat{\tau}^2_{(-i)}(\lambda))}{\sum_{j \ne i} 1 / (s_{j}^2 + \hat{\tau}^2_{(-i)}(\lambda))}.
The above procedure focuses on tuning the parameter, \lambda
, to control the penalty strength for the between-study variance. Alternatively, for the purpose of shrinking the potentially overestimated heterogeneity, we can directly treat the between-study standard deviation (SD), \tau
, as the tuning parameter. A set of candidate values of \tau
are considered, and the value that produces the minimum loss function is selected. Compared with tuning \lambda
from the perspective of penalized likelihood, tuning \tau
is more straightforward and intuitive from the practical perspective. The candidate values of \tau
can be naturally chosen from [0,\hat{\tau}_{RE}]
, with the lower and upper bounds corresponding to the CE and RE models, respectively. Denoted the candidate SDs as \tau_{t}
, the loss function with respect to \tau_{t}
is similarly defined as
\hat{L}(\tau_{t}) = \left[\frac{1}{n} \sum_{i=1}^{n} \frac{(y_i - \hat{\mu}_{(-i)}(\tau_{t}))^2}{s_{i}^2 + \hat{\tau}_{RE(-i)} + \frac{\sum_{j \ne i}(s_{j}^2 + \hat{\tau}^2_{RE(-i)}) / (s_{j}^2 + \tau_{t}^2)^2}{[\sum_{j \ne i} 1 / (s_{j}^2 + \tau_{t}^2)]^2}}\right]^{1/2},
where the overall effect size estimate (excluding study i
) is
\hat{\mu}_{(-i)}(\tau_{t}) = \frac{\sum_{j \ne i} y_j / (s_{j}^2 + \tau_{t}^2)}{\sum_{j \ne i} 1 / (s_{j}^2 + \tau_{t}^2)}.
Value
This function returns a list containing estimates of the overall effect size and their 95% confidence intervals. Specifically, the components include:
n |
the number of studies in the meta-analysis. |
I2 |
the |
tau2.re |
the maximum likelihood estimate of the between-study variance. |
mu.fe |
the estimated overall effect size of the |
se.fe |
the standard deviation of the overall effect size estimate of the |
mu.re |
the estimated overall effect size of the |
se.re |
the standard deviation of the overall effect size estimate of the |
loss |
the values of the loss function for candidate tuning parameters. |
tau.cand |
the candidate values of the tuning parameter |
lambda.cand |
the candidate values of the tuning parameter |
tau.opt |
the estimated between-study standard deviation of the penalization method. |
mu.opt |
the estimated overall effect size of the penalization method. |
se.opt |
the standard error estimate of the estimated overall effect size of the penalization method. |
Author(s)
Yipeng Wang yipeng.wang1@ufl.edu
References
Wang Y, Lin L, Thompson CG, Chu H (2022). "A penalization approach to random-effects meta-analysis." Statistics in Medicine, 41(3), 500–516. <doi: 10.1002/sim.9261>
Examples
data("dat.bohren") ## log odds ratio
## perform the penaliztion method by tuning tau
out11 <- meta.pen(y, s2, dat.bohren)
## plot the loss function and candidate taus
plot(out11$tau.cand, out11$loss, xlab = NA, ylab = NA,
lwd = 1.5, type = "l", cex.lab = 0.8, cex.axis = 0.8)
title(xlab = expression(paste(tau[t])),
ylab = expression(paste("Loss function ", ~hat(L)(tau[t]))),
cex.lab = 0.8, line = 2)
idx <- which(out11$loss == min(out11$loss))
abline(v = out11$tau.cand[idx], col = "gray", lwd = 1.5, lty = 2)
## perform the penaliztion method by tuning lambda
out12 <- meta.pen(y, s2, dat.bohren, tuning.para = "lambda")
## plot the loss function and candidate lambdas
plot(log(out12$lambda.cand + 1), out12$loss, xlab = NA, ylab = NA,
lwd=1.5, type = "l", cex.lab = 0.8, cex.axis = 0.8)
title(xlab = expression(log(lambda + 1)),
ylab = expression(paste("Loss function ", ~hat(L)(lambda))),
cex.lab = 0.8, line = 2)
idx <- which(out12$loss == min(out12$loss))
abline(v = log(out12$lambda.cand[idx] + 1), col = "gray", lwd = 1.5, lty = 2)
data("dat.bjelakovic") ## log odds ratio
## perform the penaliztion method by tuning tau
out21 <- meta.pen(y, s2, dat.bjelakovic)
## plot the loss function and candidate taus
plot(out21$tau.cand, out21$loss, xlab = NA, ylab = NA,
lwd=1.5, type = "l", cex.lab = 0.8, cex.axis = 0.8)
title(xlab = expression(paste(tau[t])),
ylab = expression(paste("Loss function ", ~hat(L)(tau[t]))),
cex.lab = 0.8, line = 2)
idx <- which(out21$loss == min(out21$loss))
abline(v = out21$tau.cand[idx], col = "gray", lwd = 1.5, lty = 2)
out22 <- meta.pen(y, s2, dat.bjelakovic, tuning.para = "lambda")
data("dat.carless") ## log odds ratio
## perform the penaliztion method by tuning tau
out31 <- meta.pen(y, s2, dat.carless)
## plot the loss function and candidate taus
plot(out31$tau.cand, out31$loss, xlab = NA, ylab = NA,
lwd=1.5, type = "l", cex.lab = 0.8, cex.axis = 0.8)
title(xlab = expression(paste(tau[t])),
ylab = expression(paste("Loss function ", ~hat(L)(tau[t]))),
cex.lab = 0.8, line = 2)
idx <- which(out31$loss == min(out31$loss))
abline(v = out31$tau.cand[idx], col = "gray", lwd = 1.5, lty = 2)
out32 <- meta.pen(y, s2, dat.carless, tuning.para = "lambda")
data("dat.adams") ## mean difference
out41 <- meta.pen(y, s2, dat.adams)
## plot the loss function and candidate taus
plot(out41$tau.cand, out41$loss, xlab = NA, ylab = NA,
lwd=1.5, type = "l", cex.lab = 0.8, cex.axis = 0.8)
title(xlab = expression(paste(tau[t])),
ylab = expression(paste("Loss function ", ~hat(L)(tau[t]))),
cex.lab = 0.8, line = 2)
idx <- which(out41$loss == min(out41$loss))
abline(v = out41$tau.cand[idx], col = "gray", lwd = 1.5, lty = 2)
out42 <- meta.pen(y, s2, dat.adams, tuning.para = "lambda")