Single- and Multiple-Group Penalized Factor Analysis: A Trust-Region Algorithm Approach with Integrated Automatic Multiple Tuning Parameter Selection

Elena Geminiani; Giampiero Marra; Irini Moustaki

doi:10.1007/s11336-021-09751-8

Single- and Multiple-Group Penalized Factor Analysis: A Trust-Region Algorithm Approach with Integrated Automatic Multiple Tuning Parameter Selection

Published online by Cambridge University Press: 01 January 2025

Elena Geminiani

Giampiero Marra and

Irini Moustaki

Show author details

Elena Geminiani*: Affiliation:
University of Bologna
Giampiero Marra: Affiliation:
University College London
Irini Moustaki: Affiliation:
London School of Economics and Political Science
*: Correspondence should be made to Elena Geminiani, Department of Statistical Sciences, University of Bologna, Via Delle Belle Arti 41, 40126 Bologna, Italy. Email: elena.geminiani4@unibo.it; geminianielena@gmail.com

Article contents

Abstract
Introduction
The Normal Linear Factor Analysis Model
Sparsity-Inducing Penalties
The Multiple-Group Factor Analysis Model
Sparsity and Invariance-Inducing Penalties
Generalized Information Criterion
Penalized Maximum Likelihood Estimation
Simulation Studies
Empirical Application
Discussion
Funding
Footnotes
References

Rights & Permissions

Abstract

Penalized factor analysis is an efficient technique that produces a factor loading matrix with many zero elements thanks to the introduction of sparsity-inducing penalties within the estimation process. However, sparse solutions and stable model selection procedures are only possible if the employed penalty is non-differentiable, which poses certain theoretical and computational challenges. This article proposes a general penalized likelihood-based estimation approach for single- and multiple-group factor analysis models. The framework builds upon differentiable approximations of non-differentiable penalties, a theoretically founded definition of degrees of freedom, and an algorithm with integrated automatic multiple tuning parameter selection that exploits second-order analytical derivative information. The proposed approach is evaluated in two simulation studies and illustrated using a real data set. All the necessary routines are integrated into the R package penfa.

Keywords

effective degrees of freedom generalized information criterion measurement invariance penalized likelihood simple structure

Information

Type: Theory and Methods
Information: Psychometrika , Volume 86 , Issue 1 , March 2021 , pp. 65 - 95

DOI: https://doi.org/10.1007/s11336-021-09751-8 [Opens in a new window]
Creative Commons: This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Copyright: Copyright © 2021 The Author(s)

1. Introduction

Factor analysis has been extensively applied in the social, behavioral and natural sciences as a data reduction method. For a given set of observed variables $x_{1}, \dots, x_{p}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_1, \ldots , x_p$$\end{document} one would like to find a set of latent factors $f_{1}, \dots, f_{r}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_1, \ldots , f_r$$\end{document} , fewer in number than the observed variables ( $r < p$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r < p$$\end{document} ), that contain essentially the same information. Factor analysis can be conducted in an exploratory (EFA; Mulaik, Reference Mulaik2009) or confirmatory (CFA; Jöreskog, Reference Jöreskog, Jöreskog, Sörbom and Magidson1979) way. EFA analyzes a set of correlated observed variables without knowing in advance either the number of factors that are required to explain their interrelationships or their meaning. Depending on the r-factor model finally chosen (based on goodness-of-fit criteria and fit measures) as well as the rotation applied, an interpretation and labelling of the factors are given. CFA postulates certain relationships among the observed and latent variables by assuming a pre-specified pattern for the model parameters (factor loadings, structural parameters, unique variances). It is mainly concerned with testing hypotheses about the values of the factor loadings (usually, that some of them are zero).

In data reduction techniques such as factor analysis, the interest is in obtaining factor solutions that exhibit a “simple structure” (Thurstone, Reference Thurstone1947), that is, with many zero loadings and pure measures (i.e., each variable loads only on a single factor). In EFA this is accomplished with orthogonal or oblique factor rotations. However, rotations often do not generate loadings precisely equal to zero, so users have to manually set to zero those loadings that are smaller than a threshold (e.g., 0.30; Hair et al., Reference Hair, Black, Babin and Anderson2010). Secondly, because each rotation is based on a specific optimization criterion, different rotations often lead to different factor structures which may all be far from “simple”. In CFA, one usually resorts to modification indices (Chou & Huh, Reference Chou, Huh and Hoyle2012) instead, but, if used extensively, they can lead to higher risks of capitalization on chance (MacCallum et al., Reference MacCallum, Roznowski and Necowitz1992), and a lower probability of finding the best model specification (Chou & Bentler, Reference Chou and Bentler1990).

Penalized factor analysis is an alternative technique that produces parsimonious models using largely an automated procedure. The resulting models are less prone to instability in the estimation process and are easier to interpret and generalize than their unpenalized counterparts. It is based on the use of penalty functions that allow a subset of the model parameters (typically the factor loadings) to be automatically set to zero. The penalty is usually non-differentiable (Fan & Li, Reference Fan and Li2001), so that it produces a sparse factor structure, that is, a loading matrix where the number of nonzero entries is much smaller than the total number of its elements. This definition does not impose any pattern on the nonzero entries, so a simple structure is not enforced if it is not supported by the data. These sparsity-inducing penalties can reduce model complexity, enhance the interpretability of the results, and produce more stable parameter estimates. These benefits come, however, with a loss in model fit (i.e., a nonzero bias), so it is crucial to balance goodness of fit and sparsity appropriately. This can be achieved via the selection of a tuning parameter, which controls the amount of sparsity introduced in the model. A grid-search over a range of tuning values is generally conducted, and the optimal model is picked on the basis of information criteria or cross-validation.

In the last few years, several works have applied penalized estimation and regularization methods to models with latent variables. Choi, Oehlert and Zou (Reference Choi, Oehlert and Zou2010) used lasso (“least absolute shrinkage and selection operator”; Tibshirani, Reference Tibshirani1996) and adaptive lasso penalties in EFA. Since the lasso leads to biased estimates and overly dense factor structures, Hirose and Yamamoto (Reference Hirose and Yamamoto2014a; Reference Hirose and Yamamotob) employed non-convex penalties, such as the scad (“smoothly clipped absolute deviation”) and the mcp (“minimax concave penalty”). Trendafilov, Fontanella and Adachi (Reference Trendafilov, Fontanella and Adachi2017) penalized a reparametrized loading matrix, whereas Jin, Moustaki and Yang-Wallentin (Reference Jin, Moustaki and Yang-Wallentin2018) considered a quadratic approximation of the objective function. Regularized methods have also been applied to structural equation models (SEM) for which CFA is a special case. Jacobucci, Grimm and McArdle (Reference Jacobucci, Grimm and McArdle2016) developed the regularized SEM (RegSEM) using a reticular action model formulation and coordinate descent or general optimization routines. Huang, Chen and Weng (Reference Huang, Chen and Weng2017) and Huang (Reference Huang2020) examined the same problem of penalizing a SEM but employed a modification of the quasi-Newton algorithm.

Penalized estimation can be also extended to multiple-group analyses, such as cross-national surveys or cross-cultural assessments in psychological or educational testing. Recently, Huang (Reference Huang2018) and Lindstrøm and Dahl (Reference Lindstrøm and Dahl2020) developed a penalized approach for multiple-group SEM, showing the benefits of using regularization techniques as alternatives to factorial invariance testing procedures (Meredith, Reference Meredith1993) to ascertain the differences and similarities of the parameter estimates across groups (see Bauer, Belzak & Cole, Reference Bauer, Belzak and Cole2020 for a regularized approach for moderated non-linear factor analysis).

In this paper, we propose a penalized-estimation strategy for single- and multiple-group factor analysis models based on a carefully structured trust-region algorithm. The penalized optimization problem requires the availability of second-order analytical derivative information and thus twice-continuously differentiable functions. Because a sparse solution can be only achieved with non-differentiable penalties, we employ differentiable approximations of them. In particular, we locally approximate several convex and non-convex penalties, including lasso, adaptive lasso, scad and mcp. We also provide a theoretically founded definition of degrees of freedom (required when performing model selection) and present an efficient automatic procedure for the estimation of the tuning parameters, hence eliminating the need for computationally intensive grid searches as done in the literature. The proposed methodology is integrated into the R package penfa (a short form for PENalized Factor Analysis).

The paper is organized as follows. Sect. 2 briefly discusses the classical linear factor analysis model. In Sect. 3 we review and develop penalized likelihood estimation via locally approximated penalties. The extension of the model and the penalized approach for the case of multiple groups are described in Sects. 4 and 5, respectively. The derivation of the model degrees of freedom is presented in Sect. 6. Parameter estimation and the automatic selection of the tuning parameters are detailed in Sect. 7. The performance of the model is evaluated in two simulation studies (Sect. 8) and an empirical application (Sect. 9). Lastly, Sect. 10 concludes the paper and gives directions for future research. Additional details can be found in the Online Resources.

2. The Normal Linear Factor Analysis Model

The classical linear factor analysis model takes the formFootnote 1:

(1)

\begin{matrix} x = Λ f + ε, \end{matrix}

where $x$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{x}}$$\end{document} is the $p \times 1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p \times 1$$\end{document} vector of observed variables, $Λ$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\Lambda }}$$\end{document} is the $p \times r$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p \times r$$\end{document} factor loading matrix, $f$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{f}}$$\end{document} is the $r \times 1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r \times 1$$\end{document} vector of common factors, and $ε$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\varepsilon }}$$\end{document} is the $p \times 1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p \times 1$$\end{document} vector of unique factors. It is assumed that $f \sim N (0, Φ)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{f}} \sim {\mathcal {N}}({\varvec{0}}, {\varvec{\Phi }})$$\end{document} , $ε \sim N (0, Ψ)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\varepsilon }} \sim {\mathcal {N}}({\varvec{0}}, {\varvec{\Psi }})$$\end{document} , and $f$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{f}}$$\end{document} is independent of $ε$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\varepsilon }}$$\end{document} . The observed variables are assumed to be conditionally independent (i.e., $Ψ$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\Psi }}$$\end{document} is a diagonal matrix), although this assumption can be relaxed if required. It then follows that $x \sim N (0, Σ)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{x}} \sim {\mathcal {N}}({\varvec{0}}, {\varvec{\Sigma }})$$\end{document} , where the model-implied covariance matrix is $Σ = Λ Φ Λ^{T} + Ψ$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\Sigma }} = {\varvec{\Lambda }} {\varvec{\Phi }} {\varvec{\Lambda }}^T + {\varvec{\Psi }}$$\end{document} .

It is possible to fix certain elements in $Λ, Φ$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\Lambda }}, {\varvec{\Phi }}$$\end{document} and $Ψ$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\Psi }}$$\end{document} to zero based on a data generating hypothesis. The remaining $m \leq min (N, \frac{p (p + 1)}{2})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m \le \min \left( N, \frac{p(p+1)}{2}\right) $$\end{document} elements, with N the total sample size, constitute the free parameters in $vec (Λ)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {vec}({\varvec{\Lambda }})$$\end{document} , $diag (Ψ)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {diag}({\varvec{\Psi }})$$\end{document} , and $vech (Φ)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {vech}({\varvec{\Phi }})$$\end{document} , and are collected in the vector $θ$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\theta }}$$\end{document} , where the vec( $\cdot$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\cdot $$\end{document} ) operator converts the enclosed matrix into a vector by stacking its columns, $diag (\cdot)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {diag}(\cdot )$$\end{document} extracts the diagonal elements of the enclosed square matrix, and vech( $\cdot$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\cdot $$\end{document} ) vectorizes the lower-diagonal part of the enclosed symmetric matrix. As it is common practice in these cases, we assume that the observed variables are measured as deviations from their means, so that the parameters only strive to reproduce the covariance matrix. As in Jöreskog (Reference Jöreskog, Jöreskog, Sörbom and Magidson1979), we fix the variances of the common factors to unity for scale setting, and $r - 1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r-1$$\end{document} elements of $Λ$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\Lambda }}$$\end{document} , in each column, to zero for uniqueness under factor rotation.

For a random sample of size N the log-likelihood is written as

(2)

\begin{matrix} ℓ (θ) = - \frac{N}{2} (log | Σ | + tr (S Σ^{- 1}) + p log (2 π)), \end{matrix}

where $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{S}}$$\end{document} is the sample covariance matrix. Since we are interested in introducing sparsity in the factor loading matrix, the estimation of the factor model will involve penalized likelihood procedures. The next section illustrates how such sparsity-inducing penalty functions can be specified and suitably approximated.

3. Sparsity-Inducing Penalties

Since the primary interest of factor analysis is a sparse loading matrix, penalization is imposed on the factor loading matrix $Λ$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\Lambda }}$$\end{document} . Let us write the parameter vector as $θ = {(θ_{1}, \dots, θ_{q^{⋆}}, θ_{q^{⋆} + 1}, \dots, θ_{m})}^{T}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\theta }} = (\theta _1, \ldots , \theta _{q^\star }, \theta _{q^\star +1}, \ldots , \theta _m)^T$$\end{document} , where the sub-vector ${(θ_{1}, \dots, θ_{q^{⋆}})}^{T}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\theta _1, \ldots , \theta _{q^\star })^T$$\end{document} collects the penalized parameters (i.e., the factor loadings), whereas ${(θ_{q^{⋆} + 1}, \dots, θ_{m})}^{T}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\theta _{q^\star +1}, \ldots , \theta _m)^T$$\end{document} the unpenalized parameters (i.e., the free elements in $Ψ$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\Psi }}$$\end{document} and $Φ$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\Phi }}$$\end{document} ). Because of the presence of fixed elements in $Λ$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\Lambda }}$$\end{document} (Sect. 2), the number of penalized factor loadings $q^{⋆}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q^\star $$\end{document} is smaller than $p \times r$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p\times r$$\end{document} . Define $R_{q} = diag (0, 0, \dots, 0, 1, 0, \dots, 0)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{R}}_q = \text {diag}(0, 0, \ldots , 0, 1, 0, \ldots , 0)$$\end{document} a diagonal matrix where the 1 on the ${(q, q)}^{th}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(q, q)^{\text {th}}$$\end{document} entry of the matrix corresponds to the $q^{th}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q^{\text {th}}$$\end{document} parameter in $θ$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\theta }}$$\end{document} , for $q = 1, \dots, q^{⋆}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q=1, \ldots , q^\star $$\end{document} , and $R_{q} = O_{m \times m}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{R}}_q = {\varvec{O}}_{m \times m}$$\end{document} for $q = q^{⋆} + 1, \dots, m$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q = q^\star +1 , \ldots , m$$\end{document} . Let $P_{η} (θ)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {P}}_\eta ({\varvec{\theta }})$$\end{document} be a penalty function on the parameter vector $θ$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\theta }}$$\end{document} , where $η \in [0, \infty)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta \in [0, \infty ) $$\end{document} is a positive tuning parameter which determines the amount of shrinkage or penalization. The overall penalty is then given by the sum of the penalty terms for each parameter, that is, $P_{η} (θ) = \sum_{q = 1}^{m} P_{η, q} (| | R_{q} {θ | |}_{1})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {P}}_\eta ({\varvec{\theta }}) = \sum _{q=1}^{m} {\mathcal {P}}_{\eta ,q} (||{\varvec{R}}_q {\varvec{\theta }}||_1)$$\end{document} , where $| | R_{q} {θ | |}_{1} = | θ_{q} |$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$||{\varvec{R}}_q {\varvec{\theta }}||_1 = |\theta _q |$$\end{document} if $q = 1, \dots, q^{⋆}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q=1, \ldots , q^\star $$\end{document} , and zero otherwise. An example clarifying the formulation of this penalty is provided in Section B.1.1. One of the best-known penalties is the lasso (Tibshirani, Reference Tibshirani1996), which is defined as

(3)

\begin{matrix} P_{η}^{L} (θ) = η \sum_{q = 1}^{q^{⋆}} | θ_{q} | . \end{matrix}

The potential problem with this penalty is that it penalizes all parameters equally, and thus can either select an overly complicated model or over-shrink large parameters. An ideal penalty should induce weak shrinkage on large effects and strong shrinkage on irrelevant effects (Tang, Shen, Zhang & Yi, Reference Tang, Shen, Zhang and Yi2017). To address this issue, alternative penalties have been developed, the most common ones being the adaptive lasso (alasso; Zou, Reference Zou2006), scad (Fan & Li, Reference Fan and Li2001) and mcp (Zhang, Reference Zhang2010). These penalties give different amounts of shrinkage to each parameter, so each factor loading is weighted differently. Because of this, they lead to sparser solutions and enjoy the so-called “oracle” property, that is, when the true parameters have some zero loadings, they are estimated as zero with probability tending to one, and the nonzero loadings are estimated as well as when the correct submodel is known (Fan & Li, Reference Fan and Li2001). The alasso is defined as

(4)

\begin{matrix} P_{η}^{A} (θ) = η \sum_{q = 1}^{q^{⋆}} w_{q} | θ_{q} | = η \sum_{q = 1}^{q^{⋆}} \frac{| θ_{q} |}{| {\hat{θ}}_{q} |^{a}} for a > 0 . \end{matrix}

It uses an adaptive weighting scheme based on a set of available weights $w_{q} = \frac{1}{| {\hat{θ}}_{q} |^{a}} (q = 1, \dots, q^{⋆})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w_q = \dfrac{1}{|{\hat{\theta }} _q |^{a}} \, (q = 1, \ldots , q^\star )$$\end{document} , which are often taken to be the maximum likelihood estimates, that is, $w_{q} = \frac{1}{| {\hat{θ}}_{q}^{MLE} |^{a}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w_q = \dfrac{1}{|{\hat{\theta }}_q^{\text {MLE}}|^{a}}$$\end{document} . As the exponent a gets larger, the relative strength of the penalization increases for smaller maximum likelihood estimates compared to larger maximum likelihood estimates.

Similarly, the scad and mcp use a varying weighting scheme. The scad is defined as

(5)

and the mcp as

(6)

where a is an additional tuning parameter. The superscripts L, A, S, M in equations (3)-(6) refer to the lasso, alasso, scad and mcp, respectively. The derivations of expressions (3)-(6) can be found in Section B.1.2. While the lasso and alasso are convex penalties, the scad and mcp are non-convex and can, therefore, make the optimization problem non-convex. In fact, a challenge with non-convex penalties is to find a good balance between sparsity and stability. To this end, both scad and mcp have an extra tuning parameter (a) which regulates their concavity so that, when it exceeds a threshold, the optimization problem becomes convex.

The above penalties help to obtain sparse solutions, however, they are non-differentiable, which is problematic for developing a coherent computational and theoretical inferential framework. The next section addresses this issue by replacing the non-differentiable penalties with their differentiable counterparts obtained via local approximations.

3.1. Locally Approximated Penalties

Ulbricht (Reference Ulbricht2010) pointed out that a good penalty function should satisfy the following properties, for $q = 1, \dots, m$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q = 1, \ldots , m$$\end{document} : (P.1) $P_{η, q} : R^{+} \to R^{+}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {P}}_{\eta ,q}: {\mathbb {R}}^+ \rightarrow {\mathbb {R}}^+$$\end{document} and $P_{η, q} (0) = 0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {P}}_{\eta ,q}(0) = 0$$\end{document} ; (P.2) $P_{η, q} (| | R_{q} {θ | |}_{1})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {P}}_{\eta ,q}(||{\varvec{R}}_q {\varvec{\theta }}||_1)$$\end{document} continuous and strictly monotone in $| | R_{q} θ {| |}_{1}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$||{\varvec{R}}_q {\varvec{\theta }}||_1$$\end{document} ; (P.3) $P_{η, q} (| | R_{q} {θ | |}_{1})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {P}}_{\eta ,q}(||{\varvec{R}}_q {\varvec{\theta }}||_1)$$\end{document} continuously differentiable $\forall | | R_{q} θ {| |}_{1} \neq 0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\forall \, ||{\varvec{R}}_q {\varvec{\theta }}||_1 \ne 0$$\end{document} , such that $\frac{\partial P_{η, q} (| | R_{q} {θ | |}_{1})}{\partial | | R_{q} θ {| |}_{1}} > 0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dfrac{\partial {\mathcal {P}}_{\eta ,q} (||{\varvec{R}}_q {\varvec{\theta }}||_1)}{\partial ||{\varvec{R}}_q {\varvec{\theta }}||_1} > 0$$\end{document} . We develop differentiable approximations of the above penalties that satisfy these properties. These approximations make the objective function differentiable, which is an indispensable prerequisite for the theoretical derivation of the degrees of freedom of the model, and a computationally and theoretically founded estimation framework (Sects. 6–7). In the same spirit, as for instance, Filippou, Marra and Radice (Reference Filippou, Marra and Radice2017), we locally approximate the non-differentiable $L_{1}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_1$$\end{document} -norms in (3)-(6) at $| | R_{q} θ {| |}_{1} = 0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$||{\varvec{R}}_q {\varvec{\theta }}||_1 = 0$$\end{document} and combine this with ideas by Fan and Li (Reference Fan and Li2001) and Ulbricht (Reference Ulbricht2010). Let $| | R_{q} {θ | |}_{1} = | | ξ_{q} {| |}_{1}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$||{\varvec{R}}_q {\varvec{\theta }}||_1 = ||{\varvec{\xi }}_q||_1 $$\end{document} , where the $q^{th}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q^{\text {th}}$$\end{document} element in $ξ_{q} = {(0, \dots, 0, θ_{q}, 0, \dots, 0)}^{T}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\xi }}_q = (0, \ldots , 0, \theta _q, 0, \ldots , 0)^T$$\end{document} corresponds to the $q^{th}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q^{\text {th}}$$\end{document} parameter in $θ$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\theta }}$$\end{document} . Assume that an approximation $K_{1} (ξ_{q}, A)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {K}}_1({\varvec{\xi }}_q, {\mathcal {A}})$$\end{document} of the $L_{1}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_1$$\end{document} -norm $| | \cdot {| |}_{1}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$||\cdot ||_1$$\end{document} exists such that

\begin{matrix} | | ξ_{q} {| |}_{1} = K_{1} (ξ_{q}, B) = lim_{A \to B} K_{1} (ξ_{q}, A), \end{matrix}

where $A$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {A}}$$\end{document} represents a set of possible tuning parameters, $B$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {B}}$$\end{document} is the set of boundary values for $| | ξ_{q} {| |}_{1}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$||{\varvec{\xi }}_q||_1$$\end{document} and $K_{1} (ξ_{q}, A)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {K}}_1({\varvec{\xi }}_q, {\mathcal {A}})$$\end{document} is at least twice differentiable. We use $| | ξ_{q} {| |}_{1} = K_{1} (ξ_{q}, A) = {(ξ_{q}^{T} ξ_{q} + \bar{c})}^{\frac{1}{2}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$||{\varvec{\xi }}_q||_1 = {\mathcal {K}}_1({\varvec{\xi }}_q, {\mathcal {A}}) = ({\varvec{\xi }}_q^T {\varvec{\xi }}_q + \bar{c})^{\frac{1}{2}}$$\end{document} (Koch, Reference Koch1996), with $\bar{c}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{c}$$\end{document} a small positive real number (e.g., $10^{- 8}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-8}$$\end{document} ) which controls the closeness between the approximation and the exact function. For all $ξ_{q}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\xi }}_q$$\end{document} for which the derivative $\frac{\partial | | ξ_{q} {| |}_{1}}{\partial ξ_{q}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dfrac{\partial ||{\varvec{\xi }}_q||_1}{\partial {\varvec{\xi }}_q}$$\end{document} is defined, we assume that

\begin{matrix} \frac{\partial | | ξ_{q} {| |}_{1}}{\partial ξ_{q}} = \frac{\partial K_{1} (ξ_{q}, B)}{\partial ξ_{q}} = lim_{A \to B} D_{1} (ξ_{q}, A), \end{matrix}

where $D_{1} (ξ_{q}, A) = \frac{\partial K_{1} (ξ_{q}, A)}{\partial ξ_{q}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {D}}_1({\varvec{\xi }}_q, {\mathcal {A}}) = \dfrac{\partial {\mathcal {K}}_1({\varvec{\xi }}_q, {\mathcal {A}})}{\partial {\varvec{\xi }}_q}$$\end{document} , and that $D_{1} (0, A) = 0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {D}}_1 ({\varvec{0}}, {\mathcal {A}}) = {\varvec{0}}$$\end{document} . Then, the first derivative $D_{1} (ξ_{q}, A) = {(ξ_{q}^{T} ξ_{q} + \bar{c})}^{- \frac{1}{2}} ξ_{q}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\mathcal {D}}_1({\varvec{\xi }}_q, {\mathcal {A}}) = ({\varvec{\xi }}_q^T {\varvec{\xi }}_q + \bar{c})^{-\frac{1}{2}} {\varvec{\xi }}_q$$\end{document} is a continuous approximation of the first-order derivative of the $L_{1}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_1$$\end{document} -norm. Notice that $K_{1} (ξ_{q}, A)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {K}}_1({\varvec{\xi }}_q, {\mathcal {A}})$$\end{document} deviates only slightly from $K_{1} (ξ_{q}, B)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {K}}_1({\varvec{\xi }}_q, {\mathcal {B}})$$\end{document} : when $ξ_{q} = 0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\xi }}_q = {\varvec{0}}$$\end{document} the deviation is $\sqrt{\bar{c}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{\bar{c}}$$\end{document} , whereas for any other value of $ξ_{q}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\xi }}_q$$\end{document} the deviation is less than $\bar{c}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{c}$$\end{document} .

Penalty $P_{η}^{T} (θ)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {P}}^{{\mathcal {T}}}_\eta ({\varvec{\theta }})$$\end{document} for $T = {L, A, S, M}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {T}} = \{L, A, S, M\}$$\end{document} can be locally approximated by a quadratic function as follows. Suppose that $\tilde{θ}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\varvec{\theta }}}$$\end{document} is an initial value close to the true value of $θ$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\theta }}$$\end{document} . Then, we approximate $P_{η}^{T} (θ)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {P}}^{{\mathcal {T}}}_\eta ({\varvec{\theta }})$$\end{document} by a Taylor expansion of order one at $\tilde{θ}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\varvec{\theta }}}$$\end{document} , that is,

(7)

\begin{matrix} P_{η}^{T} (θ) \approx P_{η}^{T} (\tilde{θ}) + \nabla_{\tilde{θ}} P_{η}^{T} {(\tilde{θ})}^{T} (θ - \tilde{θ}), \end{matrix}

where $\nabla_{\tilde{θ}} P_{η}^{T} (\tilde{θ}) = \frac{\partial P_{η}^{T} (\tilde{θ})}{\partial \tilde{θ}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nabla _{{\tilde{\varvec{\theta }}}} {\mathcal {P}}^{{\mathcal {T}}}_\eta ({\tilde{\varvec{\theta }}}) = \dfrac{\partial {\mathcal {P}}^{{\mathcal {T}}}_\eta ({\tilde{\varvec{\theta }}})}{\partial {\tilde{\varvec{\theta }}}}$$\end{document} . As proved in Section B.1.3, $P_{η}^{T} (θ)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {P}}^{{\mathcal {T}}}_\eta ({\varvec{\theta }})$$\end{document} is approximated as

The penalty matrix

is an $m \times m$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m \times m$$\end{document} block diagonal matrix of the form:

(8)

The first block is composed of the $q^{⋆} \times q^{⋆}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q^\star \times q^\star $$\end{document} diagonal matrix

and corresponds to the parameters to penalize, whereas the second block is an $(m - q^{⋆})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(m-q^\star )$$\end{document} -dimensional null matrix relative to the parameters unaffected by the penalization. The matrix

is in turn a diagonal matrix whose entries $m_{q}^{T} = \frac{\partial P_{η, q}^{T} (| | R_{q} \tilde{θ} {| |}_{1})}{\partial | | R_{q} \tilde{θ} {| |}_{1}} \frac{1}{\sqrt{{(R_{q} \tilde{θ})}^{T} R_{q} \tilde{θ} + \bar{c}}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_q^{{\mathcal {T}}} = \dfrac{\partial {\mathcal {P}}^{{\mathcal {T}}}_{\eta ,q} (||{\varvec{R}}_q {\tilde{\varvec{\theta }}}||_1)}{\partial ||{\varvec{R}}_q {\tilde{\varvec{\theta }}}||_1}\dfrac{1}{\sqrt{({\varvec{R}}_q {\tilde{\varvec{\theta }}})^T {\varvec{R}}_q {\tilde{\varvec{\theta }}} + \bar{c}}}$$\end{document} (for $q = 1, \dots, q^{⋆}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q = 1, \ldots , q^\star $$\end{document} ) determine the amount of shrinkage on ${\tilde{θ}}_{q}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\theta }}_q$$\end{document} controlled by the tuning $η$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} and required by penalty $T$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {T}}$$\end{document} . Their expressions for the lasso, alasso, scad and mcp are (see Section B.1.3.1)

(9)

(10)

(11)

(12)

4. The Multiple-Group Factor Analysis Model

In studies of multiple groups of respondents, such as cross-national surveys and cross-cultural assessments in psychological or educational testing, the interest often lies in the comparisons of the groups with respect to their factor structures. In this case, the model becomes

(13)

\begin{matrix} x_{g} = τ_{g} + Λ_{g} f_{g} + ε_{g} for g = 1, \dots, G, \end{matrix}

where the subscript g denotes the group, and $τ_{g}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\tau }}_g$$\end{document} the intercept terms. It is assumed that $f_{g} \sim N (κ_{g}, Φ_{g})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{f}}_g \sim {\mathcal {N}}({\varvec{\kappa }}_g, {\varvec{\Phi }}_g)$$\end{document} , $ε_{g} \sim N (0, Ψ_{g})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\varepsilon }}_g \sim {\mathcal {N}}({\varvec{0}}, {\varvec{\Psi }}_g)$$\end{document} , $f_{g}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{f}}_g$$\end{document} is independent of $ε_{g}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\varepsilon }}_g$$\end{document} , and $Ψ_{g}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\Psi }}_g$$\end{document} is a diagonal matrix. Then, it follows that $x_{g} \sim N (μ_{g}, Σ_{g})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{x}}_g \sim {\mathcal {N}}({\varvec{\mu }}_g, {\varvec{\Sigma }}_g)$$\end{document} , where the model-implied moments are $μ_{g} = τ_{g} + Λ_{g} κ_{g}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\mu }}_g = {\varvec{\tau }}_g + {\varvec{\Lambda }}_g {\varvec{\kappa }}_g$$\end{document} and $Σ_{g} = Λ_{g} Φ_{g} Λ_{g}^{T} + Ψ_{g}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\Sigma }}_g = {\varvec{\Lambda }}_g {\varvec{\Phi }}_g {\varvec{\Lambda }}_g^T + {\varvec{\Psi }}_g$$\end{document} .

We set the metric of the factors and the necessary identification restrictions through the “marker-variable” approach (Little, Slegers & Card, Reference Little, Slegers and Card2006), which relies on the selection of a representative variable (marker) for each factor in each group. Then, we fix the intercepts of the markers to zero, the loadings on the “marked” factors to 1.0, and those on the remaining factors to zero. All of the other parameters are estimated. The choice of the markers is crucial and should be an accurate one (Millsap, Reference Millsap2001). Alternative identification methods are discussed in Millsap (Reference Millsap2012).

The free parameters of each group appearing in $vec (Λ_{g})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {vec}({\varvec{\Lambda }}_g)$$\end{document} , $τ_{g}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\tau }}_g$$\end{document} , $diag (Ψ_{g})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {diag}({\varvec{\Psi }}_g)$$\end{document} , $vech (Φ_{g})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {vech}({\varvec{\Phi }}_g)$$\end{document} , and $κ_{g}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\varvec{\kappa }}_g$$\end{document} are collected in the $m_{g}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_g$$\end{document} -dimensional vector $θ_{g}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\theta }}_g$$\end{document} , for $g = 1, \dots, G$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g=1, \ldots , G$$\end{document} . Each group parameter vector is collected in the overall m-dimensional vector $θ = {(θ_{1}^{T}, \dots, θ_{g}^{T}, \dots, θ_{G}^{T})}^{T}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\theta }} = ({\varvec{\theta }}_1^T, \ldots , {\varvec{\theta }}_g^T, \ldots , {\varvec{\theta }}_G^T)^T$$\end{document} , where $m = \sum_{g = 1}^{G} m_{g}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m = \sum _{g=1}^{G} m_g$$\end{document} . Assume for convenience that the same set of parameters is estimated in every group, which implies that the number of observed variables p and factors r is the same across groups, the fixed elements required for identification are placed in the same positions across groups, and that $m_{1} = \dots = m_{G}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_1 = \ldots = m_G$$\end{document} , so that $m = m_{1} G$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m = m_1 G$$\end{document} . Given random samples of sizes $N_{1}, \dots, N_{G}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_1, \ldots , N_G$$\end{document} , with $N = \sum_{g = 1}^{G} N_{g}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N = \sum _{g=1}^{G} N_g$$\end{document} the total sample size across groups, the log-likelihood of the multiple-group factor model is

(14)

\begin{matrix} ℓ (θ) = - \sum_{g = 1}^{G} \frac{N_{g}}{2} {log | Σ_{g} | + tr (W_{g} Σ_{g}^{- 1}) + p log (2 π)}, \end{matrix}

where $W_{g} = S_{g} + ({\bar{x}}_{g} - μ_{g}) {({\bar{x}}_{g} - μ_{g})}^{T}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{W}}_g = {\varvec{S}}_g + (\bar{{\varvec{x}}}_g - {\varvec{\mu }}_g)(\bar{{\varvec{x}}}_g - {\varvec{\mu }}_g)^T$$\end{document} .

In multiple-group analyses, an important methodological consideration is the establishment of the comparability or “equivalence” of measurement across the groups (e.g., countries, socio-economical groups). Measurement (or factorial) invariance occurs when the factors have the same meaning in each group, which translates into equal measurement models (i.e., factor loadings, intercepts and unique variances) across groups (Millsap Reference Millsap2012). If non-equivalence of measurement exists, substantively interesting group comparisons may become distorted. Testing for measurement invariance in the parameters is, however, an intensive process. A sequence of nested tests is progressively conducted to establish the equivalence in the factor loadings, the intercepts, and optionally the unique variances (Vandenberg & Lance, Reference Vandenberg and Lance2000). The next section describes the penalty functions that can be incorporated into the multiple-group model to obtain a technique that automatically detects parameter equivalence across groups.

5. Sparsity and Invariance-Inducing Penalties

As in the single-group factor model, we can penalize the factor loadings to automatically obtain a sparse loading matrix in each of the groups. Define the diagonal matrix $R_{q} = diag (0, \dots, 0, 1, 0,$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{R}}_q = \text {diag}(0, \ldots , 0, 1, 0,$$\end{document} $\dots, 0)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ldots , 0)$$\end{document} , where the 1 on the ${(q, q)}^{th}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(q,q)^{\text {th}}$$\end{document} entry of the matrix corresponds to the $q^{th}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q^{\text {th}}$$\end{document} factor loading in $θ$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\theta }}$$\end{document} , for $q = (g - 1) m_{1} + 1, \dots, (g - 1) m_{1} + q^{⋆}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q = (g-1)m_1 + 1, \ldots , (g-1)m_1 + q^\star $$\end{document} and $g = 1, \dots, G$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g=1, \ldots , G$$\end{document} , and $R_{q} = O_{m \times m}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{R}}_q = {\varvec{O}}_{m \times m}$$\end{document} for the remaining parameters. The quantity $q^{⋆}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q^\star $$\end{document} represents the number of penalized loadings in each group. Then, the sparsity-inducing penalty on the factor loadings is $P_{η_{1}}^{T} (θ) = \sum_{q = 1}^{m} P_{η_{1}, q}^{T} (| | R_{q} {θ | |}_{1})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {P}}^{{\mathcal {T}}}_{\eta _1}({\varvec{\theta }}) = \sum _{q=1}^{m} {\mathcal {P}}^{{\mathcal {T}}}_{\eta _1,q}(||{\varvec{R}}_q {\varvec{\theta }}||_1)$$\end{document} , where $η_{1} \in [0, \infty)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _1 \in [0, \infty )$$\end{document} controls the overall amount of shrinkage.

In the same spirit as factorial invariance, we can specify a penalty encouraging the equality of the loadings across groups. Conveniently, this can be achieved by shrinking the pairwise absolute differences of every factor loading across groups. Let $D_{q}^{Λ}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{D}}_q^{{\varvec{\Lambda }}}$$\end{document} , for $q = 1, \dots, q^{⋆}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q = 1, \ldots , q^\star $$\end{document} , be the matrix computing the differences of the factor loading pairs $(θ_{(g - 1) m_{1} + q}, θ_{(g^{'} - 1) m_{1} + q})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\theta _{(g-1)m_1 +q}, \theta _{(g'-1)m_1+q})$$\end{document} for $g < g^{'}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g < g'$$\end{document} , whereas for the other parameters $D_{q}^{Λ} = O_{m_{1} (\begin{matrix} G \\ 2 \end{matrix}) \times m}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{D}}_q^{{\varvec{\Lambda }}} = {\varvec{O}}_{m_1 \left( {\begin{array}{c}G\\ 2\end{array}}\right) \times m}$$\end{document} . Then, the penalty inducing equal loadings across groups can be written as $P_{η_{2}}^{T} (θ) = \sum_{q = 1}^{m} P_{η_{2}, q}^{T} (| | D_{q}^{Λ} {θ | |}_{1})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {P}}^{{\mathcal {T}}}_{\eta _2} ({\varvec{\theta }}) = \sum _{q=1}^{m} {\mathcal {P}}^{{\mathcal {T}}}_{\eta _2,q}(||{\varvec{D}}_q^{{\varvec{\Lambda }}} {\varvec{\theta }}||_1)$$\end{document} , where $| | D_{q}^{Λ} {θ | |}_{1} = \sum_{g < g^{'}} | θ_{(g - 1) m_{1} + q} - θ_{(g^{'} - 1) m_{1} + q} |$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$||{\varvec{D}}_q^{{\varvec{\Lambda }}} {\varvec{\theta }}||_1 = \sum _{g < g'} |\theta _{(g-1)m_1 + q} - \theta _{(g'-1)m_1 + q} |$$\end{document} for $q = 1, \dots, q^{⋆}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q = 1, \ldots , q^\star $$\end{document} , and zero otherwise. If $G = 2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G = 2$$\end{document} , the absolute difference of the $q^{th}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q^{\text {th}}$$\end{document} loading across the two groups is expressed as $| | D_{q}^{Λ} {θ | |}_{1} = | θ_{q} - θ_{m_{1} + q} |$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$||{\varvec{D}}_q^{{\varvec{\Lambda }}} {\varvec{\theta }}||_1 = |\theta _q - \theta _{m_1 + q} |$$\end{document} , where $D_{q}^{Λ} = [R_{q}, - R_{q}]$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{D}}_q^{{\varvec{\Lambda }}} = [{\varvec{R}}_q, \, \, -{\varvec{R}}_q]$$\end{document} . The expression of $P_{η_{2}}^{T} (θ)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {P}}^{{\mathcal {T}}}_{\eta _2} ({\varvec{\theta }})$$\end{document} for lasso, alasso, scad and mcp is given in Section B.2.1. The tuning parameter $η_{2} \in [0, \infty)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _2 \in [0, \infty )$$\end{document} controls the amount of loading equality across groups. When the loadings are truly invariant and $η_{2}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _2$$\end{document} is properly chosen, the penalized group loading matrices “fuse”, and share the same values.

Lastly, we can encourage the equality of the intercepts across groups by specifying a penalty shrinking their pairwise absolute group differences. Let $k^{⋆}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k^\star $$\end{document} be the number of estimated intercepts in each group. Due to the presence of fixed elements in $τ_{g}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\tau }}_g$$\end{document} for identification, $k^{⋆}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k^\star $$\end{document} is smaller than p. Let $D_{q}^{τ}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{D}}_q^{{\varvec{\tau }}}$$\end{document} , for $q = (g - 1) m_{1} + q^{⋆} + 1, \dots, (g - 1) m_{1} + q^{⋆} + k^{⋆}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q=(g-1)m_1 + q^\star + 1, \ldots , (g-1)m_1 + q^\star + k^\star $$\end{document} , be a matrix of known constants computing the differences of the intercepts across groups, whereas for all of the other parameters (i.e., the loadings, the unique variances and the structural parameters) $D_{q}^{τ} = O_{m_{1} (\begin{matrix} G \\ 2 \end{matrix}) \times m}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{D}}_q^{{\varvec{\tau }}} = {\varvec{O}}_{m_1 \left( {\begin{array}{c}G\\ 2\end{array}}\right) \times m}$$\end{document} . The penalty inducing equal intercepts across groups is then written as $P_{η_{3}}^{T} (θ) = \sum_{q = 1}^{m} P_{η_{3}, q}^{T} (| | D_{q}^{τ} {θ | |}_{1})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {P}}^{{\mathcal {T}}}_{\eta _3}({\varvec{\theta }}) = \sum _{q=1}^{m} {\mathcal {P}}^{{\mathcal {T}}}_{\eta _3,q}(||{\varvec{D}}_q^{{\varvec{\tau }}} {\varvec{\theta }}||_1)$$\end{document} , where $η_{3} \in [0, \infty)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _3 \in [0, \infty )$$\end{document} governs the amount of intercept invariance.

Optionally, one can encourage the invariance of the unique variances. However, as argued by Little et al. (Reference Little, Card, Slegers, Ledford, Little, Bovaird and Card2012), these quantities contain both random sources of errors, for which there is no theoretical reason to expect equality across groups, and item-specific components, which can vary as a function of various measurement factors. In light of this, we do not introduce a penalty on the unique variances, as their cross-group equivalence would not provide any additional evidence of comparability of the constructs because the important measurement parameters (i.e., the factor loadings and the intercepts) are already encouraged to be invariant by penalties $P_{η_{2}}^{T}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {P}}^{{\mathcal {T}}}_{\eta _2}$$\end{document} and $P_{η_{3}}^{T}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {P}}^{{\mathcal {T}}}_{\eta _3}$$\end{document} .

The three penalties can be easily combined into a single penalty that simultaneously generates sparsity on the factor loading matrices and equivalent loadings and intercepts

(15)

\begin{matrix} P_{η}^{T} (θ) = & P_{η_{1}}^{T} (θ) + P_{η_{2}}^{T} (θ) + P_{η_{3}}^{T} (θ) \\ = & \sum_{q = 1}^{m} (P_{η_{1}, q}^{T}, (| |, R_{q}, {θ | |}_{1}, ) +, P_{η_{2}, q}^{T}, (| |, D_{q}^{Λ}, {θ | |}_{1}, ) +, P_{η_{3}, q}^{T}, (| |, D_{q}^{τ}, {θ | |}_{1}, )), \end{matrix}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\mathcal {P}}^{{\mathcal {T}}}_{{\varvec{\eta }}}({\varvec{\theta }})= & {} {\mathcal {P}}^{{\mathcal {T}}}_{\eta _1}({\varvec{\theta }}) + {\mathcal {P}}^{{\mathcal {T}}}_{\eta _2} ({\varvec{\theta }}) + {\mathcal {P}}_{\eta _3}^{{\mathcal {T}}} ({\varvec{\theta }}) \nonumber \\= & {} \sum _{q=1}^{m} \left\{ {\mathcal {P}}^{{\mathcal {T}}}_{\eta _1,q}(||{\varvec{R}}_q {\varvec{\theta }}||_1) + {\mathcal {P}}^{{\mathcal {T}}}_{\eta _2,q}(||{\varvec{D}}_q^{{\varvec{\Lambda }}} {\varvec{\theta }}||_1) + {\mathcal {P}}^{{\mathcal {T}}}_{\eta _3,q}(||{\varvec{D}}_q^{{\varvec{\tau }}} {\varvec{\theta }}||_1) \right\} , \end{aligned}$$\end{document}

where $η = {(η_{1}, η_{2}, η_{3})}^{T}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\eta }} = (\eta _1, \eta _2, \eta _3)^T$$\end{document} is the vector of the tuning parameters. Each penalty is controlled by its own tuning parameter, as we do not a priori expect these values to be equal. The penalties in (15) can be any of the functions illustrated in Sect. 3, including lasso, alasso, scad and mcp, and different penalty functions can be in principle combined. Suppose that the adaptive weights are available for the intercepts but not for the full loading matrices, possibly due to some inadmissible loading values. In this case, one can combine the alasso penalty for intercept similarity with the lasso (which also supports the automatic procedure, contrarily to the scad and mcp) for sparsity and loading equivalence. By following the rationale described in Sect. 3.1, we replace each non-differentiable penalty in (15) with its differentiable local approximation:

which leads to the following differentiable form of the combined penalty:

(16)

Additional details on the structure of the matrix

are given in Section B.2.1. For an example clarifying the formulation of these matrices for the multiple-group model, the reader is referred to Section B.2.2.

6. Generalized Information Criterion

The previously illustrated penalties can be directly introduced within the estimation process by means of penalized maximum likelihood estimation procedures. The penalized log-likelihood is

(17)

\begin{matrix} ℓ_{p} (θ) : = \sum_{α = 1}^{N} {ℓ (x_{α} | θ) - P_{η}^{T} (θ)} = ℓ (θ) - N P_{η}^{T} (θ) . \end{matrix}

For the normal linear factor model, $ℓ (θ)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell ({\varvec{\theta }})$$\end{document} is given in equation (2), $P_{η}^{T} (θ)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {P}}^{{\mathcal {T}}}_{{\varvec{\eta }}} ({\varvec{\theta }})$$\end{document} is one of the penalties of Sect. 3 generating a sparse factor solution, and the vector $η$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\eta }}$$\end{document} reduces to the scalar $η$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} ; for the multiple-group model, $ℓ (θ)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell ({\varvec{\theta }})$$\end{document} is given in equation (14), $P_{η}^{T} (θ)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {P}}^{{\mathcal {T}}}_{{\varvec{\eta }}} ({\varvec{\theta }})$$\end{document} is one of the penalties of Sect. 5 inducing sparsity and invariant loadings and intercepts, and $η$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\eta }}$$\end{document} is equal to the triplet ${(η_{1}, η_{2}, η_{3})}^{T}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\eta _1, \eta _2, \eta _3)^T$$\end{document} .

Simultaneous estimation of all parameters is achieved by maximizing the penalized log-likelihood in (17) and using a local approximation of $P_{η}^{T} (θ)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {P}}^{{\mathcal {T}}}_{{\varvec{\eta }}} ({\varvec{\theta }})$$\end{document} , that is,

(18)

where the function in brackets is now twice-continuously differentiable. The penalized maximum likelihood estimator (PMLE) is then defined as $\hat{θ} = arg {max}_{θ} ℓ_{p} (θ)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{{\varvec{\theta }}} = \arg \max _{{\varvec{\theta }}} \ell _p ({\varvec{\theta }})$$\end{document} . Conveniently, the gradient of the penalized log-likelihood can be written as , where $g (θ) = \frac{\partial ℓ (θ)}{\partial θ}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{g}}({\varvec{\theta }}) = \dfrac{\partial \ell ({\varvec{\theta }})}{\partial {\varvec{\theta }}}$$\end{document} , the Hessian matrix of the second-order derivatives , where , and the expected Fisher information , where .

A crucial aspect lies in the selection of $η$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\eta }}$$\end{document} , which controls the amount of penalization introduced in the model. To select $η$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\eta }}$$\end{document} , we elect to use the Generalized Information Criterion (GIC; Konishi & Kitagawa, Reference Konishi and Kitagawa1996), which is based on a theoretically founded definition of degrees of freedom. Notice that this choice is possible because the quantities we are dealing with are twice-continuously differentiable. We resort to the general approach of the GIC because the penalized maximum likelihood estimator cannot be ascribed to the ordinary maximum likelihood framework postulated by the AIC, and not for relaxing the assumption $E (- \frac{\partial^{2} ℓ (θ)}{\partial θ \partial θ^{T}}) = E (\frac{\partial ℓ (θ)}{\partial θ}, \frac{\partial ℓ (θ)}{\partial θ^{T}})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbb {E}}\left[ -\dfrac{\partial ^2 \ell ({\varvec{\theta }})}{\partial {\varvec{\theta }} \partial {\varvec{\theta }}^T}\right] = {\mathbb {E}}\left[ \dfrac{\partial \ell ({\varvec{\theta }})}{\partial {\varvec{\theta }}} \dfrac{\partial \ell ({\varvec{\theta }})}{\partial {\varvec{\theta }}^T}\right] $$\end{document} , which does hold true for the normal linear factor models considered in this paper. Let G be the true distribution function that generated the data $x_{N} = {x_{1}, \dots, x_{N}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pmb {{\mathscr {x}}}_N = \{{\varvec{x}}_1, \ldots , {\varvec{x}}_N\}$$\end{document} , which are realizations of the random vector

. Let us express the parameter vector as $θ = T (G)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\theta }} = {\varvec{T}}(G)$$\end{document} , where $T (G)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{T}}(G)$$\end{document} is the m-dimensional functional vector of G defined as the solution of the implicit equations $\int ψ (x, T (G)) d G (x) = 0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\int {\varvec{\psi }}({\varvec{x}}, {\varvec{T}}(G)) d G({\varvec{x}}) = {\varvec{0}}$$\end{document} , with

. The log-likelihood and the penalty matrix take different forms depending on whether we deal with a single- or multiple-group factor model. The GIC evaluating the goodness of fit of the penalized model, when used to predict independent future data $z$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{z}}$$\end{document} generated from the unknown distribution G, is (see Online Resource C)

where

and $η$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\eta }}$$\end{document} enters the penalty matrix

. By considering the PMLE $\hat{θ}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{{\varvec{\theta }}}$$\end{document} for $θ$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\theta }}$$\end{document} , and replacing the unknown distribution G with its empirical counterpart $\hat{G}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{G}}$$\end{document} based on the data, we have

where

The effective number or estimated degrees of freedom (edf) of the model is thus equal to

. The formula for the edf is thus readily obtained by adapting the existing results for general likelihoods (of which the differentiable function in (18) is an example) to the penalized framework and assuming the usual regularity conditions. The GIC is an extension of the Akaike Information Criterion (AIC; Akaike, Reference Akaike1974), and as such, it may inherit the tendency of the latter to select overly complex models. To avoid this issue, we can change the constant 2 of the bias term to $log (N)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\log (N)$$\end{document} (used in the Bayesian Information Criterion; Schwarz, Reference Schwarz1978). Then, given grid(s) of values, the optimal $\hat{η}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{{\varvec{\eta }}}$$\end{document} can be chosen using the following Generalized Bayesian Information Criterion (GBIC)

(19)

The optimal penalized factor model is hence chosen to be the one with the lowest BIC, as this is the information criterion routinely employed in sparse settings. However, if researchers are more interested in accuracy and achieving minimum prediction error, then the AIC is to be preferred. In the presence of moderate sample size and many variables, the extended BIC (EBIC; Chen & Chen, Reference Chen and Chen2008) may be more suitable.

The edf of an unpenalized model ( ) coincide with the dimension of the parameter vector $θ$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }$$\end{document} , since , where $I_{m}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{I}}_m$$\end{document} is the $m \times m$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m \times m$$\end{document} identity matrix. For a penalized model . This shows that $edf \to m$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {edf} \rightarrow m$$\end{document} as $η \to 0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\eta }} \rightarrow {\varvec{0}}$$\end{document} , and $edf \to m - r^{⋆}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {edf} \rightarrow m - r^\star $$\end{document} as $η \to \infty$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\eta }} \rightarrow {\varvec{\infty }}$$\end{document} , where $r^{⋆}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r^\star $$\end{document} is the number of penalized elements (equal to $q^{⋆}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q^\star $$\end{document} for the factor model, and $G (q^{⋆} + k^{⋆})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G(q^\star + k^\star )$$\end{document} for the multiple-group extension). When $0 < η < \infty$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{0}}< {\varvec{\eta }} < {\varvec{\infty }}$$\end{document} , the $edf \in [m - r^{⋆}, m]$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {edf} \in [m - r^\star , m]$$\end{document} . The overall edf of a fitted model is given by the sum of the edf for each parameter; each single edf takes a value in the range [0, 1] and quantifies precisely the extent to which each coefficient is penalized.

The existing penalized factor models (Choi et al., Reference Choi, Oehlert and Zou2010; Hirose & Yamamoto, Reference Hirose and Yamamoto2014a; Jacobucci et al., Reference Jacobucci, Grimm and McArdle2016; Huang et al., Reference Huang, Chen and Weng2017; Huang, Reference Huang2018; Jin et al., Reference Jin, Moustaki and Yang-Wallentin2018) compute the degrees of freedom as the number of nonzero parameters (referred in the following as dof), by advocating the fact that the number of nonzero coefficients in a lasso-penalized linear model gives an unbiased estimate of the total degrees of freedom (Zou et al., Reference Zou, Hastie and Tibshirani2007). This way of estimating the degrees of freedom implies that each dof can be either 0 if its parameter has been shrunken to zero, or 1 otherwise. On the contrary, the edf can take any value in [0, 1]. This suggests that, while the definitions of dof and edf may produce equivalent results (for penalties enjoying the oracle property, as the alasso, scad and mcp), in practical situations using edf is expected to yield better-calibrated degrees of freedom. The proposed method also treats the estimated edf as they are. Importantly, the definition of edf directly stems from the estimated bias term of the GIC, which gives it a theoretically founded basis.

7. Penalized Maximum Likelihood Estimation

For any given set of values of $η$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\eta }}$$\end{document} in the penalty matrix, which is hence denoted in the following as

, we minimize $- ℓ_{p} (θ)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$- \ell _p({\varvec{\theta }})$$\end{document} via a trust-region algorithm (Nocedal & Wright, Reference Nocedal and Wright2006). At iteration $t$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathfrak {t}}$$\end{document} , a “model function” $Q_{p}^{[t]}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {Q}}_p^{[{\mathfrak {t}}]}$$\end{document} is constructed, whose behavior near the current point $θ^{[t]}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\theta }}^{[{\mathfrak {t}}]}$$\end{document} is similar to that of the actual objective function. The model function is usually a quadratic approximation of $- ℓ_{p}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-\ell _p$$\end{document} at $θ^{[t]}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\theta }}^{[{\mathfrak {t}}]}$$\end{document} :

where $s$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{s}}$$\end{document} is the trial step vector aiming at reducing the model function, $g_{p} (θ^{[t]})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{g}}_p({\varvec{\theta }}^{[{\mathfrak {t}}]})$$\end{document} the penalized score function, and

the penalized Hessian matrix. For the normal linear factor model, the derivation of the second-order derivatives is a tedious and lengthy process; however, the availability of these quantities guarantees a better accuracy of the algorithm since no numerical approximation is employed. Because the Hessian requires computing many elements, we also considered the Fisher information matrix. If the elements of $(\hat{Σ} - S)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\hat{{\varvec{\Sigma }}} - {\varvec{S}})$$\end{document} are small and the second derivatives not too large, which is often the case, the information matrix is very close to the true Hessian. For the multiple-group model, due to the presence of the parameters for the mean structure besides those for the covariance structure, we only considered the information matrix as it exhibited similar numerical performances to the Hessian at a reduced computational cost. The analytical expressions of these derivatives for the single- and multiple-group model are given in Geminiani (Reference Geminiani2020, Appendices A, F, respectively).

The search for a minimizer of $Q_{p}^{[t]}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {Q}}_p^{[{\mathfrak {t}}]}$$\end{document} is restricted to some region around $θ^{[t]}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\theta }}^{[{\mathfrak {t}}]}$$\end{document} , which is usually the ball ${| | s | |}_{2} < Δ^{[t]}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$||{\varvec{s}}||_2 < \Delta ^{[{\mathfrak {t}}]}$$\end{document} , where $| | \cdot {| |}_{2}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$||\cdot ||_2$$\end{document} is the Euclidean norm, and $Δ^{[t]} > 0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta ^{[{\mathfrak {t}}]} > 0$$\end{document} the trust-region radius at iteration $t$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathfrak {t}}$$\end{document} . The size of the trust region is critical to the effectiveness of each step: if it is too small, the algorithm may miss the opportunity to take a step that moves it closer to the minimizer of the objective function; if it is too large, the minimizer of the model may be far from the one of the objective function in the region, so it may be necessary to reduce the region size and repeat the process. Each iteration of the trust-region algorithm solves the subproblem:

(20)

\begin{matrix} s^{[t]} = & arg min_{s \in R^{m}} Q_{p}^{[t]} {(s) subject to | | s | |}_{2} \leq Δ^{[t]}, \end{matrix}

(21)

\begin{matrix} θ^{[t + 1]} = & θ^{[t]} + s^{[t]}, \end{matrix}

where the current iteration $θ^{[t]}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\theta }}^{[{\mathfrak {t}}]}$$\end{document} is updated with $s^{[t]}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{s}}^{[{\mathfrak {t}}]}$$\end{document} if this step produces an improvement over the objective function. The size of the region is chosen by measuring the agreement between the model function and the objective function at previous iterations through the ratio:

(22)

\begin{matrix} r^{[t]} = \frac{- (ℓ_{p} (θ^{[t]}) - ℓ_{p} (θ^{[t]} + s^{[t]}))}{Q_{p}^{[t]} (0) - Q_{p}^{[t]} (s^{[t]})} . \end{matrix}

The numerator quantifies the actual reduction, and the denominator the predicted reduction. If $r^{[t]}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r^{[{\mathfrak {t}}]}$$\end{document} is negative, the model is an inadequate representation of the objective function over the current trust region, so the step $s^{[t]}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{s}}^{[{\mathfrak {t}}]}$$\end{document} is rejected, and the new problem is solved with a smaller region. If $r^{[t]}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r^{[{\mathfrak {t}}]}$$\end{document} is close to 1, there is good agreement between $Q_{p}^{[t]}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {Q}}_p^{[{\mathfrak {t}}]}$$\end{document} and $- ℓ_{p}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-\ell _p$$\end{document} over $s^{[t]}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{s}}^{[{\mathfrak {t}}]}$$\end{document} , so the model can accurately predict the behavior of the objective function along that step, and the trust region is enlarged for the next iteration. If $r^{[t]}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r^{[{\mathfrak {t}}]}$$\end{document} is positive, but not close to 1, the trust region is not altered, unless it is close to zero or negative, in which case it is shrunken. Algorithm 1 describes the process. The term $Δ_{max}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta _{\text {max}}$$\end{document} represents an overall bound on the step lengths. Trust-region algorithms never run too far from the current iteration as the points outside the trust region are not considered. For this reason, they were shown to be more stable and faster than line search methods, particularly for functions that are non-concave and/or exhibit regions close to flat (Radice, Marra & Wojtyś, Reference Radice, Marra and Wojtyś2016).

An alternative proposal to using a grid-search combined with GBIC is to estimate $η$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\eta }}$$\end{document} automatically and in a data-driven fashion, a development that has not been so far considered in penalized factor analysis. To this end, we propose adapting to the current context the automatic multiple tuning (a.k.a smoothing) parameter selection of Marra and Radice (Reference Marra and Radice2020, see also references therein), which is based on an approximate AIC.

Assume that, near the solution, the trust-region method behaves like a classic unconstrained Newton-Raphson algorithm (Nocedal & Wright, Reference Nocedal and Wright2006). Suppose also that $θ^{[t + 1]}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\theta }}^{[{\mathfrak {t}}+1]}$$\end{document} is the “true” parameter value, and thus $g_{p} (θ^{[t + 1]}) = 0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{g}}_p({\varvec{\theta }}^{[{\mathfrak {t}}+1]}) = {\varvec{0}}$$\end{document} . By using a first-order Taylor expansion of $g_{p} (θ^{[t + 1]})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{g}}_p({\varvec{\theta }}^{[{\mathfrak {t}}+1]})$$\end{document} at $θ^{[t]}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\theta }}^{[{\mathfrak {t}}]}$$\end{document} it follows that

Solving for $θ^{[t]}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\theta }}^{[{\mathfrak {t}}]}$$\end{document} yields, after some manipulation (see Section D.1),

(23)

where

, $K^{[t]} = μ_{K}^{[t]} + ϑ^{[t]}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{K}}^{[{\mathfrak {t}}]} = {\varvec{\mu }}_{{\varvec{K}}}^{[{\mathfrak {t}}]} + {\varvec{\vartheta }}^{[{\mathfrak {t}}]}$$\end{document} with

and

. The square root of

and its inverse are obtained by eigenvalue decomposition. If they are not positive-definite, they are corrected as described in Section D.2. From standard likelihood theory, we have that $ϑ \sim N (0, I_{m})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\vartheta }} \sim {\mathcal {N}}\left( {\varvec{0}}, {\varvec{I}}_m \right) $$\end{document} and $K \sim N (μ_{K}, I_{m})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{K}} \sim {\mathcal {N}}\left( {\varvec{\mu }}_{{\varvec{K}}}, {\varvec{I}}_m \right) $$\end{document} , where

, and $θ_{0}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\theta }}_0$$\end{document} the true parameter vector. Let ${\hat{μ}}_{K}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{{\varvec{\mu }}}_{{\varvec{K}}}$$\end{document} be the predicted value vector for $K$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{K}}$$\end{document} defined as

where

is the influence (or hat) matrix of the fitting problem and depends on the tuning parameters. The quantity

denotes the PMLE. Ideally, the estimation of the tuning parameters should suppress the model complexity unsupported by the data. This can be achieved by minimizing the expected mean squared error of ${\hat{μ}}_{K}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{{\varvec{\mu }}}_{{\varvec{K}}}$$\end{document} from its expectation $μ_{K}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\mu }}_{{\varvec{K}}}$$\end{document} (Section D.3):

(24)

\begin{matrix} E (\frac{1}{N} | | μ_{K} - {\hat{μ}}_{K} {| |}_{2}^{2}) = \frac{1}{N} E (| | K - A_{η}^{T} K {| |}_{2}^{2}) + \frac{2}{N} tr (A_{η}^{T}) - 1 . \end{matrix}

The quantity

can be interpreted as the edf of the penalized model, and is equivalent to the expression of the bias term of the GBIC. The right-hand side of (24) depends on the tuning parameters through $A_{η}^{T}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{A}}_{{\varvec{\eta }}}^{{\mathcal {T}}}$$\end{document} , whereas $K$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{K}}$$\end{document} is linked to the unpenalized part of the model. The tuning parameters are estimated by minimizing an estimate of (24):

(25)

\begin{matrix} V (η) = \frac{1}{N} | | \hat{μ_{K} - {\hat{μ}}_{K}} {| |}_{2}^{2} = \frac{1}{N} | | K - A_{η}^{T} K {| |}_{2}^{2} + \frac{2}{N} tr (A_{η}^{T}) - 1 . \end{matrix}

This is equivalent to the Un-Biased Risk Estimator (UBRE; Wood, Reference Wood2017, Ch. 6) and an approximate AIC (Section D.4), which means that $η$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\eta }}$$\end{document} is estimated by minimizing what is effectively the AIC with number of parameters given by $tr (A_{η}^{T})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {tr}({\varvec{A}}_{{\varvec{\eta }}}^{{\mathcal {T}}})$$\end{document} . In practice, given $θ^{[t + 1]}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\theta }}^{[{\mathfrak {t}}+1]}$$\end{document} , the estimation problem is expressed as

\begin{matrix} η^{[t + 1]} = arg min_{η} V^{[t + 1]} (η) = arg min_{η} (\frac{1}{N} | | K^{[t + 1]} - A_{η}^{T^{[t + 1]}} K^{[t + 1]} {| |}_{2}^{2} + \frac{2}{N} tr (A_{η}^{T^{[t + 1]}}) - 1), \end{matrix}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\varvec{\eta }}^{[{\mathfrak {t}}+1]} = \arg \min _{{\varvec{\eta }}} {\mathcal {V}}^{[{\mathfrak {t}}+1]}({\varvec{\eta }}) = \arg \min _{{\varvec{\eta }}} \left\{ \frac{1}{N} ||{\varvec{K}}^{[{\mathfrak {t}}+1]} - {\varvec{A}}^{{{\mathcal {T}}}^{[{\mathfrak {t}}+1]}}_{{\varvec{\eta }}} {\varvec{K}}^{[{\mathfrak {t}}+1]}||^2_2 + \frac{2}{N} \text {tr} ({\varvec{A}}^{{{\mathcal {T}}}^{[{\mathfrak {t}}+1]}}_{{\varvec{\eta }}}) - 1 \right\} , \end{aligned}$$\end{document}

and solved by adapting the approach by Wood (Reference Wood2004) to the current context. This approach is based on Newton’s method and can evaluate in a stable and efficient way the components in $V (η)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {V}}({\varvec{\eta }})$$\end{document} and their derivatives with respect to $log (η)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\log ({\varvec{\eta }})$$\end{document} (since the tuning parameters can only take positive values). The two steps, one for the estimation of $θ$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\theta }}$$\end{document} and the other for $η$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\eta }}$$\end{document} , are iterated until the algorithm satisfies the stopping criterion $\frac{| ℓ (θ^{[t + 1]}) - ℓ (θ^{[t]}) |}{0.1 + | ℓ (θ^{[t + 1]}) |} < 10^{- 7}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dfrac{|\ell ({\varvec{\theta }}^{[{\mathfrak {t}}+1]}) - \ell ({\varvec{\theta }}^{[{\mathfrak {t}}]}) |}{0.1 + |\ell ({\varvec{\theta }}^{[{\mathfrak {t}}+1]}) |} < 10^{-7}$$\end{document} .

Sometimes the final model could be overly dense and sparser solutions may be desired. One way to achieve this systematically is to increase the amount that each model edf counts, in the UBRE score, by a factor $γ \geq 1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma \ge 1$$\end{document} , called “influence factor” (Wood, Reference Wood2017). The slightly modified tuning criterion then is

(26)

\begin{matrix} V (η) = \frac{1}{N} | | K - A_{η}^{T} K {| |}_{2}^{2} + \frac{2}{N} γ tr (A_{η}^{T}) - 1 . \end{matrix}

For smoothing spline regression models, Kim and Gu (Reference Kim and Gu2004) found that $γ = 1.4$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = 1.4$$\end{document} can correct the tendency to over-fitting of prediction error criteria. However, this work deals with different models, and our focus is not only on fit but also on the recovery of sparse structures, thus higher values may be more appropriate.

The automatic procedure described above is general and easy to implement, but it may occasionally suffer at small sample sizes since it finds its justification asymptotically when the dependence of the Hessian on the tuning parameter(s) vanishes. As argued by Wood (Reference Wood2017), at small sample sizes, it would in principle be more reliable to select the tuning parameter(s) based on a non-approximate function, such as the GBIC and grid-searches, although implementing such an approach in the multiple-group context would introduce further complications and possibly new computational problems and instabilities. Notice also that the automatic procedure relies on the separability of the penalty matrix from the tuning parameter(s). This requirement is satisfied by the lasso and alasso (thus, $T = {L, A}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {T}} = \{L, A\}$$\end{document} ), but not by the scad and mcp which are therefore confined to the grid-search approach. However, this is not problematic because in our simulation experiments and empirical application the alasso generally represented the most convenient choice of penalty based on a number of criteria.

At convergence, the covariance matrix of $\hat{θ}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{{\varvec{\theta }}}$$\end{document} is . However, instead of $V_{\hat{θ}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{V}}_{\hat{{\varvec{\theta }}}}$$\end{document} , for practical purposes, it is more convenient to employ at convergence the alternative Bayesian result (Marra & Wood, Reference Marra and Wood2012). The goodness of fit of the penalized model can then be evaluated through confidence intervals, which are available for each model parameter, obtained from the posterior distribution $θ | {x_{1}, \dots, x_{N}}, η \sim N (\hat{θ}, V_{θ})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\theta }} |\{{\varvec{x}}_1, \ldots , {\varvec{x}}_N\}, {\varvec{\eta }} \sim {\mathcal {N}}(\hat{{\varvec{\theta }}}, {\varvec{V}}_{{\varvec{\theta }}})$$\end{document} (Section D.5). Notice that the proposed approach can be regarded as a Bayesian method with the exponential prior on the penalty function. The process of determining the optimal loading pattern can indeed be formulated as a Bayesian variable selection problem (Lu, Chow & Loken, Reference Lu, Chow and Loken2016). For instance, Bayesian Structural Equation Modeling (BSEM; Muthén & Asparouhov, Reference Muthén and Asparouhov2012)—in which the elements that would be fixed to zero in a confirmatory analysis (usually the cross-loadings) are replaced with approximate zeros based on informative, small-variance priors—is a particular case where the shrinkage is achieved through an informative ridge prior. With the proposed method, users can rely on the automatic procedure for recovering optimally sparse factor solutions without manually specifying the variance of the Bayesian prior employed in BSEM.

The presented modeling framework has been implemented in the freely available R package penfa and we refer the reader to Online Resource F for a brief description of the software and practical illustrations.

8. Simulation Studies

The performances of the proposed PMLE were evaluated and compared to the penalized methods by Jacobucci et al. (Reference Jacobucci, Grimm and McArdle2016, R package regsem) and Huang (Reference Huang2018, R package lslx) in two extensive simulation studies, one for the normal linear factor model and the other for its multiple-group extension. Despite the presence of other penalized factor analysis techniques (Choi et al., Reference Choi, Oehlert and Zou2010; Hirose & Yamamoto Reference Hirose and Yamamoto2014b, Reference Hirose and Yamamotoa; Trendafilov et al., Reference Trendafilov, Fontanella and Adachi2017; Jin et al., Reference Jin, Moustaki and Yang-Wallentin2018), our choice fell on regsem and lslx because they allow the specification of fixed, free and penalized parameters, as well the estimation of the structural model.

8.1. Simulation Study I

The first simulation evaluates the performances of the proposed technique in a single-group factor analysis model. We evaluate the impact of several conditions, including the sample size, the penalty function, the type of second-order derivative information used in the trust-region algorithm, the strategy for the choice of the tuning parameter, the magnitude of the influence factor and—for some of the penalties—the value of the additional tuning parameter. The simulation was partly inspired by the empirical application (Sect. 9), therefore the number of variables ( $p = 9$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p = 9$$\end{document} ) and of factors ( $r = 3$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r = 3$$\end{document} ) exactly match those of the real data analysis. The conditions that were varied are:

• Sample size: 300, 500, and 1000 observations. These values are in line with those investigated in similar simulation studies (Huang et al., Reference Huang, Chen and Weng2017; Jacobucci et al., Reference Jacobucci, Grimm and McArdle2016; Jin et al., Reference Jin, Moustaki and Yang-Wallentin2018; Hirose & Yamamoto, Reference Hirose and Yamamoto2014b) and include two moderate sample sizes (which are commonly found in psychometric applications) and a large one (to mimic asymptotic behavior). Note that 300 is close to the number of observations in the empirical example;
• Penalty function: lasso, alasso, scad, and mcp were examined in their ability to shrink to zero small loadings without possibly affecting the remaining ones;
• Information matrix: either the Hessian or the Fisher information matrix was used in the optimization process (see Sect. 7);
• Shrinkage parameter selection: this was achieved either by a grid-search or through the automatic procedure. The grid-search was conducted over 200 distinct values of $η$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} and for all four penalty types, with the optimal model being the one with the lowest GBIC. The elements of the grid were adapted based on the specific combination of penalty type and sample size. The automatic procedure was used with lasso and alasso;
• Influence factor: informed by the values that performed well in the application, we investigated different values for the influence factor, namely, $γ = {1, 1.4, 2, 2.5, 3, 3.5, 4, 4.5}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = \{1, 1.4, 2, 2.5, 3, 3.5, 4, 4.5 \}$$\end{document} ;
• Additional tuning parameter: we tested different values of the additional tuning parameter of the alasso, scad and mcp. For the alasso $a = {1, 2}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = \{1, 2\}$$\end{document} , for the scad $a = {2.5, 3, 3.7, 4.5}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = \{2.5, 3, 3.7, 4.5\}$$\end{document} (with 3.7 being the conventional level employed in the literature and suggested by Fan & Li, Reference Fan and Li2001), and for the mcp $a = {2.5, 3, 3.5}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = \{2.5, 3, 3.5\}$$\end{document} .

The population parameters complied to the following structure:

\begin{matrix} Λ^{T} = & (\begin{matrix} 0.85 & 0.75 & 0.65 & \underset{̲}{0} & 0 & 0 & \underset{̲}{0} & 0 & 0.30 \\ \underset{̲}{0} & 0 & 0.30 & 0.85 & 0.75 & 0.65 & \underset{̲}{0} & 0 & 0 \\ \underset{̲}{0} & 0 & 0 & \underset{̲}{0} & 0 & 0.30 & 0.85 & 0.75 & 0.65 \end{matrix}) \\ Φ = & (\begin{matrix} \underset{̲}{1} & 0.3 & 0.3 \\ \underset{̲}{1} & 0.3 \\ \underset{̲}{1} \end{matrix}) \end{matrix}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\varvec{\Lambda }}^T= & {} \left[ \begin{array}{ccccccccc} 0.85 &{} 0.75 &{} 0.65 &{} {{{\underline{0}}}} &{} 0 &{} 0 &{} {{{\underline{0}}}} &{} 0 &{} 0.30 \\ {{{\underline{0}}}} &{} 0 &{} 0.30 &{} 0.85 &{} 0.75 &{} 0.65 &{} {{{\underline{0}}}} &{} 0 &{} 0 \\ {{{\underline{0}}}} &{} 0 &{} 0 &{} {{{\underline{0}}}} &{} 0 &{} 0.30 &{} 0.85 &{} 0.75 &{} 0.65 \\ \end{array} \right] \\ {\varvec{\Phi }}= & {} \left[ \begin{array}{ccc} {{{\underline{1}}}} &{} 0.3 &{} 0.3 \\ &{} {{{\underline{1}}}} &{} 0.3 \\ &{} &{} {{{\underline{1}}}} \\ \end{array} \right] \end{aligned}$$\end{document}

with $Ψ = I_{p} - Λ Φ Λ^{T}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\Psi }} = {\varvec{I}}_p - {\varvec{\Lambda }} {\varvec{\Phi }} {\varvec{\Lambda }}^T$$\end{document} , where $I_{p}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{I}}_p$$\end{document} is the $p \times p$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p \times p$$\end{document} identity matrix. Elements in italic and underlined were fixed for scale setting and identification purposes. The specific values of the factor loadings were inspired by the numerical example in Huang et al. (Reference Huang, Chen and Weng2017). As it is common in many factor analysis applications, a subset of the observed variables does not load only on one factor but also presents a cross-loading.

All of the factor loadings were penalized for assessing the effectiveness of the proposed method in recovering the underlying factor structure and not erroneously shrinking the small cross-loadings to zero. Based on results from previous studies (see for instance Choi et al., Reference Choi, Oehlert and Zou2010 for the alasso, and Hirose & Yamamoto, Reference Hirose and Yamamoto2014b and Huang et al., Reference Huang, Chen and Weng2017 for the mcp), the alasso and the non-convex penalties are expected to outperform the lasso, which is known to be biased due to its tendency to overly shrink nonzero parameters. Concerning the influence factor, higher values favor sparsity at the expense of an increase in bias, whereas lower values favor goodness of fit.

Data were simulated in R (R Core Team, 2018) according to the population parameters. The resulting data matrix was then analyzed in penfa, regsem (Jacobucci et al., Reference Jacobucci, Grimm, Brandmaier, Serang, Kievit and Scharf2019) and lslx (Huang & Hu, Reference Huang and Hu2019) by estimating a factor model with the correct number of factors, the specified fixed elements, and all of the free loadings were penalized. Common factors were estimated to be correlated and with fixed unit variance. Whenever present, sign reversal of the factors was accounted for to ensure that the sign of the primary loadings corresponded to that of the corresponding population parameters. Based on the availability of the respective software implementations,Footnote 2 lasso, alasso, scad and mcp were tried for regsem, and lasso and mcp for lslx. For each scenario, we generated $L = 1000$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L =1000$$\end{document} replications for which the unpenalized factor model produced admissibleFootnote 3 solutions.

We evaluated the performance of the methods according to the criteria illustrated in Huang et al. (Reference Huang, Chen and Weng2017), which are briefly mentioned here. The overall estimation quality was assessed using the estimated mean squared error (MSE):

(27)

\begin{matrix} \hat{MSE} (\hat{θ}) = \frac{1}{L} \sum_{l = 1}^{L} {({\hat{θ}}^{(l)} - θ_{0})}^{T} ({\hat{θ}}^{(l)} - θ_{0}), \end{matrix}

where ${\hat{θ}}^{(l)} = {({\hat{θ}}_{1}^{(l)}, \dots, {\hat{θ}}_{m}^{(l)})}^{T}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{{\varvec{\theta }}}^{(l)} = ({\hat{\theta }}^{(l)}_1, \ldots , {\hat{\theta }}^{(l)}_m)^T$$\end{document} denotes the vector of estimated parameters in replicate l, $θ_{0}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\theta }}_0$$\end{document} the true parameter vector, and L the number of replications. The degree of bias of each estimator was evaluated by the estimated squared bias (SB):

(28)

\begin{matrix} \hat{SB} (\hat{θ}) = {(\bar{\hat{θ}} - θ_{0})}^{T} (\bar{\hat{θ}} - θ_{0}), \end{matrix}

where $\bar{\hat{θ}} = \frac{1}{L} \sum_{l = 1}^{L} {\hat{θ}}^{(l)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{\hat{{\varvec{\theta }}}} = \frac{1}{L} \sum _{l=1}^{L} \hat{{\varvec{\theta }}}^{(l)}$$\end{document} represents the empirical mean of $\hat{θ}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{{\varvec{\theta }}}$$\end{document} . Let $F = {q | θ_{0 q} \neq 0 & {\hat{θ}}_{q} penalized}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\mathcal {F}} = \{q \, |\, \theta _{0 q} \ne 0 \, \& \, {\hat{\theta }}_{q} \text { penalized} \}$$\end{document} indicate the set of indices associated to the true nonzero parameters that have been penalized (i.e., the penalized nonzero factor loadings) and $| F |$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|{\mathcal {F}} |$$\end{document} the cardinality of $F$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {F}}$$\end{document} , which in the simulation is equal to 12. The chance of correctly identifying the true nonzero parameters was evaluated via the estimated true positive rate (TPR):

(29)

Denote as $F^{c} = {q | θ_{0 q} = 0 & {\hat{θ}}_{q} penalized}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\mathcal {F}}^c = \{q \, |\, \theta _{0 q} = 0 \, \& \, {\hat{\theta }}_{q} \text { penalized} \}$$\end{document} the set collecting the indices of the true zero parameters that have been penalized (i.e., the penalized zero factor loadings), with $| F^{c} |$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|{\mathcal {F}}^c |$$\end{document} equal to 9. The estimated false positive rate (FPR) examined the degree to which the true zero parameters were incorrectly identified as nonzero:

(30)

Lastly, selection consistency was assessed via the proportion of times the true model—for which all the true zero and nonzero factor loadings were correctly identified as equal to zero and different from zero, respectively—was chosen over the replicates (proportion choosing the true model; PCTM):

(31)

where $| F | + | F^{c} | = q^{⋆}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|{\mathcal {F}} |+ |{\mathcal {F}}^c |= q^\star $$\end{document} . For the computation of PCTM and FPR,Footnote 4 the parameter estimates were rounded to one decimal digitFootnote 5 for all models. For the sake of clarity, we report a selection of the most relevant results for the configurations of penfa-alasso ( $a = 2, γ = 4.5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a=2, \gamma = 4.5$$\end{document} ), penfa-scad ( $a = 3$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a=3$$\end{document} ) and penfa-mcp ( $a = 3$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a=3$$\end{document} ) that produced the best models in terms of the aforementioned performance criteria. Due to its generally higher numerical stability in comparison to the Hessian, only penfa models estimated with the Fisher information matrix are presented. The results for penfa-lasso are described in Online Resource A. In the same spirit, the results of regsem and lslx are presented for their best performing models (i.e., with the mcp for both of them). All other results can be requested from the corresponding author.

Overall, the low values for MSE, the bias, and FPR which are very close to zero, together with high PCTM and excellent TPR show that the examined penalized techniques possess very good empirical performances and outperform the unpenalized (MLE) model (Table 1). The MSE of all penalized methods are very similar to each other and improve as the sample size increased. The results with the lower bias were associated with the use of non-convex penalties, although the bias of penfa-alasso very quickly converged to zero when the sample size increased, and hence the impact of the penalty decreased. The true positive rates were always equal to 1.0, which showed that the inspected methods never suppressed the nonzero penalized parameters (i.e., the primary loadings and the cross-loadings). In terms of both false positive rates and selection consistency, penfa-alasso with automatic tuning parameter selection presented by far the best performances for all the sample sizes. The coverage probabilities for the parameters of all fitted models (Section A.1) were generally close to their true nominal level for all penalty functions.

Table 1. Performance measures of the examined models in simulation study I by varying the sample size N.

MSE mean-squared error, SB squared bias, TPR true positive rate, FPR false positive rate, PCTM proportion choosing the true model. In brackets, the ranges of MSE, TPR, and FPR across replicates.

The mean squared error and bias of penfa-alasso with automatic tuning parameter selection were similar to those obtained with the same penalty and grid-search, but the false positives and PCTM were markedly lower and higher, respectively. This is due to the way the optimal penalized model is picked. With the automatic procedure, the optimal model is the one whose tuning parameter minimizes the criterion in (26); with the grid-search, the optimal model minimizes the GBIC in (19). However refined, a grid-search cannot compete with an approach that looks for the optimal tuning parameter on the positive real line. In addition, the presence of a sparsity-inducing quantity (the influence factor) in the optimization criterion helped the model obtain a nicer tradeoff between goodness of fit and model complexity. With reference to the exponent a in the expression of the alasso, as this quantity increased the weights became more influential, and we observed a general improvement in all the performance measures. The best results were obtained for $a = 2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = 2$$\end{document} .

By comparing the quality measures of the three methods for the same penalty function (i.e., the mcp), we notice that penfa outperformed lslx and was generally close to regsem for MSE and SB and superior for FPR and PCTM. This might be due to several aspects, e.g., the optimization algorithm, the internal software package implementations, the formulation of the degrees of freedom, and possibly the approximation of the penalty.

The examined performance criteria explored different conflicting objectives. Ideally, one desires a model with low bias and little complexity (i.e., a sparse solution), but the two measures cannot be minimized simultaneously. This can be seen by looking at the performances of the penfa-alasso model for extreme values of the influence factor (i.e., $γ = 4.5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = 4.5$$\end{document} in Table 1 and $γ = 1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = 1$$\end{document} in Table A.2 in Section A.2). The higher value of $γ$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} produced sparser solutions (i.e., smaller FPR and larger PCTM), at the cost of a larger bias. As the sample size increased, the discrepancies in the performances of the models with different values of $γ$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} diminished though.

The models fitted through the automatic tuning parameter procedure exhibited markedly shorter computational timesFootnote 6 than grid-search methods. Specifically, the average median elapsed times were 17 sec for penfa-alasso with grid (1-dim. grid for $η$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} ; $a = 2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = 2$$\end{document} ) and 0.3 sec with automatic procedure ( $a = 2; γ = 4.5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = 2; \gamma = 4.5$$\end{document} ), 21.1 sec for penfa-scad (1-dim. grid for $η$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} ; $a = 3$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = 3$$\end{document} ), 20.7 sec for penfa-mcp (1-dim. grid for $η$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} ; $a = 3$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a=3$$\end{document} ), 6.6 sec for lslx-mcp (2-dim. grid for $η$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} and a) and 42.2 sec for regsem-mcp (1-dim grid for $η$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} , $a = 3.7$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = 3.7$$\end{document} as per default software implementations). The penfa models with the automatic procedure exhibited the lowest computational times, which is also merit of the stability of the trust-region optimizer, whose parameter updates only involve the points within a proper trust-region. The computational times of lslx are lower than those of the other grid-search techniques because the underlying optimizer is implemented in C++, which significantly boosts the computations with respect to base R routines.

8.2. Simulation Study II

The second simulation evaluates the ability of the proposed technique in identifying the pattern of partial invariance in a multiple-group factor model as a function of the sample size, the size of the generated difference in the group-specific loadings and intercepts, the magnitude of the influence factor and the value of the additional tuning parameter. Since the current implementation of regsem does not allow for multiple-group analyses, our method is only compared with lslx.

Table 2. The factor loading matrices and intercepts of the two groups under each difference scenario of simulation study II.

Elements fixed for origin and scale setting and identification purposes are italic and underlined. Under the null condition, the parameters of Group 2 coincide with those of Group 1.

We consider a population multiple-group factor model with $p = 12$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p = 12$$\end{document} variables, $r = 3$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r=3$$\end{document} factors and $G = 2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G=2$$\end{document} groups. We explore a range of conditions, under which the factor loading matrices and intercepts are either invariant or non-invariant, with the level of non-invariance becoming progressively larger. Based on the findings from Simulation study I, we employ the alasso penalty for inducing sparsity and invariant loadings and intercepts, that is, . The three tuning parameters ${(η_{1}, η_{2}, η_{3})}^{T}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\eta _1, \eta _2, \eta _3)^T$$\end{document} in $η$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\eta }}$$\end{document} are estimated alongside the model parameters through the automatic multiple tuning parameter procedure. For lslx we used the mcp penalty, which had better performances than the lasso. The optimization technique currently employed in lslx makes use of a single penalty for both shrinking the parameters and their differences across groups. Therefore, there is only one shrinkage parameter $η$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} , whose optimal value is determined through a grid-search. For lslx-mcp, we carried out a grid-search over 200 values of the shrinkage parameter $η$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} and 4 of the shape parameter a. The conditions that were varied are:

• Sample size: 300, 500, and 1000 observations evenly split between the two groups, with 300 being close to the number of observations in the empirical example;
• Difference size: either null, small, medium or large group differences in the primary loadings and the intercepts of two variables were created (details are given below). This condition was partly inspired by the simulation conducted by Huang (Reference Huang2018);
• Influence factor: informed by the values that performed well in Simulation study I, we investigated three values of the influence factor, namely, $γ = {3.5, 4, 4.5}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = \{3.5, 4, 4.5\}$$\end{document} ;
• Additional tuning parameter: two values were tested for the exponent in the expression of the alasso, namely $a = {1, 2}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = \{1,2\}$$\end{document} .

The factor loading matrix and the vector of intercepts of Group 1 are reported on the left-hand side of Table 2, and are the same under every difference scenario. Elements in italic and underlined are fixed for metric setting and identification purposes. The factor loadings and intercepts of Group 2 are presented by difference scenario on the right-hand side of Table 2. In case of a null difference, the two groups share the same parameter matrices. Under the small, medium and large scenarios, the primary loadings and the intercepts of two variables (i.e., $x_{6}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_6$$\end{document} and $x_{12}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_{12}$$\end{document} ) in Group 1 differ from the corresponding parameters in Group 2 by a size of 0.1, 0.2, and 0.3, respectively. Under all conditions, the structural parameters are assumed to be invariant across groups, that is, $vech (Φ_{1}) = vech (Φ_{2}) = vech (Φ) = {(1, 0.3, 1)}^{T}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {vech}({\varvec{\Phi }}_1) = \text {vech}({\varvec{\Phi }}_2) = \text {vech}({\varvec{\Phi }}) = (1, 0.3, 1)^T$$\end{document} and $κ_{1} = κ_{2} = {(0, 0)}^{T}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\kappa }}_1 = {\varvec{\kappa }}_2 = (0, 0)^T$$\end{document} , whereas $Ψ_{g} = I_{p} - Λ_{g} Φ Λ_{g}^{T}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\Psi }}_g = {\varvec{I}}_p - {\varvec{\Lambda }}_g {\varvec{\Phi }} {\varvec{\Lambda }}_g^T$$\end{document} , for $g = 1, 2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g = 1, 2$$\end{document} . The factor loadings and the intercepts are penalized in the way described in Sect. 5 (i.e., shrinkage of the loadings and of the pairwise group differences of loadings and intercepts), whereas the remaining model parameters are estimated without penalization. For each scenario, we generated $L = 1000$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L =1000$$\end{document} replications for which the unpenalized multiple-group model produced admissible solutions, and analyzed them as in simulation study I.

The performances of the penalized models are evaluated through the criteria (27)-(31) used in simulation study I. For the sake of conciseness, we report the results for the penfa-alasso model ( $a = 2, γ = 4.5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = 2, \gamma = 4.5$$\end{document} ) that produced the best solution in terms of these performance criteria. All other results can be requested from the corresponding author. Overall, the low values of MSE, SB, FPR, high PCTM and excellent TPR show that the penalized techniques possess very good empirical performances, with all measures improving as the sample size increased (Table 3). Higher difference sizes were associated with higher MSE and squared bias, with the lower values generally occurring for penfa-alasso. We separately computed these measures for each parameter matrix (that is, $Λ_{g}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\Lambda }}_g$$\end{document} , $τ_{g}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\tau }}_g$$\end{document} , $Ψ_{g}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\Psi }}_g$$\end{document} , $Φ_{g}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\Phi }}_g$$\end{document} , $κ_{g}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\kappa }}_g$$\end{document} , for $g = 1, 2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g = 1, 2$$\end{document} ) produced by penfa-alasso. The largest MSE were observed for the factor variances and covariances, followed by the factor loadings. The bias tended to increase for the penalized parameters (factor loadings and intercepts) across the difference conditions, while remaining almost unaltered for the unique variances and the structural parameters. The squared bias quickly converged towards zero in all difference scenarios as the sample size increased. The TPR were always equal to 1.0, which showed that the examined methods never suppressed the nonzero penalized parameters.

Table 3. Performance measures of penfa-alasso ( $a = 2, γ = 4.5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = 2, \gamma = 4.5$$\end{document} ) and lslx-mcp models by sample size and difference scenario.

MSE mean-squared error, SB squared bias, TPR true positive rate, FPR false positive rate, PCTM proportion choosing the true model. In brackets, the ranges of MSE and FPR across replicates; TPR were always equal to one.

Whereas under the null and small scenarios the two methods produced similar measures, penfa-alasso markedly outperformed lslx-mcp under the medium and large conditions, especially in terms of selection consistency at the smallest sample size. On top of that, whereas these performance measures for lslx noticeably degraded as the difference size increased, they remained fairly stable for penfa-alasso; even with the smallest sample size, penfa-alasso identified the true heterogeneity pattern more than 90% of the times. Thanks to the use of the automatic multiple tuning parameter procedure, the average median computational time to fit a penfa-alasso model with 3 tuning parameters (3.2 seconds) was much lower than the one necessary to fit an lslx-mcp model with a single shrinkage parameter $η$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} and the associated shape parameter a selected through a grid-search (45 seconds).

9. Empirical Application

The Holzinger & Swineford data set (Holzinger & Swineford, Reference Holzinger and Swineford1939; Kelley, Reference Kelley2019) is a classical psychometric application containing the responses of $N = 301$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N = 301$$\end{document} students on some psychological tests. This data set (or subsets of it) has been often used to demonstrate CFA (Jöreskog, Reference Jöreskog, Jöreskog, Sörbom and Magidson1979), EFA (Browne, Reference Browne2001; Jöreskog & Sörbom, Reference Jöreskog and Sörbom1993) and various penalized factor analysis techniques (Trendafilov et al., Reference Trendafilov, Fontanella and Adachi2017; Jacobucci et al., Reference Jacobucci, Grimm and McArdle2016; Huang et al., Reference Huang, Chen and Weng2017; Jin et al., Reference Jin, Moustaki and Yang-Wallentin2018). For space constraints, the description of the data set is reported in Online Resource E.

9.1. Normal Linear Factor Model

Following Jacobucci et al. (Reference Jacobucci, Grimm and McArdle2016) and Huang et al. (Reference Huang, Chen and Weng2017), to illustrate the proposed method in the normal linear factor model, we use a subset of nine mental tests (VISUAL, CUBES, FLAGS, PARAGRAP, SENTENCE, WORDM, ADDITION, COUNTING, and STRAIGHT) underlying three latent factors. The data set was column-wise centered since the model in equation (1) assumes that the observed variables have zero means and scaled as described in Yuan and Bentler (Reference Yuan and Bentler2006). The inspection at the covariance matrix of the observed variables revealed the presence of relationships between tests designed to measure distinct mental abilities. The CFA model assuming a simple structure showed a poor fit to the data (p-value of the chi-square goodness of fit test < 0.001), which suggested the multi-dimensionality of some of the tests. In these circumstances where it may be difficult to specify the correct sparsity pattern of the loading matrix in advance, it is beneficial to resort to penalized techniques to explore and unveil the underlying loading pattern. We hence penalize all of the factor loadings and freely estimate the remaining model parameters. Factor variances are fixed to one for scale setting and some elements of the loading matrix to zero for identification purposes. Even if the proposed method does allow us to obtain sparsity, we should acknowledge that its achievement also depends on the features of the statistical model under investigation and the amount of information carried by the data. Concerning the former, as pointed out by Trendafilov et al. (Reference Trendafilov, Fontanella and Adachi2017), inducing sparsity in a factor model, and even more so one with correlated factors, is more complicated than for other types of models (e.g. principal component analysis) due to the presence of other parameters (unique variances and factor variances and covariances) affecting the overall model fit. As a result, if too large a value for the tuning parameter is chosen, an excessive number of loadings is shrunken, and the remaining parameters are forced to explode to compensate for this lack of fit. This issue can be avoided if the appropriate amount of sparsity is introduced into the model, which in turn is only possible if the tuning parameter governing the amount of sparsity is selected according to a valid procedure, such as the one introduced in the paper.

We fitted a large number of models involving all four penalties. For grid-search, 200 models corresponding to varying levels of the tuning parameter were fitted. We also tried a sequence of values for the additional tuning parameter of the alasso ( $a = {1, 1.5, 2}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a =\{ 1, 1.5, 2\}$$\end{document} ), scad ( $a = {2.5, 3.7, 4.5}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a =\{ 2.5, 3.7, 4.5\}$$\end{document} ), mcp ( $a = {1.5, 2, 2.5, 3, 3.5}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = \{1.5, 2, 2.5, 3, 3.5\}$$\end{document} ), and for the influence factor ( $γ = {1, 1.4, 2, 2.5, 3, 3.5, 4, 4.5}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = \{1, 1.4, 2, 2.5, 3, 3.5, 4, 4.5\}$$\end{document} ). The GBICFootnote 7 values were calculated for each of the fitted penfa models and are ranked in Table 4 for the best model configurations. In particular, the alasso (automatic procedure, $a = 1, γ = 4.5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = 1, \gamma = 4.5$$\end{document} ) presented the lowest BIC, closely followed by the mcp ( $a = 1.5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = 1.5$$\end{document} ) and scad ( $a = 4.5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = 4.5$$\end{document} ). Interestingly, the BIC of penfa-lasso with grid-search (7567.62) markedly decreased when the model was fitted through the automatic procedure with an influence factor of 4.5 (7562.94). Notice that both the CFA and the unpenalized solution (corresponding to the factor analysis model in equation (1) with the minimum identification restrictions) resulted in worse fits than the ones of the penalized models, probably because of the strict assumption of no cross-loadings of the former, and the unnecessary complexity of the latter. This indicates that the analysis benefited from the introduction of sparsity.

Table 4. BIC of the best configurations of the fitted models.

For penfa-alasso (automatic procedure) $a = 1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = 1$$\end{document} and $γ = 4.5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = 4.5$$\end{document} , for penfa-scad $a = 4.5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = 4.5$$\end{document} , for penfa-mcp $a = 1.5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = 1.5$$\end{document} , and for penfa-lasso (automatic procedure) $γ = 4.5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = 4.5$$\end{document} . For all models the Fisher information was used.

Table 5 (left-hand side) reports the parameter estimates of the unpenalized model and the best penfa-alasso model. A blank cell in the factor loading matrix indicates that the corresponding estimate was zero after one decimal digit rounding.Footnote 8 The unpenalized model presented various cross-loadings, which resulted in a complex model. For penfa, only four secondary loadings ( ${\hat{λ}}_{51}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{\lambda }}_{51}$$\end{document} , ${\hat{λ}}_{81}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{\lambda }}_{81}$$\end{document} , ${\hat{λ}}_{91}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{\lambda }}_{91}$$\end{document} , ${\hat{λ}}_{32}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{\lambda }}_{32}$$\end{document} ) were identified as nonzero. If a sparser loading matrix is desired, users can increase the value of the exponent a of the alasso and/or the influence factor $γ$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} in the automatic procedure. For instance, a penfa-alasso model (BIC = 7565.39) with $a = 2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = 2$$\end{document} and $γ = 5.5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = 5.5$$\end{document} (Table 5, right-hand side) produced a sparser factor solution with a single cross-loading ( ${\hat{λ}}_{91}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{\lambda }}_{91}$$\end{document} ). The data analysis was also conducted for regsem and lslx using the available penalties (i.e., lasso, alasso, scad, and mcp for the former, and lasso and mcp for the latter) and is presented in Online Resource E. The factor structures of the penalized models looked similar, but the proposed method reported the lowest BIC values, showing the potential of the presented procedure. As argued by Huang et al. (Reference Huang, Chen and Weng2017), this example shows that complex models do not necessarily outperform simpler ones when model complexity is also taken into account in the model selection criterion.

Table 5. Parameter estimates of the nine mental tests from the Holzinger & Swineford data set for the unpenalized model, and penfa-alasso with automatic procedure (on the left-hand side, $\hat{η} = 0.017, a = 1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{\eta }} = 0.017, a = 1$$\end{document} and $γ = 4.5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = 4.5$$\end{document} ; on the right-hand side, $\hat{η} = 0.011, a = 2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{\eta }} = 0.011, a = 2$$\end{document} and $γ = 5.5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = 5.5$$\end{document} ).

Fixed parameters are italic and underlined. A blank cell indicates that the corresponding estimate is zero.

9.2. Multiple-Group Factor Model

Besides considering the sample of the students as a whole, we divided it into two groups ( $N_{1} = 156, N_{2} = 145$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_1 = 156, N_2 = 145$$\end{document} ) based on the attended school, and then conducted a multiple-group analysis. One school (Pasteur) included students with parents who immigrated from Europe, whereas the other (Grant-White) was composed of students coming from middle-income American white families. Following Huang (Reference Huang2018), we considered the 19 mental tests and standardized the data to handle the scaling effect.

The traditional approach consists of the estimation of an unpenalized multiple-group CFA in which the tests are assumed to be pure measures, followed by factorial invariance testing procedures. The model assuming equal loadings across groups shows an adequate fit to the data (p-value of the chi-square goodness of fit test = 0.266), which, however, significantly worsens when the intercepts are also equated across groups (p-value of the likelihood ratio test comparing the model with invariant loadings and intercepts versus the one with only invariant loadings < 0.001). Model modifications are typically conducted to determine and freely estimate the non-invariant elements.

Alternatively, the invariance pattern can be explored via penalized techniques employing penalties that combine sparsity and cross-group equivalence of loadings and intercepts. In light of its superior performance in the single-group analysis and simulation, we employed the alasso with the automatic multiple tuning parameter procedure, and tested various values of the influence factor $(γ = {1, 2, 3, 3.5, 4, 4.5})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\gamma = \{1, 2, 3, 3.5, 4, 4.5\})$$\end{document} and the exponent $(a = {1, 2})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(a =\{ 1, 2\})$$\end{document} . The tests VISUAL, WORDM, COUNTING and NUMBERR are assumed to be the markers, and thus have fixed loadings and intercepts. The data analysis was also conducted in lslx with the mcp (see Table E.4), but not in regsem as its current implementation does not allow for multiple-group analyses. Note that lslx uses only one penalty for shrinking both the parameters and their differences, hence it has a single tuning parameter $η$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} .

The parameter estimates of penfa-alasso are reported in Table 6. The better fit of penfa-alasso ( $BIC = 14658$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {BIC} = 14658$$\end{document} ) as compared to lslx-mcp ( $BIC = 14697.75$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {BIC} = 14697.75$$\end{document} ) is also merit of the greater flexibility of the former, which employs three distinct penalties having their own tuning parameters, with respect to the latter, where a single tuning has to take care of the shrinkage of the parameters as well as their cross-group differences. penfa-alasso produces sparse loading matrices with many zero-entries, but the presence of a couple of nonzero cross-loadings demonstrates that the structure hypothesized by a multiple-group CFA is too restrictive. The factor loading matrices of penfa-alasso are also fully equivalent, in agreement to the results of invariance testing. Conversely, the intercepts are not fully invariant, which is again in line with the findings from factorial invariance testing. This example clearly shows the benefits of using properly designed penalized techniques to explore the non-equivalence pattern of the parameter matrices in a multiple-group factor model.

Table 6. Parameter estimates of the 19 mental tests from the Holzinger & Swineford data set for penfa-alasso (automatic procedure, $\hat{η} = {(0.006, 16221.852, 0.013)}^{T}, a = 1, γ = 4$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{{\varvec{\eta }}} = (0.006, 16221.852, 0.013)^T, a = 1, \gamma = 4$$\end{document} )

Fixed parameters are italic and underlined. A blank cell in the factor loading matrix indicates that the corresponding estimate is zero. Non-invariant parameters across groups are starred ( $^{⋆}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^\star $$\end{document} ).

10. Discussion

Penalized factor analysis is an efficient estimation technique that produces a factor loading matrix with many zero elements thanks to the introduction of sparsity-inducing penalty functions within the estimation process. In order to achieve sparse solutions and stable model selection procedures, the penalty functions must be non-differentiable. In this work, we adopted suitable local approximations of them. In this way, it was possible to employ in the optimization process a trust-region algorithm, which required analytical information on the score vector and the Hessian matrix (or a good approximation thereof). The use of differentiable penalties allowed us to recast the problem in a theoretically founded framework, where a precise definition of effective degrees of freedom was obtained, based on the bias term of the Generalized Information Criterion, or equivalently, the influence matrix of the model. This represents a novelty, as the existing proposals compute the degrees of freedom of a penalized factor model as the number of nonzero parameters. As an alternative to the usually time-consuming grid-searches, we also illustrated an efficient automatic technique for the estimation of the tuning parameter alongside the parameters of the factor model. The asymptotic properties of the penalized estimator can be established along the lines of Filippou et al. (Reference Filippou, Marra and Radice2017) and Fan and Li (Reference Fan and Li2001).

The simulations showed that the proposed approach produced trustworthy models with high accuracy, selection consistency, low bias and false positives. This indicates that the method is a valuable alternative to the existing techniques. Furthermore, it often generated the best tradeoff between goodness of fit and model complexity when compared to such models, as in the empirical application. As a result of this delicate balance, the proposed method may not necessarily provide the sparsest factor solution. Numerical experiments, however, confirm that the proposed method can produce very good results even if the penalized parameters are estimated just close enough to zero. This is because the edf are also being estimated close to zero, and we would actually need a considerable number of coefficients to see a substantive impact on the total edf and the GBIC. Still, if researchers desire more sparsity, they can manually and subjectively increase the value of the tuning parameter or the influence factor for the automatic procedure.

Notably, we extended the illustrated framework to multiple-group factor models by employing a penalty that simultaneously induced sparsity and cross-group equality of loadings and intercepts. As such, it revealed as a worthy alternative to invariance testing procedures. In this context, the automatic procedure proved particularly useful as it allowed for the estimation of the multiple tuning parameters composing the penalty term in a fast, stable and efficient way.

The presented framework allows one to easily and efficiently combine multiple penalty terms (like in the multiple-group model), as the automatic procedure scales well with the number of tuning parameters. In the empirical application, the alasso penalty was considered for all three penalty terms, but different penalty functions can also be combined if desired.

Another interesting modification pertains to the type of parameters that are penalized. Given the general estimation framework proposed in this work, also residual covariances (i.e., the off-diagonal elements of the covariance matrix of the unique factors) can be penalized to examine the assumption of conditional independence (that is, detect which pairs of variables are conditionally dependent). This model is known in the econometric literature as “sparse approximate factor model” (Bai & Liao, Reference Bai and Liao2016).

We envisage several interesting lines of future research. Firstly, the proposed approach can be applied to structural equation models in which, in addition to the measurement model, a structural model (usually a mediation model for the factors) is tested. Secondly, the results described in this work were derived under the $N > p$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N > p$$\end{document} scenario, as it is the case for many applications from the social and behavioral sciences. However, penalized techniques can also be extremely useful in the high-dimensional case, where maximum likelihood estimation is not feasible. It would hence be interesting to review the presented methodology in this demanding set-up. Future research may also evaluate the impact of messy data and larger model sizes on the penalized estimation framework. Finally, the observed variables were assumed to follow a multivariate normal distribution. When this is not reasonable, one can resort to pseudo maximum likelihood (Arminger & Schoenberg, Reference Arminger and Schoenberg1989) or, for categorical data, pairwise maximum likelihood (Katsikatsou et al., Reference Katsikatsou, Moustaki, Yang-Wallentin and Jöreskog2012). Further studies are needed to extend this work to the non-normal case, as this setting poses additional challenges since the asymptotic covariance matrix of the PMLE is no longer consistently estimated by the inverse Fisher information but by a “sandwich-type” covariance matrix (Yuan & Bentler Reference Yuan and Bentler1997).

Funding

Open access funding provided by Alma Mater Studiorum - Universitá di Bologna within the CRUI-CARE Agreement.

Footnotes

Supplementary Information The online version contains supplementary material available at https://doi.org/10.1007/s11336-021-09751-8.

1 An alternative factor model formulation would include the intercepts of the observed variables and the factor means. See the multiple-group extension in Sect. 4 for an example of a mean and covariance structure model.

2 For regsem, we used the default optimizer Rsolnp for lasso and alasso, and coordinate descent for scad and mcp. The additional tuning parameters of the penalties were kept to the values specified in the package, that is, a = 1 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = 1$$\end{document} for alasso, and a = 3.7 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = 3.7$$\end{document} for scad and mcp. For lslx, as per software implementations, the shape parameter of the mcp was internally selected over a varying three-dimensional grid of values adapted for each replicate. Because of the specificities of each package implementation, a customized and sensible grid of η \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} values was considered for each technique.

3 A solution is considered admissible if it does not present Heywood cases (negative unique variances), the covariance matrices of the unique factors and common factors are positive-definite, the factor loading matrix is of full column rank and does not contain any null rows (Jöreskog & Sörbom, Reference Jöreskog and Sörbom1996).

4 No rounding was necessary for TPR because, based on the simulation design, the estimates were never mistakenly estimated close to zero.

5 This choice was made on practical grounds, deeming estimates in absolute value below 0.05 to be “practically” and “substantively” equal to zero. The use of tighter rounding thresholds may worsen FPR and PCTM. For penfa, numerical experiments showed that this was the case when using the grid-search, whereas the models with the automatic procedure exhibited fairly comparable FPR and PCTM even after two or three decimal digits roundings.

6 All computations were carried out on a machine with Intel(R) Core(TM) i7-5600U 2.60GHz (quad-core) processor and 16GB of RAM.

7 We used the BIC as a criterion for model comparisons due to its widespread use in sparse settings, but different evaluation measures can be employed depending on the research question.

8 This was done for a neater visual illustration; in the penfa package, no internal rounding is implemented, and the estimates are shown as returned by the trust-region optimizer.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

Akaike, H. (1974). A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19(6), 716–723.CrossRef Google Scholar

Arminger, G., & Schoenberg, R. J. (1989). Pseudo maximum likelihood estimation and a test for misspecification in mean and covariance structure models. Psychometrika, 54(3), 409–425.CrossRef Google Scholar

Bai, J., & Liao, Y. (2016). Efficient estimation of approximate factor models via penalized maximum likelihood. Journal of Econometrics, 191(1), 1–18.CrossRef Google Scholar

Bauer, D. J., Belzak, WCM, & Cole, V. T. (2020). Simplifying the assessment of measurement invariance over multiple background variables: Using regularized moderated nonlinear factor analysis to detect differential item functioning. Structural Equation Modeling: A Multidisciplinary Journal, 27(1), 43–55.CrossRef Google Scholar PubMed

Browne, M. W. (2001). An overviewof analytic rotation in exploratory factor analysis. Multivariate Behavioral Research, 36(1), 111–150.CrossRef Google Scholar

Chen, J., & Chen, Z. (2008). Extended Bayesian information criteria for model selection with large model spaces. Biometrika, 95(3), 759–771.CrossRef Google Scholar

Choi, J., Oehlert, G., & Zou, H. (2010). A penalized maximum likelihood approach to sparse factor analysis. Statistics and its Interface, 3(4), 429–436.CrossRef Google Scholar

Chou, C., & Bentler, P. M. (1990). Model modification in covariance structure modeling: A comparison among likelihood ratio, Lagrange multiplier, and Wald tests. Multivariate Behavioral Research, 25(1), 115–136.CrossRef Google Scholar PubMed

Chou, C., Huh, J., & Hoyle, R. H. (2012). Model modification in structural equation modeling. Handbook of structural equation modeling, New York, NY:The Guilford Press. 232–246.Google Scholar

Fan, J., & Li, R. (2001). Variable selection via nonconcave penalized likelihood and its oracle properties. Journal of the American Statistical Association, 96(456), 1348–1360.CrossRef Google Scholar

Filippou, P., Marra, G., & Radice, R. (2017). Penalized likelihood estimation of a trivariate additive probit model. Biostatistics, 18(3), 569–585.CrossRef Google Scholar PubMed

Geminiani, E. (2020). A penalized likelihood-based framework for single and multiple-group factor analysis models (Doctoral dissertation, University of Bologna). Retrieved from http://amsdottorato.unibo.it/9355/.Google Scholar

Hair, J. F.,Black, W. C.,Babin, B. J., & Anderson, R. E. (2010). Multivariate data analysis, 7. Upper Saddle River, NJ: Prentice Hall.Google Scholar

Hirose, K., & Yamamoto, M. (2014a). Estimation of an oblique structure via penalized likelihood factor analysis. Computational Statistics & Data Analysis, 79, 120–132.CrossRef Google Scholar

Hirose, K., & Yamamoto, M. (2014b). Sparse estimation via nonconcave penalized likelihood in factor analysis model. Statistics and Computing, 25(5), 863–875.CrossRef Google Scholar

Holzinger, K. J. & Swineford, F. (1939). A study in factor analysis: The stability of a bi-factor solution. Supplementary Educational Monographs, 48, University of Chicago.Google Scholar

Huang, P. (2018). A penalized likelihood method for multi-group structural equation modelling. British Journal of Mathematical and Statistical Psychology, 71(3), 499–522.CrossRef Google Scholar PubMed

Huang, P. (2020). lslx: Semi-confirmatory structural equation modeling via penalized likelihood. Journal of Statistical Software, 93(7), 1–37.CrossRef Google Scholar

Huang, P. & Hu, W. (2019). lslx: Semi-confirmatory structural equation modeling via penalized likelihood [Computer software manual]. Retrieved from https://CRAN.R-project.org/package=lslx (R package version 0.6.8).Google Scholar

Huang, P., Chen, H., & Weng, L. (2017). A penalized likelihood method for structural equation modeling. Psychometrika, 82(2), 329–354.CrossRef Google Scholar PubMed

Jacobucci, R., Grimm, K. J., & McArdle, J. J. (2016). Regularized structural equation modeling. Structural Equation Modeling: A Multidisciplinary Journal, 23(4), 555–566.CrossRef Google Scholar PubMed

Jacobucci, R., Grimm, K. J., Brandmaier, A. M., Serang, S., Kievit, R. A. & Scharf, F. (2019). regsem: Regularized Structural Equation Modeling [Computer software manual]. Retrieved from https://CRAN.R-project.org/package=regsem (R package version 1.3.2).Google Scholar

Jin, S., Moustaki, I., & Yang-Wallentin, F. (2018). Approximated penalized maximum likelihood for exploratory factor analysis: An orthogonal case. Psychometrika, 83(3), 628–649.CrossRef Google Scholar

Jöreskog, K. G. Jöreskog, K. G., Sörbom, D., & Magidson, J. (1979). A general approach to confirmatory maximum likelihood factor analysis with addendum. Advances in factor analysis and structural equation models, Cambridge, MA:Abt Books. 21–43.Google Scholar

Jöreskog, K. G. & Sörbom, D. (1993). LISREL 8: Structural equation modeling with the SIMPLIS command language. Scientific Software International.Google Scholar

Jöreskog, K. G., & Sörbom, D. (1996). LISREL 8: User’s reference guide. Scientific Software. International.Google Scholar

Katsikatsou, M., Moustaki, I., Yang-Wallentin, F., & Jöreskog, K. G. (2012). Pairwise likelihood estimation for factor analysis models with ordinal data. Computational Statistics & Data Analysis, 56(12), 4243–4258.CrossRef Google Scholar

Kelley, K. (2019). MBESS: The MBESS R package [Computer software manual]. Retrieved from https://CRAN.R-project.org/package=MBESS.Google Scholar

Kim, Y., & Gu, C. (2004). Smoothing spline gaussian regression: More scalable computation via efficient approximation. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 66(2), 337–356.CrossRef Google Scholar

Koch, I. (1996). On the asymptotic performance of median smoothers in image analysis and nonparametric regression. The Annals of Statistics, 24(4), 1648–1666.CrossRef Google Scholar

Konishi, S., & Kitagawa, G. (1996). Generalised information criteria in model selection. Biometrika, 83(4), 875–890.CrossRef Google Scholar

Lindstrøm, J. C., & Dahl, F. A. (2020). Model selection with lasso in multi-group structural equation models. Structural Equation Modeling: A Multidisciplinary Journal, 27(1), 33–42.CrossRef Google Scholar

Little, T. D., Slegers, D. W., & Card, N. A. (2006). A non-arbitrary method of identifying and scaling latent variables in SEM and MACS models. Structural Equation Modeling: A Multidisciplinary Journal, 13(1), 59–72.CrossRef Google Scholar

Little, T. D., Card, N. A., Slegers, D. W., & Ledford, E. C.etalLittle, T. D., Bovaird, J. A., Card, N. A.etal (2012). Representing contextual effects in multiple-group MACS models. Modeling contextual effects in longitudinal studies, NY:Routledge New York. 121–147.Google Scholar

Lu, Z., Chow, S., & Loken, E. (2016). Bayesian factor analysis as a variable-selection problem: Alternative priors and consequences. Multivariate Behavioral Research, 51(4), 519–539.CrossRef Google Scholar PubMed

MacCallum, R. C., Roznowski, M., & Necowitz, L. B. (1992). Model modifications in covariance structure analysis: the problem of capitalization on chance. Psychological Bulletin, 111(3), 490–504.CrossRef Google Scholar PubMed

Marra, G., Radice, R. (2020). Copula link-based additive models for right-censored event time data. Journal of the American Statistical Association, 115(530), 886–895.CrossRef Google Scholar

Marra, G., & Wood, S. N. (2012). Coverage properties of confidence intervals for generalized additive model components. Scandinavian Journal of Statistics, 39(1), 53–74.CrossRef Google Scholar

Meredith, W. (1993). Measurement invariance, factor analysis and factorial invariance. Psychometrika, 58(4), 525–543.CrossRef Google Scholar

Millsap, R. E. (2001). When trivial constraints are not trivial: The choice of uniqueness constraints in confirmatory factor analysis. Structural Equation Modeling: A Multidisciplinary Journal, 8(1), 1–17.CrossRef Google Scholar

Millsap, R. E. (2012). Statistical approaches to measurement invariance, Abingdon:Routledge.CrossRef Google Scholar

Mulaik, S. A. (2009). Foundations of factor analysis, Boca Raton:Chapman and Hall/CRC.CrossRef Google Scholar

Muthén, B., Asparouhov, T. (2012). Bayesian structural equation modeling: A more flexible representation of substantive theory. Psychological Methods, 17(3), 313–335.CrossRef Google Scholar PubMed

Nocedal, J., & Wright, S. (2006). Numerical optimization, Berlin:Springer Science & Business Media.Google Scholar

R Core Team. (2018). R:Alanguage and environment for statistical computing [Computer software manual]. Vienna, Austria. Retrieved from https://www.R-project.org/.Google Scholar

Radice, R.,Marra, G., & Wojtyś, M. (2016). Copula regression spline models for binary outcomes. Statistics and Computing, 26(5), 981–995.CrossRef Google Scholar

Schwarz, G. (1978). Estimating the dimension of a model. The Annals of Statistics, 6(2), 461–464.CrossRef Google Scholar

Tang, Z., Shen, Y., Zhang, X., & Yi, N. (2017). The spike-and-slab lasso generalized linear models for prediction and associated genes detection. Genetics, 205(1), 77–88.CrossRef Google Scholar PubMed

Thurstone, L. L. (1947). Multiple-factor analysis: A development and expansion of the vectors of mind, Chicago, IL:University of Chicago Press.Google Scholar

Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B (Methodological), 58(1), 267–288.CrossRef Google Scholar

Trendafilov, N. T., Fontanella, S., & Adachi, K. (2017). Sparse exploratory factor analysis. Psychometrika, 82(3), 778–794.CrossRef Google Scholar

Ulbricht, J. (2010). Variable selection in generalized linear models. (Doctoral dissertation, Ludwig-Maximilians-Universität München). Verlag Dr. Hut.Google Scholar

Vandenberg, R. J., & Lance, C. E. (2000). A review and synthesis of the measurement invariance literature: Suggestions, practices, and recommendations for organizational research. Organizational Research Methods, 3(1), 4–70.CrossRef Google Scholar

Wood, S. N. (2004). Stable and efficient multiple smoothing parameter estimation for generalized additive models. Journal of the American Statistical Association, 99(467), 673–686.CrossRef Google Scholar

Wood, S. N. (2017). Generalized additive models: An introduction with R, Boca Raton: Chapman and Hall/CRC.CrossRef Google Scholar

Yuan, K., & Bentler, P. M. (1997). Improving parameter tests in covariance structure analysis. Computational Statistics & Data Analysis, 26(2), 177–198.CrossRef Google Scholar

Yuan, K. & Bentler, P. M. (2006). Structural equation modeling. In C. Rao & S. Sinharay (Eds.), Handbook of Statistics (Vol. 10, pp. 297-358). Elsevier.CrossRef Google Scholar

Zhang, C. (2010). Nearly unbiased variable selection under minimax concave penalty. The Annals of Statistics, 38(2), 894–942.CrossRef Google Scholar

Zou, H. (2006). The adaptive lasso and its oracle properties. Journal of the American Statistical Association, 101(476), 1418–1429.CrossRef Google Scholar

Zou, H., Hastie, T., & Tibshirani, R. (2007). On the “degrees of freedom” of the lasso. The Annals of Statistics, 35(5), 2173–2192.CrossRef Google Scholar

Table 1. Performance measures of the examined models in simulation study I by varying the sample size N.

Table 2. The factor loading matrices and intercepts of the two groups under each difference scenario of simulation study II.

Table 3. Performance measures of penfa-alasso (a=2,γ=4.5\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$a = 2, \gamma = 4.5$$\end{document}) and lslx-mcp models by sample size and difference scenario.

Table 4. BIC of the best configurations of the fitted models.

Table 5. Parameter estimates of the nine mental tests from the Holzinger & Swineford data set for the unpenalized model, and penfa-alasso with automatic procedure (on the left-hand side, η^=0.017,a=1\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{\eta }} = 0.017, a = 1$$\end{document} and γ=4.5\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma = 4.5$$\end{document}; on the right-hand side, η^=0.011,a=2\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{\eta }} = 0.011, a = 2$$\end{document} and γ=5.5\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma = 5.5$$\end{document}).

Table 6. Parameter estimates of the 19 mental tests from the Holzinger & Swineford data set for penfa-alasso (automatic procedure, η^=(0.006,16221.852,0.013)T,a=1,γ=4\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\hat{{\varvec{\eta }}} = (0.006, 16221.852, 0.013)^T, a = 1, \gamma = 4$$\end{document})

Geminiani et al. supplementary material

Online Resources for “Single and multiple-group penalized factor analysis: a trust-region algorithm approach with integrated automatic multiple tuning parameter selection”

File 1 MB

Article contents

Single- and Multiple-Group Penalized Factor Analysis: A Trust-Region Algorithm Approach with Integrated Automatic Multiple Tuning Parameter Selection

Abstract

Keywords

Information

1. Introduction

2. The Normal Linear Factor Analysis Model

3. Sparsity-Inducing Penalties

3.1. Locally Approximated Penalties

4. The Multiple-Group Factor Analysis Model

5. Sparsity and Invariance-Inducing Penalties

6. Generalized Information Criterion

7. Penalized Maximum Likelihood Estimation

8. Simulation Studies

8.1. Simulation Study I

8.2. Simulation Study II

9. Empirical Application

9.1. Normal Linear Factor Model

9.2. Multiple-Group Factor Model

10. Discussion

Funding

Footnotes

References

Geminiani et al. supplementary material

Save article to Kindle

Save article to Dropbox

Save article to Google Drive

Reply to: Submit a response

Your details

You have entered the maximum number of contributors

Conflicting interests