からっぽのしょこ

読んだら書く!書いたら読む!同じ事は二度調べ(たく)ない

【R】第1講:度数分布表とヒストグラムの作成【統計学入門(小島)のノート】

はじめに

 この記事は「統計学 Advent Calendar 2022」の13日目の記事です。

 『完全独習 統計学入門』の学習ノートです。本に載っている計算や表、グラフをR言語で再現します。本とあわせて読んでください。
 この記事では、度数分布表とヒストグラムを作成します。

【他の節の内容】

www.anarchive-beta.com

【この節の内容】

第1項 度数分布表とヒストグラムで、データの特徴を浮き彫りにする

 データセットから度数分布表とヒストグラムを作成します。

 利用するパッケージを読み込みます。

# 利用パッケージ
library(tidyverse)
library(magick)

 この記事では、基本的にパッケージ名::関数名()の記法を使うので、パッケージを読み込む必要はありません。ただし、作図コードがごちゃごちゃしないようにパッケージ名を省略しているためggplot2を読み込む必要があります。
 また、ネイティブパイプ演算子|>を使っています。magrittrパッケージのパイプ演算子%>%に置き換えても処理できますが、その場合はmagrittrも読み込む必要があります。

1-1 生データでは何もわからないから統計を使う

 まずは、データを読み込みます。

 図表1-1のデータセットをベクトルとして作成します。

# データを作成:(身長)
data_x <- c(
  151, 154, 158, 162, 
  154, 152, 151, 167, 
  160, 161, 155, 159, 
  160, 160, 155, 153, 
  163, 160, 165, 146, 
  156, 153, 165, 156, 
  158, 155, 154, 160, 
  156, 163, 148, 151, 
  154, 160, 169, 151, 
  160, 159, 158, 157, 
  154, 164, 146, 151, 
  162, 158, 166, 156, 
  156, 150, 161, 166, 
  162, 155, 143, 159, 
  157, 157, 156, 157, 
  162, 161, 156, 156, 
  162, 168, 149, 159, 
  169, 162, 162, 156, 
  150, 153, 159, 156, 
  162, 154, 164, 161
)
length(data_x)
## [1] 80

 データ数は80です。

1-2 ヒストグラムを作る

 次に、度数分布表とヒストグラムを作成します。

 データの最小値と最大値を抽出します。

# 最小値・最大値を抽出
min_x <- min(data_x)
max_x <- max(data_x)
min_x; max_x
## [1] 143
## [1] 169

 min()で最小値、max()で最大値を得られます。

 階級の上限と下限として用いる値を作成します。

# 階級の幅を指定
class_width <- 5

# 階級用の値を作成
class_vec <- seq(
  from = floor(min_x/10) * 10, 
  to = ceiling(max_x/10) * 10, 
  by = class_width
)
class_vec
## [1] 140 145 150 155 160 165 170

 階級の幅をclass_widthとして値を指定します。
 seq()by引数にclass_widthを指定して、class_width間隔の数値ベクトルを作成します。
 最小値の1桁をfloor()で切り捨てて、最大値の1桁をceiling()で切り上げて使います(もっといいやり方があれば教えて下さい)。あるいは、from, to(第1・第2)引数に値を直接指定します。

 図表1-2(離散値の場合?)の度数分布表を作成します。

# 度数分布表を作成:(図表1-2)
freq_df <- tibble::tibble(
  class_lower = class_vec[-length(class_vec)] + 1, # 階級の下限
  class_upper = class_vec[-1] # 階級の上限
) |> 
  dplyr::group_by(class_lower, class_upper) |> # 階級ごとの計算用
  dplyr::mutate(
    class_value = median(c(class_lower, class_upper)), # 階級値
    freq = sum((data_x >= class_lower) == (data_x <= class_upper)) # 度数
  ) |> 
  dplyr::ungroup() |> 
  dplyr::mutate(
    relative_freq = freq / length(data_x), # 相対度数
    cumulative_freq = cumsum(freq) # 累積度数
  )
freq_df
## # A tibble: 6 × 6
##   class_lower class_upper class_value  freq relative_freq cumulative_freq
##         <dbl>       <dbl>       <dbl> <int>         <dbl>           <int>
## 1         141         145         143     1        0.0125               1
## 2         146         150         148     6        0.075                7
## 3         151         155         153    19        0.238               26
## 4         156         160         158    30        0.375               56
## 5         161         165         163    18        0.225               74
## 6         166         170         168     6        0.075               80

 階級の下限をclass_lower列、上限をclass_upper列として、値を設定に合わせてデータフレームに格納します。
 group_by()で階級ごとにグループ化して、階級値と度数を計算します。
 階級値は、行(階級)ごとのclass_lower, class_upperの中央値をmedian()を計算します。
 度数は、class_lower以上(data_x >= class_lower)でclass_upper以下(data_x <= class_upper)であるdata_xの要素数をカウントします。2つの条件がどちらもTRUEとなる数をsum()を使ってカウントします。(もっといい方法があると思います。)
 相対度数は、度数をデータ数で割ります。
 累積度数は、度数の累積和をcumsum()で計算します。

 ただしこの方法では、連続値の場合に対応できません。
 連続値の場合の度数・相対度数・累積度数を計算します。

# 度数分布表を作成:(階級が決まっている場合)
freq_df <- tibble::tibble(
  class_lower = class_vec[-length(class_vec)], # 階級の下限
  class_upper = class_vec[-1] # 階級の上限
) |> 
  dplyr::group_by(class_lower, class_upper) |> 
  dplyr::mutate(
    class_value = median(c(class_lower, class_upper)), # 階級値
    freq = sum((data_x >= class_lower) == (data_x < class_upper)) # 度数
  ) |> 
  dplyr::ungroup() |> 
  dplyr::mutate(
    relative_freq = freq / length(data_x), # 相対度数
    cumulative_freq = cumsum(freq) # 累積度数
  )
freq_df
## # A tibble: 6 × 6
##   class_lower class_upper class_value  freq relative_freq cumulative_freq
##         <dbl>       <dbl>       <dbl> <int>         <dbl>           <int>
## 1         140         145        142.     1        0.0125               1
## 2         145         150        148.     4        0.05                 5
## 3         150         155        152.    17        0.212               22
## 4         155         160        158.    27        0.338               49
## 5         160         165        162.    23        0.288               72
## 6         165         170        168.     8        0.1                 80

 class_lower列以上でclass_upper列以下(data_x < class_upper)のデータをカウントします。

 階級の上限と下限を指定して、度数分布表を作成しました。ただし、階級の範囲を図1-2の設定と変更したため、階級値が変わってしまいました。
 階級値(階級の中央値)を指定して、度数分布表を作成します。その前に、階級値を作成します。

# 階級値を作成
class_value_vec <- seq(from = 143, to = ceiling(max_x/10)*10, by = class_width)
class_value_vec
## [1] 143 148 153 158 163 168

 データから階級値を作成します。

 度数・相対度数・累積度数を計算します。

# 度数分布表を作成:(階級値が決まっている場合)
freq_df <- tibble::tibble(
  class_lower = class_value_vec - 0.5*class_width, # 階級の下限
  class_upper = class_value_vec + 0.5*class_width, # 階級の上限
  class_value = class_value_vec # 階級値
) |> 
  dplyr::group_by(class_lower, class_upper) |> 
  dplyr::mutate(
    freq = sum((data_x >= class_lower) == (data_x < class_upper)) # 度数
  ) |> 
  dplyr::ungroup() |> 
  dplyr::mutate(
    relative_freq = freq / length(data_x), # 相対度数
    cumulative_freq = cumsum(freq) # 累積度数
  )
freq_df
## # A tibble: 6 × 6
##   class_lower class_upper class_value  freq relative_freq cumulative_freq
##         <dbl>       <dbl>       <dbl> <int>         <dbl>           <int>
## 1        140.        146.         143     1        0.0125               1
## 2        146.        150.         148     6        0.075                7
## 3        150.        156.         153    19        0.238               26
## 4        156.        160.         158    30        0.375               56
## 5        160.        166.         163    18        0.225               74
## 6        166.        170.         168     6        0.075               80

 階級値から「階級の幅の半分」を引いた値を下限、足した値を上限として、階級値と共に格納して計算します。

 図1-3の(に近い)グラフを作成するには、階級値が一致している必要があるので、このデータフレーム(度数分布表)を使います。
 ヒストグラムを作成します。

# ヒストグラムを作成
ggplot() + 
  geom_bar(data = freq_df, mapping = aes(x = class_value, y = freq), 
           stat = "identity", width = class_width, fill = "#00A968") + 
  scale_x_continuous(breaks = freq_df[["class_value"]]) + # x軸目盛
  labs(title = "Histogram", 
       subtitle = paste0("N=", length(data_x)), 
       x = "class value", y = "frequency")

geom_bar関数によるヒストグラム

 横軸を階級値(class_value列)、縦軸を度数(freq列)とする棒グラフをgeom_bar()で描画します。
 横軸の体裁をscale_x_continuous()で設定できます。

# 軸目盛ラベルを作成
x_label_vec <- paste0(freq_df[["class_lower"]], "~", freq_df[["class_upper"]])

# ヒストグラムを作成
ggplot() + 
  geom_bar(data = freq_df, mapping = aes(x = class_value, y = freq), 
           stat = "identity", width = class_width, fill = "#00A968") + 
  scale_x_continuous(breaks = freq_df[["class_value"]], 
                     labels = x_label_vec) + # x軸目盛
  labs(title = "Histogram", 
       subtitle = paste0("N=", length(data_x)), 
       x = "class value", y = "frequency")

ヒストグラム


 度数分布表が不要であれば、もう少し簡単にヒストグラムを作成できます。
 データをそのままデータフレームに格納します。

# データを格納
data_df <- tibble::tibble(
  x = data_x
)
data_df
## # A tibble: 80 × 1
##        x
##    <dbl>
##  1   151
##  2   154
##  3   158
##  4   162
##  5   154
##  6   152
##  7   151
##  8   167
##  9   160
## 10   161
## # … with 70 more rows


 階級値を使ってヒストグラムを作成します。

# 階級値を作成
class_value_vec <- seq(from = 143, to = 175, by = class_width)

# ヒストグラムを作成:(階級値が決まっている場合)
ggplot() + 
  geom_histogram(data = data_df, mapping = aes(x = x, y = ..count..), 
                 breaks = class_value_vec-0.5*class_width, fill = "#00A968") + 
  scale_x_continuous(breaks = class_value_vec) + # x軸目盛
  labs(title = "Histogram", 
       subtitle = paste0("N=", length(data_x)), 
       x = "class value", y = "frequency")

geom_histogram関数によるヒストグラム

 geom_histogram()でヒストグラムを描画できます。breaks引数に階級の下限を指定します。

 同様に、階級を使ってヒストグラムを作成します。

# 階級値を計算
class_value_vec <- cbind(class_vec[-1], class_vec[-length(class_vec)]) |> 
  apply(1, median)

# ヒストグラムを作成:(階級が決まっている場合)
ggplot() + 
  geom_histogram(data = data_df, mapping = aes(x = x, y = ..count..), 
                 breaks = class_vec, fill = "#00A968") + # ヒストグラム
  scale_x_continuous(breaks = class_value_vec) + # x軸
  labs(title = "Histogram", 
       subtitle = paste0("N=", length(data_x)), 
       x = "class value", y = "frequency")

geom_histogram関数によるヒストグラム

 class_lower, class_upper列に対応する2列のマトリクスをcbind()で作成し、apply()で行ごとに中央値(階級値)を計算して、軸目盛ラベルとして使っています。

・階級とヒストグラムの関係

 本では、最小値から最大値の範囲を5~8の小範囲(小区間)に区切るとあります。階級(階級値)とヒストグラムの関係を見てみましょう。

作図コード(クリックで展開)

 階級の幅を1ずつ大きくするヒストグラムのアニメーションを作成します。

# 画像の保存先を指定
dir_path <- "tmp_folder"

# 階級幅の最大値(フレーム数)を指定
class_width_max <- 30

# 階級を変更して作図
for(class_width in 1:class_width_max) {
  
  # 階級用の値を作成
  class_vec <- seq(from = floor(min_x/10)*10, to = ceiling(max_x/10)*10, by = class_width)
  class_vec <- c(class_vec, class_vec[length(class_vec)]+class_width) # 集計用に上限以上の値を1つ追加
  
  # 階級値を計算
  class_value_vec <- cbind(class_vec[-length(class_vec)], class_vec[-1]) |> 
    apply(1, median)
  
  # 軸目盛ラベルを作成
  x_label_vec <- paste0(
    class_vec[-length(class_vec)], "~", class_vec[-1], 
    "\n(", class_value_vec, ")"
  )
  
  # ヒストグラムを作成
  g <- ggplot() + 
    geom_histogram(data = data_df, mapping = aes(x = x), 
                   breaks = class_vec, fill = "#00A968") + 
    scale_x_continuous(breaks = class_value_vec, labels = x_label_vec) + # x軸目盛
    coord_cartesian(xlim = c(min(data_x), max(data_x)), ylim = c(0, length(data_x))) + # 表示範囲
    labs(title = "Histogram", 
         subtitle = paste0("N=", length(data_x), ", size:", class_width), 
         x = "class value", y = "frequency")
  
  # グラフを書き出し
  file_name <- sprintf(paste0("%0", nchar(class_width_max), "d"), class_width)
  ggplot2::ggsave(
    filename = paste0(dir_path, "/", file_name, ".png"), plot = g, 
    width = 800, height = 600, units = "px", dpi = 100
  )
}

# ファイル名を取得
file_name_vec <- list.files(dir_path)

# gif画像を作成
paste0(dir_path, "/", file_name_vec) |> # ファイルパスを作成
  magick::image_read() |> # 画像ファイルを読み込み
  magick::image_animate(fps = 1, dispose = "previous") |> # gif画像を作成
  magick::image_write_gif(path = "histgram.gif", delay = 1/2) -> tmp_path # gifファイルを書き出し

 階級の幅の最大値をclass_width_maxとして値を指定します。for()を使って、1からclass_width_maxの整数を順番にclass_vecとします。
 同様に作図してそれぞれグラフを保存します。画像の書き出し先フォルダdir_pathは空である必要があります。また、画像ファイルの読み込み時に、書き出した順番と(文字列基準で)同じ並びになるようなファイル名を付ける必要があります。
 保存した画像ファイルをimage_read()で読み込んで、image_animate()image_write_gif()でgif画像に変換します。


階級とヒストグラムの関係

 階級の取り方によってヒストグラムの印象が変わるのが分かります。

練習問題

 図表1-4のデータセットを作成します。

# データを作成:(体重)
data_x <- c(
  48, 54, 47, 50, 53, 43, 45, 43, 
  44, 47, 58, 46, 46, 63, 49, 50, 
  48, 43, 46, 45, 50, 53, 51, 58, 
  52, 53, 47, 49, 45, 42, 51, 49, 
  58, 54, 45, 53, 50, 69, 44, 50, 
  58, 64, 40, 57, 51, 69, 58, 47, 
  62, 47, 40, 60, 48, 47, 53, 47, 
  52, 61, 55, 55, 48, 48, 46, 52, 
  45, 38, 62, 47, 55, 50, 46, 47, 
  55, 48, 50, 50, 54, 55, 48, 50
)
length(data_x)
## [1] 80


 度数分布表とヒストグラムを作成します。階級値class_value_vecの作成時の下限(from引数)の設定を変更すれば、他は同じコードで処理できます。

## # A tibble: 7 × 6
##   class_lower class_upper class_value  freq relative_freq cumulative_freq
##         <dbl>       <dbl>       <dbl> <int>         <dbl>           <int>
## 1        35.5        40.5          38     3        0.0375               3
## 2        40.5        45.5          43    11        0.138               14
## 3        45.5        50.5          48    33        0.412               47
## 4        50.5        55.5          53    19        0.238               66
## 5        55.5        60.5          58     7        0.0875              73
## 6        60.5        65.5          63     5        0.0625              78
## 7        65.5        70.5          68     2        0.025               80


 この講では、度数分布表とヒストグラムを作成しました。次講では、平均値を計算します。

参考書籍

  • 小島寛之『完全独習 統計学入門』ダイヤモンド社,2006年.

おわりに

 ずっと放置していた本なのですが、読んでみると易しく噛み砕いた内容だったんですね。

 半年だか1年くらい前に、Rまたはdplyrで度数分布表を作る方法が分からなかったのですが、こんな感じで一応できました。なのでこの記事の私的書いておきたかったことが度数の集計処理です。
 でもtidyverseの関数でもっと簡単に作れたりしませんか?上限と下限の両方を満たすTRUE==で探してるのが何とも言えない気持ちになります。

  • 2023.01.19

 geom_histogram()の使い方をミスってたので書き直しました。

【次の節の内容】

www.anarchive-beta.com