## ----setup, include=FALSE-----------------------------------------------------
library(ggplot2)
library(metaselection)
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)

## ----one-step, echo = FALSE, fig.cap = "One-step selection model with $\\lambda_1 = 0.4$"----
lambda1 <- 0.4
lambda2 <- 0.2
pvals <- seq(0,1,0.005)
PSM3 <- step_fun(cut_vals = 0.025, weights = lambda1)
PSM4 <- step_fun(cut_vals = c(0.025, 0.500), weights = c(lambda1, lambda2))

dat <- data.frame(p = pvals, PSM3 = PSM3(pvals), PSM4 = PSM4(pvals))

ggplot(dat, aes(x = pvals)) + 
  scale_y_continuous(limits = c(0,1.1), expand = expansion(0,0)) + 
  scale_x_continuous(breaks = seq(0,1,0.2), expand = expansion(0,0)) + 
  geom_vline(xintercept = 0.025, linetype = "dashed") + 
  geom_hline(yintercept = 0) + 
  geom_area(aes(y = PSM3), fill = "green", alpha = 0.6) +   
  theme_minimal() + 
  labs(x = "p-value (one-sided)", y = "Selection probability")

## ----two-step, echo = FALSE, fig.cap = "Two-step selection model with $\\lambda_1 = 0.4, \\lambda_2 = 0.2$"----
ggplot(dat, aes(x = pvals)) + 
  scale_y_continuous(limits = c(0,1.1), expand = expansion(0,0)) + 
  scale_x_continuous(breaks = seq(0,1,0.2), expand = expansion(0,0)) + 
  geom_vline(xintercept = c(0.025, 0.500), linetype = "dashed") + 
  geom_hline(yintercept = 0) + 
  geom_area(aes(y = PSM4), fill = "purple", alpha = 0.6) +   
  theme_minimal() + 
  labs(x = "p-value (one-sided)", y = "Selection probability")

## ----beta-one, echo = FALSE, fig.cap = "Beta-density selection model with $\\lambda_1 = 0.1, \\lambda_2 = 0.9$, using truncation points $\\alpha_1 = 0.025, \\alpha_2 = 0.975$"----

pvals <- seq(0,1,0.005)
beta_strong <- beta_fun(delta_1 = 0.1, delta_2 = 0.9, trunc_1 = 0.025, trunc_2 = 0.975)
beta_mild <- beta_fun(delta_1 = 0.7, delta_2 = 1, trunc_1 = 0.025, trunc_2 = 0.500)

dat <- data.frame(p = pvals, strong = beta_strong(pvals), mild = beta_mild(pvals))

ggplot(dat, aes(x = pvals)) + 
  scale_y_continuous(limits = c(0,1.1), expand = expansion(0,0)) + 
  scale_x_continuous(breaks = seq(0,1,0.2), expand = expansion(0,0)) + 
  geom_vline(xintercept = c(0.025, 0.975), linetype = "dashed") + 
  geom_hline(yintercept = 0) + 
  geom_area(aes(y = strong), fill = "red", alpha = 0.6) +   
  theme_minimal() + 
  labs(x = "p-value (one-sided)", y = "Selection probability")


## ----beta-two, echo = FALSE, fig.cap = "Beta-density selection model with $\\lambda_1 = 0.7, \\lambda_2 = 1$,  using truncation points $\\alpha_1 = 0.025, \\alpha_2 = 0.500$"----

ggplot(dat, aes(x = pvals)) + 
  scale_y_continuous(limits = c(0,1.1), expand = expansion(0,0)) + 
  scale_x_continuous(breaks = seq(0,1,0.2), expand = expansion(0,0)) + 
  geom_vline(xintercept = c(0.025, 0.500), linetype = "dashed") + 
  geom_hline(yintercept = 0) + 
  geom_area(aes(y = mild), fill = "yellow", alpha = 0.6) +   
  theme_minimal() + 
  labs(x = "p-value (one-sided)", y = "Selection probability")

## ----echo = FALSE-------------------------------------------------------------
data("dat.lehmann2018", package = "metadat")
n_ES <- nrow(dat.lehmann2018)
n_studies <- length(table(dat.lehmann2018$Full_Citation))

## ----lehmann------------------------------------------------------------------
data("dat.lehmann2018", package = "metadat")
dat.lehmann2018$study <- dat.lehmann2018$Full_Citation
dat.lehmann2018$sei <- sqrt(dat.lehmann2018$vi)
dat.lehmann2018$esid <- 1:nrow(dat.lehmann2018) 

## ----CHE----------------------------------------------------------------------
library(metafor)
library(clubSandwich)

# Create sampling variance-covariance matrix
V_mat <- vcalc(
  vi = vi, 
  cluster = study,
  obs = esid, 
  data = dat.lehmann2018,
  rho = 0.8,
  sparse = TRUE
)

# Fit CHE working model
CHE_mod <- rma.mv(
  yi = yi, V = V_mat,
  random = ~ 1 | study / esid,
  data = dat.lehmann2018,
  sparse = TRUE
) |>
# Apply CRVE with small-sample corrections, clustering by study
  robust(cluster = study, clubSandwich = TRUE)

CHE_mod

## ----3PSM---------------------------------------------------------------------
library(metaselection)

SMD_prior <- define_priors(
  beta_mean = 0, 
  beta_precision = 0.5, 
  beta_L = 2, 
  tau_mode = 0.2, 
  tau_alpha = 1, 
  lambda_mode = 0.8, 
  lambda_precision = 0.5, 
  lambda_L = 2
)
  
mod_3PSM <- selection_model(
  data = dat.lehmann2018, 
  yi = yi,
  sei = sei,
  cluster = study,
  selection_type = "step",
  steps = 0.025, 
  priors = SMD_prior
)

mod_3PSM

## ----echo = FALSE-------------------------------------------------------------
pct_reduction_3PSM <- 100 * (1 - mod_3PSM$est["beta","Est"] / as.numeric(CHE_mod$beta))

## -----------------------------------------------------------------------------
print(mod_3PSM, transf_gamma = FALSE, transf_zeta = FALSE)

## -----------------------------------------------------------------------------
summary(mod_3PSM)

## ----3PSM-plot, fig.cap = "P-value selection probability based on estimated one-step selection model."----
selection_plot(mod_3PSM)

## -----------------------------------------------------------------------------
mod_4PSM <- selection_model(
  data = dat.lehmann2018, 
  yi = yi,
  sei = sei,
  cluster = study,
  selection_type = "step",
  steps = c(0.025, 0.500),
  priors = SMD_prior
)

print(mod_4PSM, transf_gamma = TRUE, transf_zeta = TRUE)

## ----echo = FALSE-------------------------------------------------------------
pct_4PSM <- 100 * mod_4PSM$est["beta","Est"] / as.numeric(CHE_mod$beta)

## ----4PSM-plot, fig.cap = "P-value selection probability based on estimated two-step selection model."----
selection_plot(mod_4PSM)

## ----3PSM-mod-----------------------------------------------------------------
mod_3PSM_mod <- selection_model(
  data = dat.lehmann2018, 
  yi = yi,
  sei = sei,
  cluster = study,
  selection_type = "step",
  steps = 0.025,
  mean_mods = ~ Design,
  priors = SMD_prior
)

mod_3PSM_mod

## ----beta---------------------------------------------------------------------
mod_beta <- selection_model(
  data = dat.lehmann2018, 
  yi = yi,
  sei = sei,
  cluster = study,
  selection_type = "beta",
  steps = c(0.025, 0.975),
  priors = NULL
)

print(mod_beta)

## ----beta-plot, fig.cap = "P-value selection probability based on estimated two-step selection model."----
selection_plot(mod_beta)

## ----3PSM-bootstrap-----------------------------------------------------------
set.seed(20240916)

system.time(
  mod_3PSM_boot <- selection_model(
    data = dat.lehmann2018, 
    yi = yi,
    sei = sei,
    cluster = study,
    selection_type = "step",
    steps = 0.025,
    priors = SMD_prior,
    bootstrap = "two-stage",
    CI_type = "percentile",
    R = 199
  )
)

print(mod_3PSM_boot, transf_gamma = TRUE, transf_zeta = TRUE)

## ----parallel-boot------------------------------------------------------------
library(future)
plan(multisession, workers = 4L)

system.time(
  selection_model(
    data = dat.lehmann2018, 
    yi = yi,
    sei = sei,
    cluster = study,
    selection_type = "step",
    steps = 0.025,
    bootstrap = "two-stage",
    CI_type = "percentile",
    R = 199
  )  
)


## ----sequential-boot----------------------------------------------------------
plan(sequential)

## ----eval = FALSE-------------------------------------------------------------
# library(progressr)
# 
# with_progress(
#   sel_fit <- selection_model(
#     data = dat.lehmann2018,
#     yi = yi,
#     sei = sei,
#     cluster = study,
#     selection_type = "step",
#     steps = 0.025,
#     bootstrap = "two-stage",
#     CI_type = "percentile",
#     R = 199
#   )
# )
# 

## ----eval = FALSE-------------------------------------------------------------
# progressr::handlers(global = TRUE)