Pythonでワインデータを分析してみた

#目次
[1.目的](#1目的)
[2.実行環境](#2-実行環境)
[3.データ取得元](#3-データ取得元)
[4.世の中にある分析手法の例](#4-世の中にある分析手法の例)
[5.今回実施した分析およびその結果考察](#5-今回実施した分析およびその結果考察)
[6.まとめ](#6-まとめ)
[7.今後の課題](#7-今後の課題)
[8.参考資料](#8-参考資料)



# 1.目的
機械学習の様々な分析手法の特性を理解したうえで、ワインデータに対して適用し、結果(精度)の比較を行う。

# 2. 実行環境
Google Colaboratory

# 3. データ取得元
ネットで検索すると様々なワインデータが存在するが、今回はKaggleのデータを使用する。
https://www.kaggle.com/datasets/sgus1318/winedata?resource=download
以下の赤ワインと白ワインの成分(物理化学的性質)から、Quality(ワインの味(グレード:0から10のスコア))を予測する課題である。
・ワインの色(赤・白)
・fixed acidity(酒石酸濃度(g/dm3))
・volatile acidity(酢酸濃(g/dm3))
・citric acid(クエン酸濃度(g/dm3))
・residual sugar(残留糖分濃(g/dm3))
・chlorides(塩化ナトリウム濃度(g/dm3))
・free sulfur dioxide(遊離亜硫酸濃度(mg/dm3))
・total sulfur dioxide(総亜硫酸濃度(mg/dm3))
・density(密度(g/dm3))
・pH(pH)
・sulphates(硫酸カリウム濃度(g/dm3))
・alcohol(アルコール度数(% vol.))



# 4. 世の中にある分析手法の例
※ディープラーニングの手法は、別途考察
※今回のデータの特性を踏まえ、自然言語処理や時系列は、除外

・教師なし学習
データ構造の理解のため、以下の手法で、データを確認を行う
  クラスタリング
  主成分分析
  カーネル主成分分析

・分類
  ロジスティック回帰
  決定木
  ランダムフォレスト
  線形SVM
  非線形SVM
  k-NN
※ハイパーパラメーターの調整やパラメーターのランダムサーチも行う

・回帰(正則化)
  線形回帰
  ラッソ回帰
  リッジ回帰
  ElasticNet回帰



# 5. 今回実施した分析およびその結果考察

・まず、データの読み込み～データ構造の理解まで確認します。
<details><summary>データの読み込み～データ構造の理解まで</summary>

```python
import numpy as np 
import pandas as pd 

from sklearn import preprocessing
import matplotlib.pyplot as plt 

import statsmodels.api as sm
import statsmodels.formula.api as smf

from sklearn import datasets, linear_model
from sklearn.metrics import mean_squared_error

%matplotlib inline

plt.rc("font", size=14)
import seaborn as sns
sns.set(style="white") #white background style for seaborn plots
sns.set(style="whitegrid", color_codes=True)

# from ggplot import *




### Load wine quality data into Pandas
df_red = pd.read_csv("/content/drive/MyDrive/Colab Notebooks/winequality_red.csv")
df_white = pd.read_csv("/content/drive/MyDrive/Colab Notebooks/winequality_white.csv")

# df_red["color"] = "R"
# df_white["color"] = "W"
df_red["color"] = 0
df_white["color"] = 1

df_all=pd.concat([df_red,df_white],axis=0)

# df_all.head(5000)

df = pd.get_dummies(df_all, columns=["color"])
df_all.isnull().sum()
df_all.describe()



plt.subplots(figsize=(20,15))
ax = plt.axes()
ax.set_title("Wine Characteristic Correlation Heatmap (Reds)")
corr = df_red.corr()
sns.heatmap(corr, 
            xticklabels=corr.columns.values,
            yticklabels=corr.columns.values,
           cmap="Reds")
plt.show()


plt.subplots(figsize=(20,15))
ax = plt.axes()
ax.set_title("Wine Characteristic Correlation Heatmap (Reds)")
corr = df_red.corr()
sns.heatmap(corr, 
            xticklabels=corr.columns.values,
            yticklabels=corr.columns.values,
           cmap="Purples")
```

</details>





・回帰
線形回帰・ラッソ回帰・リッジ回帰・ElasticNet回帰で実施。
ElasticNetではl1_ratioとalphaを0～1まで変えて、分析を実施した。

<details><summary>回帰関連のコード</summary>

```python
import pandas as pd
from sklearn.model_selection import train_test_split


from sklearn.linear_model import LinearRegression
from sklearn.datasets import make_regression
from sklearn.linear_model import Lasso
from sklearn.linear_model import Ridge
from sklearn.linear_model import ElasticNet



# データの読み込み
wine_train_X, wine_test_X, wine_train_y, wine_test_y = train_test_split(df_all.drop('quality', axis=1),df_all['quality'], test_size=0.2, random_state=42)


# データを生成
# X, y = make_regression(n_samples=5000, n_features=50, n_informative=50, n_targets=1, noise=100.0, random_state=42)
# train_X, test_X, train_y, test_y = train_test_split(X, y, random_state=42)

#------------------------------
# 線形回帰
model = LinearRegression()
model.fit(wine_train_X, wine_train_y)

# test_X, test_yに対する決定係数を出力してください
print("Linear regression:{}".format(model.score(wine_test_X, wine_test_y)))


#------------------------------
# ラッソ回帰
model = Lasso()
model.fit(wine_train_X, wine_train_y)

# test_X, test_yに対する決定係数を出力してください
print("Lasso regression:{}".format(model.score(wine_test_X, wine_test_y)))

##------------------------------
# リッジ回帰
model = Ridge()
model.fit(wine_train_X, wine_train_y)

# test_X, test_yに対する決定係数を出力してください
print("Ridge regression:{}".format(model.score(wine_test_X, wine_test_y)))



#------------------------------
# ElasticNet回帰
l1_ratio_range = []
alpha_range = []
score_plot = []

for l1_ratio in np.arange(0.0, 1.1 , 0.1):
  for alpha in np.arange(0.0, 1.1 , 0.1):
      model = ElasticNet()

      model = ElasticNet(l1_ratio=l1_ratio , alpha =alpha)
      model.fit(wine_train_X, wine_train_y)

      # test_X, test_yに対する決定係数を出力してください
#      print( l1_ratio , alpha , "ElasticNet regression:{}". format(model.score(wine_test_X, wine_test_y)))

      l1_ratio_range.append(l1_ratio)
      alpha_range.append(alpha)
      score_plot.append(model.score(wine_test_X, wine_test_y))



import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

arr_l1_ratio_range = np.array(l1_ratio_range )
arr_alpha_range = np.array(alpha_range)
arr_score_plot = np.array(score_plot)

arr_l1_ratio_range  = np.reshape(arr_l1_ratio_range ,( 121 , 1))
arr_alpha_range  = np.reshape(arr_alpha_range ,( 121 , 1))
arr_score_plot  = np.reshape(arr_score_plot ,( 121 , 1))


print(x.shape)
print(y.shape)
print(z.shape)
print(arr_l1_ratio_range.shape)
print(arr_alpha_range.shape)
print(arr_score_plot.shape)


#描画エリアの作成
fig2 = plt.figure()
ax2 = fig2.add_subplot(projection='3d')


#散布図の作成
ax2.scatter(arr_l1_ratio_range , arr_alpha_range , arr_score_plot ,s=40,c="red")

ax2.view_init(elev=20, azim=50)

ax2.set_xlabel('l1_ratio')
ax2.set_ylabel('alpha')
ax2.set_zlabel('score')

#描画
plt.show()

```

</details>




決定係数は以下のとおりで
Linear regression:0.26715748512455806
Lasso regression:0.0017416836745601794
Ridge regression:0.2613865721139721

また、Elastic Netで解析したところ、alphaが小さいほど決定係数が大きくなる傾向が確認でき、上記の傾向と一致する。
一般的にはLassoの方が良い結果となることが多そうだが、今回はそれとは逆の結果となった。
但し、今回は決定係数が非常に低いため、そもそも、回帰で分析すること自体に限界があると考察するのが妥当である。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/152861/5745e3cf-cee6-eeff-b030-fcb45e04b049.png)


・分類
ランダムフォレストで分析した。
scoreが65%～70%程度の精度が出た。
n_estimatorsを1,10,100,1000と変化させたが、正解率が増加する結果となった。
※今回はmax_depth=20で固定した

<details><summary>ランダムフォレスト関連</summary>

```python
# モジュールのインポート
import matplotlib.pyplot as plt
from sklearn.datasets import make_classification
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split

# データの読み込み
wine_train_X, wine_test_X, wine_train_y, wine_test_y = train_test_split(df_all.drop('quality', axis=1),df_all['quality'], test_size=0.2, random_state=42)

# データを生成
# X, y = make_classification(n_samples=1000, n_features=4, n_informative=3, n_redundant=0, random_state=42)
# train_X, test_X, train_y, test_y = train_test_split(X, y, random_state=42)



# r_seedsの値の範囲(0から99)
r_seeds = [i for i in range(100)]

# 正解率を格納する空リストを作成
accuracy = []


# random_stateを変えながらモデルを学習
for seed in r_seeds:
    model = RandomForestClassifier(n_estimators = 100, max_depth=20 ,random_state=seed)
    model.fit(wine_train_X, wine_train_y)
    accuracy.append(model.score(wine_test_X, wine_test_y))

# グラフのプロット
plt.plot(r_seeds, accuracy)
plt.xlabel("seed")
plt.ylabel("accuracy")
plt.title("accuracy by changing seed")
plt.show()



# モジュールのインポート
import matplotlib.pyplot as plt
from sklearn.datasets import make_classification
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split

# データの読み込み
wine_train_X, wine_test_X, wine_train_y, wine_test_y = train_test_split(df_all.drop('quality', axis=1),df_all['quality'], test_size=0.2, random_state=42)

# r_seedsの値の範囲(0から99)
r_seeds = [i for i in range(20)]

# 正解率を格納する空リストを作成
accuracy = []
mean = []

import statistics

n_estimators_seeds = [i for i in [1, 10 , 100 , 1000]]
# random_stateを変えながらモデルを学習
for n_estimators in n_estimators_seeds:
  accuracy = []
  for seed in r_seeds:
    model = RandomForestClassifier(n_estimators , max_depth=20 ,random_state=seed)
    model.fit(wine_train_X, wine_train_y)
    accuracy.append(model.score(wine_test_X, wine_test_y))
  mean.append(statistics.mean(accuracy))
  print(accuracy)
  print(n_estimators , statistics.mean(accuracy))
print(mean)

# グラフのプロット
plt.plot(n_estimators_seeds, mean)
plt.xlabel("seed")
plt.ylabel("accuracy")
plt.title("accuracy by changing seed")
plt.show()
```

</details>

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/152861/15b0982d-3e57-2aea-6076-fd46e0ffaff9.png)

※n_estimators = 100, max_depth=20でseedを変えて実行した際のaccuracyの変化
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/152861/645937f8-29ac-b988-902b-3b8b4eaa5203.png)

また、主成分分析を行い、ロジスティクス分析を行った結果、55%程度の精度が出た。
そして、12個の軸のうち、6～8個程度のみ採用すると、精度が収束した。
以下3つのグラフからも、6～8個が重要度の上限であることが示唆される。
※固有値の大きい方から純に採用する仕組みになっている。

<details><summary>ロジスティク分析・主成分分析など</summary>

```python
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.linear_model import LogisticRegression

# データの読み込み
wine_train_X, wine_test_X, wine_train_y, wine_test_y = train_test_split(df_all.drop('quality', axis=1),df_all['quality'], test_size=0.2, random_state=42)

# 標準化のためのインスタンスを生成
sc =StandardScaler()
# トレーニングデータから変換モデルを学習し、テストデータに適用
X_train_std = sc.fit_transform(wine_train_X)
X_test_std = sc.transform(wine_test_X)


r_components = [i for i in range(12)]

accuracy_train = []
accuracy_test = []

for components in range(1, 13):
#for components in r_components:
  # 主成分分析のインスタンスを生成
  pca =PCA(n_components = components)
  # トレーニングデータから変換モデルを学習し、テストデータに適用
  X_train_pca = pca.fit_transform(X_train_std)
  X_test_pca = pca.transform(X_test_std)

  # ロジスティック回帰のインスタンスを生成
  lr = LogisticRegression()
  # 次元削減後のトレーニングデータで分類モデルを学習
  lr.fit(X_train_pca, wine_train_y)

  # スコアの表示
  accuracy_train.append(lr.score(X_train_pca, wine_train_y))
  accuracy_test.append(lr.score(X_test_pca, wine_test_y))



# グラフのプロット
plt.plot(r_components, accuracy_train)
plt.plot(r_components, accuracy_test)
plt.xlabel("components")
plt.ylabel("accuracy")
plt.title("Logistic Regression")
plt.show()





import pandas as pd
import matplotlib.pyplot as plt

from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.linear_model import LogisticRegression

# データの読み込み
wine_train_X, wine_test_X, wine_train_y, wine_test_y = train_test_split(df_all.drop('quality', axis=1),df_all['quality'], test_size=0.2, random_state=42)

# 標準化のためのインスタンスを生成
sc =StandardScaler()
# トレーニングデータから変換モデルを学習し、テストデータに適用
X_train_std = sc.fit_transform(wine_train_X)
X_test_std = sc.transform(wine_test_X)



import numpy as np
cov_mat = np.cov(X_train_std.T)
eigen_vals, eigen_vecs = np.linalg.eig(cov_mat)
print(eigen_vals)
print(eigen_vecs)


tot = sum(eigen_vals)
var_exp = [(i/tot) for i in sorted(eigen_vals, reverse=True)]
cum_var_exp = np.cumsum(var_exp)
import matplotlib.pyplot as plt
plt.bar(range(1,13), var_exp, alpha=0.5, align='center', label='individual exp var')
plt.step(range(1,13), cum_var_exp, where='mid', label='cum exp var')
plt.ylabel('ex')
plt.xlabel('prin')
plt.legend(loc='best')
plt.show()

# 固有値の大きいものから固有対を並び替える
eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:,i]) for i in range(len(eigen_vals))]
eigen_pairs.sort(key=lambda k: k[0], reverse=True)
# 固有値の大きい順に、対応する固有ベクトルを集める
w = np.hstack((eigen_pairs[0][1][:, np.newaxis],
               eigen_pairs[1][1][:, np.newaxis],
               eigen_pairs[2][1][:, np.newaxis],
               eigen_pairs[3][1][:, np.newaxis],
               eigen_pairs[4][1][:, np.newaxis],
               eigen_pairs[5][1][:, np.newaxis],
               eigen_pairs[6][1][:, np.newaxis],
               eigen_pairs[7][1][:, np.newaxis],
               eigen_pairs[8][1][:, np.newaxis],
               eigen_pairs[9][1][:, np.newaxis],
               eigen_pairs[10][1][:, np.newaxis],
               eigen_pairs[11][1][:, np.newaxis]))
print('Matrix W:\n',w)

r_components = [i for i in range(12)]






accuracy_train = []
accuracy_test = []

for components in range(1, 13):
#for components in r_components:
  '''
  # 主成分分析のインスタンスを生成
  pca =PCA(n_components = components)
  '''
  # 124 x 13 の特徴量を変換
  X_train_pca = X_train_std.dot(w[0:12 , 0:components])
  X_test_pca = X_test_std.dot(w[0:12 , 0:components])
  print(X_train_std.shape)
  print(w.shape)
  print(X_train_pca.shape)

  '''
  # トレーニングデータから変換モデルを学習し、テストデータに適用
  X_train_pca = pca.fit_transform(X_train_std)
  X_test_pca = pca.transform(X_test_std)
  '''
  # ロジスティック回帰のインスタンスを生成
  lr = LogisticRegression()
  # 次元削減後のトレーニングデータで分類モデルを学習
  lr.fit(X_train_pca, wine_train_y)

  # スコアの表示
  print(components ,
        lr.score(X_train_pca, wine_train_y) ,
        lr.score(X_test_pca, wine_test_y))
  accuracy_train.append(lr.score(X_train_pca, wine_train_y))
  accuracy_test.append(lr.score(X_test_pca, wine_test_y))



# グラフのプロット
plt.plot(r_components, accuracy_train)
plt.plot(r_components, accuracy_test)
plt.xlabel("components")
plt.ylabel("accuracy")
plt.title("Logistic Regression")
plt.show()
```

</details>

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/152861/3b9e0e16-9600-2252-4c58-7639d8c5f51a.png)
※主成分数とaccuracyの関係

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/152861/2b32848c-e0cb-7cd2-b3a5-78af68938eb0.png)
※固有値とその累積値の推移

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/152861/f11af4b3-ed6d-3fdd-8383-038b8cdff028.png)
※エルボー法でクラスター分類した際のクラスター数とDistributionの関係

<details><summary>クラスタリング</summary>

```python
import matplotlib.pyplot as plt
import numpy as np
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs

# データの読み込み
# X, Y, wine_train_y, wine_test_y = train_test_split(df_all.drop('quality', axis=1),df_all['quality'], test_size=0.01, random_state=42)


X = df_all.drop('quality', axis=1)
Y = df_all['quality']


NUM_cluster = 15

distortions = []
for i in range(1, NUM_cluster):                # クラスター数1~10を一気に計算
    km = KMeans(n_clusters=i,
                init="k-means++",     # k-means++法によりクラスタ中心を選択
                n_init=10,
                max_iter=300,
                random_state=0)
    km.fit(X)                         # クラスタリングを実行
    distortions.append(km.inertia_)   # km.fitするとkm.inertia_が得られる

# グラフのプロット
plt.plot(range(1, NUM_cluster), distortions, marker="o")
plt.xticks(np.arange(1, NUM_cluster + 1, 1))
plt.xlabel("Number of clusters")
plt.ylabel("Distortion")
plt.show()
```

</details>





以上の結果からは、ランダムフォレストの精度が非常に良いことが分かった



# 6. まとめ
今回は、代表的な分析手法のいくつかを試し、パラメータを変えてみて、分析手法の特性を確認した。
今後、他の分析手法を実装することで、データが変わっても、Auto-ML的な分析が可能になる。



# 7. 今後の課題
・LightGBMの実装、および、精度の確認
・カーネル主成分分析の考察
・クラスタリングの検証の幅の拡大


# 8. 参考資料
[scikit-learn] 8. make_regressionによる回帰用のランダムなデータの生成
https://sabopy.com/py/scikit-learn-8/

主成分分析(PCA)による教師なし次元削減
https://notebook.community/hide-tono/python-training/python-machine-learning/ch05/pca

LightGBMを超わかりやすく解説(理論+実装)【機械学習入門33】
https://datawokagaku.com/lightgbm/

機械学習で「分からん！」となりがちな正則化の図を分かりやすく解説
https://qiita.com/c60evaporator/items/784f0640004be4eefc51