---
jupyter:
  jupytext:
    notebook_metadata_filter: all,-language_info
    split_at_heading: true
    text_representation:
      extension: .Rmd
      format_name: rmarkdown
      format_version: '1.2'
      jupytext_version: 1.16.0
  kernelspec:
    display_name: Python 3 (ipykernel)
    language: python
    name: python3
---

# Dimension reduction with Principal Component Analysis

See the [PCA introduction](pca_introduction) for background.

Start by loading the libraries we need, and doing some configuration:

```{python}
import numpy as np
import numpy.linalg as npl
# Display array values to 6 digits of precision
np.set_printoptions(precision=6, suppress=True)

import pandas as pd
pd.set_option('mode.copy_on_write', True)

import matplotlib.pyplot as plt

from sklearn.decomposition import PCA
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
```

As for the PCA page above, we will use the English Premier League (EPL) table
data for this example. See [the EPL dataset
page](https://github.com/odsti/datasets/tree/main/premier_league) for more on
this dataset.

```{python}
df = pd.read_csv('data/premier_league_2021.csv')
df.head()
```

We want to predict goal difference from spending on the various position types.


We have spending data on `keeper` `defense`, `midfield` and `forward`.

```{python}
y = df['goal_difference']
# Divide spending by 10,000 GBP to make numbers easier to read.
X = df[['keeper', 'defense', 'midfield', 'forward']] / 10_000
X.head()
```

```{python}
basic_reg = LinearRegression().fit(X, y)
basic_reg.coef_
```

How good is the model fit?

```{python}
fitted = basic_reg.predict(X)
fitted
```

```{python}
full_r2 = r2_score(y, fitted)
full_r2
```

It turns out the variables (features) are highly correlated, one with another.

```{python}
X.corr()
```

We want to reduce these four variables down to fewer variables, while still capturing the same information.

Let us first run PCA on the four columns:

```{python}
pca = PCA()
pca = pca.fit(X)

# Scikit-learn's components.
pca.components_
```

Review the scree-plot of the component weights:

```{python}
plt.plot(pca.singular_values_)
```

We can calculate the amount of each row explained by each component, by using `fit_transform` of the PCA object:

```{python}
comp_vals = pca.fit_transform(X)
comp_vals
```

If you have followed the [PCA introduction page](pca_introduction), you will know you can also get these scores by matrix multiplication:

```{python}
# Getting component values by using matrix multiplication.
demeaned_X = X - np.mean(X, axis=0)
demeaned_X.dot(pca.components_.T)
```

We can use these four columns instead of the original four columns as
regressors. Think of the component scores (values) as being a re-expression of
the original four columns, that keeps all the same information.  When we fit this model, the component scores will get their own coefficients, but because those columns contain the same information as the original columns, the resulting fitted values will be the same as for the original model (within the precision of the calculations).

```{python}
comp_reg = LinearRegression().fit(comp_vals, y)
comp_fitted = comp_reg.predict(comp_vals)
# The columns and parameters are different, but the fitted values
# are the same.
np.allclose(comp_fitted, fitted)
```

Now we've seen the scree-plot, we might conclude that only the first component
is important, and we can discard the others.  We then just include the first
component values a new model.

```{python}
comp1_vals = comp_vals[:, :1]  # A 2D array with one column
comp1_vals
```

```{python}
comp1_reg = LinearRegression().fit(comp1_vals, y)
comp1_reg.coef_
```

We don't have quite the same quality of model fit — as measured by $R^2$, for
example, but we have recovered most of the $R^2$ using this single regressor
instead of the original four.

```{python}
comp1_fitted = comp1_reg.predict(comp1_vals)
r2_score(y, comp1_fitted)
```

We could improve the fit a little by adding the second component:

```{python}
comp12_reg = LinearRegression().fit(comp_vals[:, :2], y)
r2_score(y, comp12_reg.predict(comp_vals[:, :2]))
```

And a little more with the third:

```{python}
comp123_reg = LinearRegression().fit(comp_vals[:, :3], y)
r2_score(y, comp123_reg.predict(comp_vals[:, :3]))
```

Until we get the full $R^2$ with all four components (which gives the same
underling model as using the original four regressors):

```{python}
comp_reg = LinearRegression().fit(comp_vals, y)
r2_score(y, comp_reg.predict(comp_vals))
```
