<jupyter_start><jupyter_text>Two flexible ways to visualize PCA - **Explained Variance Cumulative Plot:** This is simple but powerful. It immediately tells you how much of the variance in the data is explained by each component and how the combinations of the different components add up to the total variance.- **Principal Components Overlayed with the Original Data:** This one is my absolute favorite. You can see the progression of how each principal component brings in slightly more information, and in turn, it becomes hard to distinguish between the different components. This plot is a perfect companion to the Explained Variance Cumulative Plot.  Data PreparationThe dataset we’ll be using is the Wine Data Set from UC Irvine’s Machine Learning repository. The data set contains data about wine quality. This dataset is licensed under a Creative Commons Attribution 4.0 International (CC BY 4.0) license.<jupyter_code># Load libraries
import numpy as np
import pandas as pd
import warnings 
warnings.filterwarnings('ignore')

# Load libraries for machine learning
from sklearn import metrics
from sklearn.metrics import ConfusionMatrixDisplay
from sklearn.pipeline import Pipeline
from sklearn.pipeline import FeatureUnion
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.decomposition import IncrementalPCA
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression

# Load libraries for plotting
import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline
# Read the data
df = pd.read_csv('/content/wine.data')

# Add colunm names to the df
df.columns = ['Class','Alcohol','Malic acid','Ash','Alcalinity o ash','Magnesium','Total phenols',
              'Flavanoids','Nonflavanoid phenols','Proanthocyanins','Color intensity',
              'Hue','OD280/OD315 of diluted wines','Proline']
df.head()<jupyter_output><empty_output><jupyter_text>This particular dataset has 13 features and a target variable named Class we can use for classification. Each of the values is continuous, and therefore we can directly apply PCA to all variables. Normally we would use PCA on a much higher dimensional dataset, but this dataset will show the concepts.<jupyter_code>df.info()
# check for imbalanced data
df['Class'].value_counts()<jupyter_output><empty_output><jupyter_text>This dataset, unlike most, is pretty balanced. We'll leave it as is for this demonstration.<jupyter_code># Create X and y variable 
X = df.drop(columns = ['Class']).copy().values
y = df['Class'].copy().values
X.shape
y<jupyter_output><empty_output><jupyter_text>Find Explained Variance- The first step is to find the **explained variance** for each principal component. We calculate explained variance by first scaling our data with the *StandardScalar* and fitting a PCA model to the scaled data.- The *StandardScalar* is a simple transformation that normalizes the data to have a mean of 0 and 0 unit variance, or a standard deviation of 1.- You will also notice that the total number of *Principal Components equals the number of features in your dataset*. However, one important thing to note is that the Principal Components are not the dataset’s columns; rather, PCA is constructing new features that best explain the data.<jupyter_code>def get_variance(X,n):
  scaler = StandardScaler()
  pca = PCA(n_components = n)
  pca.fit(scaler.fit_transform(X))

  return pca.explained_variance_ratio_.cumsum()[-1:] # return the value except the last one

for i in range(1,14):
  print('Components:\t', i, '=\t', get_variance(X, i), '\tCumulative Variance')<jupyter_output>Components:	 1 =	 [0.35983071] 	Cumulative Variance
Components:	 2 =	 [0.5522435] 	Cumulative Variance
Components:	 3 =	 [0.6640381] 	Cumulative Variance
Components:	 4 =	 [0.73514919] 	Cumulative Variance
Components:	 5 =	 [0.80143663] 	Cumulative Variance
Components:	 6 =	 [0.8510403] 	Cumulative Variance
Components:	 7 =	 [0.89350044] 	Cumulative Variance
Components:	 8 =	 [0.92048035] 	Cumulative Variance
Components:	 9 =	 [0.94272497] 	Cumulative Variance
Components:	 10 =	 [0.96171025] 	Cumulative Variance
Components:	 11 =	 [0.9790518] 	Cumulative Variance
Components:	 12 =	 [0.99200387] 	Cumulative Variance
Components:	 13 =	 [1.] 	Cumulative Variance<jupyter_text>Additionally, each of the principal components is summed together, and the total of all components will equal 1. As you look through the list of cumulative variance, you will see that **the first component is the most important and the last component is the least important**; in other words, the first component contributes the most to the variance. As we include more and more of the components, the contribution amount per component starts to decrease. Let's visualize this by plotting the cumulative variance.  Plot the Threshold for Explained VarianceWhile the printout of the values is a good start, we can improve by plotting each of these values. Additionally, we will plot a line that represents `95%` of the explained variance. While there is no rule as to how much explained variance you need to include in a model, `95%` is a good threshold to start. Later, we will perform a grid search to find the optimal number of components to use in a model.<jupyter_code>scaler = StandardScaler()
data_rescaled = scaler.fit_transform(X)

pca = PCA().fit(data_rescaled)

plt.rcParams['figure.figsize'] = (8,5)

fig,ax = plt.subplots()
xi = np.arange(1, 14, step = 1) # 13 components in total
y = np.cumsum(pca.explained_variance_ratio_)

plt.ylim(0.0,1.1)
plt.plot(xi, y, marker = 'o', linestyle = '-', color = 'black')

plt.xlabel('Number of Components')
plt.xticks(np.arange(1,14,step = 1))
plt.ylabel('Cumulative variance (%)')
plt.title('The number of components needed to explain variance')

# Add a threshold line at 95%
plt.axhline(y=0.95, color = 'grey',linestyle ='--')
plt.text(1.1, 1, '95% cut-off threshold', color = 'black', fontsize = 16)

ax.grid(axis = 'x')
plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>After plotting the cumulative explained variance, we can see how the curve flattens slightly around `6` or `7` components. And where our line is drawn for 95%, the total explained variance is at approximately `9` components. And that is explained variance visualized! Let's move on to plotting the actual components.  Plotting Each Component vs. Original DataThe next step is to visualize the data resulting from each component: rather than the explained variance. We will use the **inverse_transform** method of the PCA model, and this will take each of the components and transform them back into the original data scale. We will plot the entire range of components from 1-13 and overlay them with the original data.<jupyter_code># Function to perform PCA with n-components
def transform_pca(X, n):
  pca = PCA(n_components=n)
  pca.fit(X)
  X_new = pca.inverse_transform(pca.transform(X))
  return X_new
rows = 4
cols = 4
comps = 1

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

fig, axes = plt.subplots(rows, 
                         cols, 
                         figsize=(12,12), 
                         sharex=True, 
                         sharey=True)


for row in range(rows):
    for col in range(cols):
        try:
            X_new = transform_pca(X_scaled, comps) # the orginal data
            ax = sns.scatterplot(x=X_scaled[:, 0], 
                                 y=X_scaled[:, 1], 
                                 ax=axes[row, col], 
                                 color='blue', 
                                 alpha=.3)
            ax = sns.scatterplot(x=X_new[:, 0],  # each PCA component
                                 y=X_new[:, 1], 
                                 ax=axes[row, col], 
                                 color='red',
                                 alpha = .5)
            ax.set_title(f'PCA Components: {comps}');

            comps += 1
        except:
            pass
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>With just one component displayed, it takes the form of a set of points projected on a line that goes is along the axis with the most variance in the original data. As more and more components are added, we can see how the data starts to resemble the original data even though there is less information being displayed. By the time we get into the upper values, the data starts to match the original data where finally, at all 13 components, it is the same as the original.  Comparing PCA AND Non-PCA Classification ModelsLet’s run three classifiers against our dataset, both with and without PCA applied, and see how they perform. We’ll compare them utilizing only the first two components ( n_components=2) to make it interesting. We'll compare **KNeighborsClassifier**, **RandomForestClassifier**, and **LogisticRegression** as our classifiers and see how they perform.As always, we’ll use a **pipeline** to combine our data preprocessing and PCA into a single step. We’ll scale our data with the same StandardScaler we used previously and then fit the data with a certain number of components. Since we'll use this function a few times in our workflow, we've set it up to be flexible.<jupyter_code>def create_pipe(clf, do_pca=False, n=2):

    scaler = StandardScaler()
    pca = PCA(n_components=n)

    if do_pca == True:
        combined_features = FeatureUnion([("scaler", scaler), 
                                          ("pca", pca)])
    else:
        combined_features = FeatureUnion([("scaler", scaler)])

    pipeline = Pipeline([("features", combined_features), 
                         ("clf", clf)])


    return pipeline<jupyter_output><empty_output><jupyter_text>Next is looping over the different classifiers and performing **Cross-Validation**. As mentioned, we’ll try three different classifiers and test them each with and without PCA applied.<jupyter_code>models = {'KNeighbors' : KNeighborsClassifier(),
          'RandomForest' : RandomForestClassifier(random_state=42),
          'LogisticReg' : LogisticRegression(random_state=42),
          }


def run_models(with_pca):
    for name, model, in models.items():
        clf = model
        pipeline = create_pipe(clf, do_pca = with_pca, n=2)
        scores = cross_val_score(pipeline, X, 
                                 y, 
                                 scoring='accuracy', 
                                 cv=3, n_jobs=1, 
                                 error_score='raise')
        print(name, ': Mean Accuracy: %.3f and Standard Deviation: \
             (%.3f)' % (np.mean(scores), np.std(scores)))

print(68 * '-')
print('Without PCA')
print(68 * '-')
run_models(False)
print(68 * '-')
print('With PCA')
print(68 * '-')
run_models(True)
print(68 * '-')<jupyter_output>--------------------------------------------------------------------
Without PCA
--------------------------------------------------------------------
KNeighbors : Mean Accuracy: 0.949 and Standard Deviation:              (0.024)
RandomForest : Mean Accuracy: 0.972 and Standard Deviation:              (0.021)
LogisticReg : Mean Accuracy: 0.972 and Standard Deviation:              (0.021)
--------------------------------------------------------------------
With PCA
--------------------------------------------------------------------
KNeighbors : Mean Accuracy: 0.678 and Standard Deviation:              (0.048)
RandomForest : Mean Accuracy: 0.955 and Standard Deviation:              (0.021)
LogisticReg : Mean Accuracy: 0.972 and Standard Deviation:              (0.029)
--------------------------------------------------------------------<jupyter_text>Find the Optimal Number of ComponentsFor this portion, we’ll use the train_test_split function to split the dataset into 70/30% partitions and then use the GridSearch feature to find optimal parameters. The most critical parameter we want to validate is the performance of various **PCA component numbers**.As we saw from the scatterplots above, when we approach higher numbers like 9, the dataset after transformation looks a lot like the original with 95% of the variation explained by the transformed dataset. We also saw how only 2 components performed well, so we can test the entire range from 1-13 and find the best setting.<jupyter_code># Make training and test sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.3, random_state = 53)<jupyter_output><empty_output><jupyter_text>Hyperparameter TuningAfter selecting the model based on cross-validation utilizing default parameters, we perform a **Grid Search** to fine-tune the model and select the best parameters. We'll loop over three things.<jupyter_code>def get_params(parameters, X, y, pipeline):

    grid = GridSearchCV(pipeline, 
                        parameters, 
                        scoring='accuracy', 
                        n_jobs=1, 
                        cv=3, 
                        error_score='raise')
    grid.fit(X, y)

    return grid
clf = LogisticRegression(random_state=41)
pipeline = create_pipe(clf, do_pca=True)

param_grid = dict(features__pca__n_components = list(range(2,14)),
                 clf__C = [0.1, 1.0, 10, 100],
                 clf__solver = ['liblinear', 'saga'],
                 clf__penalty = ['l2', 'l1'])

grid = get_params(param_grid, X_train, y_train, pipeline)

print("Best cross-validation accuracy: {:.3f}".format(grid.best_score_))
print("Test set score: {:.3f}".format(grid.score(X_test, y_test))) 
print("Best parameters: {}".format(grid.best_params_))<jupyter_output>Best cross-validation accuracy: 0.976
Test set score: 1.000
Best parameters: {'clf__C': 10, 'clf__penalty': 'l1', 'clf__solver': 'liblinear', 'features__pca__n_components': 6}<jupyter_text>- We can see that the best parameters are *C=10, penalty='l1', solver='liblinear' and n_components=2*. We can also see that the test set score is close to the best cross-validation score. We can also see that the cross-validation accuracy improved a little from 0.972 to 0.976.- As a side note, this is pretty typical of Hyperparameter Tuning. You don’t get massive improvements; it’s small, incremental improvements, and the biggest gains are through **feature engineering** and **model selection**.<jupyter_code><jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1055738605161999, 'fraction_numerical': 0.022377814387699067, 'mean_word_length': 3.2586962876351944, 'pattern_counts': {'":': 0, '<': 42, '<?xml version=': 0, '>': 42, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '312658', 'license': 'Apache-2.0', 'n_tokens_mistral': 4308, 'n_tokens_neox': 3963, 'n_words': 1585, 'path': 'PCA/Visualization/PCA_plotting.ipynb', 'repo_name': 'les1smore/ML-Notes', 'size': '14568'}