simrを使って一般化線形混合モデル（GLMM）の検定力分析をする

参考

Package ‘simr’

一般化線形混合モデル（GLMM）の検定力分析がしたい

　個人の主観だが，さすがに割合データは変換+分散分析ではなく，GLMMで分布に二項分布を指定して分析したいと思っている。また勝手な偏見だが，検定力分析をしないと怒られそうと思っている。そこでGLMMの検定力分析を行いたいのだが，人気のG*powerはGLMMに対応していないようなので，Rパッケージのrsimを用いて（『事後』の）GLMMの検定力分析を行った。

分析対象

　文章を読んだ後に表示される質問文の内容が正しいかどうか答える再認課題を行ってもらった時の正答率が上図である。再認課題は文章1つにつき3問で，参加者は2つの文章について全6問への回答が求められた。この時の再認課題の正答率について，参加者間条件間で有意な差があるかどうかとその検定力を見た。サンプルサイズは各条件24人ずつであった。

※条件1の人がほぼほぼ全員満点をとっているので，本当は単にGLMMするだけだと適切ではないと思う。

とりあえず分析

library(lme4)
library(lmerTest)
fit  <-  glmer(Resp ~ A + (1|Contents) , data=data, family=binomial)

　以上のモデルを用いて分析した。再認課題の各設問での回答を従属変数，参加者間条件の主効果を固定効果として，再認課題の対象となった文章をランダム切片に投入し，分布は二項分布（対数リンク）を指定した。

> summary(fit)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: binomial  ( logit )
Formula: Resp ~ A + (1 | Contents)
  Data: data

    AIC      BIC   logLik deviance df.resid 
  203.8    214.8    -98.9    197.8      285 

Scaled residuals: 
   Min      1Q  Median      3Q     Max 
-6.7885  0.2067  0.2711  0.4168  0.7885 

Random effects:
Groups   Name        Variance Std.Dev.
Contents (Intercept) 0.7501   0.8661  
Number of obs: 288, groups:  Contents, 4

Fixed effects:
           Estimate Std. Error z value Pr(>|z|)    
(Intercept)   2.2977     0.4934   4.656 3.22e-06 ***
A            -1.2202     0.4025  -3.031  0.00243 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
 (Intr)
A -0.177

　条件間の差は有意で，条件1の人の正答率が条件2の人の正答率よりも高いという結果が得られた。

検定力分析

library(simr)
A_Power <- powerSim(fit, test=fixed("A", "z"), nsim=1000, alpha=0.05)

　simrパッケージよりpowerSim関数を用いた。引数testでは，検定力分析の対象を指定している。単一の固定効果についての検定力を見たかったので，対応するものを指定した。powerSim(モデル名, ("検定力分析をしたい固定効果", "適用する検定方法"), シミュレーション回数, 有意水準)となっている。

　検定方法は，用いたフィッティングの方法によって異なるものを指定する。glm，glmer，lmerによるフィッティングを行った場合にはz検定（z）を指定する。その他，回帰分析（lm + lmerTestパッケージをインストールしていれば先に挙げた混合モデルでも可とのこと）を行った場合にはt検定（t），分散分析（anova）を行った場合には尤度比検定（lr）を指定する。これらの指定を間違えると検定力が0% or 100%でしか出てこなくなったりしたので，極端な値が出た場合には一度確認した方がいいかもしれない。

　デフォルトだと，シミュレーションの回数は1000回，有意水準は0.05に設定されている。なお，t値の絶対値が2以上であることを有意水準の代わりに用いる場合，simrパッケージ製作者（Green et al., 2016）はt=2に対応する閾値としてα=0.045を設定すれば良いとしている。

　最後に，引数testで指定出来る検定力分析の対象は，単一の固定効果（fixed）の他にも，モデル式を用いた比較（compare，fcompare，rcompare），単一のランダム効果（random）がある。

検定力分析の結果

> A_power
Power for predictor 'A', (95% confidence interval):
     85.90% (83.59, 88.00)

Test: z-test
     Effect size for A is -1.2

Based on 1000 simulations, (3 warnings, 0 errors)
alpha = 0.05, nrow = 288

Time elapsed: 0 h 3 m 2 s

nb: result might be an observed power calculation

　要因Aの検定力は85.90%と十分なものであった。なお，ランダム切片に参加者のみを指定したモデルでは66.50% [83.59, 88.00]，ランダム切片に再認課題の対象となった文章と参加者を指定したモデルでは70.70% [67.77, 73.51]であった。

　今回の分析では，警告（boundary (singular) fit: see ?isSingular）を出された。isSingularを見ろと言われたのでヘルプを覗いたところ，モデルが影響のない効果を含んでいる≒過剰適合していると警告されたっぽい。isSingular関数はモデルがそのような特異モデルかTRUE or FALSEで返してくれる。実行したところ，以下の通り特異モデルではないと帰って来たし，そもそも省ける効果も無いので警告は無視した。

> isSingular(fit, tol = 1e-4)
[1] FALSE

事前に検定力分析したい

　事故を避けるために，普通は検定力分析を事前に行うと思われるし，自分もそうしたい。その際はpowerCurve関数を用いて分析を行う。powerCurve関数は，効果量，有意水準とデータセットを指定すると，様々なサンプルサイズに対してpowerSim関数を実行し，そのサンプルサイズでの検定力を返してくれるので，検定力分析が80％以上になる現実的なサンプルサイズで実験すればいいと思う。（参考サイト）