Zhou, Haibo*; Chen, Jianwei†; Rissanen, Tiina H.‡; Korrick, Susan A.§; Hu, Howard§; Salonen, Jukka T.‡; Longnecker, Matthew P.¶

ArticlePlus

Click on the links below to access all the ArticlePlus for this article.

Please note that ArticlePlus files may launch a viewer application outside of your web browser.

* http://links.lww.com/EDE/A238

* http://links.lww.com/EDE/A239

The case–control study is a simple and familiar example of an outcome-dependent sampling design. The case–control design is logistically and economically appealing because observations made on a judiciously chosen subset of the population base provide nearly the statistical efficiency of observing the entire population base.¹ The principal idea of the case–control design and its subsequent extensions, such as case–cohort and two-stage designs,^2–7 is to concentrate resources on observations carrying the greatest amount of information. Related ideas of response-based sampling have also been developed in economics and survey sampling.^8–10

A unique property of the case–control design is that, under the logistic regression model with a binary Y, the estimated regression parameter (the odds ratio) for the exposure is the same under retrospective sampling as it is under a cohort study, with the effect of sampling confined to the intercept.^2,4 Thus, analyzing a data set from a case–control study as if it were from a cohort study will affect only the intercept. This property does not hold for a continuous response variable, or for regression models that are not logistic.

While designs for studying dichotomous outcomes have continued to develop, analogous work for studies of continuous outcomes has lagged. As the scope of epidemiologic inquiry grows, so does the need for efficient approaches to studying the determinants of the level of a continuous outcome, especially when the measurement of exposure is expensive.

As an example, we wanted to study in utero exposure to background-levels of the neurodevelopmental toxicant polychlorinated biphenyls (PCBs), and its relation to performance on the Bayley Scale of Infant Development. Maternal serum collected during pregnancy was available from a previously completed cohort study in which Bayley Scale had been measured. PCB concentration in the maternal serum provides a good surrogate measure of in utero exposure.

In studies of the relation between a relatively expensive exposure measure and the level of a continuous outcome, one approach has been to dichotomize the outcome and conduct a nested case–control study.^11,12 Dichotomizing a continuous outcome, however, will distort the estimand and usually will result in a loss of information because it uses a lower-order scale for the response.¹³

To reap the benefits of a reduced sample size, one can employ outcome-dependent sampling, which is a case–control-like sampling adapted for continuous outcomes. Zhou et al¹⁴ proposed a sampling scheme with 2 components. First, an overall random sample of the population base is taken. Second, one or more supplementary random samples are taken, in which the probability of selection depends on the level of the outcome variable (as with case–control sampling). For example, supplementary random samples might be taken only from the tails of the outcome distribution. In most settings, oversampling subjects in the tails achieves greater efficiency than a completely random sample of the same sample size, because the “tail” observations can provide greater influence on the parameter under study. Unlike the case–control study, the analysis must account for the biased sampling scheme with a continuous response to avoid biased estimates by the regression parameters (ie, to consistently estimate the same estimand as would have been obtained by studying the entire population). One commonly used method that can be adapted for this purpose is to conduct a weighted analysis (weighted estimating equation approach), with weights inversely proportional to the probability of being sampled (IPW).¹⁵ Another alternative is the weighted pseudo-likelihood method,¹⁰ which requires that one correctly specify all the underlying distributions. Misspecification of these distributions will lead to biased and erroneous conclusions. These methods, however, account for the sampling scheme in only an approximate way. A more efficient analysis, based on a likelihood function that precisely reflects the biased sampling design¹⁴ has been recently developed. The likelihood function used in our earlier paper¹⁴ reflects all observed data and characteristics of the outcome-dependent sampling design. No additional distributional assumptions about the exposure variable are needed, nor is enumeration of the base population required (as is the case with the inverse weighting method).

The present report builds on our previous, technically oriented piece¹⁴ in several ways. First, we provide a more intuitive explanation for why estimator based on outcome-dependent sampling is statistically more efficient than the alternatives. Our simulation study shows this efficiency under a wide range of exposure distributions, and translates the gain in efficiency into the reduction in sample sizes that yield equivalent statistical power. The simulation also explores the impact of various options in sampling, and expresses results in terms of statistical power. The real data examples demonstrate the wide applicability and special advantages of the outcome-dependent sampling design and estimator.

METHODS

A Semiparametric Inference Procedure for Outcome-Dependent Sampling With a Continuous Outcome

In this section, we give a brief overview of the outcome-dependent sampling design,¹⁴ and statistical inference under such a design. Let Y denote the continuous outcome variable and X the exposure. Assume that each Y falls into 1 of 3 mutually exclusive intervals: a lower tail strata, a middle section strata, and an upper tail strata. The general structure of the proposed design consists of 2 components: an overall random sample, and a supplement random sample from each of the 3 strata of Y. Let Ck, k = 1, 2, 3, denote the strata in Y. The data structure in the above design is as follows: one observes the supplement random samples conditional on Y being in strata Ck, ie, {Yki,Xki¦Y∈Ck}, where i = 1, 2, …, nk; one also observes an overall simple random sample whose individuals are denoted by {Y_0i, X_0i}, where i = 1, 2, …, n₀. The total sample size in the outcome-dependent sampling sample is therefore n = n₀ + n₁ + n₂ + n₃, where any of the nk, k = 0, 1, 2, 3 can be zero. The above general sampling strategy encompasses several special cases; for example, when n₁ = n₂ = n₃ = 0, then the outcome-dependent sampling design reduces to the simple random sample or cohort design; alternatively, when Y is binary and n₀ = 0, the outcome-dependent sampling design reduces to the usual case–control design.

Denote by fβ(Y¦X) the conditional density function for the population, where β is the vector of regression coefficients that links exposure X and the outcome Y. Let G and g denote the cumulative distribution and density functions of X, respectively. The joint likelihood function, L(β), of the observed outcome-dependent sampling data is

The component in the first bracket is data contribution to the likelihood from the simple random sample, the second bracket is the contribution from each of the supplement samples. An important feature of this likelihood is easier to appreciate if it is re-expressed. From Bayes formula,

where I is an indicator function for stratum membership, Pr(Yki∈Ck) involves both g and β through Pr(Yki∈Ck)=∫fβ(Yki¦x)g(x)dx. Plugging equation (2) into equation (1), the likelihood function we began with, denoted as L(β, G) now to reflect the dependence on the unknown distribution of X, can be rewritten as

Equation 1

Equation 2

L(β, G) now has 3 components: the specified regression model fβ(Y|X) in the first bracket, the unspecified g(X) in the second bracket, and the outcome-dependent sampling-induced probability P(Yki∈Ck) that ties fβ(Y¦X) and g(X) together in the third bracket. The first component would be the usual likelihood function for observed data, had the sampling been simple random sampling. The last component reflects the biased sampling nature of the outcome-dependent sampling design; ignoring it in the analysis would result in biased estimates of β. Hence, g(X), or G, in the second bracket cannot be simply factored out as would be the case with a simple random sample design. Statistical inference about β, using the standard maximum likelihood estimation method, will depend on a known or a parameterized G. In practice, however, G is rarely known. Misspecification of the distribution could lead to an erroneous conclusion and bias the parameter estimation. Consequently, statistical approaches that do not rely on the extra parameterization of G are desirable.

To estimate β without specifying G(X), Zhou et al¹⁴ developed a maximum likelihood based approach that maximizes L(β, G) by modeling G nonparametrically. This approach used the profile likelihood idea where it fixes β in equation (3) and solves for an empirical likelihood estimate Ĝ(β) from a constrained likelihood function, constraints placed on Ĝ that reflect its properties of being a discrete distribution function, using the Lagrange multiplier technique. An explicit solution for Ĝ(β) can be obtained. Plugging Ĝ(β) into equation (3), the Zhou et al estimator β̂Z can be obtained, using the Newton–Raphson procedure, by maximizing the resulting likelihood. An explicit standard error formula based on an asymptotic distribution has been provided.¹⁴ The statistical program for this analysis can be obtained from the web page (www.bios.unc.edu/∼zhou) or from the authors.

Equation 3

The inverse-probability weighted approach of Horvitz and Thompson¹⁵ can also be adapted in this situation by crudely treating all observed data, including the simple random sample, as if it were sampled from 3 strata, each with a given selection probability. Like the Zhou et al estimator, the inverse probability weighted approach also yields a consistent estimate for β. This method is commonly used with data from a 2-stage study.^16,17 If all N individuals were fully observed in the population, the log likelihood function would be ∑_i=1^N log P(yi¦xi;β). An estimate of this quantity is obtained if we use the completely observed individuals and weight their contributions inversely according to their selection probability into the second stage. The inverse probability weighted estimator β̂IPW is the solution to the following weighted score equation

Equation (Uncited)

where pk can be estimated by

Equation (Uncited)

if there is a complete information. Note that this method needs more information than the likelihood approach employed by Zhou et al method since it requires the sampling probabilities to be known. However, since the inverse-probability weighted approach is based on crudely accounting for the sampling scheme and on an estimating-equations approach, it may not be as efficient as a likelihood-based estimator. When the number of categories of Y is small, the estimating-equation method is not efficient.¹⁸ The realistic settings for the outcome-dependent sampling design we considered had k between 2 to 4. Generally speaking, if the same amount of information is used, among the various statistical approaches to computing any estimator, the maximum likelihood approach is always the most efficient.

A SIMULATION STUDY

We designed a simulation to study the efficiency of different methods under a variety of conditions that one might face in real applications. The basic simulation setting is modeled after a study by Daniels et al¹⁹ of prenatal exposure to low levels of PCB in relation to mental and motor development. An outcome-dependent sampling design was used in the data collection. The data were generated according to the following model,

Equation (Uncited)

where X is the exposure variable that takes on several distributions, p = 1 indicates a linear dose-response relationship and p = 2 represents a nonlinear relationship for E[Y¦X], the Zs are independent covariate variables, and e is a standard normal random error. We generate X from several distributions that include normal, exponential, lognormal, and Bernoulli. These selections reflect possible real situations where X could be a rare-binary variable, a continuous variable, or a skewed variable. We generated Z₁ from a binary distribution (Bernoulli (0.45)), Z₂ from a log-normal distribution (LN (0, 0.25)), and Z₃ from a 3 level polynomial distribution (P(n, 0.3, 0.7)).

In our simulation, we first generated a cohort with 100,000 individuals according to the above model and drew outcome-dependent samples from this cohort. We drew an overall random sample of size n₀, where we observed {Y, X, Z₁, Z₂, Z₃}. We also drew a supplemental random sample of size n₁ from the lower tail of Y defined by {Y<Ȳ−a*σY} and a supplemental random sample of size n₃ from the upper tail of Y defined by {Y>Ȳ+a*σY}, where σ is the standard deviation of Y and a is a known constant. In addition to various configurations for the parameter values, we investigated the effect of varying the value of cut point a on the performance of the methods. We also investigated the impact on statistical power of varying the contribution to the total sample size from the overall random sample and the supplement random samples (ρ = n₀/(n₀ + n₁ + n₃)). Note n₂ = 0 here does not reduce the generality.

Under each setting, we compared the Zhou et al estimator, denoted by β̂Z, with 4 other estimators: (i) the naive maximum likelihood estimator, β̂N, based on the observed ODS data but ignoring the sampling scheme; (ii) an inverse probability weighted estimator (β̂IPW); (iii) the maximum likelihood estimator based on a plain random sample of the same size as the ODS sample (β̂P); and (iv) 2 logistic regression estimators (β̂Lk) based on dichotomizing a continuous Y by defining the outcome D as D = 1 if Y > mean (Y) + kσY and D = 0 otherwise, where k = 0, 1. The weight used in calculating β̂IPW is the inverse of the observed probability of being sampled in the respective strata of Y. The βN and βP estimates are the same as the ordinary least square estimates in our simulation setting. Each set of simulations generated 1000 data sets.

SIMULATION STUDY RESULTS

Table 1 shows the simulation results for a = 1, and (n₀, n₁, n₃) = (200, 100, 100) (hence ρ = 0.5) for various exposure effects. The mean estimate given by β̂N is biased for estimating β₁ in the simulation. Thus, ignoring the sampling scheme (β̂N) leads to a biased estimate for the exposure effect. It is also clear from Table 1 that different dichotomization of a continuous Y will lead to various inconsistent estimates (β̂L0 and β̂L1) of the β₁. Perhaps, more importantly, the logistic estimators can be less able to detect the true underlying relationship, as reflected by corresponding P values of 0.08 for β̂L1, compared with P < 0.05 for all other methods using continuous response. We do not present results for these 2 estimators in future comparisons. The other 3 methods all yielded consistent estimates of β₁. The actual coverage of their nominal 95% confidence intervals (CIs) coverage are all close to 95%, indicating that a good approximation to the asymptotic normality is achieved with this sample size, and the estimated standard errors (SE) are close to the true standard deviations. Under the setting considered, β̂Z has the smallest SE while β̂P has the largest SE. Because β̂N is a biased estimator and its SE underestimates the true variation, we excluded it from the further studies of sample size and power below. The above observations are consistent across different exposure effects listed in Table 1. Under the same linear model but when the X term is quadratic, β̂Z is again more efficient.

Table 1

Results in the lower portion of Table 1 provide the contrast for the normally distributed X to extreme X, namely a skewed exposure (lognormal) and a rare binary exposure (Bernoulli (0.05)). When compared with the Normal distributed exposure, the SE for βZ is even smaller SE in the skewed exposure situations. For the rare binary exposure case, results in Table 1 demonstrate that β̂Z is still the most efficient overall, though the sample size considered was not sufficiently large enough for any of them. This is reflected in the fact that the estimated standard errors are bigger than the estimates of the slope. This is not surprising because with the distribution of a rare binary X, there may not be enough information in the data set as X = 1 could be sparse. Future development of a modified outcome-dependent sampling design for this situation is certainly warranted.

Figure 1 shows the power for testing H₀: β₁ = 0 versus H₁: β₁ = true value, for n = 400 and type I error fixed at 5%. The points corresponding to the true value of β₁ = 0 shows the empirical type I error rate for each test. All 3 methods yield close to 5% type I error. For all 3 estimators, as β₁ increases, so does the statistical power, and the power of β̂Z>β̂IPW>β̂P.

Figure 1

Table 2 shows the sample sizes required to achieve a given statistical power for 3 values of β₁ according to type of estimator under the same settings used in the top panel of Table 1. Use of outcome-dependent sampling with an appropriate estimator requires a smaller sample size. In this setting, the β̂Z method on average needs about 60% of the subjects who would be needed if the study were conducted with a simple random sampling scheme; the β̂IPW method needs about 83%. Further, for a given power, as the true value of β₁ is farther away from 0, relatively fewer subjects are needed to achieve the same power with β̂Z as compared with β̂P. That is, efficiency increased as β₁ is farther away from 0.

Table 2

Figure 2 shows the impact of 2 factors in a given setting of outcome-dependent sampling. The impact of varying a, the cut-point that determined the strata of Y where the supplement random samples were drawn, is shown in Figure 2A. When a = 0, there was only plain random sampling, and in this instance all 3 methods had the same power. As a increases, however, β̂Z (solid line) has better power than the other 2. With β̂Z, the increase in power appears to be monotonic in a. Figure 2B shows the impact of ρ, the fraction of the overall random sample in the total outcome-dependent sample (ρ = n₀/(n₀ + n₁ + n₃)). At ρ = 1, there is no supplemental sample and the outcome-dependent samples are plain random samples. However, as ρ decreases to below 0.7, β̂Z (solid curve) is again the most powerful of the 3 estimators.

Figure 2

Table 3 presents different allocations of outcome-dependent sampling when the exposure variable is skewed (X ∼ lognormal). Compared with the results in the lognormal panel of Table 1, we see that allocating more of the sample to the upper tail of Y improves the efficiency. The standard error of the slope estimator decreases as more of the sample is shifted from the lower tail to the upper tail, ie, 0.0222 → 0.0201 → 0.0192 as the allocation changes (200, 150, 50) → (200, 100, 100) → (200, 50, 150).

Table 3

Real Data Example 1

In the motivating example noted in the introductory text, a nested case–control study would have been implemented if the outcome had been dichotomous. However, the outcome of interest, (the score on the Bayley Scale of Infant Development), was a continuous variable, and, treating it as a dichotomous variable would have resulted in a loss of statistical power. Measurements of polychlorinated biphenyls (the exposure of interest here) are expensive; thus, minimizing sample size while maintaining power was especially important. The authors drew a random sample of cohort members and 2 additional random samples, one from each tail of the outcome distribution.¹⁹ Using the inverse probability-weighted estimator, the estimated β̂IPW was 0.47 BSID units/μg/L PCB with estimated SE as 0.32 (P = 0.14).¹⁴ Using the Zhou et al estimator, the estimated β̂Z was 0.44 with SE = 0.22 (P = 0.02). Using the simple random sample data alone, the β̂N was 0.29 (SE = 0.29, P = 0.32). Daniels et al¹⁹ also examined and confirmed the shape of the dose-response relation in both the outcome-dependent and the simple random samples. The example demonstrates the improved efficiency obtained by the Zhou et al estimator.

Real Data Example 2

Rissanen et al²⁰ examined the relation of serum lycopene concentration to the thickness of carotid arteries among 1028 men in Finland. Using their data, we selected 400 samples in 2 ways. First, we selected a random sample of 400. Second, we selected a random sample of 200, and, plus a random sample of 100 from those with carotid artery thickness above the 90th percentile, and a random sample of 100 from those with carotid artery thickness below the 10th percentile. We then analyzed the data using ordinary regression model with all 1028 samples, and with the 400 selected at random. In addition, we analyzed the 200-100-100 sample using inverse probability-weighted method and then with the Zhou et al estimator. The lycopene results were adjusted for the same covariates (age, year, and sonographer) in all models. Results are given in Table 4. With the full sample of 1028, the estimated coefficient for lycopene was −0.14 with an estimated SE at 0.04 and P = 0.0011. With the random sample of 400, the estimated effect β̂P was −0.10, with SE = 0.06 and P = 0.096. With the 200-100-100 sample, β̂IPW was −0.19 with SE = 0.08, P = 0.017 and β̂Z was −0.24, SE = 0.07 and P = 0.0009. This example suggests that the outcome-dependent sampling scheme and the Zhou et al estimator provided nearly as much power as analysis of the full dataset. Furthermore, the Zhou et al estimator clearly had greater efficiency compared with the inverse-probability weighted approach. The larger βs obtained using outcome-dependent sampling and estimators reflect the shape of the dose-response curve, which had larger negative slopes near the tails of the outcome. While we focused on analyzing the lycopene association as linear (trend-test), a curvilinear approach would better describe the relation and could be easily accommodated with either the inverse-probability weighted or the Zhou et al estimators.

Table 4

Real Data Example 3

Korrick et al²¹ conducted a case–control study of hypertension and bone lead level. For logistic reasons the sampling probabilities for cases and controls could not be determined. Measurements of blood pressure were available. The example, presented as Appendix 1 in the online supplementary material, available with the online version of the article, shows that the proposed outcome-dependent sampling approach could be used to estimate the coefficient for bone lead in a model of blood pressure.

DISCUSSION

When there is a fixed, limited number of observations for examining the linear relation between a continuous exposure and a continuous outcome, outcome-dependent sampling is more efficient than a random sample. The simulations provided an intuitive and simple expression of this benefit: the estimator yields the same estimand as an analysis of the underlying complete population. Stated another way, for a given level of statistical power, the number of observations can be reduced by about 40% compared with a random sample. Thus, the benefits of outcome-dependent sampling apply to continuous outcomes, as well as to dichotomous ones through case–control designs. For a binary response variable, our approach, implemented with logistic regression, is equivalent to the case–control analysis.²²

To implement the weighted method, one needs to know or estimate the weights, this requires at least empirical data about the distribution of Y. Such information may not be a problem for nested studies; however, it can be difficult to calculate the weight for studies that are not nested. The difficultly arises because good quality data on the distribution of Y, and an enumeration of potential subjects, may not be available for the base population. The inverse-probability weighting suffers from the fact that in the typical setting of the outcome-dependent sampling, the natural choice of number of categories of Y is not large enough to yield the variability in weights that would make it efficient.¹⁸ Some recently developed methods may help to identify even more optimal weights.^18,23 The Zhou et al method, on the other hand, does not use selection probability and hence does not need to enumerate the base population.

Our results showed that with a up to about 1 and ρ near 0.5 and a total sample size of 400, outcome-dependent sampling using the estimator of Zhou et al increased power about twice as much as the weighted estimating equation approach (Fig. 2). With larger values of a or smaller values of ρ, the advantage of the Zhou et al estimator was greater. This reflects a greater influence on the regression parameters for data from the tail areas than the middle areas. In general, the efficiency gain of the inverse-probability weighted method over the simple random sample analysis is notable; thus we would suggest in practice that one should at least do the weighted analysis if one cannot implement the Zhou et al method.

If an outcome-dependent sampling procedure is to be used, questions will arise regarding the optimal choice of a and ρ. For example, consider an examination of the serum level of contaminant X among pregnant women in relation to a continuous outcome in offspring. Subject matter considerations might support large values of a (eg, >1) so that a corresponds to a clinically abnormal value. Values of a >1 might seem appealing because of the resulting increase in power, especially with the Zhou et al estimator (Fig. 2). But this reward depends on an assumption about fβ(Y, X) across the range of Y (see the online Appendix 2 for discussion of an example). Similarly, choice of ρ > 0 has several advantages over ρ = 0,¹⁴ since including an overall random sample provides the flexibility of a cohort study and allows for model checking. In general, choosing a ρ from the range of [0.2, 0.5] provides much improved power with the Zhou et al estimator compared with the weighted-estimated-equation estimator, while still allocating enough observations to the simple random sample.

Similarly, the relative size of n₁ and n₃ might be affected by several factors that will vary across studies. For example, if the exposure variable is known to be skewed with a long tail to the right and β is known to be either zero or positive, then increasing the size of n₃ relative to n₁ would be sensible. A large n₃ relative to n₁ would also make sense when there is little interest in the determinants of a low value of the outcome variable.

Prospective designs coupled with relatively expensive measures of exposure are being used with increasing frequency in epidemiologic research. Furthermore, the scope of epidemiologic research increasingly includes outcomes best measured on a continuous scale. Given these trends, methods that allow cost-cutting while maintaining statistical efficiency are likely to see greater use. Recently, similar ideas using outcome-dependent sampling have been extended to the situation in which information other than the exposure variable are also available for both the sample and the rest of the base population.²⁴ Methods have been developed that account for an ordinal outcome variable in a generalized linear model setting; other methods incorporate auxiliary information about the exposure variable that is available for the entire base population.²⁵ Much work, however, remains to be done; for example, how to use outcome-dependent sampling with longitudinal data is still an open question. A survey of the recent statistical research on outcome-dependent sampling design can be found in Zhou and You.²⁶

ACKNOWLEDGMENTS

We thank Clare Weinberg and Beth Gladen for their careful reading of the paper and helpful suggestions. We also thank the reviewers for their helpful suggestions that lead to a much more complete version of the manuscript.

REFERENCES

Supplemental Digital Content

• 00001648-920070700-00001.pdf; [PDF] (25 KB)
• 00001648-920070700-00002.pdf; [PDF] (27 KB)