ベイズ統計によるアイオワ・ギャンブリング課題のモデリング方法

こんにちは。データ分析部の長尾です。

今回は心理学において有名な実験である「アイオワ・ギャンブリング課題」を題材として,ベイズ統計によるモデリング方法を紹介します。

突然ですが、パチンコのホールを思い浮かべてください。パチンコではまず,出玉の多そうな台を野生の勘で選び,その台でしばらく様子を見ます。出玉が少ないようであれば,他の台に移り,またしばらく打ちます。個人差はありますが,上述の探索行動を繰り返し,最終的には,一番期待できそうな台に落ち着いて,利益を最大化するよう努めます。このような意思決定のプロセスを模倣した実験のひとつにアイオワ・ギャンブリング課題(Iowa Gambling Tasks, IGT; Bechara, Damasio, Damasio & Anderson, 1994)があります。

IGTとは

IGTにおいて,実験参加者には,あらかじめ2000ドルの所持金がヴァーチャルに付与されます。参加者は図1で示したような4つのカードのデッキから,どれか1つのデッキを自由に選び,1枚カードを引く試行を繰り返し行います。

図1: 実験イメージ

カードの裏には報酬額もしくは損失額が記載されており,デッキ毎にそれらの金額と出現頻度は違います。表1に示したように,AとBのデッキからは毎回100ドルの報酬が得られるのに対し,CとDのデッキからは毎回50ドルと控えめです。しかしながら,10枚当たりの損失額はAとBのデッキでは1250ドルと大きいのに対し,CとDのデッキでは250ドルと抑えられています。結局,10枚当たりの純利益は,AとBのデッキで-250ドル,CとDのデッキで250ドルとなり,CもしくはDのデッキを選ぶ方が明らかに得です。

表1: デッキの性質(単位:ドル)

参加者は実験を通して,所持金を最大化するよう教示されます。このため,試行を繰り返し,デッキを探索しながら,それぞれのデッキの性質を学習することが参加者には求められます。最終的には,AとBのデッキを避け,CもしくはDのデッキを選ぶことが望ましい方略となります。しかしながら,健常群と比較して,麻薬中毒や強迫性障害といった臨床群ではこの方略を獲得する能力が低いことが指摘されており(Wetzels, Vandekerckhove, Tuerlinckx & Wagenmakers, 2010),IGTは様々な臨床群における意思決定能力の欠如の度合いを評価するために用いられます。

PVL-Deltaモデル

PVL-Deltaモデル(prospect valence learning delta model)はIGTにおける意思決定を評価するモデルとして, Ahn et al.(2008)によって提案されました。PVL-Deltaモデルは過去の試行に対する強化学習をルール化したモデルであり,4つのパラメータによって,参加者の心理プロセスを表現します。

まず,参加者は試行tで得られた純利益x^{(t)}から主観的効用u^{(t)}を評価すると仮定し,行動経済学で有名なプロスペクト理論を用いて,u^{(t)}を以下のように表現します。



ここで,\alpha(0 < \alpha < 1)は効用関数の形状母数であり,\lambda(0 < \lambda < 5)は損失の感度を表す母数です。\lambda > 1のとき,参加者は報酬よりも損失に敏感であるといえます。

次に,試行tにおけるデッキの選択をy^{(t)}とおき,参加者は得られた効用u^{(t)}から,次の試行t+1でデッキkから得られる期待効用E_{u_k}^{(t+1)}を以下のように評価すると仮定します。



母数a(0 < a < 1)は1つ前で得られた効用に対する注意の度合いを表し,デッキkが選ばれた場合のみ,期待効用を更新することから,更新比率とも呼ばれます。aの値が高い参加者は1つ前の試行で得られた効用に影響されやすく,それまでの経験を忘却する傾向があると解釈します。

さらに,試行tにおけるデッキの選択確率p_k^{(t)}をsoftmax関数を用いて,以下のように表現します。

p_k^{(t)}=\frac{\exp { (\theta\times E_{u_k}^{(t)}) }}{\sum_{j=1}^4 \exp { (\theta\times E_{u_j}^{(t)}) } }

ここで,\thetaはデッキの選択に対する感度を表し,0に近づくほどランダムな選択が行われ,大きい値をとるほど,期待効用に基づく選択が行われることを意味します。PVL_deltaモデルでは,選択の一貫性を表すパラメータc(0 < c < 5)を用いて,\thetaに以下のような変換が施されています。

\theta=3^c-1

すなわち,一貫性cの値が0に近いほどランダムな選択が行われ,逆にcの値が大きくなるにつれ,より確定的な選択が行われます。

最後に,観察される試行tにおけるデッキの選択y^{(t)}をカテゴリカル分布を用いて,以下のようにモデル化します。

y^{(t)} \sim Categorical(p_1^{(t)}, p_2^{(t)}, p_3^{(t)}, p_4^{(t)})

ひとりの参加者を分析

当社の社員7名(A~Gさん)から得たデータを用いて,ここでは参加者の1人であるDさんの分析を行います。PVL-Deltaモデルについて,母数\alphaaには区間(0,1)の一様分布,\lambdacには区間(0,5)の一様分布を仮定し,ベイズ推定を行います。

分析結果

Dさんの母数の事後分布の事後平均(EAP),事後標準偏差(post.sd),95%確信区間は表2のようになりました。

表2: Dさんの母数の推定結果

表2より,Dさんの更新比率aのEAP推定値は0.166と低い値を示しています。これは,Dさんは1つ前の試行で得られた効用にさほど影響されず,デッキの期待値を更新していることを示唆します。デッキの選択に対する一貫性はc=0.247であり,やや探索行動が多かったといえます。母数\alpha\lambdaのEAPを用いて描いたプロスペクト曲線は図2のようになりました。

図2: Dさんのプロスペクト曲線

図2より,Dさんは報酬より損失にやや敏感な傾向があることがみてとれます。

ここまで,ひとりの参加者に的を絞った分析を行ってきました。ベイズ推定の利点は母数に対する確信を区間で表現できる点ですが,表2における\lambdaの95%確信区間はやや広く,EAPへの十分な確信を保障しません。また,母数間での相関関係もひとりの参加者の分析からは明らかにできません。次節ではこれらの問題に対処する方法を紹介します。

階層ベイズモデルを用いた分析

階層ベイズモデルでは,参加者を表す添え字iを用いて,\alpha_i, \lambda_i, a_i, c_iと表現し,これら母数がさらに,それぞれ未知である超母数(hyper parameter)をもつ事前分布から発生するようにモデリングします。

心理学的な特性値である\alpha_i, \lambda_i, a_i, c_iは参加者間で正規分布していると仮定します。しかしながら,これら特性値の定義域は異なっているため,未知の平均ベクトルと共分散行列を持つ4変量正規分布から\bold{\delta}を発生させ,以下の変換を行います。

\alpha_i=\Phi(\delta_1|\mu_1, \sigma_1)
\lambda_i=5\times\Phi(\delta_2|\mu_2, \sigma_2)
a_i=\Phi(\delta_3|\mu_3, \sigma_3)
c_i=5\times\Phi(\delta_4|\mu_4, \sigma_4)

ここで,\Phi(|)は正規累積変換を表し,その値は0から1の間に収まります。\lambda_ic_iに関しては,さらに5を乗じることにより区間(0, 5)を実現しています。

超母数である4変量正規分布の平均ベクトルと共分散行列にはそれぞれ以下の正規分布と逆ウィシャート分布を事前分布として仮定します。

\mu_1, \mu_2, \mu_3, \mu_4 \sim Normal(0, 1)
\bold{\Sigma} \sim Wishart^{-1}(4, I_4)

最後に,\bold{\Sigma}から\alpha_i, \lambda_i, a_i, c_iに関する相関行列\bold{\Sigma}_\rhoを生成し,相関関係を検証します。

\bold{\Sigma}_\rho=diag(\bold{\Sigma})^{-1/2}\times\bold{\Sigma}\times diag(\bold{\Sigma})^{-1/2}
分析結果

当社の社員7名から得たデータ全てを用いた場合の階層ベイズモデルによる推定結果を表3に示します。

表3: 社員A~Gさんの推定結果

表3より,更新比率aが最も高い参加者はDさんであり,最も低い参加者はAさんです。過去の経験を忘れっぽいDさんとは対照的に,Aさんは過去の経験を鑑みた選択が出来ていたといえます。また,Dさんの\lambdaの事後標準偏差も0.439と抑えられ,前述の分析の問題点が改善されていることがわかります。

母数間の相関関係は表4のようになりました。

表4: 母数間の相関行列

表4より,\lambdaaの間には0.605と中程度の正の相関がみられます。この結果から,損失に対して敏感な参加者ほど,1つ前の試行で受けた損失に影響されやすいため,更新比率も高い値となるメカニズムが予想されます。また,acの間には,-0.603と中程度の負の相関が確認できます。これは過去の経験を鑑みた選択を行っている参加者ほど,同一のデッキを選択していたことを示唆します。

最後に,7名の社員それぞれのプロスペクト曲線を図3に示します。

図3: 社員A~Gさんのプロスペクト曲線

図3より,Aさん(赤)やEさん(紫)は純利益が正であれば,効用が極端に高くなり,純利益が負であれば,効用が極端に低くなる傾向があります。個人的には,AさんとEさんは,実験に対して並々ならぬモチベーションを有したギャンブラー気質の持ち主であるという印象を受けました。また,Bさん(緑)とCさん(青)は正の純利益に関しては,同じように効用を得ていますが,負の純利益に関しては,Cさんの方がさらに敏感に反応しており,プロスペクト理論の標準的な結果に近くなっています。

Stan・Rコード

まず,階層モデルのMCMCモデリングを行ったStanコードを以下に示します(‘PVL_delta_multi.stan’として保存)。

data {
    int<lower=1> T;
    int<lower=1> K;
    int<lower=1> N;
    matrix[4,4] R;
    int<lower=1,upper=K> y[T,N];
    int<lower=-1150,upper=100> x[T,N];
}
parameters {
    vector[4] Eta[N];
    vector[4] Mu;
    cov_matrix[4] Sigma;
}
transformed parameters {
    real<lower=0,upper=242> theta[N];
    real u[T,N];
    vector[K] Eu[T,N];
    simplex[K] p[T,N];
    real<lower=0,upper=1> alpha[N];
    real<lower=0, upper=5> lambda[N];
    real<lower=0,upper=1> a[N];
    real<lower=0, upper=5> c[N];

    for(n in 1:N){
        alpha[n] = normal_cdf(Eta[n,1],Mu[1],sqrt(Sigma[1,1]));
        lambda[n] = 5*normal_cdf(Eta[n,2],Mu[2],sqrt(Sigma[2,2]));
        a[n] = normal_cdf(Eta[n,3],Mu[3],sqrt(Sigma[3,3]));
        c[n] = 5*normal_cdf(Eta[n,4],Mu[4],sqrt(Sigma[4,4]));
    }
    for(n in 1:N){
        theta[n] = pow(3,c[n])-1;
        for(t in 1:T){
            if(x[t,n]>=0) u[t,n] = pow(x[t,n],alpha[n]);
            else u[t,n] = -lambda[n]*pow(fabs(x[t,n]),alpha[n]);
        }
    }
    for(n in 1:N){
        for(k in 1:K){
            Eu[1,n,k] = 0;    // 1回目の試行での期待効用を0に設定
            p[1,n,k] = 0.25;  // 1回目の試行での各デッキの選択確率を等確率に設定
        }
    }
    for(t in 2:T){
        for(n in 1:N){
            for(k in 1:K){
                if(y[t-1,n]==k)
                    Eu[t,n,k] = (1-a[n])*Eu[t-1,n,k] + a[n]*u[t-1,n];
                else
                    Eu[t,n,k] = Eu[t-1,n,k];
            }
        p[t,n] = softmax(theta[n]*Eu[t,n]);
        }
    }
}
model {
    Mu~normal(0,1);
    Sigma~inv_wishart(4,R);
    Eta~multi_normal(Mu,Sigma);

    for(n in 1:N){
        for(t in 1:T){
            y[t,n]~categorical(p[t,n]);
        }
    }
}
generated quantities{
    matrix[4,4] Rho;

    for(i in 1:4){
        for(j in 1:4){
            Rho[i,j] = Sigma[i,j]/sqrt(Sigma[i,i]*Sigma[j,j]);
        }
    }
}

MCMC分析を行うRコードは以下のとおりです。

# ライブラリ読み込み
library(rstan)

# データの読み込み
dat1<-read.csv("IGT_data/igtlog-A.csv",header=T)
dat2<-read.csv("IGT_data/igtlog-B.csv",header=T)
dat3<-read.csv("IGT_data/igtlog-C.csv",header=T)
dat4<-read.csv("IGT_data/igtlog-D.csv",header=T)
dat5<-read.csv("IGT_data/igtlog-E.csv",header=T)
dat6<-read.csv("IGT_data/igtlog-F.csv",header=T)
dat7<-read.csv("IGT_data/igtlog-G.csv",header=T)

# Stanコードでの変数を設定
T<-max(dat1$trialnum) # 試行数(150回)
K<-4  # デッキ数
N<-7  # 参加者数
y<-structure(c(dat1$deck,dat2$deck,dat3$deck,dat4$deck,dat5$deck,dat6$deck,dat7$deck),.Dim=c(T,N))
x<-structure(c(dat1$net,dat2$net,dat3$net,dat4$net,dat5$net,dat6$deck,dat7$deck),.Dim=c(T,N))
R<-diag(1,4)

dat<-list(T=T, K=K, N=N, y=y, x=x, R=R)

# モデルのコンパイル
stanmodel <- stan_model(file='PVL_delta_multi.stan')

# MCMCサンプリング
fit <- sampling(
    stanmodel,
    data=dat,
    seed=1234,
    chains=3, iter=10000, warmup=5000)

#結果の表示
print(fit)
参考文献
  • Ahn, W.-Y., Busemeyer, J. R., Wagenmakers, E.-J., & Stout, J. C. (2008, December). Comparison of decision learning models using the generalization criterion method. Cognitive science,32(8),1376–402.
  • Bechara, A, Damasio, A. R., Demasio, H., & Anderson, S. (1994). Insensitivity to future consequences following damage to human prefrontal cortex. Cognition, 57, 7-15.
  • Steingroever, H., Wetzels, R., & Wagenmakers, E. J. (2013). Validating the PVL-Delta model for the Iowa gambling task. Frontiers in psychology, 4.
  • Wetzels, R., Vandekerckhove, J., Tuerlinckx, F., & Wagenmakers, E. J. (2010). Bayesian parameter estimation in the Expectancy Valence model of the Iowa gambling task. Journal of Mathematical Psychology, 54(1), 14-27.