乗法モデル

単純なGLMにおいて、目的変数・説明変数（もしくは片方のみ）を対数変換したほうがより解釈しやすいケースがあります。
一般的に、目的変数・説明変数ともに対数変換したモデルは乗法モデル（弾力化モデル）と言われるようです。

（参考）対数変換を行う意味について。回帰分析において対数変換する背景にある前提とは？｜アタリマエ！

計量経済学分野では比較的ポピュラーな方法ですが、
使おうと思ったときに対数変換によってパラメータの解釈が代わる証明が不安だったので備忘録的に書きます。

なお、ここでは目的変数のみを対数変換する log-level-model になります。（※log-logがいいのか、level-logがいいのかはモデルの当てはまりや課題に応じて決めたほうが良さそうです）
log-level-modelのパラメータ解釈は、説明変数１単位の増加による目的変数の増加率（%）になるはずです。

（参考） Data Science Simplified Part 7: Log-Log Regression Models | by Pradeep Menon | Towards Data Science

証明

$$ \begin{eqnarray} \log y &=& \alpha + \beta x + \epsilon \ \end{eqnarray} $$

まず両辺を微分し、

$$ \begin{eqnarray} \dfrac{d}{dx} \log y &=& \dfrac{d}{dx} (\alpha + \beta x + \epsilon) \\ \beta &=& \dfrac{d \log y}{dx} \end{eqnarray} $$

ここで、合成関数の微分公式より

$$ \begin{eqnarray} \dfrac{d \log y}{dx} &=& \dfrac{d \log y}{dy} \cdot \dfrac{dy}{dx} \tag{1} \end{eqnarray} $$

$(1)$右辺の $d \log y / dy$ は対数関数の微分公式より

$$ \begin{eqnarray} \dfrac{d \log y}{dy} &=& (\log y)' &=& \dfrac{1}{y} \end{eqnarray} $$

したがって $(1)$ は以下のとおり

$$ \begin{eqnarray} \dfrac{d \log y}{dy} \cdot \dfrac{dy}{dx} &=& \dfrac{1}{y} \cdot \dfrac{dy}{dx} \\ &=& \dfrac{dy/y}{dx} \end{eqnarray} $$

よって、 $$ \begin{eqnarray} \beta &=& \dfrac{dy/y}{dx} \end{eqnarray} $$

$dy/y$ は $y$ の増加率
$dx$ は xの増加量

以上により、パラメータ$\beta$ は$x$の増加量に対する$y$の増加率であることが証明できました。
なお、log-log-modelも同じ手順で証明できます。

無事、log-level-modelのパラメータの解釈が理解できたところで、実データで解釈性を見ていきます。

利用データ

seabornのtipsデータセットを利用します。参考

df =  sns.load_dataset("tips")

f:id:kopaprin:20220201005511p:plain

支払額や性別、人数などに応じて、チップがいくら払われたかのデータセットのようです。
これを題材にしていきます。

なお、モデリングには以下のようにダミーデータ化し、簡単のためdayは除外しています。 f:id:kopaprin:20220201010105p:plain

データの概観

目的変数とするtipは以下のような分布となっており、やや右すそが長い対数正規分布に近い形です。 f:id:kopaprin:20220201005925p:plain

データ上0以下にはならないため、ここだけ見ても通常の正規分布を仮定した一般線形モデルはあまり適さなそうであることが分かります。
が、今回の趣旨はパラメータ解釈の違いを見たいのでこの辺は論点にしません。

level-level-model

まずはlevel-level-model、つまり単純な一般線形モデルでのモデリング結果です。

level_level_model = """

data {
  int N ; // データ数
  int C ; //説明変数カラム数
  matrix[N,C] X ; //説明変数
  vector<lower=0>[N] Y ;  // tip 目的変数
}

parameters{
  vector[C] b ; //パラメータ
  real<lower=0> s_Y ; //標準誤差
}

model {
  Y ~ normal(X*b, s_Y) ;
}

generated quantities {
  real y_pred[N];
  y_pred = normal_rng(X*b, s_Y) ;
}

"""

log-level-model

次いで、log-level-model

log_level_model = """

data {
  int N ; // データ数
  int C ; //説明変数カラム数
  matrix[N,C] X ; //説明変数
  vector<lower=0>[N] Y ;  // tip 目的変数
}

parameters{
  vector[C] b ; //パラメータ
  real<lower=0> s_Y ; //標準誤差
}

model {
  log(Y) ~ normal(X*b, s_Y) ;
}

generated quantities {
  real y_pred[N];
  y_pred = normal_rng(X*b, s_Y) ;
}

"""

本来はWAICやWBICで評価したり、バリデーションデータセットで評価すべきなのでしょうが、こちらも簡単のため（~~めんどくさい~~）とりあえず目的変数Yと事後予測分布のEAP推定量のMAEで評価してみました。

print("MAE")
print("level-level-model : ",mean_absolute_error(df_fix["tip"],fit_level_level.extract()["y_pred"].mean(axis=0)))
print("log-level-model : ",mean_absolute_error(df_fix["tip"],np.exp(fit_log_level.extract()["y_pred"].mean(axis=0))))

MAE
level-level-model : 0.736165905967111
log-level-model : 0.7633751430630242

うーん、log-level-modelのほうが悪いですね（爆）
この辺は自分のモデリングや評価が適当だったりすることにも起因する気がするので、あまり気にしないようにしましょう（白目）

パラメータ評価

level-level-modelとlog-level-modelのパラメータのEAP推定量は以下のとおりです。

f:id:kopaprin:20220201012233p:plain

例えばsexを見ると、（ここではFemale = 1, Male = 0）
女性の場合level-level-modelでは0.026ドル(?) 男性より多くチップを払うことが期待されます。一方log-level-modelでは、1.3%多く払うことが期待できます。

解釈の違いは上記の言葉のとおりですが、モデルの前提としては大きく異なります。
level-level-modelでは「どんなに高額な食事をして、大人数でも女性なら0.026ドル多く払うだけ」と考えています。
一方、log-level-modelは「総支払額が高ければ高いほど、女性の場合チップ支払い額が1.3%多くなる」と考えています。

なんとなくですが、10ドルの食事をしたケースと100ドルの食事をしたケースでは、チップの支払額も相対的に高くなりそうです。
このように、変数の影響の解釈とモデルの考え方が、対数変換するだけでかなり変わってくることがわかります。
うまく取り扱えば、マーケティングのための顧客理解や現象理解に役立てることができると思います。

おまけ

チップ支払い額には個人差が大きく関与しそうな気がします。
(これはダミーデータかもしれませんが・・・) 個人差 $mu$ をモデルに入れた階層モデルでもやってみました。

multi_log_level_model = """

data {
  int N ; // データ数
  int C ; //説明変数カラム数
  matrix[N,C] X ; //説明変数
  vector<lower=0>[N] Y ;  // tip 目的変数
}

parameters{
  vector[C] b ; //パラメータ
  vector[N] mu ; //個人差
  real<lower=0> s_Y ; //標準誤差
  real<lower=0> s_mu ; //個人差誤差
}

model {
  mu ~ normal(0, s_mu) ;
  log(Y) ~ normal(X*b + mu, s_Y) ;
}

generated quantities {
  real y_pred[N];
  y_pred = normal_rng(X*b + mu, s_Y) ;
}

"""