『Stataによる計量政治学』（浅野正彦, 矢内勇生)：Rによる分析
Yuki Yanai
September 18, 2014 （rev. June 4, 2015)
第5章 記述統計
5.0 準備
R では、様々な機能が R パッケージとして提供される。 R と共に最初からインストールされるパッケージもあるが、大部分のパッケージは 自分でインストールする必要がある。
今回は、データセットの扱いを容易にするする dplyr パッケージを使うので、インストールして それをR で呼び出す。
まず、インストールには、install.packages() 関数を使う。この関数の引数は、インストールするパッケージの名前である。また、このパッケージが依存するパッケージが他にある場合、一緒にインストールするように、dependencies = TRUE を指定する。
install.packages("dplyr", dependencies = TRUE)インストールが完了したら、このRセッションで使えるように、library() で呼び出す。
library("dplyr")##
## Attaching package: 'dplyr'
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union5.1 変数の種類と記述統計
準備として、作業ディレクトリを指定する。 作業ディレクトリの指定にはsetwd()関数を利用する。 ディレクトリへのパス（引用符の中身）は、各自で設定する。
setwd("~/projects/Stata-R/")続いて、衆議院議員総選挙のデータ（hr96-09.csv）を読み込む。 読み込むデータセットは、作業ディレクトリ内のdata というフォルダ（ディレクトリ）にあるものとする。 CSV形式のデータの読み込みには、read.csv()を使う。データの名前はHR（House of Representatives: 衆議院）にする。
HR <- read.csv("data/hr96-09.csv", na.st = ".")上のコードにある na.st = "." は、元のデータセット（CSV形式）内で欠測値が「.」（ピリオド）で表されていることをRに伝えている。 Rにおける欠測値の既定値は NA なので、それ以外の方法で欠測値を記録したデータを読むときは、na.st で欠測値の形式を指定する。
次に、データがきちんと読み込めたかどうか確認する。
names(HR) ## データセットに含まれる変数名を表示する## [1] "year" "ku" "kun" "party" "name"
## [6] "age" "status" "nocand" "wl" "rank"
## [11] "previous" "vote" "voteshare" "eligible" "turnout"
## [16] "exp"head(HR) ## データセットの一部（最初の6行）を表示する## year ku kun party name age status nocand wl rank
## 1 1996 aichi 1 1000 KAWAMURA, TAKASHI 47 2 7 1 1
## 2 1996 aichi 1 800 IMAEDA, NORIO 72 3 7 0 2
## 3 1996 aichi 1 1001 SATO, TAISUKE 53 2 7 0 3
## 4 1996 aichi 1 305 IWANAKA, MIHOKO 43 1 7 0 4
## 5 1996 aichi 1 1014 ITO, MASAKO 51 1 7 0 5
## 6 1996 aichi 1 1038 YAMADA, HIROSHIB 51 1 7 0 6
## previous vote voteshare eligible turnout exp
## 1 2 66876 40.0 346774 49.22 9828097
## 2 3 42969 25.7 346774 49.22 9311555
## 3 2 33503 20.1 346774 49.22 9231284
## 4 0 22209 13.3 346774 49.22 2177203
## 5 0 616 0.4 346774 49.22 NA
## 6 0 566 0.3 346774 49.22 NAdim(HR) ## データセットに含まれる観測数と変数の数を確認する (dimはdimension)## [1] 5614 16などで確認する。
データセットに含まれるすべての変数の基本的な統計量を確認するには、summary() 関数を使う。
summary(HR)## year ku kun party
## Min. :1996 tokyo : 529 Min. : 1.000 Min. : 1.0
## 1st Qu.:2000 kanagawa: 373 1st Qu.: 2.000 1st Qu.: 305.0
## Median :2003 osaka : 360 Median : 4.000 Median : 800.0
## Mean :2002 aichi : 314 Mean : 5.695 Mean : 727.4
## 3rd Qu.:2005 saitama : 280 3rd Qu.: 8.000 3rd Qu.:1001.0
## Max. :2009 chiba : 247 Max. :25.000 Max. :9998.0
## (Other) :3511
## name age status nocand
## ABE, SHINZO : 5 Min. :25.00 Min. :1.000 Min. :2.000
## AISAWA, ICHIRO : 5 1st Qu.:42.00 1st Qu.:1.000 1st Qu.:3.000
## AKABA, KAZUYOSHI : 5 Median :51.00 Median :1.000 Median :4.000
## AKAGI, NORIHIKO : 5 Mean :50.52 Mean :1.445 Mean :3.987
## AKAMATSU, HIROTAKA: 5 3rd Qu.:58.00 3rd Qu.:2.000 3rd Qu.:5.000
## AKUTSU, YUKIHIKO : 5 Max. :89.00 Max. :3.000 Max. :9.000
## (Other) :5584
## wl rank previous vote
## Min. :0.0000 Min. :1.000 Min. : 0.000 Min. : 177
## 1st Qu.:0.0000 1st Qu.:1.000 1st Qu.: 0.000 1st Qu.: 16235
## Median :0.0000 Median :2.000 Median : 0.000 Median : 50757
## Mean :0.4442 Mean :2.492 Mean : 1.702 Mean : 56204
## 3rd Qu.:1.0000 3rd Qu.:3.000 3rd Qu.: 3.000 3rd Qu.: 89807
## Max. :2.0000 Max. :9.000 Max. :20.000 Max. :201461
##
## voteshare eligible turnout exp
## Min. : 0.10 Min. :191241 Min. :48.89 Min. : 5540
## 1st Qu.: 7.90 1st Qu.:294013 1st Qu.:58.23 1st Qu.: 2913701
## Median :26.40 Median :343654 Median :64.31 Median : 7032390
## Mean :26.72 Mean :341226 Mean :63.78 Mean : 7958713
## 3rd Qu.:43.00 3rd Qu.:393677 3rd Qu.:68.54 3rd Qu.:11754054
## Max. :95.30 Max. :487837 Max. :83.75 Max. :27179308
## NA's :134このように、数字で表されるデータについては最小値 (Min.)、第1四分位数 (1st Qu.)、中央値 (Median)、算術平均 (Mean)、第3四分位数 (3rd Qu.)、最大値 (Max.) が表示され、カテゴリー変数についてはそれぞれのカテゴリーにいくつの観測値があるか（の一部）を表示する。NA’s はその変数の欠測値の数を示す。
変数 wlの度数分布表（表5.6）を作るには、table() 関数を使う。データセット（Rでは「データフレーム (data frame)」と呼ぶ）に含まれる変数を使うときは データフレーム名$変数名 として変数にアクセスする。
wl.freq <- table(HR$wl) ## wl の度数分布表をwl.freqとして保存
wl.freq ## wl.freq を画面に表示##
## 0 1 2
## 3617 1500 497この方法だと度数しか表示されないので、相対度数（%）も表示されるように改良する。
cbind(wl.freq, 100 * wl.freq / sum(wl.freq)) ## 度数と相対度数 (%) を並べて表示## wl.freq
## 0 3617 64.428215
## 1 1500 26.718917
## 2 497 8.852868ここで、cbind()は複数のベクトルを縦に（各ベクトルを列として）並べる関数である。 横に（各ベクトルを行として）並べるときはrbind()を使う。 また、相対度数の計算に登場するsum()はベクトルの要素を合計する関数である。 ここでは、総度数を計算している。
以上の方法で変数wlの度数分布表を作ることができたが、wl の値0, 1, 2がそれぞれ何を表しているか明らかではないので、少々不便である。教科書5.1.5 (pp.68-) と同じように、値の意味を定義しよう。 Rでは変数の種類（クラス[class] と呼ばれる）をfactorに変更した上で、Stataと同じようにラベルを添付することができる。
まず、変数wl の現在のクラスを確認してみよう。
class(HR$wl)## [1] "integer"変数wl のクラスは integer すなわち整数値であることがわかる。 これをfactor()によってfactor に変換し、値にラベルを付ける。
HR$wl <- factor(HR$wl, levels = 0:2, labels = c("lose", "win", "zombie"))
class(HR$wl)## [1] "factor"変数wl のクラスがfactorに変換されていることがわかる。 これを利用して、もう一度度数分布表を作ってみよう。 今回は、addmargins()を使って総度数 (Sum) も求める。
wl.freq <- table(HR$wl)
## margin=1 で縦方向の合計を出す; margin=2だと横方向
addmargins(cbind(wl.freq, 100 * wl.freq / sum(wl.freq)), margin = 1)## wl.freq
## lose 3617 64.428215
## win 1500 26.718917
## zombie 497 8.852868
## Sum 5614 100.000000これで、表5.6 (p.69) に似た度数分布表ができる。
同様の方法で、変数status の度数分布表も作ってみよう。 まず、status のクラスを確認する。
class(HR$status)## [1] "integer"これもinteger なので、factorに変換して値にラベルを付ける。
HR$status <- factor(HR$status, levels = 1:3,
labels = c("challenger", "incumbent", "moto"))
status.freq <- table(HR$status)
addmargins(cbind(status.freq, 100 * status.freq / sum(status.freq)), margin = 1)## status.freq
## challenger 3411 60.75882
## incumbent 1908 33.98646
## moto 295 5.25472
## Sum 5614 100.00000これで、表5.8 (p.70) と同様の度数分布表ができた。
さらに、party（政党名）にもラベルを付ける。
party.names <- c("msz", "JCP", "LDP", "CGP", "oki",
"tai", "saki", "NFP", "DPJ", "SDP",
"LF", "NJSP", "DRF", "kobe", "nii",
"sei", "JNFP", "bunka", "green", "LP",
"RC", "muk", "CP", "NCP", "ND",
"son", "sek", "NP", "NNP", "NPJ",
"NPD", "minna", "R", "H", "sho")
HR$party <- factor(HR$party, labels = party.names)
table(HR$party)##
## msz JCP LDP CGP oki tai saki NFP DPJ SDP LF NJSP
## 396 1326 1417 43 1 1 13 235 1212 247 212 38
## DRF kobe nii sei JNFP bunka green LP RC muk CP NCP
## 2 2 1 1 1 10 1 61 4 9 16 11
## ND son sek NP NNP NPJ NPD minna R H sho
## 1 1 2 11 19 8 1 14 1 292 4これを利用して、表5.11 (p.72) と同様に、2009年衆院選での所属政党の度数分布表を作ってみよう。 ここでは、with() 関数を利用する。この関数はデータフレームの中に含まれる変数に対して関数を当てはめるときに便利なもので、with(データフレーム名, 関数名(変数名))として利用する。 また、2009年のデータだけを取り出すために dplyr パッケージの filter() を使う。 使い方は、filter(データフレーム名, 条件)である。 ここではyear が2009のものだけを使いたいので、条件はyear==2009とする（year=2009 [等号がひとつ]ではダメ）。
HR09 <- filter(HR, year == 2009) ## 2009年のデータをHR09として保存
party.freq <- with(HR09, table(party))
addmargins(cbind(party.freq, 100 * party.freq / sum(party.freq)), margin = 1)## party.freq
## msz 70 6.14574188
## JCP 152 13.34503951
## LDP 291 25.54872695
## CGP 6 0.52677788
## oki 0 0.00000000
## tai 0 0.00000000
## saki 0 0.00000000
## NFP 0 0.00000000
## DPJ 271 23.79280070
## SDP 31 2.72168569
## LF 0 0.00000000
## NJSP 0 0.00000000
## DRF 0 0.00000000
## kobe 0 0.00000000
## nii 0 0.00000000
## sei 0 0.00000000
## JNFP 0 0.00000000
## bunka 0 0.00000000
## green 0 0.00000000
## LP 0 0.00000000
## RC 0 0.00000000
## muk 0 0.00000000
## CP 0 0.00000000
## NCP 0 0.00000000
## ND 0 0.00000000
## son 0 0.00000000
## sek 0 0.00000000
## NP 0 0.00000000
## NNP 9 0.79016681
## NPJ 2 0.17559263
## NPD 0 0.00000000
## minna 14 1.22914838
## R 1 0.08779631
## H 292 25.63652327
## sho 0 0.00000000
## Sum 1139 100.00000000この表と教科書の表5.11 との大きな違いは、表5.11 では度数が0の政党（つまり2009年に候補者を立てなかった政党）が除外されているのに対し、ここで作った表にはそれらの政党が含まれている点にある。 度数が0でも表にカウントするというのがRのfactorクラスの特徴である（factorではない場合、度数0は表示されない）。度数0を表示させないようにするには、次のようにする。
party.freq2 <- party.freq[party.freq > 0] ## party.freq の度数が0より大きいものだけ残す
addmargins(cbind(party.freq2, 100 * party.freq2 / sum(party.freq2)), margin = 1)## party.freq2
## msz 70 6.14574188
## JCP 152 13.34503951
## LDP 291 25.54872695
## CGP 6 0.52677788
## DPJ 271 23.79280070
## SDP 31 2.72168569
## NNP 9 0.79016681
## NPJ 2 0.17559263
## minna 14 1.22914838
## R 1 0.08779631
## H 292 25.63652327
## Sum 1139 100.00000000このように、Rではベクトル名[条件]とすることで、ベクトルの一部を取り出すことができる。 （as.character(party) として、一時的にcharacter [文字列] クラスとして扱う方法もある。）
次に、2009年衆院選における政党別の当選者数を調べよう（教科書 p.73, 表5.12）。 ここでは、%>% （パイプ）を使う。これは、左側で処理した内容を、右側の関数の第1引数として引き渡すための表現である。
## filter() でHR09のうちwlが"win"のものだけ取り出し、
## %>% で次の関数 with() の第1引数 として渡す
pty.win.freq <- HR09 %>%
filter(wl == "win") %>%
with(table(party))
## この代わりに次の行でも同じ
# pty.win.freq <- with(filter(HR09, wl == "win"), table(party))
pty.win.freq <- pty.win.freq[pty.win.freq > 0] ## 度数0を除外
tbl.pty.win.freq <- cbind(pty.win.freq, 100 * pty.win.freq / sum(pty.win.freq))
addmargins(tbl.pty.win.freq, margin = 1)## pty.win.freq
## msz 6 2.0000000
## LDP 64 21.3333333
## DPJ 221 73.6666667
## SDP 3 1.0000000
## NNP 3 1.0000000
## NPJ 1 0.3333333
## minna 2 0.6666667
## Sum 300 100.0000000これで、表5.12と同じような度数分布表ができた。
この表を見ると、みんなの党 (minna) は2人の当選者を出していることがわかる。 この2人についての詳細（各変数の値）を知りたい場合は、次のようにする。
filter(HR09, wl == "win", party == "minna")## year ku kun party name age status nocand wl rank
## 1 2009 kanagawa 8 minna EDA, KENJI 53 incumbent 4 win 1
## 2 2009 tochigi 3 minna WATANABE, YOSHIMI 57 incumbent 2 win 1
## previous vote voteshare eligible turnout exp
## 1 3 128753 49.1 376393 70.82 6851426
## 2 5 142482 95.3 385445 68.45 10621917特定の変数についてのみ表示したいときは、select() を使う。
HR09 %>%
filter(wl == "win", party == "minna") %>%
select(ku, kun, name, age, nocand, rank, vote)## ku kun name age nocand rank vote
## 1 kanagawa 8 EDA, KENJI 53 4 1 128753
## 2 tochigi 3 WATANABE, YOSHIMI 57 2 1 142482これで、表5.13 (p.74) ができる。 データセットには行（各観測値）と列（各変数）があるが、特定の観測値にアクセスするときには データフレーム名[観測値の指定, ], 特定の変数にアクセスしたい場合には データフレーム名[, 変数の指定] とする。両者を同時に指定することもできる。 指定の方法は複数ある。ここでは、filter()でデータフレーム自体を縮小し、残ったすべての観測対象について、select() で複数の変数を指定した。
表5.14 (p.74) と同じように江田けんじ氏の選挙履歴を知りたければ、
HR %>%
filter(name == "EDA, KENJI") %>%
select(year, party, age, nocand, rank, status, wl, previous, vote)## year party age nocand rank status wl previous vote
## 1 2000 LDP 44 6 2 challenger lose 0 58787
## 2 2003 msz 47 4 2 incumbent lose 1 78782
## 3 2005 msz 49 4 1 moto win 2 88098
## 4 2009 minna 53 4 1 incumbent win 3 128753とすればよい。
ラベルをつけ終わったら、データを保存しよう。 Rには様々な保存方法があるが、Rのデータ構造をそのまま残すには save()を使う。
## data フォルダの中にhr96-09labeled.dat というファイル名で保存する
save(HR, file = "data/hr96-09labeled.dat")この形式で保存したデータを呼び出すときは、load()を使う。
load("data/hr96-09labeled.dat")
## HR <- load("data/hr96-09labeled.dat") ではない！ 5.2 量的変数の視覚化
5.2.1 幹葉図 (stem-and-leaf plot)
Rでは、stem() で幹葉図を作ることができる。 教科書と同じように、2009年衆院選で小選挙区から立候補した社民党候補者の年齢を幹葉図にしてみよう。 幹葉図の階級の数は scale で調整する（既定値は1。好みの図にするには試行錯誤が必要）。
SDP09 <- filter(HR09, party == "SDP") ## 2009年社民党候補を取り出して保存
with(SDP09, stem(age, scale = 0.5))
with(SDP09, stem(age))##
## The decimal point is 1 digit(s) to the right of the |
##
## 3 | 1237
## 4 | 299
## 5 | 011233577
## 6 | 0001111245679
## 7 | 22
##
##
## The decimal point is 1 digit(s) to the right of the |
##
## 3 | 123
## 3 | 7
## 4 | 2
## 4 | 99
## 5 | 011233
## 5 | 577
## 6 | 000111124
## 6 | 5679
## 7 | 22このように図5.2, 図5.3 (p.77) と同じ図ができる。
5.2.3 ヒストグラム (histogram)
Rでヒストグラムを作るには、hist() を利用する（ggplot2 のgeom_histogram() を使うとより見栄えの良いものができるが、ここでは解説しない。ggplot2のページ などを参照されたい）。 教科書の図5.4, 図5.5 (p.79) と同様に、2009年衆院選における社民党候補の年齢のヒストグラムを作ってみよう。
## breaks で階級の境界値を設定
## right = FALSE で階級の下限値は含み、上限値は含まないようにする（Stataに合わせるため）
## right の既定値はTRUE （下限は含まず、上限は含む）
## xlab, ylab で横軸、縦軸のラベルをつける
## main で図のキャプションをつける
hist(SDP09$age, breaks = seq(30, 80, by = 10), right = FALSE,
xlab = "年齢", ylab = "度数",
main = "2009年社民党候補の年齢のヒストグラム")
## probability = TRUE で縦軸を確率密度にする
## freq = FALSE としてもまったく同じ
hist(SDP09$age, breaks = seq(30, 80, by = 10), right = FALSE,
probability = TRUE,
xlab = "年齢", ylab = "確率密度",
main = "2009年社民党候補の年齢のヒストグラム")
P.80で説明している通り、ヒストグラムは作り方によって見た目が変化するので、設定を色々変えて試してみるとよい。
試しに、図5.7 (p.80) を再現してみよう。
hist(SDP09$age, right = FALSE,
breaks = seq(min(SDP09$age), max(SDP09$age), length = 8),
xlab = "年齢", ylab = "度数", main = "2009年社民党候補の年齢のヒストグラム")
5.2.3 箱ひげ図 (box-and-whisker plot)
箱ひげ図は、boxplot() で作れる（ggplot2 でも作れる）。 2009年衆院選における社民党候補の年齢を箱ひげ図にしてみよう。
boxplot(SDP09$age, ylab = "年齢", main = "2009年社民党候補の年齢の箱ひげ図")
これで図5.9 (p.82) と同じ図ができる。
次に、2009年衆院選の政党別小選挙区当選者に関する年齢の箱ひげ図（図5.10, p.84）を再現してみよう。
## party を使うと2009年に参加していない政党も表示してしまう (factorなので)
## --> 一時的にfactorクラスをやめるため、as.character() で文字列として扱う
boxplot(age ~ as.character(party), data = subset(HR09, wl == "win"),
ylab = "年齢",
main = "2009年衆院選当選政党別候補者年齢の箱ひげ図")
政党の整列順序を除けば、図5.10 (p.84) と同じものができた。
箱ひげ図だけでは具体的な統計量の値がわからないので、自民党当選者の年齢の最小値、第1四分位数、中央値、平均値、第3四分位数、最大値と標準偏差を求めてみよう。 Rではsd() で標準偏差を、summary(変数) で変数の五数と平均値を求めることができる。
## c() を使ってwith()に複数の関数を評価させる: list() も可
HR09 %>%
filter(party == "LDP", wl == "win") %>%
with(c(summary(age), sd(age)))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 28.00000 49.75000 56.00000 56.03000 64.00000 74.00000 10.07585これで表5.16 (p.85) と同等の結果が得られる。 多少のズレがあるのは、StataとRで四分位数の計算方法が異なるためである。 今回の分析では、結果に影響が出るほど大きな差はない。
同様に、民主党、社民党の場合は、
HR09 %>%
filter(party == "DPJ", wl == "win") %>%
with(c(summary(age), sd(age)))
HR09 %>%
filter(party == "SDP", wl == "win") %>%
with(c(summary(age), sd(age)))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 28.00000 40.00000 47.00000 48.96000 59.00000 77.00000 10.46523
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 49.000000 56.500000 64.000000 60.000000 65.500000 67.000000 9.643651で表5.17, 表5.18と同じ結果が得られる。
最後に、2009年衆院選小選挙区における社民党当選者のプロフィールを確認してみよう。
SDP09 %>%
filter(wl == "win") %>%
select(name, ku, kun, party, age)## name ku kun party age
## 1 SHIGENO, YASUMASA oita 2 SDP 67
## 2 TERUYA, KANTOKU okinawa 2 SDP 64
## 3 TSUJIMOTO, KIYOMI osaka 10 SDP 49これで表5.19 (p.86) と同じ結果が得られる。