3群間の検出力検定（パワーアナライシス）のための方法

論文投稿すると、”検出力検定はしたのか？”　と聞かれることが増えてきました

普段はEZRでほとんどの検定を行っていますが、最近は統計の複雑化によりEZRの限界を感じることもあります。

今回はその一つで、”3群間での検出力比較”　です。

Rスクリプトを記載して、Rコンソールに入れてもよいのですが、RのUIは使いにくく、こういったブログから見つけたものを貼り付けるとエラーが頻発します。

Rstudioにすればよいのですが、なんとなくなじみがないので敬遠しています。

最近、ふとしたことでGoogle Colaboを使うことがあり、その使いやすさに驚愕。

”Google Colabo”と検索して、Google アカウントでログインするだけで使用可能

ここのファイルから”新しいノートブック”を開く

コーディングを開始する。コピペすればOK

Kruskal Wallis検定の検出力検定のコード

41例、28例、44例という不均等な3群間で、データが非正規分布のデータでの検討です

from scipy.stats import kruskal, norm, logistic
import numpy as np  # NumPyをインポート

def simulate_kw_power(n_sim=1000, group_sizes=[41,28,44], effect_shifts=[0,0.5,1.0], dist="normal"): 
   """ 
   Kruskal-Wallis検定の検出力をシミュレーションで計算する関数

   Parameters:
    - n_sim: シミュレーション回数 (default:1000)   
　  - group_sizes: 各群のサンプルサイズ (例: [41,28,44])
    - effect_shifts: 各群の位置パラメータシフト (基準群からの差)
    - dist: データ分布 ("normal" or "logistic")

    Returns:    - 検出力 (0~1)   
    """
   np.random.seed(42)  # 再現性確保    
   rejections = 0        

   for _ in range(n_sim):
       # データ生成
       if dist == "normal":
            groups = [ 
               norm.rvs(loc=effect_shifts[i], scale=1, size=size)
                for i, size in enumerate(group_sizes) 
           ] 
       elif dist == "logistic":
            groups = [ 
               logistic.rvs(loc=effect_shifts[i], scale=1, size=size) 
               for i, size in enumerate(group_sizes) 
           ] 
       else:  # distが"normal"または"logistic"以外の場合のエラーハンドリングを追加
           raise ValueError("dist must be either 'normal' or 'logistic'")

       # Kruskal-Wallis検定 
       try:
           stat, p_value = kruskal(*groups) 
           if p_value < 0.05:
               rejections += 1
       except ValueError: # エラータイプをValueErrorに変更
           # ランク付けに失敗する稀なケースをスキップ 
           continue

   return rejections / n_sim

# シミュレーション実行例
power = simulate_kw_power(
    n_sim=1000,
    group_sizes=[41, 28, 44],
    effect_shifts=[0, 0.5, 1.0],
    dist="normal"
)

print(f"Estimated Power: {power:.3f}")

上記のコードを貼り付けて、15秒ほど待つと

と結果が得られます。

もしよければ使って感想を教えてください
ちなみに、ANOVAは明日以降に出します。