“A two‑group t‑test is rarely enough. Real MS experiments have batches, paired samples, and repeated measurements. Ignoring these structures inflates false positives.”
Real studies are messy: the same patient sampled before and after treatment, cell lines processed in batches, age and sex riding along with the condition you actually care about. Treat these correlated, structured measurements as if they were independent and you inflate your confidence and manufacture false positives. This chapter models the structure instead of ignoring it — often the difference between a result that replicates and one that does not.
Treating repeated measures as independent samples. Paired or longitudinal data analyzed as if independent inflates your effective n and your false-positive rate. Block on the subject.
By the end of this chapter, you will be able to:
limma design matrix for MS feature tableslme4 for repeated‑measures data with random subject effects (time courses)variancePartition to identify dominant factorsMS‑based omics experiments rarely have a simple two‑group design. Real data often include:
| Structure | Example in MS | Ignoring It Causes |
|---|---|---|
| Batch effects | Samples run on different days or columns | False positives (batch confounded with treatment) |
| Paired design | Same patient before/after drug | Reduced power (between‑subject variance not removed) |
| Repeated measures | Time course from same subjects | Inflated significance (pseudoreplication) |
| Confounders | Age, sex, BMI in clinical cohorts | Biased estimates of treatment effect |
| Factorial design | Treatment × time interaction | Missed interaction effects (e.g., treatment only works at late time) |
The principle: Your statistical model must match your experimental design. For MS data, the typical workflow after preprocessing (e.g., with xcms or DEP) results in a feature table (rows = metabolites/peptides, columns = samples). This chapter assumes you already have such a table.
When the same subject is measured multiple times, observations are correlated. The intraclass correlation coefficient (ICC) quantifies this:
\rho = \frac{\sigma^2_b}{\sigma^2_b + \sigma^2_e}
where \sigma^2_b is the between-subject variance and \sigma^2_e is the within-subject (residual) variance. The effective sample size for a study with T repeated measures per subject is:
n_{\text{eff}} = \frac{nT}{1 + (T-1)\rho}
At \rho = 0.5 and T = 4, only 57 % of your nominal sample contributes to inference. This is not a small correction — it means that ignoring repeated measures can inflate your effective sample size by nearly a factor of 2, producing falsely significant results.
The covariance structure models how observations are correlated. The choice affects standard errors but not (for linear models) the fixed-effect estimates:
| Structure | Formula | When to Use |
|---|---|---|
| Compound symmetry | All pairs equally correlated | Short time series, balanced designs |
| Random intercept + slope | Correlation varies with time distance | Growth trajectories (default choice) |
| AR(1) | Correlation decays exponentially with time separation | Equally spaced time points |
| Unstructured | Each pair has its own correlation | Small T, no prior structure assumed |
When in doubt, use random intercept + slope. It captures the most important correlation pattern — that subjects differ in both baseline level and rate of change — with only 2 extra parameters.
| Dataset | Source | Best for |
|---|---|---|
msdata::metabolomics() |
msdata |
Real LC‑MS metabolomics (8 samples) |
| Simulated (included) | this chapter | Learning concepts with known truth |
sacurine |
ropls |
Covariate adjustment (sex, age) |
We will start with a simulated feature table so you can see exactly how each model recovers true biological signals.
BiocManager::install(c(
"limma", # Linear models for MS data
"lme4", # Mixed‑effects models
"variancePartition", # Variance decomposition
"broom.mixed", # Tidy mixed model summaries
"msdata", # Example MS data
"ggplot2",
"dplyr"
))# Load libraries
library(limma)
library(lme4)
library(variancePartition)
library(ggplot2)
library(dplyr)
library(tidyr)
library(broom.mixed)
library(msdata) # for real MS files
# Set seed for reproducibility
set.seed(42)Before fitting, identify your experimental structure and choose the corresponding formula.
model_reference <- data.frame(
Design = c(
"Two‑group (simple)",
"Covariate‑adjusted",
"Paired / Blocked",
"Two‑way factorial",
"Time course (fixed effects)",
"Repeated measures (mixed)"
),
Formula = c(
"~ condition",
"~ condition + batch + sex",
"~ condition + subject_id",
"~ condition * time",
"~ condition + time",
"~ condition + (1 | subject)"
),
Function = c(
"lmFit() + eBayes()",
"lmFit() + eBayes()",
"lmFit() + eBayes()",
"lmFit() + eBayes()",
"lmFit() + eBayes()",
"lmer()"
)
)
knitr::kable(model_reference, caption = "Model formulas for common MS experimental designs")Usage: Match your design to a row, adapt the formula to your column names.
We’ll create a feature table with 200 metabolites and 18 samples.
True effects: treatment upregulates first 30 metabolites; batch 2 has a systematic shift; sex has a small effect.
n_metabolites <- 200
n_samples <- 18
# Experimental factors
condition <- factor(rep(c("Control", "Treatment"), each = 9))
batch <- factor(rep(c("B1", "B2", "B3"), times = 6))
sex <- factor(sample(c("M", "F"), n_samples, replace = TRUE))
# Feature table: rows = metabolites, columns = samples
intensity_mat <- matrix(rnorm(n_metabolites * n_samples, mean = 10, sd = 1),
nrow = n_metabolites, ncol = n_samples)
# Add batch effect (batch 2 is +1.5)
intensity_mat[, batch == "B2"] <- intensity_mat[, batch == "B2"] + 1.5
# Add sex effect (males +0.2)
intensity_mat[, sex == "M"] <- intensity_mat[, sex == "M"] + 0.2
# Add treatment effect for first 30 metabolites
intensity_mat[1:30, condition == "Treatment"] <-
intensity_mat[1:30, condition == "Treatment"] + 2
rownames(intensity_mat) <- paste0("Metab_", 1:n_metabolites)
colnames(intensity_mat) <- paste0("Sample_", 1:n_samples)
# Metadata (samples as rows)
metadata <- data.frame(
sample = colnames(intensity_mat),
condition = condition,
batch = batch,
sex = sex
)
cat(sprintf("Simulated MS feature table: %d metabolites × %d samples\n",
n_metabolites, n_samples))
head(metadata)Note: In real MS data, you would replace
intensity_matwith the output fromfeatureValues()fromxcmsor a normalized expression matrix fromDEP.
limma (Batch, Sex, etc.)If batch differs between control and treatment groups, an unadjusted model will find batch differences masquerading as treatment effects.
design_adj <- model.matrix(~ condition + batch + sex, data = metadata)
colnames(design_adj)fit_adj <- lmFit(intensity_mat, design_adj)
fit_adj <- eBayes(fit_adj)
# Extract treatment effect (coefficient "conditionTreatment")
results_adj <- topTable(fit_adj, coef = "conditionTreatment",
number = Inf, sort.by = "P")
head(results_adj[, c("logFC", "AveExpr", "P.Value", "adj.P.Val")])design_unadj <- model.matrix(~ condition, data = metadata)
fit_unadj <- lmFit(intensity_mat, design_unadj) |> eBayes()
results_unadj <- topTable(fit_unadj, coef = "conditionTreatment", number = Inf)
comparison_df <- data.frame(
Metabolite = rownames(results_adj)[1:30],
Adjusted_P = results_adj$P.Value[1:30],
Unadjusted_P = results_unadj$P.Value[1:30]
)
comparison_df |>
mutate(ratio = Unadjusted_P / Adjusted_P) |>
head(10) |>
knitr::kable(digits = 4)Expected finding: Adjusted P‑values are often lower (more significant) because batch variance is removed. In MS data, failing to adjust for batch can completely reverse conclusions.
n_patients <- 9
patient_id <- factor(rep(paste0("Patient_", 1:n_patients), times = 2))
condition_paired <- factor(rep(c("Before", "After"), each = n_patients))
# Feature table: each patient contributes two samples
mat_paired <- matrix(nrow = n_metabolites, ncol = 18)
for (i in 1:n_metabolites) {
baseline <- rnorm(n_patients, mean = 10, sd = 1)
effect <- ifelse(i <= 30, 1, 0) # true after effect for first 30
mat_paired[i, ] <- c(baseline, baseline + effect + rnorm(n_patients, 0, 0.5))
}
colnames(mat_paired) <- paste0(patient_id, "_", condition_paired)
rownames(mat_paired) <- paste0("Metab_", 1:n_metabolites)
paired_metadata <- data.frame(
sample = colnames(mat_paired),
patient = patient_id,
condition = condition_paired
)# Include patient ID to absorb subject‑to‑subject variation
design_paired <- model.matrix(~ condition + patient, data = paired_metadata)
fit_paired <- lmFit(mat_paired, design_paired) |> eBayes()
results_paired <- topTable(fit_paired, coef = "conditionAfter",
number = Inf, sort.by = "P")
cat("Significant metabolites (adj.P.Val < 0.05):",
sum(results_paired$adj.P.Val < 0.05), "\n")
cat("True positives detected:",
sum(results_paired$adj.P.Val[1:30] < 0.05), "/ 30\n")Key insight: Including patient as a fixed effect removes between‑patient variability, increasing power to detect the treatment effect. This is the linear model equivalent of a paired t‑test.
n_subjects <- 20
n_timepoints <- 4
long_df <- expand.grid(
subject = paste0("S", 1:n_subjects),
time = 1:n_timepoints,
condition = c("Control", "Treatment")
) |>
mutate(
subject_effect = rep(rnorm(n_subjects, 0, 0.5), each = n_timepoints * 2),
treatment_effect = ifelse(condition == "Treatment", 0.3 * time, 0),
abundance = 10 + subject_effect + treatment_effect + rnorm(n(), 0, 0.3)
)
head(long_df)fit_lmer <- lmer(abundance ~ condition * time + (1 | subject),
data = long_df)
tidy(fit_lmer, effects = "fixed") |> knitr::kable(digits = 3)Coefficient interpretation for MS time courses:
| Coefficient | Meaning |
|---|---|
| (Intercept) | Baseline abundance for Control at time 0 |
| conditionTreatment | Difference between Treatment and Control at time 0 |
| time | Change over time in Control group (slope) |
| conditionTreatment:time | Interaction – does treatment change the slope over time? |
For real MS data, you can loop over metabolites:
# Function to extract interaction p‑value
get_interaction_p <- function(metabolite_data) {
# metabolite_data must have columns: abundance, condition, time, subject
fit <- lmer(abundance ~ condition * time + (1 | subject),
data = metabolite_data)
tidy(fit, effects = "fixed") |>
filter(term == "conditionTreatment:time") |>
pull(p.value)
}
# Example for one metabolite
p_val <- get_interaction_p(long_df)
cat("Interaction p‑value:", p_val, "\n")Note: For many metabolites, consider using
limmawithduplicateCorrelation(faster). Mixed models are more flexible but computationally heavier.
Understanding variance sources helps you decide whether to include factors in your model and whether your experiment is well‑designed.
# Prepare data: samples as rows, metabolites as columns
intensity_t <- t(intensity_mat)
colnames(intensity_t) <- paste0("Metab_", 1:n_metabolites)
# Formula for variance components
vp_form <- ~ condition + batch + sex
# Fit variance partition
vp <- fitExtractVarPartModel(intensity_t, vp_form, metadata)
# Plot across all metabolites
plotVarPart(vp, main = "Variance decomposition in MS feature table")Interpreting the plot for MS data:
| Component | If high (>30%) | Implication |
|---|---|---|
| condition | Good – biological signal is strong | Your model will detect differences |
| batch | Problematic – technical artifact dominates | Need batch correction (including batch in model) |
| sex | Interesting – sex‑specific effects exist | Include sex as covariate or stratify |
| Residual | >50% – high unexplained noise | Need more replicates or better MS measurement |
Interaction answers: Does the effect of treatment depend on time (or another factor)?
Example: A drug reduces metabolite abundance only after 24 hours.
time2 <- factor(rep(c("T0", "T24"), each = n_samples / 2))
condition2 <- condition
# Add interaction effect for first 30 metabolites
mat_interact <- intensity_mat
mat_interact[1:30, condition2 == "Treatment" & time2 == "T24"] <-
mat_interact[1:30, condition2 == "Treatment" & time2 == "T24"] + 1.5
interact_metadata <- data.frame(
sample = colnames(mat_interact),
condition = condition2,
time = time2
)design_interact <- model.matrix(~ condition * time, data = interact_metadata)
fit_interact <- lmFit(mat_interact, design_interact) |> eBayes()
# Interaction coefficient (conditionTreatment:timeT24)
results_interact <- topTable(fit_interact, coef = "conditionTreatment:timeT24",
number = Inf)
head(results_interact[, c("logFC", "P.Value", "adj.P.Val")])top_metab <- rownames(results_interact)[1]
plot_df <- data.frame(
intensity = mat_interact[top_metab, ],
condition = condition2,
time = time2
)
ggplot(plot_df, aes(x = time, y = intensity, color = condition, group = condition)) +
stat_summary(fun = mean, geom = "line", size = 1) +
stat_summary(fun = mean, geom = "point", size = 3) +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
labs(
title = paste("Interaction plot for", top_metab),
subtitle = "Non‑parallel lines indicate a significant interaction",
x = "Time point", y = "log2 intensity"
) +
theme_minimal()Interpretation for MS data: If the lines cross or diverge, the treatment effect changes over time – that’s the interaction you are testing.
flowchart TD
A[Start: Your MS feature table & metadata] --> B{Repeated measures on<br>same subjects?}
B -->|Yes – same subject measured multiple times| C{How many time points?}
C -->|2 time points| D[limma paired model<br>~ condition + subject]
C -->|3+ time points| E[lmer mixed model<br>~ condition * time + (1|subject)]
B -->|No – independent samples| F{Any known confounders<br>batch, sex, age?}
F -->|Yes| G[limma covariate‑adjusted<br>~ condition + batch + sex]
F -->|No| H{More than one factor?}
H -->|Yes – factorial| I[limma interaction model<br>~ factor1 * factor2]
H -->|No – simple| J[limma two‑group<br>~ condition]
D --> K[Extract coefficient of interest]
E --> K
G --> K
I --> K
J --> K
K --> L[Apply eBayes shrinkage]
L --> M[FDR correction (BH)]
M --> N[Export results table]Replace your_intensity_matrix with your actual feature table (samples as columns).
# Complete covariate modeling pipeline for MS data
library(limma)
library(ggplot2)
library(dplyr)
# 1. Prepare metadata (must match column order of intensity matrix)
metadata <- data.frame(
condition = factor(rep(c("Ctrl","Treat"), each = 9)),
batch = factor(rep(c("B1","B2","B3"), times = 6)),
sex = factor(sample(c("M","F"), 18, replace = TRUE))
)
# 2. Build design matrix (include all known confounders)
design <- model.matrix(~ condition + batch + sex, data = metadata)
# 3. Fit model
fit <- lmFit(your_intensity_matrix, design) |> eBayes()
# 4. Extract results for condition effect
results <- topTable(fit, coef = "conditionTreatment",
number = Inf, sort.by = "P")
# 5. Filter significant features (FDR < 0.05)
sig_results <- results |> filter(adj.P.Val < 0.05)
# 6. Export
write.csv(sig_results, "MS_differential_results_covariate_adjusted.csv")
# 7. Quick variance partition diagnostic (samples as rows)
library(variancePartition)
vp <- fitExtractVarPartModel(t(your_intensity_matrix),
~ condition + batch + sex,
metadata)
plotVarPart(vp)| Pitfall | Why It’s a Problem | Solution |
|---|---|---|
| Adding too many covariates | Overfitting, reduced degrees of freedom | Include only known confounders (batch, sex, age) – not all possible variables |
| Ignoring batch when it’s confounded | Batch effect masquerades as treatment effect | Always check batch association with condition (e.g., χ² test) |
| Using mixed models for <5 subjects | Random effects poorly estimated | Use paired design with fixed subject effect instead |
| Forgetting interactions | Main effect may not tell full story (e.g., treatment only works at late time) | Plot data; test interaction if biologically plausible |
| Treating time as categorical with many levels | Loss of power, multiple testing burden | Consider continuous time with polynomial or spline terms |
| Not checking model assumptions | Residuals may not be normal, heteroscedasticity | Use plot(fit) for limma; plot(fit_lmer) for mixed models |
| Matrix orientation confusion | lmFit expects rows = features, columns = samples |
Verify dim(your_intensity_matrix) – transpose if needed |
Using the simulated data from Step 2, add age as a continuous covariate (simulate random ages 20‑80). Fit a model with ~ condition + batch + sex + age. How do the logFC and p‑values for the treatment effect change compared to the model without age?
# Your code hereLoad the msdata::metabolomics() files, process them with xcms (Chapter 7), and extract a feature table. If your experiment had paired samples (e.g., same subject before/after), run both a paired and an unpaired limma model. How many features are significant only in the paired model?
# Your code hereLoad the sacurine dataset from ropls. Use variancePartition to quantify variance explained by sample_type (urine), age, and gender. Which factor dominates?
# Your code hereUsing the interaction‑simulated data, test all metabolites for condition × time interaction. Report: - How many true positives (first 30 metabolites) are significant at FDR < 0.05? - How many false positives (metabolites 31‑200) are significant?
# Your code here| Output | Purpose |
|---|---|
| Covariate‑adjusted results | Differential abundance controlling for batch, sex, age |
| Paired model results | Increased power for matched designs (before/after) |
| Mixed model summary | Fixed (treatment, time) and random (subject) effects |
| Variance partition plot | Which factors (condition, batch, sex) drive MS variation? |
| Interaction plot | Visualise condition × time effects |
| Model decision table | Choose the right formula for your MS experiment design |
| Your MS Design | Use This Formula |
|---|---|
| 2 groups, no confounders | ~ condition |
| 2 groups + batch/sex/age | ~ condition + covariate |
| Matched pairs (before/after) | ~ condition + subject_id |
| Time course, 2 time points | ~ condition * time (interaction) |
| Time course, 3+ time points | lmer(value ~ condition * time + (1\|subject)) |
| Multiple variance sources unknown | variancePartition first, then model |
sessionInfo()