『Stataによる計量政治学』（浅野正彦, 矢内勇生)：Rによる分析
Yuki Yanai
September 19, 2014 (rev. June 5, 2015)
第8章 平均値の比較
準備として、作業ディレクトリを指定する。 作業ディレクトリの指定にはsetwd()関数を利用する。 ディレクトリへのパス（引用符の中身）は、各自で設定する。
setwd("~/projects/Stata-R/")8.1 グループ平均の比較方法
8.1.1 Welch のt 検定（対応なし）
まず、表8.1 (p.140) を元に、データセット（Rのデータフレーム）を作る。 データフレームを作るには、data.frame() を使う。
## まず、必要な2つ変数score, starbucksを作る
score <- c(95, 95, 85, 90, 85, 75, 85, 85, 75, 65,
90, 85, 85, 80, 85, 70, 85, 75, 80, 60)
starbucks <- rep(c(0, 1), c(10, 10))
## 2つの変数（ベクトル）をもとに、データフレームを作る
coffee.unpaired <- data.frame(score = score,
starbucks = starbucks)次に、箱ひげ図を使ってデータを図示する (p.143, 図8.3)。
## ylim で縦軸の範囲を決める
boxplot(score ~ starbucks, data = coffee.unpaired, ylim = c(60, 100),
names = c("タリーズ", "スタバ"),
xlab = "コーヒー店", ylab = "味の評価（点）",
main = "コーヒーの味調査 箱ひげ図")
2つのグループの評価には統計的に意味のある違いがあるのだろうか。 Welchのt 検定で調べてみよう。 Rでは、t.test() でt 検定を行う。
t.test(score ~ starbucks, data = coffee.unpaired)##
## Welch Two Sample t-test
##
## data: score by starbucks
## t = 0.9717, df = 17.951, p-value = 0.3441
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -4.649864 12.649864
## sample estimates:
## mean in group 0 mean in group 1
## 83.5 79.5結果の読み方は、Stataの場合と同様である（pp.143-144, 表8.3の解説を参照）。 Stataの結果にある Pr(|T| > |t|) に対応するのは、p-value である。
8.1.2 対応のあるt 検定
ここでは、表8.2 の調査データを入力する。対応のあるt 検定を行う場合、データフレームを作る必要はない。
tul <- c(95, 95, 85, 90, 85, 75, 85, 85, 75, 65)
stb <- c(90, 85, 85, 80, 85, 70, 85, 75, 80, 60)対応のないときと同じように、t.test() を使って検定を行うが、使い方が異なるので、よく見比べてほしい。
t.test(tul, stb, paired = TRUE)##
## Paired t-test
##
## data: tul and stb
## t = 2.4495, df = 9, p-value = 0.03679
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.3059128 7.6940872
## sample estimates:
## mean of the differences
## 4このように、表8.4 (p.146) と同じ結果が得られる。 結果の読み方はStataと同様である（pp.145-146 を参照）。
8.1.3 サンプルサイズと標本平均
サンプルサイズの影響
表8.5のデータを入力し、検定を行う。
tul2 <- c(90, 95, 85, 80, 75)
stb2 <- c(95, 80, 80, 75, 75)
t.test(tul2, stb2, paired = TRUE)##
## Paired t-test
##
## data: tul2 and stb2
## t = 1.206, df = 4, p-value = 0.2943
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -5.208427 13.208427
## sample estimates:
## mean of the differences
## 4表8.6 と同じ結果が得られた（表の読み方はp.148 を参照）。
標本の平均値の影響
表8.7のデータを入力し、検定を行う。
tul3 <- tul2 ## タリーズは表8.5と同じ
stb3 <- c(80, 70, 75, 70, 70)
t.test(tul3, stb3, paired = TRUE)##
## Paired t-test
##
## data: tul3 and stb3
## t = 3.5386, df = 4, p-value = 0.02404
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 2.584617 21.415383
## sample estimates:
## mean of the differences
## 12表8.8 と同じ結果が得られた（表の読み方はpp.149-150 を参照）。
8.2 衆議院選挙データを使った平均値の比較
8.2.1 分析対象データの準備
分析用のデータ (hr96-09.dta) を読み込み、内容を確認する （第5章のRによる分析で作った hr96-09labeled.dat を使ってもよい。）。 RでStata形式のデータセット (拡張子が .dta のもの) を呼び出すには、foreign パッケージの read.dta() を使う。 データセットは、作業ディレクトリ内のdataフォルダ（ディレクトリ）にあるものとする。
## haven をインストール済みでないならインストールする
# install.packages("haven", dependencies = TRUE)
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) ## データの確認## year ku kun party name age status nocand wl rank
## 1 1996 aichi 1 NFP KAWAMURA, TAKASHI 47 incumbent 7 win 1
## 2 1996 aichi 1 LDP IMAEDA, NORIO 72 moto 7 lose 2
## 3 1996 aichi 1 DPJ SATO, TAISUKE 53 incumbent 7 lose 3
## 4 1996 aichi 1 JCP IWANAKA, MIHOKO 43 challenger 7 lose 4
## 5 1996 aichi 1 bunka ITO, MASAKO 51 challenger 7 lose 5
## 6 1996 aichi 1 NP YAMADA, HIROSHIB 51 challenger 7 lose 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 NAまず、2009年の衆議院議員総選挙小選挙区における獲得投票数 (vote) の記述統計を表示してみよう。 ここでは、表8.9 (p.151) に示されている、観測数、平均値（算術平均）、標準偏差、最小値、最大値をひとつずつ求めてみる。
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年のデータを取り出して保存
vote.obs <- length(HR09$vote) ## 観測値の数 (vote ベクトルの長さ)
vote.mean <- mean(HR09$vote) ## 算術平均
vote.sd <- sd(HR09$vote) ## 標準偏差
vote.min <- min(HR09$vote) ## 最小値
vote.max <- max(HR09$vote) ## 最大値
cbind(vote.obs, vote.mean, vote.sd, vote.min, vote.max) ## 並べて表示## vote.obs vote.mean vote.sd vote.min vote.max
## [1,] 1139 61939.59 53880.34 177 201461これで、表8.9と同じものが得られた。
次に、自民党所属を表すダミー変数ldp を作る。 Rでは、as.numeric() を使い、party変数からダミー変数を作ることができる。 より具体的には、as.numeric(条件) とすることで、条件に適合する観測値は1、適合しないものは0とするダミー変数ができる。 ここでは、party変数の値（ラベル）がLDPのものだけ1にする。
## LDP ダミーの作成
HR <- mutate(HR, ldp = as.numeric(party == "LDP"))これで、ldpダミーができた。 このldpダミーは候補者が自民党 (LDP) なら1を、それ以外の場合は0となっている。
この変数を元に、教科書と同じように候補者が自民党なら1、民主党 (DPJ) なら0、それ以外の場合は欠測値（Stata では「.」、Rでは NA で表す）になるように変更を加える。 また、これは実際にはLDPかDPJかを表す変数なので、factorに変えてラベルを付けよう。
## 欠測値 NA に変更すべき対象を指定する
## ここでは、party 変数が LDP でもDPJ でもないものを選ぶ
## %in% は左がわの変数の値が右側の値の集合（ここでは{LDP, DPJ})
## に含まれるかどうかを確かめる
## 私たちは「含まれない」ものを求めたいので、「== FALSE」とする
target <- (HR$party %in% c("LDP", "DPJ")) == FALSE
HR$ldp[target] <- NA ## 対象となる観測値をNAに置き換える
HR <- mutate(HR,
ldp = factor(ldp, levels = 0:1, labels = c("DPJ", "LDP")))
HR09 <- filter(HR, year == 2009) ## 元データに変数を加えたので、2009年分もアップデートこれで教科書と同じダミー変数ができる。
確認のため、表8.10 (p.152) を再現してみよう。
## 民主党候補の獲得票数
tbl.dpj <- HR09 %>%
filter(ldp == "DPJ") %>%
with(c(mean(vote), sd(vote), length(vote)))
## 自民党候補の獲得票数
tbl.ldp <- HR09 %>%
filter(ldp == "LDP") %>%
with(c(mean(vote), sd(vote), length(vote)))
## 民主党と自民党の合計：「!is.na(ldp)」はldpが欠測値NAではないということ
tbl.ttl <- HR09 %>%
filter(!is.na(ldp)) %>%
with(c(mean(vote), sd(vote), length(vote)))
rbind(c("mean", "sd", "obs."), tbl.dpj, tbl.ldp, tbl.ttl)## [,1] [,2] [,3]
## "mean" "sd" "obs."
## tbl.dpj "123558.409594096" "30792.5791130804" "271"
## tbl.ldp "94490.0515463918" "23604.136364441" "291"
## tbl.ttl "108507" "30914.5296941477" "562"これで表8.10 と同じものが得られた。
8.2.2 集計データを使った平均値比較
2009年衆院選における自民党候補者と民主党候補者の獲得票数のデータを図示してみよう。
boxplot(vote ~ ldp, data = HR09, names = c("民主党", "自民党"),
xlab = "政党", ylab = "獲得票数",
main = "民主党・自民党候補者の票 箱ひげ図")
これで図8.5 (p.152) が得られる。
Welchのt 検定
上の図に見られる民主党と自民党の得票の差が統計的に有意なものか検定してみよう。
t.test(vote ~ ldp, data = HR09)##
## Welch Two Sample t-test
##
## data: vote by ldp
## t = 12.4935, df = 505.435, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 24497.20 33639.52
## sample estimates:
## mean in group DPJ mean in group LDP
## 123558.41 94490.05これで、表8.11 (p.153) と同じ結果が得られた。