8  Effect size methods: SMD/Cohen’s d ✎ Draft

8.1 Influence of scale bounds on SMDs

The plot above assumes normal distributions, which have no minimum or maximum. Real scales do, and this limits which effect sizes are possible. On a bounded scale, the closer a mean is to the end of the scale, the smaller its SD has to be. A large Cohen’s d requires two means that are many SDs apart, but on a bounded scale the means can only move so far apart before one of them is too close to the end of the scale to have the SD that was reported.

The plot below shows this. Set the logical minimum and maximum of the scale, the SD of both conditions in scale units, and Cohen’s d. The two conditions start at the middle of the scale and move apart as d increases, so that their means are always d × SD scale points apart. Each of the two buttons switches between two settings, and is labelled with the setting that clicking it will change to. The first chooses how the conditions move apart: at first the control condition stays in the middle of the scale and only the intervention condition moves, and the alternative is that both conditions move away from the middle by the same amount. The second chooses the direction of the effect: at first the intervention increases scores, so that Cohen’s d is positive, and the alternative is that it decreases them, so that d is negative. Each condition is drawn as a beta distribution, rescaled to the range of the scale, whose parameters are calculated from its mean and SD by the method of moments. The message under the plot says whether the combination you have chosen is possible.

viewof bBothMove = toggleButton([{label: "Only intervention condition moves", value: false}, {label: "Both conditions move", value: true}])
viewof bIncrease = toggleButton([{label: "Intervention increases scores", value: true}, {label: "Intervention decreases scores", value: false}])

viewof bMin = Inputs.number({value: 1, step: 1, label: "Scale minimum"})
viewof bMax = Inputs.number({value: 7, step: 1, label: "Scale maximum"})

bBounds = Number.isFinite(bMin) && Number.isFinite(bMax) && bMin < bMax ? [bMin, bMax] : null
viewof bSd = Inputs.range([0.1, ((bBounds ?? [1, 7])[1] - (bBounds ?? [1, 7])[0]) / 2], {value: Math.round(((bBounds ?? [1, 7])[1] - (bBounds ?? [1, 7])[0]) / 5 * 10) / 10, step: 0.05, label: "SD (scale units)"})
viewof bD = Inputs.range([0, 6], {value: 0.5, step: 0.05, label: "Size of Cohen's d"})

bState = {
  if (!bBounds) return null
  const [lo, hi] = bBounds, range = hi - lo, mid = (lo + hi) / 2
  // the two conditions for a given size of d. the two means start at the middle of the scale and stay d SDs from
  // each other. either both move away from the middle by the same amount, or the control condition stays in the
  // middle and only the intervention moves. the intervention moves up the scale if it increases scores (a
  // positive d) and down the scale if it decreases them (a negative d)
  const groupsFor = d => {
    const shift = (bIncrease ? 1 : -1) * d * bSd
    return [{group: "Control", mean: bBothMove ? mid - shift / 2 : mid}, {group: "Intervention", mean: bBothMove ? mid + shift / 2 : mid + shift}].map(g => {
      // method of moments for a beta distribution in its mean-precision form, as used in beta regression:
      // rescale the mean and variance to the 0 to 1 interval, then precision = mean * (1 - mean) / variance - 1,
      // alpha = mean * precision, and beta = (1 - mean) * precision. a valid distribution needs precision > 0
      const m = (g.mean - lo) / range, v = (bSd / range) ** 2
      const precision = m * (1 - m) / v - 1
      const alpha = m * precision, beta = (1 - m) * precision
      const valid = m > 0 && m < 1 && precision > 0
      // shape: with alpha and beta above 1 the distribution is hump-shaped. with one of them at or below 1,
      // scores pile up against one end of the scale (J-shaped). with both below 1, they pile up against both
      // ends (U-shaped)
      const shape = !valid ? null : alpha > 1 && beta > 1 ? "hump" : alpha < 1 && beta < 1 ? "U" : "J"
      return {...g, m, alpha, beta, shape}
    })
  }
  // a combination is treated as possible only if neither condition would have to be U-shaped. a U-shaped
  // distribution is not strictly impossible, but it is not a plausible description of data from a scale, and the
  // plot does not rely on the beta distribution that far
  const possibleFor = d => groupsFor(d).every(g => g.shape === "hump" || g.shape === "J")
  // the largest possible d, by bisection: the conditions only become less plausible as they move apart
  let dMax = 0
  if (possibleFor(0)) {
    let lower = 0, upper = 2 * range / bSd
    for (let i = 0; i < 60; i++) {
      const middle = (lower + upper) / 2
      if (possibleFor(middle)) lower = middle; else upper = middle
    }
    dMax = lower
  }
  const possible = possibleFor(bD)
  // share of each condition's participants with the lowest and the highest possible score: on a scale with
  // whole-number response options, these are the scores within half a point of each end of the scale
  const edge = (range >= 2 ? 0.5 : range / 20) / range
  const groups = groupsFor(bD).map(g => possible ? {...g, bottom: betaI(edge, g.alpha, g.beta), top: 1 - betaI(1 - edge, g.alpha, g.beta)} : g)
  return {lo, hi, range, groups, dMax, possible, piled: possible && groups.some(g => g.shape === "J")}
}

{
  if (!bState) return html`<div style="color: var(--bs-secondary-color);">Set the minimum and maximum to numbers, with the minimum below the maximum.</div>`
  const {lo, hi, range, groups, possible} = bState
  const betaDensity = (x, g) => {
    const u = (x - lo) / range
    return Math.exp((g.alpha - 1) * Math.log(u) + (g.beta - 1) * Math.log(1 - u) + logGamma(g.alpha + g.beta) - logGamma(g.alpha) - logGamma(g.beta)) / range
  }
  // the ends of the scale are left out, because the density can be infinite there
  const xs = d3.range(1, 400).map(i => lo + range * i / 400)
  const curves = possible ? groups.flatMap(g => xs.map(x => ({group: g.group, x, y: betaDensity(x, g)}))) : []
  // the vertical axis is capped, so that a spike against the end of the scale does not flatten everything else
  const yMax = Math.min(1.3 * (d3.max(curves, d => d.y) ?? 1), 1.6 / bSd)
  // a pile of scores against the end of the scale is a spike that runs off the top of the plot, so it is easy
  // to miss how much of the distribution is in it. label each pile with the share of participants in it
  const pct = p => p > 0.99 ? "more than 99%" : p < 0.01 ? "less than 1%" : `${Math.round(100 * p)}%`
  const piles = !possible ? [] : groups.flatMap(g => [
    g.alpha <= 1 ? {group: g.group, end: "bottom", x: lo, share: g.bottom} : null,
    g.beta <= 1 ? {group: g.group, end: "top", x: hi, share: g.top} : null
  ]).filter(d => d && d.share >= 0.01)
  const pileLabel = (end, anchor, dx, dy) => Plot.text(piles.filter(d => d.end === end), {
    x: "x", y: yMax, text: d => `${pct(d.share)} of ${d.group.toLowerCase()}\nparticipants have ${d.end === "top" ? "highest" : "lowest"}\npossible score`, fill: "group", fontWeight: 600,
    textAnchor: anchor, lineAnchor: "top", dx, dy, lineHeight: 1.15
  })
  return Plot.plot({
    width, height: Math.min(360, Math.max(260, width * 0.5)), marginLeft: 20, marginRight: 20, marginTop: 10, marginBottom: 45,
    x: {label: "Scale units", labelAnchor: "center", labelOffset: 38, domain: [lo, hi], ticks: range <= 20 ? d3.range(Math.ceil(lo), Math.floor(hi) + 1) : undefined, tickFormat: range <= 20 ? "d" : undefined},
    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, clip: true}),
      Plot.lineY(curves, {x: "x", y: "y", stroke: "group", strokeWidth: 2, clip: true}),
      Plot.ruleY([0], {stroke: "currentColor"}),
      // each condition's mean, where it falls on the scale
      Plot.ruleX(groups.filter(g => g.mean >= lo && g.mean <= hi), {x: "mean", stroke: "group", strokeDasharray: "4,3"}),
      // when there is a pile at each end, the second label sits below the first, so that they cannot overlap
      pileLabel("bottom", "start", 6, 4),
      pileLabel("top", "end", -6, piles.some(d => d.end === "bottom") ? 54 : 4),
      possible ? null : Plot.text(["Not possible on this scale"], {frameAnchor: "middle", fontSize: 16, fontWeight: 600, fillOpacity: 0.6})
    ]
  })
}

{
  if (!bState) return html``
  const {lo, hi, range, groups, dMax, possible, piled} = bState
  const means = `The means would be ${groups[0].mean.toFixed(2)} (control) and ${groups[1].mean.toFixed(2)} (intervention), ${(bD * bSd).toFixed(2)} scale points apart.`
  const centred = bBothMove ? "" : " while the control condition stays in the middle of the scale"
  // round the largest possible d down, so that the value shown is itself possible
  const signed = x => (bIncrease || x === 0 ? "" : "−") + x.toFixed(2)
  const limit = `${bIncrease ? "at most" : "no lower than"} ${signed(Math.floor(dMax * 100) / 100)}${centred}`
  const pct = p => p > 0.99 ? "more than 99%" : p < 0.01 ? "less than 1%" : `about ${Math.round(100 * p)}%`
  // what each condition whose scores pile up against an end of the scale would have to look like
  const piles = groups.filter(g => g.shape === "J").map(g => g.beta <= 1
    ? `${pct(g.top)} of ${g.group.toLowerCase()} participants would have to have the highest possible score`
    : `${pct(g.bottom)} of ${g.group.toLowerCase()} participants would have to have the lowest possible score`).join(", and ")
  const sentence = s => s.charAt(0).toUpperCase() + s.slice(1) + "."
  const [kind, title, detail] =
    possible && !piled ? ["ok", "Possible.", `${means} With an SD of ${bSd.toFixed(2)} on a ${lo} to ${hi} scale, Cohen's d can be ${limit}.`]
    : possible ? ["warn", "Possible, but only with scores piled up against the end of the scale.", `${means} ${sentence(piles)} With an SD of ${bSd.toFixed(2)} on a ${lo} to ${hi} scale, Cohen's d can be ${limit}.`]
    : dMax === 0 ? ["bad", `Not possible: an SD of ${bSd.toFixed(2)} is too large for a ${lo} to ${hi} scale.`, `Even with no difference between the conditions, scores would have to be split between the two ends of the scale. Scores spread evenly across the whole scale have an SD of ${(range / Math.sqrt(12)).toFixed(2)}, and any larger SD requires this.`]
    : ["bad", `Not possible: Cohen's d = ${signed(bD)} cannot occur with an SD of ${bSd.toFixed(2)} on a ${lo} to ${hi} scale${centred}.`, `${means} ${groups.some(g => g.mean <= lo || g.mean >= hi) ? "A mean cannot be at or beyond the end of the scale." : "A mean that close to the end of the scale can only have an SD that large if scores are split between the two ends of the scale."} With this SD, Cohen's d can be ${limit}.`]
  return html`<div class="sb-feedback sb-feedback-${kind}"><strong>${title}</strong> ${detail}</div>`
}

function toggleButton(options) {
  // a button that switches between two options ({label, value}) each time it is clicked. it is labelled with what
  // clicking it will change to, which is the option that is not currently in effect. the value is stored on the
  // wrapping div, so that it can be used with viewof. starts on the first option
  let current = 0
  const label = () => { const next = options[(current + 1) % options.length].label; return `Change to ${next.charAt(0).toLowerCase()}${next.slice(1)}` }
  const button = html`<button class="btn btn-quarto">${label()}</button>`
  const view = html`<div class="sb-buttons">${button}</div>`
  view.value = options[0].value
  button.onclick = () => {
    current = (current + 1) % options.length
    button.textContent = label()
    view.value = options[current].value
    view.dispatchEvent(new Event("input", {bubbles: true}))
  }
  return view
}

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
}

boundedPlotStyle = html`<style>
  .sb-buttons { display: flex; gap: 0.5rem; flex-wrap: wrap; }
  .sb-feedback { margin: 0.5rem 0 1.5rem; padding: 0.6rem 0.9rem; border-left: 4px solid; background: rgba(128, 128, 128, 0.08); }
  .sb-feedback-ok { border-color: #3b9c6e; }
  .sb-feedback-warn { border-color: #d9a02b; }
  .sb-feedback-bad { border-color: #c0392b; }
</style>`

The message under the plot distinguishes three situations. The combination is possible with ordinary hump-shaped distributions. Or it is possible only if at least one condition’s scores are piled up against one end of the scale, as happens with floor and ceiling effects; because such a pile is a spike that runs off the top of the plot, it is labelled with the share of that condition’s participants who would have the highest or lowest possible score. Or it is not possible.

The plot treats a combination as not possible once either condition’s scores would have to be split between the two ends of the scale, with most at one end and a few at the other. Strictly, such data could exist. The absolute limit, which applies to any distribution and not only to beta distributions, is that a sample with a given mean on a scale from min to max cannot have a variance greater than (mean − min) × (max − mean), and this is only reached when every score is at one end of the scale or the other. The limit used by the plot is somewhat stricter than this, because data split between the two ends of a scale are not a plausible description of responses to that scale. One consequence is that the plot treats any SD larger than that of scores spread evenly across the whole scale, which is the range divided by √12, as not possible at all.

Try increasing the SD: the larger the SD relative to the range of the scale, the smaller the largest possible d. The largest possible d is also only half as large when the control condition’s mean stays in the middle of the scale, because only the intervention condition’s mean can then move towards the end of the scale. When an article reports a large standardised effect on a bounded scale, it is worth checking whether the reported means and SDs are compatible with the scale at all.

8.2 Statistical distribution of SMD effect sizes

Whether a reported effect size is plausible depends on more than the size of the true effect. The effect size observed in any one study also depends on sampling error, which is larger in small studies, and on how much the true effect varies between studies. The dance of the p-values (Section 4.4.0.1) showed how much the effect size and p-value of a study vary from one sample to the next, even when nothing about the true effect changes. Taken together, these determine the range within which the observed effect size of a single study can be expected to fall.

The plot below shows this range for Cohen’s d from a two-group design, as a function of the total sample size. For each sample size, it simulates 40,000 studies. Each study’s true standardised mean difference is drawn from a normal distribution with mean μ and standard deviation τ (the between-study heterogeneity). The study’s observed d is then drawn from its sampling distribution, with the total sample split equally between the two groups, with the usual correction for small-sample bias applied. The dashed line shows the mean true effect μ, and the bands show the 95%, 99%, and 99.9% prediction intervals: the ranges within which 95%, 99%, and 99.9% of studies’ observed effect sizes fall.

viewof piMu = Inputs.range([0, 2.5], {value: 0.5, step: 0.1, label: "True (population) SMD"})
viewof piTau = Inputs.range([0, 1], {value: 0.05, step: 0.05, label: "Between-study heterogeneity τ"})

piData = simulatePiD(piMu, piTau)

Plot.plot({
  width, height: Math.min(480, Math.max(320, width * 0.75)), marginLeft: 55, marginRight: 20, marginBottom: 45,
  x: {label: "Total N (both arms)", labelAnchor: "center", labelOffset: 38, ticks: [20, 50, 100, 250, 500]},
  y: {label: "Observed Cohen's d", labelAnchor: "center", grid: true,
      domain: [Math.min(-0.5, d3.min(piData, d => d.lo999)), Math.max(3.2, d3.max(piData, d => d.hi999))]},
  style: {fontSize: "13px"},
  marks: [
    Plot.frame({stroke: "currentColor"}),
    Plot.ruleY([0], {stroke: "currentColor", strokeOpacity: 0.5}),
    Plot.areaY(piData, {x: "n", y1: "lo999", y2: "hi999", fill: "currentColor", fillOpacity: 0.1}),
    Plot.areaY(piData, {x: "n", y1: "lo99", y2: "hi99", fill: "currentColor", fillOpacity: 0.15}),
    Plot.areaY(piData, {x: "n", y1: "lo95", y2: "hi95", fill: "currentColor", fillOpacity: 0.2}),
    // the mean true SMD, drawn over the bands
    Plot.ruleY([piMu], {stroke: "#2780e3", strokeDasharray: "6,4", strokeWidth: 2.5})
  ]
})

95%, 99%, and 99.9% prediction intervals (darkest to lightest bands) of an observed d.

function simulatePiD(mu, tau, draws = 40000) {
  // quantiles of the observed, bias-corrected Cohen's d of a two-group study, at each of a log-spaced grid of total sample sizes.
  // each simulated study's true SMD, delta, is drawn from N(mu, tau). its observed d is
  // t * 2 / sqrt(N), where t follows a noncentral t distribution with df = N - 2 and ncp = delta * sqrt(N) / 2,
  // multiplied by the small-sample bias correction J = 1 - 3 / (4 * df - 1). a noncentral t is drawn as (Z + ncp) / sqrt(chi-squared / df).
  // the random numbers are seeded and reused at every sample size, so that the result does not change
  // between redraws and the curves are smooth
  const nSeq = [...new Set(d3.range(16).map(i => Math.round(Math.exp(Math.log(20) + i * (Math.log(500) - Math.log(20)) / 15))))]
  return nSeq.map(N => {
    // mulberry32 pseudo-random number generator
    let state = 1
    const uniform = () => {
      state = (state + 0x6D2B79F5) | 0
      let x = Math.imul(state ^ (state >>> 15), 1 | state)
      x = (x + Math.imul(x ^ (x >>> 7), 61 | x)) ^ x
      return ((x ^ (x >>> 14)) >>> 0) / 4294967296
    }
    // standard normal, by the Box-Muller method
    const normal = () => Math.sqrt(-2 * Math.log(1 - uniform())) * Math.cos(2 * Math.PI * uniform())
    // gamma(shape, 1) for shape >= 1, by the Marsaglia-Tsang method. chi-squared(df) is 2 * gamma(df / 2)
    const gamma = shape => {
      const a = shape - 1 / 3, c = 1 / Math.sqrt(9 * a)
      for (;;) {
        const x = normal(), v = (1 + c * x) ** 3
        if (v <= 0) continue
        const u = 1 - uniform()
        if (u < 1 - 0.0331 * x ** 4 || Math.log(u) < 0.5 * x * x + a * (1 - v + Math.log(v))) return a * v
      }
    }
    const df = N - 2, J = 1 - 3 / (4 * df - 1)
    const d = new Float64Array(draws)
    for (let i = 0; i < draws; i++) {
      const delta = mu + tau * normal()
      const t = (normal() + delta * Math.sqrt(N) / 2) / Math.sqrt(2 * gamma(df / 2) / df)
      d[i] = J * t * 2 / Math.sqrt(N)
    }
    d.sort()
    const q = p => d3.quantileSorted(d, p)
    return {n: N, lo999: q(0.0005), lo99: q(0.005), lo95: q(0.025), hi95: q(0.975), hi99: q(0.995), hi999: q(0.9995)}
  })
}

Several things are worth trying:

• At small sample sizes, the intervals are wide: a study with a total N of 20 can easily observe an effect twice the size of the true effect, or an effect in the opposite direction. Large observed effects in small studies are therefore not surprising on their own.
• The intervals narrow quickly as the sample size increases. An observed d of 2 or more in a study with several hundred participants falls far outside the 99% interval unless the true effect is itself very large. Increase μ to find out how large the true effect would have to be before such a result became plausible, and consider whether an effect of that size is credible for the phenomenon in question.
• Increasing the heterogeneity τ widens the intervals at every sample size, including large ones, because it represents real differences between studies’ true effects and not sampling error.

An observed effect size that falls outside the 99% interval for its sample size, under generous assumptions about the true effect and its heterogeneity, is unlikely to have arisen from a correctly calculated result. It is then worth checking for the common errors described below.

8.3 Empirical distributions of SMD effect sizes

TODO

• Distributions of d, r, partial eta-squared, etc.
  • For different types of phenomena/tests (e.g., interventions vs manipulation checks)
  • For different fields
  • Avoiding ever saying “domain knowledge”

8.4 Cohen’s d explorer

The app below lets you build two samples on the same scale: an intervention condition and a control condition. Set the logical minimum and maximum of the scale, and then set the number of observations at each response option in either of two ways. First, you can draw each distribution: click anywhere in a column to set its bar to that height, or hold the mouse button down and drag across the plot to draw the whole histogram in one movement. The vertical axis rescales when you release, so dragging a bar to the top of the plot makes room for larger counts. Second, you can type the exact number of observations for each response option into the boxes under the axis. The scale can have up to 100 points; with more than 25, the boxes are not shown, because the columns are too narrow for them. The estimates below the two conditions update as you go, each with its 95% confidence interval. They are plotted in three groups, each with its own axis: the two means and the difference between them (intervention minus control), the two SDs and the pooled SD, and the standardised mean difference, Cohen’s d. The first two groups are in the units of the scale, whereas d is in standard deviations. The interval for the mean difference uses the t distribution, the interval for the pooled SD uses the chi-squared distribution, and the interval for d uses the noncentral t distribution, as in the worked example in Section 7.1. Try changing only the spread of the two samples while keeping their means the same, to see how much d depends on the standard deviation, and try to find the largest d that the scale allows.

// the plot below attaches the button's action
dButtons = html`<div class="sb-buttons">
  <button class="btn btn-quarto sb-clear">Clear</button>
</div>`

viewof dMin = Inputs.number({value: 1, step: 1, label: "Scale minimum"})
viewof dMax = Inputs.number({value: 7, step: 1, label: "Scale maximum"})

dPlot = {
  const scores = scaleScores(dMin, dMax)
  if (!scores) return scaleMessage()
  const groups = [
    {label: "Intervention condition", short: "Intervention", store: dStores.intervention, colour: "#2780e3"},
    {label: "Control condition", short: "Control", store: dStores.control, colour: "#d9622b"}
  ]
  // both panels share a vertical axis, so that the heights of their bars are comparable
  const tallest = () => d3.max(groups, g => d3.max(scores, s => g.store.get(s) ?? 0))
  const estimatesPlot = html`<div></div>`
  const grey = "#8a8a8a"
  // the estimates are plotted in three subplots, by type, so that each type has its own horizontal axis
  const subplots = [
   {title: "Standardized Mean Difference", names: ["Cohen's d"], unit: "Standardized units", include: [-1, 1]},
    {title: "Means", names: [...groups.map(g => `${g.short} mean`), "Mean difference"], unit: "Scale units", include: [0, 1]},
    {title: "Standard deviations", names: [...groups.map(g => `${g.short} SD`), "Pooled SD"], unit: "Scale units", include: [0, 1]}
  ]
  const update = () => {
    // rows of the estimates subplots: each condition's mean and SD, then the three effect size statistics
    const rows = []
    for (const g of groups) {
      const stats = g.panel.stats()
      g.title.textContent = `${g.label} (n = ${stats.n})`
      const [lower, upper] = meanCi(stats) ?? [stats.mean, stats.mean]
      if (stats.n > 0) rows.push({name: `${g.short} mean`, est: stats.mean, lower, upper, colour: g.colour})
      const [sdLower, sdUpper] = sdCi(stats) ?? [stats.sd, stats.sd]
      if (stats.n > 1) rows.push({name: `${g.short} SD`, est: stats.sd, lower: sdLower, upper: sdUpper, colour: g.colour})
    }
    const a = groups[0].panel.stats(), b = groups[1].panel.stats()
    const df = a.n + b.n - 2
    const pooled = Math.sqrt((a.ss + b.ss) / df)
    if (a.n > 0 && b.n > 0 && df > 0 && pooled > 0) {
      const diff = a.mean - b.mean
      const scale = Math.sqrt(1 / a.n + 1 / b.n)
      const tCrit = tQuantile(0.975, df)
      // 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
      const t = diff / (pooled * scale)
      const ncp = p => bisect(x => 1 - nctCdf(t, df, x), p, -10 - 2 * Math.abs(t), 10 + 2 * Math.abs(t))
      rows.push(
        // mean difference: t distribution
        {name: "Mean difference", est: diff, lower: diff - tCrit * pooled * scale, upper: diff + tCrit * pooled * scale, colour: grey},
        // pooled SD: chi-squared distribution
        {name: "Pooled SD", est: pooled, lower: pooled * Math.sqrt(df / chisqQuantile(0.975, df)), upper: pooled * Math.sqrt(df / chisqQuantile(0.025, df)), colour: grey},
        {name: "Cohen's d", est: diff / pooled, lower: ncp(0.025) * scale, upper: ncp(0.975) * scale, colour: grey}
      )
    }
    estimatesPlot.replaceChildren(...estimateSubplots(subplots, rows, width))
  }
  for (const g of groups) {
    g.title = html`<div class="sb-title" style="color: ${g.colour};"></div>`
    g.panel = frequencyPanel(scores, g.store, {tallest, onChange: ({rescaled}) => {
      if (rescaled) groups.forEach(other => { if (other !== g) other.panel.refresh() })
      update()
    }})
    g.panel.style.setProperty("--sb", g.colour)
  }
  dButtons.querySelector(".sb-clear").onclick = () => {
    groups.forEach(g => g.store.clear())
    groups.forEach(g => g.panel.refresh())
    update()
  }
  update()
  return html`<div>
    ${groups.map(g => html`<div>${g.title}${g.panel}</div>`)}
    <div class="sb-title">Estimates with 95% CIs</div>
    ${estimatesPlot}
  </div>`
}

// starting samples; observations outside the current scale bounds are kept but not shown or counted
dStores = ({
  intervention: new Map([[3, 2], [4, 4], [5, 7], [6, 5], [7, 2]]),
  control: new Map([[2, 2], [3, 5], [4, 7], [5, 4], [6, 2]])
})

function frequencyPanel(scores, store, {onChange = () => {}, tallest} = {}) {
  // a histogram of the counts in store (a Map of score -> count) that can be edited by drawing in the plot or by
  // typing counts under the axis. onChange is called after every edit, with rescaled: true if the vertical axis
  // may have changed. tallest optionally supplies the count that the vertical axis has to fit, for panels that
  // share an axis. the returned element has two extra methods: stats(), and refresh() to redraw after the store
  // has been changed from outside
  const count = s => store.get(s) ?? 0
  const clamp = (x, min, max) => Math.max(min, Math.min(max, x))

  // the plot is built once and then updated in place, so that it is not replaced in the middle of a drag
  const bars = scores.map(() => html`<div class="sb-bar"></div>`)
  const ghosts = scores.map(() => html`<div class="sb-ghost"><span></span></div>`)
  const grid = html`<div></div>`
  const plot = html`<div class="sb-plot">${grid}${scores.map((s, i) => html`<div class="sb-barcol">${ghosts[i]}${bars[i]}</div>`)}</div>`
  const inputs = scores.map(s => html`<input type="number" min="0" step="1" title="Number of observations at ${s}" aria-label="Number of observations at ${s}">`)
  // with more than 25 scale points, the columns are too narrow for the count boxes, so they are left out and
  // only the scale points that are multiples of 5 or 10 are labelled
  const showCounts = scores.length <= 25
  const labelStep = showCounts ? 1 : scores.length <= 50 ? 5 : 10
  const node = html`<div class="sample-builder">
    ${plot}
    <div class="sb-axis">${scores.map(s => html`<span>${s % labelStep === 0 ? s : ""}</span>`)}</div>
    ${showCounts ? html`<div class="sb-counts">${inputs}</div>` : ""}
  </div>`
  node.stats = () => sampleStats(new Map(scores.map(s => [s, count(s)])))

  // top of the vertical axis: a round number with some headroom above the tallest bar
  let yMax
  const rescale = () => {
    yMax = Math.max(10, d3.nice(0, (tallest ? tallest() : d3.max(scores, count)) * 1.25, 5)[1])
    grid.replaceChildren(...d3.ticks(0, yMax, 5).filter(d => d > 0).map(d => html`<div class="sb-gridline" style="bottom: ${d / yMax * 100}%;"><span>${d}</span></div>`))
  }
  const paint = () => {
    scores.forEach((s, i) => {
      bars[i].style.height = `${Math.min(1, count(s) / yMax) * 100}%`
      // leave the box the user is typing in alone
      if (inputs[i] !== document.activeElement) inputs[i].value = count(s)
    })
  }
  node.refresh = () => { rescale(); paint() }

  // drawing: the pointer's column gives the response option and its height gives the count
  const at = e => {
    const r = plot.getBoundingClientRect()
    return {
      i: clamp(Math.floor((e.clientX - r.left) / r.width * scores.length), 0, scores.length - 1),
      k: clamp(Math.round((r.bottom - e.clientY) / r.height * yMax), 0, yMax)
    }
  }
  const showGhost = (cur) => ghosts.forEach((g, i) => {
    g.style.display = cur && i === cur.i ? "block" : "none"
    if (cur && i === cur.i) { g.style.height = `${cur.k / yMax * 100}%`; g.firstChild.textContent = cur.k }
  })
  let last = null
  plot.onpointerdown = e => {
    if (e.button !== 0) return
    e.preventDefault()
    plot.setPointerCapture(e.pointerId)
    last = at(e)
    store.set(scores[last.i], last.k)
    showGhost(last)
    paint()
    onChange({rescaled: false})
  }
  plot.onpointermove = e => {
    const cur = at(e)
    showGhost(cur)
    if (!last) return
    // also set any columns that the pointer skipped over since the last event, interpolating their heights
    const step = Math.sign(cur.i - last.i)
    for (let i = last.i; i !== cur.i; i += step) {
      store.set(scores[i + step], Math.round(last.k + (cur.k - last.k) * (i + step - last.i) / (cur.i - last.i)))
    }
    store.set(scores[cur.i], cur.k)
    last = cur
    paint()
    onChange({rescaled: false})
  }
  plot.onpointerup = plot.onpointercancel = () => {
    if (!last) return
    last = null
    showGhost(null)
    node.refresh()
    onChange({rescaled: true})
  }
  plot.onpointerleave = () => { if (!last) showGhost(null) }

  // typing: exact counts
  inputs.forEach((input, i) => {
    input.oninput = () => {
      store.set(scores[i], clamp(Math.floor(+input.value) || 0, 0, 9999))
      node.refresh()
      onChange({rescaled: true})
    }
    input.onblur = paint
  })

  node.refresh()
  return node
}

function scaleScores(min, max) {
  // the response options of a scale, or null if the bounds are not usable
  return Number.isInteger(min) && Number.isInteger(max) && min < max && max - min <= 99 ? d3.range(min, max + 1) : null
}

function scaleMessage() {
  return html`<div style="color: var(--bs-secondary-color);">Set the minimum and maximum to whole numbers, with the minimum below the maximum and no more than 100 scale points.</div>`
}

function estimateSubplots(subplots, rows, width) {
  // one small plot of point estimates and 95% CIs per subplot ({title, names, unit, include}), drawing the rows
  // ({name, est, lower, upper, colour}) whose names it lists, or those assigned to it by a subplot property that
  // matches its title. include gives values that the axis always covers
  return subplots.map(s => {
    const data = rows.filter(d => d.subplot ? d.subplot === s.title : s.names.includes(d.name))
    return html`<div><div class="sb-subtitle">${s.title}</div>${Plot.plot({
      width, height: 55 + 32 * s.names.length, marginTop: 5, marginLeft: 130, marginRight: 20, marginBottom: 45,
      x: {label: s.unit, labelAnchor: "center", labelOffset: 38, grid: true, nice: true, ticks: 5,
          domain: d3.extent([...s.include, ...data.flatMap(d => [d.lower, d.upper])])},
      y: {label: null, type: "point", domain: s.names, padding: 0.5, tickSize: 0},
      style: {fontSize: "13px"},
      marks: [
        Plot.ruleX([0], {stroke: "currentColor", strokeDasharray: "4,3", strokeOpacity: 0.6}),
        Plot.link(data, {x1: "lower", x2: "upper", y1: "name", y2: "name", stroke: "colour", strokeWidth: 2.5}),
        Plot.dot(data, {x: "est", y: "name", fill: "colour", r: 6, symbol: "square"})
      ]
    })}</div>`
  })
}

function meanCi({n, mean, sd}) {
  // 95% CI for a mean, or null if it cannot be calculated
  if (!(n > 1 && sd > 0)) return null
  const margin = tQuantile(0.975, n - 1) * sd / Math.sqrt(n)
  return [mean - margin, mean + margin]
}

function sdCi({n, sd}) {
  // 95% CI for an SD, or null if it cannot be calculated
  if (!(n > 1 && sd > 0)) return null
  return [sd * Math.sqrt((n - 1) / chisqQuantile(0.975, n - 1)), sd * Math.sqrt((n - 1) / chisqQuantile(0.025, n - 1))]
}

function sampleStats(counts) {
  // counts: Map of score -> number of observations. ss is the sum of squared deviations from the mean,
  // and sd is the sample SD (n - 1 denominator)
  const entries = [...counts]
  const n = d3.sum(entries, ([, k]) => k)
  const mean = d3.sum(entries, ([s, k]) => s * k) / n
  const ss = n > 0 ? d3.sum(entries, ([s, k]) => k * (s - mean) ** 2) : 0
  return {n, mean, ss, sd: Math.sqrt(ss / (n - 1))}
}

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 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
}

explorerStyle = html`<style>
  .sample-builder { --sb: #2780e3; margin: 1rem 0; padding-left: 2rem; }
  .sample-builder .sb-plot { position: relative; display: flex; height: 240px; border-bottom: 1px solid currentColor; cursor: crosshair; touch-action: none; user-select: none; -webkit-user-select: none; }
  .sample-builder .sb-gridline { position: absolute; left: 0; right: 0; border-top: 1px solid rgba(128, 128, 128, 0.25); pointer-events: none; }
  .sample-builder .sb-gridline span { position: absolute; left: -2rem; width: 1.6rem; top: -0.6em; text-align: right; font-size: 11px; line-height: 1.2; }
  .sample-builder .sb-barcol { flex: 1 1 0; min-width: 0; position: relative; }
  .sample-builder .sb-bar, .sample-builder .sb-ghost { position: absolute; bottom: 0; left: 15%; right: 15%; border-radius: 2px 2px 0 0; pointer-events: none; }
  .sample-builder .sb-bar { background: var(--sb); }
  .sample-builder .sb-ghost { border: 1px dashed var(--sb); border-bottom: none; display: none; }
  .sample-builder .sb-ghost span { position: absolute; top: -1.4em; left: -1rem; right: -1rem; text-align: center; font-size: 11px; color: var(--sb); }
  .sample-builder .sb-axis { display: flex; font-size: 12px; }
  .sample-builder .sb-axis span { flex: 1 1 0; min-width: 0; text-align: center; }
  .sample-builder .sb-counts { display: flex; margin-top: 2px; }
  .sample-builder .sb-counts input { flex: 1 1 0; min-width: 0; margin: 0 2px; padding: 1px 0; text-align: center; font-size: 12px; color: inherit; background: transparent; border: 1px solid rgba(128, 128, 128, 0.4); border-radius: 3px; appearance: textfield; -moz-appearance: textfield; }
  .sample-builder .sb-counts input:focus { outline: none; border-color: var(--sb); }
  .sample-builder .sb-counts input::-webkit-inner-spin-button, .sample-builder .sb-counts input::-webkit-outer-spin-button { -webkit-appearance: none; margin: 0; }
  .sb-title { font-weight: 600; margin-top: 1.5rem; }
  .sb-subtitle { font-size: 0.9rem; font-weight: 600; margin: 0.75rem 0 0.25rem; color: var(--bs-secondary-color); }
</style>`

function chisqQuantile(p, df) {
  return bisect(x => gammaP(df / 2, x / 2), p, 0, df + 20 * Math.sqrt(2 * df) + 50)
}
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 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 gammaP(a, x) {
  // regularised lower incomplete gamma function: series for small x, continued fraction otherwise
  if (x <= 0) return 0
  const front = Math.exp(-x + a * Math.log(x) - logGamma(a))
  if (x < a + 1) {
    let term = 1 / a, sum = term
    for (let n = 1; n < 500; n++) { term *= x / (a + n); sum += term; if (Math.abs(term) < Math.abs(sum) * 1e-14) break }
    return sum * front
  }
  let b = x + 1 - a, c = 1e300, d = 1 / b, h = d
  for (let i = 1; i < 500; i++) {
    const an = -i * (i - a)
    b += 2
    d = an * d + b; if (Math.abs(d) < 1e-300) d = 1e-300
    c = b + an / 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 1 - front * h
}

8.5 Maximum positive controls and intuition benchmarks

• Using known effects as anchors for plausibility assessment