データを入手したら，分析を行う前にまずデータの概観を把握することが重要。データセットにどのような変数が含まれているかを確認し，それぞれの変数の記述統計量 (平均値，分位点，分散など)を求め，必要であれば変数間の相関関係を調べてみる。これにより，データに異常値が含まれていたり，分析するには欠損値が多かったりといったことがわかる可能性がある。
ここでは，Rのパッケージ，tidyverseとstargazerを利用するのであらかじめインストールしておく。
まず，インストールしたパッケージを読み込む。
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.2 ✔ tibble 3.2.1
✔ lubridate 1.9.4 ✔ tidyr 1.3.1
✔ purrr 1.0.4
── 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
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を読み込む。
irisというデータフレームが読み込まれる。irisデータには，植物のアヤメ150個体について，種類 (3分類)，萼片 (がくへん)の長さと幅，花弁 (かべん)の長さと幅の5つの変数が記録されている。このデータは，機械学習による分類などの練習用データとしてよく用いられているデータである。
Sepal.Length萼片の長さ
Sepal.Width萼片の幅
Petal.Length花弁の長さ
Petal.Width花弁の幅
Speciesアヤメの種類
データにどのような変数が含まれているかは，小規模データであればEnvironmentタブからスプレッド・シートを表示させれば把握できる。大規模なデータは，スプレッド・シートを表示させるだけでも時間がかかることがあるので，glimpse()関数を使うのが良い。
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()関数を用いれば，データフレームに含まれるすべての変数の要約統計量を表示させることができる。
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列目までの変数を取り出したデータフレームとなる。
この処理は，パイプを使って以下のように書くこともできる。
|> 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を指定して相関行列を求める。
<- 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
同じことをパイプを使って書くと次のようになる。
|> 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種類に分類されている。そこで，アヤメの種類ごとに萼片と花弁の長さ，幅の平均を求めてみよう。そのための方法はいくつか考えられる。
フィルター
まず，データフレームからアヤメの種類ごとにデータを取り出して分割する方法。単純でわかりやすいが，分類が多くなると大変なので，要約統計量を求めるためだけであれば効率的な方法ではない。ただし，データフレームにフィルターをかけて特定の条件を満たす観測 (行)を取り出すことはよくあるので，やり方を覚えておこう。
データフレームから条件を満たす行を取り出すには，filter()関数を用いる。filter()の第1引数にはデータフレーム名，第2引数には抽出条件を指定する。
抽出条件は，TRUEかFALSEの論理値を返す式で指定する。TRUEのデータだけが抽出される。具体的な書き方は次の通り。
x == 10xが10と等しい
x != 10xが10と等しくない !(x==10)とも書ける
x > 10xが10より大きい
x < 10xが10より小さい
x >= 10xが10以上
x <= 10xが10以下
x %in% AxがベクトルAに含まれる
このうち，とくにAとBが等しいという条件を (A == B)と書くことに注意。=が1つだと，条件ではなく代入を表すことになるため，エラーとなる。
<- filter (iris, Species == "setosa" )<- filter (iris, Species == "versicolor" )<- filter (iris, Species == "virginica" )
これでアヤメの種類ごとのデータフレームができるので，それぞれのデータフレームで要約統計表を作成する。
stargazer (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 (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 (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
-------------------------------------------------------------
group_by()を使う
たぶんこれが標準的な方法。
group_by()関数を使うと，(Rの内部で)グループに分類されたデータフレームを作成することができる。第1引数には分類したいデータフレーム名，第2引数以降は分類を行うためのファクター型変数を指定する。複数のファクター型変数を指定して，階層的な分類を行うことも可能。グループに分類されたデータフレームにsummarize()関数を適用すると，グループごとの要約統計量が計算される。
<- group_by (iris, 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
同じことを，パイプを使うと以下のように書ける。
|> 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ごとのデータ数の表を作成できる。
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