Simulation.md

Simulation

Cohorts where we chose the answer in advance, run through the identical workflow. On real data you cannot tell whether a null means nothing is there or that the design could not see it. Here you can.

Why it exists

Two questions the real cohorts cannot answer:

• when the pipeline reports nothing, is that the data or the method?
• when it reports something, would it have reported the same thing from noise?

Both need a dataset whose truth is known.

What gets planted

causal taxa	8 species that genuinely move the outcome. The right answer.
confounded taxa	12 species that track the outcome and cause none of it
conditions	effect, where the outcome responds to the 8, and no_effect, where it does not
cohorts	10 per condition, same seeds, so the two arms share a community

The confounded taxa are the point. Without them a method could succeed by flagging whatever correlates with the outcome, so the task is discrimination rather than detection.

no_effect is the important arm: every parameter and seed matches its effect twin and only the outcome's response is switched off, so anything reported there is a false positive by construction.

Running it

mbcausal sim list                                  # the 14 experiments
mbcausal sim generate --root <dir> --n-seeds 10    # write cohorts, four stages each
mbcausal sim run rtm_demo --n-rep 300              # one experiment, flags pass through

generate refuses to overwrite a populated directory. Pass --force if you mean it, or --root to write somewhere new.

Each cohort is written in pipeline order, so you can enter at whichever stage you want:

<root>/<condition>/seedNN/
  1_raw/          per sample abundance, before any QC
  2_cohort/       abundance.csv + sample_table.csv, what preprocess reads
  3_preprocess/   the CLR feature matrix
  4_discovery/    our own discovery output, for comparison
  truth/          the answer key: which species were planted, and the true effect size

Stages 3 and 4 are ours. Regenerate them yourself from 2_cohort if you would rather not trust them. truth/ can be withheld for a blind evaluation.

How a run is scored

Each recovered factor gets a causal taxa share: the fraction of its loading mass sitting on the 8 planted species. A factor above 0.30 is the planted cause recovered. This is computed after the fact, so the pipeline never learns which species were planted.

Where output goes

SIM_ROOT and PROJECT_SIM_ROOT in src/mbcausal/sim/paths.py. Both derive from settings rather than being hardcoded: SIM_ROOT follows config.preprocess_dir, and PROJECT_SIM_ROOT follows the repo location. Point preprocess_dir somewhere else in mbcausal.yaml and the simulation output moves with it.

How it is built

Six steps, in the order they happen. The equations, the calibration and which parameters were fitted are written up separately and ship with the exported cohorts they describe.

1. A community. 106 patients, 26 visits two weeks apart, 2000 species, of which 1600 to 1740 clear the detection threshold. Each abundance is a species mean, a persistent offset for that patient, and a fluctuation at that visit.
2. Patient offsets. Built from 30 shared axes rather than independently per species, because real people differ along a few community wide gradients. Against an independent background NMF never assembles the planted species into a factor at all.
3. The cause. The 8 species share a trajectory that carries over week to week, and the outcome responds at a lag.
4. The confounded species. The 12 track a non microbial driver of the outcome.
5. The assay schedule. Measured every fifth visit, then degraded by missed appointments, scheduling slop and early dropout, giving a ten week median gap that matches the real cohort.
6. Calibration. Parameters set by grid search so seven summary statistics of the simulated tables match the iHMP MetaPhlAn4 ones to about 4 percent. This is the step that ties the whole generator to one cohort.

The rule that makes it worth anything

sim/generate/ does not import the rest of the package. If the generator called the code under test, a shared defect could cancel itself out and the simulation would validate nothing. Everything downstream of generation, sim/evaluate and sim/experiments, does import mbcausal, deliberately, because it is running the real pipeline.

Limits worth knowing before quoting a number

• Every parameter is calibrated to iHMP and to calprotectin. The generator itself is not specific to a cohort, but its defaults were fitted to the iHMP MetaPhlAn4 tables, so a number measured here says nothing about a cohort of a different size, sampling density or outcome distribution until those parameters are refitted.
• Species are conditionally independent, so there are no ecological interactions.
• No count sampling step, so sparsity is a threshold rather than a consequence of depth.
• Patient offsets do not change over time.
• A cohort yields 250 to 308 usable lagged pairs, median 283, against the real 156, so any detection rate measured here is an upper bound.

See also

• Configuration for where output lands
• src/mbcausal/sim/README.md for the package layout