5  Understanding Standardized Mean Differences (Cohen’s d) ✎ Draft

5.1 Means and Standard Deviations

Standardised effect sizes are built from means and standard deviations, so it helps to have a feel for what each of them describes. The plot below shows a single sample as a normal distribution. Use the sliders to change its mean and its SD. The mean sets where the distribution sits, and is marked by the dashed line. The SD sets how spread out the distribution is: the bar marks the scores within one SD of the mean.

viewof oneMean = Inputs.range([1, 7], {value: 4, step: 0.05, label: "Mean"})
viewof oneSd = Inputs.range([0.5, 2], {value: 1, step: 0.05, label: "SD"})

{
  // one normal distribution
  const density = x => Math.exp(-(((x - oneMean) / oneSd) ** 2) / 2) / (oneSd * Math.sqrt(2 * Math.PI))
  // both axes are fixed, so that the distribution changes when the sliders move and the axes do not
  const xMin = -4, xMax = 12, yMax = 1.15 / (0.5 * Math.sqrt(2 * Math.PI))
  const peak = density(oneMean), shoulder = density(oneMean + oneSd)
  const curve = d3.range(xMin, xMax + 0.001, 0.05).map(x => ({x, y: density(x)}))
  const blue = "#2780e3"
  return Plot.plot({
    width, height: Math.min(360, Math.max(260, width * 0.5)), marginLeft: 20, marginRight: 20, marginTop: 10, marginBottom: 45,
    x: {label: "Score", labelAnchor: "center", labelOffset: 38, domain: [xMin, xMax], ticks: d3.range(-4, 13, 2)},
    y: {axis: null, domain: [0, yMax]},
    style: {fontSize: "13px"},
    marks: [
      Plot.areaY(curve, {x: "x", y1: 0, y2: "y", fill: blue, fillOpacity: 0.2}),
      Plot.lineY(curve, {x: "x", y: "y", stroke: blue, strokeWidth: 2}),
      Plot.ruleY([0], {stroke: "currentColor"}),
      // the mean
      Plot.link([{}], {x1: oneMean, x2: oneMean, y1: 0, y2: peak, stroke: blue, strokeDasharray: "4,3"}),
      Plot.text([{}], {x: oneMean, y: peak, text: () => `M = ${oneMean.toFixed(2)}`, dy: -10, fontWeight: 600}),
      // one SD either side of the mean, drawn where the curve is one SD from the mean
      Plot.link([{}], {x1: oneMean - oneSd, x2: oneMean + oneSd, y1: shoulder, y2: shoulder, stroke: "currentColor", strokeWidth: 1.5}),
      Plot.text([{}], {x: oneMean + oneSd, y: shoulder, text: () => `SD = ${oneSd.toFixed(2)}`, textAnchor: "start", dx: 8, fontWeight: 600})
    ]
  })
}

A normal distribution has no minimum or maximum, but most scales in psychology do: a 1–7 Likert item cannot produce a response of 0 or of 8. Set the mean to 6 and the SD to 2, and a large part of the distribution lies above 7. A sample on a 1–7 scale could not look like this. On a bounded scale, the closer the mean is to either end, the smaller the largest possible SD becomes, which is taken up in Section 8.1.

5.2 Standardized Mean Differences (Cohen’s d)

Cohen’s d for two independent groups is the difference between their means divided by their pooled standard deviation:

\[ d = \frac{M_1 - M_2}{SD_{pooled}}, \qquad SD_{pooled} = \sqrt{\frac{(n_1 - 1)SD_1^2 + (n_2 - 1)SD_2^2}{n_1 + n_2 - 2}} \]

Because Cohen’s d is a difference between two means in SD units, it can be pictured as the distance between two distributions relative to how wide they are. The plot below shows a control condition and an intervention condition as two normal distributions with the same SD. The three sliders are linked by d = mean difference / SD. Change Cohen’s d and the mean difference changes to match; change the mean difference or the SD and Cohen’s d changes to match. See how far apart the two distributions are, and how much they overlap, for effect sizes of different magnitudes, and notice that the same d can come from a small difference between narrow distributions or a large difference between wide ones. Two other ways of expressing the same effect size are shown under the plot: the percentage of the two distributions that overlaps, and the probability of superiority, which is the probability that a randomly chosen participant from the intervention condition scores higher than a randomly chosen participant from the control condition. Both assume normal distributions with equal SDs. This plot is a simplified version of the interactive visualisation by Kristoffer Magnusson at rpsychologist.com/cohend.

viewof dLinked = {
  // three linked sliders: Cohen's d = mean difference / SD. changing d changes the mean difference, and changing
  // the mean difference or the SD changes d. the value, {d, difference, sd}, is stored on the wrapping div
  const dMax = 3, sdMax = 2
  const round = x => Math.round(x * 100) / 100
  const dInput = Inputs.range([0, dMax], {value: 0.5, step: 0.01, label: "Cohen's d"})
  const differenceInput = Inputs.range([0, dMax * sdMax], {value: 0.5, step: 0.01, label: "Mean difference"})
  const sdInput = Inputs.range([0.5, sdMax], {value: 1, step: 0.01, label: "SD (both conditions)"})
  const view = html`<div>${dInput}${differenceInput}${sdInput}</div>`
  const store = () => { view.value = {d: dInput.value, difference: differenceInput.value, sd: sdInput.value} }
  // each slider's own listener runs before its input event reaches the wrapping div, so the other sliders are
  // already up to date when the plot reads the value
  dInput.addEventListener("input", () => { differenceInput.value = round(dInput.value * sdInput.value); store() })
  const updateD = () => {
    // the mean difference cannot be more than the largest d on the slider allows for this SD
    if (differenceInput.value > dMax * sdInput.value) differenceInput.value = round(dMax * sdInput.value)
    dInput.value = round(differenceInput.value / sdInput.value)
    store()
  }
  differenceInput.addEventListener("input", updateD)
  sdInput.addEventListener("input", updateD)
  store()
  return view
}

{
  // two normal distributions with the same SD, whose means are the mean difference apart
  const {d, difference, sd} = dLinked
  const groups = [{group: "Control", mean: 0}, {group: "Intervention", mean: difference}]
  const density = (x, mean, sd) => Math.exp(-(((x - mean) / sd) ** 2) / 2) / (sd * Math.sqrt(2 * Math.PI))
  // both axes are fixed, so that the distributions change when the sliders move and the axes do not
  const xMin = -5, xMax = 11, yMax = 1.15 * density(0, 0, 0.5)
  const peak = density(0, 0, sd)
  const curves = groups.flatMap(g => d3.range(xMin, xMax + 0.001, 0.05).map(x => ({...g, x, y: density(x, g.mean, sd)})))
  return Plot.plot({
    width, height: Math.min(360, Math.max(260, width * 0.5)), marginLeft: 20, marginRight: 20, marginTop: 10, marginBottom: 45,
    x: {label: "Score, relative to the control condition's mean", labelAnchor: "center", labelOffset: 38, domain: [xMin, xMax], ticks: d3.range(-4, 11, 2)},
    y: {axis: null, domain: [0, yMax]},
    color: {domain: groups.map(g => g.group), range: ["#d9622b", "#2780e3"], legend: true},
    style: {fontSize: "13px"},
    marks: [
      Plot.areaY(curves, {x: "x", y1: 0, y2: "y", fill: "group", fillOpacity: 0.2}),
      Plot.lineY(curves, {x: "x", y: "y", stroke: "group", strokeWidth: 2}),
      Plot.ruleY([0], {stroke: "currentColor"}),
      // each condition's mean
      Plot.link(groups, {x1: "mean", x2: "mean", y1: 0, y2: peak, stroke: "group", strokeDasharray: "4,3"}),
      // a bracket between the two means, labelled with d
      Plot.link([{}], {x1: 0, x2: difference, y1: peak + 0.06 * yMax, y2: peak + 0.06 * yMax, stroke: "currentColor", strokeWidth: 1.5}),
      Plot.text([{}], {x: difference / 2, y: peak + 0.06 * yMax, text: () => `d = ${d.toFixed(2)}`, dy: -10, fontWeight: 600})
    ]
  })
}

{
  // overlap: the area shared by the two distributions. probability of superiority: the probability that a
  // randomly chosen intervention participant scores higher than a randomly chosen control participant
  const d = dLinked.difference / dLinked.sd
  const overlap = 2 * normCdf(-Math.abs(d) / 2)
  const superiority = normCdf(d / Math.SQRT2)
  const pct = p => `${(100 * p).toFixed(1)}%`
  return html`<div class="sb-tiles">
    <div><div class="sb-tile-value">${pct(overlap)}</div><div class="sb-tile-label">Overlap</div><div class="sb-tile-note">of the two distributions is shared</div></div>
    <div><div class="sb-tile-value">${pct(superiority)}</div><div class="sb-tile-label">Probability of superiority</div><div class="sb-tile-note">chance that a randomly chosen participant from the intervention condition scores higher than a randomly chosen participant from the control condition</div></div>
  </div>`
}

Even effect sizes that are conventionally described as large leave the two distributions overlapping substantially: at d = 0.8, about 21% of participants in the control condition score higher than the average participant in the intervention condition. The picture changes quickly as d increases. At d = 2, only about 2% of the control condition scores above the intervention condition’s mean, and at d = 3 the two distributions overlap only in their tails. This is one reason why very large reported effect sizes deserve scrutiny: they describe a situation in which almost every participant in one condition scored higher than the average participant in the other.

function unitsPlot({scales, groups, captions, n, fixedViews, scaleUnits = "Scale units", width}) {
  // each scale's logical range, with the distributions of two groups drawn on it as beta curves, and buttons
  // that re-express everything in scale units or SD units. a second panel shows the difference
  // between the two groups' means (second minus first) in the current units, with its 95% CI.
  // scales: [{name, min, max, sd, means}], with one mean per group and a common SD. groups: [{name, colour}].
  // captions: one sentence per view. n: the sample size of each group, used for the CIs. fixedViews: optionally,
  // the views that the main panel and the difference panel always show, as [main, difference] (0 for scale
  // units, 1 for SD units), in which case there are no buttons and the first caption is used. scaleUnits: the
  // label of the axis in scale units. a scale can also have a sub, a short label shown under its name.
  // each view re-expresses a score v on scale s. aligned names the one quantity that the view makes the same
  // on every scale, which is highlighted in the labels
  const views = [
    {button: "(Unstandarized) Scale units", axis: scaleUnits, digits: 0, aligned: null, f: (s, v) => v,
     difference: "Difference between means with 95% CI, in scale points",
     // mean difference: t distribution
     ci: s => { const margin = tQuantile(0.975, df) * s.sd * Math.sqrt(2 / n); return [s.means[1] - s.means[0] - margin, s.means[1] - s.means[0] + margin] }},
    {button: "Standardized (SD) units", axis: "Standardized Mean Difference (Cohen's d)", digits: 2, aligned: "sd", f: (s, v) => v / s.sd,
     difference: "Difference between means with 95% CI, in SD units (Cohen's d)",
     // Cohen's d: the bounds are the noncentrality parameters under which the observed t would sit at the 97.5th
     // and 2.5th percentiles of the noncentral t distribution, converted from t to d
     ci: s => {
       const t = (s.means[1] - s.means[0]) / (s.sd * Math.sqrt(2 / n))
       const ncp = p => bisect(x => 1 - nctCdf(t, df, x), p, -10 - 2 * Math.abs(t), 10 + 2 * Math.abs(t))
       return [ncp(0.025) * Math.sqrt(2 / n), ncp(0.975) * Math.sqrt(2 / n)]
     }}
  ]
  const df = 2 * n - 2
  const blue = "#2780e3", grey = "#8a8a8a"
  const margin = {top: 10, right: 20, bottom: 45, left: 70}

  // main panel. each row's baseline is the scale's range, and the curves rise from it
  const rowHeight = 100, curveHeight = 52, height = margin.top + rowHeight * scales.length + margin.bottom
  const svg = d3.create("svg").attr("class", "sb-std").attr("width", width).attr("height", height).attr("viewBox", [0, 0, width, height])
  const axisG = svg.append("g").attr("transform", `translate(0, ${height - margin.bottom})`)
  const axisLabel = svg.append("text").attr("x", (margin.left + width - margin.right) / 2).attr("y", height - 8).attr("text-anchor", "middle").attr("font-size", 13).attr("fill", "currentColor")
  const rows = svg.selectAll("g.sb-row").data(scales).join("g").attr("class", "sb-row").attr("transform", (s, i) => `translate(0, ${margin.top + rowHeight * i + 84})`)
  // a row's name, with the optional sub-label on a second line beneath it
  const rowName = s => s.sub
    ? `<tspan x="${margin.left - 12}" dy="-0.3em">${s.name}</tspan><tspan x="${margin.left - 12}" dy="1.25em" font-size="11" font-weight="400">${s.sub}</tspan>`
    : `<tspan x="${margin.left - 12}" dy="0.35em">${s.name}</tspan>`
  rows.append("text").attr("y", -curveHeight / 2).attr("text-anchor", "end").attr("font-size", 13).attr("font-weight", 600).attr("fill", "currentColor").html(rowName)
  const labels = rows.append("text").attr("x", margin.left).attr("y", -curveHeight - 14).attr("font-size", 12).attr("fill", "currentColor")
  // each group's distribution: a filled beta curve on the range of the scale, with a dashed line at its mean
  const marks = rows.selectAll("g.sb-group").data(s => s.means.map((mean, j) => ({s, mean, colour: groups[j].colour, density: boundedDensity(mean, s.sd, s.min, s.max)}))).join("g").attr("class", "sb-group")
  const fill = marks.append("path").attr("fill", d => d.colour).attr("fill-opacity", 0.2)
  const outline = marks.append("path").attr("fill", "none").attr("stroke", d => d.colour).attr("stroke-width", 2)
  const meanLine = marks.append("line").attr("y1", 0).attr("stroke", d => d.colour).attr("stroke-dasharray", "4,3")
  // logical range of the scale, with a cap at each end
  const range = rows.append("line").attr("stroke", "currentColor").attr("stroke-opacity", 0.5).attr("stroke-width", 2)
  const capLo = rows.append("line").attr("y1", -7).attr("y2", 7).attr("stroke", "currentColor").attr("stroke-opacity", 0.5).attr("stroke-width", 2)
  const capHi = rows.append("line").attr("y1", -7).attr("y2", 7).attr("stroke", "currentColor").attr("stroke-opacity", 0.5).attr("stroke-width", 2)
  // the curves are drawn at the same positions along each scale in every view, so that they can slide between views
  const us = d3.range(0, 1.0001, 0.01)
  const marksData = marks.data()

  // difference panel
  // dBottom is the space under the difference panel's plotting area, for its axis and axis label
  const dRowHeight = 40, dBottom = 48, dHeight = 18 + dRowHeight * scales.length + dBottom
  const dSvg = d3.create("svg").attr("class", "sb-std").attr("width", width).attr("height", dHeight).attr("viewBox", [0, 0, width, dHeight])
  const dAxisG = dSvg.append("g").attr("transform", `translate(0, ${dHeight - dBottom})`)
  const dAxisLabel = dSvg.append("text").attr("x", (margin.left + width - margin.right) / 2).attr("y", dHeight - 8).attr("text-anchor", "middle").attr("font-size", 13).attr("fill", "currentColor")
  const dRows = dSvg.selectAll("g.sb-row").data(scales).join("g").attr("class", "sb-row").attr("transform", (s, i) => `translate(0, ${18 + dRowHeight * i + 22})`)
  dRows.append("text").attr("text-anchor", "end").attr("font-size", 13).attr("font-weight", 600).attr("fill", "currentColor").html(rowName)
  const dBar = dRows.append("line").attr("stroke", grey).attr("stroke-width", 3)
  const dDot = dRows.append("rect").attr("y", -6).attr("width", 12).attr("height", 12).attr("fill", grey)
  const dValue = dRows.append("text").attr("y", -12).attr("text-anchor", "middle").attr("font-size", 12).attr("font-weight", 600).attr("fill", "currentColor")
  const dZero = dSvg.append("line").attr("y1", 10).attr("y2", dHeight - dBottom).attr("stroke", "currentColor").attr("stroke-dasharray", "4,3").attr("stroke-opacity", 0.6)
  const dTitle = html`<div class="sb-subtitle"></div>`

  const buttons = views.map(v => Object.assign(html`<button class="btn btn-quarto">${v.button}</button>`, {onclick: () => show(v, v, true)}))
  const caption = html`<div class="sb-caption"></div>`
  // view is the units of the main panel, and diffView the units of the difference panel
  const show = (view, diffView, animate) => {
    views.forEach((v, i) => buttons[i].classList.toggle("sb-off", v !== view))
    caption.textContent = fixedViews ? captions[0] : captions[views.indexOf(view)]
    axisLabel.text(view.axis)
    dAxisLabel.text(diffView.axis)
    const x = d3.scaleLinear()
      .domain([Math.min(0, d3.min(scales, s => view.f(s, s.min))), d3.max(scales, s => view.f(s, s.max))]).nice()
      .range([margin.left, width - margin.right])
    const at = (s, v) => x(view.f(s, v))
    // scale points are shown without trailing zeros (28, 12.5)
    const fmt = v => view.digits === 0 ? String(+v.toFixed(1)) : v.toFixed(view.digits)
    const highlight = (text, on) => on ? `<tspan font-weight="700" fill="${blue}">${text}</tspan>` : text
    labels.html(s => [
      `range ${+fmt(view.f(s, s.min))} to ${+fmt(view.f(s, s.max))}`,
      `M = ${s.means.map((mean, j) => `<tspan font-weight="700" fill="${groups[j].colour}">${fmt(view.f(s, mean))}</tspan>`).join(" vs ")}`,
      highlight(`SD = ${fmt(view.f(s, s.min + s.sd) - view.f(s, s.min))}`, view.aligned === "sd")
    ].join(" · "))
    // every curve has the same area, so a distribution that is narrower on the page is also taller. a curve's
    // height in pixels is its density divided by the width of one scale point in pixels, and the tallest curve in
    // the current view is given the full height
    const pointWidth = s => Math.abs(at(s, s.min + 1) - at(s, s.min))
    const score = (s, u) => s.min + u * (s.max - s.min)
    const tallest = d3.max(marksData, d => d3.max(us, u => d.density(score(d.s, u))) / pointWidth(d.s))
    const curveY = (d, x) => -curveHeight * d.density(x) / pointWidth(d.s) / tallest
    const points = d => us.map(u => `${at(d.s, score(d.s, u)).toFixed(2)},${curveY(d, score(d.s, u)).toFixed(2)}`)
    // the elements slide to their new positions, unless the reader has asked for reduced motion
    // otherwise they are placed directly: a transition, even one of no duration, only starts on the next frame,
    // which would leave a newly drawn plot unpositioned for a moment
    const slide = animate && !matchMedia("(prefers-reduced-motion: reduce)").matches
    const move = selection => slide ? selection.transition().duration(800) : selection
    move(axisG).call(d3.axisBottom(x).ticks(6).tickSizeInner(-(height - margin.top - margin.bottom)).tickPadding(8))
    move(range).attr("x1", s => at(s, s.min)).attr("x2", s => at(s, s.max))
    move(capLo).attr("x1", s => at(s, s.min)).attr("x2", s => at(s, s.min))
    move(capHi).attr("x1", s => at(s, s.max)).attr("x2", s => at(s, s.max))
    move(outline).attr("d", d => `M${points(d).join("L")}`)
    move(fill).attr("d", d => `M${points(d).join("L")}Z`)
    move(meanLine).attr("x1", d => at(d.s, d.mean)).attr("x2", d => at(d.s, d.mean)).attr("y2", d => curveY(d, d.mean))
    // difference between the two groups' means, in the difference panel's units
    const diff = s => diffView.f(s, s.means[1]) - diffView.f(s, s.means[0])
    const cis = new Map(scales.map(s => [s, diffView.ci(s)]))
    const dx = d3.scaleLinear().domain(d3.extent([0, ...scales.flatMap(s => cis.get(s))])).nice().range([margin.left, width - margin.right])
    dTitle.textContent = `${diffView.difference} (n = ${n} per condition)`
    move(dAxisG).call(d3.axisBottom(dx).ticks(6).tickSizeInner(-(dHeight - dBottom - 10)).tickPadding(8))
    move(dZero).attr("x1", dx(0)).attr("x2", dx(0))
    const number = v => v.toFixed(Math.max(1, diffView.digits)).replace("-", "−")
    move(dBar).attr("x1", s => dx(cis.get(s)[0])).attr("x2", s => dx(cis.get(s)[1]))
    move(dDot).attr("x", s => dx(diff(s)) - 6)
    move(dValue).attr("x", s => dx(diff(s))).text(s => `${number(diff(s))} [${cis.get(s).map(number).join(", ")}]`)
  }
  if (fixedViews) show(views[fixedViews[0]], views[fixedViews[1]], false); else show(views[0], views[0], false)
  return html`<div>
    ${fixedViews ? "" : html`<div class="sb-buttons">${buttons}</div>`}
    <div class="sb-legend">${groups.map(g => html`<span><i style="background: ${g.colour};"></i>${g.name}</span>`)}</div>
    ${svg.node()}
    <div>${dTitle}${dSvg.node()}</div>
    ${caption}
  </div>`
}

function boundedDensity(mean, sd, lo, hi) {
  // the density of a beta distribution rescaled to run from lo to hi, with the given mean and SD. its parameters
  // come from the method of moments in the mean-precision form used in beta regression: with the mean and
  // variance rescaled to the 0 to 1 interval, precision = mean * (1 - mean) / variance - 1,
  // alpha = mean * precision, and beta = (1 - mean) * precision. unlike a normal distribution, it cannot run
  // past the ends of the scale. returns a function of the score
  const range = hi - lo, m = (mean - lo) / range, v = (sd / range) ** 2
  const precision = m * (1 - m) / v - 1
  const alpha = m * precision, beta = (1 - m) * precision
  const constant = logGamma(alpha + beta) - logGamma(alpha) - logGamma(beta)
  return x => {
    const u = (x - lo) / range
    if (u <= 0 || u >= 1) return 0
    return Math.exp((alpha - 1) * Math.log(u) + (beta - 1) * Math.log(1 - u) + constant) / range
  }
}




sampleBuilderStyle = html`<style>
  .sb-subtitle { font-size: 0.9rem; font-weight: 600; margin: 0.75rem 0 0.25rem; color: var(--bs-secondary-color); }
  .sb-buttons { display: flex; gap: 0.5rem; flex-wrap: wrap; }
  .sb-tiles { display: flex; gap: 2rem; flex-wrap: wrap; margin: 0.5rem 0 1.5rem; }
  .sb-tiles > div { flex: 1 1 12rem; }
  .sb-tile-value { font-size: 2rem; font-weight: 700; line-height: 1.1; color: #2780e3; }
  .sb-tile-label { font-weight: 600; }
  .sb-tile-note { font-size: 0.9rem; color: var(--bs-secondary-color); }
  /* standardisation plot: the buttons other than the current view are outlined */
  .btn.btn-quarto.sb-off { background: transparent; color: #2780e3; }
  .btn.btn-quarto.sb-off:hover { background: #1f6ab8; color: #fff; }
  .sb-std { max-width: 100%; overflow: visible; }
  .sb-std .domain { display: none; }
  .sb-std .tick line { stroke-opacity: 0.12; }
  .sb-std .tick text { font-size: 12px; }
  .sb-legend { display: flex; gap: 1.25rem; font-size: 0.9rem; margin-top: 0.25rem; }
  .sb-legend i { display: inline-block; width: 10px; height: 10px; margin-right: 0.4rem; }
  .sb-caption { min-height: 4.5em; margin: 0.25rem 0 1rem; font-weight: 600; }
</style>`

// quantiles of the t distribution, found by bisection on its cumulative distribution function
function tQuantile(p, df) {
  const cdf = x => x > 0 ? 1 - 0.5 * betaI(df / (df + x * x), df / 2, 0.5) : 0.5 * betaI(df / (df + x * x), df / 2, 0.5)
  return bisect(cdf, p, -1e4, 1e4)
}
function bisect(cdf, p, lower, upper) {
  for (let i = 0; i < 100; i++) {
    const mid = (lower + upper) / 2
    if (cdf(mid) < p) lower = mid; else upper = mid
  }
  return (lower + upper) / 2
}
function nctCdf(t, df, ncp) {
  // cumulative distribution function of the noncentral t distribution. T = (Z + ncp) / (U / sqrt(df)), where Z is
  // standard normal and U follows a chi distribution, so P(T <= t) is the average of normCdf(t * u / sqrt(df) - ncp)
  // over the distribution of U. the average is calculated by numerical integration (Simpson's rule)
  const logConstant = -(df / 2 - 1) * Math.LN2 - logGamma(df / 2)
  const density = u => u > 0 ? Math.exp((df - 1) * Math.log(u) - u * u / 2 + logConstant) : (df === 1 ? Math.exp(logConstant) : 0)
  const from = Math.max(0, Math.sqrt(df) - 8), to = Math.sqrt(df) + 8, steps = 400, h = (to - from) / steps
  let sum = 0
  for (let i = 0; i <= steps; i++) {
    const u = from + i * h
    sum += (i === 0 || i === steps ? 1 : i % 2 ? 4 : 2) * density(u) * normCdf(t * u / Math.sqrt(df) - ncp)
  }
  return sum * h / 3
}
function normCdf(x) {
  // Abramowitz & Stegun 7.1.26
  const t = 1 / (1 + 0.3275911 * Math.abs(x) / Math.SQRT2)
  const erf = 1 - (((((1.061405429 * t - 1.453152027) * t) + 1.421413741) * t - 0.284496736) * t + 0.254829592) * t * Math.exp(-x * x / 2)
  return x >= 0 ? 0.5 * (1 + erf) : 0.5 * (1 - erf)
}
function logGamma(x) {
  // Lanczos approximation
  const c = [76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 0.001208650973866179, -0.000005395239384953]
  let y = x, sum = 1.000000000190015
  for (const cj of c) sum += cj / ++y
  return -(x + 5.5) + (x + 0.5) * Math.log(x + 5.5) + Math.log(2.5066282746310005 * sum / x)
}
function betaI(x, a, b) {
  // regularised incomplete beta function, by continued fraction
  if (x <= 0) return 0
  if (x >= 1) return 1
  const front = Math.exp(logGamma(a + b) - logGamma(a) - logGamma(b) + a * Math.log(x) + b * Math.log(1 - x))
  const cf = (x, a, b) => {
    let c = 1, d = 1 - (a + b) * x / (a + 1)
    if (Math.abs(d) < 1e-300) d = 1e-300
    d = 1 / d
    let h = d
    for (let m = 1; m < 500; m++) {
      let aa = m * (b - m) * x / ((a + 2 * m - 1) * (a + 2 * m))
      d = 1 + aa * d; if (Math.abs(d) < 1e-300) d = 1e-300
      c = 1 + aa / c; if (Math.abs(c) < 1e-300) c = 1e-300
      d = 1 / d
      h *= d * c
      aa = -(a + m) * (a + b + m) * x / ((a + 2 * m) * (a + 2 * m + 1))
      d = 1 + aa * d; if (Math.abs(d) < 1e-300) d = 1e-300
      c = 1 + aa / c; if (Math.abs(c) < 1e-300) c = 1e-300
      d = 1 / d
      const del = d * c
      h *= del
      if (Math.abs(del - 1) < 1e-14) break
    }
    return h
  }
  return x < (a + 1) / (a + b + 2) ? front * cf(x, a, b) / a : 1 - front * cf(1 - x, b, a) / b
}

5.3 When standarization is useful

The same construct is often measured with different scales, and different scales have different units. Depression, for example, is commonly measured with the BDI-II, which ranges from 0 to 63, the PHQ-9, which ranges from 0 to 27, and the HAM-D, which ranges from 0 to 52. A difference of 5 points between two groups means something different on each of them, so results from studies that used different scales cannot be compared or combined in their original units.

The plot below shows these three scales. Each row shows a scale’s logical range, with the distributions of a control condition and an intervention condition from a hypothetical trial of a treatment for depression drawn on it with the same SD. The distributions are beta distributions, rescaled to the range of each scale, with parameters calculated from the mean and SD by the method of moments; unlike normal distributions, they cannot run past the ends of the scale. The values are illustrative, not taken from particular studies. Lower scores mean less severe depression, so a negative difference favours the intervention. The panel underneath shows the difference between the two means on each scale, in whichever units are selected, with its 95% confidence interval. Each trial is assumed to have 60 participants per condition. The buttons re-express the same three trials in different units:

• Scale units are the original units of each scale.
• SD units divide every score by that sample’s standard deviation. This is the standardisation used by Cohen’s d and other standardised mean differences.

Switch between them and pay attention to what lines up and what does not.

standardisedDifferencePlot = unitsPlot({
  // illustrative values for a hypothetical trial: means are [control, intervention], with a common SD
  scales: [
    {name: "BDI-II", min: 0, max: 63, sd: 10, means: [28, 22]},
    {name: "PHQ-9", min: 0, max: 27, sd: 5, means: [15, 12.5]},
    {name: "HAM-D", min: 0, max: 52, sd: 6, means: [20, 16.7]}
  ],
  groups: [{name: "Control", colour: "#d9622b"}, {name: "Intervention", colour: "#2780e3"}],
  captions: [
    "In scale units, the three differences cannot be compared, because a point on one scale is not a point on another.",
    "In SD units, each difference is a Cohen's d, and the three can be compared. They are similar but not identical. The means themselves still do not line up, and neither do the ranges."
  ],
  n: 60,
  width
})

In scale points, the three trials appear to have found very different effects: 6 points, 2.5 points, and 3.3 points. In SD units, they are d = −0.60, −0.50, and −0.55, which is what we might expect from three trials of the same treatment. Standardising has not made the three results identical, and it should not: it has only put them in the same units so that they can be compared.

Standardising by the SD makes one thing comparable across scales: the unit. One SD is the same length on every scale, so in SD units the distributions on all three scales have the same width, which is what allows a mean difference of 6 points on the BDI-II to be compared with a mean difference of 2.5 points on the PHQ-9. It does not make anything else comparable. The three trials still have different means in SD units, and the scales still have different ranges. Note also that the unit is a property of the sample, not of the scale: a study that recruited a more homogeneous sample would have a smaller SD, and the same difference in scale points would become a larger standardised difference (see Section 5.4).

5.4 When standardization obscures: Range restriction

Standardisation solves one problem and creates another. Because Cohen’s d divides the mean difference by the SD of the sample, anything that changes the SD changes d, even when the effect of the intervention in scale points is exactly the same. The SD of a sample depends on who was recruited. A trial that accepts anyone with a diagnosis of depression will include a wide range of severity. A trial with narrow inclusion criteria, for example one that only accepts patients whose baseline scores fall within a particular band, will include a much narrower range. This is called range restriction, and it shrinks the SD.

The plot below shows two hypothetical trials that used the same scale, the BDI-II, and found exactly the same difference between the intervention and control conditions: 5 points. Study A recruited a sample with a wide range of severity, and Study B recruited a sample with a restricted range. Use the slider to change how restricted Study B’s sample is. The top panel shows the two studies’ distributions in scale units, and the panel underneath shows each study’s standardised mean difference, Cohen’s d, with its 95% confidence interval. As in the previous plot, each trial is assumed to have 60 participants per condition.

viewof rrSd = Inputs.range([3, 10], {value: 5, step: 0.5, label: "SD in Study B"})

rrPlot = unitsPlot({
  // two hypothetical trials on the BDI-II with the same means, [control, intervention], and so the same mean
  // difference, but different SDs
  scales: [
    {name: "Study A", sub: "BDI-II", min: 0, max: 63, sd: 10, means: [28, 23]},
    {name: "Study B", sub: "BDI-II", min: 0, max: 63, sd: rrSd, means: [28, 23]}
  ],
  scaleUnits: "Scale units: BDI-II score, 0 to 63, in both studies",
  groups: [{name: "Control", colour: "#d9622b"}, {name: "Intervention", colour: "#2780e3"}],
  captions: [
    "In scale units, the two studies found exactly the same difference: 5 points on the BDI-II. But the same 5-point difference is a different Cohen's d in each study, because each study's difference is divided by its own SD."
  ],
  n: 60,
  // the distributions are always shown in scale units, and the difference in SD units (Cohen's d)
  fixedViews: [0, 1],
  width
})

The intervention did exactly the same thing in both studies: it lowered scores by 5 points. In the default setting, Study B’s d is nevertheless twice the size of Study A’s, and a reader comparing the two standardised effect sizes would conclude that the intervention worked twice as well in Study B. The difference between the two values of d comes entirely from who was recruited.

This has several implications for trustworthiness assessment:

• A standardised effect size cannot be interpreted without knowing the SD it was standardised by. When an effect size seems implausibly large, check whether the SDs are unusually small for that scale, compared with other studies of similar populations.
• Standardised effect sizes from studies with different inclusion criteria are not directly comparable, even when the studies used the same scale. When the studies did use the same scale, the mean difference in scale points is often the more informative comparison.
• Range restriction is only one of the things that change the SD. Others include the reliability of the measure, floor and ceiling effects, and which SD was used as the standardiser (for example, the baseline SD, the pooled SD at post-treatment, or the SD of change scores).