中古マンション売買データを分析してみた(3.2) - 駅が複数の路線に所属する階層モデル
引き続き，階層モデルを実際のデータに即した形で組んでみようという試みを続けていきたいと思います．昨日の記事では，1つの駅は1つの路線にしか所属しない*1という制約を持ったモデルを作成しました．で，今日の試みは，1つの駅が複数の路線に所属できるようにモデルを拡張しましょうか，というお話です．
1つの子が複数の親に所属可能な階層モデル
さっそくですが，今回作成したのは以下のようなモデルです．1行目は昨日のものと同じですが，2行目の後半部分が総和記号を使った形になっています．これは駅の効果が固定項と，当該駅が所属する駅全てのランダム効果の総和になっているという意味です．この形をとることにより，1つの駅が複数の路線に所属することが可能になります（とりあえずここでは，路線がたくさん通っているほど，駅の魅力は線形和で増加するというシンプルなモデルを想定しています*2）．
モデルの記述とシミュレーション
モデル式
ということで，Stanでこれをモデリングしてみましょう．Stanコードはこちらです．前回との違いは，dataブロックの8行目がstation-train matrixとなっているように，路線とその所属駅の関係を1/0の行列で表したものになっている点です．それに合わせてparameterブロックの5行目がvector[N_T]に，modelブロックの11行目の回帰モデルの内部が行列演算の形になっています．
data { int<lower=1> N; # sample num int<lower=1> M; # independents' num int<lower=1> N_T; # train num int<lower=1> N_S; # station num matrix[N, M] X; # independents vector[N] Y; # dependent matrix[N_S, N_T] ST; # station-train matrix int<lower=1, upper=N_S> S[N]; # station } parameters { real a; vector[M] b; real as; vector[N_T] r_t; real r_s[N_S]; real<lower=0> s; real<lower=0> s_rs; real<lower=0> s_rt; } model { # regresion model with random effect for (i in 1:N) Y[i] ~ normal(a+X[i]*b+r_s[S[i]], s); # prior distributions s ~ uniform(0, 1.0e+4); a ~ normal(39, 1.0e+4); for (i in 1:M) b[i] ~ normal(0, 1.0e+4); for (i in 1:N_S) r_s[i] ~ normal(as+ST[i]*r_t, s_rs); # hierarchical prior distribution s_rs ~ uniform(0, 1.0e+4); as ~ normal(0, 1.0e+4); for (i in 1:N_T) r_t[i] ~ normal(0, s_rt); # 2 hierarchical prior distibution s_rt ~ uniform(0, 1.0e+4); }
そしてこれをキックする側のRコードは，以下のようになります*3．
# load data d = read.delim('data/mantions.csv', header=T, sep=',') d = na.omit(d) attach(d) st = read.delim('data/station_train.csv', header=F, sep=',') # package data for stan X = t(rbind(distance, from, room, space)) ST = st S = station Y = price d.stan = list(N=nrow(X), N_T=length(unique(train)), N_S=length(unique(station)), M=ncol(X), X=X, ST=ST, S=S, Y=Y)
5行目で駅-路線行列のstを読み込んでますが，データの中身は以下のような感じです．行が駅，列が路線で駅が路線に所属していれば1に，そうでなければ0になる行列ですね．これはPythonで前処理をして別に作っています．
前処理のPythonスクリプトはこちらです．メインのモデル作成に使っているデータフレームに対して，軽く前処理して行列を作成しています．ところどころ数字がハードコーディングしてあるのは，使い捨てスクリプトだからなので，大目にみてください．
#!/usr/bin/env python #-*-coding:utf-8-*- import argparse import numpy as np # args parser = argparse.ArgumentParser() parser.add_argument('in_path') parser.add_argument('out_path') if __name__ == '__main__': # setup args = parser.parse_args() in_file = open(args.in_path, 'r') # create station-train matrix st = np.zeros((133, 23)) is_header = True for line in in_file: if is_header is True: is_header = False continue tokens = line.split(',') station, train = int(tokens[11])-1, int(tokens[12])-1 st[station][train] = 1 # write out np.savetxt(args.out_path, st, fmt="%.0f", delimiter=',') # tear down in_file.close()
シミュレーションの結果
結果からいうと，これでもやはりダメでした...推定値のサマリについて，抜粋を以下に載せますが，やはり定数項と，駅の項が基本的に収束しない感じですね．まぁそれでも前回よりはだいぶマシなわけですが...
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat a -1549.051 68.070 137.039 -1684.644 -1636.850 -1600.390 -1539.373 -1150.424 4 2.295 b[1] -4.784 0.001 0.098 -4.977 -4.848 -4.784 -4.718 -4.588 5001 1.000 b[2] 0.850 0.000 0.011 0.829 0.843 0.850 0.857 0.871 4491 1.000 b[3] -0.422 0.003 0.190 -0.782 -0.554 -0.421 -0.292 -0.047 4878 1.002 b[4] 1.935 0.002 0.149 1.642 1.835 1.936 2.032 2.232 4580 1.001 as -109.097 67.803 135.354 -510.361 -108.223 -59.499 -22.132 11.961 4 2.421 r_t[1] 2.841 0.062 4.104 -5.014 0.006 2.818 5.660 10.866 4438 1.000 r_t[2] -0.770 0.058 3.782 -8.394 -3.210 -0.696 1.698 6.654 4208 1.000 r_t[3] 2.725 0.027 1.692 -0.631 1.602 2.733 3.893 5.989 3803 1.000 ... r_t[21] -4.143 0.043 2.913 -9.847 -6.072 -4.144 -2.154 1.537 4646 1.000 r_t[22] -4.474 0.037 2.474 -9.423 -6.114 -4.483 -2.814 0.415 4384 1.000 r_t[23] 2.929 0.032 2.022 -0.984 1.543 2.883 4.320 6.887 3952 1.002 r_s[1] -111.459 67.855 135.379 -511.946 -110.823 -61.673 -24.665 9.527 4 2.424 r_s[2] -113.088 67.878 135.432 -513.094 -112.537 -63.455 -25.973 8.144 4 2.423 r_s[3] -103.502 67.866 135.411 -504.068 -102.333 -53.603 -16.753 17.481 4 2.423 ... r_s[131] -107.190 67.873 135.434 -507.730 -106.498 -57.521 -20.353 13.884 4 2.423 r_s[132] -107.515 67.853 135.364 -507.870 -106.415 -58.097 -20.576 13.855 4 2.423 r_s[133] -111.187 67.771 135.336 -510.871 -111.277 -61.242 -24.794 10.734 4 2.419 s 7.120 0.001 0.075 6.978 7.070 7.120 7.170 7.269 4594 1.000 s_rs 4.514 0.005 0.371 3.841 4.256 4.493 4.746 5.281 4638 1.001 s_rt 6.977 0.019 1.300 4.885 6.062 6.839 7.736 9.934 4554 1.000 lp__ -12201.878 0.158 9.823 -12222.260 -12208.061 -12201.493 -12195.220 -12183.965 3888 1.000
traceplotを眺めてみても，各シミュレーションごとにめちゃくちゃな動きをしているわけで，これだと収束しないよなぁ...と思わざるを得ません．
これからどうするか
現状は，ここで詰まっているわけです．路線の効果を単純なランダム効果の形でモデルに突っ込むのが問題なのかなぁという気はしているので，違う形のモデルにできたらしたいのですが，さりとてどういうモデルが望ましいのかがわからず...そもそも子と親が多対多の形になる階層モデリングが間違っているんじゃないかといわれたら，それまでのような気もしますし...StanのマニュアルとかBUGS Bookとかに，こういう場合のモデル化の例とか載ってたりするんですかね...
ということで，現状暗礁に乗り上げているので，まだ今後も続くかもしれませんが，何かしらの解決策が見つからなければ次のエントリはしばらく上がらないと思います．そんなこんなでボチボチやっていきますです，はい．
