Normal-linear factor model

Data

\(x_1,x_2,\ldots,x_n \in \Re^p\)

Model

\[X_i = \mu_0 + \Lambda_0 Y_i + E_i\]

Normal-linear factor model

Log-likelihood

\(X_i \sim \mathrm{N}(\mu_0, \Sigma_0)\), for \(\Sigma_0 = \Lambda_0 {\Lambda_0}^\top + \Psi_0\) so that

\[ \ell(\Sigma; S) = - \frac{n}{2} \left \{\log \mathrm{det}(\Sigma) + \mathrm{tr}\left(\Sigma^{-1}S \right) \right\} + C \,, \] where \(S = 1/n \sum_{i = 1}^n (x_i - \bar{x})(x_i - \bar{x})^\top\), \(\bar{x} = 1/n \sum_{i = 1}^n x_i\)

Maximum likelihood (ML) estimator

\(\hat{\Lambda}, \hat{\Psi} = \arg \max \,\ell(\Lambda \Lambda^\top + \Psi; S)\)

for \(\Lambda \in \Re^{p \times q}\), \(\Psi \in \Re^{p \times p}\) diagonal and positive definite

Job application data

Data from Kendall (1980) which scores job applicants \((n = 48)\) according to a range of criteria \((p = 15)\)

Want to use exploratory factor analysis to find common factors underlying the data

Exploratory Factor Analysis with R

Most popular R packages for ML estimation of EFA on CRAN^*:

Exploratory Factor Analysis with R

stats::factanal

factanal_fit <- stats::factanal(x = job_data, factors = 4)

factanal_uniquenesses

Exploratory Factor Analysis with R

lavaan::efa

lavaan_fit <- efa(data = job_data, nfactors = 4)

lavaan_uniquenesses

Exploratory Factor Analysis with R

psych::fa

fa_fit <- fa(job_data, nfactors = 4)

fa_fit$uniquenesses

Exploratory Factor Analysis with R

psych::fa

fa_fit <- fa(job_data, nfactors = 4)

fa_fit$uniquenesses

The ML estimate does not exist

This phenomenon is termed a Heywood (1931) case

Nonexistence of ML estimates

Dominating parameter sequence

\(\{ (\Lambda_n, \Psi_n) \}_{n \in \mathbb{N}}\):

• for all \((\Lambda, \Psi)\), there is a \(N \in \mathbb{N}\) such that for all \(n > N\), \(\ell(\Lambda_n, \Psi_n) > \ell(\Lambda, \Psi)\)

Convergence to boundary

• \(\lim_\limits{n \to \infty} \min_{i} [\Psi_n]_{ii} = 0\)

Nonexistence of ML estimates

Dominating parameter sequence

Profiled boundary path \((\Lambda_{\psi},\Psi_{\psi})\) from maximising \(\ell(\Lambda,\Psi)\) at fixed \(\psi_{\mathrm{Keenness}}=\psi\)

Maximum penalized likelihood

Add penalty function \(P(\Lambda, \Psi)\) that diverges to \(-\infty\) whenever \(\Psi \to 0\), or \(\Psi,\Lambda \to \infty\)

Penalized log-likelihood

\[\ell^*(\Lambda, \Psi; S) = \ell(\Lambda \Lambda^\top + \Psi; S) + P(\Lambda, \Psi)\]

Maximum penalized likelihood (MPL) estimator

\[\bar{\Lambda}, \bar{\Psi} \in \arg \max \, \ell^*(\Lambda, \Psi; S) \tag{1}\]

Maximum penalized likelihood

Discouraging dominating parameter sequences

Existence

Assume:

• \(S\) is full rank
• \(P(\Lambda,\Psi)\) is continuous and bounded from above
• \(P(\Lambda_n,\Psi_n)\to-\infty\) for any convergent sequence \((\Lambda_n,\Psi_n)\) with positive uniquenesses such that

\[ \lim_{n\to\infty}\min_i[\Psi_n]_{ii}=0, \qquad \lim_{n\to\infty} \lambda_{\min}(\Lambda_n\Lambda_n^\top+\Psi_n)>0 \]

Soft penalisation

Penalisation distorts MPL estimates away from ML estimates

Scale penalty by sequence \(c_n\) to make \(c_n P(\cdot, \cdot)\) asymptotically negligible

Maximum softly-penalized likelihood (MSPL) estimator

\[\bar{\Lambda}, \bar{\Psi} \in \arg \max \, \ell(\Lambda \Lambda^\top + \Psi; S) + c_n P(\Lambda, \Psi)\]

Consistency

Assume:

• The set of MPL estimates is nonempty
• \(P(\Lambda, \Psi)\) is non-positive
• \(\Sigma_0 = \Lambda_0 \Lambda_0^\top + \Psi_0\) is strongly identifiable in that for any \(\epsilon>0\) there exists a \(\delta>0\) such that \[ ||| \Sigma_0 - \Sigma |||_{\max} < \delta \implies ||| \Lambda_0 - \Lambda Q|||_{\max} < \epsilon, \, ||| \Psi_0 - \Psi|||_{\max} < \epsilon \,, \] for some \(q \times q\) orthogonal matrix \(Q\).

Consistency

If \(S \longrightarrow \Sigma_0\) and \(P(\Lambda,\Psi)\) is an \(O(1)\) function, then \(c_n=o(n)\) is sufficient for consistency

Consistency

If \(S \longrightarrow \Sigma_0\) and \(P(\Lambda,\Psi)\) is an \(O(1)\) function, then \(c_n=o(n)\) is sufficient for consistency

Can be strengthened to \[ ||| \sqrt{n}(\hat{\Lambda}Q_n - \bar \Lambda O_n) |||_{\max} = o_p(1), \qquad |||\sqrt{n}(\hat \Psi - \bar \Psi) |||_{\max} = o_p(1), \] for rotation matrices \(O_n\) and \(Q_n\)

Choosing \(c_n\)

Independence model

\[ \mu_0 = 0_p, \qquad \Lambda_0 = 0_{p\times q}, \qquad \Psi_0 = \operatorname{diag}(\sigma_1^2,\ldots,\sigma_p^2). \]

Each coordinate reduces to a univariate normal model: \[ X_{ij} \sim N(0,\sigma_j^2) \,. \]

Rate of information accumulation

\[ I_j(\sigma_j^2) = -\mathbb E\left[ \frac{\partial^2 \ell_j}{\partial(\sigma_j^2)^2} \right] = \frac{n}{2\sigma_j^4}, \qquad \sqrt{\frac{n}{2}} \left( \frac{\hat\sigma_j^2}{\sigma_j^2}-1 \right) \longrightarrow N(0,1) \,. \]

Root-inverse-information scale is \(c_n = \sqrt{2/n}\)

Candidate penalties

Akaike (1987)

\[ P(\Lambda, \Psi) = -\frac{\rho n}{2} \mathrm{tr}\left(\Psi^{-1/2} \Lambda \Lambda^\top \Psi^{-1/2}\right) \]

Hirose et al. (2011)

\[ P(\Psi) = -\frac{\rho n}{2} \mathrm{tr}\left(\Psi^{-1/2} S \Psi^{-1/2}\right) \]

Both are scale-equivariant under response rescaling and rotation-invariant under \(\Lambda \mapsto \Lambda Q\).

Setup

Based on Cooperman & Waller (2022)

• \(q = 3\)

• \(n \in \{50, 100, 400\}\)

• item-to-factor ratio \(3:1, 5:1, 8:1\)

Methods

Setup

% Heywood cases

Frequentist performance

Model selection