政治学方法論 I：Rによる統計分析の基礎
矢内 勇生
2016年4月13日
このページでは、Rを使いながら統計学の基礎を復習する。 1回の授業で復習する範囲には限りがあるので、だいぶ大雑把な説明になっている。 基礎に不安がある場合は、シラバスに記載されている統計学の参考書で復習すること。
データの要約
データを読み込む前に、setwd()で作業ディレクトリ（フォルダ）を指定する。 たとえば、ホームフォルダの中に“methods”というフォルダを作り、それを作業ディレクトリに指定するなら、
setwd('~/methods/')とする。RStudio を使うときは、Project を指定すれば自動的に作業ディレクトリが指定される。 現在の作業ディレクトリを確認したいときはConsole（プロンプト）でgetwd() と入力する。
次にデータセット fake-data-lec02.csv をダウンロードして指定した作業ディレクトリに保存する。
RのConsole（プロンプト）にdir() と入力し、使いたいデータが現在の作業ディレクトリにあるかどうか確認する。 確認できたら、read.csv()を使い、データセットを読み込んで名前（ここでは myd）を付ける。
myd <- read.csv("fake-data-lec02.csv")読み込んだデータを確認してみよう。 まず、データセットの最上部または最下部にある6行分だけを表示する。
head(myd) ## 最初の6行分## id sex age height weight income
## 1 1 male 37 168.8 59.4 5311751
## 2 2 male 51 167.5 64.6 3819610
## 3 3 male 44 178.2 70.9 1218323
## 4 4 male 28 167.2 61.6 10984692
## 5 5 male 66 165.5 64.7 2318699
## 6 6 male 22 160.7 65.4 6581828tail(myd) ## 最後の6行分## id sex age height weight income
## 95 95 female 57 165.1 57.3 1686238
## 96 96 female 25 165.1 59.6 2629502
## 97 97 female 46 153.7 46.8 3338817
## 98 98 female 67 154.8 48.8 5997860
## 99 99 female 58 151.7 34.4 2324978
## 100 100 female 28 158.3 53.2 1607819データセットに含まれる変数名を確認したいときは、
names(myd)## [1] "id" "sex" "age" "height" "weight" "income"とする。
データセットに含まれる観測数 (\(n\)) と変数の数を知りたいときは、
dim(myd)## [1] 100 6とする。最初の数字が\(n\)の数（データセットの行数）、2番目の数字が変数の数（列数）である（RStudio を使うと、右上のウィンドウの“Environment” というタブに、この情報が表示される）。
データセットの構造を確認するには、
str(myd)## 'data.frame': 100 obs. of 6 variables:
## $ id : int 1 2 3 4 5 6 7 8 9 10 ...
## $ sex : Factor w/ 2 levels "female","male": 2 2 2 2 2 2 2 2 2 2 ...
## $ age : int 37 51 44 28 66 22 68 54 59 65 ...
## $ height: num 169 168 178 167 166 ...
## $ weight: num 59.4 64.6 70.9 61.6 64.7 65.4 63.9 42.1 58.9 56.7 ...
## $ income: int 5311751 3819610 1218323 10984692 2318699 6581828 3819509 4674633 4141410 754043 ...とする（詳しい読み方は次回以降の授業で）。
データセットに含まれるすべての変数の基本的な統計量を確認したいときは、
summary(myd)## id sex age height
## Min. : 1.00 female:50 Min. :20.00 Min. :145.9
## 1st Qu.: 25.75 male :50 1st Qu.:32.00 1st Qu.:158.1
## Median : 50.50 Median :43.00 Median :161.6
## Mean : 50.50 Mean :44.04 Mean :163.1
## 3rd Qu.: 75.25 3rd Qu.:57.25 3rd Qu.:168.2
## Max. :100.00 Max. :69.00 Max. :180.2
## weight income
## Min. :34.40 Min. : 219695
## 1st Qu.:50.15 1st Qu.: 1463506
## Median :57.25 Median : 2478078
## Mean :57.43 Mean : 3752153
## 3rd Qu.:64.85 3rd Qu.: 4289740
## Max. :85.60 Max. :20619415とする。
変数の視覚化
以下では、図を描くために ggplot2 というパッケージを使う。 インストール済みでない場合は、まず
install.packages("ggplot2")とする。.Rprofile 等で普段使うレポジトリを指定していない場合、どのミラーからファイルをダウンロードするか訊かれるので、自分に近いところを選ぶ。
インストールが済んだら、library() を使ってggplot2 パッケージを読み込む（require() でもよい）。
library("ggplot2")
## 次の行は、Mac ユーザ用。Windows ユーザは実行しないように
theme_set(theme_gray(base_size = 12, base_family = "HiraKakuProN-W3"))ggplot2 の詳しい使い方は第10回の授業で説明するので、今回はとりあえず以下のコードを真似て（そのままコピペして）図を作ってみよう（例として使われている変数だけでなく、他の変数でも試してみること）。
ヒストグラム (histogram)
まず、身長 (height) のヒストグラムを作る。
p1 <- ggplot(myd, aes(height))
p1 + geom_histogram()## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
作った図を保存したいときは、ggsave()を使う たとえば、
ggsave('height-histogram.pdf')とすると、最後に作った図が height-histogram.pdf という名前のPDFファイルとして保存される（PNGファイルとして保存したいときは“.pdf” を “.png” に変えればよい）。 RStudio では、右下にあるウィンドウのPlots タブに保存したい図を表示し、Export で保存することができる。
ちなみに、ggplot2 を使わない場合は、
hist(myd$height, xlab = "身長", ylab = "度数", main = "", col = "gray")
のようにヒストグラムが描ける （見栄え・美しさはともかく、ggplot2を使う場合と使わない場合でヒストグラムの見た目が違うのはなぜか考えよ）。
デフォルトでは縦軸が度数 (count, frequency) である。これを確率密度 (probability density) に変えてみよう。
p1 + geom_histogram(aes(y = ..density..)) ## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
ggplot2を使わない場合は、
hist(myd$height, probability = TRUE, col = "gray",
xlab = "身長", ylab = "確率密度", main ="")
とする。
ヒストグラムは階級幅によって見た目が変わるので、色々変えてみる必要がある。 幅を5に変えてみよう。
p1 + geom_histogram(binwidth = 5)
これでggplot2を使わない場合のヒストグラムと似たような形になった。
ここまで作ったヒストグラムは、男女のデータを一緒に扱っている。 男女を分けて描くとどうなるだろうか。
p1 + geom_histogram(binwidth = 5) + facet_wrap( ~ sex)
箱ひげ図 (box [and whiskers] plot)
続いて、身長の箱ひげ図を男女別に描いてみよう。
p2 <- ggplot(myd, aes(sex, height))
p2 + geom_boxplot() 
散布図 (scatter plot)
今度は、身長と体重 (weight) の関係を散布図を使って確認してみよう。
p3 <- ggplot(myd, aes(height, weight))
p3 + geom_point()
ここでも、男女を色分けして区別しよう。
p3 + geom_point(aes(color = sex))
次に、身長と体重の関係を、（男女別に）直線で示してみよう。
last_plot() + geom_smooth(method = 'lm', se = FALSE, aes(color = sex))
最後に、所得 (income) も同時に図示してみよう。
p3 + geom_point(aes(color = sex, size = income))
統計量
代表値
変数（ベクトル）\(x\)の平均値（算術平均, mean）は \(\bar{x}\) と表し（「エックスバー」と読む）、 \[ \bar{x} = \frac{\sum_{i=1}^n x_i}{n}\] である。 Rでは、
x <- rnorm(20, mean = 10, sd = 4) ## 変数xを正規分布からの無作為抽出で作る（後述）
sum(x) / length(x)## [1] 10.20578とすればよい。 しかし、平均値のような基本的な統計量を求める関数はあらかじめ用意されている。 平均値は、mean() で求められる。
mean(x)## [1] 10.20578これで、先ほど求めたものと同じ平均値が得られる。
中央値はmedian()、 最小値はmin()、 最大値はmax()で求める。
median(x)## [1] 9.133965min(x)## [1] 1.018219max(x)## [1] 18.8027分位数
また、分位数 (quantile) はquantile()で求める。 百分位数 (percentile) や四分位数 (quartile) を求めるときもこの関数を使えばよい。
## 50% quantile = 50th percentile = median
quantile(x, .5)## 50%
## 9.133965## 25th percentile = 1st quartile と 75th pecentile = 3rd quartile を同時に求める
quantile(x, c(.25, .75)) ## 25th percentile と 75th pecentile を同時に求める## 25% 75%
## 7.360155 13.858118したがって、四分位範囲 (IQR) は、
quantile(x, .75) - quantile(x, .25)## 75%
## 6.497963で求められる。用意されている関数を使い、
IQR(x)## [1] 6.497963でもよい。 五数要約は、
quantile(x, c(0, .25, .5, .75, 1)) ## 五数要約：四分位数を使うとき## 0% 25% 50% 75% 100%
## 1.018219 7.360155 9.133965 13.858118 18.802697とすればよい。 第1四分位数と第3四分位数の代わりに下側ヒンジと上側ヒンジ (Tukey 1977. Exploratory Data Analysis.) を求めたいときは、fivenum() を使う。
fivenum(x) # 五数要約：ヒンジを使うとき## [1] 1.018219 7.279039 9.133965 14.047347 18.802697分散、標準偏差
変数 \(x\) の不偏分散 (unbiased variance) は\(u_x^2\) と表し、 \[u_x^2 = \frac{\sum_{i=1}^n (x_i - \bar{x})^2}{n-1}\] である。 これをR で求めると、
sum((x - mean(x))^2) / (length(x) - 1) ## 計算の順序（カッコの付け方）に注意## [1] 21.67838とすればよい。 分散は最も重要な統計量であり、当然関数が用意されているので、var() で求めることができる。
var(x)## [1] 21.67838標準偏差 (standard deviation) は分散の平方根なので、
sqrt(var(x))## [1] 4.656005で求められるが、sd() を利用することもできる。
sd(x)## [1] 4.656005母集団を扱う場合
ベクトル\(x\) が標本ではなく母集団のとき、平均は\(\mu\)、分散は\(\sigma^2\) で表される。 母集団の平均値は、標本の場合と同様、
mean(x)## [1] 10.20578で求められる。 しかし、母集団の分散は \[\sigma^2 = \frac{\sum_{i=1}^N (x_i - \mu)^2}{N}\] なので、var() をそのまま使うことはできない。代わりに、
var(x) * (length(x) - 1) / (length(x))## [1] 20.59446とすればよい。 自分で関数を定義することもできる。
pop_var <- function(population) {
return(sum((population - mean(population))^2) / length(population))
}
pop_var(x)## [1] 20.59446確率と確率分布
様々な確率分布
R には無作為抽出を行うための関数が用意されている。
たとえば、最小値 \(a\)、最大値 \(b\) の連続一様分布 (uniform distribution) から無作為に\(n\) 個の値を取り出すには、 runif(n, a, b) とする。 試しに、\(n = 10, a = 0, b = 1\) でやってみると、
runif(10, 0, 1)## [1] 0.01616203 0.94178288 0.65463629 0.75084819 0.21746764 0.25498340
## [7] 0.42603532 0.79677889 0.99156610 0.77797770となる。 無作為抽出（正確には疑似乱数だが）なので、取り出される数字は毎回異なる。 したがって、もう一度同じコマンドを使うと、
runif(10, 0, 1)## [1] 0.3930013 0.2799974 0.4877927 0.3021058 0.1998018 0.2536563 0.8137272
## [8] 0.3773307 0.9653664 0.1986425のように、先ほどとは異なる数字が得られる。 乱数を使って論文やレポートを書くとき、同じ乱数を再現したい場合がある。そんなときは、set.seed() で乱数の種 (seed) を指定してやればよい。
set.seed(2016-04-20) # 種を指定
runif(10, 0, 1)## [1] 0.6248014 0.8214425 0.8061282 0.5817437 0.4483602 0.7123616 0.2008463
## [8] 0.6576807 0.6397562 0.8438949set.seed(2016-04-20) # もう一度同じ種を使う
runif(10, 0, 1) # 同じ数字が得られる## [1] 0.6248014 0.8214425 0.8061282 0.5817437 0.4483602 0.7123616 0.2008463
## [8] 0.6576807 0.6397562 0.8438949一様分布以外の分布から無作為抽出するときも r分布名() という関数を使う。 正規分布ならrnorm()、\(t\) 分布はrt()、カイ二乗分布はrchisq() という具合である。 ただし、指定しなければならない母数（parameters) は分布によって異なる。 既に見たように、一様分布では最小値 (min) と最大値 (max) を指定するが、正規分布では平均値 (mean) と標準偏差 (sd)、\(t\) 分布やカイ二乗分布では自由度 (df) を指定する。
たとえば、平均\(\mu\)、分散\(\sigma^2\) の正規分布から無作為kに \(n\) 個の値を取り出すには、rnorm(n, mean = mu, sd = sigma) とする。平均が10、分散が4の正規分布から100個の値を取り出したいなら、
x <- rnorm(100, 10, 2)とする。
r分布名() の“r” を“d”に変えると確率密度 \(f(x)\) が、“p” に変えると累積分布関数\(F(x)\) が得られる。