『Stataによる計量政治学』（浅野正彦, 矢内勇生)：Rによる分析
Yuki Yanai
September 19, 2014 (rev. June 5, 2015)
第9章 変数間の関連性
準備として、作業ディレクトリを指定する。 作業ディレクトリの指定にはsetwd()関数を利用する。 ディレクトリへのパス（引用符の中身）は、各自で設定する。
setwd("~/projects/Stata-R/")続いて、分析用のデータ (hr96-09.dta) を読み込み、内容を確認する （第5章のRによる分析で作った hr96-09labeled.dat を使ってもよい。）。 データセットは、作業ディレクトリ内のdataフォルダ（ディレクトリ）にあるものとする。
library("foreign")
HR <- read.dta("data/hr96-09.dta")
## hr96-09labeled.dat を使う場合は、
## load("data/hr96-09.dat")
## HR <- load("data/hr96-09.dat") ではない！！！
# head(HR) ## データの確認：今回は実行しない（各自で確認してください）9.1 カテゴリー変数間の関連
9.1.1 クロス集計表
表9.2 (p.158) を作る。
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, unionHR09 <- filter(HR, year == 2009) ## 2009年のデータを取り出して保存
tbl.st.wl <- with(HR09, table(status, wl)) ## status とwl のクロス表を作る
tbl.st.wl <- tbl.st.wl[1:2, 1:2] ## status が"moto"の行と、wlが"zombie"の列を除外
addmargins(tbl.st.wl) ## 周辺度数を加えて表示する
## marginを指定しないと、行・列両方の周辺度数が加えられる## wl
## status lose win Sum
## challenger 559 81 640
## incumbent 168 174 342
## Sum 727 255 982これで同じ表ができた。
行パーセント
比率の表を作るには、まず上の方法で度数の表を作った後、prop.table()で比率に変換するのが簡単である。まず、表9.3 (p.158) を再現してみよう。
## まず列周辺度数を加える; margin=1 で列周辺度数のみ加える
row.pr.st.wl <- addmargins(tbl.st.wl, margin = 1)
## 度数を比率に変換する; margin=1 で行パーセント
row.pr.st.wl <- prop.table(row.pr.st.wl, margin = 1)
## パーセント表示に変換する; 小数第2位までの表示にする
row.pr.st.wl <- round(100 * row.pr.st.wl, 2)
## 行周辺度数を加えて表示する; margin=2 で行周辺度数のみ加える
addmargins(row.pr.st.wl, margin = 2)## wl
## status lose win Sum
## challenger 87.34 12.66 100
## incumbent 49.12 50.88 100
## Sum 74.03 25.97 100これで表9.3 と同じ比率が得られた。
列パーセント
同様に、列パーセントを計算する。 行パーセントで使ったコードの margin の値を帰ればよい。
## まず行周辺度数を加える; margin=2 で行周辺度数のみ加える
col.pr.st.wl <- addmargins(tbl.st.wl, margin = 2)
## 度数を比率に変換する; margin=2 で行パーセント
col.pr.st.wl <- prop.table(col.pr.st.wl, margin = 2)
## パーセント表示に変換する; 小数第2位までの表示にする
col.pr.st.wl <- round(100 * col.pr.st.wl, 2)
## 列周辺度数を加えて表示する; margin=1 で列周辺度数のみ加える
addmargins(col.pr.st.wl, margin = 1)## wl
## status lose win Sum
## challenger 76.89 31.76 65.17
## incumbent 23.11 68.24 34.83
## Sum 100.00 100.00 100.00これで表9.4 (p.159) と同じ比率が得られた。
全体パーセント
全体パーセントを求めるときは、margin を指定しなければよい。
## 度数を比率に変換する; margin=NULL（既定値）で全体パーセント
ttl.pr.st.wl <- prop.table(tbl.st.wl)
## パーセント表示に変換する; 小数第2位までの表示にする
ttl.pr.st.wl <- round(100 * ttl.pr.st.wl, 2)
## 周辺度数を加えて表示する
addmargins(ttl.pr.st.wl)## wl
## status lose win Sum
## challenger 56.92 8.25 65.17
## incumbent 17.11 17.72 34.83
## Sum 74.03 25.97 100.00これで、表9.5 (p.160) と同じ比率が得られた。
9.1.2 \(\chi^2\) 検定
StataR によるカイ二乗検定
まず、表9.8 (p.166) に示されている架空のデータセットを作ろう。
## male とsupport という2つの変数をもつデータフレームを作り、Fakeとして保存
## 最初はすべての値が欠測 NA
Fake <- data.frame(male = rep(NA, 100), support = rep(NA, 100))
Fake$male <- rep(c(0, 1), c(50, 50)) ## male に値を入れる
Fake$support <- rep(c(1, 0, 1), c(20, 50, 30)) ## support に値を入れる
## データを確認したければ
## Fake
## または
## View(Fake)このデータセットを使い、\(\chi^2\) 検定を行う。 Rでは、chisq.test() を利用する。
## イェーツの連続性補正は行わないので correct = FALSE
chisq.test(Fake$male, Fake$support, correct = FALSE)##
## Pearson's Chi-squared test
##
## data: Fake$male and Fake$support
## X-squared = 4, df = 1, p-value = 0.0455これで表9.9 (p.167) と同じ結果が得られた。
フィッシャーの直接確率計算法
まず、表9.10 (p.167) を元に架空のデータを入力する。 今回は、表（行列）として保存しよう。
## 行列を作る
fake.table <- matrix(c(3, 2, 2, 1), nrow = 2, byrow = TRUE)
## 行に名前を付ける
row.names(fake.table) <- c("female", "male")
## 列に名前を付ける
colnames(fake.table) <- c("not support", "support")
## 周辺度数を加えて表示。
addmargins(fake.table)## not support support Sum
## female 3 2 5
## male 2 1 3
## Sum 5 3 8このデータにFisherの直接確率計算法を適用する。 Rでは、fisher.test()を使う。
## alternative で片側検定（女性の支持 > 男性の支持）を指定;
## 「女性 < 男性」 を検定するときは alternative="greater" とする
## 既定値は両側検定
fisher.test(fake.table, alternative = "less")##
## Fisher's Exact Test for Count Data
##
## data: fake.table
## p-value = 0.7143
## alternative hypothesis: true odds ratio is less than 1
## 95 percent confidence interval:
## 0.00000 17.73797
## sample estimates:
## odds ratio
## 0.7772203これで表9.11 (p.168) と同じ結果が得られた。
9.2 数量変数間の関連
9.2.2 相関係数を使った統計的仮説検定
ここからは、上で既に作った HR09 というデータフレームを使って分析を進める。
選挙費用 (exp) と得票数 (vote) の関係をを調べるため、まず散布図を作る。 Rではplot() 関数で散布図を描くことができる。
plot(vote ~ exp, data = HR09, xlab = "選挙費用", ylab = "得票数",
main = "選挙費用と獲得票数：2009年（小選挙区）")
このように、図9.6 (p.171) によく似た図ができる。 もう少し奇麗な図を描きたいときは、ggplot2 を使うとよい。
## ggplot2 のインストールが必要なら
## install.packages('ggplot2')
library("ggplot2")
## Windows ユーザは次の行をコメントアウトせよ
theme_set(theme_gray(base_size=12, base_family="HiraKakuProN-W3"))
scat <- ggplot(HR09, aes(x = exp, y = vote)) +
geom_point() +
labs(x = "選挙費用", y = "得票数",
title = "選挙費用と獲得票数：2009年（小選挙区）")
print(scat)
続いて、選挙費用と得票数の相関係数をcor() で求め、cor.test() で相関係数の\(p\)値を求める。
with(na.omit(HR09), cor(exp, vote)) ## na.omit() で欠測値を除く## [1] 0.6616138with(na.omit(HR09), cor.test(exp, vote))##
## Pearson's product-moment correlation
##
## data: exp and vote
## t = 29.5549, df = 1122, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6274137 0.6932663
## sample estimates:
## cor
## 0.6616138これで、表9.12 (p.171) の結果が得られた。
次に、分析対象をみんなの党 (minna) に絞って同様の分析を行う。 まず、散布図を描く。
scat.minna <- ggplot(filter(HR09, party == "minna"), aes(x = exp, y = vote)) +
geom_point(size = 5) +
labs(x = "選挙費用", y = "得票数",
title = "みんなの党（選挙費用と獲得票数（2009、小選挙区）")
print(scat.minna)
これで図9.7 (p.172) と同様の図が得られた。
次に、選挙費用と得票数の相関係数をcor() で求め、cor.test() で相関係数の\(p\)値を求める。
HR09 %>%
filter(party == "minna", !is.na(vote), !is.na(exp)) %>%
with(cor(exp, vote))## [1] 0.5317654HR09 %>%
filter(party == "minna", !is.na(vote), !is.na(exp)) %>%
with(cor.test(exp, vote))##
## Pearson's product-moment correlation
##
## data: exp and vote
## t = 2.1751, df = 12, p-value = 0.05033
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.001651972 0.828569209
## sample estimates:
## cor
## 0.5317654これで表9.13 (p.173) と同じ結果が得られた。