Nothing
### To add / change:
# Add plot support for P-FPV, D-FIV and D-FAV
# Add subsections for Details
#
#' Regression tests for detecting funnel plot asymmetry
#'
#' The presence of small-study effects is a common threat to systematic reviews and meta-analyses, especially when
#' it is due to publication bias, which occurs when small primary studies are more likely to be reported (published)
#' if their findings were positive. The presence of small-study effects can be verified by visual inspection of
#' the funnel plot, where for each included study of the meta-analysis, the estimate of the reported effect size is
#' depicted against a measure of precision or sample size.
#' The premise is that the scatter of plots should reflect a funnel shape, if small-study
#' effects do not exist. However, when small studies are predominately in one direction (usually the
#' direction of larger effect sizes), asymmetry will ensue.\cr \cr
#' The \code{\link{fat}} function implements several tests for detecting funnel plot asymmetry,
#' which can be used when the presence of between-study heterogeneity in treatment effect is relatively low.
#'
#' @param b Vector with the effect size of each study. Examples are log odds ratio, log hazards ratio,
#' log relative risk.
#' @param b.se Optional vector with the standard error of the effect size of each study
#' @param n.total Optional vector with the total sample size of each study
#' @param d.total Optional vector with the total number of observed events for each study
#' @param d1 Optional vector with the total number of observed events in the exposed groups
#' @param d2 Optional vector with the total number of observed events in the unexposed groups
#' @param method Method for testing funnel plot asymmetry, defaults to \code{"E-FIV"} (Egger's test with
#' multiplicative dispersion). Other options are \code{E-UW}, \code{M-FIV}, \code{M-FPV}, \code{D-FIV} and
#' \code{D-FAV}. More info in "Details"
#'
#' @details
#'
#' \subsection{Egger regression method}{
#' A common approach to test the presence of small-study effects is to
#' estimate a regression model where the standardized effect estimate
#' (effect/SE) is regressed on a measure of precision (1/SE),
#' (\code{method="E-UW"}, Egger 1997).
#' It is possible to allow for between-study heterogeneity by adopting a
#' multiplicative overdispersion parameter by which the variance in each
#' study is multiplied (\code{method="E-FIV"}, Sterne 2000).
#'
#' Unfortunately, it has been demonstrated that the aforementioned two tests
#' are biased because: (i) the independent variable is subject to sampling
#' variability; (ii) the standardized treatment effect is correlated with its
#' estimated precision; and (iii) for binary data, the independent regression
#' variable is a biased estimate of the true precision, with larger bias for
#' smaller sample sizes (Macaskill et al. 2001).
#' }
#'
#' \subsection{Macaskill regression method}{
#' The standard approach estimates a regression model with the effect size as
#' a function of the study size (\code{method="M-FIV"}, Macaskill et al.
#' 2001). Each study is weighted by the precision of the treatment effect
#' estimate to allow for possible heteroscedasticity.
#' An alternative approach is to weight each study by a pooled' estimate of the
#' outcome proportion (\code{method="M-FPV"})
#'
#' For studies with zero events, a continuity correction is applied by adding
#' 0.5 to all cell counts.
#' }
#'
#' \subsection{Peters regression method}{
#' This approach (\code{method="P-FPV"}) estimates a regression model with the
#' treatment effect as a function of the inverse of the total sample size
#' (Peters et al. 2006).
#'
#' For studies with zero events, a continuity correction is applied by adding
#' 0.5 to all cell counts.
#' }
#'
#' \subsection{Debray regression method}{
#' This approach was proposed for survival data, and uses the total
#' number of events as independent variable in the weighted regression model
#' (Debray et al. 2017). The study weights are based on the inverse variance
#' (\code{method="D-FIV"}) or on an approximation thereof
#' (\code{method="D-FAV"}).
#' }
#'
#' @return a list containing the following entries:
#' \describe{
##' \item{"pval"}{A two-sided P-value indicating statistical significance of the funnel plot asymettry test.
##' Values below the significance level (usually defined as 10\%) support the presence of funnel plot asymmetry,
##' and thus small-study effects. }
##' \item{"model"}{A fitted \code{glm} object, representing the estimated regression model used for testing funnel
##' plot asymmetry.}
##' }
#' @author Thomas Debray <thomas.debray@gmail.com>
#'
#' @references Debray TPA, Moons KGM, Riley RD. Detecting small-study effects and funnel plot asymmetry in meta-analysis of
#' survival data: a comparison of new and existing tests. Res Syn Meth. 2018;9(1):41--50.\cr
#' \cr
#' Egger M, Davey Smith G, Schneider M, Minder C. Bias in meta-analysis detected by a simple, graphical test.
#' \emph{BMJ}. 1997;315(7109):629--34. \cr
#' \cr
#' Macaskill P, Walter SD, Irwig L. A comparison of methods to detect publication bias in meta-analysis.
#' \emph{Stat Med}. 2001;20(4):641--54.\cr
#' \cr
#' Peters JL, Sutton AJ, Jones DR, Abrams KR, Rushton L. Comparison of two methods to detect publication bias
#' in meta-analysis. \emph{JAMA}. 2006 Feb 8;295(6):676--80.\cr
#' \cr
#' Sterne JA, Gavaghan D, Egger M. Publication and related bias in meta-analysis: power of statistical tests
#' and prevalence in the literature. \emph{J Clin Epidemiol}. 2000;53(11):1119--29.
#'
#' @seealso \code{\link{plot.fat}}
#'
#' @examples
#' data(Fibrinogen)
#' b <- log(Fibrinogen$HR)
#' b.se <- ((log(Fibrinogen$HR.975) - log(Fibrinogen$HR.025))/(2*qnorm(0.975)))
#' n.total <- Fibrinogen$N.total
#' d.total <- Fibrinogen$N.events
#'
#' fat(b=b, b.se=b.se)
#' fat(b=b, b.se=b.se, d.total=d.total, method="D-FIV")
#'
#' # Note that many tests are also available via metafor
#' require(metafor)
#' fat(b=b, b.se=b.se, n.total=n.total, method="M-FIV")
#' regtest(x=b, sei=b.se, ni=n.total, model="lm", predictor="ni")
#'
#' @import stats
#' @importFrom stats pt qnorm
#' @importFrom metafor rma
#' @importFrom plyr round_any
#'
#' @export
fat <- function(b, b.se, n.total, d.total, d1, d2, method="E-FIV")
{
if (missing(b)) {
stop ("No values given for 'b'")
}
# Identify studies with complete information
if (method == "E-UW") {
if (missing(b.se)) {
stop ("No values given for 'b.se'")
}
if (length(b) != length(b.se)) {
stop("Incompatible vector sizes for 'b' and 'b.se'!")
}
studies.complete <- c(!is.na(b) & !is.na(b.se))
ds <- data.frame("y" = b,
"x" = b.se
)
} else if (method== "E-FIV") {
if (missing(b.se)) {
stop ("No values given for 'b.se'")
}
if (length(b) != length(b.se)) {
stop("Incompatible vector sizes for 'b' and 'b.se'!")
}
studies.complete <- c(!is.na(b) & !is.na(b.se))
ds <- data.frame("y" = b,
"x" = b.se,
"w" = (1/(b.se**2))
)
} else if (method == "M-FIV") {
if (missing(b.se)) {
stop ("No values given for 'b.se'")
}
if (missing(n.total)) {
stop ("No values given for 'n.total'")
}
if (length(b) != length(b.se)) {
stop("Incompatible vector sizes for 'b' and 'b.se'!")
}
if (length(b) != length(n.total)) {
stop("Incompatible vector sizes for 'b' and 'n.total'!")
}
studies.complete <- c(!is.na(b) & !is.na(b.se) & !is.na(n.total))
ds <- data.frame("y" = b,
"x" = n.total,
"w" = (1/(b.se**2))
)
} else if (method=="M-FPV") {
if (missing(n.total)) {
stop ("No values given for 'n.total'")
}
if (missing(d.total)) {
stop ("No values given for 'd.total'")
}
if (length(b) != length(n.total)) {
stop("Incompatible vector sizes for 'b' and 'n.total'!")
}
if (length(b) != length(d.total)) {
stop("Incompatible vector sizes for 'b' and 'd.total'!")
}
studies.complete <- c(!is.na(b) & !is.na(d.total) & !is.na(n.total))
# Consider continuity corrections
d.total.cc <- d.total
d.total.cc[d.total==0] <- 1 #0.5 event in exposed group and 0.5 event in non-exposed group
n.total[d.total==0] <- n.total[d.total==0]+2 #2*0.5 in the events, and 2*0.5 in the non-events
ds <- as.data.frame(cbind(b, n.total, (d.total.cc*(1-d.total.cc/n.total))))
colnames(ds) <- c("y","x","w")
} else if (method=="P-FPV") {
if (missing(n.total)) {
stop ("No values given for 'n.total'")
}
if (missing(d.total)) {
stop ("No values given for 'd.total'")
}
if (length(b) != length(n.total)) {
stop("Incompatible vector sizes for 'b' and 'n.total'!")
}
if (length(b) != length(d.total)) {
stop("Incompatible vector sizes for 'b' and 'd.total'!")
}
studies.complete <- c(!is.na(b) & !is.na(d.total) & !is.na(n.total))
# Consider continuity corrections
d.total.cc <- d.total
d.total.cc[d.total==0] <- 1 #0.5 event in exposed group and 0.5 event in non-exposed group
n.total[d.total==0] <- n.total[d.total==0]+2 #2*0.5 in the events, and 2*0.5 in the non-events
ds <- data.frame("y" = b,
"x" = 1/n.total,
"w" = (d.total.cc*(1-d.total.cc/n.total))
)
} else if (method=="D-FIV") {
if (missing(b.se)) {
stop ("No values given for 'b.se'")
}
if (missing(d.total)) {
stop ("No values given for 'd.total'")
}
if (length(b) != length(b.se)) {
stop("Incompatible vector sizes for 'b' and 'b.se'!")
}
if (length(b) != length(d.total)) {
stop("Incompatible vector sizes for 'b' and 'd.total'!")
}
studies.complete <- c(!is.na(b) & !is.na(b.se) & !is.na(d.total))
# Consider continuity corrections
d.total.cc <- d.total
d.total.cc[d.total==0] <- 1 #0.5 event in exposed group and 0.5 event in non-exposed group
ds <- data.frame("y" = b,
"x" = 1/d.total.cc,
"w" = (1/(b.se**2))
)
} else if (method=="D-FAV") {
if (missing(d1)) {
stop ("No values given for 'd1'")
}
if (missing(d2)) {
stop ("No values given for 'd2'")
}
if (length(b) != length(d1)) {
stop("Incompatible vector sizes for 'b' and 'd1'!")
}
if (length(b) != length(d2)) {
stop("Incompatible vector sizes for 'b' and 'd2'!")
}
if (!missing(d.total)) {
if (sum(d1+d2!=d.total) > 0)
stop("Incompatible information between 'd.total', 'd1' and 'd2'")
}
studies.complete <- c(!is.na(b) & !is.na(d1) & !is.na(d2))
# Consider continuity corrections
d1.cc <- d1
d2.cc <- d2
d1.cc[(d1==0 | d2==0)] <- d1.cc[(d1==0 | d2==0)]+0.5 #0.5 event in exposed group and 0.5 event in non-exposed group
d2.cc[(d1==0 | d2==0)] <- d2.cc[(d1==0 | d2==0)]+0.5
ds <- data.frame("y" = b,
"x" = 1/(d1.cc+d2.cc),
"w" = 1/((1/d1.cc)+(1/d2.cc))
)
}
else {
stop("Method for testing funnel plot asymmetry not supported")
}
# Identify which studies can be used
nstudies <- sum(studies.complete)
# Omit sudies with missing information
ds <- ds[studies.complete,]
if (nstudies < length(studies.complete)) {
warning("Some studies were removed due to missing data!")
}
# Get the fixed effect summary estimate
res <- NULL
if (!missing(b.se)) {
res <- rma(yi = b[studies.complete], sei = b.se[studies.complete], method = "FE")
}
if (method %in% c("E-FIV", "M-FIV", "M-FPV", "P-FPV", "D-FIV", "D-FAV")) {
suppressWarnings(m.fat <- try(glm(y~x, weights=ds$w, data=ds), silent=T))
} else if (method=="E-UW") {
suppressWarnings(m.fat <- try(glm(y~x, data=ds), silent=T))
} else {
stop("Method for testing funnel plot asymmetry currently not implemented")
}
if ("try-error" %in% attr(m.fat,"class")) {
warning("Estimation of the regression model unsuccessful, P-value omitted.")
t.fat <- NA
p.fat <- NA
} else {
t.fat <- coefficients(m.fat)[2]/sqrt(diag(vcov(m.fat))[2])
p.fat <- 2*pt(-abs(t.fat),df=(nstudies-2))
}
out <- list()
out$call <- match.call()
out$method <- method
out$tval <- t.fat
out$pval <- p.fat
out$fema <- res
out$df <- nstudies-2
out$model <- m.fat
class(out) <- "fat"
return(out)
}
#' @method print fat
#' @export
print.fat <- function(x, digits = max(3, getOption("digits") - 3), ...) {
cat("Call: ");
print(x$call);
if (!is.null(x$fema)) {
cat(c("\nFixed effect summary estimate: ", round(x$fema$b, digits = digits), " \n"))
}
cat("\n")
cat(paste("test for funnel plot asymmetry: t =", round(x$tval, digits = digits), ", df = ", x$df, ", ", sep=""))
cat(paste("p = ", round(x$pval, digits = digits), "\n", sep=""))
}
#' Display results from the funnel plot asymmetry test
#'
#' Generates a funnel plot for a fitted \code{fat} object.
#' @param x An object of class \code{fat}
#' @param ref A numeric value indicating the fixed or random effects summary estimate. If no value is provided
#' then it will be retrieved from a fixed effects meta-analysis (if possible).
#' @param confint A logical indicator. If \code{TRUE}, a confidence interval will be displayed for the estimated
#' regression model (based on a Student-T distribution)
#' @param confint.level Significance level for constructing the confidence interval.
#' @param confint.alpha A numeric value between 0 and 1 indicating the opacity for the confidence region.
#' @param confint.col The color for filling the confidence interval. Choose \code{NA} to leave polygons unfilled.
#' If \code{confint.density} is specified with a positive value this gives the color of the shading lines.
#' @param confint.density The density of shading lines, in lines per inch. The default value of \code{NULL} means
#' that no shading lines are drawn. A zero value of density means no shading nor filling whereas negative values
#' and \code{NA} suppress shading (and so allow color filling).
#' @param xlab A title for the x axis
#' @param add.pval Logical to indicate whether a P-value should be added to the plot
#' @param ... Additional arguments.
#'
#' @examples
#' data(Fibrinogen)
#' b <- log(Fibrinogen$HR)
#' b.se <- ((log(Fibrinogen$HR.975) - log(Fibrinogen$HR.025))/(2*qnorm(0.975)))
#' n.total <- Fibrinogen$N.total
#'
#' # A very simple funnel plot
#' plot(fat(b=b, b.se=b.se), xlab = "Log hazard ratio")
#'
#' # Plot the funnel for an alternative test
#' plot(fat(b=b, b.se=b.se, n.total=n.total, method="M-FIV"), xlab = "Log hazard ratio")
#'
#'
#' @import ggplot2
#' @importFrom stats qt
#' @importFrom graphics plot axis polygon points lines box abline
#'
#' @author Thomas Debray <thomas.debray@gmail.com>
#' @author Frantisek Bartos <f.bartos96@gmail.com>
#' @method plot fat
#'
#' @export
plot.fat <- function(x, ref, confint = TRUE, confint.level = 0.10, confint.col = "skyblue", confint.alpha = .50,
confint.density = NULL,
xlab = "Effect size", add.pval = TRUE, ...) {
if (!inherits(x, "fat"))
stop("Argument 'x' must be an object of class \"fat\".")
if (confint.level < 0 | confint.level > 1) {
stop("Argument 'confint.level' must be between 0 and 1.")
}
if (confint.alpha < 0 | confint.alpha > 1) {
stop("Argument 'confint.alpha' must be between 0 and 1.")
}
y <- NULL
xval <- x$model$data[, "y"]
if (x$method %in% c("E-UW", "E-FIV")) {
ylab <- "Standard error"
yval <- (x$model$data[, "x"])
ylim <- rev(c(0, max(yval, na.rm = T)))
xlim <- c(min(c(0, xval)), max(xval))
} else if (x$method %in% c("M-FIV")) {
ylab <- "Sample size"
yval <- (x$model$data[, "x"])
ylim <- (c(0, max(yval, na.rm = T)))
xlim <- c(min(c(0, xval)), max(xval))
} else if (x$method == "P-FPV") {
ylab <- "Sample size"
yval <- (x$model$data[, "x"])
ylim <- rev(c(0, max(yval, na.rm = T)))
xlim <- c(min(c(0, xval)), max(xval))
step <- ((max(yval) - min(yval))/5)
yax <- c(plyr::round_any(1/min(yval), 10^(sapply(round(1/min(yval)), nchar) - 1)),
plyr::round_any(1/seq(step, 4 * step, by = step), 10), plyr::round_any(1/max(yval), 10))
} else if (x$method == "D-FIV") {
ylab <- "Total events"
yval <- (x$model$data[, "x"])
ylim <- rev(c(0, max(yval, na.rm = T)))
xlim <- c(min(c(0, xval)), max(xval))
step <- ((max(yval) - min(yval))/4)
yax <- c(plyr::round_any(1/min(yval), 10^(sapply(round(1/min(yval)), nchar) - 1)),
plyr::round_any(1/seq(step, 4 * step, by = step), 10), plyr::round_any(1/max(yval), 10))
} else if (x$method == "D-FAV") {
ylab <- "Total events"
yval <- (x$model$data[, "x"])
ylim <- rev(c(0, max(yval, na.rm = T)))
xlim <- c(min(c(0, xval)), max(xval))
step <- ((max(yval) - min(yval))/4)
yax <- c(plyr::round_any(1/min(yval),10^(sapply(round(1/min(yval)), nchar) - 1)),
plyr::round_any(1/seq(step, 4 * step, by = step), 10), plyr::round_any(1/max(yval), 10))
} else {
stop("Plot not supported!")
}
newdata <- sort(c(-max(x$model$data[, "x"]), x$model$data[,"x"], 2 * max(x$model$data[, "x"])))
newdata <- as.data.frame(cbind(seq(min(newdata), max(newdata), length.out = 500), NA))
colnames(newdata) <- c("x", "y")
predy <- predict(x$model, newdata = newdata, se.fit = T)
predy.mean <- predy$fit
predy.lowerInt <- as.vector(predy$fit + qt(confint.level/2, df = x$df) * predy$se.fit)
predy.upperInt <- as.vector(predy$fit + qt((1 - confint.level/2), df = x$df) * predy$se.fit)
# restricting plotting range to the selection
predy.upperInt[predy.upperInt < min(pretty(range(xlim)))] <- min(pretty(range(xlim)))
predy.upperInt[predy.upperInt > max(pretty(range(xlim)))] <- max(pretty(range(xlim)))
predy.lowerInt[predy.lowerInt < min(pretty(range(xlim)))] <- min(pretty(range(xlim)))
predy.lowerInt[predy.lowerInt > max(pretty(range(xlim)))] <- max(pretty(range(xlim)))
newdata[, "x"][newdata[, "x"] < min(pretty(range(ylim)))] <- min(pretty(range(ylim)))
newdata[, "x"][newdata[, "x"] > max(pretty(range(ylim)))] <- max(pretty(range(ylim)))
p <- ggplot2::ggplot(data = data.frame(x = xval, y = yval))
if (confint) {
p <- p + ggplot2::geom_polygon(
mapping = ggplot2::aes(
x = x,
y = y
),
data = data.frame(
x = c(
predy.upperInt,
rev(predy.lowerInt)),
y = c(
newdata[, "x"],
rev(newdata[, "x"]))
),
fill = confint.col,
alpha = confint.alpha
)
}
p <- p +
ggplot2::geom_point(
mapping = ggplot2::aes(x = x, y = y),
shape = 19
) +
ggplot2::geom_line(
mapping = ggplot2::aes(x = x, y = y),
data = data.frame(
x = predy.mean[newdata[, "x"] > min(pretty(range(ylim))) & newdata[, "x"] < max(pretty(range(ylim)))], # another plotting range restriction
y = newdata[, "x"][newdata[, "x"] > min(pretty(range(ylim))) & newdata[, "x"] < max(pretty(range(ylim)))] # another plotting range restriction
),
linetype = 2)
if (missing(ref)) {
p <- p + ggplot2::geom_vline(xintercept = x$fema$b)
} else {
p <- p + ggplot2::geom_vline(xintercept = ref)
}
p <- p + ggplot2::scale_x_continuous(
name = xlab,
limits = range(pretty(range(xlim))),
breaks = pretty(range(xlim)))
if (x$method %in% c("P-FPV", "D-FAV", "D-FIV")){
p <- p + ggplot2::scale_y_reverse(name = ylab, breaks = 1/yax, labels = yax, limits = rev(range(pretty(ylim))))
} else if (x$method %in% c("E-UW", "E-FIV")){
p <- p + ggplot2::scale_y_reverse(name = ylab, limits = rev(range(pretty(ylim))), breaks = pretty(range(ylim)))
} else {
p <- p + ggplot2::scale_y_continuous(name = ylab, limits = range(pretty(ylim)), breaks = pretty(range(ylim)))
}
return(p)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.