Session 13 · One More Thing: Principal Component and Factor Analyses

After the course · for self-study · not examined

Thierry Warin, PhD

HEC Montréal

Outline · The hour ahead

  1. Introduction
  2. Principal Component Analysis
  3. The Netflix’s Recommender System through a PCA
  4. Factor analysis
  5. FAMD: the wine data
  6. Factor Analysis for Questionnaires

What session 13 is for

The mathematics

  • Explain why PCA requires standardised inputs, and what happens if you forget
  • Read a scree plot and defend the number of components retained
  • Distinguish a loading from a score
  • Say why a factor is identified only up to rotation

The Python

  • Compute components from the correlation matrix, and check against a library
  • Plot the scree and the loadings, and label the axes with variance explained
  • Rotate a factor solution and watch the fit stay identical

The international business

  • Build an index — competitiveness, openness, capability — from correlated indicators
  • Defend the name you gave a component to a sceptical reader
  • Say what a composite indicator hides as well as what it summarises

Your Python so far

Everything you need to type is on one sheet — how to start Python in VS Codium on your platform, how to make a file and run it, the DataFrame basics, and every name from sessions 1 to 13, each with a line showing it in use.

Python Cheatsheet · Session 13 (PDF) · in the repository at 13-…/PYTHON-CHEATSHEET.pdf

One sheet per package, functions explained with pictures: pandas · NumPy · matplotlib · prince · scikit-learn · statsmodels

Print it, or keep it open beside the editor. Nothing on it is examinable; all of it is assumed.

Quote of the day

"If you live 5 minutes away from Bill Gates, I bet you are rich."

Introduction

Where this session sits

Color (nm) Alcohol (%) Beer or Wine?
610 5 Beer
599 13 Wine
693 14 Wine
  • Labels

  • training / testing

  • supervised / unsupervised (or self-supervised) / reinforcement learning

  • System 1 (fast, intuitive) / System 2 (slow, analytical) (D. Kahneman)

The machine-learning workflow

  1. Gathering data

  2. Pre-processing

  3. Choosing a model

  4. Training the model / K-fold cross-validation

  5. Evaluation

  6. Hyper-parameter tuning

  7. Out-of-sample Prediction

Regression or classification?

  • A set of principles that can be applied to regression or classification

What is the difference between regression and classification models?

  • Regression predicts a number: a price, a volume, a share.

  • Classification predicts a category: default or not, beer or wine. A classifier draws a boundary between the classes.

Multivariate statistics

Multivariate Statistics is a group of statistical methods that focus on studying multiple variables together while focusing on the variation that those variables have in common.

Definition: Multivariate Statistics deals with the treatment of data sets with a large number of dimensions.

Its goals are therefore different from supervised modeling, but also different from segmentation and clustering models.

There are many models in the family of Multivariate Statistics.

  1. PCA
  1. Factor Analysis

Principal Component Analysis

Why PCA?

  • Modern datasets often have a very large number of variables. This makes it difficult to inspect each of the variables individually, due to the practical fact that the human mind cannot easily digest data on such a large scale.

  • When a dataset contains a large number of variables, there is often a serious amount of overlap between those variables.

  • The components that are found by PCA are ordered from the highest information content to the lowest information content.

PCA is a statistical method that allows you to “regroup” your variables into a smaller number of variables, called components.

This regrouping is done based on variation that is common to multiple variables.

Interesting book: https://shainarace.github.io/LinearAlgebra/pca.html

What PCA does

The goal of PCA is to regroup variables in such a way that the first (newly created) component contains a maximum of variation.

  • The second component contains the second-largest amount of variation, etc.

The last component logically contains the smallest amount of variation.

Thanks to this ordering of components, it is made possible to retain only a few of the newly created components, while still retaining a maximum amount of variation.

We can then use the components rather than the original variables for data exploration.

Mathematical Model — Maximize variance of the new components

The mathematical definition of the PCA problem is to find a linear combination of the original variables with maximum variance.

This means that we are going to create a (new) component. Let’s call it \(z\). This z is going to be computed based on our original variables \((X_1, X_2, …)\) multiplied by a weight for each of our variables \((u_1, u_2, …)\).

This can be written as \(z = Xu\).

[We apply it to cross-section data. For time-series cross-section data, we can use the superPCA package. Not exam material!]

The mathematical goal is to find the values for \(u\) that will maximize the variance of \(z\), with a constraint of unit length on \(u\).

This problem is mathematically called a constrained optimization using Lagrange Multiplier, but in practice, we use computers to do the whole PCA operation at once.

This can also be described as applying matrix decomposition to the correlation matrix of the original variables.

PCA is efficient in finding the components that maximize variance.

This is great if we are interested in reducing the number of variables while keeping a maximum of variance.

Sometimes, however, we are not purely interested in maximizing variance: we might want to give the most useful interpretations to our newly defined dimensions.

And this is not always easiest with the solution found by a PCA.

We can then apply Factor Analysis: an alternative to PCA that has a little bit more flexibility.

  • PCA only works with numeric data

  • Categorical data must be encoded as numeric data (e.g. "one-hot encoding" = creating dummies from categories, pd.get_dummies)

  • Numeric data must be scaled (otherwise your PCA will be misleading)

The Netflix's Recommender System through a PCA

Step 1: Collecting the data

These files contain metadata for all 45,000 movies listed in the Full MovieLens Dataset. The dataset consists of movies released on or before July 2017.

Data points include cast, crew, plot keywords, budget, revenue, posters, release dates, languages, production companies, countries, TMDB vote counts and vote averages.

import qmib
import numpy as np
import pandas as pd

movies = qmib.load("movies")

# Keep the numeric columns and drop rows with holes
df = movies[["budget", "popularity", "revenue", "runtime",
             "vote_average", "vote_count"]].dropna()
df.head()
budget popularity revenue runtime vote_average vote_count
0 30000000.0 21.946943 373554033.0 81.0 7.7 5415.0
1 65000000.0 17.015539 262797249.0 104.0 6.9 2413.0
2 0.0 11.712900 0.0 101.0 6.5 92.0
3 16000000.0 3.859495 81452156.0 127.0 6.1 34.0
4 0.0 8.387519 76578911.0 106.0 5.7 173.0

Step 2: Computing the principal components

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

# PCA is scale-dependent: standardise first, always.
Z = StandardScaler().fit_transform(df)
res_pca = PCA().fit(Z)

# Eigenvalues and the variance each component explains
eig = pd.DataFrame({
    "eigenvalue": res_pca.explained_variance_,
    "pct_variance": res_pca.explained_variance_ratio_ * 100,
    "cumulative_pct": np.cumsum(res_pca.explained_variance_ratio_) * 100,
}, index=[f"PC{i+1}" for i in range(len(res_pca.explained_variance_))])
print(eig.round(2))
     eigenvalue  pct_variance  cumulative_pct
PC1        2.96         49.41           49.41
PC2        1.11         18.57           67.97
PC3        0.85         14.13           82.10
PC4        0.60          9.96           92.06
PC5        0.31          5.20           97.26
PC6        0.16          2.74          100.00
import matplotlib.pyplot as plt

# Scree plot: where does the curve flatten?
plt.bar(range(1, len(res_pca.explained_variance_ratio_) + 1),
        res_pca.explained_variance_ratio_ * 100)
plt.xlabel("Component"); plt.ylabel("% of variance")
plt.show()

Step 3: Which variables contribute to the principal components?

# The loadings: how each variable relates to each component
loadings = pd.DataFrame(
    res_pca.components_.T * np.sqrt(res_pca.explained_variance_),
    index=df.columns, columns=[f"PC{i+1}" for i in range(res_pca.n_components_)])
print(loadings.round(3))
                PC1    PC2    PC3    PC4    PC5    PC6
budget        0.852 -0.133  0.074 -0.271  0.403 -0.121
popularity    0.716  0.064 -0.102  0.679  0.106  0.023
revenue       0.910 -0.160  0.013 -0.191 -0.106  0.314
runtime       0.222  0.704  0.674 -0.004 -0.028  0.006
vote_average  0.201  0.748 -0.614 -0.151  0.017  0.007
vote_count    0.899 -0.104 -0.032 -0.056 -0.356 -0.224
# Contribution of each variable to each component, in percent
contrib = (loadings ** 2)
contrib = contrib / contrib.sum() * 100
print(contrib.round(2))
                PC1    PC2    PC3    PC4    PC5    PC6
budget        24.50   1.59   0.65  12.33  51.95   8.98
popularity    17.29   0.37   1.22  77.22   3.58   0.32
revenue       27.93   2.30   0.02   6.10   3.58  60.06
runtime        1.66  44.54  53.52   0.00   0.25   0.02
vote_average   1.36  50.24  44.45   3.83   0.09   0.03
vote_count    27.27   0.97   0.12   0.52  40.54  30.59

Step 3 (cont.): The variables, drawn

# Variable plot: arrows on the unit circle
fig, ax = plt.subplots(figsize=(5.5, 5.5))
for name, (x, y) in loadings[["PC1", "PC2"]].iterrows():
    ax.arrow(0, 0, x, y, head_width=.02, color="#356B8C")
    ax.text(x * 1.08, y * 1.08, name, fontsize=9)
ax.add_patch(plt.Circle((0, 0), 1, fill=False, color="#B8C0CA"))
ax.set_xlim(-1.1, 1.1); ax.set_ylim(-1.1, 1.1)
ax.set_xlabel("PC1"); ax.set_ylabel("PC2")
plt.show()

contrib["PC1"].sort_values(ascending=False).head(10).plot.bar()
plt.ylabel("% contribution to PC1"); plt.show()

contrib["PC2"].sort_values(ascending=False).head(10).plot.bar()
plt.ylabel("% contribution to PC2"); plt.show()

Step 4: Where each film sits — the scores

# The scores: where each film sits in the new coordinate system
scores = pd.DataFrame(res_pca.transform(Z)[:, :5],
                      columns=[f"PC{i+1}" for i in range(5)])
scores.describe().round(3)
PC1 PC2 PC3 PC4 PC5
count 45203.000 45203.000 45203.000 45203.000 45203.000
mean -0.000 0.000 0.000 0.000 -0.000
std 1.722 1.055 0.921 0.773 0.559
min -1.193 -7.012 -10.076 -9.876 -13.220
25% -0.524 -0.386 -0.409 -0.260 -0.094
50% -0.379 0.139 -0.041 -0.098 -0.059
75% -0.099 0.572 0.338 0.236 0.002
max 54.969 21.086 21.366 72.961 10.687
scores.head().round(3)
PC1 PC2 PC3 PC4 PC5
0 10.716 -1.364 -1.501 -0.114 -6.208
1 7.305 -0.704 -0.341 -0.596 -0.758
2 0.450 0.592 -0.355 1.322 0.163
3 1.030 0.521 0.521 -0.412 0.371
4 0.902 0.152 0.092 0.613 -0.291

Step 4 (cont.): The films, drawn

# Individuals plotted on the first two components
plt.scatter(scores["PC1"], scores["PC2"], s=4, alpha=.3)
plt.xlabel("PC1"); plt.ylabel("PC2")
plt.show()

# Biplot: individuals as points, variables as arrows
fig, ax = plt.subplots(figsize=(6.5, 5.5))
ax.scatter(scores["PC1"], scores["PC2"], s=4, alpha=.2)
scale = 4
for name, (x, y) in loadings[["PC1", "PC2"]].iterrows():
    ax.arrow(0, 0, x * scale, y * scale, head_width=.08, color="#E3120B")
    ax.text(x * scale * 1.1, y * scale * 1.1, name, fontsize=9, color="#E3120B")
ax.set_xlabel("PC1"); ax.set_ylabel("PC2")
plt.show()

Factor analysis

What factor analysis assumes

Just like PCA, Factor Analysis is also a model that allows reducing information in a larger number of variables into a smaller number of variables. In Factor Analysis we call those latent variables.

Definition: Factor Analysis tries to find latent variables that make sense to us. We can rotate the solution until we find latent variables that have a clear interpretation and “make sense”.

Factor Analysis is based on a model called the common factor model. It starts from the principle that there a certain number of factors in a data set, and that each of the measured variables captures a part of one or more of those factors.

We could imagine having students who are overall good in languages, but bad in technical subjects. In this case, we could try to find a latent variable for language ability and a second latent variable for technical ability.

We now have latent variables that measure the general ability of a student for Language and Technical subjects. But it would still be possible that some students are great at languages overall, but that they are just bad at German.

This is why the Common Factor Model has specific factors: they measure the impact of one specific measured variable on the measured variable.

We could describe it as “Ability for learning German while taking into account the general ability for learning languages”.

The factor model

Where the PCA model is more of a pragmatic approach, in Factor Analysis we are hypothesizing that latent variables exist.

In a case with two latent variables, we can compute our original variables \(X\) by attributing one part of its variation to our first common latent variable (let’s call it \(k_1\)), part to the second common latent variable (\(k_2\)), and part to a specific factor (specific to this variable; called \(d\)).

In a case with 4 original variables, the Factor Analysis model would be as follows:

$$X_1 = c_{11} * k_1 + c_{12} * k_2 + d_1$$

$$X_2 = c_{21} * k_1 + c_{22} * k_2 + d_2$$

$$X_3 = c_{31} * k_1 + c_{32} * k_2 + d_3$$

$$X_4 = c_{41} * k_1 + c_{42} * k_2 + d_4$$

\(c\) are each coefficient of the coefficient matrix, which are the values that we need to estimate.

FAMD: the wine data

Fitting FAMD

# FAMD — factor analysis of mixed (numeric and categorical) data.
import prince

wine = qmib.load("redwines")

# FAMD wants the mix: the chemistry as numbers, the verdict as a category.
# `good` is dropped once it has become a label, or the same information enters
# twice and the first dimension is simply that duplication.
wine["quality_label"] = np.where(wine["good"] == 1, "good", "ordinary")
wine = wine.drop(columns=["good"])

famd = prince.FAMD(n_components=5, random_state=0).fit(wine)
# Eigenvalues of the FAMD dimensions
print(famd.eigenvalues_summary)
          eigenvalue % of variance % of variance (cumulative)
component                                                    
0              2.155        23.95%                     23.95%
1              1.861        20.68%                     44.63%
2              1.218        13.53%                     58.15%
3              1.140        12.67%                     70.82%
4              0.857         9.52%                     80.34%

Variable coordinates

# Variable coordinates: how much of each variable each dimension carries.
# `column_coordinates_` is a fitted attribute, not a method — no argument.
print(famd.column_coordinates_.round(3).to_string())
component                 0      1      2      3      4
variable                                               
fixed_acidity         0.371  0.267  0.122  0.068  0.002
volatile_acidity      0.002  0.454  0.001  0.000  0.308
residual_sugar        0.137  0.012  0.101  0.495  0.130
chlorides             0.217  0.003  0.490  0.002  0.106
total_sulfur_dioxide  0.058  0.086  0.030  0.480  0.229
density               0.819  0.011  0.081  0.006  0.002
sulphates             0.068  0.308  0.387  0.006  0.002
alcohol               0.397  0.256  0.004  0.050  0.057
quality_label         0.087  0.464  0.004  0.033  0.020
# Just the quantitative variables — the chemistry
num_cols = [c for c in wine.columns if c != "quality_label"]
print(famd.column_coordinates_.loc[num_cols].round(3).to_string())
component                 0      1      2      3      4
variable                                               
fixed_acidity         0.371  0.267  0.122  0.068  0.002
volatile_acidity      0.002  0.454  0.001  0.000  0.308
residual_sugar        0.137  0.012  0.101  0.495  0.130
chlorides             0.217  0.003  0.490  0.002  0.106
total_sulfur_dioxide  0.058  0.086  0.030  0.480  0.229
density               0.819  0.011  0.081  0.006  0.002
sulphates             0.068  0.308  0.387  0.006  0.002
alcohol               0.397  0.256  0.004  0.050  0.057
# Just the qualitative variable — the verdict. It takes a dimension largely of
# its own, which is the finding: the label is not reducible to the chemistry.
print(famd.column_coordinates_.drop(index=num_cols).round(3).to_string())
component          0      1      2      3     4
variable                                       
quality_label  0.087  0.464  0.004  0.033  0.02

The wines, drawn

# Individuals, coloured by the verdict
coords = famd.row_coordinates(wine)
plt.figure(figsize=(7, 4))
plt.scatter(coords[0], coords[1], s=6, alpha=.35,
            c=(wine["quality_label"] == "good").astype(int), cmap="RdYlBu_r")
plt.xlabel("Dim 1"); plt.ylabel("Dim 2")
plt.show()

Factor Analysis for Questionnaires

The cars questionnaire

# Exploratory factor analysis. statsmodels does this, so there is nothing to
# install: factor_analyzer, which this deck's first translation reached for,
# calls a scikit-learn function renamed in 1.7 and no longer runs.
from statsmodels.multivariate.factor import Factor

data = qmib.load("efa").dropna()

# Factor analysis is scale-dependent for the same reason PCA is. Standardise,
# so the factors describe correlation and not which item has the big numbers.
z = (data - data.mean()) / data.std()
data.head()
price safety exterior_looks space_comfort technology ... color maintenance test_drive product_reviews testimonials
0 4 4 5 4 3 ... 2 4 2 4 3
1 3 5 3 3 4 ... 4 3 2 2 2
2 4 4 3 4 5 ... 4 5 4 4 3
3 4 4 4 3 3 ... 4 4 2 5 3
4 5 5 4 4 5 ... 5 5 5 5 2

5 rows × 14 columns

# How many factors? The scree plot is the eigenvalues of the correlation
# matrix, largest first — no model has to be fitted to draw it.
ev = np.linalg.eigvalsh(np.corrcoef(z.to_numpy(), rowvar=False))[::-1]
plt.figure(figsize=(7, 2.6))
plt.plot(range(1, len(ev) + 1), ev, "o-")
plt.axhline(1, color="red", ls="--")   # Kaiser: keep eigenvalues above 1
plt.xlabel("Factor"); plt.ylabel("Eigenvalue")
plt.tight_layout(); plt.show()
print(f"Kaiser keeps {(ev > 1).sum()} factors; the elbow is nearer 3.")
print("The two rules disagree, which is why the number is argued and not read off.")

Kaiser keeps 5 factors; the elbow is nearer 3.
The two rules disagree, which is why the number is argued and not read off.

Three factors or four?

def efa(k):
    """Fit k factors by maximum likelihood, then rotate for readability.

    Rotation changes the loadings without changing the fit — the same solution
    in different coordinates. That is why a factor's *name* is never a result.
    """
    res = Factor(z.to_numpy(), n_factor=k, method="ml").fit()
    res.rotate("varimax")
    # statsmodels returns loadings as complex numbers with a zero imaginary
    # part; np.real drops it rather than printing `0.37+0.00j` on a slide.
    return pd.DataFrame(np.real(res.loadings), index=data.columns,
                        columns=[f"F{i + 1}" for i in range(k)])

loadings3 = efa(3)
print(loadings3.round(3).to_string())
                        F1     F2     F3
price               -0.055 -0.175 -0.495
safety               0.066 -0.221  0.225
exterior_looks       0.341 -0.186  0.014
space_comfort        0.073 -0.843  0.090
technology          -0.101 -0.346  0.008
after_sales_service -0.146 -0.443 -0.158
resale_value        -0.051  0.155 -0.817
fuel_type            0.024 -0.535 -0.024
fuel_efficiency     -0.407 -0.263 -0.417
color               -0.854  0.074 -0.083
maintenance         -0.373 -0.130 -0.526
test_drive          -0.148 -0.263 -0.064
product_reviews     -0.263 -0.315 -0.224
testimonials        -0.220 -0.183  0.254
loadings4 = efa(4)

# Blanking the small loadings is how a factor table is normally read, and it is
# a presentational choice with teeth: 0.3 is a convention with no theory behind
# it, and moving it moves which items appear to belong together.
print(loadings4.round(2).where(loadings4.abs() > 0.3, "").to_string())
                       F1    F2    F3    F4
price                           -0.51      
safety                      0.3            
exterior_looks       0.37                  
space_comfort              0.83            
technology                 0.36            
after_sales_service        0.45            
resale_value                    -0.75      
fuel_type                  0.57            
fuel_efficiency     -0.33       -0.48      
color               -0.94                  
maintenance         -0.34       -0.56      
test_drive                            -0.39
product_reviews                 -0.31 -0.39
testimonials                          -0.65

Which items load where

# The same table read as a list: each item under the factor it loads on.
# An item appearing twice is cross-loading; an item appearing nowhere is not
# measured by this solution at all. Both are findings, not untidiness.
for col in loadings4.columns:
    items = loadings4.index[loadings4[col].abs() > 0.3].tolist()
    print(f"{col}: " + (", ".join(items) if items else "(nothing above 0.3)"))

placed = {i for c in loadings4.columns
          for i in loadings4.index[loadings4[c].abs() > 0.3]}
orphans = [i for i in loadings4.index if i not in placed]
print(f"\nnot placed by any factor: {', '.join(orphans) if orphans else 'none'}")
F1: exterior_looks, fuel_efficiency, color, maintenance
F2: safety, space_comfort, technology, after_sales_service, fuel_type
F3: price, resale_value, fuel_efficiency, maintenance, product_reviews
F4: test_drive, product_reviews, testimonials

not placed by any factor: none

Naming the factors

Factor Items that load on it Name
F1 color, exterior_looks Aesthetics
F2 space_comfort, fuel_type, after_sales_service, technology, safety Functional benefits
F3 resale_value, maintenance, price, fuel_efficiency Economic value
F4 testimonials, test_drive, product_reviews Credibility
  • The sign of a factor is arbitrary: F3 and F4 load negatively on all their items, and flipping the whole column changes nothing. What matters is which items move together.

  • Within F1, color and exterior_looks have opposite signs, so F1 is a contrast between the two rather than a sum.

  • fuel_efficiency, maintenance and product_reviews load on two factors. The name is a claim you defend, not something the software returns.

When to use PCA and when to use Factor Analysis?

  • PCA allows us to find the components that contain the maximum amount of information in fewer variables. This makes it a great tool for dimension reduction.

  • PCA on the other hand is used in cases where we want to retain the largest amount of variation in the smallest number of variables possible.

PCA is much used in data preparation for Machine Learning tasks, where we want to help the Machine Learning model by already “summarising” the data in an easier to digest form.

  • Then as a consequence, the mathematics behind the two methods are, while close to each other, not exactly the same. Although both methods use decomposition, they differ in the details.

  • This also causes the fact that Factor Analysis has an additional possibility for rotation of the final solution, while PCA does not. Different Applications

  • As Factor Analysis is more flexible for interpretation, due to the possibility of rotation of the solution, it is very valuable in studies for marketing and psychology.

PCA’s advantage is that it allows for dimension reduction while still keeping a maximum amount of information in a data set. This is often used to simplify exploratory analyses or to prepare data for Machine Learning pipelines.

After the course

One more thing, and it is optional. Nothing here is examined. It is the method you will need the first time a client hands you forty indicators and asks for one index.

The QMIB Lab for this material — a knowledge check you can take at your own pace: warin.ca/qmib-labs

The lab still carries its original number, 07, from when this was session 7.