たとえば重回帰分析とかを行って、効果が「有意ではない」ということを積極的に主張に取り入れたいときは、検出力が十分であったかを確認しなければならない。というのも、わざと出来の悪いモデルを使ったり、サンプルサイズを少なくしたりすれば、有意でないような結果を得ることは簡単だからだ。まあ検出力も、必要な効果量をどれぐらいに設定するかとか、その効果量に含まれる残差の分散の仮定の置き方などによって、ある程度恣意的に引き上げることができてしまうが、やらないよりはいいでしょう。
G*Powerというソフトが便利なのだが（MacでもWindowsでも使える）、シンプルなケースならRのpwrパッケージでできる。
解説はこのあたりを見ればいいかと。
Power Analysis in R
Quick-R: Power Analysis
G*Powerの説明は、ちょっとわかりにくいけど、たとえばこの解説。
https://www.mizumot.com/method/mizumoto-takeuchi.pdf
検出力の分析は「有意水準」「効果量」「サンプルサイズ」「検出力」のどれかを求める手続きで、ネット上や教科書の解説のほとんどは「サンプルサイズの決定」に関するものだが、べつにそれに限らず、これらのうち3つについて仮定を与えれば残りの1つを決められるので、用途は意外と広い。
効果量については、Cohenが提案している、
f^2 > 0.02 ・・・小さめの効果
f^2 > 0.15 ・・・中くらいの効果
f^2 > 0.35 ・・・大き目の効果
という基準を採用する場合は難しく考えなくてよいのだが、たとえば「回帰係数が0.1」というように具体的な値を検出したいと考えている場合は、効果量の表現に平方和が必要になって、データのばらつき具合についても仮定を置かなくてはならなくなる。
ここでは回帰分析で「有意でない」という結果が得られたときに、それを積極的に主張してよいのかどうか判断する方法を考えたいのだが、正直良くわからないので、以下のような方法で問題がある場合はご指摘いただけると助かります。
データはもう手元にあるわけで、サンプルサイズは決まっている。有意水準はたいがい0.05とかだろう。そのうえで、「検出力0.8以上になる効果量を知りたい」あるいは「◯◯という効果量を有意水準0.05で検出できる確率（検出力）を知りたい」ということになるわけである。
で、効果量について上述のCohenのような基準を考えておけば悩まくていいのだが、こんなものでは含意がよくわからないので、現実的な値で語りたくなる。
たとえば回帰分析を行った結果、x1の係数が凄く小さくて有意でなかったとする。しかしこれだけでは、「効果がない」のか「検出力が低い」のか区別がつかない。
で、「変数x1の回帰係数が少なくとも0.1ぐらいあれば、現実社会へのインプリケーションとして意味あるよね」みたいなことが言える時がある。その場合は、「真の係数が0.1だった場合に、今回と同じサンプルサイズ等でその効果が検出できる確率」を計算して、それが0.8とか0.9ぐらいあれば、「十分検出可能なのに、分析の結果有意でなかったのだから、たぶん意味があるレベルの効果は無いのだろう」と積極的に言えるようになる。「本当は0.1ぐらいの効果があるのに、検出力不足で有意にならなかった」わけではなさそうだからだ。
しかしこの計算をするには、効果量を求めるのに、回帰係数だけでなく残差についても仮定を置かないといけない。そこが面倒なのだが、残差（期待値の周りのばらつき）はもとの回帰分析と同じということにしておけば、話がわかりやすい。要するに、予測値からの乖離は同じぐらいになるという前提にして、「意味のある回帰係数」と「元のx」から「yの予測値」を出し、それに「最初の回帰分析をやったときの残差」を足して、「社会的に意味のある効果があった場合の、目的変数y」のデータを仮につくって、そのデータに対して同じように回帰分析（と分散分析）を行って、効果量を求めるわけである。期待値のまわりのばらつきが同じでも、「意味のある回帰係数」が大きければ、yの変動に占める残差の割合は小さくなるので、検出力が上がることになる。
数学が得意な人が考えればわざわざ仮のデータを作ったりしなくても求められそうだけど、自分は「説明変数同士に相関があるかもしれないからタイプIII平方和を使う」というあたりでもう考える気がなくなったので↓のようにやってみた。
以下のような関数を作っておいて、
library(car) library(pwr) calculate_power <- function(model, val.number, sufficient.effect, sig.level=0.05) { # library{car}{pwr}が必要。 # modelにはlmの推定結果を与える。推定するときにlm(..., x=TRUE, y=TRUE, model=TRUE)としておく。 # val.numberには、切片から数えて何個目の変数について計算するか番号を入れる。 # sufficient.effectには、検出したい大きさの回帰係数を入れる。 # val.number, sufficient.effectはベクトルで複数指定してもよい。 coef.sufficient <- model$coefficients coef.sufficient[val.number] <- sufficient.effect # 回帰係数が指定した「社会的に意味ある値」だった場合のy y.sufficient <- t(coef.sufficient) %*% t(model$x) + model$residuals # xに切片も入ってる # データセットの形にする df.sufficient <- data.frame(y.sufficient=t(y.sufficient), as.data.frame(model$model[,-1])) # modelの1列目はyになってる # 再び回帰して分散分析結果を得る lm.sufficient <- lm(y.sufficient ~ ., data=df.sufficient) anova.sufficient <- Anova(lm.sufficient, type='III') ss.factor.sufficient <- anova.sufficient$`Sum Sq`[val.number] ss.resid.sufficient <- anova.sufficient$`Sum Sq`[length(anova.sufficient$`Sum Sq`)] df.factor.sufficient <- anova.sufficient$Df[val.number] df.resid.sufficient <- anova.sufficient$Df[length(anova.sufficient$Df)] # 効果量を計算 eta.sq <- ss.factor.sufficient / (ss.factor.sufficient + ss.resid.sufficient) f.sq <- eta.sq/(1-eta.sq) # この効果量を検出できる確率 pow <- pwr.f2.test(u = df.factor.sufficient, v = df.resid.sufficient, f2 = f.sq, sig.level = sig.level) output <- list(pow$power, f.sq) names(output) <- c('Power', 'f2') return(output) } }
まず、目的変数y、説明変数x1〜x3を作って、回帰分析をやってみる。
> set.seed(1) > x1 <- rnorm(100,0,1) > x2 <- rnorm(100,5,1) > x3 <- rnorm(100,10,1) > y <- 2*x1 + 3*x2 + 3*x3 +rnorm(100, 0, 10) > > model.1 <- lm(y~x1+x2+x3, x=T, y=T, model=T) > summary(model.1) Call: lm(formula = y ~ x1 + x2 + x3, model = T, x = T, y = T) Residuals: Min 1Q Median 3Q Max -25.8201 -5.2738 0.0281 6.5058 18.2364 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -7.1879 11.3108 -0.635 0.5266 x1 1.4204 1.1169 1.272 0.2065 x2 2.4506 1.0485 2.337 0.0215 * x3 4.0462 0.9712 4.166 0.0000677 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 9.98 on 96 degrees of freedom Multiple R-squared: 0.1981, Adjusted R-squared: 0.1731 F-statistic: 7.908 on 3 and 96 DF, p-value: 0.00009087
x1の真の係数は2だけど、1.4204と推定され、pは0.2065だから有意ではない。これが、サンプルサイズ（もしくはばらつきが大きい）のせいなのか、本当に効果がないのかを考えたいなら、「x1がどの程度大きければ意味のある効果なのか」を決める必要がある。
仮にそれが1ぐらいだとすると……
> calculate_power(model=model.1, val.number=2, sufficient.effect=1)$Power [1] 0.1457216
今回のようなサンプルサイズやばらつき具合のデータでは、「x1の係数が1」という効果が仮にあったとしても0.146ぐらいの確率でしか検出できないみたいだから、「x1に効果はない」と言うのは憚られる。
特に、あってほしい効果が検出されないなら単に「残念です、サンプルサイズを増やすか、実験計画を工夫して誤差を減らしましょう」で済むが、自分の説を主張する上で「あってほしくない効果」だった場合は、ついうっかり「どうだ、俺の言ったとおり効果はないじゃないか！」と言いたくなってしまい危険である。
ちなみに、真の値である「係数は2」は検出できるのかというと……
> calculate_power(model=model.1, val.number=2, sufficient.effect=2)$Power [1] 0.4330289
4割ぐらいしか検出できないようだ。
