Причинно-следственная корреляция, обусловленная неоднократными измерениями воздействия и исхода в продольных параметрах.

Abbreviations

directed acyclic graph

expanded disability status scale

iterative conditional expectation

inverse probability weighting

longitudinal targeted maximum likelihood estimator

multiple sclerosis

marginal structural model

Observatoire Français de la Sclérose en Plaques

primary progressive

root mean squared error

relapsing–remitting

secondary progressive

Introduction

In longitudinal settings with time‐dependent confounding between a time‐varying exposure and outcome, the causal inference framework provides an efficient approach to estimate an unbiased causal effect of the exposure on the outcome [ 1 ] using a counterfactual definition of causality in combination with g‐method estimators: g‐computation, inverse probability weighting (IPW), or doubly robust estimators such as the longitudinal targeted maximum likelihood estimator (LTMLE) [ 2 , 3 , 4 , 5 , 6 , 7 ]. However, in this longitudinal setting, repeated measures of exposure and outcome in the same patient is likely to lead to intra‐patient correlation for most of the measured variables, even after accounting for observed confounders. This intra‐patient correlation represents all the unmeasured causes common to the patient and can be represented in the causal directed acyclic graph (DAG) as latent variables that will affect exposure or outcome at each time point. Regardless of the absence of unmeasured confounding between the exposure and the outcome, ignoring this intra‐patient correlation theoretically yields a biased estimate of the causal effect of the exposure on the outcome.

This issue of intra‐patient correlation in a longitudinal context is sometimes overlooked in causal reasoning, and the impact of this omission has, to the best of our knowledge, not yet been explored in the causal inference literature. Mixed‐effects models can be seen as a relevant solution for dealing with intra‐patient correlation, using random effects to capture individual‐specific variability across repeated measures [ 8 ]. In the causal inference field, Shardell and Ferruci addressed intra‐patient correlation with random effects in joint models but primarily in the context of unmeasured confounding [ 9 ]. Other studies have explored the use of mixed‐effects models in causal analysis, mainly to answer mediation questions in longitudinal settings, but, although the importance of accounting for intra‐patient correlation is discussed, this was not a central aim of these studies [ 10 , 11 ]. However, no study has specifically focused on how intra‐patient correlation affects causal inference in longitudinal data, even when all confounders are adjusted for.

The work presented herein stems from a clinical question concerning the causal effect of pregnancy on the course of neurological disability in multiple sclerosis (MS) patients. The relationship between the occurrence of a pregnancy as an exposure and the neurological disability as an outcome is subject to time‐dependent confounding: patients with a higher disability have a lower probability of being exposed to pregnancy, thereby introducing a bias when assessing the causal effect of pregnancy on the subsequent disability course. However, each patient has a specific disability course reflecting the severity of her disease, which results in an intra‐patient correlation for successive disability measures; and the same is true for the probability of pregnancy in the same patient.

The novel contribution of this work lies in demonstrating the importance of accounting for intra‐patient correlation in the causal reasoning when dealing with longitudinal data with repeated measures, and to compare two ways of adapting g‐methods to account for it, using either fixed‐effects or mixed‐effects models. This was first evaluated on simulated datasets based on different scenarios, then on real‐life observational data from the French MS registry, the Observatoire Français de la Sclérose en Plaques (OFSEP), using as a demonstrative example the study of the causal effect of pregnancy on disability course in MS.

Notations and Assumptions

Time t 0 represents the baseline and time T the last follow‐up. For each discrete time t = t 0 , … , T , let Y t denote a quantitative outcome (disability), E t a binary exposure (pregnancy), and L t the vector of all possible confounding factors evaluated during regularly spaced visits. For time‐dependent variables, let Y ‾ t , E ‾ t and L ‾ t denote the history of, respectively, the outcome, the exposure and the time‐dependent confounding‐factors up to time t , represented by the vector of all values Y t 0 , Y t 1 , … , Y t , E t 0 , E t 1 , … , E t and L t 0 , L t 1 , … , L t taken by the variables from time t 0 to time t . The counterfactual outcome with exposure E = e is denoted Y e and corresponds to the value of the outcome Y that would have been observed if, possibly contrary to fact, the patient had an exposure e .

The aim of g‐methods is to estimate the mean counterfactual outcome E Y T e ‾ under a counterfactual exposure history. They rely on three assumptions for all values of time t :

(i) Sequential conditional exchangeability,

that is, the counterfactual outcome under exposure history e is independent of the observed exposure at time t conditional on the confounding and exposure histories:

(ii) Positivity,

that is, all exposures must have a non‐null probability conditional on the past exposure and covariates histories;

(iii) Consistency,

that is, the exposure is well defined, exists in only one version, and the counterfactual outcome under an observed exposure history corresponds to the outcome represented in the data [ 1 ].

Motivating Example: Pregnancy and Neurological Disability in Multiple Sclerosis

MS is the most common inflammatory disease of the central nervous system and the leading cause of non‐traumatic disability in young people [ 12 ]. MS patients suffer from increasing neurological disability over time, related to the sequelae of disease relapses (relapsing–remitting [RR] form) and progressive neurodegeneration (primary progressive [PP] or secondary progressive [SP] form). Neurological disability is assessed clinically by the Expanded Disability Status Scale (EDSS) that is based on neurological symptoms, clinical examination, walking ability, and loss of autonomy; scores range from 0 to 10, in increments of 1 point from 0 to 1, then 0.5 point [ 13 ]. As MS mainly affects young women, the question of the causal effect of pregnancy on the course of the disease is a major concern [ 14 , 15 , 16 ], but addressing this question is complicated by reverse causation, as patients with higher disability are less likely to have a pregnancy [ 17 ].

The causal DAG shown in Figure  1 represents a simplified hypothetical causal relationship between the outcome Y t and exposure E t from time t 0 to time t 3 : the occurrence of a pregnancy affects the subsequent disability (causal effect of interest E t → Y t ), the disability affects the probability of pregnancy at the next time point (time‐dependent confounding Y t → E t + 1 ), outcome at time t affects the subsequent outcome at time t + 1 ( Y t → Y t + 1 ), and exposure at time t affects the subsequent exposure at time t + 1 ( E t → E t + 1 ). For the sake of simplicity, the effect of the exposure on the outcome was restricted to the effect at the same time t with no sustained effect (i.e., effect of exposure E t on outcome Y t + 1 or subsequent outcomes), and no other baseline or time‐dependent covariates were included, nor were censored or missing data.

DAG representing the relationship between the outcome Y t and the exposure E t from time t 0 to time T  = 3, with intra‐patient correlation due to latent variables U Y and U E .

The study of pregnancy as an exposure presents several specific features compared to other exposures such as treatment. Pregnancies occur a limited number of times in a patient's history, unlike a treatment that may be taken at many time points, thus the counterfactual exposure history e ‾ 1 “exposed at all times” is meaningless for pregnancy, and only the counterfactual exposure history e ‾ 0 “never exposed” is relevant. The causal estimand was therefore the causal mean EDSS difference at the end of the study period between the observed situation with pregnancies as observed e ‾ obs and the counterfactual situation without any pregnancy e ‾ 0 in the entire study population:

Note that, it could also be of interest to report the difference between the observed mean outcome and the observed counterfactual outcome in the exposed sub‐population (as E Y T e ‾ obs = E Y T e ‾ 0 in the unexposed population), analogous to an average treatment effect in the treated.

Methods

Analysis of the Causal DAG in the Presence of Intra‐Patient Correlation

The presence of intra‐patient correlation between repeated measures for the outcome and the exposure theoretically introduces causal bias into the estimation of the effect of the exposure on the outcome if it is not accurately accounted for. This intra‐patient correlation is related to a set of common causes, specific to each patient and not measured, that will impact the outcome or the exposure values at each time point. It can be represented on the DAG as shown in Figure  1 by adding latent variables U Y and U E that respectively affect the outcome Y and the exposure E at each time. Importantly, these latent variables do not represent unmeasured confounding between the exposure and the outcome (which would be represented by a single latent variable that would simultaneously affect the exposure and the outcome).

In absence of latent variables U Y and U E , adjusting only on the last outcome value Y t − 1 is sufficient to provide an unbiased estimate of the effect of exposure E t on the outcome Y t at time t (by closing the non‐causal pathway E t ← Y t − 1 → Y t ). The presence of latent variables U Y and U E opens several non‐causal pathways between exposure and outcome at time t , as shown in Figure  2A : U Y is a common cause of the exposure (through its effect on past outcomes, and thus on past exposures due to time‐dependent confounding) and of the outcome (e.g., E t ← E t − 1 ← Y t − 2 ← U Y → Y t ), but adjusting on the exposure at time t − 1 is not sufficient to block all non‐causal paths as it opens a path through U E due to a collider effect (e.g., E t ← U E → E t − 1 ← Y t − 2 ← U Y → Y t ). Note that non‐causal pathways are still open even with the assumptions of no effect of exposure at time t on subsequent exposure at time t + 1 and of outcome at time t on subsequent outcome at time t + 1 .

DAG representing adjustments required to obtain an unbiased effect of the exposure E t on the outcome Y t at time t . Without latent variables U Y and U E , adjusting on the outcome at last time Y t − 1 is sufficient, but the presence of latent variables U Y and U E opens up several non‐causal pathways (A). The fixed‐effects approach requires blocking all non‐causal pathways by adjusting for all past outcome and exposure values in multiple fixed‐effects models (B). For the mixed‐effects approach, non‐causal pathways can be blocked by adjusting directly for latent variables of intra‐patient correlation, estimated through a random effect: Either in a single mixed‐effects outcome model for U Y , also adjusted for Y t − 1 (C), or in a single mixed‐effects exposure model for U E , also adjusted for Y t − 1 and E t − 1 (D).

Two Approaches to Account for Intra‐Patient Correlation

The Fixed‐Effects Approach

The first and more conventional approach to block all non‐causal pathways in the relationship between the exposure E t and the outcome Y t is to adjust for the entire past exposure history E t 1 , E … , E t − 1 and outcome history Y t 0 , Y … , Y t − 1 , as depicted in Figure  2B [ 6 , 18 ]. This approach can be done using only fixed‐effects in an outcome or an exposure model, hence we refer to it as “fixed‐effects approach” hereafter. It requires fitting a separate outcome or exposure model for each time point, incorporating additional parameters for the effect of the exposure and the outcome at the last time point.

The Mixed‐Effects Approach

The second approach consists of taking into account directly the intra‐patient correlation induced by latent variables through a patient random effect in a mixed‐effects model; we refer to it as the “mixed‐effects approach” hereafter. There are two possibilities for blocking non‐causal paths here: adjusting for U Y through a random effect in an outcome model also adjusted for Y t − 1 through a fixed effect (Figure  2C ), or adjusting for U E through a random effect in an exposure model adjusted for Y t − 1 and E t − 1 through fixed effects (Figure  2D ). In both cases, there is no need to further adjust for exposure and outcome at time t − 2 or earlier.

The mixed‐effects approach relies on the assumptions that the intra‐patient correlation is time‐invariant, exogenous, and can be estimated through a random individual effect, with all the assumptions of mixed‐effects models (normal distribution of the random effects, and independence of random effects and residual errors). This mixed‐effects approach involves that all the effects of exposure, outcome or time‐dependent covariates are accounted for at time t − 1 only, either with the assumption that effects at time t − 2 , …, t 0 pass through the effect on the exposure, outcome or covariate at the following time with no direct effect, or by summarizing the history of exposure, outcome or covariate at time t − 1 (e.g., with a cumulative sum of all past values). Instead of estimating new parameters for each time point, the mixed‐effects approach assumes that parameters are time‐invariant or allow for a time dependence through an interaction with a time covariate.

Adaptation of G‐Methods to Handle Intra‐Patient Correlation

G‐Computation

G‐computation [ 18 , 19 , 20 , 21 ] is a parametric version of the g‐formula that recovers the distribution of confounders under a counterfactual exposure to predict the counterfactual outcome:

The first part of the formula is approximated by an outcome model, called the model 𝒬 , which aims to predict counterfactual outcome values Y t e ‾ t for each patient as a function of the counterfactual exposure history e ‾ t and the confounding factors history L ‾ t − 1 , which, in the example used herein, is the past outcome history Y ‾ t − 1 . The second part of the formula requires an integration over the distribution of confounders history under the counterfactual exposure history, which, in the example used herein, corresponds to the mean of the counterfactual outcome predicted at the preceding time point. In some cases, the use of the mean counterfactual value is not sufficient and the entire counterfactual distribution of confounders or exposure must be approximated by a Monte Carlo approach; for example in the presence of a binary outcome, a binary time‐dependent confounder or a dynamic counterfactual binary exposure, when non‐collapsibility occurs [ 1 ].

To adapt g‐computation in the presence of intra‐patient correlation using the fixed‐effects approach, the outcome model 𝒬 must be adjusted for all past values of exposure and outcome back to time t 0 . This requires fitting a new outcome model 𝒬 t for each time point t , that includes all the past outcome and exposure values with a fixed‐effect, in the form:

where α 0 t is the intercept, α 1 t the parameter associated with the last exposure value, and α 2 t , k and α 3 t , k are the parameters associated one by one with each past outcome and exposure values, each new outcome model 𝒬 t having its own set of parameters.

For the mixed‐effects approach, a single mixed‐effects outcome model 𝒬 ( t ) is fitted over all time points in order to estimate intra‐patient correlation related to the latent variable U Y through a patient random effect a 0 i . The model takes the following form:

where α 3 . f ( t ) represents the effect of time t , to which a function f ( . ) can be applied to model non‐linear functional forms (e.g., smoothing function such as splines), or using a saturated model with one parameter per time point; α 4 . f ( t ) . Y t − 1 and α 5 . f ( t ) . E t represents the time‐varying effects of past outcome and past exposure history; and a 0 i is a random patient intercept.

Inverse Probability Weighting

For the IPW approach in a longitudinal setting [ 2 , 4 , 22 ], individuals are weighted by the inverse of their probability of having the observed exposure history E ‾ t conditional on the values of their confounding factor history L ‾ t , obtained by multiplying the stabilized weights at each time point SW t e :

The weighted pseudo‐population contains one sub‐population per possible exposure history value E ‾ in which the relationship between outcome and exposure is no longer biased by confounding factors. Thus, the weighted mean value of the outcome E Y t in the sub‐population with the exposure history E ‾ = e ‾ corresponds to the counterfactual outcome E Y t e ‾ t :

The weights can also be truncated to reduce the extreme values to lower values, thus decreasing the variance at the cost of introducing bias into the estimation [ 23 ].

The parametric IPW is based on an exposure model, the g model, to predict the probability of exposure at each time point given the past exposure and confounder histories. In addition, a marginal structural model (MSM) can also be fitted to summarize the association between the counterfactual outcome and the exposure in the weighted pseudo‐population, but in the motivating example herein, only one counterfactual history exposure was of interest (“without pregnancy”), so an MSM was not necessary as the weighted mean outcome in the non‐exposed population was sufficient.

To adapt IPW in the presence of intra‐patient correlation, the reasoning is analogous to that of g‐computation presented above. Using the fixed‐effects approach, the exposure model must be adjusted for all past values of outcome and exposure with fixed‐effects, and thus a new exposure model g t have to be fitted at each time, in the form:

For the mixed‐effects approach, a single exposure model g ( t ) is constructed including all time points in order to estimate a patient random intercept b 0 i that aims to approximate the intra‐patient correlation related to the latent variable U E , as well as fixed‐effects for the outcome and exposure value at time t − 1 only, in the form:

Longitudinal Targeted Maximum Likelihood Estimator

LTMLE is a doubly robust estimator that combines an outcome‐based estimator and an exposure‐based estimator into a single one, with the advantage of being unbiased if at least one of the two methods is correct [ 7 ]. The first step consists of a direct estimation using an outcome model, followed by a targeting step that combines the initial estimation with a clever covariate derived from an exposure model, evaluated using a loss function: if the exposure model improves the outcome predictions, it will be included in the counterfactual predictions; otherwise, only the initial estimation will be used.

For the outcome‐based estimator, and contrary to the classical description of the LTMLE algorithm and the ltmle package [ 24 ], we used the traditional g‐computation described above rather than the iterative conditional expectation (ICE) estimator because the implementation of mixed‐effects model in g‐computation is more straightforward. Indeed, the ICE estimator is based on a sequential algorithm: starting from the last time and moving backward to the first time, the new outcome model that must be fitted at each time t takes the last outcome model prediction at time t + 1 as its dependent variable. A single outcome model across all time points cannot be fitted with ICE, so no random effect can be estimated.

Simulation Study

We conducted a simulation study under different scenarios, first, to evaluate the impact of intra‐patient correlation on g‐method estimators in the context of longitudinal data with repeated measures, and second, to assess and compare the performance of g‐computation, IPW, and LTMLE to account for it using either a fixed‐effects or a mixed‐effects approach.

Designs of Simulations

All analyses were performed using R software, version 4.0.3 [ 25 ]. Data were simulated using the simcausal package [ 26 ], according to the DAG presented in Figure  1 , with the aim of obtaining simulations close to the data expected in real‐life data. R code for the data simulation and analysis is provided in Appendix  A .

Outcomes Y t were generated from an autoregressive linear model, depending on the value of last outcome Y t − 1 , the exposure E t , the measure time t , the latent variable U Y (normal distribution with a mean of zero and a fixed standard deviation σ U Y ), and the residual ε t (normal distribution with a mean of zero and a fixed standard deviation). Exposures E t were generated from a logistic model as a function of the last outcome Y t − 1 , last exposure E t − 1 , time t , and a latent variable U E (normal distribution with a mean of zero and fixed standard deviation σ U E ). The true counterfactual outcomes Y T e ‾ 0 were obtained using the same outcome generative model with exposure set to zero at each time point for all individuals. The model formulas and coefficient values used for the data generation are presented in Table  S1 .

Several simulation scenarios were defined, with different numbers of time points ( T  = 5, 10, and 20) and patients ( n  = 200 and 500), under the null hypothesis (no effect of the exposure E t − 1 on the subsequent outcome Y t ) or the alternative hypothesis (beneficial effect of pregnancy on disability), without latent variable U Y and U E ( σ U Y  = 0 and σ U E  = 0) and with latent variables (respectively, σ U Y  = 0.25 and 0.50, and σ U E  = 0.25 and 0.50 on a logit scale). All simulations were performed with time‐dependent confounding, introduced by setting a negative effect of the outcome Y t − 1 on the exposure probability E t .

Analyses of Simulated Data

Simulated data were analyzed using g‐computation, IPW and LTMLE with a manual implementation, and all outcome and exposure models were correctly specified. IPW and LTMLE were conducted with stabilized weights and with truncation of predicted probabilities below 0.05. Confidence intervals were obtained by bootstrapping on 500 replicates. The model formulas for the outcome 𝒬 ‐model and the exposure g model are provided in Table  S2 .

Evaluation of Estimators

The different estimators aim to evaluate the causal estimand, which was the mean outcome difference at the last time T between the observed situation and the counterfactual situation in the absence of exposure ∆ Y = E Y T − E Y T e ‾ 0 . The estimators were evaluated by their mean bias, relative bias under the alternative hypothesis, standard error and root mean squared error (RMSE) from 1000 simulations. As EDSS at last time point depended on many parameters (number of time points, causal effect of pregnancy on EDSS, pregnancy rate, and intra‐patient correlation U Y and U E ), the true value of the causal estimand ∆ Y in each simulation scenario was obtained empirically by simulation on a very large number of patients. The coverage probability of the 95% confidence interval was estimated under certain conditions only ( T  = 10, n  = 500, with latent variables σ U Y  = 0.50, and σ U E  = 0.50, under null and alternative hypotheses) due to the long computational time required.

Results of the Simulation Study

The results of g‐computation, IPW and LTMLE for n  = 500 individuals (mean estimate and 95% fluctuation interval over 1000 simulated datasets) are presented in Figure  3 , and the exhaustive results of absolute biases, standard errors, RMSE, and coverage probabilities for n  = 500 and n  = 200 are reported, respectively, in Figures  S1 and S2 .

Mean estimate and 95% fluctuation interval of g‐computation, IPW and LTMLE over 1000 simulated datasets. For n  = 500 individuals, under the null hypothesis, for T  = 5, 10, and 20 time points, without latent variables ( σ U Y  = 0 and σ U E  = 0), and with latent variables ( σ U Y  = 0.25 or 0.50 and σ U E  = 0.25 and 0.50), with UY and UE respectively the standard deviation of the random effect on the outcome U Y and on the exposure U E . The scale of the y‐axis differs according to the number of time points. GC, g‐computation; IPW, inverse probability weighting; LTMLE, longitudinal targeted maximum likelihood estimator.

Without Intra‐Patient Correlation

In the absence of intra‐patient correlation ( σ U Y  = 0 and σ U E  = 0), g‐computation provided an unbiased estimate when not corrected for latent variables or with a fixed‐effects approach (absolute bias < 0.01 and relative bias < 1% under all scenarios) and a small bias with a mixed‐effects approach (maximum absolute bias of 0.03 and relative bias of 4%). IPW and LTMLE were also unbiased when not corrected for latent variables or with a fixed‐effects approach (absolute bias < 0.01 and relative bias < 1% under all scenarios) for T  = 5 and 10 time points, whereas a small bias was found for T  = 20 (maximum absolute bias of 0.02 and relative bias of 2% for IPW, and maximum absolute bias of 0.03 and relative bias of 3% for LTMLE) or with a mixed‐effects approach, regardless of the number of time points (maximum absolute bias of 0.03 and relative bias of 4% for IPW, and maximum absolute bias of 0.04 and relative bias of 5% for LTMLE).

Effect of Intra‐Patient Correlation When Not Accounted for

When the latent variable U Y and U E were present but not taken into account, g‐computation estimator was highly biased under the alternative hypothesis (maximum absolute bias of 1.40 and relative bias of 134%), but no bias was found under the null hypothesis (absolute bias < 0.01 and relative bias < 1% under all scenarios). IPW was also biased, both under the null hypothesis (maximum absolute bias of 0.92) and the alternative hypothesis (maximum absolute bias of 1.21 and relative bias of 116%). Interestingly, even when not corrected for intra‐patient correlation, LTMLE showed a smaller bias than g‐computation and IPW under the null hypothesis (maximum absolute bias of 0.20) and the alternative hypothesis (maximum absolute bias of 0.11 and relative bias of 13%). In all cases where the estimator was biased, the bias increased with the number of time points T and with the strength of the intra‐patient correlation σ U Y and σ U E .

G‐Methods Adapted to Account for Intra‐Patient Correlation

For g‐computation, both the fixed‐effects and mixed‐effects approaches effectively corrected the systematic bias under the alternative hypothesis (absolute bias < 0.01 and relative bias < 1% under all scenarios), and remained unbiased under the null hypothesis. For IPW, the fixed‐effects approaches largely reduced the bias under the null and alternative hypotheses (maximum absolute bias of 0.03 and relative bias of 3%), except for T  = 20 and σ U Y , σ U E  = 0.50 where the bias remained high for fixed‐effects IPW (maximum absolute bias of 0.17 and relative bias of 16%), in relation with positivity violation issues, as indicated by high number of weight values requiring truncation and the increase in standard‐errors. The mixed‐effects approach with IPW was not effective in reducing the bias caused by intra‐patient correlation (maximum absolute bias of 1.14 and relative bias of 109%). This was explained by issues in the estimation of the random effect, as shown by a reduced bias when using the true U E rather than its estimate ( data not shown ). For LTMLE, the mixed‐effects approach effectively reduced the bias under all conditions (maximum absolute bias of 0.05 and relative bias of 5%), whereas fixed‐effects approach only reduced the bias for T  = 5 and 10 time points (maximum absolute bias of 0.01 and relative bias of 1%) but not for T  = 20 time points (maximum absolute bias of 0.43 and relative bias of 42%).

Standard Errors and Coverage Probability

Overall, g‐computation provided estimators with smaller standard errors than LTMLE and IPW under all scenarios, and LTMLE estimators were intermediate, with lower standard errors than IPW. When adapted to the presence of intra‐patient correlation, the mixed‐effects approach provided estimates with lower standard errors than the fixed‐effects approach for g‐computation (maximum standard errors of 0.07 for mixed‐effects and 0.18 for fixed‐effects) and LTMLE (maximum standard errors of 0.27 for mixed‐effects and 0.93 for fixed‐effects). IPW estimators showed high standard errors in scenarios with many time points, linked to problems of positivity violation and weight stability.

The coverage probabilities obtained by bootstrapping with T  = 10 and in the presence of latent variables ( σ U Y  = 0.50 and σ U E  = 0.50) were close to 95% for the fixed‐effects and the mixed‐effects approaches in g‐computation, IPW and LTMLE.

Real‐Life Data Analysis

We applied the previously developed methods to real‐world observational data to investigate the causal effect of pregnancy as the exposure on neurological disability as the outcome in MS patients.

Data Extraction and Inclusion Criteria

Data for all female patients from 1990 to 2020 were extracted from the OFSEP registry. All patients included in the OFSEP registry sign an informed consent for their medical data collected to be used for research purposes after anonymization and aggregation.

Inclusion criteria were a follow‐up starting between the ages of 18 and 45 years and lasting at least 10 years, with at least one annual EDSS measure and with an interval between disease onset and the start of follow‐up of less than 10 years. Patients with missing data for the number of children or their birthdate were excluded. EDSS measures performed during a relapse were not retained, and if more than one EDSS measure was made in the same year, only the first was included. The pregnancy onset was determined by subtracting 9 months from the reported child's birth date.

Statistical Analysis

The causal estimand to assess the effect of pregnancy on the course of disability was the causal mean difference in EDSS at 10 years between the observed situation and the counterfactual situation without pregnancy in the entire study population. Data were split into 10 one‐year periods and analyzed using the methods described above, g‐computation, IPW, and LTMLE, without correction for intra‐patient correlation and with correction for intra‐patient correlation using a fixed‐effects approach and a mixed‐effects approach. All models were adjusted for age at the start of follow‐up, year of the start of follow‐up, and time from disease onset to the start of follow‐up as baseline covariates, and for disease form (RR, SP, or PP) as an exogenous time‐dependent covariate. Confidence intervals were obtained by bootstrapping on 1000 replicates. The corresponding DAG is presented in Figure  S3 , and the description of the models used in g‐methods is in Table  S3 .

Results

Among the 21 744 patients aged between 18 and 45 years in the OFSEP database, 3080 had missing data regarding their child birthdate and 17 989 had less than 10 years of follow‐up or missing EDSS measures; 675 patients with complete data were included, of whom 160 (23.7%) had at least one pregnancy during their follow‐up. The baseline characteristics of the cohort are described in Table  1 . After adjustment for age, patients with a higher EDSS score were less likely to be pregnant (odds ratio per 1 point increase in EDSS score: 0.75 [95% CI, 0.68 to 0.83], p  < 0.001), indicating the existence of a reverse causation bias in the relationship between pregnancy and disability. In a crude univariate analysis, the mean EDSS score at 10 years was lower in patients who had at least one pregnancy during follow‐up (mean EDSS difference: ‐0.8 [95% CI, −1.2 to −0.4], p  < 0.001). An intra‐patient correlation was found in the mixed‐effects models: the standard deviation of σ U Y was 0.46 in the outcome model and of σ U E was 0.56 in the exposure model.

TABLE 1 Baseline characteristics of included patients from the OFSEP database. ≥ 1 pregnancy during the follow‐up No Yes Total
n  = 515
n  = 160
n  = 675 Age at MS onset (in years) Mean (SD) 31.82 (6.52) 25.23 (4.23) 30.26 (6.67) Disease duration at the start of the follow‐up (in years) Median [IQR] 2.04 [0.57; 6.09] 1.20 [0.32; 3.94] 1.91 [0.47; 5.30] EDSS score at baseline Mean (SD) 1.86 (1.53) 1.10 (1.27) 1.68 (1.51) MS phenotype at baseline
n (%) RR 460 (89.3%) 156 (97.5%) 616 (91.3%) SP 28 (5.4%) 2 (1.2%) 30 (4.4%) PP 27 (5.2%) 2 (1.2%) 29 (4.3%) Number of children at baseline Median [IQR] 1 [0; 2] 0 [0; 1] 1 [0; 2] Number of pregnancies during the follow‐up
n (%) 0 515 (100.0%) 0 (0.0%) 515 (76.3%) 1 0 (0.0%) 115 (71.9%) 115 (17.0%) 2 0 (0.0%) 43 (26.9%) 43 (6.4%) 3 0 (0.0%) 2 (1.2%) 2 (0.3%) Abbreviations: EDSS, expanded disability status scale; MS, multiple sclerosis.

Baseline characteristics of included patients from the OFSEP database.

Abbreviations: EDSS, expanded disability status scale; MS, multiple sclerosis.

The estimates of the causal mean difference in the EDSS at the last time point are presented in Figure  4 . Without correction for latent variables, g‐computation concluded that there was a small but significant beneficial effect of pregnancy on EDSS score (mean EDSS score difference: ‐0.04 [95% CI, −0.06 to −0.02]), whereas IPW and LTMLE did not reject the null hypothesis. Accounting for latent variables in g‐computation reduced the estimate of the causal effect of pregnancy and did no longer reject the null hypothesis. For IPW, the fixed‐effects approach provided a similar result, but the mixed‐effects approach resulted in a slightly different point estimate, though it did not reject the null hypothesis either. For LTMLE, all estimates provided point estimates close to zero, with or without correction for latent variables. G‐computation had the smallest standard deviations, followed by LTMLE and IPW, and the mixed‐effects approach provided estimators with lower standard deviations than the fixed‐effects approach in g‐computation and LTMLE.

Estimation of the causal effect of pregnancy on disability (causal mean difference in EDSS at 10 years with 95% confidence interval) on real‐life data. EDSS, expanded disability status scale; IPW, inverse probability weighting; LTMLE, longitudinal targeted maximum likelihood estimator.

Discussion

The presence of an intra‐patient correlation between repeated measures of the outcome and the exposure is plausible for most variables in longitudinal studies. In this setting, even in the absence of unmeasured confounding between the exposure and the outcome, the causal DAG analysis shows that this intra‐patient correlation could open up non‐causal pathways and theoretically bias the estimation of the causal effect of the exposure on the outcome. Herein, we demonstrated in the simulation study that g‐computation, IPW, and LTMLE were indeed biased under certain conditions when intra‐patient correlation was not accounted for. The absolute bias was of variable magnitude, according to the g‐method estimator, the strength of the intra‐patient correlation, the number of time points, and the effect of the exposure on the outcome. Interestingly, LTMLE was less affected by intra‐patient correlation than g‐computation and IPW, with smaller bias in our simulation scenarios.

We proposed and evaluated two ways of blocking all non‐causal pathways related to intra‐patient correlation: a fixed‐effects approach and a mixed‐effects approach. The mixed‐effects approach, by estimating intra‐patient correlation related to latent variables through a random effect in a single outcome model including all time points, has the advantage of being more parsimonious than the fixed‐effects approach, which requires fitting separate models for each time point. However, the mixed‐effects approach is based on stronger assumptions, notably that the latent variables are exogenous and time‐independent, and may require a large number of patients and time points to be estimated properly, particularly for binary variables. In the simulation study, both approaches effectively corrected the bias caused by intra‐patient correlation for g‐computation and LTMLE, with lower standard errors for the mixed‐effects approach. Other conditions are likely to impact the performance of fixed‐effects and mixed‐effects approaches but were not evaluated herein, such as the number of covariates to adjust for and their modeling complexity, or the presence of missing data. Particularly, censoring can lead to a significant decline in available data over time in longitudinal settings [ 27 ], in which case the mixed‐effects approach may be advantageous by providing a better estimate of missing values using the patient information available at other time points. In addition, the present study took place in a discrete‐time setting with regularly spaced visits, in which both approaches could perform well. However, a mixed‐effects approach could be more adapted in near‐continuous contexts, with a huge number of time points and numerous missing data, or with irregularly spaced visits, but this was not investigated here [ 28 , 29 , 30 ]. Lastly, another advantage of the mixed‐effects approach is that it can yield marginal or patient‐conditional effects, whereas the fixed‐effects approach can only provide marginal effects.

IPW presented issues that led it to be an inefficient estimator in two cases of the simulation study. On the one hand, the fixed‐effects approach applied to IPW presented weight stability issues linked to positivity violation when the number of time points or the intra‐patient correlation increased, despite the truncation of extreme weights. The higher risk of positivity violation in longitudinal settings with a large number of time points has already been discussed in the literature [ 23 , 31 ]. The number of possible exposure histories increases exponentially with the total number of time points T ( 2 T for a binary exposure), resulting in higher probability of groups of exposure history with small or null populations. In addition, a stronger intra‐patient correlation in exposure will also logically impact the positivity assumption: patients will have more distinct exposure histories, with some having more exposure events during their follow‐up, and others being more likely to be never exposed. This positivity violation led to an inflation of the variance (partially converted into bias by the truncation procedure) that might render the IPW estimator unusable in certain conditions. On the other hand, the mixed‐effects approach was also ineffective for IPW, but this was probably related to the difficulty of estimating an intra‐patient correlation for a binary exposure variable in our case. Properly estimating random effects in a logistic regression would require more time points, more patients, or higher probabilities of exposure (particularly in the motivating example, where there was a low probability of exposure at each time point). Given these limitations, we reiterate the importance of assessing weight stability when using IPW in longitudinal settings with numerous time points, and recommend particular caution when estimating random effects with a binary covariate.

The application to real‐life data herein demonstrates the necessity of accounting for intra‐patient correlation when studying longitudinal data with repeated measures, and the ability of fixed‐effects and mixed‐effects g‐methods to account for it. For g‐computation, the results were similar to those found under the alternative hypothesis in the simulation analysis, but the null hypothesis was not rejected after correction. We hypothesize that this result might be a manifestation of the g‐null paradox, whose magnitude was amplified by the bias related to intra‐patient correlation. It is important to note that the real‐life data analysis herein served as a demonstrative example. A significant proportion of patients were excluded due to missing data or insufficient follow‐up, whereas an application in an epidemiological study should include them all and use corrective approaches to deal with missing data and censoring. Also note that, as the exact conception date is incompletely recorded in the OFSEP database, it was assessed using childbirth's date, thus omitting pregnancies ending in stillbirth or miscarriage. A real‐life data analysis including a larger number of patients with incomplete follow‐up and missing data has been published in a clinical neurology journal with similar conclusions [ 32 ].

In the study reported by Shardell and Ferruci [ 9 ], the authors used joint models to address the issue of an unmeasured confounder influencing simultaneously the outcome and the exposure at each time point. Their results were, overall, similar to those presented herein: mixed‐effects g‐computation was efficient when the model was well specified, and LTMLE was able to partially account for a misspecification of the outcome model when the exposure model was well specified [ 9 ]. Although our problem can be encompassed by their work, the key point underlying the use of mixed‐effects models differs between the two studies: whereas Shardell and Ferruci [ 9 ] proposed mixed‐effects in joint models as a possible solution to unmeasured confounding, we highlight herein the importance of considering intra‐patient correlation in the causal reasoning as long as one is working on longitudinal data (even if all confounding factors between the exposure and the outcome are taken into account). In other terms, they proposed a tool that can be used in the presence of unmeasured confounding, while we demonstrated the importance of accounting for intra‐patient correlation as a prerequisite as soon as one is dealing with repeated measures. Other published studies investigated causal inference in longitudinal settings with repeated measures: however, they did not include intra‐patient correlation in the causal DAG (and thus in the exposure or outcome model), but took it into account by working at a marginal level, in the structural model, using generalized estimating equations [ 18 , 33 , 34 ], or joint modeling [ 35 ].

Conclusion

In causal inference with time‐dependent confounding between a repeatedly measured exposure and outcome, intra‐patient correlation due to latent variables must be taken into account when constructing and analyzing the causal DAG. The mixed‐effects approach is an interesting alternative to the more commonly used fixed‐effect approach, in particular when the number of time points is high and when the variables of interest are continuous.

Conflicts of Interest

S.V. received lecturing fees, travel grants, and research support from Biogen, BMS‐Celgene, Janssen, Merck, Novartis, Roche, Sandoz, and Sanofi‐Genzyme.

Supporting information

Data S1: Supporting Information.

Data S2: Supporting Information.

Acknowledgments

This study was conducted using data from the OFSEP, which is supported by a grant provided by the French State and handled by the Agence Nationale de la Recherche , within the framework of the France 2030 program, under the reference ANR‐10‐COHO‐002, by the Eugène Devic EDMUS Foundation against multiple sclerosis and by the ARSEP Foundation.

The authors thank Philip Robinson ( Direction de la Recherche en Santé, Hospices Civils de Lyon ) for his help in manuscript preparation.

Example of R Code for g‐Computation, IPW, and LTMLE for 3 Time Points

The code provided below is an example of g‐computation, IPW and LTMLE for 3 time points on simulated data, without correction for intra‐patient correlation, or accounting for intra‐patient correlation with a fixed‐effects approach and a mixed‐effects approach.

Data Simulation

Simulated dataset using the simcausal package with: time point t from 0 to 3, outcome Y and exposure E, outcome and exposure at last time Y.lag1 and E.lag1, outcome at baseline Y0, and cumulative exposure E.cum.

Preview of data in long format

Data Analysis

G‐Computation

Without Correction for Intra‐Patient Correlation

With Correction for Intra‐Patient Correlation  Using  Fixed‐Effects  Approach

Data must be swapped to the wide format

With  Correction  for  Intra‐Patient  Correlation  Using   Mixed‐Effects   Approach

Random  Effect  Might  Not  Be  Identified  in  This  Demonstrative  Example  With  Only  3  Time  Points

IPW

Without Correction for Intra‐Patient Correlation

With  Correction  for  Intra‐Patient  Correlation  Using  Fixed‐Effects  Approach

With  Correction  for  Intra‐Patient  Correlation  Using  Mixed‐Effects  Approach

Random  Effect  Might  Not  Be  Identified  in  This  Demonstrative  Example  With  Only  3  Time  Points

LTMLE

Without Correction for Intra‐Patient Correlation

LTMLE with fixed‐effects and mixed‐effects approaches depends of the IPW (weight.cum.ipw1, 2 or 3) and g‐computation (Y.cf.comp1, 2 or 3 and Qmodel1, 2 or 3) approaches used

Funding: The authors received no specific funding for this work.

Data Availability Statement

The R code for the simulation analysis is provided in Supporting Information . The dataset used for the real‐life analysis can be provided upon reasonable request by any qualified investigator, after approval by the OFSEP.

References

Associated Data

Supplementary Materials

Data S1: Supporting Information.

Data S2: Supporting Information.

Data Availability Statement

The R code for the simulation analysis is provided in Supporting Information . The dataset used for the real‐life analysis can be provided upon reasonable request by any qualified investigator, after approval by the OFSEP.