Rで日本語の都道府県名を英語（ローマ字のヘボン式）表記に変換する関数を作成した

Rで日本語の都道府県データを取り扱う際、「秋田県」や「千葉」といった日本語表記を「Akita」、「Chiba」などの英語（ローマ字）表記に変換したい場面がよくあると思う（私はある）。

しかし、単純な文字置換だと「都・府・県」がついているものとついていないものが混在していたり、都道府県以外の地名（「横浜市」や「合計」など）が混じっていてエラーや意図しないNAになってしまったりすることがあり困ることがある（行政が作成るデータで頻繁に遭遇する）。

そもそも、47都道府県名の文字列置換処理を毎回書くのが綴りミスなど気にかけることが多すぎて（内容は単純だが）大変すぎる。

そこで、自作Rパッケージinfunに、これらの表記揺れや非対応テキストを柔軟に処理できる軽量な関数to_eng_pref()を追加したので自分用にもメモしておく。

github.com

主な特徴

• 「都・府・県」の表記揺れを自動吸収: 「秋田県」「秋田」のどちらが入力されても正確にAkitaへ変換できる（「北海道」はそのまま保持）。
• 非対応テキストの柔軟な処理: 都道府県以外の地名が混ざっている場合、デフォルトでは元の文字列を維持する。必要に応じてNAや任意の代替文字列（例:"Other"）を返す指定も可能。

インストール方法

GitHubから{infun}パッケージをインストールすることで使える。

# remotesパッケージを使ってGitHubからインストール
remotes::install_github("indenkun/infun")

使い方・コード例

基本的な使い方は非常にシンプル。

library(infun)
# 基本的な使用例（都道府県以外の地名はそのまま維持）
test_data <- c("秋田県", "横浜市", "千葉", "合計", "東京都")

to_eng_pref(test_data)

## [1] "Akita"  "横浜市" "Chiba"  "合計"   "Tokyo"

該当しない要素の制御（default 引数）

都道府県名として認識されなかったデータに対する挙動はdefault引数で制御できる。

# 未マッチの要素を NA に置換したい場合
to_eng_pref(test_data, default = NA)

## [1] "Akita" NA      "Chiba" NA      "Tokyo"

# 未マッチの要素を特定の文字列に置換したい場合
to_eng_pref(test_data, default = "Other")

## [1] "Akita" "Other" "Chiba" "Other" "Tokyo"

おわりに

今月何も書いてないなと思って書いた。

データクレンジングや可視化の事前処理において、地名データの英語化（ローマ字化）にお悩みの方はぜひお試しください。

今どきならこういう単純で簡単な処理であれば毎回AIに書かせるのもありかもしれないなと思いながら、毎回意図した挙動をしてくれる（例えばマッチしなかったらどうするかの処理など）ようにするにはひと手間がかかることが割とある（要件定義の言語化が不十分なのが大きな要素だと思うが）ので、関数として持っておくのはいいと思う。

{modelsummary}でGOF統計量とp値の表示を調整するメモ

Rで{modelsummary}で表を作るときに何度も調べ直しているので、GOF統計量の名前とp値の表示方法についてメモしておく。

GOF統計量の内部名は抽出方法によって異なる

AIC, BIC, R², ICCなどのGOF統計量の表示名や丸め桁数は、modelsummary() の gof_map 引数で調整できる。

ただし、gof_map の raw には、実際に modelsummary が内部で保持している名前を指定する必要がある。この名前は、モデルの種類やGOF統計量の抽出に使われるバックエンドによって異なることがある。

たとえば、以下の混合モデルの例では、modelsummary_get の既定設定では AIC と BIC が aic, bic と小文字で返される。

これはget_gof関数で確認できる。

library(lmerTest)

## Loading required package: lme4

## Loading required package: Matrix

## 
## Attaching package: 'lmerTest'

## The following object is masked from 'package:lme4':
## 
##     lmer

## The following object is masked from 'package:stats':
## 
##     step

library(modelsummary)

fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
# aicが小文字
get_gof(fm1)

##        aic      bic r2.conditional r2.marginal       icc     rmse nobs
## 1 1755.628 1774.786      0.7992199   0.2786511 0.7216603 23.43805  180

一方、broom.mixed を読み込み、options(modelsummary_get = "all") として複数のバックエンドからGOF統計量を取得・結合すると、AICとBICは AIC, BIC と大文字で返される。その他の行を表示する

library(broom.mixed)
options(modelsummary_get = "all")
# AICが大文字
get_gof(fm1)

##   nobs   sigma    logLik      AIC      BIC REMLcrit df.residual r2.conditional
## 1  180 25.5918 -871.8141 1755.628 1774.786 1743.628         174      0.7992199
##   r2.marginal       icc     rmse
## 1   0.2786511 0.7216603 23.43805

したがって、この場合に gof_map でAICとBICの表示や丸めを指定するには、raw = "aic" や raw = "bic" ではなく、raw = "AIC" と raw = "BIC" を指定する必要がある。その他の行を表示する

p値を小数点以下3桁で出力し、p < 0.001も使用する

p値をアスタリスク等のマークではなく実際の値を小数点以下3桁で記載する必要が多くの場合であり、その場合に、それぞれの出力方式ごとに手法は若干異なるが、p = 0.000をp < 0.001に直接値を書き換える必要がある。

{gt}フォーマットの出力の場合の例を書く。

library(gt)
library(stringr)
modelsummary(fm1,
             output = "gt",
             # 小数点以下3に制御
             fmt = "%.3f",
             statistic = "p = {p.value}") |> 
  # 小数点以下3桁に制御したおかげでp < 0.001の場合p = 0.000と出力されるようになるので
  # p = 0.000をp < 0.001に書き換える
  text_transform(
    locations = cells_body(),
    fn = function(x) {
      x <- str_replace_all(x, "p = 0.000", "p < 0.001")
      return(x)
    }
  )

	(1)
(Intercept)	251.405
	p < 0.001
Days	10.467
	p < 0.001
sd__(Intercept)	24.741
sd__Days	5.922
cor__(Intercept).Days	0.066
sd__Observation	25.592
Num.Obs.	180
R2 Marg.	0.279
R2 Cond.	0.799
AIC	1755.6
BIC	1774.8
ICC	0.7
Log.Lik.	-871.814
RMSE	23.44
REMLcrit	1743.62827195995

結論

ついでに95％信頼区間もつけておく

library(lmerTest)
library(modelsummary)
library(performance)
library(parameters)

## 
## Attaching package: 'parameters'

## The following object is masked from 'package:modelsummary':
## 
##     supported_models

library(broom.mixed)
library(gt)
library(stringr)

options(modelsummary_get = "all")

fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)

coef_labels <- c("(Intercept)" = "Intercept", "Days" = "Days")

gof_labels <- list(
  list("raw" = "nobs", "clean" = "Observations (N)", "fmt" = "%.0f"),
  list("raw" = "r2.marginal", "clean" = "Marginal R² (Fixed Effects)", "fmt" = "%.3f"),
  list("raw" = "r2.conditional", "clean" = "Conditional R² (Total Model)", "fmt" = "%.3f"),
  list("raw" = "AIC", "clean" = "AIC", "fmt" = "%.1f"),
  list("raw" = "BIC", "clean" = "BIC", "fmt" = "%.1f"),
  list("raw" = "icc", "clean" = "ICC", "fmt" = "%.2f"),
  list("raw" = "logLik", "clean" = "Log Likelihood", "fmt" = "%.1f")
)

modelsummary(
  models = fm1,
  output = "gt",
  coef_map = coef_labels,
  fmt = "%.3f",
  gof_map = gof_labels,
  estimate = "{estimate} [{conf.low}, {conf.high}], p = {p.value}"
) |> 
  text_transform(
    locations = cells_body(),
    fn = function(x) {
      x <- str_replace_all(x, "p = 0.000", "p < 0.001")
      return(x)
    }
  )

	(1)
Intercept	251.405 [237.006, 265.804], p < 0.001
	(6.825)
Days	10.467 [7.206, 13.729], p < 0.001
	(1.546)
Observations (N)	180
Marginal R² (Fixed Effects)	0.279
Conditional R² (Total Model)	0.799
AIC	1755.6
BIC	1774.8
ICC	0.72
Log Likelihood	-871.8

ポイントは、gof_map の raw に指定する名前を事前に get_gof() で確認することである。特に options(modelsummary_get = "all") を使う場合、AIC/BICの内部名が AIC/BIC になることがあるため注意する。

上記のように raw = "AIC"、raw = "BIC" と指定すると、AICとBICの丸めも期待どおりに反映される。

Rのmodelsummaryで混合モデルのR²とlogLikが同時に出ない問題の対処メモ

Rの{modelsummary}パッケージで、{lmerTest}を用いて作成した混合モデルを表形式で出力したところ、決定係数（R²）は表示されるのに対数尤度（logLik）が出力されないという現象に遭遇したので、その対処法をメモしておく。

再現コード

library(lmerTest)
library(modelsummary)

# 簡単なモデル
m <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)

# そのまま出力
modelsummary(m, output = "markdown")

(Intercept)	251.405
	(6.825)
Days	10.467
	(1.546)
SD (Intercept Subject)	24.741
SD (Days Subject)	5.922
Cor (Intercept~Days Subject)	0.066
SD (Observations)	25.592
Num.Obs.	180
R2 Marg.	0.279
R2 Cond.	0.799
AIC	1755.6
BIC	1774.8
ICC	0.7
RMSE	23.44

この時、出力にはR²は含まれるが、logLikは表示されない。

原因

ドキュメントを読むと、modelsummaryはモデル情報の抽出に際して以下の2系統のパッケージを利用していることが分かる。

Model extraction functions modelsummary can use two sets of packages to extract information from statistical models: the easystats family (performance and parameters) and broom. By default, it uses easystats first and then falls back on broom in case of failure. You can change the order of priorities or include goodness-of-fit extracted by both packages by setting:

options(modelsummary_get = “easystats”)

options(modelsummary_get = “broom”)

options(modelsummary_get = “all”)

The “all” option (default) means easystats then broom.

とのことで、詳しいことはここには書いていないが、{broom.mixed}パッケージを呼び出して、その後options(modelsummary_get = "all")とすると良いよう。

解決方法

{broom.mixed}をロードした上で、取得元を両方使うように設定する。

# おまじない2行
library(broom.mixed)
options(modelsummary_get = "all")
modelsummary(m, output = "markdown")

(Intercept)	251.405
	(6.825)
Days	10.467
	(1.546)
sd__(Intercept)	24.741
sd__Days	5.922
cor__(Intercept).Days	0.066
sd__Observation	25.592
Num.Obs.	180
R2 Marg.	0.279
R2 Cond.	0.799
AIC	1755.6
BIC	1774.8
ICC	0.7
Log.Lik.	-871.814
RMSE	23.44
REMLcrit	1743.62827195995

パッケージを呼び出さないと行けないのがミソ。

正直なところ、AICやBICが出力されていれば、実務上は十分ではないかとも思っていた。

実際、AICやBICはlogLik（対数尤度）にペナルティ項を加えた指標であり、

• AIC = -2 × logLik + 2k
• BIC = -2 × logLik + log(n) × k

という関係にあるため、モデル比較という目的であれば、これらが出ていれば十分なケースも多い。

しかし、レビュー等で「logLikも明示的に出しておくように」と指摘を受けた。

余談

AIに聞いたところ「計算されているはずなので手動で取り出して追加すればよい」といった力技を提案されたが、さすがに違和感があったのでドキュメントを確認したところ、解決した問題だった。AIに聞いてハマって、無駄に時間を費やしてしまった。最初からドキュメント読んでおいたほうが圧倒的に早かった。

やはり公式ドキュメントを読むのは大事。

R 4.6.0からパッケージのREADME.md が同梱されるようになったが……

R 4.6.0がリリースされ、巷では新機能や挙動変更の話題で盛り上がっている……はず。

NEWSを眺めていて、次のような記述に目が留まった。

Package README.md files are now installed and featured in HTML help.

どうやら、パッケージのREADME.md がパッケージ本体に同梱され、HELPから参照できるようになったらしい。長年、READMEをローカルでも確認できたらいいと思っていただけに嬉しい新機能だと思った。

早速、RStudioから確認してみる。

……うーん。

確かに README.md は見られるのだが、生のMarkdownファイルがそのまま表示されるだけで、正直かなり読みにくいと思った。Markdown のコメントアウト部分までそのまま表示されてしまうので、情報量が多いREADMEほど視認性が悪くなる印象がある（試しにPositron（RStudio 以外の環境）でも確認してみたが、当然状況は同様だった）。

もちろん、ブラウザ側でMarkdownをきれいに表示してくれる拡張機能を使ったり、Markdownビューア系のアドオンを併用すれば、それなりに見やすくすることはできる。しかし、HELP から直接参照できる機能として考えるとそこまでユーザー側で工夫しないと読めないというのは、やや惜しい気もする。

正直なところ、VignettesのようにCRANのパッケージのページで表示できるREADMEのようにHTMLに変換された形で表示されてくれた方が嬉しかった（あるいは.mdとhtmlの2つ選べる状態）、というのが率直な感想である。

もう一つ、ついでにCRAN 配布パッケージに README.md（Markdown）として日本語を含むファイルを同梱できるのであれば、「それなら、パッケージのコード（.R）やヘルプファイル（.Rd）にも、日本語を含む非 ASCII 文字をもっと自由に書けてもよさそうでは？」とも思ってしまう。

実際、R自体はすでにUTFネイティブになっており、技術的には日本語をネイティブに扱える環境は整っている。それでも現状では、CRAN チェックとの兼ね合いで、コード中の日本語、ヘルプファイル中の日本語については、エスケープするなど今なお慎重にならざるを得ないのが実情でここらへんは少し手間を感じる部分となっている。

R 4.6.0からwilcox.test()関数の仕様が変わってタイがあっても正確なp値が計算できるようになった

R 4.6.0がリリースされ巷では新しいRの話題で持ちきりのはず。

R 4.6.0 の NEWS を読んでいたところ、次のような記述があった。

wilcox.test() can now perform exact (conditional) inference in case of
ties. Based on contributions by Torsten Hothorn.

これによると、R 4.6.0からはwilcox.test()において、タイ（同じ値）が存在する場合でも、ウィルコクソンの順位和検定・符号付順位検定の正確な p 値を計算できるようになった模様。

早速、実際に挙動を確認してみる。

# R 4.6.0 環境
# Introductory Statistics with R で使用されている
# wilcox.test にタイが含まれる例
daily.intake <- c(5260, 5470, 5640, 6180, 6390, 6515,
                  6805, 7515, 7515, 8230, 8770)
wilcox.test(daily.intake, mu = 7725)

##
##  Wilcoxon signed rank exact test
##
## data:  daily.intake
## V = 8, p-value = 0.02246
## alternative hypothesis: true location is not equal to 7725

一方、R 4.5.3 環境で同じコードを実行すると次のようになる。

# R 4.5.3 環境での実行結果
wilcox.test(daily.intake, mu = 7725)

## Warning in wilcox.test.default(daily.intake, mu = 7725):
##   cannot compute exact p-value with ties
##
##  Wilcoxon signed rank test with continuity correction
##
## data:  daily.intake
## V = 8, p-value = 0.0293
## alternative hypothesis: true location is not equal to 7725

R 4.5.3では、データにタイが含まれているため正確な p 値を計算できないという Warning が出力され、近似（連続性補正付き）の検定が実行される。その結果、R 4.6.0 と比べて p 値も異なっている。

この「タイが存在するとwilcox.test()では正確な p 値が計算できない」という挙動は、これまでの stats パッケージにおける 仕様（ここ重要）であり、従来はcoin やexactRankTests といった別パッケージを用いて正確検定を行うのが一般的だった。

例えば、stats::wilcox.test()とほぼ同じ感覚で使えるexactRankTests::wilcox.exact()を用いると次のようになる。

exactRankTests::wilcox.exact(daily.intake, mu = 7725)

##
##  Exact Wilcoxon signed rank test
##
## data:  daily.intake
## V = 8, p-value = 0.02246
## alternative hypothesis: true mu is not equal to 7725

結果は、R 4.6.0のwilcox.test()と一致している。

正確なp値がstats::wilcox.test()だけで得られるようになったのはいい変更だと思う。その一方で、Rのバージョンを上げただけで、同じコード・同じデータから得られる p 値が変わる点には、多少注意が必要だと思う。

この「タイがあるために正確検定ができない」という問題自体は、かなり以前から知られており、回避手法やそのためのパッケージも長らく存在していた。それにもかかわらず、今回、基本的な統計で用いる関数であるwilcox.test()の挙動として結果が変わるのは、R Core Team の比較的保守的なメンテナンス方針を考えると、やや意外に感じる部分だった。

再現性や過去との比較が重要な解析では、R のバージョン差による挙動の変化として、この点を意識しておく必要がありそう。

ggplot2の棒グラフの上や下に値をラベルとして書き込む方法のメモ

ggplot2の棒グラフの上や下に値をラベルとして書き込む方法のメモ。とりあえず簡単なgeom_col()の場合で書いていく。

上に書くだけならgeom_text()で位置を調整して書けばいい。

例えば次のようなデータの場合。

library(ggplot2)
df <- data.frame(X = LETTERS[1:4],
                 value = c(10, 20, 30, 40))
df

##   X value
## 1 A    10
## 2 B    20
## 3 C    30
## 4 D    40

ggplot(df, aes(x = X, y = value)) +
  geom_col()

この棒グラフの上に値を書き込みたいと思いたい場合、geom_text()を使い、そのままの位置に書くと棒グラフに被るので、棒グラフの上に描画するようにvjustやnudge_yを使ってラベルの位置を調整しながら次のようにかける。

ggplot(df, aes(x = X, y = value)) +
  geom_col() +
  geom_text(aes(label = value), nudge_y = 2)

ggplot(df, aes(x = X, y = value)) +
  geom_col() +
  geom_text(aes(label = value), vjust = -0.5) +
  ylim(0, 42)

nudge_yとvjustの違いは上のコードを見てもらえば分かる通り、nudge_yの場合はいい塩梅に描画範囲を調整してくれるが、vjustは描画範囲を調整してくれないので、自分で調整する必要がある。

じゃあ、いつもnudge_yでいいじゃないかとなるかもしれないが一部画面ではvjustのほうが都合がいい場合がある。

例えば次のようなデータで横棒グラフを描く場合を考える。

df <- data.frame(X1 = rep(LETTERS[1:4], 3),
                 X2 = rep(LETTERS[24:26], each = 4),
                 value = seq(10, by = 10, length.out = 12))
head(df)

##   X1 X2 value
## 1  A  X    10
## 2  B  X    20
## 3  C  X    30
## 4  D  X    40
## 5  A  Y    50
## 6  B  Y    60

ggplot(df, aes(x = X1, y = value, fill = X2)) +
  geom_col(position = position_dodge())

このとき、nudge_yだと次のようになる。

ggplot(df, aes(x = X1, y = value, fill = X2)) +
  geom_col(position = position_dodge()) +
  geom_text(aes(label = value), position = position_dodge() ,nudge_y = 2)

## Error in `geom_text()`:
## ! Both `position` and `nudge_x`/`nudge_y` are supplied.
## ℹ Only use one approach to alter the position.

position = position_dodge()を使ってずらしたいと思ったとおりにいかない（ここでposition = position_dodge()をしないと積み上げグラフの位置で数字が書かれてしまう）。ドキュメントにもちゃんと、positionとnudge_xは一緒に使えないと書いている。

ここで、vjustを使うと思ったとおりに行く。

ggplot(df, aes(x = X1, y = value, fill = X2)) +
  geom_col(position = position_dodge()) +
  geom_text(aes(label = value), 
            position = position_dodge(0.85), 
            vjust = -0.5) +
  ylim(c(0, 130))

こんな感じに行く。

さらに正負の値が混じった次のようなデータを考える。

df <- data.frame(X = c("A", "B", "C", "D"),
                 value = c(10, 20, -10, -20))

ggplot(df, aes(x = X, y = value)) + 
  geom_col()

このとき、正の値の場合は棒グラフの上に、負の値の場合は棒グラフの下に描画したいと思った場合を考える。単純にnudge_yなどで位置を書くと正負の場合分けができない。

ggplot(df, aes(x = X, y = value)) + 
  geom_col() +
  geom_text(aes(label = value), nudge_y = 2)

結論から言うとifelse()などで条件づけするとかける。

ggplot(df, aes(x = X, y = value)) + 
  geom_col() +
  geom_text(aes(label = value), nudge_y = ifelse(df$value >= 0, 2, -2))

これで思ったグラフがかける。

Rのggplot2でデータにない値を凡例に反映させるメモ2

ggplot2パッケージの3.5.x以降で仕様変更があり、データ中に存在しない値を凡例に反映させるためには、以前の手法ではうまく行かなくなったので、現在動く手法のメモ。

例としてmtcarsデータの気筒数を因子型にしたものを用意して、グラフ化する。

library(ggplot2)
mtcars <- mtcars |> 
  dplyr::mutate(cyl = factor(cyl, levels = 4:8))
ggplot(mtcars) +
  geom_point(aes(x = mpg, y = disp, color = cyl))

データ中に、5気筒、7気筒がないので凡例には反映されていない。でも、5気筒や7気筒の車も世の中にあるわけなので、データ中にないだけで凡例には凡例には合ったほうがいいなとおもうこともある。

ドキュメントを一見すると、scale_color_*()やscale_fill_*()などの中でdrop = FALSEと指定することで存在しないデータを凡例に反映することができそうだが、これだけでは次の通り、凡例に数字は反映されているがポイントは反映されなくなってしまう。

ggplot(mtcars) +
  geom_point(aes(x = mpg, y = disp, color = cyl)) +
  scale_color_hue(drop = FALSE)

これを解決するには、追加でgeom_*()でshow.legend = TRUEとするといい。

ggplot(mtcars) +
  geom_point(aes(x = mpg, y = disp, color = cyl), show.legend = TRUE) +
  scale_color_hue(drop = FALSE)

ドキュメントを読むとちゃんと書いてるのだが、わかりにくく、ISSUEでも時折取り上げられている。

github.com

完全に存在しない値を加えたいならlimitsで細かく指定するといく。

ggplot(mtcars) +
  geom_point(aes(x = mpg, y = disp, color = cyl), show.legend = TRUE) +
  scale_color_hue(limits = c(4:8, "ALL"))

参考までに以前の手法は次の通り

indenkun.hatenablog.com