Stan Advent Calendar 2020

24日目

このページは Stan Advent Calendar の24日目の記事です。

執筆者はStan初心者です。先日の日本感情心理学会セミナーでの清水先生の「心理学におけるベイズ統計モデリング：その可能性と研究事例の紹介」をお聞きして、やはりしっかり勉強しないといけないなと思いました（当日資料）。

なお、今回の記事は静岡理工科大学の紀ノ定先生から随時アドバイスをいただきました。ここに感謝の意を表します。

目的

今回は、bayes4psyパッケージにあるflankerデータを使い、Stanで統計モデリングを試みたいと思います。

bayes4psyパッケージのb_reaction_time関数を使うと反応時間をベイズ統計モデリングできるということらしいのですが、これをStanで構築したいと考えました（が、辿り着けませんでした）。無駄に長くなってしまいましたが、誤りやご助言がありましたらいただけますと幸いです。

まずはデータを要約し、可視化します。

データの可視化と仮説検定

library(bayes4psy)

library(dplyr)

library(openvis)


data <- flanker # データを準備する

control_rt <- data %>% filter(result == "correct" &

group == "control") # 正答試行のみを用いる。まずは統制群だけを扱う。

control_rt$subject <- control_rt$subject - 21 # 番号のふり直し


# 試行ごとのデータが準備されているが、それを条件ごとに平均を出す

data_summary <- control_rt %>%

dplyr::group_by(subject, congruency, result) %>%

dplyr::summarise(Mean_RT = mean(rt)) %>%

dplyr::ungroup()

このdata_summaryのデータ構造は以下のようになっています。ここで、各参加者の反応時間は全試行の平均値です。このような一般的な方法では捉えられない、試行ごとのバラツキを考慮に入れた統計モデリングは後半でチャレンジします。

> head(data_summary)

# A tibble: 6 x 4

subject congruency result Mean_RT

<dbl> <fct> <fct> <dbl>

1 1 congruent correct 0.486

2 1 incongruent correct 0.560

3 2 congruent correct 0.404

4 2 incongruent correct 0.500

5 3 congruent correct 0.523

6 3 incongruent correct 0.693

図示します。ここではopenvisパッケージを使ってrain cloud plotをしています（が、このような条件間の変化を線で繋いだうえでrain cloudにするもっと柔軟な方法が欲しいと思っています）。

# rain cloud plotを表示する

library(openvis)

df_1x1 <- data_1x1(

array_1 = subset(data_summary, congruency == "congruent")$Mean_RT,

array_2 = subset(data_summary, congruency == "incongruent")$Mean_RT,

jit_distance = .09,

jit_seed = 321)


g1 <- raincloud_1x1_repmes(df_1x1, align_clouds = TRUE) +

scale_x_continuous(breaks=c(1,2), labels=c("congruent", "incongruent"),

limits=c(0, 3)) +

xlab("flanker condition") +

ylab("RT [ms]") +

theme(plot.title = element_text(hjust = 0.5)) +

scale_y_continuous(breaks=seq(0, 1.5, 0.5), limits=c(0, 1.5)) +

theme_classic(base_size = 16)

g1

上の図の通り、一致条件（congruent）と不一致条件（incongruent）の間の平均に差が見られそうです。

仮説検定（対応のある t 検定）を行うと、条件間に差が認められました。

t = -8.006, df = 23, p-value = 4.234e-08

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

-0.1980531 -0.1167194

sample estimates:

mean of the differences

-0.1573863

Stanで統計モデリング：対応のある t 検定

これからこのデータ発生メカニズムを考えながら、Stanで統計モデリングを試みたいと思います。

このセクションは@shu_ONさんによるBoot Camp 6日目の記事を参考にしました。

まずは、2つの条件の平均値（一致条件：RT[1]、不一致条件：RT[2]）が、多変量正規分布から観察されたと考え、2条件それぞれの反応時間を推定する統計モデリングを行いました。

多変量正規分布の説明は松浦（2016）の97–98ページ、157–160ページに詳しいです。

stanコードを書きます（24_kawashima1.stanとして保存しました）。

data {

int<lower=0> N; // サンプルサイズ

vector<lower = 0>[2] RT[N]; // 反応時間（RT）のベクトルを2条件分用意。多変量正規分布が生成するベクトルの長さ。

}


parameters {

vector[2] mu;

vector<lower=0>[2] sigma;

real<lower = -1, upper = 1> r;

}


transformed parameters {

cov_matrix[2] cov;

cov[1,1] = sigma[1]^2;

cov[2,2] = sigma[2]^2;

cov[1,2] = sigma[1] * sigma[2] * r;

cov[2,1] = sigma[1] * sigma[2] * r;

}


model {

sigma ~ cauchy(0,5); // 半コーシー分布

r ~ uniform(-1, 1);

RT ~ multi_normal(mu, cov);

}


generated quantities{

real mu_diff = mu[2] - mu[1];


// 事後予測分布からの乱数を生成する

vector[2] pred[N];

for (i in 1:N){

pred[i] = multi_normal_rng(mu, cov);

}

}

MCMCサンプリングします。

library(rstan)


# 高速化のためのおまじない

options(mc.cores = parallel::detectCores())

rstan_options(auto_write = TRUE)


data1 = list(

N = length(unique(data_summary$subject)), # 参加者数

RT = cbind(subset(data_summary, congruency == "congruent")$Mean_RT,

subset(data_summary, congruency == "incongruent")$Mean_RT)

)


mcmc_result1 <- stan(

file = "24_kawashima1.stan",

data = data1,

seed = 1,

chains = 4,

iter = 2000,

warmup = 1000,

thin = 1

)

結果を見てみました。Rhatは1.1以下ですし、トレースプロットからも収束していると判断して良さそうです。

print(mcmc_result1,

pars = c("mu[1]", "mu[2]", "mu_diff"),

probs = c(0.025, 0.5, 0.975))


mean se_mean sd 2.5% 50% 97.5% n_eff Rhat

mu[1] 0.62 0 0.03 0.57 0.62 0.67 1475 1

mu[2] 0.78 0 0.04 0.70 0.78 0.86 1507 1

mu_diff 0.16 0 0.02 0.11 0.16 0.20 2715 1


library(bayesplot)

mcmc_sample1 <- rstan::extract(mcmc_result1, permuted = FALSE)

mcmc_combo(mcmc_sample1, pars = c("mu_diff")) # mu_diff だけプロットする

一致条件について事後予測チェックを行います。predは3次元です（4000回生成されたMCMC × 24サンプルサイズ × 2条件）。

mu_pred <- rstan::extract(mcmc_result1)$pred

bayesplot::ppc_hist(y = subset(data_summary, congruency == "congruent")$Mean_RT,

yrep = mu_pred[1:5, ,1])

濃い青のヒストグラムが観測データ（一致条件）の分布を示しています。薄い青色のヒストグラムが事後予測分布です事後予測分布から生成した乱数です。形状はおおよそ似ていると言って良いかと思います。

階層モデルの活用

以上は、参加者ごとに反応時間の平均値を算出し、それを個人の代表値として扱いました。その値が多変量正規分布から観察されたと考え、モデリングを行いました。

ここからは、反応時間は指数ー正規分布に従うという知識（井関、2020、のうこびなど）を使い、階層モデルを利用することを試みます。つまり、代表値ではなく、参加者の全試行のデータを使います。個々の試行の反応時間は、各条件、各実験参加者ごとに独自のパラメータ（μ）をもつ指数ー正規分布に従うと仮定しました。

以下を拡張すれば、個々の試行の反応時間が各条件・参加者ごとに独自のμ、σ、λのパラメータを持つようにできるはずなのですが、今回は途中で力尽きました…。特に、個人ごとのσ、λのパラメータを組み込んだシミュレーションデータを作るところから躓き、階層モデルまで展開できませんでした。

このセクションは「たのしいベイズモデリング 事例で拓く研究のフロンティア」第15章、Boot Camp 8日目の記事を参考にしました。

先に、上のデモデータについて、試行の平均値で代表値をとる前の、ローデータで反応時間のヒストグラムを書きます。反応時間はどのような分布を示すのでしょうか。

library(bayes4psy)

library(dplyr)

library(ggplot2)


data <- flanker

control_rt <- data %>% filter(result == "correct" &

group == "control")

control_rt$subject <- control_rt$subject - 21


library(ggridges)

ggplot() +

theme_bw(base_size = 20) +

geom_histogram(

data = control_rt,

aes(x = rt, fill = congruency)

) +

xlab("RT[sec]") + ylab("count") +

guides(fill = FALSE) +

facet_grid(~congruency)

やはり左右対象ではなく、右に裾の長い分布が確認できました。

まずはじめに、サンプルデータを作り、Stanでパラメータリカバリができるかどうかを確かめます（これをしなければ自分で作ったStanモデリングの正しさを評価できません…）。

library(brms)

library(MASS)

library(ggplot2)


set.seed(1234) # 乱数のタネ

n_trial <- 20 # 試行数

n_subject <- 5 # 参加者数


# 真値

mu <- c(500, 600) # 2条件（congruent, incongruent）の真値

cov <- matrix(c(4.0, 1.0, 1.0, 1.0), ncol = 2)

mu_ind <- mvrnorm(n = n_subject, mu = mu, Sigma = cov) / 1000 #秒単位にスケーリング

sigma <- 0.1 # 単純のため、参加者に共通のパラメータを設定

beta <- 0.1 # 単純のため、参加者に共通のパラメータを設定


# 指数ー正規分布に従う反応時間をサンプリング

cong_rt <- brms::rexgaussian(n_subject*n_trial, mu = mu_ind[,1], sigma = sigma, beta = beta)

incong_rt <- brms::rexgaussian(n_subject*n_trial, mu = mu_ind[,2], sigma = sigma, beta = beta)


# サンプルデータ

RT <- c(cong_rt, incong_rt)

Cond_ID <- rep(1:2, each = n_trial*n_subject) #条件を識別するID

Sub_ID <- rep(1:n_subject, times = n_trial*2) #個人を識別するID


ggplot(data = data.frame(RT, Cond_ID, Sub_ID),

mapping = aes(x = RT, fill = factor(Cond_ID))) +

geom_histogram() +

theme_bw(base_size = 16) +

guides(fill = FALSE) +

facet_grid(~Cond_ID)

サンプルデータが作成できました。このデータは、以下のようになっています。

  1. mu_indは参加者それぞれの条件ごとの平均値パラメータで、2つの変数の平均値を500および600とし、標準偏差をそれぞれ4.0および2.0、2変数間の共分散を1.0とした2変量正規分布に従って生成されました（参考リンク）

  2. 反応時間cong_rtならびにincong_rtは先ほど生成したmu_ind、ならびに参加者間で固定としたsigma、betaのパラメータをもつ指数ー正規分布に従って生成されました

  3. 条件のID、ならびに参加者のIDをそれぞれCond_ID、Sub_IDとして格納しました

それでは、Stanコードを実装します。ここでは階層性が与えられたことに注目します。具体的には、上のシミュレーションデータの生成に沿って、反応時間が個人・条件ごとに特有のμ（mu_ind）をもつ指数ー正規分布に従っていると仮定しました。個人レベルのパラメータμについては、独自の相関を持つ2変量正規分布から抽出されたものだと仮定しました。

stanコードを書きます（24_kawashima2.stanとして保存しました）。

data {

int<lower=0> n_data; //データ数

int<lower=0> n_cond; //条件数

int<lower=0> n_subject; // 実験参加者の数

real<lower=0> RT[n_data]; //目的変数


int<lower=0> Cond_ID[n_data]; //説明変数

int<lower=0> Sub_ID[n_data]; //グループを識別する変数（実験参加者ID）

}


parameters {

//実験参加者ごとのパラメータ

vector<lower=0>[n_cond] mu_ind[n_subject];

corr_matrix[n_cond] Omega; //それぞれの正規分布同士の相関行列

vector<lower=0>[n_cond] tau; //それぞれの正規分布の標準偏差

vector<lower=0>[n_cond] mu;

real<lower=0> sigma;

real<lower=0> lambda;

}


model {

mu_ind ~ multi_normal(mu, quad_form_diag(Omega, tau));

for(i in 1:n_data){

RT[i] ~ exp_mod_normal(mu_ind[Sub_ID[i], Cond_ID[i]],

sigma,

inv(lambda));

// 後述のように、乱数生成時に指定したbetaは、サンプリングされるlambdaと逆数の関係

}


sigma ~ student_t(3, 0, 2);

Omega ~ lkj_corr(2); //相関行列の弱情報事前分布

}


generated quantities{

real mu_diff = mu[2] - mu[1];

}



MCMCサンプリングします。上のコードとは違い、先にコンパイルしてから走らせています。

library(rstan)

options(mc.cores = parallel::detectCores())

rstan_options(auto_write = TRUE)


# シミュレーションデータをリスト化

sample_data = list(

n_data = length(RT),

n_cond = 2,

n_subject = n_subject,

RT = RT,

Cond_ID = Cond_ID,

Sub_ID = Sub_ID

)


sm = rstan::stan_model("24_kawashima2.stan")

mcmc_result_test <- rstan::sampling(

object = sm,

data = sample_data,

seed = 1,

chains = 4,

iter = 2000,

warmup = 1000,

cores = 4

)

結果を一部見てみます。

> mcmc_result_test

Inference for Stan model: 24_kawashima2.

4 chains, each with iter=2000; warmup=1000; thin=1;

post-warmup draws per chain=1000, total post-warmup draws=4000.


mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat

Omega[1,1] 1.00 NaN 0.00 1.00 1.00 1.00 1.00 1.00 NaN NaN

Omega[1,2] 0.02 0.02 0.44 -0.77 -0.33 0.02 0.36 0.80 383 1.02

Omega[2,1] 0.02 0.02 0.44 -0.77 -0.33 0.02 0.36 0.80 383 1.02

Omega[2,2] 1.00 0.00 0.00 1.00 1.00 1.00 1.00 1.00 3941 1.00

tau[1] 0.05 0.00 0.05 0.01 0.02 0.04 0.06 0.16 634 1.01

tau[2] 0.05 0.00 0.05 0.01 0.03 0.04 0.06 0.17 945 1.01

mu[1] 0.36 0.00 0.03 0.30 0.35 0.36 0.38 0.42 1146 1.00

mu[2] 0.47 0.00 0.03 0.41 0.46 0.47 0.49 0.53 837 1.00

sigma 0.06 0.00 0.01 0.04 0.05 0.06 0.06 0.08 470 1.01

lambda 0.14 0.00 0.02 0.11 0.13 0.14 0.15 0.18 698 1.01

関心はmu[1]とmu[2]ですが、上で設定した真値と0.1ずつずれていることに気づかれると思います。これは、rexgaussian関数の処理が以下のようになっているからとのことです［twitter］。

> rexgaussian

function (n, mu, sigma, beta)

{

if (isTRUE(any(sigma <= 0))) {

stop2("sigma must be greater than 0.")

}

if (isTRUE(any(beta <= 0))) {

stop2("beta must be greater than 0.")

}

mu <- mu - beta

rnorm(n, mean = mu, sd = sigma) + rexp(n, rate = 1/beta)

}

<bytecode: 0x113a188a0>

<environment: namespace:brms>

mu <- mu - betaとなっています。なので、今回は一致条件の場合、mu = 0.5 - 0.1で乱数生成されていました。したがって、mu[1]およびmu[2]はおおむねパラメータリカバリできていると判断できそうです。

なお、stanコードでexp_mod_normal(mu, sigma, inv(lambda))としている箇所がありますが、これは乱数生成時に指定したbetaは、サンプリングされるlambdaと逆数の関係になるからとのことです。

（Rhatとn_effがNaNとなっているところはどういうことなのか分かりません…Omegaは相関行列で、対角成分は1だから推定する必要がない、ということでNaNが返るときがあるようです）

ではflankerデータで試したいと思います。

library(bayes4psy)

library(dplyr)


data <- flanker

control_rt <- data %>% filter(result == "correct" &

group == "control")

control_rt$subject <- control_rt$subject - 21


# 条件の置き換え（もっといい方法があると思います…）

control_rt["congruency"] <- lapply(control_rt["congruency"],

gsub, pattern="incongruent",

replacement = 2)

control_rt["congruency"] <- lapply(control_rt["congruency"],

gsub, pattern="congruent",

replacement = 1)

control_rt$congruency <- as.numeric(control_rt$congruency)


library(rstan)

options(mc.cores = parallel::detectCores())

rstan_options(auto_write = TRUE)


data2 = list(

n_data = nrow(control_rt),

n_cond = 2,

n_subject = length(unique(control_rt$subject)), #参加者数

RT = control_rt$rt,

Cond_ID = as.numeric(control_rt$congruency),

Sub_ID = as.numeric(control_rt$subject)

)


sm = rstan::stan_model("24_kawashima2.stan")

mcmc_result_test <- rstan::sampling(

object = sm,

data = data2,

seed = 1,

chains = 4,

iter = 2000,

warmup = 1000,

cores = 4

)

結果を一部お示しします。

mean se_mean sd 2.5% 50% 97.5% n_eff Rhat

mu[1] 0.48 0 0.02 0.44 0.48 0.51 4116 1

mu[2] 0.57 0 0.03 0.52 0.58 0.62 4036 1

mu_diff 0.10 0 0.02 0.07 0.10 0.13 4290 1

条件間の差のmu_diffについて、事後期待値は0.10で、95%ベイズ信用区間は0.07から0.13です。これらは、試行の平均を代表値として扱っていた、はじめのt検定のときと比べ、異なっています。

mu_diff 0.16 0 0.02 0.11 0.16 0.20 2715 1

この差をどう評価するのか私には判断できないのですが、試行間のバラツキという情報を含んだ後者の階層モデルの方が、実際のデータ発生過程・収集過程に関する知識が反映されていることは確かだと思います。

最後に

今回は当初の目的を達成できないまま公開日を迎えてしまいました。最後の階層モデルでは、反応時間の他のパラメータ（σ、λ）の推定を行うことができず、条件間の反応時間の変化の特徴を詳細に記述することができませんでした。オリジナルのデータでは、条件間だけでなく群間の違いを反映させたデータセットだったのですが、ここまでモデリングできませんでした。また、20日目の記事にありましたように、最低限の数学の知識が無ければ、コードを書いて動いてもそれを評価することは困難だと言うことを改めて痛感した次第です。道のりは長く険しいです。