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)\)

Job application data

Data from Hirose et al. (2011) 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

Job application data

Fitting a normal-linear factor model with 4 factors in R using the psych package yields

fa_fit <- fa(job_data, nfactors = 4)

Estimated specific variances

fa_fit$uniquenesses

The ML estimate does not exist

Heywood cases

Dominating parameter sequence

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

Divergence to boundary

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

Depending optimisation routine and constraints, observe small or negative values of \(\hat{\Psi}\)

This phenomenon is termed a Heywood (1931) case

Maximum penalized likelihood

Add penalty function \(P(\Psi, \Lambda)\) 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} = \arg \max \, \ell^*(\Lambda, \Psi; S) \tag{1}\]

Existence

• \(S\) is full rank
• \(P(\Lambda, \Psi)\) is continuous and bounded from above
• \(\lim_\limits{n \to \infty} P(\Lambda_n,\Psi_n) = - \infty\) for any sequence \(\{\Lambda_n, \Psi_n\}_{n \in \mathbb{N}}\) such that \(\min_{i} [\Psi_n]_{ii} \to 0\) or either \(\max_{i} [\Psi_n]_{ii} \to \infty\), \(|||\Lambda_n|||_{\max} \to \infty\)

Soft penalisation

Distort MPL estimates away from ML estimates through penalisation

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

Maximum softly-penalized likelihood (MSPL) estimator

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

Consistency

• \(\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 O|||_{\max} < \epsilon, \, ||| \Psi_0 - \Psi|||_{\max} < \epsilon \,, \] for some orthogonal matrix \(O\) of order \(q\)

• \(P(\Lambda, \Psi)\) is non-positive

Asymptotic normality

• \((\hat{\Lambda}, \hat{\Psi})\), \((\bar{\Lambda}, \bar{\Psi})\) are consistent for \((\Lambda_0, \Psi_0)\)

• \({\sup}_{\theta \in \mathcal{N}_{0}} \, ||| n^{-1}\nabla \nabla^\top \ell(\theta) - J(\theta)|||_{\max} = o_{p}(1)\), where \(\theta = (\textrm{vec}(\Lambda)^\top, \textrm{diag}(\Psi)^\top)^\top\), \(J(\theta)\) is continuous and invertible in a neighbourhood \(\mathcal{N}_0\) around \(\theta_0\)

• \(P(\theta)\) is differentiable in a neighbourhood \(\mathcal{N}_0\) around \(\theta_0\) a

Candidate penalties

Huber loss

\[ P(\Psi) =- \sum_{j =1}^p g(\Psi_{ii}), \quad g(x) = \begin{cases} \frac{1}{2} \log(x)^2 & \textrm{if } |\log(x)| < 1, \\ |\log(x)| - \frac{1}{2} & \textrm{else} \end{cases} \]

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) \]

Simulation setup

Huber loss penalty with \(c_n = n^{-1/2}\) on \(\Psi\) (Huber[\(\Psi\)]) and \(\Psi,\Lambda \Lambda^\top\) (Huber[\(\Psi, \Lambda \Lambda'\)])

Akaike (1987), Hirose et al. (2011) penalties with \(\rho \in \{1, n^{-3/2} \}\)

% Used

Bias

Probabiliy of underestimation

Root mean-squared error

Model selection

Take away

Maximum softly-penalized likelihood framework is a robust alternative to ML estimation in factor analysis with desirable theoretical properties and strong finite-sample performance

• Penalisation guarantees existence of estimates

• Soft scaling preserves ML asymptotics

• Strong finite-sample & model selection performance

• Simple Huber loss penalty favourable to penalties of Hirose et al. (2011) and Akaike (1987)