Define marginal, residual, and graphical edges; assess uncertainty and computational cost
Before reading: Expression matrices, covariance, and estimands; follow influence-function links for sensitivity calculations.
Let X and Y be two measured features in a specified population, and let Z be a chosen adjustment set. Three targets recur:
| Edge target | Meaning | What zero does not generally establish |
|---|---|---|
| \operatorname{Corr}(X,Y) | Marginal association across units in that population | Independence or absence of a shared driver |
| E[(X-m_X(Z))(Y-m_Y(Z))] | Average conditional covariance after conditional-mean adjustment | Conditional independence or a direct biological interaction |
| A precision-matrix partial correlation | Linear residual association after adjusting for the other modeled features | Conditional independence outside the relevant Gaussian model, or causal regulation |
The normalized version of the second target is a residual correlation, not generally the average of conditional correlations. A differential-network estimand is a contrast between the same edge definition in two specified populations. State the scale, adjustment set, and feature universe before choosing the estimator. The expression examples below transfer only after the measurement assumptions for the new assay have been checked.
For Pearson correlation, normalization means
\rho=\frac{\operatorname{Cov}(X,Y)} {\sqrt{\operatorname{Var}(X)\operatorname{Var}(Y)}}.
Covariance retains the measurement units; correlation divides them out using each variable’s spread. Removing a conditional mean changes the variables being related, while normalizing changes their scale. Keep these operations separate. An average residual product, its normalized correlation, and a precision-matrix entry are three different edge definitions even when computed from the same assay.
A gene co-expression network has genes as nodes and, for each pair of genes, an edge weight that measures how strongly the two genes’ expression levels move together across samples. That last phrase is the one to hold onto. A population edge is a functional of a specified distribution, estimated from a cohort. Its estimator is not a direct measurement of regulation in an individual cell or molecule. Every property of the network is inherited from the cohort it was computed on, which is why the What is your estimator actually estimating? page applies to networks with full force.
Three decisions define a network: the association measure, whether edges are marginal or conditional, and how weights become a graph. Each has consequences for what the network can detect, how it behaves under heterogeneity, and what it costs to compute.
| Measure | Main use | Transformation behavior | Main caution |
|---|---|---|---|
| Pearson r | Linear association | Unchanged by positive affine rescaling of either variable; a sign reversal flips its sign | Sensitive to tails, mixtures, and mean shifts |
| Spearman \rho_S | Rank association | Unchanged by strictly increasing transformations of each variable | Subgroup shifts can change pooled ranks |
| Kendall \tau | Concordance | Unchanged by strictly increasing transformations of each variable | Ties and their handling matter |
| Biweight midcorrelation | Robust association | Uses robust centering, scaling, and outlier downweighting | Robustness depends on the definition and nondegenerate scale estimates |
| Distance correlation | General dependence | Translation and scalar-rescaling invariant; not invariant to arbitrary affine transformations | Population zero characterizes independence under moment conditions; finite-sample estimation still needs calibration |
| Mutual information | General dependence | Population quantity invariant to invertible transformations under suitable conditions | Estimated values depend on density, binning, and tuning choices |
| Partial correlation | Linear residual association | Depends on the conditioning variables and fitted covariance model | Zero implies conditional independence only in special models such as the Gaussian model |
Two properties matter most for our work.
Invariance to transformation. Strictly increasing transformations preserve pooled ranks. A sample-specific normalization, or a shift applied only within one group, need not preserve those ranks. Rank measures can therefore change when group composition changes. Invariance to a common transformation is not robustness to heterogeneity.
Influence and outliers. Pearson’s influence function has unbounded local sensitivity, although the correlation itself remains in [-1,1]. Rank-based and robust measures offer different forms of protection, but none establishes a biological interaction or removes an incorrect population target. The next example shows a large effect of one outlying pair on Pearson correlation.
set.seed(1)
n <- 50
x <- rnorm(n); y <- 0.3 * x + rnorm(n)
x[1] <- 6; y[1] <- 6 # one outlying sample
c(pearson = cor(x, y), spearman = cor(x, y, method = "spearman"),
pearson_without = cor(x[-1], y[-1])) pearson spearman pearson_without
0.5722860 0.1633133 0.2163756 A Pearson edge between genes A and B is marginal: it is nonzero if A and B co-vary for any reason, including both being driven by a third gene C, or by a cell-type proportion, or by batch. A conditional edge refers to a specified conditioning set. A full-feature graphical model conditions on the other modeled features; it does not follow that conditioning on every measured variable answers the biological question. Selection variables, mediators, and colliders can change associations. For Gaussian data the conditional edges are the partial correlations, read off the inverse covariance (precision) matrix \Omega = \Sigma^{-1}:
\rho_{AB \mid \text{rest}} = -\frac{\Omega_{AB}}{\sqrt{\Omega_{AA}\,\Omega_{BB}}}.
A zero partial correlation means conditional independence (under Gaussianity), and the graph of nonzero partial correlations is the Gaussian graphical model. Conditioning changes the estimand; it does not automatically make a network sparse or regulatory. Missing common causes, measurement error, selection, and the chosen feature universe can all change the interpretation.
“Rest” means the same specified set of the other p-2 features for this pair. If we keep only three genes, the A–B edge conditions on the third. If we add 100 genes, it conditions on a larger set and is a different population quantity. The minus sign in the precision formula comes from inverting the pair’s conditional covariance matrix: for a matrix with entries a,c,c,b, the inverse has off-diagonal entry -c/(ab-c^2). Dividing that precision entry by the square root of its two diagonal entries and changing the sign recovers c/\sqrt{ab}.
set.seed(2)
n <- 500
c_gene <- rnorm(n)
a <- 0.8 * c_gene + rnorm(n, sd = 0.6)
b <- 0.8 * c_gene + rnorm(n, sd = 0.6) # A and B share a driver C, no direct link
X <- cbind(a, b, c_gene)
cor(a, b) # marginal: strongly correlated[1] 0.661491Omega <- solve(cov(X))
-Omega[1, 2] / sqrt(Omega[1, 1] * Omega[2, 2]) # partial, given C: about zero[1] 0.002569492par(mfrow = c(2, 1), mar = c(4, 4, 2, 1))
plot(a, b, pch = 16, cex = 0.5, col = "#2F80C866", xlab = "Feature A", ylab = "Feature B",
main = sprintf("Marginal correlation: %.2f", cor(a, b)))
plot(a - 0.8*c_gene, b - 0.8*c_gene, pch = 16, cex = 0.5, col = "#D55E0066",
xlab = "A minus its true conditional mean", ylab = "B minus its true conditional mean",
main = "Shared driver removed")
In this shared-driver generator, A and B each have variance 0.8^2+0.6^2=1, and their covariance is 0.8^2\operatorname{Var}(C)=0.64. Their population marginal correlation is therefore 0.64. Subtracting 0.8C leaves the two independent noise terms, each with SD 0.6, giving residual covariance zero. The plotted sample correlation fluctuates around these population values. The conditional-independence interpretation here follows from the generator’s independent errors and Gaussian model, rather than from the label “adjusted.”
If X is symmetric about zero and Y=X^2, then E[X]=E[X^3]=0 and \operatorname{Cov}(X,Y)=0, despite a deterministic relationship. Taking Z constant makes this a zero residual-covariance example as well. A zero average conditional covariance can also arise by cancellation of positive and negative conditional covariances across strata. Gaussian conditional independence uses the joint Gaussian model, not zero covariance alone.
x_grid <- seq(-2, 2, length.out = 201)
plot(x_grid, x_grid^2, pch = 16, cex = 0.55, col = "#2F80C8",
xlab = "X", ylab = "Y = X squared", main = "Zero linear covariance; complete nonlinear dependence")
Let X_c be the centered N\times p data matrix. The sample covariance is S=X_c^\top X_c/(N-1). Centering makes the rows sum to zero, so its rank is at most N-1. Thus S is singular when p\ge N (and can be singular earlier through collinearity): there is no unique ordinary inverse to plug into the partial-correlation formula. Three standard approaches add structure:
glasso, or huge for a wrapper with selection.These methods regularize the estimated covariance, precision matrix, or nodewise regressions. Their tuning affects both weights and selected edges; individual transformed edge weights need not all move toward zero. Comparing two such networks estimated with the same \lambda but different sample sizes is not a like-for-like comparison. This is a trap in differential network analysis.
Regularization makes estimation possible by introducing a target, sparsity, or another structural assumption; it does not recover an arbitrary population precision matrix from insufficient information. A nonzero penalized edge is a model-selection result, not automatically a rejection of a null with a calibrated p-value. Missing values create another issue: computing each pair on a different subset can mix populations and produce a matrix that is not positive semidefinite. State how the common population and covariance estimate are defined before interpreting the graph.
An association matrix already defines a weighted network. Three common ways to transform it or select an unweighted graph are:
The significance route brings the multiple-testing problem into focus. With p = 10{,}000 genes there are about 5 \times 10^7 edges. BH needs independence or a suitable PRDS condition for the p-values. Shared-edge dependence need not meet those conditions; its presence alone does not prove they fail either. A permutation calibration can account for joint dependence only when its relabeling is justified by the tested null and study design. Between-group label permutation addresses a group-comparison null; it does not generally test that an individual population correlation is zero. See the Multiple testing when the tests are not independent page for alternatives and conditions.
First separate an edge difference from a difference in two significance labels. Suppose the estimated Fisher-transformed correlations are 0.30 (SE 0.14) and 0.28 (SE 0.15) in independent groups. Their individual normal-reference p-values are about 0.032 and 0.062: one crosses 0.05 and the other does not. But the estimated difference is only 0.02, with SE \sqrt{0.14^2+0.15^2}\approx0.205. There is little evidence of a difference. Testing the contrast directly answers the differential-edge question; subtracting two thresholded graphs does not.
The question “do the networks differ between patients and controls” has three main formulations, in increasing order of what they demand.
Edge by edge. For two independent samples, Fisher’s z difference,
\begin{aligned} d_z&=\operatorname{atanh}(r_1)-\operatorname{atanh}(r_2),\\ v_z&=\frac1{N_1-3}+\frac1{N_2-3},\\ T_z&=d_z/\sqrt{v_z}\\ &\approx\ N(0,1), \end{aligned}
tests one edge. The displayed reference is for independent i.i.d. samples from nondegenerate bivariate-normal populations under the equal-correlation null. Under other distributions its variance approximation can fail; equal group means or Gaussian marginals alone do not establish joint bivariate normality. Use an appropriate influence-function or bootstrap variance, or a permutation scheme justified for the actual null. Equal correlations alone do not make group labels exchangeable.
FisherZTest <- function(r1, n1, r2, n2) {
stopifnot(n1 > 3, n2 > 3, abs(r1) < 1, abs(r2) < 1)
z <- (atanh(r1) - atanh(r2)) / sqrt(1 / (n1 - 3) + 1 / (n2 - 3))
c(z = z, p = 2 * pnorm(-abs(z)))
}
FisherZTest(r1 = 0.6, n1 = 40, r2 = 0.2, n2 = 60) z p
2.32293885 0.02018244 Whole precision matrix. Methods that estimate the difference \Omega_1 - \Omega_2 directly (Zhao, Cai, and Li 2014) or test it entry-wise with a valid multiple-testing procedure (Xia, Cai, and Cai 2015). Direct difference estimation changes the estimation problem, but still requires tuning, sparsity/model assumptions, and an evaluation design appropriate to unequal sample sizes.
Sample-specific networks (the Jackknife pseudo-values page). Construct observation-indexed summaries and define the population contrast they are intended to estimate. Cohort-derived pseudo-values share estimated quantities, so treating them as independent raw measurements in a t-test or limma is not automatically valid. An inferential argument must handle their expansion, centering, and shared estimation, in addition to the mean/variance artifacts discussed in the linked guide.
Within-group sample correlations already center and scale within each group, so different group means alone do not create a difference between those correlations. Pooled pseudo-value constructions can inherit between-group shifts, as the Jackknife pseudo-values page demonstrates. Precision-matrix contrasts, correlation contrasts, and contrasts of observation-level network summaries are distinct targets; assess each on its own definition and assumptions.
The numbers you should know before planning an analysis.
Edges scale as p^2. A network on p genes has p(p-1)/2 edges: 499,500 for p = 1000, 2 \times 10^8 for p = 20{,}000. One dense network of 20,000 genes in double precision is 3.2 decimal GB (about 2.98 GiB), before temporary arrays; the unique upper-triangle entries require about 1.6 GB. Retaining one per person rapidly multiplies that cost. Restricting the feature universe changes the target graph and tested family; outcome-dependent filtering needs an inferential justification. Prespecified pathways or streaming edge summaries can reduce computation.
One correlation matrix is cheap. cor() on an N \times p matrix costs O(N p^2) and is dominated by matrix cross-products and associated centering/scaling operations.
set.seed(3)
N <- 100; p <- 1000
X <- matrix(rnorm(N * p), N, p)
system.time(R_all <- cor(X)) user system elapsed
0.029 0.001 0.030 Sample-specific networks the naive way are N times that, and the leave-one-out recomputation is wasteful: each cor(X[-i, ]) recomputes a cross-product that differs from the full one by a single rank-one term. Update the sufficient statistics instead.
# Leave-one-out correlation matrices for the rows in `idx`, from one cross-product,
# by rank-one updates. Assumes N > 2 and no constant columns.
LooCorAll <- function(X, idx = seq_len(nrow(X))) {
stopifnot(is.matrix(X), all(is.finite(X)), nrow(X) > 2,
length(idx) >= 1, all(idx == as.integer(idx)), all(idx >= 1), all(idx <= nrow(X)))
n <- nrow(X)
X <- sweep(X, 2, colMeans(X)) # translation preserves LOO correlations; reduces cancellation
s <- colSums(X); C <- crossprod(X)
lapply(idx, function(i) {
s_i <- s - X[i, ]
C_i <- C - tcrossprod(X[i, ])
V <- (C_i - tcrossprod(s_i) / (n - 1)) / (n - 2) # leave-one-out covariance
if (any(!is.finite(V)) || any(diag(V) <= 0)) stop("Degenerate leave-one-out variance")
V / sqrt(tcrossprod(diag(V))) # to correlation
})
}
i <- 7
all.equal(LooCorAll(X[, 1:50], idx = i)[[1]], cor(X[-i, 1:50]), check.attributes = FALSE)[1] TRUE# Keep this teaching benchmark small; storing many full large-feature matrices costs substantial memory.
X_bench <- X[, 1:250]
system.time(loo_naive <- lapply(1:20, function(i) cor(X_bench[-i, ]))) user system elapsed
0.038 0.001 0.039 system.time(loo_update <- LooCorAll(X_bench, idx = 1:20)) user system elapsed
0.011 0.006 0.014 all.equal(loo_naive[[1]], loo_update[[1]], check.attributes = FALSE)[1] TRUEFor the same 20 matrices the update version avoids repeated full-data cross-products; compare the timings on your machine rather than assuming a fixed speedup. The LIONESS network for sample i is then N * R_all - (N - 1) * loo[[i]], computed and summarized one sample at a time (call LooCorAll() with one idx at a time) so the number of retained dense matrices does not grow with the number of people. Cross-products, full estimates, updates, and temporary arrays still require several dense matrices; streaming is O(p^2) working memory, not literally one matrix. The update version is still O(p^2) per sample, which is the floor; it just avoids the factor of N from redoing the cross-product.
Permutation references multiply everything. A permutation test with B permutations costs B times the analysis. For a whole-network statistic on filtered genes this is routine; for per-edge permutation p-values across 10^7 edges it is not, and you need either a parametric reference you have checked by simulation or another calibrated approximation. Pooling edges requires assumptions about their null distributions and does not automatically solve the computational or inferential problem.
Follow the complete calculation in An end-to-end residual-network analysis, from a prespecified edge family and nuisance fits through uncertainty, multiplicity, and simulation evaluation.
Checkpoints are completion targets for this page. Additional practice below is optional. Each checkpoint gives a completion standard, a hint, and a solution outline.
Task type: Derivation and interpretation.
Let Z,U,V be independent standard normals, with X=Z+U and Y=Z+V. Calculate the marginal correlation and the covariance after subtracting E[X\mid Z] and E[Y\mid Z]. Explain the shared driver and why a marginal edge does not establish direct regulation. Give a dependent nonlinear pair with zero covariance. Is zero average residual covariance alone a general conditional-independence criterion?
Ready to continue when: you can distinguish marginal and residual targets and describe what the edge does and does not establish.
Transfer: What changes if U and V have covariance 0.3 instead of being independent?
Both marginal variances are 2. The conditional-mean residuals are U and V.
The covariance of X and Y is 1 and each variance is 2, so the marginal correlation is 1/2. After subtracting Z, the residual covariance is E[UV]=0. Shared variation in Z explains the marginal edge. In this constructed model the independent errors also imply conditional independence, but a zero average residual covariance does not establish it in general; conditional covariances can cancel, or dependence can be nonlinear. For example, a centered symmetric normal X and Y=X^2 are dependent with E[X^3]=0, hence zero covariance; take Z constant for a residual-covariance example.
These exercises extend the ideas through further derivations, simulations, or applications. Choose those useful for your work.
energy). Then add a mean shift to half the samples and see which measures move.huge has generators). Estimate the partial correlation network with shrinkage and with the graphical lasso at two values of \lambda. Count true and false edges.X. Start with equal within-group correlations and equal within-group variances and check the group contrast. Then change the variance of one gene in one group while keeping its correlation fixed. Does centering alone still remove the artifact?The original WGCNA paper describes network construction choices.