---
title: "Differential expression and abundance"
subtitle: "Define the contrast, build the design, and understand variance moderation; limma as a worked expression-data method."
aliases:
- "/guides/statistical-background/09-linear-models-limma.html"
- "/guides/statistical-background/08-linear-models-limma.html"
guide-status: "Available"
---
**Before reading:** The relevant assay primer in [Omics data and measurement](../data-and-measurement/index.qmd), basic regression, and the distinction between an estimand and a model coefficient. The worked examples use gene expression; transferring the design and contrast ideas to another assay requires its own measurement and variance assumptions.
For a population-standardized contrast with covariate-varying group differences and missing assays, see [An end-to-end expression analysis](adjusted-expression-case-study.qmd). Its target and assumptions differ from a constant linear-model contrast; the examples are not a blanket ranking of procedures.
## Start with a feature-level contrast
For a feature $g$ measured on a stated scale, an unadjusted group contrast is $\Delta_g=E[Y_g\mid A=1]-E[Y_g\mid A=0]$. If adjustment is needed, define the conditional contrast or specify a common covariate distribution $Q$ over which to average it. The corresponding standardized contrast is $\int\{E[Y_g\mid A=1,Z=z]-E[Y_g\mid A=0,Z=z]\}\,dQ(z)$. These targets need not equal a coefficient in an arbitrary regression.
Decide which units are independent, which design supports the contrast, and which preprocessing created $Y_g$. Here the worked method is limma for log-expression, with a common linear contrast across the modeled covariates. Applying a similar workflow to proteins, metabolites, or relative microbial abundances requires assay-appropriate measurement and observation models. A difference on the log scale is not a difference in mean abundance on the original scale.
The [What is your estimator actually estimating?](../questions-and-design/estimands-and-misspecification.qmd) page explains interpretation under misspecification; the [Omics data and measurement](../data-and-measurement/index.qmd) overview compares the available assay primers and the assumptions that differ across them.
## One regression per gene
The standard analysis of an expression matrix fits the same linear model to every gene. For gene $g$ with expression $y_{gi}$ in sample $i$ (on the log scale),
$$
\begin{aligned}
y_{gi}&= \mathbf{x}_i^\top \boldsymbol\beta_g + \varepsilon_{gi},\\
\varepsilon_{gi}&\sim N(0, \sigma_g^2),
\end{aligned}
$$
where independent-person errors are assumed in this unweighted teaching model. The vector $\mathbf{x}_i$ records sample $i$'s group and covariates. It is the same row of the **design matrix** for every gene; the coefficient vector $\boldsymbol\beta_g$ and residual variance $\sigma_g^2$ are gene-specific. With 15,000 genes this is 15,000 regressions sharing one design, not 15,000 independent people.
Let $X$ have $N$ rows and $r$ linearly independent columns. Choose a fixed contrast vector $c$: for example, a vector that selects the patient coefficient, or subtracts the control mean coefficient from the patient mean coefficient. The model contrast and ordinary least-squares uncertainty are
$$
\begin{aligned}
\theta_g&=c^\top\boldsymbol\beta_g,\\
\hat\theta_g&=c^\top\hat{\boldsymbol\beta}_g,\\
\operatorname{Var}(\hat\theta_g\mid X)
&=\sigma_g^2\,c^\top(X^\top X)^{-1}c,\\
\widehat{\operatorname{SE}}_g
&=s_g\sqrt{c^\top(X^\top X)^{-1}c}.
\end{aligned}
$$
The design factor describes how much information this study provides about the chosen contrast. The residual SD $s_g$ describes feature-specific noise after fitting the mean. Their product is the contrast's SE; neither alone is the SE. Dividing $\hat\theta_g-\theta_{g,0}$ by that SE gives an ordinary t-statistic for the specified null value. Under this Gaussian, constant-variance, independent-error model its reference has $N-r$ residual degrees of freedom.
Everything you know about linear regression applies, per gene. What is new is that the number of samples is small (often 10 to 50), the number of genes is large, and the questions are asked across genes. Those three facts are what limma (Smyth 2004; Ritchie et al. 2015) is built around.
## The design matrix is the analysis
Students arriving from a regression course know `lm(y ~ group)`. In expression work you write the design matrix explicitly, once, and hand it to every gene.
```{r}
group <- factor(rep(c("control", "patient"), c(6, 10)))
batch <- factor(rep(c("b1", "b2"), 8))
sex <- factor(replace(rep("F", 16), c(4, 11), "M")) # two males, fixed so the page renders the same every time
X <- model.matrix(~ group + batch + sex)
head(X)
colnames(X)
```
Reading it: the intercept is the control, batch 1, female baseline; `grouppatient` is the patient-minus-control difference *at fixed batch and sex under the stated common-contrast model*; `batchb2` and `sexM` are adjustments. Adding a covariate to the design is what "adjusting for" means in practice, and the [What is your estimator actually estimating?](../questions-and-design/estimands-and-misspecification.qmd) page says what that buys: the coefficient on `group` is the group difference within levels of the covariates, which is the estimand you want only if the covariates are the right ones and enter correctly (a continuous composition score entered linearly is a linear adjustment; the group's research on cell-composition adjustment is about what that leaves behind).
Two design decisions that come up constantly:
- **Group-means parameterization.** `model.matrix(~ 0 + group)` gives one column per group and you then form contrasts (`patient - control`). Equivalent fits, often clearer when there are three or more groups.
- **Interactions.** `~ group * sex` lets the group difference differ by sex. Use interactions when the target permits group differences to vary across the covariate. Without them, a common-contrast assumption is being imposed; their value is not decided by available degrees of freedom alone.
For a continuous covariate $Z$, a model containing $A$, $Z$, and $AZ$ has conditional group contrast $\beta_{A,g}+\beta_{AZ,g}z$. The coefficient on $A$ alone is the contrast at $z=0$. Averaging over a specified covariate population gives $\beta_{A,g}+\beta_{AZ,g}E[Z]$. These coincide only in special cases. Always translate the proposed contrast into words before fitting it.
The inverse above also needs full column rank. If every patient is in batch 2 and every control is in batch 1, separate group and batch effects cannot be learned from this design. Software cannot recover the missing comparison by adding both columns. Even with full rank, limited overlap can make the contrast imprecise or dependent on extrapolation of the mean model.
```{r}
#| label: fig-expression-design
#| fig-height: 5
#| fig-cap: "The shared design has 16 specimens and four columns. Group and batch overlap, and both groups contain the sex categories; the rank check matters more than the visual appearance alone. Rows are not automatically independent people."
#| fig-alt: "A heat map shows intercept, patient, batch-2, and male indicator columns for the 16 samples."
stopifnot(qr(X)$rank == ncol(X))
image(seq_len(ncol(X)), seq_len(nrow(X)), t(X[nrow(X):1, ]),
col = c("white", "#2F80C8"), zlim = c(0, 1), xaxt = "n", yaxt = "n",
xlab = "Design column", ylab = "Specimen")
axis(1, at = 1:4, labels = c("Intercept", "Patient", "Batch 2", "Male"), cex.axis = 0.8)
axis(2, at = 1:16, labels = 16:1, las = 1, cex.axis = 0.65)
```
The design matrix specifies mean contrasts; the covariance model and variance calculation specify how dependence is handled. An ordinary independent-error fit with one row per specimen is inadequate when repeated specimens from one subject remain correlated. The repeated-measures section below discusses alternatives.
## What limma adds: borrowing strength across genes
With $N=16$ and four linearly independent design columns, fitting the mean uses four dimensions, leaving $16-4=12$ residual degrees of freedom. The estimate is $s_g^2=\mathrm{RSS}_g/12$. Across hypothetical repetitions, a gene's variance estimate can be too small or too large relative to its true variance. An accidentally small denominator can inflate a t-statistic, although a large statistic can also reflect a large numerator; the top-ranked genes need not have the smallest absolute variances.
limma's empirical Bayes step assumes the true variances across genes are drawn from a common prior, a scaled inverse chi-square with $d_0$ degrees of freedom and scale $s_0^2$, and forms a moderated variance estimate
$$
\tilde s_g^2 = \frac{d_0\, s_0^2 + d_g\, s_g^2}{d_0 + d_g}.
$$
Each gene's variance is pulled toward the common value $s_0^2$, by an amount that depends on how informative the prior is ($d_0$) relative to the data ($d_g$). In the ideal known-prior model, the **moderated t** uses $\tilde s_g$ in the denominator and has a marginal $t_{d_0+d_g}$ reference. This is calibration over the stated hierarchical variance distribution, not an exact conditional-on-every-fixed-variance theorem. In the implementation the hyperparameters are estimated, residual degrees of freedom can vary, and the total degrees of freedom are also constrained by available information. The formula is the simple constant-prior teaching case. $d_0$ and $s_0^2$ are estimated from the spread of the observed $s_g^2$ across genes (Smyth 2004, Section 6), so the shrinkage is data-driven: if gene variances really vary a lot, $d_0$ is small and little shrinkage happens.
The formula is a weighted average: the prior scale gets weight $d_0/(d_0+d_g)$ and the observed variance gets weight $d_g/(d_0+d_g)$. With $d_0=4$, $d_g=12$, and $s_0^2=1$, an observed variance of 2 becomes 1.75, while an observed variance of 0.1 becomes 0.325. Moderation can move a variance **up or down**. The added prior degrees of freedom represent information borrowed through the variance model, not four additional people.
The moderated statistic replaces $s_g$ by $\tilde s_g$ in the SE formula above; it leaves the fitted contrast in the numerator unchanged. This is why understanding the ordinary regression calculation first makes moderation easier to interpret.
Here is the effect with a known prior, which is a simplified version of what `eBayes()` estimates. This separate simulation uses only two group means and no batch/sex adjustment, so its design has rank 2 and residual degrees of freedom $16-2=14$, rather than the earlier design's 12.
```{r}
set.seed(1)
G <- 5000; n1 <- 8; n2 <- 8; d_g <- n1 + n2 - 2
sigma2 <- 1 / rgamma(G, shape = 2, rate = 2) # inverse-chi-square scale 1, df 4; mean 2
Y1 <- matrix(rnorm(G * n1, sd = sqrt(sigma2)), G) # control; no true differences
Y2 <- matrix(rnorm(G * n2, sd = sqrt(sigma2)), G) # patient
diff <- rowMeans(Y2) - rowMeans(Y1)
s2 <- (rowSums((Y1 - rowMeans(Y1))^2) + rowSums((Y2 - rowMeans(Y2))^2)) / d_g
se_factor <- sqrt(1 / n1 + 1 / n2)
t_ord <- diff / (sqrt(s2) * se_factor)
d0 <- 4; s0_2 <- 1 # true prior scale here; limma estimates it
s2_mod <- (d0 * s0_2 + d_g * s2) / (d0 + d_g)
t_mod <- diff / (sqrt(s2_mod) * se_factor)
# Assess noise relative to each gene's own known variance, not absolute s2.
top_ord <- order(-abs(t_ord))[1:20]
rbind(all_genes = quantile(s2 / sigma2), top_ordinary_t = quantile((s2 / sigma2)[top_ord]))
# The moderated reference here is marginal under this known-prior generator.
c(ordinary = mean(2 * pt(-abs(t_ord), d_g) < 0.05),
moderated = mean(2 * pt(-abs(t_mod), d0 + d_g) < 0.05))
```
Moderation borrows information about variances; it does not shrink the fitted group coefficient in this workflow. It can improve power and stability under suitable variance models, but cannot guarantee superiority for every feature or population. Compare the simulated rejection fractions with binomial MCSEs if assessing calibration across independently generated genes; that calculation is specific to this simulation.
```{r}
#| label: fig-variance-moderation
#| fig-cap: "Known-prior variance moderation toward scale 1. Smaller residual degrees of freedom give more shrinkage. The horizontal reference is a prior scale, not a known variance for every gene or the prior mean."
#| fig-alt: "Moderated variance curves lie between the observed variance and prior scale, with stronger shrinkage at fewer residual degrees of freedom."
observed_s2 <- seq(0.05, 5, length.out = 200)
residual_df <- c(4, 14, 60)
curves <- sapply(residual_df, function(d) (4 + d * observed_s2)/(4 + d))
matplot(observed_s2, curves, type = "l", lty = 1:3,
col = c("#2F80C8", "#D55E00", "#009E73"),
xlab = "Observed residual variance", ylab = "Moderated variance")
abline(0, 1, col = "gray65"); abline(h = 1, lty = 3, col = "gray65")
legend("topleft", paste("Residual df", residual_df), lty = 1:3,
col = c("#2F80C8", "#D55E00", "#009E73"), bty = "n")
```
### Run the actual log-expression workflow
This executable example uses **synthetic log2 expression** and the design above. It requires `limma` and its dependencies. Normalization is not learned in this generator. Feature and sample identifiers remain explicit.
```{r}
set.seed(81)
G_demo <- 600
beta <- rbind(rep(6, G_demo), c(rep(0.7, 40), rep(0, G_demo - 40)),
rep(0.2, G_demo), rep(-0.1, G_demo))
variance_demo <- 1 / rgamma(G_demo, shape = 3, rate = 2)
rownames(X) <- paste0("sample_", seq_len(nrow(X)))
logexpr <- t(beta) %*% t(X) +
matrix(rnorm(G_demo * nrow(X), sd = sqrt(variance_demo)), G_demo)
dimnames(logexpr) <- list(paste0("gene_", seq_len(G_demo)), rownames(X))
stopifnot(identical(colnames(logexpr), rownames(X)), qr(X)$rank == ncol(X))
fit <- limma::lmFit(logexpr, design = X)
fit <- limma::eBayes(fit, trend = TRUE, robust = TRUE)
limma::topTable(fit, coef = "grouppatient", number = 6, sort.by = "none",
adjust.method = "BH", confint = TRUE)
```
The output is a six-row display from the full 600-feature family. `number = 6` limits printing; it does not mean that only six hypotheses received the multiplicity adjustment. The model contrast is the constant patient–control difference at the same batch/sex levels. Report `logFC`, its scale, uncertainty, and the hypothesis family. The moderation and estimated hyperparameters are part of the inference model. BH still needs valid individual tests and the appropriate dependence conditions; its adjusted p-values do not make the displayed intervals simultaneous or selection-adjusted. See the [Multiple testing when the tests are not independent](multiple-testing.qmd) page.
### A separate raw-count example with voom
Do not feed raw counts into the unweighted log-expression model. This separate synthetic count example uses edgeR's filtering and normalization factors, then voom's observation-level precision weights. The resulting modeled scale and assumptions must be stated; it is not an absolute-abundance analysis or a validation of every real RNA-seq pipeline.
```{r}
#| label: fig-voom-weights
#| fig-cap: "The voom mean–variance trend in a separate synthetic count example. The fitted trend generates observation-level precision weights; it does not change the number of independent biological units."
#| fig-alt: "Voom plots estimated variability against average log-count size and a fitted smooth trend used to obtain precision weights."
set.seed(82)
expected_counts <- 2^(t(beta) %*% t(X)) *
outer(rep(1, G_demo), exp(seq(-0.4, 0.4, length.out = nrow(X))))
counts <- matrix(rnbinom(length(expected_counts), mu = as.vector(expected_counts), size = 10),
G_demo, dimnames = dimnames(logexpr))
dge <- edgeR::DGEList(counts)
keep <- edgeR::filterByExpr(dge, design = X)
dge <- dge[keep, , keep.lib.sizes = FALSE]
dge <- edgeR::calcNormFactors(dge)
v <- limma::voom(dge, design = X, plot = TRUE)
count_fit <- limma::eBayes(limma::lmFit(v, design = X))
c(tested_features = nrow(v), independent_people = ncol(v))
```
Filtering defines the tested feature set and needs the workflow's null-calibration justification. Keep the original counts and record filtering, normalization, design, and weighting decisions. `trend=TRUE` moderates variances with an intensity trend; voom models observation-level precision. They are not identical operations. The official [limma documentation](https://bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/intro.html) and [voom paper](https://doi.org/10.1186/gb-2014-15-2-r29) describe the implementations.
A voom weight is approximately an inverse variance for a log-expression observation: a noisier observation receives less weight in fitting its gene's mean model. It is not a count of additional biological replicates. Normalization factors address the modeled library scale, precision weights address variability on that scale, and variance moderation borrows information about residual variances. Understanding their separate roles helps diagnose which part of a workflow needs attention.
## What limma assumes, and when it breaks
The empirical Bayes step buys power with assumptions. Know them.
**One residual variance per gene, shared across groups.** The linear model has a single $\sigma_g^2$ for all samples, so patients and controls are assumed equally variable. When they are not, the moderated $t$ inherits the behavior of the pooled-variance $t$-test, which depends on which group is more variable and which is larger. The simple unweighted variance model has the same vulnerability as a pooled two-sample t-test. Welch addresses unequal variances for independent observations; it does not repair dependence from cohort-derived pseudo-values or network scores. See the variance issue on independently generated observations:
```{r}
set.seed(2)
Level <- function(n1, n2, sd1, sd2, R = 2000) {
p <- replicate(R, {
a <- rnorm(n1, sd = sd1); b <- rnorm(n2, sd = sd2)
c(pooled = t.test(a, b, var.equal = TRUE)$p.value, welch = t.test(a, b)$p.value)
})
rate <- rowMeans(p <= 0.05)
data.frame(n1, n2, sd1, sd2, method = rownames(p), rejection = rate,
MCSE = sqrt(rate * (1 - rate) / R), row.names = NULL)
}
rbind(`5 vs 15, small group more variable` = Level(5, 15, 2, 1),
`5 vs 15, small group less variable` = Level(5, 15, 1, 2),
`10 vs 10, unequal variances` = Level(10, 10, 2, 1))
```
Check the simulated rates and their Monte Carlo uncertainty. The pooled test is anti-conservative when the smaller group is more variable and can be conservative in the reverse configuration. Balance reduces this particular distortion but does not establish exact validity under unequal variances. Welch usually improves calibration here; very small samples still require scrutiny. The result concerns this synthetic configuration. Weighted limma models instead allow $\operatorname{Var}(\varepsilon_{gi})=\sigma_g^2/w_{gi}$ for supplied or estimated precision weights. `arrayWeights()` and `voomWithQualityWeights()` address particular sample/observation precision structures. The weights and correlation model still need justification for the assay and design; the unweighted example does not describe every weighting option.
**Genes share a variance distribution.** Borrowing variance information trades conditional bias against variability under the variance model. A feature far from the bulk can be affected differently from a typical feature. `trend=TRUE` allows an intensity-dependent prior scale, and `robust=TRUE` robustifies hyperparameter estimation against outlying variances; neither proves that every feature satisfies the model.
**Normality on the log scale, and independence across samples.** A transformation does not prove normality or control extreme tails. The exact small-sample reference uses its distributional assumptions; residual diagnostics and suitable large-sample arguments have their own scope. Dependence across samples requires a design-appropriate covariance calculation. Correlated genes do not by themselves create additional independent people.
## Repeated samples per subject
Longitudinal cohorts give several samples per subject. For a between-person contrast with positive within-person correlation, an independent-row analysis can understate uncertainty. For a within-person contrast, ignoring pairing can instead discard precision (the [Multivariate comparisons: permutation and distance-based tests](multivariate-comparisons.qmd) page shows the same failure in a different setting). limma has two tools.
**Block as a fixed effect.** Put subject in the design matrix. Useful for specified within-person contrasts under its mean/error assumptions, but it uses person-level parameters and renders a group main effect that is constant within person (patient versus control) inestimable, since subject is nested in group.
**`duplicateCorrelation()`.** Estimates a single within-subject correlation shared across genes, then fits each gene by generalized least squares with that correlation. Keeps subject-level contrasts estimable and borrows strength for the correlation the way `eBayes` does for the variance. It assumes the within-subject correlation is the same for every gene and every pair of samples from a subject, which should be assessed against the visit structure and biological heterogeneity rather than assumed suitable for every cohort.
```{r}
set.seed(83)
subject_id <- rep(seq_len(12), each = 2)
visit <- rep(0:1, 12)
group_r <- factor(rep(rep(c("control", "patient"), each = 6), each = 2),
levels = c("control", "patient"))
X_r <- model.matrix(~ group_r + visit)
U <- matrix(rnorm(200 * 12, sd = 0.6), 200, 12)
y_r <- 6 + U[, subject_id] + matrix(rnorm(200 * 24, sd = 0.6), 200, 24) +
outer(c(rep(0.7, 15), rep(0, 185)), as.numeric(group_r == "patient")) +
outer(rep(0.2, 200), visit)
corfit <- limma::duplicateCorrelation(y_r, design = X_r, block = subject_id)
repeated_fit <- limma::lmFit(y_r, design = X_r, block = subject_id,
correlation = corfit$consensus.correlation)
repeated_fit <- limma::eBayes(repeated_fit)
c(people = length(unique(subject_id)), specimens = ncol(y_r),
estimated_consensus_correlation = corfit$consensus.correlation)
```
In this generator the residual for a specimen is its person's shared random intercept plus its own noise. Both have variance $0.6^2=0.36$. Two visits from the same person therefore have covariance 0.36 and marginal residual variance $0.36+0.36=0.72$, giving within-person correlation 0.5. Visits from different people share no random intercept. This explains why compound symmetry matches this example and why 24 specimens still represent only 12 independent people. It does not establish that one consensus correlation fits every longitudinal assay.
Gene-specific mixed models, including `dream` in `variancePartition`, can model richer correlation structures under their own assumptions. Unequal visit counts also affect target weighting and selection, not only covariance. Choose the target and design before the fitting engine; see [Units of analysis, replication, and study design](../questions-and-design/units-and-study-design.qmd).
## Check the inference before interpreting the list
Check whether the design matrix is full rank and the desired contrast is estimable; perfect group–batch confounding cannot be removed by adding both columns. Examine residuals, mean–variance trends, influential samples, and whether repeated measurements were handled at the correct level. Report the estimated contrast and uncertainty on the stated scale, alongside the multiple-testing procedure and the family of features tested.
The weighted limma workflow can accommodate modeled precision differences; the simple unweighted example above has a narrower variance model. Consult the [limma introduction and user guide](https://www.bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/intro.html) for the implemented weighting and correlation options. Variance moderation improves variance estimation under its assumptions; it does not identify a confounded contrast or correct a scientifically inappropriate adjustment set.
## Where this connects
- The [What is your estimator actually estimating?](../questions-and-design/estimands-and-misspecification.qmd) page: the design matrix is where "adjusting for $Z$" happens, and the coefficient's target is the within-covariate difference only when the adjustment is correct and complete.
- The [Jackknife pseudo-values](jackknife-pseudo-values.qmd) page: comparing per-sample network edges between groups is a limma or $t$-test on the edges as if they were expression values. Check the dependence introduced by estimating those summaries from the shared cohort as well as the mean and variance models; ordinary independent-observation testing is not automatically justified.
- The [Orthogonality and cross-fitting](../semiparametric-methods/orthogonality-and-cross-fitting.qmd) page: a misspecified linear adjustment can change the coefficient's target. An orthogonal, cross-fitted construction is an alternative when its target, model, and nuisance-rate conditions fit the question; it is not an automatic repair for every design or measurement problem.
- The [Multiple testing when the tests are not independent](multiple-testing.qmd) page: the output of limma is 15,000 $p$-values, and what you do with them is the subject of that shared inference guide.
## Exercises
**Checkpoints are completion targets for this page. Additional practice below is optional.** Each checkpoint gives a completion standard, a hint, and a solution outline.
::: {#checkpoint-de-01 .checkpoint}
### Checkpoint · DE-01 — Write an estimable contrast {.unnumbered}
**Task type:** Construction and interpretation.
Use a group-means design with columns $(A,B,C,Z)$ for three groups and one continuous covariate. Write contrast vectors for B minus A and the average of B and C minus A. State the scale and common-slope assumption needed to interpret these contrasts as adjusted feature differences. Explain what must be checked before claiming the contrasts are estimable.
**Ready to continue when:** you can connect a design/contrast calculation to its population interpretation and check rank and confounding.
**Transfer:** If the group difference varies with $Z$, which additional design terms and reference values would you need?
<details class="checkpoint-help">
<summary>Hint</summary>
A contrast takes a weighted sum of the coefficients; a common covariate term cancels in a same-$Z$ difference.
</details>
<details class="checkpoint-help">
<summary>Solution outline — open after trying</summary>
In the stated column order, the vectors are $(-1,1,0,0)$ and $(-1,1/2,1/2,0)$. Interpretation is on the feature scale modeled, such as log-expression. Without interactions the design assumes a common covariate slope and a group difference constant in $Z$. Check the design rank and whether the contrast lies in its estimable space; perfect confounding or absent groups can prevent identification. More generally, a standardized group contrast can differ from a coefficient in a misspecified model.
</details>
:::
::: {#checkpoint-de-02 .checkpoint}
### Checkpoint · DE-02 — Explain what variance moderation changes {.unnumbered}
**Task type:** Derivation and interpretation.
In the known-prior teaching model, take $d_0=4$, $s_0^2=1$, $d_g=4$, and observed $s_g^2=0.2$. Compute the moderated variance and compare it with the result at $d_g=40$. State whether the fitted group coefficient is shrunk. Explain what is averaged over in the known-prior moderated-t calibration and what remains to check in an actual weighted, repeated-measures limma analysis.
**Ready to continue when:** you can separate the target contrast, variance borrowing, hierarchical calibration, and the independent-unit/weighting assumptions.
**Transfer:** Would changing a batch-confounded design become identifiable after stronger variance moderation? Why?
<details class="checkpoint-help">
<summary>Hint</summary>
Use the weighted variance formula; the hyperprior scale is not every gene’s true variance.
</details>
<details class="checkpoint-help">
<summary>Solution outline — open after trying</summary>
The moderated variance is $(4+4(0.2))/8=0.6$ at four residual df, and $(4+40(0.2))/44=12/44$ at forty df. More residual information reduces shrinkage toward the prior scale. This variance moderation does not shrink the fitted group coefficient in the displayed workflow. The ideal moderated-t reference integrates over the specified gene-variance distribution; it is not an exact conditional-on-every-fixed-variance guarantee. Estimated hyperparameters, mean specification, precision weights, within-person covariance, contrast estimability, marginal p-values, and the tested family still need checking.
</details>
:::
### Additional practice (optional) {.unnumbered}
These exercises extend the ideas through further derivations, simulations, or applications. Choose those useful for your work.
1. In the shrinkage simulation, add a true difference of 1 to the first 200 genes and compare the number of true positives at a Bonferroni threshold for the ordinary and moderated statistics. Then set `d0 <- 50` and `d0 <- 0.5` and explain what happens.
2. Extend the `Level()` simulation so that the two groups also differ in shape (one normal, one a $t_3$ scaled to the same SD). Does Welch still hold its level? At what sample size?
3. Simulate 20 subjects with 3 samples each and a within-subject correlation of 0.5, no group effect. Fit `lm(y ~ group)` ignoring subjects and record the level over 2000 simulated genes. Then fit with subject as a random effect (`lme4::lmer`) and compare.
## Reading
- Smyth, G. K. (2004). Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. *Statistical Applications in Genetics and Molecular Biology*, 3(1), Article 3. Sections 1 to 4 and 6.
- Ritchie, M. E., et al. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. *Nucleic Acids Research*, 43(7), e47.
- Law, C. W., Chen, Y., Shi, W., and Smyth, G. K. (2014). voom: precision weights unlock linear model analysis tools for RNA-seq read counts. *Genome Biology*, 15, R29.
- Law, C. W., et al. (2020). A guide to creating design matrices for gene expression experiments. *F1000Research*, 9, 1444. A worked design-matrix reference.
- Smyth, G. K., Michaud, J., and Scott, H. S. (2005). Use of within-array replicate spots for assessing differential expression in microarray experiments. *Bioinformatics*, 21(9), 2067 to 2075. The `duplicateCorrelation` paper.
<span id="continue-through-the-collection"></span>
