対数変換して求めた予測値を指数で再変換するときの下振れの補正（スミアリング）

　昨日学生に説明したことのメモ。
　世の中には、何かと何かの「掛け算」で決まっていそうなものが色々あります。たとえば、私はオフロードバイクで山を走るのが好きなのですが、オフロードの走破能力は、マシンの性能と自分の技量の「掛け算」で高まるのだという実感があります。足し算ではなく相乗効果があるということですね。

　走破能力 = マシンの性能 × 自分の技量

　他にも、経済学に出てくるコブ=ダグラス型生産関数では、技術（生産性）×資本×労働力で供給力が決まるというふうに定式化されています。

$Y = A \cdot K^{\alpha} L^{\beta}$

$Y$ ：生産力
$A$ ：全要素生産性（TFP） ※これはパラメータとして推定される
$K$ ：資本
$L$ ：労働力

　また、交通や物流の分野では、地域間のフローをモデル化するのに重力モデルというのが使われたりします。天体の間に働く重力は、それぞれの天体が重ければ重いほど強く、また距離が離れるにつれて弱くなりますが、そのアナロジーで、2つの都市の人口が大きければ大きいほどその間の交易が盛んになり、距離が遠くなると交易が減るということです。ニュートン力学だと分母は距離の2乗ですが、何乗するかはパラメータとして推定します。

$T_{ij} = G \cdot \frac{M_i^{\alpha} M_j^{\beta}}{d_{ij}^{\gamma}}$

$T_{ij}$ ：都市iから都市jへの交易フロー（人流や物流など）
$G$ ：重力定数 ※これはパラメータとして推定される
$M_i$ , $M_j$ ：都市iと都市jの規模（人口や経済力など）
$d_{ij}$ ：2都市間の距離

　こういう乗法的な関係をモデリングするとき、伝統的には、両辺の対数をとって線形結合の形に変換します。データの各変数を対数化しておけば、そのまま線形回帰でパラメータを推定できるからです。

$\log T_{ij} = \log G + \alpha \log M_i + \beta \log M_j - \gamma \log d_{ij}$

　ところで、こういう形でOLSにより推定されたモデルを使って $Y$ や $T$ の予測値を出すのであれば、まず $\log Y$ や $\log T$ の値を出してからそれを $\exp(\log T)$ というふうに指数変換して元の真数スケールに戻すわけですが、これだけだと、予測値が系統的に下振れしてしまう再変換バイアスが生じますので、何らかの補正が必要です。
　なんでかというと、まずイエンセンの不等式というのがあって、下に凸な関数 $f$ に関して

$E[f(X)] \geq f(E[X])$

もしくは分かりやすくいうと

$\mathrm{mean}(f(X)) \geq f(\mathrm{mean}(X))$

ということですが、これは要するに、ある変数群の平均をとってから関数を適用すると、先に関数を適用してから平均をとるのに比べて、小さくなるということです（等しい場合もある）。
　さっきの対数線形回帰においては、誤差項の扱いでこれが問題になります。要は、

$\log T_{ij} = \log G + \alpha \log M_i + \beta \log M_j - \gamma \log d_{ij} + \varepsilon_{ij}$

　という形でOLS回帰をすると、「誤差の平均をゼロ」にするような線を引くことになるので（線ではないがｗ）、推定されたパラメータを使って「対数スケールでの予測値」を求める段階では、誤差項はどうせ平均的にはゼロになるものとして無視できます。ところが、対数化した世界で誤差項の平均 $\mathrm{mean}(\varepsilon)$ がゼロであっても、それを $\exp()$ した真数のスケールでは $\mathrm{mean}(\exp(\varepsilon))$ が1（対数スケールでの0に対応）より大きくなってしまうので、予測値をつくるときに誤差項の影響を復活させる必要があります。
　具体的には、対数線形モデルを $\exp()$ すると乗法モデルの形になるので、 $\mathrm{mean}(\exp(\varepsilon))$ を補正係数として予測値にかければよい。これにより、予測値が上方修正されます。これがスミアリング法です。
　ちなみに、GLMで「logリンク+ガンマ分布」を使ったりするのであれば、平均構造を直接モデル化するので、こういう下振れは生じません。

　以下、試しに、3つの説明変数からなる乗法モデルで乱数データセットを生成し、

単なる線形回帰でパラメータ推定して予測値を求める
対数変換してから線形回帰で推定して予測値を指数変換で戻す
対数変換してから線形回帰で推定して予測値を指数変換で戻してスミアリング補正
GLMのlogリンク+ガンマ分布（対数正規分布が直接指定できないので）で推定して予測値を求める

の4パターンでどれぐらいズレるかを試してみました。
　結果はこうなりました。RMSEは誤差を絶対値で評価するもの、sMAPEは誤差を割合で評価するもので、いずれも小さければ小さいほどフィットが良いです。Biasは、正負の方向つきで、モデル予測値が実測値（乱数）から平均的にどれぐらいズレてるかを表すものです。

> print(summary_table)
                   Model Mean_RMSE Mean_sMAPE Mean_Bias
1              naive OLS   43.2946     0.7949   -0.0530
2 Log-linear + naive exp   29.0960     0.3118   -2.8425
3  Log-linear + smearing   28.6847     0.3177    0.0006
4         GLM Gamma(log)   28.6884     0.3177   -0.0048

　まず、単なるOLS(1)だとRMSEもsMAPEも一番大きくなっていて、要するに当てはまりが一番悪い。指数的なモデルを線形でモデリングしてるのだから、当然です。ただ、バイアスは小さくて、これは真数のスケールで誤差の平均がゼロになるような推定をしてるからでしょう。
　対数線形モデル+スミアリング補正(3)と、GLM(4)については、RMSEもsMAPEも小さくて、バイアスも小さいです。微妙にGLMのほうがバイアスが大きいのは、対数正規分布ではなくガンマ分布を使ってるからでしょうかね。
　で、ポイントは、対数線形回帰の予測値を単純に指数変換で戻しただけ(2)だと、バイアスが負で大きくなっていて、要するにめちゃめちゃ下振れしています。
　GLMが手っ取り早くていいと思いますが、線形回帰モデルのほうが分かりやすいからか、私がいる分野（交通関係）ではよく使われている印象で、予測値をつかってシミュレーションをする（かつその水準が問題になる）のでなければ、補正の話は必要ないので出てこないです。そして、出てこないので、下振れのリスクを意識してなくて計算をミスったりします（笑）

（参考：このブログのおすすめ記事一覧はコチラ）

　乱数実験のRのコードを置いておきます。ChatGPTに書いてもらったものを手動で調整しました。

# ============================================================
# 乱数実験：
# 真のデータ生成過程を「弾力性モデル」にする
#
# 乗法（弾力性）モデル：
# y = exp(beta0) * x1^beta1 * x2^beta2 * x3^beta3 * u
# ただし u = exp(eps), eps ~ N(0, sigma_true^2)（対数正規誤差）
#
# その上で
# 1) 線形回帰（y ~ x1 + x2 + x3）
# 2) 対数線形 + 単なるexp（補正なし）
# 3) 対数線形 + exp + スミアリング補正
# 4) GLM (logリンク、Gamma分布)
# を比較する
#
# 出力する指標：
# RMSE, sMAPE, Bias（ズレ）
# ============================================================

set.seed(123)

# -------------------------------
# シミュレーション設定
# -------------------------------

n_sim <- 1000     # シミュレーション回数
n_train <- 1000   # 学習データ数
n_test  <- 10000  # テストデータ数

# 真のパラメータ
beta0 <- 1.0
beta1 <- 2.0
beta2 <- 0.8
beta3 <- -0.6

# 対数スケールでの誤差の標準偏差
# 大きいほど元のスケールでのばらつきが大きくなる
sigma_true <- 0.4

# 結果保存用
results <- data.frame(
  sim = 1:n_sim,
  
  RMSE_linear     = NA,
  RMSE_log_naive  = NA,
  RMSE_log_smear  = NA,
  RMSE_glm        = NA,
  
  sMAPE_linear    = NA,
  sMAPE_log_naive = NA,
  sMAPE_log_smear = NA,
  sMAPE_glm       = NA,
  
  Bias_linear     = NA,
  Bias_log_naive  = NA,
  Bias_log_smear  = NA,
  Bias_glm        = NA
)

# ============================================================
# シミュレーション開始
# ============================================================

for (s in 1:n_sim) {
  
  # -------------------------------
  # 1. 学習データを作る
  # -------------------------------
  # 説明変数は正の値が必要なのでexpでつくる
  x1_train <- exp(rnorm(n_train, mean = 1.0, sd = 0.5))
  x2_train <- exp(rnorm(n_train, mean = 0.5, sd = 0.4))
  x3_train <- exp(rnorm(n_train, mean = 0.8, sd = 0.6))
  
  # 真のモデル（弾力性モデル）
  eps_train <- rnorm(n_train, 0, sigma_true)
  
  y_train <- exp(beta0) *
    (x1_train ^ beta1) *
    (x2_train ^ beta2) *
    (x3_train ^ beta3) *
    exp(eps_train)
  
  train_dat <- data.frame(
    y  = y_train,
    x1 = x1_train,
    x2 = x2_train,
    x3 = x3_train
  )
  
  # -------------------------------
  # 2. テストデータを作る
  # -------------------------------
  x1_test <- exp(rnorm(n_test, mean = 1.0, sd = 0.5))
  x2_test <- exp(rnorm(n_test, mean = 0.5, sd = 0.4))
  x3_test <- exp(rnorm(n_test, mean = 0.8, sd = 0.6))
  
  eps_test <- rnorm(n_test, 0, sigma_true)
  
  y_test <- exp(beta0) *
    (x1_test ^ beta1) *
    (x2_test ^ beta2) *
    (x3_test ^ beta3) *
    exp(eps_test)
  
  test_dat <- data.frame(
    y  = y_test,
    x1 = x1_test,
    x2 = x2_test,
    x3 = x3_test
  )
  
  # ============================================================
  # 3. モデル1：何も考えない線形回帰
  #    y ~ x1 + x2 + x3
  # ============================================================
  
  fit_linear <- lm(y ~ x1 + x2 + x3, data = train_dat)
  
  pred_y_linear <- predict(fit_linear, newdata = test_dat)
  
  # ============================================================
  # 4. モデル2：対数線形モデル + naive exp
  #    log(y) ~ log(x1) + log(x2) + log(x3)
  # ============================================================
  
  fit_log <- lm(log(y) ~ log(x1) + log(x2) + log(x3), data = train_dat)
  
  pred_log_test <- predict(fit_log, newdata = test_dat)
  
  # 補正なしでそのまま exp で戻す
  pred_y_naive <- exp(pred_log_test)
  
  # ============================================================
  # 5. モデル3：対数線形モデル + smearing 補正
  # ============================================================
  
  resid_log <- residuals(fit_log)
  
  # スミアリング係数
  smear_factor <- mean(exp(resid_log))
  
  pred_y_smear <- exp(pred_log_test) * smear_factor
  
  # ============================================================
  # 6. モデル4：GLM (Gamma, log link)
  #    y ~ log(x1) + log(x2) + log(x3)
  # ============================================================
  
  fit_glm <- glm(
    y ~ log(x1) + log(x2) + log(x3),
    family = Gamma(link = "log"),
    data = train_dat
  )
  
  pred_y_glm <- predict(fit_glm, newdata = test_dat, type = "response")
  
  # ============================================================
  # 7. 予測精度を計算
  # ============================================================
  
  # ---- RMSE ----
  rmse_linear <- sqrt(mean((y_test - pred_y_linear)^2))
  rmse_naive  <- sqrt(mean((y_test - pred_y_naive)^2))
  rmse_smear  <- sqrt(mean((y_test - pred_y_smear)^2))
  rmse_glm    <- sqrt(mean((y_test - pred_y_glm)^2))
  
  # ---- sMAPE ----
  # 0〜2の範囲
  smape_linear <- mean(abs(y_test - pred_y_linear) / ((abs(y_test) + abs(pred_y_linear)) / 2))
  smape_naive  <- mean(abs(y_test - pred_y_naive)  / ((abs(y_test) + abs(pred_y_naive))  / 2))
  smape_smear  <- mean(abs(y_test - pred_y_smear)  / ((abs(y_test) + abs(pred_y_smear))  / 2))
  smape_glm    <- mean(abs(y_test - pred_y_glm)    / ((abs(y_test) + abs(pred_y_glm))    / 2))
  
  # ---- Bias ----
  # 平均予測誤差：予測値 - 実測値
  bias_linear <- mean(pred_y_linear - y_test)
  bias_naive  <- mean(pred_y_naive  - y_test)
  bias_smear  <- mean(pred_y_smear  - y_test)
  bias_glm    <- mean(pred_y_glm    - y_test)
  
  # ============================================================
  # 8. 保存
  # ============================================================
  
  results$RMSE_linear[s]     <- rmse_linear
  results$RMSE_log_naive[s]  <- rmse_naive
  results$RMSE_log_smear[s]  <- rmse_smear
  results$RMSE_glm[s]        <- rmse_glm
  
  results$sMAPE_linear[s]    <- smape_linear
  results$sMAPE_log_naive[s] <- smape_naive
  results$sMAPE_log_smear[s] <- smape_smear
  results$sMAPE_glm[s]       <- smape_glm
  
  results$Bias_linear[s]     <- bias_linear
  results$Bias_log_naive[s]  <- bias_naive
  results$Bias_log_smear[s]  <- bias_smear
  results$Bias_glm[s]        <- bias_glm
}

# ============================================================
# 9. 平均結果を表にまとめる
# ============================================================

summary_table <- data.frame(
  Model = c(
    "naive OLS",
    "Log-linear + naive exp",
    "Log-linear + smearing",
    "GLM Gamma(log)"
  ),
  Mean_RMSE = c(
    mean(results$RMSE_linear),
    mean(results$RMSE_log_naive),
    mean(results$RMSE_log_smear),
    mean(results$RMSE_glm)
  ),
  Mean_sMAPE = c(
    mean(results$sMAPE_linear),
    mean(results$sMAPE_log_naive),
    mean(results$sMAPE_log_smear),
    mean(results$sMAPE_glm)
  ),
  Mean_Bias = c(
    mean(results$Bias_linear),
    mean(results$Bias_log_naive),
    mean(results$Bias_log_smear),
    mean(results$Bias_glm)
  )
)

summary_table[, -1] <- round(summary_table[, -1], 4)

print(summary_table)