---
title: "Effect size computations: dingeszakkou2019muchstakeknowledge"
format:
  html:
    toc: true
execute:
  echo: true
  warning: true
  message: false
---

Computes standardized mean differences (`d`) and sampling variances (`v`) for the
extraction YAML `papers/dingeszakkou2019muchstakeknowledge/dingeszakkou2019muchstakeknowledge.yaml`.

For Studies 1-2, the meta-analytic effect sizes are recomputed from the published
two-group descriptives rather than taken from the paper's Tukey HSD `d`, because
the reported values align with the omnibus three-group residual SD rather than the
pairwise pooled SD used for standard between-groups SMDs. These pairwise SMDs are
computed with `esc::esc_mean_sd`, following the current project template.

## Inputs and methods

```{r}
paper_key <- "dingeszakkou2019muchstakeknowledge"
sign_convention <- "d = mean(low) - mean(high)"

effects <- list(
  list(
    study_id = 1,
    effect_id = "s1_e1",
    method_used = "between_groups",
    n_low = 51,
    n_high = 50,
    mean_low = 5.43,
    mean_high = 3.32,
    sd_low = 2.82,
    sd_high = 4.40,
    reported_d = 0.50,
    notes_on_assumptions = "Study 1 (typos): pairwise SMD recomputed from published Neutral/Stakes group descriptives. The paper's Tukey HSD d=.50 is retained under reported_test only."
  ),
  list(
    study_id = 2,
    effect_id = "s2_e1",
    method_used = "between_groups",
    n_low = 51,
    n_high = 50,
    mean_low = 5.06,
    mean_high = 1.10,
    sd_low = 3.27,
    sd_high = 5.08,
    reported_d = 1.06,
    notes_on_assumptions = "Study 2 (bank): pairwise SMD recomputed from published Neutral/Stakes group descriptives. The paper's Tukey HSD d=1.06 is retained under reported_test only."
  ),
  list(
    study_id = 3,
    effect_id = "s3_e1",
    method_used = "between_reported_d_t_df",
    t_value = 5.74,
    df = 168,
    reported_d = 0.88,
    sign_d = +1,
    notes_on_assumptions = "Study 3 (bank): d and t(df) reported; v computed from reported d + t(df). sign_d set to +1 to match the paper's claim that retraction is higher in STAKES despite the likely condition-label swap in the printed means."
  )
)

study12_denominator_checks <- list(
  list(
    study_id = 1,
    contrast = "Study 1: Neutral vs Stakes",
    n_low = 51,
    n_high = 50,
    n_other = 50,
    mean_low = 5.43,
    mean_high = 3.32,
    sd_low = 2.82,
    sd_high = 4.40,
    sd_other = 5.19,
    reported_d = 0.50
  ),
  list(
    study_id = 2,
    contrast = "Study 2: Neutral vs Stakes",
    n_low = 51,
    n_high = 50,
    n_other = 51,
    mean_low = 5.06,
    mean_high = 1.10,
    sd_low = 3.27,
    sd_high = 5.08,
    sd_other = 2.46,
    reported_d = 1.06
  )
)
```

## Shared helpers

```{r}
if (!requireNamespace("esc", quietly = TRUE)) {
  stop("Package 'esc' is required for this analysis. Install with install.packages('esc').", call. = FALSE)
}
suppressPackageStartupMessages(library(esc))

default_or <- function(x, default) {
  if (is.null(x)) default else x
}

stop_if_missing <- function(x, name) {
  if (is.na(x)) stop(sprintf("Missing required input: %s", name), call. = FALSE)
}

# Exact small-sample correction factor used by metafor (.cmicalc).
hedges_correction <- function(df) {
  ifelse(df <= 1, NA_real_, exp(lgamma(df / 2) - log(sqrt(df / 2)) - lgamma((df - 1) / 2)))
}

pooled_sd <- function(n_high, n_low, sd_high, sd_low) {
  sqrt(((n_high - 1) * sd_high^2 + (n_low - 1) * sd_low^2) / (n_high + n_low - 2))
}

d_from_groups_independent <- function(n_high, n_low, mean_high, mean_low, sd_high, sd_low) {
  s <- pooled_sd(n_high, n_low, sd_high, sd_low)
  (mean_low - mean_high) / s
}

omnibus_residual_sd_three_group <- function(n_low, n_high, n_other, sd_low, sd_high, sd_other) {
  num <- (n_low - 1) * sd_low^2 + (n_high - 1) * sd_high^2 + (n_other - 1) * sd_other^2
  den <- n_low + n_high + n_other - 3
  sqrt(num / den)
}

var_d_independent <- function(d, n_high, n_low) {
  n <- n_high + n_low
  (n / (n_high * n_low)) + (d^2 / (2 * (n - 2)))
}

var_d_between_from_d_t_df <- function(d, t_value, df) {
  (d / t_value)^2 + (d^2) / (2 * df)
}

extract_esc <- function(x) {
  list(
    d = as.numeric(x$es),
    v = as.numeric(x$var)
  )
}

compute_with_esc <- function(fun, ...) {
  d_obj <- fun(..., es.type = "d")
  g_obj <- fun(..., es.type = "g")
  d_out <- extract_esc(d_obj)
  g_out <- extract_esc(g_obj)
  list(d = d_out$d, v = d_out$v, g = g_out$d, v_g = g_out$v)
}

infer_sign <- function(mean_low, mean_high, sign_d) {
  if (!is.na(mean_low) && !is.na(mean_high) && mean_low != mean_high) {
    return(sign(mean_low - mean_high))
  }
  stop_if_missing(sign_d, "sign_d (cannot infer sign from means)")
  if (!is.element(sign_d, c(-1, 1))) stop("sign_d must be +1 or -1", call. = FALSE)
  sign_d
}
```

## Denominator check for reported d in Studies 1-2

```{r}
check_reported_d <- function(x) {
  pairwise_d <- d_from_groups_independent(
    n_high = x$n_high,
    n_low = x$n_low,
    mean_high = x$mean_high,
    mean_low = x$mean_low,
    sd_high = x$sd_high,
    sd_low = x$sd_low
  )

  omnibus_sd <- omnibus_residual_sd_three_group(
    n_low = x$n_low,
    n_high = x$n_high,
    n_other = x$n_other,
    sd_low = x$sd_low,
    sd_high = x$sd_high,
    sd_other = x$sd_other
  )

  omnibus_d <- (x$mean_low - x$mean_high) / omnibus_sd

  data.frame(
    study_id = x$study_id,
    contrast = x$contrast,
    reported_d = x$reported_d,
    pairwise_pooled_sd_d = pairwise_d,
    omnibus_residual_sd_d = omnibus_d,
    stringsAsFactors = FALSE
  )
}

denominator_check <- do.call(rbind, lapply(study12_denominator_checks, check_reported_d))
denominator_check
```

## Computation

```{r}
compute_effect <- function(effect_inputs) {
  study_id <- effect_inputs$study_id
  effect_id <- effect_inputs$effect_id
  method_used <- effect_inputs$method_used

  n_high <- default_or(effect_inputs$n_high, NA_integer_)
  n_low <- default_or(effect_inputs$n_low, NA_integer_)
  mean_high <- default_or(effect_inputs$mean_high, NA_real_)
  mean_low <- default_or(effect_inputs$mean_low, NA_real_)
  sd_high <- default_or(effect_inputs$sd_high, NA_real_)
  sd_low <- default_or(effect_inputs$sd_low, NA_real_)
  t_value <- default_or(effect_inputs$t_value, NA_real_)
  df <- default_or(effect_inputs$df, NA_real_)
  reported_d <- default_or(effect_inputs$reported_d, NA_real_)
  sign_d <- default_or(effect_inputs$sign_d, NA_real_)
  notes_on_assumptions <- default_or(effect_inputs$notes_on_assumptions, "")

  d <- NA_real_
  v <- NA_real_
  g <- NA_real_
  v_g <- NA_real_
  computed_from_suggested <- NA_character_

  if (method_used == "between_groups") {
    computed_from_suggested <- "groups"
    stop_if_missing(n_high, "n_high")
    stop_if_missing(n_low, "n_low")
    stop_if_missing(mean_high, "mean_high")
    stop_if_missing(mean_low, "mean_low")
    stop_if_missing(sd_high, "sd_high")
    stop_if_missing(sd_low, "sd_low")

    res <- compute_with_esc(
      esc::esc_mean_sd,
      grp1m = mean_low,
      grp1sd = sd_low,
      grp1n = n_low,
      grp2m = mean_high,
      grp2sd = sd_high,
      grp2n = n_high
    )
    d <- res$d
    v <- res$v
    g <- res$g
    v_g <- res$v_g
  } else if (method_used == "between_reported_d_t_df") {
    computed_from_suggested <- "reported_d"
    stop_if_missing(reported_d, "reported_d")
    stop_if_missing(t_value, "t_value")
    stop_if_missing(df, "df")
    sign_used <- infer_sign(mean_low, mean_high, sign_d)
    d <- sign_used * abs(reported_d)
    v <- var_d_between_from_d_t_df(d = d, t_value = abs(t_value), df = df)
    J <- hedges_correction(df)
    g <- J * d
    v_g <- (J^2) * v
  } else {
    stop(sprintf("Unknown method_used: %s", method_used), call. = FALSE)
  }

  inputs_used <- paste(
    c(
      sprintf("method=%s", method_used),
      sprintf("sign_convention=%s", sign_convention),
      if (!is.na(n_low)) sprintf("n_low=%s", n_low) else NULL,
      if (!is.na(n_high)) sprintf("n_high=%s", n_high) else NULL,
      if (!is.na(mean_low)) sprintf("mean_low=%s", mean_low) else NULL,
      if (!is.na(mean_high)) sprintf("mean_high=%s", mean_high) else NULL,
      if (!is.na(sd_low)) sprintf("sd_low=%s", sd_low) else NULL,
      if (!is.na(sd_high)) sprintf("sd_high=%s", sd_high) else NULL,
      if (!is.na(df)) sprintf("df=%s", df) else NULL,
      if (!is.na(t_value)) sprintf("t=%s", t_value) else NULL,
      if (!is.na(reported_d)) sprintf("reported_d=%s", reported_d) else NULL,
      if (!is.na(sign_d)) sprintf("sign_d=%s", sign_d) else NULL
    ),
    collapse = ", "
  )

  data.frame(
    paper_key = paper_key,
    study_id = study_id,
    effect_id = effect_id,
    design = "Between-Subjects",
    method_used = method_used,
    computed_from_suggested = computed_from_suggested,
    inputs_used = inputs_used,
    d = d,
    v = v,
    g = g,
    v_g = v_g,
    notes_on_assumptions = notes_on_assumptions,
    stringsAsFactors = FALSE
  )
}

audit <- do.call(rbind, lapply(effects, compute_effect))
audit
```

## Paste-ready YAML snippets

```{r}
for (i in seq_len(nrow(audit))) {
  row <- audit[i, ]
  cat(sprintf("\n# %s (study_id=%s)\n", row$effect_id, row$study_id))
  cat(sprintf(
    "effect_size:\n  metric: SMD\n  d: %.12f\n  v: %.12f\n  computed_from: %s\n  needs_review: false\n  notes: \"%s\"\n",
    row$d,
    row$v,
    row$computed_from_suggested,
    gsub("\"", "'", row$inputs_used)
  ))
}
```