Nuances of Information Criteria for Bayesian Psychometric Models

DIC

DIC is the simplest of the criteria considered here. Assuming a generic parameter vector θ and S posterior samples, it can be written as

where p_D is the effective number of parameters, ${\mathit{\theta }}^s$ is the sth posterior sample, and $\bar{\mathit{\theta }}$ is the posterior mean (or another posterior measure of central tendency). Notice that the in-sample predictive accuracy log p(y | $\bar{\mathit{\theta }}$) is multiplied by −2 and that p_D is generally positive. So smaller DIC values are favored over larger values.

Celeux et al. (2006) considered many variations of DIC, some of which replace the 2log p(y | $\bar{\mathit{\theta }}$) term with other measures of in-sample predictive accuracy. Du et al. (2023) studied some of these variations as applied to multilevel models and found evidence that they can lead to improved model selection decisions over the original DIC. While we do not describe all the variations here, we do note that they can lead to additional confusion about what DIC metric is being reported in a given application.

WAIC

WAIC is a “fully Bayesian” information criterion that utilizes pointwise predictive accuracy terms, averaging them over the entire posterior distribution. This is in contrast to DIC, which uses the overall model likelihood to approximate predictive accuracy only at the posterior mean. WAIC assumes that each datapoint is independent of the others, which immediately holds for simple linear models but becomes more complicated for mixed models and for psychometric models. Assuming N independent datapoints, we can write WAIC as

PSIS-LOO

Watanabe (2010) showed that WAIC is asymptotically equivalent to leave-one-out cross-validation, whereby we hold out each datapoint, estimate the model, and compute the likelihood of the held-out datapoint given the estimated parameters. Vehtari et al. (2017) propose a Pareto-smoothed importance sampling (PSIS) method for directly approximating leave-one-out cross-validation, as opposed to indirectly approximating leave-one-out cross-validation via WAIC. This method, PSIS-LOO, is implemented in the R package loo (Vehtari et al., 2024) and is commonly obtained alongside WAIC. Just like WAIC, the PSIS-LOO metric is based on the pointwise likelihoods evaluated at each posterior sample. The algorithm underlying PSIS-LOO produces diagnostic metrics (“Pareto k” metrics) that can signal problems with the approximation as well as influential observations in the dataset. The availability of diagnostics leads us to prefer PSIS-LOO over WAIC, though the two metrics are usually very similar to one another for the psychometric models considered in this paper.

Choice of Likelihood

The above criteria involve either the overall model likelihood p(y | θ) or the pointwise likelihoods p(y_i | θ). As seen by Equations (5) and (6), a psychometric model with latent variables has multiple likelihoods. The specific choice of likelihood influences the values of the information criteria and the type of predictive generalization that the information criteria approximate (see Merkle et al., 2019).

MCMC estimation of Bayesian psychometric models typically involves sampling latent variables for each person. By conditioning on these latent variables, we usually attain independence between individual observations within a person. The likelihood p(y_i | θ) is then based on Equation (5), where we use conditional independence to multiply across the univariate likelihoods p(y_ij | θ). When we compute information criteria using this conditional likelihood, we approximate the model’s ability to predict new data from the same people in the original dataset because we have conditioned on the latent variables η_i. Merkle et al. (2019) referred to this as “leave one unit out” cross-validation.

Frequentists typically use a marginal likelihood for model estimation, where the latent variables have been integrated out of the model. In this case, the likelihood p(y_i | θ) is based on Equation (6). For IRT models and some others, the marginal criteria require quadrature or related numerical methods because an analytic expression like Equation (6) cannot be obtained. When we compute information criteria using the marginal likelihood, we approximate the model’s ability to predict new data from new people who were not in the original dataset. Merkle et al. (2019) referred to this as “leave one cluster out” cross-validation.

In typical psychometric applications, there are many reasons to prefer marginal information criteria. First, we typically view individuals as exchangeable, and we wish to generalize beyond the people that were observed in our data (see Merkle, 2026b). Marginal information criteria approximate such generalization. Second, Millar (2018) discusses the fact that conditional forms of WAIC may not be asymptotically equivalent to leave-one-out cross-validation because the number of model parameters grows with the sample size. So there is less theoretical justification for conditional information criteria, as compared to marginal information criteria. Finally, Merkle et al. (2019) showed that marginal information criteria have less Monte Carlo error than conditional information criteria because the marginal likelihoods involve many fewer parameters. This means that we can obtain precise values of marginal information criteria using fewer posterior samples, as compared to conditional information criteria.

Despite these advantages, marginal information criteria are not always used in practice. This is especially true because traditional MCMC estimation methods sample the latent variables and rely on a conditional likelihood. So researchers will immediately think of the conditional likelihood as the one to use for information criteria, and some software like JAGS will automatically compute DIC for a specified model (which usually relies on the conditional likelihood). Additionally, the marginal likelihood is often more difficult to compute than the conditional likelihood. This can be frustrating to researchers who estimate a Bayesian model with a conditional likelihood, only to find out that they need further machinery to obtain marginal information criteria.

Comparison of Marginal and Conditional Criteria

If the conditional and marginal information criteria always led to the same conclusions, then their distinction could be ignored. On this point, Merkle et al. (2019) presented theoretical results that the conditional and marginal criteria differ, along with two examples where the criteria selected different psychometric models. The conditional criteria generally favored more complex models than marginal criteria, which makes sense when we consider the types of cross-validation that they represent. For example, we already stated that conditional criteria are related to predicting new data from people that we already observed. From this standpoint, we can afford additional model complexity that tunes parameters to the people in the observed data. In contrast, a marginal criterion needs to guard against overfitting to the people in the observed data because it is related to predicting unseen people.

On Use of Joint Likelihoods

Other researchers sometimes use information criteria that are not fully marginal or conditional. Most related to the current paper, Zhang et al. (2019) recommend a “joint” form of DIC in a two-level IRT model of, e.g., students nested in schools. They do not consider WAIC or PSIS-LOO, and we have not seen joint versions of WAIC or PSIS-LOO used in practice.

Using our notation, the joint DIC involves the likelihood p(y_i, η_i | θ) = p(y_i | η_i, θ)p(η_i | θ), where θ would contain all model parameters except the latent variables η_i. This joint likelihood is typically easier to compute than the marginal likelihood, and it appears to have been proposed because it is convenient (see Celeux et al., 2006). It resembles a one-node quadrature approximation of the marginal likelihood (e.g., Rabe-Hesketh et al., 2005), where the specific node varies for each posterior sample and where the quadrature weights do not sum to 1. Zhang et al. (2019) describe some simulations where a joint DIC selects the true model more often than a marginal DIC.

From our point of view, use of the joint likelihood for information criteria is theoretically problematic because it does not have a clear relationship to cross-validation. That is, it is unclear what it means to hold out a model parameter (the η_i) in addition to y_i when considering a model’s predictive accuracy. Inclusion of the η_i in the likelihood will also lead to a large amount of noise in the resulting information criteria, a point that Zhang et al. (2019) mention in their paper. These considerations lead us to maintain the default recommendation of fully marginal information criteria, because these criteria match the type of generalization in which researchers are usually interested. We further consider these recommendations in the General Discussion.

While our information criterion recommendations differ from those of Zhang et al. (2019), we definitely agree with them that marginal forms of information criteria can be computationally problematic. In addition to long computation times, we describe here two examples where likelihood approximations and model misspecifications can lead to unexpected behavior of the marginal criteria. In the first example, we consider a model of ordinal variables that requires numerical approximation of the marginal likelihood. We show how crude numerical approximations lead to imprecise information criteria. In the second example, we consider two-level SEM and show how information criteria are impacted by heterogeneity. Beyond raising awareness for these specific problems, we hope that these examples make it clear that use of information criteria requires some care.

Model Likelihoods

Traditional MCMC methods for estimating this model involve sampling either the ${\mathit{y}}_i^{\ast }$ or the η_i, both of which are a form of data augmentation. In the former case, we can condition on the sampled ${\mathit{y}}_i^{\ast }$ and pretend like we have continuous data when sampling the other model parameters. In the latter case, we achieve independence of the entries of y_i conditioned on η_i. Then we can evaluate the (conditional) likelihood of the y_i via univariate normal CDFs, which is computationally easy.

Neither of these estimation methods is optimal from an information criterion point of view. This is because the marginal likelihood of the y_i is marginal over both the ${\mathit{y}}_i^{\ast }$ and the η_i. This marginal likelihood can be computed with the CDF of the multivariate normal distribution, but its evaluation is generally too slow to be useful for MCMC estimation. We can instead estimate the model (draw posterior samples) by sampling the ${\mathit{y}}_i^{\ast }$ or the η_i, then compute the fully marginal likelihood of the y_i for each posterior sample after estimation.

Even if the fully marginal likelihood is computed after estimation, it is still computationally expensive. This is especially the case for metrics like WAIC and PSIS-LOO that require pointwise likelihoods. For example, say that we have a dataset of 300 individuals, and we save 1,000 posterior samples for each of three chains. If it takes one second to evaluate the pointwise log-likelihoods for one posterior sample, then it will take 50 minutes to evaluate the log-likelihoods of all 3,000 posterior samples. It is clear that we need to be as economical as possible, while still doing accurate computations. There exist multiple options for approximating the marginal likelihood, and we consider two popular options below.

Illustration

Here, we study how crude approximations of the marginal log-likelihood influence the resulting information criterion values. We fit a model to an ordinal dataset, approximate the pointwise contributions to the marginal log-likelihood in various ways, and compare the resulting information criteria to a gold standard.

Results

It is difficult to make exact statements about the computation times of the marginal likelihood approximations because the code that we used (see Merkle, 2026a) was not optimized for parallelization or for all available computational shortcuts. But to provide a rough idea of timings, each set of marginal likelihood computations took at least 5 minutes per chain on a mid-range laptop, with some taking closer to an hour per chain. GHK sampling was faster than quadrature, likely because the GHK code had better optimizations. Cruder approximations (fewer quadrature points or importance samples) provided speedups of up to 30%.

Figure 1 shows our primary results, with information criteria on the y-axes and number of quadrature points or importance samples on the x-axes. The columns show the adapted quadrature and GHK methods, and the rows show DIC and PSIS-LOO. The dashed line is the gold standard value with 8 quadrature points, and shaded areas represent 1 standard error around the criterion value (note that standard errors are not available for DIC). We see that the quadrature method generally yields information criterion values that are close to the gold standard. For GHK sampling, the criterion values are slightly larger when the approximation is crude. These differences look small, with the dashed line being well within one standard error of the PSIS-LOO values. But information criteria are generally correlated across models, so that differences between model criteria will have much smaller standard errors. For example, it is possible to compare the gold standard PSIS-LOO to other PSIS-LOO metrics, pretending like the metrics were coming from different models. In doing so, we find PSIS-LOO differences of 10 standard errors or more.

Click to enlarge

Figure 1

Information Criteria of Ordinal Factor Analysis Model

Note. Criteria are computed using marginal likelihoods approximated by adapted quadrature and by GHK sampling with differing numbers of nodes/samples. Shaded regions represent ±1 standard error.

Figure 2 is arranged similarly to Figure 1, except it shows the estimated effective number of parameters instead of the information criteria. We see here that, under crude approximations of the marginal likelihood, the effective number of parameters is usually estimated to be too large. The exception is the upper right panel, which shows DIC with GHK sampling. There, effective number of parameters is too low until we obtain 40 samples or more for the approximation.

Click to enlarge

Figure 2

Estimated Effective Numbers of Parameters Associated with DIC and PSIS-LOO

Note. Effective numbers of parameters are computed using marginal likelihoods approximated by adapted quadrature and by GHK sampling with differing numbers of nodes/samples. Shaded regions represent ±1 standard error.

Model

We now define the vector y_ik to be the data of Level 1 unit i within Level 2 unit k. In traditional applications, i would be student and k would be school, and we generally refer to students and schools below because it helps with intuition. Mixed modeling researchers may view this as a three-level dataset where individual observations are nested within students, and students are nested within schools. See Skrondal and Rabe-Hesketh (2004) for further discussion of these different viewpoints.

Specification of a two-level SEM with random intercepts can look very similar to a one-level SEM. We start out writing

which is nearly the same as the traditional model from Equations (1) and (2). The only differences are that, (i) our intercept ν has a k subscript, and (ii) some other vectors have ik subscripts because we are designating the model for student i in school k.

The ν_k term implies that the intercept is unique to each school k, which is what leads us to call this a “two-level” SEM. We complete the model by considering how to handle the ν_k. If we come from a mixed modeling background, the natural thing to do is place a multivariate normal hyperdistribution on the ν_k, i.e.,

This can add many parameters to the model, especially when each individual has many observed variables. For example, if each individual has p = 9 observed variables, then Σ_ν contains 45 free parameters.

Because an unrestricted multivariate normal can lead to many extra parameters, and also because we are already using latent variables, we could elect to place a second SEM on the ν_k. This typically involves specifying school latent variables that are predictive of individual entries of ν_k. We basically repeat our SEM, modeling ν_k instead of y_ik:

where the (c) superscripts distinguish the school-level matrices from the student-level matrices. The mean intercept across schools, ν₀, does not have a superscript because it is similar to the intercept in a traditional, one-level SEM. And while we do not consider it here, the model may include school-level observed variables that do not exist for individual students. In that case, the ν_k vector is augmented with the observed variables from school k (see Rosseel, 2021).

Conditional and Marginal Distributions

The model has student-level latent variables as well as school-level latent variables. Traditional model estimation via MCMC would proceed by sampling all the latent variables, which typically provides conditional independence between observed variables. This leads to a univariate conditional model likelihood that is relatively simple to evaluate, but it also has implications for information criteria. Depending on what latent variables are (and are not) marginalized out of the model, our information criteria approximate out-of-sample prediction at different levels (also see Zhang et al., 2019):

Note that our terminology refers only to marginalization of latent variables and should not be confused with the marginal likelihood used for Bayes factor computations. The latter marginalizes over all model parameters, including the parameters that would be viewed as fixed in a frequentist context.

In most psychometric modeling situations, researchers do not care about the specific students and/or schools in the data. So, for computing information criteria, the fully marginal likelihood would be the recommended default. Exceptions may be studies where researchers are assessing the specific students in the dataset, or school accountability studies where researchers care about the specific schools in the dataset. The former may warrant a fully conditional likelihood, and the latter may warrant a partially marginal likelihood.

Unfortunately, the two-level, fully marginal likelihoods involved in the default recommendation can be computationally difficult. The fully marginal likelihood accounts for the fact that students within a school are correlated, which typically leads to a high-dimensional multivariate normal distribution. As a specific example, imagine that we observe 30 students in each school, and each student provides 5 observed variables. Then our fully marginal likelihood above will involve a 150-dimensional normal distribution. It is difficult to evaluate this high-dimensional likelihood, which leads to excessively slow computations.

Fortunately, a good deal of attention has been devoted to this problem, and a large amount of progress has been made (e.g., du Toit & du Toit, 2008; McDonald, 1993; Rosseel, 2021). This is because the fully marginal likelihood is used for frequentist model estimation, and improved computations can lead to faster frequentist estimation. In the previous work, researchers took advantage of the structure of the model covariance matrix to evaluate the likelihood using matrices of low dimension. Let n_k be the number of students in school k, each of whom provides p observed variables. The previous authors have derived expressions that allow us to evaluate the likelihood using matrices of dimension p, instead of using matrices of dimension n_kp. The appendix includes technical detail on computing the fully marginal, pointwise log-likelihood of each Level-2 unit in a two-level SEM with random intercepts.

Illustration

In the two-level SEM framework, it is possible to observe multiple types of heterogeneity. Here, we consider a situation where the within-cluster covariance matrix differs across clusters. We generate data from a two-level model while manipulating the magnitude of heterogeneity. We then fit a model that ignores this heterogeneity and examine the resulting information criteria.

Method

We generated data from a model with 4 observed variables and 200 clusters of size 50. The data-generating model can be written as:

17

$$\begin{array}{ll}{\mathit{y}}_{ik}\mid {\mathit{\nu }}_k& \sim \{\begin{array}{ll}\text{N}({\mathit{\nu }}_k,{\mathbf{\Sigma }}_w)& i=1,\dots ,50;k=1,\dots ,100\\ \text{N}({\nu }_k,d\times {\mathbf{\Sigma }}_w)& i=1,\dots ,50;k=101,\dots ,200\end{array}\end{array}$$

18

$${\mathit{\nu }}_k\sim \text{Ν}(\mathbf{0},{\mathbf{\Sigma }}_b)k=1,\mathrm{...},200$$

19

$${\mathbf{\Sigma }}_w=\left(\begin{array}{cccc}1.5& 0.3& 0.3& 0.3\\ 0.3& 1.5& 0.3& 0.3\\ 0.3& 0.3& 1.5& 0.3\\ 0.3& 0.3& 0.3& 1.5\end{array}\right)$$

20

$${\mathbf{\Sigma }}_b=\left(\begin{array}{cccc}0.6& 0.2& 0.2& 0.2\\ 0.2& 0.6& 0.2& 0.2\\ 0.2& 0.2& 0.6& 0.2\\ 0.2& 0.2& 0.2& 0.6\end{array}\right)$$

where we manipulated the scalar d to introduce heterogeneity in the within-cluster covariance matrix. The d parameter assumed values from 1 to 4 in steps of 0.5, where a new dataset was generated for each value of d. We then fit a traditional two-level model that assumes a single within-cluster covariance matrix for all clusters. In (b)lavaan, this model could be specified as

model <-'
    level: within
        fw =~ y1 + y2 + y3 + y4
    level: between
        fb =~ y1 + y2 + y3 + y4
'

Click to enlarge

Figure 3

Magnitude of Heterogeneity in Within-Cluster Covariance Matrix (X-Axis) vs. Estimated Effective Number of Parameters (Y-Axis) of DIC and PSIS-LOO.

The prior distributions that we used for model estimation were

$$\begin{array}{ll}{\mathit{\nu }}_0& \sim \text{Normal}(0,32)\\ {\lambda }_j& \sim \text{Normal}(1,.4)j=2,\dots ,4\\ {\lambda }_j^{(c)}& \sim \text{Normal}(1,.4)j=2,\dots ,4\\ \psi & \sim \text{Gamma}(1,.5)\\ {\psi }^{(c)}& \sim \text{Gamma}(1,.5)\\ {\theta }_j& \sim \text{Gamma}(1,.5)j=1,\dots ,4\\ {\theta }_j^{(c)}& \sim \text{Gamma}(1,.5)j=1,\dots ,4.\end{array}$$

where λ₁ = ${\lambda }_1^{(c)}$ = 1 and the ψ and θ parameters are standard deviations, as opposed to variances. The priors on loadings are intended to be mildly informative about the direction of correlation between observed variables (i.e., we expect all observed variables to be positively related). The remaining priors are relatively uninformative, considering the data generating model. For each of three chains, we used blavaan to obtain 2,000 samples following 500 burnin samples.

Computational Issues

Marginal information criteria are generally more difficult to compute, as compared to conditional information criteria, because they often require extra computations after model estimation that may not be needed otherwise. Software tools can facilitate some of these extra computations. These include the blavaan package, which automatically computes marginal information criteria for the models that it is able to estimate, and the bleval package (Luo & Dong, 2025; Luo et al., 2026), which contains a quadrature method that can be flexibly applied to the posterior samples of many Bayesian models. The latter package requires users to specify a function defining the model’s conditional likelihood.

While the above packages make marginal likelihood evaluations easier, the computations will still be intensive and/or infeasible for some models. While this may lead researchers to use conditional or joint likelihoods, we caution that those criteria do not necessarily approximate the generalization that would often be expected of psychometric models. Researchers may instead consider alternative model evaluation methods and computations. For example, the loo package authors and their colleagues recently developed a subsampling approach to computing PSIS-LOO for large datasets (Magnusson et al., 2020; Magnusson et al., 2025). This approach does not require us to compute the likelihood for every observation in the dataset so can lead to faster computations, while also providing estimates of additional uncertainty due to subsampling. Alternatively, it is possible to directly carry out a cross-validation analysis by sequentially holding out data, fitting models, and examining predictions of the held-out data via metrics other than log-likelihoods. This may be advantageous when model estimation is fast and marginal likelihood computation is slow. As a third alternative, there might be an encompassing model that contains all candidate models as special cases. We can use the encompassing model to account for our uncertainty in candidate models, and the encompassing model also facilitates other model comparison metrics such as the Bayes factor (e.g., Wagenmakers et al., 2010). Finally, Bayesian model averaging, stacking, and related methods can also help us address model uncertainty without the need to select a single model (e.g., Kaplan, 2021; Yao et al., 2018). We acknowledge that none of these solutions provides a fast and simple replacement for conditional information criteria that are automatically computed by some software.

Metric Selection

Throughout the paper, we highlighted that a specific information criterion (marginal or otherwise) should be chosen based on the type of cross-validation that each criterion approximates. Other researchers sometimes focus on a criterion’s ability to select the data generating model. In our view, the model’s ability to select the true model is not necessarily relevant in situations where the true model does not exist. This topic was discussed at length in a special issue of Computational Brain & Behavior by authors with diverse points of view (e.g., Gronau & Wagenmakers, 2019; Vehtari et al., 2019). Piironen and Vehtari (2017) relatedly considered the use of different metrics depending on whether or not we expect the true model to be in the set of candidate models (see their Table 1). We also recall the earlier remarks of Browne (2000): “Any attempt to use cross-validation to detect which of a set of models is correct disregards the fundamental intent of the approach. This is to find a model that yields as small an overall discrepancy as is possible given a specific sample size, not to find a model with no error of approximation. In other words, the intent is to find a model that yields a calibration that is worth further examination despite a possibly small sample, not to seek a correct model with an unknown calibration.” (p. 130).

More generally, as pointed out by Navarro (2019), all the information criteria and other metrics mentioned in this paper emphasize quantitative measures of model performance at the neglect of whether any of the models help solve scientific problems. The overall goal of solving scientific problems should not be lost among model selection metrics.

Recommendations

We offer a variety of practical recommendations for researchers who estimate Bayesian models via MCMC and report Bayesian information criteria. First, verify that there are random model parameters (that is, parameters that would be called “random” in a frequentist context), which will often be called random effects, latent variables, or factors. If there are no random parameters, then the issues described in this paper are not applicable. Second, ensure that a large number of posterior samples are taken to reduce noise in the information criteria. We recommend starting with twice the number of samples that one would use for posterior means of individual parameters, as well as monitoring the effective sample size of the log-likelihood values that are used to compute information criteria. One may encounter a speed-accuracy tradeoff in this step, and conditional information criteria will have more noise than marginal.

Following estimation, researchers should verify the type of information criteria that their software supplies. To our knowledge, Mplus and blavaan are the only programs that supply marginal information criteria of psychometric models by default. Most other software will produce conditional information criteria because the computations are easier. The effective number of parameters associated with a criterion often provides a hint about what likelihood is being used: if the effective number of parameters is close to the frequentist count of model parameters, then a marginal likelihood is being used. Fourth, if one has conditional information criteria but does not have a good reason to prefer conditional information criteria, consider whether it is possible to compute marginal information criteria (with subsampling, if necessary), to compute a related cross-validation metric, or to estimate an encompassing model; see the Computational Issues section above. Finally, be transparent in reporting: provide details about the specific types of metrics that were reported, and be careful about using the metrics to select one model in light of the nuances described in this paper.

Conclusions

The popularity of Bayesian information criteria partly stems from the fact that they can be byproducts of MCMC estimation, allowing researchers to compare models without needing to consider post-estimation checks or computations. This simplicity becomes problematic as models increase in complexity, because many researchers do not expect the information criteria to require extra consideration or computations. In this paper, we highlighted some of these issues and made some recommendations for computing information criteria and related cross-validation metrics. The issues might be disappointing to some researchers because they imply that the information criteria require more work and attention than is desired. This could lead researchers to consider whether they are using information criteria because they are simple or because they are the best tools for their work. Because computation and effort are not infinite, researchers might sometimes find that their resources are better spent on posterior checks, assessment of qualitative differences between models, and encompassing models.