割合の差のスコア検定
以下のツイートをしましたところ、情報提供をいただきましたので、まとめた文章を作成しました。
リスク差の検定の方の式を見るとWald型の統計量で、カイ二乗検定（これはスコア型の検定に一致する）と違うものになっているので、同じにはならないと思いました。https://t.co/Aj3Gd3EoVg
— Triad sou. (@triadsou) 2024年2月28日
非劣性仮説の割合の差のWald検定でマージンを0にした時の結果と、分割表のカイ二乗検定（スコア検定に一致）の結果がずれるという状況だったので、それならば非劣性仮説の割合の差のスコア検定でマージンを0にした時の結果を出せば、分割表のカイ二乗検定と一致すると思い、Farrington–Manning検定の結果が同じになることを確認したという内容です。
Farrington–Manning検定 (Farrington & Manning, 1990) は、非劣性仮説の割合の差のスコア検定と（割合の比の方もですが）として知られていて、SASやRのパッケージではこのスコア検定はFarrington–Manning検定と称されています。
しかし、実際にはこれと同値な割合の差のスコア検定統計量がMee (1984) で最初に与えられており、佐藤（1988）ではMeeの方法がスコア検定の式に一致することが示されています。
Farrington & Manning（1990）は非劣性仮説の検定とサンプルサイズ設計を主題にしたもので、なんらかの理由でこっちの方が有名になり（Mee（1984）はレターだったことも関連しているでしょうか？）、Farrington–Manning検定という呼び方が定着してしまったのかと思います。
Meeの方法と佐藤俊哉先生の方法の詳細については、佐藤（1988）で導出も含めて非常に丁寧に書かれていますので、そちらをご参照ください。
割合の差のスコア検定統計量
佐藤（1988）の(3)式
\[
\frac{x - n(\delta + \hat{p}_2)}{(\delta + \hat{p}_2)(1 - \delta - \hat{p}_2)} + \frac{y - m \hat{p}_2}{\hat{p}_2(1 - \hat{p}_2)} = 0
\] の下の
\[
T(\delta)=\left\{ \frac{x - n \hat{p}_1}{\hat{p}_1(1 - \hat{p}_1)} \right\}^2
\left\{ \frac{\hat{p}_1(1 - \hat{p}_1)}{n} + \frac{\hat{p}_2(1 - \hat{p}_2)}{m} \right\}
\] がスコア検定統計量で、$\hat{p}_1 = \delta + \hat{p}_2$ となっています。
(3)式は3次方程式となるので $\hat{p}_2$ を解けば $P$ 値が簡単に計算できます（という話も教えていただきました）。
ということなので、(3)式を整理してみると
\[
a \hat{p}_2^3 + b \hat{p}_2^2 + c \hat{p}_2 + d = 0
\] となり、
\[
\begin{aligned}
a &= 1 + \theta, \\
b &= -\{1 + \theta + \theta\tilde{p}_1 + \tilde{p}_2 - \delta(2 + \theta) \}, \\
c &= \delta^2 - \delta(\theta + 2\tilde{p}_2 + 1) + \theta\tilde{p}_1 + \tilde{p}_2, \\
d &= \delta \tilde{p}_2 (1 - \delta),
\end{aligned}
\] となります、ただし、$\theta=n / {m}$, $\tilde{p}_1=x/n$, $\tilde{p}_2 = y / {m}$ です。
これを解くと
\[
\begin{aligned}
p &= \frac{b^2}{3a^2} - \frac{c}{a}, \\
q &= \frac{2b^3}{27a^3} - \frac{bc}{3a^2} + \frac{d}{a}, \\
\hat{p}_2 &= 2\sqrt{\frac{p}{3}} \cos\left\{ \frac{1}{3}\cos^{-1}\left( -\frac{3q}{2p}\sqrt{\frac{3}{p}} \right)+\frac{4\pi}{3} \right\}-\frac{b}{3a},
\end{aligned}
\] となります（導出は飛ばしていますが、3次方程式の解の公式を使って普通に出しています）。
あとは、$\hat{p}_2$ を $T(\delta)$ に代入し、$T(\delta)$ が帰無仮説のもとで漸近的に自由度1の $\chi^2$ 分布に従うことを利用して $P$ 値を得ることができます。
なお、しまりすくんの本で使用しているのは計算が簡単なWald検定です（あそこの説明で上の式を使うのはちょっと厳しいですね）。
Wald統計量は
\[
T_W(\delta)=\frac{(\tilde{p}_1-\tilde{p}_2-\delta)^2}{\frac{\tilde{p}_1(1-\tilde{p}_1)}{n}+\frac{\tilde{p}_2(1-\tilde{p}_2)}{m}}
\] となります。
【新刊案内】
— 朝倉書店 (@AsakuraPub) 2024年2月26日
『宇宙怪人しまりす統計よりも重要なことを学ぶ』2024/2/26発売
あの宇宙怪人が帰ってきた！医療統計を学んできたしまりすに，統計よりも重要なことはあるのか。https://t.co/fnqbMEdpRn
数値例
ツイートの数値例を使って検算してみます。
仮説$H: p_1-p_2= 0$を検定すると、
x <- 58 n <- 58 + 22 y <- 62 m <- 62 + 38 d <- 0 a3 <- 1 + n/m a2 <- -1*(1 + n/m + x/m + y/m - d*(2 + n/m)) a1 <- d^2 - d*(n/m + 2*y/m + 1) + x/m + y/m a0 <- d*y/m*(1 - d) p <- a2^2/(3*a3^2) - a1/a3 q <- 2*a2^3/(27*a3^3) - a2*a1/(3*a3^2) + a0/a3 p2h <- 2*sqrt(p/3)*cos(1/3*acos(-3*q/(2*p)*sqrt(3/p)) + 4*pi/3) - a2/(3*a3) p1h <- d + p2h Td <- ((x - n*p1h)/(p1h*(1 - p1h)))^2*(p1h*(1 - p1h)/n + p2h*(1 - p2h)/m) 1 - pchisq(Td, 1)
実行結果は
> 1 - pchisq(Td, 1) [1] 0.1375639
でした。
Farrington–Manning検定の結果も出して確認してみます。
Farrington–Manning検定はMeeタイプの統計量を二乗しないものを使っていて非劣性の検定として実装されているので、片側検定の結果が出てきます。
そのため、$P$ 値を2倍した値で確認する必要があります。
x <- c(rep(TRUE, 58), rep(FALSE, 22)) y <- c(rep(TRUE, 62), rep(FALSE, 38)) fm <- DescrTab2::farrington.manning(x, y, 0) fm$p.value*2
実行すると
> fm$p.value*2 [1] 0.1375639
となってめでたく一致しました。
$\chi^2$ 検定の方の $P$ 値も確認してみます（載せ忘れていたので追記しました [2024/3/1 10:07]）。
X <- matrix(c(58, 62, 22, 38), ncol = 2) s <- c(colSums(X)[1]/sum(X), colSums(X)[2]/sum(X)) E <- rbind(rowSums(X)[1]*s, rowSums(X)[2]*s) 1 - pchisq(sum((X-E)^2/E), 1)
実行すると
> 1 - pchisq(sum((X-E)^2/E), 1) [1] 0.1375639
こちらも一致します。
各方法と仮説で計算した $P$ 値がどうなるかをまとめると以下になります。
| 検定 | $\delta=0$ | $\delta=0.2$ |
|---|---|---|
| Wald型 | 0.132 | 0.172 |
| スコア型 | 0.138 | 0.164 |
Meeの方法とスコア統計量
佐藤（1988）ではMeeの方法がスコア検定の式に一致することが示されていますが、導出が省略されているので検算した結果をメモしておきます。
Meeの統計量は
\[
T_M(\delta)=\frac{(\tilde{p}_1-\tilde{p}_2-\delta)^2}{\frac{\hat{p}_1(1-\hat{p}_1)}{n}+\frac{\hat{p}_2(1-\hat{p}_2)}{m}}
\] で、これが $T(\delta)$ に等しいことを示します。
$T_M(\delta)$ の分子は、
\[
(\tilde{p}_1-\tilde{p}_2-\delta)^2=
\left( \frac{x}{n} - \hat{p}_1 + \hat{p}_1 - \frac{y}{m} - \delta \right)^2=
\left\{ \left(\frac{x}{n} - \hat{p}_1\right) - \left( \frac{y}{m} - \hat{p}_2 \right) \right\}^2=
\left(\frac{x}{n} - \hat{p}_1\right)^2 -2\left(\frac{x}{n} - \hat{p}_1\right) \left(\frac{y}{m} - \hat{p}_2\right) + \left( \frac{y}{m} - \hat{p}_2 \right)^2
\] と変形できます。
また、(3)式から得られる
\[
\frac{y}{m}-\hat{p}_2=-\frac{\hat{p}_2(1-\hat{p}_2)}{m}\frac{x-n\hat{p}_1}{\hat{p}_1(1-\hat{p}_1)}
\] という関係式を用いれば、$T_M(\delta)$ の分子は、
\[
\left\{ \frac{\hat{p}_1(1-\hat{p}_1)}{n} \right\}^2 \left\{ \frac{x-n\hat{p}_1}{\hat{p}_1(1-\hat{p}_1)} \right\}^2+
2\frac{\hat{p}_1(1-\hat{p}_1)}{n} \frac{\hat{p}_2(1-\hat{p}_2)}{m} \left\{ \frac{x-n\hat{p}_1}{\hat{p}_1(1-\hat{p}_1)} \right\}^2+
\left\{ \frac{\hat{p}_2(1-\hat{p}_2)}{m} \right\}^2 \left\{ \frac{x-n\hat{p}_1}{\hat{p}_1(1-\hat{p}_1)} \right\}^2 =
\left\{ \frac{x-n\hat{p}_1}{\hat{p}_1(1-\hat{p}_1)} \right\}^2 \left\{ \frac{\hat{p}_1(1-\hat{p}_1)}{n} + \frac{\hat{p}_2(1-\hat{p}_2)}{m}\right\}^2
\] となることが分かるので、Meeの統計量とスコア統計量が同じになることが分かります。
参考文献
- Farrington CP, Manning G. Test statistics and sample size formulae for comparative binomial trials with null hypothesis of non-zero risk difference or non-unity relative risk. Statistics in Medicine 1990; 9(12): 1447–1454. DOI: 10.1002/sim.4780091208.
- Mee RW. Confidence bounds for the difference between two probabilities. Biometrics 1984; 40(4): 1175–1176. [Link].
- 佐藤俊哉．コホート研究で用いる疫学指標のefficient scoreにもとづく信頼区間．応用統計学 1988; 17(1): 43–54. DOI: 10.5023/jappstat.17.43.
- 佐藤俊哉．宇宙怪人しまりす 統計よりも重要なことを学ぶ．朝倉書店 2024．[Link].
おわりに
手抜かりがありました、正確な情報をお伝えできるよう精進致します。
