PCAによる特徴量自動生成

# はじめに

　みなさん、データ分析・解析するにあたって、データの背景となる業界について下調べすると思います。
そんな時に、ぱっと与えられた特徴量の名前が専門用語だらけ...ということがよくあるのではないでしょうか。
もちろんこの下調べもデータサイエンスの醍醐味であるし、とても楽しい部分であることは間違いないです！！（笑）


でも、ちょっと思ったことありませんか？


「データサイエンスといっても結局、目的変数と説明変数の数的関係なんだからどうにか自動化できそうだよね....」

って！！！

はい、私もその一人です（笑）
もし背景知識なしでも、データの数値や変数のみで成果が出せれば...なんて夢があります。
そんなわけで実験してみました！

# データについて  



　まず仕組み作りはシンプルな例から！！と言うことでことで、kaggleのplay groundシーズン3.5のデータを使用したいと思います。
このデータは、ワインの品質を段階評価で予測する分類コンペに使われたものです。

また、構成は説明変数12の目的変数１と少なめの構成です。


# 1:  baselineの作成

何はともあれ、まずbaselineを作りましょう。

・目的変数はquality(ワインの品質)
・損失関数は二次の重み付きkappa係数
・modelは分類の勾配ブースティング
・validation設計はshuffleあり。(予測データと訓練データの個体の被りがなく、shuffle無しで再現性を高めたいところですけれど、今回、説明変数が少ないのと、ワインなので。多様性は少ないと判断します。)
・StratifiedKFoldでバランスよく分割(目的変数に偏りがまあまああるため。)

で作成します。

```python:baseline.py

from sklearn.model_selection import KFold,StratifiedKFold
import lightgbm as lgb
import numpy as np

def baseline(df,verbose=10):
    
    #trainとtestの結合データからtrainを取り出す作業
    
    train_x=df.iloc[:len(train)]
    labels=train_x['quality']
    train_x=train_x.drop('quality',axis=1)
    
    all_score=[]
    
    #validation設計

    kf=StratifiedKFold(n_splits=4,shuffle=True,random_state=123)
    
    for tr_index,va_index in list(kf.split(train_x,labels)):
        

        tr_x,va_x=train_x.iloc[tr_index],train_x.iloc[va_index]
        tr_y,va_y=labels.iloc[tr_index],labels.iloc[va_index]
        
        #model作成


        params={'objective':'multiclass',     
                'num_class':6,
                'learning_rate':0.1,
                'boosting_type':'gbdt',
                'importance_type':'gain',
                'num_iteration':10000,
                'metric':'multi_logloss',
                'num_leaves':16}

        model=lgb.LGBMClassifier(**params)
        
        #学習

        model.fit(tr_x,tr_y,eval_set=[(tr_x,tr_y),(va_x,va_y)],early_stopping_rounds=100,verbose=verbose)
        
        #推論/カッパ損失

        pred=model.predict(va_x)
        score=cohen_kappa_score(va_y,pred,weights='quadratic')
        all_score.append(score)
        
    #feature_importance_の可視化

    imp=pd.DataFrame({'name':train_x.columns,'imp':model.feature_importances_})
    display(imp.sort_values('imp',ascending=False)[0:50])
    
    #kappa損失
    
    print(np.mean(all_score))

    return model,np.mean(all_score)


```

実行結果は、

mean_kappa_score:0.48455655989665125

となりました。データを何もいじっていないのでこれくらいですね。
さあ、ここからどれだけスコアをあげられるのか。特徴量を追加していきます。

# 2: 特徴量生成の自動化

　今回の特徴量の自動生成の核としたのは、主成分分析(PCA)です。
と言ってもカラム全体にかけるのではなく、２つ、３つと特徴量を選択して成分ベクトルを生成していきます。これを特徴量として加え、ラッパー法で特徴量を選択します。
発想としましては、

「小さな集合に対する主成分分析なら、特徴量の相互関係をそれぞれあぶりだせるのではないか。また、その成分を追加すること、すなわち良質な特徴量生成となるのではないか。」

といった具合です。

自動化した関数がこちら。

```python:pca_spec.py

def pca_spec(labels1,labels2,df,m):

    #pca定義
    
    from sklearn.decomposition import PCA
    
    n_components=2
    whiten=False
    random_state=124
    
    pca=PCA(n_components=n_components,whiten=whiten,random_state=random_state)

    #kappaの初期値
    
    _,out_o=baseline(df)
    
    spects=[]

    kappa=[]

    #label総当たり式で特徴量生成
    #scoreが上がった場合のみ特徴量を残している
    #今回は第一成分のみを採用

    for i,label1 in enumerate(labels1):
        
        
        for j,label2 in enumerate(labels2):
            
            y=pca.fit_transform(df[[label1,label2]])
            y=pd.DataFrame(data=y[:,0],index=df.index)

            y.columns=[ label1+'_pca1_'+label2+'_'+str(m) for label in y.columns]
            df=pd.concat([df,y],axis=1)

            _,out=baseline(df,verbose=0)
            print(out)
            print(label2)
            
            
            if out>out_o:
                
                spects.append(label1+'_'+label2)
                kappa.append(out)
                
                print(label1+'_'+label2)
                
                out_o=out
                
                
            else:

                df=df.drop(str(label1+'_pca1_'+label2+'_'+str(m)),axis=1)
                
    return df,spects,kappa

```

　目的変数以外のデータのlabelを総当たり形式でpcaにかけて、第一主成分を訓練データに追加しました。その後modelに学習させ、kaapa_scoreが上がれば採用、上がらなければ不採用としました。
また、今回の方法の弱点として始めの方に生成されたpca成分を起点として相互特徴量の山を構築していくため、採用できない良質な特徴が埋もれてしまう可能性があることが挙げられます。
そのため、labelの順番をシャッフルして何度か実行し、最良の結果を選択するのがよいでしょう。

labelをランダムに並び替え10回ほど実行した結果、一番スコアがよかったものがこちらです。

|`特徴名の組み合わせ`    |*kappa_score*                  |
|----------|--------------------------|
|**初期値**|**0.48455655989665125**|
|**density---alcohol**  |***0.5024389483571687***|
|**density---free sulfur dioxide**|**0.5077476687421567**|
|**density---total sulfur dioxide**|**0.5122237000528017**|
|**free sulfur dioxide---volatile acidity**|**0.5177481328357583**|
|**alcohol---free sulfur dioxide**|**0.5194261932407493**|

ちなみに、二番目によかった結果がこちらです。

|`特徴名の組み合わせ`    |*kappa_score*                  |
|----------|--------------------------|
|**初期値**|**0.48455655989665125**|
|**density---free sulfur dioxide**|**0.5008467622707116**|
|**density---alcohol**|**0.5077476687421567**|
|**density---citric acid**|**0.5176992131296844**|



もとの値が0.4845であったのに対して、最終スコアが0.5194としっかり精度向上が見られました。データに対する知識なしとは思えませんね。
また、「density_free sulfur dioxide」、「density---alcohol」は二つに共通していることから、重要な関係であること。density(密度)の出現度から、相互関係において一つの注目すべき点であることが分かります。
今回はPCAの第一成分ですが、四則演算など、他の関係性も試してみると良いかもしれません。
このように、自力での特徴量生成においてのヒントも得られます。


また、学習率によって初期値が変動するので、0.01,0.05でも実験します。
条件をそろえるため、特徴量は0.1の最終結果を使用します。
結果は、

・0.01は初期値0.4930から0.5048240197927456
・0.05は初期値0.5019から0.5065597250161287

となりました。
なんと、初期値は0.01,0.05が高いものの、最終スコアは0.1が勝る結果になりました。
もしかすると、それぞれの学習率によって適切な特徴量が違うのではと思い、0.01,0.05でそれぞれ10回ずつ検証してみました。

0.01の結果がこちらです。初期値は0.4930です。

|`特徴名の組み合わせ`    |*kappa_score*                  |
|----------|--------------------------|
|**fixed acidity---chlorides**  |**0.5017524770773141**         |
|**chlorides---pH**  |**0.5066335122974438**         |
|**alcohol---density**  |**0.5148401350969285**         |
|**residual sugar---chlorides**  |**0.5177131987520922**         |
|**residual sugar---total sulfur dioxide**  |**0.5183609091743724**         |





0.05の結果がこちらです。初期値は0.5019です。

|`特徴名の組み合わせ`    |*kappa_score*                  |
|----------|--------------------------|
|**free sulfur dioxide---fixed acidity**|**0.5040698824927279**|
|**residual sugar---chlorides**|**0.5059535913659908**|
|**density---sulphates**|**0.5080298579996517**|
|**citric acid---alcohol**|**0.50843855183951**|
|**citric acid---volatile acidity**|**0.511783128707927**|
|**pH---free sulfur dioxide**|**0.5141857067441694**|

思った通りになりました。学習率ごとに最適な特徴量の組み合わせがあるように思えます。
(詳細な理由についてはもう少し探ってみます。)


さて、ここまで分類モデルで実験してきましたが、回帰モデルでも実験してみましょう。
回帰モデルなので、正解値と正解値の間を閾値で区切って、クラス予測します。また、閾値も学習することで精度向上を目指します。

閾値を学習するクラスはこちらです。

```python:coef_train.py

import scipy as sp
from functools import partial
from sklearn.metrics import cohen_kappa_score

class coef_train:
            
    def __init__(self,tr_pred,tr_y,va_pred,va_y):
                
        self.tr_pred=tr_pred
                
        self.tr_y=tr_y
                
        self.va_pred=va_pred
                
        self.va_y=va_y
                
        self.init=list(np.array([3.5, 4.5, 5.5, 6.5, 7.5]))
                
        self.coefs=0
                
    def kappa_score(self,coef,pred,y=0,p_s=False):

        preds=np.copy(pred)
            

        for p in range(len(pred)):

            if pred[p]>coef[4]:

                preds[p]=8

            elif (pred[p]<=coef[4])&(pred[p]>coef[3]):

                preds[p]=7

            elif (pred[p]<=coef[3])&(pred[p]>coef[2]):

                preds[p]=6

            elif (pred[p]<=coef[2])&(pred[p]>coef[1]):

                preds[p]=5

            elif (pred[p]<=coef[1])&(pred[p]>coef[0]):

                preds[p]=4

            elif pred[p]<=coef[0]:

                preds[p]=3
                
                
        if p_s==False:
                    
            score=cohen_kappa_score(preds,y,weights='quadratic')
                    
            #最小化すなわちscoreを-にすることでscoreの最大化をはかる。

            return -score
                
        else:

            return preds
            
    def fit(self):
                
        partial_kappa= partial(self.kappa_score,pred=self.tr_pred,y=self.tr_y)
           
        self.coefs= sp.optimize.minimize(partial_kappa, self.init, method='nelder-mead')

        score=-self.kappa_score(self.coefs['x'],pred=self.va_pred,y=self.va_y)
        
        return score
        
    def predict(self,x):
                
        pred=kappa_score(self.coefs['x'],pred=x,p_s=True)
                
        return pred


```
　初めは関数で書いていましたが、複雑になってしまったのでクラスにまとめました。
関数fitの説明を軽くすると、kohen_kappa_scoreをpartialでまとめ、'nelder-mead'(面積をだんだん小さくするような方法)で最小化し、学習した閾値で予測するといったものになります。
実は閾値の学習をするのは精度向上にかなり効果的なんです...。

結果はこちらになります。

|`特徴名の組み合わせ`    |*kappa_score*                  |
|----------|--------------------------|
|**初期値**|**0.5353487428014896**|
|**sulphates---citric acid**|**0.5562033171412317**|
|**chlorides---free sulfur dioxide**|**0.5568037938374937**|
|**chlorides---sulphates**|**0.5572875138427849**|
|**total sulfur dioxide---sulphates**|**0.55787904598868**|
|**fixed acidity---total sulfur dioxide**|**0.5584205606497032**|
|**fixed acidity---pH**|**0.5615331070595615**|
|**fixed acidity---density**|**0.563726212164812**|
|**fixed acidity---alcohol**|**0.5666017728250222**|
|**density---chlorides**|**0.5677486871307498**|


もう初期値から全然違いますね。（笑）
初めから回帰モデルでやれよっ！！ってくらいのスコア差ですが、考察を積み重ねていくのが大切です。

ちなみに、閾値を学習してない場合([3.5,4.5,5.5,6.5,7.5]のまま)の最良の結果はこちらです。
初期値は0.4860です。

|`特徴名の組み合わせ`    |*kappa_score*                  |
|----------|--------------------------|
|**citric acid---density**|**0.49025056268687794**|
|**alcohol---sulphates**|**0.49786286551520004**|
|**alcohol---density**|**0.5008424973171361**|
|**alcohol---fixed acidity**|**0.5040742968823384**|
|**pH---volatile acidity**|**0.5057734990823405**|
|**pH---sulphates**|**0.5078495719700824**|
|**pH---total sulfur dioxide**|**0.510364682036608**|
|**chlorides---citric acid**|**0.5135015590274816**|
|**chlorides---sulphates**|**0.5176273561261607**|
|**chlorides---fixed acidity**|**0.5179760903398334**|
|**residual sugar---sulphates**|**0.5191942588755812**|
|**residual sugar---density**|**0.5230455292652981**|
|**total sulfur dioxide---chlorides**|**0.5232497255750336**|

やはり、閾値の学習はかなり効果的ですね...。

ちなみに、新たに生成された特徴どうしをPCAにかけたものを追加することで、まだまだスコアが伸びます。
結果はこちらです。

|`特徴名の組み合わせ`    |*kappa_score* |
|----------|--------------------------|
|**fixed acidity_pca1_density---density_pca1_chlorides**|**0.571310719354684**|
|**fixed acidity_pca1_alcohol---chlorides_pca1_sulphates**|**0.5731771777064754**|


0.5731と、はじめの0.4846と比べるとかなりの精度向上を成し遂げられたのではないでしょうか。

また、impの値を見ると生成した特徴量が上位に入っていることもわかります。(下図)

|特徴名	|imp|
|----------|--------------------------|
alcohol	|2734.794672
sulphates|	1342.943081
sulphates_pca1_citric acid_0	|262.600824
chlorides_pca1_sulphates_0	|166.035159
total sulfur dioxide_pca1_sulphates_0	|158.866692
fixed acidity_pca1_total sulfur dioxide_0	|110.880541
fixed acidity_pca1_density_0	|91.980550
residual sugar	|91.977805
fixed acidity_pca1_alcohol_0|	85.059841
fixed acidity_pca1_pH_0	|82.564019
fixed acidity_pca1_density_0_density_pca1_chlorides|80.807850
fixed acidity_pca1_alcohol_0_chlorides_pca1_sulphates|79.634610
density	|78.664934
chlorides|	76.473363
fixed acidity	|74.028700
pH	|67.837130
free sulfur dioxide|	65.720670
volatile acidity|	60.174443
chlorides_pca1_free sulfur dioxide_0|	58.323490
density_pca1_chlorides_0|	55.214166
total sulfur dioxide|	30.461660
citric acid	|29.461067











# まとめ

　当初の予想通り、特徴量生成においてのPCAの威力を示せた実験になりました。ワインについての知識がなくても、予測精度の向上を図ることができました。また、この解析を行うことで通常のデータ解析のみでは見つけにくいデータの構造や、モデルの精度向上のための着眼点が分かるのではないでしょうか。
他にもICA(独立成分分析)やKernelPCA(非線形的な関係)など、同じような仕組みで追加できる生成法はいくつもあります。
次回はこれらでも実験してみようと思います。

# 追加考察

　今回、精度のよかった回帰モデル搭載の仕組みを10周させ、採用した組み合わせ特徴量の集計を取りました。
3回以上の出現についてまとめます。


|`特徴名の組み合わせ`    |*採用回数*                  |
|----------|--------------------------|
|**citric acidとfixed acidity**|**6**|
|**free sulfur dioxideとresidual sugar**|**4**|
|**fixed acidityとdensity**|**3**|
|**pHとdensity**|**3**|
|**volatile acidityとtotal sulfur dioxide**|**3**|
|**total sulfur dioxideとdensity**|**3**|

　一番採用されているのは、クエン酸と固定酸度の組み合わせですね。固定酸度ってブドウ由来の酸でその成分の一つがクエン酸なんですよね...。
とっても良い分散が見つけられそうです。
また、クエン酸はワインの味に直結しますし、添加したクエン酸か、自然由来か、のように品質に影響しそうな気がします。

　二番目に採用されているのは、遊離二酸化硫黄と残留糖の組み合わせです。
ワイン全ての亜硫酸のうち、一部が糖などと結合した結合型、その他は遊離亜硫酸(遊離二酸化硫黄)となります。
...いい分散が見つけられそうですね。
ちなみに遊離亜硫酸はso2とhso3-に分けられるのですが、このso2がワインのツンとくる匂いの原因なので、食した感想に直結しそうですよね。
また、phも絡んだ話にはなるのですが、phが高い(酸性が弱い)と遊離亜硫酸のうちso2でいられる割合が低くなります。このso2は細菌の殺菌や増殖の阻止に役立つのですが、これが減ってしまえば微生物汚染が進んでしまいます。(=腐敗)
また、phが高いほど微生物は増殖しやすいので二重に汚染が進みます。
ただ、あんまりphを上げても今度は酸味や匂いがきつくなってしまうので、適切なバランスが大事なようです。(調べたら出てきます、この数値は結構指標になります。)

次に揮発性酸性度と総二酸化硫黄ですが、揮発性酸性度の中身はほぼ酢酸です。so2と酢酸なので、匂いに関する分散がありそうですね。

他には「固定酸度と密度」、「総二酸化硫黄と密度」、「phと密度」、とやはり密度の関係は洗い直してみてもいいくらい重要そうですね。

　少し長くなりましたが、予測に関係ありそうな内容がピックアップされているのが分かります。
結果を言うと、この仕組みと自力の組み合わせが一番効果発揮しそうですね。
また、今回は回帰モデルについて集計を取りましたが、分類モデルではまた違う結果が得られるかもしれません。
もう少し色々試してみようと思います。

※（補足）
corr()で相関を出して、高いもの同士をPCAして特徴を生成してみましたが、相関が高いからといってPCAが予測にいい影響を与えられるわけではありませんでした。
あくまで、品質も考慮した相関関係であることが重要だなと感じました。また、いくつか集った段階で初めて効果が発揮される特徴もあるので、もう少し深く相関については探る必要がありそうです。