第2講 - データの扱い・可視化・確率シミュレーション
複数の個体について,いくつかの属性を集計した表
実データの最も一般的な形式
データフレームの例
ある小学校の1年生の身長・体重・性別・血液型のデータ
| 名前 | 身長 [cm] | 体重 [kg] | 性別 | 血液型 |
|---|---|---|---|---|
| 太郎 | 108 | 19 | 男 | B |
| 花子 | 116 | 21 | 女 | O |
| 次郎 | 130 | 25 | 男 | AB |
| … | … | … | … | … |
データ操作とグラフィクスの枠組
本講義では tidyverse を中心に説明
パッケージ集の利用には以下が必要
作成方法はいくつか用意されている
tibble::tibble())dplyr::bind_cols())データフレームの作成の例
関数 readr::write_csv() : ファイルの書き出し
関数 readr::read_csv() : ファイルの読み込み
base::data.frame() とは行名の扱いに違いがあるので注意選択方法はいくつか用意されている
データフレームの要素の選択
z <- tibble(one = c(1,2,3),
two = c("AB","CD","EF"),
three = 6:8)
z[1,2] # 1行2列の要素を選択
z[-c(1,3),] # 1,3行を除外
z[c(TRUE,FALSE,TRUE),] # 1,3行を選択
z[,"two"] # 列名"two"を選択(1列のデータフレームになる)
z["two"] # 上記と同様の結果
z[,c("one","three")] # 列名"one"と"three"を選択(データフレームになる)
z[c("one","three")] # 上記と同様の結果
z[["two"]] # 列名"two"のベクトルを選択(1列の場合しか使えない)
z$two # 上記と同様の結果関数 dplyr::filter() : 条件を指定して行を選択
行に関する条件指定には以下を用いることができる
== (否定は !=)<,>,<=,>=& (かつ), | (または)関数 dplyr::select() : 条件を指定して列を選択
条件指定には例えば以下のような方法がある
!列名, !c(列名,列名,...), !(列名:列名)starts_with("文字列")ends_with("文字列"),& (かつ), | (または)処理を順次結合する演算子 (いくつか定義がある)
|> (base R で定義; この講義ではこちらで記述する)%>% (package::magrittr)データフレームの部分集合の取得
関数 dplyr::pivot_longer() : 同じ属性の列をまとめる
pivot_longer(
data,
cols,
...,
cols_vary = "fastest",
names_to = "name", names_prefix = NULL, names_sep = NULL, names_pattern = NULL,
names_ptypes = NULL, names_transform = NULL, names_repair = "check_unique",
values_to = "value", values_drop_na = FALSE, values_ptypes = NULL,
values_transform = NULL
)
#' data: データフレーム
#' cols: 操作の対象とする列(列の番号,名前,名前に関する条件式など)
#' names_to: 対象の列名をラベルとする新しい列の名前(既定値は"name")
#' values_to: 対象の列の値を保存する新しい列の名前(既定値は"value")
#' 詳細は '?dplyr::pivot_longer' を参照列ごとのグラフを視覚化する際に多用する
成績表の形式の変更
base::sum() : 総和を計算するbase::mean() : 平均base::max() : 最大値base::min() : 最小値stats::median() : 中央値stats::quantile() : 分位点関数 dplyr::summarise() : 集計を行う
平均の算出
関数 dplyr::group_by() : グループ化を行う
グループごとに集計
package::graphics (標準で読み込まれる)package::ggplot2package::ggplot2 ではさまざまな作図関数を追加しながら描画する
関数 ggplot2::ggplot() : 初期化
関数 ggplot2::geom_line() : 線の描画
geom_line(
mapping = NULL,
data = NULL,
stat = "identity",
position = "identity",
...,
na.rm = FALSE, orientation = NA, show.legend = NA, inherit.aes = TRUE
)
#' mapping: "審美的"マップの設定
#' data: データフレーム
#' stat: 統計的な処理の指定
#' position: 描画位置の調整
#' ...: その他の描画オプション
#' na.rm: NA(欠損値)の削除(既定値は削除しない)
#' show.legend: 凡例の表示(既定値は表示)
#' 詳細は '?ggplot2::geom_line' を参照pcr_data |> select(!c(sub,total)) |>
pivot_longer(!date, names_to = "organ", values_to = "nums") |>
ggplot(aes(x = date, y = nums, colour = organ)) +
labs(title = "PCR Tests in Various Organizatios", x = "Date", y = "Number of Tests") +
geom_line(show.legend = FALSE) + # 凡例を消す
facet_grid(vars(organ)) # "organ" ごとに異なる図を並べる関数 ggplot2::geom_point() : 点の描画
geom_point(
mapping = NULL,
data = NULL,
stat = "identity",
position = "identity",
...,
na.rm = FALSE, show.legend = NA, inherit.aes = TRUE
)
#' mapping: 審美的マップの設定
#' data: データフレーム
#' stat: 統計的な処理の指定
#' position: 描画位置の調整
#' ...: その他の描画オプション
#' na.rm: NA(欠損値)の削除(既定値は削除しない)
#' show.legend: 凡例の表示(既定値は表示)
#' 詳細は '?ggplot2::geom_point' を参照if(Sys.info()["sysname"] == "Darwin") { # MacOSか調べて日本語フォントを指定
theme_update(text = element_text(family = "HiraginoSans-W4"))}
pcr_data |>
ggplot(aes(x = niid, y = mi)) + # x軸を niid,y軸を mi に設定
geom_point(colour = "blue", shape = 19) + # 色と形を指定(点の形は '?points' を参照)
labs(x = pcr_colnames["niid"], y = pcr_colnames["mi"]) # 軸の名前を指定散布図行列は複数の散布図を行列状に配置したもの
関数 GGally::ggpairs() : 散布図行列の描画
#' 必要であれば 'install.packages("GGally")' を実行
library(GGally) # パッケージのロード
ggpairs(
data, mapping = NULL,
columns = 1:ncol(data),
upper = list(continuous = "cor", # 連続量の扱い
combo = "box_no_facet", # 連続・離散の扱い
discrete = "count", # 離散量の扱い
na = "na"), # 欠損の扱い
lower = list(continuous = "points",
combo = "facethist",
discrete = "facetbar",
na = "na"),
diag = list(continuous = "densityDiag",
discrete = "barDiag",
na = "naDiag"),
...,
axisLabels = c("show", "internal", "none"),
columnLabels = colnames(data[columns]),
legend = NULL
)
#' columns: 表示するデータフレームの列を指定
#' upper/lower/diag: 行列の上三角・下三角・対角の表示内容を設定
#' axisLabels: 各グラフの軸名の扱い方を指定
#' columnLabels: 表示する列のラベルを設定(既定値はデータフレームの列名)
#' legend: 凡例の設定(どの成分を使うか指定)
#' 詳細は '?GGally::ggpairs' を参照pcr_data |> select(!c(sub,total)) |> # 日付から四半期の因子を作成
mutate(quarter = as_factor(quarter(date, with_year = TRUE))) |>
ggpairs(columns = 2:8, columnLabels = pcr_colnames[-c(1,8,10)], axisLabels = "none",
aes(colour = quarter), legend = c(2,1), # 四半期ごとに色づけて(1,1)の凡例を使用
upper = "blank", diag = list(continuous = "barDiag")) +
theme(legend.position = "top") # 凡例を上に表示RStudioの機能を使う (少数の場合はこちらが簡便)
関数 ggsave() : 図の保存
データの値の範囲をいくつかの区間に分割し, 各区間に含まれるデータの個数を棒グラフにした図
データ散らばり具合を考察するための図
geom_boxplot(
mapping = NULL, data = NULL, stat = "boxplot", position = "dodge2",
...,
outlier.colour = NULL, outlier.color = NULL, outlier.fill = NULL,
outlier.shape = 19, outlier.size = 1.5, outlier.stroke = 0.5, outlier.alpha = NULL,
notch = FALSE, notchwidth = 0.5, varwidth = FALSE,
na.rm = FALSE, orientation = NA, show.legend = NA, inherit.aes = TRUE
)
#' ourlier.*: 外れ値の描画方法の指定
#' notch*: ボックスの切れ込みの設定
#' varwidth: ボックスの幅でデータ数を表示
#' 詳細は '?ggplot2::geom_boxplot' を参照項目ごとの量を並べて表示した図
#' 機関ごとの月の検査件数の推移 (2021年分)
pcr_data |>
filter(year(date) == 2021) |>
mutate(month = as_factor(month(date))) |> # 月を作成
select(!c(date,sub,total)) |> # 機関に限定
group_by(month) |> # 月でグループ化
summarize(across(everything(), sum)) |> # 全て(月以外)を集計
pivot_longer(!month, names_to = "organ", values_to = "nums",
names_transform = list(organ = as_factor)) |>
## 最後のオプションは organ 列のラベルを出てきた順で因子化して元の列の並びにしている
ggplot(aes(x = organ, y = nums, fill = month)) +
geom_bar(stat = "identity", position = "dodge", na.rm = TRUE) +
theme(legend.position = "top") + guides(fill = guide_legend(nrow = 1))?Random 参照)set.seed())base::sample() : ランダムサンプリングstats::rbinom() : 二項乱数stats::runif() : 一様乱数stats::rnorm() : 正規乱数定理
\(X_1,X_2,\dotsc\) を独立同分布な確率変数列とし, その平均を \(\mu\) ,標準偏差を \(\sigma\) とする. このとき,すべての実数 \(a< b\) に対して
\[\begin{equation} P\Bigl(a\leq\frac{\sqrt{n}(\bar{X}_n-\mu)}{\sigma}\leq b \Bigr) \to\frac{1}{\sqrt{2\pi}}\int_a^be^{-\frac{x^2}{2}}dx\quad (n\to\infty) \end{equation}\]
が成り立つ.
#' 確率変数の分布の設定 (例 : 区間[-1,1]の一様乱数)
mc_rand <- function(n) { # n個の乱数を生成
return(runif(n, min = -1, max = 1))
}
#' 標本平均の計算
mc_mean <- function(n) { # n個のデータで計算
return(mean(mc_rand(n)))
}
#' Monte-Carlo実験
set.seed(123) # 実験を再現したい場合はシード値を指定する
mu <- 0; sigma <- sqrt(1/3) # 理論平均と標準偏差
mc_num <- 5000 # 実験の繰り返し回数
for(n in c(1,2,4,8,16)){ # nを変えて実験
p <- tibble(x = replicate(mc_num, mc_mean(n))) |> # 繰り返し実験し標本平均を記録
ggplot(aes(x)) +
geom_histogram(aes(y = after_stat(density)), # 密度表示
fill = "orchid", alpha = 0.5, # 塗り潰しの色
colour = "purple") + # 境界線の色
geom_function(fun = \(x) dnorm(x, mean = mu, sd = sigma/sqrt(n)),
colour = "orange", linewidth = 1.5) + # 理論曲線を重ねる
labs(x = expression(bar(X)), # x軸の表示
title = paste0("n=", n)) # タイトルにnを記載
print(p) # for 文の中では明示的に print する必要がある
}AとBの二人で交互にコインを投げる. 最初に表が出た方を勝ちとするとき, AとBそれぞれの勝率はいくつとなるか?
コイン投げは関数 sample(), rbinom() などを用いて模擬できる
#' コイン投げの試行 (いろいろな書き方があるので以下は一例)
mc_trial <- function(){
while(TRUE){ # 永久に回るループ
if(rbinom(1, size = 1, prob = 0.5)==1){return("A")} # Aが表で終了
if(rbinom(1, size = 1, prob = 0.5)==1){return("B")} # Bが表で終了
#' どちらも裏ならもう一度ループ
}
}
#' Monte-Carlo実験
set.seed(8888) # 実験を再現したい場合はシード値を指定する
mc_num <- 10000 # 実験回数を設定
mc_data <- replicate(mc_num, mc_trial())
#' 簡単な集計
table(mc_data) # 頻度
table(mc_data)/mc_num # 確率(推定値)R言語の基礎 第5章 (pp71-82)
などの実装例がある