This section explains every mathematical concept used in the pipeline, in plain language with formulas. Written for a math/CS undergraduate who knows linear algebra and basic statistics but has no biology background.
Microbiome data is compositional: the numbers for all species in one sample always add up to 1 (or 100%). This is called the constant sum constraint.
Imagine you have three species: A, B, C with abundances:
Sample 1: A=0.33, B=0.33, C=0.33 (all equal)
Sample 2: A=0.50, B=0.25, C=0.25 (A doubled)In sample 2, only A changed biologically. But the numbers for B and C also changed — not because B and C changed, but because they had to, since everything must sum to 1. This creates spurious (fake) correlations between species that have nothing to do with each other biologically.
Standard Pearson correlation between two compositional variables is always biased:
\[\text{Cov}(x_i, x_j) = -\frac{1}{p^2} \sum_{k \neq i,j} \text{Var}(x_k) \neq 0 \text{ even when truly independent}\]
In plain terms: even if species A and species B have no biological relationship, their correlation will be non-zero just because they’re both measured as fractions of a whole.
CLR moves data from the simplex space (where everything sums to 1) into real Euclidean space (where standard statistics work normally).
For a sample with species abundances x₁, x₂, …, xₚ:
\[\text{clr}(x_i) = \log(x_i) - \frac{1}{p}\sum_{j=1}^{p} \log(x_j)\]
The second term is the geometric mean of all abundances in log space.
After CLR, the row means are always exactly zero:
\[\sum_{i=1}^{p} \text{clr}(x_i) = 0\]
This is the mathematical proof it’s working correctly — our step 3 code checks this.
We cannot take log(0). So before the CLR we add a tiny number (1e-6) to all values:
X_plus = X + 1e-6
clr = log(X_plus) - mean(log(X_plus), axis=1)Removes systematic differences between studies (batch effects) while preserving disease signal.
ComBat models each feature (taxon) value as:
\[Y_{ijg} = \alpha_g + X\beta_g + \gamma_{ig} + \delta_{ig}\epsilon_{ijg}\]
Where: - \(Y_{ijg}\) = observed CLR value for sample j, in batch i, for taxon g - \(\alpha_g\) = overall mean for taxon g - \(X\beta_g\) = biological signal to preserve (disease label) - \(\gamma_{ig}\) = additive batch effect to remove (study-level mean shift) - \(\delta_{ig}\) = multiplicative batch effect to remove (study-level variance scaling) - \(\epsilon_{ijg}\) = residual noise
Think of it as: observed = biology + lab artifacts + noise. ComBat estimates the lab artifact terms (γ and δ) using an empirical Bayes approach (borrows information across all taxa to get stable estimates even when n is small), then subtracts them out, leaving only biology + noise.
Instead of estimating γᵢ separately for each taxon (which is noisy), it assumes all the γᵢ values across taxa come from a common distribution, and uses that distribution to regularise each individual estimate. This is why it works even when n < 20 samples per batch.
XGBoost builds many simple decision trees one after another. Each tree tries to correct the mistakes of all the trees before it. The final prediction is a weighted sum of all trees.
Step 1 — Start with a constant prediction
For multi-class, initialise all predictions to equal probability: 1/K
for K classes.
Step 2 — Compute residuals (what we got wrong)
For each sample, compute the difference between what we predicted and
what was actually true. These are called
pseudo-residuals or gradients:
\[g_i = \frac{\partial L(y_i, \hat{y}_i)}{\partial \hat{y}_i}\]
where L is the loss function (log-loss for classification).
Step 3 — Fit a tree to the residuals
Build a decision tree that predicts the residuals, not the original
labels. The tree splits are chosen to maximise the
gain:
\[\text{Gain} = \frac{1}{2}\left[\frac{G_L^2}{H_L + \lambda} + \frac{G_R^2}{H_R + \lambda} - \frac{(G_L+G_R)^2}{H_L+H_R+\lambda}\right] - \gamma\]
Where: - G = sum of gradients in a leaf - H = sum of second-order
gradients (hessians) in a leaf - λ = L2 regularisation (our
reg_lambda=1.0) - γ = minimum gain required to make a split
(our min_child_weight)
Step 4 — Add the tree to the ensemble (with shrinkage)
\[\hat{y}_i^{(t)} = \hat{y}_i^{(t-1)} + \eta \cdot f_t(x_i)\]
Where η is the learning rate (our learning_rate=0.1). A
smaller learning rate means each tree contributes less — the model
learns more slowly but generalises better.
Step 5 — Repeat for t = 1 to n_estimators
reg_alpha = 0.1 → L1: pushes some feature weights to exactly 0
(like LASSO — automatic feature selection)
reg_lambda = 1.0 → L2: keeps all weights small
(like Ridge — prevents any one feature dominating)
min_child_weight=5 → minimum sum of sample weights in a leaf
(prevents overfitting on small subgroups)
subsample = 0.8 → use 80% of samples per tree (like bagging)
colsample_bytree=0.6 → use 60% of features per tree
(reduces feature correlation between trees)Our dataset has:
IBS: 305 samples (27.7%)
Healthy_Control: 353 samples (32.1%)
Crohn Disease: 265 samples (24.1%)
Ulcerative Colitis: 177 samples (16.1%) ← smallest classA naive model could achieve 32% accuracy by always predicting “Healthy” — without learning anything. This is the class imbalance problem.
\[w_k = \frac{N}{K \cdot n_k}\]
Where: - N = total samples (880 training) - K = number of classes (4) - nₖ = samples in class k
For Ulcerative Colitis (the smallest class): \[w_{UC} = \frac{880}{4 \times 142} = 1.549\]
For Healthy_Control (the largest class): \[w_{HC} = \frac{880}{4 \times 282} = 0.780\]
UC samples count 1.549/0.780 = 2× as much as Healthy samples during training. This forces the model to pay equal attention to all classes regardless of their frequency.
We need an honest estimate of how the model will perform on new patients it has never seen. If we evaluate on the same data we trained on, the model looks better than it really is (overfitting).
A simple random split might, by chance, put all UC samples in training and none in test. Stratified splitting guarantees each class appears in both sets at the same proportion:
Before split: UC = 16.1% of all 1100 samples
After split: UC = 16.1% of train AND 16.1% of testWe went further and stratified on disease × dataset jointly, because IBS only comes from one study — we needed MARS samples in both train and test.
Instead of one train/test split, we split the training data into K=5 equal parts (folds):
Fold 1: [val] [train] [train] [train] [train]
Fold 2: [train] [val] [train] [train] [train]
Fold 3: [train] [train] [val] [train] [train]
Fold 4: [train] [train] [train] [val] [train]
Fold 5: [train] [train] [train] [train] [val]Train 5 models, each evaluated on a different validation fold. Average the 5 accuracy scores. This gives a more reliable performance estimate than one split and uses all the data for training.
Critical rule: the held-out test set (220 samples) is never used during cross-validation. It is touched exactly once — at the very end, after all model decisions are finalised.
If 90% of patients are healthy and 10% have Crohn’s, a model that always predicts “healthy” gets 90% accuracy — without learning anything useful.
For each class, compute:
\[\text{Precision}_k = \frac{\text{True Positives}_k}{\text{True Positives}_k + \text{False Positives}_k}\]
\[\text{Recall}_k = \frac{\text{True Positives}_k}{\text{True Positives}_k + \text{False Negatives}_k}\]
\[F1_k = 2 \times \frac{\text{Precision}_k \times \text{Recall}_k}{\text{Precision}_k + \text{Recall}_k}\]
Macro F1 = simple average of F1 across all classes. Equal weight to each class regardless of size. Our result: 0.864.
A single number summarising a multi-class confusion matrix. Ranges from -1 (always wrong) to +1 (perfect) with 0 = no better than random:
\[\text{MCC} = \frac{\sum_{k}\sum_{l}\sum_{m} C_{kk}C_{lm} - C_{kl}C_{mk}}{\sqrt{\sum_k(\sum_l C_{kl})(\sum_{k'\neq k}\sum_l C_{k'l})} \times \sqrt{\sum_k(\sum_l C_{lk})(\sum_{k'\neq k}\sum_l C_{lk'})}}\]
In plain terms: MCC measures the correlation between predicted labels and true labels across all classes simultaneously. Our result: 0.822 — substantially better than random (0.0).
For each class, plot the True Positive Rate vs False Positive Rate at every possible probability threshold. AUC = area under this curve.
SHAP (SHapley Additive exPlanations) assigns each feature a contribution score to each prediction, based on game theory.
The SHAP value φᵢ for feature i is:
\[\phi_i = \sum_{S \subseteq F \setminus \{i\}} \frac{|S|!(|F|-|S|-1)!}{|F|!} \left[f(S \cup \{i\}) - f(S)\right]\]
Where: - F = set of all features - S = a subset of features not including i - f(S) = model prediction using only features in S
Think of each feature as a player in a game. The game’s prize is the model’s prediction. SHAP calculates each player’s fair share of the prize by averaging their contribution across all possible orderings of players joining the game.
A positive SHAP value for “Escherichia coli” in the Crohn’s Disease class means: “having high E. coli abundance pushed this patient’s Crohn’s risk score upward.”
For each feature, take the average of the absolute SHAP values across all test samples:
\[\text{Importance}_i = \frac{1}{n}\sum_{j=1}^{n} |\phi_i^{(j)}|\]
The higher this number, the more that feature consistently influenced predictions — in either direction.
The VAE trains by minimising two terms simultaneously:
\[\mathcal{L} = \underbrace{\mathbb{E}_{q(z|x)}[\log p(x|z)]}_{\text{Reconstruction loss}} - \underbrace{\beta \cdot D_{KL}(q(z|x) \| p(z))}_{\text{KL divergence}}\]
\[\mathcal{L}_{\text{recon}} = \frac{1}{n}\sum_{i=1}^{n}(x_i - \hat{x}_i)^2\]
Measures how well the decoder rebuilt the original CLR values. Minimising this forces the network to preserve information through the bottleneck.
\[D_{KL}(q \| p) = -\frac{1}{2}\sum_{j=1}^{d}\left(1 + \log\sigma_j^2 - \mu_j^2 - \sigma_j^2\right)\]
Where μ and σ are the encoder’s output for each latent dimension. This term pulls the latent space toward a standard normal distribution N(0,1). Without it, the encoder would memorise training data rather than learning a smooth, generalisable representation.
We cannot backpropagate through a random sampling operation. So instead of sampling Z ~ N(μ, σ²) directly, we write:
\[Z = \mu + \sigma \cdot \epsilon, \quad \epsilon \sim \mathcal{N}(0, I)\]
The randomness is in ε (which has no parameters), so gradients flow through μ and σ normally.
| Concept | Formula/Method | Pipeline step | Undergrad owns |
|---|---|---|---|
| Compositionality | Simplex constraint | Background | ✓ Explain in report |
| CLR transform | log(xᵢ) − mean(log(x)) | Step 3 | ✓ Implement + verify |
| Pseudocount | x + 1e-6 before log | Step 3 | ✓ Implement |
| ComBat | Empirical Bayes batch model | Step 4 | ✓ Run + validate PCA |
| Missingness flags | Binary indicator per study | Step 5 | ✓ Build mask logic |
| Stratified split | Joint disease×dataset key | Step 6 | ✓ Implement + verify |
| Class weights | wₖ = N/(K·nₖ) | Step 7 | ✓ Compute + apply |
| XGBoost gain | G²/(H+λ) split criterion | Step 7 | ✓ Tune hyperparameters |
| F1, MCC, AUC | Precision/recall curves | Step 7 | ✓ Compute + report |
| SHAP | Shapley value attribution | Step 7 | ✓ Visualise + interpret |
| VAE loss | Recon + β·KL | Step 3b | ✓ Implement + train |
| Reparameterisation | Z = μ + σ·ε | Step 3b | ✓ Implement |