2.1 The original RA-L model

Let x be a p × 1 vector for p predictor variables and y be a q × 1 vector for q criterion variables. With p predictor variables, one can construct up to p redundancy variates. Let $\boldsymbol{\xi} ={\left({\xi}_1\kern0.5em {\xi}_2\kern0.5em \cdots \kern0.5em {\xi}_p\right)}^{\prime }$ be the vector that includes all p redundancy variates. According to Van Den Wollenberg (Reference Van Den Wollenberg1977), ξ_i (i = 1, 2, …, p) must satisfy two restrictions. First, ξ_i is uncorrelated with ξ_j (i ≠ j). Second, ξ_i has unit variance (i = 1, 2, …, p). With these restrictions, Gu et al. (Reference Gu, Yung, Cheung, Joo and Nimon2023) specified the covariance structure of the original RA-L model as

(1) $$\begin{align}\boldsymbol{\Sigma} &=\boldsymbol{\Sigma} \left({\mathbf{D}}_x,{\mathbf{D}}_y,{\mathbf{L}}_{x\xi},{\mathbf{L}}_{y\xi},{\mathbf{R}}_{y y}\right)\nonumber\\ &=\left(\begin{array}{cc}{\mathbf{D}}_x& \mathbf{0}\\ {}\mathbf{0}& {\mathbf{D}}_y\end{array}\right)\left(\begin{array}{cc}{\mathbf{L}}_{x\xi}& \mathbf{0}\\ {}\mathbf{0}& {\mathbf{I}}_q\end{array}\right)\left(\begin{array}{cc}{\mathbf{I}}_p& {\mathbf{L}}_{y\xi}^{\prime}\\ {}{\mathbf{L}}_{y\xi}& {\mathbf{R}}_{y y}\end{array}\right)\left(\begin{array}{cc}{\mathbf{L}}_{x\xi}^{\prime }& \mathbf{0}\\ {}\mathbf{0}& {\mathbf{I}}_q\end{array}\right)\left(\begin{array}{cc}{\mathbf{D}}_x& \mathbf{0}\\ {}\mathbf{0}& {\mathbf{D}}_y\end{array}\right),\end{align}$$

where I _p and I _q are identity matrices of orders p and q, separately, D _x is a p × p diagonal matrix whose diagonal elements are the standard deviations of p predictor variables, D _y is a q × q diagonal matrix whose diagonal elements are the standard deviations of q criterion variables, L _xξ is a p × p square matrix that includes the redundancy loadings (i.e., the correlations between p predictor variables and p redundancy variates), L _yξ is a q × p matrix that includes the cross-loadings (i.e., the correlations between q criterion variables and p redundancy variates), and R _yy is a q × q correlation matrix whose off-diagonal elements are the correlations of q criterion variables.

To identify the original RA-L model, three types of constraints must be imposed. The first type of constraints is applicable only when the number of predictor variables exceeds that of criterion variables by two or more (i.e., p - q ≥ 2). Specifically, let d = p − q be a positive integer. When d ≥ 2, the first type of constraints requires one to arbitrarily fix d(d − 1)/2 elements in the last d columns of L _xξ . When d = 1 or p ≤ q, the first type of constraints is not applicable. The second type of constraints is

(2) $$\begin{align}\mathrm{vecdiag}\left({\mathbf{L}}_{x\xi}{\mathbf{L}}_{x\xi}^{\prime}\right)-{\mathbf{1}}_p={\mathbf{0}}_p,\end{align}$$

where vecdiag(M) denotes a column vector created with the diagonal elements of M, and 1 _p denotes a unit vector of order p, and 0 _p denotes a null vector of order p. Finally, the third type of constraints is

(3) $$\begin{align}\mathrm{vecb}\left({\mathbf{L}}_{y\xi}^{\prime }{\mathbf{L}}_{y\xi}\right)=\mathbf{0},\end{align}$$

where vecb(M) denotes a column vector created with the off-diagonal elements below the main diagonal of M, and 0 denotes a null vector of appropriate orderFootnote ⁴. The third type of constraints indicate that ${\mathbf{L}}_{y\xi}^{\prime }{\mathbf{L}}_{y\xi}$ must be a diagonal matrix, but the number of constraints required by equation (3) depends on the relative magnitude of p and q. When p ≤ q, all p columns of L _yξ include non-zero cross-loadings. In this situation, ${\mathbf{L}}_{y\xi}^{\prime }{\mathbf{L}}_{y\xi}$ has p(p − 1)/2 unique off-diagonal elements that must be 0. When p > q, only the first q columns of L _yξ include non-zero cross-loadings, while the last d = p - q columns of L _yξ are null vectors (see Appendix A of Gu et al. Reference Gu, Yung, Cheung, Joo and Nimon2023). In this situation, the first q × q submatrix of ${\mathbf{L}}_{y\xi}^{\prime }{\mathbf{L}}_{y\xi}$ has q(q − 1)/2 unique off-diagonal elements that must be 0. This completes the three types of constraints for the original RA-L model.

To count the number of parameters of the RA-L model, it is obvious that D _x has p standard deviations, D _y has q standard deviations, and R _yy has q(q − 1)/2 correlations. For L _xξ and L _yξ , however, the number of parameters in these two matrices also depends on the relative magnitude of p and q. For p ≤ q, L _xξ has p ² redundancy loadings, and L _yξ has pq cross-loadings. For p > q, L _xξ has p ² − d(d − 1)/2 = (p ² + 2pq − q ² + p − q)/2 redundancy loadings, and L _yξ has q ² cross-loadings in the first q columns because the last d columns of L _yξ are null vectors. Finally, given the number of constraints for identification and the number of parameters, we can verify that the RA-L model is a saturated model regardless of the relative magnitude of p and q (see Appendix B of Gu et al. Reference Gu, Yung, Cheung, Joo and Nimon2023).

2.2 Matrix partitions

To specify the two modified RA-L models in the next two subsections, it is necessary to partition some matrices of the original RA-L model. Let m be a positive integer that indicates the number of redundancy variates to be rotated. When p ≤ q, m must be equal to or less than p. When p > q, m must be equal to or less than q, because there is no need to rotate the last d = p − q redundancy variates.

With these settings, we first partition L _xξ as

(4) $$\begin{align}{\mathbf{L}}_{x\xi}=(\!\begin{array}{cc}{\mathbf{L}}_{x\xi \mid m}& {\mathbf{L}}_{x\xi \mid u}\end{array}\!),\end{align}$$

where ${\mathbf{L}}_{x\xi \mid m}$ is a p × m matrix, ${\mathbf{L}}_{x\xi \mid u}$ is a p × u matrix, and u = p - m. Correspondingly, the submatrices I _p and L _yξ in $\left(\!\begin{array}{cc}{\mathbf{I}}_p& {\mathbf{L}}_{y\xi}^{\prime}\\ {}{\mathbf{L}}_{y\xi}& {\mathbf{R}}_{y y}\end{array}\!\right)$ of equation (1) should be partitioned as

(5) $$\begin{align}{\mathbf{I}}_p=\left(\!\begin{array}{cc}{\mathbf{I}}_m& \mathbf{0}\\ {}\mathbf{0}& {\mathbf{I}}_u\end{array}\!\right)\kern0.36em \mathrm{and}\kern0.36em {\mathbf{L}}_{y\xi}=\left(\!\begin{array}{cc}{\mathbf{L}}_{y\xi \mid m}& {\mathbf{L}}_{y\xi \mid u}\end{array}\!\right),\end{align}$$

where ${\mathbf{L}}_{y\xi \mid m}$ is a q × m matrix and ${\mathbf{L}}_{y\xi \mid u}$ is a q × u matrix.

Based on the partitions in equations (4) and (5), the covariance structure of the original RA-L model can be re-written as

(6) $$\begin{align}\boldsymbol{\Sigma}&=\boldsymbol{\Sigma} \left({\mathbf{D}}_x,{\mathbf{D}}_y,{\mathbf{L}}_{x\xi \mid m},{\mathbf{L}}_{x\xi \mid u},{\mathbf{L}}_{y\xi \mid m},{\mathbf{L}}_{y\xi \mid u},{\mathbf{R}}_{y y}\right)\nonumber\\ &=\left(\!\begin{array}{cc}{\mathbf{D}}_x& \mathbf{0}\\ {}\mathbf{0}& {\mathbf{D}}_y\end{array}\!\right)\left(\!\begin{array}{cc}\left(\!\begin{array}{cc}{\mathbf{L}}_{x\xi \mid m}& {\mathbf{L}}_{x\xi \mid u}\end{array}\!\right)& \mathbf{0}\\ {}\mathbf{0}& {\mathbf{I}}_q\end{array}\!\right)\left(\!\begin{array}{cc}\left(\!\begin{array}{cc}{\mathbf{I}}_m& \mathbf{0}\\ {}\mathbf{0}& {\mathbf{I}}_u\end{array}\!\right)& \left(\!\begin{array}{c}{\mathbf{L}}_{y\xi \mid m}^{\prime}\\ {}{\mathbf{L}}_{y\xi \mid u}^{\prime}\end{array}\!\right)\\ {}\left(\!\begin{array}{cc}{\mathbf{L}}_{y\xi \mid m}& {\mathbf{L}}_{y\xi \mid u}\end{array}\!\right)& {\mathbf{R}}_{y y}\end{array}\!\right)\left(\!\begin{array}{cc}\left(\!\begin{array}{c}{\mathbf{L}}_{x\xi \mid m}^{\prime}\\ {}{\mathbf{L}}_{x\xi \mid u}^{\prime}\end{array}\!\right)& \mathbf{0}\\ {}\mathbf{0}& {\mathbf{I}}_q\end{array}\!\right)\left(\!\begin{array}{cc}{\mathbf{D}}_x& \mathbf{0}\\ {}\mathbf{0}& {\mathbf{D}}_y\end{array}\!\right).\end{align}$$

In the next two subsections, we will show the effect of orthogonal and oblique rotations on ${\mathbf{L}}_{x\xi \mid m}$ , I _m , and ${\mathbf{L}}_{y\xi \mid m}$ in equation (6) and define the two modified RA-L models for orthogonal and oblique rotations, separately.

2.3 The modified RA-L model for orthogonal rotations

When the first m redundancy variates are rotated with an orthogonal rotation method, ${\mathbf{L}}_{x\xi \mid m}$ is transformed by an m × m orthogonal matrix T^orth to produce ${\mathbf{L}}_{x\xi \mid m}^{\mathrm{orth}}$ , which is a p × m matrix that includes the rotated redundancy loadings. That is,

(7) $$\begin{align}{\mathbf{L}}_{x\xi \mid m}{\mathbf{T}}^{\mathrm{orth}}={\mathbf{L}}_{x\xi \mid m}^{\mathrm{orth}}.\end{align}$$

At the same time, I _m and ${\mathbf{L}}_{y\xi \mid m}$ are also transformed by T^orth. For I _m , the transformation is

(8) $$\begin{align}{\left({\mathbf{T}}^{\mathrm{orth}}\right)}^{-1}{\mathbf{I}}_m{\left({\mathbf{T}}^{\mathrm{orth}}\right)}^{\prime -1}={\left({\mathbf{T}}^{\mathrm{orth}}\right)}^{-1}\left({\mathbf{T}}^{\mathrm{orth}}\right)={\mathbf{I}}_m.\end{align}$$

For ${\mathbf{L}}_{y\xi \mid m}$ , the transformation is

(9) $$\begin{align}{\mathbf{L}}_{y\xi \mid m}{\left({\mathbf{T}}^{\mathrm{orth}}\right)}^{\prime -1}={\mathbf{L}}_{y\xi \mid m}{\mathbf{T}}^{\mathrm{orth}}={\mathbf{L}}_{y\xi \mid m}^{\mathrm{orth}}.\end{align}$$

Obviously, ${\mathbf{L}}_{y\xi \mid m}^{\mathrm{orth}}$ is a q × m matrix that includes the rotated cross-loadings. Given equations (7)–(9), the covariance structure of the modified RA-L model for orthogonal rotations is defined as

(10) $$\begin{align}\boldsymbol{\Sigma} &=\boldsymbol{\Sigma} \left({\mathbf{D}}_x,{\mathbf{D}}_y,{\mathbf{L}}_{x\xi \mid m}^{\mathrm{orth}},{\mathbf{L}}_{x\xi \mid u},{\mathbf{L}}_{y\xi \mid m}^{\mathrm{orth}},{\mathbf{L}}_{y\xi \mid u},{\mathbf{R}}_{y y}\right)\nonumber\\ &=\left(\!\begin{array}{cc}{\mathbf{D}}_x& \mathbf{0}\\ {}\mathbf{0}& {\mathbf{D}}_y\end{array}\!\right)\left(\!\begin{array}{cc}\left(\!\begin{array}{cc}{\mathbf{L}}_{x\xi \mid m}^{\mathrm{orth}}& {\mathbf{L}}_{x\xi \mid u}\end{array}\!\right)& \mathbf{0}\\ {}\mathbf{0}& {\mathbf{I}}_q\end{array}\!\right)\left(\!\begin{array}{cc}\left(\!\begin{array}{cc}{\mathbf{I}}_m& \mathbf{0}\\ {}\mathbf{0}& {\mathbf{I}}_u\end{array}\!\right)& \left(\!\begin{array}{c}{\left({\mathbf{L}}_{y\xi \mid m}^{\mathrm{orth}}\right)}^{\prime}\\ {}{\mathbf{L}}_{y\xi \mid u}^{\prime}\end{array}\!\right)\\ {}\left(\!\begin{array}{cc}{\mathbf{L}}_{y\xi \mid m}^{\mathrm{orth}}& {\mathbf{L}}_{y\xi \mid u}\end{array}\!\right)& {\mathbf{R}}_{y y}\end{array}\!\right)\left(\!\begin{array}{cc}\left(\!\begin{array}{c}{\left({\mathbf{L}}_{x\xi \mid m}^{\mathrm{orth}}\right)}^{\prime}\\ {}{\mathbf{L}}_{x\xi \mid u}^{\prime}\end{array}\!\right)& \mathbf{0}\\ {}\mathbf{0}& {\mathbf{I}}_q\end{array}\!\right)\left(\!\begin{array}{cc}{\mathbf{D}}_x& \mathbf{0}\\ {}\mathbf{0}& {\mathbf{D}}_y\end{array}\!\right).\end{align}$$

To identify the modified RA-L model for orthogonal rotations, we must impose four types of constraints. The first three types of constraints are inherited with or without changes from the three types of constraints for the original RA-L model, whereas the fourth type of constraints is introduced to remove rotational indeterminacy. The first type of constraints is identical to that for the original RA-L model. That is, when x has 2 or more variables than y, one should arbitrarily fix d(d − 1)/2 elements in the last d columns of ${\mathbf{L}}_{x\xi \mid u}$ .

The second type of constraints involves both rotated and unrotated redundancy loadings. That is,

(11) $$\begin{align}\mathrm{vecdiag}\left[\left(\!\begin{array}{cc}{\mathbf{L}}_{x\xi \mid m}^{\mathrm{orth}}& {\mathbf{L}}_{x\xi \mid u}\end{array}\!\right)\left(\!\begin{array}{c}{\left({\mathbf{L}}_{x\xi \mid m}^{\mathrm{orth}}\right)}^{\prime}\\ {}{\mathbf{L}}_{x\xi \mid u}^{\prime}\end{array}\!\right)\right]-{\mathbf{1}}_p={\mathbf{0}}_p.\end{align}$$

Compared to the p constraints in equation (2), the first m constraints in equation (11) are different, because these constraints are imposed on the rotated redundancy loadings in ${\mathbf{L}}_{x\xi \mid m}^{\mathrm{orth}}$ .

To derive the third type of constraints, we must express ${\mathbf{L}}_{y\xi}^{\prime }{\mathbf{L}}_{y\xi}$ in equation (3) with the partitioned matrix ${\mathbf{L}}_{y\xi}=\left(\!\begin{array}{cc}{\mathbf{L}}_{y\xi \mid m}& {\mathbf{L}}_{y\xi \mid u}\end{array}\!\right)$ . That is,

$$\begin{align*}{\mathbf{L}}_{y\xi}^{\prime }{\mathbf{L}}_{y\xi}&=\left(\!\begin{array}{c}{\mathbf{L}}_{y\xi \mid m}^{\prime}\\ {}{\mathbf{L}}_{y\xi \mid u}^{\prime}\end{array}\!\right)\left(\!\begin{array}{cc}{\mathbf{L}}_{y\xi \mid m}& {\mathbf{L}}_{y\xi \mid u}\end{array}\!\right)\\ &=\left(\!\begin{array}{cc}{\mathbf{L}}_{y\xi \mid m}^{\prime }{\mathbf{L}}_{y\xi \mid m}& {\mathbf{L}}_{y\xi \mid m}^{\prime }{\mathbf{L}}_{y\xi \mid u}\\ {}{\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid m}& {\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid u}\end{array}\!\right).\end{align*}$$

Given the constraints required by equation (3), we can see that ${\mathbf{L}}_{y\xi \mid m}^{\prime }{\mathbf{L}}_{y\xi \mid m}$ and ${\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid u}$ must be diagonal matrices and ${\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid m}$ must be a null matrix. Thus, we can re-write equation (3) as

$$\begin{align*}\left[\begin{array}{c}\mathrm{vec}\mathrm{b}\left({\mathbf{L}}_{y\xi \mid m}^{\prime }{\mathbf{L}}_{y\xi \mid m}\right)\\ {}\mathrm{vec}\mathrm{b}\left({\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid u}\right)\\ {}\mathrm{vec}\left({\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid m}\right)\end{array}\right]=\mathbf{0},\end{align*}$$

where vec(M) denotes a column vector created with all elements of M. With orthogonal rotations, ${\mathbf{L}}_{y\xi \mid m}$ should be substituted with ${\mathbf{L}}_{y\xi \mid m}^{\mathrm{orth}}={\mathbf{L}}_{y\xi \mid m}{\mathbf{T}}^{\mathrm{orth}}$ so that the first and last components in the above expression must be changed as follows:

$$\begin{align*}\left\{\begin{array}{c}\mathrm{vec}\mathrm{b}\left[{\left({\mathbf{L}}_{y\xi \mid m}^{\mathrm{orth}}\right)}^{\prime }{\mathbf{L}}_{y\xi \mid m}^{\mathrm{orth}}\right]\\ {}\mathrm{vec}\mathrm{b}\left({\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid u}\right)\\ {}\mathrm{vec}\left({\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid m}^{\mathrm{orth}}\right)\end{array}\right\}=\left\{\begin{array}{c}\mathrm{vec}\mathrm{b}\left[{\left({\mathbf{T}}^{\mathrm{orth}}\right)}^{\prime }{\mathbf{L}}_{y\xi \mid m}^{\prime }{\mathbf{L}}_{y\xi \mid m}{\mathbf{T}}^{\mathrm{orth}}\right]\\ {}\mathrm{vec}\mathrm{b}\left({\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid u}\right)\\ {}\mathrm{vec}\left({\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid m}{\mathbf{T}}^{\mathrm{orth}}\right)\end{array}\right\}.\end{align*}$$

It is easy to verify that $\mathrm{vecb}\left({\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid u}\right)$ and $\mathrm{vec}\left({\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid m}{\mathbf{T}}^{\mathrm{orth}}\right)$ remain to be null vectors after orthogonal rotations, but $\mathrm{vecb}\left[{\left({\mathbf{T}}^{\mathrm{orth}}\right)}^{\prime }{\mathbf{L}}_{y\xi \mid m}^{\prime }{\mathbf{L}}_{y\xi \mid m}{\mathbf{T}}^{\mathrm{orth}}\right]$ may not be a null vector, because ${\left({\mathbf{T}}^{\mathrm{orth}}\right)}^{\prime }{\mathbf{L}}_{y\xi \mid m}^{\prime }{\mathbf{L}}_{y\xi \mid m}{\mathbf{T}}^{\mathrm{orth}}$ in general is an m × m symmetric matrix. It means that rotation violates the first m(m − 1)/2 constraints required by equation (3). Therefore, the third type of constraints for the modified RA-L model for orthogonal rotations is

(12) $$\begin{align}\left[\begin{array}{c}\mathrm{vec}\mathrm{b}\left({\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid u}\right)\\ {}\mathrm{vec}\left({\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid m}^{\mathrm{orth}}\right)\end{array}\right]=\mathbf{0}.\end{align}$$

In the fourth type of constraints, the results derived by Archer and Jennrich (Reference Archer and Jennrich1973) are adapted to remove rotational indeterminacy for orthogonal rotations. That is, the fourth type of constraints requires ${\left({\mathbf{L}}_{x\xi \mid m}^{\mathrm{orth}}\right)}^{\prime}\frac{\partial {h}^{\mathrm{orth}}}{\partial {\mathbf{L}}_{x\xi \mid m}^{\mathrm{orth}}}$ to be a symmetric matrix, where ${h}^{\mathrm{orth}}={h}^{\mathrm{orth}}\left({\mathbf{L}}_{x\xi \mid m}^{\mathrm{orth}}\right)$ denotes the simplicity function of ${\mathbf{L}}_{x\xi \mid m}^{\mathrm{orth}}$ for a particular orthogonal rotation criterion, and this type of constraints includes m(m − 1)/2 constraints. Formally, we can write the fourth type of constraints as

(13) $$\begin{align}\mathrm{vecb}\left[{\left({\mathbf{L}}_{x\xi \mid m}^{\mathrm{orth}}\right)}^{\prime}\frac{\partial {h}^{\mathrm{orth}}}{\partial {\mathbf{L}}_{x\xi \mid m}^{\mathrm{orth}}}-\frac{\partial {h}^{\mathrm{orth}}}{\partial {\left({\mathbf{L}}_{x\xi \mid m}^{\mathrm{orth}}\right)}^{\prime }}{\mathbf{L}}_{x\xi \mid m}^{\mathrm{orth}}\right]={\mathbf{0}}_{m\left(m-1\right)/2}.\end{align}$$

This completes the four types of constraints for the modified RA-L model for orthogonal rotations.

It can be seen that the number of parameters of the modified RA-L model for orthogonal rotations is the same as that of the original RA-L model, because orthogonal rotations do not increase the number of parameters. As for the number of constraints, equation (12) has m(m − 1)/2 fewer constraints than equation (3), while equation (13) introduces m(m − 1)/2 new constraints. Therefore, the modified RA-L model for orthogonal rotations is still a saturated model.

2.4 The modified RA-L model for oblique rotations

When the first m redundancy variates are rotated with an oblique rotation method, ${\mathbf{L}}_{x\xi \mid m}$ is transformed by an m × m nonsingular matrix T^obli that must satisfy the restriction $\operatorname{diag}{\left[{\left({\mathbf{T}}^{\mathrm{obli}}\right)}^{\prime }{\mathbf{T}}^{\mathrm{obli}}\right]}^{-1}={\mathbf{I}}_m$ to produce ${\mathbf{L}}_{x\xi \mid m}^{\mathrm{obli}}$ , which is a p × m matrix that includes the rotated redundancy loadings. That is,

(14) $$\begin{align}{\mathbf{L}}_{x\xi \mid m}{\mathbf{T}}^{\mathrm{obli}}={\mathbf{L}}_{x\xi \mid m}^{\mathrm{obli}}.\end{align}$$

At the same time, I _m and ${\mathbf{L}}_{y\xi \mid m}$ are also transformed by T^obli. For I _m , the transformation is

(15) $$\begin{align}{\left({\mathbf{T}}^{\mathrm{obli}}\right)}^{-1}{\mathbf{I}}_m{\left({\mathbf{T}}^{\mathrm{obli}}\right)}^{\prime -1}={\left[{\left({\mathbf{T}}^{\mathrm{obli}}\right)}^{\prime }{\mathbf{T}}^{\mathrm{obli}}\right]}^{-1}=\boldsymbol{\Phi},\end{align}$$

where is Φ a m × m correlation matrixFootnote ⁵ of the rotated redundancy variates. For ${\mathbf{L}}_{y\xi \mid m}$ , the transformation is

(16) $$\begin{align}{\mathbf{L}}_{y\xi \mid m}{\left({\mathbf{T}}^{\mathrm{obli}}\right)}^{\prime -1}={\mathbf{L}}_{y\xi \mid m}^{\mathrm{obli}},\end{align}$$

where ${\mathbf{L}}_{y\xi \mid m}^{\mathrm{obli}}$ is a q × m matrix that includes the rotated cross-loadings. Based on equations (14)–(16), the covariance structure of the modified RA-L model for oblique rotations is defined as

(17) $$\begin{align}\boldsymbol{\Sigma}&=\boldsymbol{\Sigma} \left({\mathbf{D}}_x,{\mathbf{D}}_y,{\mathbf{L}}_{x\xi \mid m}^{\mathrm{obli}},{\mathbf{L}}_{x\xi \mid u},\boldsymbol{\Phi}, {\mathbf{L}}_{y\xi \mid m}^{\mathrm{obli}},{\mathbf{L}}_{y\xi \mid u},{\mathbf{R}}_{y y}\right)\nonumber\\ &=\left(\!\begin{array}{cc}{\mathbf{D}}_x& \mathbf{0}\\ {}\mathbf{0}& {\mathbf{D}}_y\end{array}\!\right)\left(\!\begin{array}{cc}\left(\!\begin{array}{cc}{\mathbf{L}}_{x\xi \mid m}^{\mathrm{obli}}& {\mathbf{L}}_{x\xi \mid u}\end{array}\!\right)& \mathbf{0}\\ {}\mathbf{0}& {\mathbf{I}}_q\end{array}\!\right)\left(\!\begin{array}{cc}\left(\!\begin{array}{cc}\boldsymbol{\Phi} & \mathbf{0}\\ {}\mathbf{0}& {\mathbf{I}}_u\end{array}\!\right)& \left(\!\begin{array}{c}{\left({\mathbf{L}}_{y\xi \mid m}^{\mathrm{obli}}\right)}^{\prime}\\ {}{\mathbf{L}}_{y\xi \mid u}^{\prime}\end{array}\!\right)\\ {}\left(\!\begin{array}{cc}{\mathbf{L}}_{y\xi \mid m}^{\mathrm{obli}}& {\mathbf{L}}_{y\xi \mid u}\end{array}\!\right)& {\mathbf{R}}_{y y}\end{array}\!\right)\left(\!\begin{array}{cc}\left(\!\begin{array}{c}{\left({\mathbf{L}}_{x\xi \mid m}^{\mathrm{obli}}\right)}^{\prime}\\ {}{\mathbf{L}}_{x\xi \mid u}^{\prime}\end{array}\!\right)& \mathbf{0}\\ {}\mathbf{0}& {\mathbf{I}}_q\end{array}\!\right)\left(\!\begin{array}{cc}{\mathbf{D}}_x& \mathbf{0}\\ {}\mathbf{0}& {\mathbf{D}}_y\end{array}\!\right).\end{align}$$

Note that equation (17) has m(m − 1)/2 more parameters than equations (6) due to the off-diagonal elements of Φ.

To identify the modified RA-L model for oblique rotations, we also need to impose four types of constraints. The first type of constraints is that when x has 2 or more variables than y, one should arbitrarily fix d(d − 1)/2 elements in the last d columns of ${\mathbf{L}}_{x\xi \mid u}$ in equation (17).

The second type of constraints involves not only the rotated and unrotated redundancy loadings but also the correlations of the rotated redundancy variates. That is,

(18) $$\begin{align}\mathrm{vecdiag}\left[\left(\!\begin{array}{cc}{\mathbf{L}}_{x\xi \mid m}^{\mathrm{obli}}& {\mathbf{L}}_{x\xi \mid u}\end{array}\!\right)\left(\!\begin{array}{cc}\boldsymbol{\Phi} & \mathbf{0}\\ {}\mathbf{0}& {\mathbf{I}}_u\end{array}\!\right)\left(\!\begin{array}{c}{\left({\mathbf{L}}_{x\xi \mid m}^{\mathrm{obli}}\right)}^{\prime}\\ {}{\mathbf{L}}_{x\xi \mid u}^{\prime}\end{array}\!\right)\right]-{\mathbf{1}}_p={\mathbf{0}}_p.\end{align}$$

Compared to the p constraints in equation (2), the first m constraints in equation (18) are different, because these m constraints involve the rotated redundancy loadings in ${\mathbf{L}}_{x\xi \mid m}^{\mathrm{obli}}$ and the correlations in Φ.

The derivation of the third type of constraints for the modified RA-L model for oblique rotations is similar to that for the orthogonal rotations. Recall that equation (3) requires

$$\begin{align*}\left[\begin{array}{c}\mathrm{vec}\mathrm{b}\left({\mathbf{L}}_{y\xi \mid m}^{\prime }{\mathbf{L}}_{y\xi \mid m}\right)\\ {}\mathrm{vec}\mathrm{b}\left({\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid u}\right)\\ {}\mathrm{vec}\left({\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid m}\right)\end{array}\right]=\mathbf{0}.\end{align*}$$

With oblique rotations, ${\mathbf{L}}_{y\xi \mid m}$ should be substituted with ${\mathbf{L}}_{y\xi \mid m}^{\mathrm{obli}}={\mathbf{L}}_{y\xi \mid m}{\left({\mathbf{T}}^{\mathrm{obli}}\right)}^{\prime -1}$ so that the first and last components in the above expression must be changed as follows:

$$\begin{align*}\left\{\begin{array}{c}\mathrm{vec}\mathrm{b}\left[{\left({\mathbf{L}}_{y\xi \mid m}^{\mathrm{obli}}\right)}^{\prime }{\mathbf{L}}_{y\xi \mid m}^{\mathrm{obli}}\right]\\ {}\mathrm{vec}\mathrm{b}\left({\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid u}\right)\\ {}\mathrm{vec}\left({\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid m}^{\mathrm{obli}}\right)\end{array}\right\}=\left\{\begin{array}{c}\mathrm{vec}\mathrm{b}\left[{\left({\mathbf{T}}^{\mathrm{obli}}\right)}^{\prime }{\mathbf{L}}_{y\xi \mid m}^{\prime }{\mathbf{L}}_{y\xi \mid m}{\mathbf{T}}^{\mathrm{obli}}\right]\\ {}\mathrm{vec}\mathrm{b}\left({\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid u}\right)\\ {}\mathrm{vec}\left({\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid m}{\mathbf{T}}^{\mathrm{obli}}\right)\end{array}\right\}.\end{align*}$$

It is easy to verify that $\mathrm{vecb}\left({\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid u}\right)$ and $\mathrm{vec}\left({\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid m}{\mathbf{T}}^{\mathrm{obli}}\right)$ remain to be null vectors after oblique rotations, but $\mathrm{vecb}\left[{\left({\mathbf{T}}^{\mathrm{obli}}\right)}^{\prime }{\mathbf{L}}_{y\xi \mid m}^{\prime }{\mathbf{L}}_{y\xi \mid m}{\mathbf{T}}^{\mathrm{obli}}\right]$ may not be a null vector, because ${\left({\mathbf{T}}^{\mathrm{obli}}\right)}^{\prime }{\mathbf{L}}_{y\xi \mid m}^{\prime }{\mathbf{L}}_{y\xi \mid m}{\mathbf{T}}^{\mathrm{obli}}$ in general is an m × m symmetric matrix. Therefore, the third type of constraints for the modified RA-L model for oblique rotations is

(19) $$\begin{align}\left[\begin{array}{c}\mathrm{vec}\mathrm{b}\left({\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid u}\right)\\ {}\mathrm{vec}\left({\mathbf{L}}_{y\xi \mid u}^{\prime }{\mathbf{L}}_{y\xi \mid m}^{\mathrm{obli}}\right)\end{array}\right]=\mathbf{0}.\end{align}$$

In the fourth type of constraints, the results derived by Jennrich (Reference Jennrich1973) are adapted to remove rotational indeterminacy for oblique rotations. That is, the fourth type of constraints requires ${\left({\mathbf{L}}_{x\xi \mid m}^{\mathrm{obli}}\right)}^{\prime}\frac{\partial {h}^{\mathrm{obli}}}{\partial {\mathbf{L}}_{x\xi \mid m}^{\mathrm{obli}}}{\boldsymbol{\Phi}}^{-1}$ to be a diagonal matrix, where ${h}^{\mathrm{obli}}={h}^{\mathrm{obli}}\left({\mathbf{L}}_{x\xi \mid m}^{\mathrm{obli}}\right)$ denotes the simplicity function of ${\mathbf{L}}_{x\xi \mid m}^{\mathrm{obli}}$ for a particular oblique rotation criterion, and this type of constraints includes m(m − 1) constraints. Formally, we can write the fourth type of constraints as

(20) $$\begin{align}\mathrm{veco}\left[{\left({\mathbf{L}}_{x\xi \mid m}^{\mathrm{obli}}\right)}^{\prime}\frac{\partial {h}^{\mathrm{obli}}}{\partial {\mathbf{L}}_{x\xi \mid m}^{\mathrm{obli}}}{\boldsymbol{\Phi}}^{-1}\right]={\mathbf{0}}_{m\left(m-1\right)},\end{align}$$

where veco(M) denotes a column vector created with all off-diagonal elements of M. This completes the four types of constraints for the modified RA-L model for oblique rotations.

It can be seen that the modified RA-L model for oblique rotations has m(m − 1)/2 more parameters (i.e., the off-diagonal elements of Φ) than the original RA-L model, equation (19) has m(m − 1)/2 fewer constraints than equation (3), and equation (20) introduces m(m − 1) new constraints. Therefore, the modified RA-L model for oblique rotations is still a saturated model.

3.1 Notations of the parameter vectors

Strictly speaking, we should use θ^orth and θ^obli to denote the parameter vectors for the two modified RA-L models, separately. With these notations, we have $\boldsymbol{\Sigma} \left({\boldsymbol{\unicode{x3b8}}}^{\mathrm{orth}}\right)=\boldsymbol{\Sigma} \left({\mathbf{D}}_x,{\mathbf{D}}_y,{\mathbf{L}}_{x\xi \mid m}^{\mathrm{orth}},{\mathbf{L}}_{x\xi \mid u},{\mathbf{L}}_{y\xi \mid m}^{\mathrm{orth}},{\mathbf{L}}_{y\xi \mid u},{\mathbf{R}}_{y y}\right)$ and $\boldsymbol{\Sigma} \left({\boldsymbol{\unicode{x3b8}}}^{\mathrm{obli}}\right)=\boldsymbol{\Sigma} \left({\mathbf{D}}_x,{\mathbf{D}}_y,{\mathbf{L}}_{x\xi \mid m}^{\mathrm{obli}},{\mathbf{L}}_{x\xi \mid u},\boldsymbol{\Phi}, {\mathbf{L}}_{y\xi \mid m}^{\mathrm{obli}},{\mathbf{L}}_{y\xi \mid u},{\mathbf{R}}_{y y}\right)$ . However, to avoid repetitive descriptions in this section, we use θ as a generic symbol to denote the parameter vector for both modified RA-L models. As such, $\boldsymbol{\Sigma} \left(\boldsymbol{\theta} \right)$ is used to refer to either $\boldsymbol{\Sigma} \left({\boldsymbol{\theta}}^{\mathrm{orth}}\right)$ or $\boldsymbol{\Sigma} \left({\boldsymbol{\theta}}^{\mathrm{obli}}\right)$ .

3.2 Jacobian matrix and partial differentials

For both modified RA-L models, the ULS fitting function is defined as

(21) $$\begin{align}F=0.5\mathrm{tr}{\left[\mathbf{S}-\boldsymbol{\Sigma} \left(\boldsymbol{\theta} \right)\right]}^2.\end{align}$$

Then, the estimating equations have the following form

(22) $$\begin{align}\mathbf{g}\left(\boldsymbol{\theta}, \mathbf{S}\right)=\left[\begin{array}{c}\frac{\partial F}{\partial \boldsymbol{\theta}}\\ {}{\boldsymbol{\varphi}}_1\left(\boldsymbol{\theta} \right)\\ {}{\boldsymbol{\varphi}}_2\left(\boldsymbol{\theta} \right)\\ {}{\boldsymbol{\varphi}}_3\left(\boldsymbol{\theta} \right)\end{array}\right]=\mathbf{0},\end{align}$$

where ${\boldsymbol{\varphi}}_1\left(\boldsymbol{\theta} \right)$ , ${\boldsymbol{\varphi}}_2\left(\boldsymbol{\theta} \right)$ , and ${\boldsymbol{\varphi}}_3\left(\boldsymbol{\theta} \right)$ represent the second, third, and fourth type of constraints for either modified RA-L model. Specifically, ${\boldsymbol{\varphi}}_1\left(\boldsymbol{\theta} \right)$ includes p constraints from either equation (11) for orthogonal rotations or equation (18) for oblique rotations, ${\boldsymbol{\varphi}}_2\left(\boldsymbol{\theta} \right)$ includes p(p − 1)/2 − m(m − 1)/2 or q(q − 1)/2 − m(m − 1)/2 constraints, depending on the relative magnitude of p and q, from either equation (12) for orthogonal rotations or equation (19) for oblique rotations, and ${\boldsymbol{\varphi}}_3\left(\boldsymbol{\theta} \right)$ includes either m(m − 1)/2 constraints from equation (13) for orthogonal rotations or m(m − 1) constraints from equation (20) for oblique rotations.

Given equation (22), the Jacobian matrix of $\mathbf{g}\left(\boldsymbol{\theta}, \mathbf{S}\right)$ with respect to θ is

(23) $$\begin{align}\mathbf{J}\left(\boldsymbol{\theta}, \mathbf{S}\right)=\frac{\partial \mathbf{g}\left(\boldsymbol{\theta}, \mathbf{S}\right)}{\partial {\boldsymbol{\theta}}^{\prime }}=\left[\begin{array}{c}\frac{\partial^2F}{\partial \boldsymbol{\theta} \partial {\boldsymbol{\theta}}^{\prime }}\\[1pt] {}\frac{\partial {\boldsymbol{\varphi}}_1\left(\boldsymbol{\theta} \right)}{\partial {\boldsymbol{\theta}}^{\prime }}\\[1pt] {}\frac{\partial {\boldsymbol{\varphi}}_2\left(\boldsymbol{\theta} \right)}{\partial {\boldsymbol{\theta}}^{\prime }}\\[1pt] {}\frac{\partial {\boldsymbol{\varphi}}_3\left(\boldsymbol{\theta} \right)}{\partial {\boldsymbol{\theta}}^{\prime }}\end{array}\right],\end{align}$$

where $\frac{\partial^2F}{\partial \boldsymbol{\theta} \partial {\boldsymbol{\theta}}^{\prime }}$ is the Hessian matrix of the ULS fitting function, and the remaining components are the partial derivatives of the constraints with respect to θ.

Let ${\partial}_2{\mathbf{g}}_{\left(\boldsymbol{\theta}, \mathbf{S}\right)}\left(d\mathbf{S}\right)$ be the partial differential of $\mathbf{g}\left(\boldsymbol{\theta}, \mathbf{S}\right)$ with respect to S evaluated at $\left(\boldsymbol{\theta}, \mathbf{S}\right)$ , and we define k _n as

(24) $$\begin{align}{\mathbf{k}}_n&={\partial}_2{\mathbf{g}}_{\left(\boldsymbol{\theta}, \mathbf{S}\right)}\left[\left({\mathbf{z}}_n-\overline{\mathbf{z}}\right){\left({\mathbf{z}}_n-\overline{\mathbf{z}}\right)}^{\prime}\right]\nonumber\\ &=\left(\!\begin{array}{c}-{\frac{\partial \left\{\mathrm{vec}\left[\boldsymbol{\Sigma} \left(\boldsymbol{\theta} \right)\right]\right\}}{\partial \boldsymbol{\theta}}}^{\prime}\mathrm{vec}\left[\left({\mathbf{z}}_n-\overline{\mathbf{z}}\right){\left({\mathbf{z}}_n-\overline{\mathbf{z}}\right)}^{\prime}\right]\\ {}\mathbf{0}\\ {}\mathbf{0}\\ {}\mathbf{0}\end{array}\!\right),\end{align}$$

where n = 1, 2, …, N, N is the sample size, z _n is a column vector for the nth observation of all predictor and criterion variables, and $\overline{\mathbf{z}}$ is a column vector of the sample means of all predictor and criterion variables. The last three components in equation (24) are null vectors, because ${\boldsymbol{\varphi}}_1\left(\boldsymbol{\theta} \right)$ , ${\boldsymbol{\varphi}}_2\left(\boldsymbol{\theta} \right)$ , and ${\boldsymbol{\varphi}}_3\left(\boldsymbol{\theta} \right)$ are not functions of S.

3.3 Pseudo values and asymptotic covariance matrix of parameter estimates

Given the Jacobian matrix and the partial differentials, the pseudo values for each observation can be computed. Let λ _n (n = 1, …, N) be a column vector collecting the pseudo values for the nth observation, and it can be solved from

(25) $$\begin{align}\mathbf{J}\left(\boldsymbol{\theta}, \mathbf{S}\right){\boldsymbol{\lambda}}_n=-{\mathbf{k}}_n.\end{align}$$

Note that $\mathbf{J}\left(\boldsymbol{\theta}, \mathbf{S}\right)$ defined in equation (23) has more rows than columns so that the system of equations in equation (25) appears to be over-determined. Thus, we apply the QR decomposition to $\mathbf{J}\left(\boldsymbol{\theta}, \mathbf{S}\right)$ to solve for λ _n .

After λ _n is obtained for all observations, the IJ estimate of the asymptotic covariance matrix of $\hat{\boldsymbol{\theta}}$ is

(26) $$\begin{align}{\mathrm{acov}}^{\mathrm{IJ}}\left(\hat{\boldsymbol{\theta}}\right)=\mathrm{scov}\left({\boldsymbol{\lambda}}_n\right),\end{align}$$

where $\mathrm{scov}\left({\boldsymbol{\lambda}}_n\right)$ is the sample covariance matrix of all λ _n . Finally, the standard error estimates for $\hat{\boldsymbol{\theta}}$ are obtained from dividing the square roots of the diagonal elements of ${\mathrm{acov}}^{\mathrm{IJ}}\left(\hat{\boldsymbol{\theta}}\right)$ by $\sqrt{N}$ .

4.1 Data generation

Two factors are manipulated in this simulation study. The first factor is the data distribution, including 1) multivariate normality and 2) multivariate nonnormality. The second factor is the sample size, including 1) 200, 2) 400, and 3) 600. In total, there are 6 combinations of data distribution and sample size. At each combination, we use the following population covariance matrix to generate 1000 random data sets:

$$\begin{align*}{\boldsymbol{\Sigma}}_0=\left(\!\begin{array}{cc}{\boldsymbol{\Sigma}}_{xx}& {\boldsymbol{\Sigma}}_{yx}^{\prime}\\ {}{\boldsymbol{\Sigma}}_{yx}& {\boldsymbol{\Sigma}}_{yy}\end{array}\!\right),\end{align*}$$

where the first eight variables are the predictor variables and the last eight variables are the criterion variablesFootnote ⁶. The submatrices of Σ ₀ are

$$\begin{align*}&{\boldsymbol{\Sigma}}_{xx}=\left(\begin{array}{rrrrrrrr}1.00& & & & & & & \\ {}.71& 1.00& & & & & & \\ {}.72& .72& 1.00& & & & & \\ {}.73& .73& .73& 1.00& & & & \\ {}.74& .74& .74& .74& 1.00& & & \\ {}.20& .10& .10& .10& .20& 1.00& & \\ {}.10& .20& .20& .20& .10& .52& 1.00& \\ {}.20& .10& .10& .10& .20& .53& .53& 1.00\end{array}\right),\\&{\boldsymbol{\Sigma}}_{yx}=\left(\begin{array}{cccccccc}.40& .50& .35& .50& .40& .05& .04& .03\\ {}.35& .35& .40& .40& .35& .04& .02& .01\\ {}.50& .40& .50& .35& .50& .03& .01& .04\\ {}.35& .35& .40& .40& .35& .02& .03& .02\\ {}.40& .50& .35& .50& .40& .01& .05& .05\\ {}.01& .01& .02& .02& .03& .40& .30& .35\\ {}.02& .03& .01& .03& .02& .35& .40& .30\\ {}.03& .02& .03& .01& .01& .30& .35& .40\end{array}\right),\\&{\boldsymbol{\Sigma}}_{yy}=\left(\begin{array}{rrrrrrrr}1.00& & & & & & & \\ {}.51& 1.00& & & & & & \\ {}.52& .52& 1.00& & & & & \\ {}.53& .53& .53& 1.00& & & & \\ {}.54& .54& .54& .54& 1.00& & & \\ {}.20& .00& .20& .00& .20& 1.00& & \\ {}.00& .20& .00& .20& .00& .52& 1.00& \\ {}.20& .00& .20& .00& .20& .53& .53& 1.00\end{array}\right).\end{align*}$$

To generate the multivariate normal data, the RANDNORMAL function in SAS PROC IML is used. To generate the multivariate non-normal data, we use the procedure developed by Qu et al. (Reference Qu, Liu and Zhang2020). This procedure is implemented by the MNONR package in R, which requires the user to specify the population values of multivariate skewness and multivariate kurtosis. In this simulation study, we set the values of multivariate skewness and multivariate kurtosis to 10 and 400, respectivelyFootnote ⁷.

4.2 Data analysis and evaluation criteria

By applying RA to Σ ₀, we obtain the population values of the unrotated redundancy loadings and unrotated cross-loadings:

$$\begin{align*}{\boldsymbol{L}}_{x\xi}=\left(\!\begin{array}{rrrrrrrr}.8401& .1341& -.2650& -.1836& .0465& .0428& .2225& .3447\\ {}.8901& .1721& .0814& -.2285& .0159& -.0920& .1257& -.3077\\ {}.8316& .1220& -.3780& .2809& -.0830& .1485& .0850& -.1886\\ {}.9066& .1797& .1418& .2203& -.1025& -.0981& -.1493& .1872\\ {}.8402& .1346& -.2597& -.1343& .2253& .0643& -.3674& -.0258\\ {}.0103& .8084& -.1790& .1924& .5022& .0500& .1456& .0374\\ {}.0121& .8119& -.1419& -.0416& -.2345& -.5073& -.0606& -.0527\\ {}.0123& .8081& -.1880& -.1739& -.2621& .4449& -.1104& .0502\end{array}\!\right),\end{align*}$$

$$\begin{align*}{\mathbf{L}}_{y\xi}=\left(\!\begin{array}{rrrrrrrr}.5159& .0848& .2053& -.0324& .0356& -.0033& .0099& -.0000\\ {}.4209& .0139& -.0475& .1029& -.0095& .0014& .0053& .0001\\ {}.4729& -.0492& -.3565& -.0723& .0294& .0048& -.0001& .0000\\ {}.4203& .0119& -.0564& .0895& -.0406& -.0023& -.0053& -.0001\\ {}.5159& .0779& .1938& -.0580& -.0221& .0004& -.0098& .0000\\ {}-.0568& .4327& -.0045& .0378& .0801& .0614& -.0040& -.0000\\ {}-.0536& .4222& -.0495& .0030& .0102& -.1050& -.0006& .0000\\ {}-.0577& .4204& -.0616& -.0390& -.0908& .0434& .0045& -.0000\end{array}\!\right),\end{align*}$$

and the first two population redundancy indices are .1399 and .0698, while the subsequent population redundancy indices are less than .03. Thus, for each random data set, we only rotate the first two columns of redundancy loadings. In terms of the rotation method, we use a widely accepted oblique rotation method: QUARTIMIN (Browne, Reference Browne2001; Carroll, Reference Carroll1953) with Kaiser’s normalization (Reference Kaiser1958). In general, oblique rotations are more flexible than orthogonal rotations in the sense that oblique rotations can accommodate correlations among rotated factors/variates. If the rotated factors/variates are indeed uncorrelated, the resulting correlations from oblique rotations would be small and negligible. By applying QUARTIMIN to the first two columns of unrotated redundancy loadings, we obtain the population values of rotated redundancy loadings, rotated cross-loadings, and correlation of rotated redundancy variates:

$$\begin{align*}{\mathbf{L}}_{x\xi \mid m}^{\mathrm{obli}}=\left(\!\begin{array}{rr}.8525& -.0097\\ {}.9028& .0199\\ {}.8440& -.0203\\ {}.9194& .0240\\ {}.8525& -.0092\\ {}-.0011& .8087\\ {}.0006& .8119\\ {}.0008& 8080\end{array}\!\right),\kern0.36em \boldsymbol{\Phi} =\left(\!\begin{array}{cc}1.0000& \\ {}.1826& 1.0000\end{array}\!\right),\kern0.36em {\mathbf{L}}_{y\xi \mid m}^{\mathrm{obli}}=\left(\!\begin{array}{rr}.5228& .0921\\ {}.4172& .0199\\ {}4578& -.0425\\ {}.4162& .0178\\ {}.5216& .0852\\ {}.0170& .4318\\ {}.0184& .4214\\ {}.0140& .4195\end{array}\!\right).\end{align*}$$

The normalized QUARTIMIN rotation is implemented by SAS PROC FACTOR, and the IJ method is implemented by customized code written in SAS PROC IML.

After the analyses are completed, we compute the means, standard deviations, and average standard error estimates across 1000 replications at each combination of data distribution and sample size. The standard deviations are used as the true standard errors to evaluate the performance of the IJ method. The first evaluation criterion we use is the relative bias of the average standard error estimate, which is calculated as

$$\begin{align*}\mathrm{Relative}\kern0.17em \mathrm{bias}=\frac{\mathrm{Avg}\;\mathrm{SE}-\mathrm{SD}}{\mathrm{SD}}.\end{align*}$$

According to Hoogland and Boomsma (Reference Hoogland and Boomsma1998), the standard error estimate is acceptable when the absolute value of relative bias is less than .1. Additionally, we use the estimate and the associated standard error estimate to construct a symmetric 95% confidence interval (CI) and evaluate if the population value is included in the symmetric 95% CI. Thus, the second evaluation criterion is the coverage rate for each parameter across 1000 replications.

Table 1 Results from simulations under multivariate normality

Note: Parm = parameter, SD = standard deviation, Avg SE = average standard error, lx denotes the element of ${\mathbf{L}}_{x\xi \mid m}^{\mathrm{obli}}$ , ϕ denotes the element of Φ, ly denotes the element of ${\mathbf{L}}_{y\xi \mid m}^{\mathrm{obli}}$ , and the subscript after lx, ϕ, and ly refers to the location of the element in the corresponding matrix.

Table 2 Results from simulations under multivariate nonnormality

Note: Parm = parameter, SD = standard deviation, SE = standard error, lx denotes the element of ${\mathbf{L}}_{x\xi \mid m}^{\mathrm{obli}}$ , ϕ denotes the element of Φ, ly denotes the element of ${\mathbf{L}}_{y\xi \mid m}^{\mathrm{obli}}$ , and the subscript after lx, ϕ, and ly refers to the location of the element in the corresponding matrix.

4.3 Results

Because our purpose is to validate the standard error estimates from the IJ method, the means of rotated estimates are omitted in this section but can be found from the Supplementary Materials. Instead, we show the standard deviations, average standard errors, relative biases, and coverage rates in Tables 1 and 2 under multivariate normality and multivariate nonnormality, separately. It is observed that 1) the means are getting closer to their population values as the sample size increases, 2) all the absolute values of relative biases are less than 0.1, and 3) all the coverage rates are close to 95%. Therefore, we conclude that the IJ method performs well under both multivariate normality and multivariate nonnormality.

5.1 Example 1

In the first example, we use the data from van Dam and van Trijp (Reference Van Dam and Van Trijp2011), who collected 851 survey responses from the light users of sustainable products and applied RA to predict 10 variables measuring the motivational structure of sustainability by 15 variables that include psychographic variables and purchase behavior. The 10 motivational structure variables are healthiness (y ₁), price (y ₂), convenience (y ₃), naturalness (y ₄), taste (y ₅), local production (y ₆), environment friendliness (y ₇), fair trade (y ₈), animal friendliness (y ₉), and waste (y ₁₀). The 15 predictor variables are concern for future consequences (x ₁), prevention focus (x ₂), promotion focus (x ₃), altruistic value (x ₄), biospheric value (x ₅), egoistic value (x ₆), NEP Footnote ⁸ scale (x ₇), connectedness to nature (x ₈), environment affect (x ₉), ethical orientation (x ₁₀), health prevention (x ₁₁), health promotion (x ₁₂), social SVO Footnote ⁹ (x ₁₃), individual SVO (x ₁₄), and competitive SVO (x ₁₅). More details of these variables can be found from van Dam and van Trijp (Reference Van Dam and Van Trijp2011).

By applying RA, we find that the first three redundancy indices are .2503, .0357, and .0074, which are exactly the same as those reported by van Dam and van Trijp (Reference Van Dam and Van Trijp2011, p. 736), and all subsequent redundancy indices are smaller than .005. According to van Dam and van Trijp (Reference Van Dam and Van Trijp2011), the first two redundancy indices are meaningful, and the third and subsequent redundancy indices can be ignored. Thus, we focus on the first two columns of the redundancy loadings and the cross-loadings.

To obtain the standard error estimates for unrotated RA estimates, we fit the original RA-L model. The estimation method we use include maximum likelihood (ML), which requires the multivariate normality assumption of the data, and ML with the Satorra–Bentler correction (referred to as MLSB hereafter), which does not require any distribution assumptions of the data. Table 3 shows the first two columns of L _xξ and L _yξ and the associated standard error estimates from ML and MLSB, separately.

Table 3 The first two columns of unrotated redundancy loadings and unrotated cross-loadings and the associated standard error estimates from ML and MLSB

Note: SE = standard error estimate, ML = maximum likelihood, MLSB = maximum likelihood with the Satorra–Bentler correction.

By applying the normalized VARIMAX, we obtain ${\mathbf{L}}_{x\xi \mid m}^{\mathrm{orth}}$ and ${\mathbf{L}}_{y\xi \mid m}^{\mathrm{orth}}$ . To obtain the standard error estimates for rotated RA estimates, we fit the modified RA-L model for orthogonal rotations estimated by ULS, and apply the IJ method described in this paper. Table 4 shows ${\mathbf{L}}_{x\xi \mid m}^{\mathrm{orth}}$ , ${\mathbf{L}}_{y\xi \mid m}^{\mathrm{orth}}$ , and the associated standard error estimates from the IJ method.

Table 4 Rotated redundancy loadings, rotated cross-loadings, and the associated standard error estimates from the IJ method

Note: SE = standard error estimate, IJ = infinitesimal jackknife. The rotated redundancy loadings whose absolute values are significantly larger than .3 are in boldface.

Using the standard error estimates, we can test if the absolute value of a rotated redundancy loading in ${\mathbf{L}}_{x\xi \mid m}^{\mathrm{orth}}$ is larger than some cutoff value. Because the rotated redundancy loadings are correlations, we take .3 as the cutoff value, which means that at least 9% of the variance of a predictor variable must be shared with a rotated redundancy variate. Because we need to test the statistical significance of 30 rotated redundancy loadings simultaneously, it is necessary to adjust the typical significance level of .05. For convenience, we use the Bonferroni adjustment so that the adjusted significance level is .00167. It means that we will select a rotated redundancy loading if the associated p-value is smaller than .00167.

Based on the selected rotated redundancy loadings, we use the corresponding predictor variables to interpret the rotated redundancy variates. Specifically, the first rotated redundancy variate should be interpreted in terms of biospheric value (x ₅), NEP scale (x ₇), connectedness to nature (x ₈), environment affect (x ₉), and ethical orientation (x ₁₀); however, all the rotated redundancy loadings are smaller than .3 in the second column of ${\mathbf{L}}_{x\xi \mid m}^{\mathrm{orth}}$ . Accordingly, the first rotated redundancy variate can be interpreted as people’s concern for environmental sustainability.

It is worth noting that if we only compared the absolute values of rotated redundancy variates against .3 but did not consider the sampling variability, we would select two more rotated redundancy loadings in the first column of ${\mathbf{L}}_{x\xi \mid m}^{\mathrm{orth}}$ (i.e., .5698 and .5495) that correspond with altruistic value (x ₄) and health prevention (x ₁₁). Nevertheless, the significance tests indicate that the rotated redundancy loadings on these two variables are not really larger than .3, and their magnitude observed in this example just appears to be larger than .3 due to randomness. If these two variables would be used to interpret the first rotated redundancy variate, it would totally change the current interpretation of the first rotated redundancy variate. This reflects the advantage of the use of standard error estimates in selecting the rotated redundancy loadings.

Table 5 Results of the individual redundancy indices and cumulative redundancy for the real example

Note: CI = confidence interval, ML = maximum likelihood, MLSB = maximum likelihood with the Satorra–Bentler correction.

5.2 Example 2

In the second example, we use the data from Jurukasemthawee et al. (Reference Jurukasemthawee, Pisitsungkagarn, Taephant and Sittiwong2021) that collected responses from 424 young adults (mean age = 19.97, standard deviation of age = 1.64) on 9 psychological variables, serving as the predictor variables, and 7 spiritual well-being variables, serving as the outcome variables. The 9 predictor variables are family and environment background (x ₁), crisis in life that contributed to self-development (x ₂), positive personal predisposition (x ₃), good role models (x ₄), faith activities (x ₅), mindfulness and self-regulation (x ₆), voluntary activities (x ₇), self-reflection (x ₈), and listening to positive experience (x ₉). The 7 spiritual well-being variables are: inner peace (y ₁), acceptance in diversity (y ₂), compassion (y ₃), self-transcendence (y ₄), value in self (y ₅), meaning in life (y ₆), and insight in learnings (y ₇). Each of the predictor and outcome variables is computed from the sum of item scores that are measured on a Likert scale ranging from 0 to 6. The number of items used for each of the predictor and outcome variables is from 5 to 12 items. More details of these items can be found from Jurukasemthawee et al. (Reference Jurukasemthawee, Pisitsungkagarn, Taephant and Sittiwong2021).

To determine the dimensionality in this example, we apply a new criterion proposed by Gu et al. (Reference Gu, Yung, Cheung, Joo and Nimon2023), which relies on the inferential information of redundancy indices. Specifically, we need to compare the lower limit of the 95% confidence interval (CI) for cumulative redundancy with some cutoff value. As a result, the smallest cumulative redundancy, of which the lower limit is larger than the specified cutoff value, can be identified. The identified cumulative redundancy determines the dimensionality in RA. In other words, we should retain the individual redundancy indices that constitute the identified cumulative redundancy. As for the cutoff value, we choose .3, meaning that at least 30% of the variance of criterion variables must be explained. To apply this new criterion, we need to fit the original RA-L model. As for the estimation method, we still use ML and MLSB.

Table 5 shows the results of the individual redundancy indices and cumulative redundancy for this example. By examining the lower limit of the 95% CI of cumulative redundancy, we find that the second cumulative redundancy is the smallest cumulative redundancy whose lower limit is larger than .3. It means that we should retain the first two individual redundancy indices. In addition, we notice that the second and third redundancy indices have comparable magnitude and both of them are distinctively larger than the fourth and subsequent redundancy indices, all of which are smaller than .01. Thus, we further study the difference between the second and third redundancy indices and their sumFootnote ¹⁰. The results in Table 6 show that the 95% CI for the difference includes 0, indicating that the second and third redundancy indices are not significantly different; simultaneously, the lower limit of the 95% CI for their sum is larger than .06 and the upper limit is nearly .10, indicating that the second and third redundancy indices can explain about 6–10% of the variance of criterion variables. Based on these results, we decide to retain the first three redundancy variates. The unrotated redundancy loadings and unrotated cross-loadings of the first three redundancy variates are shown in Table 7.

Table 6 Difference between the 2nd and 3rd individual redundancy indices and sum of the 2nd and 3rd individual redundancy indices

Note: Difference = the 2nd individual redundancy index − the 3rd individual redundancy index, Sum = the 2nd individual redundancy index + the 3rd individual redundancy index, CI = confidence interval, ML = maximum likelihood, MLSB = maximum likelihood with the Satorra–Bentler correction.

Table 7 The unrotated redundancy loadings and unrotated cross-loadings for the first three redundancy variates

Note: SE = standard error estimate, ML = maximum likelihood, MLSB = maximum likelihood with the Satorra–Bentler correction.

By applying the normalized QUARTIMIN, we obtain ${\mathbf{L}}_{x\xi \mid m}^{\mathrm{obli}}$ , ${\mathbf{L}}_{y\xi \mid m}^{\mathrm{obli}}$ , and Φ. To obtain the standard error estimates for rotated RA estimates, we fit the modified RA-L model for oblique rotations estimated by ULS, and apply the IJ method described in this paper. Table 8 shows ${\mathbf{L}}_{x\xi \mid m}^{\mathrm{obli}}$ , ${\mathbf{L}}_{y\xi \mid m}^{\mathrm{obli}}$ , and Φ, and the associated standard error estimates from the IJ method.

Table 8 Results of the rotated redundancy loadings, the rotated cross-loadings, and the correlations of the three rotated redundancy variates

Note: SE = standard error estimate, IJ = infinitesimal jackknife. The rotated redundancy loadings whose absolute values are significantly larger than .3 are in boldface.

Using the standard error estimates, we can test if the absolute value of a rotated redundancy loading in ${\mathbf{L}}_{x\xi \mid m}^{\mathrm{obli}}$ is larger than some cutoff value. Again, we take .3 as the cutoff value. Because we need to test the statistical significance of 27 rotated redundancy loadings simultaneously, it is necessary to adjust the typical significance level of .05. We use the Bonferroni adjustment again so that the adjusted significance level is .00185. It means that we will select a rotated redundancy loading if the associated p-value is smaller than .00185.

Based on the selected rotated redundancy loadings, we use the corresponding predictor variables to interpret the three rotated redundancy variates. Specifically, the first rotated redundancy variate should be interpreted in terms of positive personal predisposition (x ₃), voluntary activities (x ₇), self-reflection (x ₈), and listening to positive experience (x ₉); the second rotated redundancy variate should be interpreted in terms of family and environment background (x ₁), crisis in life that contributed to self-development (x ₂), and Mindfulness and Self-Regulation (x ₆); and the third rotated redundancy variate should be interpreted in terms of faith activities (x ₅). Accordingly, the first rotated redundancy variate can be interpreted as positive personal predispositions that facilitated attention to positive experiences, self-reflection, and voluntary activities; the second rotated redundancy variate can be interpreted as safe family and environmental backgrounds that facilitated the use of mindfulness and self-regulation in transforming crisis into self-development; and the third rotated redundancy variate can be interpreted as engagement in activities that were related to own faiths. Also, we found that the correlation between the first and second rotated redundancy variates is .7029 (with standard error estimate = .0265), suggesting that the first and second rotated redundancy variates share almost 50% of their variance. It implies that positive personal predispositions and safe family and environmental backgrounds are closely and significantly related. It should be noted that only oblique rotations can produce correlated rotated redundancy variates and the resulting correlations may bring more meaningful interpretations and insights to the study than the orthogonal rotations.

It is worth noting that if we only compared the absolute values of rotated redundancy variates against .3 but did not consider the sampling variability, we would select one more rotated redundancy loading in the third column of ${\mathbf{L}}_{x\xi \mid m}^{\mathrm{obli}}$ (i.e., .3102) that corresponds with voluntary activities (x ₇). Nevertheless, the significance test indicates that the rotated redundancy loading on this variable is not really larger than .3. If this variable would be used to interpret both the first and third rotated redundancy variates, it would cause some inconvenience in the interpretation, which in turn reflects the advantage of the use of standard error estimates in selecting the rotated redundancy loadings.