Do the Mega and Titan Tests Yield Accurate Results? An Investigation into Two Experimental Intelligence Tests

: The Mega and Titan Tests were designed by Ronald K. Hoeﬂin to make ﬁne distinctions in the intellectual stratosphere. The Mega Test purported to measure above-average adult IQ up to and including scores with a rarity of one in a million of the general population. The Titan Test was billed as being even more di ﬃ cult than the Mega Test. In this article, these claims are subjected to scrutiny. Both tests are renormed using the normal curve of distribution. It is found that the Mega Test has a higher ceiling and a lower ﬂoor than the Titan Test. While the Mega Test may thus seem preferable as a psychometric instrument, it is somewhat marred by a number of easy items in its verbal section. Although o ﬃ cial scores reported to test-takers are too high, it is likely that the Mega Test does stretch to the one in a million level. The Titan Test does not. Testees who had previously taken standard intelligence tests achieved average scores of 135–145 IQ on those. Since the mean of all scores on the Mega and Titan Tests was found to be IQ 137 and IQ 138, respectively, testees had considerable scope to ﬁnd their true level without ceiling e ﬀ ects. Both are unusual and non-standard tests which require a great deal of e ﬀ ort to complete. Nevertheless, they deserve consideration as they represent an inventive experimental method of measuring the very highest levels of human intelligence and have been taken by enough subjects to allow norming.


Introduction
Intelligence tests were invented by Alfred Binet and his student Théodore Simon in 1905 with the purpose of identifying pupils in need of remedial help in French public education.Within a few years, they had been translated into English and were to reach their apogee in the United States where Lewis M. Terman, a young professor of education at Stanford University, made his reputation as the foremost authority on all matters connected with intelligence.Terman's first book on the topic, The Measurement of Intelligence, featured examples of individuals within the various classifications [1].By the time his The Intelligence of School Children was published three years later, it was clear that Terman's primary interest was in subjects scoring at the highest levels [2].He had already begun a study of exceptional children, which became the basis for longitudinal research into the lives and careers of the gifted.
This study, published in five volumes as Genetic Studies of Genius between 1926 and 1959, required the construction of a special instrument to accommodate Terman's subjects as adults, called the Concept Mastery Test [3].This marked the beginning of experimental research on adults in the intellectual stratosphere using psychological techniques.Due to the rarity of the individuals concerned, it was fraught with practical difficulties.One possible method was to give adolescents achievement tests designed for adults.That was the approach chosen by the Study of Mathematically Precocious Youth, which began in 1971 at Johns Hopkins University and which, despite its name, also considered verbal

Introduction
Structural equation models (SEMs) and confirmatory factor analysis (CFA) are important statistical methods for analyzing multivariate data in the social sciences [1][2][3][4][5].In these models, a multivariate vector X = (X 1 , . . ., X I ) of I continuous observed variables (also referred to as items or indicators) is modeled as a function of a vector of latent variables (i.e., factors or traits) η.SEMs represent the mean vector µ and the covariance matrix Σ of the random variable X as a function of an unknown parameter vector θ.In SEMs, constrained estimation of the moment structure of the multivariate normal distribution is applied [6].
The measurement model in an SEM is given as We denote the covariance matrix Var( ) = Ψ.The vectors η and are multivariate normally distributed.In addition, η and are uncorrelated random vectors.In CFA, the multivariate normal (MVN) distribution is represented as η ∼ MVN(α, Φ) and ∼ MVN(0, Ψ).Hence, one can represent the mean and the covariance matrix in CFA as µ(θ) = ν + Λα and Σ(θ) = ΛΦΛ + Ψ . ( In SEM, a matrix B of regression coefficients can additionally be specified such that η = Bη + ξ with E(ξ) = α and Var(ξ) = Φ .
where I is the identity matrix.
Researchers often parsimoniously parameterize the mean vector and the covariance matrix using a parameter θ as a summary in an SEM.The model assumptions in SEMs are, at best, merely an approximation of a true data-generating model.In SEMs, model deviations (i.e., model errors) in covariances emerge as a difference between a population covariance matrix Σ and a model-implied covariance matrix Σ(θ) (see [7][8][9]).Simultaneously, model errors in the mean vector cause a difference between the population mean vector µ and the model-implied mean vector µ(θ).As a result, the SEM is misspecified at the population level.It should be noted that the model errors are defined at the population level in infinite sample sizes.In real-data applications with limited sample sizes, the empirical covariance matrix S estimates the population covariance matrix Σ, while the mean vector x estimates the population mean vector µ.
In this work, estimators with some resistance to model deviations are investigated.In more depth, the presence of some amount of model errors should have no impact on the parameter estimate θ.This robustness property is denoted as model robustness, and it adheres to robust statistics principles [10][11][12].Model errors in SEMs appear as residuals in the modeled mean vector and the modeled covariance matrix, whereas in traditional robust statistics, observations (i.e., cases or subjects) do not obey an imposed statistical model and should be considered as outliers.That is, an estimator in an SEM should automatically recognize large deviations in µ − µ(θ) and Σ − Σ(θ) as outliers that should not significantly damage the estimated parameter θ.
In previous research, the L p loss function has been used for model-robust estimation based on moments [13,14].Non-robust estimators such as maximum likelihood or unweighted and weighted least square estimates will typically result in biased estimates in the presence of model error [14].In this article, we more thoroughly discuss the choice of the tuning parameter ε in differentiable approximations of the nondifferentiable L p loss functions.Furthermore, we compare the L p loss function with a recently proposed differentiable approximation of the L 0 loss function and a direct minimization of a smoothed version of the Bayesian information criterion [15] in regularized estimation.Notable, the L 0 loss function minimizes the number of model deviations in a fitted model.If only a few entries in the modeled mean vector (or mean vectors) or covariance matrix (or covariance matrices) deviate from zero at the population level, while all other entries equal zero, the L 0 loss function would be the most appropriate fit function.In contrast, if all model deviations differ from zero and unsystematically fluctuate around zero, the L p or L 0 loss functions with p ≤ 1 would be less appropriate.Finally, standard errors for the proposed modelrobust estimators based on the delta method are derived in this article.Their performance is assessed by evaluating coverage rates.
To sum up, this article focuses on implementation details of SEM estimation based on the L p (0 < p ≤ 1) and the newly proposed L 0 loss functions, while [16] was devoted to regularized SEM estimation, which can also be utilized for model-robust estimation.A comparison of regularized estimation and robust loss functions can be found in [14].
The remainder of the article is organized as follows.Model-robust SEM estimation based on the robust L 0 and L p loss functions is treated in Section 2. Section 3 introduces direct BIC minimization as a special approach to regularized maximum likelihood estimation.Section 4 is devoted to details of standard error computation.In Section 5, research questions are formulated that are addressed in two subsequent simulation studies.In Section 6, bias and root mean square error of model-robust SEM estimators are of interest.Section 7 reports findings regarding the standard error estimation regarding coverage rates.Finally, the article closes with a discussion in Section 8.

L 0 and L p Loss Functions in SEM Estimation
We now describe model-robust moment estimation of multiple-group SEMs.The treatment closely follows previous work in Refs.[14,16].
The empirical mean vector x and the empirical covariance matrix S are sufficient statistics for estimating µ and Σ when modeling multivariate normally distributed data with no missing values.In particular, they are also sufficient statistics of µ(θ) and Σ(θ) that are constrained functions of a parameter vector θ.Hence, x and S are also sufficient statistics for the parameter vector θ = (θ 1 , . . ., θ K ) that contains K elements.Now, assume that there are G groups with sample sizes N g , mean vectors x g , and covariance matrices S g (g = 1, . . ., G).Let ξ g = x g , vech(S g ) be the vector of sufficient statistics in group g, where vech denotes the operator that stacks all nonredundant matrix entries on top of one another.Furthermore, the vector ξ = (ξ 1 , . . ., ξ G ) contains the sufficient statistics of all G groups.
The population mean vectors and covariance matrices are denoted by µ g and Σ g , respectively.The model-implied mean vectors and covariance matrices are denoted by µ g (θ) and Σ g (θ), respectively.It is worth noting that the parameter vector θ lacks an index g, indicating that there can be common and unique parameters across groups.Equal factor loadings and item intercepts across groups are frequently imposed in a multiple-group CFA (i.e., measurement invariance is specified [17,18]).
In model-robust SEM estimation discussed in this article, discrepancies x g − µ g (θ) and vech(S g ) − vech(Σ g (θ)) are minimized according to a loss function ρ.There are two kinds of errors that are only, for simplicity, discussed for the mean structure.We can express the discrepancy in the mean structure as The first term x g − µ g describes a discrepancy due to sampling variation (i.e., with respect to the sampling of subjects).This term can typically be reduced when larger samples are drawn.The second term µ g − µ g (θ) indicates a model error.This term exists at the population level and, therefore, does not vanish with increasing sample sizes.In modelrobust estimation, a few entries in the model error are allowed to differ from zero (i.e., a sparsity assumption), corresponding to model misspecification.Note that the sparsity assumption is vital for the performance of model-robust estimators.
In robust moment estimation, the following fit function F rob is minimized: In the first term at the right side of the equation in (6), discrepancies in the sample means x g,i from the model-implied mean µ g,i (θ) for item i in group g are considered.In the second term at the right side of the equation in (6), discrepancies in the sample covariances s g,ij from the model-implied covariance σ g,ij (θ) for items i and j in group g are considered.The weights w 1g,i (i = 1, . . ., I) and w 2g,ij (i, j = 1, . . ., I) are known but can be set to one if all variables have (approximately) the same standard deviation in the sample comprising all groups or the original scaling of the variables reflects the intended weighing of sampling and model errors.The loss function ρ in (6) should be chosen such that it is resistant to outlying effects in the mean and the covariance structure.
The robust mean absolute deviation (MAD) loss function ρ(x) = |x| was examined in [19,20].When compared to usually employed SEM estimation approaches, this fit function is more robust to a few model violations, such as unmodeled item intercepts or unmodeled residual correlations of residuals (see [7]).
In this article, we investigate the L p loss function It has been shown that p < 1 provides more efficient model-robust estimates than p = 1 (see [7,13]).The L p loss function with p = 2 is the square loss function ρ(x) = x 2 and corresponds to unweighted least squares (ULSs) estimation.However, this loss function does not possess the model robustness property [7].
The L p loss function ρ(x) = |x| = |x| 0.5 (i.e., p = 0.5) is implemented in invariance alignment [21][22][23] and penalized structural equation modeling [24] in the popular Mplus software (Version 8.10., https://www.statmodel.com/support/index.shtml(accessed on 25 September 2023)).The critical aspect of the L p loss function ρ defined in (7) is that it is a nondifferentiable function.Consequently, the fit function F rob in (6) is also nondifferentiable, which does not allow the application of general-purpose optimizers that rely on the differentiable optimization functions.As a remedy, the nondifferentiable L p loss function ρ can be replaced by a differentiable approximation ρ ε , which is close to ρ, but differentiable on the entire real line.The approximating function ρ ε is defined as where ε > 0 is a tuning parameter that should be small enough such that ρ ε is close to ρ but large enough to ensure estimation stability.Replacing the nondifferentiable ρ by the differentiable approximation ρ ε has been previously recommended by [13,21,25,26].
The loss function ρ and its differentiable approximation ρ ε are displayed for six different values of p in Figure 1.It can be seen (and shown) that ρ(x) ≤ ρ ε (x) for all p > 0 and ε > 0. Furthermore, the loss function ρ is much steeper at x = 0 for a smaller p.With a larger ε, the approximation ρ ε becomes smoother.Choosing an appropriate tuning parameter ε > 0 is therefore important when applying model-robust moment estimation based on the L p loss function.It might be tempting to use a very small p close to zero.Such a loss function is close to the L 0 loss function, which takes the value of 1 for all arguments differently from x = 0 and 0 for x = 0.If the sparsity assumption of model errors holds, the L 0 would theoretically be the most desirable loss function [27,28].However, as shown in Figure 1, ρ ε for p = 0.01 does not have a clear minimum zero, making this differentiable loss function difficult to apply in practical optimization.
O'Neill and Burke [15,29] proposed the following differentiable approximation χ ε of the L 0 loss function in a recent work related to regularized estimation: where ε > 0 is again a tuning parameter.The differentiable approximation χ ε is displayed for different ε values in Figure 2. It can be seen that the functional form of χ ε seems much nicer to use in optimization than ρ ε with a p close to 0. Hence, χ ε might be a useful alternative robust loss function whose performance in the presence of model errors has to be evaluated.
In practical minimization of F rob in (6), when ρ is replaced with ρ ε (using an appropriate p) or χ ε , it is advisable to use reasonable starting values and to minimize F rob using a sequence of differentiable approximations ρ ε with decreasing ε values (i.e., subsequently fitting ρε with ε = 10 −1 , 10 −2 , 10 −3 , 10 −4 , while using the previously obtained parameter estimate as the initial value for the subsequent minimization problem).
Loss function χ ε (see ( 9)) as a function of the tuning parameter ε.

A Direct BIC Minimization in Regularized Maximum Likelihood Estimation
Most frequently, SEMs are estimated with maximum likelihood (ML) estimation.This estimation method provides the most efficient estimates for correctly specified models.However, the efficiency properties are lost in the case of misspecified SEMs.
As an alternative, regularized ML estimation can be used that introduces an overidentified SEM by allowing free group-specific item intercepts and residual covariances.To identify the model, a penalty function on the overidentified parameters is imposed that is targeted at the sparsity structure of model errors.The methodological literature has extensively documented the regularized estimation of single-group and multiple-group SEMs [30][31][32][33].Cross-loadings, residual covariances, or item intercepts are regularized in these applications.Regularized SEM estimation enables flexible yet parsimonious model specifications.
In ML estimation, the fit function F ML is the negative log-likelihood function based on the multivariate normal distribution that is defined as [2,4] In empirical applications, the model-implied mean vectors µ g and covariance matrices Σ g will often be misspecified [34][35][36], and θ can be understood as a pseudo-true parameter that is defined as the minimizer of F ML in (10).
In regularized SEM estimation, a penalty function P is added to the log-likelihood fit function F ML that imposes some sparsity assumption on a subset of model parameters [31,33].In order to enforce sparsity, the penalty function P is often chosen to be nondifferentiable.Define a known parameter ι k ∈ {0, 1} for all parameters θ k , where ι k = 1 indicates that for the kth entry θ k in θ a penalty function is applied.The penalized log-likelihood function is defined as where λ > 0 is a regularization parameter, and N * is a scaling factor that frequently equals the total sample size N = ∑ G g=1 N g .The regularized (or penalized) ML estimate is defined as the minimizer of F pen (θ; ξ).
The least absolute shrinkage and selection operator (LASSO; ref. [37]) or the smoothly clipped absolute deviation (SCAD; ref. [38]) penalty functions have been frequently used in regularized SEM estimation.For a fixed value of λ, a subset of θ k parameters for which the penalty function is applied (i.e., ι k = 1) will result in estimates of zero.That is, a sparse θ vector is obtained as the result of regularized estimation.However, the estimate of θ depends on the fixed regularization parameter λ; that is, As a result, the parameter estimate θ(λ) of θ depends on the unknown parameter λ.To avoid this problem, the regularized SEM can be repeatedly estimated on a finite grid of regularization parameters λ (e.g., on an equidistant grid between 0.01 and 1.00 with increments of 0.01).The Bayesian information criterion (BIC), defined by BIC = 2F ML (θ, ξ) + log(N)H may be used to choose an optimal regularization parameter λ, where H denotes the number of parameters.Because the minimization of BIC is equivalent to the minimization of BIC/2, the final parameter estimate θ is determined as where the function χ is an indicator whether |x| is different from 0: In particular, the quantity 13) counts the number of parameter estimates θ k (λ) for k = 1, . . ., K for which the penalty function is applied (i.e., ι k = 1) and differs from 0.
As becomes clear, regularization SEM estimation necessitates fitting an SEM on a grid of the regularization parameter λ.This approach is computationally intensive, especially for SEMs with a large number of parameters.The final parameter estimate is obtained by minimizing the BIC across all estimated regularized SEMs.A naïve idea might be directly minimizing the BIC to avoid repeated estimation that involves regularization parameter selection.It should be noted that only a subset of parameters for which sparsity should be imposed is relevant in the BIC computation.Hence, a parameter estimate θ by minimizing the BIC is provided by The optimization function in (15) employs an L 0 penalty function [39,40] with a fixed regularization parameter log(N)/2.This optimization function contains the nondifferen-tiable indicator function χ that counts the number of regularized parameters that differ from 0. The ingenious idea of O'Neill and Burke [15] was to replace the nondifferentiable L 0 loss function χ with its differentiable approximation χ ε (see (8) and Ref. [16] for a more comprehensive treatment).Therefore, the parameter θ can be estimated as The estimation approach from ( 16) is referred to as the smoothed direct BIC minimization (DBIC) approach.This method has been used to estimate regularized distributional regression models [15].
It has been shown in SEMs that the DBIC approach performs similarly to regularized estimation based on the indirect approach by minimizing the BIC on a finite grid of regularization parameters [16].Hence, we confine ourselves in this article to compare the model-robust moment estimation methods with different values of the power p with the DBIC estimator.

Computation of Standard Errors
In this section, the computation of the variance matrix of parameter estimates θ from model-robust moment estimation and DBIC estimation using the fit functions in ( 6) and ( 16), respectively, is described (see also [14,16] for a similar treatment).Both methods minimize a differentiable (approximation) function F(θ, ξ) with respect to θ as a function of sufficient statistics ξ (see also [41]).The vector of estimated sufficient statistics ξ is approximately normally distributed (see [3]); that is, for a true population parameter ξ 0 of sufficient statistics.Let F θ = (∂F)/(∂θ) be the vector of partial derivatives with respect to θ.The parameter estimate θ fulfills the nonlinear equation F θ ( θ, ξ) = 0.The delta method [34] can be employed to derive the variance matrix of θ.Assume that there exists a (pseudo-)true parameter θ 0 such that F θ (θ 0 , ξ 0 ) = 0. Now, we conduct a Taylor expansion of F θ (see [3,5,42]).Denote by F θθ and F θξ the matrices of second-order partial derivatives of F θ with respect to θ and ξ, respectively.The Taylor expansion can be written as By solving (18) for θ, we get the approximation By defining A = −F θθ ( θ, ξ) −1 F θξ ( θ, ξ) when substituting θ 0 and ξ 0 with θ and ξ, respectively, we get by using the multivariate delta method [34] The square root of diagonal elements of Var( θ) computed from (20) may be used to calculate standard errors for elements in θ.

Statistical Inference for Parameter Differences of Different Models Based on the Same Dataset
In this section, statistical inference for differences in parameters from different models based on the same dataset is discussed.For example, researchers could use the L p loss function with p = 2, p = 0.5, and p = 0.It should be evaluated whether the estimated factor means from different models that employ the different loss functions are statistically significant.Importantly, the different models rely on the same dataset and its vector of sufficient statistics ξ.Hence, the standard error of a parameter difference from parameters of different models can be smaller than the standard error of a single model because the data are used twice.The M-estimation framework can also be utilized to derive the variance estimate of a parameter difference [43,44].The different loss functions provide different estimates θ m for models m = 1, . . ., M. At the population level, the parameters are denoted as θ m .Note that the population parameters differ in the case of misspecified SEMs.In the following, we discuss the case M = 2 to reduce notation.
Following the lines of the variance derivation in the previous section, we can approximate the estimate of model m by using ( 19) where A m is defined as in (19).Note that A 1 = A 2 due to choosing different loss functions.
Researchers can now ask whether the parameter difference ∆ = θ 1 − θ 2 (or some entries of it) significantly differs from 0. From ( 21), we obtain where ∆ = θ 1 − θ 2 .We then obtain the variance estimate of ∆ as The unknown matrix A 1 − A 2 can be estimated by its sample analog A 1 − A 2 .

Research Purpose
In this article, several research questions are addressed connected to model-robust estimation.First, simulations should clarify which tuning parameter ε > 0 in the differentiable approximation should be chosen to minimize bias and variance in estimated structural SEM parameters.Second, according to our knowledge, no research compares the performance of the L p loss function (approximation) with the newly proposed L 0 loss function approximation of O'Neill and Burke in terms of bias and root mean square error (RMSE).Third, it should be examined whether the standard error computation based on the delta method provides valid standard error estimates regarding coverage rates even though a differentiable approximation of the involved model-robust loss function is used.The first two research questions are addressed in Simulation Study 1 (see Section 6), while the third is investigated in Simulation Study 2 (see Section 7).

Simulation Study 1: Bias and RMSE
This Simulation Study 1 examined the impact of group-specific item intercepts in a multiple-group one-dimensional factor model on bias and RMSE of factor means and factor variances.In the data-generating model (DGM), measurement invariance was violated.That is, differential item functioning (DIF; refs.[45,46]) occurred, and hence, DIF effects in item intercepts were simulated.

Method
The DGM in the simulation study was identical to Simulation Study 2 in [14] and mimicked [21].The data were simulated from a one-dimensional factor model involving five items and three groups.The factor variable η 1 was normally distributed with group means α 1,1 = 0, α 2,1 = 0.3, and α 3,1 = 0.8 and group variances φ 1,11 = 1, φ 2,11 = 1.5, and φ 3,11 = 1.2, respectively.All five factor loadings were set to 1, and all measurement error variances were set to 1 in all groups and uncorrelated with each other.The factor variable and residual variables were normally distributed.
Only a subset of group-specific item intercepts was simulated differently from zero.These nonzero item intercepts indicated measurement noninvariance (i.e., the presence of DIF effects).One of the five items in each group had a DIF effect.However, different items across the three groups were affected by the DIF.In the first group, the fourth item intercepts had a DIF effect δ.In the second group, the first item had a DIF effect −δ, while the second item had a DIF effect −δ in the third group.The DIF effect δ was chosen as 0, 0.3, or 0.6.The value of δ = 0 represented the situation of measurement invariance.The sample size per group was chosen as N = 250, N = 500, N = 1000, or N = 2000.
All analysis models were a multiple-group one-factor model.For identification reasons, the mean of the factor variable in the first group was fixed at 0, and the standard deviation in the first group was fixed at 1. Invariant factor loadings and residual variances were specified across groups.In model-robust moment estimation (ME), we also assumed invariant item intercepts.We utilized the powers p = 0.5, p = 0.25, and p = 0.1 for the loss function ρ ε defined in (8) combined with values of the tuning parameter ε chosen as ε = 10 −2 (=0.01), ε = 10 −3 (=0.001), and ε = 10 −4 (=0.0001).The resulting estimators will be denoted by ME0.5, ME0.25, and ME0.1, in the following Results sections, respectively.Moreover, we used the loss function χ ε defined in (9) with the same ε tuning parameter values as for ρ ε to approximate the L 0 loss function (denoted by ME0 in the following Section 6.2).Furthermore, we used DBIC estimation in which all group-specific item intercepts were allowed to differ across groups.The indicator variables ι k involved in the DBIC approach (see ( 16)) of the model parameter vector θ take only the value 1 for item intercepts.For all other elements of θ, they are set to 0. Therefore, the DBIC approach effectively minimizes the number of estimated item intercepts in its penalty.We also chose the tuning parameter values as ε chosen as 10 −2 , 10 −3 , and 10 −4 .
We did not use the power p = 1 in model-robust moment estimation because it resulted in biased estimates in the presence of DIF [14].However, we included the non-robust ML estimation method for a more comprehensive comparison of the estimation methods.
In total, R = 5000 replications were conducted for all 3 (DIF effect size δ) × 4 (sample size N) = 12 conditions of the simulation study.We investigated the estimation quality of factor means (i.e., α 2,1 and α 3,1 ) and factor variances (i.e., φ 2,11 and φ 3,11 ) in the second and third groups.Bias, RMSE, and relative RMSE were computed to assess the performance of the different estimators.Let θr be a model parameter estimate in replication r = 1, . . ., R.

The bias was estimated by
where θ denotes the true parameter value.The RMSE was estimated by Note that RMSE( θ) ≥ |Bias( θ)| holds because the mean square error (i.e., the square of the RMSE) is the sum of the square of the bias and the variance of an estimator.A relative RMSE can be defined by dividing the RMSE of an estimator by the RMSE of a chosen reference model.To ease the reading of the numeric values of the relative RMSE, the values were multiplied by 100.This quantity can then easily be converted as a percentage gain or loss of a particular estimator compared to a reference model.
The entire simulation study was carried out in the R [47] software (Version 4.3.1).The SEMs were estimated using the sirt::mgsem() function in the R package sirt (Version 4.0-19; ref. [48]).Information about model specification can be found in the material located at https://osf.io/ng6s3(accessed on 25 September 2023).
Researchers interested in fitting a particular model without interest in studying the entire simulation code are referred to the manual help site of the mgsem function in the R package sirt [48].They can type ?sirt::mgsem in the R console (or look at https://alexanderrobitzsch.r-universe.dev/sirt/doc/manual.html#mgsem (accessed on 25 September 2023)) and find an example for applying the ME estimator for an existing dataset.

Results
Figure 3 displays the RMSE of factor mean α 2,1 in the second group for the different estimators for DIF effect sizes δ = 0.3 and δ = 0.6 and sample sizes N = 500 and N = 1000 as a function of the tuning parameter ε.It can be seen that ε = 10 −3 was optimal with respect to the RMSE for ME0.5, ME0.25, and ME0.1, while ε = 10 −2 resulted in the smallest RMSE for ME0 and DBIC.The findings were very similar for the factor mean α 3,1 in the third group and the sample sizes N = 250 and N = 2000.Table 1 displays the bias and relative RMSE of factor means and factor variances as a function of the DIF effect size and the sample size for the different estimation methods.We chose the tuning parameter ε = 10 −3 for the estimators ME0.5, ME0.25, and ME0.1, and ε = 10 −2 for the estimators ME0 and DBIC because they resulted in the smallest RMSE of factor means according to the findings of Figure 3.
Overall, all estimators except for ML were approximately unbiased for factor variances in all conditions and for factor means in the absence of DIF effects (i.e., δ = 0 and measurement invariance holds) except for the sample size N = 250.Slightly biased estimates were obtained for ME with a larger p, such as p = 0.5 in ME0.5.However, the bias decreased with increasing sample size.The ME0 and DBIC were unbiased across all conditions, followed by the estimators ME0.1, ME0.25, and ME0.5, with increasing absolute bias.There was substantial bias for estimated factor means for ML estimation, while the estimates of factor variances were approximately unbiased for ML.Note.Par = parameter; α g,1 = factor mean in group g = 2, 3; φ g,11 = factor variance in group g = 2, 3; ML = maximum likelihood estimation; MEp = robust moment estimation with p = 0.5 (with ε = 0.001), p = 0.25 (with ε = 0.001), p = 0.1 (with ε = 0.001), or p = 0 (with ε = 0.01); DBIC = direct BIC minimization (with ε = 0.01).Absolute biases larger than 0.015 are shown with a gray background.† ME0.5 is the reference method in computation of the relative RMSE.Relative RMSE values smaller than 98.0 are printed in bold font.Relative RMSE values larger than 102.0 are show with a gray background.
To compare the estimation accuracy in terms of the relative RMSE, we chose ME0.5 (with ε = 10 −3 ) as the reference model.Hence, the relative RMSE was 100 for this method in Table 1.It turned out that ME0 and DBIC were superior to the other estimators with non-negligible efficiency gains.For example, for the factor mean α 2,1 in the second group, δ = 0.3, and N = 1000, the efficiency gain in terms of the RMSE was 6.4% (=100 − 93.6) for ME0 and 6.5% for DBIC.Across all conditions, no noteworthy differences between ME0 and DBIC estimators were found.
ML estimates were slightly more efficient than ME estimates in the absence of DIF.However, the efficiency gains of ML decrease with increasing sample size.
To conclude, it seems promising to replace the L p loss function using the power p = 0.5 (i.e., ME0.5) with the L 0 loss function implemented in ME0 or regularized estimation in DBIC for large sample sizes N = 1000 or N = 2000.Similar performance of the different estimators was found for N = 500, while the power p = 0.5 was the frontrunner for N = 250.

Simulation Study 2: Coverage
The second Simulation Study 2 investigated the assessment of coverage rates of the model-robust moment estimators and the DBIC method.

Method
The same DGM as in Simulation Study 1 was employed to simulate data (see Section 6.1).We evaluated the standard error computation described in Section 4 for the five estimators using the tuning parameter ε that resulted in the smallest RMSE of factor means.In detail, we used ε = 10 −3 for the estimators ME0.5, ME0.25, and ME0.1; and we used ε = 10 −2 for the estimators ME0 and DBIC (see Figure 3).Moreover, we included ML estimation for reasons of comparison.
The same analysis models as in Simulation Study 1 were specified.Confidence intervals at the confidence level of 95% were computed using a normal distribution approximation (i.e., the estimated confidence interval was θ ± 1.96 × SE( θ)).The coverage rate at the confidence level of 95% was computed as the percentage of the events that a computed confidence interval covers the true parameter value.

Discussion
In this article, we compared model-robust moment estimation with a recently proposed variant of regularized ML estimation by O'Neill and Burke [15] that directly maximizes the BIC (i.e., the DBIC estimator).In the DBIC estimation, these authors suggested a differentiable approximation of the L 0 loss function, which was also used in model-robust moment estimation.Interestingly, the L 0 loss function outperformed L p loss functions for p > 0 regarding bias and RMSE.Furthermore, model-robust moment estimation with the L 0 loss function performed very similarly to the DBIC estimator.Moreover, the estimation of standard errors was successfully implemented for all estimators because coverage rates were acceptable for all parameters in all simulation conditions.
In line with previous studies, we anticipate that sample sizes must be sufficiently large in order to achieve the model-robustness properties of the L p and L 0 loss functions.If the sample size is too small (e.g., N = 100 subjects in a multiple-group SEM analysis), the sampling error in moments that are used as sufficient statistics in the SEM can exceed model errors (i.e., unmodeled group-specific item intercepts).In this case, model-robust methods are not expected to perform well.In fact, Simulation Study 1 revealed that p = 0.5 is preferable to other L p loss functions with p < 0.5 or the L 0 loss function for N = 250, while the situation changes with larger sample sizes such as N = 1000.
Simulation Study 1 indicated that the tuning parameter ε = 0.001 in the differentiable approximation of the nondifferentiable L p loss function should be used for the L p loss function for p = 0.5, p = 0.25, or p = 0.1.In contrast, ε = 0.01 was found optimal for the L 0 loss function and the direct BIC minimization estimation approach in regularized estimation.We expect that these findings will transfer to other models that involve standardized variables.In our experience from regularized estimation and the invariance alignment approach [21], using a value of ε (such as ε = 10 −5 and smaller values) tuning parameter that is too small is prohibitive because it more likely results in convergence issues in the local optima of the fit function.
In our simulation studies, we only considered model errors in item intercepts.Future simulation studies could also investigate model errors in the covariance structure (i.e., unmodeled residual correlations).Of course, model errors in the mean and the covariance structure can also be simultaneously examined.
A requirement to achieve model robustness of SEM estimators is that model errors in the mean or covariance structure are sparsely distributed.That is, only a few entries are allowed to differ from zero, while the majority of model errors must be (approximately) zero.If factor loadings were not invariant across groups, more densely distributed model errors would result.As a consequence, model-robust moment estimation will likely not work, and regularized maximum likelihood estimation might be preferred.However, the DBIC minimization method for regularized maximum likelihood estimation can also be utilized in this case if the number of group-specific factor loadings is counted in the BIC penalty term.
In this article, we only applied the L 0 and L p loss functions to continuous items.However, the principle directly transfers to SEMs of ordinal data that are based on fitting thresholds and polychoric correlations [50] instead of means and covariances for continuous data, respectively.Moreover, the model-robust estimators could also be applied to two-step estimation methods of multilevel structural equation models [51,52].
To sum up, model-robust SEM estimators based on the L p (for p < 1) and L 0 loss functions are attractive to researchers who do not want model estimates being influenced by the presence of a few model deviations (i.e., model errors).In contrast, usually employed (non-robust) SEM estimators such as maximum likelihood estimation are impacted by model errors.In this sense, misfitting models do not necessarily result in biased estimates when using model-robust estimation.

Table 1 .
Simulation Study 1: Bias and relative root mean square error (RMSE) of factor means and factor variances as a function of DIF effect size δ and sample size N.

Table 2 .
Simulation Study 2: Coverage rates (in percentages) of factor means and factor variances as a function of DIF effect size δ and sample size N.