Factor Analysis

Factor analysis posits a generative model: observed data arise from latent factors through the equation x = Λz + ε. But how do we estimate the parameters—the factor loadings Λ and uniquenesses Ψ—from observed data?

Maximum Likelihood (ML) estimation is the gold standard approach. It finds parameter values that maximize the probability of observing the data we actually observed, given the factor model. ML estimation provides:

• Asymptotically efficient estimates: Among consistent estimators, ML has the smallest variance as sample size grows
• Standard errors: Allowing inference about whether loadings differ from zero
• Goodness-of-fit tests: Chi-square statistics to assess model adequacy
• Model comparison criteria: AIC, BIC for selecting number of factors

This page develops the ML estimation framework from first principles, covering the likelihood function, optimization procedures, inference, and model assessment. We'll also compare ML with alternative estimation methods.

The Multivariate Normal Likelihood

Under the factor model with Gaussian assumptions, the marginal distribution of observed x is:

$$\mathbf{x} \sim \mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Sigma})$$

where the model-implied covariance is:

$$\boldsymbol{\Sigma} = \mathbf{\Lambda \Lambda}' + \boldsymbol{\Psi}$$

Given n independent observations x₁, x₂, ..., xₙ, the likelihood is the product of individual densities:

$$L(\boldsymbol{\mu}, \mathbf{\Lambda}, \boldsymbol{\Psi} | \mathbf{X}) = \prod_{i=1}^{n} (2\pi)^{-p/2} |\boldsymbol{\Sigma}|^{-1/2} \exp\left( -\frac{1}{2} (\mathbf{x}_i - \boldsymbol{\mu})' \boldsymbol{\Sigma}^{-1} (\mathbf{x}_i - \boldsymbol{\mu}) \right)$$

The Log-Likelihood

Taking the natural logarithm (and ignoring constants):

$$\ell(\boldsymbol{\mu}, \mathbf{\Lambda}, \boldsymbol{\Psi}) = -\frac{n}{2} \left[ \ln |\boldsymbol{\Sigma}| + \text{tr}(\boldsymbol{\Sigma}^{-1} \mathbf{S}) + (\bar{\mathbf{x}} - \boldsymbol{\mu})' \boldsymbol{\Sigma}^{-1} (\bar{\mathbf{x}} - \boldsymbol{\mu}) \right]$$

where:

• S is the sample covariance matrix: $\mathbf{S} = \frac{1}{n} \sum_{i=1}^{n} (\mathbf{x}_i - \bar{\mathbf{x}})(\mathbf{x}_i - \bar{\mathbf{x}})'$
• $\bar{\mathbf{x}}$ is the sample mean

The ML estimate of μ is simply $\hat{\boldsymbol{\mu}} = \bar{\mathbf{x}}$ (the sample mean). For covariance structure, we can work with centered data.

Concentrating Out the Mean

Since $\hat{\boldsymbol{\mu}} = \bar{\mathbf{x}}$ regardless of Λ and Ψ, we can work with the concentrated log-likelihood that depends only on the covariance parameters:

$$\ell^*(\mathbf{\Lambda}, \boldsymbol{\Psi}) = -\frac{n}{2} \left[ \ln |\boldsymbol{\Sigma}| + \text{tr}(\boldsymbol{\Sigma}^{-1} \mathbf{S}) \right]$$

This focuses on matching the model-implied covariance Σ = ΛΛ' + Ψ to the sample covariance S.

Parameter Count and Identification

Free parameters:

• Λ: p × k loadings
• Ψ: p uniquenesses (diagonal elements)
• Total: pk + p parameters

However, the model is not identified without constraints because of rotational indeterminacy. Standard constraint: require Λ'Ψ⁻¹Λ to be diagonal. This fixes k(k-1)/2 parameters, giving:

Net free parameters: pk + p - k(k-1)/2

Degrees of freedom for the model:

$$df = \frac{(p)(p+1)}{2} - \left[ pk + p - \frac{k(k-1)}{2} \right] = \frac{(p-k)^2 - (p+k)}{2}$$

For the model to be testable, we need df > 0.

Iterative Optimization

The ML estimates $\hat{\mathbf{\Lambda}}$ and $\hat{\boldsymbol{\Psi}}$ are found by maximizing the log-likelihood iteratively. There is no closed-form solution; numerical optimization is required.

Common algorithms:

1. Newton-Raphson: Uses first (gradient) and second (Hessian) derivatives

   • Fast convergence near the optimum
   • Requires Hessian computation (expensive)

2. Fisher Scoring: Similar to Newton-Raphson but uses expected Hessian

   • More stable in some cases
   • Standard in many FA implementations

3. EM Algorithm: Treats factors as missing data

   • E-step: Compute expected factor scores
   • M-step: Update Λ and Ψ
   • Guaranteed to increase likelihood each iteration
   • May converge slowly

4. Quasi-Newton (BFGS, L-BFGS): Approximates Hessian

   • Good for high-dimensional problems
   • Less per-iteration computation

Initialization and Convergence

Starting values matter. Poor initialization can lead to:

• Slow convergence
• Convergence to a local (not global) maximum
• Non-convergence

Common initialization strategies:

• Principal components (use PCA loadings scaled by eigenvalue)
• Prior reliable communality estimates
• Random starts with selection of best solution

The EM Algorithm for Factor Analysis

The EM algorithm treats the latent factors z as "missing data":

E-step: Given current Λ⁽ᵗ⁾, Ψ⁽ᵗ⁾, compute:

• Expected factor scores: E[z|x] = Λ'(ΛΛ' + Ψ)⁻¹(x - μ)
• Expected factor covariance: E[zz'|x]

M-step: Given expected sufficient statistics, update:

$$\mathbf{\Lambda}^{(t+1)} = \mathbf{S} \cdot \mathbf{\Beta}' \cdot (\mathbf{\Beta} \mathbf{S} \mathbf{\Beta}' + \mathbf{I} - \mathbf{\Beta} \mathbf{\Lambda}^{(t)})^{-1}$$

where β = Λ'Σ⁻¹ (factor score regression weights).

$$\psi_j^{(t+1)} = s_{jj} - \sum_m \lambda_{jm}^{(t+1)} \beta_{mj}$$

Iterate until convergence (change in log-likelihood < tolerance).

EM advantages: Monotonic likelihood increase, handles some missing data naturally. EM disadvantages: Can be slow; may converge to saddle points or local maxima.

Asymptotic Theory

Under regularity conditions, ML estimates are:

• Consistent: $\hat{\theta} \xrightarrow{p} \theta$ as n → ∞
• Asymptotically normal: $\sqrt{n}(\hat{\theta} - \theta) \xrightarrow{d} \mathcal{N}(0, \mathbf{I}^{-1})$
• Asymptotically efficient: Achieve the Cramér-Rao lower bound

The expected information matrix I is: $$\mathbf{I} = -E\left[ \frac{\partial^2 \ell}{\partial \theta \partial \theta'} \right]$$

The asymptotic covariance of estimates is I⁻¹/n.

Standard Errors for Loadings

For each loading λᵢⱼ, the standard error allows construction of confidence intervals and z-tests:

$$SE(\hat{\lambda}{ij}) = \sqrt{[\mathbf{I}^{-1}]{(i,j),(i,j)} / n}$$

Practical computation: Modern software computes the observed or expected information matrix and extracts standard errors from its inverse.

Testing Individual Loadings

To test H₀: λᵢⱼ = 0:

z-statistic: $$z = \frac{\hat{\lambda}{ij}}{SE(\hat{\lambda}{ij})}$$

Under H₀, z ~ N(0,1) asymptotically. A loading is "significant" if |z| > 1.96 (α = 0.05).

Robust Standard Errors

The normality assumption often fails. Robust (sandwich) standard errors account for non-normality:

$$\hat{\mathbf{V}}_{robust} = \mathbf{I}^{-1} \mathbf{B} \mathbf{I}^{-1}$$

where B is estimated from the outer product of score vectors, capturing the actual variability of the data.

When to use robust SEs:

• Multivariate non-normality (skewness, kurtosis)
• Violation of distributional assumptions
• Cautious practice in applied research

Software support: Most modern SEM software (Mplus, lavaan) offers robust SEs as an option (e.g., MLR, MLM estimators).

Bootstrap Standard Errors

An alternative to analytic SEs is bootstrapping:

1. Resample n observations with replacement
2. Re-estimate the factor model on the bootstrap sample
3. Repeat B times (e.g., B = 1000)
4. Standard error = Standard deviation of bootstrap estimates

Advantages:

• No distributional assumptions
• Works for complex statistics (factor scores, factor correlations after rotation)

Disadvantages:

• Computationally intensive
• Requires convergence on each bootstrap sample
• Confidence interval coverage may not be exact

Unlike PCA, factor analysis is a model that can be right or wrong. Fit assessment asks: does the hypothesized factor structure adequately explain the observed covariances?

The Likelihood Ratio Test

The fundamental test compares the factor model to a saturated model (which exactly reproduces S):

Test statistic: $$\chi^2 = n \cdot F_{ML}$$

where F_ML is the minimized fitting function.

Under H₀ (k factors are sufficient): $\chi^2 \sim \chi^2_{df}$

Degrees of freedom: df = [(p - k)² - (p + k)] / 2

Interpretation:

• Significant χ² (p < 0.05): Reject the k-factor model
• Non-significant χ² (p > 0.05): Do not reject; k factors may be sufficient

The Problem with χ²

The χ² test has well-known limitations:

• Sensitive to sample size: With large n, even trivial misfits are significant
• Sensitive to normality: Non-normality inflates χ²
• All models are wrong: Rejecting a useful model because it's not perfect is counterproductive

Thus, we supplement χ² with descriptive fit indices.

Information Criteria

For comparing models with different numbers of factors:

AIC (Akaike Information Criterion): $$AIC = -2 \ell + 2q$$

where q is the number of free parameters and ℓ is the maximized log-likelihood.

BIC (Bayesian Information Criterion): $$BIC = -2 \ell + q \ln(n)$$

Using AIC/BIC:

• Fit models with k = 1, 2, 3, ... factors
• Select the model with lowest AIC or BIC
• BIC penalizes complexity more heavily, tends to favor simpler models

Residual Analysis

Residual covariance matrix: $$\mathbf{R} = \mathbf{S} - \hat{\boldsymbol{\Sigma}}$$

where $\hat{\boldsymbol{\Sigma}} = \hat{\mathbf{\Lambda}}\hat{\mathbf{\Lambda}}' + \hat{\boldsymbol{\Psi}}$.

Standardized residuals: $r_{ij}^* = r_{ij} / \sqrt{var(r_{ij})}$

• Large residuals (|r*| > 2-3) indicate specific misfits
• Patterns in residuals suggest missing factors or local dependencies

This is analogous to regression residual analysis.

The Heywood Case Problem

A Heywood case occurs when:

• An estimated uniqueness ψᵢ ≤ 0 (impossible—variance cannot be negative)
• Equivalently, communality h²ᵢ ≥ 1 (impossible—cannot explain more than 100%)

This is an improper solution—parameter estimates that violate constraints inherent to the model.

Why Heywood Cases Occur

1. Model misspecification:

• Wrong number of factors (often too few)
• Missing factors that would account for the variance

2. Sampling error:

• Small samples produce unstable estimates
• Random fluctuation can push estimates to boundaries

3. Insufficient indicator variables:

• Factors defined by too few variables (1-2 indicators per factor)
• Identification problems

4. Near-zero true uniqueness:

• Some variables may genuinely have near-zero unique variance
• Slight sampling error pushes the estimate negative

Remedies for Heywood Cases

1. Constrain uniquenesses to be positive: Most software allows fixing ψᵢ ≥ ε (small positive number). This yields a proper solution but may mask underlying problems.

2. Try different number of factors: Add more factors: the Heywood case may disappear when the model has sufficient factors. Remove factors: if you've overfactored, simpler models may converge properly.

3. Examine the offending variable:

• Is the variable well-measured?
• Does it correlate extremely highly with others?
• Consider removing it from the analysis

4. Increase sample size: If sampling error is the cause, more data stabilizes estimates.

5. Use ridge FA or regularization: Adding a small constant to the diagonal of S can prevent ill-conditioning.

Ultra-Heywood Cases

Sometimes loadings become extremely large (implausible values like λ = 5). This is a sign of severe model problems:

• Complete model misspecification
• Collinearity issues
• Numerical instability

Action: Do not interpret such solutions. Revisit model specification fundamentally.

While ML is the gold standard, other estimation methods are available, each with advantages and drawbacks.

Principal Axis Factoring (PAF)

Also called Principal Factor Analysis or Iterated Principal Factor.

Procedure:

1. Replace diagonal of R with estimated communalities (initially 1 or SMC)
2. Extract eigenvalues and eigenvectors of the reduced matrix
3. Update communalities as sum of squared loadings
4. Repeat until communalities converge

Advantages:

• No distributional assumptions
• Computationally simpler than ML
• Sometimes converges when ML fails

Disadvantages:

• No standard errors
• No formal fit test
• May converge to different solution than ML

When to use: Exploratory analysis when normality is doubtful and formal testing is not needed.

Minimum Residual (MINRES)

Minimizes the sum of squared off-diagonal residuals:

$$\min_{\Lambda, \Psi} \sum_{i \neq j} (s_{ij} - \sigma_{ij})^2$$

Advantages:

• No distributional assumptions
• Focuses on reproducing correlations (off-diagonal)
• Robust to some types of misspecification

Disadvantages:

• Less efficient than ML under correct model
• No formal inference

Comparison of Estimation Methods

Which to Choose?

Use ML when:

• Distributional assumptions are approximately met
• You need standard errors and fit tests
• Sample size is adequate (n > 200 or so)

Use PAF/MINRES when:

• Data are clearly non-normal
• Exploratory analysis without formal testing
• ML doesn't converge
• Quick approximate solution is sufficient

Use robust ML (MLR, MLM) when:

• Some non-normality but you want inferential results
• Available in Mplus, lavaan, and other SEM software

Step-by-Step Procedure

1. Data Preparation:

• Check for missing data (consider imputation or FIML)
• Examine univariate distributions (skewness, kurtosis)
• Check multivariate normality (Mardia's test)
• Standardize if using correlation matrix

2. Determine Number of Factors:

• Eigenvalue analysis (Kaiser rule as starting point)
• Scree plot
• Parallel analysis (more reliable)
• Theory and interpretability

3. Run Initial ML Estimation:

• Extract k factors
• Check convergence
• Check for Heywood cases

4. Rotate Solution:

• Apply varimax or oblimin rotation
• Examine loading patterns

5. Assess Model Fit:

• χ² test (with caution)
• RMSEA, CFI, TLI, SRMR
• Examine residuals

6. Refine If Needed:

• Try k-1 or k+1 factors
• Consider removing poorly fitting items
• Re-estimate and compare fit

Common R Code (lavaan)

library(lavaan)

# Specify 3-factor EFA model
efa_model <- efa(
  data = mydata,
  nfactors = 3,
  rotation = "oblimin"
)

# Summary with fit indices
summary(efa_model, fit.measures = TRUE, standardized = TRUE)

# Get loadings
inspect(efa_model, what = "std")$lambda

# Get factor correlations
lavInspect(efa_model, what = "cov.lv")

Common Python Code (factor_analyzer)

from factor_analyzer import FactorAnalyzer

# Fit model
fa = FactorAnalyzer(n_factors=3, rotation='oblimin', method='ml')
fa.fit(data)

# Loadings
print(fa.loadings_)

# Communalities
print(fa.get_communalities())

# Eigenvalues
print(fa.get_eigenvalues())

Maximum likelihood estimation provides a principled, inferentially rich framework for factor analysis. Let's consolidate the key concepts:

What's Next

With estimation and model assessment covered, the final page addresses Interpretation—how to make sense of factor solutions, name factors, communicate findings, and connect factor analysis to broader research contexts.