Two-Phase Sampling For GLMs

Why Two-Phase Sampling Changes the Estimating Equation

In a two-phase study, researchers first screen a relatively large set of units and record inexpensive information on them. They then collect a response, detailed covariates, or both for a selected subset. The phase-two sample is therefore not simply a smaller copy of the population: selection may depend on screening information, and some units may have much higher chances of receiving the detailed measurements than others.

A generalized linear model for the population is defined using the variables of scientific interest, not merely the units that happen to receive phase-two measurements. If the phase-two inclusion chances differ, the ordinary unweighted score generally targets the selected sample's distribution. Inverse-probability weighting instead gives each observed unit a contribution scaled to represent the population units it could stand for.

Two-Phase Inclusion Probabilities

Let $U=\{1,\ldots,N\}$ be a finite target population. Let $I_i$ equal 1 if unit $i$ is included at phase one, and 0 otherwise. Phase one records screening information $Z_i$ for each selected unit. Let $J_i$ equal 1 if unit $i$ is included at phase two. The phases are nested, so $J_i\leq I_i$. Phase two collects the response and covariates needed for the GLM.

Definition: The phase-one inclusion probability is $\pi_{i1}=P(I_i=1)$. Given the phase-one sample and its observed screening information, write $q_i(S_1)=P(J_i=1\mid S_1)$ for the conditional phase-two inclusion probability of a phase-one unit. The overall phase-two inclusion probability is $\pi_i=P(J_i=1)$. For distinct units $i,j$, the pairwise phase-two inclusion probability is $\pi_{ij}=P(J_i=1,J_j=1)$; use $\pi_{ii}=\pi_i$.

If $q_i(S_1)$ is a fixed probability for unit $i$ whenever it enters phase one, then the stages combine simply:

\pi_i=\pi_{i1}q_i.

More generally, the conditional phase-two probability can depend on the realized phase-one sample. Then the marginal probability is

\pi_i = \pi_{i1} E\!\left[q_i(S_1)\mid I_i=1\right].

This distinction matters when calculating weights: a conditional probability for the realized phase-one sample is not automatically the same thing as a marginal inclusion probability. In either representation, units that can contribute to the target score need positive inclusion probability. In particular, assume $\pi_i>0$ for every population unit whose score is to be represented.

The pairwise probabilities describe how units enter together. In general, $\pi_{ij}$ is not $\pi_i\pi_j$: fixed sample sizes, stratification, and other design features can make inclusion indicators dependent. The joint probabilities are important not only for describing the design but also for quantifying the sampling variability of weighted estimates.

From a Binary-Response Model to a Weighted Score

Suppose the population response $Y_i$ is binary and the covariate vector $x_i$ includes an intercept. For a logistic GLM, specify

P(Y_i=1\mid x_i)=\mu_i(\beta), \qquad \operatorname{logit}\{\mu_i(\beta)\}=x_i^{T}\beta, \qquad \mu_i(\beta)=\frac{\exp(x_i^{T}\beta)}{1+\exp(x_i^{T}\beta)}.

The Bernoulli log-likelihood contribution for a population unit is $Y_i\log\{\mu_i(\beta)\}+(1-Y_i)\log\{1-\mu_i(\beta)\}$. Differentiating with respect to $\beta$ gives the unit score

u_i(\beta)=x_i\{Y_i-\mu_i(\beta)\}.

Indeed, the derivative of $\mu_i$ with respect to $\beta$ is $\mu_i(1-\mu_i)x_i$, while the derivative of the Bernoulli log-likelihood with respect to $\mu_i$ is $(Y_i-\mu_i)/\{\mu_i(1-\mu_i)\}$; their product is the displayed score.

If all population outcomes and covariates were measured, the likelihood score would be $U_N(\beta)=\sum_{i\in U}u_i(\beta)$, and the usual score equation would set it to zero. Under two-phase sampling, $u_i(\beta)$ is only available for phase-two units. The inverse-probability-weighted score is

U_w(\beta) = \sum_{i\in U:J_i=1}\frac{u_i(\beta)}{\pi_i} = \sum_{i\in U:J_i=1}\frac{x_i\{Y_i-\mu_i(\beta)\}}{\pi_i}.

The weighted estimator solves $U_w(\beta)=0$. The weight $1/\pi_i$ makes a rarely observed unit represent more population units than a commonly observed one. If the two stages have fixed probabilities with $\pi_i=\pi_{i1}q_i$, the weight is $1/(\pi_{i1}q_i)$, not merely $1/q_i$: phase-one selection also affects representation in the population.

Definition: The inverse-probability-weighted logistic estimator is a solution $\hat\beta_w$ of $$ \sum_{i:J_i=1}\frac{x_i\{Y_i-\mu_i(\beta)\}}{\pi_i}=0. $$ It is an estimating-equation estimator. Equivalently, it maximizes the Bernoulli log-likelihood with each observed contribution multiplied by $1/\pi_i$; that weighted objective is not, in general, the likelihood of the phase-two sample under its sampling design.

Worked Example: Combining the Two Sampling Phases

Suppose every unit has phase-one inclusion probability $\pi_{i1}=0.5$. The screening variable $Z$ has two values. Conditional on phase one, units with $Z=0$ have phase-two inclusion probability $q_i=0.5$, while those with $Z=1$ have $q_i=0.25$. For a phase-two unit with $Z=0$, the overall inclusion probability is $0.5(0.5)=0.25$, so its weight is $1/0.25=4$. For a phase-two unit with $Z=1$, it is $0.5(0.25)=0.125$, so its weight is $1/0.125=8$.

Consider four observed phase-two records, with $Z$ used as the single covariate and an intercept in the model: two records have $Z=0$ and responses 0 and 1; two have $Z=1$ and responses 1 and 1. Evaluate the weighted score at $\beta=(0,0)^T$. Then $\mu_i=1/2$ for all four records, so each residual $Y_i-\mu_i$ is $-1/2$ for a response of 0 and $1/2$ for a response of 1.

The two $Z=0$ records each have weight 4. Their combined intercept-score contribution is $4(-1/2)+4(1/2)=0$, and their slope-score contribution is 0 because $Z=0$. The two $Z=1$ records each have weight 8. Their combined intercept contribution is $8(1/2)+8(1/2)=8$, and their slope contribution is $1$ times that amount, also 8. Thus $U_w((0,0)^T)=(8,8)^T$. This calculation is a weighted score evaluation, not a claim that the unknown population score must equal this realized value.

If only the phase-two probability had been used, the weights would have been 2 for $Z=0$ and 4 for $Z=1$. Those weights would represent the phase-one sample, not the full target population, because they omit the factor for phase-one selection.

Why the Weighted Score Targets the Population Score

The key property is design unbiasedness for the score at any fixed $\beta$. Here, design expectation means expectation over the sampling design while holding the finite-population responses and covariates fixed. It does not require the GLM to be correctly specified.

Theorem: Suppose $u_i(\beta)$ is fixed for every $i\in U$, and the phase-two inclusion probability $\pi_i$ is known and positive. Then $$ E_d\!\left[ \sum_{i:J_i=1}\frac{u_i(\beta)}{\pi_i} \right] = \sum_{i\in U}u_i(\beta). $$

Proof. Write $J_i$ for the phase-two inclusion indicator. Since $E_d[J_i]=P(J_i=1)=\pi_i$, and $u_i(\beta)$ and $\pi_i$ are fixed under the design expectation,

\begin{aligned} E_d\!\left[\sum_{i\in U}\frac{J_i u_i(\beta)}{\pi_i}\right] &=\sum_{i\in U}\frac{E_d[J_i]u_i(\beta)}{\pi_i}\\ &=\sum_{i\in U}\frac{\pi_i u_i(\beta)}{\pi_i}\\ &=\sum_{i\in U}u_i(\beta). \end{aligned}

The sum is finite, so taking expectation term by term requires no interchange-of-limit argument. The result applies component by component when $u_i(\beta)$ is a vector. $\square$

This theorem says the weighted score recovers the full-population score in design expectation. It does not say that a root of the weighted score is exactly unbiased for a finite-population logistic coefficient: solving a random nonlinear equation is a separate operation. Under suitable regularity conditions, including positive inclusion probabilities, a well-defined population target, and design conditions that make the normalized weighted score concentrate around its population counterpart, the resulting estimator can be consistent for that target.

Pairwise Inclusion and Score Variance

Although the weighted score is unbiased, its variance depends on the dependence between inclusion indicators. Define $\Delta_{ij}=\pi_{ij}-\pi_i\pi_j$, where $\pi_{ii}=\pi_i$. Treating each score as a column vector, the design covariance matrix is

\operatorname{Var}_d\{U_w(\beta)\} = \sum_{i\in U}\sum_{j\in U} \frac{\Delta_{ij}}{\pi_i\pi_j} u_i(\beta)u_j(\beta)^T.

Derivation. By the theorem, $E_d\{U_w(\beta)\}=\sum_i u_i(\beta)$. Subtract this mean from the random score to obtain

U_w(\beta)-E_d\{U_w(\beta)\} = \sum_{i\in U}\frac{J_i-\pi_i}{\pi_i}u_i(\beta).

For any two units, $E_d[(J_i-\pi_i)(J_j-\pi_j)]=\operatorname{Cov}_d(J_i,J_j)=\pi_{ij}-\pi_i\pi_j=\Delta_{ij}$. Expanding the covariance of the displayed sum therefore gives

\begin{aligned} \operatorname{Var}_d\{U_w(\beta)\} &=\sum_{i\in U}\sum_{j\in U} \frac{E_d[(J_i-\pi_i)(J_j-\pi_j)]}{\pi_i\pi_j} u_i(\beta)u_j(\beta)^T\\ &=\sum_{i\in U}\sum_{j\in U} \frac{\pi_{ij}-\pi_i\pi_j}{\pi_i\pi_j} u_i(\beta)u_j(\beta)^T. \end{aligned}

Thus marginal weights suffice for the weighted score, but marginal probabilities alone do not describe its design variance. Estimating uncertainty for a fitted coefficient also involves how the score changes with $\beta$, as well as the sampling variation of the score; the pairwise inclusion structure is one essential part of that calculation. If the design makes some relevant pairwise inclusion probabilities zero, a variance estimator based on observing those pairs needs additional care.

Worked Example: Pairwise Probabilities Under Fixed-Size Sampling

Suppose phase one is a census, and phase two selects exactly 2 of $N=4$ units by simple random sampling without replacement. Each of the $\binom{4}{2}=6$ possible pairs is equally likely. A particular unit appears in 3 of the 6 pairs, so $\pi_i=3/6=1/2$. Two distinct units appear together in exactly 1 pair, so $\pi_{ij}=1/6$. Their joint inclusion probability is not the product $1/4$.

For example, let the four fixed scalar score contributions at a chosen $\beta$ be $u_1=1$, $u_2=2$, $u_3=0$, and $u_4=-1$. The full score total is $1+2+0-1=2$. If the selected pair is $\{1,4\}$, the weighted total is $1/(1/2)+(-1)/(1/2)=2-2=0$. Across all six equally likely pairs, the weighted totals are $6,2,4,4,6,0$, which average to $22/6$—that list would be incorrect if calculated without checking the pair contributions. Calculating carefully, the six pair totals are: $\{1,2\}: 2(1+2)=6$; $\{1,3\}:2(1+0)=2$; $\{1,4\}:2(1-1)=0$; $\{2,3\}:2(2+0)=4$; $\{2,4\}:2(2-1)=2$; and $\{3,4\}:2(0-1)=-2$. Their average is $(6+2+0+4+2-2)/6=2$, as the theorem requires.

The pairwise probability $1/6$, rather than $1/4$, records the dependence induced by selecting a fixed-size sample. Replacing it with a product would give the wrong covariance terms for this design.

Worked Example: A Weighted Intercept-Only Logistic Fit

Suppose phase one is a census of 100 people. Phase two samples 20 people from screening group $Z=1$, where the group contains 40 people, and 6 from $Z=0$, where the group contains 60 people. Assume simple random sampling within each group. The phase-two inclusion probabilities are $20/40=0.5$ for $Z=1$ and $6/60=0.1$ for $Z=0$, so the weights are 2 and 10, respectively.

Suppose the phase-two sample contains 8 cases among the 20 sampled from $Z=1$, and 1 case among the 6 sampled from $Z=0$. For an intercept-only logistic model, the weighted score is $\sum_i w_i(Y_i-\mu)$. Its weighted total weight is $20(2)+6(10)=100$, and its weighted case total is $8(2)+1(10)=26$. Therefore the score equation is $26-100\mu=0$, giving $\hat\mu_w=0.26$. The fitted intercept is

\hat\beta_{0,w} = \log\!\left(\frac{0.26}{1-0.26}\right) = \log\!\left(\frac{13}{37}\right) \approx -1.046.

The unweighted sample prevalence is $9/26\approx0.346$, which differs because the two screening groups were sampled at different rates. The weighted fit corrects this imbalance for the intercept-only population response model. It does not, by itself, solve every issue of model specification or uncertainty estimation.

Common Mistakes and Practical Interpretation

Using only the phase-two probability. If phase one sampled from a larger population, the overall probability must include both stages. With fixed stage probabilities, use $\pi_i=\pi_{i1}q_i$.
Assuming pairwise inclusion probabilities factor. Fixed sample sizes and other dependent designs generally give $\pi_{ij}\ne\pi_i\pi_j$. Use the design's joint probabilities when assessing score variability.
Calling a weighted score an ordinary sample likelihood. The weighted Bernoulli objective is useful for estimation, but it is generally a pseudo-likelihood or estimating criterion, not the likelihood of the observed sample under the sampling design.
Claiming the coefficient estimate is unbiased because the score is unbiased. The result proved here concerns the score at a fixed $\beta$. A root of a random score equation need not be exactly unbiased in finite samples.
Ignoring positivity. A population unit with zero chance of reaching phase two cannot be represented by a finite inverse-probability weight. The weighted score cannot recover contributions from such units without additional assumptions or a changed target.

The central design principle is straightforward: describe how a unit reaches phase two, use its overall inclusion probability in the score, and retain the joint probabilities when considering uncertainty. This lets the GLM estimating equation refer to the target population even when detailed measurements are collected only on a selected subset.

Check Your Understanding

Work through each question by identifying the sampling probability or score contribution that is relevant.

A unit has phase-one inclusion probability 0.4 and conditional phase-two inclusion probability 0.5. What is its overall inclusion probability and inverse-probability weight?
Why does the weighted-score unbiasedness result not imply that the fitted coefficient is exactly unbiased in a finite sample?
In a fixed-size phase-two sample, why might $\pi_{ij}$ differ from $\pi_i\pi_j$?
For a logistic model, write the contribution of a phase-two unit with covariate vector $x_i$, response $Y_i$, and overall inclusion probability $\pi_i$ to the weighted score.
A stratum has no chance of being selected at phase two. Can inverse-probability weighting recover its score contribution? Explain.

What You'll Learn