比較政治学B

RStudio 入門
データの要約
- 変数の視覚化
- 統計量
確率と確率分布
- 様々な確率分布

このページでは、RStudio を使いながら統計学の基礎を復習する。 1回の授業で復習する範囲には限りがあるので、だいぶ大雑把な説明になっている。基礎に不安がある場合は、統計学の教科書を参照してほしい。

RStudio 入門

プロジェクトの作成

RStudio にはプロジェクト機能がある。この機能を使うと、プロジェクトの管理が容易になる。ここでは、「計量経済学応用」をプロジェクトの1つと考え、新規プロジェクトを作成してみよう。

以下のステップを踏めば、プロジェクトが作れる。

コンピュータの自分のファイルが保存できる場所に、比較政治学Bのためのディレクトリ（directory、フォルダ）を作る

ディレクトリ名はアルファベットと数字のみを利用する
ディレクトリ名の最初の文字はアルファベットにする
ディレクトリ名にスペースを使わない
例: comparativeB

上部のメニューで、File -> New Project を選ぶ
Existing Directory（既存のディレクトリ）を選ぶ
browse（閲覧）を押して、自分が作ったディレクトリ（フォルダ）を選び、右下の Create Project をクリックする

これで新しいプロジェクトができる。プロジェクト名（自分で作ったフォルダの名前がそのまま使われる）は RStudioの右上に表示される。

次回以降、このプロジェクトを開くには、File -> Open Project でこのプロジェクトを選べばよい。

Rスクリプトの作成と利用

Rのコマンド（命令文）は、RStudio 左側（あるいは左下）の Console に直接打ち込むこともできる。しかし、通常はそのような使い方はしない。代わりに、Rの命令が書かれたファイルを別に作り、その中に命令を記入する。

RStudio で新しいRスクリプトを作るには、RStudio 上部のメニューで、File -> New File -> R Script の順番で選ぶ。そうすると、RStudio の左側のウィンドウが上下に2分割されるはずである。このとき、左下にはConsole が見える。左上に新たに開くのがRスクリプトである。

Rスクリプトができたら、名前をつけて保存しよう（名前の付け方はディレクトリ名の付け方と同じルールで）。このファイルにRの命令を書き込んでいく。基本的には、1つの行には1つの命令しか書かない。

このファイルに書いた命令を実行したいときは、実行したい部分を選択し、Command と Return (Windows の場合は、Control と Enter) を同時に押す。すると、命令がConsole に送られ、実行される。

コメントの利用

Rスクリプトには、Rに送る命令以外に、自分（あるいは他の人間）用のコメントを書き込むことができる。Rでコメントを書くときは、# という記号を使う（注意：# は半角！）。Rは、その行で # より後にあるものを無視する。

例えば、以下の3行をRで実行すると、2行目と3行目は無視される。

a <- c(4, 5, 3, 4, 6, 7, 1, 2, 9)
# a の標準偏差を求めたい
# 標準偏差は英語では standard deviation: sd() という関数を使う
sd(a)

## [1] 2.505549

スクリプトに命令だけ書いても、命令の意味を忘れてしまったり、自分がなぜその命令を書いたのか後でわからなくなったりするので、どんどんコメントを書き込もう。

データの要約

まず、今回の実習で利用するデータをダウンロードしよう。準備として、現在利用しているプロジェクト（上で作ったディレクトリ）の中に、data という名前のディレクトリを作ろう。ディレクトリはdir.crate() で作れる。このコマンドは一度しか実行しないので、Consoleに直接入力する。

dir.create("data")

次にデータセット fake-data-01.csv をダウンロードし、今作った data ディレクトリの中に保存する。

download.file(url = "https://git.io/fAy9a",
              destfile = "data/fake-data-01.csv")

このデータは CSVと呼ばれる形式なので、Rのreadr パッケージに含まれている read_csv()という関数を使ってデータセットを読み込むことができる。このデータセットを myd という名前で利用することにしよう（これはデータセットのファイル名を変えるのではなく、Rで使う上での呼び名を決めている）。

myd <- readr::read_csv("data/fake-data-01.csv")

## Parsed with column specification:
## cols(
##   id = col_integer(),
##   sex = col_character(),
##   age = col_integer(),
##   height = col_double(),
##   weight = col_double(),
##   income = col_integer()
## )

読み込んだデータを確認してみよう。まず、データセットの最上部または最下部にある6行分だけを表示する。

head(myd)   ## 最初の6行分

## # A tibble: 6 x 6
##      id sex     age height weight  income
##   <int> <chr> <int>  <dbl>  <dbl>   <int>
## 1     1 male     52   174    63.1 3475810
## 2     2 male     33   175.   70.2  457018
## 3     3 male     22   175    82.6 1627793
## 4     4 male     33   170.   81.8 6070642
## 5     5 male     26   167.   51.2 1083052
## 6     6 male     37   159.   57.8 2984929

tail(myd)   ## 最後の6行分

## # A tibble: 6 x 6
##      id sex      age height weight   income
##   <int> <chr>  <int>  <dbl>  <dbl>    <int>
## 1    95 female    21   165.   56.3  1339138
## 2    96 female    65   161.   46.8  6127136
## 3    97 female    45   161.   48.7  1062663
## 4    98 female    53   166.   64.2 10154200
## 5    99 female    43   158.   48.5  8287163
## 6   100 female    48   154.   42    1125390

データセットに含まれる変数名を確認したいときは、

names(myd)

## [1] "id"     "sex"    "age"    "height" "weight" "income"

とする。

データセットに含まれる観測数 (\(n\)) と変数の数を知りたいときは、

dim(myd)

## [1] 100   6

とする。最初の数字が\(n\)の数（データセットの行数）、2番目の数字が変数の数（列数）である（RStudio を使うと、右上のウィンドウの“Environment” というタブに、この情報が表示される）。

また、データの中身の確認には、dplyrパッケージのglimpse()も便利である。

dplyr::glimpse(myd)

## Observations: 100
## Variables: 6
## $ id     <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, ...
## $ sex    <chr> "male", "male", "male", "male", "male", "male", "male",...
## $ age    <int> 52, 33, 22, 33, 26, 37, 50, 30, 62, 51, 55, 36, 66, 42,...
## $ height <dbl> 174.0, 175.3, 175.0, 170.1, 167.4, 159.3, 173.3, 162.5,...
## $ weight <dbl> 63.1, 70.2, 82.6, 81.8, 51.2, 57.8, 68.6, 47.2, 68.2, 5...
## $ income <int> 3475810, 457018, 1627793, 6070642, 1083052, 2984929, 14...

データセットに含まれるすべての変数の基本的な統計量を確認したいときは、

summary(myd)

##        id             sex                 age            height     
##  Min.   :  1.00   Length:100         Min.   :20.00   Min.   :148.0  
##  1st Qu.: 25.75   Class :character   1st Qu.:36.00   1st Qu.:158.1  
##  Median : 50.50   Mode  :character   Median :45.00   Median :162.9  
##  Mean   : 50.50                      Mean   :45.96   Mean   :163.7  
##  3rd Qu.: 75.25                      3rd Qu.:57.25   3rd Qu.:170.2  
##  Max.   :100.00                      Max.   :70.00   Max.   :180.5  
##      weight          income        
##  Min.   :28.30   Min.   :  240184  
##  1st Qu.:48.95   1st Qu.: 1343679  
##  Median :59.95   Median : 2987818  
##  Mean   :59.18   Mean   : 4343425  
##  3rd Qu.:67.33   3rd Qu.: 6072696  
##  Max.   :85.60   Max.   :23505035

とする。

変数の視覚化

以下では、図を描くために ggplot2 というパッケージを使う。インストール済みでない場合（大学のコンピュータには既にインストール済み（のはず）なので、この部分は実行しなくてよい。）は、まず

install.packages("ggplot2")

とする。1度インストールすれば、次回からはこの過程をスキップしてよい。

インストールが済んだら、library() を使ってggplot2 パッケージを読み込む。

library("ggplot2")

ggplot2 の詳しい使い方は今後で説明するので、今回はとりあえず以下のコードを真似て（そのままコピペして）図を作ってみよう（例として使われている変数だけでなく、他の変数でも試してみること）。

ヒストグラム (histogram)

まず、身長 (height) のヒストグラムを作る。

p1 <- ggplot(myd, aes(height)) 
p1 + geom_histogram(color = "black")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

作った図を保存したいときは、ggsave()を使うたとえば、

ggsave('height-histogram.pdf')

とすると、最後に作った図が height-histogram.pdf という名前のPDFファイルとして保存される。 RStudio では、右下にあるウィンドウのPlots タブに保存したい図を表示し、Export で保存することもできる。

ちなみに、ggplot2 を使わない場合は、

hist(myd$height, xlab = "身長", ylab = "度数", main = "", col = "gray")

のようにヒストグラムが描ける（見栄え・美しさはともかく、ggplot2を使う場合と使わない場合でヒストグラムの見た目が違うのはなぜでしょう？）。

デフォルトでは縦軸が度数 (count, frequency) である。これを確率密度 (probability density) に変えてみよう。

p1 + geom_histogram(aes(y = ..density..), color = "black")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot2を使わない場合は、

hist(myd$height, probability = TRUE, col = "gray",
     xlab = "身長", ylab = "確率密度", main ="")

とする。

ヒストグラムは階級幅によって見た目が変わるので、色々変えてみる必要がある。幅を5に変えてみよう。

p1 + geom_histogram(binwidth = 5, color = "black")

これでggplot2を使わない場合のヒストグラムと似たような形になった。

ここまで作ったヒストグラムは、男女のデータを一緒に扱っている。男女を分けて描くとどうなるだろうか。

p1 + geom_histogram(binwidth = 5, color = "black") + 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] 11.48006

とすればよい。しかし、平均値のような基本的な統計量を求める関数はあらかじめ用意されている。平均値は、mean() で求められる。

mean(x)

## [1] 11.48006

これで、先ほど求めたものと同じ平均値が得られる。

中央値はmedian()、最小値はmin()、最大値はmax()で求める。

median(x)

## [1] 11.96396

min(x)

## [1] 3.933282

max(x)

## [1] 16.19784

分位数

また、分位数 (quantile) はquantile()で求める。百分位数 (percentile) や四分位数 (quartile) を求めるときもこの関数を使えばよい。

## 50% quantile = 50th percentile = median
quantile(x, .5)

##      50% 
## 11.96396

## 25th percentile = 1st quartile と 75th pecentile = 3rd quartile を同時に求める
quantile(x, c(.25, .75))  ## 25th percentile と 75th pecentile を同時に求める

##      25%      75% 
## 10.02478 13.65789

したがって、四分位範囲 (IQR) は、

quantile(x, .75) - quantile(x, .25)

##      75% 
## 3.633108

で求められる。用意されている関数を使い、

IQR(x)

## [1] 3.633108

でもよい。五数要約は、

quantile(x, c(0, .25, .5, .75, 1))   ## 五数要約：四分位数を使うとき

##        0%       25%       50%       75%      100% 
##  3.933282 10.024781 11.963964 13.657889 16.197839

とすればよい。第1四分位数と第3四分位数の代わりに下側ヒンジと上側ヒンジ (Tukey 1977. Exploratory Data Analysis.) を求めたいときは、fivenum() を使う。

fivenum(x)   # 五数要約：ヒンジを使うとき

## [1]  3.933282  9.340644 11.963964 13.814015 16.197839

分散、標準偏差

変数 \(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] 11.84859

とすればよい。 分散は最も重要な統計量であり、当然関数が用意されているので、var() で求めることができる。

var(x)

## [1] 11.84859

標準偏差 (standard deviation) は分散の平方根なので、

sqrt(var(x))

## [1] 3.442179

で求められるが、sd() を利用することもできる。

sd(x)

## [1] 3.442179

母集団を扱う場合

ベクトル\(x\) が標本ではなく母集団のとき、平均は\(\mu\)、分散は\(\sigma^2\) で表される。母集団の平均値は、標本の場合と同様、

mean(x)

## [1] 11.48006

で求められる。しかし、母集団の分散は \[\sigma^2 = \frac{\sum_{i=1}^N (x_i - \mu)^2}{N}\] なので、var() をそのまま使うことはできない。代わりに、

var(x) * (length(x) - 1) / (length(x))

## [1] 11.25616

とすればよい。自分で関数を定義することもできる。

pop_var <- function(population) {
    return(sum((population - mean(population))^2) / length(population))
}
pop_var(x)

## [1] 11.25616

確率と確率分布

様々な確率分布

R には無作為抽出を行うための関数が用意されている。

たとえば、最小値 \(a\)、最大値 \(b\) の連続一様分布 (uniform distribution) から無作為に\(n\) 個の値を取り出すには、 runif(n, a, b) とする。試しに、\(n = 10, a = 0, b = 1\) でやってみると、

runif(10, 0, 1)

##  [1] 0.98915274 0.01465155 0.48844010 0.78348320 0.56397424 0.74716996
##  [7] 0.65689909 0.31081817 0.72014599 0.52662381

となる。無作為抽出（正確には疑似乱数だが）なので、取り出される数字は毎回異なる。したがって、もう一度同じコマンドを使うと、

runif(10, 0, 1)

##  [1] 0.6514845 0.4081328 0.2644889 0.6752236 0.4468969 0.9418136 0.6389834
##  [8] 0.7637832 0.2623205 0.4273976

のように、先ほどとは異なる数字が得られる。乱数を使って論文やレポートを書くとき、同じ乱数を再現したい場合がある。そんなときは、set.seed() で乱数の種 (seed) を指定してやればよい。

set.seed(2018-04-12) # 種を指定
runif(10, 0, 1)

##  [1] 0.45755673 0.21033244 0.94678455 0.02105657 0.91747148 0.79369380
##  [7] 0.75890287 0.98959130 0.98801158 0.90748110

set.seed(2018-04-12) # もう一度同じ種を使う
runif(10, 0, 1)      # 同じ数字が得られる

##  [1] 0.45755673 0.21033244 0.94678455 0.02105657 0.91747148 0.79369380
##  [7] 0.75890287 0.98959130 0.98801158 0.90748110

おすすめのseed（種）は、ここで使ったようにその日の日付 (年-月-日)である。

一様分布以外の分布から無作為抽出するときも 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)\) が得られる。

授業の内容に戻る

比較政治学B

Rによる統計分析の基礎

矢内　勇生

2018-09-18

RStudio 入門

プロジェクトの作成

Rスクリプトの作成と利用

コメントの利用

データの要約

変数の視覚化

ヒストグラム (histogram)

箱ひげ図 (box [and whiskers] plot)

散布図 (scatter plot)

統計量

代表値

分位数

分散、標準偏差

母集団を扱う場合

確率と確率分布

様々な確率分布

比較政治学B

Rによる統計分析の基礎

矢内 勇生

2018-09-18

RStudio 入門

プロジェクトの作成

Rスクリプトの作成と利用

コメントの利用

データの要約

変数の視覚化

ヒストグラム (histogram)

箱ひげ図 (box [and whiskers] plot)

散布図 (scatter plot)

統計量

代表値

分位数

分散、標準偏差

母集団を扱う場合

確率と確率分布

様々な確率分布

矢内　勇生