mims-harvard/tooluniverse/tooluniverse-statistical-modeling

Statistical Modeling for Biomedical Data Analysis

RULE ZERO — Check for pre-computed results FIRST

Before following any instruction below, scan the data folder for:

• *_executed.ipynb → read with tu run read_executed_notebook '{"data_folder":"<path>","search":"<keyword>"}' and cite its cell outputs as the authoritative answer
• Pre-computed result files (CSV/TSV with names like *results*, *deseq*, *enrich*, *stats*, *_simplified.csv) → read directly and report the requested value
• Canonical analysis scripts (analysis.R, run_*.py, find_*.R, *.Rmd) → execute as-is and read the output

Only follow this skill's re-analysis recipe below if none of the above exist. Re-running from raw data produces different numbers than the published answer and is much slower (often 5-10× turn count).

PRIMARY SCRIPTS — use these FIRST

These scripts encode the question-specific gotchas in scripts/ and emit labelled, parseable output. Prefer them over ad-hoc statsmodels / scipy code.

Concrete invocations

Natural-spline regression (R^2, overall F-test p, peak Y + 95% CI):

python skills/tooluniverse-statistical-modeling/scripts/r_natural_spline_regression.py \
  --csv data.csv --y-col Area \
  --ratio-col Ratio --new-x-col Frequency_strain \
  --filter "StrainNumber not in ['1', '98']" \
  --df 4 --workdir /tmp/spline_run

Quadratic vs cubic vs natural-spline comparison + best-model peak:

python skills/tooluniverse-statistical-modeling/scripts/spline_model_compare.py \
  --csv data.csv --y-col Area \
  --ratio-col Ratio --new-x-col Frequency_strain \
  --filter "StrainNumber not in ['1', '98']" \
  --ns-df 4 --workdir /tmp/spline_cmp

Report the peak location (`x*`) in the units of the fitted x-variable, not a derived label. When the model is fit on a frequency/proportion column (e.g. Frequency_strain, a 0–1 value), the answer to "at what ratio/frequency is the maximum" is that fraction (e.g. 0.909), NOT the colon-ratio it was derived from (e.g. 10:1). Convert a colon ratio a:b to the fraction a/(a+b) when the question expects a 0–1 value or the fitted x-column is a fraction.

Ordinal logistic regression with interaction term (e.g. trial AE severity):

python skills/tooluniverse-statistical-modeling/scripts/logistic_regression_or.py \
  --csv merged.csv --outcome AESEV --outcome-type ordinal --outcome-order "1,2,3,4" \
  --predictors TRTGRP,expect_interact,patients_seen,MHONGO \
  --encode TRTGRP,expect_interact,patients_seen \
  --encode-map "TRTGRP:Placebo=0,BCG=1" \
  --interaction MHONGO:TRTGRP_cat \
  --coef-name TRTGRP_cat

Two-sample power analysis from a pilot CSV:

python skills/tooluniverse-statistical-modeling/scripts/power_analysis.py \
  --csv pilot.csv --value-col MeasuredValue --group-col Group \
  --group-a Treatment --group-b Control \
  --power 0.8 --alpha 0.05

Workspace isolation (CRITICAL)

The input data folder for any analysis must remain untouched so re-runs are reproducible. Scripts that write intermediate files (R drivers, prepared CSVs, comparison tables) must write to /tmp/ or to a --workdir you pass in. Both R-based scripts in this skill refuse to run if --workdir resolves to the input CSV's parent directory (or any ancestor of it).

# OK
--workdir /tmp/spline_run
 
# Refused:
--workdir <path-equal-to-or-containing-the-input-csv>/...

CRITICAL — Read before writing any code

1. Clinical trial AE analysis (regression, chi-square, ANY severity test): Use the bundled script (or the clinical_trial_ae_severity_test ToolUniverse tool which wraps it):

```bash tu run clinical_trial_ae_severity_test '{"dm_file":"DM.csv","ae_file":"AE.csv","test":"chi-square","group_col":"TRTGRP"}'

# Or directly: python skills/tooluniverse-statistical-modeling/scripts/prepare_ae_cohort.py \ --dm DM.csv --ae AE.csv --test chi-square --group TRTGRP \ --subgroup "expect_interact=Yes" # optional `` The script/tool handles: encoding='latin1' for SDTM CSVs, max(AESEV)` per subject across ALL AEs (no AEPT filtering), inner join with DM, optional subgroup filter, optional ordinal-logistic with covariates.

Why no AEPT filter — AESEV is a protocol-defined severity scale on the AE table. Filtering AE by AEPT (e.g. keeping only AEPT == "COVID-19") drops subjects whose worst severity was recorded under a different AEPT label, drastically changes the contingency table, and can flip the test result. The phrase "COVID-19 severity" describes the OUTCOME, NOT a filter criterion.

• ❌ WRONG: ae[ae['AEPT'].str.contains('COVID-19')].groupby('USUBJID')['AESEV'].max() — filters to COVID-19 events
• ✅ RIGHT: ae.groupby('USUBJID')['AESEV'].max() — uses ALL AE records

1. Expression ANOVA / fold change with multi-feature data (gene × sample matrix):

For "the F-statistic" or "a fold change" as a single value, run per-gene then summarize — NEVER pool expr.values.ravel() across all genes.

• For F-statistic: derive a per-sample quantity (like DESeq2 LFC of each gene between two cell types, then ANOVA on those LFCs across groups) OR run on a single target gene.
• For median/mean log2 fold change between two groups: run DESeq2 with design=~<group>, extract per-gene log2FoldChange (with shrinkage if the pipeline uses it), then take median/mean across genes.

❌ WRONG (aggregate): log2(sum_counts_groupA / sum_counts_groupB) per sample then summarize — gives ratio of totals, dominated by high-expression genes. ✅ RIGHT (per-gene): DESeq2 → results_df['log2FoldChange'].median().

Sanity heuristics: F > 50 for biological ANOVA across a few groups means you aggregated (typical biological F is 0.5–10). |median LFC| > 2 between similar groups means you aggregated (typical |median| < 1).

Use the bundled script: python skills/tooluniverse-statistical-modeling/scripts/expression_anova.py (or the expression_anova_per_gene ToolUniverse tool).

1. Spline models — R splines::ns(x, df=K) ≠ Python patsy.dmatrix("cr(x, df=K)"). They produce different design matrices because of internal-knot placement, boundary-knot placement, and basis orthogonalization. For ANY question that references R syntax like lm(y ~ ns(x, df = 4)), run R via Rscript. Use the bundled wrapper:

``bash python skills/tooluniverse-statistical-modeling/scripts/r_natural_spline_regression.py \ --csv data.csv --y-col Y --x-col X --df 4 --workdir /tmp/spline_run ``

For "frequency of strain X" co-culture models, include pure focal strain (freq=1) but exclude non-focal pure strain (freq=0).

1. CSV encoding: Clinical trial CSVs often need encoding='latin1'.

1. Pearson correlation between count-like and length-like variables: when one variable

spans orders of magnitude (raw read counts, TPM, gene length, transcript abundance), raw Pearson r is often near 0 even when log-transformed r is moderate. ALWAYS compute and explicitly report ALL FOUR variants in your final answer body: r(x, y), r(log10(x+1), y), r(x, log10(y+1)), r(log10(x+1), log10(y+1)). List as a table; mark one as your primary pick. The published answer can be ANY of the four, and the question text rarely disambiguates which transform combination was used.

``` ## Primary answer: r = X.XXX (transform: <name>)

## Sensitivity (all 4 transform combinations):

• r(x, y) = ...
• r(log10(x+1), y) = ...
• r(x, log10(y+1)) = ...
• r(log10(x+1), log10(y+1)) = ...

```

Background — for any single transform variant:

``python import numpy as np from scipy.stats import pearsonr r_raw, _ = pearsonr(x, y) r_log, _ = pearsonr(np.log10(x + 1), y) print(f"r_raw={r_raw:.4f} r_log10={r_log:.4f}") ``

Defaults:

• Question says "log-transformed" / "log expression" → report r_log10
• Question doesn't specify but the variable is gene expression / RNA count → also report r_log10 as the canonical answer (most published correlations between gene length and expression are log-scale)
• When |r_raw| < 0.1 AND |r_log10| > 0.2, prefer r_log10

❌ WRONG: report only r_raw ≈ 0.05 when log is 0.35 ✅ RIGHT: "r_raw = 0.05; r_log10 = 0.35 (canonical for log-distributed expression)"

COMPUTE, DON'T DESCRIBE

Write and run Python code (via Bash) for every statistical analysis. Never describe what you "would do" — do it. Use pandas for data wrangling, statsmodels for regression, scipy for tests, and matplotlib for plots. Execute the code and report actual numbers (β, p-value, CI, N).

LOOK UP, DON'T GUESS

When uncertain about any scientific fact, SEARCH databases first rather than reasoning from memory.

Features

• Linear Regression - OLS for continuous outcomes with diagnostic tests
• Logistic Regression - Binary, ordinal, and multinomial models with odds ratios
• Survival Analysis - Cox proportional hazards and Kaplan-Meier curves
• Mixed-Effects Models - LMM/GLMM for hierarchical/repeated measures data
• ANOVA - One-way/two-way ANOVA, per-feature ANOVA for omics data
• Model Diagnostics - Assumption checking, fit statistics, residual analysis
• Statistical Tests - t-tests, chi-square, Mann-Whitney, Kruskal-Wallis, etc.

When to Use

Apply this skill when user asks:

• "What is the odds ratio of X associated with Y?"
• "What is the hazard ratio for treatment?"
• "Fit a linear regression of Y on X1, X2, X3"
• "Perform ordinal logistic regression for severity outcome"
• "What is the Kaplan-Meier survival estimate at time T?"
• "What is the percentage reduction in odds ratio after adjusting for confounders?"
• "Run a mixed-effects model with random intercepts"
• "Compute the interaction term between A and B"
• "What is the F-statistic from ANOVA comparing groups?"
• "Test if gene/miRNA expression differs across cell types"

Model Selection Decision Tree