都議選当選への影響要因をロジスティック回帰で推定してみた！【Rで実践編】

前回は、都議選の当選・落選という結果が候補者のどんな属性に影響を受けているのか、Pythonで推定を行いました。
また必要なデータは、PythonのBeautiful Soupというライブラリを使ってスクレイピングして取得しました。
randpy.hatenablog.com

ロジスティック回帰に関心が高いのか、扱ったデータが都議選のデータだったからか分かりませんが、これまで公開している記事の中だと最もアクセスがありました！
皆さまありがとうございます！

さて、「Python、Pythonばっか言ってるけどさー。こちとらRで分析したいわけよ！ちょっとそこら辺考慮してくれませんかねー？」
というR~~おじさん~~ユーザーからの声が聞こえてきそうだったので、
今回は上記のスクレイピング＆ロジスティック回帰をRを使って実装してみます！

是非参考にしてみてください。
なお、これから解説するコードはGithub上で公開しています。基本的にはそのままコピペすれば実行が完了するようになっています。（あまり綺麗なコードじゃないですが、怒らないでください…。）
https://github.com/chan-ume/npurProject/blob/master/glm/logit_for_tokyo_assembly_elections.R

あ、実装の話に入る前に...。ロジスティック回帰って何？という方は、下の理論記事ををまずはご覧になっていただけると！
randpy.hatenablog.com

できるだけ直観的な理解ができるように頑張って書いたので、あまり統計に馴染みが無い方にとっても分かりやすいのではと思います。

まずは必要なデータをスクレイピング

都議選の候補者の属性による、当選確率への影響を分析したいので、まずは候補者の属性データをスクレイピングしてきます。

Python編の記事でも書きましたが、できるだけ多くの属性データを取得したかったのですが...。
あまりWeb上には載っていなかったりして、下のリンクから年齢・所属・当選回数を取得することにしました。
党派別立候補者数:都議選２０１７:東京新聞(TOKYO Web)

変数としてはかなり少ないですが、ロジスティック回帰を含む一般化線形モデルの分析をRで実践してみることが目的なので、まあ許してください...。

まずは必要なライブラリを使えるように宣言しましょう。もしまだインストールしていない場合は、2行目から6行目のコメントアウトを外してインストールしてください。

次に、先ほどのURLから、各選挙区のURLを取得してきます。
今回はxpathにて指定しました。

そして、取得した各選挙区のURLからfor文を使って、それぞれの候補者情報を取得します。

掲載されている情報が下の画像のようなので、名前と年齢が一緒に取得されてしまいます。
f:id:Np-Ur:20170716141903p:plain

そこで44行目で名前と年齢を分解して、行列に加えています。名字と名前の間にスペースが入っているので削除の処理も加えています。

この辺りの正規表現を使った文字列処理、毎回忘れて調べながらやっています…。
正規表現の実力はレベル1ぐらいなので、「もっと良い書き方あるよ！」とかあったら共有してもらえると助かります。

さて、head関数を使ってどんな行列が作られたのか確認してみましょう。

> head(candidateInfo)
         nameInfoList       partyInfoList newOrOldInfoList winCountInfoList                
[1,] "0" "樋口　高顕（34）" "都"          "新"             "0"              "樋口高顕" "34"
[2,] "0" "須賀　和男（61）" "無"          "新"             "0"              "須賀和男" "61"
[3,] "0" "中村　彩（27）"   "自"          "新"             "0"              "中村彩"   "27"
[4,] "0" "後藤　輝樹（34）" "諸"          "新"             "0"              "後藤輝樹" "34"
[5,] "0" "森山　高至（51）" "無"          "新"             "0"              "森山高至" "51"
[6,] "0" "立石　晴康（75）" "無"          "現"             "8"              "立石晴康" "75"

こんな感じです。
一番左の列は全候補者に「0」を代入しています。このあと、当選した候補者に関しては「1」を代入して、当選者と落選者を区別する処理を加えます。
そのための準備と思ってください。

左から

当選 or 落選を代入する列
名前と年齢
政党
現職か新人か
過去の当選回数
名前
年齢

となっています。

続いて当選 or 落選のデータをスクレイピング

次に、各候補者が当選したか落選したかのデータを取得してきます。

色々探しましたが、下のサイトがわりとスクレイピングしやすかったです。
開票結果 - 2017都議選：朝日新聞デジタル

ここから当選と落選を取得するスクレイピングコードですが、以下のように書いてみました。

48行目で各選挙区の当選情報が掲載されているURlを取得しています。先ほどの処理と一緒ですね。先ほどは2行に分けたのを1行で今度はまとめました。

下の画像のように当選したか落選したかが表現されています。目印は左の画像ですね！
f:id:Np-Ur:20170716143854p:plain

58行目で画像が貼ってある部分の情報を取得しています。
その後66行目で画像が貼ってある候補者のみ抽出する処理を加えています。

その後、下のコードで、先ほど作った候補者情報の一番左の列に、当選している場合は「1」を代入しています。

これにて、候補者の属性データと当選しかた落選したかの情報をすべて取得することに成功しました！
皆さまお疲れ様です！

どんなデータが取得できているのか、下を実行して確認してみましょう！

>candidateInfo

ロジスティック回帰を実行

さて、スクレイピング処理が終わったので、ロジスティック回帰を行いましょうか。

Rでロジスティック回帰を行うのはとても簡単です。
randpy.hatenablog.com
こちらのコンテンツでも少し触れましたが、「glm」という関数を用います。

glmは一般化線形モデルを扱うための関数で、目的変数と説明変数、そしてどんな確率分布・リンク関数を用いるかを指定してあげればOKです。
今回は、当選したか落選したか（0 or 1）を目的変数、候補者の年齢・政党・当選回数・現職か新人かが説明変数ですね。

次にロジスティック回帰なので、確率分布はベルヌーイ分布、リンク関数はロジットを指定しましょう。

resultという変数にglmの結果が入りました。
次に、summaryで結果を確認しましょう。

> summary(result)

Call:
glm(formula = p ~ age + numberOfWinning + candidateInfo[, 3] + 
    candidateInfo[, 4], family = binomial(link = "logit"))

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.76213  -0.49756  -0.00018   0.24564   2.38595  

Coefficients:
                       Estimate Std. Error z value Pr(>|z|)    
(Intercept)            -0.12160    2.02110  -0.060 0.952025    
age                    -0.02722    0.01984  -1.372 0.170003    
numberOfWinning        -0.02610    0.16762  -0.156 0.876278    
candidateInfo[, 3]維    0.21076    1.91554   0.110 0.912387    
candidateInfo[, 3]共    2.02264    1.41251   1.432 0.152158    
candidateInfo[, 3]公    4.22457    1.71752   2.460 0.013906 *  
candidateInfo[, 3]自   -0.10528    1.34462  -0.078 0.937591    
candidateInfo[, 3]社  -14.80310 3956.18057  -0.004 0.997015    
candidateInfo[, 3]諸  -15.04794  950.25507  -0.016 0.987365    
candidateInfo[, 3]都    6.19448    1.72790   3.585 0.000337 ***
candidateInfo[, 3]民   -0.20480    1.47764  -0.139 0.889767    
candidateInfo[, 3]無   -0.39719    1.43023  -0.278 0.781235    
candidateInfo[, 4]現    1.64379    0.92147   1.784 0.074442 .  
candidateInfo[, 4]新   -1.19858    0.98608  -1.215 0.224177

はいー！やっぱり都民ファースト強いー！！ということが分かりました。
あとは現職の方が有利ということも分かりました。

なるほど、まあそうか......。という結果ですね。

係数に関しての理解についてはPython編をご覧くださいませ。
Python編と結果が異なるのは、最尤法のアルゴリズムの違いでしょうかね。

この辺りも考察するべきなんでしょうけど、今回は割愛します。

AICを基に変数選択

結果を求めることができましたが、せっかくなのでAICを用いた変数選択の手法について簡単に触れておきます。

AICとは、尤度を$L$、説明変数の数を$k$としたときに
$$AIC = -2\log L + 2k$$
で求めることができる値です。

これは式を見てもらうと分かるように、尤度が大きいほど、そして説明変数の数が少ないほどAICは小さくなります。

説明変数の数が多くなると、尤度は上がりますが、その分新しいデータが追加されたときに当てはまりが悪くなる危険があります。
極端な話、サンプルサイズが100のデータに対して説明変数を100個用意してしまえば、当てはまりは完璧です。

しかし、このパラメータの推定値は、用いたサンプルサイズ100のデータにのみ当てはまりがよいもので、この推定値で別のデータを予測することができるかいうと、恐らくほとんど使えない結果になるでしょう。

例えばとても旦那さん思いの奥さんがいたとします。その旦那さんはとても激辛好き&お肉好きのため、通常のカレー（にんじん・じゃがいも・たまねぎ・豚肉）に唐辛子を沢山入れて、更に肉をドロドロ追加で入れて専用カレーを作ったとします。このカレーはその旦那さんにとっては100点ですが、お隣さんにおすそ分けして食べてもらったら、「辛すぎる！肉多すぎ！もう0点！」となるかもしれません。

そんなカレーよりは、にんじん・じゃがいも・たまねぎ・豚肉のみが入っていて誰が食べてもまあまあ美味しいカレーの方がいいですよね？（この例分かりやすいのか？笑）

そのあたりをうまいこと表現したのがAICというわけです。このAICが最も小さくなるように変数を選べば、それが最も良いモデルだろうと考えます。

ちなみにAICは、偉大なる日本の研究者である赤池先生が考案しました。

RでAICを実践

RでAICを用いた変数選択を実践するには、stepという関数を使います。

>step(result)

はい、これだけです。色々AICの値を計算しながら、最もAICが小さいモデルが何なのかを教えてくれます。

出力された結果の最後の方を見てみると

Call:  glm(formula = p ~ age + candidateInfo[, 3] + candidateInfo[, 
    4], family = binomial(link = "logit"))

こんな風に最適なモデルはこちらですよ！と教えてくれています。

では、今度はこちらのモデルのsummaryを見てみましょう。

> summary(glm(p ~ age + candidateInfo[,3] + candidateInfo[,4], family=binomial(link="logit")))

Call:
glm(formula = p ~ age + candidateInfo[, 3] + candidateInfo[, 
    4], family = binomial(link = "logit"))

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.77792  -0.49464  -0.00018   0.24663   2.38587  

Coefficients:
                       Estimate Std. Error z value Pr(>|z|)    
(Intercept)            -0.06537    1.98740  -0.033 0.973761    
age                    -0.02841    0.01833  -1.550 0.121100    
candidateInfo[, 3]維    0.18138    1.90515   0.095 0.924151    
candidateInfo[, 3]共    1.99148    1.39607   1.426 0.153728    
candidateInfo[, 3]公    4.18779    1.70029   2.463 0.013779 *  
candidateInfo[, 3]自   -0.13377    1.33058  -0.101 0.919917    
candidateInfo[, 3]社  -14.82824 3956.18057  -0.004 0.997009    
candidateInfo[, 3]諸  -15.08263  949.44107  -0.016 0.987326    
candidateInfo[, 3]都    6.15882    1.71044   3.601 0.000317 ***
candidateInfo[, 3]民   -0.24014    1.45861  -0.165 0.869228    
candidateInfo[, 3]無   -0.44055    1.40164  -0.314 0.753286    
candidateInfo[, 4]現    1.61823    0.90571   1.787 0.073987 .  
candidateInfo[, 4]新   -1.16666    0.96399  -1.210 0.226189    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

まあそんな変わらないですね。
元々の変数が少ないので、あまりやる価値はないかもしれません。

今回の場合は、どの変数を選択するか？というよりも、もっと変数を増やすことを考えた方が良いでしょう。