R/estimate_dispersion.R
estimate_dispersion-methods.RdPort of HPCell's `se_estimate_dispersion()` (<https://github.com/MangiolaLaboratory/HPCell/blob/ed763c6f53cfedf63d3a3353ec5ffcfd09bfe0e6/R/differential_expression.R#L36>): tagwise dispersion when there are fewer than 1000 samples, otherwise trended dispersion on a random subset of up to 2000 samples. The design is a fixed-effect formula: edgeR has no random effects, so write the analogue yourself (`~ dex + cell` for `~ dex + (1 | cell)`, `~ dex * cell` for `~ dex + (dex | cell)`). [estimate_gene()] does not compute this; [estimate()] calls this helper once on the full object before iterating genes.
estimate_dispersion(
.data,
formula_abundance,
abundance = "counts",
dispersion_column = "dispersion",
dispersion_degrees_freedom_column = "dispersion_degrees_freedom"
)
# S4 method for class 'SummarizedExperiment'
estimate_dispersion(
.data,
formula_abundance,
abundance = "counts",
dispersion_column = "dispersion",
dispersion_degrees_freedom_column = "dispersion_degrees_freedom"
)
# S4 method for class 'RangedSummarizedExperiment'
estimate_dispersion(
.data,
formula_abundance,
abundance = "counts",
dispersion_column = "dispersion",
dispersion_degrees_freedom_column = "dispersion_degrees_freedom"
)A `SummarizedExperiment`
Fixed-effect model for the mean, used to build the edgeR design. This is not the mixed model you pass to [estimate_gene()]: random-effect terms (`|`) are rejected. Translate them yourself, e.g. `(1 | donor)` as `donor` and `(1 + treatment | donor)` as `donor + treatment:donor`.
Assay name (default `"counts"`).
Name of the `rowData` column for tagwise (or trended) dispersion \(\phi_g\) (default `"dispersion"`).
Name of the `rowData` column for the effective degrees of freedom \(d_{eff}\) behind that estimate (default `"dispersion_degrees_freedom"`).
`.data` with `rowData(.data)[[dispersion_column]]` and `rowData(.data)[[dispersion_degrees_freedom_column]]` filled. The raw edgeR object is in `metadata(.)$tidybulk$estimateDisp`.
A `SummarizedExperiment` object
A `SummarizedExperiment` object
Two columns are written, named by `dispersion_column` and `dispersion_degrees_freedom_column`. The first holds \(\phi_g\); the second holds the effective degrees of freedom behind that estimate. [estimate_gene()] turns the pair into a **prior** on the negative binomial shape, not a plug-in: \(\phi_g\) locates the prior and \(d_{eff}\) sets its tightness, so the gene-wise likelihood can still pull the posterior away from the edgeR value when the counts disagree.
The effective degrees of freedom are $$d_{eff} = (n - \mathrm{ncol}(design)) + \mathrm{prior.df}$$ The first term is the residual degrees of freedom of the fixed-effect design; `estimateDisp()` does not return it, so it is computed here. The second is the prior degrees of freedom that `estimateDisp()` does return, quantifying how far each gene-wise estimate is shrunk toward the mean-dispersion trend by empirical Bayes. `estimateTrendedDisp()` (the branch used above 1000 samples) performs no such shrinkage and returns no `prior.df`, so there \(d_{eff}\) is the residual degrees of freedom alone.
A grouping factor written as a fixed effect with close to one level per sample will exhaust the design and leave no residual degrees of freedom. That is an error rather than a column of `NA`s: pass a simpler `formula_abundance` for this step if it happens.
The link from degrees of freedom to a standard deviation runs through the chi-square distribution of the estimator. For \(s^2 \sim \sigma^2 \chi^2_d / d\), the log of the estimate has variance $$\mathrm{Var}(\log s^2) = \psi'(d/2)$$ with \(\psi'\) the trigamma function, so \(\mathrm{SD}(\log \hat\phi_g) \approx \sqrt{\psi'(d_{eff}/2)}\). The familiar \(\sqrt{2/d}\) is the large-\(d\) approximation to this and is 12 is the same empirical-Bayes variance model that limma and edgeR use to moderate gene-wise variances; it is exact for a scaled chi-square and is applied to the negative binomial dispersion by analogy.
The same \(d_{eff}\) supports a second, conjugate parameterisation, which [estimate_gene()] offers as `shape_prior = "gamma"`. Inverting the scaled inverse chi-square gives a gamma, so the precision \(1/\phi_g\) – the quantity brms calls `shape` – has prior $$\mathrm{Gamma}(d_{eff}/2,\ \mathrm{rate} = d_{eff}\phi_g/2)$$ with mean \(1/\phi_g\). Because \(\mathrm{Var}(\log X) = \psi'(a)\) for a gamma of shape \(a\), this reproduces \(\psi'(d_{eff}/2)\) exactly: the two routes are the same calculation, since \(\chi^2_d\) is itself \(\mathrm{Gamma}(d/2, \mathrm{scale} = 2)\). Its coefficient of variation is \(\sqrt{2/d_{eff}}\), which is where that familiar approximation comes from.
The two forms are not reparameterisations of one another, and differ in two ways worth knowing.
They centre different summaries of the shape on edgeR's estimate. brms gives the shape submodel a log link, so a Student-t on the intercept is symmetric in \(\log(\mathrm{shape})\) and places the *median* of the shape at \(1/\phi_g\). The gamma places its *mean* there, and a gamma's median lies below its mean. On the log scale the centres differ by $$\psi(d_{eff}/2) - \log(d_{eff}/2) \approx -1/d_{eff}$$ with \(\psi\) the digamma function: about -0.11, or 0.22 prior standard deviations, at \(d_{eff} = 9.8\), and vanishing as samples accumulate. Neither is wrong. Mean-centring is what the conjugate hierarchy dictates, since a scaled inverse chi-square on \(\phi_g\) implies \(E[1/\phi_g] = 1/\phi_g\); median-centring is what an additive offset on a log link implies.
They also differ in tail weight, which is why the Student-t is the default. Exponentiating a Student-t leaves a prior on the shape with polynomial tails and no finite mean, proper but very permissive, so a gene whose true dispersion is far from the shrunken edgeR estimate can still escape. The gamma decays exponentially and holds such a gene closer to the trend. The gamma is also left-skewed in \(\log(\mathrm{shape})\) (skewness \(\psi''(a) / \psi'(a)^{3/2}\), about -0.47 here) where the Student-t is symmetric.
Note that this standard deviation describes how precisely edgeR estimated its own \(\phi_g\). It does not account for the fact that the design here fits grouping factors as fixed rather than partially pooled, nor for a zero-inflated likelihood in which the `zi` component absorbs part of the overdispersion. Treating them as fixed is the conservative choice: it spends the full degrees of freedom that shrinkage would have given back, so \(d_{eff}\) understates rather than overstates the information behind \(\phi_g\).
Smyth GK (2004). Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. *Statistical Applications in Genetics and Molecular Biology* 3(1). [PDF](https://gksmyth.github.io/pubs/ebayes.pdf) — derives the \(\mathrm{Var}(\log s^2) = \psi'(d/2)\) result and the prior degrees of freedom used to moderate it.
McCarthy DJ, Chen Y, Smyth GK (2012). Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. *Nucleic Acids Research* 40(10):4288-4297. [PMC3378882](https://pmc.ncbi.nlm.nih.gov/articles/PMC3378882/) — Cox-Reid adjusted profile likelihood dispersion conditional on a design.
Robinson MD, McCarthy DJ, Smyth GK (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. *Bioinformatics* 26(1):139-140. [PMC2796818](https://pmc.ncbi.nlm.nih.gov/articles/PMC2796818/)
Phipson B, Lee S, Majewski IJ, Alexander WS, Smyth GK (2016). Robust hyperparameter estimation protects against hypervariable genes and improves power to detect differential expression. *Annals of Applied Statistics* 10(2):946-963. doi:10.1214/16-AOAS920 — `robust = TRUE`, which makes `prior.df` gene-specific and possibly infinite.
[`estimateDisp()` reference manual](https://rdrr.io/bioc/edgeR/man/estimateDisp.html), [edgeR on Bioconductor](https://bioconductor.org/packages/release/bioc/html/edgeR.html)
if (FALSE) { # \dontrun{
data("airway", package = "airway")
se <- estimate_dispersion(airway[1:150, ], ~ dex + cell)
} # }