Reports

Twitter: @mtknnktm.

未知の非線形な交互作用をマルチレベルモデリングでモデル化してみる

はじめに

なんか変なタイトルですが、この記事では次のような現象について考えます。

y \sim Normal(E(x_1, x_2), \sigma)
E(x_1, x_2) = \beta_0 + \beta_1 x_1 + \beta_2 (x_1) x_2

これはぱっと見は通常の重回帰モデルですが、よく見るとx_2の回帰係数\beta_2はx_1の関数になっています。本記事ではこういったx_1, x_2 間の未知の（非線形かもしれない）交互作用について考えることを目的とします。

どういった状況でこんなよくわからないモデルを考えないといけないのでしょうか？
例えば、とあるシステムの開発プロジェクトの開発にかかる工数を分析するモデルだとして、目的変数 y が開発期間、説明変数 x_1 が開発規模（人数）、x_2が開発者のコミュ力だとしましょう（※ 単なる例なので現象の妥当性は気にしてはいけません）。

小規模（x_1が小）ならばコミュニケーションを取る相手が少なく文脈も常に共有しているため、コミュ力x_2の影響は小さいと考えられます。つまり、\beta_2(x_1)は小さくなります。
一方で規模が大きくなる（x_1が大きくなる）ほどコミュニケーションを取る相手が多く、かつ、文脈も共有されていないことが増えるため、コミュ力x_2の影響は非常に大きくなると考えられます。つまり、\beta_2(x_1)は非常に大きいと考えられます。

従って、\beta_2(x_1)は線形関数ではなく、例えば\beta_2(x_1)=\exp(a x_1)みたいな感じでx_1について指数的に増大したりしそうです。また、N次の多項式で表される極大値や極小値があるような複雑な形をしているかもしれません。ここでの目的はコミュ力x_2の開発期間への影響が規模によってどのように変わってくるか？ （つまり\beta_2(x_1)の形）について知ることです。

モデル

では、上記の\beta_2(x_1)が未知の関数だとして、どうにかして\beta_2(x_1)の形を推定する方法を考えます。
この記事では\beta_2(x_1)をx_1で適当に区切って、その区間内を線形関数で近似することを考えます。
すなわち、
E(x_1, x_2, c) = \beta_0 + \beta_1 x_1 + (\beta_{21c} (x_1 - m_c) + \beta_{20c}) x_2
という形で近似することを考えます。ここでcはx_1が所属する区間、m_cはその区間の中央の値を表します。この場合、\beta_{20c}, \beta_{21c}は区間数L個分だけ推定する必要があります。x_1で階層分けをしたマルチレベルモデリングと言ってもいいと思います。

例えば、\beta_2(x_1) = (x_1 - 2.5) ^2で、区間数L=4のときは[x_2]の回帰係数である関数\beta_2は次のように近似されるイメージです。
f:id:swarm_of_trials:20150320005934p:plain:h200
二次曲線が4つの線形関数で近似されていることがわかります。

実験

それでは実際に未知の非線形関数beta_2(x_1)を線形関数で近似できるかやってみましょう。
ここでは\beta_2(x_1) = (x_1 - 2.5) ^2/5、x_1 \in [0, 5]とし、\beta_0=1, \beta_1=2, \sigma=0.1とします。

以下にこのモデルを表すStanのコードを示します。

data {
  int<lower=0> N;   //データの数
  int<lower=0> C;   //区間の数
  real<lower=0> m; //x_1の最大値（5）

  int c[N];         //データが所属する区間id（1〜C）
  vector[N] x1; //説明変数x_1
  vector[N] x2; //説明変数x_2
  vector[N] y;   //目的変数
}
parameters {
  real beta0;
  real beta1;
  vector[C] beta2_1;
  vector[C] beta2_0;
  real s;
}
model {
  for(i in 1:N)
    y[i]~normal(beta0 + beta1 * x1[i] + (beta2_1[c[i]] * (x1[i] - m*(c[i]-0.5)/C) + beta2_0[c[i]]) * x2[i], s);  
}
generated quantities{
  vector[N] log_lik;  //対数尤度。WAICを求めるため
  vector[N] pred;     //予測値
  
  for (i in 1:N){
    log_lik[i] <- normal_log(y[i], beta0 + beta1 * x1[i] + (beta2_1[c[i]] * (x1[i] - m*(c[i]-0.5)/C) + beta2_0[c[i]]) * x2[i], s);
    pred[i] <- beta0 + beta1 * x1[i] + (beta2_1[c[i]] * (x1[i] - m*(c[i]-0.5)/C) + beta2_0[c[i]]) * x2[i];
  }
}

以下にテスト用のデータを生成して上記のStanを実行するRのコードを示します。ここでは区間数を5としています。※ 2015/03/23 使用ライブラリが書いてなかったので追記しました。

library(ggplot2)
library(dplyr)
library(rstan)
library(matrixStats)

#パラメータ
N <- 1000
beta1 <- 2
beta0 <- 1

#データ生成
d <- tbl_df(data.frame(x1=runif(N, 0, 5), x2=runif(N, 1, 5)))
beta2 <- (d$x1 - 5/2)^2 / 5
d <- d %>% mutate(y=rnorm(N, beta0 + beta1 * x1 + beta2 * x2, 0.1))

#区間作成
L <- 5
d <- d %>% mutate(cls=ceiling(x1 / (5/L)))

dat <- list(y=d$y, x1=d$x1, x2=d$x2, c=d$cls, N=N, C=L, m=5)

#MCMC
fit <- stan(file='ml.stan', data=dat, iter=3000)

上記のコードを実行して、結果を見てみましょう。Rhatは1.00付近でいい感じだったのでMCMCは収束したと判断します。

まず、モデルが正しい結果を推定できているか否か？ について見てます。以下を実行することで予測値（pred）と実際の値（y）の散布図を書きます。

d %>% mutate(pred=colMeans(extract(fit, 'pred')$pred)) %>% ggplot(aes(y, pred)) + geom_point()

以下が、その結果です。モデルはフィットしてそうであることがわかります。
f:id:swarm_of_trials:20150320115104p:plain

それでは、真の関数\beta_2(x_1) = (x_1 - 2.5) ^2/5と求めた\beta_{21c} (x_1- m_c) + \beta_{20c}を比較してみましょう。\beta_{21c}, \beta_{20c}は各区間ごとに中央値で評価しています。
f:id:swarm_of_trials:20150320115435p:plain
黒い線が真の関数、カラフルな直線が各区間で求めた線形関数です。二次関数が線形関数である程度近似されていることがわかります。

それでは区間数Lを変えてどんな感じになるか見てみましょう。
以下に区間数L=1, 2, 3, 4, 5, 10, 15, 20, 100としたときの(x_1, \beta_{21c} (x_1- m_c) + \beta_{20c})を示します。
L=1では\beta_2(x_1)の特徴が全く表現できていないこと、Lが大きくなるに伴い、\beta_2(x_1)によくフィットしたことがわかります。ただし、L=100のように区間数が大きくなりすぎると過学習が発生してしまうこともわかります。
f:id:swarm_of_trials:20150320134512p:plain

以下に区間数LとそれぞれのモデルのWAICの値を示します。Lは1〜20と100について求めました。L=1から6までは急激に減少し、その後は似たりよったりで、L=100になるとL=6〜20よりも値が大きくなりました。ここから、L\in[6, 100)の間にWAICを最小にする値があることがわかります。そのため、WAICによってLを適切に決めることでモデルを決めることができそうです。
ただし、ここでのWAICが支持するモデルは必ずしも解釈しやすいとは限らないので、結果を見て適当にLを決めてしまってもいいんじゃないかなと思います（今回の件で言えばL=5ぐらいでも十分に傾向を見て取ることができます）。

f:id:swarm_of_trials:20150320134601p:plain:h250

ちなみにWAICは以下のようにして求める事ができます（多分合ってるはず(´・ω・`)）。

log_lik <- extract(fit, 'log_lik')$log_lik
lppd <- sum(log(colMeans(exp(log_lik))))
p_waic <- sum(colVars(log_lik))
(waic.pp <- -2*lppd + 2*p_waic)

まとめ

ちょっと無理矢理感がありますがやってみたらできました。

データサイエンティスト養成読本 [ビッグデータ時代のビジネスを支えるデータ分析力が身につく! ] (Software Design plus)

データサイエンティスト養成読本 [ビッグデータ時代のビジネスを支えるデータ分析力が身につく! ] (Software Design plus)

ベイズ統計モデリング (統計ライブラリー)

ベイズ統計モデリング (統計ライブラリー)