Variance Estimation Under Case-Control Sampling

Why Variance Estimation Needs the Sampling Design

Inverse-probability weighting can make a case-control estimating equation target a population quantity, but that does not by itself determine a valid standard error. The uncertainty depends both on how the weighted equation changes with the coefficient and on how the study selected cases and controls. Treating the weighted equation as an ordinary unweighted likelihood can therefore produce an uncertainty estimate that does not match the sampling process.

This tutorial focuses on weighted logistic estimating equations when cases and controls are sampled separately. The key distinction is between the bread, which describes the equation's local sensitivity to the coefficient, and the meat, which describes variability in its score under the relevant sampling scheme. A design-aware calculation uses the actual sampling strata and, for sampling without replacement from finite strata, a finite-population correction.

Definition: Let $Y_i\in\{0,1\}$, let $x_i$ be a column vector of covariates (including an intercept if required), and let $\mu_i(\beta)=\operatorname{expit}(x_i^\mathsf{T}\beta)$. The logistic score contribution is $g_i(\beta)=x_i\{Y_i-\mu_i(\beta)\}$. For sampled units with inverse-probability weights $w_i$, the weighted estimating equation is $U(\beta)=\sum_{i\in s}w_i g_i(\beta)=0$. Its negative derivative is $A(\beta)=\sum_{i\in s}w_i\mu_i(\beta)\{1-\mu_i(\beta)\}x_ix_i^\mathsf{T}$.

Here $s$ denotes the observed sample. In an outcome-stratified case-control design, the sampling strata are often $h=1$ for cases and $h=0$ for controls. More generally, there may be additional sampling strata. The formulas below first assume that, within each stratum $h$, a simple random sample without replacement of size $n_h$ is drawn from a finite set of $N_h$ units. The inclusion probability is then $n_h/N_h$, and the common inverse-probability weight is $d_h=N_h/n_h$.

Model-Based and Sandwich Uncertainty

The inverse-Hessian, or model-based, covariance estimate uses $A^{-1}$ as the coefficient covariance. It is appropriate only when the likelihood and its information calculation match the data-generating and sampling assumptions being used. In a weighted estimating equation, the weights alter the score's variability as well as its derivative; using the inverse weighted Hessian alone does not generally account for that variability.

A sandwich estimate instead has the form $A^{-1}BA^{-\mathsf{T}}$, where $B$ is the covariance of the estimating equation. If independent sampled units are treated as the relevant sampling units, a common empirical meat is $\sum_{i\in s}w_i^2g_i(\hat\beta)g_i(\hat\beta)^\mathsf{T}$. This is a robust, independent-unit estimate, but it does not automatically account for finite-population sampling without replacement. When the design is known, the meat should reflect that design.

Key distinction: The weighted Hessian $A$ describes sensitivity of the estimating equation. The meat $B$ describes score variability. The inverse-Hessian covariance $A^{-1}$ effectively substitutes a particular information identity for $B$; a sandwich or design-based covariance estimates $B$ separately.

A design-based formulation starts with the finite population of units in each sampling stratum. At a fixed coefficient $\beta$, define the finite-population score total $T_h(\beta)=\sum_{i\in h}g_i(\beta)$. The inverse-probability-weighted sample estimates it by $\hat T_h(\beta)=d_h\sum_{i\in s_h}g_i(\beta)$. Because the case and control sample sizes are fixed by design, these totals are estimated separately within the strata.

Theorem: Under independent simple random sampling without replacement within strata, the design covariance of the weighted score total at a fixed $\beta$ is $$ \operatorname{Var}_d\{\hat T(\beta)\} =\sum_h \frac{N_h^2(1-f_h)}{n_h}S_{g,h}(\beta), \qquad f_h=\frac{n_h}{N_h}, $$ where $\hat T=\sum_h\hat T_h$, $S_{g,h}=(N_h-1)^{-1}\sum_{i\in h}\{g_i-\bar g_h\}\{g_i-\bar g_h\}^\mathsf{T}$, and $\operatorname{Var}_d$ is variance over the sampling design, conditional on the finite population.

Proof. Within stratum $h$, write $I_i=1$ if unit $i$ is sampled and $I_i=0$ otherwise. Under simple random sampling without replacement, $E_d[I_i]=f_h$, $\operatorname{Var}_d(I_i)=f_h(1-f_h)$, and for distinct units $i,j$, $\operatorname{Cov}_d(I_i,I_j)=-f_h(1-f_h)/(N_h-1)$. Since $d_h=1/f_h$,

\operatorname{Var}_d\{\hat T_h\} =d_h^2\operatorname{Var}_d\left(\sum_{i\in h}I_i g_i\right) =\frac{N_h^2(1-f_h)}{n_h}S_{g,h}.

The last equality follows by expanding the covariance of the sum using the stated variances and covariances of the indicators, then collecting the centered outer products $(g_i-\bar g_h)(g_i-\bar g_h)^\mathsf{T}$. Samples in different strata are independent, so their covariance contributions add. This proves the result. $\square$

In practice, the finite-population score values are not all observed. Estimate $S_{g,h}$ by the sample covariance $s_{g,h}$ of the fitted score contributions in stratum $h$. The resulting estimated meat is $$ \hat B_d=\sum_h \frac{N_h^2(1-f_h)}{n_h}s_{g,h}(\hat\beta). $$ If $A$ estimates the negative derivative of the target equation, a first-order Taylor expansion gives the coefficient covariance estimate $A^{-1}\hat B_d A^{-\mathsf{T}}$. This linearization requires a nonsingular limiting derivative and an estimator close enough to the target root for the first-order approximation to be informative. The formula also presumes the stated stratified simple-random-sampling design; unequal-probability or clustered designs require their own design variance calculation.

Worked Design-Based Calculations

Worked Example: Two Outcome Strata

Suppose a case-control study samples $n_1=10$ of $N_1=100$ cases and $n_0=40$ of $N_0=400$ controls, by simple random sampling without replacement within each outcome stratum. At a fixed fitted coefficient, suppose the sample score variances for one scalar component are $s_{g,1}^2=0.04$ among cases and $s_{g,0}^2=0.01$ among controls. Calculate the estimated design variance of the weighted score total.

For cases, $f_1=10/100=0.10$, so the variance contribution is $$ \frac{100^2(1-0.10)}{10}(0.04) =\frac{10{,}000(0.90)}{10}(0.04) =900(0.04)=36. $$

For controls, $f_0=40/400=0.10$, so the contribution is $$ \frac{400^2(1-0.10)}{40}(0.01) =\frac{160{,}000(0.90)}{40}(0.01) =3{,}600(0.01)=36. $$

The two independently sampled strata contribute additively, giving $\hat B_d=36+36=72$ for this score component. If the scalar derivative magnitude is $A=12$, the linearized coefficient variance is $72/12^2=72/144=0.5$, and the standard error is $\sqrt{0.5}\approx0.707$. The score variance and the coefficient variance are different quantities: the latter also depends on the equation's derivative.

Worked Example: How the Sampling Fraction Changes the Correction

Now suppose a case stratum has $N_1=100$, $n_1=20$, and score variance $0.04$, while the control stratum has $N_0=400$, $n_0=80$, and score variance $0.01$. Both sampling fractions are $0.20$. The case contribution is $$ \frac{100^2(1-0.20)}{20}(0.04) =\frac{10{,}000(0.80)}{20}(0.04) =400(0.04)=16. $$

The control contribution is $$ \frac{400^2(1-0.20)}{80}(0.01) =\frac{160{,}000(0.80)}{80}(0.01) =1{,}600(0.01)=16. $$

Thus the estimated design variance of this score component is $16+16=32$. For comparison, omitting the finite-population correction would give case contribution $500(0.04)=20$ and control contribution $2{,}000(0.01)=20$, totaling $40$. Sampling a larger fraction of each finite stratum reduces design variance; the correction is $1-f_h=0.80$ in both strata here.

From Score Variability to Coefficient Variability

Theorem: Suppose the weighted estimating equation has a target root $\beta_0$, its negative derivative at that root is nonsingular, and the score-total estimation error is small enough for a first-order Taylor expansion. Then $$ \hat\beta-\beta_0\ \approx\ A^{-1}\{\hat T(\beta_0)-T(\beta_0)\}, \qquad \operatorname{Var}_d(\hat\beta)\ \approx\ A^{-1}\operatorname{Var}_d\{\hat T(\beta_0)\}A^{-\mathsf{T}}. $$

Derivation. Define the estimated equation $\hat T(\beta)=0$ and let $T(\beta_0)=0$ be its target-population counterpart. A first-order expansion of the estimated equation around $\beta_0$ gives $$ 0=\hat T(\hat\beta)\approx \hat T(\beta_0)-A(\hat\beta-\beta_0), $$ where $A$ is the negative derivative of the target equation, estimated in applications by the corresponding weighted derivative. Rearranging gives $\hat\beta-\beta_0\approx A^{-1}\hat T(\beta_0)$, since $T(\beta_0)=0$. Taking design covariance on both sides yields the displayed sandwich expression. This is a linearization, not an exact finite-sample identity; its quality depends on the regularity of the estimating equation and the amount of information in the sample. $\square$

Worked Example: Propagating a Two-Parameter Meat

Suppose a fitted two-coefficient model has estimated bread $$ A=\begin{pmatrix}10&0\\0&5\end{pmatrix}, \qquad \hat B_d=\begin{pmatrix}72&12\\12&48\end{pmatrix}. $$ Because $A^{-1}=\begin{pmatrix}0.1&0\\0&0.2\end{pmatrix}$, the sandwich covariance is $$ A^{-1}\hat B_d A^{-\mathsf{T}} = \begin{pmatrix}0.1&0\\0&0.2\end{pmatrix} \begin{pmatrix}72&12\\12&48\end{pmatrix} \begin{pmatrix}0.1&0\\0&0.2\end{pmatrix} = \begin{pmatrix}0.72&0.24\\0.24&1.92\end{pmatrix}. $$

The standard errors are $\sqrt{0.72}\approx0.849$ and $\sqrt{1.92}\approx1.386$. The off-diagonal covariance is $0.24$. Applying only the inverse bread would instead report diagonal entries $0.1$ and $0.2$, which are not the design-based variances in this example.

Bootstrap Resampling Within Sampling Strata

A bootstrap for a case-control sample must respect the way the observations entered the study. To preserve the fixed case and control sample sizes, resample cases from observed cases and controls from observed controls separately. In each replicate, draw $n_1$ cases with replacement from the observed case records and $n_0$ controls with replacement from the observed control records, retain the appropriate weights, and refit the weighted equation. The empirical covariance of the replicate coefficient estimates estimates the variability represented by that bootstrap scheme.

That basic stratified bootstrap is a with-replacement approximation to sampling within the observed strata. It does not, by itself, reproduce the finite-population correction for simple random sampling without replacement. When that correction matters, use a survey bootstrap or replicate-weight procedure calibrated to the design, or use the design-based sandwich calculation directly. Resampling all observations together is also inappropriate when it allows the case-control composition to fluctuate despite a design that fixed the stratum sample sizes.

Worked Example: A Stratified Bootstrap Plan

Suppose the observed data contain 20 sampled cases and 80 sampled controls, with inverse-probability weights known for both strata. To generate one ordinary stratified bootstrap replicate, draw 20 case records with replacement from the 20 observed cases, then draw 80 control records with replacement from the 80 observed controls. Refit the weighted logistic estimating equation using the case and control weights associated with the resampled records. Repeat this procedure $B$ times, obtaining estimates $\hat\beta^{*(1)},\ldots,\hat\beta^{*(B)}$.

For a coefficient component $j$, its bootstrap variance is $$ \widehat{\operatorname{Var}}_{\mathrm{boot}}(\hat\beta_j) =\frac{1}{B-1}\sum_{b=1}^B \left(\hat\beta_j^{*(b)}-\bar\beta_j^*\right)^2, \qquad \bar\beta_j^*=\frac{1}{B}\sum_{b=1}^B\hat\beta_j^{*(b)}. $$

This calculation preserves the observed stratum sizes in every replicate. If the study sampled substantial fractions of known finite case and control populations, however, this ordinary bootstrap variance should not be interpreted as including the finite-population reduction. The resampling method must be chosen to match the intended source of uncertainty.

Common Pitfalls and Practical Checks

Using the inverse Hessian as if weights were ordinary likelihood counts. The derivative supplies the bread, not necessarily the score covariance. Check whether the assumed information identity is justified for the weighted sampling procedure.
Confusing $N_h$, $n_h$, and the weight. Under equal-probability sampling in stratum $h$, $d_h=N_h/n_h$. The design variance factor is $N_h^2(1-f_h)/n_h$, not a factor obtained by substituting the sample size for the population size.
Dropping the finite-population correction automatically. The factor $1-f_h$ applies to simple random sampling without replacement from a finite stratum. It is not a universal multiplier for every sampling scheme.
Bootstrapping cases and controls as one pool. If the design fixed their sample sizes, resampling separately within the outcome strata preserves that feature. A basic with-replacement bootstrap still does not automatically reproduce finite-population sampling.
Reporting score variance as coefficient variance. The meat is on the scale of the estimating equation. It must be propagated through the bread on both sides.

Key takeaway: For weighted case-control fits, uncertainty depends on both the weighted equation's derivative and the variability induced by sampling. Build the meat from the actual design—separately within sampling strata—and then use the sandwich to translate score variability into coefficient variability.

Check Your Understanding

Answer these questions using the definitions and calculations above.

Why does the inverse weighted Hessian alone not generally give the design-based covariance of an inverse-probability-weighted fit?
Under stratified simple random sampling without replacement, what does the factor $1-f_h$ represent?
A stratum has $N_h=200$, $n_h=20$, and a scalar score variance of $0.03$. Calculate its estimated design variance contribution.
Why should a bootstrap preserve the case and control sample sizes when those sizes were fixed by design?
Does an ordinary with-replacement bootstrap within outcome strata automatically include a finite-population correction? Explain.