タイタニックの質的変数で数量化理論

Kaggleの有名なコンペの一つにタイタニックで生き残る人を見つけるのがあります。
そこには質的変数が非常に多いので質的変数を使って数量化理論をやってみようと思います。
# ライブラリのインポート
```Python3
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import pandas as pd
import statsmodels.api as sm
import numpy as np
import matplotlib.pyplot as plt
```

# データの読み込み
欠損値は下手に代表値で補完すると実データと異なる可能性があるので消します。
```Python3
df = pd.read_csv("titanic.csv")
df = df.dropna()
df.head()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/b002e33d-c5f7-dfc1-80b4-563a4584047b.png)



# 説明変数選びとダミー変数
質的変数のみ選び必要なものはダミー変数にして量的変数は説明変数から削除します。
また、多重共線性の観点からダミー変数は各一つずつ削除します。
```Python3
x1 = pd.get_dummies(df.drop(["PassengerId", "Survived","Pclass","Name", "Age", "Fare", "Ticket", "Cabin"], axis=1))
x1 = x1.drop(["Sex_female", "Embarked_C"], axis=1)
x1.head()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/88a304d7-55a6-c544-47d6-03556d277e9e.png)

# 数量化理論
## 数量化理論Ⅰ類
重回帰分析を使います。
### 年齢を分析
```Python3
y1 = df["Age"]
model1 = sm.OLS(y1, x1).fit()
model1.summary()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/4d61e581-1220-634b-3a94-946c3af3eb84.png)

### 運賃を分析
```Python3
y2 = df["Fare"]
model2 = sm.OLS(y2, x1).fit()
model2.summary()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/5fefed97-6112-7ab9-1d17-85679942ef36.png)

## 数量化理論Ⅱ類
ロジスティック回帰を使います。
生存したかどうかを分析します。
```Python3
y3 = df["Survived"]
x3 = x1
model3 = sm.Logit(y3, x3).fit_regularized()
model3.summary()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/af7826d5-ce44-0476-21ac-62889cc061f6.png)
女性が生き残りやすくEmbarkedがSの方が生き残りやすくSibSpが1の人が生き残りやすいそうですね。
## 数量化理論Ⅲ類
解釈しやすいようにbiplotにします。
```Python3
x4 = pd.get_dummies(df.drop(["PassengerId", "Survived","Pclass","Name", "Age", "Fare", "Ticket", "Cabin"], axis=1))
col = x4.columns
ss = StandardScaler()
ss.fit(x4)
x4 = ss.transform(x4)
model4 = PCA()
model4.fit(x4)
tx = model4.transform(x4)
exp = model4.explained_variance_ratio_
com = model4.components_
fac = []
for i in range(len(exp)):
    fac.append(np.sqrt(exp[i])*com[i])
fig, ax = plt.subplots()
ax1 = ax.twinx()
ax1.scatter(tx[:, 0], tx[:, 1], c=y3, cmap="brg")
ylim = [abs(max(fac[1])), abs(min(fac[1]))]
xlim = [abs(max(fac[0])), abs(min(fac[0]))]
ax2 = ax.twiny()
for i in range(len(col)):
    ax2.plot([0, fac[0][i]], [0, fac[1][i]], color="#FF0000")
    ax2.text(fac[0][i], fac[1][i], col[i])
ax2.set_xlim(-max(xlim), max(xlim))
ax2.set_ylim(-max(ylim), max(ylim))
plt.show()
df_fac = pd.DataFrame(fac)
df_fac.columns = col
ind = []
for i in range(len(fac)):
    ind.append("PC%d"%(i+1))
df_fac.index = ind
df_fac
```
![Untitled.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/4aa25ebd-d6a7-6700-2374-0b2e011b1e66.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/2b2da7f6-3ec8-0f49-385e-e320efb7b8fd.png)

青が生き残れず、緑が生き残れた人です。
Ⅱ類の分析と同様に生き残れた人が女性でEmbarkedがSの方が生き残れやすそうですね。

# 余談
## 数量化理論Ⅱ類の精度
```Python3
from sklearn.metrics import classification_report

pred3 = model3.predict(x3)
y_pred3 = np.where(pred3 >= 0.5, 1, 0)
print(classification_report(y3, y_pred3))
```
```
              precision    recall  f1-score   support

           0       0.59      0.70      0.64        60
           1       0.84      0.76      0.80       123

    accuracy                           0.74       183
   macro avg       0.72      0.73      0.72       183
weighted avg       0.76      0.74      0.75       183
```
意外と74%当たってますね、まあ訓練データとテストデータが同じなので何とも言えませんけど。

## LightGBMで予測精度分布
```Python3
from sklearn.model_selection import train_test_split
from lightgbm import LGBMClassifier

score = []
for i in range(200):
    x_train, x_test, y_train, y_test = train_test_split(x3, y3, random_state=i, test_size=0.3)
    model5 = LGBMClassifier()
    model5.fit(x_train, y_train)
    score.append(model5.score(x_test, y_test))
df_scr = pd.DataFrame(score).describe()
df_scr.columns = ["Score"]
df_scr
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/de22e1ea-091c-2d59-3b8d-45bdc95fff3f.png)

```Python3
plt.boxplot(score)
plt.scatter([1],[df_scr.loc["mean", "Score"]])
plt.show()
```
![Untitled-1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/7389e9fa-5c30-e1c0-23b9-5781f0947ec4.png)

最大値は外れ値ではないみたいですね。むしろ最小値が外れ値のようです。
いずれにせよ70%の精度が質的変数だけでもあるみたいです。

## 普通に予測
数量化理論を使わず量的変数も使って予測してみました
```Python3
score = []
x = pd.get_dummies(df.drop(["PassengerId", "Survived", "Name", "Ticket", "Cabin"], axis=1))
x = x.drop(["Sex_female", "Embarked_C"], axis=1)
y = df["Survived"]
for i in range(200):
    x_train, x_test, y_train, y_test = train_test_split(x, y, random_state=i, test_size=0.3)
    model6 = LGBMClassifier()
    model6.fit(x_train, y_train)
    score.append(model6.score(x_test, y_test))
df_scr = pd.DataFrame(score).describe()
df_scr.columns = ["Score"]
df_scr
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/63226f3b-b092-d0bb-0a69-8351527dc01d.png)

```Python3
plt.boxplot(score)
plt.scatter([1],[df_scr.loc["mean", "Score"]], marker="x", color="#000000")
plt.show()
```
![Untitled.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/3cadc5a2-36a6-ce2f-ec1a-3fa912d129fb.png)

あまり安定していないみたいですが中央値までは精度が向上したみたいです。
なので低い精度は出にくくなったと言えます。ただし精度に限界はあります。

(裏を返せば質的変数が割と目的変数に関連しているのかも)