An end-to-end residual-network analysis

From a conditional covariance estimand to tested edges and a reproducible evaluation

Before reading: Network inference from omics data, Residual covariance and nuisance-estimation error, and Multiple testing when the tests are not independent. The residual-covariance page supplies the derivative, nuisance remainder, and inference conditions used here.

1 A question about co-variation after adjustment

Do two measured features co-vary differently between groups after removing their conditional means given recorded covariates? A marginal correlation can change because a shared predictor varies or predicts the features differently. Here we define the adjusted association before choosing a network algorithm.

This complete synthetic analysis uses 600 independent people, one fully observed assay per person, six continuous features, a binary group A, and covariate Z. Feature values represent a fixed log2-intensity scale. All IDs and measurements are invented. There is no molecular-interaction ground truth: the generator supplies statistical covariance truth for checking an estimator. The case study does not model assay detection, annotation, normalization estimation, or biological causality.

The downloadable base-R script contains data generation, alignment, cross-fitting, all edge tests, a failure-aware simulation, and output writing. It runs without additional packages. An end-to-end expression analysis is the complementary example with missingness and a standardized mean contrast.

2 Define the population edge and its contrast

For feature j and group g, let m_{jg}(z)=E[Y_j\mid Z=z,A=g]. Define

\begin{aligned} U_{jg}&=Y_j-m_{jg}(Z),\\ \psi_{jk,g}&=E[U_{jg}U_{kg}\mid A=g],\\ \Delta_{jk}&=\psi_{jk,1}-\psi_{jk,0}. \end{aligned}

The edge is an average conditional covariance, averaging over the covariate distribution within group g. Its units are the product of the two outcome units. It is not a correlation, a coefficient from a precision matrix, conditional independence, or a causal interaction. Rescaling a feature rescales its incident edges.

If conditional covariances vary with Z, the contrast can reflect both their differences and the groups’ different covariate distributions. A common-population standardized covariance contrast is a different estimand and requires a different estimation/identification argument. In this generator the conditional covariance is constant within group, making the truth simple. Both groups have the same covariate distribution, but this convenience is not an assumption silently imposed on arbitrary data.

Six features yield 15 unordered pairs. Prespecify all 15 as the testing family; do not select promising pairs from the same outcomes and then report ordinary p-values as if prespecified. Each null is H_{0,jk}:\Delta_{jk}=0, not a null that either group has zero covariance.

3 Derive the contribution before implementing it

Write W=(A,Z,Y_1,\ldots,Y_6) for one person’s complete record. Let \pi_g=P(A=g)>0 and R_{jk,g}=U_{jg}U_{kg}. For the nonparametric observed-data model, a mean-zero influence function for \psi_{jk,g} is

D_{jk,g}(W)=\frac{1\{A=g\}}{\pi_g} \{R_{jk,g}-\psi_{jk,g}\}.

To see the nuisance cancellation, differentiate the residual product along a smooth distributional path. Its mean-function derivative has terms involving E[U_{jg}\mid Z,A=g] or E[U_{kg}\mid Z,A=g], both zero. The remaining distribution derivative is the group-conditional mean derivative, including its group-probability denominator. Thus D_{jk}=D_{jk,1}-D_{jk,0} represents the contrast derivative. In the unrestricted observed-data model this is the efficient gradient when square integrable. Efficiency here refers to a regular-estimation asymptotic bound for this target and model.

The group factor can also be understood as the influence function of a ratio: \psi_{jk,g}=E[1\{A=g\}R_{jk,g}]/\pi_g. Changing the numerator alone would miss the change in the group probability. Differentiating both gives the centered residual product multiplied by 1\{A=g\}/\pi_g. People outside the group contribute zero to that group’s first-order term; people inside it represent the group through the inverse probability factor.

There is also an exact population remainder. For fixed trained mean estimates, put e_{jg}=\widehat m_{jg}-m_{jg}, and let \widehat R_{jk,g} denote the product of their two fitted residuals. Conditional mean-zero residuals give

\begin{aligned} E[\widehat R_{jk,g}\mid A=g]&=\psi_{jk,g}+E[e_{jg}e_{kg}\mid A=g],\\ |E[e_{jg}e_{kg}\mid A=g]|&\leq\|e_{jg}\|_{2,g}\|e_{kg}\|_{2,g}. \end{aligned}

The contrast bias is the group-1 product remainder minus the group-0 remainder. A small held-out prediction error is useful evidence, but does not alone establish the little-o product condition needed at the inference scale. The Orthogonality and cross-fitting page gives a counterexample where splitting fails to remove population bias.

4 Data, alignment, and held-out residual products

The generator has group-specific quadratic conditional means. Residual variances are one. Three disjoint feature pairs have covariance 0.2 in group 0. In group 1, those covariances are 0.65, −0.35, and 0.2. All other conditional cross-covariances are zero. The block matrices are positive definite. Hence only F1–F2 and F3–F4 have nonzero contrasts, 0.45 and −0.55.

Code

source("../examples/residual-network.R")
study <- GenerateResidualNetwork()
fit <- EstimateResidualNetwork(study$assay, study$metadata)
fit$group_counts


  0   1 
301 299

Code

fit$family_size

[1] 15

Code

head(fit$fold)

           id group fold
1 person_0001     0    1
2 person_0002     1    1
3 person_0003     0    1
4 person_0004     0    5
5 person_0005     1    1
6 person_0006     0    5

The assay has features in rows and person IDs in columns. The metadata are deliberately shuffled. The estimator matches IDs, rejects duplicate/mismatched IDs or nonfinite values, and checks group counts and training-design rank. Simulation truth is a separate object and is never passed to the estimator.

Five folds are formed within groups. For each held-out fold and group, the script fits all six means using only the other people in that group, with basis (1,Z,Z^2). It then multiplies pairs of held-out residuals. The group average of these products is \widehat\psi_{jk,g}, and their difference is \widehat\Delta_{jk}. Each person has one held-out product per edge. Training sets overlap; this does not make the fitted residual products independent finite-sample data.

For one edge, follow the calculation in this order: predict both features from other people; subtract those predictions from the held-out values; multiply the two residuals; average the products separately by group; subtract group 0 from group 1. The product can be negative even though both residual variances are positive: it records whether the two features depart from their fitted means in the same or opposite directions. The estimator averages these signed products, rather than fitting a separate biological network to each person.

The influence calculation justifies their first-order variance under further conditions. Require nonvanishing group probabilities, independent people, fixed feature family and fold count, finite/nonzero influence variances, negligible fold-weighted product remainders, and held-out score differences converging in L2. Suitable L4/moment control on residuals and fits can ensure the latter; nuisance L2 consistency by itself need not control every product. In this fixed-dimensional, correctly specified quadratic generator, the regular least-squares fits have the requisite moments/rates asymptotically. This is a teaching case, not a proof for an arbitrary machine learner.

5 Estimate uncertainty for the actual target

The asymptotic variance separates by group:

\operatorname{Var}(\widehat\Delta_{jk})\approx\frac{\operatorname{Var}(R_{jk,1}\mid A=1)}{n_1} +\frac{\operatorname{Var}(R_{jk,0}\mid A=0)}{n_0}.

The implementation substitutes sample variances of held-out products and observed group counts. These finite-sample degrees-of-freedom conventions are asymptotically equivalent to the full influence-function variance. They do not establish an exact t distribution for learned residuals. Normal p-values and 95% intervals are asymptotic and pointwise.

To see the scale, a bivariate-normal residual pair with unit variances and covariance c has product variance 1+c^2. For F1–F2 the oracle product variances are therefore 1.04 and 1.4225 in groups 0 and 1. With 300 people per group, the oracle SE of the difference is about \sqrt{1.04/300+1.4225/300}=0.091. The actual group counts and fitted residuals vary, so the implementation estimates this quantity. The spread of the individual residual products is roughly one; averaging over independent people reduces estimation uncertainty to roughly a tenth in this example.

Code

shown <- fit$results[c("edge", "delta", "SE", "lower", "upper", "BY", "discovery")]
knitr::kable(shown, digits = 3)

edge	delta	SE	lower	upper	BY	discovery
F1–F2	0.546	0.094	0.361	0.731	0.000	TRUE
F1–F3	0.126	0.080	-0.031	0.283	1.000	FALSE
F1–F4	-0.029	0.081	-0.189	0.130	1.000	FALSE
F1–F5	0.078	0.073	-0.066	0.221	1.000	FALSE
F1–F6	0.003	0.079	-0.153	0.158	1.000	FALSE
F2–F3	0.214	0.077	0.064	0.365	0.087	FALSE
F2–F4	0.068	0.083	-0.094	0.231	1.000	FALSE
F2–F5	-0.003	0.084	-0.168	0.162	1.000	FALSE
F2–F6	-0.144	0.082	-0.305	0.016	0.972	FALSE
F3–F4	-0.697	0.092	-0.878	-0.516	0.000	TRUE
F3–F5	0.007	0.081	-0.152	0.165	1.000	FALSE
F3–F6	0.046	0.091	-0.133	0.225	1.000	FALSE
F4–F5	-0.018	0.081	-0.177	0.140	1.000	FALSE
F4–F6	-0.058	0.086	-0.227	0.112	1.000	FALSE
F5–F6	-0.045	0.091	-0.223	0.133	1.000	FALSE

BY adjusts over the 15 prespecified tests. It handles arbitrary dependence among valid marginal p-values; our normal p-values have asymptotic justification under the conditions above. With a fixed finite family, that supports an asymptotic error guarantee, not exact finite-sample control. Pointwise 95% intervals are neither simultaneous intervals nor intervals with automatic post-selection coverage.

A raw 95% interval can exclude zero while BY does not declare a discovery, because these are different decision rules. BY applies the step-up procedure at 0.05/H_{15}, with H_{15}\approx3.32, to the entire family. The table retains all edges so the reader can distinguish effect size, pointwise uncertainty, and the family-level reporting rule. In real data the number of false discoveries is unknown; only the generator makes the realized FDP available for the evaluation below.

Code

r <- fit$results
truth <- study$truth$delta[match(r$edge, study$truth$edge)]
par(mar = c(4, 6, 1, 1))
plot(r$delta, seq_len(nrow(r)), xlim = range(c(r$lower, r$upper, truth)),
     yaxt = "n", ylab = "", xlab = "Group 1 minus group 0 residual covariance",
     pch = ifelse(r$discovery, 17, 16), col = ifelse(r$discovery, "#D55E00", "#2F80C8"))
segments(r$lower, seq_len(nrow(r)), r$upper, seq_len(nrow(r)))
points(truth, seq_len(nrow(r)), pch = 4, cex = 1.2)
axis(2, seq_len(nrow(r)), r$edge, las = 1)
abline(v = 0, lty = 2)
legend("topright", c("BY discovery", "Other edge", "Truth"), pch = c(17, 16, 4),
       col = c("#D55E00", "#2F80C8", "black"), bty = "n", cex = 0.8)

A forest plot shows 15 edge contrasts and their pointwise intervals, with truth crosses and a reference at zero. — Figure 1: The synthetic reference dataset: point estimates and pointwise 95% normal intervals for all 15 residual-covariance contrasts. Orange triangles mark BY discoveries; blue circles mark other edges. Crosses are generator truth, used only for evaluation. Interval exclusion and BY discovery are different criteria.

6 Check why marginal and residual networks disagree

The figure below uses F1–F3, whose true residual covariance contrast is zero. Their conditional mean patterns differ across groups, so a raw covariance comparison concerns another target. The residual plot removes estimated means; it does not discover an underlying molecular interaction graph.

Code

d <- study$metadata[match(colnames(study$assay), study$metadata$id), ]
cols <- c("#2F80C899", "#D55E0099")[d$a + 1]
par(mfrow = c(2, 1), mar = c(4, 4, 1, 1))
plot(study$assay['F1',], study$assay['F3',], col = cols, pch = 16,
     xlab = "Raw F1", ylab = "Raw F3")
legend("topleft", c("Group 0", "Group 1"), pch = 16,
       col = c("#2F80C8", "#D55E00"), bty = "n")
plot(fit$residual[, 'F1'], fit$residual[, 'F3'], col = cols, pch = 16,
     xlab = "Held-out F1 residual", ylab = "Held-out F3 residual")

Two stacked scatterplots compare F1 and F3 raw values with their held-out residuals, coloring the two groups. — Figure 2: One synthetic null residual edge (F1–F3), before and after group-specific covariate adjustment. The raw cloud includes shared mean structure. The fitted residual cloud illustrates the chosen target, not an exact finite-sample independence claim.

7 Challenge the estimator with a paired simulation

Use 200 independently generated datasets at each of n=300 and 900. Within each dataset, compare the correct quadratic nuisance basis with an intercept-only basis, using the same folds. The misspecified fits leave covariate mean structure in the residuals; cross-fitting still separates training people but does not repair that remainder.

The script retains every attempted dataset/model fit and its error message. Summaries report successes and failures. If any fail, performance conditional on success is not the whole procedure’s performance; investigate the failure pattern before making a claim. The reference run below has no failures.

Code

performance <- EvaluateResidualNetwork(replicates = 200)
performance$summary

    n     model attempted success failures sentinel_bias   bias_MCSE  coverage
1 300 quadratic       200     200        0  -0.002730365 0.008971570 0.9466667
2 300 intercept       200     200        0   0.363797730 0.012437164 0.7280000
3 900 quadratic       200     200        0   0.004981226 0.005135492 0.9523333
4 900 intercept       200     200        0   0.374269563 0.006605971 0.4770000
  coverage_MCSE         FDR    FDR_MCSE  power
1   0.004357041 0.015833333 0.005230251 0.8250
2   0.009467804 0.226154762 0.018778409 0.7875
3   0.004289301 0.009166667 0.004126270 1.0000
4   0.008534065 0.603067100 0.012053399 1.0000

Coverage is averaged over the 15 pointwise intervals within each dataset, then across datasets. FDR is the mean realized FDP, using zero FDP for no discoveries. Their MCSEs use variation across datasets, preserving within-dataset edge dependence. The sentinel F1–F3 is a prespecified null edge for showing persistent target bias. Power is averaged over the two nonzero edges; it is a separate criterion from FDR.

Code

s <- performance$summary
par(mar = c(4, 4, 1, 1))
plot(range(s$n), range(c(s$sentinel_bias - 2*s$bias_MCSE,
                        s$sentinel_bias + 2*s$bias_MCSE, 0)), type = "n",
     xlab = "Independent people", ylab = "Mean estimated contrast minus truth")
for (model in c("quadratic", "intercept")) {
  t <- s[s$model == model, ]; color <- if (model == "quadratic") "#2F80C8" else "#D55E00"
  lines(t$n, t$sentinel_bias, type = "b", pch = 16, col = color)
  arrows(t$n, t$sentinel_bias - 2*t$bias_MCSE, t$n, t$sentinel_bias + 2*t$bias_MCSE,
         angle = 90, code = 3, length = 0.05, col = color)
}
abline(h = 0, lty = 2)
legend("right", c("Quadratic means", "Intercept only"),
       col = c("#2F80C8", "#D55E00"), lty = 1, pch = 16, bty = "n")

Bias points and error bars compare correct quadratic fits with intercept-only fits at sample sizes 300 and 900. — Figure 3: Estimated bias for the prespecified null residual edge F1–F3, plus/minus two Monte Carlo SEs across 200 datasets at each size. The intercept-only fit targets a different covariance even though it uses cross-fitting. This finite experiment supports diagnosis, not a general coverage theorem.

8 Reproduce the analysis and audit the claim

Download residual-network.R, save it in a project, and run the command below in a terminal. It writes the synthetic assay, shuffled metadata, separate truth, folds, all 15 edge results, simulation attempts, evaluation edge records/summaries, and R session information. The fixed data and fold seeds are visible in the script. The default evaluation uses 200 replicates per sample size/model; replace the final argument with 3 for a workflow smoke test, not a precision evaluation.

Rscript residual-network.R results 200

A defensible report says: “We compared average group-specific conditional covariances of the defined log2-intensity features, adjusting their means for Z, over all 15 prespecified pairs. We used person-level cross-fitting, influence-based pointwise uncertainty, and BY-adjusted asymptotic p-values.” It also states the measurement limitations and simulation scope. It does not call the discoveries physical interactions or claim that nonsignificance proves equality.

For repeated specimens, first define whether the target averages people, visits, scheduled times, or within-person changes. The contribution and weights must match that choice. Keep all records from a person in a fold, and use whole-person contributions under independent-person sampling. Simply plugging repeated row products into this script’s SE does not implement that extension; the Dependent observations and cluster-aware inference page provides the necessary framework.

9 Exercises

Try each checkpoint before opening its hint or solution.

Checkpoint · NETCASE-01 — Derive the residual-edge contrast and its contribution

Task type: Derivation and interpretation.

For one prespecified pair, define ψ_g and Δ, derive D_g including its group-probability factor, and derive the exact nuisance product remainder. Explain why this is neither a partial-correlation edge nor necessarily a common-covariate-population comparison.

Ready to continue when: you can connect the edge definition, centering, nuisance cancellation, and remainder to the estimator actually implemented.

Transfer: How would multiplying F1 by ten change the edge and its SE?

Hint

Differentiate the conditional group average, then use the zero conditional means of the true residuals.

Solution outline — open after trying

ψ_g=E[(Y_j−m_jg)(Y_k−m_kg)|A=g] and Δ=ψ_1−ψ_0. D_g=1(A=g)(R_g−ψ_g)/π_g, and D=D_1−D_0. Expanding trained residual products leaves E[e_jg e_kg|g] after the two cross terms vanish. Average conditional covariance is not precision-matrix partial correlation. The outer covariate law is group-specific; standardizing both groups to a common law changes the target.

Checkpoint · NETCASE-02 — Audit uncertainty and reproduce a failure case

Task type: Implementation and diagnosis.

Rerun the script, check ID alignment and all 15 family members, and locate the pointwise intervals, BY decisions, failure log, and MCSEs. Compare the quadratic and intercept-only fits for the sentinel null edge. State the independent unit and what must change for repeated specimens.

Ready to continue when: you can trace each conclusion to its target, inputs, conditions, and evaluation record rather than relying on a plot of discoveries.

Transfer: If no edge is significant after BY, what evidence would you need before claiming the two covariance networks are equivalent?

Hint

Compare a study SE with a simulation MCSE, and compare training separation with population remainder control.

Solution outline — open after trying

The implementation matches metadata by unique person IDs and evaluates all 15 pairs. Pointwise intervals do not give simultaneous or post-selection coverage. BY needs valid marginal p-values; normal approximations here are asymptotic. Monte Carlo SEs quantify estimated simulation metrics. Intercept-only means leave a nonvanishing product remainder despite cross-fitting. Repeated specimens need a newly specified person/visit target, cluster contributions and person-level folds, rather than treating extra rows as independent.

Additional practice (optional)

Derive the sentinel F1–F3 marginal covariance contrast under the generator and compare it with the simulation bias.
Change the random seed and report changes in point estimates, intervals, and BY discoveries without changing the tested family.
Design a separate estimator for a common-population standardized covariance contrast; identify its additional nuisance functions before coding.
Prespecify an equivalence margin for one edge and explain how its inferential question differs from testing Δ=0.

10 Reading

Kennedy’s semiparametric review: influence functions and estimator construction.
R p.adjust documentation: multiplicity procedures and conditions.

11 Continue through the guide

Compare An end-to-end expression analysis, extend the unit argument with Dependent observations and cluster-aware inference, and use Reporting a statistical methods study to match claims to evidence.

--- title: "An end-to-end residual-network analysis" subtitle: "From a conditional covariance estimand to tested edges and a reproducible evaluation" guide-status: "Available" --- **Before reading:** [Network inference from omics data](network-inference.qmd), [Residual covariance and nuisance-estimation error](residual-covariance-and-nuisance-error.qmd), and [Multiple testing when the tests are not independent](multiple-testing.qmd). The residual-covariance page supplies the derivative, nuisance remainder, and inference conditions used here. ## A question about co-variation after adjustment Do two measured features co-vary differently between groups after removing their conditional means given recorded covariates? A marginal correlation can change because a shared predictor varies or predicts the features differently. Here we define the adjusted association before choosing a network algorithm. This complete **synthetic** analysis uses 600 independent people, one fully observed assay per person, six continuous features, a binary group $A$, and covariate $Z$. Feature values represent a fixed log2-intensity scale. All IDs and measurements are invented. There is no molecular-interaction ground truth: the generator supplies statistical covariance truth for checking an estimator. The case study does not model assay detection, annotation, normalization estimation, or biological causality. The [downloadable base-R script](../examples/residual-network.R) contains data generation, alignment, cross-fitting, all edge tests, a failure-aware simulation, and output writing. It runs without additional packages. [An end-to-end expression analysis](adjusted-expression-case-study.qmd) is the complementary example with missingness and a standardized mean contrast. ## Define the population edge and its contrast For feature $j$ and group $g$, let $m_{jg}(z)=E[Y_j\mid Z=z,A=g]$. Define $$ \begin{aligned} U_{jg}&=Y_j-m_{jg}(Z),\\ \psi_{jk,g}&=E[U_{jg}U_{kg}\mid A=g],\\ \Delta_{jk}&=\psi_{jk,1}-\psi_{jk,0}. \end{aligned} $$ The edge is an **average conditional covariance**, averaging over the covariate distribution within group $g$. Its units are the product of the two outcome units. It is not a correlation, a coefficient from a precision matrix, conditional independence, or a causal interaction. Rescaling a feature rescales its incident edges. If conditional covariances vary with $Z$, the contrast can reflect both their differences and the groups' different covariate distributions. A common-population standardized covariance contrast is a different estimand and requires a different estimation/identification argument. In this generator the conditional covariance is constant within group, making the truth simple. Both groups have the same covariate distribution, but this convenience is not an assumption silently imposed on arbitrary data. Six features yield 15 unordered pairs. Prespecify all 15 as the testing family; do not select promising pairs from the same outcomes and then report ordinary p-values as if prespecified. Each null is $H_{0,jk}:\Delta_{jk}=0$, not a null that either group has zero covariance. ## Derive the contribution before implementing it Write $W=(A,Z,Y_1,\ldots,Y_6)$ for one person's complete record. Let $\pi_g=P(A=g)>0$ and $R_{jk,g}=U_{jg}U_{kg}$. For the nonparametric observed-data model, a mean-zero influence function for $\psi_{jk,g}$ is $$ D_{jk,g}(W)=\frac{1\{A=g\}}{\pi_g} \{R_{jk,g}-\psi_{jk,g}\}. $$ To see the nuisance cancellation, differentiate the residual product along a smooth distributional path. Its mean-function derivative has terms involving $E[U_{jg}\mid Z,A=g]$ or $E[U_{kg}\mid Z,A=g]$, both zero. The remaining distribution derivative is the group-conditional mean derivative, including its group-probability denominator. Thus $D_{jk}=D_{jk,1}-D_{jk,0}$ represents the contrast derivative. In the unrestricted observed-data model this is the efficient gradient when square integrable. Efficiency here refers to a regular-estimation asymptotic bound for this target and model. The group factor can also be understood as the influence function of a ratio: $\psi_{jk,g}=E[1\{A=g\}R_{jk,g}]/\pi_g$. Changing the numerator alone would miss the change in the group probability. Differentiating both gives the centered residual product multiplied by $1\{A=g\}/\pi_g$. People outside the group contribute zero to that group's first-order term; people inside it represent the group through the inverse probability factor. There is also an exact population remainder. For fixed trained mean estimates, put $e_{jg}=\widehat m_{jg}-m_{jg}$, and let $\widehat R_{jk,g}$ denote the product of their two fitted residuals. Conditional mean-zero residuals give $$ \begin{aligned} E[\widehat R_{jk,g}\mid A=g]&=\psi_{jk,g}+E[e_{jg}e_{kg}\mid A=g],\\ |E[e_{jg}e_{kg}\mid A=g]|&\leq\|e_{jg}\|_{2,g}\|e_{kg}\|_{2,g}. \end{aligned} $$ The contrast bias is the group-1 product remainder minus the group-0 remainder. A small held-out prediction error is useful evidence, but does not alone establish the little-o product condition needed at the inference scale. The [Orthogonality and cross-fitting](../semiparametric-methods/orthogonality-and-cross-fitting.qmd) page gives a counterexample where splitting fails to remove population bias. ## Data, alignment, and held-out residual products The generator has group-specific quadratic conditional means. Residual variances are one. Three disjoint feature pairs have covariance 0.2 in group 0. In group 1, those covariances are 0.65, −0.35, and 0.2. All other conditional cross-covariances are zero. The block matrices are positive definite. Hence only F1–F2 and F3–F4 have nonzero contrasts, 0.45 and −0.55. ```{r} source("../examples/residual-network.R") study <- GenerateResidualNetwork() fit <- EstimateResidualNetwork(study$assay, study$metadata) fit$group_counts fit$family_size head(fit$fold) ``` The assay has features in rows and person IDs in columns. The metadata are deliberately shuffled. The estimator matches IDs, rejects duplicate/mismatched IDs or nonfinite values, and checks group counts and training-design rank. Simulation truth is a separate object and is never passed to the estimator. Five folds are formed within groups. For each held-out fold and group, the script fits all six means using only the other people in that group, with basis $(1,Z,Z^2)$. It then multiplies pairs of held-out residuals. The group average of these products is $\widehat\psi_{jk,g}$, and their difference is $\widehat\Delta_{jk}$. Each person has one held-out product per edge. Training sets overlap; this does not make the fitted residual products independent finite-sample data. For one edge, follow the calculation in this order: predict both features from other people; subtract those predictions from the held-out values; multiply the two residuals; average the products separately by group; subtract group 0 from group 1. The product can be negative even though both residual variances are positive: it records whether the two features depart from their fitted means in the same or opposite directions. The estimator averages these signed products, rather than fitting a separate biological network to each person. The influence calculation justifies their first-order variance under further conditions. Require nonvanishing group probabilities, independent people, fixed feature family and fold count, finite/nonzero influence variances, negligible fold-weighted product remainders, and held-out score differences converging in L2. Suitable L4/moment control on residuals and fits can ensure the latter; nuisance L2 consistency by itself need not control every product. In this fixed-dimensional, correctly specified quadratic generator, the regular least-squares fits have the requisite moments/rates asymptotically. This is a teaching case, not a proof for an arbitrary machine learner. ## Estimate uncertainty for the actual target The asymptotic variance separates by group: $$ \operatorname{Var}(\widehat\Delta_{jk})\approx\frac{\operatorname{Var}(R_{jk,1}\mid A=1)}{n_1} +\frac{\operatorname{Var}(R_{jk,0}\mid A=0)}{n_0}. $$ The implementation substitutes sample variances of held-out products and observed group counts. These finite-sample degrees-of-freedom conventions are asymptotically equivalent to the full influence-function variance. They do not establish an exact t distribution for learned residuals. Normal p-values and 95% intervals are asymptotic and pointwise. To see the scale, a bivariate-normal residual pair with unit variances and covariance c has product variance $1+c^2$. For F1–F2 the oracle product variances are therefore 1.04 and 1.4225 in groups 0 and 1. With 300 people per group, the oracle SE of the difference is about $\sqrt{1.04/300+1.4225/300}=0.091$. The actual group counts and fitted residuals vary, so the implementation estimates this quantity. The spread of the individual residual products is roughly one; averaging over independent people reduces estimation uncertainty to roughly a tenth in this example. ::: {.table-responsive} ```{r} shown <- fit$results[c("edge", "delta", "SE", "lower", "upper", "BY", "discovery")] knitr::kable(shown, digits = 3) ``` ::: BY adjusts over the 15 prespecified tests. It handles arbitrary dependence **among valid marginal p-values**; our normal p-values have asymptotic justification under the conditions above. With a fixed finite family, that supports an asymptotic error guarantee, not exact finite-sample control. Pointwise 95% intervals are neither simultaneous intervals nor intervals with automatic post-selection coverage. A raw 95% interval can exclude zero while BY does not declare a discovery, because these are different decision rules. BY applies the step-up procedure at $0.05/H_{15}$, with $H_{15}\approx3.32$, to the entire family. The table retains all edges so the reader can distinguish effect size, pointwise uncertainty, and the family-level reporting rule. In real data the number of false discoveries is unknown; only the generator makes the realized FDP available for the evaluation below. ```{r} #| label: fig-network-case-intervals #| fig-height: 6 #| fig-cap: "The synthetic reference dataset: point estimates and pointwise 95% normal intervals for all 15 residual-covariance contrasts. Orange triangles mark BY discoveries; blue circles mark other edges. Crosses are generator truth, used only for evaluation. Interval exclusion and BY discovery are different criteria." #| fig-alt: "A forest plot shows 15 edge contrasts and their pointwise intervals, with truth crosses and a reference at zero." r <- fit$results truth <- study$truth$delta[match(r$edge, study$truth$edge)] par(mar = c(4, 6, 1, 1)) plot(r$delta, seq_len(nrow(r)), xlim = range(c(r$lower, r$upper, truth)), yaxt = "n", ylab = "", xlab = "Group 1 minus group 0 residual covariance", pch = ifelse(r$discovery, 17, 16), col = ifelse(r$discovery, "#D55E00", "#2F80C8")) segments(r$lower, seq_len(nrow(r)), r$upper, seq_len(nrow(r))) points(truth, seq_len(nrow(r)), pch = 4, cex = 1.2) axis(2, seq_len(nrow(r)), r$edge, las = 1) abline(v = 0, lty = 2) legend("topright", c("BY discovery", "Other edge", "Truth"), pch = c(17, 16, 4), col = c("#D55E00", "#2F80C8", "black"), bty = "n", cex = 0.8) ``` ## Check why marginal and residual networks disagree The figure below uses F1–F3, whose true residual covariance contrast is zero. Their conditional mean patterns differ across groups, so a raw covariance comparison concerns another target. The residual plot removes estimated means; it does not discover an underlying molecular interaction graph. ```{r} #| label: fig-network-case-residuals #| fig-height: 6.5 #| fig-cap: "One synthetic null residual edge (F1–F3), before and after group-specific covariate adjustment. The raw cloud includes shared mean structure. The fitted residual cloud illustrates the chosen target, not an exact finite-sample independence claim." #| fig-alt: "Two stacked scatterplots compare F1 and F3 raw values with their held-out residuals, coloring the two groups." d <- study$metadata[match(colnames(study$assay), study$metadata$id), ] cols <- c("#2F80C899", "#D55E0099")[d$a + 1] par(mfrow = c(2, 1), mar = c(4, 4, 1, 1)) plot(study$assay['F1',], study$assay['F3',], col = cols, pch = 16, xlab = "Raw F1", ylab = "Raw F3") legend("topleft", c("Group 0", "Group 1"), pch = 16, col = c("#2F80C8", "#D55E00"), bty = "n") plot(fit$residual[, 'F1'], fit$residual[, 'F3'], col = cols, pch = 16, xlab = "Held-out F1 residual", ylab = "Held-out F3 residual") ``` ## Challenge the estimator with a paired simulation Use 200 independently generated datasets at each of $n=300$ and $900$. Within each dataset, compare the correct quadratic nuisance basis with an intercept-only basis, using the same folds. The misspecified fits leave covariate mean structure in the residuals; cross-fitting still separates training people but does not repair that remainder. The script retains every attempted dataset/model fit and its error message. Summaries report successes and failures. If any fail, performance conditional on success is not the whole procedure's performance; investigate the failure pattern before making a claim. The reference run below has no failures. ```{r} #| label: network-case-evaluation #| cache: true performance <- EvaluateResidualNetwork(replicates = 200) performance$summary ``` Coverage is averaged over the 15 pointwise intervals **within each dataset**, then across datasets. FDR is the mean realized FDP, using zero FDP for no discoveries. Their MCSEs use variation across datasets, preserving within-dataset edge dependence. The sentinel F1–F3 is a prespecified null edge for showing persistent target bias. Power is averaged over the two nonzero edges; it is a separate criterion from FDR. ```{r} #| label: fig-network-case-evaluation #| fig-height: 5 #| fig-cap: "Estimated bias for the prespecified null residual edge F1–F3, plus/minus two Monte Carlo SEs across 200 datasets at each size. The intercept-only fit targets a different covariance even though it uses cross-fitting. This finite experiment supports diagnosis, not a general coverage theorem." #| fig-alt: "Bias points and error bars compare correct quadratic fits with intercept-only fits at sample sizes 300 and 900." s <- performance$summary par(mar = c(4, 4, 1, 1)) plot(range(s$n), range(c(s$sentinel_bias - 2*s$bias_MCSE, s$sentinel_bias + 2*s$bias_MCSE, 0)), type = "n", xlab = "Independent people", ylab = "Mean estimated contrast minus truth") for (model in c("quadratic", "intercept")) { t <- s[s$model == model, ]; color <- if (model == "quadratic") "#2F80C8" else "#D55E00" lines(t$n, t$sentinel_bias, type = "b", pch = 16, col = color) arrows(t$n, t$sentinel_bias - 2*t$bias_MCSE, t$n, t$sentinel_bias + 2*t$bias_MCSE, angle = 90, code = 3, length = 0.05, col = color) } abline(h = 0, lty = 2) legend("right", c("Quadratic means", "Intercept only"), col = c("#2F80C8", "#D55E00"), lty = 1, pch = 16, bty = "n") ``` ## Reproduce the analysis and audit the claim Download [residual-network.R](../examples/residual-network.R), save it in a project, and run the command below in a terminal. It writes the synthetic assay, shuffled metadata, separate truth, folds, all 15 edge results, simulation attempts, evaluation edge records/summaries, and R session information. The fixed data and fold seeds are visible in the script. The default evaluation uses 200 replicates per sample size/model; replace the final argument with 3 for a workflow smoke test, not a precision evaluation. ```{.bash} Rscript residual-network.R results 200 ``` A defensible report says: “We compared average group-specific conditional covariances of the defined log2-intensity features, adjusting their means for Z, over all 15 prespecified pairs. We used person-level cross-fitting, influence-based pointwise uncertainty, and BY-adjusted asymptotic p-values.” It also states the measurement limitations and simulation scope. It does not call the discoveries physical interactions or claim that nonsignificance proves equality. For repeated specimens, first define whether the target averages people, visits, scheduled times, or within-person changes. The contribution and weights must match that choice. Keep all records from a person in a fold, and use whole-person contributions under independent-person sampling. Simply plugging repeated row products into this script's SE does not implement that extension; the [Dependent observations and cluster-aware inference](../estimation-and-inference/dependent-observations.qmd) page provides the necessary framework. ## Exercises Try each checkpoint before opening its hint or solution. ::: {#checkpoint-netcase-01 .checkpoint} ### Checkpoint · NETCASE-01 — Derive the residual-edge contrast and its contribution {.unnumbered} **Task type:** Derivation and interpretation. For one prespecified pair, define ψ_g and Δ, derive D_g including its group-probability factor, and derive the exact nuisance product remainder. Explain why this is neither a partial-correlation edge nor necessarily a common-covariate-population comparison. **Ready to continue when:** you can connect the edge definition, centering, nuisance cancellation, and remainder to the estimator actually implemented. **Transfer:** How would multiplying F1 by ten change the edge and its SE? <details class="checkpoint-help"> <summary>Hint</summary> Differentiate the conditional group average, then use the zero conditional means of the true residuals. </details> <details class="checkpoint-help"> <summary>Solution outline — open after trying</summary> ψ_g=E[(Y_j−m_jg)(Y_k−m_kg)|A=g] and Δ=ψ_1−ψ_0. D_g=1(A=g)(R_g−ψ_g)/π_g, and D=D_1−D_0. Expanding trained residual products leaves E[e_jg e_kg|g] after the two cross terms vanish. Average conditional covariance is not precision-matrix partial correlation. The outer covariate law is group-specific; standardizing both groups to a common law changes the target. </details> ::: ::: {#checkpoint-netcase-02 .checkpoint} ### Checkpoint · NETCASE-02 — Audit uncertainty and reproduce a failure case {.unnumbered} **Task type:** Implementation and diagnosis. Rerun the script, check ID alignment and all 15 family members, and locate the pointwise intervals, BY decisions, failure log, and MCSEs. Compare the quadratic and intercept-only fits for the sentinel null edge. State the independent unit and what must change for repeated specimens. **Ready to continue when:** you can trace each conclusion to its target, inputs, conditions, and evaluation record rather than relying on a plot of discoveries. **Transfer:** If no edge is significant after BY, what evidence would you need before claiming the two covariance networks are equivalent? <details class="checkpoint-help"> <summary>Hint</summary> Compare a study SE with a simulation MCSE, and compare training separation with population remainder control. </details> <details class="checkpoint-help"> <summary>Solution outline — open after trying</summary> The implementation matches metadata by unique person IDs and evaluates all 15 pairs. Pointwise intervals do not give simultaneous or post-selection coverage. BY needs valid marginal p-values; normal approximations here are asymptotic. Monte Carlo SEs quantify estimated simulation metrics. Intercept-only means leave a nonvanishing product remainder despite cross-fitting. Repeated specimens need a newly specified person/visit target, cluster contributions and person-level folds, rather than treating extra rows as independent. </details> ::: ### Additional practice (optional) {.unnumbered} 1. Derive the sentinel F1–F3 marginal covariance contrast under the generator and compare it with the simulation bias. 2. Change the random seed and report changes in point estimates, intervals, and BY discoveries without changing the tested family. 3. Design a separate estimator for a common-population standardized covariance contrast; identify its additional nuisance functions before coding. 4. Prespecify an equivalence margin for one edge and explain how its inferential question differs from testing Δ=0. ## Reading - [Kennedy’s semiparametric review](https://arxiv.org/abs/2203.06469): influence functions and estimator construction. - [R p.adjust documentation](https://stat.ethz.ch/R-manual/R-devel/library/stats/html/p.adjust.html): multiplicity procedures and conditions. ## Continue through the guide {#continue-through-the-collection} Compare [An end-to-end expression analysis](adjusted-expression-case-study.qmd), extend the unit argument with [Dependent observations and cluster-aware inference](../estimation-and-inference/dependent-observations.qmd), and use [Reporting a statistical methods study](../research-practice/reporting-methods-studies.qmd) to match claims to evidence.