Association targets, pooling, and imperfect adjustment

Estimands, parameters, misspecification, and pulling out the bias term.

Before reading: The target/estimator distinction in What is your estimator actually estimating?, covariance, and basic regression. This is a worked network-analysis extension; it is not an entry requirement for the theory guides.

1 Specify the research question

A useful target statement names the population, the independent sampling unit, the measurement scale, the contrast, and any adjustment. For example, let A indicate group and let Y_g be a gene’s measured log-expression. The descriptive contrast E[Y_g\mid A=1]-E[Y_g\mid A=0] is not automatically a treatment effect. Averaging conditional contrasts over a specified distribution of covariates defines a further, distinct target.

For an edge, \operatorname{Corr}(X,Y) in the pooled population, \operatorname{Corr}(X,Y\mid A=a) within a group, and a residual association after adjustment for Z are different estimands. A pooled correlation can consistently estimate the first while failing to estimate the second. That is a mismatch between question and procedure, even when the procedure has no difficulty estimating its own target.

Defining a target is also different from identifying it from observed data. Sampling, missingness, measurement, and any causal interpretation require their own assumptions. Use the Research questions, study design, and estimands overview before fitting the final model.

2 The one idea

2.1 Working model, population, and misspecification

A model is a set of distributions, usually described by an equation or a list of assumptions with a few free numbers. A working model supplies assumptions used to justify a particular interpretation or uncertainty calculation. A procedure may remain well-defined outside that model. Sample correlation estimates the pooled population correlation under suitable moment conditions, including in a heterogeneous mixture. Least squares estimates a linear projection without requiring a correct linear conditional mean. Interpreting either as a within-group association, a structural coefficient, or a causal effect requires additional assumptions. State which claim the working model supports.

The data come from the population, a single distribution we write F (for two variables, their joint distribution). Nobody knows F in real life. In theory and in simulation we choose one. We write down a richer model, a bigger set of distributions that contains the working model as a special case (for example, two groups with different means, or an outcome that also depends on a variable we left out), and we pick one member of it to play the role of F. So the working model and the richer model are both things we write down, and F is one distribution. They differ in their role. The working model is what the procedure assumes, and F is what we ask the procedure to face.

We call it model misspecification when F lies outside the working model, so that the assumptions behind the procedure are false. Research mostly lives here: a procedure is designed for one setting and used in another. The structure the working model leaves out may be observed but ignored (group labels we have but the procedure never uses) or truly latent (hidden subtypes, unmeasured cell composition). The bias formulas below are the same either way. What changes is what you can do about the bias: known labels let you adjust for it directly, while latent structure can usually only be diagnosed or estimated. The question of this page is what a procedure is actually estimating when the population departs from its working model in one specific way.

2.2 Three objects that are easy to confuse

Answering that question means keeping three objects apart:

• The estimand is the quantity you want. It is a feature of a population, not of your data: a mean, a correlation, a regression slope, a network edge. Write it as \theta(F), a functional of a distribution F.
• The estimator is the recipe you apply to data: \hat\theta = T(\text{data}).
• The limit of the estimator is what the recipe converges to as the sample grows (its probability limit). Call it \theta^\ast. When a model is misspecified, you will also see it called the pseudo-true value.

The picture below (Figure 1) puts the models and the three objects together. In this example the richer model is two groups whose means differ by \delta, and \delta = 0 is the working model. The population F is the member of the richer model at \delta = 2, and turning \delta up asks what the procedure does as the departure from its working model grows.

Figure 1: The estimand, the estimator, and the limit, and where the models come in. Left: the richer model is a family of worlds F_\delta indexed by a departure \delta (here, a shift in group means). The working model, the one the procedure assumes, sits inside it as the special case \delta = 0 (drawn as a point, with the other parameters held fixed). The population F = F_\delta at \delta = 2 is the world we ask about. The question asked of F gives the estimand \theta(F). Drawing data from F and applying the recipe T gives an estimate \hat\theta_n, which converges to the limit \theta^\ast as n grows. Right: as the departure is turned up from zero, the estimand stays put and the limit moves, and the gap between them is the bias term. The curves are for the pooled-correlation example below (within-group correlation 0.5, 20% controls).

A fourth word is probably more familiar than any of these: the parameter. It is easy to blur with the estimand and the limit, so it helps to be exact. A parameter belongs to a model: it is a number that indexes a family of distributions, such as \beta in Y = \beta X + \varepsilon, and its interpretation depends on the model and domain on which it is defined. A model parameter may also define a useful functional beyond that model, such as a linear projection coefficient. An estimand belongs to the question and the population, and must be defined on a stated set of distributions. Some scientific targets require identifying assumptions before they can be expressed using observed data. A limit belongs to an estimator and the population, and is whatever that estimator converges to, again whether or not the model behind it is right. The first example below pins this down, and a table after it summarizes.

Statistics courses spend most of their time on the case where the working model is correct, so the parameter, the estimand, and the limit are all the same number and you only have to worry about variance. Under misspecification the limit quietly moves. The question “what is this estimator actually estimating?” means: compute \theta^\ast under the conditions you actually have, and compare it to \theta(F).

2.3 The bias term

The difference \theta^\ast - \theta(F) is the bias term (strictly, the asymptotic bias, since \theta^\ast is a limit). Writing it down explicitly is the core of a great deal of methodological work:

\text{procedure's limit}=\text{target}+\text{asymptotic bias}.

Once the bias term is a formula, two things become possible. If it depends on quantities you can estimate, you can subtract it or redesign the procedure so it vanishes: a correction. If it depends on quantities you can test for, you can check whether the conditions that create it are present: a diagnostic. A bias you have not written down can do neither.

3 A first example: omitted variable bias

The classic case. We posit a richer model for the population, in which the outcome depends on two variables:

Y = \beta X + \gamma Z + \varepsilon,

with mean-zero \varepsilon independent of (X,Z), finite second moments, and \operatorname{Var}(X)>0. Here \beta describes the conditional mean slope at fixed Z. Reading it as an intervention effect would require a causal model and additional identifying assumptions. An equation like this, together with distributions for X and Z, pins down a joint distribution F of (X,Y,Z). The analyst leaves out Z and fits

Y = a + bX + e,

but wants its slope b to describe the conditional slope \beta. Setting \gamma=0 is one sufficient case in which those coefficients agree. More generally, the least-squares slope of Y on X alone converges to its linear projection coefficient:

\begin{aligned} \beta^\ast&= \frac{\operatorname{Cov}(X, Y)}{\operatorname{Var}(X)}\\ &= \beta + \gamma\,\frac{\operatorname{Cov}(X, Z)}{\operatorname{Var}(X)}. \end{aligned}

To derive this, expand \operatorname{Cov}(X,Y) using the richer equation. The three terms are \beta\operatorname{Var}(X), \gamma\operatorname{Cov}(X,Z), and \operatorname{Cov}(X,\varepsilon)=0. Dividing by \operatorname{Var}(X) gives the result. No assumption that the reduced conditional mean is linear is needed for this projection limit.

The estimand is \beta; the limit is \beta^\ast. The discrepancy is the product of how much Z contributes to the conditional mean and how much Z co-moves with X per unit of X variation. It vanishes if \gamma=0 or \operatorname{Cov}(X,Z)=0. The latter guarantees equality of these slopes; uncorrelatedness alone does not guarantee that the entire reduced conditional-mean or error-distribution model is correct. When Z is recorded, its contribution can be investigated directly. When it is unmeasured, this decomposition exposes what is missing from identification (Figure 2).

Figure 2: Omitted variable bias as a picture. Left: the richer conditional-mean model contains an X term (slope \beta) and a Z term (slope \gamma), and Z moves with X. Arrows show these model terms; they do not identify a causal graph. The figure’s “effect” labels refer here to conditional mean slopes. Right: regressing Y on X alone mixes the Z association into the X coefficient. The reported slope is the conditional slope plus \gamma\operatorname{Cov}(X,Z)/\operatorname{Var}(X). The numbers are those of the simulation below, with both bars on the same scale.

set.seed(1)
n <- 1e5
z <- rnorm(n)
x <- 0.6 * z + rnorm(n)            # Cov(X, Z) = 0.6, Var(X) = 1.36
y <- 1 * x + 2 * z + rnorm(n)      # beta = 1, gamma = 2
coef(lm(y ~ x))["x"]               # estimates the limit, about 1.88

       x 
1.885318 

1 + 2 * 0.6 / 1.36                 # the formula, from the true Cov(X, Z) and Var(X)

[1] 1.882353

The regression recovers the limit, about 1.88, which is well above the estimand \beta = 1.

3.1 Estimand, parameter, and limit in this example

This example shows how the parameter, the estimand, and the limit line up. In the richer model, \beta and \gamma are parameters, and the estimand we care about is how much the mean of Y changes per unit change in X at a fixed value of Z. That estimand is \beta. In the working model, the parameters are a and b. Fitting it by least squares gives \hat b, which converges to the limit \beta^\ast found above, so the question is how b relates to \beta.

The working model asserts that its error is unrelated to X. Writing the error as e = Y - a - bX, being uncorrelated with X means \operatorname{Cov}(X, Y) - b\operatorname{Var}(X) = 0, which gives b = \operatorname{Cov}(X, Y)/\operatorname{Var}(X). This is the expression for \beta^\ast at the start of the section, so the working model’s slope is

\begin{aligned} b&= \beta^\ast\\ &= \beta + \gamma\,\frac{\operatorname{Cov}(X, Z)}{\operatorname{Var}(X)}. \end{aligned}

The ratio \operatorname{Cov}(X,Z)/\operatorname{Var}(X) is the slope of Z’s linear projection on X. The reduced regression combines the conditional X slope \beta with the contribution from the associated Z. All of these are associational statements in this example.

In the jointly Gaussian simulation, the regression residual is jointly normal with X, so being uncorrelated also makes it independent, so the working model describes the (X, Y) data exactly, and \hat b estimates b without difficulty. The gap between b and \beta arises because b answers a different question, namely how much Y differs per unit of X if we ignore Z. The conditional slope at a fixed Z is a feature of the richer model, and the (X, Y) data cannot reveal it by themselves. Different richer models, with different values of \gamma and \operatorname{Cov}(X, Z), can produce the same (X, Y) distribution and therefore the same b, while having different values of \beta.

A working model can also be wrong about the data themselves. In the nonlinear example in optional exercise 4, no value of b makes the straight-line conditional-mean model true. Its least-squares limit is still a well-defined linear projection, the line minimizing population mean squared error. That is a pseudo-true value for the restricted conditional-mean model.

	Belongs to	Still meaningful if the model is wrong?	In this example
Parameter	a model	not necessarily	\beta in the richer model, b in the working model
Estimand	the population and your question	yes	\beta, the conditional mean slope at fixed Z
Limit	an estimator and the population	yes	\beta^\ast = \beta + \gamma\operatorname{Cov}(X,Z)/\operatorname{Var}(X)

The textbook order is to start with a parameter and ask whether \hat\theta is consistent for it. The order here is the reverse: start with the estimand, work out the limit, and compare. A parameter enters only if you have a model, and the model behind a procedure is often not the model that generated your data.

4 A running example: the pooled correlation

Now a version that comes up constantly in expression studies. Two groups, controls and patients, with proportions \pi and 1-\pi. Within each group the two genes have the same correlation \rho. The patients have higher mean expression of both genes, by \delta standard deviations. Here \pi is the control fraction in the population being pooled. In a deliberately sampled case–control cohort it can be the sampling mixture, rather than the natural control prevalence. Changing that mixture changes a pooled target even when neither within-group distribution changes.

The picture below (Figure 3) shows one sample of 20 controls and 80 patients, first without a shift and then with it.

Figure 3: The pooled-correlation example, drawn with one illustrative sample of 20 controls (open circles) and 80 patients (filled circles). The points are constructed so that the correlation within each group is exactly 0.5. Left: with no shift between the groups (\delta = 0), the line through everyone has the same slope as the line within each group, and the correlation is 0.5 either way. Right: the same sample with the patients shifted by \delta = 2 in both genes. Each group still has correlation 0.5 (blue lines), but the cloud as a whole is stretched along the diagonal, and the line through everyone (dashed) is steeper, with correlation \rho^\ast = 0.695. Both genes have equal variance, so the slope of each line equals its correlation, and the bars compare the two correlations. The patients outnumber the controls four to one, so the line through everyone passes much closer to the patients’ center than to the controls’.

The estimand most people have in mind is \rho, the within-group correlation. It is also the parameter of a simple model: each group is bivariate normal with correlation \rho. The procedure most people run is cor(x, y) on everyone. What is its limit?

Write the pooled distribution as a mixture F = \pi F_0 + (1 - \pi) F_1, where F_0 is the distribution in the controls and F_1 in the patients. Each of these is a distribution for the pair of genes (gene A, gene B), so the group means \mu_0 and \mu_1 are vectors with one entry per gene, and the group covariance matrices \Sigma_0 and \Sigma_1 are 2 \times 2. The mean of a mixture is the weighted mean, \mu = \pi\mu_0 + (1-\pi)\mu_1. The covariance of a mixture is the weighted within-group covariance plus a between-group term:

\begin{aligned} \Delta_\mu&=\mu_1-\mu_0,\\ \Sigma&=\pi\Sigma_0+(1-\pi)\Sigma_1 +\pi(1-\pi)\Delta_\mu\Delta_\mu^\top. \end{aligned}

This is the law of total covariance. To derive it, write the feature vector minus its pooled mean as a within-group residual plus its group’s displacement from the pooled mean. The residual has conditional mean zero, so the two cross products vanish after averaging. The residual outer product contributes the weighted within-group covariances; the mean-displacement outer product contributes the between-group term. Both kinds of variation are real components of the pooled population, but only the first answers the within-group question.

In our example both groups have unit variances and correlation \rho. The mean difference is (\delta,\delta), so its outer product has every entry equal to \delta^2. Set c=\pi(1-\pi)\delta^2 to keep the matrix readable:

\Sigma = \begin{pmatrix}1+c&\rho+c\\\rho+c&1+c\end{pmatrix}.

and dividing its off-diagonal entry (the covariance) by its diagonal entry (the common variance) gives the pooled correlation:

\begin{aligned} \rho^\ast&= \frac{\rho+c}{1+c},\\ c&=\pi(1-\pi)\delta^2. \end{aligned}

Using the target/limit distinction, the bias term is the difference between the limit and the estimand, \rho^\ast - \rho. To write it out, put \rho over the denominator of the formula above and subtract. The numerator becomes \rho + \pi(1-\pi)\delta^2 - \rho\bigl(1 + \pi(1-\pi)\delta^2\bigr), which simplifies to (1-\rho)\,\pi(1-\pi)\,\delta^2, so

\rho^\ast - \rho = \frac{(1-\rho)\,\pi(1-\pi)\,\delta^2}{1 + \pi(1-\pi)\delta^2}.

The limit is therefore the estimand plus the bias term:

\rho^\ast = \rho + \underbrace{\frac{(1-\rho)\,\pi(1-\pi)\,\delta^2}{1 + \pi(1-\pi)\delta^2}}_{\text{bias term}}.

For \rho < 1, the bias term is positive whenever \delta \neq 0 and both groups are present, it grows with \delta^2, and it is largest at \pi = 1/2. With \rho = 0.5, \delta = 2, and \pi = 0.2, the formula gives \rho^\ast = (0.5 + 0.64)/(1 + 0.64) = 0.695, so the bias term is 0.695 - 0.5 = 0.195. The limit \rho^\ast is the correlation of the mixture. It is also the pseudo-true value of the correlation if you fit a single bivariate normal to everyone, and it is a derived functional of the two-group model parameters, rather than the common within-group correlation parameter.

pooled_cor <- function(rho, delta, pi, delta_b = delta) {   # delta_b: the second gene's shift (default: same as the first)
  v <- pi * (1 - pi)
  (rho + v * delta * delta_b) / sqrt((1 + v * delta^2) * (1 + v * delta_b^2))
}
pooled_cor(rho = 0.5, delta = 2, pi = 0.2)

[1] 0.695122

set.seed(1)
SimulateCohort <- function(n_ctrl, n_pat, shift, rho = 0.5) {
  group <- rep(c("control", "patient"), c(n_ctrl, n_pat))
  n <- length(group)
  z1 <- rnorm(n)
  z2 <- rho * z1 + sqrt(1 - rho^2) * rnorm(n)
  data.frame(group = group,
             gene_a = z1 + shift * (group == "patient"),
             gene_b = z2 + shift * (group == "patient"))
}
d <- SimulateCohort(n_ctrl = 20000, n_pat = 80000, shift = 2)
cor(d$gene_a, d$gene_b)                                   # the limit

[1] 0.6936609

with(subset(d, group == "control"), cor(gene_a, gene_b))  # the estimand

[1] 0.4978741

The simulation agrees with the formula. That agreement is the point of the simulation: it checks that the algebra is right, and the algebra tells you what the simulation would never tell you on its own, namely exactly which quantities drive the bias.

Shifting both genes by the same amount is a special case. If the shifts are \delta_a and \delta_b (in within-group standard deviations), so that \mu_1-\mu_0=(\delta_a,\delta_b), put v=\pi(1-\pi). The pooled correlation is

\rho^\ast = \frac{\rho+v\delta_a\delta_b} {\sqrt{(1+v\delta_a^2)(1+v\delta_b^2)}}.

The numerator contains the product \delta_a\delta_b, so the direction of the distortion depends on whether the two genes move in the same direction between groups. Genes that shift by similar amounts in the same direction look more correlated than they are. When one shift is much larger than the other, the denominator grows faster than the numerator, and the pooled correlation can fall below \rho even though both shifts are positive. Genes that move in opposite directions look less correlated, and can look negatively correlated:

pooled_cor(rho = 0.5, delta = 2, delta_b = -2, pi = 0.2)   # the within-group correlation is +0.5

[1] -0.08536585

4.1 What the bias term buys you

Because the bias term is a formula, it supports three uses.

The first is a correction. The mean difference \mu_1 - \mu_0 enters the pooled covariance only through the last term of the law of total covariance, \pi(1-\pi)(\mu_1 - \mu_0)(\mu_1 - \mu_0)^\top. Without that term the pooled covariance would equal the within-group covariance, and the pooled correlation would equal \rho. So the bias term comes entirely from the mean difference. The formula for it shows the same thing: it increases with \delta^2 and is zero exactly when \delta = 0. If each group’s mean is subtracted from its own members before the correlation is computed, the last term is zero by construction and the bias term disappears. Within-group centering is a correction that follows directly from the bias formula.

d$a_c <- d$gene_a - ave(d$gene_a, d$group)
d$b_c <- d$gene_b - ave(d$gene_b, d$group)
cor(d$a_c, d$b_c)

[1] 0.4988346

The second is detection. The derivation of the bias term also tells us which hypothesis to test. In this section’s setting the bias term is zero exactly when the two groups have the same mean vector, so the hypothesis

\begin{aligned} H_0&:\ \mu_1-\mu_0=(0,0),\\ H_1&:\ \mu_1-\mu_0\ne(0,0) \end{aligned}

is the one that matters for whether a pooled correlation can be trusted as an estimate of \rho. Rejecting H_0 is a warning that the bias term is generally nonzero, and the pooled correlation should then be compared with the within-group correlation. The choice of test depends on the model we assume. Here each group is bivariate normal with the same covariance matrix, and under that model Hotelling’s T^2 is an appropriate test of H_0. In R, manova() with test = "Hotelling-Lawley" gives it for two groups, and testing each gene separately with a two-sample t test, the usual differential expression test, is a simpler alternative. The size of the distortion is estimated by the pooled sample correlation minus the correlation after within-group centering.

# the diagnostic: test H0 that the two mean vectors are equal
summary(manova(cbind(gene_a, gene_b) ~ group, data = d), test = "Hotelling-Lawley")

             Df Hotelling-Lawley approx F num Df den Df    Pr(>F)    
group         1          0.84865    42431      2  99997 < 2.2e-16 ***
Residuals 99998                                                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# the estimated bias term: pooled correlation minus the correlation after centering
cor(d$gene_a, d$gene_b) - cor(d$a_c, d$b_c)

[1] 0.1948263

In the simulated cohort the test rejects decisively, and the estimated bias term is 0.195, in line with the formula.

The third is design. The bias term vanishes at \delta = 0 and is largest at \pi = 1/2. In a different problem the dependence on \pi is reversed and the bias vanishes at balance (see the Jackknife pseudo-values page). Either way, the formula shows how the cohort’s makeup, meaning the group sizes and the size of the mean shift, controls the artifact.

5 The general recipe

When a procedure is applied outside the setting it was built for, these steps show what it is doing.

1. Write down the working model the procedure assumes, and a richer model that departs from it in one specific way. Pick one member of the richer model to play the role of F.
2. Name the estimand, \theta(F), as a functional of that distribution. Be concrete about which distribution: the control group’s, the patients’, or the conditional distribution given a covariate value.
3. Work out the limit, \theta^\ast, which is what the procedure converges to when the data come from F, with heterogeneity, unbalanced groups, estimated nuisances, or whatever the richer model contains. Often it is the same functional applied to a different distribution, such as the pooled one or an adjusted one.
4. Subtract the estimand from the limit. The difference \theta^\ast - \theta(F) is the bias term, and the aim is to write it in closed form. Mixtures, projections, and first-order expansions are the three tools that work most often.
5. Check the formula by simulation. Simulate from F, compute the procedure’s output at large n, and compare it with the limit. If they disagree, either the formula or the simulation code is wrong, so check both (the Simulation studies for methods research page says to include a cell whose answer you already know).
6. Read the formula. When is the bias term zero? What makes it large? When its ingredients are identified from observed data, a zero-bias condition can sometimes motivate a diagnostic or correction. Restrictions involving unobserved outcomes or confounders generally cannot be tested this way; a non-rejection does not certify absence of bias. Its dependence on the design, such as the group sizes, gives advice on how to collect data.

6 A second pattern: adjustment that leaves something behind

The limit can also move because an adjustment is incomplete. Suppose you want the association between X and Y after removing a nuisance Z, and you do this by regressing both on an estimate \hat Z and correlating the residuals. The working model treats \hat Z as good enough to stand in for Z. In the richer model, \hat Z captures only part of Z, and since only the estimate is available, the rest of Z is unobserved even though a proxy for Z is recorded.

To see what the residuals contain, write X = bZ + e_X and Y = cZ + e_Y, where the errors e_X and e_Y are unrelated to Z and to \hat Z, and \operatorname{Cov}(e_X, e_Y) is the association given Z that we want. Write Z = \kappa\hat Z + U, where U is what is left of Z after projecting it on \hat Z, so that U is uncorrelated with \hat Z. If the estimate is calibrated, \kappa = 1 and U = Z - \hat Z. Let r_X and r_Y be the residuals from regressing X and Y on \hat Z. They still contain U, and their covariance is

\operatorname{Cov}(r_X,r_Y)=\operatorname{Cov}(e_X,e_Y)+bc\operatorname{Var}(U).

The bias term, bc\operatorname{Var}(U), is the part of the nuisance the adjustment missed, weighted by how strongly the nuisance drives each variable (the product bc). If the missed part differs between two groups, the adjusted correlations differ between groups even when the true associations given Z are identical. This is the mechanism behind cell-composition adjustment in bulk expression data, where Z is composition and \hat Z comes from a deconvolution method.

This is a projection calculation with centered variables (or intercepts in both regressions). Substituting Z=\kappa\hat Z+U gives r_X=bU+e_X and r_Y=cU+e_Y. Their covariance has the desired error covariance, the bc\operatorname{Var}(U) term, and zero cross terms under the stated uncorrelatedness assumptions. For a residual correlation, the remaining U also enters both variances:

\begin{aligned} v_X^r&=\operatorname{Var}(e_X)+b^2\operatorname{Var}(U),\\ v_Y^r&=\operatorname{Var}(e_Y)+c^2\operatorname{Var}(U). \end{aligned}

Divide the residual covariance by \sqrt{v_X^rv_Y^r} to get the residual correlation. A noisy proxy \hat Z is different from an estimated regression that becomes consistent as the training sample grows. The proxy’s error can remain a property of the measurement process at arbitrarily large sample sizes. Splitting the sample does not remove that measurement error.

set.seed(2)
n <- 20000
z <- rnorm(n)                      # the nuisance
z_hat <- 0.8 * z + rnorm(n, sd = 0.6)   # an imperfect estimate of it
x <- 1.5 * z + rnorm(n)
y <- 1.5 * z + rnorm(n)            # X and Y are independent given Z
cor(resid(lm(x ~ z)), resid(lm(y ~ z)))        # adjust on the truth: about 0

[1] 0.004645802

cor(resid(lm(x ~ z_hat)), resid(lm(y ~ z_hat)))  # adjust on the estimate: not 0

[1] 0.4469891

In this example b = c = 1.5, \operatorname{Var}(Z) = 1, and \operatorname{Cor}(Z, \hat Z) = 0.8, so \operatorname{Var}(U) = 1 - 0.8^2 = 0.36. The bias term for the covariance is 2.25 \times 0.36 = 0.81. The association we want, \operatorname{Cov}(e_X, e_Y), is zero, so the residual covariance is 0.81. Each residual has variance 1 + 1.5^2 \times 0.36 = 1.81, so the residual correlation is 0.81 / 1.81 = 0.45, which is what the second line above returns. That is the limit of the adjusted correlation, whose estimand is 0. Here \hat Z is not calibrated, with \kappa = 0.8, so U is not simply Z - \hat Z, and using that would give a covariance of 0.90 instead of 0.81.

The same formula explains the claim about groups. Give two groups the same truth, but let the estimate of Z be noisier in the second:

set.seed(3)
Adj <- function(n, sd_hat) {      # X and Y are independent given Z; the estimate of Z has noise SD sd_hat
  z <- rnorm(n)
  z_hat <- 0.8 * z + rnorm(n, sd = sd_hat)
  x <- 1.5 * z + rnorm(n)
  y <- 1.5 * z + rnorm(n)
  cor(resid(lm(x ~ z_hat)), resid(lm(y ~ z_hat)))
}
Pred <- function(sd_hat) {        # b * c * Var(U) / (1 + b^2 * Var(U)), with Var(U) = 1 - Cov(Z, Zhat)^2 / Var(Zhat)
  var_u <- 1 - 0.8^2 / (0.8^2 + sd_hat^2)
  2.25 * var_u / (1 + 2.25 * var_u)
}
round(rbind(simulated = c(group_1 = Adj(1e5, 0.3), group_2 = Adj(1e5, 0.9)),
            predicted = c(group_1 = Pred(0.3), group_2 = Pred(0.9))), 3)

          group_1 group_2
simulated   0.219   0.557
predicted   0.217   0.557

Both groups have a true association given Z of zero. The adjusted correlations are about 0.22 and 0.56, matching the formula, because the part of Z that the estimate misses, \operatorname{Var}(U), differs between the groups.

7 Why this matters in practice

The three examples in this page fit one pattern:

	Omitted variable	Pooled correlation	Imperfect adjustment
Working model	Y = a + bX + e, with b read as the conditional X slope at fixed Z	everyone comes from one homogeneous population	adjusting for \hat Z removes the nuisance Z
Richer model	Y = \beta X + \gamma Z + \varepsilon, with Z related to X	two groups whose means differ by \delta	\hat Z is a noisy estimate of Z
Estimand	\beta	the within-group correlation \rho	\operatorname{Cov}(e_X, e_Y), the association given Z (zero in the example)
Limit	\beta^\ast, the slope of Y on X alone	\rho^\ast, the correlation of everyone pooled	the covariance of the adjusted residuals
Bias term	\gamma\operatorname{Cov}(X,Z)/\operatorname{Var}(X)	\frac{(1-\rho)\pi(1-\pi)\delta^2}{1 + \pi(1-\pi)\delta^2}	bc\operatorname{Var}(U)
Bias term is zero when	\gamma = 0 or \operatorname{Cov}(X, Z) = 0	\delta = 0	bc\operatorname{Var}(U) = 0, for example when \hat Z captures all of Z

Each column identifies conditions under which the displayed discrepancy vanishes. Whether those conditions are testable depends on which quantities are observed and identified.

Much methodological work starts by working out the limit. Someone uses a procedure (a sample-specific network, a composition adjustment, a permutation test) in a setting it was not built for (unbalanced groups, repeated measures, estimated nuisances). The paper is then: here is the working model and the departure from it, here is the estimand, here is the limit, here is the difference in closed form, here is the simulation that confirms it, here is the correction or the diagnostic, here is what it changes in real data.

When you read a methods paper in this area, look for these seven pieces. In well-built papers they are almost always there.

8 The expansion applied to the wrong distribution

For a differentiable functional, write T(G)=T(F)+G\phi_F+R_2(G,F), with F\phi_F=0. This applies the local expansion from Influence functions to a different subgroup distribution; a larger cohort does not force that distributional change to be small.

SdN <- function(v) sqrt(mean((v - mean(v))^2))

Now the part that matters for pseudo-value regression. Suppose the data are a mixture F = \pi F_0 + (1-\pi) F_1 of controls (a fraction \pi) and patients, and the estimand is T(F_1), the patients’ correlation. The pseudo-values are computed at the pooled F_N \approx F. By the expansion with G = F_1,

T(F_1)=T(F) +\int\operatorname{IF}(x;T,F)\,dF_1(x) +R_2(F_1,F).

The average pseudo-value over the patients converges to the limit T(F) + \int \operatorname{IF}(x; T, F)\,dF_1(x), which is the first-order approximation to the estimand T(F_1) around the pooled distribution. The gap between the estimand and the limit is the remainder R_2(F_1, F). It does not shrink as N grows, because it depends on how far F_1 is from F, a property of the population that a larger sample cannot change.

There are two different approximations here. Estimating a fixed pooled F more accurately makes sampling error smaller. Moving from F to a fixed subgroup distribution F_1 is a population change whose size does not decrease with N. Local differentiability gives a tangent at F; it does not say that tangent must be accurate all the way to F_1. This differs from a nuisance remainder that vanishes because the fitted regression approaches its truth.

The picture below (Figure 4) shows the expansion for the two groups of this example. The pooled sample is the starting point, and the paths lead toward each group’s own distribution.

Figure 4: Where the first-order approximation goes wrong. Both paths start at the pooled sample (20% controls and 80% patients, a shift of 2 in both genes, within-group correlation 0.5) and move toward one group’s own distribution. The blue curves are the correlation along each path, which is the truth at every point, and they end at the estimand, 0.5. The dashed vermilion lines are the tangents at the start, the first-order approximation. A tangent’s slope is the group’s average influence function, and its value at the far end is the limit of the group’s average pseudo-value. The vertical gap at the far end is the bias term, which is minus the remainder R_2. It is larger for the controls because their own distribution is farther from the pooled distribution. The curves are exact for the population mixture, and the limits are 1.05 for the controls and 0.61 for the patients.

For a linear functional, T(F) = \int h\,dF, the remainder is zero, so the limit equals the estimand. For correlation the remainder can be nonzero when group means or variances differ; the discrepancy depends on their actual moment configuration and can cancel in particular cases. In the What is your estimator actually estimating? page the bias term is the limit minus the estimand, so it is minus this remainder: that bias term is now identified as the second-order remainder of a von Mises expansion. Pulling it out tells you the correction (center or standardize within group so that F_1 and F agree in the moments that drive R_2) and the diagnostic (check whether the group difference in pseudo-values tracks the mean difference).

In the setting of this section, with equal unit variances, within-group correlation \rho in both groups, and means that differ by \delta in both genes, the difference between the two remainders has a closed form. Again put c=\pi(1-\pi)\delta^2:

R_2(F_1,F)-R_2(F_0,F)=-(1-\rho)\frac{\delta^2(2\pi-1)}{(1+c)^2}.

It has the sign of 1 - 2\pi, so it changes sign when the group sizes are swapped, and it vanishes at balance. A variance difference between the groups adds a separate issue that centering does not remove (optional exercise 1).

set.seed(4)
n_ctrl <- 2000; n_pat <- 8000; rho <- 0.5; shift <- 2
grp <- rep(0:1, c(n_ctrl, n_pat))
z1 <- rnorm(n_ctrl + n_pat); z2 <- rho * z1 + sqrt(1 - rho^2) * rnorm(n_ctrl + n_pat)
a <- z1 + shift * grp; b <- z2 + shift * grp
r_pool <- cor(a, b)
IFcor <- function(x, y) { r <- cor(x, y); xt <- (x - mean(x)) / SdN(x); yt <- (y - mean(y)) / SdN(y); xt * yt - r / 2 * (xt^2 + yt^2) }
first_order <- r_pool + tapply(IFcor(a, b), grp, mean)   # what pseudo-value means estimate
truth <- c(cor(a[grp == 0], b[grp == 0]), cor(a[grp == 1], b[grp == 1]))
rbind(first_order = first_order, truth = truth, remainder_R2 = truth - first_order)

                      0          1
first_order   1.0444850  0.5977956
truth         0.4669724  0.4922634
remainder_R2 -0.5775127 -0.1055322

In this equal-variance example, the population remainder is zero when the mean shift is zero; the displayed empirical remainder still fluctuates. With a fixed nonzero shift, a larger sample estimates the persistent discrepancy more precisely rather than removing it. It differs between the two groups and is larger for the smaller group, here -0.58 for the controls and -0.11 for the patients. A between-group comparison of average pseudo-values converges to the difference of the two limits, which misses the comparison of estimands, T(F_1) - T(F_0), by minus the difference between the two remainders. At balance each group’s remainder is still far from zero (try n_ctrl <- n_pat <- 5000, where both are about -0.25), so a between-group comparison can have the right limit while each group’s own limit differs from its estimand. That distinction, contrast versus level, matters whenever the quantity of interest is a comparison between groups.

9 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 · DES-01 — Separate the scientific target from a pooled association

Task type: Derivation and interpretation.

Suppose group G is 0 or 1 with P(G=0)=\pi, retaining the control-fraction convention from the worked example. Within each group, X and Y have variance 1 and correlation \rho. Their means are (0,0) in group 0 and (\delta_X,\delta_Y) in group 1.

Use total covariance and total variance to derive the pooled correlation. Identify the within-group target and the probability limit of a pooled sample correlation. Explain when their difference represents a mismatch with the research question.

Ready to continue when: you can name the population target and derive the procedure’s limit without treating a pooled quantity as intrinsically wrong.

Transfer: Let the shifts have opposite signs. Can the pooled association have the opposite sign from the within-group association?

Hint

The between-group covariance is \operatorname{Cov}(\delta_XG,\delta_YG), and \operatorname{Var}(G)=\pi(1-\pi).

Solution outline — open after trying

The pooled covariance is \rho+\pi(1-\pi)\delta_X\delta_Y and its marginal variances are 1+\pi(1-\pi)\delta_X^2 and 1+\pi(1-\pi)\delta_Y^2. Their normalized ratio is the pooled correlation. A pooled sample correlation consistently targets that quantity under suitable sampling/moment conditions. If the scientific question is within-group correlation, the target is \rho and the pooled limit may differ. If the pooled population correlation is the intended target, heterogeneity alone does not make that estimator inconsistent.

Additional practice (optional)

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

1. Extend pooled_cor() to the case where the groups have different within-group variances, and check it by simulation. Try a variance difference in one gene only, and in both genes equally, with no mean shift. Which of the two artifacts, the mean shift or the variance difference, does within-group centering remove, and what would remove the other?
2. In the omitted-variable example, what happens to the bias if Z is measured with error and you adjust for the noisy version? Write the limit and check it by simulation. This is the “adjustment leaves something behind” pattern in its simplest form.
3. Use pooled_cor() with delta_b = -delta to plot the pooled correlation against \pi for a few values of \delta. At what \pi is it most negative, and for which \delta does it change sign? Check one point by simulation.
4. Let X \sim N(1, 1) and Y = X^2 + \varepsilon with E[\varepsilon\mid X]=0. A linear conditional-mean model is wrong, but a linear projection coefficient is still a valid population target. Work out the limit of the least-squares slope by hand (it is \operatorname{Cov}(X,Y)/\operatorname{Var}(X)), check it by simulation, and say what question it answers.

10 Reading

• Any regression text’s treatment of omitted variable bias; Wooldridge, Introductory Econometrics: A Modern Approach, does it cleanly.
• Hernán and Robins, Causal Inference: What If, Chapters 1 to 3, for a careful treatment of “what are we estimating” that applies well beyond causal inference.
• Vansteelandt, S. and Dukes, O. (2022). Assumption-lean inference for generalised linear model parameters. Journal of the Royal Statistical Society B, 84(3), 657 to 685. Sections 1 and 2 ask, of the ordinary regression coefficient, exactly the question this page asks, and answer it by defining an estimand that stays meaningful when the model is wrong. The rest of the paper constructs an estimator and its inferential justification.
• White, H. (1982). Maximum likelihood estimation of misspecified models. Econometrica, 50(1), 1 to 25. The standard reference for what a maximum likelihood estimator converges to when its model is wrong, which is the pseudo-true value.