このページでは、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  6581828
tail(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.133965
min(x)
## [1] 1.018219
max(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.8438949
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.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)\) が得られる。



授業の内容に戻る