Distribution of significant effects

math = require("mathjs")
norm = import('https://unpkg.com/norm-dist@3.1.0/index.js?module')

quad_points = JSON.parse(all_quad_points).at(qp - 1)

sigma = 2 / math.sqrt(ESS)

zeta_sd = math.sqrt(tau**2 + rho * sigma**2)

crit = norm.icdf(1 - alpha)

binomial_coefs = Array(m+1).fill(null).map((x,index) => {
  return math.combinations(m, index);
})

probs = quad_points.map(zeta => {
  let Z = (zeta[0] * zeta_sd + mu - crit * sigma) / math.sqrt(omega**2 + (1 - rho) * sigma**2);
  return [norm.cdf(Z), zeta[1]];
})

p_binom_norm = binomial_coefs.map((coef, a) => {
  let p = probs.map((x) => {
    return (x[0]**a) * ((1 - x[0])**(m - a)) * x[1];
  });
  return coef * math.sum(p);
})

Plot.plot({
  x: {
    label: "Number of significant effect sizes"
  },
  y: {
    domain: [0, 1],
    label: "Probability"
  },
  marks: [
    Plot.ruleY(0),
    Plot.barY(p_binom_norm, {
      fill: "steelblue"
    }),
  ]
})

viewof m = Inputs.range(
  [1, 30], 
  {value: 6, step: 1, label: "Number of effect sizes (m):"}
)

viewof ESS = Inputs.range(
  [4, 300], 
  {value: 80, step: 1, label: "Effective sample size:"}
)

viewof mu = Inputs.range(
  [-2, 2], 
  {value: 0.3, step: 0.01, label: "Average effect size (mu):"}
)

viewof tau = Inputs.range(
  [0, 1], 
  {value: 0.1, step: 0.01, label: "Between-study SD (tau):"}
)

viewof omega = Inputs.range(
  [0, 1], 
  {value: 0.1, step: 0.01, label: "Within-study SD (omega):"}
)

viewof rho = Inputs.range(
  [0, 1], 
  {value: 0.6, step: 0.01, label: "Correlation of sampling errors (rho):"}
)

viewof alpha = Inputs.range(
  [0.005, 0.995], 
  {value: 0.025, step: .005, label: "One-sided significance threshold (alpha):"}
)

viewof qp = Inputs.range(
  [1, 30], 
  {value: 21, step: 1, label: "Number of quadrature points:"}
)