ベイズマックスへようこそ。
統計に関するあれやこれやを自由気ままに投稿します。
全記事一覧はこちらから。

2017年12月11日

StanとRで折れ線回帰──空間的視点取得課題の反応時間データを説明する階層ベイズモデルを例に──[スライド紹介]

 こんにちは。mutopsyです。この記事は,『StanとRでベイズ統計モデリング』読書会最終回 (Osaka.Stan #6) のLT(ライトニングトーク)セッションで使用したスライドの紹介記事です。この記事では,このスライドで用いたStanコードの紹介と補足説明を行います。

 このスライドでは,空間的視点取得課題の反応時間データに対し,Stanを使って(階層)折れ線回帰モデル(分割点が1つの場合)を適用した例を紹介しています。その背景や詳しい方法については上記スライドをご覧ください。

 分析で用いたStanコードは以下の通りです。データとしては,各行がある参加者のある条件に対応したデータを用いています。使用している変数(列)は,「実験参加者ID」「角度(0度・40度・80度・120度・160度)」「平均反応時間」の3つです。したがって,行数は,角度の水準数×実験参加者の人数となります。もちろん,以下のコードは他のデータセットにも適用可能です(事前分布やパラメーターの範囲などを修正する必要はありますが)。

Stanコード01: 階層折れ線回帰モデルのStanコード例 (スライド12枚目)
data {
  //分析するデータ
  int N;  //合計試行数 (データの行数)
  int Np;  //実験参加者の人数
  real Angle[N];  //角度の配列
  real RT[N];  //平均反応時間の配列
  int Par[N];  //実験参加者IDの配列(整数)
  
  //以下は,平均値の事後分布を求めるために使用。
 //分析そのものとは無関係なので,不要なら消す。
  int N_new;  //事後分布用のN
  real Angle_new[N_new];  //事後分布用の角度の配列
}

parameters {
  //4つのパラメータの全体平均
  vector[4] mu0;  //BP, b0, b1, b2の全体平均を格納するベクトル

  //実験参加者ごとのパラメータ
  real<lower=1, upper=159> BP[Np];  //分割点:適切な範囲を指定しないとうまく収束しない。
  real<lower=0> b0[Np];  //切片:反応時間なので非負の実数。
  real b1[Np];  //分割点未満の角度における傾き
  real b2[Np];  //分割点以上の角度における傾き
  real<lower = 0> resid[Np];  //残差の標準偏差

  //パラメータ間の分散・共分散構造(階層化のために使用)
  cov_matrix[4] cov;  //パラメータ間の分散・共分散行列
}

transformed parameters{
  //階層化しやすいように,実験参加者ごとのパラメータベクトルを作る。
  vector[4] mu_p[Np];  //1人分のパラメータを格納するベクトル × 人数分
  mu_p[,1] = BP;
  mu_p[,2] = b0;
  mu_p[,3] = b1;
  mu_p[,4] = b2;
}

model {
  //階層化(こうすると短く書けて便利)
  for (np in 1:Np){
    mu_p[np,] ~ multi_normal(mu0, cov);
    //↑4つのパラメータの個人差が正規分布に従うと仮定。
    // 多変量正規分布を指定することで,パラメータ間の相関も推定できる。
  }

  //折れ線回帰モデル
  for (n in 1:N){
    if (Angle[n] < BP[Par[n]]){
      //角度が分割点未満の時のモデル
      RT[n] ~ normal(b0[Par[n]] + (b1[Par[n]] * Angle[n]),
                     resid[Par[n]]);
    } else {
      //角度が分割点以上の時のモデル
      RT[n] ~ normal(b0[Par[n]] + (b1[Par[n]] * BP[Par[n]]) 
                      + (b2[Par[n]] * (Angle[n] - BP[Par[n]])),
                     resid[Par[n]]);
    }
  }

  //事前分布
  for (np in 1:Np){
    resid[np] ~ cauchy(0,1000);  //残差の標準偏差の事前分布(半コーシー分布)
    BP[np] ~ normal(90,180);  //分割点の事前分布として弱情報を与えておく
  }
}

generated quantities{
  corr_matrix[4] rho;  //4つのパラメータの相関行列
  real RT_new[N_new];  //平均値の事後分布用

  //パラメータ間の相関係数を求める
  for (i in 1:4){
    for (j in 1:4){
      rho[i,j] = cov[i,j]/sqrt(cov[i,i]*cov[j,j]);
    }
  }

  //平均値の事後分布を求める(不要なら削除)
  for (n in 1:N_new){
    if (Angle_new[n] < mu0[1]){
      RT_new[n] = mu0[2] + (mu0[3] * Angle_new[n]);
    } else {
      RT_new[n] = mu0[2] + (mu0[3] * mu0[1]) + (mu0[4] * (Angle_new[n] - mu0[1]));
    }
  }
}


 なお,スライドでも触れている通り,if文を含むモデルはStanでは収束効率が悪くなります。モデルの書き方を変えることで収束効率が改善する可能性があるので,後日,また比較してみて記事にしようと思います。

 ちなみに余談ですが,この記事およびスライドのタイトル(「StanとRで〜」)は,アヒル本の書名(『StanとRでベイズ統計モデリング』)のオマージュです。Osaka.Stanでのアヒル本読書会最終回ということで,それにちなんだタイトルにしてみました。アヒル本読書会は終わりましたが,Osaka.Stan自体はまだまだ続くそうですので,乞うご期待。

posted by mutopsy at 04:00 | 統計