はじめに

DeBruine & Barr (2021) に基づいてpower analysisの基本を説明します
線形混合モデル自体についてもごく簡単に説明します
- 主目的はpower analysisですが、ここで使うデータのシミュレーションは線形混合モデル自体が何をしているのかを理解するのにも役立ちます

Power Analysisとは

ある実験をするときに、どのくらいの量のデータを集めれば充分か？
データ量が少なければ、第二種の過誤を犯す（帰無仮説が偽であるのにそれを棄却しない）可能性（\(\beta\)で表される）が高まる
したがって、検定力＝第二種の過誤を起こさない確率（\(1-\beta\)）を高める必要がある
検定力が低いと、効果量が誇張されたり正負が逆転したりする問題もある（Nicenboim et al., 2018）
検定力を計算するのがPower Analysis

考え方

データ生成用の線形混合モデル（実際の実験では我々が知り得ない「真実」）を用意し、そこからランダムに生成した\(m\)個のデータからモデルを推定する、という手続き（＝机上実験）を\(n\)回繰り返す
\(n\)回のうち\(k\)回、生成元のモデルに存在する効果が、推定されたモデルでも有意となれば、実験でデータを\(m\)個とった場合の検定力は\(k/n\)となる。検定力が低ければ（80%が目安）、\(m\)を増やせばよい

準備

パッケージを読み込み

library(lme4)
library(lmerTest)
library(tidyverse)

データ表示用の設定

kable <- function(df) {
  knitr::kable(df, "html")
}
options(digits=4)

1要因2水準の実験の場合

より複雑なデザインになっても基本的にやることは同じなので、ここでは最も単純なケースを例に説明する ## 題材とする実験
被験者に顔の写真を見てその感情を判定させる
- 説明変数：写真の人物が被験者の身内in-groupか他人out-groupか
- 従属変数：判定までの時間（RT)
つまり、1条件2水準、従属変数は連続データ、within-subject / between-item

モデル

線形混合モデルの下で、被験者\(s\)、アイテム\(i\)の反応時間\(RT_{si}\)は次の式で与えられると仮定される \[RT_{si}=\beta_0+T_{0s}+O_{0i}+(\beta_1+T_{1s})X_i+e_{si}\]
- \(\beta_0\)：切片（intercept）。コーディングによるが、ここでは全試行の平均RTに一致する
- \(T_{0s}\)：被験者ごとのランダム切片。被験者によって反応の速い遅いがあるので、それを調整する。全体的に反応が遅めの被験者なら正の値に、速めの被験者なら負の値になる
- \(O_{0i}\)：アイテムごとのランダム切片。上記と同様。
- \(\beta_1\)：固定効果の傾斜（slope）。下記\(X_i\)と掛け算されることで、説明変数（in-groupかout-groupか）の影響を従属変数（\(RT_{si}\)）にもたらす。写真の人物がin-groupのほうがRTが遅いなら正の値、逆なら負の値をとる（ように\(X_i\)をコーディングする）
- \(T_{1s}\)：被験者ごとのランダム傾斜。写真の人物がin-groupかどうかがRTに影響を与える程度もまた、被験者によって異なると考えられるので、\(\beta_1\)の値を調節する
- \(X_i\)：説明変数。アイテム\(i\)がin-groupなら0.5、out-groupなら-0.5の値をとる（このような値の設定をコーディングという）
- \(e_{si}\)：残差。上記の全ての要因を考慮しても残る個別のデータのブレ
実験結果を分析する際は、観察されたデータから、固定効果がないという帰無仮説（\(\beta_1=0\)）をある有意水準（慣習的に\(\alpha=0.05\)）のもとで棄却できるかを問う（頻度主義）

パラメータの設定

Power analysisでは、データ生成用にモデルを与える
- パラメータは先行研究やパイロット実験の結果をもとにそれらしい値を決めるしかない
- パイロットや過去の実験データがあれば、lmer関数の出力から値を読み取れる（後述）
固定効果は直に設定する

beta_0 <- 800 #全体平均
beta_1 <- 50 #固定効果の傾斜：in-groupのほうが50ms遅い

被験者・アイテムごとのランダム切片はそれぞれ、平均0で、ある標準偏差\(\tau_0,\omega_0\)を持った正規分布をなすと考え、この標準偏差をパラメータとする
- たとえば、被験者1のランダム切片\(T_{0,1}\)は、標準偏差\(\tau_0\)によって与えられる正規分布から採ったひとつのサンプルと考える

tau_0 <- 100 #被験者ごとのランダム切片の標準偏差
omega_0 <- 80 #アイテムごとのランダム切片の標準偏差

被験者ごとのランダム傾斜も、平均0、標準偏差\(\tau_1\)の正規分布をなすと考える
- ただし、ランダム傾斜とランダム切片は全く無関係とは考えにくい。反応が全体的に速め（ランダム切片が小さい）被験者は、刺激の種類によって受ける影響の幅（ランダム傾斜）も小さいと考えられる。そこで両者に相関を想定し、その相関係数を\(\rho\)とする

tau_1 <- 40 #被験者ごとのランダム傾斜の標準偏差
rho <- .2 #被験者ごとのランダム切片とランダム傾斜の相関係数

残差も平均0、標準偏差\(\sigma\)の正規分布をなすと考える

sigma <- 200 #残差の標準偏差

データの生成

被験者とアイテムの数を決める

n_subj <- 100
n_ingroup <- 25
n_outgroup <- 25

アイテムの一覧表を生成する

items <- data.frame(
  #アイテムのIDを1,2,3,...,50とする
  item_id = seq_len(n_ingroup + n_outgroup),
  #1〜25をingroup, 26〜50をoutgroupとする
  category = rep(c("ingroup","outgroup"), c(n_ingroup, n_outgroup)),
  #平均0、標準偏差omega_0の正規分布から生成した50個の値を、アイテムごとのランダム切片とする
  O_0i = rnorm(n = n_ingroup + n_outgroup, mean = 0, sd = omega_0)
)
#固定効果をコーディング
items$X_i <- recode(items$category, "ingroup" = -0.5, "outgroup" = +0.5)

生成したデータ（抜粋）

items[c(1:3,48:50),] %>% kable()

	item_id	category	O_0i	X_i
1	1	ingroup	42.129	-0.5
2	2	ingroup	-30.798	-0.5
3	3	ingroup	5.336	-0.5
48	48	outgroup	36.004	0.5
49	49	outgroup	-60.000	0.5
50	50	outgroup	106.865	0.5

被験者の一覧表も生成する
- 先述の通りランダム切片ランダム傾斜に相関を想定する
- 相関のあるデータの生成をやってくれるfaux::rnorm_multi関数を使う
- 標準インストールされているMASSパッケージを使うやり方も書かれているがわかりにくいのでこちらがおすすめ

subjects <- faux::rnorm_multi(
  n = n_subj,
  mu = 0, #ランダム切片・傾斜の平均はいずれも0
  sd = c(tau_0, tau_1), # ランダム切片・傾斜の標準偏差（これによりサンプルも2個組で生成される）
  r = rho, #相関係数
  varnames = c("T_0s","T_1s") #上記の通り2個組のサンプルが生成されるので、それぞれT_0s（切片）、T_1s（傾斜）と名付ける
)
subjects$subj_id <- seq_len(n_subj) #被験者IDを追加

生成したデータ（抜粋）
- なんとなく相関してる

head(subjects) %>% kable()

T_0s	T_1s	subj_id
-28.60	-35.66	1
37.05	54.66	2
-7.57	22.81	3
131.64	15.08	4
-109.21	-31.36	5
-109.46	-57.66	6

アイテム・被験者のデータをかけ合わせて全試行のデータフレームを作成

trials <- crossing(subjects, items) %>% #被験者・アイテムの表をかけ合わせて
  mutate(e_si = rnorm(nrow(.), mean=0, sd=sigma)) %>% #残差の行を追加（平均0、標準偏差sigmaの正規分布）
  select(subj_id, item_id, category, X_i, everything()) #列の順番を整える

生成したデータ（抜粋）

trials[order(trials$subj_id),][c(1:3,4998:5000),] %>% kable()

subj_id	item_id	category	X_i	T_0s	T_1s	O_0i	e_si
1	1	ingroup	-0.5	-28.6	-35.66	42.129	25.65
1	2	ingroup	-0.5	-28.6	-35.66	-30.798	105.80
1	3	ingroup	-0.5	-28.6	-35.66	5.336	408.16
100	48	outgroup	0.5	-108.8	-22.45	36.004	-165.84
100	49	outgroup	0.5	-108.8	-22.45	-60.000	-145.05
100	50	outgroup	0.5	-108.8	-22.45	106.865	-28.09

これで従属変数（\(RT_{si}\)）は自動的に決まる（冒頭の式を参照）

dat_sim <- trials %>% #先ほどのtrialsデータフレームをもとに
  mutate(RT = beta_0 + T_0s + O_0i + (beta_1 + T_1s) * X_i + e_si) %>% #モデルに基づいてRTを計算し
  select(subj_id, item_id, category, X_i, RT) #必要な行だけを抜き出す

生成したデータ（抜粋）

dat_sim[order(dat_sim$subj_id),][c(1:3,4998:5000),] %>% kable()

subj_id	item_id	category	X_i	RT
1	1	ingroup	-0.5	832.0
1	2	ingroup	-0.5	839.2
1	3	ingroup	-0.5	1177.7
100	48	outgroup	0.5	575.2
100	49	outgroup	0.5	500.0
100	50	outgroup	0.5	783.8

ここまで見てきたデータ生成過程は、線形混合モデルの仮定そのもの
- 従属変数の分布は、線形モデルによって与えられる
- 我々が観察できるのは、そこから採った一定数のサンプル
- 普段の線形混合モデルによる分析は、そのサンプルから、見えないモデルを推定している
- 今回はモデルを初めに与えているが、以下では、そこから採ったサンプルから、元のモデルが適切に推定できるかを問う
ここまでのデータ生成過程が実験1回分なので、繰り返せるように1つの関数としてまとめておく

my_sim_data <- function(
  #パラメータの値は変更できるように、関数の引数（デフォルト値付き）としておく
    n_subj = 100, n_ingroup = 25, n_outgroup = 25,
    beta_0 = 800, beta_1 = 50,
    omega_0 = 80, tau_0 = 100, tau_1 = 40, rho = 0.2,
    sigma = 200
  ) {
  items <- data.frame(
    item_id = seq_len(n_ingroup + n_outgroup),
    category = rep(c("ingroup","outgroup"), c(n_ingroup, n_outgroup)),
    X_i = rep(c(-0.5,+0.5), c(n_ingroup, n_outgroup)),
    O_0i = rnorm(n = n_ingroup + n_outgroup, mean = 0, sd = omega_0)
  )
  subjects <- data.frame(
    subj_id = seq_len(n_subj),
    faux::rnorm_multi(n = n_subj, mu = 0, sd = c(tau_0, tau_1), r = rho, varnames = c("T_0s","T_1s"))
  )
  crossing(subjects, items) %>%
    mutate(e_si = rnorm(nrow(.), mean=0, sd=sigma),
           RT = beta_0 + T_0s + O_0i + (beta_1 + T_1s) * X_i + e_si) %>%
    select(subj_id, item_id, category, X_i, RT)
}

生成したデータの分析

生成したデータをlmer関数で分析する

mod_sim <- lmer(RT ~ 1 + X_i + (1|item_id) + (1+X_i|subj_id), data = dat_sim)

見慣れた出力
- ちゃんと設定したパラメータに近い値が推定されている

summary(mod_sim)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RT ~ 1 + X_i + (1 | item_id) + (1 + X_i | subj_id)
##    Data: dat_sim
## 
## REML criterion at convergence: 67441
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -3.610 -0.656 -0.010  0.660  3.601 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  subj_id  (Intercept)  9688     98.4         
##           X_i          1311     36.2    -0.20
##  item_id  (Intercept)  4235     65.1         
##  Residual             38977    197.4         
## Number of obs: 5000, groups:  subj_id, 100; item_id, 50
## 
## Fixed effects:
##             Estimate Std. Error    df t value Pr(>|t|)    
## (Intercept)    793.1       13.8 124.1   57.64   <2e-16 ***
## X_i             64.8       19.6  51.1    3.31   0.0017 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## X_i -0.027

今回は少し違った整形をする

broom.mixed::tidy(mod_sim) %>% kable()

effect	group	term	estimate	std.error	statistic	df	p.value
fixed	NA	(Intercept)	793.1486	13.76	57.64	124.13	0.0000
fixed	NA	X_i	64.7771	19.57	3.31	51.11	0.0017
ran_pars	subj_id	sd__(Intercept)	98.4269	NA	NA	NA	NA
ran_pars	subj_id	cor__(Intercept).X_i	-0.2025	NA	NA	NA	NA
ran_pars	subj_id	sd__X_i	36.2142	NA	NA	NA	NA
ran_pars	item_id	sd__(Intercept)	65.0752	NA	NA	NA	NA
ran_pars	Residual	sd__Observation	197.4256	NA	NA	NA	NA

1回実験を行い、分析して結果を出力するまでを関数化

single_run <- function(...) {
  #データを生成（引数は...によりsingle_runからmy_sim_dataに渡される）
  dat_sim <- my_sim_data(...)
  #モデルを推定
  mod_sim <- lmer(RT ~ 1 + X_i + (1|item_id) + (1+X_i|subj_id), data = dat_sim)
  #モデルを返す
  broom.mixed::tidy(mod_sim)
}

purrr関数で実験を繰り返す
- boundary (singular) fitという警告が出ることがあるが無視してよい
- failed to convergeは、数回ならよいがあまり何度も出るなら、設定しているパラメータ数に対してデータ量が少なすぎるので増やしてみる
- 安定した結果を得るには1000回必要だそうだが時間かかるので注意

n_runs <- 100 #繰り返す回数
sims <- purrr::map_df(1:n_runs, ~ single_run(n_subj=100))  #被験者数・アイテム数などを変えたい場合はここの引数をいじる

結果を一旦ファイルに保存

write_csv(sims, "sims.csv")

分析結果から固定効果を取り出す

sims <- read.csv("sims.csv") %>%
  filter(effect == "fixed") %>%
  select(term, estimate, std.error, p.value)

分析結果（抜粋）

head(sims) %>% kable()

term	estimate	std.error	p.value
(Intercept)	804.348	15.58	0.0000
X_i	7.279	25.82	0.7792
(Intercept)	792.402	16.35	0.0000
X_i	48.595	26.26	0.0701
(Intercept)	836.912	16.87	0.0000
X_i	36.124	25.60	0.1644

このうち、当該の効果が有意になっていた割合が検定力！

alpha <- 0.05 #有意水準
sims %>%
  group_by(term) %>% #効果ごとに結果をまとめる
  summarize(
    mean_estimate = mean(estimate),
    mean_se = mean(std.error),
    power = mean(p.value < alpha), #p値がalpha未満なら1, そうでないなら0と評価されるので、平均をとれば割合となる
    .groups = "drop"
  ) %>%
  kable()

term	mean_estimate	mean_se	power
(Intercept)	800.93	15.41	1.0
X_i	48.42	23.67	0.5

より複雑な実験デザインの場合

たとえば、読文実験でよくある2x2、within subjectかつwithin itemの場合
基本的な考え方は同じで、モデルが複雑になるのでその分の変数を増やす
- within itemなのでアイテムのランダム傾斜（\(O_{1i}\)）が加わる
- 要因ごとに固定効果（\(\beta\)）、被験者のランダム効果（\(T\)）、アイテムのランダム効果（\(O\)）がある
- 相互作用の固定効果・ランダム効果も独立して設定する必要があることに注意 \[RT_{si}=\beta_0+T_{0s}+O_{0i}+(\beta_1+T_{1s}+O_{1i})X_{1si}+(\beta_2+T_{2s}+O_{2i})X_{2si}+(\beta_{12}+T_{12s}+O_{12i})X_{1si}X_{2si}+e_{si}\]
コードはDeBruine & Barr (2021) のAppendix 2（オンライン）を参照
ラテン方格をシミュレートする小技
- subj_idとitem_idを足して条件数で割った余りによって条件を割り当てればよい
- たとえばこんな感じ

trials <-  crossing(subjects, items) %>%
    mutate(e_si = rnorm(nrow(.), mean=0, sd=sigma))
trials$cond <- (trials$subj_id+trials$item_id)%%4 #Rでは%%で剰余を求める
trials$X_1 <- c(-1/2,-1/2, 1/2, 1/2)[trials$cond+1]
trials$X_2 <- c(-1/2, 1/2,-1/2, 1/2)[trials$cond+1]

ベイズ統計の場合

ベイズ統計には有意水準という概念がないので、検定力という概念もない
ただし、power analysisに準じた方法で、(95%) CIが0を含まず予測した符号の側にある確率を求めることはできる
事前分布のパラメータを、データ生成用のパラメータにも使うことができる
シミュレーションのやり方自体は同じ
詳細は Quant (2021)のAppendix Cを参照

謝辞

きょうの発表は小川さん・陳さんにご紹介いただいた資料を元に作成しました。

参考文献

DeBruine, L. M., & Barr, D. J. (2021). Understanding mixed-effects models through data simulation. Advances in Methods and Practices in Psychological Science, 4(1). https://doi.org/10.1177/2515245920965119
Nicenboim, B., Vasishth, S., Engelmann, F., & Suckow, K. (2018). Exploratory and confirmatory analyses in sentence processing: A case study of number interference in German. Cognitive Science 42(suppl. 4), 1075–1100. https://doi.org/10.1111/cogs.12589
Quandt, J. Power analysis by data simulation. https://julianquandt.com/categories/power-analysis/

Power Analysis

磯野真之介

1/28/2022

はじめに

Power Analysisとは

考え方

準備

1要因2水準の実験の場合

モデル

パラメータの設定

データの生成

生成したデータの分析

より複雑な実験デザインの場合

ベイズ統計の場合

謝辞

参考文献