Lab 3 - Dimension Reduction via PCA¶
This lab covers the use of principal component analysis (PCA) as a dimension reduction approach in sklearn.
Directions: Please read through the contents of this lab with your partner and try the examples. After you're both confident that you understand a topic you should attempt the associated exercise and record your answer in your own Jupyter notebook that you will submit for credit. The notebook you submit should only contain answers to the lab's exercises (so you should remove any code you ran for the examples, or use a separate notebook to test out the examples).
We'll use our standard set of libraries:
## Libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import sklearn
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
We will look at two types of data in this lab:
- A collection of synthetic data sets with increasing numbers of non-informative features.
- A random sample of $n=6000$ grayscale images of handwritten digits from the MNIST database.
The synthetic data use a binary outcome with an expected 50-50 distribution of each class. A specified number of features (two in our examples) are positively associated with the outcome $Y=1$, while other features are non-informative, as their values are independent of the outcome.
## Function used to make the synthetic data
def make_synth_data(n_sig, n_noise):
n = 400
# Binary outcome (50-50 class distribution)
np.random.seed(7)
y = np.random.binomial(1, 0.5, n)
# Create the signal features
X_signal = np.zeros((n, n_sig))
X_signal[y == 0] = np.random.normal(-1, 1, size=(sum(y == 0), n_sig))
X_signal[y == 1] = np.random.normal(+1, 1, size=(sum(y == 1), n_sig))
# Noise features
X_noise = np.random.normal(0, 1, size=(n, n_noise))
# Combine signal and noise features
X = np.hstack([X_signal, X_noise])
return X, y
## Make 3 versions of the synthetic data w/ increasing n_noise
X2, y2 = make_synth_data(2, 2)
X48, y48 = make_synth_data(2, 48)
X198, y198 = make_synth_data(2, 198)
The MNIST data were originally recorded as 28 by 28 px images with each feature being the grayscale color intensity of a pixel. They are given in "flattened" form, or as a 6000 by 784 matrix (with an extra column for the class label). To understand what the data look like we can reshape and plot a few examples:
### Read flattened version of MNIST
mnist = pd.read_csv("https://remiller1450.github.io/data/mnist_small.csv")
### Separate the target variable (label column)
label = mnist['label']
mnist = mnist.drop(['label'], axis=1)
### Convert to numpy array and reshape to 28 by 28
mnist_unflattened = mnist.to_numpy()
mnist_unflattened = mnist_unflattened.reshape(6000,28,28)
## Plot some examples
import matplotlib.cm as cm
fig, axs = plt.subplots(ncols=7)
for i in range(7):
axs[i].imshow(mnist_unflattened[i], cmap=cm.Greys)
axs[i].title.set_text(f'label={label[i]}')
plt.show()
Part 1 - When to Consider PCA¶
Dimension reduction via PCA has the greatest potential to improve model performance when:
- The modeling approach is sensitive to high-dimensional data (for example, KNN) and we suspect that many features contain little useful signal or are primarily noise.
- Many features are highly correlated, suggesting that the data may contain a smaller number of underlying latent factors. In this case, a few principal components can often capture most of the information contained in the original features.
However, it is also important to note that PCA does not know which variables are signal and which are noise, so we can never be fully certain (without an empirical investigation) that dimension reduction will yield better performance when working with complex, real-world data.
Assessing the dimensionality of our data is straightforward, but understanding the correlation structure of the data can be trickier. Below we visualize the correlation matrix of the features in the MNIST data set:
## Separate training and testing sets
from sklearn.model_selection import train_test_split
mnist_train, mnist_test, label_train, label_test = train_test_split(mnist,label, test_size=0.2,random_state=123)
## Correlation plot
plt.matshow(mnist_train.corr())
plt.gcf().set_size_inches(12, 12)
plt.show()
There are a lot of features in the MNIST data, so interpretting this plot can be challenging, but you should notice that bands of adjacent pixels do tend to have high positive correlations, suggesting PCA might be a useful pre-processing step.
Question #1:
- Part A: Split the synthetic data with 198 noise features (
X198andy198) into a training set containing $n=300$ observations and a test set containing $n=100$ observations usingrandom_state=7. - Part B: Using the training data from Part A, plot the correlation matrix and briefly describe what you see. Hint: the
.corr()method only exists forpandasdataFrames. - Part C: Considering all that you currently know about the synthetic data used in Parts A and B, do you believe that PCA might improve the performance of KNN classifier? Briefly explain your reasoning.
Part 2 - PCA in sklearn¶
The implementation of PCA in sklearn is similar to the scaling and transformation pre-processing operations covered in our previous lab. More specifically, we can use the fit() method to fit the transformation (determine the principal component loadings) and the transform() method to apply the transformation (convert the original data values to coordinates on principal component axes).
We demonstrate this below the MNIST data, projecting the original coordinates onto 10 retained principal components:
## Fit the PCA
from sklearn.decomposition import PCA
pca10_reducer = PCA(n_components=10).fit(mnist_train)
## Apply the PCA transformation to reduce the dim of mnist_train
mnist_train_pca10 = pca10_reducer.transform(mnist_train)
# Examine the transformed data
print(mnist_train_pca10.shape)
(4800, 10)
Before we move on, you should note that PCA is sensitive to the scale of the input features, so most applications require a scaling step before dimension reduction via PCA. We did not need to rescale these data because each pixel is already measured on the same color intensity scale.
Rather than blindly guessing that 10 principal components provide a sufficient representation of the interesting axes of variation in our data, we should be more deliberate and consider variance explained by component. This can be done using the explained_variance_ratio_ attribute of a fitted PCA object:
## Full PCA
full_pca_reducer = PCA(n_components=784).fit(mnist_train)
full_pca_reducer.explained_variance_ratio_[0:20]
array([0.09986798, 0.07159082, 0.06129261, 0.0533461 , 0.04809149,
0.04324373, 0.0332702 , 0.0279454 , 0.02751206, 0.02404282,
0.02149484, 0.02055013, 0.01738087, 0.01666684, 0.01588959,
0.01436303, 0.01299478, 0.01251487, 0.01172244, 0.01149775])
This example shows the variance explained by the first 20 components. While a component explaining ~1% of the variance in the data might seem small, for a data set with 784 there might still be predictive information in this feature.
A common practice is retain however many components are need to large fraction of the total variance in the data (such as 80%, 90%, 99%, etc.), which can be done automatically by providing a decimal value to the n_components argument of PCA():
## PCA to keep 95% of the total variation
v95_pca_reducer = PCA(n_components=0.95).fit(mnist_train)
v95_pca_reducer.explained_variance_ratio_.shape
(149,)
Here we see that 149 components are needed to retain 95% of the variation in the original 784, which aligns with what we saw in the correlation matrix earlier.
Question #2:
- Part A: Create a scree plot showing the component number on the x-axis and the variance explained on the y-axis for all principal components in the MNIST example. Do you notice an "elbow", or a location where the variance explained by additional components begins to level off? Approximately how many principal components appear to capture most of the important variation in the data according to your analysis of the scree plot?
- Part B: Perform PCA on the training set for synthetic data with 198 noise features (
X198andy198). Retain all components and create a scree plot showing the component number of the x-axis and the variance explained on the y-axis. Considering what you know about the correlation matrix of each data set (the synthetic data and the MNIST data), why might their scree plots exhibit different shapes? - Part C: Using the data from Part B, find the number of components needed to retain 90% of the variance in the original data. Does this appear to correspond with an "elbow" on the scree plot? Briefly explain.
- Part D: Considering the sensitivity of PCA to the scale of the input variables, is it possible that we made a mistake by not re-scaling the synthetic data prior to using PCA? Consider how these data were generated when explaining your answer.
Part 3 - Dimension Reduction and Performance¶
We've now covered the essential functions needed to use PCA as a pre-processing tool. Let's now see if it can improve the performance of KNN and decision tree classifiers on the data we've been working with.
Question #3:
- Part A: For each sythentic data set (X2, X48, and X198):
- Create a 75%-25% train-test split using
random_state = 7 - Fit a PCA transformation to the data, keeping enough principal components to contain 80% of the total variation in the data.
- Report the number of retained principal components for each data set.
- Create a 75%-25% train-test split using
- Part B: Using the training and testing sets from Part A, fit a KNN classifier with $k=10$ to both the original data and the PCA-transformed data for each synthetic data set. Record the accuracy of both approaches on the test set.
- Part C: Does PCA seem to improve KNN classifier performance when the number of noise variables is large? Briefly describe any relationships you see between the number of noise variables, the dimensionality of the data given to KNN, and classification performance.
Question #4:
- Part A: Using the training-testing split performed on the MNIST dataset earlier in the lab, fit a KNN classifier with $k=10$ using the original 784 pixel features. Report the accuracy of this model on the test set.
- Part B: Now apply a PCA-transformation to the MNIST training data to retain 75% of the total variation in the data. How many principal components are retained?
- Part C: Now fit a KNN classifier with $k=10$ to the PCA-transformed data from Part B. Report the accuracy of this model on the test set, being sure to first apply the PCA-transformation you fit in Part B to the testing data.
- Part D: Briefly compare and contrast your findings in this question with those in Question 3 involving the synthetic data. Did PCA seem to work better for one of these scenarios? How might the correlation structures of these datasets help explain the results.
Part 4 - PCA for Inference and Feature Engineering¶
Part of being a skillful machine learning practitioner is understanding the contextual nuances of your data and knowing how to leverage them when making modeling choices. In the context of PCA, this might entail recognizing when the principal components themselves provide useful information about underlying structures in the data.
For example, in our lecture we saw an example where responses to personality survey items were highly correlated because each were measuring a common latent trait such as extroversion, agreeableness, etc. In this setting we might only only consider PCA as a pre-processing tool, but also as a way to explore our data and create new features.
To facilitate this, we need to inspect the most influential loadings in the most important principal components. The example below uses the Big 5 personality data set from our lecture, printing the ten most influential features in the first principal component:
## Big 5 data
bf = pd.read_csv("https://remiller1450.github.io/data/big5data.csv", sep='\t')
## Split the personality questions from the demographics
bf_q = bf.drop(['race','age','engnat','gender','hand','source','country'], axis=1)
bf_demo = bf[['race','age','engnat','gender','hand','source','country']]
## Loadings for a particular principal component (sorted by abs magnitude)
pca_bfq = PCA().fit(bf_q)
PC1_loadings = pd.DataFrame({'Question': bf_q.columns, 'PC1': pca_bfq.components_[0]})
print(PC1_loadings.sort_values(by='PC1', ascending=False, key=abs).head(10))
Question PC1 6 E7 0.263503 2 E3 0.252589 4 E5 0.240867 9 E10 -0.223513 19 N10 -0.222202 3 E4 -0.211569 17 N8 -0.205233 18 N9 -0.197345 5 E6 -0.196707 1 E2 -0.195601
The survey items for the four most influential features in PC1 are:
- E7 - I talk to a lot of different people at parties.
- E3 - I feel comfortable around people.
- E5 - I start conversations.
- E10 - I am quiet around strangers.
We might notice that these seem to measuring extroversion, so perhaps engineering a feature to capture a person's level of extroversion (such as a total extroversion score) might yield better performance than giving each feature to model individually, as this allows for lower dimensionality and cuts down on noise as each individual survey item might be subject to random measurement error, but a composite score will be more stable by averaging out possible misunderstanding or inconsistencies in a subject's responses.
Question #5: In this question you'll use data from an online questionnaire aimed at capturing "dark triad" personality traits, which are machiavellianism (a manipulative attitude), narcissism (excessive self-love), and psychopathy (lack of empathy). For more information you can view the data source at this link: http://openpsychometrics.org/tests/SD3/
- Part A: Separate 80% of the data into a training set using
random_state=7, then visualize the correlation matrix of the Dark Triad questionaire items (found in the image below). Briefly describe what you seen in this visualization. - Part B: Perform PCA on the Dark Triad questionaire items and print the nine survey items with the largest loadings for the first three principal components.
- Part C: Using your results from Part A, assign each of the dark triad traits, Machiavellianism, Narcissism, and Psychopathy, to one of the first three principal components.
- Part D: Create a scree plot displaying the variance explained by the first ten principal components. Based upon this plot, does it seem like three components provides a reasonable representation of the data?
- Part E: Suppose you are building a machine learning model that can use these questionaire responses. Explain why replacing groups of highly correlated survey items with a smaller number of engineered features (such as principal component scores or manually-created trait scores) might be beneficial.
## Dark Triad data
dt = pd.read_csv("https://remiller1450.github.io/data/dark_triad.csv", sep="\t")
## Screenshot of the Dark Triad survey items
from IPython.display import Image
Image("C:\\Users\\millerry\\OneDrive - Grinnell College\\Documents\\STA-395_Intro_ML\\Spring24\\Labs\\dark.PNG")
