大阪都構想の投票結果を区ごとに分析してみた

# はじめに
私は現在大阪市に住んでおり、一昨日の都構想投票はテレビにかじりつきながら見ていました。

経過を見ていて思ったのが、**区ごとの結果の差が顕著**に出ており、分析対象として適したデータが得られそうだと感じたため、詳しく分析してみました。

※下図はおおさか維新の会HP掲載の、都構想における新旧区分け
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/f4d0d6a7-2ee8-3fc3-f4f2-671c8d456aae.png)

[GitHubに、使用したスクリプトやクレンジング後のデータをアップロードしています](https://github.com/c60evaporator/osaka_metropolis)

[Qiitaのガイドライン](https://help.qiita.com/ja/articles/qiita-community-guideline)にあるように、あくまで技術記事としての領分を超えないよう、政治的な深い考察は避け、得られた事実のみを列挙していこうと思います。

また、私は因果推論のような高度な分析のスキルは持ち合わせていないので、「さらに深い知見を得るためにはこうしたらいい」
というような**手法に関するアドバイスがございましたら、コメント頂けると大変ありがたいです！**

# 結論
結論に至るまでの手順は次章以降で述べますが、以下のような結論が得られました。

**「30～50代割合が高い区ほど、賛成率が高い」**
**「平均世帯人数が少ない(ファミリー層が少ない)区ほど、賛成率が高い」**
**「緯度が高い(北側の)区ほど、賛成率が高い」**
**「男性比率が高い区ほど、賛成率が高い」**

※影響度の高い順に並べています

# 分析の手順
データ分析の分野ではおなじみの、[TJO氏のブログ](https://tjo.hatenablog.com/entry/2013/11/14/202610)に記載されているように、

1. 目的の明確化（目的変数の設定）
2. 特徴量選択（説明変数の選択）
3. データの収集とクレンジング
4. 分析と結果の可視化
5. 考察とまとめ(今回は深入りしない)

という手順で進めていきます

## 1. 目的の明確化
今回、私が気になったのは、下記の2点です
①賛成が多い区と少ない区の偏り
②前回投票との差

①に関しては、下図（毎日新聞のサイトより引用）のように、賛成率に明確な地理的偏りが見られます（市の中央部北寄りの賛成率が高い？）
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/881d375f-0a87-6198-d413-b86cad9f6df8.png)
この偏りの陰にある要因（特徴量）を、数値的な根拠をもってあぶり出したいと思いました。

ですので、
**目的：賛成率の高低に対し、根拠のある要因候補を抽出する**
**目的変数：賛成率**
**説明変数：要因候補となる特徴量**
として分析を進めようと思います

②に関しては、前回と比べ賛成率が大きく変化した区の要因を抽出したいので、
**目的変数：今回の賛成率**
**説明変数：前回の賛成率**
として、外れ値をグルーピングしたいと思います。

## 2. 特徴量選択（説明変数の選択）
前述のように、目的変数である賛成率に影響を与えそうな指標を、特徴量として抽出していきます。

### 年齢構成
下図NHKの集計のように、出口調査では30～50代の働き盛りの賛成率が高く、
60代以上の高齢者、および20代の若者の賛成率が比較的低いです。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/1b7dfada-9e15-3354-1f8f-3dc30b6dde6c.png)
今回は、20代のみ傾向が逆転していることが少し特異だと感じたので、この影響の有無を考慮し、
**特徴量1「60歳以上割合」**(20代の逆転現象を考慮しない)
**特徴量2「30～50代割合」**(20代の逆転現象を考慮)
を特徴量に加えます

### 男女比率
年齢ほどの差はありませんが、男女比も影響がありそうなので、
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/3a02cfb6-fe43-80e5-c6c9-4a70c04c246c.png)
**特徴量3「男性比率」**
を特徴量に加えます

### 新区役所（本庁舎）からの距離
下図のように、新区役所（本庁舎）は市の中心部に集中して置かれる計画となっていました。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/acd8aa0d-4798-1083-a039-53858827adb9.png)

私は本庁舎が置かれる区とは異なる区に住んでいますが、「自分が住んでいる地域が発展から置いていかれるのでは？」という不安を少し感じたので、この感覚（ドメイン知識）を信じ、
**特徴量4「本庁舎が置かれる駅から、旧区中心駅までの所要時間」**
を、特徴量として加えました

### 平均世帯人数、平均年収
これらは明確な根拠はありませんが、何となく投票行動に影響を与えそうで、かつデータも入手しやすいので、特徴量に加えました
**特徴量5「平均世帯人数」**
**特徴量6「平均1人あたり年収」**

## 3. データの収集とクレンジング
目的変数となる賛成率はたくさんソースがあるので割愛しますが、
各種特徴量は下記の場所・クレンジング方法で収集しました
### 具体的な特徴量
下記6特徴量を使用します
#### 「1.60歳以上割合」「2.30～50代割合」「3.男性比率」
データ入手方法：下記大阪市HPより2020年の年齢・男女別推計人口をダウンロード
https://www.city.osaka.lg.jp/toshikeikaku/page/0000015211.html
クレンジング方法：全体に対する対象年齢層の割合、男性の割合を区ごとに算出しフィールド化

#### 「4.本庁舎までの時間」
データ入手・クレンジング方法：下記乗換サイトで、旧区中心駅～新区役所最寄駅までの所要時間を、土曜8, 12, 15, 18時で求め平均
https://ekitan.com/

#### 「5.平均世帯人数」
データ入手方法：下記大阪市HPより2020年の1世帯当たり人員をダウンロード
https://www.city.osaka.lg.jp/toshikeikaku/page/0000068035.html
クレンジング方法：区ごとににフィールド化

#### 「6.平均1人あたり所得」
データ入手方法：下記大阪市HPより世帯別収入階層データをダウンロード
https://www.city.osaka.lg.jp/shimin/cmsfiles/contents/0000180/180789/20.pdf
クレンジング方法：階層ごとに平均年収と割合を積算して平均世帯年収を求め、5.で求めた平均世帯人員で割る

### 得られたデータ
下図のように、フィールド化されたデータが得られました。これでクレンジング完了です
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/10914876-1d9d-fdb2-0662-ff7e06205fff.png)

日本語だとPythonで扱った際に良くない事が色々起こるので、英語化しておきます
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/7322a82b-5f3a-dea1-7c75-6fe179681633.png)

## 4. データの分析と可視化
今回は「目的変数：賛成率」、「説明変数：前章の6種類の特徴量」として回帰をしたいので、
下記の手順で進めていきます
4-1) データ全体の可視化
4-2) 特徴量選択（相関の強い特徴量の除外）
4-3) 回帰モデルの作成
4-4) 性能評価
4-5) モデルの改善

### 4-1. データ全体の可視化
[過去に作成したこちらのツール](https://qiita.com/c60evaporator/items/20f11b6ee965cec48570)を使用し、散布図と相関係数を求めます。

```python
#%%読込
import pandas as pd
from custom_pair_plot import CustomPairPlot

#使用するフィールド
KEY_VALUE = 'ward_before'#キー列
OBJECTIVE_VARIALBLE = 'approval_rate'#目的変数
EXPLANATORY_VALIABLES = ['1_over60','2_between_30to60','3_male_ratio','4_required_time','5_household_member','6_income']#説明変数

#pair_analyzerで相関を可視化
df = pd.read_csv(f'./osaka_metropolis_english.csv')
use_cols = [OBJECTIVE_VARIALBLE] + EXPLANATORY_VALIABLES
gp = CustomPairPlot()
gp.pairanalyzer(df[use_cols])
```
※特徴量名を英語に変えていますが、番号から類推してください
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/390eb1c2-461e-e0ea-8711-506ea9bb6b3e.png)
赤枠で囲った部分が、目的変数（賛成率）との相関係数になります。

#### 可視化により得られた知見
下記のような知見が得られそうです
①20代の逆転を考慮しない特徴量「1. 60歳以上割合」よりも、考慮した「2. 30～50代比率」の方が賛成率との相関係数が高い（もくろみ通り）
②男性比率は賛成率との相関係数が低い（外れ値(西成区)の影響がありそうだが‥）
③説明変数「2. 30～50代比率」と「6. 平均年収」の相関係数が0.89と高い

### 4-2. 特徴量選択
天下り的ですが、今回は最近の回帰分析の定番手法である、勾配ブースティング（XGBoost）を使用します。
この手法は多重共線性（相関の強い説明変数を両方使用した際の不具合）に比較的強いと言われていますが、[それでもこちらのサイトのように相関の強い変数をそのまま使用すると、無視できない影響が出てくるようです](https://note.com/hetiyaborake/n/n818e4af6d297)

そこで今回は、他の説明変数と相関係数が高い説明変数を、分析対象から除外しました。
多重共線性の基準としてはVIF=10(R=0.95に相当)、あるいは5(R=0.9に相当)がよく使われるようですが、
今回は厳しめに見て、相関係数0.9以上の説明変数を除外します。
具体的には
**・「1. 60歳以上割合」は除外**‥「2. 30～50代比率」と相関係数-0.96で、かつ賛成率との相関が2.より低い
**・「6. 平均1人あたり年収は除外」**「2. 30～50代比率」との相関係数0.89と基準をやや下回っているが、「働き盛りが多いと年収が高い」という明確な因果関係が特定できるため、除外する

というわけで、今後の分析には下記の4特徴量を使用することとします。
**特徴量2「30～50代割合」**(20代の逆転現象を考慮)
**特徴量3「男性比率」**
**特徴量4「本庁舎が置かれる駅から、旧区中心駅までの所要時間」**
**特徴量5「平均世帯人数」**

特徴量除外後の散布図は下図のようになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/43ed34c5-24c7-2e2a-bb6a-42dec2283988.png)
この説明変数を使用し、回帰分析を進めていきます。
また、以下は回帰分析であることを強調するため、「特徴量」の用語を「説明変数」に統一させて頂きます。

### 4-3. 回帰モデルの作成
前述のように、最近の回帰分析の主流である勾配ブースティング（XGBoost）を使用します。
便利なPythonのライブラリが存在するので、カンタンに使用することができますが、
ハイパーパラメータがいくつか存在するので、グリッドサーチにより最適化します

#### グリッドサーチによるパラメータの最適化
グリッドサーチとクロスバリデーションを組み合わせて、最適なハイパーパラメータを探します。

クロスバリデーションとは、端的に言うと
「データをN等分して、性能と過学習防止のバランスが取れたパラメータ設定を探す」
という操作になります（データの等分方法は下図Wikipedia画像参照）
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/fe534e49-0559-d395-5a60-e1ee0ebb071c.png)

またグリッドサーチとは、あらかじめ決め打ちしたパラメータの組み合わせを総当たりで探索する方法です。時間が掛かることがデメリット、実装が簡単であることがメリットとなります。
ネットを検索し、XGBoostで実際によく使われているパラメータの範囲を候補として使用しました。

```python
import xgboost as xgb
from sklearn import metrics as met
import sklearn as skl
from sklearn.model_selection import LeaveOneOut
from sklearn.model_selection import train_test_split
from sklearn.model_selection import GridSearchCV
import numpy as np
import matplotlib.pyplot as plt
import os
from datetime import datetime
import seaborn as sns

#使用するフィールド
KEY_VALUE = 'ward_before'#キー列
OBJECTIVE_VARIALBLE = 'approval_rate'#目的変数
USE_EXPLANATORY = ['2_between_30to60','3_male_ratio','4_required_time','5_household_member']#使用する説明変数
#データ読込
df = pd.read_csv(f'./osaka_metropolis_english.csv')

#目的変数と説明変数を取得（pandasではなくndarrayに変換）
y = df[[OBJECTIVE_VARIALBLE]].values
X = df[USE_EXPLANATORY].values
#グリッドサーチと性能評価の共通パラメータ
num_round=10000#最大学習回数
early_stopping_rounds=50#評価指標がこの回数連続で改善しなくなった時点で学習をストップ
seed = 42#乱数シード

#グリッドサーチ用パラメータ(詳細はhttps://qiita.com/R1ck29/items/50ba7fa5afa49e334a8f)
cv_params = {'eval_metric':['rmse'],#データの評価指標
             'objective':['reg:squarederror'],#最小化させるべき損失関数
             'random_state':[seed],#乱数シード
             'booster': ['gbtree'],
             'learning_rate':[0.1,0.3,0.5],
             'min_child_weight':[1,5,15],
             'max_depth':[3,5,7],
             'colsample_bytree':[0.5,0.8,1.0],
             'subsample':[0.5,0.8,1.0]
            }

#XGBoostのインスタンス作成
cv_model = xgb.XGBRegressor()
#グリッドサーチのインスタンス作成
# n_jobs=-1にするとCPU100%で全コア並列計算。とても速い。
cv = GridSearchCV(cv_model, cv_params, cv = 5, scoring= 'r2', n_jobs =-1)

#グリッドサーチ実行
evallist = [(X, y)]
cv.fit(X,
        y,
        eval_set=evallist,
        early_stopping_rounds=early_stopping_rounds
        )

#最適パラメータの表示
print('最適パラメータ ' + str(cv.best_params_))
print('変数重要度' + str(cv.best_estimator_.feature_importances_))
```
なお、ここでグリッドサーチに使用するデータとテストデータを分けていない事に違和感を覚えられる方もいるかもしれませんが、本節ではあくまでパラメータの最適化が目的で、性能評価は次節で別途行うため、ご容赦いただければと思います。
（Leave_One_Outで1区ごとにグリッドサーチを回すと所要時間が‥）

### 4-4. 性能評価
テストデータと学習データを分けて性能評価します。

指標としては、下記を使用します
**RMSE平均**：全体としての予測誤差の大きさ(小さいほどGood)
**予測と実際の相関係数**：予測と実際の相関の大きさ(大きいほどGood)
**予測誤差の最大値**：予測誤差が極端に大きい区が存在しないか評価(小さいほどGood)

#### テストデータの分割方法
分け方には何種類かありますが

・データ数が少ないので、学習データ数を確保する必要がある
・1区ごとに性能指標を算出し、モデルに適合しない区をあぶりだしたい

という観点から、
1個ずつテストデータを取り出し、残りを学習データとする「**Leave_One_Out**」で性能評価を行います
（他の区全てのデータで学習したモデルで性能を評価する、Wikipediaの下図参照）
![450px-LOOCV.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/5077fc02-f12d-7859-7092-8a27b2bd66b6.gif)

### 評価に使用したコード
前記パラメータ最適化用のコードに、下記のコードを追記して性能評価します。

```python
#%%3. 性能評価(Leave-One-Out)
#パラメータにグリッドサーチでの最適パラメータを使用
params = cv.best_params_
#結果保持用のDataFrame
df_result = pd.DataFrame(columns=['test_index','eval_rmse_min','train_rmse_min','num_train'])

#Leave-One-Outでデータ分割して性能評価
loo = LeaveOneOut()
for train_index, test_index in loo.split(X):#全データに対して分割ループ
    X_train, X_test = X[train_index], X[test_index]
    y_train, y_test = y[train_index], y[test_index]
    dtrain = xgb.DMatrix(X_train, label=y_train)#学習用データ
    dtest = xgb.DMatrix(X_test, label=y_test)#テストデータ
    evals = [(dtest, 'eval'),(dtrain, 'train')]#結果表示用の学習データとテストデータを指定
    evals_result = {}#結果保持用

    #学習実行
    model = xgb.train(params,
                    dtrain,#訓練データ
                    num_boost_round=num_round,
                    early_stopping_rounds=early_stopping_rounds,
                    evals=evals,
                    evals_result=evals_result
                    )
    
    #モデルの性能評価
    test_pred = model.predict(dtest, ntree_limit=model.best_ntree_limit)
    num_train = len(evals_result['eval']['rmse'])
    df_result = df_result.append({'test_index': test_index[0],
                    'key_value': df[[KEY_VALUE]].iloc[test_index[0],0],
                    'pred_value': test_pred[0],
                    'real_value': df[[OBJECTIVE_VARIALBLE]].iloc[test_index[0],0],
                    'eval_rmse_min': evals_result['eval']['rmse'][num_train - 1],
                    'train_rmse_min': evals_result['train']['rmse'][num_train - 1],
                    'num_train': num_train},
                    ignore_index=True)

#性能評価結果の表示
print('RMSE平均' + str(df_result['eval_rmse_min'].mean()))
print('相関係数' + str(df_result[['pred_value','real_value']].corr().iloc[1,0]))
print('予測誤差の最大値' + str(max((df_result['pred_value'] - df_result['real_value']).abs())))

#結果を出力
dt_now = datetime.now().strftime('%Y%m%d%H%M%S')
feat_use = 'feat' + '-'.join([ex.split('_')[0] for ex in USE_EXPLANATORY])
#評価結果
df_result.to_csv(f"{os.getenv('HOMEDRIVE')}{os.getenv('HOMEPATH')}\Desktop\{feat_use}_{dt_now}_result.csv")

path = f"{os.getenv('HOMEDRIVE')}{os.getenv('HOMEPATH')}\Desktop\{feat_use}_{dt_now}_result.txt"
with open(path, mode='w') as f:
        f.write('特徴量' + str(USE_EXPLANATORY))
        f.write('\n最適パラメータ' + str(cv.best_params_))
        f.write('\nグリッドサーチ対象' + str(cv_params))
        f.write('\n変数重要度' + str(cv.best_estimator_.feature_importances_))
        f.write('\nRMSE平均' + str(df_result['eval_rmse_min'].mean()))
        f.write('\n相関係数' + str(df_result[['pred_value','real_value']].corr().iloc[1,0]))
        f.write('\n予測誤差の最大値' + str(max((df_result['pred_value'] - df_result['real_value']).abs())))

#散布図表示
sns.regplot(x="pred_value", y="real_value", data=df_result, ci=0)
```

#### 評価結果
今回求めた性能評価指標は下記のようになりました。

RMSE平均(小さいほどGood)：0.0226
予測と実際の相関係数(大きいほどGood)：0.619
予測誤差の最大値(小さいほどGood)：0.0709

この指標が良いのか悪いのかは次節で評価するので、
**まずは予測値(横軸)と実際の値(縦軸)とを比較した散布図を見てみます**
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/c512aca2-5d94-4b15-e11f-482a78f9a7d8.png)

上記散布図を見て、予測値と実測値の差の大きい福島区の結果に、少し違和感を覚えました。
例えば、同じ合併後の北区に所属する都島区と比較すると、
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/63adb054-bc43-44d7-04ee-10f8690deb8d.png)

全ての特徴量の値が近いにも関わらず福島区の予測値が不自然に低くなっており(むしろ30～50代割合が大きい分、感覚的には都島区よりも予測値が大きくなりそう)、
Leave-One-Out時に分割した学習データに過学習してしまっていると思われます。
次節以降で過学習防止法を考えてみます。

### 4-5. モデルの改善
前節で過学習と思われる現象が発生していたため、改善策を適用してモデルの性能向上を図ります
#### 改善策1 特徴量を絞る
過学習防止のためグリッドサーチ以上の方法が思いつかなかったので、使用する説明変数を絞って対処してみました。

|使用説明変数|RMSE平均<br>(小さいほどGood)|予測と実際の相関係数<br>(大きいほどGood)|予測誤差の最大値<br> (小さいほどGood)|
|---|---|---|---|
|2,3,4,5全て使用|0.0226|0.619|0.0709|
|2,3,4|0.0213|0.648|0.0531|
|**2,3,5**|**0.0191**|**0.744**|**0.0474**|
|2,4,5|0.0222|0.614|0.0739|
|3,4,5|0.0221|0.642|0.0593|
|2,3|0.0197|0.707|0.0519|
|**2,4**|**0.0181**|**0.733**|**0.0604**|
|2,5|0.0248|0.572|0.0590|
|3,4|0.0229|0.568|0.0610|
|3,5|0.0232|0.596|0.0503|
|4,5|0.0304|0.460|0.0750|
|2|0.0209|0.646|0.0559|
|3|0.0303|-0.996|0.0726|
|4|0.0246|0.541|0.0532|
|5|0.0292|0.407|0.0686|

特に性能が高かった組み合わせを太字にしました。
先ほどと同様に散布図も下記します

##### 説明変数2,3,5を使用した場合の予測値vs実測値散布図
![pred_real_2-3-5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/7eedacf6-7ad8-ae43-3b0c-85a1a233144e.png)

##### 説明変数2,4を使用した場合の予測値vs実測値散布図
![pred_real_2-4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/91ec791c-b266-b419-99ab-c79add8c20da.png)

#### なぜ指標が改善した？
使用した4説明変数間には3者間以上での相関・因果関係（例：「本庁舎に近い中心街ほど、30～50代の働き盛り世代が多く、ファミリー層が少ないので平均世帯人員が少ない」など）が含まれるため、全てを使用すると悪影響を与えるのではないかと仮説を立てています。

また、正直私の勾配ブースティング活用スキルが足りないため、4変数では性能を出し切れなかった感も否めないため、
**「こうすると過学習を防げる！」というような知見をお持ちの方がいれば**、コメント頂けるとありがたいです。

#### 改善策2　地域の位置関係の考慮
性能が良かった説明変数2,3,5を使用したモデルの、予測値と実測値の散布図を見ると、予測より実測値が低い区には、南側の区(特に赤で示した新天王寺区に相当する区)が目立つことが分かります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/cc3ca275-e8f4-146e-e068-4038fb251366.png)
また逆に北区や淀川区など緑で示した北側の区では、予測よりも実測値が高い区が多いです。

大阪市は北側と南側で地域制が大きく異なると言われており、この地域性の差が投票行動に影響を与えたと仮説が立てられます。

そこで今回、南北を表す**説明変数として、「緯度(latitude)」を追加しました**
(「北緯〇度△分□秒」の〇の部分は全て"34度"なので、今回は分以下の値を使用します)

|使用説明変数|RMSE平均<br>(小さいほどGood)|予測と実際の相関係数<br>(大きいほどGood)|予測誤差の最大値<br> (小さいほどGood)|
|---|---|---|---|
|2,3,5|0.0191|0.744|0.0474|
|**2,3,5+緯度**|**0.0174**|**0.767**|**0.0444**|
|2,4|0.0181|0.733|0.0604|
|**2,4+緯度**|**0.0182**|**0.755**|**0.0586**|

説明変数「2,3,5」を使用した場合も、「2,4」を使用した場合も、軒並み指標が改善しております。
やはり南北の差が賛成率の要因となっているようです。

##### 説明変数2,3,5,latitudeを使用した場合の予測値vs実測値散布図
![pred_real_2-3-5-latitude.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/efb97412-7cc3-6900-6f38-8119085cc080.png)

##### 説明変数2,4,latitudeを使用した場合の予測値vs実測値散布図
![pred_real_2-4-latitude.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/810b2e84-37da-015b-e422-2eaa0a028164.png)

##### 説明変数2,3,5,latitudeを使用した場合の特徴量重要度(Feature Importance)
![importance_2-3-5-latitude.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/d558219c-299b-234c-607c-1d74f7db5d84.png)

## 5. 考察とまとめ
最も性能が良いモデル(説明変数2,3,5,latitudeを使用)は、重要度の大きい順に、下記の特徴量を使用しています
特徴量2 **「30～50代割合」**
特徴量5 **「平均世帯人数」**
追加特徴量 **「緯度」**
特徴量3 **「男性比率」**

ここで使用した特徴量が、目的変数である賛成率に影響を与える要因とみなせます。

深い原因考察は避けますが、
**「働き盛り世代の男性比率が大きく、平均世帯人数が少なく、北側にある区」で賛成が多く**
**「働き盛り世代の男性比率が小さく、ファミリー層が多く、南側にある区」で反対が多い**
ことが事実として言えそうです。
（前節の特徴量重要度を見る限り、年齢構成の影響が最も大きく、男性比率の影響が最も小さそうです）

また、4-1での相関分析からもわかるように
**20代の反対が比較的多い(30代との逆転現象)ことも、区による賛成率の差に影響を与えている**
とも言えるかと思います。

## 6.感じたこと
今回の分析を通じて、下記のような感触を得ることができました
### ドメイン知識の重要さ
4.4で過学習を疑うきっかけとして、「福島区が外れ値となる理由が思い浮かばない」というものがあります。
福島区はグルメで有名な特色のある地域ですが、
投票に影響を与えるような要素としては市内でも平均的なイメージ（ドメイン知識）があったので、「ここが外れ値になるということは、作ったモデルの性能が悪いに違いない」という考えに至ることができました。
### ドメイン知識を過信しすぎない
2章では、市民の感覚として感じていた「本庁舎から遠いほど不便さを感じて賛成率が低そう」という考えに基づいて、
特徴量「4.本庁舎までの所要時間」を追加しました。

乗換案内をポチポチしながらデータ収集しているときは、
**「性能アップに繋がる特徴量になるに違いない！ドメイン知識最強！」**
と思いながらハイテンションになっていましたが、
実際に分析を進めてみると、他の特徴量に埋もれてしまう程度の影響度しかありませんでした。

そして頑張って収集したデータには思い入れも生まれるので、
**「何としてもこの特徴量は捨てたくない！無理やりでも使いたい」**
という、分析者として客観性を欠いた心情が生じてきます。

このようなドメイン知識から生まれた先入観や思い入れを封印して、得られた事実のみから判断する事の難しさを、身をもって感じることができました。

### 生きたデータを分析する面白さ
タイムリーかつ話題性の高いデータを扱うと、モチベーションが上がることをひしひしと感じました。
今後も勉強した結果をアウトプットする場として、時事ネタを分析していこうと思います！

# おまけ：前回投票との差の分析
おまけ程度の分析となりますが、前回投票との差が大きい区についても分析してみます。

### 前回投票と今回投票の賛成率比較
前回2015年投票の賛成率と、今回賛成率を比較してみました(横軸：前回、縦軸：今回)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/b4e95ca1-e109-4385-983c-4b2262039e2b.png)
概ね前回と同様の傾向(R2=0.88)ですが、いくつか前回から変化のあった区があります（特に顕著な区が赤色）

### 区分けの比較
前述の前回からの変化の一因に、合併後の区分けの変化があるかもしれないと思い、
前回と今回で区分けを比較してみました。

#### 前回の区分け
下図は産経新聞HPより
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/c514a52c-0aeb-288e-29f3-a3ee89c8bc34.png)

#### 今回の区分け
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/f4d0d6a7-2ee8-3fc3-f4f2-671c8d456aae.png)

中央区、北区の中心部の区分けには変化がなく、
周辺部の区（前回の東区、南区、湾岸区が相当）の区分けが大きく変化が出ている（本庁舎が市中心部に近い位置にある区に編入）ようです。

### 分析方針
前回から2%以上の増減があった区を抽出してみました。
前節の分類（中央区・北区 or その他）に基づき、合併後の区分けの変化も併記します。

|区名|合併後の区分けの変化|増減|
|---|---|---|
|大正区|あり|+2.54%|
|西成区|なし|+2.49%|
|北区|なし|-2.73%|
|中央区|あり|-3.12%|
|阿倍野区|あり|-3.32%|
|西区|なし|-3.89%|
|港区|あり|-4.82%|

合併後の区分けの変化があった区と、なかった区に分けて分析したいと思います。

### 1. 区分けに変化があった区
港区、阿倍野区、大正区が相当します。
天下り的にはなりますが、港区とそれ以外で異なる要因がありそうだと感じています。

#### 港区
減少幅がダントツで大きい港区ですが、[各種報道でも触れられているように](https://www.news-postseven.com/archives/20201008_1601456.html?DETAIL)、投票前から地理的な遠さを指摘する声があったようです。
合併後の地図を見ると、港区は本庁舎のある十三駅から、非常に離れた場所にあります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/39c875b5-4592-3f42-2889-e9c058723a0b.png)
ただ離れた場所にあるだけでなく、間に淀川という大きな川を挟んでいて地域的な繋がりも小さいですし、
電車で行こうとしても途中の梅田駅での乗換時間が非常に長く、とても不便です。
（東京に例えると、お台場(江東区)と新小岩(葛飾区)くらいの位置関係のイメージ）

前回構想では港区に本庁舎が置かれる予定であったことや、
似たような位置関係にある此花区でも1.5%賛成率が低下していることからも、
**本庁舎までの地理的な遠さ**が影響しているのではと推察されます。

#### 阿倍野区、大正区
（多少の主観が入ってしまうことをご容赦ください）
阿倍野区に関しては、前回の構想では阿倍野区には新"南区"の本庁舎が置かれる予定だったのが、今回構想では旧天王寺区に本庁舎が置かれ、新区の名称も"天王寺区"となります。
[「あべのハルカス」という名前がついたいきさつが、「阿倍野という名前を全国区にしたい」](https://ja.wikipedia.org/wiki/%E3%81%82%E3%81%B9%E3%81%AE%E3%83%8F%E3%83%AB%E3%82%AB%E3%82%B9#%E6%A6%82%E8%AA%AC)
というものだったことから考えても、ある種のライバル関係にある隣の天王寺区に区庁舎も区名も移ってしまうことは、区民にとってデメリットに映ってしまったのでは、と感じます。

逆に大正区は、前回の本庁舎が隣で規模の近い港区であったのに対し、今回はより規模の大きい中央区に移るので、
**隣接した規模の近い区に本庁舎が移る**という状況を脱したことで、賛成率が上がったのではと推察されます

### 2. 区分けに変化がなかった区
西成区(賛成率増)、北区(賛成率減)、中央区(減)、西区(減)が相当します。

あまりはっきりとしたことは分かりませんが、区分けの変更とは無関係に、この5年間で何らかの政治的変化が発生したと思われます。
（賛成率の減った北区、中央区、西区が平均1人あたり年収のトップ4に入る区であるのに対し、
賛成率の増えた西成区が1人あたり年収の最も低い区であることから、
収入に関係した変化なのでは？という仮説が立てられます）

## 前回と今回で賛成率が変化した区　まとめ
下記の3分類に分けられそうです
**区役所との地理的隔絶**
港区、(此花区)

**近隣の規模の近い区に本庁舎**
阿倍野区、大正区、(旭区)

**区分けの変更とは無関係に、何らかの変化が発生（収入と関係ありそう？）**
西成区、中央区、西区、北区