2020-12-15

複数の.csvファイル（じゃなくてもいいけど）を一気に結合する

本記事について

Online Psychological Experiment Advent Calendar 2020 第17日目の記事です。

最初は気軽にポエムでも書こうかと思っていたものの、他の記事がどれも真面目だったので、「だから、オレだってなんかしなくっちゃあな…カッコ悪くて2021年に行けねーぜ…」ということで急きょtipsを書くことにしました。

オンライン実験のデータがどのように保存されるか

もちろん、どんな実験をするかや、どのサービスを利用するかによって、話は変わります。しかし一般的には、「協力してくださった実験参加者の数だけ、ファイルが蓄積される」と思います。オフラインの（実験室で行うような）実験ならなおさらですね。

例えば僕自身は、こちらの記事を参考にして、lab.jsで作成した.jsonファイルを、Open Labにアップロードし、単語記憶課題のオンライン実験を行いました。その結果、実験参加者の数だけ個別の.csvファイルがサーバ上に蓄積されました。

複数の.csvファイル（じゃなくてもいいけど）を一気に結合したい

もちろん.csvファイルに限らず、全てのファイルが.xlsxでも良いのですが、ともかく拡張子が同じ複数のファイルを、1つのファイルに結合しなければ、前処理も分析もできないわけです。

実験参加者数が少なければ（例えば10人など）、1つのファイルを開いてデータをコピーして、別のファイルにペーストするのを繰り返す...という手作業でも、さほど手間にはならないでしょう（※言うまでもなく、これはヒューマンエラーが混入する恐れがあるので、避けるべきですが）。

しかしオンライン実験となると、一般的にサンプルサイズが大きくなりやすいと思うので、手作業は労力の観点から現実的ではありません（※くどいようですが、労力以外の観点からも、データをコピペでまとめるのは避けるべきです）。

purrrパッケージを使おう

R環境限定の話になりますが、purrrパッケージを用いれば、同じディレクトリ内に存在する、拡張子が同じ全てのファイルを、1つのデータフレームに結合することが、容易にできます。

その方法は...これだっ！（他の方が作成された資料を共有していくツェペリ魂）

上の画像のコードは、tidyverseパッケージをロードする前提で書かれているので、最低限必要なpurrrパッケージとreadrパッケージだけをロードした場合には、以下の書き方になります。
（tidyverseパッケージをロードしたら、これら2つのパッケージも自動的にロードされますが）

library(purrr)
library(readr)


target_files = list.files("data/",
                          pattern = ".csv$",
                          full.names = TRUE)

dat = purrr::map_df(target_files, readr::read_csv)

まず、特定のディレクトリ内に存在する、拡張子が.csvのファイル名を全て、target_filesという名前のリストにまとめています。この場合は、現在の作業ディレクトリの中に、「data」という名前のフォルダがあり、その中に全ファイルが存在する想定です。

次に、 readr::read_csv()関数で個々の.csvファイルを読み込むたび、purrr::map_df()で1つのデータフレームとして結合していきます。結合されたデータフレームは、datという名前のオブジェクトに格納されます。

自分用にカスタマイズも可能

データによっては、

ファイルを読み込む際に、空欄のセルをNAにしたい
ファイル名も、データフレーム内の変数として保存したい

などの、細かいカスタムが必要となることもあるでしょう。そのような場合には、データを読み込む関数（上のコードではreadr::read_csv()）をカスタムした自作関数を定義すればよいです。

ここで新たにdplyrパッケージをロードしているので、パイプ演算子%>%を使用しています（それだったらtidyverseパッケージを1つロードすればよい気もしなくもない）。

library(purrr)
library(readr)
library(dplyr)

# 自作関数
read_func = function(x){readr::read_csv(file = x, na = "") %>%
    dplyr::mutate(filename = x)}


target_files = list.files("data/",
                          pattern = ".csv$",
                          full.names = TRUE)

dat = purrr::map_df(target_files, read_func)

このコードでは、read_func()という関数を自作しています。
まずreadr::read_csv()でファイルを読み込む際に、空白のセルをNAにするよう、引数を指定しています。さらに、1つのファイルをそのようにして読み込んだあとで、dplyr::mutate()によって、ファイル名をfilenameという名前の列に追記しています。
あとはpurrr::map_df()のなかで、自作した関数read_func()を使用することを明記するだけです。

おわり

全てのデータが、1つのデータフレームにまとめられれば、第15日目の記事のように前処理を行うことも、分析することもできます。

Enjoy !!

2020-12-11

Google Colaboratory + PyStanを使ってみた話

Stan Python

この記事について

Stan Advent Calendar 2020 第14日目の記事です。
14日目のカレンダーが空いていたので、何か埋めなきゃということで、「そういえば使ったことないから、Python環境でStanを使ってみよう」という思い付きで書いています。

普段はR環境でStanを使っているので、全てが手探りです。そんなわけで、この記事は自分のための備忘録です。世の中にはもっと体系的にまとまった記事がたくさんあるので、ぜひそれらを参考にしてください。

今回はGoogle Colaboratory上で、PyStanを用いることにします。なお今回はGoogle Driveのマウントはしませんが、Google ColaboratoryをGoogle Driveと連携させることで、より使いやすくなると思います。

実演

.stanファイルを用意する

これはRやPythonとは関係のない作業なので*1、好きなエディタを用いて書きます。ここでは説明変数が1つの、単回帰モデルを例とします。

$y_i=\beta_0+\beta_1x_i+\varepsilon_i \qquad i = 1, 2, ..., n \\ \varepsilon_i \sim {\rm Normal}(0, \sigma)$

$n$ 個のデータがあり、目的変数 $y$ を、説明変数 $x_1$ で予測しています。 $\beta_0$ が切片、 $\beta_1$ が回帰係数、 $\varepsilon$ が残差です。
これは以下の確率モデルと同義です。

$y_i \sim {\rm Normal}(\beta_0+\beta_1x_i, \sigma) \qquad i = 1, 2, ..., n \\$

つまり、

目的変数 $y_i$ は、正規分布に従う。その平均は $\beta_0+\beta_1x_i$ で、標準偏差が $\sigma$ である

というモデルを立てたことになります。以下のStanコードを、regression.stanという名前で保存しておきます。

data {
  int<lower=0> n_data;   //サンプルサイズ
  real y[n_data];        //目的変数
  vector[n_data]  x;     //説明変数
}

parameters {
  real beta_0;           //切片
  real beta_1;           //回帰係数
  real<lower=0> sigma;   //標準偏差
}

model {
  y ~ normal(beta_0 + beta_1*x, sigma);
}

Google Colaboratory上に.stanファイルをアップロードする

自分のPCにPython環境がなくても、インターネットブラウザ上でPythonコードが実行できる、Google Colaboratoryを使用します。ここにアクセスして、[ファイル　→　ノートブックを新規作成]します。

f:id:das_Kino:20201211103021p:plain

ローカルの.stanファイルは、以下の手順で容易にアップロードできます。 f:id:das_Kino:20201211103808p:plain

「注: アップロードしたファイルはランタイムのリサイクル時に削除されます。」というメッセージが現れると思いますが、これはGoogle Colaboratoryを使う以上仕方のないことです。

アップロードされた場所は、対象のファイル名にカーソルを合わせて、右端の「･･･」をクリックすると現れるメニューのうち、「パスをコピー」を選択するとクリップボードに格納されるので、知ることができます。実際にやってみると、/content/regression.stanであることが分かります。

f:id:das_Kino:20201211104119p:plain

ライブラリの読み込み

最低限、これらのライブラリを読み込めば良いと思います。PyStanのサンプリング結果を、トレースプロットやチェインごとの事後分布といった形で可視化するためには、arvizというライブラリが適しているようですが、Google Colaboratoryではまず!pip install arvizでインストールする必要があります。

import numpy as np #サンプルデータを作成する際に、randomモジュールで乱数生成するため
import pystan #stanのPython用インタフェース

# 可視化用ライブラリ ----------------------
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns

# PyStanのサンプリング結果の可視化用ライブラリ ------------------
!pip install arviz # arvizというライブラリをまずインストールする
import arviz

サンプルデータの作成

単回帰モデルを真のモデルとして、パラメータの真値を指定してサンプルデータを作成します。

n = 50 #サンプルサイズ

np.random.seed(seed = 1234) #乱数のシード
x = np.random.randint(low = 1, high = 30, size = n) #説明変数。1～29の一様乱数（整数）をn個生成
y = np.random.normal(loc = 20 + 1.5 * x,
                     scale = 3.0,
                     size = n) #目的変数。正規分布に従い、その平均は20 + 1.5x。標準偏差は3.0

可視化してみます。

# サンプルデータの可視化
sns.jointplot(x = x, y = y)

f:id:das_Kino:20201211110044p:plain

.stanファイルのコンパイルとサンプリング

まずはアップロード済みの.stanファイルをコンパイルします。

sm = pystan.StanModel(file = '/content/regression.stan')

次に、サンプリングです。Rだと、rstan::sampling()関数の中で、object = smとコンパイルしたオブジェクトを指定しますが、Pythonだとちょっと書き方が違うことに注意です。

また、Rでは、rstan::sampling()関数の中で、list型でデータを渡しますが、PyStanでは辞書型でデータを渡す必要があるところも違いますね。辞書型のデータはdict()などで作成できます。

その他の引数はRStanでもPyStanでもほとんど同じですね。

fit = sm.sampling(
    data = dict(
        n_data = n,
        y = y,
        x = x
        ),
    seed = 1234,
    iter = 2000,
    warmup = 1000,
    chains = 4)

収束診断とサンプリング結果の出力

まずはうまく収束したかどうか、視覚的に確認してみましょう。ここではarvizライブラリのplot_trace()を用いて、トレースプロットやチェインごとの事後分布を可視化します。

arviz.plot_trace(fit)

f:id:das_Kino:20201211111419p:plain

収束していると考えてよさそうです。
要約統計量は以下の通りです。うまくパラメータリカバリできていますね。

print(fit)

f:id:das_Kino:20201211111202p:plain

おわり

PyStan、初めて使ってみたんですが、Google Colaboratoryを併用することで、環境構築の手間もほとんどいらず、手軽に利用することが出来ました。

Enjoy !!

*1:もっとも、RStudioの文法チェック機能が優れているので、慣れないうちはRStudio上で書くのが良いような気はします

2020-12-09

Stan Advent Boot Camp 第9日目: 階層線形モデルをやってみよう

R Stan

この記事について

Stan Advent Calendar 2020 9日目の記事です。本記事で紹介する階層線形モデルは、第4日目に紹介された重回帰モデルや、第5日目に紹介されたロジスティック回帰モデルなどの、一般線形モデルや一般化線形モデルを拡張したモデルになります。Stanで階層線形モデルの確率モデルを書く方法を紹介します。

なおこの記事では階層線形モデルという名称を使用しますが、線形混合モデルやマルチレベルモデルとも本質的には同じです。

一般線形モデルの復習

ここでは説明変数が1つの、単回帰モデルを例とします。

$y_i=\beta_0+\beta_1x_i+\varepsilon_i \qquad i = 1, 2, ..., n \\ \varepsilon_i \sim {\rm Normal}(0, \sigma)$

$y_i \sim {\rm Normal}(\beta_0+\beta_1x_i, \sigma) \qquad i = 1, 2, ..., n \\$

つまり、

目的変数 $y_i$ は、正規分布に従う。その平均は $\beta_0+\beta_1x_i$ で、標準偏差が $\sigma$ である

というモデルを立てたことになります。

階層線形モデルへの拡張

データの階層性を考慮すべき状況

StanとRでベイズ統計モデリングの121ページから引用すると、階層モデルとは

説明変数だけでは説明がつかない、グループに由来する差異（グループ差）を上手く扱うための一手法

です。
「グループ」というのは例えば、

実験参加者ID	グループ
1	◇◇大学
2	◇◇大学
3	◇◇大学
4	▲▲大学
5	▲▲大学
6	▲▲大学
7	☆☆大学
8	☆☆大学
9	☆☆大学

などです。グループは上表のような「所属集団」に限らず、以下のように個人をグループと見なすことも可能です。下表では、「実験参加者ID」がグループに相当しています。

実験参加者ID	試行番号
1	1
1	2
1	3
2	1
2	2
2	3
3	1
3	2
3	3

上の2つの表ではいずれも、9つのデータが得られていることになっています。では、目的変数 $y_i$ と説明変数 $x_i$ も同様に9つずつデータが得られているとき、

$y_i \sim {\rm Normal}(\beta_0+\beta_1x_{i}, \sigma) \qquad i = 1, 2, ..., n \\$

という回帰モデルを考えてよいかというと...そうとは限りません。ここで、データの階層性を考慮しなければならない場合があります。その理由は、以下のスライドの12枚目に詳しいです。

マルチレベルモデル講習会理論編 from Hiroshi Shimizu

階層線形モデルの確率モデル

上記の回帰モデルを階層線形モデルに拡張したい場合、色々な拡張の方法がありますが、典型的には以下の3通りが考えられます。

切片のみに、グループ差を考慮する
回帰係数のみに、グループ差を考慮する
切片と回帰係数に、グループ差を考慮する

ここでは一例として、「3. 切片と回帰係数に、グループ差を考慮する」場合の階層線形モデルの確率モデルを紹介します。

$y_{i} \sim {\rm Normal}(\beta_{0j\lbrack i \rbrack}+\beta_{1j\lbrack i \rbrack}x_{i}, \sigma) \\ \beta_{0j} \sim {\rm Normal}(\mu_0, \sigma_0) \\ \beta_{1j} \sim {\rm Normal}(\mu_1, \sigma_1)$

添え字 $i$ は一つひとつのデータを、添え字 $j$ は一つひとつのグループを表します。上の表でいえば（一つ目の表）、添え字 $i$ は一人一人の実験参加者を、添え字 $j$ は各大学を表しています。切片パラメータ $\beta_0$ と回帰係数パラメータ $\beta_1$ に添え字 $j$ が付いて、 $\beta_{0j}$ と $\beta_{1j}$ になっているので、グループごとにこれらのパラメータが異なることを仮定しています。

では、グループごとに異なると仮定した切片と回帰係数は、どのように、どれくらいグループごとに異なるのでしょうか。それを反映しているのが、

$\beta_{0j} \sim {\rm Normal}(\mu_0, \sigma_0) \\ \beta_{1j} \sim {\rm Normal}(\mu_1, \sigma_1)$

です。つまりグループごとの切片と回帰係数の違いは、正規分布に従うというモデルを立てたことになります。

実演

まずは、パラメータの真値と真のモデルが分かっている状況からサンプルデータを生成します。

library(ggplot2)
library(rstan)

set.seed(1234)

n_data = 100 #総データ数
n_group = 5 #グループ数
n_each_group = n_data / n_group　#各グループ内の人数 = 20

mu_0 = rnorm(n = n_group, mean = 30, sd = 10) #切片は、平均30、標準偏差10の正規分布に従う
mu_1 = rnorm(n = n_group, mean = 3, sd = 1.5) #回帰係数、平均3、標準偏差1.5の正規分布に従う

x = runif(n = n_data, min = 1, max = 30) #説明変数
y = rnorm(n = n_data, mean = (mu_0 + mu_1*x), sd = 15) #目的変数
group = rep(1:n_group, time = n_each_group) #グループを識別する変数

このようなサンプルデータの作り方は、StanとRでベイズ統計モデリングに載っています。

グループごとに回帰直線を描くと、こんな感じです。確かに切片も傾きも、グループごとに異なっています。

ggplot(data = data.frame(x, y, group), mapping = aes(x = x, y = y, color = factor(group))) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)

f:id:das_Kino:20201208170634p:plain

この確率モデルを指定してパラメータを推定したら、真値に近くなることが期待できるはずです。やってみましょう。Stanコードは以下です。

方法1：回帰モデルを拡張した書き方

data {
  int<lower=0> n_data;      //データ数
  int<lower=0> n_group;     //グループ数
  
  real y[n_data];           //目的変数
  vector[n_data] x;         //説明変数
  
  int group[n_data];        //グループを識別する変数
}

parameters {
  vector[n_group] beta[2];  //beta[1]が各グループの切片
                            //beta[2]が各グループの回帰係数
  
  real mu[2];               //mu[1]が各グループの切片の平均
                            //mu[2]が各グループの回帰係数の平均
                            
  real<lower=0> sigma[3];   //sigma[1]が各グループの切片の標準偏差
                            //sigma[2]が各グループの回帰係数の標準偏差
                            //sigma[3]が目的変数が従う正規分布の標準偏差
}

model {
  beta[1] ~ normal(mu[1], sigma[1]);
  beta[2] ~ normal(mu[2], sigma[2]);
  
  y ~ normal(beta[1, group] + beta[2, group].*x, sigma[3]);
}

.*という演算子が用いられていることに注意してください。これは同じ長さのベクトルを、要素ごとに掛算するための演算子です。

これをmixed_1.stanという名前で保存して、以下のRコードからサンプリングを行います。

sm_1 = rstan::stan_model("mixed_1.stan")

fit_1 = rstan::sampling(object = sm_1,
                        data = list(
                          n_data = n_data,
                          n_group = n_group,
                          y = y,
                          x = x,
                          group = group
                        ),
                        seed = 1234,
                        iter = 4000,
                        warmup = 1000)

print(fit_1)

トレースプロットやchainごとの事後分布を視覚的にチェックした限りでは、収束していると判断して良さそうです。

f:id:das_Kino:20201209102951p:plain

事後分布の要約統計量は以下の通りです。

f:id:das_Kino:20201209102236p:plain

グループごとの切片が従う正規分布の平均 $\mu[1$ ]は、真値が30のところ、事後中央値は28.27
グループごとの切片が従う正規分布の標準偏差 $\sigma[1$ ]は、真値が10のところ、事後中央値は18.26
グループごとの回帰係数が従う正規分布の平均 $\mu[2$ ]は、真値が3のところ、事後中央値は2.86
グループごとの回帰係数が従う正規分布の標準偏差 $\sigma[2$ ]は、真値が1.5のところ、事後中央値は1.62
目的変数が従う正規分布の標準偏差 $\sigma[3$ ]は、真値が15のところ、事後中央値は15.22

グループごとの切片が従う正規分布の標準偏差 $\sigma[1$ ]だけちょっと真値と事後中央値が離れていますが、おおむね良くパラメータリカバリできているように思います。

方法2：多変量正規分布を用いた書き方

上記の階層モデルは、多変量正規分布を用いて、以下のように書くこともできます。多変量正規分布は、Stan Advent Calendar 2020 第6日目の記事でも紹介されています。

※前述の確率モデルに対して、違うStanコードの書き方ができるということではなく、確率モデルも多変量正規分布を用いた別のものに変わっていることに注意してください

data {
  int<lower=0> n_data;             //データ数
  int<lower=0> n_group;            //グループ数
  int<lower=0> n_slope;            //グループ差を考慮する回帰係数の数（切片含む）
  
  real y[n_data];                  //目的変数
  matrix[n_data, n_slope] X;       //説明変数（切片含む）の行列
  
  int group[n_data];               //グループを識別する変数
}

parameters {
  vector[n_slope] beta[n_group];   //グループごとの、回帰係数（切片含む）のベクトル
  
  real<lower=0> sigma;             //目的変数が従う正規分布の標準偏差
  vector[n_slope] mu;              //各回帰係数（切片含む）の平均ベクトル

  cov_matrix[n_slope] Cov;         //多変量正規分布の分散共分散行列
}

model {
  for(n in 1:n_data){
    y[n] ~ normal(X[n]*beta[group[n]], sigma);
  }
  
  beta ~ multi_normal(mu, Cov);
}

modelブロックのy[n] ~ normal(X[n]*beta[group[n]], sigma)は、行列の掛算によって、目的変数 $y$ が従う正規分布の平均を計算しています。重回帰モデルにおける同様の書き方は、Stan Advent Calendar2020の4日目の記事でも紹介されています。

ポイントは、各回帰係数（切片を含む）がvector[n_slope] beta[n_group]で宣言されているところです。回帰係数の数（切片を含む）の要素を持つbetaというパラメータが、グループの数（n_group）だけあると宣言しています。

modelブロックのbeta ~ multi_normal(mu, Cov)は、各回帰係数（切片を含む）betaが、多変量正規分布に従うことを仮定しています。多変量正規分布のパラメータは、以下の2つです。

平均ベクトルmu
分散共分散行列Cov

要するに、それぞれの回帰係数は同時に正規分布に従っており、かつ互いに相関があるというモデルになっています。

ただし多変量正規分布の分散共分散行列を、quad_form_diag()というStanの関数を用いて以下のように定義すると、サンプリング効率がよくなるようです。

data {
  int<lower=0> n_data;             //データ数
  int<lower=0> n_group;            //グループ数
  int<lower=0> n_slope;            //グループ差を考慮する回帰係数の数（切片含む）
  
  real y[n_data];                  //目的変数
  matrix[n_data, n_slope] X;       //説明変数（切片含む）の行列
  
  int group[n_data];               //グループを識別する変数
}

parameters {
  vector[n_slope] beta[n_group];   //グループごとの、回帰係数（切片含む）のベクトル
  
  real<lower=0> sigma;             //目的変数が従う正規分布の標準偏差
  vector[n_slope] mu;              //各回帰係数（切片含む）の平均ベクトル
  
  
  corr_matrix[n_slope] Omega;      //それぞれの正規分布同士の相関行列
  vector<lower=0>[n_slope] tau;    //それぞれの正規分布の標準偏差
}

model {
  for(n in 1:n_data){
    y[n] ~ normal(X[n]*beta[group[n]], sigma);
  }
  
  beta ~ multi_normal(mu, quad_form_diag(Omega, tau));
  
  Omega ~ lkj_corr(2);             //相関行列の弱情報事前分布。
}

Stan モデリング言語: ユーザーガイド・リファレンスマニュアル（日本語翻訳版）によればquad_form_diag()という関数は、

quad_form_diag(Sigma,tau)がdiag_matrix(tau) * Sigma * diag_matrix(tau)と等価になるように定義されています. ここで, diag_matrix_(tau)は, 対角成分がtauとなり, それ以外が0の行列を返します

という働きを持ちます（上記のStanコードでは、SigmaというパラメータをOmegaという名前で宣言しています）。

これで本当に分散共分散行列が作れることを確かめてみましょう。今回のように、グループ差を仮定している回帰係数が2つの場合（切片と、説明変数の回帰係数の2つ）、以下のような計算になります。

f:id:das_Kino:20201209121459p:plain

対角成分は、標準偏差に相当する $\tau$ の2乗なので、分散を表します。非対角成分は、標準偏差 $\tau$ の積に、相関係数 $r$ を掛けているので、共分散を表します。確かに、分散共分散行列が作れていますね。

それでは上記のStanコードをmixed_2.stanという名前で保存し、以下のRコードを実行してみましょう。使用しているサンプルデータは、mixed_1.stanに渡したものと同じです。
model.matrix()は、引数にformula（回帰式）を指定すると、説明変数を行列形式でまとめてくれる関数です。今回は切片も回帰係数とみなしているので、全ての値が1の説明変数を追加する必要があるため、X = model.matrix(y ~ 1 + x)という書き方にしています（ただしX = model.matrix(y ~ x)でも同じです）。

X = model.matrix(y ~ 1 + x)  # = の左側は大文字のX、右側は小文字のxであることに注意

fit_2 = rstan::sampling(object = sm_2,
                        data = list(
                          n_data = n_data,
                          n_group = n_group,
                          n_slope = ncol(X),
                          y = y,
                          X = X, # = の左右ともに、大文字のXであることに注意
                          group = group
                        ),
                        seed = 1234,
                        iter = 4000,
                        warmup = 1000)

パラメータ数が多いので、トレースプロットやチェインごとの事後分布は掲載しませんが、収束していると判断して良さそうです。事後分布の要約統計量は以下の通りです。

f:id:das_Kino:20201209122814p:plain

グループごとの切片が従う正規分布の平均 $\mu[1$ ]は、真値が30のところ、事後中央値は28.15
グループごとの切片が従う正規分布の標準偏差 $\tau[1$ ]は、真値が10のところ、事後中央値は19.41
グループごとの回帰係数が従う正規分布の平均 $\mu[2$ ]は、真値が3のところ、事後中央値は2.85
グループごとの回帰係数が従う正規分布の標準偏差 $\tau[2$ ]は、真値が1.5のところ、事後中央値は1.67
目的変数が従う正規分布の標準偏差 $\sigma$ は、真値が15のところ、事後中央値は15.20

相変わらず、グループごとの切片が従う正規分布の標準偏差 $\tau[1$ ]だけちょっと真値と事後中央値が離れていますが、おおむね良くパラメータリカバリできていますね。

まとめ

~~いかがだったでしょうか~~
パラメータに階層性を仮定したモデルは、ちょっとだけ書き方が難しくなりますが、基本的には回帰モデルの拡張と考えられます。実際のデータ分析では、階層性を考慮したほうがよい状況も多いので、参考になれば幸いです。

Enjoy!

2020-12-04

昔はHarukara.Rというのがあったんだが、膝に矢を受けてしまってな

雑談

この記事について

ベイズ塾 Advent Calendar 2020 5日目の記事です。このブログは技術ブログ・書評が中心なんですが、今回の記事はポエムです。

広島ベイズ塾は、創設メンバーの一人、某氏（ここでは仮にA氏とする）が広島にいらっしゃったときにできました。詳しい話は、もしかしたら最終日（2020/12/25）に紹介されるのかもしれません。

磯野ー、カラオケしようぜー

しばらくして、A氏が関西に異動されました。カラオケに誘われました。

A氏「カラオケの前に寿司でも食うか」
僕「おっ、いいですねー。あ、でもどこも混んでますねー。ググったら、近くに1件だけ寿司屋ありました、もうそこにしますか」
A氏「せやな」

繁華街の片隅にひっそりと佇むその店内には、ショーケースに色んなネタの切り身が陳列されていた。ネタの名前も、価格も表示されていない。名もなき詩魚（うお）。

僕「（ひそひそ、Aさん、これ、何のネタやと思います？）」
A氏「（ひそひそ、わからへん。赤いな）」

ええい、残さずに全部食べてやる oh darlin 君は誰

そんなこんなで僕らの間には、「寿司を食べてからカラオケに行く」という謎のルーティーンが出来あがったのでした。僕らがよく行っていた寿司屋の名前（はる〇〇）と、カラオケを合わせて、Harukaraの誕生です。二人組音楽ユニットじゃありません~~（いやそれHALCALI）~~

Rの勉強会でもする？

時を同じくして広島では、R旋風が巻き起こっていました。その中心はHijiyama.RというRの勉強会（僕らも準レギュラーでしたけど）。

関西でもやる？　どちらからともなくそう切り出したのは、偶然ではなく必然だったのかもしれない。

そんなわけで、Rの勉強会をしてから寿司を食べてカラオケに行く、Harukara.Rが誕生しました。

まあRの勉強会といいつつ、その実態は統計の勉強会でした。Rを使って計算をしているんだから間違ってない。ここで階層線形モデルとかベイズ（入門編）とか色々学んでゆきました。広島から参加者が来ることもあって、しかもみんな寿司が好きで歌が上手いので、いつしか常連メンバーが増えていきました。

いいぞ、Harukara.R、どうしてそんなに大きくなってしまったんだい？　真面目にやってきたからよ、ね？~~（アリさんマークのひっこｓ）~~

アヒル本の読書会でもする？

そんなこんなで2016年秋。いつものように寿司を食べていた僕らは、「次のHarukara.Rでは何の勉強する？」と雑談していました。

聞くところによると、StanとRでベイズ統計モデリングという本（のちに、アヒル本という通称が生まれる）が発売になるらしい。しかも著者は、A氏が愛読しているブログの著者らしい。

じゃあもう、この本の読書会するっきゃないでしょ。というわけで当初はHarukara.Rのいち企画として読書会するつもりだったんですが、せっかくだから詳しい人や興味を持っている人もたくさん巻き込んでやろう、ということで、より大きな組織に企画をお任せすることになり、そこで生まれたのがOsaka.Stanでした。

あれよあれよという間に、色んな方にご参集いただくことができ、ついにはアヒル本の著者にまでゲストトークしていただくことになり、ここで多くの出会いが生まれ、色んなことを学びました。

なんかもう、どこに行っても同じ人いるよね

気づけば、色んな勉強会が乱立していたわけですが、どこに行っても同じ人と会うわけです。そんなこんなでいつしか、広島ベイズ塾に色んな人が集まるようになって、今に至るという。

色んな人とつながり、色んな方向から学ばせていただき、一人では知り得なかったことを知ってゆく。人と人を繋ぐ、キューピットの矢を膝に受けて、走り出していくんだ、遠い地平へーーー

（完）

2020-08-18

【書評】PsychoPyでつくる心理学実験

書評

本記事について

朝倉書店様と訳者の先生方より、2020年9月1日出版の、PsychoPyでつくる心理学実験をご恵投いただきました！

www.asakura.co.jp

PsychoPyの作者自身が筆者であるため、本当に素晴らしい教科書でした。まずPsychoPyについて簡単に紹介したのち、書評を記します。

※ただ書いてから気づきましたが、書評というより要約みたいになってしまいました。まあいいか。

そもそもPsychoPyとは？

University of NottinghamのJonathan Peirce氏が作成された、「心理学、神経科学、言語学のための多彩でダイナミックな実験を作成するための、オープンソース（無料）のソフトウェアパッケージ」です（本書より直接引用）。

端的にいえば、「心理学、神経科学、言語学でよく用いられる実験課題を無料で作成できるソフト」です。Pythonをベースとしているので、名前に「Py」が付いています。

PsychoPyのダウンロードはこちらから。

PsychoPyの作者は海外の研究者ですが、愛媛大学・十河宏行先生が日本語化や、日本語ドキュメントの執筆を行ってくださっています。

インストール方法や基本的な使用方法は、関西学院大学・小川洋和先生の解説が最も分かりやすいと思います。

これまでにも同様の目的を実現する商用ソフトはありました。しかしPsychoPyは、それらと遜色ない、いやむしろ今やそれら以上に充実した機能を備えたソフトであり、しかもタダ。作者自身が実験心理学者であるため、痒い所にちゃんと手が届く仕様。

特徴的なのは、以下の3通りの方法で実験課題を作成できること。

Builder（プログラミングのコードを一切書く必要がないGUIで、マウス操作が主体）
Coder（Pythonのコードを書く）
BuilderとCoderのハイブリッド（大半はBuilderで作成し、一部の機能をCoderで実現する）

このうち、Builderの操作が本当に分かりやすく、それでいて多機能なので、学部生向けの実験実習から、プロ研究者の実験まで、幅広く用いることができます。僕も学部生向けの実験実習で長年使用していますが、上記の関西学院大学・小川洋和先生の解説に沿って実演するだけで、30分程度で基本的な機能は習得してもらえます。

僕自身でも、長年にわたりBuilderで実験課題を作成してきましたが、最近はちょっと複雑な実験課題を組むことがあるため、柔軟性の高いCoderに移行しています。PsychoPyはPythonベースのソフトと書きましたが、上述のように「パッケージ」なので、心理実験等でよく組み込まれる処理（例：刺激を200msec呈示）を行うためのPythonコードが、関数として多数用意されています。

十河先生が執筆された「心理学実験プログラミング: Python/PsychoPyによる実験作成・データ処理」を丁寧に写経すれば、Python初心者でも十分にCoderで実験課題を作成できます。

www.asakura.co.jp

この記事も非常におすすめです。

qiita.com

書評

ではここから書評を記します。僕が本書の特長だと思った点は、以下の4点です。

Builder操作に関する解説が主体
いくつかの代表的な実験課題（例：視覚探索課題）を例に挙げながら、機能を解説
実験心理学の教科書（←これが良い！）
外部機器との連携

1つずつ説明していきます。

1. Builder操作に関する解説

上で、「Coderのほうが柔軟性が高いから、僕自身では、最近はCoderを使っている」と書きました。え、じゃあBuilderって、GUIで簡単に操作できるとはいえ、できることが限られるんじゃないの？　Coderの解説書の方が需要が高いんじゃないの？　と思うかもしれません。

その疑問に対する答えは、「Yes、でもそれでいいっ！　それがBEST！」です（リサリサ先生の声で）。

Builderの設計思想は、GUIで任意の実験課題を作れるようにすることではありません。Builderは、GUIで「標準的な実験課題」を作れるように用意されたもの。Builderには多くのメニューが用意されていますが、意図的に機能を網羅していません。なぜなら、多すぎるメニューは利便性を低下させるし、何より初心者を混乱させるからです。幅広いユーザーが、簡便に標準的な実験課題を作れるように調整されたインタフェースが、Builderなのです。

したがって、初心者は難なく使用方法を習得できるし、プロ研究者であっても標準的な実験課題を用いる場合にはBuilderを使って何の問題もありません。むしろ素早く実験を作成できるのだから、その方がいいじゃないか。デバグも楽になるし。

あ、熱く語りましたが、これは本書に記載されていることです。僕はここを読んだとき、ハートが震えてヒートが燃えつきました（いや燃えつきたらあかん）。PsychoPyの作者自身が執筆した書籍なので、こういう作者の思いが垣間見えるのがいいですね。

ちなみに、後述もしますが、PsychoPyはオンライン実験にも対応しています。PsychoPyで作成した実験課題を専用のサーバにアップロードすることで、インターネットブラウザを介して実験を受けられるようになるということです。本記事の執筆時点（2020/8/18）では、オンライン実験用のプログラムは、原則としてCoderを使用せず、Builderで作成しなければならないと思います。そういう意味でも、Builderで課題を作成する方法の解説には意義があります。

※なお本書では、オンライン実験のための手続きは掲載されていません

2. 代表的な実験課題（例：視覚探索課題）を例に挙げながら、機能を解説

PsychoPyが見据えている研究領域（心理学、神経科学、言語学）の特徴は、「こういうデータを測定したかったら、こういう実験がよく用いられるよね」という実験パラダイムが存在することかと思います。

例えば抑制能力（特定の情報処理を我慢して、別の情報処理を行う）を測定したかったら、ストループ課題が例に挙がると思います。ストループ課題とは、「あか」「あお」のように、単語の意味とインクの色が一致する刺激（あか）と、一致しない刺激（あお）を呈示し、インクの色を回答させるという実験です（「あか」でも「あお」でも、インクの色は赤）。

あるいは、特定の情報を探し出す速さを測定したかったら、視覚探索課題が用いられると思います。例えば画面上に多くの「L」が散りばめられており、その中に「T」が1つだけ紛れ込んでいるか否かを判断するという実験です。

本書では、このような実験課題を例に挙げながら、Builderの機能を広く解説しています。汎用的な機能も学びつつ、特定の実験課題を作り上げられるので、一石二鳥ですね。

しかも、冒頭で書いたような、「作者自身が実験心理学者であるため、痒い所にちゃんと手が届く仕様」の紹介も随所に存在しています。

実験では、刺激を呈示して、それに対する被験者の反応を測定することが基本ですが、それだけすれば十分というわけではありません。ときに、測定精度を高めるための”工夫”を交えることがあります。

例えば、通常は本番の実験の前に練習を行い、被験者に実験手続きを理解していただくことが必要になります。被験者が実験手続きを理解していなかったら、本番中に測定されるデータ（例：反応の速度）は、どのような情報処理を反映しているのか分かりにくくなるからです（例：反応時間が長くても、それは刺激が処理しにくかったのか、反応方法が分からなかったのか、区別できない）。

そこで、例えば練習試行の最後の5試行ぶんだけ抽出し、反応時間や反応の正確性を確認したいというモチベーションが生まれます。そこで極端に反応が遅かったり、正答率が低かったりしなければ、被験者が実験手続きを理解したと推測することができます（もちろん、これだけで十分という意味ではありませんが）。

じゃあPsychoPyでどう実装するのか。ちゃんと実例が紹介されているんですね、これが。

このように、基本的な機能の紹介にとどまらず、「実践的テクニック」も随所に盛り込まれているのが、素晴らしい点だと思います。

3. 実験心理学の教科書

さて、僕が何より感動したのが、この点です。

ここまで、「PsychoPyを使えば簡単に実験課題が作れるよ！　実例も載ってるよ！」と、”楽ができる”ことを推してきました。しかし一般的に、「実験の準備」という手続き自体は、楽ではありません。

実験を実施する前に、丁寧に時間をかけて確認しなければならないこと、調整しなければならないことが、たくさんあります。そのために知らなければならないことが、たくさんあります。

PsychoPyは、「作る」という作業においてユーザーに楽を与えてくれます。しかし、それは諸手続きを軽視してよいことを意味しません。

例えば、現在ではPCの液晶ディスプレイを用いて実験刺激の呈示を行うことが多いと思います。一般的には、画面は60Hzで更新されることが多いでしょう（= 1秒間に60回、画面が更新される）。いわゆる、リフレッシュレートというやつです。

画面に、0.2秒だけ刺激を呈示したいとします。これは可能です。1回の更新が1/60秒なのだから、12*1/60 = 0.2であり、12回画面を更新すればよいからです。しかし同様の理屈により、画面に0.22秒だけ正確に刺激を呈示することは困難です。13.2*1/60 = 0.22なので、画面を13.2回更新すればよいわけですが、そんな半端なことはできません。

PsychoPy上で、刺激を0.22秒呈示すると設定することはできます。しかし、それは実際に刺激が0.22秒呈示されていることを意味しません。

というわけで、環境（例：使用するディスプレイ）によって実験手続きに制約が実は存在しているわけです。え...そんなの知らんかった...言うといてや...。安心してください、書いてますよ。

そもそも、液晶ディスプレイって実験で使っていいのだろうか。あれ、そういえば昔受けた実験だと、ブラウン管（CRT）モニターを使用していた気がするぞ。そもそも液晶ディスプレイとCRTモニターの違いはどうなっているのだろう？　安心してください、書いてますよ。

他にも、気を払わなければならない（かもしれない）ポイントはたくさんあります。刺激の色や輝度を厳密にコントロールしたい場合もあるでしょう。特定の情報を探す速さを測定する視覚探索課題の場合、当たり前ですが、物理的に目立つ刺激は見つけられやすくなります。こういった事情を気にする場合、キャリブレーションやガンマ補正などのキーワードが頭に浮かびます。安心してください、書いてますよ。

刺激の呈示ルールはどうすればいいだろう。とりあえずランダムな順序で呈示されるようにしておけばいいかな。でもランダム呈示って本当に万能なのかなあ。安心してください、書いてますよ。

被験者が特定の情報を認識できる／できない境界線（閾値）を知りたい。テクニックがあるのかなあ。そういえば階段法とか聞いたことがあるな...。安心してください、書いてますよ。

と、このように、実は本書は実験心理学の教科書としても非常に有用です。PsychoPyの作者自身が実験心理学者である強みが、まさにここに表れています。巻末には、三角関数の簡単な紹介も載っていて（刺激を特定の位置に呈示するためには、sinやcosの計算が必要だから）、本当に...痒い所に手が届く...。

4. 外部機器との連携

僕自身がこの用途を体験したことがないので、ここはさらっとした紹介になってしまいますが、脳活動関係（fMRI、EEG）の測定をする場合や、眼球運動の計測をする場合における、外部機器との連携が求められる発展的な機能も解説されています。これらの測定やデータの解析と相性の良い商用ソフトもありますが、フリーのPsychoPyでも実現できるのは嬉しいところだと思います。

結び

ここまでレビューしてきたように、本書は、初心者から熟練者まで、幅広いユーザーが待ち望んでいた「実験心理学の教科書」なのではないかと思います。

PsychoPyの作者自身が執筆した書籍であることによる充実度はもちろんのこと、訳者も全員、実験心理学を専門とし精力的に活躍する研究者であるため、翻訳の精度も高いと思います（まだ隅々まで精読したわけではありませんが、ざっと読んだ限りでは、非常に読みやすいと思いました）。

皆様もぜひ！　Enjoy！

2020-05-06

Windows上でzipファイルを解凍するときの、「ファイル名文字化け」への対処

Tips

本記事について

自分用メモです。

僕は論文ファイルを保存するとき、ファイル名を論文の概要にする習慣があります（この記事に感銘を受けました）。

たまった論文をフォルダごと移動させようと思って、ZIPにした後で解凍したところ、ことごとく日本語のファイル名が文字化けしました。どうやらWindowsだとこういうことが起きやすいみたいですね（理由）。

楽な解決方法が見つかったので、記録しておきます。

前提

（最近の）Windows10ではデフォルトで用意されている、Windows Subsystem for Linux（WSL）の機能を用いるなどして、Windowsマシン上でUbuntuを動かせる環境にあることを前提とします。この環境構築方法については、以下の記事が参考になります。

kscscr.com

Ubuntu上でzipファイルの解凍

あとはもうこっちのもんです。

Ubuntu上でunarコマンドを実行すれば、解凍時に、文字コードを自動検出して良きに計らってくれると。

qiita.com

まずはunarをインストールして、次にunarコマンドで解凍するだけ。

sudo apt-get install unar

unar ファイル名.zip

a-zumi.net

補足

WSLとWindows間でファイルをやり取りする場合は、この方法で。

一番楽なのは、解凍対象のzipファイルを

\\wsl$\Ubuntu-(バージョン)\home\(ユーザー名)\

のディレクトリに入れた後で、上記のunarコマンドを実行することかな。

qiita.com

2019-12-22

実験参加者内要因の条件間でパラメータを比較する（パラメータ間の相関を考慮した比較）

Stan

はじめに

本記事は、Stan Advent Calendar 2019　22日目の記事です。

本記事の内容はたのしいベイズモデリング　事例で拓く研究のフロンティア第15章の内容とほぼ同じですが、ちょっとだけ発展的な内容を含んでいます。

次のような手続きで実験を行ったとしましょう。PCのディスプレイ上に、「O」と「D」がランダムに呈示されます。実験参加者は、

「O」が出たと思ったら右のボタンを押す
「D」が出たと思ったら左のボタンを押す

という風に反応することを求められます。刺激が呈示されてから、実験参加者がボタンを押すまでの反応時間が記録されます。一見簡単に見えますが、これをずーーーっと続けていると、時折間違えた反応をしてしまうこともあるでしょう。また、反応が徐々に遅くなっていくと予想されます。すなわち反応の正確性や反応時間は、「情報処理が円滑に行われている程度」を反映する指標と考えられます。

モデリングの方針

このように、実験参加者に二通りの反応が許されているとき（今回は、左右いずれかのボタンを押す）、その結果として得られる反応時間のデータを、drift diffusion modelによってモデリングしてみましょう。そのようなモデリングを行う利点については、Stan Advent Calendar 2019　16日目の記事を参照してください。

本記事の目的

複数の条件を、異なる実験参加者グループが経験している場合（実験参加者「間」計画）と、同じ実験参加者が経験している場合（実験参加者「内」計画）で、モデリングの仕方が変わります。

実験参加者「内」計画の場合、同じ実験参加者がすべての条件を経験するので、条件同士で反応に相関関係が認められると予想されます。例えば反応時間を測定する実験の場合、慎重な人はすべての条件で反応が遅いかもしれないし、せっかちな人はすべての条件で反応が早いかもしれません。

本記事ではこのような実験参加者内計画における反応時間のデータを、条件間の相関を考慮したうえでモデリングする方法を紹介します。

実践

サンプルデータの作成

drift diffusion modelでは、反応時間がWiener分布に従うと考えます。そこで、Wiener分布に従う乱数を発生させましょう。brmsパッケージには、brms::rwiener()という関数が用意されています。これは内部ではRWienerというパッケージを利用しています。もし以下のコードを実行したとき「RWienerというパッケージは無い」というエラーが出た場合には、追加でRWienerパッケージをインストールしておいてください。

# load packages ----------------------------
library(tidyverse)
library(brms)
library(rstan)

# sample data ------------------------------
trial = 1000 # 一つの条件における試行数
N_subject = 5 # 実験参加者数
tau = 0.5 # Wiener分布の非決定時間パラメータ
beta = 0.5 # Wiener分布の開始点パラメータ

set.seed(123)

# 以下は順に、実験参加者1の条件1、条件2、実験参加者2の条件1、...、実験参加者5の条件2
dat_all = rbind(brms::rwiener(n = trial, alpha = 0.6, tau = tau, beta = beta, delta = 2),
                brms::rwiener(n = trial, alpha = 0.75, tau = tau, beta = beta, delta = 2.5),
                brms::rwiener(n = trial, alpha = 0.5, tau = tau, beta = beta, delta = 1.7),
                brms::rwiener(n = trial, alpha = 0.55, tau = tau, beta = beta, delta = 2.4),
                brms::rwiener(n = trial, alpha = 0.65, tau = tau, beta = beta, delta = 2.2),
                brms::rwiener(n = trial, alpha = 0.8, tau = tau, beta = beta, delta = 2.5),
                brms::rwiener(n = trial, alpha = 0.5, tau = tau, beta = beta, delta = 1.8),
                brms::rwiener(n = trial, alpha = 0.6, tau = tau, beta = beta, delta = 2.2),
                brms::rwiener(n = trial, alpha = 0.7, tau = tau, beta = beta, delta = 1.6),
                brms::rwiener(n = trial, alpha = 0.75, tau = tau, beta = beta, delta = 2.1)) %>% 
  dplyr::mutate(ID = rep(x = 1:N_subject, each = 2 * trial),
                condition = rep(1:2, each = trial, time = N_subject))

Wiener分布は4つのパラメータを持ちます。このうち、非決定時間tauと、情報蓄積の開始点betaについては、今回はすべて共通としています。

反応が起こるために必要な情報蓄積の閾値alphaと、情報蓄積の速さdeltaが、実験参加者や条件によって異なるサンプルデータです。このサンプルデータには、以下のように4つの変数があります。

q
- 反応時間（秒）
resp
- 右のボタンを押したか（1）、左のボタンを押したか（0）
ID
- 実験参加者のID。5人いるので、1~5。
condition
- 「O」条件（1）、「D」条件（2）

head(dat_all)

          q resp ID condition
1 0.5569206    1  1         1
2 0.6966415    0  1         1
3 0.5378940    1  1         1
4 0.5503652    1  1         1
5 0.6447536    1  1         1
6 0.7017261    0  1         1

tail(dat_all)
              q resp ID condition
9995  0.6605005    0  5         2
9996  0.6017658    1  5         2
9997  0.6058292    1  5         2
9998  0.5311541    1  5         2
9999  0.7406046    1  5         2
10000 0.5725912    0  5         2

Stanコード

Stanにはwiener_lpdf()という分布関数が用意されているので、基本的にはここに目的変数とパラメータを書き入れるだけです。なお16日目の記事とはやや違う書き方をしています。

パラメータ名は以下のように名付けなおしています。

閾値alpha　→　boundary parameter a
非決定時間tau　→　t_raw及びt_new
- これら二つの違いは後述
情報蓄積の速さdelta　→　drift rate parameter d

data {
  int<lower=0> N;                         //全データ数
  int<lower=0> N_ID;                      //被験者数
  int<lower=0> N_cond;                    //被験者内要因の条件数
  
  int<lower=0> ID[N];                     //被験者を識別するIDが格納されたベクトル
  int<lower=0> cond[N];                   //被験者内条件を識別するラベルが格納されたベクトル
  
  vector<lower=0>[N] RT;                  //反応時間ベクトル
  real<lower=0> minRT;                    //全被験者のなかでの最小RT（推定されたtがこれを超えるとエラーになるため）
  
  int Resp[N];                            //反応のベクトル
  
  real<lower=0, upper = 1> beta;         //starting point (0.5に固定)
}

parameters {
  vector[N_cond] a[N_ID];                 //各条件のboundary parameterを計算する材料。被験者内条件数の長さで、そのベクトルが被験者数ある。
  vector[N_cond] mean_mu_a;               //各条件のboundary parameterが従う、多変量正規分布の平均ベクトル
  corr_matrix[N_cond] omega_mu_a;         //各条件のboundary parameterが従う、多変量正規分布の、分散共分散行列を作るための、相関行列
  vector<lower=0>[N_cond] tau_mu_a;       //各条件のboundary parameterが従う、多変量正規分布の、分散共分散行列を作るための、標準偏差ベクトル
  
  real t_raw;                             //全被験者に共通するnon-decision time parameterを求めるための材料
  
  vector[N_cond] d[N_ID];                 //各条件のdrift rate parameterを計算する材料。被験者内条件数の長さで、そのベクトルが被験者数ある。
  vector[N_cond] mean_mu_d;               //各条件のdrift rate parameterが従う、多変量正規分布の平均ベクトル
  corr_matrix[N_cond] omega_mu_d;         //各条件のdrift rate parameterが従う、多変量正規分布の、分散共分散行列を作るための、相関行列
  vector<lower=0>[N_cond] tau_mu_d;       //各条件のdrift rate parameterが従う、多変量正規分布の、分散共分散行列を作るための、標準偏差ベクトル
}          

transformed parameters{
  real<lower=0, upper=1> t_new;           //全被験者に共通するnon-decision time parameter
  cov_matrix[N_cond] mu_a_covmatrix;      //各条件のboundary parameterが従う、多変量正規分布の、分散共分散行列
  cov_matrix[N_cond] mu_d_covmatrix;      //各条件のdrift rate parameterが従う、多変量正規分布の、分散共分散行列
  
  t_new = inv_logit(t_raw) * minRT;       //ロジスティック変換で0～1に変換し、それに最小RTを掛けているため、
                                          //Wiener分布から推定するt_newは最小RTを超えないように工夫
  
  mu_a_covmatrix = quad_form_diag(omega_mu_a, tau_mu_a); //各条件のboundary parameterが従う、多変量正規分布の、分散共分散行列
  mu_d_covmatrix = quad_form_diag(omega_mu_d, tau_mu_d); //各条件のdrift rate parameterが従う、多変量正規分布の、分散共分散行列
}

model {
  // model ----------------------------------------------------
  for(n in 1:N){
    if(Resp[n] == 1){
      target += wiener_lpdf(RT[n] | exp(a[ID[n], cond[n]]), t_new, beta, exp(d[ID[n], cond[n]])); //exp(a_raw), exp(v_raw)により、両parameterが多変量「対数」正規分布から生成されたことを表現
    }
    else{
      target += wiener_lpdf(RT[n] | exp(a[ID[n], cond[n]]), t_new, beta, -exp(d[ID[n], cond[n]]));
    }
  }
  
  // boundary parameter
  target += multi_normal_lpdf(a | mean_mu_a, mu_a_covmatrix);
  
  // drift rate parameter
  target += multi_normal_lpdf(d | mean_mu_d, mu_d_covmatrix);
  
  
  // priors -----------------------------------------------------
  target += normal_lpdf(mean_mu_a | 0, 5);
  target += lkj_corr_lpdf(omega_mu_a | 2);
  target += student_t_lpdf(tau_mu_a | 3, 0, 3) - student_t_lccdf(0 | 3, 0, 3);
  
  target += normal_lpdf(t_raw | 0, 5);
  
  target += normal_lpdf(mean_mu_d | 0, 5);
  target += lkj_corr_lpdf(omega_mu_d | 2);
  target += student_t_lpdf(tau_mu_d | 3, 0, 3) - student_t_lccdf(0 | 3, 0, 3);
}

ポイントは以下の点です。

wiener_lpdf()のなかで、aとdのパラメータに指数変換を施している
- これらのパラメータは正の値を取ります。ということは、両パラメータにlower=0と下限を指定するだけでもいいような気がしますが、そうすると収束が不安定になりがちでした。
- aやd自体には下限を指定せず、exp()によって正の値に変換しています。
非決定時間パラメータ
- 反応時間の最小値を超えてしまうと推定できなくなるので、それを抑える工夫です。詳しくはコード内のコメントを見てください。
各条件のaとdのパラメータが、多変量正規分布に従うと仮定している
- 今回は2条件しかないので、実質的に2変量正規分布。
- 多変量正規分布は、分散共分散行列をパラメータに持つので、ここで条件間の相関を表現できる
- なお今回は、aとdを最終的に指数変換しているため、結果的にこれらのパラメータは多変量対数正規分布に従うというモデルになっている

1.と2.は推定効率を高めるための工夫ですが、3.が本記事の目的である、相関のあるパラメータを推定するために必要な部分です。

実験参加者「内」条件間でパラメータを比較したいとき、この多変量（対数）正規分布の平均ベクトルを比較することが出来ると思います。今回のモデルではmean_mu_aやmean_mu_dです。

もっとも、そういったやり方ではなく、条件間でパラメータが等しいというモデルと、パラメータが異なることを仮定するモデルを、例えばベイズファクター等を用いてモデル比較するという方法もあると思います。

推定

Rコードは以下の通りです。自己相関が高かったので、thin = 10としています。もうちょっとモデルを工夫できるのかもしれません。

# model -------------------
sm = rstan::stan_model("ddm.stan")

# data --------------------
sdata = list(N = nrow(dat_all),
             N_ID = dat_all %>% dplyr::distinct(ID) %>% nrow(),
             N_cond = dat_all %>% dplyr::distinct(condition) %>% nrow(),
             ID = dat_all$ID,
             cond = dat_all$condition,
             RT = dat_all$q,
             Resp = dat_all$resp,
             minRT = min(dat_all$q),
             z = 0.5)

# sampling --------------------
fit = rstan::sampling(object = sm,
                      data = sdata,
                      seed = 1234,
                      iter = 8000,
                      warmup = 3000,
                      thin = 10)

結果

パラメータが多いので、ここでは結果の要約は載せず、パラメータリカバリできたかどうかに注目します。推定されたaやdを指数変換したものが、Wiener分布のパラメータであることに注意してください。

f:id:das_Kino:20191222215445p:plain

直線上に点が載っていたら、パラメータの真値と推定結果（パラメータの事後期待値）が一致していることを意味します。閾値aについてはうまくリカバリできたようです。

しかし情報蓄積の速さdについては、真値に近い値は得られているものの、ちょっと誤差があります。このパラメータの推定はもうちょっと工夫が要るのかもしれません。

おわりに

今回は2肢選択事態における反応時間データを、drift diffusion modelによってモデリングしました。さらにその際に、実験参加者内要因、つまり同じ実験参加者が全ての実験条件を経験したことにより、条件間でパラメータに相関が仮定される状況を取り上げました。

実験心理学ではこのような実験計画を組むことが多いので、参考になれば幸いです。

Enjoy!