## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----echo=FALSE, message=FALSE, warning=FALSE, results='hide'-----------------
library(metaConvert)
library(metafor)
library(DT)

## -----------------------------------------------------------------------------
data(df.psychom)

## ----eval=FALSE---------------------------------------------------------------
# str(df.psychom)

## ----echo=FALSE---------------------------------------------------------------
DT::datatable(df.psychom, options = list(
    scrollX = TRUE,
    dom = c('t'),
    scrollY = "300px",
    pageLength = 50,
    ordering = FALSE,
    rownames = FALSE,
    columnDefs = list(
                  list(width = '130px',
                       targets = "_all"),
                  list(className = 'dt-center',
                                     targets = "_all"))))

## ----eval=FALSE---------------------------------------------------------------
# # Extraction sheet for Cronbach's alpha
# data_extraction_sheet(measure = "alpha", extension = "data.frame")
# 
# # Extraction sheet for ICC
# data_extraction_sheet(measure = "icc", extension = "data.frame")
# 
# # Extraction sheet for correlations (Pearson or Spearman)
# data_extraction_sheet(measure = "r", extension = "data.frame")
# 
# # Extraction sheet for proportions
# data_extraction_sheet(measure = "prop", extension = "data.frame")

## -----------------------------------------------------------------------------
# Subset internal consistency studies
dat_alpha <- df.psychom[df.psychom$outcome == "internal_consistency", ]

# Compute effect sizes with Bonett transformation (default)
res_alpha <- convert_df(dat_alpha, measure = "alpha",
                        verbose = FALSE, split_adjusted = FALSE)
summary_alpha <- summary(res_alpha, flags = TRUE, guidance = TRUE)

## ----eval=FALSE---------------------------------------------------------------
# summary_alpha[, c("author", "es", "se", "es_ci_lo", "es_ci_up", "info_used")]

## ----echo=FALSE---------------------------------------------------------------
DT::datatable(summary_alpha[, c("author", "es", "se", "es_ci_lo", "es_ci_up", "info_used")],
              options = list(scrollX = TRUE, dom = 't', ordering = FALSE,
                             rownames = FALSE, pageLength = 50,
                             columnDefs = list(
                               list(width = '130px', targets = "_all"),
                               list(className = 'dt-center', targets = "_all"))))

## ----fig.width=7, fig.height=4------------------------------------------------
# Random-effects meta-analysis on Bonett-transformed scale
ma_alpha <- metafor::rma(yi = es, sei = se, data = summary_alpha,
                         method = "REML")

# Forest plot
metafor::forest(ma_alpha, slab = paste0(summary_alpha$author, " (", dat_alpha$year, ")"),
                xlab = "Bonett-transformed alpha: ln(1 - alpha)",
                header = TRUE)

## -----------------------------------------------------------------------------
# Back-transform pooled estimate and CI
pooled_alpha    <- 1 - exp(ma_alpha$beta[1])
pooled_alpha_lo <- 1 - exp(ma_alpha$ci.ub)
pooled_alpha_up <- 1 - exp(ma_alpha$ci.lb)

cat(sprintf("Pooled alpha = %.3f [%.3f, %.3f]\n",
            pooled_alpha, pooled_alpha_lo, pooled_alpha_up))

## ----eval=FALSE---------------------------------------------------------------
# # Raw alpha (not recommended for meta-analysis)
# res_alpha_raw <- convert_df(dat_alpha, measure = "alpha", alpha_to_es = "raw",
#                             verbose = FALSE, split_adjusted = FALSE)
# summary(res_alpha_raw)

## -----------------------------------------------------------------------------
# Subset test-retest studies
dat_icc <- df.psychom[df.psychom$outcome == "test_retest", ]

# Compute effect sizes with Bonett transformation (default)
res_icc <- convert_df(dat_icc, measure = "icc",
                      verbose = FALSE, split_adjusted = FALSE)
summary_icc <- summary(res_icc, flags = TRUE)

## ----eval=FALSE---------------------------------------------------------------
# summary_icc[, c("author", "es", "se", "es_ci_lo", "es_ci_up",
#                  "icc_type", "info_used")]

## ----echo=FALSE---------------------------------------------------------------
DT::datatable(summary_icc[, c("author", "es", "se", "es_ci_lo", "es_ci_up",
                                "icc_type", "info_used")],
              options = list(scrollX = TRUE, dom = 't', ordering = FALSE,
                             rownames = FALSE, pageLength = 50,
                             columnDefs = list(
                               list(width = '130px', targets = "_all"),
                               list(className = 'dt-center', targets = "_all"))))

## ----fig.width=7, fig.height=4------------------------------------------------
# Random-effects meta-analysis
ma_icc <- metafor::rma(yi = es, sei = se, data = summary_icc,
                       method = "REML")

# Forest plot
metafor::forest(ma_icc, slab = paste0(summary_icc$author, " (", dat_icc$year, ")"),
                xlab = "Bonett-transformed ICC: ln(1 - ICC)",
                header = TRUE)

## -----------------------------------------------------------------------------
# Back-transform pooled estimate
pooled_icc    <- 1 - exp(ma_icc$beta[1])
pooled_icc_lo <- 1 - exp(ma_icc$ci.ub)
pooled_icc_up <- 1 - exp(ma_icc$ci.lb)

cat(sprintf("Pooled ICC = %.3f [%.3f, %.3f]\n",
            pooled_icc, pooled_icc_lo, pooled_icc_up))

## -----------------------------------------------------------------------------
# Use test-retest studies that report both ICC and SD
dat_me <- df.psychom[df.psychom$outcome == "test_retest" &
                     !is.na(df.psychom$sd_scores), ]

# Compute SEM for each study (delta-method SE included)
sem_results <- compute_sem(
  sd       = dat_me$sd_scores,
  icc      = dat_me$icc,
  n_sample = dat_me$n_sample
)

# Chain to SDC
sdc_results <- compute_sdc(
  sem    = sem_results$sem,
  sem_se = sem_results$sem_se
)

## ----eval=FALSE---------------------------------------------------------------
# # Combined measurement error table
# measurement_error <- data.frame(
#   study  = dat_me$study_id,
#   n      = dat_me$n_sample,
#   sd     = dat_me$sd_scores,
#   icc    = dat_me$icc,
#   sem    = round(sem_results$sem, 3),
#   sem_se = round(sem_results$sem_se, 3),
#   sdc    = round(sdc_results$sdc, 3),
#   sdc_se = round(sdc_results$sdc_se, 3)
# )
# measurement_error

## ----echo=FALSE---------------------------------------------------------------
measurement_error <- data.frame(
  study  = dat_me$study_id,
  n      = dat_me$n_sample,
  sd     = dat_me$sd_scores,
  icc    = dat_me$icc,
  sem    = round(sem_results$sem, 3),
  sem_se = round(sem_results$sem_se, 3),
  sdc    = round(sdc_results$sdc, 3),
  sdc_se = round(sdc_results$sdc_se, 3)
)
DT::datatable(measurement_error,
              options = list(scrollX = TRUE, dom = 't', ordering = FALSE,
                             rownames = FALSE, pageLength = 50,
                             columnDefs = list(
                               list(width = '100px', targets = "_all"),
                               list(className = 'dt-center', targets = "_all"))))

## ----fig.width=7, fig.height=4------------------------------------------------
ma_sem <- metafor::rma(yi = sem_results$sem, sei = sem_results$sem_se,
                       method = "REML")

metafor::forest(ma_sem,
                slab = paste0(dat_me$author, " (", dat_me$year, ")"),
                xlab = "Standard Error of Measurement (SEM)",
                header = TRUE)

## -----------------------------------------------------------------------------
cat(sprintf("Pooled SEM = %.3f [%.3f, %.3f]\n",
            ma_sem$beta[1], ma_sem$ci.lb, ma_sem$ci.ub))

# Derive pooled SDC from pooled SEM
pooled_sdc <- 1.96 * sqrt(2) * ma_sem$beta[1]
cat(sprintf("Pooled SDC = %.3f\n", pooled_sdc))

## -----------------------------------------------------------------------------
# Subset criterion and construct validity studies
dat_validity <- df.psychom[df.psychom$outcome %in%
                           c("criterion_validity", "construct_validity"), ]

# convert_df automatically handles both Pearson r and Spearman rho
res_r <- convert_df(dat_validity, measure = "r",
                    verbose = FALSE, split_adjusted = FALSE)
summary_r <- summary(res_r, flags = TRUE)

## ----eval=FALSE---------------------------------------------------------------
# summary_r[, c("author", "es", "se", "info_used")]

## ----echo=FALSE---------------------------------------------------------------
DT::datatable(summary_r[, c("author", "es", "se", "es_ci_lo", "es_ci_up", "info_used")],
              options = list(scrollX = TRUE, dom = 't', ordering = FALSE,
                             rownames = FALSE, pageLength = 50,
                             columnDefs = list(
                               list(width = '130px', targets = "_all"),
                               list(className = 'dt-center', targets = "_all"))))

## ----fig.width=7, fig.height=5------------------------------------------------
# Pool observed (uncorrected) correlations on Fisher's z scale
z_obs <- atanh(summary_r$es)
r_bounded <- pmin(pmax(summary_r$es, -0.9999), 0.9999)
z_se  <- summary_r$se / (1 - r_bounded^2)

ma_r_obs <- metafor::rma(yi = z_obs, sei = z_se, method = "REML")

# Forest plot on correlation scale
metafor::forest(ma_r_obs,
                slab = paste0(summary_r$author, " (", dat_validity$year, ")"),
                xlab = "Fisher's z (observed correlation)",
                header = TRUE)

## -----------------------------------------------------------------------------
# Back-transform to correlation scale
pooled_r_obs <- tanh(ma_r_obs$beta[1])
pooled_r_lo  <- tanh(ma_r_obs$ci.lb)
pooled_r_up  <- tanh(ma_r_obs$ci.ub)

cat(sprintf("(a) Pooled uncorrected r = %.3f [%.3f, %.3f]\n",
            pooled_r_obs, pooled_r_lo, pooled_r_up))

## -----------------------------------------------------------------------------
# Disattenuate using study-specific reliabilities
corrected <- es_disattenuate(
  r             = summary_r$es,
  r_se          = summary_r$se,
  reliability_x = dat_validity$reliability_target,
  reliability_y = dat_validity$reliability_comparator,
  n_sample      = dat_validity$n_sample
)

## ----eval=FALSE---------------------------------------------------------------
# # Comparison table
# data.frame(
#   study              = dat_validity$study_id,
#   info_used          = summary_r$info_used,
#   r_observed         = round(summary_r$es, 3),
#   r_corrected        = round(corrected$r_corrected, 3),
#   attenuation_factor = round(corrected$attenuation_factor, 3)
# )

## ----echo=FALSE---------------------------------------------------------------
comp <- data.frame(
  study              = dat_validity$study_id,
  info_used          = summary_r$info_used,
  r_observed         = round(summary_r$es, 3),
  r_corrected        = round(corrected$r_corrected, 3),
  attenuation_factor = round(corrected$attenuation_factor, 3)
)
DT::datatable(comp,
              options = list(scrollX = TRUE, dom = 't', ordering = FALSE,
                             rownames = FALSE, pageLength = 50,
                             columnDefs = list(
                               list(width = '130px', targets = "_all"),
                               list(className = 'dt-center', targets = "_all"))))

## ----fig.width=7, fig.height=5------------------------------------------------
# Pool corrected correlations on Fisher's z scale
ma_r_corr <- metafor::rma(yi = corrected$z_corrected,
                           sei = corrected$z_corrected_se,
                           method = "REML")

metafor::forest(ma_r_corr,
                slab = paste0(summary_r$author, " (", dat_validity$year, ")"),
                xlab = "Fisher's z (disattenuated correlation)",
                header = TRUE)

## -----------------------------------------------------------------------------
pooled_r_corr <- tanh(ma_r_corr$beta[1])
pooled_r_corr_lo <- tanh(ma_r_corr$ci.lb)
pooled_r_corr_up <- tanh(ma_r_corr$ci.ub)

cat(sprintf("(b) Pooled disattenuated r = %.3f [%.3f, %.3f]\n",
            pooled_r_corr, pooled_r_corr_lo, pooled_r_corr_up))

## ----eval=FALSE---------------------------------------------------------------
# # Use pooled alpha and ICC from earlier sections as proxies
# es_disattenuate(
#   r = 0.65,
#   r_se = 0.08,
#   reliability_x = pooled_alpha,  # from Section 2
#   reliability_y = 0.85,          # from external published data
#   n_sample = 150
# )

## -----------------------------------------------------------------------------
# Subset responsiveness studies
dat_resp <- df.psychom[df.psychom$outcome == "responsiveness", ]

# Step 1: Compute change-score reliability
change_rel <- reliability_change_score(
  reliability = dat_resp$single_occasion_reliability,
  r_pre_post  = dat_resp$r_pre_post
)

# Step 2: Get responsiveness correlations via convert_df
res_resp <- convert_df(dat_resp, measure = "r",
                       verbose = FALSE, split_adjusted = FALSE)
summary_resp <- summary(res_resp)

# Step 3: Disattenuate using CHANGE-SCORE reliability
corrected_resp <- es_disattenuate(
  r             = summary_resp$es,
  r_se          = summary_resp$se,
  reliability_x = change_rel$rel_change,
  reliability_y = dat_resp$reliability_comparator,
  n_sample      = dat_resp$n_sample
)

## ----eval=FALSE---------------------------------------------------------------
# data.frame(
#   study        = dat_resp$study_id,
#   r_observed   = round(summary_resp$es, 3),
#   rel_single   = dat_resp$single_occasion_reliability,
#   r_pre_post   = dat_resp$r_pre_post,
#   rel_change   = round(change_rel$rel_change, 3),
#   r_corrected  = round(corrected_resp$r_corrected, 3)
# )

## ----echo=FALSE---------------------------------------------------------------
resp_table <- data.frame(
  study        = dat_resp$study_id,
  r_observed   = round(summary_resp$es, 3),
  rel_single   = dat_resp$single_occasion_reliability,
  r_pre_post   = dat_resp$r_pre_post,
  rel_change   = round(change_rel$rel_change, 3),
  r_corrected  = round(corrected_resp$r_corrected, 3)
)
DT::datatable(resp_table,
              options = list(scrollX = TRUE, dom = 't', ordering = FALSE,
                             rownames = FALSE, pageLength = 50,
                             columnDefs = list(
                               list(width = '110px', targets = "_all"),
                               list(className = 'dt-center', targets = "_all"))))

## ----fig.width=7, fig.height=4------------------------------------------------
# Pool observed change-score correlations
z_resp     <- atanh(summary_resp$es)
r_resp_bounded <- pmin(pmax(summary_resp$es, -0.9999), 0.9999)
z_resp_se  <- summary_resp$se / (1 - r_resp_bounded^2)

ma_resp_obs <- metafor::rma(yi = z_resp, sei = z_resp_se, method = "REML")

cat(sprintf("Pooled observed r = %.3f [%.3f, %.3f]\n",
            tanh(ma_resp_obs$beta[1]),
            tanh(ma_resp_obs$ci.lb),
            tanh(ma_resp_obs$ci.ub)))

# Pool disattenuated change-score correlations
ma_resp_corr <- metafor::rma(yi = corrected_resp$z_corrected,
                              sei = corrected_resp$z_corrected_se,
                              method = "REML")

cat(sprintf("Pooled disattenuated r = %.3f [%.3f, %.3f]\n",
            tanh(ma_resp_corr$beta[1]),
            tanh(ma_resp_corr$ci.lb),
            tanh(ma_resp_corr$ci.ub)))

## -----------------------------------------------------------------------------
# Example: one study with r_observed = 0.50, rel_single = 0.85
r_pre_post_grid <- c(0.20, 0.40, 0.60, 0.80)

sensitivity <- data.frame(
  r_pre_post = r_pre_post_grid,
  rel_change = sapply(r_pre_post_grid, function(rpp) {
    reliability_change_score(reliability = 0.85, r_pre_post = rpp)$rel_change
  }),
  r_corrected = sapply(r_pre_post_grid, function(rpp) {
    rc <- reliability_change_score(reliability = 0.85, r_pre_post = rpp)$rel_change
    es_disattenuate(r = 0.50, r_se = 0.08, reliability_x = rc,
                    reliability_y = 0.80, n_sample = 100)$r_corrected
  })
)
sensitivity

## -----------------------------------------------------------------------------
# Subset floor/ceiling studies
dat_prop <- df.psychom[df.psychom$outcome == "floor_ceiling", ]

# Freeman-Tukey transformation (recommended)
res_prop <- convert_df(dat_prop, measure = "prop", prop_to_es = "freeman_tukey",
                       verbose = FALSE, split_adjusted = FALSE)
summary_prop <- summary(res_prop)

## ----fig.width=7, fig.height=3.5----------------------------------------------
ma_prop <- metafor::rma(yi = es, sei = se, data = summary_prop,
                        method = "REML")

metafor::forest(ma_prop,
                slab = paste0(summary_prop$author, " (", dat_prop$year, ")"),
                xlab = "Freeman-Tukey transformed proportion",
                header = TRUE)

## -----------------------------------------------------------------------------
# Approximate back-transformation (for the pooled estimate)
pooled_prop <- sin(ma_prop$beta[1])^2
cat(sprintf("Pooled proportion (approx.) = %.3f\n", pooled_prop))

