Nothing

```
#' 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)[2])
# 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)[2], 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)[1], 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", ...)
```

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

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.