Using Data-Dependent Priors to Mitigate Small Sample Bias in Latent Growth Models

Abstract

Mixed-effects models (MEMs) and latent growth models (LGMs) are often considered interchangeable save the discipline-specific nomenclature. Software implementations of these models, however, are not interchangeable, particularly with small sample sizes. Restricted maximum likelihood estimation that mitigates small sample bias in MEMs has not been widely developed for LGMs, and fully Bayesian methods, while not dependent on asymptotics, can encounter issues because the choice for the factor covariance matrix prior distribution has substantial influence with small samples. This tutorial discusses differences between LGMs and MEMs and demonstrates how data-dependent priors, an established class of methods that blend frequentist and Bayesian paradigms, can be implemented within Mplus 7.1 to abate the small sample bias that is prevalent with LGM software while keeping additional programming to the bare minimum.

Keywords

latent growth models small samples data-dependent prior second-order growth model latent basis model

In the social and behavioral sciences, research questions concerning the growth of a trait or construct over time are quite common. To model growth over time, mixed-effects models (MEMs; aka hierarchical linear models, random-effects models, linear mixed models, and multilevel models) and latent growth models (LGMs; aka latent curve models) are two widely used methods. MEMs and LGMs share many common features such as estimating fixed effects that capture the growth trajectory across the entire sample and also random effects that allow for subject-specific growth trajectories that vary from person to person. Curran (2003) delves into great detail about the similarity of MEMs and LGMs and provides examples that demonstrate the near identical results provided if both methods are applied to the same data under fairly idealistic conditions (e.g., linear growth and large sample sizes). Despite the similarity of MEMs and LGMs, Curran (2003) notes that “[LGMs] may provide advantages over the current [MEM] under certain empirical conditions” (p. 541). To elaborate, some types of models can only be fit as an LGM or are much more straightforward to fit as an LGM rather than an MEM, especially in software (e.g., latent outcomes, nonlinear growth trajectories; more detail on these differences are provided shortly).

As will be discussed in detail in the subsequent sections of this article, the estimates of the variability of the subject-specific growth trajectories can be quite different between MEMs and LGMs when sample sizes are small. Specifically, LGMs most commonly implement full information maximum likelihood (FIML) estimation, which is known to yield downwardly biased estimates of the variance components of the factor covariance matrix with small sample sizes (e.g., Browne & Draper, 2006; Cheung, 2013). On the other hand, MEMs most commonly implement restricted maximum likelihood (REML) estimation, which has been shown to greatly reduce bias in small sample size contexts (e.g., McNeish & Stapleton, 2014). However, due to the nature in which LGMs are specified, an REML-like correction is more complex and therefore has not yet been developed broadly for LGMs. Because some types of growth models are impossible or at the very least rather difficult to fit as an MEM, the known small sample bias inherent with FIML and LGM is troublesome. In this article, our goal is to outline one method (and details for implementation in Mplus) in which this bias in the variability of the subject-specific growth trajectories can be reduced.

To outline the remaining sections of this article, first we will briefly overview the differences between FIML and REML. Then, we will compare and contrast MEMs and LGMs and discuss the types of models that may be difficult to fit as an MEM. We will overview the software available to fit these models and also introduce Bayesian estimation of LGMs with a focus on data-dependent priors (DDPs). Illustrative simulations are provided to demonstrate the utility of this approach to reduce bias in the elements of the factor covariance matrix with LGMs and small samples. The procedure for fitting a model with DDPs is then detailed with an applied example.

Brief Comparison of FIML and REML

Consider a set of data that repeatedly measures N individuals over m time points and there is a research interest in subject-specific growth trajectories rather than merely the average growth of the whole sample. For this scenario, an MEM for continuous outcomes can be fit such that,¹

Y = X β + Z b + ∊,

where Y is a vector of responses, X is a design matrix for the fixed effects, β is a vector of fixed effects, Z is a design matrix for the random effects, b is a vector of random effects, E(b) = 0, Cov(b) = G, ∊ is a vector of residuals where E(∊) = 0 and Cov(∊) = R, both b and ∊ are normally distributed, and Cov(b, ∊) = 0. If the model were estimated with FIML, the objective log-likelihood function for the variance components with FIML is:

l_{FIML} (G, R) = - \frac{1}{2} \log | V | - \frac{1}{2} {(Y - X {\hat{β}}_{GLS})}^{T} V^{- 1} (Y - X {\hat{β}}_{GLS}) + c,

where c is a constant that does not rely on model parameters, V = Var(Y) = ZGZ ^T + R and ${\hat{β}}_{GLS} = {(X^{T} V^{- 1} X)}^{- 1} X^{T} V^{- 1} Y$ , the generalized least squares estimator for the fixed effects. Because the estimates of G, R, and β are dependent upon each other, the estimation procedure alternates between estimating the variance components (G and R) and β. That is, the initial estimates of G and R from the first iteration are used to estimate β, which is then used to update the estimates of G and R and so forth until convergence is reached. For this reason, the estimation process is sometimes referred to as ML via iterative generalized least squares (Browne & Draper, 2006; de Leeuw, Meijer, & Goldstein, 2008 p. 40; Goldstein, 1989). In Equation 2, no adjustment is made to account for β being estimated in the model. REML takes the loss of degrees of freedom from estimating β into account. To do this, REML uses the residuals from an ordinary least squares fit, $r = Y - X {\hat{β}}_{OLS} = Y - X {(X^{T} X)}^{- 1} X^{T} Y$ , and maximum likelihood is then applied to these residuals, which are independent of the fixed effects by definition. The objective log-likelihood function for the variance components with REML is thus:

l_{REML} (G, R) = - \frac{1}{2} \log | V | - \frac{1}{2} \log | X^{T} V^{- 1} X | - \frac{1}{2} {(Y - X {\hat{β}}_{GLS})}^{T} V^{- 1} (Y - X {\hat{β}}_{GLS}) + c .

The $\frac{1}{2} \log | X^{T} V^{- 1} X |$ term in Equation 3 accounts for the degrees of freedom lost in estimating β and reduces the bias from estimates of the variance components with small samples. Instead of estimating β and the variance components simultaneously, β is estimated separately from the variance components with REML.

Conceptually, the difference between FIML and REML is akin to dividing a sample variance estimate by n or n − 1. With large samples, the difference is hardly noticeable; however, for decreasing sample sizes, the difference becomes increasingly discernable. Currently, software for MEMs (e.g., SAS Proc Mixed 9.3, HLM 7, R lme4 package) estimate models with REML by default or feature REML as an option (Stata xtmixed and MLwiN do not use REML by default but can implement it).

Brief Introduction to LGMs

The general linear LGM with k time-invariant covariates can be written as a confirmatory factor analysis model with an imposed mean structure such that:

\begin{matrix} Y_{i j} = η_{0 i} + η_{1 i} t_{i j} + ∊_{i j} \\ η_{0 i} = α_{β 0} + γ_{01} X_{1 i} + \dots + γ_{0 k} X_{k i} + ζ_{0 i} \\ η_{1 i} = α_{β 1} + γ_{11} X_{1 i} + \dots + γ_{1 k} X_{k i} + ζ_{1 i}, \end{matrix}

where Y_ij is the response for the ith individual at the jth time, η_0i is the latent intercept for the ith individual, η_1i is the latent slope for the ith individual, $ζ_{i} = {(ζ_{0 i}, ζ_{1 i})}^{T}$ is a vector of random effects for the ith individual, t_ij is the jth time point for the ith individual, and ∊ _ij is the residual for the ith individual at the jth time. In matrix notation, Equation 4 becomes:

Y = Λ η + ∊

and

η = α + Γ X + ζ .

When translating between Equation 4 and Equation 5, readers may note that Equation 4 does not explicitly include Λ whereas Equation 5 does not explicitly include time. The equivalence is not directly explicit because the Λ in Equation 5 contains time as one of its columns in addition to information about other growth parameters. For instance, for the case of linear growth with four equally spaced time points, Λ would equal ${[\begin{matrix} 1 & 1 & 1 & 1 \\ 0 & 1 & 2 & 3 \end{matrix}]}^{T}$ . The first row (before transposing) corresponds to the intercept factor (whose effect is constant across measures), and the second row contains the time information that is captured by t in Equation 4.

The model-implied mean and covariance structures of the repeated measures based on Equations 5 and 6 are thus,

μ = Λ (α + Γ κ)

and

Σ = Λ (Γ Φ Γ^{T} + Ψ) Λ^{T} + Θ,

where μ is a vector of model-implied means of the outcome variables, Λ is a matrix of loadings that can, but are not always, be prespecified to fit a specific type of growth trajectory; α is a vector of latent factor means; Σ is the model-implied covariance matrix of the outcome variables; Ψ is the covariance matrix of the latent factors; Θ is a matrix of residual variances and covariances among the repeated measures; Γ is a matrix of coefficients for the predicted effect of time-invariant covariates on the latent growth trajectory factors; κ is a vector of covariate means; and Φ is a covariance matrix of the covariates X_ki (Bauer & Curran, 2003; Biesanz, Deeb-Sossa, Papadakis, Bollen, & Curran, 2004). Comparing Equation 1 to Equations 5 through 8, α and Γ would be related to β, Ψ is related to G, and Θ is related to R. Λ would be contained within X (often in the first two columns), and the subsequent columns of X would consist of observed values for variables used to compute κ.² For a full comparison of MEMs and LGMs, readers are referred to Curran (2003).

Unlike MEMs, the structural equation modeling framework and the associated software used to fit LGMs (e.g., Mplus, LISREL, and EQS 6.2) does not have an REML analog and instead uses FIML. As a result, researchers who prefer LGMs to MEMs due to familiarity or their increased flexibility (discussed in more detail shortly) may be at a disadvantage when modeling data with smaller sample sizes that are common in behavioral sciences. This tutorial demonstrates a method by which the flexibility of LGMs can be maintained without succumbing to small sample bias incurred through the use of FIML in current software.

MEMs Versus LGMs

Although the difference between MEMs and LGMs is often a matter of disciplinary semantics because models can be formulated interchangeably (Curran, 2003), software implementation of MEMs and LGMs is not always equally interchangeable. Figure 1 shows a conceptual Venn diagram of the software differences between programs for MEMs and LGMs. Although many models can be formulated in either framework (as indicated by the large overlap in the center), there are a handful of nontrivial differences between methods.³ This section will focus in particular on the ability of LGMs to accommodate the latent basis model and to accommodate latent outcomes.

Figure 1.

Venn diagram contrasting software implementations of mixed-effects models (MEMs) and latent growth models (LGMs).

Latent Basis Model

The LGM framework in structural equation modeling programs such as Mplus, LISREL, or EQS provides much greater flexibility for growth modeling compared to analyses that can be done straightforwardly in MEM software. Within software for the MEM framework, time must be prespecified as a variable in the data. However, the LGM framework treats time as a parameter than can be estimated from the data to allow for a flexible and straightforward model for nonlinear growth (commonly referred to as a latent basis model; Grimm, Ram, & Hamagami, 2011; McArdle & Epstein, 1987; Meredith & Tisak, 1990). To demonstrate, using the notation in Equation 1, MEM software require time to be prespecified such that, for instance, $X = {[\begin{matrix} 1 & 1 & 1 & 1 \\ 0 & 1 & 2 & 3 \end{matrix}]}^{T}$ in a hypothetical unconditional, four repeated measure model. In LGM software, prespecified time can similarly be accommodated; however, these elements (the loadings from the slope latent variable to the observed indicators) can also be estimated such that, for instance, $Λ = {[\begin{matrix} 1 & 1 & 1 & 1 \\ 0 & 1 & λ_{22} & λ_{32} \end{matrix}]}^{T}$ where λ₂₂ and λ₃₂ are estimated from the data and the first two loadings are constrained to give the estimates scale. For data collected at four equally spaced intervals, λ₂₂ and λ₃₂ may be estimated to be, for instance, 5 and 7 (rather than constrained to 2 and 3), which would be indicative of nonlinear growth because the estimated loadings are both unequally spaced and larger than would be anticipated, given the time points at which the data were collected.

Latent Variable Outcomes

LGM software can easily incorporate latent outcome variables that have their own particular measurement model (i.e., second-order growth models and curve of factor models), which has particular appeal in behavioral science research where outcomes may not be well-defined or are not represented by a single, observable variable (Curran, 2003; Hancock, Kuo, & Lawrence, 2001; McArdle, 1988; Meredith & Tisak, 1990). That is, rather than directly modeling the change in the observed variables Y_ij , one can model the change in a latent variable η _ij where the first-order measurement model is defined by Y = τ + Λ_ηη + ∊, for τ a vector of intercepts, Λ_η a matrix of loadings, η a matrix of latent constructs, and ∊ a vector of residuals. The second-order model then models the growth over time in the first-order latent variables. MEM software cannot accommodate second-order growth models as a simultaneous system using REML.

Seeing as latent variable outcomes are commonplace in many behavioral sciences and that nonlinear growth is ubiquitous, particularly in areas such as developmental and educational psychology (Grimm et al., 2011), researchers resorting to LGM software to accommodate these popular models may be at a disadvantage with smaller samples because LGM theory and software have yet to advance an REML equivalent and thus models are estimated with FIML. Alternatively, as will be discussed in the next section, a Bayesian analysis using typical “noninformative” prior distributions for the variance components is unlikely to be an adequate solution with small sample sizes because these priors can actually be quite informative under small sample conditions and lead to considerable bias in latent factor variance estimates (Carlin & Louis, 2000; Gelman, 2006).

Brief Discussion of Additional Differences

Although additional differences will not be discussed in detail, a brief description will be provided as will references for readers interested in further detail.

MEMs

Most MEM software have a prepackaged command for modeling the residual covariance matrix (the R matrix in Equation 1) including common structures such as autoregressive, Toeplitz, or compound symmetry. Different programs will have different structures included (e.g., SPSS MIXED has 17 structures, SAS Proc Mixed has 37 structures, and Stata xtmixed has 8 structures; Peugh & Enders, 2005; SAS Institute, 2008, pp. 3952–3962; StataCorp, 2013). These structures can be specified in LGM software; however, it must be done manually and doing so (e.g., by setting appropriate constraints) is not always intuitive.

MEM software often allow for models to be easily extended to multiple levels (i.e., nested random effects). For example, in Proc Mixed, one can simply include an additional RANDOM statement to go from a two-level to a three-level model. For researchers in the behavioral sciences, this is a common occurrence with a natural example being measuring students’ growth over time with students being meaningfully nested within schools. In Proc Mixed, generic RANDOM statements for this example would look similar to:

RANDOM intercept / SUBJECT=School_ID;

RANDOM intercept / SUBJECT=Student_ID (School_ID);

Three-level analysis can be done in Mplus beginning in Version 7 but can be more cumbersome to program.

Given the small sample interest of this tutorial, it is important to note that MEMs will inferentially test regression coefficients using t-tests or F-tests which are known to have more desirable finite sample properties compared to the asymptotic Z-tests and χ² tests that are output by LGM programs such as LISREL and Mplus. The denominator degrees of freedom for finite sample inferential tests is not always clear, and there is debate about what is the appropriate value (see, e.g., Ender, 2011; Schaalje, McBride, & Fellingham, 2002). SPSS is able to calculate Satterthwaite degrees of freedom (Satterthwaite, 1946) and SAS allows users to select from five different options including the oft recommended Kenward–Roger (Kenward & Roger, 1997, 2009). The recently released Stata 14 also now allows users to implement Kenward–Roger among its five denominator degree of freedom options.

MEM software also allows observations to be taken at different time points (e.g., Person 1 is measured at 2 and 4 weeks after baseline and Person 2 is measured 1 and 3 weeks after baseline). Hox (2010) discusses this point in some detail and states that LGMs are inherently a fixed occasion model: If participants are measured at different time points, all possible time points would need to be modeled within the LGM framework.

LGMs

As elaborated upon above, latent basis and second-order growth models cannot be parameterized directly within MEM software. As additional advantages of modeling with LGMs, measurement error in predictor variables can be modeled (Blozis & Cudeck, 1999; Curran, 2003) and categorical latent variables can be added as a natural extension to identify latent subgroups (B. O. Muthén, 2001; B. Muthén & Shedden, 1999; Nagin, 1999), provided that distributional assumptions of the latent variables are met. LGM software (AMOS 23, EQS, LISREL, and Mplus) all can perform FIML to accommodate missing values on both the outcome variable and the predictor variables (under the assumption that data are missing at random), whereas MEMs (SAS Proc Mixed and Stata xtmixed) can use likelihood methods to accommodate missing dependent variables but will listwise delete missing values for predictors (Allison, 2012; Enders, 2010). The listwise deletion of missing predictors in MEM software can be especially problematic with the small samples of interest here because each observation is crucial to retain to maximize power. Lastly, Wu, West, and Taylor (2009) and Hox (2010) point out that LGMs produce fit statistics and indices that assess the fit of the model from a global perspective compared to MEM software which focuses primarily on inferentially testing individual parameters. LGMs also produce modification indices that can suggest which paths to add, delete, or edit to atheoretically improve model fit (Hox, 2010).

Bayesian Markov Chain Monte Carlo (MCMC) Estimation for Small Samples

A source of the small sample bias prevalent with FIML estimation lies in which maximum likelihood estimation is an asymptotic procedure, meaning that it is unbiased for large sample sizes but not necessarily so with small or even moderate samples. On the contrary, Bayesian MCMC does not operate under the same asymptotic property that has led to MCMC being a suggested method for challenging data structures, smaller sample sizes in particular (Dunson, 2000; Gelman, 2006; Hox, van de Schoot, & Matthijsse, 2012; Lee & Song, 2004; Price, 2012; Scheines, Hoijtink, & Boomsma, 1999; Soares, Gonçalves, & Gamerman, 2009; van de Schoot, Broere, Perryck, Zondervan-Zwijnenburg, & van Loey, 2015).⁴ Previous studies have addressed the small sample properties of MCMC in the context of confirmatory factor analytic models (Lee & Song, 2004), multilevel SEM models for survey data (Hox et al., 2012), and for autoregressive time series models (Price, 2012) with the general conclusion that MCMC is preferable over FIML with small sample sizes. Our goal is not to introduce Bayesian MCMC at the most fundamental level, so readers unfamiliar with the basics of Bayesian MCMC are referred to Kruschke, Aguinis, and Joo (2012); van de Schoot et al. (2014); or Zyphur, Oswald, and Rupp (2015) for highly readable introductions to Bayesian MCMC.

However, with small samples, the choice of prior distributions for the factor covariance matrix with MCMC is vitally important, especially considering the increased importance of the factor covariance matrix in LGMs compared to other SEMs (i.e., if the factor variances are not of interest, random effects can be treated as a nuisance and averaged over with design-based methods, eliminating the necessity of LGMs). In the case of small sample sizes, the prior distribution has a much larger influence on the posterior distribution compared to the data, relative to when sample sizes are larger (Depaoli, 2014; Ghosh & Mukerjee, 1992; Lambert, Sutton, Burton, Abrams, & Jones, 2005, van de Schoot et al., 2015). Combined with the common constraint that disallows variances from taking on negative values (see Savalei & Kolenikov, 2008, for a detailed discussion of constrained and unconstrained variances), a sufficiently noninformative prior for the factor covariance matrix (the equivalent of the individual-level variance components in MEMs) has still been a point of contention (Browne & Draper, 2006; Gelman, 2006; Stegmueller, 2013).

In typical implementations of MCMC, univariate variances are assigned an inverse gamma prior when uninformative priors are desired. Gelman (2006) showed, however, that under small sample size conditions, the inverse gamma prior was no longer noninformative and caused the estimates to be biased toward 0. Baldwin and Fellingham (2013) compared the inverse gamma prior to a uniform prior for a univariate random effect in Bayesian MEMs and found intercept variance estimates to be highly underestimated with 16 or fewer clusters.

Although the advantages of MCMC are often discussed, implementing models within a Bayesian framework can at times be difficult because much of the literature and related software is of a fairly high technical level (B. Muthén, 2010). For researchers without a strong background in Bayesian statistics, frequentist models in Mplus can be extended to a Bayesian framework with relative ease compared to specifying the model in more flexible but complex Bayesian specific software such as WinBUGS 1.4.3, Stan 2.6.0, JAGS 4.0.0, SAS Proc MCMC, or certain R packages (e.g., R2Bugs and BRugs). The MCMC capabilities within Mplus can simplify analyses for researchers looking to capitalize on the advantages of MCMC but who may not be as well-versed with software and/or syntax routinely utilized by fully Bayesian statisticians. Because of the broad appeal of modeling Bayesian LGMs in Mplus, this article will emphasize software implementation in Mplus (see Online Appendix A for an overview of fitting Bayesian MCMC models in Mplus).

Demonstrating the Suboptimal Small Sample Performance of Noninformative Priors

A simulation was conducted to demonstrate the poor small sample performance of FIML and traditional noninformative priors and to motivate the need for an alternative. This simulation was meant to be illustrative, not comprehensive, so it featured only three conditions for the number of individuals (20, 30, and 50) and four repeated measures. The generation model and population values were inspired by previous simulation studies by Bauer and Curran (2003) and L. K. Muthén and Muthén (2002). The generated model featured positive linear growth over time, and the latent intercept and latent slope had nonzero variance but their covariance was null in the population.⁵ Two binary individual-level predictors were included for both the latent intercept and the latent slope, one with a 50:50 split reminiscent of biological sex and the other with a 60:40 split more commonly seen in an ethnic minority status indicator. Predictors accounted for approximately 10% of the total variance in the latent intercept and latent slope. Following Bauer and Curran (2003), the residual variances at each time point were selected such that the variance accounted for at each time point was approximately 50%. Figure 2 shows the full path diagram for the model along with the unstandardized population values.

Figure 2.

Path diagram and population values for generation model. ϕ are variances of observed time-invariant predictors, ψ are disturbance variances of the latent growth parameters, γ are effects of observed predictors on latent growth parameters, α are latent factor means, and θ are residual variances.

Models were estimated in Mplus using two different estimation methods: FIML (Estimator = ML) and MCMC (Estimator = Bayes). Two inverse Wishart prior distributions were implemented for the factor covariance matrix in the MCMC conditions. An advantage of the inverse Wishart prior is that the diagonal terms of the covariance matrix follow an inverse gamma distribution, but each draw from the prior is guaranteed to be positive definite provided that the degrees of freedom are greater than or equal to the dimension of the covariance matrix (Lunn, Jackson, Best, Thomas, & Spiegelhalter, 2013). The inverse Wishart is also available in Mplus.

The first prior considered was an improper inverse Wishart W⁻¹ (0, −p −1; the Mplus default)⁶ where the first argument is the scale matrix, the second argument is the degrees of freedom (larger values for degrees of freedom increase the informativeness of the prior), and p is the dimension of the factor covariance matrix. The second was a proper inverse Wishart W⁻¹ (I, p) where I is an identity matrix of dimension p as recommended by Chen (2011, pp. 14–16) and Lunn, Jackson, Best, Thomas, and Spiegelhalter (2013, p. 239).⁷ A related prior that is often recommended is W⁻¹ (I, p + 1; e.g., Asparouhov & Muthén, 2010); this prior was considered, but the results were nearly identical to W⁻¹ (I, p). Therefore, only W⁻¹ (I, p) will be discussed and results reported for the remainder of the article although readers may note that similar results would be obtained with W⁻¹ (I, p + 1).

Mplus defaults for MCMC were implemented that include no thinning, the Gibbs (PX1) algorithm with a maximum of 50,000 iterations, using the Potential Scale Reduction (PSR) convergence criteria (Gelman & Rubin, 1992), and running two chains. The posterior is summarized with the median in accordance with recommendations in Browne and Draper (2006). To briefly explicate on these modeling choices with MCMC, the maximum number of iterations control how long MCMC will continue to update parameter estimates (the Gibbs algorithm is the process used to update parameters) unless convergence is reached prior to the maximum number, which will be assessed with the PSR criteria. With the PSR criteria, multiple chains are run and then the within-chain and between-chain variance are calculated at each iteration, for each parameter. From the within-chain and between-chain variance, the PSR factor is calculated that assesses the proportion of variance attributable between chains. Once there is very little variance between chains for all parameters (indicated by a PSR factor near 1), one can conclude that both chains have arrived at roughly the same location, and convergence has been reached. Thinning results in estimation that is based on every kth iteration rather than each iteration, which can be useful in the presence of high autocorrelation (e.g., if the posterior distribution is thinned by 10, every 10th iteration will be kept for use in estimation and all other replications will be discarded).⁸ The number of simulation study replications within each cell of the design was 2,500.

Figure 3 shows the relative bias for the latent intercept variance (top panel) and the latent slope variance (bottom panel) with the dashed gray lines representing the ±10% cutoff typically used to denote significant bias (Flora & Curran, 2004; Hoogland & Boomsma, 1998). As one would expect, the magnitude of the bias decreases as sample size increases, but bias is rather large at smaller samples, particularly for the variance of the latent slope (note the different scales of the panels in Figure 3).

Figure 3.

Relative bias for latent intercept variance (top) and latent slope variance (bottom).

As a partial explanation for the bias observed with the MCMC conditions, the priors that are generally considered to be noninformative often have large ranges to avoid concentrating too heavily on a specific area of the distribution. However, with small samples where the likelihood contributes relatively less (compared to larger samples) to the posterior distribution, the prior (and thus some larger values that may not be supported by the data) has a larger influence on the posterior, resulting in upward bias. As sample size increases, the influence of the data on the posterior becomes larger and bias diminishes.

DDP Approach for the Factor Covariance Matrix

Casella (2001) mentions three general methods for setting hyperparameters (e.g., the elements of the inverse Wishart scale matrix, S, and the degrees of freedom p) in prior distributions. First, one may select specific fixed values based on prior beliefs, subject matter experts’ opinions, or recommended noninformative priors (as was done in the previous simulation). Second, the hyperparameters can be modeled through an additional level of hierarchy such that the hyperparameters are themselves given a prior distribution (although the hierarchy must terminate with fixed values at some point). Setting a hyperprior for the elements of the scale matrix in an inverse Wishart prior would present some difficulty because the elements of the scale matrix are rather difficult to interpret themselves, never mind deciding on an appropriate prior distribution for them. Third, hyperparameters can be estimated from the data with what have been generally classified as DDPs (see, e.g., Darnieder, 2011).

Because specifying specific values that are commonly stated to be noninformative is difficult with small samples and did not perform well in the previous illustrative simulation and introducing another level of hierarchy would be trying, the following sections will discuss a data-dependent method for estimating values for the hyperparameters. That is, rather than selecting arbitrary values for the elements of the scale matrix (e.g., a matrix of zeroes or an identity matrix) or giving these hyperparameters a prior distribution, values for these elements will be estimated from the data.

DDPs

DDPs and the related and popular empirical Bayes approach (classified under the broader umbrella of DDPs in Darnieder, 2011)⁹ have existed in the literature for decades, going as far back as work on nonparametric empirical Bayes in Robbins (1955).

As implied by the name, the general idea behind DDPs is to use the data to inform values for the hyperparameters that are typically estimated either through frequentist methods such as maximum likelihood or a method of moments (Casella, 2001). In essence, the process of selecting the prior is taken out of the researchers’ hands and delegated to the data. Once the hyperparameters are estimated from the data, the estimated values are substituted into the prior distribution(s) as constants and the Bayesian analysis proceeds as usual. The combination of frequentist methods to obtain hyperparameters for the prior distribution with Bayesian machinery to estimate model parameters has led Efron (2010) to state that “Empirical Bayes blurs the line between testing and estimation as well as between frequentism and Bayesianism,” (p. X of the prologue). In the viewpoint of Berger (2006), DDPs may be necessitated when noninformative priors are poorly behaved and Berger asserts that noninformative priors are only approximations to useful informative priors in the first place. In a fully Bayesian analysis with small sample sizes, no prior distribution can be completely noninformative and selecting hyperparameters for a prior on more-or-less arbitrary grounds can lead estimates astray as was seen in the previous illustrative simulation.

DDPs are not without their critics. The goal of this tutorial is to outline the use of DDPs rather than present the differential philosophical views concerning their use. Readers interested in the detractions of DDPs are referred to Deely and Lindley (1981), Lindley (1969), or Robert (2007). Despite reservations regarding DDPs, constructing a reasonable prior for the factor covariance matrix with small samples is quite difficult and the resulting estimates when using traditional noninformative priors or FIML can be quite biased as seen in the previous simulation. Thus, following Berger (2006), the subsequent sections of this tutorial will demonstrate the utility of DDPs in the context of growth models with small samples, particularly in Mplus.

To construct a DDP in Mplus, first, the model is estimated using FIML and the estimates for the factor variances and covariances are recorded. Then, the estimates are substituted into the scale matrix of the inverse Wishart prior and scaled by the degrees of freedom (in accordance with Mplus’ parameterization of the inverse Wishart distribution). Then, the model is estimated with straightforward MCMC. Parameter estimates are obtained in the same way as they would were the model estimated with fully Bayesian estimation.

As a conceptual explanation of this strategy, with small samples, the prior distribution has a stronger impact relative to the data than with large samples. With small samples, the contribution from the data does not separate the posterior very well from the prior due to the limited information provided by the small sample. The FIML estimates are known to be biased with small samples, but, assuming that FIML estimation converges, the estimates will be in a general vicinity of the population values. By estimating the hyperparameters with FIML, the prior is generally close to the population value, so the limited information from the data (limited in the sense that the sample size is small) will not need to distance itself much from the prior to arrive at the appropriate posterior distribution. Despite the small sample size, the data will update the posterior distribution enough so that the bias is removed from the FIML estimates. Although Darnieder (2011) mentions the ability of DDPs to drive the posterior distribution as a negative aspect of DDPs, with small samples, the MCMC machinery can be manipulated to take advantage of the fact that FIML is known to be close to population values but consistently downwardly biased.

DDPs With LGM Software-Specific Models

Recall that software implementations of MEM are less flexible than LGM and some models can only be fit in LGM software where REML is not available. Therefore, a second simulation will show how DDPs perform for LGM software-specific models (a latent basis and second-order growth model) with small samples. Again, since the general DDP methodology has been widely studied and supported (although not without detractors, as noted above), this simulation is intended to be demonstrative rather than comprehensive so it will only feature three sample size conditions (20, 30, and 50) and four repeated measures. The latent basis model is very similar to Figure 2 except that the loadings from the slope factor to Y₃ and Y₄ will no longer be constrained but rather estimated from the data. The population values for these parameters are 5 and 7, respectively. The second-order growth model is also similar to Figure 2 except that the four observed Y variables are replaced with four latent variables that each have a three variable measurement model that is identical across time as required by the model (Hancock et al., 2001). The growth of these latent variables was linear and constrained across time points as in Figure 2. Because the models were more complex and required more computational time, 1,000 simulation study replications were used rather than 2,500.

Figure 4 shows the relative bias for the latent intercept variance (top left panel), latent slope variance (top right panel), and the estimated loading from the slope to Y₃ (bottom panel; bias for loadings between Y₃ and Y₄ was very similar so only Y₃ is shown) across estimation methods and sample size for the latent basis model. FIML did not perform too poorly when estimating the intercept and slope variances; however, the estimates of the loadings were highly upwardly biased. The noninformative priors often performed poorly (although the scaled identity inverse Wishart performed well for the intercept variance). The DDP was not pristine (slope variance relative bias exceeded 15% with 20 observations and exceeded 10% with 30 observations), but it performed much better than FIML or the noninformative priors for estimates of all three parameters. Again, note that this model cannot be fit as an MEM in programs that use REML, so very few alternatives exist to accommodate the small sample size in LGM software. In addition to abating the bias in the variance of the growth parameters, DDP provided the best coverage intervals and the best statistical power for all parameters in the model. In particular, FIML consistently had coverage intervals that were too short (i.e., standard error estimates were too small which would inflate type I error rates), and noninformative Bayesian methods estimated the standard errors of the loadings from the slope to the observed variables very poorly. Table B1 in Online Appendix B shows specific values for coverage intervals and statistical power for each parameter in each sample size condition for the latent basis model.

Figure 4.

Relative bias for intercept variance (top left), slope variance (top right), and slope loading (bottom) for the latent basis model.

Figure 5 shows the relative bias for the latent intercept variance (top panel) and latent slope variance (bottom panel) across estimation methods and sample sizes for the second-order growth model. The estimates of the factor loadings in the first-order measurement model were reasonable across methods and sample sizes, so they are not reported. Similar to results in Figures 3 and 4, FIML underestimated the factor variances by a fair amount and both inverse Wishart priors performed particularly poorly for the slope variance. Again, DDPs did not remove all concerns about small sample bias (−24% bias for the slope variance with 20 observations), but it did appreciably mitigate the bias with other options (−52% for FIML, 238% for improper inverse Wishart, and 177% for identity Wishart for the slope variance with 20 observations). As with the latent basis model, FIML consistently had coverage intervals that were far too short. All MCMC and DDP conditions have reasonable converge intervals that tended to be too wide for the latent intercept variance and latent slope variance (although the standard errors for these parameters may be of questionable utility anyway because the variance components are unlikely to be symmetric). Table B2 in Online Appendix B shows specific values for coverage intervals and statistical power for each parameter in each sample size condition for the latent basis model.

Figure 5.

Relative bias for intercept variance (top) and slope variance (bottom) for the second-order growth model.

To provide evidence for similarity between REML and the DDP approach for linear growth, the data from the first simulation whose results were presented in Figure 3 were modeled as an MEM since the model could be equivalently parameterized and because the population values were known. Table 1 shows the relative bias of the factor variances for FIML, REML, and DDP and full simulation results can be viewed in Table B3 in Online Appendix B. As can be seen, the REML and DDP values are less biased even for small sample sizes, especially compared to FIML.

Table 1.

Factor Variance Relative Bias Comparison of FIML, REML, and DDP

	Sample Size
	20			30			50
	FIML	REML	DDP	FIML	REML	DDP	FIML	REML-	DDP
Var(I)	−25	2	9	−17	0	6	−11	−2	3
Var(S)	−55	6	−2	−34	5	0	−24	2	0
Corr(I, S)	1	0	2	1	0	0	1	1	−1

Note. Bolded values indicate the considerable bias based on Hoogland and Boomsma (1998). FIML = full maximum likelihood; REML = restricted maximum likelihood; DDP = data-dependent prior.

Demonstration of DDP for Linear Growth Model

To elucidate the process of constructing a DDP with an example, consider the classical Potthoff and Roy (1964) data set that models the growth of 27 children’s dental measurements at four occasions. The model has one binary exogenous predictor (biological sex) that predicts both the latent intercept and the latent slope, the intercept and slope covary, and the error variances have a homogeneous diagonal structure. A path diagram for this model is shown in Figure 6.

Figure 6.

Path diagram for Potthoff and Roy (1964) model. ψ are disturbance co/variances of the latent growth parameters, γ are effects of observed predictors on latent growth parameters, α are latent factor means, and θ are residual variances.

The REML random effect covariance matrix for the Pothoff and Roy data is $[\begin{matrix} 3.23 & - 0.06 \\ - 0.06 & 0.13 \end{matrix}]$ . When the model is run as an LGM in Mplus using standard FIML, the factor covariance matrix is $[\begin{matrix} 2.91 & - 0.02 \\ - 0.02 & 0.10 \end{matrix}]$ , showing the downward small sample bias with FIML (10% for the intercept variance and 23% for the slope variance compared to REML).

Although the FIML estimates are known to be biased, they can still be useful to help assuage the bias through a DPP. Because the factor covariance matrix is multivariate, an inverse Wishart prior of the form $W^{- 1} (S, d)$ where d = p, the dimension of the factor covariance matrix, and S is the scale matrix is appropriate. The FIML estimates will be used to set the hyperparameters in S. In Mplus and many other Bayesian software programs, the scale matrix for the Wishart or inverse Wishart distributions represents sums of squares and cross-products rather than precisions or variances (see Note 6). As Lunn et al. (2013) state, “it is advisable, whether attempting to be informative or not, to set [S] equal to [d] multiplied by some prior guess at the value of the between-unit covariance matrix” (p. 226). As noted previously, the FIML estimates will be used as the “prior guess.” Thus, in Mplus, S in this example would be equal to the FIML estimates multiplied by the degrees of freedom,

(\begin{matrix} d (2.91) & d (- 0.02) \\ d (- 0.02) & d (0.10) \end{matrix}) = (\begin{matrix} 5.82 & - 0.04 \\ - 0.04 & 0.20 \end{matrix})

where d = p = 2.¹⁰ After running an MCMC analysis with the DDP, the factor covariance matrix is estimated to be $[\begin{matrix} 3.31 & - 0.04 \\ - 0.04 & 0.11 \end{matrix}]$ which is much more in-line with the REML estimates than the FIML estimates.

DDP Compared to Other Priors Example for Linear Growth

The Potthoff and Roy (1964) model will be estimated in Mplus with each of the other priors included in the previous simulation. Mplus code for each model is provided in Table 2 as is code that could be used in MEM software to obtain REML estimates. Of particular interest to users of statistical methods who may not necessarily be statisticians, in Mplus, the code for conducting an MCMC analysis is largely equivalent to a frequentist analysis. The model is specified identically and one only needs to list the estimator as BAYES in the ANALYSIS step and list the desired priors after the model should they be different from the Mplus defaults.

Table 2.

Mplus and General MEM Software Code for Potthoff and Roy (1964) Model

FIML	Improper Inverse Wishart^a	Identity Inverse Wishart	DDP
ANALYSIS: estimator = ML;	ANALYSIS: estimator = Bayes; fbiterations = 1e5; convergence = 1e-8;	ANALYSIS: estimator = Bayes; fbiterations = 1e5; convergence = 1e-8;	ANALYSIS: estimator = Bayes; fbiterations = 1e5; convergence = 1e-8;
MODEL: i s \| y1@0 y2@1 y3@2 y4@ 3; y1-y4 (a); [y1-y4@0]; [i s]; i (b) ; s (c); S WITH I (d); I on sex; S on sex;	MODEL: i s \| y1@0 y2@1 y3@2 y4@ 3; y1-y4 (a); [y1-y4@0]; [i s]; i (b) ; s (c); S WITH I (d) ; I on sex; S on sex;	MODEL: i s \| y1@0 y2@1 y3@2 y4@ 3; y1-y4 (a); [y1-y4@0]; [i s]; i (b) ; s (c); S WITH I (d); I on sex; S on sex;	MODEL: i s \| y1@0 y2@1 y3@2 y4@ 3; y1-y4 (a); [y1-y4@0]; [i s]; i (b) ; s (c); S WITH I (d); I on sex; S on sex;
		model priors: b∼IW(2,2); c∼IW(2,2); d∼IW(0,2);	model priors: b∼IW(5.82,2); c∼IW(0.20,2); d∼IW (-.04,2);
SAS Proc Mixed	R lme4	SPSS Mixed	Stata xtmixed
proc mixed data=PR64; model y=time\|sex/ solution; random int time/ sub=person type=un; run;	require(lme4) PR<-lmer(y∼Time+ sex+ Time*sex+ (1+Time\|Person), data=PR64) summary(PR)	mixed y with time sex /print = solution testcov /method = reml /fixed = intercept time sex time*sex /random intercept time \| subject(person) covtype(un).	xtmixed y time sex time_sex \|\| Person: Time, reml cov(unstruct)

Note. The fbiterations and convergence options are not required but are included because the model is fairly simple and requiring more iterations and stricter convergence criteria than the default in Mplus was not computationally demanding. These options also help reproducibility should readers want to replicate the results. MEM = mixed-effects model; DDP = data-dependent prior.

^aThe improper inverse Wishart Prior is the Mplus default, so the model priors statement is not needed so long as the estimator is set to Bayes.

Table 3 shows the model estimates and standard errors (posterior distribution standard deviations for MCMC models) for the two noninformative prior distributions, FIML, DDP, and REML (as implemented in SAS Proc Mixed). Although no population values are available, the data will be compared to the REML estimates. The intercept and slope factor means (equivalent to the fixed effect for the intercept and slope in an MEM) were nearly identical across methods and are consequently not reported for brevity.

Table 3.

Model Results for Varying Estimation Methods With Potthoff and Roy (1964) Data

	Parameter Estimate (Standard Error/Posterior Standard Deviation)
Method	Var(I)	Var(S)	Corr(I, S)	γ₀₁	γ₁₁
Improper IW	4.27 (2.02)	.25 (.21)	−.20 (.49)	−1.41 (.95)	−.61 (.32)
Identity IW	3.09 (1.37)	.35 (.16)	−.21 (.37)	−1.41 (.84)	−.61 (.34)
FIML	2.91 (1.14)	.10 (.14)	−.02 (.30)	−1.41 (.79)	−.61 (.26)
DDP	3.31 (1.36)	.11 (.10)	−.02 (.27)	−1.41 (.86)	−.61 (.28)
REML	3.23 (1.28)	.13 (.15)	−.03 (.29)	−1.41 (.83)	−.61 (.27)

Note. Correlations between the intercept and slope and provided instead of the covariance because correlations are often more readily interpretable. FIML = full maximum likelihood; REML = restricted maximum likelihood; DDP = data-dependent prior; IW= inverse Wishart.

The Mplus default improper inverse Wishart prior yielded elements of the factor covariance matrix much larger than REML (4.27 vs. 3.23 for the latent intercept; 0.25 vs 0.13 for the latent slope) and resulted in similarly large posterior standard deviations. The identity inverse Wishart prior estimated the latent intercept variance fairly close to REML but estimated the variance of the latent slope at close to 3 times the REML estimate. DDP estimated the elements of the factor covariance matrix closest to REML of the competing methods. The estimates of the variance of the latent slope and factor correlation were almost exactly the same and were both far closer to the REML estimate than any other method. The DDP latent intercept variance estimate was also within 3% of REML and was also the closest to REML of any competing method. As an aside, it should be noted that the REML estimates are unbiased on expectation but may not necessarily reflect the population values exactly for this particular data set, so the REML estimates should not be considered exact truth but only the estimate from the best small sample estimator developed to date for small sample correlated data models fitting in the MEM framework.

Discussion

Mathematically, MEMs and LGMs are often considered to be one in the same with the nomenclature changing only across disciplinary lines. The types of models that MEM software and LGM software can fit largely overlap; however, the mapping is not one to one. The requirements of MEM software that time is a known quantity located in the fixed effect design matrix, that all predictors and outcomes are manifest, and that missing predictors are listwise deleted in software restrict the flexibility of models that can be implemented in MEM software. Models that introduce latent predictors or latent outcomes, some models for nonlinear growth, and treating missing predictors with likelihood methods fall solely within the realm of LGM software, which have less developed methods to accommodate small samples.

For as long as priors that are dependent on data have existed, critics have aptly pointed out that these methods “double dip” and use the data twice. This point is certainly valid with the DDP method demonstrated in this tutorial and the goal was not to debate the merits or faults of the philosophy of this approach. Methodologists who subscribe to a more pure statistical and philosophical ideology may eschew the ideas presented here; however, random-effects models in general (which subsumes both MEMs and LGMs) most commonly estimate the values of the random effects with empirical Bayes, a type of DDP (Laird & Ware, 1982; Raudenbush & Bryk, 2002, p. 47). Therefore, opposition to DDPs in the prior for the random effect variance but acceptance of DDPs in estimating the values of the random effects is a rather nuanced criticism.

It is important to reiterate that the emphasis of this tutorial was on a practical and accessible method to address the small sample problem in LGM software that allows researchers to navigate within familiar software. The noninformative prior distributions that were included in the simulations are more on the tame end of the spectrum of possible prior distributions and were selected with the express purpose of being available in the user-friendly Mplus software or easily programmable in other Bayesian-specific software. For researchers who are well versed in Bayesian methods, prior distributions of increasing complexity could also be fit that may perform better than the prior distributions that were investigated here. For instance, for models with only a random intercept, Gelman (2006) suggested either a half-t or a half-Cauchy prior distribution over the inverse gamma or uniform distribution for the univariate variance components. Barnard, McColloch, and Meng (2000) have alternatively proposed decomposing the covariance matrix into a correlation matrix and standard deviations, with each having their own priors and O’Malley and Zaslavsky (2005) have proposed a scaled inverse Wishart that allows the off-diagonal terms in the scale matrix to be estimated without affecting the diagonal elements by pre- and postmultiplying the scale matrix S by Diag(ξ) where ξ is a p × p matrix of scale parameters that are given their own prior. These distributions are not easily programmable in Mplus but may potentially outperform the priors used in the simulation if implemented in more specialized Bayesian software such as WinBUGS, Stan, or JAGS. Nonetheless, there is fair evidence that, given the ease of implementation, DDPs do a fairly reasonable job reducing the sizable small sample bias that was observed with FIML or traditional noninformative priors for samples sizes of 50 or less. As models grow more complex, the required sample size to obtain unbiased estimates will increase likewise and DDPs may have potential utility for sample sizes into the low triple digits.

Seeing as REML estimation has not been well developed broadly (see Cheung, 2013, for some instances where REML can be applied to linear LGMs), DDPs may suffice to provide a viable, relatively unbiased alternative for LGMs with small samples, especially since not all LGMs can be conventionally reparameterized into an MEM to take advantage of REML. However, if one has a small sample and can parameterize the model in the MEM framework, it would be recommended to do so and take advantage of REML. As such, research into continued development of analogs of REML and Kenward–Roger techniques for LGMs could be a valuable area of future research, especially given the heightened flexibility of LGMs.

As additional considerations not discussed in-depth previously, Bayesian methods considered in this tutorial still require the normality assumption to be met (nonparametric Bayesian methods not requiring normality also exist); data in the simulation studies were generated such that this assumption was upheld but data in applied studies may not always meet this requirement and this assumption can be difficult to assess in general with small samples. Although touched upon only briefly, the power of DDPs almost unanimously surpassed that of other methods in conditions where comparisons were warranted (e.g., well-behaved type I errors rates), albeit by only a few percentage points (see tables in Online Appendix B for exact values). In the illustrative simulations provided in this study, DDP provided reasonable standard error estimates. It should be noted, however, that DDP may underestimate standard errors because the variability of the estimates of the hyperparameters in the random effect covariance matrix is ignored (e.g., Carlin & Louis, 2000). This shortcoming has been widely researched and remedies have been advanced in the statistical literature (e.g., Carlin & Gelfand, 1990, 1991; Laird & Louis, 1987).

With regard to data model fit, with small sample sizes and FIML, the minimum fit function χ² produces vastly inflated type I error rates, meaning that well-fitting models are rejected too often (e.g., Kenny & McCoach, 2003). Small sample corrections such as those by Bartlett (1950), Swain (1975), and Yuan (2005) have been investigated in previous studies (see Fouladi, 2000; Herzog & Boomsma, 2009; Nevitt & Hancock, 2004; Savalei, 2010). From the somewhat limited selection of prior distributions in Mplus, other priors may be potentially useful to be noninformative and/or beneficial with small samples. Asparouhov and Muthén (2010) mention Γ⁻¹ (1, 2), W⁻¹ (0, 0), and W⁻¹ (I, p + 1) as other potential options. This study used prior distributions that have been most strongly advocated in the literature but that is not to say that these other distributions are not without merit, particularly with small samples.

As a final remark, as Bayesian software becomes more widely accessible and user-friendly, researchers must still be cognizant of the ramifications of their modeling decisions, especially with small samples. More specifically, as seen in this tutorial, although software such as Mplus provide helpful default priors, blindly proceeding with the defaults can lead to seriously biased results and unwarranted conclusions. In fact, the Mplus default prior consistently produced the most seriously biased estimates (although do keep in mind the illustrative nature of the simulations and the small number of conditions). Truly noninformative priors are more myth than reality with small samples and Bayesian methods are useful tools for these scenarios but are by no means a magic bullet.

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

Notes

Supplemental Material

The online appendices are available at

References

Allison

P. D.

(2012, 4). Handling Missing data by maximum likelihood. Statistics and Data Analysis Keynote Presentation at the SAS Global Forum, Orlando, FL.

Asparouhov

Muthén

(2010). Bayesian analysis of latent variable models using Mplus, version 2. Unpublished manuscript. Retrieved from http://www.statmodel.com

Baldwin

S. A.

Fellingham

G. W.

(2013). Bayesian methods for the analysis of small sample multilevel data with a complex variance structure. Psychological Methods, 18, 151–164.

Barnard

McCulloch

Meng

X. L.

(2000). Modeling covariance matrices in terms of standard deviations and correlations, with application to shrinkage. Statistica Sinica, 10, 1281–1312.

Bartlett

M. S.

(1950). Tests of significance in factor analysis. British Journal of Statistical Psychology, 3, 77–85.

Bauer

D. J.

Curran

P. J.

(2003). Distributional assumptions of growth mixture models: Implications for overextraction of latent trajectory classes. Psychological Methods, 8, 338–363.

Berger

(2006). The case for objective Bayesian analysis. Bayesian Analysis, 1, 385–402.

Biesanz

J. C.

Deeb-Sossa

Papadakis

A. A.

Bollen

Curran

(2004). The role of coding time in estimating and interpreting growth curve models. Psychological Methods, 9, 30–52.

Blozis

S. A.

Cudeck

(1999). Conditionally linear mixed-effects models with latent variable covariates. Journal of Educational and Behavioral Statistics, 24, 245–270.

10.

Browne

W. J.

Draper

(2006). A comparison of Bayesian and likelihood-based methods for fitting multilevel models. Bayesian Analysis, 1, 473–514.

11.

Carlin

B. P.

Gelfand

A. E.

(1990). Approaches for empirical Bayes confidence intervals. Journal of the American Statistical Association, 85, 105–114.

12.

Carlin

B. P.

Gelfand

A. E.

(1991). A sample reuse method for accurate parametric empirical Bayes confidence intervals. Journal of the Royal Statistical Society. Series B, 53, 189–200.

13.

Carlin

B. P.

Louis

T. A.

(2000). Empirical Bayes: Past, present and future. Journal of the American Statistical Association, 95, 1286–1289.

14.

Casella

(2001). Empirical Bayes Gibbs sampling. Biostatistics, 2, 485–500.

15.

Chen

(2011). The RANDOM statement and more: Moving on with PROC MCMC. In Proceedings of the SAS Global Forum 2011 Conference. Cary, NC: SAS Institute, Inc.

16.

Cheung

M. W. L.

(2013). Implementing restricted maximum likelihood estimation in structural equation models. Structural Equation Modeling: A Multidisciplinary Journal, 20, 157–167.

17.

Curran

P. J.

(2003). Have multilevel models been structural equation models all along? Multivariate Behavioral Research, 38, 529–569.

18.

Darnieder

W. F.

(2011). Bayesian methods for data-dependent priors (Doctoral dissertation). The Ohio State University, Columbus, OH.

19.

Depaoli

(2014). The impact of inaccurate “informative” priors for growth parameters in Bayesian growth mixture modeling. Structural Equation Modeling: A Multidisciplinary Journal, 21, 239–252.

20.

de Leeuw

Meijer

Goldstein

(2008). Handbook of multilevel analysis. New York, NY: Springer.

21.

Deely

J. J.

Lindley

D. V.

(1981). Bayes empirical Bayes. Journal of the American Statistical Association, 76, 833–841.

22.

Dunson

D. B.

(2000). Bayesian latent variable models for clustered mixed outcomes. Journal of the Royal Statistical Society: Series B, 62, 355–366.

23.

Efron

(2010). Large-scale inference: empirical Bayes methods for estimation, testing, and prediction. New York, NY: Cambridge University Press.

24.

Ender

(2011, 7). xtmixed and denominator degrees of freedom: Myth or magic. Stata Conference, Chicago, IL.

25.

Enders

C. K.

(2010). Applied missing data analysis. New York, NY: Guilford Publications.

26.

Flora

D. B.

Curran

P. J.

(2004). An empirical evaluation of alternative methods of estimation for confirmatory factor analysis with ordinal data. Psychological Methods, 9, 466–491.

27.

Fouladi

R. T.

(2000). Performance of modified test statistics in covariance and correlation structure analysis under conditions of multivariate nonnormality. Structural Equation Modeling, 7, 356–410.

28.

Gelman

(2006). Prior distributions for variance parameters in hierarchical models. Bayesian Analysis, 1, 515–534.

29.

Gelman

Rubin

D. B.

(1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7, 457–472.

30.

Ghosh

J. K.

Mukerjee

(1992). Non-informative priors (with discussion). In Bernardo

J. M.

Berger

J. O.

Dawid

A. P.

Smith

A. F. M.

(Eds.), Bayesian statistics 4 (pp. 195–2010). Oxford, England: University Press.

31.

Goldstein

(1989). Restricted unbiased iterative generalized least-squares estimation. Biometrika, 76, 622–623.

32.

Grimm

K. J.

Ram

Hamagami

(2011). Nonlinear growth curves in developmental research. Child Development, 82, 1357–1371.

33.

Hancock

G. R.

Kuo

W. L.

Lawrence

F. R.

(2001). An illustration of second-order latent growth models. Structural Equation Modeling, 8, 470–489.

34.

Herzog

Boomsma

(2009). Small-sample robust estimators of noncentrality-based and incremental model fit. Structural Equation Modeling: A Multidisciplinary Journal, 16, 1–27.

35.

Hoogland

J. J.

Boomsma

(1998). Robustness studies in covariance structure modeling: An overview and a meta-analysis. Sociological Methods & Research, 26, 329–367.

36.

Hox

(2010). Multilevel analysis: Techniques and applications. Mahwah, NJ: Routledge.

37.

Hox

van de Schoot

Matthijsse

(2012). How few countries will do? Comparative survey analysis from a Bayesian perspective. Survey Research Methods, 6, 87–93.

38.

Kenny

D. A.

McCoach

D. B.

(2003). Effect of the number of variables on measures of fit in structural equation modeling. Structural Equation Modeling, 10, 333–351.

39.

Kenward

M. G.

Roger

J. H.

(1997). Small sample inference for fixed effects from restricted maximum likelihood. Biometrics, 53, 983–997.

40.

Kenward

M. G.

Roger

J. H.

(2009). An improved approximation to the precision of fixed effects from restricted maximum likelihood. Computational Statistics & Data Analysis, 53, 2583–2595.

41.

Kruschke

J. K.

Aguinis

Joo

(2012). The time has come Bayesian methods for data analysis in the organizational sciences. Organizational Research Methods, 15, 722–752.

42.

Laird

N. M.

Louis

T. A.

(1987). Empirical Bayes confidence intervals based on bootstrap samples. Journal of the American Statistical Association, 82, 739–750.

43.

Laird

N. M.

Ware

J. H.

(1982). Random-effects models for longitudinal data. Biometrics, 38, 963–974.

44.

Lambert

P. C.

Sutton

A. J.

Burton

P. R.

Abrams

K. R.

Jones

D. R.

(2005). How vague is vague? A simulation study of the impact of the use of vague prior distributions in MCMC using WinBUGS. Statistics in Medicine, 24, 2401–2428.

45.

Lee

S. Y.

Song

X. Y.

(2004). Evaluation of the Bayesian and maximum likelihood approaches in analyzing structural equation models with small sample sizes. Multivariate Behavioral Research, 39, 653–686.

46.

Lindley

D. V.

(1969). Discussion of compound decisions and empirical Bayes, Copas

J. B.

. Journal of the Royal Statistical Society, Series B, 31, 419–421.

47.

Lunn

Jackson

Best

Thomas

Spiegelhalter

(2013). The BUGS book. Boca Raton, FL: CRC Press.

48.

McArdle

J. J.

(1988). Dynamic but structural equation modeling of repeated measures data. In Nesselroade

J. R.

Cattell

R. B.

(Eds.), The handbook of multivariate experimental psychology (pp. 561–614). New York, NY: Plenum Press.

49.

McArdle

J. J.

Epstein

(1987). Latent growth curves within developmental structural equation models. Child Development, 58, 110–133.

50.

McNeish

D. M.

Stapleton

L. M.

(2014). The effect of small sample size on two level model estimates: A review and illustration. Educational Psychology Review. doi:10.1007/s10648-014-9287-x

51.

Meredith

Tisak

(1990). Latent curve analysis. Psychometrika, 55, 107–122.

52.

Muthén

(2010). Bayesian analysis in Mplus: A brief introduction. Unpublished manuscript. Retrieved from http://www.statmodel.com

53.

Muthén

B. O.

(2001). Second-generation structural equation modeling with a combination for categorical and continuous latent variables: New opportunities for latent class/latent curve modeling. In Sayers

Collins

(Eds.), New methods for the analysis of change (pp. 291–322). Washington, DC: American Psychological Association.

54.

Muthén

L. K.

Muthén

B. O.

(2002). How to use a Monte Carlo study to decide on sample size and determine power. Structural Equation Modeling: A Multidisciplinary Journal, 9, 599–620.

55.

Muthén

Shedden

(1999). Finite mixture modeling with mixture outcomes using the EM algorithm. Biometrics, 55, 463–469.

56.

Nagin

D. S.

(1999). Analyzing developmental trajectories: a semiparametric, group-based approach. Psychological Methods, 4, 139–157.

57.

Nevitt

Hancock

G. R.

(2004). Evaluating small sample approaches for model test statistics in structural equation modeling. Multivariate Behavioral Research, 39, 439–478.

58.

Malley

A. J.

Zaslavsky

A. M.

(2005). Variance-covariance functions for domain means of ordinal survey items. Survey Methodology, 31, 169–182.

59.

Peugh

J. L.

Enders

C. K.

(2005). Using the SPSS mixed procedure to fit cross-sectional and longitudinal multilevel models. Educational and Psychological Measurement, 65, 717–741.

60.

Potthoff

R. F.

Roy

S. N.

(1964). A generalized multivariate analysis of variance model useful especially for growth curve problems. Biometrika, 51, 313–326.

61.

Price

L. R.

(2012). Small sample properties of Bayesian multivariate autoregressive time series models. Structural Equation Modeling: A Multidisciplinary Journal, 19, 51–64.

62.

Raudenbush

S. W.

Bryk

A. S.

(2002). Hierarchical linear models: Applications and data analysis methods. Thousand Oaks, CA: Sage.

63.

Robbins

(1955). An empirical Bayes approach to statistics. In Neyman

. (Ed.), Proceedings of the third Berkeley symposium on mathematical statistics and probability (pp. 157–164). Berkeley: University of California Press.

64.

Robert

(2007). The Bayesian choice: A decision-theoretic motivation. New York, NY: Springer.

65.

SAS Institute, Inc. (2008). SAS/STAT^® 9.2 user’s guide. Cary, NC: Author.

66.

Satterthwaite

F. E.

(1946). An approximate distribution of estimates of variance components. Biometrics Bulletin, 2, 110–114.

67.

Savalei

(2010). Small sample statistics for incomplete nonnormal data: Extensions of complete data formulae and a Monte Carlo comparison. Structural Equation Modeling, 17, 241–264.

68.

Savalei

Kolenikov

(2008). Constrained versus unconstrained estimation in structural equation modeling. Psychological Methods, 13, 150–170.

69.

Schaalje

G. B.

McBride

J. B.

Fellingham

G. W.

(2002). Adequacy of approximations to distributions of test statistics in complex mixed linear models. Journal of Agricultural, Biological, and Environmental Statistics, 7, 512–524.

70.

Scheines

Hoijtink

Boomsma

(1999). Bayesian estimation and testing of structural equation models. Psychometrika, 64, 37–52.

71.

Soares

T. M.

Gonçalves

F. B.

Gamerman

(2009). An integrated Bayesian model for DIF analysis. Journal of Educational and Behavioral Statistics, 34, 348–377.

72.

StataCorp. (2013). Stata 13 base reference manual. College Station, TX: Stata Press.

73.

Stegmueller

(2013). How many countries for multilevel modeling? A comparison of frequentist and Bayesian approaches. American Journal of Political Science, 57, 748–761.

74.

Swain

A. J.

(1975). Analysis of parametric structures for variance matrices (Unpublished doctoral dissertation). Department of Statistics, University of Adelaide, Australia.

75.

van de Schoot

Broere

J. J.

Perryck

K. H.

Zondervan-Zwijnenburg

Van Loey

N. E.

(2015). Analyzing small data sets using Bayesian estimation: The case of posttraumatic stress symptoms following mechanical ventilation in burn survivors. European Journal of Psychotraumatology, 6, 25216.

76.

van de Schoot

Kaplan

Denissen

Asendorpf

J. B.

Neyer

F. J.

Aken

M. A.

(2014). A gentle introduction to Bayesian analysis: Applications to developmental research. Child Development, 85, 842–860.

77.

Wakefield

(2013). Bayesian and frequentist regression methods. New York, NY: Springer.

78.

West

S. G.

Taylor

A. B.

(2009). Evaluating model fit for growth curve models: Integration of fit indices from SEM and MLM frameworks. Psychological Methods, 14, 183–201.

79.

Yuan

K.-H.

(2005). Fit indices versus test statistics. Multivariate Behavioral Research, 40, 115–148.

80.

Zyphur

M. J.

Oswald

F. L.

Rupp

D. E.

(2015). Rendezvous Overdue Bayes analysis meets organizational research. Journal of Management, 41, 387–389.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.20 MB