データの不備を統計的に見抜く (Gelman's Secret Weapon)

リクルートの高柳さん、Yahooの簑田さんと共同で翻訳した本が出版されます。
「みんなのR」(原題：R for Everyone)です。

みんなのR -データ分析と統計解析の新しい教科書-

作者: Jared P. Lander,Tokyo.R（協力）,高柳慎一,牧山幸史,簑田高志
出版社/メーカー: マイナビ
発売日: 2015/06/30
メディア: 単行本（ソフトカバー）
この商品を含むブログ (7件) を見る

この本は、統計言語 R のインストール・基本的な使い方から始まり、統計解析の基礎からちょっと高度な話題まで、幅広く取り扱っています。

特徴としては、

RStudio の使用を推奨
グラフィクスはすべて ggplot2 を使用
plyr, data.table, stringr といった、モダンな便利パッケージを使用*1
説明に使用されるデータはすべて Web からダウンロード可能
R の個々の分析関数に対して、使う際の細かい TIPS が数多くちりばめられている

といったことが挙げられます。

原著者の Jared P. Lander さんは、coefplot パッケージ、useful パッケージの作者です。
この coefplot パッケージは、作成した統計モデルの係数を可視化するのに非常に便利なパッケージです。

今日は「みんなのR」19章の話題から、coefplot パッケージを使って、取得したデータに問題があることを統計的に見抜くという話を紹介しようと思います。

データ

今回使用するデータは、アメリカ大統領選挙のデータです。
このデータは、4 年に一度の大統領選挙に対して、1948 年から 2000 年に渡って、のべ 15240 人の有権者から取ったアンケート結果です。
データは著者の Web サイトから取得することができます。

connection <- url("http://jaredlander.com/data/ideo.rdata")
load(connection)
close(connection)
head(ideo)

  Year       Vote Age Gender  Race                             Education               Income                  Religion
1 1948   democrat  NA   male white     grade school of less (0-8 grades)  34 to 67 percentile                protestant
2 1948 republican  NA female white high school (12 grades or fewer, incl 96 to 100 percentile                protestant
3 1948   democrat  NA female white high school (12 grades or fewer, incl  68 to 95 percentile catholic (roman catholic)
4 1948 republican  NA female white some college(13 grades or more,but no 96 to 100 percentile                protestant
5 1948   democrat  NA   male white some college(13 grades or more,but no  68 to 95 percentile catholic (roman catholic)
6 1948 republican  NA female white high school (12 grades or fewer, incl 96 to 100 percentile                protestant

選挙年(Year)
民主党(democrat)と共和党(republican)のどちらに投票したか(Vote)
年齢(Age)
性別(Gender)
人種(Race)
教育(Education)
収入(Income)
宗教(Religion)

という変数が格納されています。

coefplot による係数データの可視化

民主党と共和党のどちらに投票したか(Vote)を出力変数、人種(Race)・収入(Income)・性別(Gender)・教育(Education)を入力変数として、ロジスティック回帰モデルにあてはめてみます。

model <- glm(Vote ~ Race + Income + Gender + Education, data=ideo, family=binomial(link="logit"))
coef(summary(model))

                                                  Estimate Std. Error     z value      Pr(>|z|)
(Intercept)                                     0.15929845 0.06960771   2.2885172  2.210742e-02
Raceasian                                      -0.28332904 0.20446278  -1.3857243  1.658311e-01
Raceblack                                      -2.48060238 0.09841886 -25.2045424 3.571728e-140
Racehispanic                                   -0.85340290 0.11786821  -7.2403147  4.476449e-13
Racenative american                            -0.39195681 0.14576603  -2.6889448  7.167827e-03
Raceother                                      -0.58168979 0.48829616  -1.1912643  2.335498e-01
RaceUnknown                                    -0.37908755 0.24233631  -1.5643036  1.177463e-01
Income17 to 33 percentile                      -0.06018027 0.06884955  -0.8740837  3.820727e-01
Income34 to 67 percentile                       0.02307343 0.06238488   0.3698561  7.114897e-01
Income68 to 95 percentile                       0.15514729 0.06327773   2.4518467  1.421252e-02
Income96 to 100 percentile                      0.73730304 0.09115114   8.0887969  6.025693e-16
IncomeUnknown                                   0.17534821 0.08894610   1.9713985  4.867831e-02
Gendermale                                      0.09542550 0.03472874   2.7477386  6.000783e-03
GenderUnknown                                   0.34912397 0.41827387   0.8346779  4.038991e-01
Educationgrade school of less (0-8 grades)     -0.32777601 0.06203652  -5.2835977  1.266713e-07
Educationhigh school (12 grades or fewer, incl -0.10646145 0.04702443  -2.2639605  2.357654e-02
Educationsome college(13 grades or more,but no  0.18023660 0.05430502   3.3189678  9.035085e-04
EducationUnknown                               -0.12275738 0.22466669  -0.5463978  5.847925e-01

係数データを表示すると、大量に数値が出てきて何が何やら分からなくなってしまいます。
そこで、coefplot を使って可視化してみます。

library(coefplot)
coefplot(model)

このグラフは、各変数に対する係数の推定値とその標準誤差を表しています。*2
青い太線が 1 標準誤差範囲、細線が 2 標準誤差範囲になります。
このグラフを見ると、

高収入(Income96 to 100 percentile)だと共和党に投票しやすい
黒人(Raceblack)は民主党に投票しやすい
Unknown や Others などの係数は誤差範囲が広い

といったことが見て取れるようになります。
便利ですね！

Gelman の「秘密兵器」

ここで趣向を変えて、各年ごとにモデルを作成し、係数がどう変化するかを見てみることにします。
まずは各年ごとのモデルを作成します。*3

library(pforeach)
theYears <- unique(ideo$Year)
npforeach(year=theYears, .c=list)({
  glm(Vote ~ Race + Income + Gender + Education,
      data=ideo[ideo$Year==year, ], family=binomial(link="logit"))
}) -> results
names(results) <- theYears

これらのモデルに対して、係数の時系列プロットを作成するには、coefplot パッケージの multiplot 関数を使えば簡単です。
例えば、高所得者についての係数の時系列プロットは次のようになります。

multiplot(results, coefficients="Income96 to 100 percentile", secret.weapon=TRUE) +
  coord_flip(xlim=c(-3, 3))

1972 年(ニクソン二期目)から、共和党に対する高所得者層の支持が上がっていることが分かります。
1971年のニクソン・ショックに関係がありそうですね。

このように、一連のモデルに対して係数の推定値を並べてプロットするというアイデアは、単純でありながら非常に便利です。
応用統計学の巨人 Andrew Gelman はこのようなプロットを、その有用さから「Secret Weapon(秘密兵器)」と呼びました。*4

データの不備を見抜く

同様に、黒人についての係数の時系列プロットを見てみましょう。

multiplot(results, coefficients="Raceblack", secret.weapon=TRUE) +
  coord_flip(xlim=c(-20, 10))

これはどうしたことでしょうか。
1964 年についてだけ、係数の誤差範囲が異常に広いです。
multiplot 関数の plot 引数を FALSE にして、データだけを表示させてみましょう。

multiplot(results, coefficients="Raceblack", secret.weapon=TRUE, plot = FALSE)

          Value Coefficient   HighInner     LowInner   HighOuter    LowOuter Model
1    0.07119541   Raceblack   0.6297813   -0.4873905   1.1883673   -1.045976  1948
2   -1.68490828   Raceblack  -1.3175506   -2.0522659  -0.9501930   -2.419624  1952
3   -0.89178359   Raceblack  -0.5857195   -1.1978476  -0.2796555   -1.503912  1956
4   -1.07674848   Raceblack  -0.7099648   -1.4435322  -0.3431811   -1.810316  1960
5  -16.85751152   Raceblack 382.1171424 -415.8321655 781.0917963 -814.806819  1964
6   -3.65505395   Raceblack  -3.0580572   -4.2520507  -2.4610605   -4.849047  1968
7   -2.68154861   Raceblack  -2.4175364   -2.9455608  -2.1535242   -3.209573  1972
8   -2.94158722   Raceblack  -2.4761518   -3.4070226  -2.0107164   -3.872458  1976
9   -3.03095296   Raceblack  -2.6276580   -3.4342479  -2.2243631   -3.837543  1980
10  -2.47703741   Raceblack  -2.1907106   -2.7633642  -1.9043839   -3.049691  1984
11  -2.79340230   Raceblack  -2.4509285   -3.1358761  -2.1084548   -3.478350  1988
12  -2.82977980   Raceblack  -2.4609290   -3.1986306  -2.0920782   -3.567481  1992
13  -4.45297040   Raceblack  -3.4433048   -5.4626360  -2.4336392   -6.472302  1996
14  -2.67827784   Raceblack  -2.2777557   -3.0788000  -1.8772336   -3.479322  2000

なんと、1964 年だけ、誤差範囲が 100 倍以上広くなっています。

これは一体なぜでしょうか？
各年のアンケートに答えた黒人の数を確認してみます。

library(dplyr)
ideo %>% 
  filter(Race == "black") %>%
  count(Year)

   Year   n
1  1948  15
2  1952  48
3  1956  49
4  1960  42
5  1964  92
6  1968  85
7  1972 134
8  1976 102
9  1980 104
10 1984 127
11 1988 119
12 1992 172
13 1996 109
14 2000 107

1964 年だけサンプルサイズが特に小さい、ということはないようです。

本書「みんなのR」では、この原因は "underrepresented"、すなわち、サンプルが母集団をうまく代表できていないためと結論づけています。

このように、coefplot を使って係数の推定値の比較を行うことで、データの不備を見つけることができました。
Gelman の「秘密兵器」は分析者にとって非常に強力な武器であることがわかります。

さらなる分析へ

本書「みんなのR」では、このデータに対して、ベイズ統計学の視点から立ち向かいます。
1964 年の係数が他の年の係数と比べて 100 倍になることは考えにくいと思います。
このことは、1964 年の係数は少なくともこの程度の範囲内にあるのではないか(例えば 9 割がた -15 〜 15 の範囲内に収まるだろう)ということを事前知識として持っていることになります。
この事前知識は、ベイズ統計モデルにおける、事前分布としてモデルに取り込むことができます。
このような分析も、R のパッケージを使えば簡単にできます。