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