Worked Regularized NB Normalization Example

From a Count Model to an Adjusted Expression Value

A negative-binomial (NB) normalization model describes how a gene’s expected count changes with cell depth. Its fitted mean can be used in more than one way: one can rescale an observed count to a reference depth, or compare the observed count with the fitted mean using a residual. Those quantities answer different questions, so the model and its assumptions must be clear before interpreting either.

This worked example uses four hypothetical cells and two genes. It is intentionally small enough to calculate by hand. The depth factors are known inputs, not estimates from these four cells; the example illustrates the fitting and adjustment arithmetic rather than recommending a normalization strategy for a real data set.

Definition: For gene $g$ and cell $i$, let $Y_{gi}$ be the observed count, $d_i$ the positive depth factor, and $\mu_{gi}$ the expected count. The gene-wise NB depth model is $$ Y_{gi}\sim \operatorname{NB}(\mu_{gi},\phi_g),\qquad \log(\mu_{gi})=\alpha_g+\beta_g\log(d_i), $$ where $\operatorname{Var}(Y_{gi})=\mu_{gi}+\phi_g\mu_{gi}^{2}$. Here $\alpha_g$ is the log expected count at $d_i=1$, $\beta_g$ is the gene-specific depth coefficient, and $\phi_g\geq 0$ is the dispersion. The model assumes conditional independence across cells for a gene, with the specified mean-variance relationship and depth factors treated as fixed.

The factors in the example are $d=(1/4,1/2,2,4)$. Their product is $1$, so their geometric mean is $1$; equivalently, the four values of $x_i=\log(d_i)$ sum to zero. This centering makes $d=1$ a natural reference depth. It does not make the factors independent of biology or guarantee that they capture every technical effect.

Cell	Depth factor $d_i$	Gene A count	Gene B count
1	0.25	4	1
2	0.50	4	2
3	2	4	8
4	4	4	16

For this hypothetical exercise, fix both dispersions at $\phi_g=0.20$. Fixing dispersion isolates the depth-coefficient calculations; in practice, dispersion also needs to be estimated or otherwise specified, and uncertainty in it can affect inference. The table has an especially transparent pattern: Gene A is constant across cells, while Gene B’s counts are proportional to depth.

Fitting the Gene-Wise Coefficients

For each gene, the likelihood is the product of the NB probabilities for its four counts, with $\mu_{gi}=\exp(\alpha_g+\beta_g x_i)$. In general, one numerically maximizes this likelihood over $\alpha_g$ and $\beta_g$, holding the specified dispersion fixed. The fit is gene-wise: each gene has its own intercept and depth coefficient.

In this constructed table, each gene’s observed counts equal the means at a particular coefficient pair. Since the NB probability for an individual observation, as a function of its mean, is maximized at $\mu=y$, these pairs give the maximum likelihood fits. Gene A has $\hat{\alpha}_A=\log(4)$ and $\hat{\beta}_A=0$, so all fitted means equal 4. Gene B has $\hat{\alpha}_B=\log(4)$ and $\hat{\beta}_B=1$, so its fitted means are $4d_i=(1,2,8,16)$. The model reproduces the counts exactly in this deliberately simple example; that will not generally happen in real data.

Worked Example: Verifying the Maximum-Likelihood Fits

The NB log-mean score for one observation is $(y-\mu)/(1+\phi\mu)$. Thus, if a proposed coefficient pair gives $\mu_i=y_i$ in every cell, both the intercept and slope scores are zero. More strongly, each fitted mean maximizes that cell’s NB likelihood, so the pair is a global maximum for the joint likelihood.

For Gene A, $\alpha=\log(4)$ and $\beta=0$ give $\mu_i=\exp(\log(4))d_i^0=4$ in each cell, matching $(4,4,4,4)$. For Gene B, $\alpha=\log(4)$ and $\beta=1$ give $\mu_i=4d_i=(1,2,8,16)$, matching its observed counts. These are maximum-likelihood estimates for these specific data and fixed dispersion, not assumed fitting outputs.

A model-based standard error describes local uncertainty under the fitted NB model. For log mean $\eta=\log(\mu)$, the score is $(Y-\mu)/(1+\phi\mu)$. Its expected information is the variance of that score: $$ \mathcal I_{\eta\eta} =\frac{\operatorname{Var}(Y)}{(1+\phi\mu)^2} =\frac{\mu}{1+\phi\mu}. $$

For an intercept and slope, each cell contributes this weight times the outer product of $(1,x_i)$. After allowing for the estimated intercept, the information for the slope is $\sum_i w_i(x_i-\bar{x}_w)^2$, where $w_i=\mu_i/(1+\phi\mu_i)$ and $\bar{x}_w=\sum_iw_ix_i/\sum_iw_i$. The usual model-based standard error is the reciprocal square root of that information. This is a local, large-sample approximation; four cells are far too few to make it a reliable inferential analysis.

Worked Example: Calculating a Slope Standard Error

For Gene A, each fitted mean is 4, so each information weight is $4/(1+0.20\times4)=2.2222$. The log depths are $-1.3863,-0.6931,0.6931,1.3863$. Since the weights are equal and these log depths sum to zero, their weighted mean is zero. Therefore,

\mathcal I_{\beta_A\beta_A} =2.2222\left[(-1.3863)^2+(-0.6931)^2+(0.6931)^2+(1.3863)^2\right] \approx 10.68, \qquad SE(\hat{\beta}_A)\approx\frac{1}{\sqrt{10.68}}=0.3060.

The estimate is $\hat{\beta}_A=0$, with an approximate standard error of 0.3060 on the log-depth slope scale. For Gene B, the fitted-mean weights are $0.8333,1.4286,3.0769,3.8095$. Using their weighted log-depth variance gives slope information of approximately 8.05 and $SE(\hat{\beta}_B)\approx0.3524$. The weighted centering matters here: Gene B’s weights are unequal, so its weighted mean log depth is not zero.

Regularizing the Depth Coefficients

A gene with few counts may have a noisy depth estimate. Regularization pulls gene-specific estimates toward a shared target, trading some possible bias for reduced variability. The precise target and strength of shrinkage depend on the method. To make the arithmetic explicit, use a normal approximation to the gene-wise estimate and a normal prior centered at a common target $\beta_0=1$, with prior variance $\tau^2=0.25$. This is a stated illustrative shrinkage rule, not a claim that every NB normalization procedure uses this prior.

Definition: If $\hat{\beta}_g$ has approximate sampling variance $v_g=SE(\hat{\beta}_g)^2$, and the prior is $\beta_g\sim N(\beta_0,\tau^2)$, the normal-normal posterior mean is $$ \tilde{\beta}_g =\frac{\hat{\beta}_g/v_g+\beta_0/\tau^2}{1/v_g+1/\tau^2} =\beta_0+\frac{\tau^2}{\tau^2+v_g}(\hat{\beta}_g-\beta_0). $$

Worked Example: Shrinking the Two Depth Slopes

For Gene A, $v_A=0.3060^2\approx0.09364$. Its shrinkage weight on the difference from the target is $0.25/(0.25+0.09364)\approx0.7275$. Hence

\tilde{\beta}_A =1+0.7275(0-1) \approx0.2725.

Gene A’s slope is pulled from 0 toward the shared target 1. Gene B’s unregularized estimate already equals the target, so the same rule gives $\tilde{\beta}_B=1+0.6681(1-1)=1$, using $v_B\approx0.3524^2=0.1242$. This simple example illustrates that regularization does not necessarily change every gene’s coefficient.

For the calculations below, retain each gene’s fitted intercept $\hat{\alpha}_g=\log(4)$ and replace only its slope by the regularized slope. That defines the post-fit adjustment being illustrated. Other procedures may jointly refit the intercept after regularizing slopes; results should not be mixed between these procedures.

Adjusted Counts and Pearson Residuals

With the stated post-fit rule, the regularized fitted mean is $\tilde{\mu}_{gi}=4d_i^{\tilde{\beta}_g}$. An observed count can be rescaled to the reference depth by dividing by the depth multiplier $d_i^{\tilde{\beta}_g}$. This yields a count-equivalent adjusted observation, which need not be an integer. The model’s expected count at the reference depth is 4 for both genes because the intercept is held fixed. An adjusted observation is not the same thing as that expected value.

A Pearson residual instead expresses the difference between the observed count and the fitted mean in estimated standard-deviation units under the NB variance model:

r_{gi} =\frac{Y_{gi}-\tilde{\mu}_{gi}} {\sqrt{\tilde{\mu}_{gi}+\phi_g\tilde{\mu}_{gi}^{2}}}.

Worked Example: Adjusting a Gene A Count

For Gene A in the cell with $d=4$, the regularized slope is 0.2725. The depth multiplier, fitted mean, and adjusted observed count are

4^{0.2725}=\exp(0.2725\log 4)\approx1.459,\qquad \tilde{\mu}_{A,4}=4(1.459)\approx5.836,\qquad \frac{Y_{A,4}}{4^{0.2725}}=\frac{4}{1.459}\approx2.742.

Thus the observed count of 4 becomes approximately 2.742 count-equivalent units at the reference depth under this slope adjustment. It is not a new observed count and should not be interpreted as an exact estimate of the cell’s biological expression.

Worked Example: Computing a Pearson Residual

For the same Gene A cell, the NB variance at the regularized fitted mean is

\tilde{\mu}+\phi\tilde{\mu}^{2} =5.836+0.20(5.836)^2 \approx12.65.

The standard deviation is $\sqrt{12.65}\approx3.556$, and the Pearson residual is

r_{A,4} =\frac{4-5.836}{3.556} \approx-0.516.

The negative sign indicates that this observed count is below its fitted mean; the magnitude is about 0.516 fitted standard deviations. It is not a probability or a direct measure of biological effect.

For comparison, Gene B in the cell with $d=4$ has $\tilde{\beta}_B=1$, fitted mean $4(4^1)=16$, and observed count 16. Its adjusted observed count is $16/4=4$, and its Pearson residual is $(16-16)/\sqrt{16+0.20(16)^2}=0$. This zero follows from the deliberately exact fit in this hypothetical table.

Assumptions, Checks, and Common Pitfalls

The calculations require positive depth factors, a correctly specified NB mean-variance relationship, and counts that are conditionally independent across cells given the model quantities. They also treat depth factors and dispersion as known for this illustration. In real analyses, estimated factors, uncertain dispersion, cell dependence, and biological differences can all affect the fitted coefficients and their uncertainty. The normal approximation used for shrinkage is another assumption; it may be poor for weakly identified slopes.

Before using adjusted values or residuals, check whether the fit is plausible: inspect observed and fitted counts across the range of depth, look for extreme coefficients or unstable estimates, and check whether genes with very few nonzero counts are driving the result. A coefficient at an implausible extreme, a failure to converge, or a pattern of residuals that changes systematically with depth indicates that the model or its inputs may be inadequate. Residual diagnostics can assess remaining depth structure more directly; the calculations here do not establish that such structure has been removed.

Confusing a depth factor with a raw library size. The model uses the specified $d_i$; changing its scale changes the intercept reference. State whether factors are centered and how.
Assuming regularization is a universal formula. The shrinkage rule and target must be specified. Different choices produce different adjusted means.
Treating adjusted counts as residuals. Rescaling divides by a depth multiplier; residualizing subtracts a fitted mean and scales by model variability. They are different outputs for different purposes.
Reading a small residual as proof of a good model. A residual near zero for one cell-gene pair does not show that the model fits across cells or genes.

A regularized NB fit is useful when gene-wise depth responses vary but individual estimates are noisy. Its benefit is stability, not a guarantee of unbiased expression recovery. If depth is entangled with biology, shrinking or adjusting can also alter biological signal; model fit and the study’s assumptions must guide interpretation.

Check Your Understanding

Answer each question, then compare with the solutions below.

Why do the factors $(1/4,1/2,2,4)$ have geometric mean 1?
If $\hat{\beta}=0.5$, $SE=0.4$, $\beta_0=1$, and $\tau^2=0.25$, what is the regularized coefficient?
At $d=2$, with $\tilde{\beta}=0.5$, what is the depth multiplier?
What does a negative Pearson residual say about an observed count relative to its fitted mean?

Solutions. (1) Their product is $1/4\times1/2\times2\times4=1$, so the fourth root of the product is 1. (2) The sampling variance is $0.4^2=0.16$, and the weight is $0.25/(0.25+0.16)=0.6098$. Thus $\tilde{\beta}=1+0.6098(0.5-1)\approx0.6951$. (3) The multiplier is $2^{0.5}=\sqrt{2}\approx1.414$. (4) The observed count is below the fitted mean; the residual’s magnitude gives the difference in fitted standard-deviation units, not a probability.