一般(化)線形混合モデルの例数設計<固定効果編>

AnsNumSim.jpeg SubNumSim.jpeg AnsSubNumSim.jpeg
Figure:  一般化線形混合モデルの検定力シミュレーション
左からグループ内人数、グループ数、グループ内人数とグループ数を定数倍した場合の検定力の推移を示す。
author: Unadon (見習い飯炊き兵) 動作環境:Mac OS Sierra 10.12.1; R version3.3.1

 

本稿の目的

  1. 一般化線形混合モデルの例数設計ができるようになる
  2. なお、効果量には固定効果の傾きを使用する(ランダム効果を用いた例数設計は今回は扱わない)

主要参考文献および出典はこちらです

・ 本稿のメイン(コードの出典)

http://onlinelibrary.wiley.com/doi/10.1111/2041-210X.12504/full

・ GLMM検定力のシミュレーションをした論文(この論文から、色々な効果量を使った場合のシミュレーション研究論文に出会えます)

http://journals.sagepub.com/doi/pdf/10.1177/1094428107308906

・ パイロットデータを取らずに、シミュレーションデータを全て一から作る方法(本稿では使わない)

http://onlinelibrary.wiley.com/store/10.1111/2041-210X.12504/asset/supinfo/mee312504-sup-0002-AppendixS2.html?v=1&s=111a69aebb9f8086bd5684d0da84cc1959b52e0f

はじめに

最近、心理学領域で経験サンプリング法という調査方法が定着してきています。スマートフォンを介したアンケート調査を行い、個人に対して複数回の回答を求めるというものです。 反復測定データが得られます。

変数間の関係を(非)線形モデルで分析していくケースがほとんどで、 個人レベルと集団レベルの両方のトレンドを考慮したマルチレベルな分析(一般(化)線形混合モデル:GLMM)が頻繁に使用されています。

その時、問題になるのが検定力、サンプルサイズの話。 検定力(power)とは、帰無仮説を正しく棄却できる確率のことです。power = 1-β,と書かれることもあります。 必要な水準として0.8を設定するのが一般的です。

検定力は効果量とサンプルサイズ、危険率(有意水準,一般に0.05)から計算することが可能であり、 逆に、必要な検定力、危険率と効果量が定まれば、それらを使って必要サンプルサイズを求めることができます。

ですが、相関分析やt検定の場合とは違って、GLMMでこれをやろうとすると、結構難しいんですよね。根本的にやり方が違います。

今回は、そんな経験サンプリング(GLMMの、といったほうがいいかもしれません)の例数設計をやってみたいと思います。

例数設計とは、つまるところデータを何人分とればいいのか計算する手続きのことです。 必要な効果量、必要な検定力、必要な危険率のもとで、必要最小となるサンプルサイズを求める手続きというと精確でしょうか。 事前の検定力分析(apriori power alnalysis)と呼ばれたりします。

なお、GLMMそのものの解説、例数設計そのものの基本的な考え方については、ブログ程度では書ききれませんので、別途お調べ願いたいと思います。

一般(化)線形混合モデルの例数設計(例)

GLMMの例数設計にはいくつかのやり方があります。

1. 方程式で計算する

古くには提案されていたりしますが、現時点ではシミュレーションを使った方法が主流だと思います。今回は扱いません。

2. シミューレーションで求める

シミュレーションによる方法のなかにも、2パターンのやり方があります。

2.1 全データシミュレーション

ハードルが高いので今回は扱いません。線形代数学の知識と複雑な乱数生成法が必要です。

2.2 パイロットデータを用いたシミュレーション

こちらは実際の研究においてもメリットが大きいように思います。そして、方法もある程度整っています。 「パイロットデータを取って例数設計をする」ということは、必然的に統計モデルを確定させ、実際に分析まで行ってみる手続きが含まれます。 これで致命的なミスを防ぐことができます。”このデータの取り方間違ってた!”とか、”データってこういうふうな表にまとめるのか!”ということにも気がつくことができます。 また、仮定の上の仮定の上の…全部仮定のシミュレーションデータで例数設計するよりは、本解析時と同じモデルで実際の小サンプルデータで例数設計ができるので、間違った検定力分析にはなりません。

以上より、今回は「パイロットデータを取得した上で、シミュレーションによる例数設計を行う」ことを実践していきます。

なお、使用する効果量には、固定効果(e.g. グループレベル)の傾きを今回は使用します。およそ、多くの方の関心は固定効果の傾きだと思いますので。

一般化線系混合モデルの例数設計の手続き

1. 研究デザインと目的を定め、主たる統計モデルを事前に決めておく

ここが一番重要だと思います。

今回は従属変数zを標準化せず、標準化しない独立変数xで予測する。ランダム切片のみのモデルを使います。 参照元に合わせました。

が、実際の研究では、

「経験サンプリング法を使用してデータを取得。従属変数は怒り感情、独立変数は腰痛強度。従属変数はbetweenで標準化して、独立変数はwithinで標準化。グループごとに平均値を求めてbetweenで標準化した変数を統制変数に用いる。一時点前の腰痛を統制するために、within標準化した腰痛得点を1時点ずらした変数を作る。欠損値は補完しない。その他、時間と日時、曜日、性別と年齢を統制変数に加え、ランダム切片とランダム傾きを用いたモデルをメインに据える。モデル比較のため、切片だけを仮定したモデル、ランダム切片のみのモデル...etcを用意し、最終的にAICを用いてモデル比較する」 」

くらいまで決めておきましょう。

2. 必要な効果量(何の変数の傾きか)とその値を決める

一番関心のある(標準化した)独立変数の、グループレベルの(固定効果)傾きを”中程度の効果量”0.5”と定めるのが一般的でしょうか。

今回は参照元にあわせて非標準化効果量-0.05を期待する効果量にします。

3. パイロットデータの取得

実際の研究と同じ構造で少サンプルのデータを取得します。基準は特にありません。

今回はデフォルトで用意されたサンプルデータを使用します。

3. 回答回数と参加者数を変動させた際の必要サンプルサイズをシミュレーションする

Rを使ってシミュレーションしていきます。結構時間がかかります。

手順としては、

3.1 一回GLMMで分析する

3.2 分析結果の固定効果の傾き(効果量)を設定したい値に変えてしまう

3.3 人数、グループ数、それらの定数倍、を変化させたときの検定力の推移をシミュレーションする

という流れになります。

例数設計の実施

では、実際にやっていきましょう。

ライブラリの読み込みとデータの確認

今回使用するデータを確認します。使用するのは、従属変数にz、独立変数にx、グループの識別変数gだけです。

経験サンプリングをした、という場合は、gをグループではなく”個人のID番号”と読み替えましょう。

############################################################
#  Generalized Linear Mixed Model: A priori Power simulation
############################################################
invisible({rm(list=ls());gc();gc()})

#Install
  #install.packages("simr")

#loading library
library(simr) #power simulation for linear mixed model
library(lme4) #linear model, (Genaralized) Linear Mixed Model Analysis
library(sjPlot) #plot
library(arm) 
library(ggplot2) 

print(simdata)

SimData.jpeg

パイロットデータで分析

続いて、一回GLMMで分析を行います。従属変数にポアソン分布を仮定した一般化線形混合モデルによる分析です。

独立変数はxのみ。ランダム切片のみの仮定をおく、シンプルなGLMMです。 パッケージは必ずlme4を使用してください。

#pilot Study fit----------------------------------------
model1 <- glmer(z ~ x + (1|g), family="poisson", data=simdata)

#SimpleSummary
summary(model1)

次のような結果が出てくると思います。xの傾きは有意らしい。

Generalized linear mixed model fit by maximum
  likelihood (Laplace Approximation) [glmerMod]
 Family: poisson  ( log )
Formula: z ~ x + (1 | g)
   Data: simdata

     AIC      BIC   logLik deviance df.resid 
   109.0    113.2    -51.5    103.0       27 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.28918 -0.41836 -0.03916  0.57284  1.29631 

Random effects:
 Groups Name        Variance Std.Dev.
 g      (Intercept) 0.08345  0.2889  
Number of obs: 30, groups:  g, 3

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.54079    0.27173   5.670 1.43e-08 ***
x           -0.11481    0.03955  -2.903   0.0037 ** 
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
  (Intr)
x -0.666

パイロットデータの分析結果を可視化

加工しにくいため使いにくい部分はありますが、GLMMの結果を一瞬で可視化できるパッケージsjplotを使用します。

便利ですので使ってみてください。lme4による分析結果しか扱えません。


#plot and table-------------------------------------
#以下 全ggplotに適用(theme_set)するもの
#"グラフタイトル14pt",”日本語表示”
ggplot()+theme_set(theme_bw(base_size = 14,base_family="HiraKakuProN-W3"))

#ランダム切片の値と95%信頼区間 (と ランダム傾き) : 今回はランダム傾き使ってないのででない
sjp.glmer(model1)

#固定効果の傾きと信頼区間
sjp.glmer(model1,type="fe")

#固定効果の線をプロット
sjp.glmer(model1,type="fe.slope",y.offset = .4)

#固定効果の傾きで、ランダム切片。グループごと。
sjp.glmer(model1,type="fe.ri")

#結果を表で出力
sjt.glmer(model1)

sjplot1.jpeg sjplot2.jpeg sjplot3.jpeg sjplot4.jpeg sjt.jpeg

シミュレーションの実行

1時間近くかかるので、お覚悟を。

まずは、変数xの固定効果の傾き(fixef()[“x”])を、-0.05に変更します。

標準化されたデータを使用した場合は、cohenの基準に従い中程度の効果量0.5を使用する場合が一般的です(冒頭の出典論文を参照) ここでは、コードの参照元にあわせて、-0.05の傾きを期待する効果量とします。

#power simulation using pilot data--------------------------------------
#expected fixed effect 
fixef(model1)["x"] <- -0.05

続いて、今のサンプルサイズだと検定力がどれくらいあるのか(事後の検定力分析です)を確認しましょう。

#posteriori power(in the case of expected fe)
powerSim(model1)

Power for predictor 'x', (95% confidence interval):
      32.80% (29.89, 35.81)

Test: z-test
      Effect size for x is -0.050

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

Time elapsed: 0 h 2 m 28 s

32.8%しかありません。これは困った。

サンプルサイズを変化させて検定力シミュレーション

続いて、サンプルサイズ(ESMなら回答回数)を増やしたときの検定力の推移をシミュレーションしていきます。 n=10から、n=20まで増やしてみましょう。グループは3グループあるので、合計60人まで参加することを想定したものになります。 なお、経験サンプリングならば、「3人が20回ずつ回答した場合」 まで、回答回数を拡張してみます。

#sample size(with-in group) extend
  #in other words, "extend the number of answers in each subjects participated ESM"
  model2 <- extend(model1, along="x", n=20)
  powerSim(model2)

#power curve
pc2<-powerCurve(model2)
print(pc2)
plot(pc2)

結果をprintします。

Power for predictor 'x', (95% confidence interval),
by largest value of x:
      3:  5.50% ( 4.17,  7.10) - 9 rows
      5:  9.30% ( 7.57, 11.27) - 15 rows
      7: 15.30% (13.12, 17.68) - 21 rows
      9: 24.00% (21.38, 26.77) - 27 rows
     11: 39.80% (36.75, 42.91) - 33 rows
     12: 46.80% (43.67, 49.95) - 36 rows
     14: 66.90% (63.89, 69.81) - 42 rows
     16: 80.50% (77.91, 82.91) - 48 rows
     18: 89.40% (87.32, 91.24) - 54 rows
     20: 95.40% (93.91, 96.61) - 60 rows

Time elapsed: 0 h 16 m 16 s

n=18とれば、十分そうですね。

AnsNumSim.jpeg

グループ数を変化させて検定力シミュレーション

次は、グループ数を変化させたときのシミュレーションです。 n=10(サンプルサイズ、ESMなら回答回数10回)で固定し、グループ数(ESMなら参加者数)を3から15まで増やしていきます。

経験サンプリングをするなら、現実的にはこのシミュレーションを使うことになるでしょうね。

回答回数を増やすと欠損値が増えますし、ドロップアウトも出てきます。

回答回数は研究上必要最低限で、参加者数を増やすのが、一番セーフティーな気がします。

(一人あたりに高額な参加報酬を渡すなら、話はべつかもしれませんが)

#sample size(number of group or subjects) extend
#in other words, "extend the number of subjecs  participated ESM"
model3 <- extend(model1, along="g", n=15)
pc3 <- powerCurve(model3, along="g")
plot(pc3)

print(pc3)による結果は省略します。図を見ると、11グループ取ればOKというかんじでしょうか。

SubNumSim.jpeg

グループ数と参加者数をパイロットデータ✕a倍したときのシミュレーション

パイロットデータが各グループ10人で3グループですので、定数倍すると、

2倍: 20人ずつ6グループ(トータル120行)

3倍: 30人ずつ9グループ(トータル270行)…

と言った具合に増えていきます。

#multiplication each n of anser and subjects in ESM constantly
model4 <- extend(model1, within="x+g", n=5)
pc4 <- powerCurve(model4, within="x+g", breaks=1:5)
print(pc4)
plot(pc4)

結果を見ると、4倍(40人ずつ12グループ)のところで検定力が80%を超えていますね。

AnsSubNumSim.jpeg

結果の保存

.Rデータとして、結果の全てを保存しておきましょう。保存して、再度読み込めば、シミュレーションの結果が全て復元されます。

#----------------------------------------------
#saving Rdata(シミュレーション1回につき時間がかなりかかるので、Rデータとして結果を全て保存しておく)
#Rデータの保存
save.image("AprioriPowerAnalysisForGLMM.Rdata",compress="xz")

#Rdataの読み込み
load("AprioriPowerAnalysisForGLMM.Rdata")
#----------------------------------------------

以上、GLMMの固定効果を使った事前の検定力分析のシミュレーションでした。

パイロットデータを取って、一回分析し、ご自身のモデルで例数設計を行い、いざ本調査にのぞまれてください。

Enjoy!!

Written on January 19, 2017