Statistical Methods to Study Timing of Vulnerability with Sparsely Sampled Data on Environmental Toxicants

1Identifying critical windows of vulnerability to environmental toxicants is an important area in children’s health research. Critical windows of vulnerability are defined as periods during life when an exposure causes a stronger deficit in health later in life compared with other periods when exposure (could have) occurred. The research question is often formulated as “Is the exposure more strongly associated with the health outcome if it occurred in time window 1, 2, …, or n?” Formulating the question in terms of discrete time windows is advantageous from a clinical and practical perspective. However, considering timing of exposure in a continuous fashion may be more advantageous from an analytical standpoint.

2We discuss three statistical approaches that may be useful in studies of windows, or timing, of vulnerability; two of these use time of exposure as a continuous variable. We compare the proposed methods with commonly used statistical approaches such as fitting separate regression models for each potential window and fitting simultaneously adjusted multiple regressions that include all exposure measures (time windows) in one model.

3As exposition, we use data from a study of lead exposure during pregnancy and mental development. We refer readers to published studies of lead exposure for information regarding the toxicologic effects of lead on development and the relevance of time windows of vulnerability in that context. We focus on statistical modeling issues and apply all methods to the same data so that the inferences and interpretations across methods can be compared.

4We use data from an ELEMENT (Early Life Exposures in Mexico to Environmental Toxicants) study cohort of mother–child pairs recruited during pregnancy or before conception. We use maternal blood lead concentrations ascertained during study visits scheduled within each trimester. Low birth weight (< 2,500 g) and preterm (< 37 weeks) children were excluded. Children’s mental development was measured at 24 months using Bayley’s Mental Development Index (MDI24). We also collected other participant characteristics (Table 1); data collection procedures are reported elsewhere. This analysis was restricted to n = 169 participants with complete covariates (mother’s age and IQ, breast-feeding duration, and child’s sex, height, weight, and blood lead level at 24 months) and at least one measure of prenatal lead exposure.

5Participants gave written informed consent before data collection. The study was approved by the institutional review boards of the hospitals where participants were recruited, the National Institute of Public Health of Mexico, Brigham and Women’s Hospital, Harvard School of Public Health, and the University of Michigan.

6For each subject i, assume X1i, X2i, … XKi are measures of exposure taken during Ki time points t1i, t2i, …, tKi; Yi is the health outcome, and Zi are mean-centered covariates. For the lead study, Yi is MDI24, X1i, X2i, … XKi are loge(maternal blood lead levels, micrograms per deciliter), and Zi are maternal age, maternal IQ, child’s sex, and child’s 24-month weight and height-for-age z-score. For methods 1 and 2, we assume the Ki measures of exposure fall within K time windows (i.e., K = three trimesters).

7A commonly used statistical approach in studies of prenatal windows of vulnerability is to fit regression models for each potential window. The magnitude and significance of the regression coefficients are compared to draw conclusions about which window of exposure may be more important.

8This approach is intuitive, but it has two assumptions that are not always discussed. First, the exposure is assumed constant within the window. For example, a single maternal blood lead measure during the first trimester is assumed to be representative of average exposure during the entire trimester. Second, although it is recognized that “the closer we look, the more evident it is that often there is not a uniform response within a given window”, the exposure effect (e.g., the effect of prenatal lead exposure on mental development) is assumed to be constant within a time window. When these assumptions are violated, the estimated effects within a predefined window may be biased toward or away from the null, depending on the true nature of the exposure effect and the time at which exposures are measured. Specific examples are shown in the Supplemental Material, Figures 1 and 2 (doi:10.1289/ehp.1002453).

9The approach is also subject to other statistical limitations. First, estimating separate regressions for each time window precludes formal testing of the differences in effects across time windows. Interpreting nonoverlapping confidence intervals (CIs) as proof that effects differ is not valid, because the estimated coefficients are not statistically independent. This issue is analogous to the incorrect choice of a two-sample t-test over a (correct) paired t-test for testing differences between repeated measures within an individual. Second, it is possible that not all regressions will be based on the same group of subjects because of missing data. For example, a popular approach to maximize the sample size for each regression is to include all observations with data for a given time window in the regression for that time window, even though some may be excluded from models of other time windows because of missing data. However, the available case approach is valid only if data are missing completely at random (MCAR), and it is not possible to verify this assumption based on available (nonmissing) data (for example, not finding significant differences in bivariate analyses comparing subjects with complete versus incomplete data is not sufficient proof of MCAR). Third, because the same outcome is used in all regressions, the hypotheses tested are not independent, and precision may be over- or underestimated because the correlations among the residual errors from the regressions are ignored.

10A straightforward alternative to fitting separate regressions for each time window is to fit a single simultaneously adjusted regression model that includes all time windows. An advantage of this approach is that estimates of the independent effect of exposure during each window may be obtained. However, the approach is not always feasible because of collinearity issues (e.g., inflated standard errors, unstable regression coefficients). Further, models are restricted to observations that have complete data for all variables.

11We estimated separate and simultaneously adjusted multiple linear regression models to estimate trimester-specific associations between MDI24 and loge-transformed maternal blood lead levels.

12Multiple informant data refers to information gathered from different individuals or sources used to measure the same construct. One classic example is where mother and child both respond to questions regarding the physical activity of the child, so that the mother and child are the informants. Methods for multiple informant data also can be applied when exposure information is obtained for the same individual at different time points (windows) by treating the exposure windows as informants. Multiple informant methods can be used to test whether the information relayed by different informants (in this case, whether exposure is measured during different time windows) relates in the same manner to an outcome of interest. Although this approach retains the interpretation (and assumptions) of a set of separate multiple regressions (by providing a single estimate of effect for exposure in each time window), it also provides a way to test differences in associations between the exposure and the outcome across time windows.

13The associations of primary interest are estimated by β1k from window-specific multiple regressions Yi = β0k + β1kXki + β2kZi + ɛki, for window k = 1,2, …, K. The multiple informants approach jointly estimates the regression models.

14Joint estimation enables us to impose and test constraints on the regression coefficients across exposure windows. That is, test whether the exposure coefficients are equal across time windows, Ho: β11 = β12 = … = β1K versus Ha: at least one β1K differs. Parameter estimation can be conducted using generalized estimating equations (GEE) or maximum likelihood (ML) estimation. Details regarding data structure, model fitting, hypothesis testing, and macros to analyze the data using this method are provided in the Supplemental Material (doi:10.1289/ehp.1002453).

15ML estimation proceeds by modeling the exposure at different windows and the outcome as having a K + 1 multivariate normal distribution (i.e., exposure is considered a random variable) and subsequently obtaining the distribution of the outcome conditional on the exposure. The test of equal association between exposure and outcome across time windows can also be performed on standardized regression coefficients (i.e., adjusted correlation coefficient denoted ρ1k = σXkβ1k/σY, k = 1, …, K, where σXk is the SD of the exposure at window k, and σY is the SD of the outcome. Testing adjusted correlation coefficients may be desirable when the variability of the exposure changes over time or if protocols for measuring exposure change across time windows. ML can account for exposure and outcome data missing at random (MAR) and is robust to distributional assumptions when missing data are MCAR. MAR means that the fact the data are missing is independent of the actual missing value, after accounting for other observed participant characteristics. MCAR means that the fact the data are missing is independent of the missing value as well as other observed data and hence is more restrictive than MAR.

16The GEE approach embeds separate linear regression models for each time window into a unified set of estimating equations. In contrast to the ML approach, the exposures are not considered random (dependent) variables. However, unless more sophisticated methods are used, the GEE method retains the MCAR assumption for the missing data.

17We applied the multiple informant approach to the ELEMENT study data. For the GEE approach we used a score type test of the hypothesis Ho: β11 = β12 = β13 [degrees of freedeom (df) = 2]. Because the variance of the exposure during trimester 2 was smaller than during the other two trimesters, in the ML approach we tested for homogeneity of standardized exposure effects across time windows (Ho: ρ11 = ρ12 = ρ13) using a likelihood ratio test (df = 2). The p-values for these tests are denoted pint, alluding to an interaction between exposure level and timing of exposure.

18When assumptions regarding the timing of exposure in methods 1 and 2 are violated (i.e., when timing of exposure measurement or the effect of exposure on the outcome varies within time windows), when a large number of time windows are considered, or when the number of measurement occasions for the exposure varies across participants, the multiple informant approach may no longer be feasible. An alternative is to model the pattern of exposure for each individual over time and then relate exposure features to the outcome.

19The underlying idea in this approach is to reduce the number of exposure measures for each individual from Ki, which may vary across participants, to a smaller number equal across participants. A simple example where the exposures are summarized into two exposure features is to model the exposures Xik using a random intercepts and random slopes model, Xik = θ0i + θ1i tik + ɛik, where k = 1, …, Ki represents the kth measurement occasion. The random effects θ0i, and θ1i are person-specific intercepts and slopes that jointly describe the pattern of exposure for individual i and have a population mean θ = (θ0,θ1) and variance Φ. Next, the outcome is modeled in relation to the exposure features θ0i, and θ1i, for example, Yi = β0 + β11 θ0i + β12 θ1i + β2 Zi + ɛi. [Details regarding data structure, model fitting, hypothesis testing, and macros to analyze the data using this method are provided in the Supplemental Material (doi:10.1289/ehp.1002453).]

20For the ELEMENT study, where a few measures of exposure were available, we employed a model where a random intercept and slope were the only subject-specific exposure features. The time variable was centered at 7 weeks, the middle of the first trimester, such that θ0i represents blood lead levels during gestational week 7. First, the exposure was modeled as Xik = θ0i + θ1i (tik −7) + ɛik, with tik being gestational time in weeks. Because week 7 is the midpoint of the first trimester, θ0i can also be interpreted as the average exposure for individual i during the first trimester, so that study participants with a random intercept θ0i higher than the mean θ0 would have a higher than average lead exposure in early pregnancy. The slope θ1i represents the subject’s rate of change in exposure across the pregnancy. For example, if the population slope θ1 was negative (indicating declining lead levels across pregnancy), a study participant with a random slope θ1i larger than the mean θ1 would have a slower than average rate of decline in lead levels over the course of pregnancy. We modeled the outcome (MDI24) as Yi = β0 + β11 θ0i + β12 θ1i + β2 Zi+ ɛi, where β11 is interpreted as the association between exposure during the first trimester and MDI24, and β12 is the association between changes in exposure over the course of pregnancy and MDI24.

21We fit the model for the exposures and the outcome model jointly using ML, which yields consistent and efficient estimates of model parameters. Heuristically, a two-step approach also could be used where empirical Bayes estimates of the exposure features are obtained, for example, θ̂0i, and θ̂1i, and then substituted in the outcome model. However, the two-step approach gives inconsistent (likely attenuated) estimates of the regression coefficients.

22Method 3 does not require that participants be measured at the same time points, whereas methods 1 and 2 assume the timing of measurements is similar across participants. Furthermore, because the timing of exposure is not discretized into windows, this approach is more useful than the multiple informant approach when windows are not well defined or when little prior information about the etiologically relevant windows are is available. It requires that each participant has sufficient information such that the features can be reliably estimated (i.e., at least one more measurement than features in most participants). In the lead example, at least three measurements are required for most participants, because two features are estimated. Although the interpretation of model parameters may not be as straightforward because the reference to discrete windows is lost, conclusions about the relative importance of broad time periods can still be drawn. In the lead example, this method helps answer the question “After accounting for first-trimester exposure, does changing the exposure in subsequent trimesters matter?” In the lead example, this approach assumes a constant rate of decline in exposure, and the association of exposure with outcome over the course of pregnancy diminishes (or increases) at a linear rate.

23The fourth method consists of describing the population-average exposure pattern for levels of the outcome. For example, in ELEMENT, this method can be used to compare the pattern of exposure for children who achieve high MDI24 scores with the pattern in children with low scores. This approach differs from methods 1–3 because it models the exposure given the health outcome.

24The exposure is modeled as Xik = f0(tik) + f1(tik)Yci + δik, where Yci is the outcome for subject i centered at the sample mean (Yci = Yi − Ȳ) or centered at its predicted value (Yci = Yi − Ŷi, where Ŷi is the predicted outcome given Zi) given factors other than exposure (e.g., suspected confounders and possibly other independent predictors of the outcome). The residuals δik are assumed to have mean zero and covariance Δ within individuals but are independent across individuals. The term f0(tik) represents the exposure pattern over time for those with an average outcome, which can be modeled as a parametric function [e.g., f0(tik) = α00 + α01tik + α02t2ik] or a semiparametric curve (e.g., penalized splines). The term f0(tik)Yci quantifies the differences in exposure over time across levels of the outcome. Both f1(tik) and f0(tik) are modeled in the same fashion [e.g., f1(tik) = α10 + α11tik + α12t2ik if f0(tik) is modeled as a quadratic function]. The coefficients α10, α11, α12 jointly describe the pattern of exposure curve over time.

25An overall test of association between exposure and outcome involves testing Ho: f1(t) = 0 vs. Ha: f1(t) ≠ 0. When f1(t) is given by a quadratic function, then the null hypothesis is Ho: α10 = α11 = α12 = 0 versus Ha: at least one α differs from zero. Testing whether the exposure pattern over time is associated with outcome amounts to testing Ho: α11 = α12 = 0 versus Ha: at least one α differs from zero. For timing of susceptibility (i.e., whether exposure effects vary depending upon the timing of exposure), the more relevant hypothesis is the latter.

26The estimation of this model can be conducted in two steps. First, a model for the outcome is estimated where only covariates Zi are predictors. The residuals Yci = Yi − Ŷi are then constructed, and the model for Xik (i.e., a model of the exposure conditional on the outcome) is estimated. Details regarding data structure, model fitting, hypothesis testing, and macros to analyze the data using method 4 are provided in the Supplemental Material (doi:10.1289/ehp.1002453).

27In the ELEMENT data, we obtained Yci = Yi − Ŷi, and interpreted it as the deviation of the individual from the expected MDI24 score given mother’s age, mother’s IQ, duration of breastfeeding, sex, and weight and height z-score at 24 months. For example, Yci will be positive when the MDI24 score of a child is higher than predicted given the characteristics listed above. We then estimated the exposure model using tik as gestational time in weeks and nonparametric and quadratic models for f0(t) and f1(t). Finally, we constructed the predicted exposure pattern (back-transformed to natural units of blood lead) for children in the 10th and 90th percentile of the covariate-adjusted MDI distribution and determined the relative difference in exposure between the 10th and 90th percentiles of the covariate-adjusted MDI distribution across time by exponentiating the difference in the predictions of loge-transformed blood lead.

28The study visits occurred, on average, toward the end of each trimester, with considerable variability in their timing (Table 1). Notably, although the earliest visit among all study participants was at 3.8 weeks of gestation, the average first-trimester visit occurred at 13.7 weeks. The variability in timing of measurements raises concerns about interpretation of regression coefficients as the association between exposure during each trimester and mental development and was the primary motivation for seeking alternative analytical methods for this type of data. Nevertheless, we keep the language of trimesters in the subsequent descriptions of the results for methods 1 and 2.

29The estimated associations from the separate and simultaneously adjusted regression models were as follows: −5.42 (95% CI, −10.2 to −0.64) and −2.74 (95% CI, −5.78 to 0.29) MDI24 points per 1-loge-unit increase in blood lead in the first trimester, respectively (Table 2). Although not significant, the coefficients for trimesters 2 and 3 from the simultaneously adjusted regression were positive (suggesting higher MDI24 scores with higher lead exposure) and imprecise, which would be consistent with collinearity. The coefficients from separate regressions were all negative and decreased in strength across trimesters, with the strongest association observed at trimester 1. However, none of the estimated effects were statistically significant (p > 0.05).

30The estimated associations between first-trimester exposure and MDI24 from the multiple informant approach were −2.74 (95% CI, −5.82 to 0.33) and −4.13 (95% CI, −7.54 to −0.72) MDI24 points per 1-loge-unit increase in blood lead from GEE and ML estimation, respectively (Table 2). As is always the case, the GEE estimates were equal to the point estimates obtained from separate regressions, but the CIs varied slightly because the GEE method takes into account within-individual correlation across the time windows. Because there were missing data, the ML and GEE approaches do not give the same estimates. The ML approach estimated a stronger association for trimester 1 (e.g., βGEE = −2.74 vs. βML = −4.13). Although a trend of increasing effect of lead in earlier times during pregnancy is suggested, tests of a varying exposure effect were not significant (GEE pint = 0.56, score test of equal regression coefficients; MLE pint = 0.23, likelihood ratio test of equal standardized coefficients). The ML approach detected a significant association at trimester 1 and exhibited a smaller p-value than GEE for the test of differences in association across trimesters. This reflects greater efficiency of ML estimation relative to GEE when data are missing and when differences in standardized regression coefficients are tested, in addition to using score versus likelihood ratio tests.

31The exposure model parameters from method 3 indicate that the average exposure at 7 weeks of gestation was 1.90 loge(blood lead) units (approximately = 6.69 μg/dL), SD 0.49 (Table 3). This is consistent with the average and SDs of the observed measures during trimester 1 (Table 1). There was a significant average decline across pregnancy (average linear decline is 0.04 ln(blood lead) units per 12 weeks) (Table 3).

32The regression coefficients for blood lead level at 7 weeks and the average change in blood lead level are expressed in MDI points associated with a 1-SD change in these predictors. A 1-SD increase in loge-transformed blood lead levels at 7 weeks was significantly associated with a 2.11-point decrease in MDI24 (95% CI, −0.01 to 4.23). Although the changes in exposure over the course of pregnancy were not significantly associated with MDI24 (p = 0.73), the positive association, indicating that a slower rate of decline in blood lead (i.e., sustained exposure) was associated with a 0.58 point increase in MDI24, was unexpected. However, this relation may have resulted from collinearity; although the random intercept and slope were less correlated than the raw lead concentrations (e.g., correlation between trimester 1 and 2 measurements is 0.66), they were still significantly and negatively correlated (correlation is −0.56).

33Notably, the correlation between the predicted random intercept and the observed blood lead level at trimester 1 was 0.97, indicating that the random intercept was closely related to the measured exposure at the first visit. The advantage of the random intercept is that it is interpreted as exposure level at 7 weeks for all study participants and can be estimated for participants who did not have a measured blood lead level at trimester 1. This advantage translates to stronger standardized estimates for the association between exposure and MDI [−2.11/0.49 = −4.31 points (95% CI, −8.66 to 0.02) for a 1-log-unit increase in exposure] compared with corresponding estimates based on methods 1 and 2 (2.74 points via separate regressions or GEE, −4.13 points via MLE).

34The exposure patterns from method 4 were estimated with quadratic functions (the shapes of f0(t) and f1(t) did not deviate significantly from quadratic); the coefficients for these functions are given in Table 4. Lead exposure and MDI24 were significantly associated (Ho: α10 = α11 = α12 = 0, p = 0.03) based on this model. The exposure patterns over time comparing high with low achievers marginally differed in their shape over time (Ho: α11 = α12 = 0, p = 0.086). Figure 1 portrays the prenatal exposure patterns for children in the 90th and 10th percentile of the covariate-adjusted MDI distribution (top) and their relative difference (bottom). The earliest time at which we estimated exposures was at 4 weeks, the first observed time point among all mothers. As can be observed, the exposure patterns (i.e., the shapes of the curves over time) were similar toward the end of pregnancy but had less overlap earlier in pregnancy. More specifically, children in the 10th percentile of the distribution had significantly higher exposures during the first 17 weeks of pregnancy compared with children in the lower 90th percentile of the distribution. At 7 weeks of pregnancy, those in the 10th percentile had exposure 1.63 (95% CI, 1.09 to 2.45) times higher compared with those in the 90th percentile. The significant relative exposure difference persisted until approximately 17 weeks (relative difference = 1.24; 95% CI, 1.00–1.53).

35Identifying critical windows of vulnerability is an emerging field in children’s health research. Existing statistical approaches used in this area, primarily based on multiple regression, have important limitations: Formal tests for the difference across windows cannot be performed, missing data cannot be easily incorporated, and variation in the timing of exposure cannot be easily accommodated. We presented three alternative approaches that mitigate these limitations. Table 5 summarizes their assumptions, and compares them with multiple regression. To our knowledge, this is the first application of a multiple informants approach (method 2) and population exposure pattern approach (method 4) to studies of timing of vulnerability in environmental epidemiology.

36Modifications of multiple regression have been discussed previously in the literature (e.g.,;;, ). These methods involve a) including as predictors the residuals from models of the exposure during each window regressed on a reference window; b) using ratios of exposure during each window to a reference window; and c) including exposures outside the windows of interests (e.g., postnatal exposure) as a control variable. However, these methods are limited primarily by their lack of ability to incorporate missing data.

37Methods beyond multiple regression that address windows of vulnerability have been discussed in the literature under the rubric of exposure–time response models. These models focus on estimating a weight function, w(t), that measures the relative effect of an exposure increment at time t compared with other times. The time metric t can be time at exposure as in our case and other examples in cancer research or time since exposure, as in studies quantifying latency. These methods, and the ones we propose, are geared toward settings of protracted exposures. Specifically, our method 3 can be seen as a special case of exposure–time response models, because the coefficients corresponding to features θ0i (level at week 7) and θ1i (rate of change through pregnancy) can be transformed to obtain an estimate of w(t). The estimated w(t) in our case would be a linear function [see Supplemental Material, Figure 4 (doi:10.1289/ehp.1002453)]. Furthermore, related approaches have also been discussed in life course epidemiology, where one of the approaches is to develop a taxonomy of exposure patterns based on dichotomized exposures at each exposure window and relate the patterns to health outcomes. All of these methods are related, but a distinguishing factor is the time scale on which they focus and whether time of exposure is discretized or treated as continuous.

38Although the proposed approaches circumvent some drawbacks of fitting multiple regression models, they are not without limitations. For example, the tests offered by method 2 are low-powered. In principle, method 3 allowed us to flexibly model individual exposure patterns over time; however, relating outcome to many exposure features may become cumbersome to interpret (e.g., individual-level quadratic trends). Method 4 could be improved by incorporating confidence bands instead of pointwise CIs when presenting the results graphically. However, the test of whether population exposure patterns vary over time between outcome groups is valid.

39When the windows of vulnerability refer to prenatal windows, inaccurate ascertainment of gestational age may limit the applicability/interpretation of the above methods. Inaccurate measurement of gestational time may arise from incorrect reporting of last menstrual period, for example. This introduces measurement error and may attenuate effect estimates. When gestational time (e.g., weeks) is used to define windows of vulnerability (e.g., trimesters), then even a few days’ error in gestational age could lead to misclassifying the exposure as occurring in one trimester versus another. Furthermore, because of field issues (e.g., rescheduled appointments), some participants may be seen during a later trimester but their visit recorded as occurring within the previous time window. Misclassification of the timing of exposure would lead to biased exposure–outcome associations. Method 2 (like method 1) is not robust to misclassification of the timing of exposure. Methods 3 and 4 may be less affected by this issue, however, because the time is treated as continuous, such that the extent of misclassification is of a lesser magnitude. Other issues of how inaccurate ascertainment of gestational age may introduce bias in exposure–birth outcome studies, even when using ultrasound to ascertain gestational age, have been discussed in the literature.

40We described several approaches to study timing of vulnerability using an example with few exposure windows and few covariates. Extensions of the approach to incorporate more covariate information on exposure models for methods 3 and 4 or using categorical instead of continuous outcomes are discussed in the Supplemental Material (doi:10.1289/ehp.1002453). Sensitivity analyses are also included in the Supplemental Material, where blood lead measured at 24 months is evaluated both as an extra window using method 2 and as a confounder using all methods.

41Uncovering windows of vulnerability to environmental pollutants is a complex question that requires sophisticated data analysis tools. Methods that test exposure effect differences across time of exposure should be employed before concluding which window is most important.

42Of the methods presented, method 2 is preferred over method 1, because it enables formal testing of differences in effects across a priori–defined windows. Methods 3 and 4 are preferred over method 2 when there is large variability in the timing of exposure across participants, as was the case in the example presented.