2024-02-26

ロジスティック回帰に最尤推定量が存在するか判定する

引き続きロジスティック回帰について調べている。

前回、ロジスティック回帰の最尤推定量にはバイアスがあることを調べた。

ロジスティック回帰の最尤推定量にはバイアスがある - ほくそ笑む

このバイアスは、サンプルサイズが大きい場合は無視できるが、入力変数の数に対してサンプルサイズが小さい場合には無視できないほど大きくなる。

今回はサンプルサイズをさらに小さくすると起こる問題について考える。

前回の最初の例を少しだけ変えて実行してみよう。この例では、入力変数の数 $p = 300$ として、パラメータ 300個の真の値を、最初の 100個は $\beta = 10$ 、次の 100個は $\beta = -10$ 、残りの 100個は $\beta = 0$ と設定した。ここまでは同じだが、変更点として、前回はサンプルサイズ $n = 1500$ だったのに対して、今回は $n = 1300$ で実行する。

n <- 1300
p <- 300

# データの生成
set.seed(314)
x <- rnorm(n * p, mean = 0, sd = sqrt(1/n))
X <- matrix(x, nrow = n, ncol = p)
beta <- matrix(rep(c(10, -10, 0), each = p/3))
prob <- plogis(X %*% beta)
y <- rbinom(n, 1, prob)

# ロジスティック回帰モデルの適用
fit <- glm(y ~ X, family = binomial, control = list(maxit = 50))

警告メッセージ: 
glm.fit: 数値的に 0 か 1 である確率が生じました

df <- data.frame(index = seq_len(p), coef = fit$coefficients[-1])
library(ggplot2)
theme_set(theme_bw())

ggplot(df, aes(index, coef)) +
  geom_point(color = "blue") +
  annotate("segment", x = c(1, 101, 201), xend = c(100, 200, 300),
           y = c(10, -10, 0), yend = c(10, -10, 0), size = 1.5) +
  xlab("Index") + ylab("Coefficients (true and fitted)")

実行結果は見てのとおり、推定されたパラメータが 1e+16 (= 10¹⁶) などと、とんでもない値になっていることがわかる。これは最尤推定量のバイアスとは異なる原因で起こっている。

何が起こっているかというと、実は、このデータに対して、ロジスティック回帰モデルには最尤推定量が存在しない。ロジスティック回帰モデルには最尤推定量が存在しない場合があるのである。しかし、R言語の最尤推定値を求めるアルゴリズムは、最尤推定量が存在しない場合にも、なんらかの値を返してしまう。

まとめると、入力変数の数に対してサンプルサイズが小さい場合、推定値には2つのケースが存在する。

最尤推定量が存在する場合、バイアスのある推定値
最尤推定量が存在しない場合、なんの意味もない値

この記事では、この2つのケースを判別する方法、すなわち、ロジスティック回帰に最尤推定量が存在するかどうかを判定する方法を紹介する。

判定材料

ロジスティック回帰に最尤推定量が存在しない場合、パラメータの推定値が異常な値を取るため、すぐにわかるのではないか？と思ってしまうが、そうでもない。

上記の例でサンプルサイズ $n$ を変えて推定されたパラメータの最初の100個 ( $\beta=10$ ) がどのような値を取るかを見てみよう。

$n = 1300$ でだけ、推定値が異常な値を取ることがわかる。

$n$ の小さい範囲を拡大してみると、パラメータが大きめに推定されていることがわかるが、これは真の値が 10 であることを知っているため大きいとわかるのであって、真の値を知らない場合はこの推定値が大きいかどうかの判断は難しい。

判定方法

ロジスティック回帰に最尤推定量が存在するかどうかを判定する方法は、40年ほど前に研究された [Albert & Anderson, 1984][Santer & Duffy, 1986]。

文献 [Albert & Anderson, 1984] によれば、観測データが完全分離 (complete separation) または準完全分離 (quasi-complete separation) の場合に限り、最尤推定量は存在しない。

完全分離とは、入力変数が作る空間に、出力変数を完全に分離できる超平面が存在することを言う。すなわち、入力変数 $X$ 、出力変数 $y$ に対して、

\begin{cases} Xb > 0 & y = 1 \\ Xb < 0 & y = 0 \\ \end{cases}

を満たす $b$ が存在することを言う。ここで、

\begin{align} 2y - 1 = \begin{cases} 1 & y = 1 \\ -1 & y = 0 \\ \end{cases} \end{align}

であるので、

\begin{align} Xb \times (2y - 1) > 0 \end{align}

が成り立つ $b$ が存在すればよい。

これは、線形計画問題に実行可能解があるかどうかで判定できる。

ROI パッケージを使えば線形計画問題を簡単に解くことができる。

# ROI パッケージで線形計画問題を解く
library(ROI)

A <- cbind(X, 1) * (2 * y - 1)
b <- rep(0.001, n)
c <- rep(1, p + 1)

lp <- ROI::OP(objective = ROI::L_objective(c), 
              constraints = ROI::L_constraint(A, rep(">=", n), b),
              bounds = ROI::V_bound(lb = rep(-Inf, p + 1)))
solution <- ROI::ROI_solve(lp)

# 実行可能解 (feasible solution) があれば完全分離可能
status_code <- solution$status$msg[["code"]]
STATUS_CODE_NO_FEASIBLE_SOLUTION_EXISTS <- 4L
STATUS_CODE_SOLUTION_IS_OPTIMAL <- 5L
STATUS_CODE_SOLUTION_IS_UNBOUNDED <- 6L
is_separable <- 
  if (status_code == STATUS_CODE_NO_FEASIBLE_SOLUTION_EXISTS) {
    FALSE
  } else if (status_code == STATUS_CODE_SOLUTION_IS_OPTIMAL) {
    TRUE
  } else if (status_code == STATUS_CODE_SOLUTION_IS_UNBOUNDED) {
    TRUE
  }

ここでは簡易的に完全分離を判定する方法を書いたが、完全分離、準完全分離、そしてオーバーラップ（最尤推定量が存在する）のいずれであるかを判別する線形計画問題の定式化が文献 [Santer & Duffy, 1986] に載っているので参照してほしい。

文献 [Konis, 2007] では、より洗練された方法が提案されている。これを実装したのが detectseparation パッケージである。このパッケージでは完全分離または準完全分離かどうかを簡単に判定できる。

library(detectseparation)
result <- detect_separation(X, y, family = binomial())
is_separable <- result$outcome

完全分離または準完全分離であれば最尤推定量は存在しないので、これでロジスティック回帰に最尤推定量が存在するかどうかを判定できる。

参考文献

Albert, Adelin & Anderson, J.A.. (1984). On the Existence of Maximum Likelihood Estimates in Logistic Regression Models.
Santner, Thomas & Duffy, Diane. (1986). A Note on A. Albert and J. A. Anderson's Conditions for the Existence of Maximum Likelihood Estimates in Logistic Regression Models.
Konis, K.P. (2007). Linear programming algorithms for detecting separated data in binary logistic regression models.

2024-01-23

ロジスティック回帰の最尤推定量にはバイアスがある

ロジスティック回帰について調べている。

ロジスティック回帰モデルのパラメータの最尤推定量は、不偏推定量ではなく、バイアスがある。

例として、サンプルサイズ $n = 1500$ 、入力変数の数 $p = 300$ のときを考える。パラメータ 300個の真の値を、最初の 100個は $\beta = 10$ 、次の 100個は $\beta = -10$ 、残りの 100個は $\beta = 0$ に設定して推定してみよう。

n <- 1500
p <- 300

# データの生成
set.seed(314)
x <- rnorm(n * p, mean = 0, sd = sqrt(1/n))
X <- matrix(x, nrow = n, ncol = p)
beta <- matrix(rep(c(10, -10, 0), each = p/3))
prob <- plogis(X %*% beta)
y <- rbinom(n, 1, prob)

# ロジスティック回帰モデルの適用
fit <- glm(y ~ X, family = binomial)

# パラメータの最尤推定値の抽出
df <- data.frame(index = seq_len(p), coef = fit$coefficients[-1])

library(ggplot2)
theme_set(theme_bw())

ggplot(df, aes(index, coef)) +
  geom_point(color = "blue") +
  annotate("segment", x = c(1, 101, 201), xend = c(100, 200, 300),
           y = c(10, -10, 0), yend = c(10, -10, 0), size = 1.5) +
  xlab("Index") + ylab("Coefficients (true and fitted)")

最初の 100個のパラメータは、真の値 $\beta = 10$ よりも大きめに推定されている。次の 100個のパラメータは、真の値 $\beta = -10$ よりも小さめに推定されている。最後の 100個のパラメータは真の値 ( $\beta = 0$ ) の周りに均等にばらついている。

つまり、ロジスティック回帰のパラメータの最尤推定量には、絶対値を大きめに推定するようなバイアスが乗っている。

最尤推定量には一致性と漸近正規性があるので、サンプルサイズ $n$ が十分大きければ、このバイアスは無視できるようになる。上記の例でサンプルサイズを変えて同じシミュレーションをしてみよう。

入力変数の数 $p = 300$ に対して $n = 6000$ まで増やせば、最尤推定量のバイアスは無視できるほど小さくなるようだ。

しかし、逆に言うと、入力変数の数に対してサンプルサイズが小さい場合は無視できないほどのバイアスが乗ってしまう。

このバイアスは、統計解析において次の問題を引き起こす。

効果量を過大に見積もってしまう
出力の確率を極端に見積もってしまう

上記の例はパラメータの真の値が作為的すぎるので、もう少し現実的な例で説明しよう。

真のパラメータ $\beta$ を正規分布 $N(3, 4)$ から抽出し、パラメータの真の値と最尤推定値 (MLE) の散布図を描いてみる。

n <- 4000
p <- 800

# データの生成
set.seed(314)
x <- rnorm(n * p, mean = 0, sd = sqrt(1/n))
X <- matrix(x, nrow = n, ncol = p)
beta <- rnorm(p, mean = 3, sd = 4)
prob <- plogis(X %*% beta)
y <- rbinom(n, 1, prob)

# ロジスティック回帰モデルの適用
fit <- glm(y ~ X, family = binomial)

# パラメータの最尤推定値の抽出
df <- data.frame(true = beta, mle = fit$coefficients[-1])

ggplot(df, aes(true, mle)) +
  geom_point(color = "blue") +
  geom_abline(slope = 1, size = 1) +
  xlab("True signal") + ylab("MLE")

黒線は傾き 1 の直線であり、もしバイアスがなければ推定値はこの線の周りに分布するはずである。しかし、明らかに分布がずれていることから、この例でも最尤推定量にバイアスがあることがわかる。

また、散布図の傾きが大きくなっていることから、絶対値が大きくなる方向にバイアスが乗っていることがわかる。

絶対値が大きくなる方向にバイアスの乗った推定値を、入力変数の効果量として解釈すると、その入力変数は過大評価され、誤った結論を導いてしまう恐れがある。

もう一つの問題は、このようなバイアスの乗った推定値で算出された出力の確率は、極端な値になりやすいということである。

出力 = 1 となる確率をパラメータの推定値から算出してみる。

x <- X %*% beta
pred <- predict(fit, type = "response")

df <- data.frame(x = x, prob = prob, pred = pred, pred2 = pred2)
ggplot(df, aes(x, pred)) +
  geom_point(color = "blue", alpha = 0.1) +
  geom_point(y = prob) +
  scale_x_continuous(breaks = NULL, limits = c(-6, 6)) +
  scale_y_continuous(breaks = 0:5 * 1/5, minor_breaks = NULL) +
  xlab(NULL) + ylab("Probabilities (True and predicted)")

左図はバイアスがある場合、右図はバイアスがない場合である。黒線が真の確率、青点が予測確率である。

左図は右図にくらべて、予測確率は 0 または 1 に寄っている。つまり、推定値にバイアスがある場合、予測確率は両極端に偏ってしまう。この予測確率をもとに「出力はほぼ確実に 1 である」と予測しても、それは推定値のバイアスのせいであって、実際はそれほど確実ではないかもしれない。

続き

hoxo-m.hatenablog.com

2024-01-08

ロジスティック回帰モデルの負の対数尤度

ロジスティック回帰モデルの負の対数尤度は次のように書ける。

\begin{align} NLL(\beta) = \sum_{i=1}^n F(x_i \beta) - y_i (x_i \beta) \end{align}

ただし、 $F(t)$ はロジスティック関数の積分

\begin{align} F(t) = \int \frac{1}{1 + e^{-t}} dt = \log(1 + e^{t}) \end{align}

である。

導出

入力変数 $x$
出力変数 $y$
回帰係数 $\beta$
逆リンク関数 $f(t) = \frac{1}{1 + e^{-t}} = \frac{e^{t}}{1 + e^{t}}$
モデル $P(y \mid x) = f(x\beta)$

尤度関数

\begin{align} Lik(\beta) = f(x\beta)^{y} (1-f(x\beta))^{(1-y)} \end{align}

$t = x\beta$ とおく

\begin{align} Lik(\beta) = f(t)^{y} (1-f(t))^{(1-y)} \end{align}

対数尤度関数

\begin{align} LogLik(\beta) &= \log(Lik(\beta)) \\ &= \log \left( f(t)^{y} (1-f(t))^{(1-y)} \right) \\ &= \log f(t)^{y} + \log (1-f(t))^{(1-y)} \\ &= y \log f(t) + (1-y) \log (1-f(t)) \\ &= y \log \frac{e^{t}}{1 + e^{t}} + (1-y) \log (1-\frac{e^{t}}{1 + e^{t}}) \\ &= y \log e^{t} - y \log (1 + e^{t}) + (1-y) \log (\frac{1 + e^{t} - e^{t}}{1 + e^{t}}) \\ &= y t - y \log (1 + e^{t}) + (1-y) \log (\frac{1}{1 + e^{t}}) \\ &= y t - y \log (1 + e^{t}) - (1-y) \log (1 + e^{t}) \\ &= y t - y \log (1 + e^{t}) - \log (1 + e^{t}) + y \log (1 + e^{t}) \\ &= y t - \log (1 + e^{t}) \\ \end{align}

負の対数尤度関数

\begin{align} NegLogLik(\beta) &= -LogLik(\beta) \\ &= \log (1 + e^{t}) - y t \\ &= F(t) - y t \\ &= F(x\beta) - y (x\beta) \\ \end{align}

2023-10-28

状態を持つループ処理を accumulate() でシンプルに書く

R言語のコミュニティ https://r-wakalang.slack.com で回答したのでメモ。

質問はこんな感じ（意訳しています）。

次のようなデータを以下のルールで処理したい。データを上から下に見ていき、 (1) before に TRUE が出たら、それ以降は after を TRUE にする。 (2) ただし、condition が FALSE になったら after を FALSE にして状態をリセットする。これを、for を使わないやり方で書きたい。（データにすでにある after は答えあわせ用）

   before condition after
 1 FALSE  FALSE     FALSE
 2 FALSE  TRUE      FALSE
 3 TRUE   TRUE      TRUE 
 4 FALSE  TRUE      TRUE 
 5 FALSE  TRUE      TRUE 
 6 FALSE  FALSE     FALSE
 7 FALSE  TRUE      FALSE
 8 FALSE  TRUE      FALSE
 9 TRUE   TRUE      TRUE 
10 FALSE  TRUE      TRUE

以下、データは次のコードで生成したものを使う（元の質問より簡略化しています）。

library(tidyverse)

# データの準備
d <- tibble(
  before    = c(F,F,T,F,F,F,F,F,T,F),
  condition = c(F,T,T,T,T,F,T,T,T,T), 
  after     = c(F,F,T,T,T,F,F,F,T,T),
)

for を使った書き方

これは、状態 (state) を持つ繰り返し処理であり、for を使って書ける。

データを1行ずつ上から処理する。初期値は before = FALSE, condition = FALSE でこのときの出力 after は FALSE
初期状態 (state = 0) で before が TRUE になると after を TRUE にして次の状態 (state = 1) に移る。before が FALSE ならば after は FALSE のままで状態も変わらない
次の状態 (state = 1) で condition が FALSE になると after を FALSE にして初期状態に戻る。condition が TRUE ならば after は TRUE のままで状態も変わらない

state <- 0  # 状態を初期化
for (i in seq_len(nrow(d))) {
  row <- d[i, ]  # 一行ずつ取得
  if (state == 0) {
    if (row$before == TRUE) {
      d[i, "after2"] <- TRUE
      state <- 1  # 次の状態に移る
    } else {
      d[i, "after2"] <- FALSE
    }
  } else {  # state == 1
    if (row$condition == FALSE) {
      d[i, "after2"] <- FALSE
      state <- 0  # 初期状態に戻る
    } else {
      d[i, "after2"] <- TRUE
    }
  }
}

結果は d の after2 に格納される。期待する結果である after と一致する。

   before condition after after2
 1 FALSE  FALSE     FALSE FALSE 
 2 FALSE  TRUE      FALSE FALSE 
 3 TRUE   TRUE      TRUE  TRUE  
 4 FALSE  TRUE      TRUE  TRUE  
 5 FALSE  TRUE      TRUE  TRUE  
 6 FALSE  FALSE     FALSE FALSE 
 7 FALSE  TRUE      FALSE FALSE 
 8 FALSE  TRUE      FALSE FALSE 
 9 TRUE   TRUE      TRUE  TRUE  
10 FALSE  TRUE      TRUE  TRUE

accumulate() の使い方

これを for などのループを使わずに書くにはどうすればよいだろうか？

R言語の通常のベクトル化された関数では、状態を持つ処理を行うのは困難である。

そこで、accumulate() 関数を使う。この関数は、ベクトルを1つずつ処理していく関数であり、その際に前の処理の結果を使うことができる（基本関数の Reduce() と同様）。

状態の扱いが問題だったので、状態をいったん出力し、その状態を使って after を求めるという方針を取ろう。まずは、前の状態を入力にとり、次の状態を結果として返す関数を作成する（上記の for の中身から state に関する部分を抜き出せば簡単にできる）。

compute_state <- function(state, before, condition) {
  if (state == 0) {
    if (before == TRUE) {
      state <- 1
    }
  } else {  # state == 1
    if (condition == FALSE) {
      state <- 0
    }
  }
  state
}

accumulate2() 関数を使って状態を計算するには次のように書く（accumulate() は1変数を受け取り、accumulate2() は2変数を受け取る）。

state_vec <- accumulate2(d$before, d$condition, compute_state, .init = 0)
# 0 0 0 1 1 1 0 0 0 1 1

length(state_vec)
# 11

ただし、得られる結果はデータの行数より1つ多い。これは初期値 (.init = 0) が先頭に含まれるためである。

そのため、現在の状態を得るには最後の値を取り除けばよい。

current_state_vec <- state_vec |> head(-1)
# 0 0 0 1 1 1 0 0 0 1

ちなみに、次の状態を得たい場合は先頭を取り除く。

next_state_vec <- state_vec |> tail(-1)
# 0 0 1 1 1 0 0 0 1 1

accumulate() を使った書き方

現在の状態 state を得ることができた。出力 after は state, before, condition を使って通常のベクトル化された関数で計算できる。

state が 0 のとき、before が TRUE の場合のみ after = TRUE
state が 1 のとき、condition が FALSE の場合のみ after = FALSE、すなわち condition が TRUE の場合のみ after = TRUE

d |>
  mutate(state = accumulate2(before, condition, compute_state, .init = 0) |> head(-1)) |>
  mutate(after3 = (state == 0 & before) | (state == 1 & condition))

   before condition after after2 state after3
 1 FALSE  FALSE     FALSE FALSE      0 FALSE 
 2 FALSE  TRUE      FALSE FALSE      0 FALSE 
 3 TRUE   TRUE      TRUE  TRUE       0 TRUE  
 4 FALSE  TRUE      TRUE  TRUE       1 TRUE  
 5 FALSE  TRUE      TRUE  TRUE       1 TRUE  
 6 FALSE  FALSE     FALSE FALSE      1 FALSE 
 7 FALSE  TRUE      FALSE FALSE      0 FALSE 
 8 FALSE  TRUE      FALSE FALSE      0 FALSE 
 9 TRUE   TRUE      TRUE  TRUE       0 TRUE  
10 FALSE  TRUE      TRUE  TRUE       1 TRUE

非常にシンプルに書くことができた。

余談

ここからは余談であるが、この処理はもっとシンプルにできる。

まず、状態を計算する関数 compute_state() は次のように書ける。

compute_state <- function(state, before, condition) {
  if (state == 0) {
    before
  } else {  # state == 1
    condition
  }
}

さらにシンプルにすると次のようになる。

compute_state <- function(state, before, condition) {
  (state == 0 && before) || (state == 1 && condition)
}

これをよく見ると、after3 を計算したときの式と全く同一であることがわかる。

すなわち、今回の問題では、状態の計算結果と出力が一致する。

この場合、状態を計算しなくても、直接出力を計算することができる。

compute_after <- function(after, before, condition) {
  (!after && before) || (after && condition)
}

d |>
  mutate(after4 = accumulate2(before, condition, compute_after, .init = FALSE) |> tail(-1))

   before condition after after2 after4
 1 FALSE  FALSE     FALSE FALSE  FALSE 
 2 FALSE  TRUE      FALSE FALSE  FALSE 
 3 TRUE   TRUE      TRUE  TRUE   TRUE  
 4 FALSE  TRUE      TRUE  TRUE   TRUE  
 5 FALSE  TRUE      TRUE  TRUE   TRUE  
 6 FALSE  FALSE     FALSE FALSE  FALSE 
 7 FALSE  TRUE      FALSE FALSE  FALSE 
 8 FALSE  TRUE      FALSE FALSE  FALSE 
 9 TRUE   TRUE      TRUE  TRUE   TRUE  
10 FALSE  TRUE      TRUE  TRUE   TRUE

Enjoy!

2023-09-10

LighitGBM で変数重要度を出すと列名が文字化けするときの解決策

LighitGBM で変数重要度を出すと列名が文字化けするという相談を受けたときの調査メモ。

解決策は3つある。

現象を再現

まずは文字化け現象を再現してみる。データとしてはこんな感じ。

Sys.setlocale(locale = "Japanese_Japan.932")

df_all <- iris
colnames(df_all) <- c("あ", "い", "う", "え", "お")
colnames(df_all)[1] <- iconv(colnames(df_all)[1], from = "SJIS", to = "UTF-8")
head(df_all)

   あ  い  う  え     お
1 5.1 3.5 1.4 0.2 setosa
2 4.9 3.0 1.4 0.2 setosa
3 4.7 3.2 1.3 0.2 setosa
4 4.6 3.1 1.5 0.2 setosa
5 5.0 3.6 1.4 0.2 setosa
6 5.4 3.9 1.7 0.4 setosa

データを表示したときには列名は文字化けしていないのがポイント。

これを LightGBM にかけて変数重要度を出す。

library(lightgbm)

train_data <- as.matrix(df_all[, 1:4])
label <- as.integer(df_all[, 5]) - 1
lgb_train <- lgb.Dataset(data = train_data, label = label)

params <- list(objective = "multiclass",
               num_class = 3,
               metric = "multi_logloss")

model <- lgb.train(params = params, data = lgb_train)

lgb.importance(model)

   Feature       Gain     Cover Frequency
1:  縺<86> 0.62786358 0.3758313 0.3721847
2:  縺<88> 0.33334783 0.2347444 0.2139640
3:  縺<84> 0.02101790 0.2002376 0.2246622
4:  縺<82> 0.01777069 0.1891867 0.1891892

列名 (Feature) が文字化けしてしまっている。

この原因は、列名の文字エンコーディングに UTF-8 と Shift_JIS が混在しているためである。

encoding <- colnames(df_all) |> stringi::stri_enc_detect() |>
  purrr::map_chr(~ .x$Encoding[1])
encoding

[1] "UTF-8"     "Shift_JIS" "Shift_JIS" "Shift_JIS" "Shift_JIS"

現実のデータ分析では、データソースの異なるデータを結合して使用することも多いため、このような現象が起こる場合がある。

ひとつでも混じり物があると、すべての列名が文字化けしてしまうというのは興味深い。

解決策1. ロケールを UTF-8 にする

最もおすすめしたい解決策は、ロケールの文字エンコーディングを UTF-8 に変更することだ。これができるなら一番早い。

Sys.setlocale(locale = "Japanese_Japan.utf8")

lgb.importance(model)

   Feature       Gain     Cover Frequency
1:      う 0.62786358 0.3758313 0.3721847
2:      え 0.33334783 0.2347444 0.2139640
3:      い 0.02101790 0.2002376 0.2246622
4:      あ 0.01777069 0.1891867 0.1891892

しかし、様々な事情があって、これが常に可能とは限らない。

解決策2. 学習時に列名を指定する

2番目におすすめなのは、学習時に列名を指定する方法だ。 lgb.train() 関数には colnames という引数がある。この引数に文字エンコーディングが統一された列名を入力する。列名をエディタ上で定義してしまえば文字エンコーディングは自然に統一される。

Sys.setlocale(locale = "Japanese_Japan.932")

col_names <- c("あ", "い", "う", "え")

model <- lgb.train(params = params, data = lgb_train,
                   colnames = col_names)

lgb.importance(model)

   Feature       Gain     Cover Frequency
1:      う 0.62786358 0.3758313 0.3721847
2:      え 0.33334783 0.2347444 0.2139640
3:      い 0.02101790 0.2002376 0.2246622
4:      あ 0.01777069 0.1891867 0.1891892

この方法のデメリットは、列数が多いときになかなか面倒なことである。

解決策3. 列名の文字エンコーディングを統一する

最もおすすめしない解決策は、列名の文字エンコーディングを統一する方法だ。

まずは、それぞれの列名の文字エンコーディングを調べる。

colnames(df_all) |> stringi::stri_enc_detect() |> 
  purrr::set_names(colnames(df_all)) |> head(4)

$あ
  Encoding Language Confidence
1    UTF-8                 0.8

$い
   Encoding Language Confidence
1 Shift_JIS       ja        0.1
2   GB18030       zh        0.1
3      Big5       zh        0.1

$う
   Encoding Language Confidence
1 Shift_JIS       ja        0.1
2   GB18030       zh        0.1
3      Big5       zh        0.1

$え
   Encoding Language Confidence
1 Shift_JIS       ja        0.1
2   GB18030       zh        0.1
3      Big5       zh        0.1

1番目の列名「あ」が UTF-8 であるのが文字化けの原因なので、これを Shift_JIS に変換すればよい。

Sys.setlocale(locale = "Japanese_Japan.932")

colnames(df_all)[1] <- iconv(colnames(df_all)[1], from = "UTF-8", to = "SJIS")

model <- lgb.train(params = params, data = lgb_train)

lgb.importance(model)

   Feature       Gain     Cover Frequency
1:      う 0.62786358 0.3758313 0.3721847
2:      え 0.33334783 0.2347444 0.2139640
3:      い 0.02101790 0.2002376 0.2246622
4:      あ 0.01777069 0.1891867 0.1891892

一見するとスマートな方法に思える。

しかし、列名の文字エンコーディングを見抜くのは割と難しい。ちょっとしたクイズを出してみよう。次の列名のうち、文字エンコーディングが異なるのは何番目だろうか？

colnames(df_all) |> stringi::stri_enc_detect() |> 
  purrr::set_names(colnames(df_all))

$スペード
      Encoding Language Confidence
1 windows-1252       es       0.42
2     UTF-16BE                0.10
3     UTF-16LE                0.10
4    Shift_JIS       ja       0.10
5      GB18030       zh       0.10
6         Big5       zh       0.10

$ハート
   Encoding Language Confidence
1  UTF-16BE                 0.1
2  UTF-16LE                 0.1
3 Shift_JIS       ja        0.1
4   GB18030       zh        0.1
5      Big5       zh        0.1

$ダイヤ
   Encoding Language Confidence
1     UTF-8                 0.8
2  UTF-16BE                 0.1
3  UTF-16LE                 0.1
4 Shift_JIS       ja        0.1
5   GB18030       zh        0.1

$クラブ
   Encoding Language Confidence
1  UTF-16BE                 0.1
2  UTF-16LE                 0.1
3 Shift_JIS       ja        0.1
4   GB18030       zh        0.1
5      Big5       zh        0.1

$ジョーカー
      Encoding Language Confidence
1 windows-1250       pl        0.5
2     UTF-16BE                 0.1
3     UTF-16LE                 0.1
4    Shift_JIS       ja        0.1
5      GB18030       zh        0.1
6       EUC-JP       ja        0.1
7       EUC-KR       ko        0.1
8         Big5       zh        0.1

正解は3番目の「ダイヤ」である。このデータは次のコードで作成した。

df_all <- iris
colnames(df_all) <- c("スペード", "ハート", "ダイヤ", "クラブ", "ジョーカー")
colnames(df_all)[3] <- iconv(colnames(df_all)[3], from = "SJIS", to = "UTF-8")
head(df_all)

  スペード ハート ダイヤ クラブ ジョーカー
1      5.1    3.5    1.4    0.2     setosa
2      4.9    3.0    1.4    0.2     setosa
3      4.7    3.2    1.3    0.2     setosa
4      4.6    3.1    1.5    0.2     setosa
5      5.0    3.6    1.4    0.2     setosa
6      5.4    3.9    1.7    0.4     setosa

このクイズが解けた人でも、現実のデータのぐちゃぐちゃな文字エンコーディングを見ると腰が引けるだろう。

また、最初に見たように、「ひとつでも異なる文字エンコーディングが混在していると、すべての列名が文字化けする」というのも、この方法での解決を困難にする。

おわりに

LighitGBM で変数重要度を出すと列名が文字化けするときの解決策を3つ紹介した。

ここまで読んでくれた人に、ハドリー・ウイッカムの次の言葉を送ろう。

Why are you using SJIS ?

https://github.com/tidyverse/dplyr/issues/339#issuecomment-38159109

参考文献

Rにおける文字列処理について詳しく書かれている（第8章）

Rによる自動データ収集: Webスクレイピングとテキストマイニングの実践ガイド

作者:Munzert,Simon,Rubba,Christian,Meissner,Peter,Nyhuis,Dominic
共立出版

Amazon

2022-12-27

Clustered ICE でサンプルごとのクラスタ割付を取得する

最近、『機械学習を解釈する技術』という本を読んでいます。

機械学習を解釈する技術〜予測力と説明力を両立する実践テクニック

作者:森下光之助
技術評論社

Amazon

非常に分かりやすくて良い本なのでおすすめです。

この本には出てこないのですが、著者のブログで紹介されている Clustered ICE という手法がとても便利です。

Clustered ICE は、データサンプルごとに得られる ICE を傾向ごとにクラスタ化してくれる手法です。この手法により、ICE の傾向が異なるサンプル群があるかどうかを簡単に調べることができます。

しかし、これを実装している ingredients パッケージの cluster_profiles() 関数は、どのサンプルがどのクラスタに属するのかという情報を返してくれません。

つまり、クラスタ化すると傾向の異なるサンプル群が発見できても、どのサンプルがどのクラスタに割り付けられたのかが分からないという状況が起きます。

そこで、この記事では、cluster_profiles() 関数から、サンプルごとのクラスタ割付を取得する方法を説明します。

Clustered ICE

データとモデルは下記のブログ記事の「シミュレーション2」をそのまま使います。

dropout009.hatenablog.com

df |> head()

       Y     X1     X2    X3
1 -5.13  -0.802  0.889     0
2 -2.23  -0.457 -0.394     1
3 -0.843  0.533  0.292     0
4  4.20  -0.551  0.934     1
5 -1.29  -0.595 -0.174     1
6  3.62  -0.393  0.790     1

これを使って explainer を作成します。

library(DALEX)
explainer <- explain(fitted, data = df |> select(-Y), y = df |> pull(Y))

サンプルをランダムに取り出して、特徴量 X2 に対する ICE を作成します。

library(ingredients)
df_ice <- df |> sample_n(100)
ice <- ceteris_paribus(explainer, variables = "X2", new_observation = df_ice)

この ICE について、 Clustered ICE を適用します。

clustered_ice <- cluster_profiles(ice, k = 2)
plot(clustered_ice)

特徴量 X2 が増加したとき、出力変数が増加するサンプル群と減少するサンプル群の2つのクラスタに分かれることが分かります。

この結果を受けて、それぞれのクラスタに属するサンプルの違いを調べたいのですが、cluster_profiles() 関数は、どのサンプルがどのクラスタに属するのかという情報を返しません。

そこで、次のようにして無理やり取得します。

クラスタ割付の取得

まずは cluster_profiles() 関数の中身を確認します。

> cluster_profiles
function (x, ..., aggregate_function = mean, variable_type = "numerical", 
    center = FALSE, k = 3, variables = NULL) 
{
    （省略）
    clus <- cutree(hclust(as.dist(dist_mat), method = "ward.D2"), 
        k = k)
    （省略）
}

クラスタ割付の情報は関数内の clus 変数に格納されていることが分かります（ついでにクラスタの手法は階層クラスタリング (hclust) であることも分かります）。

この clus 変数を無理やり取得します。それには trace() 関数を使用します。trace() 関数は、任意の関数内で任意のコードを実行できるようにする関数です。

ここでは、cluster_profiles() の終了時点での clus 変数をグローバル環境に代入します。

trace(cluster_profiles, exit = expression(clus <<- clus))
cluster_profiles(ice, k = 2)
untrace(cluster_profiles)

これにより、グローバル環境に clus という変数ができます。あとはこれを df_ice の順番に合うように並び替えて新しいカラムとして加えるだけです。

names(clus) <- str_pad(names(clus), width = 3, pad = "0")
clus <- clus[sort(names(clus))]

d <- df_ice |> cbind(clus = clus)
d |> head()

             Y         X1          X2 X3 clus
001  0.4637253  0.3238136 -0.04910185  0    1
002  1.6953002  0.9427162  0.14749948  1    2
003 -3.1748771 -0.6070075 -0.50735357  1    2
004  0.5291305 -0.8114684  0.26476546  1    2
005 -0.1732673  0.4395402 -0.12539468  1    2
006  0.7308131  0.5727512 -0.04552853  0    1

サンプルごとのクラスタ割付が分かったので、クラスタ間でサンプルがどのような違いを持つのかを調査できます。試しに、クラスタごとに各特徴量の平均を集計してみます。

d |> group_by(clus) |> summarise_at(vars(X1, X2, X3), mean)

   clus      X1      X2    X3
1     1  0.0807  0.0582     0
2     2 -0.0493 -0.0523     1

特徴量 X3 が 0 と 1 に完全に分かれていることが分かります。これにより、クラスタ 1 と 2 の違いが、特徴量 X3 による違いだという知見を得ることができました。

`model_profile()` の場合

さて、Clustered ICE は、DALEX パッケージでは model_profile() 関数に引数 k を指定することで実行できます。この場合についてもサンプルごとのクラスタ割付を取得する方法を書いておきます。

model_profile() 関数は、内部で cluster_profiles() 関数を使っているので、trace() の実行は同じコードになります。異なるのは、explainer を作るときに df_ice を使う必要があるという点です。

library(DALEX)
N <- 100
k <- 2

df_ice <- df |> sample_n(N)

explainer_ice <- explain(fitted, data = df_ice |> select(-Y), y = df_ice |> pull(Y))

trace(cluster_profiles, exit = expression(clus <<- clus))
clustered_ice <- model_profile(explainer_ice, variables = "X2", N = N, k = k)
untrace(cluster_profiles)

names(clus) <- str_pad(names(clus), width = floor(log10(N)+1), pad = "0")
clus <- clus[sort(names(clus))]

d <- df_ice |> cbind(clus = clus)
d |> head()

             Y         X1          X2 X3 clus
001  0.4637253  0.3238136 -0.04910185  0    1
002  1.6953002  0.9427162  0.14749948  1    2
003 -3.1748771 -0.6070075 -0.50735357  1    2
004  0.5291305 -0.8114684  0.26476546  1    2
005 -0.1732673  0.4395402 -0.12539468  1    2
006  0.7308131  0.5727512 -0.04552853  0    1

model_profile() でもサンプルごとのクラスタ割付を取得することができました。

参考

機械学習を解釈する技術〜予測力と説明力を両立する実践テクニック

作者:森下光之助
技術評論社

Amazon

dropout009.hatenablog.com

2022-12-07

R言語こぼれ話(2) tryCatch はなんでもつかめる

R言語を使う上で知っていても知らなくても特に問題にならないようなこぼれ話をしていくシリーズ、第2回は tryCatch() の話です。

tryCatch() は主にエラー処理に使われる関数です。例えば、ファイルを読み込もうとしたときに、そのファイルが存在しない場合は警告が発生します。警告が発生したとき、ファイルの読み込みを中断し、その警告が持つメッセージを表示するには次のようにします。

tryCatch({
  read.csv("hogehoge.csv")
}, warning = function(w) {
  print(w$message)
})

[1] "cannot open file 'hogehoge.csv': No such file or directory"

一般に、なんらかの処理中にエラーが発生した場合に、別の処理を実行するには次のように書きます。

tryCatch({
  # なんらかの処理
}, warning = function(w) {
  # 処理中に警告が発生したときの処理
}, error = function(e) {
  # 処理中にエラーが発生したときの処理
}, finally = {
  # 後処理
})

つまり、処理の実行を試してみて (try) 、警告やエラーなどが発生した場合にはそれをつかんで (catch)、別の処理を実行する、というわけですね。

この「つかむ部分」に書けるものにはどんなものがあるでしょうか？ Rがデフォルトで提供しているのは次の5つです。

error（エラー）
warning（警告）
message（メッセージ）
interrupt（割り込み）
condition（上記すべてをキャッチできる特別なやつ）

ふつうのRユーザーであれば、error, warning, message の3つと、後処理を書くために finally を知っていれば十分だと思います。

しかし、この「つかむ部分」、実は「なんでもつかめる」ことはあまり知られていないように思います。

`tryCatch()` はなんでもつかめる

tryCatch() は実はなんでもつかめます。

tryCatch() に自分の好きなものをつかませるためには、カスタムエラーを作成します。ここでは例として hoge をつかませるためのカスタムエラーを作成しましょう。カスタムエラーの作成方法は、message を持つ list オブジェクトを作成し、hoge, error, condition をクラスとして設定するだけです。

# カスタムエラーを作成する
my_error <- list(message = "ほげほげ")
class(my_error) <- c("hoge", "error", "condition")

エラーを発生させる関数 stop() を使って、このカスタムエラーを発生させると、tryCatch() で hoge をつかむことができます。

tryCatch({
  stop(my_error)
}, hoge = function(e) {
  "hogeをキャッチしたよ"
})

[1] "hogeをキャッチしたよ"

このように、hoge でも fuga でも piyo でも、なんでも思い通りに tryCatch() につかませることができるというわけです。

ちなみに、rlang パッケージにはカスタムエラーを作成する便利な関数があります。次のコードは上記と同じ動きをします。

my_error <- rlang::error_cnd("hoge", message = "ほげほげ")

tryCatch({
  rlang::cnd_signal(my_error)
}, hoge = function(e) {
  "hogeをキャッチしたよ"
})

何に使えるのか？

tryCatch() はなんでも好きなものをつかめることが分かりました。では、これが何の役に立つのか？というのは難しいところです。なかなか使いどころが思いつきません。けっこうな無駄知識だと思います。

ただ、1つ思いつくのは、エラー処理の細分化に使えそうです。

例えば、ファイルが存在しない場合とそれ以外の場合で、エラー処理を変えたいとします。 Rのエラーはメッセージしか持たないので、メッセージの内容によってエラー処理を分けます。

file_name <- "piyopiyo.txt"
tryCatch({
  if (!file.exists(file_name)) {
    stop("File not found")
  } else {
    message("ファイル見つけてえらい")
  }
}, error = function(e) {
  if (e$message == "File not found") {
    message(file_name, "は見つかりませんでした")
  } else {
    message(e)
  }
})

これだと、メッセージの文言を変えるたびにエラー処理のコードも変える必要があります。

そこで、ファイルが存在しない場合に対応したカスタムエラーを作成すると次のように書けます。

fnf_error <- rlang::error_cnd("file_not_found", message = "File not found")

file_name <- "piyopiyo.txt"
tryCatch({
  if (!file.exists(file_name)) {
    rlang::cnd_signal(fnf_error)
  } else {
    message("ファイル見つけてえらい")
  }
}, file_not_found = function(e) {
  message(file_name, "は見つかりませんでした")
}, error = function(e) {
  message(e)
})