## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)

## ----message=FALSE, warning=FALSE---------------------------------------------
library(confMeta)
library(ggplot2)

set.seed(42)

n <- 5
conf_level <- 0.95
estimates <- rnorm(n, mean = 0.5, sd = 0.2)
SEs <- rgamma(n, shape = 5, rate = 20)
study_names <- c("Study A", "Study B", "Study C", "Study D", "Study E")

## ----eval = FALSE-------------------------------------------------------------
# ?confMeta
# ?p_value_functions

## -----------------------------------------------------------------------------
cm_edgington <- confMeta(
  estimates = estimates,
  SEs = SEs,
  study_names = study_names,
  conf_level = conf_level,
  fun = p_edgington,
  fun_name = "Edgington"
)

## -----------------------------------------------------------------------------

# Use the specific print argument
print(cm_edgington)

# See what elements it has
names(cm_edgington)


# Check out the combined confidence interval(s)
cm_edgington$joint_cis

## -----------------------------------------------------------------------------
# Test the null hypothesis that mu = 0 using Edgington's method
p_val_0 <- p_edgington(
  estimates = estimates, 
  SEs = SEs, 
  mu = 0, 
  input_p = "two.sided", 
  output_p = "two.sided"
)

p_val_0

## ----eval = FALSE-------------------------------------------------------------
# ?autoplot.confMeta

## -----------------------------------------------------------------------------
# show the p-value function
autoplot(cm_edgington, type = "p")

## -----------------------------------------------------------------------------
# show the forest plot
autoplot(cm_edgington, type = "forest")

## ----fig.height = 6-----------------------------------------------------------
# show both
autoplot(cm_edgington)

## -----------------------------------------------------------------------------
cm_fisher <- confMeta(
  estimates = estimates, SEs = SEs, study_names = study_names,
  fun = p_fisher, fun_name = "Fisher", input_p = "two.sided"
)

cm_hmean <- confMeta(
  estimates = estimates, SEs = SEs, study_names = study_names,
  fun = p_hmean, fun_name = "Harmonic Mean", distr = "chisq"
)

## ----fig.height=8-------------------------------------------------------------
autoplot(cm_edgington, cm_fisher, cm_hmean)

## -----------------------------------------------------------------------------
weights <- 1 / SEs

cm_weighted <- confMeta(
  estimates = estimates, SEs = SEs,
  fun = p_edgington_w, fun_name = "Edgington (Weighted)",
  w = weights, input_p = "two.sided"
)

## -----------------------------------------------------------------------------
# Estimate tau^2
tau2_est <- estimate_tau2(estimates, SEs, method.tau = "PM")

# Create a random-effects confMeta object using additive heterogeneity
cm_additive <- confMeta(
  estimates = estimates, 
  SEs = SEs,
  study_names = study_names,
  fun = p_edgington, 
  fun_name = "Edgington (Additive RE)",
  heterogeneity = "additive", 
  tau2 = tau2_est,
  input_p = "two.sided"
)

## -----------------------------------------------------------------------------
# Estimate phi
phi_est <- estimate_phi(estimates, SEs)

# Create a confMeta object using multiplicative heterogeneity
cm_multiplicative <- confMeta(
  estimates = estimates, 
  SEs = SEs,
  study_names = study_names,
  fun = p_edgington, 
  fun_name = "Edgington (Multiplicative)",
  heterogeneity = "multiplicative", 
  phi = phi_est,
  input_p = "two.sided"
)

## ----fig.height=6-------------------------------------------------------------
autoplot(cm_edgington, cm_additive, cm_multiplicative)

## ----message=FALSE, warning=FALSE---------------------------------------------

if (requireNamespace("bayesmeta", quietly = TRUE)) {
  
  # Run a standard bayesmeta model
  bm <- bayesmeta::bayesmeta(y = estimates, sigma = SEs)
  
  # Add the bayesmeta posterior diamond to the forest plot
  autoplot(cm_edgington, type = "forest", bayesmeta = bm)
}

## -----------------------------------------------------------------------------
# Simulated 2x2 table data for 3 clinical trials
clinical_data <- data.frame(
  ai = c(12, 5, 22),   # Events in experimental group
  bi = c(88, 45, 78),  # Non-events in experimental group
  ci = c(18, 10, 30),  # Events in control group
  di = c(82, 40, 70),  # Non-events in control group
  n1i = c(100, 50, 100),
  n2i = c(100, 50, 100)
)

# Use metafor::escalc to obtain log Odds Ratios (yi) and variances (vi)
if (requireNamespace("metafor", quietly = TRUE)) {
  
  es <- metafor::escalc(
    measure = "OR", 
    ai = ai, bi = bi, ci = ci, di = di, 
    data = clinical_data
  )
  
  # Create a confMeta object using MH pooling for the Fixed Effect reference
  cm_clinical <- confMeta(
    estimates = es$yi, 
    SEs = sqrt(es$vi),         # Remember to take the square root of the variance!
    study_names = c("Trial 1", "Trial 2", "Trial 3"),
    fun = p_fisher, 
    fun_name = "Fisher (Clinical)",
    MH = TRUE,                 # Turn on Mantel-Haenszel pooling
    table_2x2 = clinical_data, # Provide the raw table
    measure = "OR"             # Specify the effect measure
  )

  # Notice how the "Fixed effect" comparison now uses the MH method
  cm_clinical$comparison_cis
}