Latent growth curve modeling (LGCM) is widely used in psychological sciences to understand changes over time. It uses latent constructs to estimate the average coefficients of trajectories, e.g., latent intercept and slopes, accounting for individual variations in these coefficients. Although the multilevel modeling (MLM) approach also investigates the growth trajectory of repeated-measures outcomes of interests, LGCM has some advantages over it. For example, LGCM is a structural equation modeling-based approach, and therefore, it explicitly accounts for measurement errors to test the theoretical models. In addition, LGCM enables researchers to investigate relationship among latent growth components (i.e., latent intercept and slopes) and other variables (e.g., covariates and distal outcomes).

Since missing data has been an inevitable issue in longitudinal studies, researchers who employ LGCM must select appropriate methodological approaches to address this challenge. Typical issues caused by missing data are biased parameter estimates, distorted model fit, and power loss (Enders, 2022). Traditional missing data handling approaches, such as listwise deletion and single imputations, have been considered ineffective and problematic. The state-of-arts missing data techniques, such as full information maximum likelihood (FIML) and multiple imputation (MI), have gained increasing popularity in recent decades. FIML is wildly recognized in the fields of latent variable analysis, including LGCM. Previous research revealed when the multivariate normality assumption holds, FIML could produce unbiased parameter estimates and appropriate test statistics, taking into account missing data simultaneously (Collins et al., 2001; Enders and Bandalos, 2001; Rubin, 1987). In contrast, MI is a family of approaches that addresses missing data and perform data analysis in separate steps. While FIML is often the default approach for handling missing data in the LGCM, MI offers several important advantages in many applied research settings. For example, MI can readily accommodate a large number of auxiliary variables that are not part of the analysis model but improve missing data recovery, whereas incorporating numerous auxiliary variables with FIML can be cumbersome and computationally challenging (Graham, 2003; Savalei and Bentler, 2009). In addition, with complex models or complicated data structures, FIML may fail to converge or become computationally infeasible (Savalei and Bentler, 2009). In contrary, MI can often overcome these limitations because the estimation is conducted on completed data sets following imputation (Enders and Mansolf, 2018; Gottschall et al., 2012). MI also offers greater flexibility in specifying imputation models, allowing it to better accommodate nonnormal or nonlinear relationships (Asparouhov and Muthén, 2010; White et al., 2011; Jia and Wu, 2022). This article therefore focuses on conditions under which both FIML and MI are practically applicable and directly comparable, in order to examine their relative performance in a controlled setting. Although MI’s primary advantages often arise in scenarios where FIML is difficult or impossible to apply, establishing their comparative behavior under shared feasible conditions provides an important foundation for understanding their respective strengths and limitations.

Despite its advantages, MI is less straightforward than FIML. It generally requires three steps: imputation (i.e., generating multiple complete data sets with missing values filled in), analysis (i.e., fitting multiple imputed data sets into the target analysis model), and pooling (i.e., combining the analysis results from these multiple imputed data sets). To ensure unbiased parameter estimates, each step requires to implement appropriate approaches. Generally speaking, the imputation model should align with the research context, the nature of the missing data, and the underlying assumptions of missing data mechanism (we will describe it in the next section). For longitudinal data, the process for selecting an appropriate imputation strategy can be even more intricate. Specific to longitudinal studies, data can be organized in either a wide format or long format. In the wide format, observations across measurement occasions are presented in separate variables. In such data structure, each row stores all observational data for one individual but with a greater number of columns. In the LGCM framework, the most common practice is to perform both imputation and a subsequent data analysis in the wide format.

The long format involves vertically stacking the observations from multiple measurement occasions for each individual, with an additional variable indicating the measurement occasion. This data format is commonly used within the MLM framework. For data organized in long format, several imputation methods suitable for the multilevel data structure have been proposed, especially for cross-sectional analysis (e.g., students clustered within classroom; Audigier and Resche-Rigon, 2018; Carpenter et al., 2011; Enders et al., 2020; Quartagno and Carpenter, 2016, 2019; Resche-Rigon and White, 2018; van Buuren and Groothuis-Oudshoorn, 2011; Zhao and Schafer, 2015). Conceptually, the multilevel imputation approach is appropriate for the longitudinal data for a LGCM model, i.e., repeated measures inherently nested within individuals. However, limited research has examined their performance with such longitudinal data structures. One complication is that such investigation requires data transformation between wide and long formats. Specifically, one needs to convert original data to long format for the multilevel imputation process, then revert imputed data back to wide format to accommodate the requirements of the latent growth analysis model. Given the growing awareness in multilevel MI methods for longitudinal designs and recent software advancements in the MLM framework (e.g., Enders et al., 2020), it is timely and relevant to evaluate how these conceptually appropriate multilevel MI methods (a.k.a., long-format MI) perform when applied to repeated-measures data in LGCM, compared to standard, single-level wide-format MI.

Another recent methodological advancement in MI further motivated this study. While parameter estimates and standard errors can be directly pooled using Rubin’s (1987) rules, strategies for pooling likelihood ratio test (LRT) statistics and practical fit indices have received comparatively little attention. Enders (2022) summarizes two strategies for pooling LRT: D₂ (Li et al., 1991) and D₃ (Meng and Rubin, 1992), which have been evaluated in a few simulation studies (e.g., Enders and Mansolf, 2018; Liu and Sriutaisuk, 2020) under different scenarios. In more recent work, Chan and Meng (2022) proposed a new LRT statistic pooling strategy, which has been denoted as D₄ and proven at least as effective as D₃ while offering greater efficiency (Grund et al., 2023). To our knowledge, none of these pooling approaches has been evaluated within the LGCM framework.

This current article aims to provide researchers with guidance on the best practice in handling missing data using MI when performing a LGCM, including recommendations for selecting appropriate imputation methods (e.g., wide-format vs. long-format) and optimal pooling strategies for model fit evaluation. The remaining of the article is organized as follows. We first introduce the basics and principals of LGCM, and then discuss missing data mechanisms and highlight the difference among commonly used missing data techniques, with a focus on MI. We then delve into the imputation phase in MI, introducing imputation algorithms and imputation models that are suitable for scenarios in the LGCM framework. Additionally, the pooling phase is reviewed, with an emphasis on the pooling strategies for the LRT statistic. Moving forward, we present a simulation study that compared the combinations of these imputation and pooling strategies with FIML. Finally, we provided an empirical example and conclude with practical guidance on implementing these strategies in MI in the context of LGCM.

A classic application of LGCM is to analyze continuous repeated-measures data collected from a group of individuals and estimate intraindividual (within-person) growth over time [i.e., latent intercept and latent slope(s)], while accounting for variability at interindividual (between-person) levels (Kohli and Harring, 2013; Bollen and Curran, 2006). A latent intercept captures the initial level; and latent slope(s) usually represent a growth pattern as a polynomial function of time, such as linear, quadratic, cubic, etc. These two types of latent growth factors are measured by their mean levels, individual variances, and covariances between them (Depaoli et al., 2023). The repeated measures are indicators of latent growth factors, linked by factor loadings. A typical latent quadratic model, for example, can be presented in matrix form in Equation 1.

where t j is time code for the j^th measurement occasion; Y j represents a vector of outcome for all individuals at the j^th measurement occasion; and ε j is a residual vector at the j^th measurement occasion (j = 1, 2, … T). The latent intercept, latent linear slope, and latent quadratic slope are denoted by η I , η S 1 , and η S 2 , respectively. For a linear model, the model is a reduced form of Equation 1, with the Λ y matrix only containing the first two columns, and η matrix only containing η I and η S 1 .

The model fit evaluation is commonly based on the likelihood ratio test (LRT) statistics, which assesses the difference between the model implied and observed covariance matrices and mean vectors. The LRT statistic, denoted T, follows a χ 2 distribution assuming data are normally distributed and can be written as follows.

where l ( θ ̂ ) is the maximized log-likelihood under the hypothesized model, and l ( β ˜ ) is the corresponding log-likelihood of the saturated model (Bollen, 1989).

Data can be missing during the data collection process through various ways. Rubin (1976) summarized three primary missing data mechanisms, including missing completely at random (MCAR), missing at random (MAR), and missing not at random (MNAR). MAR occurs when probability of data being missing on a variable is only related to the observed variables. MCAR is sometimes considered as a special case of MAR; the former refers to the situation where the probability of data being missing on a variable is unrelated to its values or any other observed variables in the dataset. These two mechanisms are commonly referred to as “ignorable missingness” when the condition of “distinct parameters” holds. That means, the parameters governing the distribution of the data are distinct from those governing the missing data process. In an LGCM, this implies that the joint likelihood can be separated into two parts: a part for the LGCM and a part for the missingness mechanism. Because the missingness part does not contain information about the growth parameters, it does not need to be modeled in order to obtain valid estimates of the growth parameters. In contrast, MNAR, represents a “non-ignorable” missingness scenario where probability of data being missing on a variable is determined by unobserved values themselves, demanding more intricate methodological approaches for an effective solution. A common type of missingness in longitudinal studies, called attrition, happens when participants drop out from a study after initially contributing data in the early measurement occasions. Attrition at later occasions is likely to be informed by the data available from the preceding occasions (Allison, 2001), in which case, the missing data due to attrition can reasonably be assumed to be MAR.

In situations where an “ignorable” missing data mechanism is present (i.e., MCAR or MAR), modern missing data techniques, such as full information maximum likelihood (FIML) or multiple imputation (MI), can effectively recover the missing information. FIML is an approach that estimates the parameters of a statistical model and handle missing data simultaneously. It is based on an iterative process that aims to find the values maximizing the overall log-likelihood, while accounting for the patterns of missing data (Enders, 2022). MI, in contrast, typically handle missing data through three steps. In MI, missing values are first imputed multiple times based on a certain parametric or nonparametric model, creating multiple complete datasets (the imputation step). These imputed datasets are subsequently used for the target statistical analysis individually (the analysis step). Then, the results of all the multiple imputed datasets are pooled to generate the final result (the pooling step; Rubin, 1987). This multi-step approach could avoid potential convergence problems FIML may have with complex data or models and enables flexibility in incorporating auxiliary variable and handling various types of missing data. However, the implementation of MI could be more challenging than FIML as it requires more complex methodological considerations. For example, researchers need to make decisions on the choices of (1) imputation algorithm, (2) imputation model, and (3) pooling strategies for different types of statistics, especially for the LRT statistic and model fit indices.

In this section, we applied the aforementioned imputation approaches and pooling strategies to a real-world data set for demonstration. The data were obtained from the Early Childhood Longitudinal Studies Program (ECLS-K: 2011 Kindergarten–Fifth Grade; Mulligan et al., 2019; Tourangeau et al., 2019). This national data set contains multisource information that offers a broad and multifaceted view of young children’s early childhood experiences at home and in school. The sample consists of children who either entered kindergarten for the first time or repeated kindergarten during the 2010–11 school year. We selected data collected during the fall semesters from kindergarten through second grade, and the spring semester from second grade through third grade. Among all the outcome measures related to children’s home and school experiences, our primary focus was on their trajectory of mathematics achievement. For demonstration purposes, we randomly selected 600 children who had complete data on their math item response theory (IRT) scores at all measurement points. We then imposed attrition under MAR at all time points except for the first one, following Equation 5, with b₁ set at 0.4, and missing data proportions of five measurement occasions set at (0, 0.12, 0.24. 0.30, 0.36), respectively. At each time point, the probability of a value being missing is impacted by the values observed at the immediately preceding time point. We fit a quadratic latent growth curve model to the data. Then we intended to center time around the middle point and wrote loading matrix ( Λ ) as follows:

where the three columns represent the loadings for intercept, the liner slope and the quadratic slope, respectively. The analysis results are shown in the Table 2. The parameter estimates obtained from all MI methods were comparable to those from FIML.

In the empirical analysis, the LRT statistic evaluates the fit of the LGCM (null hypothesis) against the corresponding saturated model (alternative hypothesis). FIML yields the conventional model LRT statistic, whereas each multiple imputation method pools the LRT statistics generated from regular maximum likelihood across imputations using the D₂, D₃, and D₄ pooling procedures. The D₂ test statistics with MI-wEMB were most similar to that from FIML, while the D₃ and D₄ methods with MI-wEMB and all pooling methods with MI-wFCS yielded lower LRT statistics values. The lowest LRT statistic values were found from all pooling approaches with MI-lqMLM. These results were consistent with the simulation results from Study 2.

The current study evaluated the efficacy of various MI approaches applied to wide- and long-format data in estimating parameters and assessing model fit in linear and quadratic latent growth curve models. Our aim was to offer researchers practical guidance for handling missing data within this specific modeling framework.

Our analysis revealed that, similar to FIML, the MI approaches consistently yielded unbiased parameter estimates and comparable MSEs cross all examined parameters. Specifically, all methods performed well in accurately estimating the means, variances, and covariances of the growth factors. However, divergent patterns emerged in CIC. For the linear growth model, CICs for the latent variance of the linear slope occasionally fell below 0.925 when using MI-wEMB and MI-lqMLM, particularly under higher levels of missingness in MAR conditions. For the quadratic model, larger proportions of missing data also led to overly liberal CICs (i.e., < 0.925) for MI-wEMB and MI-lqMLM, particularly for parameters related to the linear slope.

Model fit evaluation was one of the two central objects of our study. To assess model fit, we examined three pooled test statistics (D₂, D₃, and D₄) obtained from the MI methods, and compared them with the FIML LRT statistic. When the analysis model was correctly specified, the Type I error rates of the pooled test statistics from the wide-format MI methods (MI-wEMB and MI-wFCS) were generally comparable to those produced by FIML. In contrast, the long-format MI methods (MI-llMLM and MI-lqMLM) yielded near-zero Type I error rates, indicating an overly conservative tendency. Under the quadratic model misspecification, where a latent covariance was incorrectly fixed to zero, MI-wEMB and MI-wFCS demonstrated the second- and third-highest power, respectively, to detect misspecification with smaller sample sizes, ranking just below FIML. As expected from the correct-model results, the long-format MI-lqMLM showed the lowest power. Further comparison under model misspecification revealed that the D₂ test statistics from MI-wEMB produced power values comparable to FIML, whereas D₃ and D₄ consistently yielded lower power values.

Overall, unlike the wide-format imputation methods (MI-wEMB and MI-wFCS), the long-format methods (MI-llMLM and MI-lqMLM) performed poorly in model fit evaluation. This shortfall likely stems from the long-format methods’ transformation of data into long format prior to imputation, followed by the specification of a parametric linear (MI-llMLM) or quadratic (MI-lqMLM) mixed-effects model. These impose greater constraints than less parametric wide-format approaches (no transformation needed), which permits an unstructured covariance matrix across time points without assuming specific random-effects distributions (as in MI-wEMB), or imputes missing values in time points sequentially based on Bayesian regression (as in MI-wFCS). Consequently, the constrained imputation model in long-format methods compels imputed values to adhere closely to the assumed trajectory and clustering structure, thereby reducing unexplained variability in the completed datasets more than a less restrictive model would.

This study provides practical guidance for researchers to select the appropriate approach dealing with missing data in the context of the LGCM framework. Given the growing attention to multilevel long-format imputation methods in recent literature (e.g., Quartagno and Carpenter, 2019; Resche-Rigon and White, 2018; Audigier and Resche-Rigon, 2018; Enders et al., 2020; Grund et al., 2018), it is important to explore how these approaches can be effectively extended to the LGCM, particularly in addressing the inherent clustered data structure. Building on our findings, we highlight key insights as follows to support better methodological decision-making and identifying promising directions for future studies.

While our results indicate that FIML performs best under the specific conditions examined, this does not diminish the importance of MI in many real-world applications where FIML cannot be reliably implemented. MI is preferable in several key scenarios: when researchers require complete datasets for secondary analyses, encounter convergence issues with FIML (which are common in repeat measures designs), or need to incorporate auxiliary variables not included in the final analysis model. While MI’s stepwise approach yielded comparable results in most cases and offers greater flexibility, it requires careful specification of the imputation model and thoughtful selection of pooling strategies.

Despite differences observed in parameter estimates and model fit statistics, MI methods generally served as viable alternatives to FIML. All examined MI approaches generated unbiased parameter estimates, with MI-wFCS (implemented via Blimp) standing out for its robust performance, producing accurate CICs across all parameters and demonstrating minimal sensitivity to both the mechanism and proportion of missing data. MI-wEMB exhibited greater power in LRT overall than other MI methods. For model fit evaluation, the D₂ test statistic paired with MI-wEMB most closely approximated FIML’s ability to detect misspecification in latent factor covariances, whereas D₃ and D₄ required larger sample sizes to perform adequately. Notably, as sample size increased, the difference in power among imputation and pooling methods diminished substantially.

The long-format multilevel imputation is particularly appealing for complex LGCMs with a large number of time points or nonlinear trajectories, as it can, in theory, account for within-person dependencies and between-person heterogeneity through the specification of random effects during the imputation phase. However, our findings under the examined conditions indicate that these theoretical advantages, primarily derived from the multilevel modeling framework, did not translate into clear practical benefits over standard wide-format imputation approaches in LGCM. The long-format imputation methods examined in this study showed outstanding performance in parameter recovery relative to their wide-format counterparts, however, they could be more sensitive to missing data proportion when estimating CIs. More critically, the long-format imputation methods associated with the D₂, D₃, or D₄ pooling strategies produced near-zero Type I error rates. This result likely arose because the imputed values were generated deterministically and were overly consistent with the analysis model (e.g., linear or quadratic). Without valid between-imputation variability, these pooling strategies tended to underestimate the true uncertainty associated with the missing data. As a result, the pooled LRT statistics became overly conservative and were unlikely to reject the null hypothesis, even when they should have. This conservativeness not only deflated the Type I error rate but also severely reduced statistical power, limiting the ability to detect model misspecifications when they existed.

It is important to acknowledge the limitations of our findings within the context of the specific LGCMs and conditions examined. First, the conclusions may not generalize to alternative model configurations, such as conditional growth curve models with covariates, multivariate growth curve models, or multilevel latent growth curve models with additional hierarchical levels. Further research is needed to investigate these extensions. Second, the observed performance of MI methods may not hold across different types of model misspecifications. While our study focused on misspecification in the analysis model, misspecification can potentially arise in both the imputation and analysis models. The performance of long-format multilevel imputation under various misspecification scenarios presents an intriguing avenue for future investigation. Third, in our simulation, we focus exclusively on conditions allowing a straightforward head-to-head comparison of FIML and MI (e.g., no auxiliary variables under MCAR and few auxiliary variables under MAR), where both approaches are easily implemented and should be equally efficient in theory. However, as documented by Savalei and Bentler (2009), MI generally provides more straightforward inclusion of a large number of auxiliary variables, while FIML often requires saturated-correlates modeling with auxiliary variables that could substantially increase complexity and compromise convergence. This important practical distinction suggests avenues for future research, such as targeted comparisons of MI and saturated-correlates FIML in LGCMs with many auxiliaries or under conditions that challenge FIML convergence. Fourth, our simulation was conducted under the assumption of data normality, utilizing ML or FIML estimators. However, when normality is violated, methodological adjustments become critical. Nonnormality can introduce significant challenges, including biased standard errors, inaccurate test statistics, and potentially misleading fit indices. An important direction for future research involves exploring how different imputation and pooling strategies address nonnormality in LGCM. Fifth, our evaluation on long-format multilevel imputation was limited to the approaches available in Blimp (Enders et al., 2020; Keller and Enders, 2023) and did not include more recent or alternative substantive-model-compatible (SMC) imputation developments for multilevel data, such as the SMC-JM approach in the R package jomo (Quartagno and Carpenter, 2022; 2023), and the sequential modeling approach in the R package mdmb (Grund et al., 2021). Future research could incorporate direct comparisons of these methods under conditions similar to or more extreme than those examined in our study, to further inform methodological choices for handling complex missing data. Finally, the pooling strategies employed in this study did not appear to be well-suited for long-format multilevel imputation approaches. The observed deflated Type I error rates and reduced statistical power highlight the need for further methodological investigation and the development of more targeted post-imputation techniques.