Power Analysis – MetaLab

How many infants does your study need? Ground your planning in everything the literature already knows: pick a phenomenon, optionally condition on age and method variables, and get the meta-analytically estimated effect size with the sample size needed to detect it at 80% power (α = .05, two-sided, within-subjects normal approximation).

Loading data…

pnorm = (x) => {
  const t = 1 / (1 + 0.3275911 * Math.abs(x) / Math.SQRT2);
  const y = 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 + y) : 0.5 * (1 - y);
}

qnorm975 = 1.959963984540054

// Lanczos log-gamma (self-recursion wrapped: OJS cells cannot reference
// their own name)
lgamma = {
  const g = [676.5203681218851, -1259.1392167224028, 771.32342877765313,
    -176.61502916214059, 12.507343278686905, -0.13857109526572012,
    9.9843695780195716e-6, 1.5056327351493116e-7];
  const f = (z) => {
    if (z < 0.5) return Math.log(Math.PI / Math.sin(Math.PI * z)) - f(1 - z);
    z -= 1;
    let x = 0.99999999999980993;
    for (let i = 0; i < 8; i++) x += g[i] / (z + i + 1);
    const t = z + 7.5;
    return 0.5 * Math.log(2 * Math.PI) + (z + 0.5) * Math.log(t) - t + Math.log(x);
  };
  return f;
}

// regularized incomplete beta I_x(a, b) by continued fraction (Lentz)
ibeta = {
  const f = (x, a, b) => {
    if (x <= 0) return 0;
    if (x >= 1) return 1;
    if (x > (a + 1) / (a + b + 2)) return 1 - f(1 - x, b, a);
    const lbeta = lgamma(a) + lgamma(b) - lgamma(a + b);
    const front = Math.exp(a * Math.log(x) + b * Math.log(1 - x) - lbeta) / a;
    let h = 1, c = 1, d = 0;
    for (let i = 0; i <= 300; i++) {
      const m = Math.floor(i / 2);
      let numerator;
      if (i === 0) numerator = 1;
      else if (i % 2 === 0) numerator = (m * (b - m) * x) / ((a + 2 * m - 1) * (a + 2 * m));
      else numerator = -((a + m) * (a + b + m) * x) / ((a + 2 * m) * (a + 2 * m + 1));
      d = 1 + numerator * d;
      if (Math.abs(d) < 1e-30) d = 1e-30;
      d = 1 / d;
      c = 1 + numerator / c;
      if (Math.abs(c) < 1e-30) c = 1e-30;
      h *= c * d;
      if (Math.abs(1 - c * d) < 1e-9) break;
    }
    return front * (h - 1);
  };
  return f;
}

// two-sided p for a t statistic with df degrees of freedom
t_pvalue = (t, df) => {
  const x = df / (df + t * t);
  return ibeta(x, df / 2, 0.5);
}

// upper-tail p for an F statistic
f_pvalue = (F, df1, df2) => {
  if (F <= 0) return 1;
  return ibeta(df2 / (df2 + df1 * F), df2 / 2, df1 / 2);
}

// standard normal draws (Box-Muller)
rnorm = (n, mean = 0, sd = 1) => {
  const out = new Array(n);
  for (let i = 0; i < n; i += 2) {
    const u = Math.random() || 1e-12, v = Math.random();
    const r = Math.sqrt(-2 * Math.log(u));
    out[i] = mean + sd * r * Math.cos(2 * Math.PI * v);
    if (i + 1 < n) out[i + 1] = mean + sd * r * Math.sin(2 * Math.PI * v);
  }
  return out;
}

mean_ = (xs) => xs.reduce((a, b) => a + b, 0) / xs.length
sd_ = (xs) => {
  const m = mean_(xs);
  return Math.sqrt(xs.reduce((a, b) => a + (b - m) * (b - m), 0) / (xs.length - 1));
}

// paired t test -> two-sided p
paired_t_p = (x, y) => {
  const d = x.map((v, i) => v - y[i]);
  const t = mean_(d) / (sd_(d) / Math.sqrt(d.length));
  return t_pvalue(t, d.length - 1);
}

// power of a two-sided one-sample z test with effect h at size n
// (the exact formula behind pwr::pwr.p.test, as the legacy app used)
power_z = (h, n) =>
  pnorm(Math.sqrt(n) * Math.abs(h) - qnorm975) +
  pnorm(-Math.sqrt(n) * Math.abs(h) - qnorm975)

// smallest n achieving target power (legacy "N for 80% power")
n_for_power = (h, target = 0.8, nmax = 100000) => {
  if (Math.abs(h) < 1e-9) return Infinity;
  let lo = 2, hi = 4;
  while (power_z(h, hi) < target && hi < nmax) { lo = hi; hi *= 2; }
  if (hi >= nmax) return Infinity;
  for (let i = 0; i < 60; i++) {
    const mid = (lo + hi) / 2;
    if (power_z(h, mid) < target) lo = mid; else hi = mid;
  }
  return hi;
}

// critical t (quantile) via bisection on t_pvalue
qt975 = (df) => {
  let lo = 1, hi = 300;
  for (let i = 0; i < 80; i++) {
    const mid = (lo + hi) / 2;
    if (t_pvalue(mid, df) > 0.05) lo = mid; else hi = mid;
  }
  return (lo + hi) / 2;
}

html`<style>
  .${ojs_ns} div label {
    background-color: #eef2f6;
    padding: 0.2rem 0.5rem;
    border-radius: 0.4rem;
    margin-bottom: 0.25rem;
    width: auto;
  }

  .${ojs_ns} select,
  .${ojs_ns} input[type="text"],
  .${ojs_ns} input[type="file"] {
    background-color: #fff;
    border: 1px solid #d8dee4;
    border-radius: 0.4rem;
    padding: 0.2rem 0.4rem;
  }

  .${ojs_ns} input[type="checkbox"] {
    accent-color: #2780e3;
    margin-bottom: 0;
  }
</style>`

html`<style>
  .${ojs_ns} button {
    background-color: #fff;
    color: #212529;
    border: 1px solid #c6ced6;
    border-radius: 0.4rem;
    padding: 0.25rem 0.9rem;
    cursor: pointer;
  }

  .${ojs_ns} button:hover {
    background-color: #eef2f6;
  }
</style>`

unpack_vcov = (lower, p) => {
  const V = Array.from({ length: p }, () => new Array(p).fill(0));
  let idx = 0;
  for (let j = 0; j < p; j++)
    for (let i = j; i < p; i++) {
      V[i][j] = V[j][i] = lower[idx++];
    }
  return V;
}

// design vector for one data row (or synthetic row of moderator values)
design_row = (row, combo, mods) => {
  return combo.coefs.map((c) => {
    if (c.name === "intrcpt") return 1;
    if (mods.includes(c.name)) return +row[c.name];  // numeric moderator
    for (const m of Object.keys(combo.levels ?? {})) {
      if (c.name.startsWith(m))
        return String(row[m]) === c.name.slice(m.length) ? 1 : 0;
    }
    return 0;
  });
}

dot = (a, b) => a.reduce((s, v, i) => s + v * b[i], 0)

quad_form = (x, V) => {
  let s = 0;
  for (let i = 0; i < x.length; i++)
    for (let j = 0; j < x.length; j++) s += x[i] * V[i][j] * x[j];
  return s;
}

// R's make.unique: append .1, .2, ... to repeats
make_unique = (names) => {
  const seen = new Map();
  return names.map((n) => {
    if (!seen.has(n)) { seen.set(n, 0); return n; }
    const k = seen.get(n) + 1; seen.set(n, k);
    return `${n}.${k}`;
  });
}

fmt_p = (p) => p < 0.0001 ? "< .0001" : "= " + p.toFixed(4).replace(/^0/, "")

days_per_month = 365.2425 / 12

registry = d3.json("slices/datasets.json")

domain_titles = new Map([
  ["early_language", "Early language"],
  ["cognitive_development", "Cognitive development"]
])

standard_mods = ["mean_age", "response_mode", "exposure_phase"]
mod_labels = ({ mean_age: "Age", response_mode: "Response mode",
  exposure_phase: "Exposure phase" })

viewof domain = Inputs.select([...domain_titles.keys()], {
  label: "Domain", format: (d) => domain_titles.get(d)
})

viewof dataset_name = Inputs.select(
  registry.filter((d) => d.domain === domain).map((d) => d.name),
  { label: "Meta-analysis" })

dataset_meta = registry.find((d) => d.name === dataset_name)

viewof subset_choice = {
  const subsets = dataset_meta.subset ?? [];
  if (!subsets.length) {
    // no subsets for this dataset: fixed value, no control shown
    const el = html`<div></div>`;
    el.value = "All data";
    return el;
  }
  return Inputs.radio(["All data", ...subsets], {
    label: "Subset",
    value: "All data"
  });
}

slice_all = d3.json(`slices/es/${dataset_meta.short_name}.json`)

models_all = d3.json(`slices/models/${dataset_meta.short_name}.json`)

slice = subset_choice === "All data"
  ? slice_all
  : slice_all.filter((d) => d[subset_choice] === true)

pwr_rows = slice.filter((d) =>
  Number.isFinite(+d.d_calc) && Number.isFinite(+d.d_var_calc))

// standard moderators available for this dataset (>1 observed value; age
// hidden for longitudinal datasets, as in the legacy app)
pwr_pool = standard_mods.filter((m) => {
  if (m === "mean_age" && dataset_meta.longitudinal) return false;
  const vals = new Set(pwr_rows.map((d) => d[m]).filter((v) => v !== null && v !== undefined && v !== ""));
  return vals.size > 1;
})

viewof pwr_mods = Inputs.checkbox(pwr_pool, {
  label: "Condition on",
  format: (m) => mod_labels[m] ?? m
})

viewof age_months = {
  if (!pwr_mods.includes("mean_age"))
    return Inputs.input(null);
  const ages = pwr_rows.map((d) => +d.mean_age_months).filter(Number.isFinite);
  return Inputs.range([0, Math.ceil(d3.max(ages))], {
    label: "Age (months)", step: 1, value: Math.round(d3.mean(ages))
  });
}

viewof response_mode_choice = pwr_mods.includes("response_mode")
  ? Inputs.select([...new Set(pwr_rows.map((d) => d.response_mode).filter(Boolean))],
      { label: "Response mode" })
  : Inputs.input(null)

viewof exposure_phase_choice = pwr_mods.includes("exposure_phase")
  ? Inputs.select([...new Set(pwr_rows.map((d) => d.exposure_phase).filter(Boolean))],
      { label: "Exposure phase" })
  : Inputs.input(null)

The estimate comes from the same multilevel random-effects model as the visualization page (the legacy app used a single-level model here; see the changelog).

mods_key = standard_mods.filter((m) => pwr_mods.includes(m)).join(",")

combo = models_all[`${subset_choice}|d|${mods_key}`]

d_pwr = {
  if (!combo || !combo.converged) return null;
  const synth = ({
    mean_age: (age_months ?? 0) * days_per_month,
    response_mode: response_mode_choice,
    exposure_phase: exposure_phase_choice
  });
  const X = design_row(synth, combo, pwr_mods);
  return dot(X, combo.coefs.map((c) => c.est));
}

n80 = d_pwr === null ? null : n_for_power(d_pwr, 0.8)

power_boxes = {
  clear_loading();
  if (d_pwr === null)
    return html`<p style="color:#777;"><em>No converged model for this selection.</em></p>`;
  return html`<div style="display: flex; gap: 0.75rem; flex-wrap: wrap; margin-bottom: 0.75rem;">
    <div class="value-box" style="background:#dd4b39;">
      <div class="vb-value">${d_pwr.toFixed(2)}</div>
      <div class="vb-label">Effect size (d)</div></div>
    <div class="value-box" style="background:#dd4b39;">
      <div class="vb-value">${!Number.isFinite(n80) || n80 >= 200 ? "> 200" : Math.ceil(n80)}</div>
      <div class="vb-label">N for 80% power</div></div>
  </div>`;
}

Power as a function of sample size

N is the number of infants per group: for a within-participant design (most MetaLab datasets) that is the total number of infants; for a between-participant design, each condition needs N infants.

power_plot = {
  if (d_pwr === null) return html``;
  const n90 = n_for_power(d_pwr, 0.9);
  const max_n = Math.min(Math.max(60, Number.isFinite(n90) ? Math.ceil(n90) : 200), 200);
  const ns = d3.range(5, max_n + 1, 5);
  const curve = ns.map((n) => ({ n, power: power_z(d_pwr, n) }));
  return Plot.plot({
    width: 780, height: 430,
    marginLeft: 55,
    x: { label: "Number of participants per group (N)", domain: [0, max_n] },
    y: { label: "Power to reject the null at p < .05", domain: [0, 1] },
    marks: [
      Plot.ruleY([0.8], { strokeDasharray: "5,4", stroke: "#666" }),
      Number.isFinite(n80) && n80 <= max_n
        ? Plot.ruleX([n80], { strokeDasharray: "2,3", stroke: "#666" }) : null,
      Plot.line(curve, { x: "n", y: "power", stroke: "#333" }),
      Plot.dot(curve, { x: "n", y: "power", fill: "#333", r: 2.5,
        title: (d) => `N = ${d.n}: power = ${d.power.toFixed(2)}`, tip: true })
    ].filter(Boolean)
  });
}

d_pwr === null ? html`` : Inputs.table(
  d3.range(5, 205, 5).map((n) => ({
    N: n, power: +power_z(d_pwr, n).toFixed(3)
  })), { rows: 15 })

Power is computed for a two-sided test at α = .05 using the normal approximation (the same formula as pwr::pwr.p.test with h = d, as in the legacy application): a reasonable approximation for a within-subjects design. For between-subjects designs, required samples are substantially larger.