Statistical modeling — linear/logistic/ordinal/Poisson regression, ANOVA, Kruskal-Wallis, chi-square, Mann-Whitney, Cox survival, spline fits (R `ns()`), odds ratios, Cohen's d, F-statistic, p-value computation. Specializes in clinical-trial AE analysis (SDTM DM/AE), severity ordinal regression, and

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. Script When to use it r natural spline regression.py ANY question that mentions R syntax lm(y ~ ns(x, df = K)) , "natural spline", or asks for spline R²/F/peak prediction CIs. Always shells out to Rscript so splines::ns() matches. spline model compare.py "Best fitting model among quadratic, cubic and natural spline" / "max colony area at the optimal x". Fits all three in R, ranks by adj R²/AIC/BIC, and reports the BEST model's peak (x , y ) with 95% CI. logistic regression or.py Binary or ordinal logistic regression where the answer is an OR (or OR + 95% CI). Handles label encoding, explicit Placebo=0/BCG=1 maps, AND interaction terms ( interaction A:B creates A B = A B ). Prints OR + CI for every coefficient and a SCALARS block for the requested coef name . power analysis.py "Minimum sample size per group", "TTestIndPower", "given Cohen's d, what N for power=0.8". Computes pooled SD Cohen's d from a CSV (or accepts effect size ), then TTestIndPower.solve power . expression anova.py Per gene ANOVA / median LFC across cell types or sample groups (NOT pooled across genes — see warnings below). prepare ae cohort.py Clinical trial AE severity tests (chi square / ordinal) on SDTM DM/AE files ( encoding='latin1' , max(AESEV) per subject across ALL AEs — no AEPT filter). stat tests.py Stdlib only chi square goodness of fit, Fisher's exact, simple OLS. Use when scipy/statsmodels aren't available. Concrete invocations Natural spline regression (R^2, overall F test p, peak Y + 95% CI): Quadratic vs cubic vs natural spline comparison + best model peak: 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): Two sample power analysis from a pilot CSV: 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). 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): 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 2. 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). 3. 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: For "frequency of strain X" co culture models, include pure focal strain (freq=1) but exclude non focal pure strain (freq=0). 4. CSV encoding : Clinical trial CSVs often need encoding='latin1' . 5. 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. Background — for any single transform variant: 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 Workflow Phase 0: Data Validation Goal : Load data, identify variable types, check for missing values. CRITICAL: Identify the Outcome Variable First Before any analysis, verify what you're actually predicting: 1. Read the full question Look for "predict [outcome]", "model [outcome]", or "dependent variable" 2. Examine available columns List all columns in the dataset 3. Match question to data Find the column that matches the described outcome 4. Verify outcome exists Don't create outcome variables from predictors Common mistake : Question mentions "obesity" Assumed outcome = BMI = 30 (circular logic with BMI predictor). Always check data columns first: print(df.columns.tolist()) Phase 1: Model Fitting Goal : Fit appropriate model based on outcome type. Use the decision tree above to select model type, then refer to the appropriate reference file for detailed code: Linear Regression : references/linear models.md Logistic Regression (binary): references/logistic regression.md Ordinal Logistic : references/ordinal logistic.md Cox Proportional Hazards : references/cox regression.md ANOVA / Statistical Tests : anova and tests.md Quick reference for key models : Phase 1b: ANOVA for Multi Feature Data When data has multiple features (genes, miRNAs, metabolites), use per feature ANOVA (not aggregate). This is the most common pattern in genomics. See anova and tests.md for the full decision tree, both methods, and worked examples. Default for gene expression data : Per feature ANOVA (Method B). Phase 2: Model Diagnostics Goal : Check model assumptions and fit quality. Key diagnostics by model type: OLS : Shapiro Wilk (normality), Breusch Pagan (heteroscedasticity), VIF (multicollinearity) Cox : Proportional hazards test via cph.check assumptions() Logistic : Hosmer Lemeshow, ROC/AUC See references/troubleshooting.md for diagnostic code and common issues. Phase 3: Interpretation Goal : Generate publication quality summary. For every result, report: effect size (OR/HR/coefficient), 95% CI, p value, and model fit statistic. See common patterns summary.md for common question answer patterns. Common Patterns Pattern Question Type Key Steps 1 Odds ratio from ordinal regression Fit OrderedModel, exp(coef) 2 Percentage reduction in OR Compare crude vs adjusted model 3 Interaction effects Fit A B , extract A:B coef 4 Hazard ratio Cox PH model, exp(coef) 5 Multi feature ANOVA Per feature F stats (not aggregate) See common patterns summary.md for solution code for each pattern. See references/common patterns.md for 15+ detailed question patterns. Statsmodels vs Scik

tooluniverse-statistical-modeling