Jackknife pseudo-values

A number for every observation: how it is built, what it measures, and when it can be used as data.

Before reading: Sample correlation and the What is your estimator actually estimating? page; Influence functions explains the first-order expansion.

1 What is the population quantity?

A jackknife pseudo-value is constructed from a cohort-level estimator; its observation-level appearance does not by itself make it a measurement of a person’s biological network. Before using pseudo-values as outcomes, define the population quantity their group means or regression coefficients should estimate. Because the values share the full cohort and leave-one-out fits, they are not independent raw observations. Inference needs a justification for that construction and the sampling design.

The examples below use correlation to expose the effect of pooled group composition. The construction is a useful optional application of resampling and Influence functions; the main semiparametric learning path can proceed without it.

2 The construction

Let \hat\theta be any statistic computed from N observations, typically an estimate of some feature \theta(F) of the population distribution F, and let \hat\theta_{(-i)} be the same statistic computed with observation i removed. The jackknife pseudo-value for observation i is

\theta_i = N\,\hat\theta - (N - 1)\,\hat\theta_{(-i)}.

Rearranged, \theta_i = \hat\theta + (N-1)(\hat\theta - \hat\theta_{(-i)}): the overall estimate, plus N - 1 times the amount observation i moved it. An observation that pulls the estimate up gets a high pseudo-value; one that pulls it down gets a low one; one with no influence gets \hat\theta itself.

“The same statistic” means rerunning its required estimation steps on the delete-one data. Holding part of a fitted pipeline fixed defines a different deletion calculation and needs its own argument. Here an observation is an independent person with their pair of gene values; deleting one gene is not the same operation. For clustered sampling, a person-level deletion can remove several records at once.

The picture below (Figure 1) shows three observations from one sample of 20.

Figure 1: How a pseudo-value is built, for one sample of 20 observations on two genes whose true correlation is 0.5. Each panel highlights one observation (vermilion). The correlation is computed with all 20 observations and again without the highlighted one, and the pseudo-value is the first number plus 19 times the difference. An observation far out along the diagonal pulls the correlation up and gets a high pseudo-value. An observation near the center barely moves the correlation and gets about the overall value. An observation far off the diagonal pulls the correlation down and gets a low one. Pseudo-values are not correlations, and they can leave the interval from -1 to 1.

PseudoValues <- function(x, statistic, ...) {
  n <- length(x)
  theta_all <- statistic(x, ...)
  sapply(seq_len(n), function(i) n * theta_all - (n - 1) * statistic(x[-i], ...))
}

Quenouille (1949, 1956) introduced leave-one-out recomputation to reduce bias; Tukey (1958) named the jackknife and proposed the pseudo-values, suggesting they could be treated approximately as independent observations. That suggestion is why pseudo-values are used as data today, and why they can mislead.

3 Three facts you can check in a few lines of code

If \hat\theta = \bar x, then \theta_i = x_i exactly, so for the mean the pseudo-values are the data. Try it.

Indeed, (N-1)\bar x_{(-i)}=N\bar x-x_i. Substituting that identity into the definition cancels the entire cohort sum and leaves x_i. The cancellation works because a mean is a sum of fixed observation-level contributions. A correlation has estimated means and SDs in its numerator and denominator; deleting one person changes those quantities for every retained person.

set.seed(1)
x <- rnorm(8)
all.equal(PseudoValues(x, mean), x)

[1] TRUE

The average of the pseudo-values is the bias-corrected jackknife estimate, \bar\theta = N\hat\theta - (N-1)\overline{\hat\theta_{(\cdot)}}. For the mean it equals \bar x. For sufficiently smooth statistics such as the variance example below, its difference from \hat\theta is often of order 1/N under suitable moments. Nonlinearity by itself does not establish this rate. A classic example is the plug-in variance, which divides by N and is biased low by the factor (N-1)/N. Its jackknife estimate is exactly the sample variance with divisor N - 1.

PlugInVar <- function(x) mean((x - mean(x))^2)     # divisor N
c(plug_in = PlugInVar(x), mean_pseudovalue = mean(PseudoValues(x, PlugInVar)), var_N_minus_1 = var(x))

         plug_in mean_pseudovalue    var_N_minus_1 
       0.6369472        0.7279396        0.7279396 

For this divisor-N plug-in variance, the pseudo-value has the exact form

\theta_i=\frac{N}{N-1}(x_i-\bar x)^2.

Its average is the divisor-(N-1) sample variance. It still uses the pooled mean \bar x, so a person’s large value can reflect displacement of their group from that pooled center, rather than large within-group variation. This simple example foreshadows the correlation artifact below. If the starting statistic is instead R’s var() with divisor N-1, the exact pseudo-value formula changes; optional exercise 3 asks for it.

The sample variance of the pseudo-values gives the familiar jackknife variance formula. Let \bar\theta be their mean and \bar\theta_{(-\cdot)} the mean of the delete-one estimates. Centering the definition gives

\begin{aligned} \theta_i-\bar\theta &=-(N-1)(\hat\theta_{(-i)}-\bar\theta_{(-\cdot)}),\\ \frac{s_\theta^2}{N} &=\frac{N-1}{N}\sum_i (\hat\theta_{(-i)}-\bar\theta_{(-\cdot)})^2. \end{aligned}

Thus computing the usual s_\theta/\sqrt N from the pseudo-values reproduces the jackknife SE exactly. This is an algebraic identity between two formulas. Whether either consistently estimates sampling uncertainty still depends on the estimator’s smoothness, moments, and sampling design; the identity does not make the pseudo-values independent or justify an arbitrary regression SE.

pv <- PseudoValues(x, var)
loo <- sapply(seq_along(x), function(i) var(x[-i]))     # leave-one-out estimates
n <- length(x)
c(jackknife_se = sqrt((n - 1) / n * sum((loo - mean(loo))^2)),
  se_of_pseudovalue_mean = sd(pv) / sqrt(n))

          jackknife_se se_of_pseudovalue_mean 
              0.334519               0.334519 

4 Pseudo-values as data

Two literatures took Tukey’s suggestion further and use pseudo-values as a response variable.

In survival analysis, a subject’s survival status at time t is unobserved if they were censored before t. Andersen, Klein, and Rosthøj (2003) compute the Kaplan-Meier estimate of S(t) with and without each subject, form the pseudo-value, and regress the pseudo-values on covariates with generalized estimating equations. The pseudo-value stands in for the unobservable indicator 1\{T_i > t\}. The same idea extends to cumulative incidence and restricted mean survival time. The theory (Graw, Gerds, and Schumacher 2009; Overgaard, Parner, and Pedersen 2017) establishes when this works, meaning when a regression of the pseudo-values on covariates estimates how the probability of surviving past t depends on the covariates. The estimand is an average of a per-subject quantity, and, under independent-subject sampling and suitable censoring/support conditions, the pseudo-value can have, to first order,¹ the same conditional mean given the covariates as that per-subject quantity. A regression on the pseudo-values then estimates the regression of the indicator 1\{T_i > t\} on the covariates, which is what the analysis is after.

In sample-specific network estimation, LIONESS (Kuijjer et al. 2019) estimates a co-expression network for every sample, starting from a method that estimates one network from all N samples. With Pearson correlation as that method, \hat\theta is the correlation between two genes across all N samples, \hat\theta_{(-i)} is the same with sample i removed, and LIONESS defines sample i’s edge as

e_i = N\,\hat\theta - (N - 1)\,\hat\theta_{(-i)},

which is the pseudo-value of the correlation. The resulting per-sample edges are then compared between groups, for example with t-tests or limma, regressed on phenotypes, or clustered. Other single-sample network methods use related ideas. SSN (Liu et al. 2016) measures how the correlation changes when a sample is added to a reference set, and BONOBO (Saha et al. 2024) builds each sample’s coexpression from a Bayesian model whose prior is formed from all the other samples.

# LIONESS edge for one gene pair = pseudo-value of the correlation
set.seed(2)
n <- 60
g1 <- rnorm(n); g2 <- 0.5 * g1 + sqrt(0.75) * rnorm(n)
r_all <- cor(g1, g2)
edge <- sapply(seq_len(n), function(i) n * r_all - (n - 1) * cor(g1[-i], g2[-i]))
c(mean_edge = mean(edge), r_all = r_all)

mean_edge     r_all 
0.5168562 0.5153655 

The mean edge, 0.517, is close to the pooled correlation, 0.515, as the second fact says.

5 What a pseudo-value measures

The trap is in how the number is read.

In both literatures the pseudo-value \theta_i is read as a noisy measurement of the quantity of interest for people like subject i (their survival probability, or their correlation). For survival, that reading holds to first order under stated conditions,² because S(t) is an average of individual indicators, and the average pseudo-value among people like subject i estimates their survival probability. For a correlation, the reading fails in a way that can be written down: the average pseudo-value among people like subject i need not estimate their correlation.

A pseudo-value measures the influence of observation i on the statistic computed from the whole pooled sample. It is anchored at \hat\theta, the pooled estimate. If the sample is a mixture of two groups with different means, an observation from the smaller group sits far from the pooled center, has a large influence on the pooled correlation, and gets a pseudo-value that reflects its group membership. That pseudo-value says little about the correlation within its group.

The calculation below connects a pseudo-value to the pooled estimate plus the observation’s influence function, to first order. The practical version is this simulation: two groups with identical within-group correlation, different means, and a mean difference in their LIONESS edges that is entirely an artifact.

set.seed(3)
n_ctrl <- 200; n_pat <- 800; rho <- 0.5; shift <- 2
group <- rep(c("control", "patient"), c(n_ctrl, n_pat))
n <- n_ctrl + n_pat
z1 <- rnorm(n); z2 <- rho * z1 + sqrt(1 - rho^2) * rnorm(n)
a <- z1 + shift * (group == "patient"); b <- z2 + shift * (group == "patient")
LionessEdge <- function(x, y) {
  n <- length(x); r <- cor(x, y)
  sapply(seq_len(n), function(i) n * r - (n - 1) * cor(x[-i], y[-i]))
}
e <- LionessEdge(a, b)
tapply(e, group, mean)                               # groups differ...

 control  patient 
1.078159 0.595102 

c(control = cor(a[group == "control"], b[group == "control"]),
  patient = cor(a[group == "patient"], b[group == "patient"]))   # ...but the correlations do not

  control   patient 
0.5094897 0.4758453 

In this sample the smaller group, the controls, has a mean edge of 1.08 and the patients 0.60, while both groups’ correlations are about 0.5 (0.51 and 0.48, which differ only by sampling variation). In the vocabulary of the What is your estimator actually estimating? page, the estimand for a group’s edge is that group’s own correlation, and the group’s average edge converges to a different limit.

The direction is systematic: the smaller group gets the higher edges. The gap changes sign when the two sample sizes are swapped and vanishes at balance, although at balance each group’s mean edge is still far above the true correlation (about 0.75 against 0.5 in a large sample with a shift of 2). The picture below (Figure 2) shows the mechanism and the size of the effect.

Figure 2: Why the smaller group gets the higher edges. Left: the sample of 20 controls (open circles) and 80 patients (filled circles) from Association targets, pooling, and imperfect adjustment, with the patients shifted by 2 in both genes, so that the correlation within each group is 0.5. The pooled center, which averages over everyone, lies close to the patients and far from the controls, so a control has a larger influence on the pooled correlation. Right: what each group’s average edge converges to as the fraction of controls changes, with the shift fixed at 2. The dotted line marks the 20% controls of the simulation above. The estimand, the within-group correlation, is 0.5 throughout (blue). The limits (vermilion) differ between the groups, are equal at balance where both are 0.75, and stay above 0.5. Curves use the first-order formula developed below, and dots are averages over simulated samples of 20,000.

In this equal-variance mixture, the asymptotic gap vanishes if each gene is centered within its group before the edges are computed: the relevant pooled and group moments then agree. A variance difference can create an additional artifact that centering does not remove. Group standardization can align the required moments in this common-correlation example. These are repairs for the stated target mismatch, not a general proof of pseudo-value independence or regression coverage. Optional exercise 2 explores the mean and variance changes.

6 When the average pseudo-value estimates what you want

Earlier we described how pseudo-values are used: regress them on covariates, or compare their averages between groups. Behind that use is a claim that the average pseudo-value among people with covariates x estimates the quantity of interest for people like them. In the vocabulary of the What is your estimator actually estimating? page, that quantity is the estimand \theta(F_x), the same feature \theta (a correlation, a survival probability) evaluated on F_x, the distribution of the subpopulation with covariates x. The average pseudo-value is computed from data, so it is an estimator, and it converges to a limit. Pseudo-value regression works as intended when the limit equals the estimand, because then the coefficients describe how the quantity of interest changes with the covariates. This is a claim about what the regression estimates. Standard errors and p-values are a separate question.

The limit equals the estimand when the estimand is linear in the distribution, meaning that \theta(F) is an average \int h\,dF of something computed on each observation. Then the pseudo-value is exactly h(x_i), as the mean’s pseudo-values are the data, and the average pseudo-value among people with covariates x is the average of h in their subpopulation, which is \theta(F_x). Survival probabilities, cumulative incidences, and restricted means are of this kind for complete data. Their observed-data pseudo-value regressions additionally need the censoring and estimation conditions above. (With censoring, the Kaplan-Meier pseudo-value approximates that quantity to first order.)

Correlations, variances, and regression coefficients are nonlinear in F, so this exact linear-functional argument is unavailable. Their pseudo-values can carry the pooled sample’s composition into subgroup averages (see Figure 2). Whether a discrepancy occurs depends on the subgroup distributions and functional; nonlinearity alone does not prove a nonzero bias in every configuration. The relevant limit and remainder must be checked, as developed in Association targets, pooling, and imperfect adjustment.

7 Check the influence approximation (optional)

This numerical check uses a fresh bivariate-normal sample, separately from the earlier cohort example.

set.seed(1)
n <- 500
x <- rnorm(n); y <- 0.6 * x + 0.8 * rnorm(n)

# Weighted correlation of (x, y) with weights w
WCor <- function(x, y, w) {
  cv <- cov.wt(cbind(x, y), wt = w, method = "ML")$cov
  cv[1, 2] / sqrt(cv[1, 1] * cv[2, 2])
}

# Numerical IF of the correlation at the empirical distribution, evaluated at a new point
NumericIF <- function(x0, y0, eps = 1e-4) {
  w <- c(rep((1 - eps) / n, n), eps)
  (WCor(c(x, x0), c(y, y0), w) - cor(x, y)) / eps
}

# Analytic IF, standardized with divisor N (the scaling of the functional at the empirical distribution)
SdN <- function(v) sqrt(mean((v - mean(v))^2))
AnalyticIF <- function(x0, y0) {
  r <- cor(x, y)
  xt <- (x0 - mean(x)) / SdN(x); yt <- (y0 - mean(y)) / SdN(y)
  xt * yt - r / 2 * (xt^2 + yt^2)
}

pts <- rbind(c(0, 0), c(2, 2), c(2, -2), c(-1, 3))
cbind(pts, numeric = apply(pts, 1, function(p) NumericIF(p[1], p[2])),
      analytic = apply(pts, 1, function(p) AnalyticIF(p[1], p[2])))

                 numeric      analytic
[1,]  0  0 -0.0007942872 -0.0007942872
[2,]  2  2  1.7040853368  1.7047466919
[3,]  2 -2 -5.9163710249 -5.9186147672
[4,] -1  3 -5.7475178557 -5.7507604588

if_vals <- AnalyticIF(x, y)
c(if_based_se = sd(if_vals) / sqrt(n),
  fisher_se   = sqrt((1 - cor(x, y)^2)^2 / n))   # textbook SE of r under normality

if_based_se   fisher_se 
 0.02941070  0.03066164 

The two standard errors agree to within sampling noise (0.029 and 0.031 here). The influence-function version does not assume normal data, and the textbook formula does.

The second use is the link to pseudo-values. Removing observation i from F_N gives F_{N,-i} = F_N + \frac{1}{N-1}(F_N - \delta_{X_i}). Write \phi_N(x)=\operatorname{IF}(x;T,F_N) for the empirical influence function. Apply the expansion with G=F_{N,-i} and F=F_N, using \int\phi_N\,dF_N=0:

T(F_{N,-i})-T(F_N)\approx-\frac{\phi_N(X_i)}{N-1}.

Substitute into the pseudo-value \theta_i = T(F_N) + (N-1)\big(T(F_N) - T(F_{N,-i})\big):

\theta_i\approx T(F_N)+\phi_N(X_i).

For a sufficiently smooth plug-in functional, this deletion expansion can give the estimate plus its empirical influence value, with an O_p(1/N) error for a typical observation under suitable moment and local-curvature control.³ A generic root-N asymptotic expansion alone does not imply this deletion rate: the leave-one-out remainder is multiplied by N-1. This is the precise content of Tukey’s heuristic, and it is exact when T is linear (for the mean, \bar x + (x_i - \bar x) = x_i).

r <- cor(x, y)
pv <- sapply(seq_len(n), function(i) n * r - (n - 1) * cor(x[-i], y[-i]))
summary(pv - (r + AnalyticIF(x, y)))   # differences are O(1/N)

      Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
-1.690e-02 -4.204e-04 -1.185e-05  3.534e-04  2.229e-04  5.668e-02 

8 Exercises

Checkpoints are completion targets for this page. Additional practice below is optional. Each checkpoint gives a completion standard, a hint, and a solution outline.

Checkpoint · PV-01 — Verify the linear pseudo-value case

Task type: Derivation and interpretation.

Let \hat\theta=N^{-1}\sum_i h(X_i), where h is a fixed function and the required expectation exists. Derive the delete-one estimate and show that N\hat\theta-(N-1)\hat\theta_{(-i)}=h(X_i). Explain which step need not hold for a nonlinear functional such as correlation.

Ready to continue when: you can derive the identity and explain why it does not make every pseudo-value a raw biological measurement.

Transfer: For h(x)=x^2, what does the average pseudo-value target, and how does that differ from the variance?

Hint

Write (N-1)\hat\theta_{(-i)} as the full sum with term i removed.

Solution outline — open after trying

Multiplying the full estimate by N gives \sum_jh(X_j). Multiplying the delete-one estimate by N-1 gives \sum_jh(X_j)-h(X_i). Subtraction leaves h(X_i). For a nonlinear statistic, multiplication by sample size does not turn the estimator into a sum of fixed observation contributions; the cohort and leave-one-out fits remain part of each pseudo-value.

Additional practice (optional)

These exercises extend the ideas through further derivations, simulations, or applications. Choose those useful for your work.

1. Reproduce the two-group simulation above with 500 controls and 500 patients. What happens to the edge difference, and how do the two mean edges compare with the true correlation of 0.5? Then try 800 controls and 200 patients. What is the sign of the difference now?
2. Center each gene within its group (a - ave(a, group)) and recompute the edges. Then give the patients a larger within-group variance in gene a only (multiply their centered a values by 2) and check whether centering alone still fixes the problem. Use large groups (thousands) so that the bias is large compared with the noise.
3. Compute pseudo-values for the sample variance of a normal sample and plot them against (x_i - \bar x)^2. They are nearly, but not exactly, equal. Why “nearly”? In fact the relationship is an exact straight line, so find its slope and intercept as functions of N and s^2.

9 Reading

• Quenouille, M. H. (1956). Notes on bias in estimation. Biometrika, 43(3-4), 353 to 360.
• Tukey, J. W. (1958). Bias and confidence in not-quite large samples (abstract). Annals of Mathematical Statistics, 29, 614.
• Efron, B. (1982). The Jackknife, the Bootstrap and Other Resampling Plans. SIAM. Chapters 1 to 3.
• Andersen, P. K., Klein, J. P., and Rosthøj, S. (2003). Generalised linear models for correlated pseudo-observations, with applications to multi-state models. Biometrika, 90(1), 15 to 27.
• Andersen, P. K. and Pohar Perme, M. (2010). Pseudo-observations in survival analysis. Statistical Methods in Medical Research, 19(1), 71 to 99. A practical review of the survival side.
• Graw, F., Gerds, T. A., and Schumacher, M. (2009). On pseudo-values for regression analysis in competing risks models. Lifetime Data Analysis, 15(2), 241 to 255.
• Overgaard, M., Parner, E. T., and Pedersen, J. (2017). Asymptotic theory of generalized estimating equations based on jack-knife pseudo-observations. The Annals of Statistics, 45(5), 1988 to 2015. The theory behind pseudo-value regression, including the conditions under which it estimates what it should.
• Kuijjer, M. L., Tung, M. G., Yuan, G., Quackenbush, J., and Glass, K. (2019). Estimating sample-specific regulatory networks. iScience, 14, 226 to 240.
• Liu, X., Wang, Y., Ji, H., Aihara, K., and Chen, L. (2016). Personalized characterization of diseases using sample-specific networks. Nucleic Acids Research, 44(22), e164.
• Saha, E., Fanfani, V., Mandros, P., Ben Guebila, M., Fischer, J., Shutta, K. H., DeMeo, D. L., Lopes-Ramos, C. M., and Quackenbush, J. (2024). Bayesian inference of sample-specific coexpression networks. Genome Research, 34, 1397 to 1410.
• Kuijjer, M., De Marzio, M., and Glass, K. (2026). Challenges and opportunities in single-sample network modeling. bioRxiv preprint, doi:10.64898/2026.02.27.708608. Recasts several single-sample network methods with common notation and compares what each one assumes.

Footnotes

1. “To first order” is a phrase from calculus. It means keeping the largest term of an approximation and dropping the smaller ones. The standard picture is a curve near one point: the tangent line follows the curve closely there, and the gap between them grows like the square of the distance from the point. For example, \sqrt{1 + x} \approx 1 + x/2 for small x, and the dropped terms are of order x^2. A second-order approximation keeps the next term as well.

   Here the point is the pooled sample and the line is the influence function, a derivative defined in the Influence functions page. A pseudo-value is treated as the estimate computed from everyone plus one correction term for observation i, with everything smaller ignored. The ignored part is exactly zero when the quantity being estimated is linear in the distribution, as the mean is, and it can be large when the quantity is nonlinear and the group of interest differs from the pooled sample, as with the correlation in the next section. A claim “to first order” is only as good as the reason to think the ignored terms are small.

2. Results in statistics are proved by assuming conditions, such as independent observations, a correct model, or a large sample. “Under stated conditions” means the result is guaranteed when the conditions the authors list hold, and it can fail when they do not. It is a prompt to find the list and check each item against your data.

   For pseudo-value regression with survival data, the conditions in the theory papers cited above include the following. The quantity being modeled is an average of a per-subject quantity, such as the probability of surviving past a fixed time, a cumulative incidence, or a restricted mean survival time. A simple sufficient censoring setup has independent subjects, censoring time independent of the joint event-time/covariate vector (including event type for competing risks), and positive censoring survival over the time range studied. Independence from covariates alone is not enough if censoring still depends on event time. The needed smoothness, moment, and large-sample conditions must also hold. And, as in any regression, the model for how the covariates relate to the average is adequate. When censoring depends on the covariates, pseudo-values built from the ordinary Kaplan-Meier estimate can be biased.

3. Here O_p(1/N) means multiplying the random error by N produces a sequence bounded in probability. It is a rate statement, not a promise that every realized error halves when N doubles. Under the stated local controls it is smaller than the root-N sampling scale, since N^{-1}/N^{-1/2}\to0. A bound for a typical observation also need not be a uniform bound over the most extreme observations in a growing sample.

   Here the remainder is the leftover from the expansion above, evaluated at one observation. For a typical observation it is tiny, but it is larger for an observation with a very large influence function, because the leftover grows with the square of the influence. The code below shows this: most differences are within 0.001 in absolute value, and the largest, at an extreme point, is about 0.06 when N = 500.