R で分割プロット分散分析に必要なサンプルサイズを計算する方法

分割プロット分散分析のサンプルサイズ計算の方法。

• 分割プロット分散分析とは？
• 分割プロット分散分析のサンプルサイズ計算のために pwr.anova.test() の一部を抜き出す
• 分割プロット分散分析のサンプルサイズ計算関数
• 分割プロット分散分析のサンプルサイズ計算をする際に大事なこと
• まとめ
• 参考文献

分割プロット分散分析とは？

同じ対象者のある指標を反復して測定したときのデータを分析したいことがある。

対象者をいくつかの群に振り分けて比較したいことが多い。

ここでは、この解析方法を「分割プロット分散分析」と呼んでいる。

分割プロットデザインとも呼ばれる方法だ。

参考となる過去記事はこちら。

https://toukeier.hatenablog.com/entry/2019/01/17/202433toukeier.hatenablog.com

分割プロット分散分析のサンプルサイズ計算のために pwr.anova.test() の一部を抜き出す

まずは、分散分析のサンプルサイズ計算スクリプトを pwr パッケージの pwr.anova.test() から抜き出してみる。

pwr.anova.test.n <- function(k, f, sig.level=0.05, power=0.8)
{
p.body <- quote({
lambda <- k*n*f^2
pf(qf(sig.level, k-1, (n-1)*k, lower=FALSE), k-1, (n-1)*k, lambda, lower=FALSE)
})
n <- uniroot(function(n) eval(p.body)-power, c(2+1e-10, 1e+09))$root
NOTE <- "n is number in each group"
METHOD <- "Balanced one-way analysis of variance power calculation"
structure (list(k=k, n=n, f=f, sig.level=sig.level, power=power, note=NOTE, method=METHOD),
class="power.htest")
}

例えば、4グループの比較で、effect size (f) が中程度で0.25だったとすると、一群45例必要と計算される。

> pwr.anova.test.n(k=4,f=0.25)
Balanced one-way analysis of variance power calculation
k = 4
n = 44.59927
f = 0.25
sig.level = 0.05
power = 0.8
NOTE: n is number in each group

これをもとに、分割プロット分散分析のサンプルサイズ計算のスクリプトを作成してみた。

分割プロット分散分析のサンプルサイズ計算関数

比較する群の数をp、反復して測定する回数をq、群間比較のeffect sizeをf.A.CR、対象者の測定値における自己相関をrhoとすると、以下のような関数になる。

pwr.split.plot.test.n <- function(p, q, f.A.CR, rho, sig.level=0.05, power=0.8)
{
p.body <- quote({
f.A <- f.A.CR/sqrt(1+(q-1)*rho)
lambda <- p*q*n*f.A^2
pf(qf(sig.level, p-1, (n-1)*p, lower=FALSE), p-1, (n-1)*p, lambda, lower=FALSE)
})
n <- uniroot(function(n) eval(p.body)-power, c(2+1e-10, 1e+09))$root
NOTE <- "n is number in each group"
METHOD <- "Split-plot design analysis of variance power calculation"
structure (list(p=p, q=q, n=n, f.A.CR=f.A.CR, rho=rho,
sig.level=sig.level, power=power, note=NOTE, method=METHOD),
class="power.htest")
}

群の数を4、反復回数も4、f.A.CRを0.25、rhoを0.5とすると、一群29例と計算される。

> pwr.split.plot.test.n(p=4, q=4, f.A.CR=0.25, rho=0.5)
Split-plot design analysis of variance power calculation
p = 4
q = 4
n = 28.2526
f.A.CR = 0.25
rho = 0.5
sig.level = 0.05
power = 0.8
NOTE: n is number in each group

rhoが1（完璧な自己相関）の場合、一群45例と計算される。

> pwr.split.plot.test.n(p=4, q=4, f.A.CR=0.25, rho=1)
Split-plot design analysis of variance power calculation
p = 4
q = 4
n = 44.59927
f.A.CR = 0.25
rho = 1
sig.level = 0.05
power = 0.8
NOTE: n is number in each group

この場合は通常の反復測定しない分散分析と同じになる。

完璧に自己相関するなら1回の測定でも同じということだ。

rhoが0（全くの自己相関なし）の場合、一群12例と計算される。

> pwr.split.plot.test.n(p=4, q=4, f.A.CR=0.25, rho=0)
Split-plot design analysis of variance power calculation
p = 4
q = 4
n = 11.92611
f.A.CR = 0.25
rho = 0
sig.level = 0.05
power = 0.8
NOTE: n is number in each group

これは通常あり得ない。

同じ対象者は必ず少しは相関する。

例えば、血圧が高めの人は何度測っても高めのはず。

低めの人はいつも低め。

これは感覚的にわかる。

分割プロット分散分析のサンプルサイズ計算をする際に大事なこと

分割プロット分散分析のサンプルサイズ計算をする際に決めるのが大変なのは、f.A.CRとあるeffect sizeとrhoという対象者内の自己相関だろう。

f.A.CRは、いわゆる群間の標準偏差と、群内の標準偏差の比である。

パイロットスタディや先行研究などがあれば参考にできるが、そうでない場合は、困ってしまう。

行動科学分野では、以下のような目安があるので、それに従うこともできる。

• Small Effect: f=0.1
• Medium Effect: f=0.25
• Large Effect: f=0.4

自己相関 rho のほうはもっと困ってしまう。

これもパイロットスタディや先行研究のデータがある場合は、自己相関を計算してみてもいいかもしれない。

利用可能なデータがない場合は、rho=c(0.3, 0.5, 0.7)などいくつか振って計算してみるという方法もあるかもしれない。

以下、effect sizeを振ってみた計算結果、rhoを振ってみた計算結果を列挙してみる。

# effect size=0.1
> pwr.split.plot.test.n(p=4, q=4, f.A.CR=0.1, rho=0.5)
Split-plot design analysis of variance power calculation
p = 4
q = 4
n = 171.3325
f.A.CR = 0.1
rho = 0.5
sig.level = 0.05
power = 0.8
NOTE: n is number in each group
# effect size=0.25
> pwr.split.plot.test.n(p=4, q=4, f.A.CR=0.25, rho=0.5)
Split-plot design analysis of variance power calculation
p = 4
q = 4
n = 28.2526
f.A.CR = 0.25
rho = 0.5
sig.level = 0.05
power = 0.8
NOTE: n is number in each group
# effect size=0.4
> pwr.split.plot.test.n(p=4, q=4, f.A.CR=0.4, rho=0.5)
Split-plot design analysis of variance power calculation
p = 4
q = 4
n = 11.67164
f.A.CR = 0.4
rho = 0.5
sig.level = 0.05
power = 0.8
NOTE: n is number in each group
# rho=0.3
> pwr.split.plot.test.n(p=4, q=4, f.A.CR=0.25, rho=0.3)
Split-plot design analysis of variance power calculation
p = 4
q = 4
n = 21.71697
f.A.CR = 0.25
rho = 0.3
sig.level = 0.05
power = 0.8
NOTE: n is number in each group
# rho=0.5
> pwr.split.plot.test.n(p=4, q=4, f.A.CR=0.25, rho=0.5)
Split-plot design analysis of variance power calculation
p = 4
q = 4
n = 28.2526
f.A.CR = 0.25
rho = 0.5
sig.level = 0.05
power = 0.8
NOTE: n is number in each group
# rho=0.7
> pwr.split.plot.test.n(p=4, q=4, f.A.CR=0.25, rho=0.7)
Split-plot design analysis of variance power calculation
p = 4
q = 4
n = 34.79044
f.A.CR = 0.25
rho = 0.7
sig.level = 0.05
power = 0.8
NOTE: n is number in each group

まとめ

分散分析プロット分散分析のためのサンプルサイズ計算を紹介した。

pwr パッケージの pwr.anova.test() の一部を抜き出して、修正したものである。

計算に必要はeffect sizeと自己相関を決めねばならず、それが悩ましいポイントだ。

何らかの前提・仮定をおいて、進めていくしかないと思う。

必ずしも先行研究が見つかるわけでもなく、パイロットスタディを行っている時間もないことも多いからである。

参考文献

Bradley, D.R., Russell, R.L. Some cautions regarding statistical power in split-plot designs. Behavior Research Methods, Instruments, & Computers 30, 462–477 (1998).

Some cautions regarding statistical power in split-plot designs | Behavior Research Methods