# R/pooling_metrics.R In baggr: Bayesian Aggregate Treatment Effects

#### Documented in heterogeneitypoolingweights.baggr

#' Pooling metrics and related statistics for baggr
#'
#' Compute statistics relating to
#' pooling in a given [baggr] meta-analysis model returns statistics, for
#' either the entire model or individual groups, such as
#' pooling statistic by Gelman & Pardoe (2006), _I-squared_, _H-squared_, or study weights;
#' heterogeneity is a shorthand for pooling(type = "total")
#' weights is shorthand for pooling(metric = "weights")
#'
#' @param bg a [baggr] model
#' @param metric "pooling" for Gelman & Pardoe statistic _P_,
#'               "isq" for I-squared statistic (_1-P_, Higgins & Thompson, 2002)
#'               "hsq" for H squared statistic (_1/P_, ibid.);
#'               "weights" for study weights;
#'               also see _Details_
#' @param type In pooling calculation is done for each of the "groups"
#'            (default) or for "total" hypereffect(s).
#' @param summary logical; if FALSE a whole vector of pooling values is returned,
#'                otherwise only the means and intervals
#'
#' @details
#'
#' Pooling statistic (Gelman & Pardoe, 2006) describes the extent to which
#' group-level estimates of treatment
#' effect are "pooled" toward average treatment effect in the meta-analysis model.
#' If pooling = "none" or "full" (which you specify when calling [baggr]),
#' then the values are always 0 or 1, respectively.
#' If pooling = "partial", the value is somewhere between 0 and 1.
#' We can distinguish between pooling of individual groups and overall pooling in
#' the model.
#'
#' In many contexts, i.e. medical statistics, it is typical to report _1-P_, called \eqn{I^2}
#' (see Higgins and Thompson, 2002; sometimes another statistic, \eqn{H^2 = 1 / P},
#' is used).
#' Higher values of _I-squared_ indicate higher heterogeneity;
#' Von Hippel (2015) provides useful details for _I-squared_ calculations (and some
#' issues related to it, especially in frequentist models).
#' See Gelman & Pardoe (2006) Section 1.1 for a short explanation of how \eqn{R^2}
#' statistic relates to the pooling metric.
#'
#' @section Group pooling:
#'
#' This is the calculation done by pooling() if type = "groups" (default).
#' In a partial pooling model (see [baggr] and above), group _k_ (e.g. study) has
#' standard error of treatment effect estimate, \eqn{se_k}.
#' The treatment effect (across _k_ groups) is variable across groups, with
#' hyper-SD parameter \eqn{\sigma_(\tau)}.
#'
#' The quantity of interest is ratio of variation in treatment effects to the
#' total variation.
#' By convention, we subtract it from 1, to obtain a _pooling metric_ _P_.
#'
#' \deqn{p = 1 - (\sigma_(\tau)^2 / (\sigma_(\tau)^2 + se_k^2))}
#'
#' * If \eqn{p < 0.5}, the variation across studies is higher than variation within studies.
#' * Values close to 1 indicate nearly full pooling. Variation across studies dominates.
#' * Values close to 0 indicate no pooling. Variation within studies dominates.
#'
#' Note that, since \eqn{\sigma_{\tau}^2} is a Bayesian parameter (rather than a
#' single fixed value),
#' _P_ is also a parameter. It is typical for _P_ to have very high dispersion,
#' as in many cases we
#' cannot precisely estimate \eqn{\sigma_{\tau}}. To obtain samples from the distribution
#' of _P_ (rather than summarised values), set summary=FALSE.
#'
#'
#' @section Study weights:
#'
#' Contributions of each group (e.g. each study) to the mean meta-analysis estimate
#' can be calculated by calculating for each study *w_k* the inverse of sum of group-specific
#' SE squared and between-study variation.
#' To obtain weights, this vector (across all studies) has to be normalised to 1, i.e.
#' *w_k/sum(w_k)* for each _k_.
#'
#' SE is typically treated as a fixed quantity
#' (and usually reported on the reported point estimate),
#' but between-study variance is a model parameter,
#' hence the weights themselves are also random variables.
#'
#'
#'
#'
#' @section Overall pooling in the model:
#'
#' Typically researchers want to report a single measure from the model,
#' relating to heterogeneity across groups.
#' This is calculated by either pooling(mymodel, type = "total") or simply
#' heterogeneity(mymodel)
#'
#' Formulae for the calculations below are provided in main package vignette and
#' almost analogous to the group calculation above, but using mean variance across
#' all studies. In other words, pooling _P_ is simply ratio of the expected within-study
#' variance term to total variance.
#'
#' To obtain such single estimate we need to substitute average variability of group-specific
#' treatment effects and then calculate the same way we would calculate \eqn{p}.
#' By default we use the mean across _k_ \eqn{se_k^2} values. Typically, implementations of
#' \eqn{I^2} in statistical packages use a different calculation for this quantity,
#' which may make _I_'s not comparable when different studies have different SE's.
#'
#' Same as for group-specific estimates, _P_ is a Bayesian parameter and its
#' dispersion can be high.
#'
#'
#' @section Value:
#'
#' Matrix with mean and intervals for chosen pooling metric,
#' each row corresponding to one meta-analysis group.
#'
#' @references
#' Gelman, Andrew, and Iain Pardoe.
#' "Bayesian Measures of Explained Variance and Pooling in Multilevel (Hierarchical) Models."
#' _Technometrics 48, no. 2 (May 2006): 241-51_.
#'
#' Higgins, Julian P. T., and Simon G. Thompson.
#' "Quantifying Heterogeneity in a Meta-Analysis."
#' _Statistics in Medicine, vol. 21, no. 11, June 2002, pp. 1539-58_.
#'
#' Hippel, Paul T von. "The Heterogeneity Statistic I2 Can Be Biased in Small Meta-Analyses."
#' _BMC Medical Research Methodology 15 (April 14, 2015)._
#'
#'
#' @export
#'

pooling <- function(bg,
metric = c("pooling", "isq", "hsq", "weights"),
type = c("groups", "total"),
summary = TRUE) {
type <- match.arg(type, c("groups", "total"))
metric <- match.arg(metric, c("pooling", "isq", "hsq", "weights"))

# we have to rig it for no pooling cases
# because sigma_tau parameter might be meaningless then
if(bg$pooling == "none") return(array(0, c(3, bg$n_groups, bg$n_parameters))) if(bg$pooling == "full")
return(array(1, c(3, bg$n_groups, bg$n_parameters)))

# we'll replace by switch() in the future
# if(bg$model %in% c("rubin", "mutau", "logit", "rubin_full")) { if(bg$n_parameters == 1) {
# 1-dimensional vector of S values (S=N samples)
sigma_tau <- treatment_effect(bg)$sigma_tau # Grab the appropriate SE (k values) sigma_k <- switch(bg$model,
"mutau"       = bg$data$se.tau,
"rubin"       = bg$data$se,
# May need to apply rare event correction for
# the pooling models to work correctly
"logit"       = apply_cont_corr(bg$summary_data, 0.25, "single", add_or = TRUE, pooling = TRUE)$se,
# These are SDs after pooling, don't use this (here for tests)
# "rubin_full"  = group_effects(bg, summary = TRUE)[, "sd", 1],
"rubin_full"  = bg$summary_data$se.tau,
"mutau_full"  = bg$summary_data$se.tau
)

if(type == "groups")
ret <- sapply(sigma_k, function(se) se^2 / (se^2 + sigma_tau^2))
if(type == "total")
ret <- replicate(1, mean(sigma_k^2) / (mean(sigma_k^2) + sigma_tau^2))
if(metric == "weights" && type == "groups"){
precisions <- sapply(sigma_k, function(se) 1 / (se^2 + sigma_tau^2))
ret <- t(apply(precisions, 1, function(x) x/sum(x)))
}
if(metric == "weights" && type == "total")
stop("Weights can be calculated only for type = 'groups'")

ret <- replicate(1, ret) #third dim is always N parameters, by convention,
#so we set it to 1 here

} else if(bg$model == "quantiles") { # compared to individual-level, here we are dealing with 1 more dimension # which is number of quantiles; so if sigma_tau above is sigma of trt effect # here it is N effects on N quantiles etc. # Nq-dimensional sigma_tau sigma_tau <- treatment_effect(bg)$sigma_tau
# for(i in 1:dim(sigma_tau))
# sigma_tau[,i,1] <- sigma_tau[,i,i]
# sigma_tau <- sigma_tau[,,1] #rows are samples, columns are quantiles

# now SE's of study treatment effects: this is our input data!
# Sigma_y_k_1 is now vcov of trt effect (used to be control + effect, mind)
sigma_k <- t(apply(bg$inputs$Sigma_y_k_1, 1, diag)) #rows are studies, cols are quantiles

# output is an array with dim = samples, studies, parameters
# so let's convert to these dimensions:
# for sigma_k the values are fixed over samples
sigma_k <- array(rep(sigma_k, each = nrow(sigma_tau)),
dim = c(nrow(sigma_tau), nrow(sigma_k), ncol(sigma_k)))
# for sigma_tau it doesn't change with studies, so:
sigma_tau <- replicate(dim(sigma_k), sigma_tau)
# but this has (parameters, studies) rather than (studies, parameters) so
sigma_tau <- aperm(sigma_tau, c(1, 3, 2))

# now we can just do the operation on arrays and preserve dimensions:
ret <- sigma_k^2 / (sigma_k^2 + sigma_tau^2)

if(type == "total"){
warning("Total pooling not implemented for quantiles model")
return(array(NA, c(3, bg$n_groups, bg$n_parameters)))
}
} else if(bg$model == "sslab") { if(type == "total"){ warning("Total pooling not implemented for spike & slab model") return(array(NA, c(3, bg$n_groups, bg$n_parameters))) } # Must calculate pooling for each effect vector? te <- treatment_effect(bg)$sigma_tau
warning("In this version of baggr calculations of pooling for spike & slab may be wrong")
warning("Please contact package authors if youy are using this feature")
ge <- group_effects(bg, summary = T)[,"sd",]
sigma_tau <- aperm(replicate(5, te), c(1,3,2))
sigma_k   <- aperm(replicate(dim(sigma_tau), ge), c(3,1,2))
ret <- sigma_k^2 / (sigma_k^2 + sigma_tau^2)

} else {
stop("Cannot calculate pooling metrics.")
}

if(metric == "isq") ret <- 1-ret
if(metric == "hsq") ret <- 1/ret
if(metric == "h")   ret <- sqrt(1/ret)

if(summary)
ret <- apply(ret, c(2,3), mint)

# if(bg\$n_parameters == 1)
#   ret <- ret[1,,]

return(ret)

}

#' @rdname pooling
#' @export
heterogeneity <- function(
bg,
metric = c("pooling", "isq", "hsq", "weights"),
summary = TRUE)
pooling(bg, metric, type = "total", summary = summary)

#' @rdname pooling
#' @param object [baggr] model for which to calculate group (study) weights
#' @param ... Unused, please ignore.
#' @export
#' @importFrom stats weights
weights.baggr <- function(object, ...)
pooling(object, metric = "weights", type = "groups", ...)


## Try the baggr package in your browser

Any scripts or data that you put into this service are public.

baggr documentation built on March 18, 2022, 7:02 p.m.