2.1.1 LV measurement model

Let $Y_{ij} \in \mathcal {R}$ be the individual i’s response for MV j, and $\mathbf {X}_i=(X_{i1}, \dots , X_{id})^{^{\top }} \in \mathcal {R}^d$ be the d-dimensional LV for individual i. Denote the conditional density of $Y_{ij}|\mathbf {X}_{i}$ by $f_j(y_{ij}|\boldsymbol {x})$ , where $y_{ij}$ and $\boldsymbol {x}=(x_1, \dots , x_d)^{^{\top }}$ indicate the realized values of $Y_{ij}$ and $\mathbf {X}_i$ , respectively. Let $\mathbf {Y}_i=(Y_{i1}, \dots , Y_{im})^{\top }$ be a collection of m MVs from the same individual i. It is typically assumed that $Y_{i1}, \dots , Y_{im}$ are conditionally independent given $\mathbf {X}_i$ , often referred to as the local independence assumption (McDonald, Reference McDonald1999, pp. 255–257). The conditional density of $\mathbf {Y}_i=\boldsymbol {y}_i=(y_{i1}, \dots , y_{im})^{\top }$ given $\mathbf {X}_i=\boldsymbol {x}$ can then be factorized asFootnote ¹

(1) $$ \begin{align} f(\boldsymbol{y}_i|\boldsymbol{x})=\prod_{j=1}^{m}f_j(y_{ij}|\boldsymbol{x}). \end{align} $$

It follows that the marginal density/likelihood of $\mathbf {Y}_i=\boldsymbol {y}_i$ is obtained as the following integral:

(2) $$ \begin{align} f(\boldsymbol{y}_i)=\int \prod_{j=1}^{m}f_j(y_{ij}|\boldsymbol{x})\phi(\boldsymbol{x}) d\boldsymbol{x}, \end{align} $$

in which $\phi $ denotes the population density of $\mathbf {X}_i$ .

Pooling across a sample of n independent and identically distributed (i.i.d.) random vectors of MVs, the sample log-likelihood can be expressed as

(3) $$ \begin{align} l_n(\boldsymbol{\xi}; \mathbf{Y}) = \sum_{i=1}^n \log f(\boldsymbol{y}_i; \boldsymbol{\xi}), \end{align} $$

in which $\mathbf {Y}=(\mathbf {Y}_1, \dots , \mathbf {Y}_n)^{\top } \in \mathcal {R}^{n \times m}$ denotes the data matrix, and the dependency on model parameters, denoted as $\boldsymbol {\xi } \in \mathcal {R}^q$ , is now included in the notation. In parametric LV measurement models where the distribution of LVs and the conditional distribution of MVs given LVs are specified, $\boldsymbol {\xi }$ is often estimated by maximum likelihood (ML). ML estimation seeks for $\boldsymbol {\xi }$ that solves the following estimating equation:

(4) $$ \begin{align} \nabla_{\boldsymbol{\xi}}l_n(\boldsymbol{\xi}; \mathbf{Y}) = \mathbf{0}, \end{align} $$

in which $\nabla _{\boldsymbol {\xi }}l_n(\boldsymbol {\xi }; \mathbf {Y})$ is a $1 \times q$ vector of partial derivatives of $l_n(\boldsymbol {\xi }; \mathbf {Y})$ with respect to $\boldsymbol {\xi }$ . Provided that the Hessian matrix is negative definite, the solution to Equation (4), denoted $\hat {\boldsymbol {\xi }}$ , corresponds to a local maximizer of the sample log-likelihood function. Under suitable regularity conditions and given correct model specification, $\hat {\boldsymbol {\xi }}$ satisfies

(5) $$ \begin{align} \sqrt{n}(\hat{\boldsymbol{\xi}}-\boldsymbol{\xi}) \overset{d}{\to} \mathcal{N}\left(\mathbf{0},\boldsymbol{\mathcal{I}}^{-1}(\boldsymbol{\xi}) \right) \end{align} $$

as $n \to \infty $ , where $\boldsymbol{\mathcal {I}}(\boldsymbol {\xi })=\mathbb {E} [\nabla _{\boldsymbol {\xi }} \log f(\mathbf {Y}_i; \boldsymbol {\xi })^{{}^{\top }} \nabla _{\boldsymbol {\xi }} \log f(\mathbf {Y}_i; \boldsymbol {\xi })]$ denotes the $q \times q$ per-observation Fisher information matrix. By Equation (5), local identifiability of the model parameters is implicitly assumed.

2.1.2 Common factor model

In the current work, we focus on a specific LV measurement model for continuous MVs: the common factor model. We assume that $\mathbf {X}_i \sim \mathcal {N}(\mathbf {0}, \boldsymbol {\Phi })$ , with the population density of $\mathbf {X}_i=\boldsymbol {x}$ expressed by

(6) $$ \begin{align} \phi(\boldsymbol{x})=(2\pi)^{-\frac{d}{2}} \text{det}(\boldsymbol{\Phi})^{-\frac{1}{2}} \exp\left(-\frac{1}{2}\boldsymbol{x}^{^{\top}} \boldsymbol{\Phi}^{-1} \boldsymbol{x} \right). \end{align} $$

We further assume that $Y_{ij}|\mathbf {X}_i=\boldsymbol {x} \sim \mathcal {N}(\nu _j+\boldsymbol {\lambda }_j^{^{\top }} \boldsymbol {x}, \theta _j)$ for $j=1,\dots ,m$ , where $\nu _j$ denotes an intercept, $\boldsymbol {\lambda }_j$ denotes a $d \times 1$ vector of factor loadings, and $\theta _j$ denotes an error variance. With this model specification, the conditional density $f_j(y_{ij}|\boldsymbol {x})$ is given by

(7) $$ \begin{align} f_j(y_{ij}|\boldsymbol{x})=\frac{1}{\sqrt{2\pi \theta_j}}\exp\left(-\frac{1}{2\theta_j}(y_j-(\nu_j+\boldsymbol{\lambda}_j^{^{\top}} \boldsymbol{x}))^2\right). \end{align} $$

The marginal likelihood $f(\boldsymbol {y}_i)$ can then be derived from Equation (2), resulting in a normal density function for $\mathcal {N}(\boldsymbol {\nu }, \boldsymbol {\Lambda } \boldsymbol {\Phi } \boldsymbol {\Lambda }^{\top } + \boldsymbol {\Theta })$ :

(8) $$ \begin{align} f(\boldsymbol{y}_i)=(2\pi)^{-\frac{m}{2}} \text{det}(\boldsymbol{\Lambda} \boldsymbol{\Phi} \boldsymbol{\Lambda}^{\top} + \boldsymbol{\Theta})^{-\frac{1}{2}} \exp\left(-\frac{1}{2}(\boldsymbol{y}_i - \boldsymbol{\nu})^{^{\top}} (\boldsymbol{\Lambda} \boldsymbol{\Phi} \boldsymbol{\Lambda}^{\top} + \boldsymbol{\Theta})^{-1} (\boldsymbol{y}_i - \boldsymbol{\nu}) \right), \end{align} $$

in which $\boldsymbol {\nu }=(\nu _1, \dots , \nu _m)^{^{\top }}$ is an $m \times 1$ vector of intercepts, $\boldsymbol {\Lambda }=(\boldsymbol {\lambda }_1, \dots , \boldsymbol {\lambda }_m)^{^{\top }}$ is an $m \times d$ matrix of factor loadings, and $\boldsymbol {\Theta } = \text {Diag}(\theta _1,\dots ,\theta _m)$ is an $m \times m$ diagonal error covariance matrix. The model parameters $\boldsymbol {\xi }$ now consist of the free elements in $\boldsymbol {\nu }$ , $\boldsymbol {\Lambda }$ , $\boldsymbol {\Phi }$ , and $\boldsymbol {\Theta }$ . The ML estimator $\hat {\boldsymbol {\xi }}$ can be obtained by solving Equation (4) with the normal log-likelihood $l_n(\boldsymbol {\xi }; \mathbf {Y})$ .

2.1.3 Summary quantity and fit assessment

Under this common factor model, our goal now is to develop test statistics that can capture misfit beyond the mean and covariance structures. To achieve this, consider a summary quantity $H(\mathbf {Y}_i; \boldsymbol {\xi })$ for each individual i, which is a real continuous function allowed to depend on the MVs $\mathbf {Y}_i$ and the model parameters $\boldsymbol {\xi }$ . Depending on the context, this function can also incorporate fixed LV values $\boldsymbol {x}$ . Using $H(\mathbf {Y}_i; \boldsymbol {\xi })$ , we obtain its population expectation by

(9) $$ \begin{align} \eta(\boldsymbol{\xi})=\mathbb{E}[H(\mathbf{Y}_i; \boldsymbol{\xi})]=\int H( \boldsymbol{y}_i; \boldsymbol{\xi}) f(\boldsymbol{y}_i; \boldsymbol{\xi})d\boldsymbol{y}_i. \end{align} $$

An empirical estimator of Equation (9) is the sample average

(10) $$ \begin{align} \hat{\eta}(\boldsymbol{\xi})=\frac{1}{n}\sum_{i=1}^n H(\mathbf{Y}_i; \boldsymbol{\xi}). \end{align} $$

The difference $\hat {\eta }(\boldsymbol {\xi })-\eta ( \boldsymbol {\xi })$ then quantifies the discrepancy between the sample average and the model-implied population expectation, forming the basis for assessing model fit. Because model parameters are typically unknown in practice and must be estimated from data, we replace $\boldsymbol {\xi }$ by the ML estimator $\hat{\boldsymbol{\xi}}$ and rely on the residual $\hat {\eta }(\hat {\boldsymbol {\xi }})-\eta ( \hat {\boldsymbol {\xi }})$ as the final basis for constructing fit statistics.

When the same idea is applied to models for categorical data, the difference between sample and population averages of summary quantities is referred to as generalized residuals (Haberman & Sinharay, Reference Haberman and Sinharay2013). Generalized residuals are asymptotically normal under the correct model specification, which has enabled their wide application in constructing various fit statistics in IRT (e.g., Haberman et al., Reference Haberman, Liu and Lee2019; Haberman et al., Reference Haberman, Sinharay and Chon2013; Monroe, Reference Monroe2021; van Rijn et al., Reference van Rijn, Sinharay, Haberman and Johnson2016).

In this article, we extend the existing theory of generalized residuals (Haberman & Sinharay, Reference Haberman and Sinharay2013) to establish asymptotic properties of transformed residuals under more general measurement models. Our proposed framework not only accommodates continuous MVs, but also enables flexible construction of diverse fit test statistics not considered in extant work on generalized residuals. In the following sections, we begin by presenting the most general theory in Section 2.2, construct two different types of test statistics in Section 2.3, and then narrow down to specific examples of fit assessment in Section 2.4.

2.2.1 Asymptotic normality of residuals

Consider a $k \times 1$ vector of summary quantities

(11) $$ \begin{align} \mathbf{H}(\mathbf{Y}_i; \boldsymbol{\xi})=(H_1(\mathbf{Y}_i; \boldsymbol{\xi}), H_2(\mathbf{Y}_i; \boldsymbol{\xi}), \dots, H_k(\mathbf{Y}_i; \boldsymbol{\xi}))^{\top}, \end{align} $$

whose components are formed by different choices of summary quantities introduced in Section 2.1.3. Analogous to Equations (9) and (10), let $\boldsymbol {\eta }(\boldsymbol {\xi })$ and $\hat {\boldsymbol {\eta }}(\boldsymbol {\xi })$ denote the population expectation and sample average of the vector-valued $\mathbf {H}(\mathbf {Y}_i; \boldsymbol {\xi })$ , respectively:

(12) $$ \begin{align}\boldsymbol{\eta}(\boldsymbol{\xi}) &= \mathbb{E}[\mathbf{H}(\mathbf{Y}_i; \boldsymbol{\xi})], \end{align} $$

(13) $$ \begin{align} \hat{\boldsymbol{\eta}}(\boldsymbol{\xi}) &=\frac{1}{n}\sum_{i=1}^n \mathbf{H}(\mathbf{Y}_i; \boldsymbol{\xi}). \end{align} $$

With the ML estimates for model parameters plugged in, the $k \times 1$ difference vector $\hat {\boldsymbol {\eta }}(\hat {\boldsymbol {\xi }})-\boldsymbol {\eta }(\hat {\boldsymbol {\xi }})$ quantifies the discrepancy between the data and the model, evaluated based on $H_{\ell }(\mathbf {Y}_i; \hat {\boldsymbol {\xi }})$ for $\ell =1,\dots ,k$ .

To further generalize, let $\boldsymbol {\varphi }: \mathcal {R}^{k} \to \mathcal {R}^{k'}$ be a twice continuously differentiable transformation function applied to both $\hat {\boldsymbol {\eta }}$ and $\boldsymbol {\eta }$ . This transformation allows us to construct a $k' \times 1$ vector of transformed residuals:

(14) $$ \begin{align} \boldsymbol{e}_{\boldsymbol{\varphi}}=\boldsymbol{\varphi}(\hat{\boldsymbol{\eta}}(\hat{\boldsymbol{\xi}}))-\boldsymbol{\varphi}(\boldsymbol{\eta}(\hat{\boldsymbol{\xi}})), \end{align} $$

which evaluates model-data fit with respect to the population quantities of interest $\boldsymbol {\varphi }(\boldsymbol {\eta }(\boldsymbol {\xi }))$ . This general formulation enables more flexible construction of statistics, including those involving ratios of summary quantities, as will be illustrated in Sections 2.4.2 and 2.4.3.

In the supplementary material, we establish the following asymptotic normality of the transformed residuals:

(15) $$ \begin{align} \sqrt{n}\boldsymbol{e}_{\boldsymbol{\varphi}} \overset{d} \to \mathcal{N}\left(\mathbf{0},\nabla\boldsymbol{\varphi}(\boldsymbol{\eta}(\boldsymbol{\xi})) (-\mathbf{A}( \boldsymbol{\xi})\boldsymbol{\mathcal{I}}^{-1}(\boldsymbol{\xi})\mathbf{A}( \boldsymbol{\xi})^{\top} + \boldsymbol{\Sigma}_{\mathbf{H}}( \boldsymbol{\xi})) \nabla\boldsymbol{\varphi}(\boldsymbol{\eta}(\boldsymbol{\xi}))^{\top} \right), \end{align} $$

in which $\nabla \boldsymbol {\varphi }$ stands for a $k' \times k$ Jacobian matrix of $\boldsymbol {\varphi }$ , $\mathbf {A}(\boldsymbol {\xi })=\mathbb {E}\left [\mathbf {H}(\mathbf {Y}_i; \boldsymbol {\xi }) \nabla _{\boldsymbol {\xi }} \log f(\mathbf {Y}_i; \boldsymbol {\xi }) \right ]$ , and $\boldsymbol {\Sigma }_{\mathbf {H}}(\boldsymbol {\xi })=\mathrm {Cov}[\mathbf {H}( \mathbf {Y}_i; \boldsymbol {\xi })].$ In the special case where the identity transformation $\boldsymbol \varphi = \boldsymbol {\iota }$ is applied, with $\boldsymbol {\iota }(\mathbf {h})=\mathbf {h}$ for any $\mathbf {h}\in \mathcal {R}^{k}$ , Equation (14) reduces to a $k \times 1$ vector of original residuals,

(16) $$ \begin{align} \boldsymbol{e}_{\boldsymbol{\iota}}=\hat{\boldsymbol{\eta}}(\hat{\boldsymbol{\xi}})-\boldsymbol{\eta}(\hat{\boldsymbol{\xi}}). \end{align} $$

In this case, $\nabla \boldsymbol {\varphi }(\boldsymbol {\eta }(\boldsymbol {\xi }))$ in Equation (15) simplifies to the $k \times k$ identity matrix.

2.2.2 Estimation of the asymptotic covariance matrix

In practice, the asymptotic covariance matrix (ACM) of residuals (Equation (15)) needs to be estimated from the data. When the population ACM is substituted with its consistent estimator, the asymptotic normality results remain unchanged by Slutsky's theorem (Bickel & Doksum, Reference Bickel and Doksum2015, p. 512).

To obtain the consistent estimator of the ACM in Equation (15), we replace $\boldsymbol {\xi }$ with the ML estimator $\hat {\boldsymbol {\xi }}$ wherever they appear in the formula. For the population mean and (co)variance in $\mathbf {A}(\hat {\boldsymbol {\xi }})$ and $\boldsymbol {\Sigma }_{\mathbf {H}}(\hat {\boldsymbol {\xi }})$ , a straightforward approach is to use the sample mean and sample covariance matrix, as done by Haberman et al. (Reference Haberman, Sinharay and Chon2013) and Monroe (Reference Monroe2021). However, our pilot simulation study revealed that this approach could result in severely inflated Type I error rates, particularly with small sample sizes. To address this issue, we suggest obtaining more stable estimates for $\mathbf {A}( \hat {\boldsymbol {\xi }})$ and $\boldsymbol {\Sigma }_{\mathbf {H}}( \hat {\boldsymbol {\xi }})$ by computing the MC mean and covariance. This MC approach has been similarly used to approximate the expected Fisher information for IRT models (Monroe, Reference Monroe2019). The steps are outlined as follows:

1. 1. Sample $\tilde {\boldsymbol {y}}_i, i=1,\dots ,M$ , from the marginal distribution of $\mathbf {Y}_i$ with density $f(\boldsymbol {y}_i;\hat {\boldsymbol {\xi }})$ .

2. 2. Use $\tilde {\boldsymbol {y}}_i$ to calculate $\mathbf {H}( \tilde {\boldsymbol {y}}_i; \hat {\boldsymbol {\xi }})$ and $\nabla _{\boldsymbol {\xi }} \log f(\tilde {\boldsymbol {y}}_i; \hat {\boldsymbol {\xi }})$ .

3. 3. Estimate $\mathbf {A}(\hat {\boldsymbol {\xi }})$ and $\boldsymbol {\Sigma }_{\mathbf {H}}(\hat {\boldsymbol {\xi }})$ by

   (17) $$ \begin{align} \hat{\mathbf{A}}(\hat{\boldsymbol{\xi}})&=\frac{1}{M}\sum_{i=1}^M \left[\mathbf{H}(\tilde{\boldsymbol{y}}_i; \hat{\boldsymbol{\xi}}) \nabla_{\boldsymbol{\xi}} \log f(\tilde{\boldsymbol{y}}_i; \hat{\boldsymbol{\xi}})\right], \end{align} $$
   (18) $$ \begin{align} &\hat{\boldsymbol{\Sigma}}_{\mathbf{H}}(\hat{\boldsymbol{\xi}})= \frac{1}{M-1}\sum_{i=1}^M \left[\left(\mathbf{H}(\tilde{\boldsymbol{y}}_i; \hat{\boldsymbol{\xi}}) - \bar{\mathbf{H}}( \hat{\boldsymbol{\xi}})\right)\left(\mathbf{H}(\tilde{\boldsymbol{y}}_i; \hat{\boldsymbol{\xi}}) - \bar{\mathbf{H}}(\hat{\boldsymbol{\xi}})\right)^{\top}\right]. \end{align} $$

In Equation (18), $\bar {\mathbf {H}}(\hat {\boldsymbol {\xi }})=\frac {1}{M}\sum _{i=1}^M \mathbf {H}(\tilde {\boldsymbol {y}}_i; \hat {\boldsymbol {\xi }})$ . As $M \to \infty, \hat {\mathbf {A}}(\hat {\boldsymbol {\xi }}) \text{ and } \hat{\boldsymbol {\Sigma}}_{\mathbf {H}}(\hat {\boldsymbol {\xi }})$ are consistent estimators of $\mathbf {A}(\hat {\boldsymbol {\xi }}) \text{ and } \boldsymbol {\Sigma }_{\mathbf {H}}(\hat {\boldsymbol {\xi }})$ , respectively.

In the subsequent notation for ACMs, a caret “ˆ” is added to indicate that the ACMs are obtained using $\hat {\boldsymbol {\xi }}, \hat {\mathbf {A}},$ and $\hat {\boldsymbol {\Sigma }}_{\mathbf {H}}$ .

2.3.1 z-statistic

Let $e_{\boldsymbol {\varphi },r}$ be the rth component of $\boldsymbol {e}_{\boldsymbol {\varphi }}$ in Equation (14) ( $r=1, \dots , k'$ ). Similarly, let $\hat {\sigma }_{\boldsymbol {\varphi },r}$ be the square root of the rth diagonal element of $\hat {\boldsymbol {\Sigma }}_{\boldsymbol {\varphi }}$ , where the shorthand notation $\boldsymbol {\Sigma }_{\boldsymbol {\varphi }}$ is introduced to denote the ACM in Equation (15):

(19) $$ \begin{align} \boldsymbol{\Sigma}_{\boldsymbol{\varphi}}(\boldsymbol{\xi})=\nabla\boldsymbol{\varphi}(\boldsymbol{\eta}( \boldsymbol{\xi})) (-\mathbf{A}( \boldsymbol{\xi})\boldsymbol{\mathcal{I}}^{-1}(\boldsymbol{\xi})\mathbf{A}( \boldsymbol{\xi})^{\top} + \boldsymbol{\Sigma}_{\mathbf{H}}(\boldsymbol{\xi})) \nabla\boldsymbol{\varphi}(\boldsymbol{\eta}(\boldsymbol{\xi}))^{\top}. \end{align} $$

For a single residual $e_{\boldsymbol {\varphi },r}$ , we standardize it to construct the z-statistic:

(20) $$ \begin{align} z=\frac{e_{\boldsymbol{\varphi},r}}{\hat{\sigma}_{\boldsymbol{\varphi},r} / \sqrt{n}}. \end{align} $$

Under the null hypothesis that the model is correctly specified, the z-statistic converges in distribution to $\mathcal {N}(0, 1)$ . Consequently, values of $|z|$ > 1.96 suggest significant misfit at $\alpha =0.05$ based on $e_{\boldsymbol {\varphi },r}$ .

2.3.2 $\chi ^2$ -statistic

Building on the joint asymptotic multivariate normality of $\boldsymbol {e}_{\boldsymbol {\varphi }}$ established in Equation (15), we construct a quadratic form statistic $n\boldsymbol {e}_{\boldsymbol {\varphi }}^{\top } \mathbf {W} \boldsymbol {e}_{\boldsymbol {\varphi }}$ , where $\mathbf {W}$ is a $k' \times k'$ symmetric weight matrix. This statistic provides a summary of model misfit based on multiple residuals. Such summary information can help practitioners make informed judgments about the overall model misfit of interest. In the remainder of this section, we first review some existing approaches for constructing quadratic form statistics based on normal variables, followed by our newly proposed approach.

When the model is correctly specified, it is known that a quadratic form statistic of normal variates is distributed as a mixture of independent $\chi _1^2$ random variables (Box, Reference Box1954). The correct p-value for this mixture- $\chi ^2$ distribution can be computed, for example, by using the inversion formula given in Imhof (Reference Imhof1961) or by adjusting the statistic based on its mean and variance, approximating a $\chi ^2$ distribution (Satorra & Bentler, Reference Satorra and Bentler1994; Satterthwaite, Reference Satterthwaite1946). However, calculating p-values under a $\chi ^2$ -mixture reference can often be computationally challenging.

To simplify the p-value calculations, the weight matrix $\mathbf {W}$ can be chosen such that the resulting quadratic form statistic is asymptotically distributed as $\chi ^2$ . Specifically, this condition is met when

(21) $$ \begin{align} \boldsymbol{\Sigma}_{\boldsymbol{\varphi}}\mathbf{W}\boldsymbol{\Sigma}_{\boldsymbol{\varphi}}\mathbf{W}\boldsymbol{\Sigma}_{\boldsymbol{\varphi}}=\boldsymbol{\Sigma}_{\boldsymbol{\varphi}}\mathbf{W}\boldsymbol{\Sigma}_{\boldsymbol{\varphi}} \end{align} $$

and $tr(\mathbf {W}\boldsymbol {\Sigma }_{\boldsymbol {\varphi }})=s$ , where s is the degrees of freedom of the resulting $\chi ^2$ distribution (Schott, Reference Schott2016, Theorem 11.11). A common way to satisfy these conditions is to use the Moore–Penrose pseudoinverse of $\boldsymbol {\Sigma }_{\boldsymbol {\varphi }}$ , $\boldsymbol {\Sigma }_{\boldsymbol {\varphi }}^{+},$ as the weight matrix (Magnus & Neudecker, Reference Magnus and Neudecker1999, p. 36). This approach has been employed in previous studies to construct quadratic form fit statistics (e.g., Liu & Maydeu-Olivares, Reference Liu and Maydeu-Olivares2014; Liu et al., Reference Liu, Yang and Maydeu-Olivares2019; Maydeu-Olivares & Liu, Reference Maydeu-Olivares and Liu2015; Reiser, Reference Reiser1996). However, the use of $\boldsymbol {\Sigma }_{\boldsymbol {\varphi }}^{+}$ can introduce numerical instability, as the rank of $\boldsymbol {\Sigma }_{\boldsymbol {\varphi }}$ needs to be estimated numerically by counting non-zero eigenvalues in its sample estimate. This procedure relies on an arbitrarily selected tolerance for identifying zero eigenvalues, which may lead to inconsistent performance of the resulting test statistic (Maydeu-Olivares & Joe, Reference Maydeu-Olivares and Joe2008).

To address this issue, we propose defining the weight matrix as

(22) $$ \begin{align} \mathbf{W}=\boldsymbol{\Sigma}_{\boldsymbol{\varphi}}^{-}=\mathbf{U} \boldsymbol{\Omega}^{-} \mathbf{U}^{\top}, \end{align} $$

in which $\mathbf {U}$ is the $k' \times k'$ matrix of eigenvectors from the eigendecomposition of $\boldsymbol {\Sigma }_{\boldsymbol {\varphi }}$ (i.e., $\boldsymbol {\Sigma }_{\boldsymbol {\varphi }}=\mathbf {U} \boldsymbol {\Omega } \mathbf {U}^{\top }$ ), $\boldsymbol {\Omega }$ is the diagonal matrix of eigenvalues, and $\boldsymbol {\Omega }^{-}=\text {Diag}(\omega _1^{-1}, \dots , \omega _s^{-1}, 0, \dots , 0)$ . Here, $\boldsymbol {\Omega }^{-}$ stands for the inverse of $\boldsymbol {\Omega }$ only with the largest s eigenvalues, $\omega _1, \dots , \omega _s$ , where $s \leq \text {rank}(\boldsymbol {\Sigma }_{\boldsymbol {\varphi }})$ . The proposed weight matrix $\boldsymbol {\Sigma }_{\boldsymbol {\varphi }}^{-}$ generalizes the pseudoinverse $\boldsymbol {\Sigma }_{\boldsymbol {\varphi }}^{+}$ , subsuming it as a special case when $s=\text {rank}(\boldsymbol {\Sigma }_{\boldsymbol {\varphi }})$ . For $s \leq \text {rank}(\boldsymbol {\Sigma }_{\boldsymbol {\varphi }})$ , $\boldsymbol {\Sigma }_{\boldsymbol {\varphi }}$ serves as a generalized inverse of $\boldsymbol {\Sigma }_{\boldsymbol {\varphi }}^{-}$ , such that $\boldsymbol {\Sigma }_{\boldsymbol {\varphi }}^{-} \boldsymbol {\Sigma }_{\boldsymbol {\varphi }} \boldsymbol {\Sigma }_{\boldsymbol {\varphi }}^{-}=\boldsymbol {\Sigma }_{\boldsymbol {\varphi }}^{-}$ holds. With this formulation, the conditions in Equation (21) and $tr(\boldsymbol {\Sigma }_{\boldsymbol {\varphi }}^{-}\boldsymbol {\Sigma }_{\boldsymbol {\varphi }})=s$ are satisfied.

Accordingly, we define our $\chi ^2$ -statistic using $\hat {\mathbf {W}}=\boldsymbol {\hat {\Sigma }}_{\boldsymbol {\varphi }}^{-}$ as

(23) $$ \begin{align} T=n\boldsymbol{e}_{\boldsymbol{\varphi}}^{\top} \hat{\boldsymbol{\Sigma}}_{\boldsymbol{\varphi}}^{-} \boldsymbol{e}_{\boldsymbol{\varphi}}, \end{align} $$

such that T asymptotically follows a $\chi ^2$ distribution with s degrees of freedom ( $\chi ^2_s)$ under correct model specification. Since this asymptotic property holds for any $1 \leq s \leq \text {rank}(\boldsymbol {\Sigma }_{\boldsymbol {\varphi }})$ , a sufficiently small value, such as $s=1$ , can be chosen as the reference degrees of freedom to avoid numerical instability.

With the choice of $s=1$ , the reference distribution becomes $\chi ^2_1$ , where $T> 3.84$ indicates significant overall misfit at $\alpha =0.05$ based on multiple residuals in $\boldsymbol {e}_{\boldsymbol {\varphi }}$ . The performance of the proposed $\chi ^2$ -statistic with $s=1$ will be evaluated through our simulation studies in Section 3.

2.4.1 Testing normality of LV density

Let $\boldsymbol {x}_1, \dots , \boldsymbol {x}_Q$ denote the Q distinct values of the d-dimensional LV, for which we aim to assess the misfit in LV density. Given our goal of evaluating normality assumption of LVs, the dependency of summary quantities on LV values is now naturally introduced in our formulation, even when not explicitly expressed in the notation.

First, construct $\mathbf {H}$ in Equation (11) with $k=Q$ as

(24) $$ \begin{align} \mathbf{H}(\mathbf{Y}_i; \boldsymbol{\xi})=(H_1(\mathbf{Y}_i; \boldsymbol{\xi}), H_2(\mathbf{Y}_i; \boldsymbol{\xi}), \dots, H_Q(\mathbf{Y}_i; \boldsymbol{\xi}))^{\top}, \end{align} $$

in which each component $H_\ell $ ( $\ell =1, \dots , Q)$ represents the posterior density of $\boldsymbol {x}_\ell $ given $\mathbf {Y}_i$ :

(25) $$ \begin{align} H_\ell(\mathbf{Y}_i; \boldsymbol{\xi})=f(\boldsymbol{x}_\ell|\mathbf{Y}_i; \boldsymbol{\xi}). \end{align} $$

With this choice of $H_\ell $ , its population expectation, denoted by $\eta _\ell $ , corresponds to the LV density $\phi (\boldsymbol {x}_\ell; \boldsymbol {\xi })$ defined in Equation (6):

(26) $$ \begin{align} \eta_\ell(\boldsymbol{\xi}) &= \mathbb{E} [H_\ell(\mathbf{Y}_i; \boldsymbol{\xi})]=\phi(\boldsymbol{x}_\ell; \boldsymbol{\xi}). \end{align} $$

Equation (26) follows from Bayes’ theorem,

(27) $$ \begin{align} f(\boldsymbol{x}_\ell|\mathbf{Y}_i; \boldsymbol{\xi})=\frac{\phi(\boldsymbol{x}_\ell; \boldsymbol{\xi})f(\mathbf{Y}_{i}|\boldsymbol{x}_\ell; \boldsymbol{\xi})}{f(\mathbf{Y}_i; \boldsymbol{\xi})}, \end{align} $$

along with Equation (9). The empirical estimator of $\eta _\ell $ , denoted by $\hat {\eta }_\ell $ , is then constructed as

(28) $$ \begin{align} \hat{\eta}_\ell(\boldsymbol{\xi}) &= \frac{1}{n}\sum_{i=1}^n H_\ell(\mathbf{Y}_i; \boldsymbol{\xi})=\frac{1}{n}\sum_{i=1}^n f(\boldsymbol{x}_\ell|\mathbf{Y}_i; \boldsymbol{\xi}), \end{align} $$

representing the sample mean of the posterior densities. Now in the vector form, the difference between $\hat {\boldsymbol {\eta }}$ and $\boldsymbol {\eta }$ (Equations (13) and (12)) collectively represents the discrepancies between the population LV density and its sample estimate on $\boldsymbol {x}_1, \dots , \boldsymbol {x}_Q$ . Thus, with $\boldsymbol {\varphi }=\boldsymbol {\iota }$ being the identity transformation (see Equation (16)), the resulting $Q \times 1$ vector of residuals,

(29) $$ \begin{align} \boldsymbol{e}_{\boldsymbol{\iota}}=\hat{\boldsymbol{\eta}}(\hat{\boldsymbol{\xi}})-\boldsymbol{\eta}(\hat{\boldsymbol{\xi}}), \end{align} $$

reflects the nonnormality of the LVs.

With this formulation, we construct the z-statistic from Equation (20) based on each single residual $e_{\boldsymbol {\iota }, \ell }=\hat {\eta }_\ell (\hat {\boldsymbol {\xi }})-\eta _\ell (\hat {\boldsymbol {\xi }})$ , which captures the misfit in LV density at $\boldsymbol {x}_\ell $ . In addition, we construct the $\chi ^2$ -statistic from Equation (23) based on $\boldsymbol {e}_{\boldsymbol {\iota }}$ , which reflects the overall misfit in LV density across $\boldsymbol {x}_1, \dots , \boldsymbol {x}_Q$ . Given the nature of fit assessment in our example, we will also use the terms pointwise statistic and summary statistic for the z- and $\chi ^2$ -statistics, respectively. Also, we denote these statistics developed for this LV density example as $z_1(\boldsymbol {x}_\ell )$ and $T_1$ , respectively, for use in later sections. The pointwise statistic $z_1$ coincides with the one proposed in the context of IRT for detecting misfit in the LV density under unidimensional IRT models (Monroe, Reference Monroe2021).

2.4.2 Testing MV-level linearity

On the same grid of LV values $\boldsymbol {x}_1, \dots , \boldsymbol {x}_Q$ , construct $\mathbf {H}$ in Equation (11) with $k=2Q$ as follows:

(30) $$ \begin{align} \mathbf{H}(\mathbf{Y}_i; \boldsymbol{\xi})=(H_{1}(\mathbf{Y}_i; \boldsymbol{\xi}), \dots, H_{Q}(\mathbf{Y}_i; \boldsymbol{\xi}), H_{Q+1}(\mathbf{Y}_i; \boldsymbol{\xi}), \dots, H_{2Q}(\mathbf{Y}_i; \boldsymbol{\xi}))^{\top}, \end{align} $$

in which

(31) $$ \begin{align} H_{\ell}(\mathbf{Y}_i; \boldsymbol{\xi})&=Y_{ij} f(\boldsymbol{x}_\ell|\mathbf{Y}_i; \boldsymbol{\xi}), \end{align} $$

(32) $$ \begin{align} &H_{Q+\ell}(\mathbf{Y}_i; \boldsymbol{\xi})=f(\boldsymbol{x}_\ell|\mathbf{Y}_i; \boldsymbol{\xi}), \end{align} $$

for $\ell =1, \dots , Q$ .

Following the construction of $\hat {\boldsymbol {\eta }}$ and $\boldsymbol {\eta }$ as in Equations (13) and (12), define the $Q \times 1$ vector of transformed residuals (i.e., $k'=Q$ ) as

(33) $$ \begin{align} \boldsymbol{e}_{\boldsymbol{\varphi}}=\boldsymbol{\varphi}(\hat{\boldsymbol{\eta}}(\hat{\boldsymbol{\xi}}))-\boldsymbol{\varphi}(\boldsymbol{\eta}(\hat{\boldsymbol{\xi}})), \end{align} $$

in which a transformation function $\boldsymbol {\varphi }: \mathcal {R}^{2Q} \to \mathcal {R}^Q$ is applied to produce ratios:

(34) $$ \begin{align} \boldsymbol{\varphi}(\boldsymbol{\gamma})=\left(\frac{\gamma_1}{\gamma_{Q+1}}, \dots, \frac{\gamma_Q}{\gamma_{2Q}}\right)^{\top} \end{align} $$

with $\boldsymbol {\gamma }=(\gamma _1,\dots ,\gamma _Q,\gamma _{Q+1},\dots ,\gamma _{2Q})^{\top }$ . With this setup, the residuals in Equation (33) represent the misfit in the mean response for the jth MV, conditional on $\boldsymbol {x}_1,\dots ,\boldsymbol {x}_Q$ .

To elaborate, the expectations of $H_{\ell }$ and $H_{Q+\ell }$ are obtained by

(35) $$ \begin{align} &\eta_{\ell}(\boldsymbol{\xi}) = \mathbb{E} [H_{\ell}(\mathbf{Y}_i; \boldsymbol{\xi})]=\phi(\boldsymbol{x}_\ell; \boldsymbol{\xi})\mathbb{E}(Y_{ij} | \boldsymbol{x}_\ell; \boldsymbol{\xi}), \end{align} $$

(36) $$ \begin{align} &\eta_{Q+\ell}(\boldsymbol{\xi}) = \mathbb{E} [H_{Q+\ell}(\mathbf{Y}_i; \boldsymbol{\xi})]=\phi(\boldsymbol{x}_\ell; \boldsymbol{\xi}) \end{align} $$

from Equations (9) and (27), with corresponding empirical estimators

(37) $$ \begin{align} &\hat{\eta}_{\ell}(\boldsymbol{\xi}) = \frac{1}{n}\sum_{i=1}^n H_{\ell}(\mathbf{Y}_i; \boldsymbol{\xi})=\frac{1}{n}\sum_{i=1}^n Y_{ij}f(\boldsymbol{x}_\ell|\mathbf{Y}_i; \boldsymbol{\xi}), \end{align} $$

(38) $$ \begin{align} &\hat{\eta}_{Q+\ell}(\boldsymbol{\xi}) = \frac{1}{n}\sum_{i=1}^n H_{Q+\ell}(\mathbf{Y}_i; \boldsymbol{\xi})=\frac{1}{n}\sum_{i=1}^n f(\boldsymbol{x}_\ell|\mathbf{Y}_i; \boldsymbol{\xi}). \end{align} $$

After applying the ratio transformation function in Equation (34) to the vector form counterparts $\hat {\boldsymbol {\eta }}$ and $\boldsymbol {\eta }$ , the $\ell $ th component of $\boldsymbol {e}_{\boldsymbol {\varphi }}$ in Equation (33) becomes

(39) $$ \begin{align} e_{\boldsymbol{\varphi}, \ell} = \frac{\hat{\eta}_{\ell}(\hat{\boldsymbol{\xi}})}{\hat{\eta}_{Q+\ell}(\hat{\boldsymbol{\xi}})} - \frac{\eta_{\ell}(\hat{\boldsymbol{\xi}})}{\eta_{Q+\ell}(\hat{\boldsymbol{\xi}})} = \frac{\sum_{i=1}^n Y_{ij}f(\boldsymbol{x}_\ell|\mathbf{Y}_i; \hat{\boldsymbol{\xi}})}{\sum_{i=1}^n f(\boldsymbol{x}_\ell|\mathbf{Y}_i; \hat{\boldsymbol{\xi}})} - \mathbb{E}(Y_{ij} | \boldsymbol{x}_\ell; \hat{\boldsymbol{\xi}}). \end{align} $$

Note that the conditional expectation in Equation (39), $\mathbb {E}(Y_{ij} | \boldsymbol {x}_\ell; \hat {\boldsymbol {\xi }})$ , is known as the item characteristic curve under a unidimensional IRT model with dichotomous response variables. Based on this, similar form of residuals has been proposed to assess item-level fit in IRT models (e.g., Haberman et al., Reference Haberman, Sinharay and Chon2013). In the context of the common factor model, the conditional expectation is a linear function of $\boldsymbol {x}_\ell $ :

(40) $$ \begin{align} \mathbb{E}(Y_{ij} | \boldsymbol{x}_\ell; \hat{\boldsymbol{\xi}})=\hat{\nu}_j+\boldsymbol{\hat{\lambda}}_j^{^{\top}} \boldsymbol{x}_\ell. \end{align} $$

Therefore, the residual $e_{\boldsymbol {\varphi }, \ell }$ in Equation (39) serves as a basis for assessing this linearity assumption for the jth MV.

With this background, the pointwise statistic is now constructed from Equation (20) to assess the MV-level linearity at $\boldsymbol {x}_\ell $ , while the summary statistic is constructed from Equation (23) to provide an overall assessment across multiple LV grids $\boldsymbol {x}_1, \dots , \boldsymbol {x}_Q$ . We refer to these statistics as $z_2(\boldsymbol {x}_\ell )$ and $T_2$ , respectively.

Before moving on to the next example, it is worth noting that alternative ways to formulate residuals targeting the same population quantity may exist. For instance, with the target population quantity $\mathbb {E}(Y_{ij} | \boldsymbol {x}_\ell; \boldsymbol {\xi })$ as in this example, we can keep the formulation of $\mathbf {H}$ as in Equation (24), with each component defined as

(41) $$ \begin{align} H_\ell(\mathbf{Y}_i; \boldsymbol{\xi})=\frac{Y_{ij}f(\mathbf{Y}_{i}|\boldsymbol{x}_\ell; \boldsymbol{\xi})}{f(\mathbf{Y}_i; \boldsymbol{\xi})}. \end{align} $$

With this choice of $H_\ell $ , its population expectation $\eta _\ell $ directly becomes $\mathbb {E}(Y_{ij} | \boldsymbol {x}_\ell; \boldsymbol {\xi })$ , eliminating the need for a transformation and thus simplifying the formulation of residuals. The corresponding empirical estimator is

(42) $$ \begin{align} \hat{\eta}_\ell(\boldsymbol{\xi}) = \frac{1}{n}\sum_{i=1}^n \frac{Y_{ij}f(\mathbf{Y}_i|\boldsymbol{x}_\ell; \boldsymbol{\xi})}{f(\mathbf{Y}_i; \boldsymbol{\xi})}=\frac{1}{\phi(\boldsymbol{x}_\ell; \boldsymbol{\xi})} \frac{1}{n} \sum_{i=1}^n Y_{ij}f(\boldsymbol{x}_\ell|\mathbf{Y}_i; \boldsymbol{\xi}). \end{align} $$

We observed in pilot simulations that the performance of generalized residuals based on Equation (42) is very sensitive to misspecification of the LV density $\phi (\boldsymbol {x}_\ell; \boldsymbol {\xi })$ . In contrast, the ratio formulation we presented in Equation (39) also estimates the LV density in Equation (42), resulting in a more targeted fit diagnostic. Further discussion on this is provided in the supplementary material using a real-data example.

2.4.3 Testing MV-level homoscedasticity

In this example, we retain the setup from Section 2.4.2, modifying only the first $\mathcal{Q}$ elements of $\mathbf{H}(\mathbf{Y}_i; \boldsymbol{\xi})$ defined in Equation (31) as follows:

(43) $$\begin{align} H_{\ell}(\mathbf{Y}_i; \boldsymbol{\xi})&=[Y_{ij}-\mathbb{E}(Y_{ij}|\boldsymbol{x}_\ell; \boldsymbol{\xi})]^2 f(\boldsymbol{x}_\ell|\mathbf{Y}_i; \boldsymbol{\xi}), \quad \ell=1,\dots,Q. \end{align}$$

The expectation of $H_{\ell }$ then becomes

(44) $$ \begin{align} \eta_{\ell}(\boldsymbol{\xi}) &= \mathbb{E} [H_{\ell}(\mathbf{Y}_i; \boldsymbol{\xi})]=\phi(\boldsymbol{x}_\ell; \boldsymbol{\xi})\mathrm{Var}(Y_{ij} | \boldsymbol{x}_\ell; \boldsymbol{\xi}) \end{align} $$

with the empirical estimator

(45) $$ \begin{align} \hat{\eta}_{\ell}(\boldsymbol{\xi}) =\frac{1}{n}\sum_{i=1}^n [Y_{ij}-\mathbb{E}(Y_{ij}|\boldsymbol{x}_\ell; \boldsymbol{\xi})]^2 f(\boldsymbol{x}_\ell|\mathbf{Y}_i; \boldsymbol{\xi}). \end{align} $$

The definitions of $H_{Q+\ell }, \eta _{Q+\ell },$ and $\hat \eta _{Q+\ell }$ remain unchanged, as specified in Equations (32), (36), and (38). After applying the ratio transformation function in Equation (34) to the corresponding vector forms $\boldsymbol {\eta }$ and $\boldsymbol {\hat \eta }$ , the $\ell $ th component of the resulting $\boldsymbol {e}_{\boldsymbol {\varphi }}$ becomes

(46) $$ \begin{align} e_{\boldsymbol{\varphi}, \ell} = \frac{\hat{\eta}_{\ell}(\hat{\boldsymbol{\xi}})}{\hat{\eta}_{Q+\ell}(\hat{\boldsymbol{\xi}})} - \frac{\eta_{\ell}(\hat{\boldsymbol{\xi}})}{\eta_{Q+\ell}(\hat{\boldsymbol{\xi}})} = \frac{\sum_{i=1}^n [Y_{ij}-\mathbb{E}(Y_{ij}|\boldsymbol{x}_\ell;\hat{\boldsymbol{\xi}})]^2 f(\boldsymbol{x}_\ell|\mathbf{Y}_i; \hat{\boldsymbol{\xi}})}{\sum_{i=1}^n f(\boldsymbol{x}_\ell|\mathbf{Y}_i; \hat{\boldsymbol{\xi}})} - \mathrm{Var}(Y_{ij} | \boldsymbol{x}_\ell; \hat{\boldsymbol{\xi}}). \end{align} $$

Under the common factor model, the conditional variance of the jth MV in Equation (46) is a constant (i.e., $\mathrm {Var}(Y_{ij} | \boldsymbol {x}_\ell; \hat {\boldsymbol {\xi }})=\hat {\theta }_j$ ) and does not depend on LVs $\boldsymbol {x}_\ell $ . Therefore, this residual can be used to assess the constant-variance assumption for the jth MV. Similar to the previous examples, the pointwise statistic is constructed from Equation (20) to assess the MV-level homoscedasticity at $\boldsymbol {x}_\ell $ , while the summary statistic is constructed from Equation (23) to summarize results across multiple LV grids $\boldsymbol {x}_1, \dots , \boldsymbol {x}_Q$ . We refer to these statistics as $z_3(\boldsymbol {x}_\ell )$ and $T_3$ , respectively.

In summary, we have discussed three examples of fit assessments that can be developed within our framework. We constructed pointwise statistics, $z_1(\boldsymbol {x}_\ell ), z_2(\boldsymbol {x}_\ell ),$ and $z_3(\boldsymbol {x}_\ell )$ ( $\ell =1, \dots , Q$ ), and summary statistics, $T_1, T_2,$ and $T_3$ , optimized for testing the normality of LV density, MV-level linearity of the mean function, and MV-level homoscedasticity of the variance function, respectively. Table 1 provides a summary of the three examples discussed so far.

Table 1 Summary of components for constructing the three examples discussed in Section 2.4

Note: For brevity, dependency on model parameters $\boldsymbol {\xi }$ is omitted. In the first column of the second row, $[\boldsymbol {\varphi }(\boldsymbol {\eta })]_\ell $ denotes the $\ell $ th component of the transformed $\boldsymbol {\eta }$ , representing the population quantity of interest. In these examples, the index $\ell $ distinguishes the LV values under evaluation, where $\ell = 1, \dots , Q.$

3.1.1 Simulation setup

As a special case of the common factor model presented in Section 2.1.2, we considered a two-dimensional independent-cluster model with LVs: $\mathbf {X}_i=(X_{i1}, X_{i2})^{^{\top }}$ . The number of MVs was set to 20 (i.e., $m=20$ ), with the first ten MVs loaded only on $X_{i1}$ , and the remaining ten only on $X_{i2}$ . The two LVs were assumed to be correlated with each other.

For data generation, two different LV densities were specified to represent correctly specified and misspecified conditions. For the correctly specified condition, $\mathbf {X}_i, i=1, \dots , n,$ were generated from the standard bivariate normal distribution with a correlation coefficient of 0.2. For the misspecified condition, we simulated the LVs from a mixture of two bivariate normal distributions: $\mathbf {X}_i \sim 0.5\mathcal {N}(\boldsymbol {\mu }^{(1)}, \boldsymbol {\Phi }^{(1)})+0.5\mathcal {N}(\boldsymbol {\mu }^{(2)}, \boldsymbol {\Phi }^{(2)}),$ in which $\boldsymbol {\mu }^{(1)}=(-0.6, -0.6)^{\top }$ , $\boldsymbol {\mu }^{(2)}=(0.6, 0.6)^{\top }$ , and

(47) $$ \begin{align} \boldsymbol{\Phi}^{(1)}=\boldsymbol{\Phi}^{(2)}=\begin{bmatrix} 0.8^2 & -0.16 \\ -0.16 & 0.8^2 \end{bmatrix}. \end{align} $$

With this setup, the LV distributions in both conditions have means of zero and the covariance matrix of

(48) $$ \begin{align} \boldsymbol{\Phi}=\mathrm{Cov}(\mathbf{X}_i)= \begin{bmatrix} 1 & 0.2 \\ 0.2 & 1 \end{bmatrix}. \end{align} $$

The left panel of Figure 1 depicts the contour plot of the two data-generating LV densities, with one superimposed on the other. The right panel displays the conditional density of $X_{i1}=x_1$ given the mean of $X_{i2}$ , which is zero in this case. The right panel will be revisited in the discussion of results in Section 3.1.2.

Figure 1 Left panel: Contour plot of the data-generating LV densities under correctly specified (gray solid lines) and misspecified (black dashed lines) conditions. Note: Right panel: Conditional density of $X_{i1}=x_1$ given $X_{i2}=x_2=0$ under the same conditions. The conditional density of $X_{i2}=x_2$ given $X_{i1}=x_1=0$ is identical to that shown in the right panel.

In both correctly and incorrectly specified conditions, other data-generating parameters were set equal. In particular, MVs were assumed to have zero intercepts and exhibit an independent-cluster factor structure: $\boldsymbol {\nu }=\mathbf {0}$ and

(49) $$ \begin{align} \boldsymbol{\Lambda}= \begin{bmatrix} \boldsymbol{a}_1 & \textbf{0}_{10} \\ \textbf{0}_{10} & \boldsymbol{a}_2\\ \end{bmatrix}, \end{align} $$

in which $\boldsymbol {a}_1=\boldsymbol {a}_2=(\sqrt {0.3},\sqrt {0.5},\sqrt {0.7},\sqrt {0.3},\sqrt {0.5},\sqrt {0.7},\sqrt {0.3},\sqrt {0.5},\sqrt {0.7},\sqrt {0.3})^{\top }$ and $\textbf{0}_{10}$ indicates a $10 \times 1$ vector of zeros. Error variances $\theta _j, j=1, \dots , 20,$ were then determined by diagonal elements of the square matrix $\mathbf {I}-\boldsymbol {\Lambda } \boldsymbol {\Phi } \boldsymbol {\Lambda }^{\top }$ . This model specification implies that the communalities alternate among 0.3, 0.5, and 0.7, covering a range of values typically observed in practice. A value of 0.3 represents low communality, while 0.7 represents high communality (MacCallum et al., Reference MacCallum, Widaman, Zhang and Hong1999). The value of 0.5 was chosen as an intermediate point to reflect moderate communality.

Upon generating data using the specified parameter values, the proposed tests for assessing the LV density fit were conducted using the following procedures: (1) estimate the model parameters $\boldsymbol {\xi }$ by ML to get $\hat {\boldsymbol {\xi }}$ , (2) construct $\boldsymbol {\hat {\eta }}(\hat {\boldsymbol {\xi }})$ and $\boldsymbol {\eta }(\hat {\boldsymbol {\xi }})$ and calculate residuals by $\boldsymbol {e}_{\boldsymbol {\iota }}=\boldsymbol {\hat {\eta }}(\hat {\boldsymbol {\xi }})-\boldsymbol {\eta }(\hat {\boldsymbol {\xi }})$ , (3) estimate the ACM of $\boldsymbol {e}_{\boldsymbol {\iota }}$ , and (4) compute $z_1(\boldsymbol {x}_\ell )$ for $\ell =1,\dots , Q$ and $T_1$ to conduct pointwise z-tests and the overall $\chi ^2$ -test.

For parameter estimation in Step 1, the lavaan package version 0.6-14 (Rosseel, Reference Rosseel2012) was used with the default option of ML estimation with the normal likelihood. Mean structure was added to the model (meanstructure = TRUE). To check model convergence, we used the following criteria: (1) the absence of Heywood cases, (2) the maximum absolute element of the gradient vector is less than $10^{-4}$ , and (3) the minimum eigenvalue of the Hessian matrix is greater than zero. Default maximum number of iterations was used. In Step 3, the estimated ACM was obtained following the process described in Section 2.2.2 with $M=10,000$ . The inverse of the expected information matrix was obtained using the inverted.information option in the lavInspect function, and the gradient of the individual log-likelihood, $\nabla _{\boldsymbol {\xi }} \log f(\mathbf {Y}_i; \hat {\boldsymbol {\xi }})$ , was obtained using the lavScores function. In Step 4, $19$ evenly spaced values between $-3$ and $3$ were used as evaluation points per LV, resulting in a total of $19^2=361$ pointwise statistics (i.e., $Q=361$ ) for evaluating the fit at every possible combination of LVs. For the summary statistic, a subgrid of $7^2=49$ points, obtained from seven equally spaced values from $-2$ to $2$ per LV, were used for illustration purposes. In addition, we set the degree of freedom $s=1$ in Equation (23) so that the reference limiting distribution for $T_1$ becomes $\chi _1^2$ .

3.2.1 Simulation setup

In Study 2, the linear normal one-factor model with $X_i \sim \mathcal {N}(0, 1)$ was under investigation with or without the misspecification in MV-level mean/variance functions. Let $\mu _j$ and $\sigma _j^2$ denote the mean and variance of the conditional distribution $Y_{ij}|X_i$ . Under the linear normal factor model, $Y_{ij}|X_i=x \sim \mathcal {N}(\mu _j(x), \sigma _j^2(x))$ with the mean function $\mu _j(x)=\nu _j+\lambda _j x$ and variance function $\sigma _j^2(x)=\theta _j$ . MV-level fit assessments were performed by investigating the linearity of $\mu _j$ and homoscedasticity of $\sigma _j^2$ .

The data-generating model consisted of four types of MVs summarized in Table 3. The four types were determined by fully crossing linear versus quadratic mean functions and constant versus log-linear variance functions, labeled as LMCV, QMCV, LMLV, and QMLV, respectively. Regardless of the type of MVs, the intercept was set to zero, i.e., $\nu _j=0$ , and the error variance was obtained by $\theta _j=1-\lambda _j^2$ . The factor loading $\lambda _j$ was set to alternate between $\sqrt {0.3}, \sqrt {0.5},$ and $\sqrt {0.7}$ for LMCVs; $\lambda _j$ was fixed at $\sqrt {0.5}$ for the other three types of MVs. After that, coefficients $\kappa _j$ , $\gamma _{0j}$ , and $\gamma _{1j}$ were added for the misspecified MVs. The quadratic coefficient $\kappa _j$ for QMCV and QMLV was fixed at $-0.1$ ; $\gamma _{0j}$ and $\gamma _{1j}$ for LMLV and QMLV were set as $\gamma _{0j}=\log \theta _j$ and $\gamma _{1j}=0.3$ , respectively. In Figure 4, the black dashed and solid lines illustrate the shapes of $\mu _j(x)$ and $\sigma _j^2(x)$ with and without these additional coefficients, respectively, the combination of which comprises the four types of MVs. As an illustration, the values of $\nu _j=0$ , $\lambda _j=\sqrt {0.5}$ , and $\theta _j=0.5$ were used in generating the figure.

Table 3 Four types of MVs used for data generation

Note: LMCV: linear mean function, constant variance function. QMCV: quadratic mean function, constant variance function. LMLV: linear mean function, log-linear variance function. QMLV: quadratic mean function, log-linear variance function.

Figure 4 Black lines illustrate the mean function $\mu _j(x)$ (left panel) and variance function $\sigma _j^2(x)$ (right panel) under correctly specified (LM, CV) and misspecified (QM, LV) conditions in the data-generating model.

Note: Gray dotted lines (LMQ, CVL) approximate the expected functions under misspecification: a misfitting linear mean under a quadratic population (left) and constant variance under a log-linear population (right). LM: linear mean, QM: quadratic mean, CV: constant variance, LV: log-linear variance, LMQ: linear mean under quadratic population, and CVL: constant variance under log-linear population.

For both correctly and incorrectly specified conditions, the number of MVs was fixed at ten ( $m=10$ ). In the correctly specified condition, all ten MVs were LMCVs. In the misspecified condition, three out of ten MVs were misspecified—one for QMCV, LMLV, and QMLV, respectively. The first seven MVs in the misspecified condition were set to be correctly specified to allow for the evaluation of the empirical false detection rate—the proportion of cases in which the tests falsely detect these correctly specified MVs as misfitting MVs in the presence of both correctly and incorrectly specified MVs.

Upon generating the data, the proposed MV-level fit assessments were performed using the following procedures: (1) estimate the model parameters $\boldsymbol {\xi }$ by ML to get $\hat {\boldsymbol {\xi }}$ , (2) construct $\boldsymbol {\hat {\eta }}(\hat {\boldsymbol {\xi }})$ and $\boldsymbol {\eta }(\hat {\boldsymbol {\xi }})$ , and apply the ratio transformation to obtain residuals by $\boldsymbol {e}_{\boldsymbol {\varphi }}=\boldsymbol {\varphi }(\boldsymbol {\hat {\eta }}(\hat {\boldsymbol {\xi }}))-\boldsymbol {\varphi }(\boldsymbol {\eta }(\hat {\boldsymbol {\xi }}))$ , (3) estimate the ACM of $\boldsymbol {e}_{\boldsymbol {\varphi }}$ , and (4) compute $z_2(x_\ell ), z_3(x_\ell )$ for $\ell =1,\dots ,Q$ to conduct pointwise z-tests and $T_2, T_3$ to conduct the overall $\chi ^2$ -tests. Computational details follow those introduced in Section 3.1.1 for Study 1, except for the number of evaluation points used. In Step 4, we adopted $31$ evenly spaced LV values between $-3$ and 3 as evaluation points (i.e., $Q=31$ ) for $z_2(x_\ell )$ and $z_3(x_\ell )$ . Among these, 11 sub-points ranging from $-2$ to $2$ were selected to construct $T_2$ and $T_3$ for illustration.