Power function for testing τ = 0 for random effects meta-analyses

power_tau.R

# Based on methods described by Hedges and Pigott, https://doi.org/10.1037/1082-989X.9.4.426, p. 438
# Written 2019-10-25 by Brenton M. Wiernik
# Licensed GPL v3.0

power_curve <- function(object, tau, level = .95) {
  if (! inherits(object, "rma.uni")) stop("'object' must be an object of class 'rma.uni'")
  tau2 <- tau^2
  obs_tau2 <- object$tau2
  obs_tau <- sqrt(obs_tau2)
  obs_H2 <- object$H2
  s2 <- obs_tau2 / (obs_H2 - 1)
  df <- object$k - object$p
  crit <- qchisq(.95, df)
  wi <- 1 / object$vi
  a <- sum(wi) - sum(wi^2) / sum(wi)
  mu_q <- function(tau2) a * tau2 + df
  b <- function(tau2) df + 2 * a * tau2 + (sum(wi^2) - 2 * sum(wi^3) / sum(wi) + sum(wi^2)^2 / sum(wi)^2) * tau2^2
  var_q <- function(tau2) 2 * b(tau2)
  r <- function(tau2) var_q(tau2) / (2 * mu_q(tau2))
  s <- function(tau2) 2 * mu_q(tau2)^2 / var_q(tau2)
  H <- function(tau2) pchisq(crit / r(tau2), s(tau2))
  power <- 1 - H(tau2)
  out <- list(obs_tau = obs_tau, power = data.frame(tau = tau, power = power))
  return(out)
}

This function takes a metafor meta-analysis object, not a psychmeta one. tau isn't the estimated value of tau from your meta-analysis, it is a vector of population tau (SDrho, SDres) values you want to estimate power for.

If you want to get something analogous using a psychmeta object, estimate a psychmeta meta-analysis, then:

es <- get_escalc(psycmeta)[[1]][[1]] # or whichever subset is for the meta-analysis whose power you want to estimate
ma_obj <- rma.uni(yi, vi, weights = weight, method="DL", data = es)
sd_rho_values <- seq(0.05, .50, by = .05) # or whatever SDrho values you want to test
power_curve(ma_obj, sd_rho_values)

Excellent! thanks for your reply and that should fix things up properly for me on this end. Makes perfect sense it would reference a metafor object 😑 😆

Thanks for sharing the psychmeta option too!!

This is gonna be a handy tool to have. THANKS!

Quick note for anyone who finds this in the future.

If you are using the code for the psychmeta object you'll need to have library(metafor) loaded for the rma.uni() to function correctly.