【Marketing×Data Science勉強会#1】Pythonを用いたロジスティック回帰による広告クリック率の予測

# はじめに
本記事ではData Science/Machine Learningを活用することで、データドリブンにマーケティング活動を進めていく方法について、6つのテーマに分けて様々な多変量解析や機械学習モデルをPythonで構築することで理解を深めていく。
### 本企画の目的
1. 単なるデータの可視化等の簡易な分析スキルからさらに一歩進んだ分析スキルを身につけ、データエンジニアリングよりは多変量解析や機械学習の手法をつかった予測モデルや要因探索を扱えるようにすること
2. 分析手法ありきでその手法が、実際のデータコンサルティング/マーケティングにおいてどのような場面で使えるのかイメージを沸かせること

### 本企画の構成
頻度は月に一度で扱うトピックは以下の6つ。
1. ロジスティック回帰による広告クリック率の予測
2. 決定木分析による顧客のクラスタリング
3. 主成分・因子分析によるブランド知覚マップ
4. 協調フィルタリング分析によるレコメンドモデルの構築
5. テキストマイニングによるTwitter口コミ分析 
6. RNN, LSTMの時系列予測によるPV数予測

### 本記事のゴール
本記事では一つ目のテーマ**「ロジスティック回帰による広告クリック率の予測」**([Source](https://medium.com/@gouravbais08/ad-click-prediction-9478a9ecf066))をとりあげ、以下の二点をゴールに進めていく。

1. ロジスティック回帰モデルの大まかな理論が理解し、Pythonでモデルを構築できるようになること
    1. プログラミングの基本「入力」「演算」「出力」を理解する
    2. 細かなPythonのコーディングよりも、何のための操作なのかを理解できるようにする(出力→演算→入力をイメージする。ググれるようにする)


2. 一連のロジスティック回帰分析をとして、マーケティングでの応用例を理解すること
    1. 実際にどのようなデータを入力して、どのようなデータを出力するのかイメージを掴む
    2. 出力結果から実プロジェクトにおけるマーケティング施策に活かせるインサイトは何かを考えれるようにする

今回は[このスライド](https://docs.google.com/presentation/d/1D3-DvdkuCR678s5E7lx287Cr7ZCsj3HG5ZdMbra20hw/edit#slide=id.g78043fb566_0_117)に沿って進めていくので随時確認してください。
このノートでは基本的にコードベースで全体の流れとコードの説明を行う。

# ロジスティック回帰モデルの構築
## ロジスティック回帰について
ロジスティック回帰と混同行列について理解を深めたい場合以下を確認ください。
[Scikit-learn でロジスティック回帰（クラス分類編）](https://qiita.com/0NE_shoT_/items/b702ab482466df6e5569)
[混同行列について](https://qiita.com/TsutomuNakamura/items/a1a6a02cb9bb0dcbb37f)

## Pythonによるモデル構築
### Step 0 : Building Environment
今回はデータ分析ツールJupyter NotebooksかAzure Notebooksを使うので、以下の資料を参考に各自インストールして環境構築します。
[環境構築資料](https://drive.google.com/drive/folders/1yfyn8KpC9VuIGZGCh0plpYaetpUxsjGS)

### Step 1 : Importing Dependencies and Loading Data set
Pythonにはデータ分析に特化したライブラリが豊富にある。

* pandas：データフレーム
* numpy：配列
* matplotlib.pyplot, seaborn：可視化  [seaborn例](https://qiita.com/hik0107/items/3dc541158fceb3156ee0)

```python:AdClickPrediction.ipynb
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
```
データの読み込み

```python:AdClickPrediction.ipynb

amex = pd.read_csv('Ad_click_prediction_train.csv')
amex.head()
amex.info()
```
![スクリーンショット 2020-06-13 22.45.25.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/266935/065541c9-164e-c38c-0e66-9c9e50314e65.png)



### Step 2 : Exploratory Analysis
データの特徴を把握するために、データを可視化する。

```python:AdClickPrediction.ipynb

sns.countplot(amex['is_click'], hue=amex['gender'])
sns.countplot(amex['is_click'], hue=amex['age_level'])
sns.heatmap(amex.isna(), yticklabels=False, cbar=False)
```
性別ごとのクリック有無の総数
![スクリーンショット 2020-06-13 22.45.06.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/266935/7c7e98e0-f39f-b1be-5eff-3a26ce07a5a1.png)
年代ごとのクリック有無の総数
![スクリーンショット 2020-06-13 22.45.13.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/266935/88d7de71-c1a6-680a-4793-0d1c2200f60b.png)
欠損値ヒートマップ
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/266935/d78773ea-20ac-29f4-f1e3-864fbdeae3ad.png)




### Step 3 : Data Cleaning
実際に得られるデータは顧客の入力漏れやテキストなどのモデルに直接入力できない型のデータが存在するため、データの前処理を行う必要がある。この処理がデータ分析において膨大な時間と労力を要する。今回の処理もコードは複雑に見えるが、コードの内容を理解するというよりかはどのようなフローで何をしているのかを理解することに力を入れましょう。

前処理は具体的に以下のようなことを実施する。

* カテゴリーデータの処理
* 欠損値処理
* 特徴量の変換,追加
* 次元削除
* エンコーディング

[参考リンク](https://qiita.com/takuya_tsurumi/items/53b9e3f7427b631b17cf)


実際に欠損値を確認してみる。

```python:AdClickPrediction.ipynb
amex.isnull().sum()
```
まずはプロダクトのカテゴリをエンコードする必要があるが、同じディレクトリに以下のファイルを作成する。

```python:product_category.py
import pandas as pd

def productCatagory2(columns):
    product = columns[0]
    campaign = columns[1]
    webpage = columns[2]
    catagory = columns[3]
    
    if pd.isnull(catagory):
        
        if(product == 'H' and campaign == 404347  and webpage == 53587 ):
            return(235358.0)
        elif(product == 'C' and campaign == 404347  and webpage == 53587 ):
            return(234846.0)
        elif(product == 'A' and campaign == 404347  and webpage == 53587 ):
            return(234846.0)
        elif(product == 'F' and campaign == 404347  and webpage == 53587 ):
            return(301147.0)
        elif(product == 'H' and campaign == 404347  and webpage == 53587 ):
            return(99226.0)
        elif(product == 'B' and campaign == 404347  and webpage == 53587 ):
            return(408790.0)
        elif(product == 'D' and campaign == 404347  and webpage == 53587 ):
            return(146115.0)
        elif(product == 'E' and campaign == 404347  and webpage == 53587 ):
            return(146115.0)
        elif(product == 'G' and campaign == 404347  and webpage == 53587 ):
            return(327439.0)
        elif(product == 'I' and campaign == 404347  and webpage == 53587 ):
            return(419804.0)
        elif(product == 'J' and campaign == 404347  and webpage == 53587 ):
            return(450184.0)
        
        
        elif(product == 'B' and campaign == 414149  and webpage == 45962 ):
            return(143597.0)
        elif(product == 'C' and campaign == 414149  and webpage == 45962 ):
            return(254132.0)
        elif(product == 'D' and campaign == 414149  and webpage == 45962 ):
            return(254132.0)
        elif(product == 'E' and campaign == 414149  and webpage == 45962 ):
            return(254132.0)
        
        
        elif(product == 'A' and campaign == 396664  and webpage == 51181 ):
            return(143597.0)
        elif(product == 'F' and campaign == 396664  and webpage == 51181 ):
            return(18595.0)
        
        
        elif(product == 'C' and campaign == 98970  and webpage == 6970 ):
            return(270147.0)
        elif(product == 'I' and campaign == 98970  and webpage == 6970 ):
            return(327439.0)
        elif(product == 'B' and campaign == 98970  and webpage == 6970 ):
            return(202351.0)
        
        
        elif(product == 'C' and campaign == 359520  and webpage == 13787 ):
            return(408831.0)
        
        
        elif(product == 'A' and campaign == 118601  and webpage == 28529 ):
            return(82527.0)
        elif(product == 'B' and campaign == 118601  and webpage == 28529 ):
            return(82527.0)
        elif(product == 'H' and campaign == 118601  and webpage == 28529 ):
            return(82527.0)
        elif(product == 'G' and campaign == 118601  and webpage == 28529 ):
            return(82527.0)
        elif(product == 'I' and campaign == 118601  and webpage == 28529 ):
            return(82527.0)
        
        
        elif(product == 'B' and campaign == 82320  and webpage == 1734 ):
            return(235358.0)
        elif(product == 'G' and campaign == 82320  and webpage == 1734):
            return(255689.0)
        elif(product == 'H' and campaign == 82320  and webpage == 1734):
            return(255689.0)
        elif(product == 'I' and campaign == 82320  and webpage == 1734):
            return(18595.0)
        elif(product == 'C' and campaign == 82320  and webpage == 1734):
            return(18595.0)
        else:
            return(211061.500000)
            
        
    else:
        return(catagory)
```

Jupyter Notebooksに戻り、先ほどのファイルを使ってプロダクトカテゴリ2をエンコードする。

```python:AdClickPrediction.ipynb

from product_category2 import productCatagory2
amex['product_category_2'] = amex[['product','campaign_id','webpage_id','product_category_2']].apply(productCatagory2,axis = 1)
amex.head()
```
プロダクトカテゴリがエンコードされていることが確認できる。

次に、genderとuserGroupIDをエンコードし、時間等の他の特徴量も処理する。

```python:AdClickPrediction.ipynb

def gender(col):
    pro1 = col[0]
    gender = col[1]
    
    if pd.isnull(gender):
        if(pro1 <= 3):
            return('Male')
        else:
            return('Female')
    else:
        return(gender)

def userGroupId(col):
    gender = col[0]
    user = col[1]
    
    if pd.isnull(user):
        if(gender == 'Female'):
            return(8)
        else:
            return(3)
    else:
        return(user)

#gender
amex['gender'] = amex[['product_category_1','gender']].apply(gender,axis = 1)

#userID
amex['user_group_id'] = amex[['gender','user_group_id']].apply(userGroupId,axis = 1)
amex['age_level'].fillna(3,inplace = True)
amex['user_depth'].fillna(3,inplace = True)
amex['city_development_index'].fillna(3,inplace = True)

amex['Time'] = amex['DateTime'].apply(lambda x : x.split(' ')[1])
amex['Date'] = amex['DateTime'].apply(lambda x : x.split(' ')[0])
amex.drop(['DateTime'],inplace=True,axis=1)

amex['Year'] = amex['Date'].apply( lambda x : x.split('-')[0] )
amex['Month'] = amex['Date'].apply( lambda x : x.split('-')[1] )
amex['Date'] = amex['Date'].apply( lambda x : x.split('-')[2] )

amex['Hour'] = amex['Time'].apply( lambda x : x.split(':')[0] )
amex['Minutes'] = amex['Time'].apply( lambda x : x.split(':')[1] )

amex.drop(['Date'],axis = 1,inplace=True)
amex.drop(['Time'],axis = 1,inplace=True)

gender = pd.get_dummies(amex['gender'],drop_first=True)
product_type = pd.get_dummies(amex['product'])

amex.drop(['gender'],axis = 1,inplace=True)
amex.drop(['product'],axis = 1,inplace = True)
amex['gender'] = gender 

amex = pd.concat([amex,product_type],axis=1)
```

先ほどのデータテーブルから商品カテゴリがエンコードされた以下のような状態になった。

![スクリーンショット 2020-06-13 22.45.25.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/266935/eaeec854-bcf0-05bb-60b2-c7a8ed79e045.png)

このステップにおけるコードは難しいですが、実際にどのようなことをしているのかを理解できるようにしましょう！

### Step 4 : Train Test Split
さて、やっと長い前処理が終わりました。。
次にモデルの構築の前に学習データとテストデータに分ける。


```python:AdClickPrediction.ipynb

from sklearn.model_selection import train_test_split
x = amex.drop(['is_click'],axis=1)
y = amex['is_click']
X_train,x_test,Y_train,y_test = train_test_split(x,y,test_size = 0.6,random_state = 101)
```

### Step 5 : Training the Machine Learning Model
ついに念願のモデル構築！！！
先ほど作成した入力データをロジスティック回帰モデルに入れて結果を予測する。

```python:AdClickPrediction.ipynb

from sklearn.linear_model import LogisticRegression

lr = LogisticRegression() # ロジスティック回帰モデルのインスタンスを作成
lr.fit(X_train, Y_train) # ロジスティック回帰モデルの重みを学習
y_pred = lr.predict(x_test)
print(y_pred)
```

### Step 6 : Model Evaluation
さて、先ほど構築したモデルが実際にどの程度の精度なのか、どの特徴量がクリック率に大きく影響しているのか？というのを混同行列とオッズ非を用いてみる。


```python:AdClickPrediction.ipynb

from sklearn.metrics import confusion_matrix, accuracy_score

print('Confusion Matrix = \n', confusion_matrix(y_true=y_test, y_pred=y_pred))
print('\n Accuracy = ', accuracy_score(y_true=y_test, y_pred=y_pred))

print(259251/(259251+18724))

# Check the each input influence

print("coefficient = ", lr.coef_)
print("intercept = ", lr.intercept_)

lr.coef_.reshape(-1,)

odds = pd.DataFrame({"Name":amex.columns[:len(lr.coef_.reshape(-1,))],
                    "Coefficients":np.abs(lr.coef_.reshape(-1,))}).sort_values(by='Coefficients')
print(odds)
```

`Accuracy =  0.9326414245885422`
という高精度でクリック率を予測できるモデルを作成することができた。


ここまで、予測と推論に重きを置いて実際にコードを用いて説明していきましたが、次に統計解析を用いてよりモデルの解釈性を上げるためにクリック数(目的変数)と説明変数の相関を見ていく。

### Step 7 : Statistical Analysis

次にモデルの解釈性に統計モデルを用いてアプローチしていく。
それぞれの変数がクリック率に対してどの程度の影響を与えているのか？というイシューに対して重回帰モデルと異なり、係数をそのまま解釈することができないので、定量的な指標として**オッズ比**を用いる。

モデルに加える説明変数を選んでいく際の方法には以下の二つの点に注意する。

* 興味のある変数から選んでいく。
* その際に多重共線性に注意する。（説明変数同士に0.7以上の相関があると怪しい）

```python:AdClickPrediction.ipynb

# 統計解析に必要なライブラリをインポート
import statsmodels.formula.api as smf
import statsmodels.api as sm

# ヒートマップを作成
plt.figure(figsize=(20, 15))
sns.heatmap(amex.corr(),  cmap= sns.color_palette('coolwarm', 10), annot=True,fmt='.2f', vmin = -1, vmax = 1)
```

![heatmap.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/266935/5f7b61b7-8837-432b-6719-6a496d237f27.png)
**gender とuser_group_idに強い相関がある。このペアだけ気を付ければ大丈夫そう！**

今回は説明変数をage_level, city_development_index, gender,  user_id, webpage_idとしてモデルを構築してみる。

```python:AdClickPrediction.ipynb

## Modeling 
mod_glm = smf.glm(formula = "is_click ~  age_level+city_development_index + gender+user_id + webpage_id "
                  , data = amex, family = sm.families.Binomial()).fit()
mod_glm.summary()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/266935/af0a1dca-b67d-f5cf-2d34-27b9f700a1ec.png)

P値が0.05よりも小さい変数はクリック率に対して統計的に有意な影響を与えていることがわかる
**genderは統計的に有意な影響を与えていない。**


次に、モデルで推定された係数を指数変換することでオッズ比を推定して、定量的な影響度合いを測る。

```python:AdClickPrediction.ipynb

import scipy as sp
# モデルの係数を指数変換します。その結果をparam_expという名前で保存する
param_exp = sp.exp(mod_glm.params)
# 実際に中身をみてみる
param_exp
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/266935/82f7fa12-3ad8-2377-692b-8e4aac826023.png)

これだと少しみにくいので綺麗にしていく。

```python:AdClickPrediction.ipynb
# 先ほどのparam_expをデータフレームの形に
param_exp = pd.DataFrame(param_exp)
# 列の名前を指定
param_exp.columns = ['オッズ比']
# 昇順に並び替え
param_exp.sort_values('オッズ比')
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/266935/28e84f32-f427-fd21-603e-973fadaac4d5.png)


かなり微妙な結果ですね。おそらくis_clickと相関の強い変数があまりなかったことが原因か、もしかしたら、データでトラック出来ていない変数が大きな影響を持っているのかもしれない。

強いていうならば、

- **age**
    - levelが増えるにつれてクリック確率は下がっていく
    - 若い人の方がクリックする可能性が高い
- **city_development_index**
    - 発展している地域ほど、クリック確率は下がっていく

- **gender**
    - 0→1になることで、クリック確率が1.02倍になる
    - 男女間でクリックするかどうかの差はあまりみられない
    - ただし、今回のモデルの場合、統計的に有意な影響は確認出来なかった

**1つのモデルだけ構築して、結果を決めるのはオススメできない**
説明変数を入れ替えて、結果が変化するかどうかを確認する必要がある。結果がさほど変わらないのであれば、結果の頑健性が示されたことになる。

### Step8 : Challenge！！
これまで学んだことを生かして、[このKaggle](https://www.kaggle.com/swankhede/customer-classification-using-logistic-regression/data)にアクセスして、顧客が商品を購入するかどうかをロジスティック回帰で予測してみましょう！