34  Cohen’s d is a biased estimator of the population standardized mean difference

35 Introduction

This introduction is structured following the ADEMP framework (Aims, Data-generating process, Estimands, Methods, Performance) as described by Morris, White, and Crowther (2019) and operationalised for psychology by Siepe et al. (2024). Working through the framework explicitly makes the design decisions behind the simulation visible and auditable, so that a reader can judge whether the conclusions are supported by the conditions actually examined.

35.1 Project description

Cohen’s d, the standardized mean difference calculated by dividing the sample mean difference by the pooled sample standard deviation, is among the most widely reported effect sizes in psychology. It is also a biased estimator of the population standardized mean difference: in small samples it systematically overestimates the magnitude of the population effect. Hedges (1981) derived the exact multiplicative correction factor that removes this bias, and the corrected estimator is conventionally called Hedges’ g.

This simulation quantifies that bias across the range of sample sizes at which it matters, and verifies that Hedges’ correction removes it. The motivating empirical context is a two-arm between-subjects study with a continuous outcome, of the size routinely seen in psychology — including the small samples characteristic of pilot studies, within-lab replications, and the primary studies that feed meta-analyses, where the bias is largest and where standardized effect sizes are most often the quantity actually pooled.

The study serves a second, didactic purpose. It demonstrates that the experiment implemented in a Monte Carlo simulation often contains within-iteration factors, not only between-iteration ones. Here the estimator (uncorrected vs. corrected) is a within factor: both estimators are computed from the same simulated dataset inside analyze(), rather than appearing as a column in expand_grid(). Some factors of the experiment are therefore a property of the analysis function or of the piped workflow, not of the parameter grid. This is not merely a stylistic choice — see Methods and extracted quantities below for why it is the statistically preferable design here.

35.2 Prior related work

The bias of the standardized mean difference and its correction are established analytic results, not open questions. Hedges (1981) showed that for two independent samples with \(df = n_1 + n_2 - 2\),

\[\mathbb{E}[d] \;=\; \frac{\delta}{J(df)}, \qquad J(df) \;=\; \frac{\Gamma\!\left(\frac{df}{2}\right)}{\sqrt{df/2}\;\; \Gamma\!\left(\frac{df-1}{2}\right)} \;\approx\; 1 - \frac{3}{4\,df - 1}\]

so that \(g = J \cdot d\) is unbiased for \(\delta\). Because \(J < 1\) and \(J \to 1\) as \(df \to \infty\), the bias of d is upward, proportional to \(\delta\), and vanishes asymptotically. The correction is implemented in standard software, including effsize::cohen.d(hedges.correction = TRUE) as used here, and is recommended by meta-analytic texts (Hedges & Olkin, 1985; Borenstein et al., 2009).

We therefore do not expect novel methodological findings. The value of the simulation is threefold: it makes an abstract analytic result concrete and visible; it provides a check that the simulation infrastructure is implemented correctly, since the simulated bias can be compared against a known closed-form expectation; and it quantifies at which sample sizes the bias is large enough to matter in practice.

No preliminary simulations were conducted by the contributors on this specific question prior to this study.

35.3 Why d is biased: the sampling distribution of the standard deviation

The formula above says that Cohen’s d is biased, but not why. The mechanism is worth spelling out, because the usual one-line explanation is only about a third right.

Cohen’s d is a ratio, \(d = (M_{\text{intervention}} - M_{\text{control}}) / s_{\text{pooled}}\). The numerator is an unbiased estimator of the population mean difference at every sample size — the negative control in this study confirms it. So the entire bias must come from the denominator. For Normally distributed data the numerator and denominator are statistically independent, which lets the expectation be factored:

\[\mathbb{E}[d] \;=\; \mathbb{E}\!\left[M_{\text{intervention}} - M_{\text{control}}\right] \times \mathbb{E}\!\left[\frac{1}{s_{\text{pooled}}}\right] \;=\; (\mu_{\text{intervention}} - \mu_{\text{control}}) \times \mathbb{E}\!\left[\frac{1}{s_{\text{pooled}}}\right]\]

Everything therefore hinges on \(\mathbb{E}[1/s]\).

Step 1: the sample SD underestimates the population SD. The sample variance \(s^2\) (with the \(n-1\) denominator) is unbiased for \(\sigma^2\). The sample standard deviation \(s\) is not unbiased for \(\sigma\), because taking a square root is a concave transformation and Jensen’s inequality gives \(\mathbb{E}[s] = \mathbb{E}[\sqrt{s^2}] < \sqrt{\mathbb{E}[s^2]} = \sigma\). The plot below shows why: the sampling distribution of \(s\) is right-skewed at small \(df\), so its mean sits below \(\sigma\).

# this chunk is self-contained so that the mechanism can be shown here in the
# introduction, before the simulation's own dependencies are loaded under Methods
library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)
library(knitr)
library(kableExtra)

# exact sampling density of the pooled SD. for Normal data, df * s^2 / sigma^2 follows
# a chi-square distribution with df degrees of freedom, so the density of s itself is
# obtained by change of variables. using the exact density rather than simulating keeps
# the figure free of Monte Carlo noise
density_of_sd <- function(s, df, sigma = 1) dchisq(df * s^2 / sigma^2, df) * 2 * df * s / sigma^2

# E[s] and E[1/s] in closed form (both are gamma-function ratios)
expected_sd     <- function(df) sqrt(2 / df) * gamma((df + 1) / 2) / gamma(df / 2)
expected_inv_sd <- function(df) sqrt(df / 2) * gamma((df - 1) / 2) / gamma(df / 2)

ns_to_show <- c(5, 10, 25, 50)
panel_levels <- paste0("n = ", ns_to_show, " per group")

sd_densities <- expand_grid(n_per_condition = ns_to_show,
                            s = seq(0.01, 2.2, length.out = 800)) |>
  mutate(df      = 2 * n_per_condition - 2,
         density = density_of_sd(s, df),
         panel   = factor(paste0("n = ", n_per_condition, " per group"), levels = panel_levels))

sd_means <- tibble(n_per_condition = ns_to_show) |>
  mutate(df     = 2 * n_per_condition - 2,
         mean_s = expected_sd(df),
         panel  = factor(paste0("n = ", n_per_condition, " per group"), levels = panel_levels))

ggplot(sd_densities, aes(x = s, y = density)) +
  geom_area(fill = "grey70", alpha = 0.6) +
  geom_vline(xintercept = 1, linetype = "dashed") +
  geom_vline(data = sd_means, aes(xintercept = mean_s), color = "red") +
  facet_wrap(~ panel, scales = "free_y") +
  scale_x_continuous(name = "Pooled sample SD", breaks = breaks_pretty(n = 5)) +
  scale_y_continuous(name = "Density", breaks = NULL) +
  theme_linedraw() +
  theme(panel.grid.minor = element_blank()) +
  ggtitle("Sampling distribution of the pooled SD when the population SD is 1",
          subtitle = "Dashed line = population SD of 1. Red line = mean of the sampling distribution.")

At n = 5 per group the distribution is visibly right-skewed and its mean sits at 0.969 rather than 1. By n = 50 per group it is nearly symmetric and centred almost exactly on 1. This is the part of the story most people remember: the SD is underestimated in small samples.

Step 2: but d divides by s, and that makes it worse. The intuition “the SD is about 3% too small, so d is about 3% too big” gets the direction right but substantially understates the magnitude, because \(\mathbb{E}[1/s] \neq 1/\mathbb{E}[s]\). Since \(1/x\) is a convex function, Jensen’s inequality now runs the other way and gives \(\mathbb{E}[1/s] > 1/\mathbb{E}[s]\). Intuitively: \(1/s\) blows up whenever \(s\) happens to come out small, so the left tail of the SD’s sampling distribution has an outsized effect on the average of its reciprocal.

# density of 1/s, again by change of variables from the density of s
inv_densities <- expand_grid(n_per_condition = ns_to_show,
                             inv_s = seq(0.25, 3.5, length.out = 800)) |>
  mutate(df      = 2 * n_per_condition - 2,
         density = density_of_sd(1 / inv_s, df) / inv_s^2,
         panel   = factor(paste0("n = ", n_per_condition, " per group"), levels = panel_levels))

inv_means <- tibble(n_per_condition = ns_to_show) |>
  mutate(df         = 2 * n_per_condition - 2,
         mean_inv_s = expected_inv_sd(df),
         panel      = factor(paste0("n = ", n_per_condition, " per group"), levels = panel_levels))

ggplot(inv_densities, aes(x = inv_s, y = density)) +
  geom_area(fill = "grey70", alpha = 0.6) +
  geom_vline(xintercept = 1, linetype = "dashed") +
  geom_vline(data = inv_means, aes(xintercept = mean_inv_s), color = "red") +
  facet_wrap(~ panel, scales = "free_y") +
  scale_x_continuous(name = "1 / pooled sample SD", breaks = breaks_pretty(n = 5)) +
  scale_y_continuous(name = "Density", breaks = NULL) +
  theme_linedraw() +
  theme(panel.grid.minor = element_blank()) +
  ggtitle("Sampling distribution of 1 / pooled SD when the population SD is 1",
          subtitle = "Dashed line = 1. Red line = mean of the sampling distribution: it sits clearly above 1.")