"""False-positive rate when cells are counted as n, /guides/what-counts-as-n/.

No true treatment effect. 2 groups x 5 mice x 20 cells per mouse; total variance 1, split into
mouse-to-mouse and cell-to-cell parts by the intraclass correlation (ICC). Each simulated experiment is
analysed twice with Welch's t-test at alpha = 0.05: on all cells (wrong unit) and on mouse means.
100,000 experiments per ICC, seed 20261110. Python 3.13, NumPy, SciPy 1.18.1. Run: python simulation.py
"""
import json

import numpy as np
from scipy import stats

rng = np.random.default_rng(20261110)
N_SIM, MICE, CELLS, ALPHA = 100_000, 5, 20, 0.05
results = []
for icc in (0.0, 0.05, 0.1, 0.2, 0.5):
    sd_mouse, sd_cell = np.sqrt(icc), np.sqrt(1 - icc)
    hits_cells = hits_means = 0
    for start in range(0, N_SIM, 10_000):
        k = min(10_000, N_SIM - start)
        mouse = rng.normal(0, sd_mouse, size=(k, 2, MICE, 1))
        y = mouse + rng.normal(0, sd_cell, size=(k, 2, MICE, CELLS))
        cells = y.reshape(k, 2, MICE * CELLS)
        p_cells = stats.ttest_ind(cells[:, 1], cells[:, 0], axis=1, equal_var=False).pvalue
        means = y.mean(axis=3)
        p_means = stats.ttest_ind(means[:, 1], means[:, 0], axis=1, equal_var=False).pvalue
        hits_cells += int((p_cells < ALPHA).sum())
        hits_means += int((p_means < ALPHA).sum())
    results.append({"icc": icc, "design_effect": 1 + (CELLS - 1) * icc,
                    "false_positive_cells_as_n": hits_cells / N_SIM,
                    "false_positive_mouse_means": hits_means / N_SIM})

out = {"n_sim": N_SIM, "mice_per_group": MICE, "cells_per_mouse": CELLS, "alpha": ALPHA,
       "seed": 20261110, "results": results}
with open("simulation_output.json", "w") as f:
    json.dump(out, f, indent=2)
print(json.dumps(out, indent=2))