Time-varying associations between an exposure history and a subsequent health outcome: a landmark approach to identify critical windows

Abstract

Background

Long-term behavioral and health risk factors constitute a primary focus of research on the etiology of chronic diseases. Yet, identifying critical time-windows during which risk factors have the strongest impact on disease risk is challenging. To assess the trajectory of association of an exposure history with an outcome, the weighted cumulative exposure index (WCIE) has been proposed, with weights reflecting the relative importance of exposures at different times. However, WCIE is restricted to a complete observed error-free exposure whereas exposures are often measured with intermittent missingness and error. Moreover, it rarely explores exposure history that is very distant from the outcome as usually sought in life-course epidemiology.

Methods

We extend the WCIE methodology to (i) exposures that are intermittently measured with error, and (ii) contexts where the exposure time-window precedes the outcome time-window using a landmark approach. First, the individual exposure history up to the landmark time is estimated using a mixed model that handles missing data and error in exposure measurement, and the predicted complete error-free exposure history is derived. Then the WCIE methodology is applied to assess the trajectory of association between the predicted exposure history and the health outcome collected after the landmark time. In our context, the health outcome is a longitudinal marker analyzed using a mixed model.

Results

A simulation study first demonstrates the correct inference obtained with this approach. Then, applied to the Nurses’ Health Study (19,415 women) to investigate the association between body mass index history (collected from midlife) and subsequent cognitive decline (evaluated after age 70), the method identified two major critical windows of association: long before the first cognitive evaluation (roughly 24 to 12 years), higher levels of BMI were associated with poorer cognition. In contrast, adjusted for the whole history, higher levels of BMI became associated with better cognition in the last years prior to the first cognitive interview, thus reflecting reverse causation (changes in exposure due to underlying disease).

Conclusions

This approach, easy to implement, provides a flexible tool for studying complex dynamic relationships and identifying critical time windows while accounting for exposure measurement errors.

Background

Long-term lifestyle, environmental or occupational exposures may have a major impact on the occurrence of chronic diseases. Yet, identifying critical time-windows during which risk factors have the strongest impact on disease risk is challenging. Chronic diseases often result from a long and accumulating pathological process that evolves over years before diagnosis [1]. In such context, the exposure time-windows close to the clinical event may be less meaningful in terms of etiology and, importantly, may be obscured by reverse causality (which occurs when behaviors or exposures change as the disease progresses in infra-clinic stages). For example, while obesity in midlife is a risk factor for subsequent dementia, it is consistently reported as a protective factor later in life; indeed advanced neuropathology, by altering olfactory function and taste, potentially leads to malnutrition and weight loss in late-life [2–5]. Thus, the only valuable approach to evaluate causal associations linking cumulative unhealthy body weight to cognitive aging might be to estimate the relationship between the entire history of exposure that precedes and begins well upstream of the period at risk of the event.

Exposure history metrics have been developed to estimate the cumulative effect of time-varying exposure on disease endpoints [6, 7]. However, to be relevant in life-course epidemiology, such methodologies absolutely need to handle: (i) associations that can vary according to the age during exposure or the distance between the exposure and the disease endpoint, (ii) exposures that are measured only at sparse time points and with error, and (iii) exposure and disease endpoint time-windows that do not necessarily coincide. Finally, they should also apply to longitudinal disease outcomes (not only survival or binary) and their change over time such as cognitive decline that is a strong predictor of dementia. This work specifically aimed at extending the methodology of exposure history metrics to address all these important issues.

The most common exposure history metric is the cumulative index of exposure (CIE) [6, 7]. Computed as the un-weighted sum of all past exposures, CIE assumes that past values of exposure are of equivalent importance (see Fig. 1, Scenario A) which may induce etiologically inadequate conclusions when the effect of the exposure on health outcomes likely varies according to the age or timing of exposures (see Fig. 1, Scenarios B and C). To address this challenge, Breslow et al. [8] and Thomas [9] introduced the concept of weighted CIE (WCIE), that combines information about duration, intensity and timing of the exposure. WCIE represents the time-weighted sum of past exposures, with weights reflecting the relative importance of exposures at different times. Therefore, unlike CIE, the estimated effects of the exposure history on health status are time-varying, a key aspect for the identification of critical windows of exposure in life-course epidemiology. Weights may be assigned a priori by using a known parametric weight function in the presence of biologically/clinically relevant assumptions [10–12] or when assessing very specific assumptions such as different life course hypotheses (e.g., progressive accumulation, critical period, recent period) [13, 14]. Nevertheless, in most cases, the form of the weight function is not known and needs to be directly determined from the data. WCIE models are found in a wide range of applications with binary and time-to-event outcomes, such as concentration of occupational asbestos in relation to mesothelioma [15, 16] or drug use linked to fall-related injuries [17]. In environmental health research, complex associations have been explored through distributed lag models (DLM) [18], which are very similar to WCIE models in the context where long-term exposure and outcome are collected at identical times across individuals.

Fig. 1

Three examples of trajectories of association between an exposure history from time −S to time 0 and a subsequent health outcome: (A) constant association, (B) time-varying with remote association or (C) time-varying with opposing remote and recent associations. The dashed horizontal lines represent the 0 association

A major limit of CIE, WCIE and even DLM methodologies is that they always require the exposure history to be complete. Yet, intermittent missing data on individual exposures are frequently reported in case-control and cohort studies, in particular when follow-up is long, thus preventing any analysis using these metrics. Second, although measurement error is inevitable, CIE and WCIE metrics also require the exposures to be measured without error and the violation of this assumption (which is likely in most epidemiological contexts) may lead to biased exposure-risk relationships [19]. To overcome these limitations, Mauff et al. [20] introduced a weighted cumulative association structure within a joint model to assess the impact of an endogenous time-varying exposure measured intermittently and with error on the risk of a time-to-event outcome; the exposure model accounted for independent centered Gaussian measurement errors through a standard linear mixed model. However, in this approach, exposure and outcome have to be evaluated during the same time-window, while as introduced previously, in some contexts the interest is in the cumulative effect of an exposure history that precedes evaluation of disease risk. In environmental health research, Chen et al. [21] proposed a reverse DLM, which models the exposure trajectory as a function of the health outcome rather than the inverse. By doing so, the reverse DLM naturally accounts for measurement errors and missing values in exposures. However, as a counterpart, this model conditions on the health outcome and thus considers it is measured without error. Finally, all these methodologies have been developed for survival or cross-sectional continuous or binary outcomes when we are interested in a longitudinal health outcome. To date, we are aware of only one pharmaco-epidemiological study applying the WCIE in association with a longitudinal outcome [22]. Nonetheless, it considered a complete error-free exposure history measured concomitantly with the outcome, and it only assessed its association with the level of the longitudinal outcome while in etiological studies, the interest lies in the change over time of the longitudinal outcome.

In this manuscript, we propose to extend the WCIE methodology to address life-course epidemiology research questions. We rely on a landmark approach which consists in limiting the analysis of the outcome of interest after a certain time of interest, called landmark, and considering as predictors the exposure history up to the landmark. Although mainly developed in dynamic prediction context [23, 24], it also addresses the issue of non-concurrency between the exposure and the outcome time-windows central in life-long epidemiology.

In “Methods” section, we first introduce the proposed methodology, and describe the estimation procedure. In “Application to body mass index history starting at mid-life and subsequent cognitive trajectories after age 70” section, we apply the method in the Nurses’ Health Study to investigate the association between BMI history starting in midlife and repeated cognition assessed after age 70. In “Simulation study” section, we demonstrate that the estimation procedure provides correct inference in a simulation study. Finally we discuss the results and methodology, and conclude.

Methods

A landmark approach for assessing the dynamics of association between an intermittently evaluated time-varying exposure and a subsequent health outcome

General study framework

Within a prospective cohort study, we consider a sample of N individuals present in the cohort at a landmark time of interest from which the health outcome is to be studied, hereafter called time 0 (see Fig. 2). The temporal framework of the study consists in repeated measures of the exposure prior to the landmark time, that is from a negative time −S (e.g., time of entry in the cohort) up to time 0, as well as repeated measures of the outcome collected at and after the landmark time. Following the landmark framework [23], exposures data collected after time 0, if any, are not included. Both exposure and health outcome are collected at discrete and individual-specific occasions with error. The design of the Nurses’ Health Study used in the application (“Application to body mass index history starting at mid-life and subsequent cognitive trajectories after age 70” section) naturally follows the landmark framework since the exposure was assessed several decades before the first assessment of outcome. In this manuscript, we will focus mainly on a longitudinal heath outcome and an exposure history constructed based on the time preceding the first outcome assessment. Nevertheless, the proposed approach also applies to other types of event, such as time-to-event and binary endpoints, and different timescales sur as age for the exposure and landmark age for the outcome (e.g., age at menopause).

Fig. 2

Temporal representation of the non-concomitant measurement times of the exposure and the subsequent health outcome considered in our study framework. For both exposure and health outcome, measurements are collected at discrete and individual-specific occasions with error. The exposure history is constructed according to the time t_U prior to the first outcome assessment at time 0 (−S≤t_U≤0). The longitudinal health outcome Y is modeled prospectively according to time t (t≥0)

Mixed model for the exposure

Let U_il denote the exposure values collected for individual i (i = 1,..., N) at the retrospective measurement time t_Uil (l=1,...,m_i) before the landmark time so that t_Uil∈[−S,0] (see Fig. 2). The times and the total number m_i of repeated measures can differ from one individual to the other thus implicitly handling intermittent missing values. Although the methodology could apply beyond, we consider here a continuous exposure and describe its trajectory via a standard linear mixed effects model [25]:

$$ \begin{aligned} U_{{il}} &= U^{*}_{i}(t_{{Uil}}) + \varepsilon_{{il}}\\ &= X_{i}(t_{{Uil}})^{\top}\beta + Z_{i}(t_{{Uil}})^{\top} b_{i} + \varepsilon_{{il}} \end{aligned} $$

(1)

where $U^{*}_{i}(t_{{Uil}})$ is the true underlying exposure value at measurement time t_Uil,X_i(t_Uil) and Z_i(t_Uil) are two vectors of covariates at time t_Uil associated with the vectors of fixed effects β and of random effects b_i, respectively; $b_{i} \sim \mathcal {N}(0,B)$. We consider random Gaussian measurement errors ε_il with mean 0 and variance $\sigma ^{2}_{\varepsilon }$; b_i and ε_il are independent. The objective of this linear mixed model is to accurately model the underlying true exposure $U^{*}_{i}(t_{U})$ at any time t_U in [−S,0] for inclusion in the WCIE. As such, both X_i(t_U) and Z_i(t_U) should include flexible functions of time, for instance a basis of splines [26] or fractional polynomials [27].

Weighted cumulative index of exposure

WCIE is defined as the weighted sum of the underlying true exposures during the period of interest. Without loss of generality, we consider the entire window [−S,0] although any sub-window could be considered instead:

$$ {WCIE}_{i} = \sum \limits_{t_{U}=-S}^{0} w(t_{U})U^{*}_{i}(t_{U}) $$

(2)

where w(t_U) is the weight function assigned to the history of the true exposure $U_{i}^{*}(t_{U})$ during the S+1 time units (e.g., years) preceding the initial assessment of the outcome (see Fig. 1). Note that for ease of epidemiological interpretation, we chose to define WCIE as the sum of exposure levels at discrete times (e.g., every year) but the methodology also applies to WCIE defined as the integral of the exposure levels in [−S,0].

Model for the health outcome

Although transposable to cross-sectional or time-to-event outcomes, our primary interest was to estimate the time-varying effects of a past exposure history on a subsequent longitudinal health outcome (see Fig. 2). Thus, let’s denote Y_ij the measure for individual i (i=1,...,N) collected at time t_ij≥0 with j=1,...,n_i. The times and the total number n_i of repeated measures can differ from one individual to the other. Change over time of Y is modeled using a linear mixed model [25]. For ease of presentation, we consider here a linear trajectory over time, and thus introduce two time-varying effects of the exposure, one on the level at the landmark time (denoted WCIE_I) and one on the change over time of Y (denoted WCIE_S):

$$ {}\begin{aligned} Y_{{ij}} =~&\alpha_{0}~+~\tilde{X_{i}}(t_{{ij}})^{\top} \alpha_{1}~+~ {WCIE}_{{Ii}}~\gamma_{I}~+~c_{0i}\\ &\!\! +\left(\!\alpha_{2}\,+\,\tilde{X_{i}}(t_{{ij}})\!^{\top} \alpha_{3}\,+\, {WCIE}_{{Si}}~\gamma_{S}~\,+\,~c_{1i} \right)\!\times t_{{ij}}~\,+\,~\tilde{\varepsilon}_{{ij}}\\ =~\!& \alpha_{0}~\!\,+\,\!~\tilde{X_{i}}(t_{{ij}})\!^{\top} \alpha_{1}~\!\,+\,\!~ \sum \limits_{t_{U}=-S}^{0} \gamma_{I}~w_{I}(t_{U})U^{*}_{i}(t_{U})~\!\,+\,~c_{0i}\\ &\!\,+\,\left(\!\!\alpha_{2}\!\,+\,\tilde{X_{i}}(\!t_{{ij}}\!)\!^{\top} \alpha_{3}\,+\,\! \sum \limits_{t_{U}=-S}^{0} \gamma_{S}~w_{S}(\!t_{U}\!)U^{*}_{i}(t_{U}\!)\!\,+\,c_{1i} \!\!\right)\\ &\times t_{{ij}} + \tilde{\varepsilon}_{{ij}}\\ \end{aligned} $$

(3)

where $\tilde {X_{i}}(t_{{ij}})$ is a vector of covariates at time t_ij associated with the vectors of fixed effects α₁ and α₃; WCIE_Ii and WCIE_Si are the weighted cumulative exposure covariates associated to the initial level and to the slope of Y, respectively. They are defined according to Eq. 2 with different weights w_I(t_U) and w_S(t_U), and are associated with fixed effects γ_I and γ_S; c_0i and c_1i are correlated individual random intercept and slope, respectively, with $c_{i} \sim \mathcal {N}(0,\tilde {B})$; $\tilde {\varepsilon }_{{ij}}$ are the independent Gaussian measurement errors with mean zero and variance $\sigma ^{2}_{\tilde {\varepsilon }}$.

Identification of the trajectories of association

The linear mixed model of the health outcome defined in Eq. 3 provides two trajectories of association (i.e., time-varying effects) of the history of past exposure on the subsequent health outcome, one for the initial level (noted $\gamma _{I}^{*}(t_{U})$), one for the slope (noted $\gamma _{S}^{*}(t_{U})$):

$$ \begin{aligned} \gamma_{I}^{*}(t_{U}) &= \gamma_{I} w_{I}(t_{U})\\ \gamma_{S}^{*}(t_{U}) &= \gamma_{S} w_{S}(t_{U}) \end{aligned} $$

(4)

They represent the mean difference of initial level (i.e. $\gamma ^{*}_{I}(t_{U})$) or of change per unit of time (i.e. $\gamma ^{*}_{S}(t_{U})$) when the exposure increases of 1-unit at time t_U, adjusted for other covariates and when the exposure history is stricly similar at any other time in [−S,0].

Specification of weights

In some applications, the form of the weight function w_.(t_U) (where subscript. refers either to I or S in Eqs. 3 and (4)) is known. However, most often, it has to be estimated directly from the data. As others previously [22, 28, 29], we chose to estimate the weight function by regression using a basis of splines [26], which are flexible enough to capture a variety of clinically plausible shapes [30]. Thus, the weight function can be written:

$$ w_.(t_{U}) = \theta_{.0} + \sum \limits_{k=1}^{K} \theta_{.k} B_{k}(t_{U}) = \sum \limits_{k=0}^{K} \theta_{.k} B_{k}(t_{U}) $$

(5)

where (B_k)_k=1,...K refers to the K basis of splines functions and (θ_.k)_k=1,...K the coefficients to estimate. For ease of calculation, we denote the intercept θ_.0B₀(t_U) with B₀(t_U)=1,∀ t_U.

Although any type of splines could be considered, we favored natural cubic splines that limit erratic behavior at the boundary of the time window thanks to linearity constraints [31]. In that case, the basis of splines is built from K+1 knots which are to be chosen inside the time window [−S,0]. The number of knots has to be carefully determined from the data. A large number of knots implies high flexibility but it may lead to overfitting. On the contrary, a small number of knots may result in an oversmooth estimate that is prone to under-fit bias [31]. In our work, we considered between one and five inner knots (i.e., K∈[2,6]) and relied on the Akaike information criterion [32] to select the final splines basis. In addition, in the absence of prior knowledge, we considered equidistant knots.

One additional advantage of approximating the weight function using splines is that the WCIE can be rewritten as a linear combination of K+1 components:

$$ \begin{aligned} {WCIE}_{.i} &= \sum \limits_{t_{U}=-S}^{0} w.(t_{U})U^{*}_{i}(t_{U})\\ &=\sum \limits_{t_{U}=-S}^{0} \sum \limits_{k=0}^{K} \theta_{.k} B_{k}(t_{U}) U^{*}_{i}(t_{U}) \\ &= \sum \limits_{k=0}^{K} \theta_{.k} \underbrace{\sum \limits_{t_{U}=-S}^{0} B_{k}(t_{U}) U^{*}_{i}(t_{U})} \\ &= \sum \limits_{k=0}^{K} \theta_{.k} \hspace*{1.2cm} H_{{ki}} \end{aligned} $$

(6)

where H_ki (for k=0, …,K) denote intermediate summary covariates of the exposure history [22, 29, 33].

Identifiability constraints

The time-varying effects of the exposure history defined in Eq. 4 is overparameterized; coefficients γ_I and γ_S, and the weights w_I(t_U) and w_S(t_U) cannot be simultaneously estimated from the data without further constraint. Hauptmann et al. [28] proposed in a case-control study to estimate the total effect of the history of exposures γ_. and constrain w_.(t_U) with the constraint $\phantom {\dot {i}\!}\sum \limits _{t_{U}=-S}^{0} w_.(t_{U}) = S+1$. Following Sylvestre et al. [29] and Danieli et al. [22], we directly estimated the time-varying coefficients without constraining the weights by setting the coefficients γ_I and γ_S to 1. Therefore, the time-varying effects directly correspond to the weights and Eq. 4 becomes:

$$ \begin{aligned} \gamma_{I}^{*}(t_{U}) &= w_{I}(t_{U}) \\ \gamma_{S}^{*}(t_{U}) &= w_{S}(t_{U}) \end{aligned} $$

(7)

With this parameterization, the overall mean effect of the history of exposure is $\phantom {\dot {i}\!}\frac {1}{S+1}\sum \limits _{t_{U}=-S}^{0} \gamma _{I}^{*}(t_{U}) = \frac {1}{S+1}\sum \limits _{t_{U}=-S}^{0} w_{I}(t_{U})$ and $\frac {1}{S+1}\sum \limits _{t_{U}=-S}^{0} \gamma _{S}^{*}(t_{U}) = \frac {1}{S+1}\sum \limits _{t_{U}=-S}^{0} w_{S}(t_{U})$.

Finally, by considering an approximation of the weight functions by splines (Eq. 6) and by not constraining the weights (Eq. 7), the linear mixed model defined in Eq. 3 becomes:

$$ {}\begin{aligned} Y_{{ij}} &=\alpha_{0}~+~\tilde{X_{i}}(t_{{ij}})^{\top} \alpha_{1}~+~\sum \limits_{k=0}^{K} \theta_{{Ik}}H_{{ki}}~+~ c_{0i}~~\\ &\quad+\left(\alpha_{2}~+~\tilde{X_{i}}(t_{{ij}})^{\top} \alpha_{3}+\sum \limits_{k=0}^{K} \theta_{{Sk}}H_{{ki}}~+~ c_{1i} \right)\\ &\quad\times t_{{ij}}~+~\tilde{\varepsilon}_{{ij}} \end{aligned} $$

(8)

where θ_Ik and θ_Sk are K+1 unconstrained parameters to be estimated, and H_ki are the intermediate covariates defined from Eq. 6.

Maximum likelihood estimation using a two-stage procedure

By sharing the underlying true exposure level ${U}^{*}_{i}(t_{U})$, the sub-model for the exposure in Eq. 1 and the health outcome model in Eq. 3, define a shared parameter joint model, which would call for the simultaneous estimation of all parameters in order to avoid any bias [34]. However, given the peculiar temporal framework of the landmark approach with the inclusion of subjects present in the cohort at the landmark time and the distinct windows of observation for the two variables, a two-stage estimation procedure is unlikely to generate any bias – this assumption is confirmed by simulations in “Simulation study” section.

Let ϕ_U = (β,σ_ε,vec(B)) and $\phi _{Y} \,=\, (\alpha, \theta _{.1},..., \theta _{.K}, \sigma _{\tilde {\varepsilon }}, \text {vec}(\tilde {B})$) denote the vectors of parameters to be estimated in the linear mixed model for the exposure (in Eq. 1) and in the linear mixed model for the outcome (in Eq. 8), respectively.

First stage

In the first stage, we classically compute the maximum likelihood estimators of ϕ_U. Then, based on the estimated fixed effects parameters $\widehat {\beta }$ and the best linear unbiased predictors of the random effects $\widehat {b_{i}}$, the individual-specific true exposure can be predicted at any time t_U preceding the initial assessment of the outcome:

$$ {}\widehat{U}^{*}_{i}(t_{U}) = X_{i}(t_{U})^{\top}\widehat{\beta} + Z_{i}(t_{U})^{\top} \widehat{b_{i}}~~~~~~~~~~\forall~t_{U} \in [-S,0] $$

(9)

Second stage

In the second stage, the true underlying exposure level $U^{*}_{i}(t_{U})$ is replaced by its individual estimation $\widehat {U}^{*}_{i}(t_{U})$. Therefore, in the case of an approximation by splines (as introduced in “Specification of weights” section), ${WCIE}_{.i} = \sum \limits _{k=0}^{K} \theta _{.k} H_{{ki}}$ in Eq. 6 is replaced by $\widehat {WCIE}_{.i} = \sum \limits _{k=0}^{K} \theta _{.k} \widehat {H}_{{ki}}$ with $\phantom {\dot {i}\!}\widehat {H}_{{ki}} =\sum \limits _{t_{U}=-S}^{0} B_{k}(t_{U}) \widehat {U}^{*}_{i}(t_{U})$. The model for the outcome thus becomes:

$$ {}\begin{aligned} Y_{{ij}} =~&\alpha_{0}~+~\tilde{X_{i}}(t_{{ij}})^{\top} \alpha_{1}~+~\sum \limits_{k=0}^{K} \theta_{{Ik}} \widehat{H}_{{ki}}~+~c_{0i}\\ &\,+\,\left(\alpha_{2}+\tilde{X_{i}}(t_{{ij}})^{\top} \alpha_{3}+\sum \limits_{k=0}^{K} \theta_{{Sk}} \widehat{H}_{{ki}}+ c_{1i} \right)\times t_{{ij}}+\tilde{\varepsilon}_{{ij}} \end{aligned} $$

(10)

and parameters ϕ_Y can be again estimated by classical likelihood maximization.

In the end, the estimated time-varying effects of the pre-landmark-time exposure on the level of the outcome at the landmark time (i.e. $\widehat {\gamma }^{*}_{I}$) and on the change over time of the outcome after the landmark time (i.e. $\widehat {\gamma }^{*}_{S}$) are:

$$ \begin{aligned} ~~~~~~~~~~~~~\widehat{\gamma}^{*}_{I}(t_{U}) &= \widehat{w_I}(t_{U}) = \sum \limits_{k=0}^{K} \widehat{\theta}_{{Ik}} B_{k}(t_{U})\\ \widehat{\gamma}^{*}_{S}(t_{U}) &= \widehat{w_S}(t_{U}) = \sum \limits_{k=0}^{K} \widehat{\theta}_{{Sk}} B_{k}(t_{U}) \end{aligned} $$

(11)

Standard error estimation

One drawback of two-stage estimation approach is that the variances of the parameters obtained in the second stage do not account for the variability of the parameters estimated in the first stage. To properly take into account this variability of the first stage, we used parametric bootstrap [35]: instead of including in the second stage $\widehat {U^{*}}(t_{U})$ computed at the maximum likelihood estimate $\widehat {\phi _{U}}$, we generated M bootstrap replicates ϕ_Um (for m=1,...,M) from the asymptotic distribution $\mathcal {N}(\widehat {\phi _{U}},\widehat {V(\widehat {\phi _{U}}}))$ and computed the corresponding $\widehat {U^{*}}_{m}(t_{U})$ to be included in the second stage where $\widehat {\phi _{{Ym}}}$ and $\widehat {V(\widehat {\phi _{{Ym}}})}$ were then obtained. The total variance $\widehat {V_{{tot}}(\widehat {\phi _{Y}})}$ of the estimated parameters $\widehat {\phi _{Y}}$ was obtained as the sum of the intra- and inter-replicate variances:

$$ {}\widehat{V_{{tot}}(\widehat{\phi_{Y}})} \,=\,\! \frac{1}{M} \sum \limits_{m=1}^{M} \!\widehat{V(\widehat{\phi_{{Ym}}})} + \frac{1}{M}\! \sum \limits_{m=1}^{M} (\widehat{\phi_{{Ym}}}-\overline{\widehat{\phi_{{Ym}}}})(\widehat{\phi_{{Ym}}}-\overline{\widehat{\phi_{{Ym}}}})\!^{\top} $$

(12)

The variance of $\widehat {\gamma.}(t_{U})$ was then easily derived using the Delta-method [36].

Implementation

The methodology can be implemented in any standard statistical software. We chose R programming (version 3.6.0) and function hlme of the R package lcmm (version 1.7.8) [37] for estimating the linear mixed models. The codes of the simulation studies are openly available in GitHub at https://github.com/MaudeWagner/WHistory.

Application to body mass index history starting at mid-life and subsequent cognitive trajectories after age 70

To emphasize the utility of our methodology, we investigated in a prospective cohort study the relationship of BMI history collected since mid-life with subsequent cognitive function and cognitive decline in older ages, where prior data indicate changing relations over time, with the possibility of reverse causation at older age [3–5].

The Nurses’ Health Study

We relied on data from the Nurses’ Health study (NHS). NHS began in 1976, when 121,700 female registered nurses, aged 30-55 years and residing in 11 US states, returned a mailed questionnaire about their lifestyle and health, including their weight and height [38]. Thereafter, the participants continued to complete biennial questionnaires, with updated data. BMI was computed as self-reported weight in kilograms divided by height in meters squared; self-reported weight was highly correlated with measured weight in a previous validation study in 184 participants [39]. From 1995 to 2001, all nurses who had reached age 70 or older with no history of stroke were invited to participate in a telephone-based study of cognitive function; the entry in this sub-study corresponds to our landmark time. Among eligible women, 19,415 (92%) completed the initial validated Telephone Interview for Cognitive Status (TICS, score range 0 to 41), a telephone-based modified version of the Mini-Mental State Examination [40]. Cognitive assessments were repeated at three occasions approximately every 2 years, with a high follow-up rate (>90% among those remaining alive at each follow-up point).

Following our landmark framework (see Fig. 2), we considered the history of BMI from the entry in the cohort through the first assessment of TICS (at time 0); more than 95% of participants had information on BMI up to 24 years before the initial cognitive interview. For the current analyses, among the 19,415 participants of the cognitive sub-study, we only excluded 34 women who did not have at least one measure of BMI between enrollment and the first cognitive interview or with missing data for educational level (an important potential confounder), leading to a study sample of 19,381 women.

At enrollment, the mean age of women was 50.5 years-old (SD=2.5), 77.9% had the registered nurse diploma as highest educational level and 34.9% were overweight or obese whereas only 1.3% were underweight. At the landmark time of the initial cognitive interview, mean age was 73.3 years-old (SD=2.3), a majority of women was overweight or obese (53.4%), average TICS was 33.7 points (SD=2.8) and 11.2% of nurses had cognitive impairment. On average, 11.4 (SD=1.7) BMI measurements per person were collected over the 23.7 years (SD=1.2) of the window of exposure, followed by 3.2 (SD=1.1) TICS measurements collected subsequently over 4.9 years (SD=2.5). Figure 3 shows the observed individual trajectories of BMI in the 24 years of the window of exposure preceding the first cognitive interview and those of the TICS in the 8 years of the window of cognitive assessment in a subsample of 150 randomly selected women. In general, BMI tended to increase over time while TICS decreased.

Fig. 3

Observed individual trajectories of body mass index in the 24 years of the window of exposure preceding the first cognitive interview (left panel) and of Telephone Interview for Cognitive Status (TICS) score over the window of cognitive assessment (right panel) for 150 randomly selected women from Nurses’ Health Study, United States (1976-2008)

Specification of the statistical models

We considered the following linear mixed effect models for the BMI and the subsequent TICS score:

$$ \begin{aligned} \text{BMI}_{{il}} &=\beta_{0} + \beta_{1} \text{age0}_{i} + \beta_{2} \text{education}_{i} + F(t_{{Uil}})^{\top} \beta_{3} + b_{0i} + F(t_{{Uil}})^{\top} b_{1i} + \varepsilon_{{il}}\\ \text{TICS}_{{ij}} &=\alpha_{0} + \alpha_{1}\text{age0}_{i} + \alpha_{2} \text{education}_{i} + \sum \limits_{k=1}^{K} \theta_{{Ik}}H_{{ki}} + c_{0i} + \alpha_{3} \text{V0}_{{ij}} \\ &\quad +\left(\alpha_{4} + \alpha_{5} \text{age0}_{i} + \alpha_{6} \text{education}_{i} + \sum \limits_{k=1}^{K} \theta_{{Sk}}H_{{ki}} + c_{1i}\right)\times t_{{ij}} + \tilde{\varepsilon}_{{ij}} \end{aligned} $$

(13)

where age0 is the age of women at enrollment into the cognitive sub-study in years; education is the highest degree of diploma (binary variable: Registered nurse versus Bachelor’s, Master’s or Doctorate); V0 is an indicator for the first cognitive assessment which captures a possible first-passing effect [41]; F(t_Uil) is a basis of natural cubic splines with 4 inner knots located at 20th, 40th, 60th and 80th percentiles of the elapsed time of exposure (the number of knots was determined by the AIC); and H_ki are the intermediate covariate summarizing BMI history: $\phantom {\dot {i}\!}H_{{ki}} = \sum \limits _{t_{U}=-S}^{0} B_{k}(t_{U}) \text {BMI}_{i}^{*}(t_{U})$. The weight function of the WCIE was approximated using natural cubic splines with K−1=2 inner knots located at 33th and 66th percentiles (the number of knots was determined by AIC). β,α,b_i,c_i,ε_il, and ε_ij have been previously defined in “Methods” section. In order to facilitate the interpretation of the results and as introduced in the methods, we assumed linear trajectories of cognitive decline. This choice seemed reasonable since cognition was evaluated in the NHS over a short follow-up.

Results

Figure 4 represents the trajectory of association of BMI pre-landmark-time history in the 24 years with the initial level of TICS (left panel) and its change over time (right panel). Each point represents the association coefficient of BMI at a specific year adjusted for age, educational level, and all BMI history at any other year.

Fig. 4

Trajectories of associations between body mass index (BMI) history in the 24 years prior to the first cognitive interview on the initial level (left panel) or the change with time (right panel) of the Telephone Interview for Cognitive Status (TICS) score in the Nurses’ Health Study (N=19,381), United States (1976-2000). 95% confidence intervals were obtained by parametric bootstrap with 500 replicates. A negative estimate indicates that increased BMI is related to worse cognition/more cognitive decline and a positive estimate indicates better cognition/less cognitive decline

The overall mean association of BMI history over the whole 24 year period of exposure was significant with a negative relation to cognition for both the initial level (-0.0013 [95% CI: -0.0014;-0.0012]) and the annual slope of decline in TICS (-0.00016 [95% CI:-0.00019;-0.00015]). In addition, as expected, the association between BMI exposure history and cognition was non-constant, with shapes of associations on both the initial level and the slope of TICS suggesting opposite remote and recent effects. For the initial level (see Fig. 4, left panel), higher BMI from -24 to -20 years prior to cognitive assessment, corresponding to mid-life, were significantly associated with lower mean TICS at the first cognitive interview assessed after age 70. For example, for similar confounders (i.e., age and education level) and BMI trajectories at any times except at -24 years, a 1-kg/m² increase of BMI was associated with a lower initial TICS level of 0.025 point on average (95% CI = -0.039; -0.009). Between -20 and -12 years, for similar confounders and BMI trajectories at any times except at the year evaluated, BMI levels were no longer significantly associated with initial TICS level whereas, from -12 to -5 years, higher BMI levels were again associated with worse cognitive function, consistently with that observed earlier in midlife. In contrast, starting around 5 years before cognitive assessment, higher BMI levels became significantly associated with higher initial level of TICS (likely reflecting reverse causation).

For the slope of TICS (see Fig. 4, right panel), results showed no significant association with BMI levels from -24 to -13 years preceding the first cognitive interview. However, higher BMI levels between -12 and -5 years and lower BMI levels between -5 years and the first cognitive interview were associated with worse cognitive decline, similarly to what observed with the initial level of TICS (but to a lesser extent).

To help interpret these complex trajectories of association, we examined hypothetical BMI trajectories (e.g., stable BMI of 25 kg/m² over the 24 years of exposure) in relation to subsequent changes in TICS (see Fig. 5). As observed with the trajectories of association, these predictions by profile of BMI generally showed that women with greater BMI over time had higher (better) initial mean levels of TICS compared to women who had a stable or a decreasing BMI with age, reflecting as expected the major influence of reverse causation. Overall, this application provides additional evidence that relationships between BMI and cognition are complex and largely depend on careful consideration of the window of exposure when BMI is assessed.

Fig. 5

Mean change of the TICS score over the study course following 5 theoretical BMI history profiles defined all over the whole 24 years preceding the first cognitive interview: (i) linear decrease of 5 points of BMI (purple lines), (ii) linear decrease of 10 points of BMI (blue lines), (iii) stable BMI (dashed lines), (iv) linear increase of 5 points of BMI (pink lines), (v) linear increase of 10 points of BMI (green lines). For example, on the top left panel, women who had a BMI of 25kg/m² 24 years prior to the first cognitive assessment and linearly dropped to a BMI of 20kg/m² at the initial cognitive interview (i.e. legend labelled “25 >20”) have an initial average TICS level of 33.8 points that decreases over time to 32.4 points after 7 years of follow-up

When considering piecewise constant weights instead of weights approximated by natural cubic splines, the trajectories of time-varying effects remained similar (see eFigure 1) suggesting that the splines function well reflected the data. Moreover, results remained generally the same after adjusting the linear mixed model of the outcome for additional potential confounders (i.e., chronic disease history, smoking status or postmenopausal hormone therapy) categorized as binary variables and collected over the same period as BMI history. Likewise, we observed similar trajectories of associations of BMI on cognitive function when we stratified analyses for these factors (results not shown).

Simulation study

Overview of the simulation design

To validate our methodology and its two-stage estimation procedure, we generated cohort data that mimicked the NHS design and temporal framework of Fig. 2, according to a joint model with the shared true exposure as introduced in “A landmark approach for assessing the dynamics of association between an intermittently evaluated time-varying exposure and a subsequent health outcome” section. In the main simulation study, we generated individual-specific visit times for exposure, using a uniform distribution in [-6, 6] months around theoretical visits every 2 years from -24 years to landmark time 0. At each visit time, we simulated the observed value of the exposure U according to a mixed model as defined in Eq. 13 and considered a 10% proportion of missing visits completely at random (as observed in NHS). The trajectory of U was modelled according to natural cubic splines (with 4 inner knots at 20th, 40th, 60th and 80th percentiles of the observation times) at both the population and individual levels (with correlated individual random intercept and slope which modelled the within-subject correlation over time), and was adjusted for two covariates mimicking age at study entry (normal distribution: mean 51, standard deviation 3) and binary education level (0.25-probability Bernoulli distribution).

We also generated individual-specific visit times for the outcome with a uniform distribution in [-6, 6] months around theoretical visits every 2 years from 0 to 7 years. We then generated the values of the health outcome Y according to the mixed model as defined in Eq. 13 adjusted for age and education (on both the intercept and the slope) and the indicator for the first outcome assessment, and considering a 3% proportion of missing visits completely at random (as observed in NHS). In this model, the true underlying exposure level was considered.

Scenarios

Parameters in these generating models corresponded roughly to those obtained in the application data for BMI and TICS (parameters provided in eTables 1 and 2) with exception for the trajectory of association with BMI history for which we assumed three different relevant scenarios plotted in Fig. 1:

• Scenario A: constant negative association over the time window of exposure (this corresponds to a standard CIE) with γ_I(t_U)=−0.05 and γ_S(t_U)=−0.01 for all t_U∈[−24,0];

• Scenario B: negative association far upstream of the outcome assessment and null association at the approach of the initial assessment using a truncated centred normal distribution: $\gamma _{I}(t_{U}) = \left (\Phi \left (\frac {{t_{U}}+24}{6}\right) - 1\right)\times 0.1$ and $\gamma _{S}(t_{U}) = \left (\Phi \left (\frac {{t_{U}}+24}{6}\right) - 1\right)\times 0.03$ with t_U∈[−24,0];

• Scenario C: trajectories of association obtained in the application between BMI history and TICS trajectory (see Fig. 4).

In the main simulations, we considered a framework close to NHS with frequent exposure data (roughly every two years), a small proportion of missing visits (10% for the exposure and 3% for the outcome), and a small measurement error for the exposure (σ_ε = 0.9). However, to further evaluate the methodology in less favorable situations, we also considered cases where (i) the exposure was measured every 4 years, (ii) the proportion of missing data was larger for the exposure (20%), and (iii) the error of measurement was larger (σ_ε = 1.8).

Results

We systematically considered 500 replicates of samples of 1,000 participants each. For each scenario, we focused on the estimated trajectories of association with both the initial level and the slope of the repeated outcome, and on the coverage rate of its pointwise 95% confidence interval obtained by parametric bootstrap. As shown in Fig. 6 for Scenario B, the two-stage estimation approach retrieved without bias the true generated trajectory of association, and the parametric Bootstrap provided satisfying coverage rate of the 95% confidence interval around the 95% nominal value. The same conclusions were drawn for Scenario A and C with other shapes of trajectories of association (see eFigure 2 and eFigure 3, respectively).

Fig. 6

Boxplots of the trajectory of association between the exposure history over the 24 years prior to the initial health outcome assessment on the initial level (top left panel) or the slope (top right panel) of the outcome of interest across 500 simulations of 1000 subjects each, and corresponding coverage rates of the 95% pointwise confidence interval (lower panels) for Scenario B (distant negative effect)

Furthermore, when considering less repeated measurement points for exposure (see eFigures 4, 5, 6), higher rate of missing data (see eFigures 7, 8, 9), higher measurement error (see eFigures 10, 11, 12), or higher number of inner knots of cubic splines in the definition of the BMI history (see eFigure 13), parameters were again well estimated with negligible bias and no departure from the expected 95% coverage rate of the 95% confidence interval. Overall the simulation study demonstrated that in this specific temporal landmark framework, the two-stage procedure combined with parametric bootstrap provided a correct inference.

Discussion

We proposed a flexible approach for estimating the trajectory of association between a time-varying exposure of any kind (e.g., related to lifestyle, occupation, environmental toxins), and a subsequent health outcome when exposure and outcome are assessed non-concurrently. We relied for that on the WCIE methodology [8, 9] that we incorporated within a landmark approach to handle the non-concurrent windows for the exposure and outcomes. We followed previous works in WCIE and directly estimated the weights assigned to past exposures from the data using flexible natural cubic splines [31]. However, in contrast with most previous WCIE approaches that are restricted to complete error-free exposure, our approach considers that the exposure is endogenous, prone to measurement error, and assessed at discrete and individual-specific times; this naturally handles the intermittent missing visits usually encountered in cohort studies and limits the biases induced by neglecting the measurement error that is most likely present in many epidemiological contexts. Mauff et al. [20] had previously extended WCIE to handle endogenous exposure data. However, their method was limited to concomitant exposure and outcome assessments while in life-long epidemiology, exposures remote from the outcome risk-period are also of great interest. In addition, they had only considered time-to-event outcomes while we wanted to also consider repeated outcomes that are frequently encountered in longitudinal studies. Danieli et al. [22] have just proposed a WCIE methodology for repeated outcome data. However, their method is limited to exogenous complete error-free exposures, and evaluates the cumulative effect of the exposure on the current level of the outcome while in etiological studies such as ours, the interest is usually in the effect of the exposure on the rate of change [42].

Other methods were proposed in the literature to assess the association between two longitudinal processes. By modeling the association through correlated individual random effects, the bivariate mixed-effect model [43] can describe the correlation structure between the exposure and the outcome but usually fails at providing a proper evaluation of the association of the exposure on the outcome. Cross-lagged models (CLM) [44] and causal dynamic models (CDM) [45] circumvent this problem by directly analyzing the effect of the exposure assessed at a current time on the subsequent trajectory of the outcome. While CLM focuses on the outcome level and is intrinsically linked to the discrete visit process, CDM focuses on the change of the outcome in continuous time and is thus to be favored in etiological studies. However, to the best of our knowledge, none of them were extended to acount for the history of the exposure through a WCIE as we proposed, and both methods are limited to concomitant exposure and outcome assessments. Thanks to the separation between the assessment periods of the exposure and the outcome in the landmark approach, we were able to correctly estimate the time-varying effects of the exposure history using a two-stage approach, as demonstrated by simulations. This two-stage procedure makes our methodology easy to implement in standard statistical software, and easy to adapt to other contexts as long as exposure and outcome are not assessed simultaneously. First, any type of outcome could be considered such as time-to-event or binary endpoint as usually encountered in case-control studies. Note however that (i) in case-control studies (unlike prospective studies or case-control studies nested within a cohort), exposure data are collected retrospectively from the landmark time, when the groups were selected, which may result in non-differential measurement errors; (ii) with a time-to-event endpoint, our approach would still require that the exposure time- window ends at the landmark time. Second, non-ignorable dropout for the repeated outcome could be easily taken into account by estimating the outcome parameters in the second stage using a shared random-effect model [46] rather than a standard mixed model. Third, with a repeatedly-measured outcome, any relevant nonlinear shape of trajectory over time could be considered instead of a linear trajectory. Fourth, the approach could be easily adapted to handle other types of exposures, such as categorical variables, by considering a generalized linear mixed model in the first stage. Fifth, the NHS women were highly homogeneous in age and we used the same time scale (i.e., elapsed time since entry in cognitive substudy) to (i) model the underlying true exposure and (ii) define the WCIE; however, in other applications, users could consider other timescales, and even two different timescales, such as age for the exposure trajectory and years before the landmark for the weights. Finally, for ease of epidemiological interpretation we defined the WCIE as the sum of weighted exposures at each time unit over the [−S,0] period but the integral over the continuous time window can also be considered, and leads to the same trajectory of association as illustrated in supplementary eFigure 14 with the trajectory of BMI association in the NHS. Simulation studies showed that the proposed two-stage estimation model recovers various sets of shapes of the true weight function (i.e., constant over time or time-varying), and provides satisfactory estimates of the strength of the associations. In particular, splitting the estimation procedure into two parts does not seem to lead to bias when the exposure and outcome are non-concurrent, even in the case of poorer exposure information (i.e., higher missing data or higher time-interval between two visits or higher measurement error).

The method presents however several limitations that warrant consideration. First, while our estimates are capable of capturing the shape of the weight function, the splines functions may not accurately reflect rapid and sudden changes in the relationship [29]. If these effects are expected, the user can improve the model by including one or more knot(s) in regions where a strong curvature is anticipated or by considering a piecewise constant trajectory of association instead of splines. Second, our methodology assumes a linear association between the exposure at each time and the outcome. Extensions to account for nonlinear associations both over time and over the exposure range is a direction for future research. Third, we considered a standard linear mixed model for the exposure which assumes zero-mean homoscedastic Gaussian independent measurement errors. Additional simulations showed that the method remained robust to errors with heavier tails of distribution when these errors were centered around 0 (see eFigure 15) but the method was likely to provide biased estimates in the case of errors not centered around zero (see eFigure 16). In epidemiological contexts where the distribution of the errors is expected to be non-standard (e.g., when an under-reporting is likely), the linear mixed model should thus be revised accordingly. Finally, as in previous WCIE methods, we only considered one time-varying exposure while it could be of interest in the future to simultaneously examine several exposure histories that are known to be interrelated (e.g., BMI and physical activity, or multi-pollutant/chemical mixtures).

Applied to the association between BMI history starting at midlife and cognition at older ages, our methodology confirmed the complex and age-dependent relationship. Indeed, with no a priori assumptions on the shape of relations, and after controlling for BMI history, age and educational level, we found that higher BMI levels at midlife but lower BMI levels at older ages were significantly associated with poorer cognition after the age 70 in women. This bidirectional relationship is consistent with previous studies, including our own work, that showed that obesity at middle age but low BMI late in life were associated with subsequent development of cognitive impairment and dementia [3–5, 47, 48]. Note however that our interpretations are based on a specific population of generally healthy women and that these results cannot be generalized to populations with other socioeconomic or health status. Although a non-linear J-shaped association (i.e., underweight and overweight are associated with an increased risk of worst cognition) was sometimes reported for BMI in the literature [47], this is unlikely to be the case in our sample in which very few nurses were underweight during the whole exposure period. Finally, we note that the follow-up rate was high in the NHS and there was little or no attrition so that we did not need to account for death or possible informative dropout in our application.

Conclusions

This methodology offers great flexibility in capturing complex dynamic relationships over the life-course, but also allows application to the majority of epidemiological contexts when the shape of the effects is not known and even more importantly when the history of exposure is measured incompletely and with error. Overall, this method may significantly contribute to broadening the applications of WCIE for a variety epidemiological contexts.