2.1 Standard error estimates for standardized regression coefficients of Z-scores

Let z_iy be the standardized criterion variable and z_ij (j = 1, …, p) be the jth standardized predictor variable. To analyze the data of unstandardized variables (i.e., y_i , x _i1, …, x_ip ) or standardized variables (i.e., z_iy , z _i1, …, z_ip ), we define the first covariance structure model as

(4) $$\begin{align}\boldsymbol{\Sigma} \left({\mathbf{D}}_1,{\mathbf{P}}_1\right)={\mathbf{D}}_1{\mathbf{P}}_1{\mathbf{D}}_1,\end{align}$$

where ${\mathbf{D}}_1=\operatorname{diag}\left(\begin{array}{cccc}{d}_y& {d}_1& \cdots & {d}_p\end{array}\right)$ is a diagonal matrix of standard deviations, ${\mathbf{P}}_1=\left(\begin{array}{rr}1& \\ {}\boldsymbol{\unicode{x3c1}} & \mathbf{R}\end{array}\right)$ is a symmetric matrix, $\boldsymbol{\unicode{x3c1}} ={\left(\begin{array}{ccc}{\rho}_{1y}& \cdots & {\rho}_{py}\end{array}\right)}^{\prime }$ is a p × 1 vector of correlations between x _i1, …, x_ip and y_i or between z _i1, …, z_ip and z_iy , and R is the same as that in (Equation 3). When the data of y_i , x _i1, …, x_ip are analyzed, the estimates of the diagonal elements of D ₁ are equal to sample standard deviations, whereas when the data of z_iy , z _i1, …, z_ip are analyzed, the estimates of the diagonal elements of D ₁ are equal to 1. In contrast, the estimates of the elements of ρ and those of the off-diagonal elements of R are correlations, which are invariant to the scales of the input data.

In some situations, only the criterion variable and a subset of the predictor variables should be standardized. Let q be a positive integer that is smaller than p, indicating the number of predictor variables that are standardized. To analyze the data of z_iy , z _i1, …, z_iq , x _i,q+1, …, x_ip , we define the second covariance structure model as

(5) $$\begin{align}\boldsymbol{\Sigma} \left({\mathbf{D}}_2,{\mathbf{P}}_2\right)={\mathbf{D}}_2{\mathbf{P}}_2{\mathbf{D}}_2,\end{align}$$

where ${\mathbf{D}}_2=\operatorname{diag}\left(\begin{array}{rrrrrrr}{d}_y& {d}_1& \cdots & {d}_q& {1}_{q+1}& \cdots & {1}_p\end{array}\right)$ is a diagonal matrix with the first q + 1 diagonal elements being standard deviations and the last p – q elements fixed to 1, ${\mathbf{P}}_2=\left(\begin{array}{ccc}1& & \\ {}\boldsymbol{\unicode{x3c1}} & \mathbf{R}& \\ {}{\boldsymbol{\unicode{x3c3}}}_{xy}& {\boldsymbol{\Phi}}_{xz}& {\boldsymbol{\Phi}}_{xx}\end{array}\right)$ is a symmetric matrix, $\boldsymbol{\unicode{x3c1}} ={\left(\begin{array}{ccc}{\rho}_{1y}& \cdots & {\rho}_{qy}\end{array}\right)}^{\prime }$ is a q × 1 vector of correlations between z _i1, …, z_iq and z_iy , $\mathbf{R}=\left(\begin{array}{cccc}1& & & \\ {}{\rho}_{21}& 1& & \\ {}\vdots & \vdots & \ddots & \\ {}{\rho}_{q1}& {\rho}_{q2}& \cdots & 1\end{array}\right)$ is a q × q correlation matrix of z _i1, …, z_iq , ${\boldsymbol{\unicode{x3c3}}}_{xy}={\left(\begin{array}{ccc}{\sigma}_{q+1,y}& \cdots & {\sigma}_{py}\end{array}\right)}^{\prime }$ is a (p – q) × 1 vector of covariances between x _i,q+1, …, x_ip and z_iy , Φ _xz is a (p – q) × q matrix of covariances between x _i,q+1, …, x_ip and z _i1, …, z_iq , and Φ _xx is a (p – q) × (p – q) covariance matrix of x _i,q+1, …, x_ip . The estimates of the first q + 1 diagonal elements of D ₂ should be equal to 1, because the standard deviations of Z-scores always equal to 1; the estimates of the elements of ρ and those of the off-diagonal elements of R are correlations, which are invariant to the scales of the input data; and finally, the estimates of σ _xy , Φ _xz , and Φ _xx are influenced by the scales of the last p – q predictor variables.

In terms of the estimation method, generalized least squares (GLS) and maximum likelihood (ML) can be used with the normality assumption, while ML with the Satorra–Bentler correction (referred to as MLSB hereafter; Satorra & Bentler, Reference Satorra, Bentler, von Eye and Clogg1994), the sandwich estimator (known as MLR in the covariance structure modeling literature; Huber, Reference Huber1967; White, Reference White1982), and the ADF method can be used without the normality assumption. All these estimation methods provide not only the point estimates but also the asymptotic covariance matrix of the estimates. Because the parameters to be estimated in (Equations 4) and (5) can also be estimated without specifying any models, one can calculate the sample statistics first and use them as the starting values. As a result, the optimization algorithm can always reach convergence and avoid numerical problems. Essentially, the purpose of estimating the models in (Equations 4) and (5) is to obtain the asymptotic covariance matrix of the estimates. At the same time, there is no need to derive any new formulas.

Given the estimates from the first two models, one can compute the standardized regression coefficients. Let ${\boldsymbol{\unicode{x3b2}}}_1={\left(\begin{array}{ccc}{\beta}_1& \cdots & {\beta}_p\end{array}\right)}^{\prime }$ be a p × 1 vector of standardized regression coefficients of z _i1, …, z_ip . An estimate of β ₁ can be obtained from the estimates of the first model. That is, ${\widehat{\boldsymbol{\unicode{x3b2}}}}_1={\widehat{\mathbf{R}}}^{-1}\widehat{\boldsymbol{\unicode{x3c1}}}$ . Let ${\boldsymbol{\unicode{x3b2}}}_2={\left(\begin{array}{cccccc}{\beta}_1& \cdots & {\beta}_q& {b}_{q+1}& \cdots & {b}_p\end{array}\right)}^{\prime }$ be a p × 1 vector of both standardized regression coefficients of z _i1, …, z_iq and unstandardized regression coefficients of x _i,q+1, …, x_ip . An estimate of β ₂ can be obtained from the estimates of the second model. That is, ${\widehat{\boldsymbol{\unicode{x3b2}}}}_2={\left(\begin{array}{cc}\widehat{\mathbf{R}}& \\ {}{\widehat{\boldsymbol{\Phi}}}_{xz}& {\widehat{\boldsymbol{\Phi}}}_{xx}\end{array}\right)}^{-1}\left(\begin{array}{c}\widehat{\boldsymbol{\unicode{x3c1}}}\\ {}{\widehat{\boldsymbol{\unicode{x3c3}}}}_{xy}\end{array}\right)$ .

Based on the formulas of ${\widehat{\boldsymbol{\unicode{x3b2}}}}_1$ and ${\widehat{\boldsymbol{\unicode{x3b2}}}}_2$ , one can apply the delta method to compute $\mathrm{SE}\left({\widehat{\boldsymbol{\unicode{x3b2}}}}_1\right)$ and $\mathrm{SE}\left({\widehat{\boldsymbol{\unicode{x3b2}}}}_2\right)$ . Specifically, given the asymptotic covariance matrix of the estimates obtained from the first two models, one only needs the Jacobian matrices of partial derivatives of β ₁ and β ₂ with respect to the parameters of the first two models, separately. Jones and Waller (Reference Jones and Waller2015, Eqs. (25) and (26)) provided the matrix formulas of the partial derivatives of β ₁ with respect to correlations. In contrast, the partial derivatives of β ₂ involve both covariances and correlations, but no formula has been derived for such a problem. Following the recommendations by Lord (Reference Lord1975) and Browne and Du Toit (Reference Browne and Du Toit1992), we compute the partial derivatives numerically for the purpose of automated implementation of the delta method. Admittedly, computational efficiency can be further improved with the use of analytic derivatives, but any changes in the delta method will require extra efforts to derive analytic derivatives and correspondingly modify the software code.

Finally, it should be noted that our approach to compute $\mathrm{SE}\left({\widehat{\boldsymbol{\unicode{x3b2}}}}_1\right)$ is essentially the same as the correlation-based delta method as Jones and Waller (Reference Jones and Waller2015, Section 3) with the advantage that the asymptotic covariance matrix of the estimated correlations (that is a subset of the estimates) can be obtained automatically as a byproduct of parameter estimation, whereas Jones and Waller (Reference Jones and Waller2015, Section 3) referred to the formulas derived elsewhere (e.g., Nel, Reference Nel1985). Additionally, if a robust method (e.g., MLSB, MLR, ADF) is used to estimate the model in (Equation 4), we can easily get the robust estimator of the asymptotic covariance matrix of the estimated correlations. In principle, given a particular estimation method, the correlation-based delta method and the covariance-based delta method should produce the same $\mathrm{SE}\left({\widehat{\boldsymbol{\unicode{x3b2}}}}_1\right)$ .

2.2 Standard error estimates for unstandardized regression coefficients of Z-scores and products of Z-scores

In some situations, quadratic or higher-order interaction terms are involved in research questions (LeBreton et al., Reference LeBreton, Tonidandel and Krasikova2013). We show that our analytic approach can accommodate the interaction terms of any orders. With p predictor variables, one can create ${{}_2C}_p$ second-order interaction terms, ${{}_3C}_p$ third-order interaction terms, …, ${{}_{p-1}C}_p$ (p – 1)th-order interaction terms, and ${{}_pC}_p$ pth-order interaction term. In total, there are up to ${{}_2C}_p+{{}_3C}_p+\cdots +{{}_{p-1}C}_p+{{}_pC}_p$ interaction terms that can be included in the regression model. Let k be an integer between 1 and ${{}_2C}_p+{{}_3C}_p+\cdots +{{}_{p-1}C}_p+{{}_pC}_p$ and ω _i1, …, ω_ik be the interaction terms created from the Z-scores of predictor variables. To analyze the data of z_iy , z _i1, …, z_ip , ω _i1, …, ω_ik , we define the third covariance structure model as

(6) $$\begin{align}\boldsymbol{\Sigma} \left({\mathbf{D}}_3,{\mathbf{P}}_3\right)={\mathbf{D}}_3{\mathbf{P}}_3{\mathbf{D}}_3,\end{align}$$

where ${\mathbf{D}}_3=\operatorname{diag}\left(\begin{array}{rrrrrrr}{d}_y& {d}_1& \cdots & {d}_p& {1}_1& \cdots & {1}_k\end{array}\right)$ is a diagonal matrix with the first p + 1 diagonal elements being standard deviations and the last k elements fixed to 1, ${\mathbf{P}}_3=\left(\begin{array}{ccc}1& & \\ {}\boldsymbol{\unicode{x3c1}} & \mathbf{R}& \\ {}{\boldsymbol{\unicode{x3c3}}}_{\omega y}& {\boldsymbol{\Phi}}_{\omega z}& {\boldsymbol{\Phi}}_{\omega \omega}\end{array}\right)$ is a symmetric matrix, $\boldsymbol{\unicode{x3c1}} ={\left(\begin{array}{ccc}{\rho}_{1y}& \cdots & {\rho}_{py}\end{array}\right)}^{\prime }$ is a p × 1 vector of correlations between z _i1, …, z_ip and z_iy , $\mathbf{R}=\left(\begin{array}{cccc}1& & & \\ {}{\rho}_{21}& 1& & \\ {}\vdots & \vdots & \ddots & \\ {}{\rho}_{p1}& {\rho}_{p2}& \cdots & 1\end{array}\right)$ is a p × p correlation matrix of z _i1, …, z_ip , ${\boldsymbol{\unicode{x3c3}}}_{\omega y}={\left(\begin{array}{ccc}{\sigma}_{\omega_1,y}& \cdots & {\sigma}_{\omega_k,y}\end{array}\right)}^{\prime }$ is a k × 1 vector of covariances between ω _i1, …, ω_ik and z_iy , Φ _ωz is a k × p matrix of covariances between ω _i1, …, ω_ik and z _i1, …, z_ip , and Φ _ωω is a k × k covariance matrix of ω _i1, …, ω_ik .

When only q (< p) predictor variables are standardized and p – q predictor variables remain unstandardized, one can create ${{}_2C}_q+{{}_3C}_q+\cdots +{{}_{q-1}C}_q+{{}_qC}_q$ interaction terms with the Z-scores of predictor variables, ${{}_2C}_{p-q}+{{}_3C}_{p-q}+\cdots +{{}_{p-q-1}C}_{p-q}+{{}_{p-q}C}_{p-q}$ interaction terms with the unstandardized predictor variables, and $\left({{}_2C}_p+{{}_3C}_p+\cdots +{{}_{p-1}C}_p+{{}_pC}_p\right)$ - $\left({{}_2C}_q+{{}_3C}_q+\cdots +{{}_{q-1}C}_q+{{}_qC}_q\right)$ - $\left({{}_2C}_{p-q}+{{}_3C}_{p-q}+\cdots +{{}_{p-q-1}C}_{p-q}+{{}_{p-q}C}_{p-q}\right)$ interaction terms with the combinations of unstandardized and standardized predictor variables. No matter how the interaction terms are created, we still use ω _i1, …, ω_ik to denote the interaction terms. To analyze the data of z_iy , z _i1, …, z_iq , x _i,q+1, …, x_ip , ω _i1, …, ω_ik , we define the fourth covariance structure model as

(7) $$\begin{align}\boldsymbol{\Sigma} \left({\mathbf{D}}_4,{\mathbf{P}}_4\right)={\mathbf{D}}_4{\mathbf{P}}_4{\mathbf{D}}_4,\end{align}$$

where ${\mathbf{D}}_4=\operatorname{diag}\left(\begin{array}{cccccccccc}{d}_y& {d}_1& \cdots & {d}_q& {1}_{q+1}& \cdots & {1}_p& {1}_1& \cdots & {1}_k\end{array}\right)$ is a diagonal matrix with the first q + 1 diagonal elements being standard deviations and the last p – q + k elements fixed to 1, ${\mathbf{P}}_4=\left(\begin{array}{ccc}1& & \\ {}\boldsymbol{\unicode{x3c1}} & \mathbf{R}& \\ {}{\boldsymbol{\unicode{x3c3}}}_{\left(\begin{array}{cc}x& \omega \end{array}\right),y}& {\boldsymbol{\Phi}}_{\left(\begin{array}{cc}x& \omega \end{array}\right),z}& {\boldsymbol{\Phi}}_{\left(\begin{array}{cc}x& \omega \end{array}\right),\left(\begin{array}{cc}x& \omega \end{array}\right)}\end{array}\right)$ is a symmetric matrix, $\boldsymbol{\unicode{x3c1}} ={\left(\begin{array}{ccc}{\rho}_{1y}& \cdots & {\rho}_{qy}\end{array}\right)}^{\prime }$ is a q × 1 vector of correlations between z _i1, …, z_iq and z_iy , $\mathbf{R}=\left(\begin{array}{cccc}1& & & \\ {}{\rho}_{21}& 1& & \\ {}\vdots & \vdots & \ddots & \\ {}{\rho}_{q1}& {\rho}_{q2}& \cdots & 1\end{array}\right)$ is a q × q correlation matrix of z _i1, …, z_iq , ${\boldsymbol{\unicode{x3c3}}}_{\left(\begin{array}{cc}x& \omega \end{array}\right),y}={\left(\begin{array}{cccccc}{\sigma}_{q+1,y}& \cdots & {\sigma}_{py}& {\sigma}_{\omega_1,y}& \cdots & {\sigma}_{\omega_k,y}\end{array}\right)}^{\prime }$ is a (p – q + k) × 1 vector of covariances between x _i,q+1, …, x_ip , ω _i1, …, ω_ik and z_iy , Φ _{(x ω),z} is a (p – q + k) × q matrix of covariances between x _i,q+1, …, x_ip , ω _i1, …, ω_ik and z _i1, …, z_iq , and Φ _{(x ω),(x ω)} is a (p – q + k) × (p – q + k) covariance matrix of x _i,q+1, …, x_ip , ω _i1, …, ω_ik .

To estimate the third and fourth models, one should always use the robust methods (e.g., MLSB, MLR, and ADF), because the interaction terms are not normally distributed regardless of the distributions of the predictor variables. With an appropriate estimation method, the asymptotic covariance matrix of the estimates can be easily obtained from the software packages. Because the point estimates of the third and fourth models can also be obtained from the corresponding sample statistics, these point estimates can be used as the starting values to facilitate the optimization algorithm to reach convergence immediately. Essentially, the purpose of estimating the models in (Equations 6) and (7) is to obtain the asymptotic covariance matrix of the estimates. Again, there is no need to derive any new formulas.

Given the estimates from the third and fourth models, one can compute the relevant regression coefficients. Let ${\boldsymbol{\unicode{x3b2}}}_3={\left(\begin{array}{cccccc}{\beta}_1& \cdots & {\beta}_p& {b}_1& \cdots & {b}_k\end{array}\right)}^{\prime }$ be a (p + k) × 1 vector of regression coefficients of z _i1, …, z_ip , ω _i1, …, ω_ik . An estimate of β ₃ can be computed from the estimates of the third model. That is, ${\widehat{\boldsymbol{\unicode{x3b2}}}}_3={\left(\begin{array}{cc}\widehat{\mathbf{R}}& \\ {}{\widehat{\boldsymbol{\Phi}}}_{\omega, z}& {\widehat{\boldsymbol{\Phi}}}_{\omega, \omega}\end{array}\right)}^{-1}\left(\begin{array}{c}\widehat{\boldsymbol{\unicode{x3c1}}}\\ {}{\widehat{\boldsymbol{\unicode{x3c3}}}}_{\omega, y}\end{array}\right)$ . Let ${\boldsymbol{\unicode{x3b2}}}_4={\left(\begin{array}{ccccccccc}{\beta}_1& \cdots & {\beta}_q& {b}_{q+1}& \cdots & {b}_p& {b}_{p+1}& \cdots & {b}_{p+k}\end{array}\right)}^{\prime }$ be a vector of regression coefficients of z _i1, …, z_iq , x _i,q+1, …, x_ip , ω _i1, …, ω_ik . An estimate of β ₄ can be computed from the estimates of the fourth model. That is, ${\widehat{\boldsymbol{\unicode{x3b2}}}}_4={\left(\begin{array}{cc}\widehat{\mathbf{R}}& \\ {}{\widehat{\boldsymbol{\Phi}}}_{\left(\begin{array}{cc}x& \omega \end{array}\right),z}& {\widehat{\boldsymbol{\Phi}}}_{\left(\begin{array}{cc}x& \omega \end{array}\right),\left(\begin{array}{cc}x& \omega \end{array}\right)}\end{array}\right)}^{-1}\left(\begin{array}{c}\widehat{\boldsymbol{\unicode{x3c1}}}\\ {}{\widehat{\boldsymbol{\unicode{x3c3}}}}_{\left(\begin{array}{cc}x& \omega \end{array}\right),y}\end{array}\right)$ .

Based on the formulas of ${\widehat{\boldsymbol{\unicode{x3b2}}}}_3$ and ${\widehat{\boldsymbol{\unicode{x3b2}}}}_4$ , one can apply the delta method to compute $\mathrm{SE}\left({\widehat{\boldsymbol{\unicode{x3b2}}}}_3\right)$ and $\mathrm{SE}\left({\widehat{\boldsymbol{\unicode{x3b2}}}}_4\right)$ . Specifically, given the asymptotic covariance matrix of the estimates obtained from the third and fourth models, one only needs the Jacobian matrices of partial derivatives of β ₃ and β ₄ with respect to the parameters of the third and fourth models, separately. For automated implementation of the delta method, we compute the partial derivatives numerically.

3.1 Software implementations

To implement our analytic approach, the first step is to estimate one of the covariance structure models defined in this article, and the second step is to apply the delta method. In the first step, we use the lavaan (Rosseel, Reference Rosseel2012) and OpenMx (Neale et al., Reference Neale, Hunter, Pritikin, Zahery, Brick, Kirkpatrick, Estabrook, Bates, Maes and Boker2016) packages in R and SAS PROC CALIS to estimate the covariance structure models. The model specifications in the OpenMx package and SAS PROC CALIS are easy and straightforward, because the OpenMx package and SAS PROC CALIS allow the user to define and specify matrices directly. In contrast, the model specifications in the lavaan package rely on the use of latent variables, each latent variable corresponding to a manifest variable, so that the standard deviations and correlations can be estimated properly. Additionally, all the predictor variables are modeled (rather than treated as fixed variables) in the first step. Because no interaction term is involved in the first example, we can use five estimation methods: ML, GLS, MLSB, MLR, and ADF. In contrast, we can only use MLSB, MLR, and ADF due to the interaction term in the second example. Among the software packages, the lavaan package provides all five estimation methods, SAS PROC CALIS provides ML, GLS, MLSB, and ADF, and the OpenMx package only provides ML and ADF. In the second step, it is much more convenient to use the OpenMx package than the lavaan package to implement the delta method, because the OpenMx package allows users to define additional parameters with matrix formulas, whereas matrix formulas must be translated to scalar functions for the lavaan package. For the second example, the lavaan package is only used in the first step of our analytic approach due to the inconvenience of working out the scalar functions.

To implement the existing methods, we also use several software packages. Specifically, we use the lavaan package to estimate Bentler and Lee’s constrained model, the betaDelta package to implement the covariance-based delta method, and the stdmod package (Cheung, Reference Cheung2024) to construct the percentile-based CIs. Additionally, we take the bootstrap results generated from the stdmod package to calculate the bootstrap standard errors, which are the standard deviations across bootstrap replications. All the programs used in this section can be found in the Supplementary Material of this article.

Table 1

Point estimates of D ₁ and P ₁ and associated standard error estimates from the first step of our analytic approach in Example 1

Note: * indicates a fixed value.

Table 2

Standardized regression coefficients and associated standard error estimates from the second step of our analytic approach, two existing methods, and regular regression in Example 1

3.2 Example 1

The data of the first example were retrieved from the electronic supplementary material of Kwan and Chan (Reference Kwan and Chan2011). In this example, the criterion variable is child’s reading ability (y_i ), and the predictor variables are parental occupational status (x _i1), parental educational level (x _i2), and child’s home possession (x _i3). The sample size is 200.

To compute the standard error estimates for the three standardized regression coefficients, we estimate the covariance structure model defined in (Equation 4) and apply the correlation-based delta method. Table 1 shows the point estimates of D ₁ and P ₁ and the associated standard error estimates from the first step of our analytic approach, and Table 2 shows the standard error estimates for the three standardized regression coefficients from the second step of our analytic approach, two existing methods,Footnote ⁵ and regular regression. Several observations can be made. First, the standard error estimates for the estimates of D ₁ and P ₁ in Table 1 are the same between GLS and ML, and this is also the case between MLSB, MLR, and ADF. Second, given the first observation, our analytic approach produces the same standard error estimates for the three standardized regression coefficients in Table 2 between GLS and ML and between MLSB, MLR, and ADF. Third, given a particular estimation method, our analytic approach produces the same standard error estimates as those from the existing methods: (1) the covariance-based delta method, (2) Bentler and Lee’s constrained covariance structure modeling method, and (3) Kwan and Chan’s unconstrained covariance structure modeling method (using ML only; see Kwan & Chan, Reference Kwan and Chan2011, Table 3). Fourth, the standard error estimates produced by our analytic approach using MLSB, MLR, and ADF are very close to the bootstrap standard errors based on 5,000 bootstrap replications. Note that the bootstrap standard errors obtained from the stdmod package are also close to the bootstrap standard errors obtained from EQS based on 1,000 bootstrap replications (see Kwan & Chan, Reference Kwan and Chan2011, Table 3). Finally, the standard error estimates from regular regression are inflated. For parental occupational status (x _i1), the standard error estimate from regular regression is larger than the corresponding bootstrap standard error and the one from our analytic approach using MLSB, MLR, and ADF by more than 10%. For parental educational level (x _i2), the standard error estimate from regular regression is larger than the corresponding bootstrap standard error and the one from our analytic approach using MLSB, MLR, and ADF by more than 5%.

3.3 Example 2

The data of the second example were retrieved from the stdmod package in R (Cheung et al., Reference Cheung, Cheung, Lau, Hui and Vong2022). In this example, the criterion variable is sleep duration (y_i ), and the predictor variables are emotional stability (x _i1), conscientiousness (x _i2), age (x _i3), and gender (x _i4, coded as female and male). The sample size is 500.

To study the interaction effect of emotional stability and conscientiousness in standardized metric, Cheung et al. (Reference Cheung, Cheung, Lau, Hui and Vong2022) standardized all the variables to Z-scores, except the categorical variable gender (x _i4), and used the standardized emotional stability (z _i1) and the standardized conscientiousness (z _i2) to create the interaction term (ω _i1). To match the results produced by the stdmod package, we re-code gender as 0 (for female) and 1 (for male). Given these settings, we analyze the data of z_iy , z _i1, z _i2, z _i3, x _i4, ω _i1.

To compute the standard error estimates for the five unstandardized regression coefficients, we estimate the covariance structure model defined in (Equation 7) and apply the delta method. Table 3 shows the point estimates of D ₄ and P ₄ and the associated standard error estimates from the first step of our analytic approach, and Table 4 shows the standard error estimates for the five unstandardized regression coefficients and the symmetric 95% CIs from the second step of our analytic approach, the bootstrap standard errors and percentile-based CIs obtained from the stdmod package based on 5,000 bootstrap replications, and the standard error estimates from regular regression. Several observations can be made. First, MLSB, MLR, and ADF produce the same standard error estimates for the estimates of D ₄ and P ₄ in Table 3, and this is also the case for the five unstandardized regression coefficients in Table 4. Second, the delta method produces the standard error estimates that are close to the bootstrap standard errors. Third, the symmetric 95% CIs from our analytic approach are similar to the percentile-based 95% CIs from the nonparametric bootstrap procedure. Finally, the standard error estimate for standardized age (z _i3) from regular regression is larger than the corresponding bootstrap standard error by 24% and the one from our analytic approach by 28%, whereas the standard error estimates for standardized conscientiousness (z _i2), gender (x _i4), and the interaction term (ω _i1) from regular regression are smaller than the corresponding bootstrap standard errors and those from our analytic approach by 9%–12%.

Table 3

Point estimates of D ₄ and P ₄ and associated standard error estimates from the first step of our analytic approach in Example 2

Note: * indicates a fixed value.

Table 4

Unstandardized regression coefficients, standard error estimates and symmetric confidence intervals from the second step of our analytic approach, standard error estimates and percentile-based confidence intervals from the nonparametric bootstrap procedure, and standard error estimates from regular regression in Example 2

Note: Symmetric 95% CI = point estimate ± t*×SE, where t* = 1.9648 with df = 494.

4.1 First simulation study: Data generation, analysis, and evaluation criteria

In the first simulation study, we generate random data from multivariate normal distributions to mimic the second example in the previous section. First, a random dataset of size n is generated from MVN(0 ₅, Σ) for x _i1, x _i2, x _i3, x _i4, and ε_i , where Σ = HKH, $\mathbf{H}=\left(\begin{array}{@{\!}rrrrr@{\!}} {\sigma}_x& & & & \\ {}0& {\sigma}_x& & & \\ {}0& 0& \sqrt{2}& & \\ {}0& 0& 0& \sqrt{2}& \\ {}0& 0& 0& 0& {\sigma}_{\varepsilon}\end{array}\right)$ , $\mathbf{K}=\left(\begin{array}{@{\!}cc@{\!}}\mathbf{R}& \\ {}{\mathbf{0}}_4^{\prime }& 1\end{array}\right)=\left(\begin{array}{@{\!}rrrrr@{\!}}1& & & & \\ {}r& 1& & & \\ {}.3& .5& 1& & \\ {}.6& .4& .2& 1& \\ {}0& 0& 0& 0& 1\end{array}\right)$ , r is the correlation between x _i1 and x _i2, ${\sigma}_x^2$ is the variance of x _i1 and x _i2, and ${\sigma}_{\varepsilon}^2$ is the variance of ε_i . Four factors are manipulated: (1) n, (2) ${\sigma}_x^2$ , (3) r, and (4) ${\sigma}_{\varepsilon}^2$ . The levels of these four factors are as follows: n = 100, 200, 300; ${\sigma}_x^2$ = 1, 6, 12; r = .2, .4, .8; and ${\sigma}_{\varepsilon}^2$ = 10, 20, 40. Second, we create the criterion variable from ${y}_i={x}_{i1}+{x}_{i2}+0.7{x}_{i3}+0.5{x}_{i4}+0.3{x}_{i1}{x}_{i2}+{\varepsilon}_i$ , where i = 1, …, n, and gather the data of y_i , x _i1, x _i2, x _i3, and x _i4. Third, we standardize y_i , x _i1, x _i2, and x _i3 to z_iy , z _i1, z _i2, and z _i3 while leaving x _i4 unstandardized. Finally, we create a new variable ω _i1 = z _i1 z _i2 and gather the data of z_iy , z _i1, z _i2, z _i3, x _i4, and ω _i1. We repeat the above steps 1,000 times at each combination of n, ${\sigma}_x^2$ , r, and ${\sigma}_{\varepsilon}^2$ . Thus, a total of 81,000 (= 3 × 3 × 3 × 3 × 1,000) random datasets are generated.

For data analysis, we use regular regression and our analytic approach to analyze each random data of z_iy , z _i1, z _i2, z _i3, x _i4, and ω _i1, while the stdmod package in R can take the data of either y_i , x _i1, x _i2, x _i3, x _i4 or z_iy , z _i1, z _i2, z _i3, x _i4 as the input dataset to implement the nonparametric bootstrap procedure. The bootstrap replications are set to 5,000 in our simulations.

Upon completion of data analysis, we have 1,000 point estimates, from which we calculate the mean and standard deviation, for each unstandardized regression coefficient at each combination of n, ${\sigma}_x^2$ , r, and ${\sigma}_{\varepsilon}^2$ . To evaluate the point estimate of each unstandardized regression coefficient, we calculate the relative bias of mean (relative to the true value in the population, which will be discussed in the next paragraph). As for the standard deviation, it will be used as the true standard error to calculate the relative bias of the standard error estimate. In addition, we have 1,000 standard error estimates from regular regression and our analytic approach and 1,000 bootstrap standard errors for each unstandardized regression coefficient at each combination of n, ${\sigma}_x^2$ , r, and ${\sigma}_{\varepsilon}^2$ , from which we calculate the average standard error estimate from regular regression and our analytic approach and the average bootstrap standard error. Finally, we calculate the relative bias of the average standard error estimate from regular regression and our analytic approach, as well as the relative bias of average bootstrap standard error. The general formula of relative bias is

$$\begin{align*}\mathrm{relative}\kern0.17em \mathrm{bias}=\frac{\overline{\theta}-{\theta}_0}{\theta_0},\end{align*}$$

where $\overline{\theta}$ is the mean of an unstandardized regression coefficient or the average standard error estimate across replications, and θ ₀ is the corresponding true value of the unstandardized regression coefficient or the true standard error.

One challenge in our simulations is to figure out the true values of the unstandardized regression coefficients of z _i1, z _i2, z _i3, x _i4, and ω _i1. Even though the regression model is not too complicated, it is still difficult to obtain the true values, because the moments of Z-scores and products of Z-scores must be derived analytically. As an alternative, we approximate the true values by analyzing the simulated data of huge sample sizes. That is, at each combination of ${\sigma}_x^2$ , r, and ${\sigma}_{\varepsilon}^2$ , we generate a random dataset of size 1,000,000 and regress z_iy on z _i1, z _i2, z _i3, x _i4, and ω _i1. The unstandardized regression coefficients of z _i1, z _i2, z _i3, x _i4, and ω _i1 obtained at huge sample sizes are used to approximate the true values. Additionally, we record the associated R ². Table 5 shows the approximate true values of regression coefficients and R ² at each combination of ${\sigma}_x^2$ , r, and ${\sigma}_{\varepsilon}^2$ . The values of R ² range from .1530 to .8848, covering a wide range of scenarios in social and behavioral research.

Table 5

Approximate true values of the unstandardized regression coefficients and R ² at each combination of ${\sigma}_x^2$ , r and ${\sigma}_{\varepsilon}^2$ in the first simulation study

4.2 First simulation study: Results

In the first simulation study, even though the data of x _i1, x _i2, x _i3, x _i4, and ε_i are generated from multivariate normal distributions, the interaction term (i.e., ω _i1 = z _i1 z _i2, the product of Z-scores) is, in general, not normally distributed. According to the reporting practice recommended by Cain, Zhang, and Yuan (Reference Cain, Zhang and Yuan2017), we first report a variety of percentiles of the sample univariate skewness and the sample univariate kurtosis of ω _i1 at each combination of n, ${\sigma}_x^2$ , r, and ${\sigma}_{\varepsilon}^2$ (see Univariate_skew_kurt_MVN.docx in the Supplementary Material). These percentiles cover a wide range of values that are commonly encountered in practical data analysis (see Cain et al., Reference Cain, Zhang and Yuan2017, Table 1).

Table 6

Means and standard deviations of the unstandardized regression coefficients in the first simulation study

The relative bias of the mean is included in Relative_bias_of_mean_MVN.docx in the Supplementary Material, and the main results of the first simulation study are included in Tables 6–9. Specifically, Table 6 shows the means and standard deviations of the unstandardized regression coefficients; Table 7 shows the average standard error estimate from regular regression, and the associated relative bias; Table 8 shows the average standard error estimate from our analytic approach and the associated relative bias; and Table 9 shows the average bootstrap standard error and the associated relative bias.

Table 7

The average standard error estimate from regular regression and its relative bias in the first simulation study

Table 8

The average standard error estimate from our analytic approach and its relative bias in the first simulation study

Table 9

The average bootstrap standard error and its relative bias in the first simulation study

Regarding the relative bias of mean, 16 of a total of 405 (= 5 × 3 × 3 × 3 × 3) absolute values are larger than .05, of which eight absolute values are associated with the unstandardized regression coefficient of ω _i1. Regarding the relative bias of average standard error estimate, drastically different performances are observed among regular regression, our analytic approach, and the nonparametric bootstrap procedure. Based on the pattern in Table 7, we observe that the relative bias of average standard error estimate from regular regression is mainly determined by ${\sigma}_x^2$ . More specifically, 4 absolute values of relative bias of average standard error estimate are larger than .1 at ${\sigma}_x^2$ = 6 (particularly, 1 within n = 100, 2 within n = 200, and 1 within n = 300; see the last column of Table 7), and 45 absolute values of relative bias of average standard error estimate are larger than .1 at ${\sigma}_x^2$ = 12 (particularly, 15 within each level of n; see the last column of Table 7). Based on the pattern in Table 8, we observe that the relative bias of average standard error estimate from our analytic approach is mainly determined by n, and most values of relative bias are negative.Footnote ⁶ More specifically, 24 absolute values of relative bias of average standard error estimate are larger than .1 at n = 100 (particularly, 6 within ${\sigma}_x^2$ = 1, 8 within ${\sigma}_x^2$ = 6, and 10 within ${\sigma}_x^2$ = 12; see the last column of Table 8), 7 absolute values of relative bias of average standard error estimate are larger than .1 at n = 200 (particularly, 1 within ${\sigma}_x^2$ = 1, 3 within ${\sigma}_x^2$ = 10 and 20, separately; see the last column of Table 8), and 4 absolute values of relative bias of average standard error estimate are larger than .1 at n = 300 (particularly, 4 within ${\sigma}_x^2$ = 12; see the last column of Table 8). Finally, from the nonparametric bootstrap procedure, all absolute values of relative bias of average standard error estimate in Table 9 are smaller than .1. To better demonstrate the different performances of the three methods, the values of relative bias of average standard error estimate for the unstandardized regression coefficient of ω _i1 in Tables 7–9 are plotted in Figure 1.

Figure 1

Relative bias in the first simulation study.

According to Hoogland and Boomsma (Reference Hoogland and Boomsma1998), the point estimates are acceptable when the relative biases of mean are less than .05 in absolute value, whereas the standard error estimates are acceptable when the relative biases of average standard error estimates are less than .1 in absolute value. Therefore, we can reach the following conclusions. First, more than 95% of the point estimates of unstandardized regression coefficients are acceptable (because 16/405 = 3.95% are unacceptable). Second, all standard error estimates from regular regression are acceptable at ${\sigma}_x^2$ = 1, more than 95% from regular regression are acceptable at ${\sigma}_x^2$ = 6 [because 4/(5 × 3 × 3 × 3) = 2.96% are unacceptable], and between 65% and 70% from regular regression are acceptable at ${\sigma}_x^2$ = 12 [because 45/(5 × 3 × 3 × 3) = 33.33% are unacceptable]. Third, between 80% and 85% of the standard error estimates from our analytic approach are acceptable at n = 100 [because 24/(5 × 3 × 3 × 3) = 17.78% are unacceptable], about 95% from our analytic approach are acceptable at n = 200 [because 7/(5 × 3 × 3 × 3) = 5.19% are unacceptable], and more than 95% from our analytic approach are acceptable at n = 300 [because 4/(5 × 3 × 3 × 3) = 2.96% are unacceptable]. Finally, all bootstrap standard errors are acceptable in all conditions. Remember that all these conclusions are drawn from the specific regression model, the multivariate normal distributions of x _i1, x _i2, x _i3, x _i4, and ε_i , and the particular levels of n, ${\sigma}_x^2$ , r, and ${\sigma}_{\varepsilon}^2$ used in this first simulation study.

4.3 Second simulation study: Data generation, analysis, and evaluation criteria

In the second simulation study, we use the procedure developed by Qu et al. (Reference Qu, Liu and Zhang2020) to generate the data of x _i1, x _i2, x _i3, x _i4, and ε_i from multivariate non-normal distributions. This procedure is implemented by the mnonr package (Qu & Zhang, Reference Qu and Zhang2020) in R, which requires the user to specify the covariance matrix of the variables and the population values of multivariate skewness and multivariate kurtosis. The covariance matrix Σ of x _i1, x _i2, x _i3, x _i4, and ε_i used in the second simulation study is the same as the one used in the first simulation study, and the factors manipulated are also the same (i.e., n = 100, 200, 300; ${\sigma}_x^2$ = 1, 6, 12; r = .2, .4, .8; and ${\sigma}_{\varepsilon}^2$ = 10, 20, 40). The population values of multivariate skewness and multivariate kurtosis are 3 and 35, respectively. For each dataset generated from the mnonr package in R, we follow the same steps as in the first simulation study to prepare the dataset for analysis.

Table 10

Approximate true values of the unstandardized regression coefficients and R ² at each combination of ${\sigma}_x^2$ , r and ${\sigma}_{\varepsilon}^2$ in the second simulation study

To approximate the true values of the unstandardized regression coefficients of z _i1, z _i2, z _i3, x _i4, and ω _i1 at each combination of ${\sigma}_x^2$ , r, and ${\sigma}_{\varepsilon}^2$ , we generate a random dataset of size 1,000,000 and regress z_iy on z _i1, z _i2, z _i3, x _i4, and ω _i1. Table 10 shows the approximate true values of regression coefficients and R ² at each combination of ${\sigma}_x^2$ , r, and ${\sigma}_{\varepsilon}^2$ . The values of R ² range from .1570 to .9150, covering a wide range of scenarios in real applications.

4.4 Second simulation study: Results

In the second simulation study, the data of x _i1, x _i2, x _i3, x _i4, and ε_i are generated from multivariate non-normal distributions. Thus, at each combination of n, ${\sigma}_x^2$ , r, and ${\sigma}_{\varepsilon}^2$ , we first report the percentiles of the sample multivariate skewness and the sample multivariate kurtosis of x _i1, x _i2, x _i3, x _i4, and ε_i (see Multivariate_skew_kurt_ms3mk35.docx in the Supplementary Material), followed by the percentiles of the sample univariate skewness and the sample univariate kurtosis of ω _i1 (see Univariate_skew_kurt_ms3mk35.docx in the Supplementary Material). The reported percentiles cover a wide range of values that are commonly encountered in practical data analysis (see Cain et al., Reference Cain, Zhang and Yuan2017, Tables 1 and 2). Additionally, the percentiles of the sample univariate skewness and the sample univariate kurtosis in Univariate_skew_kurt_ms3mk35.docx are larger than those in Univariate_skew_kurt_MVN.docx, indicating more severe nonnormality of ω _i1 in the second simulation study.

The relative bias of mean is included in Relative_bias_of_mean_ms3mk35.docx in the Supplementary Material, and the main results of the second simulation study are included in Tables 11–14. Specifically, Table 11 shows the mean and standard deviation of the unstandardized regression coefficients; Table 12 shows the average standard error estimate from regular regression, and the associated relative bias; Table 13 shows the average standard error estimate from our analytic approach and the associated relative bias; and Table 14 shows the average bootstrap standard error and the associated relative bias.

Regarding the relative bias of mean, 18 of a total of 405 absolute values of the relative bias of mean are larger than .05, of which 5 absolute values are associated with the unstandardized regression coefficient of ω _i1. Regarding the relative bias of average standard error estimate, dramatically different performances are observed among regular regression, our analytic approach, and the nonparametric bootstrap procedure. Based on the pattern in Table 12, we observe that the relative bias of average standard error estimate from regular regression is again mainly determined by ${\sigma}_x^2$ . More specifically, there are 43 absolute values of relative bias of average standard error estimate larger than .1 at ${\sigma}_x^2$ = 12 (particularly, 12 within n = 100, 16 within n = 200, and 15 within n = 300; see the last column of Table 12). Based on the pattern in Table 13, we observe that the relative bias of average standard error estimate from our analytic approach is again mainly determined by n, and most values of relative bias are negative. More specifically, 15 absolute values of relative bias of average standard error estimate are larger than .1 at n = 100 (particularly, 4 within ${\sigma}_x^2$ = 1, 2 within ${\sigma}_x^2$ = 6, and 9 within ${\sigma}_x^2$ = 12; see the last column of Table 13), 9 absolute values of relative bias of average standard error estimate are larger than .1 at n = 200 (particularly, 9 within ${\sigma}_x^2$ = 12; see the last column of Table 13), and 7 absolute values of relative bias of average standard error estimate are larger than .1 at n = 300 (particularly, 7 within ${\sigma}_x^2$ = 12; see the last column of Table 13). Finally, 3 absolute values of the relative bias of average bootstrap standard error are larger than .1 (particularly, 1 at the combination of n = 200 and ${\sigma}_x^2$ = 12 and 2 at the combination of n = 300 and ${\sigma}_x^2$ = 12; see the last column of Table 14). To better demonstrate the different performances of the three methods, the values of relative bias of average standard error estimate for the unstandardized regression coefficient of ω _i1 in Tables 12–14 are plotted in Figure 2.

Given these observations, we have the following conclusions. First, more than 95% of the point estimates of unstandardized regression coefficients are acceptable [because 18/405 = 4.44% are unacceptable]. Second, all standard error estimates from regular regression are acceptable at ${\sigma}_x^2$ = 1 and 6, and between 65% and 70% from regular regression are acceptable at ${\sigma}_x^2$ = 12 [because 43/(5 × 3 × 3 × 3) = 31.85% are unacceptable]. Third, between 85% and 90% of the standard error estimates from our analytic approach are acceptable at n = 100 [because 15/(5 × 3 × 3 × 3) = 11.11% are unacceptable], between 90% and 95% from our analytic approach are acceptable at n = 200 [because 9/(5 × 3 × 3 × 3) = 6.67% are unacceptable], and about 95% from our analytic approach are acceptable at n = 300 [because 7/(5 × 3 × 3 × 3) = 5.19% are unacceptable]. Finally, more than 99% bootstrap standard errors are acceptable [because 3/405 = .74% are unacceptable]. Remember that all these conclusions are drawn from the specific regression model, the particular non-normal distributions of x _i1, x _i2, x _i3, x _i4, and ε_i , and the particular levels of n, ${\sigma}_x^2$ , r, and ${\sigma}_{\varepsilon}^2$ used in the second simulation study.

Table 11

Means and standard deviations of the unstandardized regression coefficients of the second simulation study

Figure 2

Relative bias in the second simulation study.

Table 12

The average standard error estimate from regular regression and its relative bias in the second simulation study

Table 13

The average standard error estimate from our analytic approach and its relative bias in the second simulation study

Table 14

The average bootstrap standard error and its relative bias in the second simulation study