こんにちは。mutopsyです。この記事は，『StanとRでベイズ統計モデリング』読書会 (Osaka.Stan #4) で使用したスライドの紹介記事です。
『StanとRでベイズ統計モデリング』，通称「アヒル本」のChapter 7「回帰分析の悩みどころ」の前半部分を解説しているスライドです。回帰分析を用いる際に留意するべき点，あるいは工夫できる点として，「交互作用」「対数をとるか否か」「非線形の関係」「多重共線性」「交絡」の5つに焦点を当てています。内容に関してはスライドをご覧頂くか，「アヒル本」を読んで頂ければ良いかと思いますので，この記事ではスライド内で紹介しているStanコードを（コピペしやすいように）改めて紹介したいと思います。
以下のStanコードは，2つの説明変数（イケメンか否か・年収）で1つの応答変数（魅力）を予測する重回帰モデルです（この例は清水先生の記事を参考にしました）。まずは交互作用項を含まないモデル（すなわち，交互作用に関わるパラメータが0であると仮定するモデル）から。
Stanコード01: 交互作用を含まないモデル (スライド11枚目)
data {
int N; //データの数
int<lower=1, upper=10> Miryoku[N]; //魅力の評定値データ
int<lower=0, upper=1> Ikemen[N]; //イケメンか否かの2値データ
real<lower=0> Nenshu[N]; //年収データ
}
parameters{
real b[3]; //回帰係数 b[1], b[2], b[3]
real sigma; //残差の標準偏差
}
model{
for (n in 1:N){
Miryoku[n] ~ normal(b[1] + b[2]*Ikemen[n] + b[3]*Nenshu[n], //重回帰モデル
sigma);
}
}
このコードのモデル部分を以下のように書き換えて，2つの説明変数を掛けたものを新たな項として加えることで，交互作用項を含むモデルになります。
Stanコード02: 交互作用を含むモデル (スライド13枚目)
model{
for (n in 1:N){
Miryoku[n] ~ normal(b[1] + b[2]*Ikemen[n] + b[3]*Nenshu[n]
+ b[4]*Ikemen[n]*Nenshu[n], //この部分が交互作用項を表す
sigma);
}
また，以下のtransformed paramatersブロックをmodelブロックの前に挿入することで，回帰係数の和についても推定することができます。これは，年収が魅力に与える効果の大きさ（傾き）を，イケメンと非イケメンのそれぞれについて推定することに相当します。
Stanコード03: パラメータの和を推定 (スライド16枚目)
transformed parameters{
real Intrcpt_I0;
real Slope_I0;
real Intrcpt_I1;
real Slope_I1;
Intrcpt_I0 = b[1];
Slope_I0 = b[3];
Intrcpt_I1 = b[1] + b[2];
Slope_I1 = b[3] + b[4];
}
非線形な関係を表すモデルも，モデル式を変えてやれば簡単に表現できます。
Stanコード04: 二次曲線の当てはめ (スライド33枚目)
data {
int N; //データの数
real X[N]; //変数X（説明変数）
real Y[N]; //変数Y（応答変数）
}
parameters{
real a;
real b;
real x0;
real<lower=0> s_Y;
}
model{
for (n in 1:N){
Y[n] ~ normal(a + b*(X[n]-x0)^2, s_Y);
}
}
Stanコード05: 時系列データへの指数曲線の当てはめ (スライド35枚目)
data {
int N;
real X[N];
real Y[N];
}
parameters {
real<lower=0, upper=100> a;
real<lower=0, upper=5> b;
real<lower=0> s_Y;
}
model {
for (n in 1:N)
Y[n] ~ normal(a*(1 - exp(-b*X[n])), s_Y);
}
最後に，多重共線性がある場合（＝説明変数同士が強く相関している場合）に何が起こるかをシミュレーションするためのRコードを紹介します。ここではベイズではなく普通にlm()を使って得られる結果を用いています。パラメータをいろいろ変えてみながら遊んでみると楽しいです。
Rコード01: 多重共線性のあるデータに重回帰分析を適用するとどうなるか (スライド42枚目)
##パラメータ等の指定##
rAB <- 0.9 #説明変数間の母相関係数。この値をいろいろ変えてみる。
rAY <- 0.5 #AとYの母相関係数
rBY <- 0.6 #BとYの母相関係数
n <- 100 #サンプルサイズ
Rep <- 1000 #サンプリング回数
##分散・共分散行列の作成##
Mat <- matrix(c(1, rAB, rAY, rAB, 1, rBY, rAY, rBY, 1), ncol=3)
##N = nのデータセットをRep回生成し，それぞれに対して重回帰分析##
Res_beta <- data.frame(b2 = 1:n, b3 = 1:n)
for (i in 1:Rep){
d <- as.data.frame(mvrnorm(n= n, mu= c(0, 0, 0), Sigma= Mat, empirical= FALSE))
colnames(d) <- c("A","B","Y")
reg <- lm(Y ~ A + B, data = d)
Res_beta[i,] <- reg$coefficients[2:3]
}
Res_beta #Rep回分の偏回帰係数 (b2とb3) が格納されたデータフレーム
思い返してみると，このOsaka.Stanでの発表が，僕にとってのベイズ（とStan）の始まりだったように思います。Osaka.Stanに参加するまでまともにStanを動かしたことのなかった僕ですが，なんやかんやでそれなりに自由にモデルが書けるようになりました。モデリングの感覚をつかむには，自分で手を動かしてみるのが一番だと思います。その点，アヒル本ではStanコードの各行の意味までしっかりと説明がなされているので，実際にコードを走らせながら読めば相当身になると思います。
