はじめに
この記事は「統計学 Advent Calendar 2022」の13日目の記事です。
『完全独習 統計学入門』の学習ノートです。本に載っている計算や表、グラフをR言語で再現します。本とあわせて読んでください。
この記事では、度数分布表とヒストグラムを作成します。
【他の節の内容】
【この節の内容】
第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")
横軸を階級値(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()
でヒストグラムを描画できます。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")
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()
の使い方をミスってたので書き直しました。
【次の節の内容】