bio-data-visualization-statistical-annotation

What this skill does

## Version Compatibility

Reference examples tested with: ggpubr 0.6+, ggsignif 0.6+, rstatix 0.7+, statannotations 0.6+ (Python), seaborn 0.13+.

Before using code patterns, verify installed versions match. If versions differ:
- R: `packageVersion('<pkg>')` then `?function_name`
- Python: `pip show <package>` then `help(module.function)` to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

# Statistical Annotation

**"Add p-values to my plot"** -> Render pairwise group comparisons as brackets with the correct statistical test (parametric vs non-parametric, paired vs unpaired, independent vs nested), adjusted for multiple testing, with rendering of significance as either numerical p OR asterisks. The choices that matter: which test is appropriate for the data, what multiple-testing adjustment applies, and whether to show n.s. (non-significant) results.

- R: `ggpubr::stat_compare_means`, `ggsignif::geom_signif`, `rstatix::t_test`/`wilcox_test`
- Python: `statannotations.Annotator`, `scipy.stats` directly

## The Single Most Important Modern Insight -- The Test Must Match the Data

Tool defaults are NOT data-appropriate. `ggpubr::stat_compare_means(method='t.test')` uses Welch's two-sample t-test assuming approximate normality and unequal variances. This is wrong when:

1. **Data are non-normal and N is small (<30):** use Mann-Whitney U (`method='wilcox.test'`).
2. **Data are paired:** use paired t-test or paired Wilcoxon (`paired = TRUE`).
3. **Comparing >2 groups:** ANOVA / Kruskal-Wallis with post-hoc, not all-pairs t-test (multiple-testing penalty).
4. **Data are nested (cells within patients, replicates within samples):** linear mixed model, NOT pairwise test.

The bracket-and-asterisk visual is the same; the underlying statistics are not. Choose the test deliberately.

## Decision Tree for Test Selection

| Question | Recommended test | Function |
|----------|------------------|----------|
| 2 unpaired groups, normal, N≥30 | Welch t-test | `t.test()`, `stat_compare_means(method='t.test')` |
| 2 unpaired groups, non-normal or small N | Mann-Whitney U (Wilcoxon rank-sum) | `wilcox.test()`, `stat_compare_means(method='wilcox.test')` |
| 2 paired groups | Paired t-test OR Wilcoxon signed-rank | `paired = TRUE` |
| 3+ groups, normal | One-way ANOVA + Tukey HSD post-hoc | `aov()`, `TukeyHSD()` |
| 3+ groups, non-normal | Kruskal-Wallis + Dunn post-hoc | `kruskal.test()`, `dunn.test()` |
| Nested data (cells in patients) | Linear mixed model | `lme4::lmer()` |
| Time-course / repeated measures | Repeated-measures ANOVA OR LMM | `nlme::lme()` |
| Two-way factorial | Two-way ANOVA + interaction term | `aov(y ~ a*b)` |
| Survival / time-to-event | Log-rank, NOT t-test | `survdiff()` |
| Categorical outcome | Chi-square OR Fisher exact | `chisq.test()`, `fisher.test()` |

## Multiple Testing Adjustment

For pairwise comparisons among K groups, there are K(K-1)/2 unadjusted p-values. Without adjustment, family-wise error rate inflates rapidly:
- 3 groups: 3 comparisons; α_FW = 14% at nominal 5%
- 4 groups: 6 comparisons; α_FW = 26%
- 6 groups: 15 comparisons; α_FW = 54%

```r
# rstatix supports per-comparison adjustment
library(rstatix)
df %>%
    pairwise_wilcox_test(value ~ group, p.adjust.method = 'bonferroni') %>%
    add_xy_position(x = 'group')
```

`p.adjust.method` options:
- `'bonferroni'` — strictest; controls FWER
- `'holm'` — stepwise Bonferroni; uniformly more powerful than Bonferroni
- `'BH'` (Benjamini-Hochberg) — FDR; less strict than FWER; standard for genomics
- `'fdr'` — alias for BH

For figure annotations, **holm** is the modern default — controls FWER and is more powerful than bonferroni. For a small number of pre-planned comparisons (≤3), Bonferroni is fine.

## ggpubr -- Standard ggplot2 Workflow

**Goal:** Add per-comparison p-value brackets between groups on a distribution plot, using a test appropriate to data shape and adjusting for multiple comparisons.

**Approach:** Build the base plot with `ggboxplot()`; add `stat_compare_means()` with explicit `method`, `comparisons`, `p.adjust.method`, and `label` arguments; render as asterisks (`p.signif`) for terse display or numeric (`p.format`) for precise display.

```r
library(ggpubr)

# Default boxplot + p-value bracket(s)
ggboxplot(df, x = 'group', y = 'value', color = 'group',
          add = 'jitter', palette = 'npg') +
    stat_compare_means(method = 'wilcox.test',           # explicit; default is t-test
                       comparisons = list(c('Control', 'Treatment'),
                                          c('Control', 'Vehicle'),
                                          c('Treatment', 'Vehicle')),
                       label = 'p.signif',               # 'p.signif' for asterisks; 'p.format' for numeric
                       p.adjust.method = 'holm',
                       method.args = list(alternative = 'two.sided'))
```

For an overall test plus pairwise:

```r
# Overall + per-comparison
ggboxplot(df, x = 'group', y = 'value', color = 'group') +
    stat_compare_means(method = 'kruskal.test',          # overall test
                       label.y = 1.05 * max(df$value)) +
    stat_compare_means(comparisons = pairs,
                       method = 'wilcox.test',
                       p.adjust.method = 'holm',
                       label = 'p.signif')
```

## ggsignif -- Lighter Alternative

```r
library(ggsignif)
ggplot(df, aes(group, value, fill = group)) +
    geom_boxplot() +
    geom_signif(comparisons = list(c('Control', 'Treatment')),
                test = 'wilcox.test',
                map_signif_level = TRUE,                  # asterisks vs numeric p
                step_increase = 0.1) +
    scale_fill_manual(values = c('#0072B2', '#D55E00'))
```

`map_signif_level = TRUE` converts p-values to asterisks per Wasserstein-Lazar 2016 convention:
- `***` p < 0.001
- `**` p < 0.01
- `*` p < 0.05
- `ns` p ≥ 0.05

For literal p-values, set `FALSE`.

## statannotations (Python)

```python
import seaborn as sns
from statannotations.Annotator import Annotator

ax = sns.boxplot(x='group', y='value', data=df, palette=['#0072B2', '#D55E00', '#009E73'])

pairs = [('Control', 'Treatment'),
         ('Control', 'Vehicle'),
         ('Treatment', 'Vehicle')]

annotator = Annotator(ax, pairs, data=df, x='group', y='value')
annotator.configure(test='Mann-Whitney',                    # 't-test_ind', 't-test_paired', 'Wilcoxon', etc
                    comparisons_correction='holm',
                    text_format='star',                     # 'star', 'simple', 'full'
                    line_height=0.02,
                    text_offset=0.5)
annotator.apply_and_annotate()
```

## Per-Method Failure Modes

### Default t-test on non-normal data

**Trigger:** `stat_compare_means(method='t.test')` (default) on log-distributed expression.

**Mechanism:** t-test assumes approximate normality; non-normal data with small N inflates Type-I.

**Symptom:** Significant p where rank test gives p > 0.05.

**Fix:** Switch to `method='wilcox.test'` for non-normal or small-N data. Verify normality with Shapiro-Wilk if borderline.

### Pairwise tests without adjustment

**Trigger:** Multiple bracket annotations with raw p-values.

**Mechanism:** K(K-1)/2 comparisons inflate FWER without adjustment.

**Symptom:** All-pairwise significant at nominal 0.05; doesn't replicate.

**Fix:** `p.adjust.method = 'holm'` (or 'bonferroni' or 'BH'). Document choice.

### Paired data tested as independent

**Trigger:** Before/after measurements in same subjects, tested with unpaired t-test.

**Mechanism:** Ignores within-subject correlation; loses power.

**Symptom:** Non-significant p where paired test gives significant.

**Fix:** `paired = TRUE` (R) or `t-test_paired` (statannotations). Verify subjects are correctly matched.

### Nested data tested with pairwise t

**Trigger:** Hu