library(tidyverse)
# tidyverse系命名規則
# 1. 関数・オブジェクト・列名: snake_case
# 例: read_data(), daily_sales, customer_id, is_active
# 2. S3メソッド: generic.class
# 例: print.my_model, predict.my_model
# 3. クラス名（S4/R6）: CamelCase
# 例: LinearModel, DataLoader
# 4. 定数（ほぼ使わない）: ALL_CAPS
#
# コツ
# 1. ブールは接頭辞：is_, has_, can_, should_（例: is_valid, has_missing）
# 2. 単位や型は接尾辞：_pct, _ms, _id, _dt, _usd
# 3. 関数名は動詞から：load_data(), fit_model(), plot_residuals()
# 4. 略語は最小限に：num_customers はOK、ncst はNG
tidyverse_frequently_used <-
iris %>%
# filter()関数。引数を識別しやすいように改行。
dplyr::filter(
Species == "setosa"
) %>%
# 短いコードなら改行しなくても構わない。
dplyr::select(Species, Sepal.Length, Petal.Length) %>%
# ライブラリ名は宣言したほうがトラブルが少ない。
# たとえライブラリをロードした後としても。
dplyr::mutate(ratio = Sepal.Length / Petal.Length) %>%
dplyr::group_by(Species) %>%
dplyr::summarise(
mean_ratio = mean(ratio, na.rm = TRUE)
)Rいろいろ
R
Rで書いたいろいろなコードから、メモを残したほうがよさげなものをよりぬき。
1 コードを書く
2 数多くシートを含むMSExcelファイルをたっぷり読み込んで結合してTable 1をつくる
賃金構造基本統計調査（都道府県パネル）について、データを読み込んでTable 1と相関行列とをつくって、多重共線性評価までやるRコード。重みは従業者数。
再現性に関する方針は、以下に示す通り。
+ 列名は `janitor::clean_names()` 後の **snake_case**を用いる（例：`N_of_employees` → `n_of_employees`）。
+ 欠損・ゼロは明示的に扱う（例：`filter(n_of_employees > 0)`）。
+ .featherにて保存・読込。.rdsや.csvなど、これまで用いたファイル形式に比べて作業速いから。Rコードをすべて列記すると、読むだけで疲れてしまう。小職がひっかかったところだけ、抜粋して解説する。完全なコードは、Githubにおいてある。それぞれ、 データ作成とTable 1など記述的処理をみてね。
2.1 データを作る
読み取るべき変数を指定する。もともとMSExcel（.csvだけど）ファイルで、変数は横に並んでいる。この並んでいる各変数を書き忘れるから、原本にあたって漏れなく書き出す。
items <-
c(
"年齢",
"勤続年数",
"所定内実労働時間数",
"超過実労働時間数",
"きまって支給する現金給与額",
"所定内給与額",
"年間賞与その他特別給与額",
"労働者数"
)stringr()を用いて、必要がない記号を削除する。
本Rコードは、年次をMSExcelファイルを収納するフォルダ名から、都道府県名と産業名とをMSExcelシート名から頂戴する。ファイル名は、バッチ処理するときに読み取りファイルリストとしてのみ使う。
いま、ファイル名は(1-8-1-1)2010pref_01-02.xlsxで、また同ファイルに含まれるシート群には01-04北海道Ｅ製造業があるとしよう。シート名については、「北海道」と「製造業」だけが必要な情報である。言い換えれば、シート名を構成する英数字と記号は不要で、かつ都道府県名と産業名とを切り分けられるならば、自動的に同情報を引っ張り出せる。
というわけで、不要な文字列を削除する関数。stringr::str_replace_all()がいい仕事をする。以下はいい仕事をするために必要な正規表現。
[]：角括弧内にあるいずれかの記号A-Za-z0-9：アルファベット大文字小文字AからZまでと数字。[:digit:]や[:alpha:]でもよさげ。[:punct:]：記号類。\\s：空白。
clean_sheet_name <-
function(x) {
x %>%
stringi::stri_trans_nfkc() %>%
stringr::str_replace_all(
"[A-Za-z0-9[:punct:]\\s]+",
""
) %>%
stringr::str_squish()
}シートを読む関数。
元データがもつ性質を利用する。すなわち、
- 読み込む範囲を指定する：賃金センサスは、対象全ファイルにおいて、必要な数値はすべて
D12:AI50に必ず入っている。これ以外は読まなくてもよい。行名や列名は後から与えればよい。やはり同じだから。 - 形推定を省略させる：これが意外に計算資源を消費する。型は後からどうせ指定し直すから、さしあたりはテキストとして読み込む。
- エラー避けをいくつか仕込む。
# A function reading a sheet
read_one_sheet <- function(file, sheet) {
year <- extract_year(file)
# Fix reading ranges of target cell
rng <- "D12:AI50"
# Read data
# If failed, the process is skipped.
# (No skips were there)
raw <-
tryCatch(
readxl::read_excel(
file,
sheet = sheet,
range = rng,
col_names = FALSE,
col_types = "text"
),
error = function(e) tibble()
)
}先ほどつくったread_onesheet関数を繰り返して、多数あるファイルに含まれる多数シートを読み込む関する。
furrr::future_シリーズを、全CPU利用して走らせるようplan()にて設定する。multisessionにてたいていは問題ない。
次に、指定したフォルダ（ディレクトリ）に含まれるMSExcelファイルを探して名前を取得する。OSが変わったり、文字コード毎にいろいろ面倒が起きた時の例外処理を考えると、3つに分けたほうがいい。まったく、OS毎になんでこんなに挙動が異なるやら。もう少しなんとかしてほしい。これがあるから、WSLはありがたいともおもえるし、計算用にLinuxマシンを1台用意したくなる。
以下は、それぞれがもっている機能。
fs::dir_ls(..., glob = "**/*.{xls,xlsx,XLS,XLSX,xlsm,XLSM}")：シェルっぽいワイルドカードとglobを使う。書きやすくて速い。主にfs::由来。fs::はC++だから。 読みやすい・書きやすい（拡張子の列挙が楽）。ただし大文字小文字はOS毎に扱いが異なることから、同じような内容を列記するハメになる。fs::dir_ls(..., regexp = "\\.(xls|xlsx|xlsm)$")：正規表現にてマッチ。1パターンで拡張子を表現できる。base::list.files(..., pattern="\\.(xls|xlsx|xlsm)$", recursive=TRUE, full.names=TRUE, ignore.case=TRUE)：やはり正規表現。ignore.case=TRUEを使えるから、大文字小文字に由来する問題を回避できる。baseパッケージだからたいてい動作する。
そしてfurrr::future_map_dfr()。これがあるおかげで速い。CPUがリストにあるファイル群から手当たり次第にシートを読みまくる。おかげで計算中はCPU利用率は100%になる。お部屋を涼しくしましょう。360ファイル36,000シートでだいたい1時間弱かかる。もうひとつ入れ子になっている方もfurrr::にする手もあるが、シートまで手当たり次第に読みまくることになり、収拾がつかなくなるたぶん。
# strategy
plan(multisession, workers = parallel::detectCores() - 1)
#
read_all <-
function(
data_dir = "data",
save_per_sheet = TRUE,
out_dir = "out_sheets",
feather_only = TRUE
)
{
data_dir <- fs::path_expand(data_dir)
if (!fs::dir_exists(data_dir)) stop("フォルダが見つかりません: ",
data_dir)
files_fs_glob <- fs::dir_ls(
data_dir, recurse = TRUE, type = "file",
glob = "**/*.{xls,xlsx,XLS,XLSX,xlsm,XLSM}"
)
files_fs_re <- fs::dir_ls(
data_dir, recurse = TRUE, type = "file",
regexp = "\\.(xls|xlsx|xlsm)$"
)
files_base <- list.files(
data_dir, pattern="\\.(xls|xlsx|xlsm)$",
recursive=TRUE, full.names=TRUE, ignore.case=TRUE
)
files <- unique(c(files_fs_glob, files_fs_re, files_base))
message("Excelファイル数: ", length(files))
if (length(files) == 0) return(tibble())
# 必殺入れ子map
furrr::future_map_dfr(
files,
function(f) {
sh <- readxl::excel_sheets(f)
# 内側は通常の map_dfr（逐次）
map_dfr(seq_along(sh), function(j) {
s <- sh[j]
d <- read_one_sheet(f, s)
if (save_per_sheet) save_one_sheet(
d,
j,
out_dir = out_dir,
feather_only = feather_only
)
d
}
)
},
.options = furrr::furrr_options(scheduling = 1)
)
}
# 2.2 Table 1
毎度おなじみTable 1。dplyr::summarise()で作って整形してもいいけれど、ここはgtsummaryライブラリにお願いしましょう。例のごとく、必要な関数群を作ってから、ぶん回す。
まずは労働者数合計。特に物珍しいところはない。
# Helper: one-row total employees by gender (within a size block)
make_block_total_employees <-
function(svy_data) {
df <- svy_data$variables %>% dplyr::as_tibble()
if (!"gender" %in% names(df)) df <- df %>% dplyr::mutate(gender =factor("(Missing)"))
gtsummary::tbl_custom_summary(
data = df,
by = gender,
include = n_of_employees,
stat_fns = everything() ~ (function(x) sum(x, na.rm = TRUE)),
statistic = everything() ~ "{stat}",
label = list(n_of_employees ~ "Employees [total]")
) %>%
gtsummary::add_overall(last = TRUE) %>%
gtsummary::modify_header(all_stat_cols() ~ "**{level}**") %>%
gtsummary::bold_labels() %>%
gtsummary::modify_fmt_fun(everything() ~ scales::comma)
}次は連続量部分。使う変数選んで、名前（ラベル）を与えて、男女毎に平均（標準偏差）計算する。ポイントは、statistic = list(gtsummary::all_continuous() ~ "{mean} ({sd})")とdigits = list(gtsummary::all_continuous() ~ 1)。前者は連続量すべてについて平均と標準偏差を計算させる。後者は有効数字桁数指定。たくさん指定したくなるが、表が横長になるだけであまり意味はない。
まぁなんというか、コード長いな。とはいえ、いままでさわった言語ではもっとも書きやすい。Perlはなんだったんだろう。あれは若いうちでないとできない。
# Helper: one-row total employees by gender (within a size block)
# Helper: continuous block (weighted means with SD)
make_block_cont <- function(svy_data) {
vars_cont <-
c(
"mean_age","length_of_service",
"actual_number_of_scheduled_hours_worked",
"actual_number_of_overtime_worked",
"contractual_cash_earnings",
"hourly_scheduled_cash_earnings",
"annual_special_cash_earnings"
) %>%
intersect(names(svy_data$variables))
labels_cont <-
list(
mean_age ~ "Mean age [years]",
length_of_service ~ "Length of service [years]",
actual_number_of_scheduled_hours_worked ~ "Scheduled hours worked [h/mo]",
actual_number_of_overtime_worked ~ "Overtime hours [h/mo]",
contractual_cash_earnings ~ "Contractual cash earnings [k JPY/mo]",
hourly_scheduled_cash_earnings ~ "Hourly scheduled cash earnings [JPY/h]",
annual_special_cash_earnings ~ "Annual special cash earnings [k JPY/yr]"
)
gtsummary::tbl_svysummary(
data = svy_data,
by = gender,
include = dplyr::all_of(vars_cont),
statistic = list(gtsummary::all_continuous() ~ "{mean} ({sd})"),
digits = list(gtsummary::all_continuous() ~ 1),
missing = "ifany",
label = labels_cont
) %>%
gtsummary::add_overall(last = TRUE) %>%
gtsummary::bold_labels()
}いよいよTable 1生成。gtsummary::tbl_stack()で上から積み上げられるあたりがすばらしい。gtsummary::デフォルトだけでこれをやろうとすると、ちょっと大変。
# Helper: one-row total employees by gender (within a size block)
# Assemble the table: spanners = company_size, columns = gender (Female/Male/Overall)
size_levels <-
svy2$variables$company_size %>%
droplevels() %>%
levels()
block_list <- size_levels %>%
purrr::map(function(sz) {
svy_sz <- svy2 %>% dplyr::filter(company_size == sz)
list(
make_block_total_employees(svy_sz),
make_block_cont(svy_sz),
make_block_cat(svy_sz, "age_class_top", "Age class (top)"),
make_block_cat(svy_sz, "industry_top", "Industry (top)")
) %>%
gtsummary::tbl_stack() %>%
gtsummary::modify_header(all_stat_cols() ~ "**{level}**")
}
)
table1a <- gtsummary::tbl_merge(tbls = block_list,
tab_spanner = paste0("**", size_levels, "**"))2.3 相関行列
Table 1とセットにする相関行列。連続量がいくつかあるなら、モデリング前に必ずつくる。後にあらわれる多重共線性を評価しないといけないから。多重共線性評価も重要だけれど、実は相関行列を眺めて、相関が高い変数同士は、あたりまえだけれど多重共線性が発生するもとになる。
ちなみに、やはり賃金同士に強い相関がある。1つだけ残してあとは除去してVIFを計算すると、十分に多重共線性を避けられそうなことがわかった。
今回は、労働者数で重み付けした。
Rコードには特に目新しいことはない。
vars_cont <- c(
"mean_age","length_of_service",
"actual_number_of_scheduled_hours_worked",
"actual_number_of_overtime_worked",
"contractual_cash_earnings",
"hourly_scheduled_cash_earnings",
"annual_special_cash_earnings"
) %>%
intersect(names(wc))
# Weighted Pearson
form <-
as.formula(paste("~", paste(vars_cont, collapse = " + ")))
cov_mat <-
survey::svyvar(
form, design = svy, na.rm = TRUE
) %>%
as.matrix()
cor_w <-
cov2cor(cov_mat) %>% round(3)
cor_w %>%
as.data.frame() %>%
tibble::rownames_to_column("Variable") %>%
gt::gt() %>%
gt::tab_caption("Table S1. Employee-weighted Pearson correlation matrix")
# Save CSV for appendix
readr::write_csv(as_tibble(cor_w, rownames = "variable"), "correlation_weighted_pearson.csv")