私 実験心理学者のための統計モデリング

私がベイジアンになってから1年近く経過し、その間にそれなりのモデリング経験を積んできました。ただ、「ベイズモデリングについて勉強する」ことと「実際にモデリングを研究で使う」こととの間に、そこそこ大きなギャップがあったように感じています。実際に学会等で皆さんのご意見を聞いていると、こういったギャップが実験心理学者をモデリングから遠ざける要因になっていると思います。そこで本セッションでは、なるべく実験心理学者の視点に立って、実験心理学の分野でベイズモデリングをどのように使っていくことができるのか一例を示します。実験心理学においてベイズモデリングの導入や利用法に関する議論の材料となれば幸いです。

本セッションの内容のうち、「統計モデリング入門」はStatistical Rethinking第2版 および、その内容をtidyverseとbrmsで再現したStatistical Rethinking with brms, ggplot2, and the tidyverseにかなり依拠していますが、内容に誤りがあった場合は私の勘違いや誤りである可能性が高いです。ご指摘ください。

研究に対する姿勢について考え直すきっかけをくれた先輩に「なんでモデリング？」と聞かれたときに、ろくな返事ができなかったことをずっと後悔していました。これはそのときの返事の代わりです。統計モデリングってのもちょっとおもしろそう、とその人に思ってもらえたら幸せだな、と妄想しながら書きました。

導入

「ベイズについて勉強はしてみたが、どう使うものかがよくわからない」という声をときどき耳にします（実際に、私が先輩方からこの質問をされたこと、それに対して満足のいく返答ができなかったことが本セッションの動機になっています）。そこで、本セッションでは実験心理学者には馴染み深い課題である視覚探索課題 (visual search) を例に、ベイズモデリングをどう使うことができるのか（あるいは使わないでおくのがよいのか）を考えたいと思います。

実験心理学者でない人のための視覚探索課題入門

私たちの周囲には多くの情報があり、それらすべてを一度に処理することはできません（今、皆さんはこの文章を読んでいると思いますが、その間にも触覚や聴覚からさまざまな情報が入ってきているはずです。たとえば椅子の感触はどうでしょう？今この瞬間まで椅子の感触には“注意”していなかったと思いますが、その間もずっと椅子の座面の感覚は入力されてきているはずです。なぜ常に椅子の感触を感じないでいられるのでしょう？）。注意研究の特徴は、豊富かつ独創的な実験パラダイムです（私個人の意見です）。

なかでも視覚探索課題は、視覚情報の選択システムである視覚的注意 (visual attention) の働きを調べるために広く用いられる課題です。

図1. T (標的) を探し、左向きか右向きかを答えてください

上に、典型的な視覚探索課題の試行画面を示しました。参加者の課題は標的刺激 (ここではT) を妨害刺激 (ここではL) の中から探すことです。

標的刺激に対して効率的に注意を誘導できるなら、Tをすぐに見つけることができます。たとえばこの画面でTのみが赤色で描かれていたとしたら、わざわざ探そうとしなくてもひとりでにTが目に飛び込んでくるような感じを覚えるでしょう。注意の誘導が非効率的なときには、各刺激に対して順番に注意を向けていく必要があり、そのために標的が見つかるまでの時間は長くなります。このような状況を非効率探索と呼ぶことにします。非効率探索の状況下では、画面上の刺激の数 (セットサイズ) が増えるほど反応時間は長くなります。反応時間は線形に増えていくことがわかっていて、この増加の様子を\[y=ax+b\]という一次関数 (探索関数) であらわし、傾き\(a\)を探索勾配と呼んで探索の効率性の指標として使うこともあります。

統計モデリング入門

上に示したような視覚探索課題のサンプルデータを使って、実際に統計モデリングをやってみます。その前に、ベイズ統計モデリングとはいったい何をやっているのか、簡単に紹介します。

n <- 1001
n_water <- 6
n_trials  <- 9

d <- tibble(p_grid     = seq(from = 0, to = 1, length.out = n),
            prior      = 1) %>% 
  dplyr::mutate(likelihood = dbinom(n_water, size = n_trials, prob = p_grid)) %>% 
  dplyr::mutate(posterior  = (likelihood * prior) / sum(likelihood * prior))

d %>%
  dplyr::rename("パラメタ" = p_grid, "事前確率" = prior,
                "尤度" = likelihood, "事後確率" = posterior) %>%
  knitr::kable() %>%
  kableExtra::kable_styling() %>%
  scroll_box(width = "100%", height = "200px")

# how many samples would you like?
n_samples <- 1e4

# make it reproducible
set.seed(3)

samples <-
  d %>% 
  sample_n(size = n_samples, weight = posterior, replace = T)

samples %>% 
  mutate(sample_number = 1:n()) %>% 
  ggplot(aes(x = sample_number, y = p_grid)) +
  geom_line(size = 1/10, color = "skyblue4") +
  labs(x = "サンプル",
       y = "海が占める割合 (p)",
       title = "図7: サンプリングのイメージ")

簡単な線形モデル

'data.frame':   10000 obs. of  3 variables:
 $ RT       : num  594 508 512 483 416 ...
 $ ID       : int  1 1 1 1 1 1 1 1 1 1 ...
 $ condition: num  5 5 5 5 5 5 5 5 5 5 ...

data {
  int<lower=1> N;
  vector[N] setsize;
  vector[N] RT;
}

parameters {
  real alpha;
  real beta;
  real<lower=0> sigma;
}

model {
  vector[N] mu = alpha + beta * setsize;
  target += normal_lpdf(RT | mu, sigma);
  target += normal_lpdf(alpha | 0, 10);
  target += normal_lpdf(beta | 0, 10);
  target += cauchy_lpdf(sigma | 0, 10) - cauchy_lccdf(0 | 0, 10);
}

install.packages("brms")
library(brms)

dat1 <- list(N = NROW(dat),
             ID = dat$ID,
            setsize = dat$condition,
            RT = dat$rt.s)

fit1 <- brm(RT ~ 1 + setsize,
            family = "gaussian",
            data = dat1,
            iter = 3000,
            warmup = 1000,
            chain = 2,
            cores = 2,
            prior = c(prior(normal(0, 10), class = Intercept),
                      prior(normal(0, 10), class = b),
                      prior(cauchy(0, 10), class = sigma))
            )

summary(fit1)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: RT ~ 1 + setsize 
   Data: dat (Number of observations: 10000) 
Samples: 2 chains, each with iter = 3000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    -1.75      0.03    -1.80    -1.70 1.00     3813     3187
setsize       0.23      0.00     0.23     0.24 1.00     3935     2979

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.81      0.01     0.80     0.82 1.00     3407     2846

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

bayesplot::color_scheme_set("teal")
plot(fit1)

mod2 <- rstan::stan_model("200211_mod_search1.stan")

fit2 <- sampling(mod2, 
                 data = dat2, 
                 iter = 3000, 
                 warmup = 1000,
                 chains = 2, 
                 cores = 2)

Inference for Stan model: 200211_mod_search1.
2 chains, each with iter=3000; warmup=1000; thin=1; 
post-warmup draws per chain=2000, total post-warmup draws=4000.

           mean se_mean   sd      2.5%       25%       50%       75%     97.5%
alpha     -1.75    0.00 0.03     -1.80     -1.77     -1.75     -1.73     -1.70
beta       0.23    0.00 0.00      0.23      0.23      0.23      0.24      0.24
sigma      0.81    0.00 0.01      0.80      0.81      0.81      0.82      0.82
lp__  -12111.99    0.03 1.19 -12115.07 -12112.52 -12111.69 -12111.12 -12110.62
      n_eff Rhat
alpha  1473    1
beta   1521    1
sigma  1793    1
lp__   1277    1

Samples were drawn using NUTS(diag_e) at Tue Feb 18 11:23:56 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

交互作用

視覚探索課題のもう一つの特徴は、非効率探索と効率探索の存在です。図1のような、刺激一つ一つに注意を向けていかないといけないような探索を非効率探索と呼びます。非効率探索の特徴が、上でモデリングしたセットサイズの効果です。刺激1つずつに注意を向けているので、刺激の数が増えると当然標的にたどり着くまでに経由する刺激の数は増えます。

図1で標的であるTだけが赤色だったら、どうでしょうか？セットサイズにかかわらず、すぐに見つかると思います。こういった探索を効率探索と呼びます。効率探索・非効率探索のサンプルデータを作って、下にプロットを描いてみました。

交互作用です。交互作用を含むモデルも、brmsで簡単に実装できます。モデル式はこういう風に書けます。\(\gamma_i search_i\)の部分が交互作用を表しています。\(\gamma_i = \beta_1 + \beta_3 x_i\)が傾きについての線形モデルになっていることがわかるでしょうか？

\[t_i\sim {\sf Normal}(\mu_i, \sigma)\]

\[\mu_i = \alpha + \gamma_i x_i + \beta_E search_i \]

\[\gamma_i = \beta_S + \beta_I search_i\]

dat2 <- list(
  N = NROW(dat_all2),
  ID = dat_all2$ID,
  RT = (dat_all2$RT - mean(dat_all2$RT)) / sd(dat_all2$RT),
  setsize = dat_all2$setsize,
  efficient = as.factor(dat_all2$efficient)
)

fit3 <- brms::brm(RT ~ 1 + setsize * efficient,
            family = "gaussian",
            data = dat2,
            iter = 3000,
            warmup = 1000,
            chain = 2,
            cores = 2,
            prior = c(prior(normal(0, 10), class = Intercept),
                      prior(normal(0, 10), class = b),
                      prior(cauchy(0, 10), class = sigma))
            )

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: RT ~ 1 + setsize * efficient 
   Data: dat2 (Number of observations: 20000) 
Samples: 2 chains, each with iter = 3000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept             -1.16      0.02    -1.20    -1.11 1.00     2386     2876
setsize                0.22      0.00     0.22     0.23 1.00     2441     2631
efficient1             0.60      0.03     0.54     0.67 1.00     1826     2166
setsize:efficient1    -0.22      0.00    -0.22    -0.21 1.00     1836     2087

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.77      0.00     0.76     0.77 1.00     3380     2388

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

練習の効果をモデリング

視覚探索課題に関わる現象はたくさんありますが、中でも私が注目しているのが練習や学習の効果です。ここでは一番基本的な学習の効果をモデリングしてみたいと思います。当然かもしれませんが、視覚探索課題を何度も何度も繰り返すと、どんどん反応時間が短くなっていきます。課題への慣れや注意誘導の効率性向上など、さまざまな理由が考えられるでしょう。だいたいこのような形で反応時間が変化します。

こういった形の関係性を、上で用いたような線形モデルでなんとかするのは難しそうです。ましてや、繰り返し回数ごとの平均反応時間を比較したところで、繰り返し回数によって反応時間がどのように変わっていくのかを調べるのは無理そうです。

この変化を表すのに使えそうなのがべき関数 (power function) です。

\[RT=a+bN^{-c}\]

こういうモデルになりそうです。

\[t_i\sim {\sf Normal}(\mu, \sigma)\]

\[\mu = a + bN^{-c} \]

\(a\)が練習前の反応時間 (ベースライン)、\(b\)が練習を通じてどれぐらい反応時間が短くなったのか、\(c\)は減り方 (大きいほど急速に減る) を示しているパラメータ、ととらえることができそうです。

今度はStanで書いてみました。

data {
  int<lower=1> N;
  vector[N] repN;
  vector[N] RT;
}

parameters {
  real a;
  real<lower=0> b;
  real<lower=0> c;
  real<lower=0> sigma;
}

model {
  vector[N] mu;
  for (n in 1:N) {
    mu[n] = a + b * repN[n]^(-c);
  }

  target += normal_lpdf(RT | mu, sigma);
  target += normal_lpdf(a | 0, 10);
  target += normal_lpdf(b | 0, 10);
  target += normal_lpdf(c | 0, 10);
  target += cauchy_lpdf(sigma | 0, 10) - cauchy_lccdf(0 | 0, 10);
}

推定結果はこんな感じになりました。

Inference for Stan model: 200219_power.
2 chains, each with iter=3000; warmup=1000; thin=1; 
post-warmup draws per chain=2000, total post-warmup draws=4000.

         mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
a       -0.67    0.00 0.17   -1.06   -0.76   -0.65   -0.56   -0.42  1159    1
b        3.03    0.00 0.20    2.67    2.89    3.02    3.15    3.41  1999    1
c        0.87    0.00 0.14    0.60    0.77    0.87    0.96    1.16  1485    1
sigma    0.77    0.00 0.03    0.71    0.75    0.77    0.78    0.82  2187    1
lp__  -472.25    0.04 1.44 -475.85 -472.98 -471.92 -471.20 -470.46  1449    1

Samples were drawn using NUTS(diag_e) at Wed Feb 19 10:27:47 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Stanコードを書いて推定した後には、うまくコードが走っているか・推定がうまくいっているかを確かめる必要があります。Stanの事後処理については、この資料が役立ちます。例としてトレースプロット、と呼ばれるプロットを載せます。今回は2本のチェーンを走らせているので、この2本の線がいい感じに重なっていたら、収束していそうだな、とわかります。

tidybayesパッケージを使って事後処理をすることもできます。詳しい説明はこちらで見ることができます。

draws <- fit_power %>%
  tidybayes::spread_draws(a, b, c)

head(draws)

# A tibble: 6 x 6
  .chain .iteration .draw      a     b     c
   <int>      <int> <int>  <dbl> <dbl> <dbl>
1      1          1     1 -0.483  2.95 1.13 
2      1          2     2 -0.645  3.14 0.850
3      1          3     3 -0.550  2.80 0.967
4      1          4     4 -0.698  2.74 0.732
5      1          5     5 -0.731  2.84 0.767
6      1          6     6 -0.780  2.92 0.745

点推定・区間推定値を出すことができます。tidyverseに慣れている人には随分使いやすいと思います。

draws %>%
  tidybayes::mean_hdci()

# A tibble: 1 x 12
       a a.lower a.upper     b b.lower b.upper     c c.lower c.upper .width
   <dbl>   <dbl>   <dbl> <dbl>   <dbl>   <dbl> <dbl>   <dbl>   <dbl>  <dbl>
1 -0.670  -0.976  -0.372  3.03    2.66    3.40 0.869   0.592    1.16   0.95
# ... with 2 more variables: .point <chr>, .interval <chr>

点推定値を使って、事後予測をプロットしてみました (曲線)。実際に得られたデータはマーカーで示してあります。

今回は1条件しか設けていないのですが、学習の速さや学習後の反応時間に影響しそうな現象は多くありますので、そういう操作を入れたときに各パラメータがどう変化するのかを見るのも楽しそうです。

似たようなことを私の研究でもやっていますので、興味のある方はご覧ください。視覚探索課題の正答率が処理時間の増加に伴って上がっていく様子をモデリングして、その増加率が学習の有無によって変化するかを調べました。

このモデルでは、入力 (繰り返し回数) を出力（反応時間）に変換する機能≒関数として、べき関数を利用しました。もちろん他の関数を試してみることもできますし、どちらの関数のあてはまりがよいのかを調べることもできます。また、実験条件を複数設けて、パラメタを比較することもできます。ただし、このようなモデルはあくまでも記述的なものであり、パラメタの解釈に注意が必要です (パラメタを心的過程に結びつけて解釈できない)。

おまけ: 階層モデル

ベイズ推定の特徴は、階層化が容易なことです。視覚探索課題のデータなどを見ていると気づくことですが、人によってベースの反応時間は結構違います。かなり速い人から、時間がかかる人までいます。非効率探索でも効率的に探索できる人もいれば、そうでない人もいます。探索の効率性に影響するような独立変数を操作したとき、その効果の出方にも個人差はあります。そこで、個人ごとの変動や試行間での変動を考慮してモデリングする必要が出てきます。そのときに登場するのが階層モデルです。階層モデルの利点を以下に示しました。

1. 反復測定の際に推定を向上する
2. 個人内/間・条件間・群間などの変動を推定できる
3. データを平均せずに利用することができる

brmsでは階層モデルが簡単にできます。Stanで書くには、少し訓練が必要かもしれません。単純な線形モデルを階層化してみます。自信がないので、コードが変だったら教えてください。

fit1_h <- brm(RT ~ 1 + setsize + (1|ID) + (1|setsize) + (1|ID:setsize),
            family = "gaussian",
            data = dat1,
            iter = 3000,
            warmup = 1000,
            chain = 2,
            cores = 2,
            prior = c(prior(normal(0, 10), class = Intercept),
                      prior(normal(0, 10), class = b),
                      prior(cauchy(0, 10), class = sigma))
            )

lmerTestなどで線形混合モデルを使ったことがある人には、見慣れた書き方だと思います。

summary(fit1_h)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: RT ~ 1 + setsize + (1 | ID) + (1 | setsize) + (1 | ID:setsize) 
   Data: dat1 (Number of observations: 10000) 
Samples: 2 chains, each with iter = 3000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~ID (Number of levels: 5) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.27      0.36     0.01     1.16 1.01       99      226

~ID:setsize (Number of levels: 10) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.55      0.18     0.31     1.00 1.01      205      274

~setsize (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     4.72      4.05     0.20    14.95 1.04      107      179

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    -4.03      8.43   -23.24    11.47 1.03      121      148
setsize       0.51      1.04    -1.47     2.79 1.03      119      217

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.72      0.01     0.71     0.73 1.00      266      749

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

個人ごとの変動をプロットすることもできます。

この資料で扱ったような他のモデルも、すべて階層化できます。