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")3 写真から位置情報を引っこ抜いて地図に投影する
タイにある工業団地にて集めた求人票写真がある。
たくさんあるが、どこで撮影したかは引用しないと。企業がある場所を記録することにもなる。 広すぎるから数日かけて回る
ということで、画像から位置情報を抜き取る。昔はtcl/tkを使って書いた記憶があるけれど、いまならexifr()という便利なライブラリがある。その他ライブラリを含めて、読み込み部分コードは省略。
はじめに、画像データがある場所をRに教える。here::here()を使うと、子ディレクトリ内も巡回してくれて便利。引数patternには、一通り画像ファイル拡張子を教える。もっとも、jpgしかないから、いらないかもしれない。(?i)を付して、大文字小文字を区別しないようにする。画像ファイル拡張子は.jpgだったり.JPGだったりしていまひとつ統一されない事がある一方、システムはしっかり大文字小文字を区別する。なので、区別しないでね、と教えることはとっても大事。
exif引数名候補もここで与える。たいていはGPSLatitudeとGPSLongitudeとがあれば足りるが、スマホ機種次第で多少違ったりする。なので、複数候補を与える。
このプロセスは、少々時間がかかる。画像1,500枚で数分程度。
#
## ---- load.data ----
imgs <- list.files(
path = here::here("img"),
recursive = TRUE, full.names = TRUE,
pattern = "(?i)\\.(jpe?g|heic|heif|png|tif?f|webp)$"
)
# obtain exif data from images
img_metadata <-
exifr::read_exif(imgs)
#
## ---- candidates.variables ----
candidate_lat <- c("GPSLatitude", "Composite.GPSLatitude", "EXIF.GPSLatitude")
candidate_lon <- c("GPSLongitude","Composite.GPSLongitude","EXIF.GPSLongitude")
candidate_dt <- c("DateTimeOriginal", "CreateDate", "ModifyDate", "FileModifyDate")
# 画像は、一度exif情報を取得後は削除する。特にOnedriveなど使う場合、容量制限があるから注意する。
次はお助け関数集。
exifが60進法でもつ緯度経度を10進法に変換する関数、時刻フォーマットを整えつつタイムゾーンをBKKにする関数、緯度経度が入っているであろう変数をもれなく抽出する関数。
#
## ---- helper.function ----
# 60deg. → decimal
dms_to_decimal <- function(x) {
if (is.numeric(x)) return(x)
if (is.na(x)) return(NA_real_)
s <- tolower(x)
nums <-
# ここがポイント。
# 数字部分のみを取り出す。
stringr::str_match(s, "([0-9\\.]+)\\D+([0-9\\.]+)?\\D*([0-9\\.]+)?")
# hhmmssを取り出す。
# この時点ではまだ60進法
deg <-
as.numeric(nums[,2]); min <- as.numeric(nums[,3]); sec <- as.numeric(nums[,4])
deg[is.na(deg)] <- 0; min[is.na(min)] <- 0; sec[is.na(sec)] <- 0
# 60進法から10進法へと変換
val <- deg + min/60 + sec/3600
# おまじない
if (str_detect(s, "[sw]")) val <- -val
val
}
# day
parse_dt_tidy <-
function(x, tz_out = "Asia/Bangkok") {
x <- as.character(x)
x <-
x %>%
stringr::str_replace_all("^([0-9]{4}):", "\\1-") %>%
stringr::str_replace_all(
"^([0-9]{4}-[0-9]{2}):",
"\\1-"
)
dt <- suppressWarnings(
lubridate::parse_date_time(
x,
orders = c("Y-m-d H:M:S", "Y/m/d H:M:S", "Ymd HMS"),
tz = "Asia/Bangkok"
)
)
with_tz(dt, tzone = tz_out)
}
# make a complete candidates' list
coalesce_from <-
function(.df, candidates) {
cols <-
# choose any of existing candidates (GPSLattitude...ect.) from .df
# and make a tibble data frame
dplyr::pick(
.df,
dplyr::any_of(candidates)
)
# When no candidates were there, return NA
# Otherwise, dplyr::coalesce() returns errors.
if (ncol(cols) == 0) return(NA)
# find the first non-NA element
# !!! is a splice operator splitting vector / list as argument
# coalesce(c(1,NA,3), c(10,20,NA)) -> c(1,20,3)
# (tibble's style is LIST. Errors will be returned if not be split)
dplyr::coalesce(!!!cols)
}dplyr::coalesce()は、ちょっとわかりにくい。日頃使うところもないし。以下がわかりやすいかも、と書いておいてなんだか実はよくわからん。
library(tidyverse)
# おためしデータ
df <-
dplyr::tibble(
prefecture = c("長崎", NA, NA),
PREF_NAME = c(NA, "佐賀", NA),
`都道府県名` = c(NA, NA, "福岡")
)
# 都道府県名が入った変数名候補
std_pref <- c("prefecture","Pref","PREF_NAME","都道府県名")
# ふさわしい都道府県名を入れた新しい変数を追加する。
df %>%
# 変数を1つ追加する
dplyr::mutate(
pref_std = {
# 中間変数(?)としてとりあえずac。
# std_prefに入ってる都道府県名中いずれかにあてはまる変数を抽出。
ac <- dplyr::across(
dplyr::any_of(std_pref)
)
# そんな変数なければNA
if (ncol(ac) == 0) NA_character_
# そんな変数あればdplyr::coalesce()にて
# pref_std内容になる観測値を返す。
# 内容は、ac中std_prefにあてはまる変数について、最初にあらわれる都道府県名。
# purrr::reduce()を使って、最低限の列（変数）のみ残るようにする。
else purrr::reduce(
ac,
dplyr::coalesce
)
}
)おつぎは本番。必要なところろ抜き出す。またもやcoalesce()を頻用しまくり。だって便利なんですもの。
#
## ---- refine.data ----
img_metadata_refined <-
img_metadata %>%
dplyr::mutate(
lat_raw = {
ac_df <- across(any_of(candidate_lat))
if (ncol(ac_df) == 0) NA else do.call(dplyr::coalesce, ac_df)
},
lon_raw = {
ac_df <- across(any_of(candidate_lon))
if (ncol(ac_df) == 0) NA else do.call(dplyr::coalesce, ac_df)
},
dt_raw = {
ac_df <- across(any_of(candidate_dt))
if (ncol(ac_df) == 0) NA else do.call(dplyr::coalesce, ac_df)
}
) %>%
dplyr::mutate(
latitude = suppressWarnings(as.numeric(lat_raw)),
longitude = suppressWarnings(as.numeric(lon_raw))
) %>%
dplyr::mutate(
latitude = ifelse(is.na(latitude), dms_to_decimal(lat_raw), latitude),
longitude = ifelse(is.na(longitude), dms_to_decimal(lon_raw), longitude)
) %>%
dplyr::mutate(
taken_at_bkk = parse_dt_tidy(dt_raw, tz_out = "Asia/Bangkok")
) %>%
dplyr::transmute(
source_file = dplyr::coalesce(SourceFile, FileName),
latitude, longitude,
altitude = GPSAltitude,
taken_at_bkk
) %>%
dplyr::filter(!is.na(latitude), !is.na(longitude)) %>%
sf::st_as_sf(coords = c("longitude", "latitude"), crs = 4326, remove = FALSE)
# save
readr::write_excel_csv(
sf::st_drop_geometry(img_metadata_refined),
"img_metadata.csv"
)