2016年の社会調査士科目 D は最終的にブートストラップへの入門を目指している。この授業で扱う分析を行うための関数をいくつか紹介しておこう。ただし、事前に以下のページの内容を理解しておくこと。
記述統計は summary() 関数でも得られるが、以下の様な関数も知っておくといろいろ便利である。
mean(ベクトル) # ベクトルの平均値を返す
median(ベクトルの) # ベクトルの中央値を返す
以下が計算例
x1 <- c(2, 4, 7, 5, 3, 22)
mean(x1)## [1] 7.166667median(x1)## [1] 4.5mean() と median() はデフォルトでは欠損値 (NA) があると、NA を返し、平均値を計算してくれないので、na.omit() であらかじめ欠損値を除外したベクトルを作ってから平均などを計算するか、
mean(ベクトル, na.rm=TRUE)
median(ベクトル, na.rm=TRUE)
とすると、うまく計算できる。例えば以下のように。
x1.1 <- c(x1, NA)
x1.1## [1] 2 4 7 5 3 22 NAmean(x1.1)## [1] NAmean(na.omit(x1.1))## [1] 7.166667mean(x1.1, na.rm = TRUE)## [1] 7.166667\(\alpha\)パーセンタイルとは、ある変数の値を昇順に並べたときに、下位からちょうど\(\alpha\)パーセントのところにある値のことである。例えば、50 人の身長を調べて、背の低い人から高い人まで順に並べたときに、下から10番目あたりの値が 20パーセンタイルであり、下から45番目あたりの値が90パーセンタイルである。
パーセンタイルをパーセント単位ではなく、0.2, 0.9 といった比率で表す場合は分位点 (quantile) ともいう。これを計算する関数は、以下のとおり。
quantile(データのベクトル, どの分位を計算するかを指定するベクトル) # 指定した分位の分位点のベクトルを返す
例えば、20 パーセンタイルと 50 パーセンタイルと 90 パーセンタイルの計算例が以下である。
height <- round(rnorm(50, mean=165, sd=10),1) # 架空の50人分のデータを作る rnorm() については後述
height## [1] 146.8 153.4 169.0 168.7 164.7 164.9 169.0 163.9 158.9 166.5 168.8
## [12] 147.5 145.1 159.1 145.4 179.7 148.4 156.7 161.6 165.1 176.9 173.5
## [23] 150.5 148.3 164.3 166.5 178.8 181.4 168.0 184.5 168.6 152.5 161.9
## [34] 174.8 136.7 151.7 175.1 167.2 166.7 147.7 154.6 168.6 166.1 167.2
## [45] 164.5 180.8 174.9 160.3 148.0 171.4sort(height) # sort() はベクトルを昇順で並べ替え## [1] 136.7 145.1 145.4 146.8 147.5 147.7 148.0 148.3 148.4 150.5 151.7
## [12] 152.5 153.4 154.6 156.7 158.9 159.1 160.3 161.6 161.9 163.9 164.3
## [23] 164.5 164.7 164.9 165.1 166.1 166.5 166.5 166.7 167.2 167.2 168.0
## [34] 168.6 168.6 168.7 168.8 169.0 169.0 171.4 173.5 174.8 174.9 175.1
## [45] 176.9 178.8 179.7 180.8 181.4 184.5quantile(height, c(0.2, 0.5, 0.9))## 20% 50% 90%
## 151.46 165.00 177.09ineq パッケージの ineq() でジニ係数は計算できる。
ineq(ベクトル) # ベクトルのジニ係数を返す
以下が計算例。
library(ineq)## Warning: package 'ineq' was built under R version 3.2.3ineq(x1)## [1] 0.4379845データからピアソンの積率相関係数（ふつう単に相関係数といえばこれのこと）を計算する場合、cor() を使う。
cor(データフレーム) # データフレーム内のすべての変数のあいだの相関係数を返す
以下は計算例。
head(USArrests) # USArrests というサンプルデータを使う## Murder Assault UrbanPop Rape
## Alabama 13.2 236 58 21.2
## Alaska 10.0 263 48 44.5
## Arizona 8.1 294 80 31.0
## Arkansas 8.8 190 50 19.5
## California 9.0 276 91 40.6
## Colorado 7.9 204 78 38.7summary(USArrests)## Murder Assault UrbanPop Rape
## Min. : 0.800 Min. : 45.0 Min. :32.00 Min. : 7.30
## 1st Qu.: 4.075 1st Qu.:109.0 1st Qu.:54.50 1st Qu.:15.07
## Median : 7.250 Median :159.0 Median :66.00 Median :20.10
## Mean : 7.788 Mean :170.8 Mean :65.54 Mean :21.23
## 3rd Qu.:11.250 3rd Qu.:249.0 3rd Qu.:77.75 3rd Qu.:26.18
## Max. :17.400 Max. :337.0 Max. :91.00 Max. :46.00cor(USArrests)## Murder Assault UrbanPop Rape
## Murder 1.00000000 0.8018733 0.06957262 0.5635788
## Assault 0.80187331 1.0000000 0.25887170 0.6652412
## UrbanPop 0.06957262 0.2588717 1.00000000 0.4113412
## Rape 0.56357883 0.6652412 0.41134124 1.0000000これもデータに NA が含まれると計算してくれないので、na.omit() でリストワイズするか、引数に use=“complete” としてリストワイズするか、use=“pairwise”、としてペアワイズすればよい。
d0 <- matrix(c(NA, -10, -30, NA, 200, NA, NA, 100), 2, 4) # NA を含む架空の2行分のデータ
colnames(d0) <- colnames((USArrests))
d1 <- rbind(USArrests, d0) # NA を含むデータを2行追加
cor(d1) # NA があるので、相関係数は計算してくれない## Murder Assault UrbanPop Rape
## Murder 1 NA NA NA
## Assault NA 1 NA NA
## UrbanPop NA NA 1 NA
## Rape NA NA NA 1cor(d1, use="complete") # リストワイズ## Murder Assault UrbanPop Rape
## Murder 1.00000000 0.8018733 0.06957262 0.5635788
## Assault 0.80187331 1.0000000 0.25887170 0.6652412
## UrbanPop 0.06957262 0.2588717 1.00000000 0.4113412
## Rape 0.56357883 0.6652412 0.41134124 1.0000000cor(na.omit(d1)) # 同じくリストワイズ## Murder Assault UrbanPop Rape
## Murder 1.00000000 0.8018733 0.06957262 0.5635788
## Assault 0.80187331 1.0000000 0.25887170 0.6652412
## UrbanPop 0.06957262 0.2588717 1.00000000 0.4113412
## Rape 0.56357883 0.6652412 0.41134124 1.0000000cor(d1, use="pairwise") # ペアワイズ## Murder Assault UrbanPop Rape
## Murder 1.00000000 0.8018733 0.06957262 -0.06886945
## Assault 0.80187331 1.0000000 -0.10826583 0.66524123
## UrbanPop 0.06957262 -0.1082658 1.00000000 0.41134124
## Rape -0.06886945 0.6652412 0.41134124 1.00000000相関係数が多すぎて見にくい場合は、以下のようにして計算したい変数の組み合わせを指定することもできる。
cor(表側に配する変数のデータフレーム（またはベクトル）, 表頭配する変数のデータフレーム（またはベクトル）)
例えば以下のように。
cor(USArrests$UrbanPop, USArrests[, -3]) # USArrests[, -3] は USArrests から 3列目 (UrbanPop) を除いたデータフレーム## Murder Assault Rape
## [1,] 0.06957262 0.2588717 0.4113412二項分布を使った比率の検定には以下の様な専用の関数がある。
binom.test(成功回数, 試行回数, p = 帰無仮説の比率) # p はデフォルトでは 0.5 で両側検定の結果を返す
例えば、ある大学の学生から 50人無作為抽出したらアルバイト経験者が 45人いたとする。全国平均の 75% （架空）と有意に異なるかどうか検定した計算が以下である。
binom.test(45, 50, p=0.75)##
## Exact binomial test
##
## data: 45 and 50
## number of successes = 45, number of trials = 50, p-value = 0.01331
## alternative hypothesis: true probability of success is not equal to 0.75
## 95 percent confidence interval:
## 0.7818646 0.9667249
## sample estimates:
## probability of success
## 0.9両側検定の結果、 5% 水準で帰無仮説を棄却できる。比率の 95% 信頼区間もいっしょに求めてくれる。
t.test(ベクトル, mu = 帰無仮説の平均値) # ベクトルの平均値の検定結果を返す。デフォルトの帰無仮説は mu = 0
例えばパーセンタイルの節で使った身長のデータに関して、帰無仮説を「平均身長 = 160」として両側検定した結果は以下のとおり。
summary(height)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 136.7 153.7 165.0 163.1 169.0 184.5t.test(height, mu = 160)##
## One Sample t-test
##
## data: height
## t = 1.9634, df = 49, p-value = 0.05529
## alternative hypothesis: true mean is not equal to 160
## 95 percent confidence interval:
## 159.927 166.281
## sample estimates:
## mean of x
## 163.104p 値が 0.001 より小さいので、0.1% 水準で帰無仮説は棄却できる。
対応がある場合、2つのベクトルの差をとって、それを上と同じように t 検定するという方法が比較的簡単である。例えば、7組の夫婦の年齢を調べて、夫の年齢と妻の年齢の平均値が等しいかどうか検定したい場合、以下のようにする。
husband <- c(20, 31, 46, 57, 72, 80, 99) # 夫の年齢
wife <- c(19, 33, 46, 50, 70, 75, 92) # 妻の年齢
dif.couple <- husband - wife # 夫婦年齢の差
dif.couple## [1] 1 -2 0 7 2 5 7t.test(dif.couple) # 帰無仮説は「差が 0」 なので mu = は省略して良い##
## One Sample t-test
##
## data: dif.couple
## t = 2.1401, df = 6, p-value = 0.07614
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.409565 6.123851
## sample estimates:
## mean of x
## 2.857143p値が 0.05 以上なので、帰無仮説は棄却できない。以下のようにしても対応のある平均値の差の検定はできる。
t.test(ベクトル1, ベクトル2, paired = TRUE)
t.test(husband, wife, paired = TRUE)##
## Paired t-test
##
## data: husband and wife
## t = 2.1401, df = 6, p-value = 0.07614
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.409565 6.123851
## sample estimates:
## mean of the differences
## 2.857143対応がない平均値の差の検定の場合は、上の “, paired = TRUE” という引数を削除すればよい。
t.test(ベクトル1, ベクトル2)
例えば A校と B校からそれぞれ、4人と 5人生徒を無作為抽出して、平均身長の差を検定したい場合、以下のようにすればよい。
set.seed(1986)
A <- round(rnorm(4, mean=140, sd=5), 1) # A校の架空のデータを作成
set.seed(777)
B <- round(rnorm(5, mean=160, sd=8), 1) # B 〃
A## [1] 139.8 141.4 141.3 135.2B## [1] 163.9 156.8 164.1 156.8 173.1t.test(A, B)##
## Welch Two Sample t-test
##
## data: A and B
## t = -7.0384, df = 5.6771, p-value = 0.0005233
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -31.80402 -15.22598
## sample estimates:
## mean of x mean of y
## 139.425 162.940上の例は、A校と B校のデータが別々のベクトルとして引数に指定されていたが、両者が一つのデータフレームにまとめられていることも多い。そのような場合、以下の様な引き数の指定の仕方が便利である。
t.test(平均値を計算する変数 ~ グループ分けする変数, data = データフレーム名)
例えば以下のように
school <- c(rep("A", 4), rep("B", 5))
schools.height <- data.frame(school, height=c(A, B))
schools.height # こういうかたちにデータがまとめられていることも多い。## school height
## 1 A 139.8
## 2 A 141.4
## 3 A 141.3
## 4 A 135.2
## 5 B 163.9
## 6 B 156.8
## 7 B 164.1
## 8 B 156.8
## 9 B 173.1t.test(height ~ school, data=schools.height)##
## Welch Two Sample t-test
##
## data: height by school
## t = -7.0384, df = 5.6771, p-value = 0.0005233
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -31.80402 -15.22598
## sample estimates:
## mean in group A mean in group B
## 139.425 162.940相関係数を検定したいときは、以下の関数を使う。
cor.test(ベクトル1, ベクトル2)
または以下のように指定しても良い
cor.test(~ 変数1 + 変数2, data = データフレーム名)
以下は使用例。
cor.test(USArrests$UrbanPop, USArrests$Assault)##
## Pearson's product-moment correlation
##
## data: USArrests$UrbanPop and USArrests$Assault
## t = 1.8568, df = 48, p-value = 0.06948
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.02098836 0.50111118
## sample estimates:
## cor
## 0.2588717cor.test( ~ UrbanPop + Assault, data = USArrests)##
## Pearson's product-moment correlation
##
## data: UrbanPop and Assault
## t = 1.8568, df = 48, p-value = 0.06948
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.02098836 0.50111118
## sample estimates:
## cor
## 0.2588717p 値が 0.05 以上なので、帰無仮説は棄却できない。相関係数の 95% 信頼区間も計算してくれている。
クロス表を行列で用意して、以下のようにすればよい。
chisq.test(行列)
以下は使用例。
UCBAdmissions # UC Berkely の 6学部の受験者の性別と合否（1973年）## , , Dept = A
##
## Gender
## Admit Male Female
## Admitted 512 89
## Rejected 313 19
##
## , , Dept = B
##
## Gender
## Admit Male Female
## Admitted 353 17
## Rejected 207 8
##
## , , Dept = C
##
## Gender
## Admit Male Female
## Admitted 120 202
## Rejected 205 391
##
## , , Dept = D
##
## Gender
## Admit Male Female
## Admitted 138 131
## Rejected 279 244
##
## , , Dept = E
##
## Gender
## Admit Male Female
## Admitted 53 94
## Rejected 138 299
##
## , , Dept = F
##
## Gender
## Admit Male Female
## Admitted 22 24
## Rejected 351 317tab1 <- apply(UCBAdmissions, c(1,2), sum) # 6学部分を合併したクロス表を作成
tab1## Gender
## Admit Male Female
## Admitted 1198 557
## Rejected 1493 1278round(prop.table(tab1, 2), digits=2) # 合格率を男女別に計算。男性の方が合格率が高く見えるが？## Gender
## Admit Male Female
## Admitted 0.45 0.3
## Rejected 0.55 0.7chisq.test(tab1)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tab1
## X-squared = 91.61, df = 1, p-value < 2.2e-16ちなみに、6学部あわせると男性の方が有意に合格率が高くなるが、これは男性の方が競争率の低い学部を受験する傾向があるからで、学部別に見ると、合格率は男女でほぼ同じか、女性のほうが高い。
round(prop.table(UCBAdmissions, c(2,3))[1,,], digits=2) # 煩雑なので不合格率を省略して学部別に男女の合格率を計算## Dept
## Gender A B C D E F
## Male 0.62 0.63 0.37 0.33 0.28 0.06
## Female 0.82 0.68 0.34 0.35 0.24 0.07度数分布表の検定も chisq.test() を使えばよい。
chisq.test(度数分布表のベクトル, p=帰無仮説の比率のベクトル) # 帰無仮説のデフォルトはすべてのセルが等比率
例えば、上の UCBAdmissions のデータで、学部によって受験生の数が有意に異なるのかどうか知りたいとしよう。その場合、以下のように計算すればよい。
tab2 <- apply(UCBAdmissions, 3, sum) # 6学部の受験者数の度数分布表
tab2## A B C D E F
## 933 585 918 792 584 714chisq.test(tab2)##
## Chi-squared test for given probabilities
##
## data: tab2
## X-squared = 158.34, df = 5, p-value < 2.2e-16以下のようにすれば、関数のヘルプが見られる。
?関数名
例えば、mean() のヘルプを見たければ “?mean” とすればよい。リンクを辿って mean() のヘルプを見て欲しい。ヘルプの構成はほぼ決まっており、以下のようになっている。
mean(x, …)
5 Default S3 method:
mean(x, trim = 0, na.rm = FALSE, …)
とあり、用法は2つある。上は x という引き数だけを指定する用法であり、下は x, trim, na.rm という3つの引き数が示されている。trim = 0 となっているのは、trim のデフォルトの値は 0 であり、 na.rm = FALSE となっているのは、na.rm のデフォルトの値は FALSE であることを示している。
Arguments: 引き数の簡単な解説。関数が計算している事柄について読み手がよく知っていることを前提に書かれていることが多いため、何を書いてあるのかわからないことも結構あるが、関数について理解する手がかりにはなる。mean() の場合、3つの引き数の意味が簡単に解説してある。x については、An R object. Currently there are methods for numeric/logical vectors and date, date-time and time interval objects. Complex vectors are allowed for trim = 0, only. とある。
Details: 関数についてもう少し詳しい解説。mean() のヘルプには Details は無い。
Value: その関数がどんな値を返すのかの簡単な解説。mean() の場合、trim にどんな値を指定するのかによって、どう返り値が変わるのか、ごく簡単に述べてある。trim = 0 （デフォルト）の場合はただの平均値を返すが、例えば trim=0.2 の場合は、x の 20 パーセンタイル以下の要素と 80 パーセンタイル以上の要素を除外したベクトルの平均値を返す、といった趣旨のことが書かれている。以下はそのことを具体的に示す計算例。
set.seed(1986)
x2 <- round(rnorm(10), 1) # 10ケースの架空のデータ
x2## [1] 0.0 0.3 0.3 -1.0 0.5 -0.7 0.8 -0.7 -1.6 -1.1mean(x2)## [1] -0.32x2s <- sort(x2) # x2 を昇順で並べ替え
x2s## [1] -1.6 -1.1 -1.0 -0.7 -0.7 0.0 0.3 0.3 0.5 0.8mean(x2s[2:9])## [1] -0.3mean(x2, trim=0.1) # 上と同じことをやってくれる ## [1] -0.3References: 関数がどんな計算をしているのか、くわしく知りたい人のための参考文献
See Also: 関連する関数のヘルプへのリンク
Examples: 関数の使用例。Arguments や Value を見てよくわからなくても Examples を良く検討してみるとどんな関数なのかわかることも多い。
ヘルプには分かりにくい記述も多いが、これを部分的にでよいので、読みこなすことが R プログラミング上達への道なので、チャレンジして欲しい。