An end-to-end expression analysis

From a standardized contrast to semiparametric inference and a reproducible report

Before reading: Read Preprocessing and the quantity being estimated, From a research question to an identified estimand, and From an influence function to an estimator. The target and data-audit sections can be read before completing all the proofs.

1 One research argument, from measurements to a reported contrast

We will analyze 24 synthetic expression features from independently sampled people. Group membership is associated with a recorded covariate and technical batch. Some whole-assay measurements are missing. The question is: how do the groups differ on the supplied log2 expression scale after standardizing to a common covariate distribution?

This extended worked application follows one argument: define the target, audit the data, construct its estimator and uncertainty, then assess the reported family of claims. Simulation checks specific predictions of that argument. The From an influence function to an estimator page supplies the proof pattern; here we combine group standardization with missing outcomes.

Everything in this example is synthetic. Feature IDs are invented. The generator supplies known truth for evaluation, while the estimator receives only the observed assay and metadata. No clinical or biological finding is being asserted.

2 Define the population, scale, and target

Let A be group (0 or 1), C=(Z,B) the recorded covariate and batch, and Y_g a feature’s calibrated log2 intensity. Define

\begin{aligned} \mu_{ag}(c)&=E[Y_g\mid A=a,C=c],\\ \psi_g&=E_C[\mu_{1g}(C)-\mu_{0g}(C)]. \end{aligned}

The outer distribution is that of all eligible people in the source population, estimated using everyone’s recorded C, including people with missing assays. We do not condition on the realized covariate values as a fixed population. Variation in the sampled covariate mix therefore contributes to the influence function.

Decision	Specification in this case
Independent unit	One person, one specimen, one assay occasion
Population	The source population sampled by the synthetic generator
Measurement	Supplied calibrated log2 intensity; no learned normalization
Contrast	Mean conditional group difference standardized over the source population’s C distribution
Missingness	One recorded indicator R for a missing whole assay; all covariates remain observed
Family	All 24 prespecified feature hypotheses H_{0g}:\psi_g=0
Interpretation	A standardized association; no causal interpretation is assumed

A log2 contrast exponentiates to a ratio of geometric means of the two standardized intensity distributions. It is not a ratio of arithmetic mean expression or an absolute molecules-per-cell contrast. The Preprocessing and the quantity being estimated page explains those distinctions. Real count data would need their own measurement and normalization argument before this analysis could be reused.

Group contrasts can vary with C. We do not impose a constant coefficient or a common regression slope across groups. Nor does adjustment alone turn this association into an intervention effect. A causal interpretation would need a well-defined intervention and additional identification assumptions.

3 Reconstruct the observation process

We observe W=(C,A,R,RY_1,\ldots,RY_G). Algebraically RY_g contributes zero when R=0; the stored assay uses NA for those outcomes. Treating those outcomes as measured zeros would change the analysis.

Assume E[|Y_g|]<\infty, outcome missingness at random given (A,C), and positive group and observation probabilities on the target support:

\begin{aligned} e(c)&=P(A=1\mid C=c),\\ \pi_a(c)&=P(R=1\mid A=a,C=c),\\ q_a(c)&=P(A=a,R=1\mid C=c). \end{aligned}

Thus q_1=e\pi_1 and q_0=(1-e)\pi_0. MAR identifies \mu_{ag} from the observed outcomes within group a. Positive q_a provides data for each required group/covariate combination. Stronger lower bounds will support the simple variance and remainder argument below.

The missingness assumption is known in this generator. In real assays, failure to detect a low-abundance feature can depend on its unobserved value; that mechanism is not automatically MAR. Feature-specific non-detection, missing covariates, repeated specimens, and perfect group–batch confounding each require a revised analysis. See Missingness, selection, and the observed-data target and Units of analysis, replication, and study design.

4 Generate and audit the study

The complete R functions use base R. Start from a fresh R session; each random operation has an explicit seed.

Code

source("../examples/adjusted-expression.R")
config <- ExpressionStudyConfig()
study <- GenerateExpressionStudy(config$n, config$genes, config$data_seed)
aligned <- AlignExpressionStudy(study$assay, study$metadata)
d <- aligned$metadata
c(features = nrow(study$assay), specimens = ncol(study$assay),
  people = length(unique(d$id)))

 features specimens    people 
       24       800       800

Code

with(d, table(group = a, assay_observed = r))

     assay_observed
group   0   1
    0 132 221
    1 236 211

The assay has features in rows and person IDs in columns. Metadata are intentionally shuffled. AlignExpressionStudy() checks unique IDs, an exact match of specimen sets, and explicit column alignment. It rejects nonfinite observed outcomes and nonmissing values on rows declared unobserved. No complete-case sample exclusion is used to define the covariate distribution. Every person remains in the analysis.

The generator uses Z\sim\operatorname{Uniform}(-1,1) and B\sim\operatorname{Bernoulli}(1/2) independently, with

\begin{aligned} \operatorname{logit}e(C)&=-0.1+Z+0.5B,\\ \operatorname{logit}\pi_A(C)&=0.6-0.6A +0.5Z-0.3B. \end{aligned}

Here \operatorname{logit}(p)=\log\{p/(1-p)\}; applying its inverse expit turns each linear predictor into a probability. This is a convenient factorization of the joint distribution, not a claim that technical batch causes group status. The conditional outcome mean is

\mu_{ag}(Z,B)=b_g+c_gZ+d_g(Z^2-1/3) +t_gB+a(\delta_g+h_gZ).

A shared bounded person-level noise component and independent bounded feature noise induce correlated features. Their conditional means are zero. Since E[Z]=0, the population standardized contrast is \psi_g=\delta_g. The first six features have nonzero \delta_g; the remaining eighteen have zero standardized contrasts, even though their conditional group contrasts vary with Z.

The generator returns these oracle values in a separate truth table. EstimateExpressionContrast() has no truth argument and does not receive the complete missing outcomes.

Code

g <- 12; z_grid <- seq(-1, 1, length.out = 201)
conditional <- study$truth$psi[g] + study$truth$heterogeneity[g] * z_grid
plot(z_grid, conditional, type = "n", xlab = "Covariate Z",
     ylab = "Conditional group contrast (log2 units)")
for (side in list(z_grid <= 0, z_grid >= 0)) {
  polygon(c(z_grid[side], rev(z_grid[side])),
          c(conditional[side], rep(0, sum(side))),
          col = if (all(z_grid[side] <= 0)) "#D55E0033" else "#2F80C833", border = NA)
}
lines(z_grid, conditional, lwd = 2, col = "#2F80C8")
abline(h = study$truth$psi[g], lty = 2)
legend("topleft", c("Conditional contrast", "Population standardized contrast"),
       col = c("#2F80C8", "black"), lty = c(1, 2), bty = "n", cex = 0.8)

A conditional contrast rises from a negative value at Z=-1 to a positive value at Z=1, crossing the horizontal standardized contrast of zero. — Figure 1: A zero standardized contrast can coexist with nonzero conditional contrasts. For synthetic gene_12, the positive and negative areas cancel under the uniform population distribution of Z. This is an association, not a causal curve.

5 Derive the contribution before fitting a model

The influence function in the unrestricted observed-data model, under identification and square integrability, is

\phi_g(W)=\mu_{1g}(C)-\mu_{0g}(C)-\psi_g +\frac{AR}{q_1(C)}\{Y_g-\mu_{1g}(C)\} -\frac{(1-A)R}{q_0(C)}\{Y_g-\mu_{0g}(C)\}.

The first line captures variation in the population’s covariate mix. The second uses observed group-1 residuals; the third subtracts observed group-0 residuals. Conditional centering gives mean zero. In an inner product with an outcome-law score, the joint observation weight cancels the probability of observing that group’s outcome. Scores changing group assignment or observation alone have zero derivative for this functional. Covariate-distribution scores pair with the first line. This follows the derivation in From an influence function to an estimator separately for the two standardized means.

The uncentered contribution starts with the two predicted means for everyone’s C. An observed group-1 person then adds their weighted prediction error; an observed group-0 person subtracts theirs. Someone with a missing assay contributes only the predicted contrast, but remains part of the target covariate population. For example, suppose a held-out person’s predictions are 6.5 and 6.0. If they are observed in group 1 with outcome 7.0 and fitted joint probability 0.25, their contribution is 0.5+(7.0-6.5)/0.25=2.5. That is a weighted statistical contribution, not their measured or causal individual contrast.

Why the joint probability? Among people with the same C, only a fraction q_1(C) both belong to group 1 and have an observed assay. Multiplying by AR/q_1(C) makes that subgroup’s residual correction represent its conditional population average. Weighting by the group probability alone would leave assay selection uncorrected; weighting by observation alone would leave the unequal group distribution of C uncorrected.

For held-out person i, estimate all nuisances on other people, construct the uncentered contribution \hat D_{ig}=\hat\phi_{ig}+\hat\psi_g, and average it.

Algorithm: cross-fitted standardized expression contrast

Inputs: the aligned observed assay, all people’s covariates and observation indicators, and prespecified fitting choices.

Randomly assign independent people to five folds. All features for a person stay together.
On each training set, fit group membership using a ~ z + batch and observation using r ~ a + z + batch.
Fit a quadratic outcome regression separately in each group, using only its observed training outcomes. All genes use the same design, but have separate coefficients.
Predict both group means, both observation probabilities, and the group probability on held-out covariates.
Add residual corrections only where the corresponding outcome is observed, then combine all held-out contributions.

Output: one contribution per person and feature, and their feature-wise means. Uncertainty additionally requires the remainder and sampling conditions below.

The residual correction uses q_a, not just the group probability or just the assay observation probability. It never multiplies an NA residual by zero. In code, the group-1 update is:

Code

observed_1 <- a == 1 & r == 1
score[observed_1, ] <- score[observed_1, , drop = FALSE] +
  sweep(y[observed_1, , drop = FALSE] - m1[observed_1, , drop = FALSE],
        1, q1[observed_1], "/")

This is the same code used in the downloadable function. For repeated specimens, splitting arbitrary columns would leak person-level information and the independent-person SE below would be inappropriate.

6 Name the remainder and the inferential conditions

Hold one training fit fixed. Let \Delta\mu_{ag}=\hat\mu_{ag}-\mu_{ag}. Its exact population bias is

P\hat D_g-\psi_g=E\!\left[\frac{\hat q_1-q_1}{\hat q_1}\Delta\mu_{1g}\right] -E\!\left[\frac{\hat q_0-q_0}{\hat q_0}\Delta\mu_{0g}\right].

If each fitted joint probability is at least a fixed \kappa>0, its magnitude is bounded by

\kappa^{-1}\big\{\|\hat q_1-q_1\|_2\|\Delta\mu_{1g}\|_2 +\|\hat q_0-q_0\|_2\|\Delta\mu_{0g}\|_2\big\}.

Both correct outcome regressions can make this bias vanish with controlled, incorrect weights. Correct joint probabilities can do so with controlled, incorrect outcome regressions. These are consistency statements. They do not automatically justify the same estimated influence-function SE under misspecification.

For the efficient first-order argument used here, require i.i.d. people; fixed fold count with nonvanishing fold proportions; MAR; true and fitted joint probabilities bounded below; bounded outcomes and fitted means; L2 consistency of all nuisances; each displayed product o_p(n^{-1/2}); and finite nonzero influence variance. The fold-wise conditional proof in From an influence function to an estimator controls the empirical remainder. The bias bound controls the other remainder. Then

\begin{aligned} \hat\psi_g-\psi_g&=(P_n-P)\phi_g+o_p(n^{-1/2}),\\ \widehat{\operatorname{SE}}_g&=\frac{\operatorname{SD}(\hat D_{ig}-\hat\psi_g)}{\sqrt n}. \end{aligned}

The default generator and correctly specified finite-dimensional fits provide a transparent setting for these conditions. This does not demonstrate that those nuisance models are adequate in a real cohort. Nor does the fixed-24-feature argument establish uniform validity when the number of features grows rapidly with sample size.

The SE formula uses variation in the constructed contributions across all sampled people, including those with missing assays. It is not the raw expression SD divided by the square root of the observed assay count. Cross-fitting contributions share training fits; their dependence becomes negligible at first order only through the displayed expansion. The proof, rather than treating fitted scores as new independent data, justifies this variance calculation.

7 Fit, diagnose, and compare the targets

Code

fit <- EstimateExpressionContrast(study$assay, study$metadata,
  folds = config$folds, seed = config$fold_seed,
  probability_bounds = config$probability_bounds,
  mean_bounds = config$mean_bounds, alpha = config$alpha)
fit$diagnostics

    n observed group0 group1 observed0 observed1     e_min     e_max    pi_min
1 800      432    353    447       221       211 0.2602384 0.8299704 0.2877082
    pi_max      q_min max_observed_weight probability_clipped mean_clipped
1 0.801157 0.08226094            11.23615                   0            0
  warning_count
1             0

The teaching bounds constrain fitted probabilities to [0.1,0.9] and fitted means to [0,16]. The true functions lie inside these prespecified bounds. Joint probabilities are products, so their lower fitted bound is 0.01. Clipping is recorded; it does not establish population overlap or MAR. Training fits stop on inadequate observations, singular outcome designs, nonconvergence, or invalid results; warnings are retained.

Code

boxplot(e ~ a, data = fit$predictions, names = c("Group 0", "Group 1"),
        ylim = c(0, 1), col = c("#D55E0033", "#2F80C833"),
        xlab = "Recorded group", ylab = "Fitted P(A=1 | C)")
abline(h = config$probability_bounds, lty = 2, col = "gray40")

Boxplots of fitted group-1 probabilities by observed group, with substantial overlap and dashed lines at 0.1 and 0.9. — Figure 2: Out-of-fold fitted group probabilities for all 800 synthetic people. The distributions overlap; dashed lines are the prespecified clipping bounds, not a test of population positivity. Observation probabilities and joint weights need separate inspection in the diagnostics.

Join truth only after fitting, for the synthetic comparison. The complete-case crude difference compares different covariate distributions in the two groups and a selected set of outcomes; it does not target \psi_g in general.

Code

results <- fit$results
results$truth <- study$truth$psi[match(results$gene, study$truth$gene)]
results[c(1:6, 12), c("gene", "truth", "complete_case_crude", "estimate", "se")]

      gene truth complete_case_crude    estimate         se
1  gene_01  0.70          0.97462630  0.66973970 0.04966932
2  gene_02 -0.65         -0.49572086 -0.64003498 0.05076546
3  gene_03  0.50          0.77493285  0.46298809 0.04820242
4  gene_04 -0.45         -0.03338629 -0.41017725 0.04694681
5  gene_05  0.35          0.46473893  0.32263596 0.04835757
6  gene_06 -0.30          0.08724268 -0.24268797 0.04836003
12 gene_12  0.00          0.40556590  0.01270766 0.05177327

8 Address the family of claims

Use the two-sided asymptotic Wald p-value for H_{0g}:\psi_g=0. All 24 tests belong to the prespecified family. Shared noise and nuisance fitting induce dependence. We use Benjamini–Yekutieli (BY) as the primary FDR adjustment because its multiplicity result permits arbitrary dependence among valid null p-values. This does not remove the need for valid individual tests; here those are asymptotic approximations.

The calculation takes all 24 raw p-values, sorts them, and applies the BH step-up rule at 0.05/H_{24}, where H_{24}\approx3.78. It then maps the decisions back to feature IDs. This is the BY algorithm in Multiple testing when the tests are not independent; selecting just the displayed or smallest-p-value features first would change the family and the guarantee. Each null concerns the standardized contrast, so a gene with varying conditional contrasts can still be a true null in this family.

With a fixed family and jointly adequate marginal approximations, the argument supports asymptotic FDR control at 0.05. It is not an exact finite-sample guarantee. BH can be less conservative under appropriate dependence conditions, but positive feature correlations alone do not establish the required condition for these two-sided tests. See the Multiple testing when the tests are not independent page and R’s adjustment documentation.

Code

results[, c("gene", "estimate", "lower", "upper", "p_BY", "reject_BY")]

      gene      estimate        lower       upper         p_BY reject_BY
1  gene_01  0.6697397005  0.572389627  0.76708977 1.761300e-39      TRUE
2  gene_02 -0.6400349785 -0.739533457 -0.54053650 8.678789e-35      TRUE
3  gene_03  0.4629880883  0.368513082  0.55746309 2.298770e-20      TRUE
4  gene_04 -0.4101772488 -0.502191296 -0.31816320 5.420311e-17      TRUE
5  gene_05  0.3226359563  0.227856868  0.41741504 4.577265e-10      TRUE
6  gene_06 -0.2426879695 -0.337471894 -0.14790404 7.871356e-06      TRUE
7  gene_07  0.0687310636 -0.020289329  0.15775146 1.000000e+00     FALSE
8  gene_08  0.0418306533 -0.052756760  0.13641807 1.000000e+00     FALSE
9  gene_09  0.0032866151 -0.082295057  0.08886829 1.000000e+00     FALSE
10 gene_10  0.0286631291 -0.069126115  0.12645237 1.000000e+00     FALSE
11 gene_11  0.0234872549 -0.072376394  0.11935090 1.000000e+00     FALSE
12 gene_12  0.0127076633 -0.088766073  0.11418140 1.000000e+00     FALSE
13 gene_13  0.0007318202 -0.096101815  0.09756545 1.000000e+00     FALSE
14 gene_14  0.0073960085 -0.088834690  0.10362671 1.000000e+00     FALSE
15 gene_15  0.0280877638 -0.068984522  0.12516005 1.000000e+00     FALSE
16 gene_16  0.0686265837 -0.024782556  0.16203572 1.000000e+00     FALSE
17 gene_17  0.0149326366 -0.084105813  0.11397109 1.000000e+00     FALSE
18 gene_18  0.0831296867 -0.015230569  0.18148994 9.830168e-01     FALSE
19 gene_19  0.0348696613 -0.063540210  0.13327953 1.000000e+00     FALSE
20 gene_20  0.0798578970 -0.010199049  0.16991484 9.312820e-01     FALSE
21 gene_21  0.0215018226 -0.068730377  0.11173402 1.000000e+00     FALSE
22 gene_22 -0.0397497000 -0.134430897  0.05493150 1.000000e+00     FALSE
23 gene_23  0.0807160265 -0.008977473  0.17040953 9.312820e-01     FALSE
24 gene_24 -0.0117349386 -0.102335702  0.07886582 1.000000e+00     FALSE

The confidence intervals remain pointwise. BY-adjusted p-values do not turn them into simultaneous or selection-adjusted intervals. A reported discovery list also has no known realized false-discovery proportion in real data.

Code

show <- c(1:8, 12, 16, 20, 24); a <- results[show, ]; position <- rev(seq_len(nrow(a)))
par(mar = c(4, 5, 1, 1))
plot(a$estimate, position, xlim = range(c(a$lower, a$upper, a$truth)),
     yaxt = "n", ylab = "", xlab = "Standardized contrast (log2 units)",
     pch = 16, col = ifelse(a$reject_BY, "#2F80C8", "gray35"))
segments(a$lower, position, a$upper, position, col = "gray50")
points(a$estimate, position, pch = 16, col = ifelse(a$reject_BY, "#2F80C8", "gray35"))
points(a$truth, position, pch = 4, col = "#D55E00", cex = 1.1)
axis(2, at = position, labels = a$gene, las = 1, cex.axis = 0.8)
abline(v = 0, lty = 2, col = "gray60")

A forest plot compares adjusted estimates and pointwise confidence intervals with known synthetic contrasts for six nonnull and six null features. — Figure 3: Twelve prespecified display features from the 24-feature family. Points and bars are estimates and pointwise approximate 95% intervals; crosses are known simulation truth. Blue indicates BY rejection. The intervals are not simultaneous or adjusted for selection.

9 Challenge the analysis in repeated samples

Evaluate all features on the same synthetic datasets under four prespecified fitting choices: both model families correct; only outcome models correct; only weight models correct; both wrong. Incorrect models are intercept-only fits. Repeat at n=400 and 1200. The two singly correct cases investigate double-robust consistency; the efficient-SE theorem above requires all nuisance functions consistent and does not promise interval calibration in those cases.

Code

records <- EvaluateExpressionStudy(B = config$simulation_replicates,
  ns = config$simulation_n, genes = config$genes, folds = config$folds,
  seed = config$simulation_seed, alpha = config$alpha,
  probability_bounds = config$probability_bounds, mean_bounds = config$mean_bounds)
performance <- SummarizeExpressionEvaluation(records)

Code

performance[, c("n", "scenario", "attempted", "successful", "sentinel_bias", "bias_mcse")]

     n             scenario attempted successful sentinel_bias   bias_mcse
1  400         Both correct       100        100  -0.005799467 0.006943310
2  400 Outcome only correct       100        100  -0.004718295 0.006827738
3  400 Weights only correct       100        100  -0.032840763 0.010586325
4  400           Both wrong       100        100   0.371145986 0.009351854
5 1200         Both correct       100        100  -0.003759116 0.003970181
6 1200 Outcome only correct       100        100  -0.003351759 0.003935657
7 1200 Weights only correct       100        100  -0.018433270 0.006295249
8 1200           Both wrong       100        100   0.372832048 0.006000246

The sentinel is prespecified null feature gene_12. Bias here is the average signed error across independent datasets; its MCSE uses the replicate-level SD divided by the square root of successful replicate count. A small bias in one finite experiment does not prove consistency.

Code

scenarios <- unique(performance$scenario)
colors <- c("#2F80C8", "#009E73", "#7B61A8", "#D55E00")
layout(matrix(1:2, ncol = 1), heights = c(1, 4))
par(mar = c(0, 0, 0, 0)); plot.new()
legend("center", scenarios, col = colors, pch = 15:18, lty = 1,
       ncol = 2, bty = "n", cex = 0.8)
par(mar = c(4, 4, 0.5, 1))
plot(NA, xlim = range(config$simulation_n) * c(0.9, 1.08),
     ylim = range(c(0, performance$sentinel_bias + 2 * performance$bias_mcse,
                    performance$sentinel_bias - 2 * performance$bias_mcse)),
     xlab = "Independent people n", ylab = "Mean error for gene_12")
abline(h = 0, lty = 2, col = "gray50")
for (j in seq_along(scenarios)) {
  a <- performance[performance$scenario == scenarios[j], ]
  x <- a$n * (1 + (j - 2.5) * 0.025)
  lines(x, a$sentinel_bias, type = "b", pch = j + 14, col = colors[j])
  arrows(x, a$sentinel_bias - 2 * a$bias_mcse, x, a$sentinel_bias + 2 * a$bias_mcse,
         angle = 90, code = 3, length = 0.04, col = colors[j])
}

Bias curves for four nuisance-model scenarios at 400 and 1200 people, with Monte Carlo error bars and a zero-bias reference line. — Figure 4: Mean estimation error for the prespecified null gene_12. Bars are plus/minus two MCSEs for the simulated bias, not SEs of an individual gene estimate. Each scenario uses the same datasets and folds. Small horizontal offsets separate scenarios at each sample size. Only two sample sizes are shown; the asymptotic claims come from the remainder analysis.

For interval and multiplicity summaries, first compare the fully specified setup with the deliberately wrong setup:

Code

performance[performance$scenario %in% c("Both correct", "Both wrong"),
  c("n", "scenario", "coverage", "coverage_mcse", "fdr", "fdr_mcse", "true_positive_fraction")]

     n     scenario   coverage coverage_mcse         fdr    fdr_mcse
1  400 Both correct 0.95833333   0.005772288 0.005595238 0.003296084
4  400   Both wrong 0.12833333   0.015036898 0.737973926 0.017038970
5 1200 Both correct 0.95125000   0.008185880 0.010714286 0.004439868
8 1200   Both wrong 0.01666667   0.003792085 0.803885440 0.002060112
  true_positive_fraction
1              0.9783333
4              0.6800000
5              1.0000000
8              0.7216667

Within each replicate, coverage is the fraction of the 24 fixed targets covered, and FDP is false rejections divided by all rejections (zero when no tests reject). We average those replicate summaries. Their MCSEs use variation across independent datasets; genes within a dataset are not independent simulation replicates. A plug-in MCSE can be zero if every observed FDP is zero; that does not prove true FDR is zero.

Each attempt is retained with seeds, success status, and any error. Metrics are conditional on successful complete-family fits, so report success counts alongside them; do not interpret conditional performance as unconditional reliability when fits fail. One hundred replicates make this a runnable teaching experiment, with visibly reported Monte Carlo uncertainty. A research benchmark should plan its replicate count from the precision needed.

10 Reproduce the inputs, run, and report

Download the complete project, which includes the R functions, editable configuration, runner, observed synthetic assay/metadata, truth stored separately for evaluation, diagnostics, contributions, all simulation attempts, manifests, software information, Quarto report source, and rendered report. You can also read the standalone report.

The archive’s results/ is the supplied reference run. From the unpacked project root, rebuild into a new directory so the reference run remains available:

Code

Rscript run.R rerun-results
quarto render report.qmd -P results_dir:rerun-results

The runner refuses to overwrite a populated output directory. It records input/output identity and code/configuration checksums for local change tracking. Analysis needs base R; rendering the report additionally needs a working Quarto R/knitr installation. The report reads explicit result files and has freezing disabled. It can be regenerated without an interactive workspace.

A methods paragraph should state the scale, population, independent unit, target distribution, missingness assumptions, nuisance models, folds, clipping, family, uncertainty calculation, and inferential limits. For example, the following is generated from this run:

Code

s <- results[1, ]; a <- fit$diagnostics
cat(sprintf(paste0(
  "In this synthetic study, %d independently sampled people supplied covariates, ",
  "and %d supplied the whole expression assay. We standardized group contrasts ",
  "to the source-population covariate distribution using all people. ",
  "For gene_01, the cross-fitted augmented estimate was %.3f log2 units ",
  "(SE %.3f; pointwise approximate 95%% interval %.3f to %.3f). ",
  "The prespecified family contained %d features; BY at 0.05 rejected %d nulls. ",
  "Fitted probabilities and means required clipping for %.1f%% and %.1f%% ",
  "of their predictions. These are standardized associations under the stated ",
  "MAR, overlap, sampling, and nuisance-rate conditions; they do not establish ",
  "causal or absolute-abundance effects.\n\n"),
  a$n, a$observed, s$estimate, s$se, s$lower, s$upper,
  nrow(results), sum(results$reject_BY), 100 * a$probability_clipped, 100 * a$mean_clipped))

In this synthetic study, 800 independently sampled people supplied covariates, and 432 supplied the whole expression assay. We standardized group contrasts to the source-population covariate distribution using all people. For gene_01, the cross-fitted augmented estimate was 0.670 log2 units (SE 0.050; pointwise approximate 95% interval 0.572 to 0.767). The prespecified family contained 24 features; BY at 0.05 rejected 6 nulls. Fitted probabilities and means required clipping for 0.0% and 0.0% of their predictions. These are standardized associations under the stated MAR, overlap, sampling, and nuisance-rate conditions; they do not establish causal or absolute-abundance effects.

The scientific interpretation is bounded by the target and assumptions, even when execution is reproducible and the synthetic answer is close to truth. Return to Reporting a statistical methods study for the full claim/evidence checklist.

11 Exercises

Try each checkpoint before opening its hint or solution.

Checkpoint · CASE-01 — Specify and identify the expression contrast

Task type: Interpretation and derivation.

Write the independent observation, measurement scale, population distribution, estimand, and tested family for this case. Derive why the generator gives \psi_g=\delta_g. Explain how a feature can have \psi_g=0 with conditional contrasts that vary with Z, and state the missingness/support assumptions needed to identify the full-population contrast.

Ready to continue when: you can distinguish the standardized target from a complete-case crude contrast, a constant regression coefficient, and a causal effect.

Transfer: If the target covariate distribution were that of group 1, would the same delta remain the target? Explain what changes in estimation and uncertainty.

Hint

Average the conditional contrast over the specified population, then identify which outcome laws are actually observed.

Solution outline — open after trying

One independent observation is a person’s covariates, group, observation indicator, and observed assay vector. The outcome is calibrated log2 intensity, and the family contains all 24 feature-wise standardized-contrast nulls. The conditional contrast is \delta_g+h_gZ; since the population has E[Z]=0, its average is \delta_g. A zero average need not imply zero contrast at every covariate value. MAR within group/covariates identifies the outcome regressions from observed outcomes; positive group and observation probabilities are needed throughout the target support, and covariates must be observed for the target population. None of these statements alone gives a causal or absolute-abundance interpretation.

Checkpoint · CASE-02 — Reconstruct the estimator and its remainder

Task type: Derivation and implementation.

Derive the conditional expectation of each residual correction, then the two-term exact bias displayed in this page. Explain why the weight uses the joint probability q_a. Reconstruct the fold procedure without accessing simulation truth or unobserved outcomes, and identify the conditions needed before using the estimated influence-function SE.

Ready to continue when: you can connect every correction and condition to an observed-data quantity and explain what double-robust consistency does and does not establish.

Transfer: If observation probabilities depend on the unobserved expression even after conditioning on group and C, which step of the argument fails?

Hint

Hold a training fit fixed. For group a, the expected observed residual is q_a times the difference between the true and fitted conditional mean.

Solution outline — open after trying

The group-a correction has fresh-observation expectation E[q_a(\mu_{ag}-\hat\mu_{ag})/\hat q_a]. Combining it with each plug-in mean yields the displayed signed product remainders. Weighting only by group or observation probability fails to cancel the joint observation process. Fit all nuisance functions on other people, evaluate both group predictions on all held-out covariates, and correct residuals only on observed rows of the matching group. For efficient inference, check nuisance consistency, lower probability bounds, product little-o rates, the held-out empirical term, finite variance and sampling/moment conditions. A wrong but controlled nuisance can permit consistency while invalidating that efficient-SE argument.

Checkpoint · CASE-03 — Interpret discoveries and uncertainty correctly

Task type: Interpretation and computation.

Explain the difference between a pointwise 95% interval, BY FDR control for the 24-feature family, a realized FDP in one simulation, and MCSE for estimated FDR. State the assumptions left after selecting BY and explain why you cannot treat 24 genes in each of 100 datasets as 2,400 independent replicates.

Ready to continue when: you can label each uncertainty quantity and error criterion, including the conditions multiplicity adjustment cannot repair.

Transfer: If you display only the twelve prespecified plot features, does the tested family shrink to twelve? What if those twelve were chosen after inspecting significance?

Hint

Ask whether the expectation is over people, repeated datasets, or the feature family.

Solution outline — open after trying

A pointwise interval targets one feature’s population contrast across samples. BY targets the expected FDP across repeated analyses of a stated family, allowing arbitrary p-value dependence but requiring valid null p-values; the Wald p-values here are asymptotically justified under the estimator’s conditions. One simulated FDP is known only because truth is supplied. Its average estimates FDR; its MCSE is the SD of replicate FDPs divided by the square root of the number of independent successful datasets. Within-dataset genes share people, noise, and fits. BY does not create simultaneous confidence intervals or repair missingness, target, or sampling errors.

Checkpoint · CASE-04 — Reproduce and challenge the complete analysis

Task type: Execution and diagnosis.

Unpack the project, regenerate it in a new results directory, render its report, and compare target estimates and simulation records with the reference run. Audit the input alignment, missing-outcome conventions, seeds, code/configuration identities, warnings, clipping, and attempted/successful counts. Before examining the four scenario summaries, predict which cases support consistency and which support the efficient-SE theorem.

Ready to continue when: you can rebuild the result, trace a claim to its inputs and assumptions, and diagnose a failed condition without silently discarding failures.

Transfer: Change only the number of folds or one model class, record a new run, and distinguish stochastic changes from a changed target or invalid theorem.

Hint

Use a fresh session and retain the reference directory. Derive the exact bias prediction before looking at the simulation plot.

Solution outline — open after trying

The documented runner sources explicit functions and configuration, writes observed inputs separately from oracle truth, and records results, diagnostics, attempts, manifests, and software information. With the tested environment and seeds the statistical outputs reproduce; session metadata can differ across environments. Correct outcomes or correct joint weights can remove population bias under the stated controls. The fully consistent, sufficiently fast nuisance case supports the efficient-SE argument; single-model correctness alone does not. Both wrong leaves a generally nonzero remainder. Inspect failures and their denominators rather than assuming missing rows succeeded.

Additional practice (optional)

Extend the generator to feature-specific missingness and derive the corresponding observation probabilities; do not silently reuse the whole-assay indicator.
Use an external target covariate sample. Derive the new target and account for its sampling variation instead of treating the current influence function as unchanged.
Compare the four scenarios across more sample sizes and choose replicate counts from a prespecified MCSE goal. Treat singly correct nuisance cases as a separate inference problem.
Construct a simultaneous interval procedure for the prespecified family and explain how its claim differs from the plotted pointwise intervals.

12 Reading

Kennedy, review of influence-function estimation: standardized means, estimator construction, and remainder analysis.
Benjamini and Yekutieli (2001), FDR under dependency, and the official R p.adjust documentation.

13 Continue through the guide

Return to the Omics analysis problems overview, compare with the Differential expression and abundance page, or assess the four case-study tasks in the Completion checkpoints page.

--- title: "An end-to-end expression analysis" subtitle: "From a standardized contrast to semiparametric inference and a reproducible report" guide-status: "Available" --- **Before reading:** Read [Preprocessing and the quantity being estimated](../data-and-measurement/preprocessing-and-estimands.qmd), [From a research question to an identified estimand](../questions-and-design/question-to-estimand.qmd), and [From an influence function to an estimator](../semiparametric-methods/influence-function-to-estimator.qmd). The target and data-audit sections can be read before completing all the proofs. ## One research argument, from measurements to a reported contrast We will analyze 24 synthetic expression features from independently sampled people. Group membership is associated with a recorded covariate and technical batch. Some whole-assay measurements are missing. The question is: **how do the groups differ on the supplied log2 expression scale after standardizing to a common covariate distribution?** This extended worked application follows one argument: define the target, audit the data, construct its estimator and uncertainty, then assess the reported family of claims. Simulation checks specific predictions of that argument. The [From an influence function to an estimator](../semiparametric-methods/influence-function-to-estimator.qmd) page supplies the proof pattern; here we combine group standardization with missing outcomes. **Everything in this example is synthetic.** Feature IDs are invented. The generator supplies known truth for evaluation, while the estimator receives only the observed assay and metadata. No clinical or biological finding is being asserted. ## Define the population, scale, and target {#target-specification} Let $A$ be group (0 or 1), $C=(Z,B)$ the recorded covariate and batch, and $Y_g$ a feature's calibrated log2 intensity. Define $$ \begin{aligned} \mu_{ag}(c)&=E[Y_g\mid A=a,C=c],\\ \psi_g&=E_C[\mu_{1g}(C)-\mu_{0g}(C)]. \end{aligned} $$ The outer distribution is that of **all eligible people in the source population**, estimated using everyone's recorded $C$, including people with missing assays. We do not condition on the realized covariate values as a fixed population. Variation in the sampled covariate mix therefore contributes to the influence function. ::: {.table-responsive} | Decision | Specification in this case | |---|---| | Independent unit | One person, one specimen, one assay occasion | | Population | The source population sampled by the synthetic generator | | Measurement | Supplied calibrated log2 intensity; no learned normalization | | Contrast | Mean conditional group difference standardized over the source population's $C$ distribution | | Missingness | One recorded indicator $R$ for a missing whole assay; all covariates remain observed | | Family | All 24 prespecified feature hypotheses $H_{0g}:\psi_g=0$ | | Interpretation | A standardized association; no causal interpretation is assumed | ::: A log2 contrast exponentiates to a ratio of geometric means of the two standardized intensity distributions. It is not a ratio of arithmetic mean expression or an absolute molecules-per-cell contrast. The [Preprocessing and the quantity being estimated](../data-and-measurement/preprocessing-and-estimands.qmd) page explains those distinctions. Real count data would need their own measurement and normalization argument before this analysis could be reused. Group contrasts can vary with $C$. We do not impose a constant coefficient or a common regression slope across groups. Nor does adjustment alone turn this association into an intervention effect. A causal interpretation would need a well-defined intervention and additional identification assumptions. ## Reconstruct the observation process We observe $W=(C,A,R,RY_1,\ldots,RY_G)$. Algebraically $RY_g$ contributes zero when $R=0$; the stored assay uses `NA` for those outcomes. Treating those outcomes as measured zeros would change the analysis. Assume $E[|Y_g|]<\infty$, outcome missingness at random given $(A,C)$, and positive group and observation probabilities on the target support: $$ \begin{aligned} e(c)&=P(A=1\mid C=c),\\ \pi_a(c)&=P(R=1\mid A=a,C=c),\\ q_a(c)&=P(A=a,R=1\mid C=c). \end{aligned} $$ Thus $q_1=e\pi_1$ and $q_0=(1-e)\pi_0$. MAR identifies $\mu_{ag}$ from the observed outcomes within group $a$. Positive $q_a$ provides data for each required group/covariate combination. Stronger lower bounds will support the simple variance and remainder argument below. The missingness assumption is known in this generator. In real assays, failure to detect a low-abundance feature can depend on its unobserved value; that mechanism is not automatically MAR. Feature-specific non-detection, missing covariates, repeated specimens, and perfect group–batch confounding each require a revised analysis. See [Missingness, selection, and the observed-data target](../questions-and-design/missingness-and-selection.qmd) and [Units of analysis, replication, and study design](../questions-and-design/units-and-study-design.qmd). ## Generate and audit the study The complete [R functions](../examples/adjusted-expression.R) use base R. Start from a fresh R session; each random operation has an explicit seed. ```{r} source("../examples/adjusted-expression.R") config <- ExpressionStudyConfig() study <- GenerateExpressionStudy(config$n, config$genes, config$data_seed) aligned <- AlignExpressionStudy(study$assay, study$metadata) d <- aligned$metadata c(features = nrow(study$assay), specimens = ncol(study$assay), people = length(unique(d$id))) with(d, table(group = a, assay_observed = r)) ``` The assay has features in rows and person IDs in columns. Metadata are intentionally shuffled. `AlignExpressionStudy()` checks unique IDs, an exact match of specimen sets, and explicit column alignment. It rejects nonfinite observed outcomes and nonmissing values on rows declared unobserved. No complete-case sample exclusion is used to define the covariate distribution. Every person remains in the analysis. The generator uses $Z\sim\operatorname{Uniform}(-1,1)$ and $B\sim\operatorname{Bernoulli}(1/2)$ independently, with $$ \begin{aligned} \operatorname{logit}e(C)&=-0.1+Z+0.5B,\\ \operatorname{logit}\pi_A(C)&=0.6-0.6A +0.5Z-0.3B. \end{aligned} $$ Here $\operatorname{logit}(p)=\log\{p/(1-p)\}$; applying its inverse expit turns each linear predictor into a probability. This is a convenient factorization of the joint distribution, not a claim that technical batch causes group status. The conditional outcome mean is $$ \mu_{ag}(Z,B)=b_g+c_gZ+d_g(Z^2-1/3) +t_gB+a(\delta_g+h_gZ). $$ A shared bounded person-level noise component and independent bounded feature noise induce correlated features. Their conditional means are zero. Since $E[Z]=0$, the population standardized contrast is $\psi_g=\delta_g$. The first six features have nonzero $\delta_g$; the remaining eighteen have zero standardized contrasts, even though their conditional group contrasts vary with $Z$. The generator returns these oracle values in a separate `truth` table. `EstimateExpressionContrast()` has no truth argument and does not receive the complete missing outcomes. ```{r} #| label: fig-case-conditional-target #| fig-cap: "A zero standardized contrast can coexist with nonzero conditional contrasts. For synthetic gene_12, the positive and negative areas cancel under the uniform population distribution of Z. This is an association, not a causal curve." #| fig-alt: "A conditional contrast rises from a negative value at Z=-1 to a positive value at Z=1, crossing the horizontal standardized contrast of zero." g <- 12; z_grid <- seq(-1, 1, length.out = 201) conditional <- study$truth$psi[g] + study$truth$heterogeneity[g] * z_grid plot(z_grid, conditional, type = "n", xlab = "Covariate Z", ylab = "Conditional group contrast (log2 units)") for (side in list(z_grid <= 0, z_grid >= 0)) { polygon(c(z_grid[side], rev(z_grid[side])), c(conditional[side], rep(0, sum(side))), col = if (all(z_grid[side] <= 0)) "#D55E0033" else "#2F80C833", border = NA) } lines(z_grid, conditional, lwd = 2, col = "#2F80C8") abline(h = study$truth$psi[g], lty = 2) legend("topleft", c("Conditional contrast", "Population standardized contrast"), col = c("#2F80C8", "black"), lty = c(1, 2), bty = "n", cex = 0.8) ``` ## Derive the contribution before fitting a model The influence function in the unrestricted observed-data model, under identification and square integrability, is $$ \phi_g(W)=\mu_{1g}(C)-\mu_{0g}(C)-\psi_g +\frac{AR}{q_1(C)}\{Y_g-\mu_{1g}(C)\} -\frac{(1-A)R}{q_0(C)}\{Y_g-\mu_{0g}(C)\}. $$ The first line captures variation in the population's covariate mix. The second uses observed group-1 residuals; the third subtracts observed group-0 residuals. Conditional centering gives mean zero. In an inner product with an outcome-law score, the joint observation weight cancels the probability of observing that group's outcome. Scores changing group assignment or observation alone have zero derivative for this functional. Covariate-distribution scores pair with the first line. This follows the derivation in [From an influence function to an estimator](../semiparametric-methods/influence-function-to-estimator.qmd#derive-the-observed-data-influence-function) separately for the two standardized means. The **uncentered contribution** starts with the two predicted means for everyone's $C$. An observed group-1 person then adds their weighted prediction error; an observed group-0 person subtracts theirs. Someone with a missing assay contributes only the predicted contrast, but remains part of the target covariate population. For example, suppose a held-out person's predictions are 6.5 and 6.0. If they are observed in group 1 with outcome 7.0 and fitted joint probability 0.25, their contribution is $0.5+(7.0-6.5)/0.25=2.5$. That is a weighted statistical contribution, not their measured or causal individual contrast. Why the joint probability? Among people with the same $C$, only a fraction $q_1(C)$ both belong to group 1 **and** have an observed assay. Multiplying by $AR/q_1(C)$ makes that subgroup's residual correction represent its conditional population average. Weighting by the group probability alone would leave assay selection uncorrected; weighting by observation alone would leave the unequal group distribution of $C$ uncorrected. For held-out person $i$, estimate all nuisances on other people, construct the uncentered contribution $\hat D_{ig}=\hat\phi_{ig}+\hat\psi_g$, and average it. ::: {#expression-crossfit-procedure .callout-note icon=false title="Algorithm: cross-fitted standardized expression contrast"} **Inputs:** the aligned observed assay, all people's covariates and observation indicators, and prespecified fitting choices. 1. Randomly assign independent people to five folds. All features for a person stay together. 2. On each training set, fit group membership using `a ~ z + batch` and observation using <code>r ~ a + z + batch</code>. 3. Fit a quadratic outcome regression separately in each group, using only its observed training outcomes. All genes use the same design, but have separate coefficients. 4. Predict both group means, both observation probabilities, and the group probability on held-out covariates. 5. Add residual corrections only where the corresponding outcome is observed, then combine all held-out contributions. **Output:** one contribution per person and feature, and their feature-wise means. Uncertainty additionally requires the remainder and sampling conditions below. ::: The residual correction uses $q_a$, not just the group probability or just the assay observation probability. It never multiplies an `NA` residual by zero. In code, the group-1 update is: ```{r} #| eval: false observed_1 <- a == 1 & r == 1 score[observed_1, ] <- score[observed_1, , drop = FALSE] + sweep(y[observed_1, , drop = FALSE] - m1[observed_1, , drop = FALSE], 1, q1[observed_1], "/") ``` This is the same code used in the downloadable function. For repeated specimens, splitting arbitrary columns would leak person-level information and the independent-person SE below would be inappropriate. ## Name the remainder and the inferential conditions {#case-remainder} Hold one training fit fixed. Let $\Delta\mu_{ag}=\hat\mu_{ag}-\mu_{ag}$. Its exact population bias is $$ P\hat D_g-\psi_g=E\!\left[\frac{\hat q_1-q_1}{\hat q_1}\Delta\mu_{1g}\right] -E\!\left[\frac{\hat q_0-q_0}{\hat q_0}\Delta\mu_{0g}\right]. $$ If each fitted joint probability is at least a fixed $\kappa>0$, its magnitude is bounded by $$ \kappa^{-1}\big\{\|\hat q_1-q_1\|_2\|\Delta\mu_{1g}\|_2 +\|\hat q_0-q_0\|_2\|\Delta\mu_{0g}\|_2\big\}. $$ Both correct outcome regressions can make this bias vanish with controlled, incorrect weights. Correct joint probabilities can do so with controlled, incorrect outcome regressions. These are **consistency** statements. They do not automatically justify the same estimated influence-function SE under misspecification. For the efficient first-order argument used here, require i.i.d. people; fixed fold count with nonvanishing fold proportions; MAR; true and fitted joint probabilities bounded below; bounded outcomes and fitted means; L2 consistency of all nuisances; each displayed product $o_p(n^{-1/2})$; and finite nonzero influence variance. The fold-wise conditional proof in [From an influence function to an estimator](../semiparametric-methods/influence-function-to-estimator.qmd#checkpoint-sp-12) controls the empirical remainder. The bias bound controls the other remainder. Then $$ \begin{aligned} \hat\psi_g-\psi_g&=(P_n-P)\phi_g+o_p(n^{-1/2}),\\ \widehat{\operatorname{SE}}_g&=\frac{\operatorname{SD}(\hat D_{ig}-\hat\psi_g)}{\sqrt n}. \end{aligned} $$ The default generator and correctly specified finite-dimensional fits provide a transparent setting for these conditions. This does not demonstrate that those nuisance models are adequate in a real cohort. Nor does the fixed-24-feature argument establish uniform validity when the number of features grows rapidly with sample size. The SE formula uses variation in the **constructed contributions across all sampled people**, including those with missing assays. It is not the raw expression SD divided by the square root of the observed assay count. Cross-fitting contributions share training fits; their dependence becomes negligible at first order only through the displayed expansion. The proof, rather than treating fitted scores as new independent data, justifies this variance calculation. ## Fit, diagnose, and compare the targets ```{r} fit <- EstimateExpressionContrast(study$assay, study$metadata, folds = config$folds, seed = config$fold_seed, probability_bounds = config$probability_bounds, mean_bounds = config$mean_bounds, alpha = config$alpha) fit$diagnostics ``` The teaching bounds constrain fitted probabilities to $[0.1,0.9]$ and fitted means to $[0,16]$. The true functions lie inside these prespecified bounds. Joint probabilities are products, so their lower fitted bound is $0.01$. Clipping is recorded; it does not establish population overlap or MAR. Training fits stop on inadequate observations, singular outcome designs, nonconvergence, or invalid results; warnings are retained. ```{r} #| label: fig-case-overlap #| fig-cap: "Out-of-fold fitted group probabilities for all 800 synthetic people. The distributions overlap; dashed lines are the prespecified clipping bounds, not a test of population positivity. Observation probabilities and joint weights need separate inspection in the diagnostics." #| fig-alt: "Boxplots of fitted group-1 probabilities by observed group, with substantial overlap and dashed lines at 0.1 and 0.9." boxplot(e ~ a, data = fit$predictions, names = c("Group 0", "Group 1"), ylim = c(0, 1), col = c("#D55E0033", "#2F80C833"), xlab = "Recorded group", ylab = "Fitted P(A=1 | C)") abline(h = config$probability_bounds, lty = 2, col = "gray40") ``` Join truth only after fitting, for the synthetic comparison. The complete-case crude difference compares different covariate distributions in the two groups and a selected set of outcomes; it does not target $\psi_g$ in general. ```{r} results <- fit$results results$truth <- study$truth$psi[match(results$gene, study$truth$gene)] results[c(1:6, 12), c("gene", "truth", "complete_case_crude", "estimate", "se")] ``` ## Address the family of claims Use the two-sided asymptotic Wald p-value for $H_{0g}:\psi_g=0$. All 24 tests belong to the prespecified family. Shared noise and nuisance fitting induce dependence. We use **Benjamini–Yekutieli (BY)** as the primary FDR adjustment because its multiplicity result permits arbitrary dependence among valid null p-values. This does not remove the need for valid individual tests; here those are asymptotic approximations. The calculation takes all 24 raw p-values, sorts them, and applies the BH step-up rule at $0.05/H_{24}$, where $H_{24}\approx3.78$. It then maps the decisions back to feature IDs. This is the BY algorithm in [Multiple testing when the tests are not independent](multiple-testing.qmd#by-procedure); selecting just the displayed or smallest-p-value features first would change the family and the guarantee. Each null concerns the standardized contrast, so a gene with varying conditional contrasts can still be a true null in this family. With a fixed family and jointly adequate marginal approximations, the argument supports asymptotic FDR control at 0.05. It is not an exact finite-sample guarantee. BH can be less conservative under appropriate dependence conditions, but positive feature correlations alone do not establish the required condition for these two-sided tests. See the [Multiple testing when the tests are not independent](multiple-testing.qmd) page and [R's adjustment documentation](https://stat.ethz.ch/R-manual/R-devel/library/stats/html/p.adjust.html). ```{r} results[, c("gene", "estimate", "lower", "upper", "p_BY", "reject_BY")] ``` The confidence intervals remain **pointwise**. BY-adjusted p-values do not turn them into simultaneous or selection-adjusted intervals. A reported discovery list also has no known realized false-discovery proportion in real data. ```{r} #| label: fig-case-contrasts #| fig-height: 5.5 #| fig-cap: "Twelve prespecified display features from the 24-feature family. Points and bars are estimates and pointwise approximate 95% intervals; crosses are known simulation truth. Blue indicates BY rejection. The intervals are not simultaneous or adjusted for selection." #| fig-alt: "A forest plot compares adjusted estimates and pointwise confidence intervals with known synthetic contrasts for six nonnull and six null features." show <- c(1:8, 12, 16, 20, 24); a <- results[show, ]; position <- rev(seq_len(nrow(a))) par(mar = c(4, 5, 1, 1)) plot(a$estimate, position, xlim = range(c(a$lower, a$upper, a$truth)), yaxt = "n", ylab = "", xlab = "Standardized contrast (log2 units)", pch = 16, col = ifelse(a$reject_BY, "#2F80C8", "gray35")) segments(a$lower, position, a$upper, position, col = "gray50") points(a$estimate, position, pch = 16, col = ifelse(a$reject_BY, "#2F80C8", "gray35")) points(a$truth, position, pch = 4, col = "#D55E00", cex = 1.1) axis(2, at = position, labels = a$gene, las = 1, cex.axis = 0.8) abline(v = 0, lty = 2, col = "gray60") ``` ## Challenge the analysis in repeated samples {#case-evaluation} Evaluate all features on the same synthetic datasets under four prespecified fitting choices: both model families correct; only outcome models correct; only weight models correct; both wrong. Incorrect models are intercept-only fits. Repeat at $n=400$ and $1200$. The two singly correct cases investigate double-robust consistency; the efficient-SE theorem above requires all nuisance functions consistent and does not promise interval calibration in those cases. ```{r} records <- EvaluateExpressionStudy(B = config$simulation_replicates, ns = config$simulation_n, genes = config$genes, folds = config$folds, seed = config$simulation_seed, alpha = config$alpha, probability_bounds = config$probability_bounds, mean_bounds = config$mean_bounds) performance <- SummarizeExpressionEvaluation(records) ``` ::: {.table-responsive} ```{r} performance[, c("n", "scenario", "attempted", "successful", "sentinel_bias", "bias_mcse")] ``` ::: The sentinel is prespecified null feature `gene_12`. Bias here is the average signed error across independent datasets; its MCSE uses the replicate-level SD divided by the square root of successful replicate count. A small bias in one finite experiment does not prove consistency. ```{r} #| label: fig-case-bias #| fig-cap: "Mean estimation error for the prespecified null gene_12. Bars are plus/minus two MCSEs for the simulated bias, not SEs of an individual gene estimate. Each scenario uses the same datasets and folds. Small horizontal offsets separate scenarios at each sample size. Only two sample sizes are shown; the asymptotic claims come from the remainder analysis." #| fig-alt: "Bias curves for four nuisance-model scenarios at 400 and 1200 people, with Monte Carlo error bars and a zero-bias reference line." scenarios <- unique(performance$scenario) colors <- c("#2F80C8", "#009E73", "#7B61A8", "#D55E00") layout(matrix(1:2, ncol = 1), heights = c(1, 4)) par(mar = c(0, 0, 0, 0)); plot.new() legend("center", scenarios, col = colors, pch = 15:18, lty = 1, ncol = 2, bty = "n", cex = 0.8) par(mar = c(4, 4, 0.5, 1)) plot(NA, xlim = range(config$simulation_n) * c(0.9, 1.08), ylim = range(c(0, performance$sentinel_bias + 2 * performance$bias_mcse, performance$sentinel_bias - 2 * performance$bias_mcse)), xlab = "Independent people n", ylab = "Mean error for gene_12") abline(h = 0, lty = 2, col = "gray50") for (j in seq_along(scenarios)) { a <- performance[performance$scenario == scenarios[j], ] x <- a$n * (1 + (j - 2.5) * 0.025) lines(x, a$sentinel_bias, type = "b", pch = j + 14, col = colors[j]) arrows(x, a$sentinel_bias - 2 * a$bias_mcse, x, a$sentinel_bias + 2 * a$bias_mcse, angle = 90, code = 3, length = 0.04, col = colors[j]) } ``` For interval and multiplicity summaries, first compare the fully specified setup with the deliberately wrong setup: ::: {.table-responsive} ```{r} performance[performance$scenario %in% c("Both correct", "Both wrong"), c("n", "scenario", "coverage", "coverage_mcse", "fdr", "fdr_mcse", "true_positive_fraction")] ``` ::: Within each replicate, coverage is the fraction of the 24 fixed targets covered, and FDP is false rejections divided by all rejections (zero when no tests reject). We average those replicate summaries. Their MCSEs use variation across independent datasets; **genes within a dataset are not independent simulation replicates**. A plug-in MCSE can be zero if every observed FDP is zero; that does not prove true FDR is zero. Each attempt is retained with seeds, success status, and any error. Metrics are conditional on successful complete-family fits, so report success counts alongside them; do not interpret conditional performance as unconditional reliability when fits fail. One hundred replicates make this a runnable teaching experiment, with visibly reported Monte Carlo uncertainty. A research benchmark should plan its replicate count from the precision needed. ## Reproduce the inputs, run, and report {#case-reproduction} Download the [complete project](../examples/adjusted-expression-project.zip), which includes the R functions, editable configuration, runner, observed synthetic assay/metadata, truth stored separately for evaluation, diagnostics, contributions, all simulation attempts, manifests, software information, Quarto report source, and rendered report. You can also [read the standalone report](../examples/adjusted-expression-report.html). The archive's `results/` is the supplied reference run. From the unpacked project root, rebuild into a new directory so the reference run remains available: ```{bash} #| eval: false Rscript run.R rerun-results quarto render report.qmd -P results_dir:rerun-results ``` The runner refuses to overwrite a populated output directory. It records input/output identity and code/configuration checksums for local change tracking. Analysis needs base R; rendering the report additionally needs a working Quarto R/knitr installation. The report reads explicit result files and has freezing disabled. It can be regenerated without an interactive workspace. A methods paragraph should state the scale, population, independent unit, target distribution, missingness assumptions, nuisance models, folds, clipping, family, uncertainty calculation, and inferential limits. For example, the following is generated from this run: ```{r} #| results: asis s <- results[1, ]; a <- fit$diagnostics cat(sprintf(paste0( "In this synthetic study, %d independently sampled people supplied covariates, ", "and %d supplied the whole expression assay. We standardized group contrasts ", "to the source-population covariate distribution using all people. ", "For gene_01, the cross-fitted augmented estimate was %.3f log2 units ", "(SE %.3f; pointwise approximate 95%% interval %.3f to %.3f). ", "The prespecified family contained %d features; BY at 0.05 rejected %d nulls. ", "Fitted probabilities and means required clipping for %.1f%% and %.1f%% ", "of their predictions. These are standardized associations under the stated ", "MAR, overlap, sampling, and nuisance-rate conditions; they do not establish ", "causal or absolute-abundance effects.\n\n"), a$n, a$observed, s$estimate, s$se, s$lower, s$upper, nrow(results), sum(results$reject_BY), 100 * a$probability_clipped, 100 * a$mean_clipped)) ``` The scientific interpretation is bounded by the target and assumptions, even when execution is reproducible and the synthetic answer is close to truth. Return to [Reporting a statistical methods study](../research-practice/reporting-methods-studies.qmd) for the full claim/evidence checklist. ## Exercises Try each checkpoint before opening its hint or solution. ::: {#checkpoint-case-01 .checkpoint} ### Checkpoint · CASE-01 — Specify and identify the expression contrast {.unnumbered} **Task type:** Interpretation and derivation. Write the independent observation, measurement scale, population distribution, estimand, and tested family for this case. Derive why the generator gives $\psi_g=\delta_g$. Explain how a feature can have $\psi_g=0$ with conditional contrasts that vary with $Z$, and state the missingness/support assumptions needed to identify the full-population contrast. **Ready to continue when:** you can distinguish the standardized target from a complete-case crude contrast, a constant regression coefficient, and a causal effect. **Transfer:** If the target covariate distribution were that of group 1, would the same delta remain the target? Explain what changes in estimation and uncertainty. <details class="checkpoint-help"> <summary>Hint</summary> Average the conditional contrast over the specified population, then identify which outcome laws are actually observed. </details> <details class="checkpoint-help"> <summary>Solution outline — open after trying</summary> One independent observation is a person’s covariates, group, observation indicator, and observed assay vector. The outcome is calibrated log2 intensity, and the family contains all 24 feature-wise standardized-contrast nulls. The conditional contrast is $\delta_g+h_gZ$; since the population has $E[Z]=0$, its average is $\delta_g$. A zero average need not imply zero contrast at every covariate value. MAR within group/covariates identifies the outcome regressions from observed outcomes; positive group and observation probabilities are needed throughout the target support, and covariates must be observed for the target population. None of these statements alone gives a causal or absolute-abundance interpretation. </details> ::: ::: {#checkpoint-case-02 .checkpoint} ### Checkpoint · CASE-02 — Reconstruct the estimator and its remainder {.unnumbered} **Task type:** Derivation and implementation. Derive the conditional expectation of each residual correction, then the two-term exact bias displayed in this page. Explain why the weight uses the joint probability $q_a$. Reconstruct the fold procedure without accessing simulation truth or unobserved outcomes, and identify the conditions needed before using the estimated influence-function SE. **Ready to continue when:** you can connect every correction and condition to an observed-data quantity and explain what double-robust consistency does and does not establish. **Transfer:** If observation probabilities depend on the unobserved expression even after conditioning on group and $C$, which step of the argument fails? <details class="checkpoint-help"> <summary>Hint</summary> Hold a training fit fixed. For group a, the expected observed residual is q_a times the difference between the true and fitted conditional mean. </details> <details class="checkpoint-help"> <summary>Solution outline — open after trying</summary> The group-a correction has fresh-observation expectation $E[q_a(\mu_{ag}-\hat\mu_{ag})/\hat q_a]$. Combining it with each plug-in mean yields the displayed signed product remainders. Weighting only by group or observation probability fails to cancel the joint observation process. Fit all nuisance functions on other people, evaluate both group predictions on all held-out covariates, and correct residuals only on observed rows of the matching group. For efficient inference, check nuisance consistency, lower probability bounds, product little-o rates, the held-out empirical term, finite variance and sampling/moment conditions. A wrong but controlled nuisance can permit consistency while invalidating that efficient-SE argument. </details> ::: ::: {#checkpoint-case-03 .checkpoint} ### Checkpoint · CASE-03 — Interpret discoveries and uncertainty correctly {.unnumbered} **Task type:** Interpretation and computation. Explain the difference between a pointwise 95% interval, BY FDR control for the 24-feature family, a realized FDP in one simulation, and MCSE for estimated FDR. State the assumptions left after selecting BY and explain why you cannot treat 24 genes in each of 100 datasets as 2,400 independent replicates. **Ready to continue when:** you can label each uncertainty quantity and error criterion, including the conditions multiplicity adjustment cannot repair. **Transfer:** If you display only the twelve prespecified plot features, does the tested family shrink to twelve? What if those twelve were chosen after inspecting significance? <details class="checkpoint-help"> <summary>Hint</summary> Ask whether the expectation is over people, repeated datasets, or the feature family. </details> <details class="checkpoint-help"> <summary>Solution outline — open after trying</summary> A pointwise interval targets one feature’s population contrast across samples. BY targets the expected FDP across repeated analyses of a stated family, allowing arbitrary p-value dependence but requiring valid null p-values; the Wald p-values here are asymptotically justified under the estimator’s conditions. One simulated FDP is known only because truth is supplied. Its average estimates FDR; its MCSE is the SD of replicate FDPs divided by the square root of the number of independent successful datasets. Within-dataset genes share people, noise, and fits. BY does not create simultaneous confidence intervals or repair missingness, target, or sampling errors. </details> ::: ::: {#checkpoint-case-04 .checkpoint} ### Checkpoint · CASE-04 — Reproduce and challenge the complete analysis {.unnumbered} **Task type:** Execution and diagnosis. Unpack the project, regenerate it in a new results directory, render its report, and compare target estimates and simulation records with the reference run. Audit the input alignment, missing-outcome conventions, seeds, code/configuration identities, warnings, clipping, and attempted/successful counts. Before examining the four scenario summaries, predict which cases support consistency and which support the efficient-SE theorem. **Ready to continue when:** you can rebuild the result, trace a claim to its inputs and assumptions, and diagnose a failed condition without silently discarding failures. **Transfer:** Change only the number of folds or one model class, record a new run, and distinguish stochastic changes from a changed target or invalid theorem. <details class="checkpoint-help"> <summary>Hint</summary> Use a fresh session and retain the reference directory. Derive the exact bias prediction before looking at the simulation plot. </details> <details class="checkpoint-help"> <summary>Solution outline — open after trying</summary> The documented runner sources explicit functions and configuration, writes observed inputs separately from oracle truth, and records results, diagnostics, attempts, manifests, and software information. With the tested environment and seeds the statistical outputs reproduce; session metadata can differ across environments. Correct outcomes or correct joint weights can remove population bias under the stated controls. The fully consistent, sufficiently fast nuisance case supports the efficient-SE argument; single-model correctness alone does not. Both wrong leaves a generally nonzero remainder. Inspect failures and their denominators rather than assuming missing rows succeeded. </details> ::: ### Additional practice (optional) {.unnumbered} 1. Extend the generator to feature-specific missingness and derive the corresponding observation probabilities; do not silently reuse the whole-assay indicator. 2. Use an external target covariate sample. Derive the new target and account for its sampling variation instead of treating the current influence function as unchanged. 3. Compare the four scenarios across more sample sizes and choose replicate counts from a prespecified MCSE goal. Treat singly correct nuisance cases as a separate inference problem. 4. Construct a simultaneous interval procedure for the prespecified family and explain how its claim differs from the plotted pointwise intervals. ## Reading - Kennedy, [review of influence-function estimation](https://arxiv.org/html/2203.06469v2): standardized means, estimator construction, and remainder analysis. - Benjamini and Yekutieli (2001), [FDR under dependency](https://doi.org/10.1214/aos/1013699998), and the official [R p.adjust documentation](https://stat.ethz.ch/R-manual/R-devel/library/stats/html/p.adjust.html). ## Continue through the guide {#continue-through-the-collection} Return to the [Omics analysis problems](index.qmd) overview, compare with the [Differential expression and abundance](differential-expression-and-abundance.qmd) page, or assess the four case-study tasks in the [Completion checkpoints](../estimation-and-inference/checkpoints.qmd) page.