ベイズマックスへようこそ。
このブログでは,心理学と統計が大好きなmutopsyが,
統計に関するあれやこれやを自由気ままに投稿します。

2017年12月10日

Stanでpsychophysics──階層ベイズモデルで恒常法データを分析する──[スライド紹介]

 こんにちは。mutopsyです。この記事は,『StanとRでベイズ統計モデリング』読書会 (Osaka.Stan #5) のLT(ライトニングトーク)セッションで使用したスライドの紹介記事です。このスライドで用いたデータの配布,および分析に用いたRコード・Stanコードの紹介もします。

 このスライドでは,恒常法で集めたデータに対して累積正規分布関数をフィッティングする時に,(1) 参加者ごとに最尤推定した場合,(2) 参加者ごとにベイズ推定した場合,(3) 階層ベイズモデルでまとめて推定した場合の3つの推定方法で結果がどう変わるのかを比較しています。

 分析対象は,Webで収集したミュラー・リヤー錯視実験のデータです(データはこちらで公開しています[csvへのリンク])。このデータ(csv形式)はこちらをクリックするとダウンロードできます。

 データは1行1試行の形式で,7つの列(変数)が格納されています。
・1列目 TrialIndex:通し番号。
・2列目 Participant:実験参加者の番号。1-25までの整数値。
・3列目 Type:標準刺激の種類。O = 外向図形,I = 内向図形。スライドでは外向図形のみを使用。
・4列目 Length:比較刺激の長さ(単位はパーセント)。
・5列目 SelectC:比較刺激が「長い」と判断された (1) or 標準刺激が「長い」と判断された (0)。
・6列目 RT:反応時間(ミリ秒)。スライドでは未使用。
・7列目 PresenPos:刺激の提示位置。cs = 比較刺激が左,sc = 標準刺激が左。スライドでは未使用。


 まずは,最尤法を用いて実験参加者ごとの平均値と標準偏差を推定するRコードを紹介します。

Rコード01: [1] 参加者ごとに最尤推定する (スライド9枚目)
d <- read.csv("mldata_raw")  #データの読み込み
d <- d[d$Type == "O",]  #ここでは外向図形条件のみを抽出
n <- 25  #実験参加者の人数
Agg <- tapply(d$SelectC,list(d$Participant,d$Length),mean) #集計して選択率データ(行列形式)にする
n <- nrow(Agg)  #実験参加者の人数
M <- 1:n  #平均値を格納するためのベクトルを用意
SD <- 1:n  #標準偏差を格納するためのベクトルを用意
for (i in 1:n){  
  di <- data.frame(  #i番目の実験参加者のデータを抽出
  x = seq(50,150,10),  #比較刺激の幅
  Rate = Agg[i,]  #比較刺激の選択率データのベクトル
  )

  fit <- glm(formula = di$Rate ~ di$x, family = binomial(probit))  #glm()で最尤推定
  c1 <- coefficients((fit))
  p1 <- c(-c1[1]/c1[2],1/c1[2])
  M[i] <- p1[1]
  SD[i] <- p1[2]
}
M  #参加者ごとの平均値ベクトル
SD  #参加者ごとの標準偏差ベクトル

 続いて,Stanを用いて参加者ごとの平均・標準偏差を推定するコードを紹介します。赤色でハイライトした箇所(17行目)は,恒常法で収集したデータに累積正規分布関数(プロビット関数)を当てはめる肝の部分です。モデル上は,この累積正規分布関数が返す確率pをパラメータとするベルヌーイ分布から反応(どちらの刺激を選択するか)が生成されると仮定しています。この累積正規分布関数形を決める平均と標準偏差が推定したいパラメータとなります。最後のgenerated quantitiesブロックでは,参加者ごとに得られた平均と標準偏差の全体平均を計算しています。

Stanコード01: [2] 参加者ごとにベイズ推定する (スライド10-11枚目)
data {
  int N;  #全部で何試行か。今回は2200 (= 88試行 × 25人)。
  int NP;  #実験参加者の人数。今回は25。
  real Length[N];  #比較刺激の水準の配列 (単位は%) = {50, 60, …, 140, 150}
  int P[N];  #参加者番号の配列 = {1, 2, …, 24, 25}
  int <lower=0, upper=1> SelectC[N]; 	#比較刺激を選んだか否かの配列 = {0, 1, 1, …}
}

parameters{
  real<lower=0, upper=200> mu[NP];  #参加者一人ひとりの平均
  real<lower=0, upper=50> sigma[NP];  #参加者一人ひとりの標準偏差
}

transformed parameters{
  real<lower=0, upper=1> p[N];
  for (n in 1:N){
    p[n] = normal_cdf(Length[n], mu[P[n]], sigma[P[n]]);  #累積正規分布の確率を取得
  }
}

model{
  for (n in 1:N){
    SelectC[n] ~ bernoulli(p[n]);  #ベルヌーイ分布から反応データ (0/1) が生じるモデル
  }
}

generated quantities{  #参加者間の平均 (mu0とsigma0) を生成
  real<lower=0, upper=200> mu0;
  real<lower=0, upper=50> sigma0;
  mu0 = 0;
  sigma0 = 0;
  for (np in 1:NP){
    mu0 = mu0 + mu[np];
    sigma0 = sigma0 + sigma[np];
  }
  mu0 = mu0 / NP;
  sigma0 = sigma0 / NP;
}


 最後に,階層ベイズでまとめて推定するためのStanコードを紹介します。上記のコードと違う箇所(階層化に関わる部分)を赤色でハイライトしています。

Stanコード02: [3] 階層ベイズでまとめて推定 (スライド12-13枚目)
data {
  int N; 
  int NP; 
  real Length[N]; 
  int P[N]; 
  int <lower=0, upper=1> SelectC[N]; 
}

parameters{
  real<lower=0, upper=200> mu[NP];
  real<lower=0, upper=50> sigma[NP];
  real<lower=0, upper=200> mu0;  #全体平均
  real<lower=0, upper=50> sigma0;  #全体標準偏差
  real<lower=0, upper=200> s_mu0;  #参加者間の平均のばらつき
  real<lower=0, upper=200> s_sigma0;  #参加者間の標準偏差のばらつき
}

transformed parameters{
  real<lower=0, upper=1> p[N];
  for (n in 1:N){
    p[n] = normal_cdf(Length[n], mu[P[n]], sigma[P[n]]);
  }
}

model{
  for (np in 1:NP){
    mu[np] ~ normal(mu0,s_mu0);  #参加者ごとの平均を得る
    sigma[np] ~ normal(sigma0,s_sigma0);  #参加者ごとの標準偏差を得る
  }
  for (n in 1:N){
    SelectC[n] ~ bernoulli(p[n]);
  }
}

 実際に上記のコードを走らせるか,あるいはスライドをご覧頂けば詳細な推定結果を確認できます。ざっくり結果をまとめると,平均の点推定値はどの推定方法でもほとんど変わらなかったものの,区間推定の幅は(見かけ上)推定方法によって異なっていました。ただし,この結果からどのモデルが「正しい」かを議論することはできません(情報量基準などを用いて比較することは一応できますが,少なくとも今回のケースではあまり意味がないような気がします)。モデルによって置かれている仮定が異なるというだけに過ぎません。今回用いたような比較的綺麗なデータに対しては,少なくとも錯視量の平均値に関しては推定方法によらずそれなりにうまく推定できるようですが,もっと微妙なデータになると挙動が変わってくるかもしれません。まだまだ遊び甲斐がありそうです。

posted by mutopsy at 00:00 | 統計

2017年12月09日

回帰分析の悩みどころ (「アヒル本」7.1-7.5) [スライド紹介]

 こんにちは。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コードの各行の意味までしっかりと説明がなされているので,実際にコードを走らせながら読めば相当身になると思います。

posted by mutopsy at 00:00 | 統計

2017年12月08日

ベイズ推定の初歩──二項分布を例に── [スライド紹介]

 こんにちは。mutopsyです。この記事では,約2年前に公開した「ベイズ推定の初歩──二項分布を例に──」というスライドを紹介します。



 この資料は,2015年7月18日に関西学院大学で開催された読書会で用いたものです。文献は“Bayesian Cognitive Modeling: A Practical Course”,通称「コワイ本」と呼ばれる本です。今年の夏,日本語訳が出版されて再び話題に上がりましたね(『ベイズ統計で実践モデリング: 認知モデルのトレーニング』)。
 このスライドでは,コワイ本の第3章「Inferences with binomials(二項分布を使った推論)」について解説をしています。内容としては,ベイズ推定の基本的な考え方と,確率的プログラミング言語の1つであるWinBUGSを使って実際にベイズ推定を行う方法について,二項分布に従う変数を例に説明しています。二項分布は整数値を返す離散分布で,コインを投げて表が出る回数であったり,n問中k問正解する回数であったりと,直感的に理解しやすい分布なので,導入としてよく用いられます。しかし,ここで紹介している基本的な内容は,正規分布のような連続分布に対しても適用可能です。そういった理由もあって,スライドのタイトルを「ベイズ推定の初歩──二項分布を例に──」としました。
 当時の僕は博士前期課程2年で,まだベイズの「ベ」の字も知らないような初心者でしたが,この本を読んでこの資料を作って,なんとなくベイズの「ヘ」ぐらいまでは分かったかなあといったところです。今であれば,『ベイズ統計モデリング: R,JAGS, Stanによるチュートリアル 原著第2版』,通称「犬4匹本」や,『StanとRでベイズ統計モデリング (Wonderful R) 』といった,日本語で読める良著も増えており,ベイズを学ぶ敷居はだいぶ低くなりつつあるように思います。インターネット上にも有志の方々による分かりやすい資料がたくさん公開されています。とてもいい時代ですね。
 ちなみに,コワイ本で用いられているWinBUGSは現在開発が終了していますが,その後継であるJAGSが利用できます。また,より効率的なアルゴリズムを採用しているStanもよく使われています(筆者もStanを愛用しています)。さらに,特定の目的に特化した,簡単に使えるソフトウェアやパッケージも増えつつあります。今やベイズはコワくない。Let's Do Bayesian Data Analysis!
posted by mutopsy at 00:00 | 統計

2017年12月07日

統計ブログを作りました。

こんにちは。mutopsyこと武藤拓之です。
大阪大学で認知心理学の研究をしています。
心理学と統計が好きな,どこにでもいる普通のベイジアンです。

勢いあまって統計ブログというものを作ってしまいました。
RやStanを使って統計に関する記事を書いたり,あるいは書かなかったりする予定です。

このブログの名前はベイズマックス。あなたの知的好奇心をくすぐります。さういふブログに私はしたい。
何卒,よろしくお願いいたします。
posted by mutopsy at 16:08 | 雑記