A comparison of multiple imputation methods for missing data in longitudinal studies

Background Multiple imputation (MI) is now widely used to handle missing data in longitudinal studies. Several MI techniques have been proposed to impute incomplete longitudinal covariates, including standard fully conditional specification (FCS-Standard) and joint multivariate normal imputation (JM-MVN), which treat repeated measurements as distinct variables, and various extensions based on generalized linear mixed models. Although these MI approaches have been implemented in various software packages, there has not been a comprehensive evaluation of the relative performance of these methods in the context of longitudinal data. Method Using both empirical data and a simulation study based on data from the six waves of the Longitudinal Study of Australian Children (N = 4661), we investigated the performance of a wide range of MI methods available in standard software packages for investigating the association between child body mass index (BMI) and quality of life using both a linear regression and a linear mixed-effects model. Results In this paper, we have identified and compared 12 different MI methods for imputing missing data in longitudinal studies. Analysis of simulated data under missing at random (MAR) mechanisms showed that the generally available MI methods provided less biased estimates with better coverage for the linear regression model and around half of these methods performed well for the estimation of regression parameters for a linear mixed model with random intercept. With the observed data, we observed an inverse association between child BMI and quality of life, with available data as well as multiple imputation. Conclusion Both FCS-Standard and JM-MVN performed well for the estimation of regression parameters in both analysis models. More complex methods that explicitly reflect the longitudinal structure for these analysis models may only be needed in specific circumstances such as irregularly spaced data. Electronic supplementary material The online version of this article (10.1186/s12874-018-0615-6) contains supplementary material, which is available to authorized users.


Background
Longitudinal studies, where information on the same participants is obtained repeatedly over time, are frequently used in clinical and population health research. Analysis of data obtained from such studies is often impeded by the presence of missing data due to item or visit non-response and loss to follow-up [1,2].
Commonly used analytic approaches exclude patients or records with missing data, which may lead to biased estimates and considerable loss of precision [3,4].
Multiple imputation (MI) has become a very popular tool for dealing with missing data in recent years [5,6]. MI involves the generation of multiple copies of the dataset in each of which missing values are replaced by imputed values sampled from their posterior predictive distribution given the observed data. Each completed dataset is analysed using the statistical model appropriate to the epidemiological question of interest, and the resulting estimates and standard errors are combined using Rubin's rules [4]. MI methods generally assume data are missing at random (MAR), which requires that the probability of data being missing, conditional on observed data, is independent of the missing data.
Two general approaches for imputing missing data in the presence of multiple incomplete variables are available in standard computer packages [7][8][9]: MI based on the joint posterior distribution of incomplete variables, often referred to as joint modelling (JM), and fully conditional specification (FCS; also known as sequential regression and MI using chained equations (MICE)) [10][11][12][13]. The JM approach commonly assumes that the incomplete variables follow a multivariate normal distribution, often referred to as multivariate normal imputation [12]. FCS, on the other hand, imputes missing values using univariate conditional distributions for each incomplete variable given all the others, cycling iteratively through the univariate imputation models [13].
Both of these MI approaches were originally proposed for cross-sectional data but can be used to impute longitudinal data collected at equal intervals by considering repeated measurements of time-dependent variables as distinct variables [14], denoted as JM-MVN and FCS-Standard, respectively. Although relatively straightforward to implement using existing software, these methods (i) cannot accommodate longitudinal time-dependent covariates that are measured at irregular time intervals and (ii) may experience non-convergence due to model over-fitting and/or collinearity with large numbers of repeated measurements. Several extensions of the standard JM and FCS approaches for imputing cluster/longitudinal data have been proposed in the literature over recent years [14][15][16][17][18][19][20]. The extensions include limiting the number of time-dependent variables in the univariate imputation models within FCS [20]; and specifying imputation models based on the generalized linear mixed-effects model (GLMM) [14][15][16][17][18][19]. The GLMM-based approaches are generally based on more restrictive assumptions about modelling the correlation structure than the JM-MVN and FCS-Standard as these approaches allow for arbitrary dependence of each variable on other variables. The GLMM-based approaches also use the observation at a given time-point as the unit of analysis rather than the individual. It is unclear how important these differences are in practice as currently available comparisons of the various MI models in the literature are limited to a few methods and are in very specific settings [19, [21][22][23][24][25][26] and no comprehensive comparison of the available methods has been conducted.
In the current paper we present a comprehensive simulation-based comparison of the MI methods available in standard software packages for imputation of incomplete longitudinal data. Specifically, we evaluated estimators of regression coefficients for both a linear regression model and a linear mixed-effects model (LMM) in the presence of incomplete binary and continuous predictors. Our primary aim was to investigate whether the GLMM-based MI approaches, which are specifically designed for imputing longitudinal data, provide more accurate estimates than the standard MI approaches. We based our simulation study on a previously conducted analysis exploring the association between the burden of overweight and quality of life (QoL) in the Longitudinal Study of Australian Children (LSAC) [27].

Methods
Longitudinal study of Australian children (LSAC) data and analysis models LSAC is a nationally representative study that examines the development and wellbeing of Australian children. Following recruitment, data have been collected every 2 years (referred to as waves of data collection) using face-to-face interviews, questionnaires and direct anthropometric measurements. The study is ongoing with six waves of data currently available. The detailed study procedure has been described elsewhere [28]. The simulation study was based on data from the kindergarten (K) cohort of LSAC (n = 4983), who were aged 4-5 years when recruited in 2004. The specific question of interest that we focused on was motivated by previous research examining whether the number of overweight occasions in waves 1-5 of data collection predicts poor QoL at wave 6, adjusted for child age, sex, English speaking background and family socio-economic position [27]. In this analysis, the child's age-and-gender specific body mass index (BMI z-score) at each wave was dichotomized into normal weight (including underweight) and overweight (including obesity) using International Obesity Task Force criteria [29].
In our study we slightly modified the original analysis model to include an additional known predictor of child QoL: family structure (which also had missing data), reflecting whether the child was living with a single parent or two parents [30]. Specifically, we were interested in exploring: 1. the association between the cumulative burden of overweight from 4 to 5 to 12-13 years (waves [1][2][3][4][5] and QoL z-scores at age 14-15 years (wave 6), and 2. the cross-sectional association between QoL z-score and BMI z-score across all 6 waves, as these two analysis models are of interest to many researchers in longitudinal data settings. Both analyses were adjusted for sex, English language background and socioeconomic position (all at baseline), age and family structure. The association between QoL z-score at wave 6 and cumulative burden of overweight (OverWtCat) were adjusted for child sex (Sex), English language background (language) and socio-economic position (SEP), all at baseline, age at wave 6 and burden of family structure(FamSt-Cat). OverWtCat and FamStCat both are derived variables representing cumulative burden of overweight and living with a single parent, respectively. These variables were derived by counting the total number of waves, from 1 to 5, on which the child was overweight or lived with a single parent (followed by categorization as none, 1, 2, 3-4, or all five waves), respectively. More precisely, the analysis model to address (1) was a linear regression model: where i = 1, 2, …, n indexes the participant, QoLz 6 is the QoL z-score at wave 6; ageyr 6 represents the age (year) of the study child at wave 6; Sex, SEP and language represent sex, socio-economic position and English language, respectively measured at wave 1. The residual error ε i was assumed to be normally distributed.
The analysis model to address (2) was a LMM applied to all 6 waves of data: where i = 1, 2, …, n and j = 1, 2, …, 6 index participants and waves respectively, BMIz is the BMI z-score, Fam-Struc represents whether the child was living with a single parent or two parents and u 0i represents a subject-specific random intercept. The other variables have the same definition as above. Note that in analysis model (1) only the QoL z-score data from wave 6 is analysed whereas analysis model (2) uses all the QoL z-score data. Both analysis models were susceptible to missing data when applied to LSAC. Data were missing for BMI z-score and QoL z-score in all six waves. Family structure was completely observed for the first wave but had missing values in the subsequent waves. Age and sex were fully observed at baseline, but the other baseline variables had very small amounts of missing information (SEP 0.2% and English-speaking background, 1.4%). Age was occasionally missing at the later waves but was completed by adding the time difference between the wave and wave 1 to the child's age at wave 1. To simplify our example, we excluded cases where QoL z-score was missing across all 6 waves (given that is the outcome of interest) and participants with missing socio-economic position, and/or English language background information at wave 1, leaving a total of 4661 participants for analysis. Except for QoL z-score, most of the missingness observed in our case study was due to dropout (Additional file 1: Table S1).

MI methods to impute longitudinal data
We have identified the following MI methods available for imputing longitudinal data in standard software packages (see Table 1 for additional details).
1. JM-MVN: Originally implemented in Schafer's NORM software [12], this approach utilizes data augmentation, a form of Markov chain Monte Carlo algorithm, to impute missing data under the assumption of (unstructured) multivariate normality. For longitudinal data, this approach treats repeated measurements as distinct variables and imputes all variables in the imputation model as continuous. The model is clearly mis-specified for imputing binary and categorical data, but JM-MVN has been shown to perform reasonably well even for such variables, unless they are severely skewed [  Mplus, R package pan [42].

JM-MLMM-LN
• Repeated measurements of time-dependent variables are imputed using hierarchical models. • All incomplete variables are imputed using a joint multivariate LMM. • Binary and categorical incomplete variables are imputed using latent normal variables. • Can be fitted assuming either a constant or a subject-specific residual error variance.
Fully conditional specification (FCS) (Imputes using a univariate conditional model for each variable with missing data)

FCS-Standard
• Repeated measurements of time-dependent variables are imputed as distinct variables. • Imputes variables using conditional univariate regression models for each incomplete variable, conditional on the time-dependent variables at all waves.

FCS -Twofold
• Repeated measurements of time-dependent variables are imputed as distinct variables. • Imputes variables using univariate regression model for each incomplete variable, conditional on a subset of all time-dependent variables in the data based on a window period. • Imputation carried out in a two-step iterative process.

FCS-MTW
• Repeated measurements of time-dependent variables are imputed as distinct variables. • Imputes variables using univariate regression models for each incomplete variable, conditional on a subset of all time-dependent variables in the data based on a window period. • Imputation carried out in a single step iterative process.

FCS-LMM
• Repeated measurements of time-dependent variables are imputed using hierarchical models.

MI approaches
Method Details Software process is repeated for all incomplete variables in turn. • Binary and categorical incomplete variables are imputed using a latent normal variable. • Can be fitted using either a constant or a subject-specific residual error variance.

FCS-LMM-LN
• Repeated measurements of time-dependent variables are imputed using hierarchical models. • Assumes a conditional LMM for incomplete variables.
• Binary and categorical incomplete variables are imputed using a latent normal variable • Can be fitted using either a constant or a subject-specific residual error variance.
Blimp [45] FCS-LMM-PMM • Repeated measurements of time-dependent variables are imputed using hierarchical models. • Imputes incomplete values using a draw from a pool of observed values who have the closest predicted mean to that of the incomplete case.
R package miceadds [46] The following abbreviations are used to denote different MI methods, e.g., allow the subject-specific error variance to vary across individuals, doing so may not produce stable results in the case of LMM with only random intercepts. We therefore considered homoscedastic error variance for all the methods where possible (all but FCS-LMM-het).

Simulation study Data simulation
To compare the behaviour of the regression coefficient estimators using the above MI approaches, we generated 1000 datasets, each with eight variables (age and sex of the child, English language background, mother's highest education level (whether completed year 12), family structure, socio-economic position of the family, child's BMI z-scores and QoL z scores) and 5000 individuals. In each dataset, variables were simulated in a sequential manner as follows: 1. Child sex, whether English is the main language spoken at home (language) and mother's education (moedu: whether or not completed year 12) at baseline were generated using binomial distributions with probabilities 0.5, 0.9 and 0.6 respectively. 2. Child age (in years) at each wave was generated according to the following model: where 48 + u i. = 48 + N(11,1.5), is the distribution of age (in months) of the participant at the recruitment and v ij = N(0, 2), j = 2,3,...,6, is the random noise in the age distribution to accommodate the random variation of the age at interview during follow-up as observed in the LSAC.
3. The family structure variable (FamStruc, whether the child was living with a single parent or two parents) in each wave was simulated using the following logistic mixed-effects model: where ξ oi = N(0,5.5) is the random intercept.
4. The socio-economic position (SEP) at baseline was generated based on mother education and family structure at baseline (wave = 1) using the following model: where, τ i = N(0,1.2) is the residual error.

5.
The main exposure variable of interest, child BMI z-score was generated based on age, sex and family structure using the following linear mixed-effects model: where φ oi = N(0,0.9) is the random intercept and φ ij = N (0,0.6) is the residual error.
6. Finally, the outcome, QoL z-score at each wave was generated using the following linear mixed-effects model: where ω 0i~N (0,0.7) and ω ij~N (0,0.7). Note that we used a single linear mixed model (7) to generate data for both analysis models (1) and (2). Thus, for analysis model (2), the regression coefficients used in the data generating model were considered as the true values. The true values of the regression coefficients for analysis model (1) were estimated by generating a simulated population of 10 million individuals and fitting the model of interest. All the parameters in the above data generating models were estimated from LSAC data to ensure that the simulated datasets were comparable to a real data example.
As both analysis models were susceptible to missing data when applied to LSAC, upon generating the complete datasets, we set BMI z-scores for all waves and family structure from waves 2 to 6 to missing according to the MAR assumption. Specifically, we used the following equation to induce missingness: i. Values of BMI z-score at j th wave (j = 1, 2,…,6) were set to missing dependent on QoL z-score and age at that wave: ii. Similarly, family structure in waves 2 to 6 was set to missing at a given wave dependent on baseline family structure and the value of the child QoL zscore of that wave: The parameters for the above model (θ 0 , θ 1j , θ 2j and ϑ 0 , ϑ 1j , ϑ 2j ) were chosen based on the LSAC data to ensure a similar proportion of missing observations for each of the variables at each wave as in LSAC. A summary of the proportions of missing data by wave in both the LSAC and simulated data is given in Table 2.
We have also provided parameter values for the above models in the Additional file 1: Table S2, see supplementary materials. Although in the LSAC data and in many longitudinal studies the outcome variable might also have missing values, in our simulation study, we considered the outcome to be fully observed, in order to ensure that the simulated missing data scenarios satisfied the MAR assumption.

Imputation strategies and comparison of methods
We applied each of the 12 imputation models described above to both the simulated and LSAC datasets. For each MI approach, we imputed continuous BMI z-score, binary family structure variable and QoL z-score (only in the case study) at each wave. The imputation model also included age and sex of the child, socio-economic position and whether English was the main language spoken at home at Wave 1, and QoL z-score at each of the 6 waves. For MI methods that impute binary variables as continuous in the imputation model, we rounded the imputed values (to either 0 or 1) using adaptive rounding [31] for analysis model (1), as rounding was required to derive the FamStCat variable (i.e., the number of waves that the participant lived with a single parent). However, for analysis model (2) we used the unrounded values in the analysis as rounding can introduce bias [11,34]. Forty imputations were generated for each approach to limit Monte Carlo (imputation-related) error for the regression coefficient of interest to approximately 5% of its standard error. Of note, for both FCS-MTW and FCS-Twofold method we used measurements within adjacent time-period in the imputation model.
Using simulated datasets, we compared the bias, standard errors (both average of the model based and empirical standard error) and coverage probability of the estimated regression coefficients among the 12 imputation approaches, complete data analysis (fully observed simulated dataset of 5000 observations, before inducing missingness) and available data analysis (that excludes patients/records with missing data). The sampling properties of the estimators are estimated from 1000 simulated datasets.

Simulation results
The sampling distribution of the estimated bias and the coverage of the regression coefficients for analysis model (1) across the 1000 simulated datasets are displayed in Figs. 1 and 2, respectively. A detailed numerical summary of the estimated bias and standard errors is provided in Additional file 1: Figure S1. It is clear from Figs. 1 and 2 that analysis based on the available data resulted in biased estimation of the regression coefficients, and inadequate coverage probabilities. All the MI approaches provided similar estimation of the regression coefficients for fully observed covariates (age at wave 6, sex, socio-economic position and language). The estimated coverage probabilities for these covariates were very  For analysis model 2 (the LMM), we observed mixed results across MI methods in performance associated with incomplete predictors (Figs. 3 and 4). Again, the available data analysis, FCS-Twofold, FCS-MTW and FCS-MLMM-LN resulted in biased estimates of the regression coefficients and under-coverage for both the incomplete predictors. In addition, FCS-GLMM, FCS-LMM-het and FCS-LMM-PMM resulted in biased estimation of regression coefficients and under-coverage for the family structure indicator. All of the MI approaches demonstrated reliable Fig. 1 Distribution of the bias in the estimated regression coefficients (i.e., mean changes in the QoL z-score associated with each covariate) for analysis model (1) across the 1000 simulated datasets following complete data, available data and 12 multiple imputation methods. Top and bottom panel show the distribution of the bias in the estimated regression coefficients for covariates with missing data whereas the middle panel shows the distribution of the bias associated with fully observed covariate estimation of the intra-cluster correlation (Additional file 1: Figure S1).
The implementation time of all of the 12 methods (implemented on a standard Windows PC intel core i-5, 3.2 GHz processor and 8GB RAM) is summarized in Fig. 5

Analysis of LSAC data
Estimated regression coefficients with 95% confidence intervals (CIs) for the two analysis models applied using different imputation approaches to the analysis of the LSAC data are given in Figs. 6 and 7. Both the figures suggest an inverse association between BMI and QoL. For both analysis models, we observed that FCS-Standard and JM-MVN resulted in very similar estimates of the regression coefficients and 95% CIs. The estimated regression coefficients and corresponding 95% CIs for other methods were also generally similar, although some variability among the approaches was observed, especially in some of the results for FCS-LMM-LN and FCS-LMM-het for analysis models (1) and (2), respectively.

Discussion
In this paper we evaluated the performance of currently available MI methods for handling incomplete variables in longitudinal studies in the context of estimating regression coefficients from a linear regression model (with a cumulative measure of exposure) and a LMM, two commonly used models in the analysis of longitudinal data. We found that FCS-Standard and JM-MVN provided reliable estimates of the regression coefficients for both analysis models, often with better coverage  The FCS-Twofold and FCS-MTW methods produced slightly more biased and less precise estimates than FCS-Standard. The observed bias was likely to be because these approaches limit the variables in the univariate imputation models, thereby potentially throwing away important information about the missing data. Similar bias in the estimated regression coefficients following FCS-Twofold was observed by De Silva and colleagues [22], who showed that this bias can be reduced by relaxing the restriction to include a wider time window. In our simulations, we observed very similar results from FCS-Twofold and FCS-MTW, casting doubt on the need for the nested within-and among-time-point iterations in the FCS-Twofold approach.
We observed relatively poor performance for several GLMM-based approaches, namely FCS-LMM-het, FCS-GLMM, FCS-MLMM-LN and FCS-LMM-PMM. Bias in the estimated regression coefficients for binary covariates Fig. 3 Distribution of the bias in the estimated regression coefficients (i.e., mean changes in the QoL z-score associated with each covariate) for analysis model (2) across the 1000 simulated datasets following complete data, available data and 12 multiple imputation methods. Top, left and bottom right panels show the distribution of the bias in the estimated regression coefficients for covariates with missing data and all other panels show the distribution of the bias associated with fully observed covariate following FCS-LMM-het was also reported in a previous simulation study [35]. The poor-performance of FCS-MLMM-LN might be due to its implementation in the FCS framework which requires many iterations to converge. Moreover, this method is computationally demanding, limiting its usefulness in practice. Finally, FCS-LMM-PMM appeared to perform relatively poorly in our simulation study, especially for binary incomplete variables. It would be important to explore the use of this approach in further datasets as it has been shown that predictive mean matching provides good results in the context of cross-sectional data [36].
In general, our results suggest that using either FCS-Standard or JM-MVN produces reliable estimates of the regression coefficients for incomplete longitudinal data. These methods are known to be equivalent in the case of normally distributed covariates [23,37], but for incomplete data in both binary and continuous covariates, the FCS-Standard approach can be more robust than a mis-specified joint modelling approach [38]. This is because for the joint incomplete binary and continuous data the FCS-Standard approach produces imputed values that are compatible with a restricted general location joint model. But for many situations, FCS-Standard may not correspond to a joint model when the conditional models are mis-specified, potentially resulting in sub-optimal imputation [39]. A detailed study of compatibility in the context of longitudinal data is beyond the scope of the present paper but warrants future study.
Both the JM-MVN and FCS-Standard approaches require study designs with a fixed number of sampling waves and fixed time intervals between successive waves. Therefore, these methods may not be appropriate if data are collected in irregular time intervals, and GLMM based approaches should be used.  The results presented here are consistent with the results of a number of previous simulation studies. Comparable results between JM-MLMM and FCS-LMM were also obtained by Zhao and Yucel [32]. Our findings are also consistent with those of Kalaycioglu et al. (2013), who showed that JM-MLMM produced better results, both in Fig. 7 Estimated regression coefficients with 95% CI for analysis model (2) applying available data and all the MI approaches to handle missing data in LSAC terms of bias and precision, than FCS-MTW and FCS-LMM-het, with FCS-LMM-het performing the worst [35]. Similarly consistent results were reported by Audigier et al. (2018), who compared FCS-LMM-het, FCS-GLMM, JM-MLMM-LN methods for imputing incomplete binary and continuous data in the context of individual patient data meta-analysis, and found that JM-MLMM-LN performed better than FCS-LMM-het and FCS-GLMM [26]. Audigier et al. also considered an additional "FCS-2stage" method that fits separate imputation models within each cluster and then combines the results using a multivariate random effects model. They reported that FCS-2stage performed worse than FCS-GLMM and FCS-LMM-het when the cluster size was small. Given that in our study, and often in longitudinal studies, we only had a few waves of data collection, we did not include this approach in our comparison. Our results are also consistent with those of Enders et al. (2016) who showed that JM-MLMM-LN provided better results than FCS-LMM-het for imputing incomplete binary data in the context of a LMM with a random intercept [21].
Although our simulation study was based on a real dataset and thus has a realistic level of complexity, it is always difficult to draw general recommendations from a single simulation study. Further extensions could be to explore the behaviours of these MI methods when data are collected at irregular time intervals, when neither JM-MVN nor the FCS-Standard will be feasible. In our simulation study we only considered covariates to be missing and the outcome to be fully observed. This is because we restricted our simulation study to MAR scenarios, in which missingness in the covariate depends only on the outcome of interest and other (fully observed) covariates in the model. Although it may be reasonable to include a case where both outcome and covariates are missing this may introduce additional complexity with the data being missing not at random (MNAR). Exploring the sensitivity of the performance of the MI methods under various missingness assumptions is beyond the scope of the present paper. Both our simulations and case study are based on simplistic imputation models with few variables. However, in many situations the multiple imputation method may use high-dimensional data with a large number of predictors, in such situations JM-MVN and FCS-Standard may incur convergence problems. To overcome convergence issues, several methods, including using principal component analysis to reduce the dimensionality of the predictors in the imputation model, have been proposed [40,41]. However, a full exploration of methods proposed to address the complexities arising in high-dimensional data, e.g., appropriate selection of predictors and convergence issues is beyond the scope of the present paper. Despite the use of a simplistic imputation model, this study has for the first time provided an overview and a systematic comparison of a growing number of approaches to MI with longitudinal data; as methods continue to develop, further evaluations will undoubtedly be needed.

Conclusion
In summary, both the cross-sectional MI methods (FCS-Standard and JM-MVN) and some GLMM-based approaches (JM-MLMM, JM-MLMM-LN, FCS-LMM and FCS-LMM-LN) performed well for the estimation of regression parameters in the case of a linear regression model and LMM. Both approaches have important strengths and limitations. No single method is appropriate for every situation. As the GLMM-based approaches are still developing and are generally more complex than the cross-sectional methods, these may only be needed in specific circumstances such as irregularly spaced data or high-dimensional data that create convergence problems.

Additional file
Additional file 1: Table S1 and S2 contains values of the average biases, average of the model standard errors and empirical standard errors under analysis model (1) and (2), respectively over 1000 simulated datasets. Figure S1 shows the distribution of estimated intra-cluster correlation coefficients for analysis model (2)