Rによるグループ別 基本統計量の出し方
勉強し始めたばかりのRについて、色々試行錯誤したので記録しておきます。
とあるレポートを作成するのに、データセットの基本統計量を出して概要を見ようとしたらR初心者には意外に難易度が高かった・・・。
基本統計量が出せても見づらかったりしたので、定番のsummary()の他、関数を色々比較してみました。
--------------------------------------
基本統計量とは
データセットに含まれる量的変数の、平均や分散、最小値／最大値、中央値、データの個数等のことです。
Excelだと、アドインのデータ分析を使うと簡単に下の図のような一覧表ができたりします。
もちろん、同じようなことはRでもsummary()関数を使うことで可能です。
summary()は、R初心者が最初に覚える関数のひとつなのではないでしょうか。使うと下の図のように出力されます。（画面はR Studioのコンソール）
ここで問題が2点。
１）グループ別に分けて基本統計量を出したい
２）基本統計量の表に、各項目名（平均とか標準偏差とかMin.とかMedianとか…）が繰り返し出てくるのがうっとおしい
グループに分けたい…！
今回利用しているのは "iris" というRのパッケージの中に標準で入っているデータセットです。これにはユリの品種、花弁の長さ・幅、萼片の長さ・幅のデーが品種ごと50個、全部で150個（N=150）入っています。
このデータセットに対して、なにも手を加えずにsummary関数を使うと、当然のことながら出てくる基本統計量は、品種問わず全てのデータをつかったものになります。とはいえ、データ分析するからには傾向が違いそうなグループごとに分けてみたくなるのは当然でしょう。
実はこのグループ分け、良いやり方を知らないと読み込むデータセットそのものを分割しないといけなかったり、意外にめんどくさかったりするのです。
項目名がうっとおしい…！
Excelのアドインと、Rのsummary、どちらもデータセットの複数の変数について1回の操作で基本統計量を出してくれます。楽なのですが、項目名を一つ一つ繰り返してくれるのは冗長だし、読みづらい。Excelなら出した後に列を削除するのも簡単ですが、それでも手間です。
その１：is.element()でデータフレーム抽出してグループ分け
Rに読み込んだcsvやExcelのデータセットは、データフレームというデータ型になっているそうです。まずはこのデータフレームをirisのSpeciesでグループ分けしてみます。
データセットの列に含まれている文字列データの種類を見るにはlebels()が使える模様。実行してみると、こんな感じで列に含まれるデータを返してくれます。データ整理から始めないといけないときは、タイプミスとか見つけるのに使えるかも・・・。
> levels(iris$Species)
[1] "setosa" "versicolor" "virginica" Speciesの中身が分かったところで、この3つについてデータフレームを分けてみます。
#Speciesごとにデータフレームを抽出
sp0<-iris[is.element(iris$Species,"setosa"),,drop=F]
sp1<-iris[is.element(iris$Species,"versicolor"),,drop=F]
sp2<-iris[is.element(iris$Species,"virginica"),,drop=F]
#とりあえず、グループ分けしたそれぞれのデータフレームにsummary関数を使う
summary(sp0[,c(2,4)])
summary(sp1[,c(2,4)])
summary(sp2[,c(2,4)])summary関数の引数では、データフレームの内、2列目、4列目のみ集計対象として指定してみています。Excelは隣り合う列だけしか集計できないので、事前に整形が必要なのですが、Rはその必要がない点は楽ですね。
出てきた結果は下図。グループ別集計ができてバンザイなのですが、やはり項目名がうっとおしい。。。
> summary(sp0[,c(2,4)])
Sepal.Width Petal.Width
Min. :2.300 Min. :0.100
1st Qu.:3.200 1st Qu.:0.200
Median :3.400 Median :0.200
Mean :3.428 Mean :0.246
3rd Qu.:3.675 3rd Qu.:0.300
Max. :4.400 Max. :0.600
> summary(sp1[,c(2,4)])
Sepal.Width Petal.Width
Min. :2.000 Min. :1.000
1st Qu.:2.525 1st Qu.:1.200
Median :2.800 Median :1.300
Mean :2.770 Mean :1.326
3rd Qu.:3.000 3rd Qu.:1.500
Max. :3.400 Max. :1.800
> summary(sp2[,c(2,4)])
Sepal.Width Petal.Width
Min. :2.200 Min. :1.400
1st Qu.:2.800 1st Qu.:1.800
Median :3.000 Median :2.000
Mean :2.974 Mean :2.026
3rd Qu.:3.175 3rd Qu.:2.300
Max. :3.800 Max. :2.500 この項目名をまとめる方法を調べて、見つけたのがapply()関数。
一つの関数を複数オブジェクトにまとめて実行できるらしいです。
#Speciesごとにデータフレームを抽出。ここまではさっきと一緒。
sp0<-iris[is.element(iris$Species,"setosa"),,drop=F]
sp1<-iris[is.element(iris$Species,"versicolor"),,drop=F]
sp2<-iris[is.element(iris$Species,"virginica"),,drop=F]
#apply()で複数列のsummaryを一括実行＆出力
tmp<-apply(sp0[,c(2,4)], 2, summary)
print(tmp)
#summary(sp1[,c(2,4)])
tmp<-apply(sp1[,c(2,4)], 2, summary)
print(tmp)
#summary(sp2[,c(2,4)])
tmp<-apply(sp2[,c(2,4)], 2, summary)
print(tmp)結果はこちら。ねらいどおりに項目名は先頭列に集約されました。これで列が増えても、シンプルに表示できるはず。
> #apply()で複数列のsummaryを一括実行＆出力
> tmp<-apply(sp0[,c(2,4)], 2, summary)
> print(tmp)
Sepal.Width Petal.Width
Min. 2.300 0.100
1st Qu. 3.200 0.200
Median 3.400 0.200
Mean 3.428 0.246
3rd Qu. 3.675 0.300
Max. 4.400 0.600
>
> #summary(sp1[,c(2,4)])
> tmp<-apply(sp1[,c(2,4)], 2, summary)
> print(tmp)
Sepal.Width Petal.Width
Min. 2.000 1.000
1st Qu. 2.525 1.200
Median 2.800 1.300
Mean 2.770 1.326
3rd Qu. 3.000 1.500
Max. 3.400 1.800
>
> #summary(sp2[,c(2,4)])
> tmp<-apply(sp2[,c(2,4)], 2, summary)
> print(tmp)
Sepal.Width Petal.Width
Min. 2.200 1.400
1st Qu. 2.800 1.800
Median 3.000 2.000
Mean 2.974 2.026
3rd Qu. 3.175 2.300
Max. 3.800 2.500その２：describeBy()関数（psychパッケージ）使用
グループ分けできたことに満足して、出力されたsummary()の結果を見ていたのですが・・・
あれ？データの個数とか分散とか、標準偏差とかがない。
summary()って、この辺りは出力してくれないんですね・・・。せめてデータの個数は欲しい。
それで調べてみたら、psychパッケージのdescribeBy()関数が良さそう。
Rをインストールした標準の状態のライブラリには入っていないので、パッケージをインストールして来て、使用する際にはlibrary()で読み込む必要があります。だけど、逆に言えばそれだけ。
library(psych)
#引数でデータセットの他、グループ分けの元になる変数の入った列を指定する
describeBy(iris,group = iris$Species)実行結果はこんな感じで、かなり情報量多めです。
グループの仕分け元になったデータを表の冒頭に示してくれるのが、とても良い感じです。
ただ、やはり、もう少しシンプルで良いのだけど…。
> library(psych)
> describeBy(iris,group = iris$Species)
Descriptive statistics by group
group: setosa
vars n mean sd median trimmed mad min max range skew kurtosis se
Sepal.Length 1 50 5.01 0.35 5.0 5.00 0.30 4.3 5.8 1.5 0.11 -0.45 0.05
Sepal.Width 2 50 3.43 0.38 3.4 3.42 0.37 2.3 4.4 2.1 0.04 0.60 0.05
Petal.Length 3 50 1.46 0.17 1.5 1.46 0.15 1.0 1.9 0.9 0.10 0.65 0.02
Petal.Width 4 50 0.25 0.11 0.2 0.24 0.00 0.1 0.6 0.5 1.18 1.26 0.01
Species* 5 50 1.00 0.00 1.0 1.00 0.00 1.0 1.0 0.0 NaN NaN 0.00
----------------------------------------------------------------
group: versicolor
vars n mean sd median trimmed mad min max range skew kurtosis se
Sepal.Length 1 50 5.94 0.52 5.90 5.94 0.52 4.9 7.0 2.1 0.10 -0.69 0.07
Sepal.Width 2 50 2.77 0.31 2.80 2.78 0.30 2.0 3.4 1.4 -0.34 -0.55 0.04
Petal.Length 3 50 4.26 0.47 4.35 4.29 0.52 3.0 5.1 2.1 -0.57 -0.19 0.07
Petal.Width 4 50 1.33 0.20 1.30 1.32 0.22 1.0 1.8 0.8 -0.03 -0.59 0.03
Species* 5 50 2.00 0.00 2.00 2.00 0.00 2.0 2.0 0.0 NaN NaN 0.00
----------------------------------------------------------------
group: virginica
vars n mean sd median trimmed mad min max range skew kurtosis se
Sepal.Length 1 50 6.59 0.64 6.50 6.57 0.59 4.9 7.9 3.0 0.11 -0.20 0.09
Sepal.Width 2 50 2.97 0.32 3.00 2.96 0.30 2.2 3.8 1.6 0.34 0.38 0.05
Petal.Length 3 50 5.55 0.55 5.55 5.51 0.67 4.5 6.9 2.4 0.52 -0.37 0.08
Petal.Width 4 50 2.03 0.27 2.00 2.03 0.30 1.4 2.5 1.1 -0.12 -0.75 0.04
Species* 5 50 3.00 0.00 3.00 3.00 0.00 3.0 3.0 0.0 NaN NaN 0.00その３：summarise_each()関数（dplyr パッケージ）使用
手ごろな関数がないのなら、自分の好きなように計算させて表っぽくまとめられればいいのでは…というわけでこちらの関数。
ただ出力結果が見づらい…。
#dplyr パッケージの関数使用
library(dplyr)
#Species列に入っているデータでグループ分け
data_grp<-group_by(iris[,c(2,4)],iris$Species)
#summarise_each()関数を使って、データ個数、最小、平均、中央値、最大、標準偏差を
#Sepal.Width,Petal.Width列のデータについて出しています。
#data_sum_eachに結果を入れているのは、R StudioのEnvironmentビューにデータとして表示させたいため。
#summarise_each()の実行だけでも結果は表示されます。
data_grp<-group_by(iris[,c(2,4)],iris$Species)
data_sum_each<-summarise_each(data_grp,list(length, min, mean, median, max,sd), Sepal.Width,Petal.Width)
print(data_sum_each)
実行結果はこちら。
悪くはないのだけど、列名が「変数名＋fn(n)」で並んでいくのは読みにくい＆スペースを取りすぎ・・・。
コンソール画面にも表示しきれないので、省略の注意が出てしまってます。
R Studioの場合、Environmentビューから結果を表として見られるので、こちらで確認するのが便利。
【2019.06.17 追記】列名表示を見やすくする方法を教えてもらいました…！
https://note.mu/kotoko_tyo/n/na2f5c9ca24ee
# A tibble: 3 x 13
`iris$Species` Sepal.Width_fn1 Petal.Width_fn1 Sepal.Width_fn2 Petal.Width_fn2 Sepal.Width_fn3 Petal.Width_fn3
<fct> <int> <int> <dbl> <dbl> <dbl> <dbl>
1 setosa 50 50 2.3 0.1 3.43 0.246
2 versicolor 50 50 2 1 2.77 1.33
3 virginica 50 50 2.2 1.4 2.97 2.03
# ... with 6 more variables: Sepal.Width_fn4 <dbl>, Petal.Width_fn4 <dbl>, Sepal.Width_fn5 <dbl>,
# Petal.Width_fn5 <dbl>, Sepal.Width_fn6 <dbl>, Petal.Width_fn6 <dbl>
その４：stargazer()関数（stargazer パッケージ）使用
ふたたび基本統計量を出してくれる関数を探してstargazer()関数。
グループ分けは別に行わないといけません。
library(stargazer)
#Speciesごとにデータフレームを抽出
sp0<-iris[is.element(iris$Species,"setosa"),,drop=F]
sp1<-iris[is.element(iris$Species,"versicolor"),,drop=F]
sp2<-iris[is.element(iris$Species,"virginica"),,drop=F]
#stargazerがほぼ理想！
stargazer(sp0[,c(1,2,3,4)],type="text",summary = TRUE)
stargazer(sp1[,c(1,2,3,4)],type="text",summary = TRUE)
stargazer(sp2[,c(1,2,3,4)],type="text",summary = TRUE)
実行結果がこちら。表の形と基本統計量の項目は、ほぼ最初に望んでいた通り！
課題があるとしたら、グループ分けの数が増えたときに、データフレーム抽出が煩雑かも…くらいでしょうか。その時は、グループ分けとデータフレーム抽出についてfor文でも書くことになるのでしょう。
> stargazer(sp0[,c(1,2,3,4)],type="text",summary = TRUE)
============================================================
Statistic N Mean St. Dev. Min Pctl(25) Pctl(75) Max
------------------------------------------------------------
Sepal.Length 50 5.006 0.352 4.300 4.800 5.200 5.800
Sepal.Width 50 3.428 0.379 2.300 3.200 3.675 4.400
Petal.Length 50 1.462 0.174 1.000 1.400 1.575 1.900
Petal.Width 50 0.246 0.105 0.100 0.200 0.300 0.600
------------------------------------------------------------
> stargazer(sp1[,c(1,2,3,4)],type="text",summary = TRUE)
============================================================
Statistic N Mean St. Dev. Min Pctl(25) Pctl(75) Max
------------------------------------------------------------
Sepal.Length 50 5.936 0.516 5 5.6 6.3 7
Sepal.Width 50 2.770 0.314 2.000 2.525 3.000 3.400
Petal.Length 50 4.260 0.470 3.000 4.000 4.600 5.100
Petal.Width 50 1.326 0.198 1.000 1.200 1.500 1.800
------------------------------------------------------------
> stargazer(sp2[,c(1,2,3,4)],type="text",summary = TRUE)
============================================================
Statistic N Mean St. Dev. Min Pctl(25) Pctl(75) Max
------------------------------------------------------------
Sepal.Length 50 6.588 0.636 4.900 6.225 6.900 7.900
Sepal.Width 50 2.974 0.322 2.200 2.800 3.175 3.800
Petal.Length 50 5.552 0.552 4 5.1 5.9 7
Petal.Width 50 2.026 0.275 1.400 1.800 2.300 2.500
------------------------------------------------------------結論
今回試してみた結果、私の感想としては次の通り。
手間をかけずにグループ別基本統計量を把握したい： describeBy()
カスタマイズして集計したい： group_by() ＋ summarise_each()
見やすい、シンプルな表が良い： stargazer()
…ここまでやってきてなんですが、Rを使う必要がなくて、データ数もそこそこの基本統計量なら、変に意地を張らずにExcelのアドインで処理して整形が楽かもとも思ったのは内緒です。
【2019.06.17 追記】
twitterで教えていただき、stby()関数とsummarise_each()関数の表示改善についてnoteに書きました！