計量経済学応用

分析内容の説明

R Markdown ファイルには、通常の文章や数式、写真やリンクなどを書き込むことができる。文章については、通常のテキストエディタやワープロソフトのように、単に文章を書く。

このとき Markdown という記法が使える。Markdown記法については、ココやココやココ を参照。

Rの命令だけ書いても、命令の意味を忘れてしまったり、自分がなぜその命令を書いたのか後でわからなくなったりするので、説明する文章を積極的に（たくさん）書き込もう。

Rの命令

R Markdownファイルの中で、Rの命令はコードチャンクと呼ばれる部分に書く。コードチャンクは、 ```{r}（コードチャンクの開始）と ```（コードチャンクの終了）に囲まれた部分である。RStudioでは、Ctrl + Alt + i （3つのキーを同時に押す）でコードチャンクが作れる。RStudio では、コードチャンクが正しくできるとその部分の背景の色が他と少し変わる。

コードチャンクに書いた命令を実行したいときは、実行したい行にカーソルを合わせ、Ctrl と Enter (Mac の場合は、Cmd と Return) を同時に押す。すると、命令が実行され、結果がR Markdownファイル上に表示される（設定によって、Consoleにのみ表示するように変更することもできる）。

このように、R Markdownファイルには文章とRのコードの両方を書き込める。最終的には、文章と分析結果をまとめて HTMLファイルやPDFファイルに出力することができる。よって、RStudioだけでレポートや論文を書くことが可能である。（So long, MS Word!）

もう少し詳しい説明は、Rマークダウン入門にある（他にも、ネット常に有益な情報がたくさんある）。

ヒストグラム (histogram)

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

p1 <- ggplot(myd, aes(x = 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 で保存することもできる。

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

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

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

ヒストグラムは階級幅によって見た目が変わるので、色々変えてみる必要がある。 幅を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(x = sex, y = height))
p2 + geom_boxplot() 

箱ひげ図だけでは、分布の形状がわからないので、バイオリンプロットを描いてみよう。

p2v <- p2 +  geom_violin()
print(p2v)

バイオリンプロットの上に、箱ひげ図を重ね描きしてみよう。

p2v + geom_boxplot(fill = "gray", width = 0.4)

箱ひげ図よりも情報量が多く、近年よく使われるようになってきた beeswarmplot （蜂群図）を描いてみよう。 追加のパッケージが必要なので、まずインストールして読み込む。

install.packages("ggbeeswarm")

library("ggbeeswarm")

geom_quasirandom() でbeeswarmplot を描く。

p2b <- p2 +  geom_quasirandom()
print(p2b)

箱ひげ図と併用してみよう。

p2 + geom_boxplot() + geom_quasirandom(aes(color = sex))

散布図 (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.88367

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

mean(x)

## [1] 10.88367

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

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

median(x)

## [1] 10.34832

min(x)

## [1] 4.225534

max(x)

## [1] 17.1569

分位数

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

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

##      50% 
## 10.34832

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

##       25%       75% 
##  8.840022 12.793482

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

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

##     75% 
## 3.95346

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

IQR(x)

## [1] 3.95346

でもよい。 五数要約は、

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

##        0%       25%       50%       75%      100% 
##  4.225534  8.840022 10.348323 12.793482 17.156903

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

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

## [1]  4.225534  8.783684 10.348323 13.492319 17.156903

分散、標準偏差

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

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

var(x)

## [1] 12.97088

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

sqrt(var(x))

## [1] 3.60151

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

sd(x)

## [1] 3.60151

母集団を扱う場合

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

mean(x)

## [1] 10.88367

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

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

## [1] 12.32233

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

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

## [1] 12.32233