Basic Linear Algebra : Dimension reduction ด้วย PCA ตอน Python Implementation

กลุ่มโรคผิวหนัง Erythemato Squamous Disease (ESD) มี 6 ชนิด แต่ละชนิดมี clinical และ histopathology ที่คล้ายคลึงกัน การทำ differential diagnosis หรือการวินิจฉัยว่าผู้ป่วยเป็น ESD ชนิดใด ต้องอาศัยลักษณะทางคลินิคและข้อมูลอื่น (dimensions) มากถึง 34 อย่างประกอบกัน (อิงจาก dataset) ตัวอย่างนี้จะนำเอา dataset จาก https://archive.ics.uci.edu/ml/datasets/dermatology มาทำการลด dimension ลงด้วยวิธีการ PCA

1. import libraries

  
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.neighbors import KNeighborsClassifier

2.Pre-process

 
# 2.1 load data from website 
df = pd.read_csv("https://archive.ics.uci.edu/ml/machine-learning-databases/dermatology/dermatology.data")

# reshuffle data
df = df.sample(frac = 1)

# preview
df.tail()

 
# 2.2 -- copy data
X_derm = df.iloc[:,:-1].values
derm_diag = df.iloc[:,-1].values
derm_diag = derm_diag - 1

การสำเนาข้อมูลมาไว้ใน X_derm เพื่อให้ทำงาน matrix operation สะดวกขึ้น สำหรับ derm_diag ใช้เก็บ class label เอา 1 ไปลบออกจากค่าเดิมเพื่อปรับ label ให้อยู่ระหว่าง [0,5] แทนค่าเดิมที่อยู่ในช่วง [1,6]

   
# 2.3 standardize data

X_std = np.empty_like(X_derm).astype(np.float64)

for i in range(X_derm.shape[1]):
    mean = X_derm[:,i].mean()
    std = X_derm[:,i].std()
    X_std[:,i] = (X_derm[:,i] - mean)/std

การทำ standardize data คือการทำให้ข้อมูลในแต่ละ dimension มี mean = 0 และ standard deviation = 1 \[ X_{\text{std}} = \frac{X - \bar{X}}{\sigma_x}\] ขั้นตอนมีความสำคัญที่ต้องทำก่อนที่จะไปทำในขั้นตอน(ดูเนื้อหาในตอน Principal Component Analysis [2])

 
# 2.4 split data
nd = int(0.8*X_std.shape[0])

X_train = X_std[:nd,:]
X_test = X_std[nd:,:]

y_train = derm_diag[:nd]
y_test = derm_diag[nd:]

แยกข้อมูลออกเป็นสองส่วนด้วยสัดส่วน 8:2 ตัวแปร X_train และ y_train เก็บ data และ label ที่ใช้สำหรับสอนตัวแบบ X_testc และ y_test เก็บ data และ label สำหรับทดสอบตัวแบบ

3. Process PCA

   
# 3.1 compute covariance matrix, eigen vectors and eigen values

def eigen(A):
    # compute and sort eigen value descendingly     
    eval, evec = np.linalg.eig(A)
    idx = np.argsort(-eval)
    return (eval[idx], evec[:,idx])

cov_mat = (X_std.T @ X_std)/(X_std.shape[0]-1)
eig_val,eig_vec = eigen(cov_mat)

ขั้นตอนนี้คำนวณหา covariance matrix ของ standardized data eigen value(ดูการคำนวณใน Covariance Matrix and Correlation Matrix [3]) แล้วนำ covariance matrix ที่ได้ไปคำนวณหา eigen vectors และ eigen values ต่อ (ดูขั้นตอนการคำนวณจาก Eigen decomposition of matrix [4])

   .
# 3.2 visualize covariance matrix with heatmap correlation

plt.figure(figsize=(15, 10))
sns.heatmap(cov_mat,center= 0,vmin=-1, vmax=1,cmap="coolwarm")

heat map correlatiion บอกเราว่าระหว่าง dimension ทั้ง 34 dimension มีบางคู่ที่มีความสัมพันธ์กันค่อนข้างสูง สังเกตุจากสีจะเข้ม (ฟ้าเข้มหรือแดงเข้ม) ก็เป็นไปได้ว่าจะสามารถหา Principal component ที่ดีได้

   
# 3.3 compute explained variance ratio and visualize with scree plot

def explained_var(eig_val):
    sum_eig_val = np.sum(eig_val)
    exp_var = [val/sum_eig_val for val in eig_val]
    cs_var = np.cumsum(exp_var)
    return (exp_var,cs_var)
    
explained_var,cumsum_var = explained_var(eig_val)    

# -- scree plot
plt.figure(figsize=(15, 8))

plt.bar(range(cumsum_var.shape[0]), 
              explained_var, 
              alpha=0.8, 
              align='center',
              label='individual explained variance')
plt.step(range(cumsum_var.shape[0]), 
               cumsum_var, 
               where='mid',
               c='r',
              label='cumulative explained variance')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal components')
plt.legend(loc='best')
plt.grid()
plt.tight_layout()
plt.show()

scree plot บอกเราเกี่ยวกับความสัมพันธ์ระหว่างจำนวน principal component กับปริมาณ information ที่จะได้ ถ้าเราเลือกจำนวน principal component จำนวน 11 component ปริมาณ information ที่จะได้คือประมาณ 0.8 หรือ 80% เป็นต้น

   
# 3.4 create vector space for principal component 

def PCA(eig_vec,c_num=2):
    # principal components
    _pca = eig_vec[0].reshape((eig_vec.shape[0],1))
    for i in range(1,c_num):
        _pca = np.hstack((_pca,eig_vec[i].reshape((eig_vec.shape[0],1)))) 
    return _pca

cn = 16
pca = PCA(eig_vec,cn) # create vector space
X_pca = X_train @ pca # project data onto new vector space

ในตัวอย่างเลือกใช้จำนวน principal component เป็น 16 component เพื่อสร้าง vector space

4. Test

การทดสอบ principal component ที่ได้มาจะยังคงรักษา information ไว้ได้จริงหรือไม่ ในที่นี้จะใช้ K-clustering algorithm [5] ที่มีอยู่ใน sklearn.neighbors package มาใช้ ทดสอบด้วยข้อมูลที่แยกออกมาก่อนหน้าและไม่ได้นำไปรวมไว้ในการคำนวณหา principal component ตั้งแต่แรก

   
# 4.1 create and train KNeighborsClassifier

# there are 6 classes of ESD
model = KNeighborsClassifier(n_neighbors=6)

# train model
model.fit(X_pca, y_train)

# test prediction
pred = model.predict(X_pca)

# preview
acc = 0
for p,e  in zip(pred,y_train):
    if p==e : 
        acc += 1
print("Accuracy = {:.2f} %".format(100* acc/derm_diag_train.shape[0]))

Accuracy = 97.20 %

สังเกตุว่าค่าความถูกต้องของ model จะไม่ใช่ 100 % นั่นเพราะเป็นธรรมชาติของ algorithm ที่ใช้ [7] จำนวน principal component หรือ dimension ใน vector space ใหม่ก็มีส่วนด้วย ถ้าใช้จำนวน principal component มากความถูกต้องก็จะมากตาม

 
# 4.2 test pca with test-data

pred_test = model.predict(X_test @ pca)

# simple check
acc = 0
for p,e  in zip(pred_test, y_test):
    if p==e : 
        acc += 1
print("# Components : {}, Accuracy = {:.2f} %".format(cn,100* acc/derm_test.shape[0]))

# Components : 16, Accuracy = 95.83 %

ค่าความถูกต้องสูงใช้ได้ทีเดียว แสดงว่า principal component ที่ได้มาสามารถเก็บ information ได้ใกล้เคียงกับข้อมูลต้นฉบับ

ทดลองเปลี่ยนจำนวน principal component แล้วดูค่าความถูกต้องในการแยกกลุ่มของ ESD ได้ผลตามตาราง ดูเหมือนว่าความถูกต้องจะหยุดอยู่ที่ 95.83 %

# Principal component	Accuracy (%)
5	80.56
10	88.89
15	95.83
20	95.83
25	95.83

PCA ช่วยทำ visualizing

ประโยชน์ของการทำ data dimension reduction มีหลายด้าน ที่เห็นได้ง่ายที่สุดก็น่าจะเป็นเรื่อง data visualizing ด้วยความสามารถของสมองมนุษย์สามารถสร้างและตีความภาพได้เพียง 3 dimensions แต่ด้วยข้อมูลของ ESD มีถึง 34 dimensions ซึ่งเป็นไม่ได้เลยที่จะนำเอาทุก dimension มาสร้าง scatter plot พร้อมกัน อาจทำได้โดยการจัดมา plot ครั้งละ 3 dimensions ซึ่งจะเป็นภาระพอสมควร แต่การ project ข้อมูล 34 dimensions ไปยัง vector space ที่มี 3 dimensions ผ่านกระบวนการ PCA จะทำให้สร้าง scatter plot เป็นไปได้ ดังรูปที่ 1