データの要約1 (記述統計量)

データを入手したら，分析を行う前にまずデータの概観を把握することが重要。データセットにどのような変数が含まれているかを確認し，それぞれの変数の記述統計量 (平均値，分位点，分散など)を求め，必要であれば変数間の相関関係を調べてみる。これにより，欠損値や異常値が含まれていないかということや，どのような仮説を検証することができそうかなど，分析の方向性を決めることができる。

ここでは，Rに含まれているirisデータを例に，データの要約統計量を求める方法を説明する。

Rのパッケージ，tidyverseとstargazerを利用するのであらかじめインストールしておく。

まず，インストールしたパッケージを読み込む。

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.6
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.1     ✔ tibble    3.3.0
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.2.0     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(stargazer)

Please cite as: 

 Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
 R package version 5.2.3. https://CRAN.R-project.org/package=stargazer 

データの概観

練習としてRに組み込まれているサンプル・データirisを概観してみよう。まず，サンプル・データirisを読み込む。

data(iris)

irisというデータフレームが読み込まれる。irisデータには，植物のアヤメ150個体について，種類 (3分類)，萼片 (がくへん)の長さと幅，花弁 (かべん)の長さと幅の5つの変数が記録されている。このデータは，機械学習による分類などの練習用データとしてよく用いられているデータである。

変数名	説明
Sepal.Length	萼片の長さ
Sepal.Width	萼片の幅
Petal.Length	花弁の長さ
Petal.Width	花弁の幅
Species	アヤメの種類

データにどのような変数が含まれているかは，小規模データであればEnvironmentタブからスプレッド・シートを表示させれば把握できる。大規模なデータは，スプレッド・シートを表示させるだけでも時間がかかることがあるので，glimpse()関数を使うのが良い。

glimpse(iris)

Rows: 150
Columns: 5
$ Sepal.Length <dbl> 5.1, 4.9, 4.7, 4.6, 5.0, 5.4, 4.6, 5.0, 4.4, 4.9, 5.4, 4.…
$ Sepal.Width  <dbl> 3.5, 3.0, 3.2, 3.1, 3.6, 3.9, 3.4, 3.4, 2.9, 3.1, 3.7, 3.…
$ Petal.Length <dbl> 1.4, 1.4, 1.3, 1.5, 1.4, 1.7, 1.4, 1.5, 1.4, 1.5, 1.5, 1.…
$ Petal.Width  <dbl> 0.2, 0.2, 0.2, 0.2, 0.2, 0.4, 0.3, 0.2, 0.2, 0.1, 0.2, 0.…
$ Species      <fct> setosa, setosa, setosa, setosa, setosa, setosa, setosa, s…

irisデータは，行数 (観測数)が150，列数 (変数の数)が5であることがわかる。また，それぞれの変数のはじめのいくつかの観測値を見ることができるので，データのイメージがつかめる。<dbl>や<fct>はそれぞれの変数の型で，<dbl>は実数型 (double, numericと同じ)，<fct>はファクター型であることを表している。

変数の型について

変数の型には，実数型 (numeric/double)，ファクター型 (factor)のほかに，文字列型 (character)<chr>や論理型 (logical)<lgl>などがある。Rでは通常，あまり厳密に変数の型を気にする必要はない。ただし，変数が数値型だと思って演算を行おうとすると，文字列型になっていてエラーになるということはよくあるので。数値型か文字型かくらいは意識しよう。

要約統計表の作成

summary()関数

データ全体のイメージがつかめたら，次は記録されている変数の要約統計表を作成してみよう。summary()関数を用いれば，データフレームに含まれるすべての変数の要約統計量を表示させることができる。

summary(iris)

  Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
 Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
 1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
 Median :5.800   Median :3.000   Median :4.350   Median :1.300  
 Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
 3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
 Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500  
       Species  
 setosa    :50  
 versicolor:50  
 virginica :50  
                
                
                

数値型の変数については，最小値(Min.)，第一四分位(1st Qu.)，中央値(Median)，平均値(Mean)，第三四分位(3rd Qu.)，最大値(Max)が表示される (summary()関数では，分散や標準偏差は表示されない)。また，ファクター型の変数については，度数分布表が表示される。

特定の変数を取り出す

特定の変数についてのみ，要約統計量を求めたい場合は，select()関数を使えば良い。

summary(select(iris, Sepal.Length, Sepal.Width))

  Sepal.Length    Sepal.Width   
 Min.   :4.300   Min.   :2.000  
 1st Qu.:5.100   1st Qu.:2.800  
 Median :5.800   Median :3.000  
 Mean   :5.843   Mean   :3.057  
 3rd Qu.:6.400   3rd Qu.:3.300  
 Max.   :7.900   Max.   :4.400  

select()関数は，データフレームから指定した変数だけを取り出したデータフレームを作成する。第1引数には変数を取り出すデータフレームを，第2引数以降は取り出したい変数名をカンマで区切って指定する。たとえば，select(dataframe, a, b)はdataframe$aとdataframe$bの2つの変数が含まれるデータフレームとなる。逆にデータフレームから指定した変数を除外したいときには，第2引数以降で変数名に-をつけて指定する。たとえば，select(dataframe, -a)はdataframeから変数aを除いたデータフレームとなる。また，変数は列番号で指定することも可能である。たとえば，select(dataframe, 2:4)はdataframeから2列目から4列目までの変数を取り出したデータフレームとなる。

この処理は，パイプを使って以下のように書くこともできる。

iris |>
  select(Sepal.Length, Sepal.Width) |>
  summary()

  Sepal.Length    Sepal.Width   
 Min.   :4.300   Min.   :2.000  
 1st Qu.:5.100   1st Qu.:2.800  
 Median :5.800   Median :3.000  
 Mean   :5.843   Mean   :3.057  
 3rd Qu.:6.400   3rd Qu.:3.300  
 Max.   :7.900   Max.   :4.400  

パイプというのは，関数を|>でつなげることにより，前の関数の結果を次の関数の第1引数として渡す。複数の処理を一度に行うことができ，無駄な中間オブジェクトを作成する必要がないため便利。この例だと，irisというデータフレームが2行目のselect()関数の第1引数として渡されるので，2行目ではirisからSepal.Length, Sepal.Widthの二つの変数を取り出すことになる。さらに，その結果が3行目のsummary()関数に渡されるため，Sepal.Length, Sepal.Widthの要約統計量が表示される。

summarize()関数

より柔軟にカスタマイズして要約統計表を作成したい場合には，summarize()関数を用いる。summarize()関数は，結果をデータフレームで返し，たとえば以下のように使う。

summarize(iris,
  萼片長平均     = mean(Sepal.Length),
  萼片長標準偏差 = sd(Sepal.Length),
  花弁長平均     = mean(Petal.Length),
  花弁長標準偏差 = sd(Petal.Length)
)

  萼片長平均 萼片長標準偏差 花弁長平均 花弁長標準偏差
1   5.843333      0.8280661      3.758       1.765298

summarize()関数の第1引数には要約統計量を求めたいデータフレーム名を指定する。第2引数以降は，要約統計量の名前 = 関数(変数名)という形式で指定する。名前は出力される要約統計表データフレームの列名なので自由に指定すれば良い。関数の部分には，例として下の表のような関数が利用できる。変数名xには，第1引数として与えたデータフレーム内の変数名を指定すればよい (データフレーム名は必要ない)。

関数	説明
n()	観測数
mean(x)	xの平均値
var(x)	xの分散
sd(x)	xの標準偏差
max(x)	xの最大値
median(x)	xの中央値
min(x)	xの最小値
quantile(x, XX)	xのXX×100パーセンタイル, e.g. xの25パーセンタイルであればquantile(x,0.25)

summarize()関数は，後で述べるようにグループ別の集計のときに役に立つ。

stargazerパッケージを使う

レポートや論文に貼り付けるための要約統計表を作成するには，stargazerパッケージが便利。

stargazer(
  iris,
  type = "text",
  title = "要約統計表",
  digits = 2,
  summary.stat = c(
    "mean",
    "sd",
    "min",
    "p25",
    "median",
    "p75",
    "max"
  )
)

要約統計表
=============================================================
Statistic    Mean St. Dev. Min  Pctl(25) Median Pctl(75) Max 
-------------------------------------------------------------
Sepal.Length 5.84   0.83   4.30   5.10    5.80    6.40   7.90
Sepal.Width  3.06   0.44   2.00   2.80    3.00    3.30   4.40
Petal.Length 3.76   1.77   1.00   1.60    4.35    5.10   6.90
Petal.Width  1.20   0.76   0.10   0.30    1.30    1.80   2.50
-------------------------------------------------------------

stargazer()関数は，要約統計表や回帰分析の結果などを，良い感じの見た目で出力してくれる。

stargazer()関数で要約統計表を作成する場合は，第1引数としてデータフレーム (この例ではiris)を指定する。そのほかの引数は必要に応じて指定すれば良い。よく使うオプションだけ説明する。まず，typeは表の出力形式で，“text”(プレーンテキスト)，“html”，“latex”のいずれかを指定する。titleは表のタイトル，digitsは表の数値の小数点以下の桁数を指定する。summary.statsは表に含める要約統計量をベクトル形式で指定する。

要約統計量	指定方法
平均	mean
標準偏差	sd
最小値	min
最大値	max
XXパーセンタイル(分位点)	p90 (90パーセンタイル)

要約統計表を作成することができるパッケージはstargazerのほかにもgtパッケージなどさまざまなものがある。用途や好みに応じて使い分けよう。stargazer()関数は，とくにLatex形式で出力する際に便利だ。

相関表

次に，変数間の相関を求めよう。相関表を求めるにはcor()関数を使う。引数にはデータフレーム名を指定する。ただし，irisにファクター型のSpeciesという変数が含まれているため，そのままではエラーになる。そこで，irisに含まれる5つの変数のうち数値型の4つを取り出して相関表を作成したいので，ファクター型の変数Speciesだけを除外したiris2を作成し，cor()関数の引数に，iris2を指定して相関行列を求める。

iris2 <- select(iris, - Species)
cor(iris2)

             Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length    1.0000000  -0.1175698    0.8717538   0.8179411
Sepal.Width    -0.1175698   1.0000000   -0.4284401  -0.3661259
Petal.Length    0.8717538  -0.4284401    1.0000000   0.9628654
Petal.Width     0.8179411  -0.3661259    0.9628654   1.0000000

同じことをパイプを使って書くと次のようになる。

iris |>
  select(- Species) |>
  cor()

             Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length    1.0000000  -0.1175698    0.8717538   0.8179411
Sepal.Width    -0.1175698   1.0000000   -0.4284401  -0.3661259
Petal.Length    0.8717538  -0.4284401    1.0000000   0.9628654
Petal.Width     0.8179411  -0.3661259    0.9628654   1.0000000

ここでは，パイプの理解のためにこのようにしたが，これくらいの作業であればパイプを使わずに以下のように書いても良い。

cor(select(iris, - Species))

グループ分け

データをグループに分けて，グループごとに要約統計を求めたいということは，実際によくある。irisデータでは，アヤメはsetosa，versicolor，virginicaの3種類に分類されている。そこで，アヤメの種類ごとに萼片と花弁の長さ，幅の平均を求めてみよう。そのための方法はいくつか考えられる。

group_by()を使う

たぶんこれが標準的な方法。

group_by()関数を使うと，(Rの内部で)グループに分類されたデータフレームを作成することができる。第1引数には分類したいデータフレーム名，第2引数以降は分類を行うためのファクター型変数を指定する。複数のファクター型変数を指定して，階層的な分類を行うことも可能。グループに分類されたデータフレームにsummarize()関数を適用すると，グループごとの要約統計量が計算される。

iris_grouped <- group_by(iris, Species)

summarize(
  iris_grouped,
  萼片長平均     = mean(Sepal.Length),
  萼片長標準偏差 = sd(Sepal.Length),
  花弁長平均     = mean(Petal.Length),
  花弁長標準偏差 = sd(Petal.Length)
)

# A tibble: 3 × 5
  Species    萼片長平均 萼片長標準偏差 花弁長平均 花弁長標準偏差
  <fct>           <dbl>          <dbl>      <dbl>          <dbl>
1 setosa           5.01          0.352       1.46          0.174
2 versicolor       5.94          0.516       4.26          0.470
3 virginica        6.59          0.636       5.55          0.552

同じことを，パイプを使うと以下のように書ける。

iris |> 
  group_by(Species) |> 
    summarize(
      萼片長平均     = mean(Sepal.Length),
      萼片長標準偏差 = sd(Sepal.Length),
      花弁長平均     = mean(Petal.Length),
      花弁長標準偏差 = sd(Petal.Length)
    )

# A tibble: 3 × 5
  Species    萼片長平均 萼片長標準偏差 花弁長平均 花弁長標準偏差
  <fct>           <dbl>          <dbl>      <dbl>          <dbl>
1 setosa           5.01          0.352       1.46          0.174
2 versicolor       5.94          0.516       4.26          0.470
3 virginica        6.59          0.636       5.55          0.552

パイプのメリットの1つは，不要な中間オブジェクト (この場合は，iris_grouped)を作成せずに済むという点。もちろん，中間オブジェクトを作成せずに同じ結果を得るのは，以下のコードでも可能。

summarize(
  group_by(iris, Species),
  萼片長平均     = mean(Sepal.Length),
  萼片長標準偏差 = sd(Sepal.Length),
  花弁長平均     = mean(Petal.Length),
  花弁長標準偏差 = sd(Petal.Length)
)

# A tibble: 3 × 5
  Species    萼片長平均 萼片長標準偏差 花弁長平均 花弁長標準偏差
  <fct>           <dbl>          <dbl>      <dbl>          <dbl>
1 setosa           5.01          0.352       1.46          0.174
2 versicolor       5.94          0.516       4.26          0.470
3 virginica        6.59          0.636       5.55          0.552

パイプを使わないと，関数の引数にまた関数を使うということになるので，コードが見づらくなる。パイプは，作業の流れがわかりやすいという点でも有用。

xtabs()関数を使う

xtabs()関数を使うと，多次元分割表を柔軟に作成することができる。第1引数には形式を，第2引数にはデータフレーム名を指定する。形式は，~ x1 + x2 + ...かy ~ x1 + x2 + ...の形で指定する。たとえば，~ x1 + x2を指定すれば，x1, x2で表を分割し，各セルの観測数が得られる。y ~ x1 + x2を指定すれば，x1, x2で表を分割し，各セルごとのyの合計が得られる。

irisでは分類に用いることができる変数はSpeciesの1つしかないが，それでもxtabs()を適用できる。以下のコードでSpeciesごとのデータ数の表を作成できる。

xtabs(~ Species, iris)

Species
    setosa versicolor  virginica 
        50         50         50 

次に，以下のコードでSpeciesごとの萼片長の合計値の表を作成できる。

xtabs(Sepal.Length ~ Species, iris)

Species
    setosa versicolor  virginica 
     250.3      296.8      329.4 

平均値は合計値をデータ数で割ったものなので，Speciesごとの萼片長の平均値を以下のように求められる。

xtabs(Sepal.Length ~ Species, iris) / xtabs(~ Species, iris)

Species
    setosa versicolor  virginica 
     5.006      5.936      6.588 

フィルター

まず，データフレームからアヤメの種類ごとにデータを取り出して分割する方法。単純でわかりやすいが，分類が多くなると大変なので，要約統計量を求めるためだけであれば効率的な方法ではない。ただし，データフレームにフィルターをかけて特定の条件を満たす観測 (行)を取り出すことはよくあるので，やり方を覚えておこう。

データフレームから条件を満たす行を取り出すには，filter()関数を用いる。filter()の第1引数にはデータフレーム名，第2引数には抽出条件を指定する。

抽出条件は，TRUEかFALSEの論理値を返す式で指定する。TRUEのデータだけが抽出される。具体的な書き方は次の通り。

条件	意味
x == 10	xが10と等しい
x != 10	xが10と等しくない !(x==10)とも書ける
x > 10	xが10より大きい
x < 10	xが10より小さい
x >= 10	xが10以上
x <= 10	xが10以下
x %in% A	xがベクトルAに含まれる

このうち，とくにAとBが等しいという条件を (A == B)と書くことに注意。=が1つだと，条件ではなく代入を表すことになるため，エラーとなる。

iris_setosa <- filter(iris, Species == "setosa")
iris_versicolor <- filter(iris, Species == "versicolor")
iris_virginica <- filter(iris, Species == "virginica")

これでアヤメの種類ごとのデータフレームができるので，それぞれのデータフレームで要約統計表を作成する。

stargazer(
  iris_setosa, type = "text", title = "setosa", digits = 2,
  summary.stat = c("mean", "sd", "min", "p25", "median", "p75", "max")
)

setosa
=============================================================
Statistic    Mean St. Dev. Min  Pctl(25) Median Pctl(75) Max 
-------------------------------------------------------------
Sepal.Length 5.01   0.35   4.30   4.80    5.00    5.20   5.80
Sepal.Width  3.43   0.38   2.30   3.20    3.40    3.68   4.40
Petal.Length 1.46   0.17   1.00   1.40    1.50    1.58   1.90
Petal.Width  0.25   0.11   0.10   0.20    0.20    0.30   0.60
-------------------------------------------------------------

stargazer(
  iris_versicolor, type = "text", title = "versicolor", digits = 2,
  summary.stat = c("mean", "sd", "min", "p25", "median", "p75", "max")
)

versicolor
=============================================================
Statistic    Mean St. Dev. Min  Pctl(25) Median Pctl(75) Max 
-------------------------------------------------------------
Sepal.Length 5.94   0.52   4.90   5.60    5.90    6.30   7.00
Sepal.Width  2.77   0.31   2.00   2.52    2.80    3.00   3.40
Petal.Length 4.26   0.47   3.00   4.00    4.35    4.60   5.10
Petal.Width  1.33   0.20   1.00   1.20    1.30    1.50   1.80
-------------------------------------------------------------

stargazer(
  iris_virginica, type = "text", title = "virginica", digits = 2,
  summary.stat = c("mean", "sd", "min", "p25", "median", "p75", "max")
)

virginica
=============================================================
Statistic    Mean St. Dev. Min  Pctl(25) Median Pctl(75) Max 
-------------------------------------------------------------
Sepal.Length 6.59   0.64   4.90   6.23    6.50    6.90   7.90
Sepal.Width  2.97   0.32   2.20   2.80    3.00    3.18   3.80
Petal.Length 5.55   0.55   4.50   5.10    5.55    5.88   6.90
Petal.Width  2.03   0.27   1.40   1.80    2.00    2.30   2.50
-------------------------------------------------------------