R/mapower_ul.R

In dsquintana/metameta: A Suite of Tools for calculating and visualising study-level statistical power for meta-analyses

#' Calculate statistical power using effect size and confidence interval data
#'
#' Calculate statistical power using effect size and confidence interval data that
#' is typically reported in meta-analyses. Power is calculated assuming
#' the true effect size and a range of effect sizes. The script to calculate power
#' was adapted from the "viz_sunset" function of the metaviz package.
#'
#' \strong{The true effect size}
#'
#' For the purposes of power analysis a "true" effect size needs to be specified. However,
#' the true effect size is unknown when testing hypotheses and reported summary
#' effect sizes are typically inflated. Here, the observed summary effect size estimate 
#' reported in the meta-analysis results is used as one possible true effect size, however, 
#' other effect sizes can be specified here instead. Additionally, statistical power for a range 
#' of true effect sizes are returned. The default is a "medium" range from 0.1 to 1, 
#' in increments of 0.1. However, a "small" and "large" range is also available 
#' (see the "size" argument below for details).
#'
#' @param dat A dataset that contains one column with observed effect sizes or outcomes
#'     labelled "yi", a column labelled "lower" with the lower confidence
#'     interval bound, and column labelled "upper" with the upper confidence
#'     interval bound. This function assumes a 95\% confidence interval.
#' @param size The effect size range. The default is a medium range (0.1 to 1, in increments of 0.1).
#'     In other words, if you do not use this argument, the medium range will be used.
#'     Other options include a "small" range (0.05 to 0.5 in increments of 0.05) and a "large"
#'     range (0.25 to 2.5 in increments of 0.25).
#' @param observed_es The observed summary effect size estimate for the meta-analysis, which is one output
#'     for the set of possible "true" effect sizes. See 'Details'
#' @param name A label with a name for the meta-analysis, which is required for using the "firepower" function.
#'     See 'Examples'
#' @return This function returns the following:
#'       \item{dat}{A dataset with results from power analyses for a range of effect sizes, including the
#'       specified observed effect size, in a column labelled "es_observed". The additional added columns include
#'       results for power analysis assuming a range of true effect sizes, with the default beginning at 0.1 ("power_es01"),
#'       then 0.2 ("power_es02"), then continuing in increments of 0.1 up to 1 ("power_es1"). A smaller
#'       or larger range of effect sizes can be specified (see the 'size' argument for details)}
#' @examples
#' 
#' ### An analysis using the default effect size range
#' 
#' keech_power <- mapower_ul(
#'     dat = dat_keech,
#'     observed_es = 0.08,
#'     name = "Keech et al 2017")
#' keech_power
#'
#' ### An analysis using the large effect size range
#'
#' keech_power_l <- mapower_ul(
#'                  dat = dat_keech,
#'                  size = "large",
#'                  observed_es = 0.08,
#'                  name = "Keech et al 2017")
#' keech_power_l
#'
#' ### An analysis using the small effect size range
#'
#' keech_power_s <- mapower_ul(
#'                  dat = dat_keech,
#'                  size = "small",
#'                  observed_es = 0.08,
#'                  name = "Keech et al 2017")
#' keech_power_s
#' 


mapower_ul <-
function(dat,
         size = "medium",
         observed_es,
         name)
#'@import dplyr
{
  dat[["sei"]] <- ((dat[["upper"]]) -
                     (dat[["lower"]])) / (2 * 1.96) # Calculate standard error if not already available
  if(size == "small"){
    dat[["power_es_observed"]] <-
      (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                        abs(observed_es), dat[["sei"]])) +
      stats::pnorm(stats::qnorm(0.05 / 2) *
                     dat[["sei"]], abs(observed_es),
                   dat[["sei"]]) # Calculate power for each study for observed effect
    dat[["power_es005"]] <-
      (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                        abs(0.05), dat[["sei"]])) +
      stats::pnorm(stats::qnorm(0.05 / 2) *
                     dat[["sei"]], abs(0.05),
                   dat[["sei"]]) # Calculate power for each study for an effect of 0.05
    dat[["power_es01"]] <-
      (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                        abs(0.1), dat[["sei"]])) +
      stats::pnorm(stats::qnorm(0.05 / 2) *
                     dat[["sei"]], abs(0.1),
                   dat[["sei"]]) # Calculate power for each study for an effect of 0.1
    dat[["power_es015"]] <-
      (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                        abs(0.15), dat[["sei"]])) +
      stats::pnorm(stats::qnorm(0.05 / 2) *
                     dat[["sei"]], abs(0.15),
                   dat[["sei"]]) # Calculate power for each study for an effect of 0.15
    dat[["power_es02"]] <-
      (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                        abs(0.2), dat[["sei"]])) +
      stats::pnorm(stats::qnorm(0.05 / 2) *
                     dat[["sei"]], abs(0.2),
                   dat[["sei"]]) # Calculate power for each study for an effect of 0.2
    dat[["power_es025"]] <-
      (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                        abs(0.25), dat[["sei"]])) +
      stats::pnorm(stats::qnorm(0.05 / 2) *
                     dat[["sei"]], abs(0.25),
                   dat[["sei"]]) # Calculate power for each study for an effect of 0.25
    dat[["power_es03"]] <-
      (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                        abs(0.3), dat[["sei"]])) +
      stats::pnorm(stats::qnorm(0.05 / 2) *
                     dat[["sei"]], abs(0.3),
                   dat[["sei"]]) # Calculate power for each study for an effect of 0.3
    dat[["power_es035"]] <-
      (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                        abs(0.35), dat[["sei"]])) +
      stats::pnorm(stats::qnorm(0.05 / 2) *
                     dat[["sei"]], abs(0.35),
                   dat[["sei"]]) # Calculate power for each study for an effect of 0.35
    dat[["power_es04"]] <-
      (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                        abs(0.4), dat[["sei"]])) +
      stats::pnorm(stats::qnorm(0.05 / 2) *
                     dat[["sei"]], abs(0.4),
                   dat[["sei"]]) # Calculate power for each study for an effect of 0.4
    dat[["power_es045"]] <-
      (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                        abs(0.45), dat[["sei"]])) +
      stats::pnorm(stats::qnorm(0.05 / 2) *
                     dat[["sei"]], abs(0.45),
                   dat[["sei"]]) # Calculate power for each study for an effect of 0.45
    dat[["power_es05"]] <-
      (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                        abs(0.5), dat[["sei"]])) +
      stats::pnorm(stats::qnorm(0.05 / 2) *
                     dat[["sei"]], abs(0.5),
                   dat[["sei"]]) # Calculate power for each study for an effect of 0.5
    dat <- as.data.frame(dat)
    power_median_dat <- data.frame(power_es_observed=numeric(1),
                                   es005=numeric(1),
                                   es01=numeric(1),
                                   es015=numeric(1),
                                   es02=numeric(1),
                                   es025=numeric(1),
                                   es03=numeric(1),
                                   es035=numeric(1),
                                   es04=numeric(1),
                                   es045=numeric(1),
                                   es05=numeric(1)
    )
    power_median_dat[["es_observed"]] <- median(dat[["power_es_observed"]])
    power_median_dat[["es005"]] <- median(dat[["power_es005"]])
    power_median_dat[["es01"]] <- median(dat[["power_es01"]])
    power_median_dat[["es015"]] <- median(dat[["power_es015"]])
    power_median_dat[["es02"]] <- median(dat[["power_es02"]])
    power_median_dat[["es025"]] <- median(dat[["power_es025"]])
    power_median_dat[["es03"]] <- median(dat[["power_es03"]])
    power_median_dat[["es035"]] <- median(dat[["power_es035"]])
    power_median_dat[["es04"]] <- median(dat[["power_es04"]])
    power_median_dat[["es045"]] <- median(dat[["power_es045"]])
    power_median_dat[["es05"]] <- median(dat[["power_es05"]])
    power_median_dat <- as.data.frame(power_median_dat)
    power_median_dat <- dplyr::mutate(power_median_dat,
                                      meta_analysis_name = name)
  } else {
    if(size == "medium"){
      dat[["power_es_observed"]] <-
        (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                          abs(observed_es), dat[["sei"]])) +
        stats::pnorm(stats::qnorm(0.05 / 2) *
                       dat[["sei"]], abs(observed_es),
                     dat[["sei"]]) # Calculate power for each study for observed effect
      dat[["power_es01"]] <-
        (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                          abs(0.1), dat[["sei"]])) +
        stats::pnorm(stats::qnorm(0.05 / 2) *
                       dat[["sei"]], abs(0.1),
                     dat[["sei"]]) # Calculate power for each study for an effect of 0.1
      dat[["power_es02"]] <-
        (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                          abs(0.2), dat[["sei"]])) +
        stats::pnorm(stats::qnorm(0.05 / 2) *
                       dat[["sei"]], abs(0.2),
                     dat[["sei"]]) # Calculate power for each study for an effect of 0.2
      dat[["power_es03"]] <-
        (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                          abs(0.30), dat[["sei"]])) +
        stats::pnorm(stats::qnorm(0.05 / 2) *
                       dat[["sei"]], abs(0.30),
                     dat[["sei"]]) # Calculate power for each study for an effect of 0.3
      dat[["power_es04"]] <-
        (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                          abs(0.40), dat[["sei"]])) +
        stats::pnorm(stats::qnorm(0.05 / 2) *
                       dat[["sei"]], abs(0.40),
                     dat[["sei"]]) # Calculate power for each study for an effect of 0.4
      dat[["power_es05"]] <-
        (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                          abs(0.50), dat[["sei"]])) +
        stats::pnorm(stats::qnorm(0.05 / 2) *
                       dat[["sei"]], abs(0.50),
                     dat[["sei"]]) # Calculate power for each study for an effect of 0.5
      dat[["power_es06"]] <-
        (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                          abs(0.60), dat[["sei"]])) +
        stats::pnorm(stats::qnorm(0.05 / 2) *
                       dat[["sei"]], abs(0.60),
                     dat[["sei"]]) # Calculate power for each study for an effect of 0.6
      dat[["power_es07"]] <-
        (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                          abs(0.70), dat[["sei"]])) +
        stats::pnorm(stats::qnorm(0.05 / 2) *
                       dat[["sei"]], abs(0.70),
                     dat[["sei"]]) # Calculate power for each study for an effect of 0.7
      dat[["power_es08"]] <-
        (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                          abs(0.80), dat[["sei"]])) +
        stats::pnorm(stats::qnorm(0.05 / 2) *
                       dat[["sei"]], abs(0.80),
                     dat[["sei"]]) # Calculate power for each study for an effect of 0.8
      dat[["power_es09"]] <-
        (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                          abs(0.90), dat[["sei"]])) +
        stats::pnorm(stats::qnorm(0.05 / 2) *
                       dat[["sei"]], abs(0.90),
                     dat[["sei"]]) # Calculate power for each study for an effect of 0.9
      dat[["power_es1"]] <-
        (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                          abs(1), dat[["sei"]])) +
        stats::pnorm(stats::qnorm(0.05 / 2) *
                       dat[["sei"]], abs(1),
                     dat[["sei"]]) # Calculate power for each study for an effect of 1
      dat <- as.data.frame(dat)
      power_median_dat <- data.frame(power_es_observed=numeric(1),
                                     es01=numeric(1),
                                     es02=numeric(1),
                                     es03=numeric(1),
                                     es04=numeric(1),
                                     es05=numeric(1),
                                     es06=numeric(1),
                                     es07=numeric(1),
                                     es08=numeric(1),
                                     es09=numeric(1),
                                     es1=numeric(1)
                                     
      )
      power_median_dat[["es_observed"]] <- median(dat[["power_es_observed"]])
      power_median_dat[["es01"]] <- median(dat[["power_es01"]])
      power_median_dat[["es02"]] <- median(dat[["power_es02"]])
      power_median_dat[["es03"]] <- median(dat[["power_es03"]])
      power_median_dat[["es04"]] <- median(dat[["power_es04"]])
      power_median_dat[["es05"]] <- median(dat[["power_es05"]])
      power_median_dat[["es06"]] <- median(dat[["power_es06"]])
      power_median_dat[["es07"]] <- median(dat[["power_es07"]])
      power_median_dat[["es08"]] <- median(dat[["power_es08"]])
      power_median_dat[["es09"]] <- median(dat[["power_es09"]])
      power_median_dat[["es1"]] <- median(dat[["power_es1"]])
      power_median_dat <- as.data.frame(power_median_dat)
      power_median_dat <- dplyr::mutate(power_median_dat,
                                        meta_analysis_name = name)
    } else {
      if(size == "large"){
        dat[["power_es_observed"]] <-
          (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                            abs(observed_es), dat[["sei"]])) +
          stats::pnorm(stats::qnorm(0.05 / 2) *
                         dat[["sei"]], abs(observed_es),
                       dat[["sei"]]) # Calculate power for each study for observed effect
        dat[["power_es025"]] <-
          (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                            abs(0.25), dat[["sei"]])) +
          stats::pnorm(stats::qnorm(0.05 / 2) *
                         dat[["sei"]], abs(0.25),
                       dat[["sei"]]) # Calculate power for each study for an effect of 0.25
        dat[["power_es05"]] <-
          (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                            abs(0.5), dat[["sei"]])) +
          stats::pnorm(stats::qnorm(0.05 / 2) *
                         dat[["sei"]], abs(0.5),
                       dat[["sei"]]) # Calculate power for each study for an effect of 0.5
        dat[["power_es075"]] <-
          (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                            abs(0.75), dat[["sei"]])) +
          stats::pnorm(stats::qnorm(0.05 / 2) *
                         dat[["sei"]], abs(0.75),
                       dat[["sei"]]) # Calculate power for each study for an effect of 0.75
        dat[["power_es1"]] <-
          (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                            abs(1), dat[["sei"]])) +
          stats::pnorm(stats::qnorm(0.05 / 2) *
                         dat[["sei"]], abs(1),
                       dat[["sei"]]) # Calculate power for each study for an effect of 1
        dat[["power_es125"]] <-
          (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                            abs(1.25), dat[["sei"]])) +
          stats::pnorm(stats::qnorm(0.05 / 2) *
                         dat[["sei"]], abs(1.25),
                       dat[["sei"]]) # Calculate power for each study for an effect of 1.25
        dat[["power_es15"]] <-
          (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                            abs(1.5), dat[["sei"]])) +
          stats::pnorm(stats::qnorm(0.05 / 2) *
                         dat[["sei"]], abs(1.5),
                       dat[["sei"]]) # Calculate power for each study for an effect of 1.5
        dat[["power_es175"]] <-
          (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                            abs(1.75), dat[["sei"]])) +
          stats::pnorm(stats::qnorm(0.05 / 2) *
                         dat[["sei"]], abs(1.75),
                       dat[["sei"]]) # Calculate power for each study for an effect of 1.75
        dat[["power_es2"]] <-
          (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                            abs(2), dat[["sei"]])) +
          stats::pnorm(stats::qnorm(0.05 / 2) *
                         dat[["sei"]], abs(2),
                       dat[["sei"]]) # Calculate power for each study for an effect of 2
        dat[["power_es225"]] <-
          (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                            abs(2.25), dat[["sei"]])) +
          stats::pnorm(stats::qnorm(0.05 / 2) *
                         dat[["sei"]], abs(2.25),
                       dat[["sei"]]) # Calculate power for each study for an effect of 2.25
        dat[["power_es25"]] <-
          (1 - stats::pnorm(stats::qnorm(1 - 0.05 / 2) * dat[["sei"]],
                            abs(2.5), dat[["sei"]])) +
          stats::pnorm(stats::qnorm(0.05 / 2) *
                         dat[["sei"]], abs(2.5),
                       dat[["sei"]]) # Calculate power for each study for an effect of 2.5
        dat <- as.data.frame(dat)
        power_median_dat <- data.frame(power_es_observed=numeric(1),
                                       es025=numeric(1),
                                       es05=numeric(1),
                                       es075=numeric(1),
                                       es1=numeric(1),
                                       es125=numeric(1),
                                       es15=numeric(1),
                                       es175=numeric(1),
                                       es2=numeric(1),
                                       es225=numeric(1),
                                       es25=numeric(1)
                                       
        )
        power_median_dat[["es_observed"]] <- median(dat[["power_es_observed"]])
        power_median_dat[["es025"]] <- median(dat[["power_es025"]])
        power_median_dat[["es05"]] <- median(dat[["power_es05"]])
        power_median_dat[["es075"]] <- median(dat[["power_es075"]])
        power_median_dat[["es1"]] <- median(dat[["power_es1"]])
        power_median_dat[["es125"]] <- median(dat[["power_es125"]])
        power_median_dat[["es15"]] <- median(dat[["power_es15"]])
        power_median_dat[["es175"]] <- median(dat[["power_es175"]])
        power_median_dat[["es2"]] <- median(dat[["power_es2"]])
        power_median_dat[["es225"]] <- median(dat[["power_es225"]])
        power_median_dat[["es25"]] <- median(dat[["power_es25"]])
        power_median_dat <- as.data.frame(power_median_dat)
        power_median_dat <- dplyr::mutate(power_median_dat,
                                          meta_analysis_name = name)
      }
      
    }
  }
  value <- list(
    dat = dat,
    power_median_dat = power_median_dat
  ) # Create a list of output objects
  attr(value, "class") <- "mapower_se_t2"
  value
}