Power Simulation

What does an infant looking-time experiment look like when the true effect is known? Simulate an experiment: pick a sample size and a true effect size, then see whether the simulated study comes out significant. Resample to feel how much results vary — especially at small N.

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

// input styling adapted from levante-datapage: pill labels, borderless
// rounded selects (scoped to the ojs namespace class)
ojs_ns = Inputs.text().classList[0]

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

sim_mu = 10
sim_sd = 5
sim_colors = ({ "Longer looking predicted": "#cb4b16", "Shorter looking predicted": "#268bd2" })

viewof sim_n = Inputs.range([4, 120], { label: "Infants per group (N)", step: 2, value: 16 })

viewof sim_d = Inputs.range([0, 2], { label: "Effect size (Cohen's d)", step: 0.1, value: 0.5 })

viewof sim_control = Inputs.radio(
  new Map([["Experimental only", false], ["Experimental & control", true]]),
  { label: "Conditions", value: false })

viewof sim_interval = Inputs.radio(
  new Map([["95% confidence interval", "ci"], ["Standard error of the mean", "sem"]]),
  { label: "Error bars", value: "ci" })

viewof sim_plot_type = Inputs.radio(["Bar graph", "Scatter plot"], {
  label: "Plot type", value: "Bar graph" })

viewof sim_go = Inputs.button("Sample again")

sim_cells = {
  sim_go;  // re-run on button press
  const groups = sim_control ? ["Experimental", "Control"] : ["Experimental"];
  const conds = ["Longer looking predicted", "Shorter looking predicted"];
  const cells = [];
  for (const g of groups) {
    for (const c of conds) {
      // true effect only in the experimental group's "longer" condition,
      // matching the legacy app (which, however, accidentally duplicated
      // every draw when the control group was on -- fixed here)
      const mu = g === "Experimental" && c === "Longer looking predicted"
        ? sim_mu + (sim_d * sim_sd) / 2
        : sim_mu - (sim_d * sim_sd) / 2;
      cells.push({ group: g, condition: c, values: rnorm(sim_n, mu, sim_sd) });
    }
  }
  return cells;
}

sim_points = sim_cells.flatMap((c) =>
  c.values.map((v) => ({ group: c.group, condition: c.condition, looking_time: v })))

sim_summary = sim_cells.map((c) => {
  const m = mean_(c.values);
  const sem = sd_(c.values) / Math.sqrt(c.values.length);
  const half = sim_interval === "sem" ? sem : qt975(c.values.length - 1) * sem;
  return { group: c.group, condition: c.condition, mean: m,
           lo: m - half, hi: m + half };
})

Simulated data

sim_plot = {
  // dodged layout on a continuous x: group centers at 1, 2; conditions
  // offset ±0.22 within group (band scales in Plot don't dodge)
  const groups = sim_control ? ["Experimental", "Control"] : ["Experimental"];
  const conds = Object.keys(sim_colors);
  const xpos = (g, c) => groups.indexOf(g) + 1 + (conds.indexOf(c) === 0 ? -0.22 : 0.22);
  const summary = sim_summary.map((d) => ({ ...d, x: xpos(d.group, d.condition) }));
  const points = sim_points.map((d) => ({
    ...d, x: xpos(d.group, d.condition) + (Math.random() - 0.5) * 0.12 }));
  const ymax = Math.max(20, d3.max(summary, (d) => d.hi),
    sim_plot_type === "Bar graph" ? 0 : d3.max(points, (d) => d.looking_time));
  const ymin = sim_plot_type === "Bar graph" ? 0 :
    Math.min(0, d3.min(points, (d) => d.looking_time));
  return Plot.plot({
    width: 720, height: 420,
    marginLeft: 55,
    x: { label: null, domain: [0.4, groups.length + 0.6],
         ticks: groups.map((g, i) => i + 1),
         tickFormat: (t) => groups[t - 1] ?? "" },
    y: { label: "Simulated looking time (s)", domain: [ymin, ymax] },
    color: { domain: conds, range: Object.values(sim_colors), legend: true },
    marks: [
      sim_plot_type === "Bar graph"
        ? Plot.rect(summary, { x1: (d) => d.x - 0.18, x2: (d) => d.x + 0.18,
            y1: 0, y2: "mean", fill: "condition" })
        : Plot.dot(points, { x: "x", y: "looking_time", fill: "condition",
            fillOpacity: 0.45, r: 3.5 }),
      Plot.ruleX(summary, { x: "x", y1: "lo", y2: "hi", stroke: "black",
        strokeWidth: 1.5 }),
      sim_plot_type === "Bar graph" ? null :
        Plot.ruleY(summary, { y: "mean", x1: (d) => d.x - 0.15,
          x2: (d) => d.x + 0.15, stroke: "black", strokeWidth: 2 })
    ].filter(Boolean)
  });
}

sim_stats = {
  const cell = (g, c) => sim_cells.find((x) => x.group === g && x.condition === c).values;
  const p_exp = paired_t_p(cell("Experimental", "Longer looking predicted"),
                           cell("Experimental", "Shorter looking predicted"));
  const fmt = (p) => p < 0.001 ? p.toExponential(2) : p.toFixed(3);
  let out = `A paired t test of the experimental condition is
    <b>${p_exp < 0.05 ? "significant" : "non-significant"}</b> at p = ${fmt(p_exp)}.`;
  if (sim_control) {
    const p_ctl = paired_t_p(cell("Control", "Longer looking predicted"),
                             cell("Control", "Shorter looking predicted"));
    out += ` The control condition is ${p_ctl < 0.05 ? "significant" : "non-significant"}
      at p = ${fmt(p_ctl)}.`;
    // 2x2 ANOVA interaction (equal cell sizes)
    const cells = sim_cells;
    const all = cells.flatMap((c) => c.values);
    const gm = mean_(all);
    const n = sim_n;
    const cellMeans = cells.map((c) => mean_(c.values));
    const gMeans = ["Experimental", "Control"].map((g) =>
      mean_(cells.filter((c) => c.group === g).flatMap((c) => c.values)));
    const cMeans = ["Longer looking predicted", "Shorter looking predicted"].map((cc) =>
      mean_(cells.filter((c) => c.condition === cc).flatMap((c) => c.values)));
    let ss_int = 0;
    cells.forEach((c) => {
      const gi = c.group === "Experimental" ? 0 : 1;
      const ci = c.condition === "Longer looking predicted" ? 0 : 1;
      const dev = mean_(c.values) - gMeans[gi] - cMeans[ci] + gm;
      ss_int += n * dev * dev;
    });
    let ss_w = 0;
    cells.forEach((c) => {
      const m = mean_(c.values);
      c.values.forEach((v) => ss_w += (v - m) * (v - m));
    });
    const df_w = 4 * n - 4;
    const F = (ss_int / 1) / (ss_w / df_w);
    const p_int = f_pvalue(F, 1, df_w);
    out += ` A group × condition ANOVA ${p_int < 0.05 ? "shows" : "does not show"}
      a significant interaction at p = ${fmt(p_int)}.`;
  }
  return html`<p style="font-size: 1.05rem;">${html([out])}</p>`;
}