【統計】コクラン-アーミテージ検定（Cochran-Armitage test）

Cochran-Armitage test（コクラン-アーミテージ検定またはCA検定）について。

【目次】

帰無仮説、対立仮設
計算式等
計算例
プログラムコード
その他
参考

CA検定は、割合の増加/減少傾向を評価する傾向性の検定。
検定の中身は、割合を目的変数とした回帰分析であり、回帰直線の傾きが０であるかどうかを見ている。
なお、目的変数が連続値の場合は、Jonckheere-Terpstra trend test（ヨンクヒール-タプストラ検定、またはヨンキー検定）等が用いられる。

【統計】ヨンキー検定（傾向性検定） - こちにぃるの日記

帰無仮説、対立仮設

帰無仮説 $H_{0}: p_{1}=p_{2}=...=p_{k}$ （直線の傾きが０）
対立仮設 $H_{1}: p_{1} \leqq p_{2} \leqq ... \leqq p_{k}$ （直線の傾きが０でない）

計算式等

CA検定では $\chi^{2}$ 統計量を用いる。
通常の $\chi^{2}$ 検定の統計量を $\chi _{Total}^{2}$ 、CA検定で扱う統計量を $\chi _{Slope}^{2}$ 、残差を $\chi _{Resid}^{2}$ とすると次の式が成り立つ。

$\quad \chi _{Total}^{2} = \chi _{Slope}^{2} + \chi _{Resid}^{2}$

これら $\chi^{2}$ 統計量の導出は下表のとおり。

	$\chi ^{2}$ 統計量	自由度(df)
Slope	$\chi _{Slope}^{2}= S_{sp}^{2}/S_{ss}$	$1$
Residual	$\chi _{Resid}^{2}=\chi _{Total}^{2} - \chi _{Slope}^{2}$	$k-1-1$
Total	$\chi _{Total}^{2}=S_{pp}$	$k-1$

表の平方和は次の計算式。

$\quad S_{ss}= \sum \{ n_{i} w \left( s_{i} - \overline{s} \right) ^{2} \}$

$\quad S_{pp}= \sum \{ n_{i} w \left( p_{i} - \overline{p} \right) ^{2} \}$

$\quad S_{sp}= \sum \{ n_{i} w \left( s_{i} - \overline{s} \right) \left( p_{i} - \overline{p} \right) \}$

$\quad \overline{p}= \dfrac{A}{N}= \dfrac{\sum a_{i}}{\sum n_{i}},\quad p_{i}= \dfrac{a_{i}}{n_{i}}$

$\quad \overline{s}= \dfrac{\sum \left( n_{i} s_{i} \right)}{N}$

$\quad w= \dfrac{1}{\overline{p} \left( 1-\overline{p} \right)}$

※ $s_{i}$ =当該水準でのスコア、 $n_{i}$ =当該水準での全体数、 $a_{i}$ =当該水準での該当数

計算例

簡単なデータで計算方法を確認。
今回は、病気Xの病状（Score）が悪いほど死亡発生確率（Yes/Total）が増加する傾向があるかを見てみる（架空データ）。
※計算の理解のためScoreにも値を指定。特に指定がない場合は1,2,3,4...と等間隔な値を仮定する。

■ 観測度数

Score	Yes	No	Total	Prop.
1	10	100	110	0.091
2	20	80	100	0.200
3.5	20	60	80	0.250
5.5	20	40	60	0.333

先ず、全体の基本情報；
例数 $N＝110+100+80+60=350$ 、死亡数 $A=10+20+20+20=70$ 、死亡割合 $\overline{p}＝70/350=0.2$

CA検定では、割合を目的変数とした回帰分析を行うが、手法として割合の標準誤差の逆数を重みとした「重み付き最小二乗法」を用いる。この重みと各スコアの重み付きN数は次のとおり。
$\quad w=\dfrac{1}{\overline{p}(1-\overline{p})}=\dfrac{1}{0.2×0.8}=6.25$
$\quad \verb|Score 1 : | n_{1} w=110 \times 6.25 = 687.5$
$\quad \verb|Score 2 : | n_{2} w=100 \times 6.25 = 625.0$
$\quad \verb|Score 3.5 : | n_{3} w=80 \times 6.25 = 500.0$
$\quad \verb|Score 5.5 : | n_{4} w=60 \times 6.25 = 375.0$

重みを使った回帰直線：

f:id:cochineal19:20210207155857p:plain

　※プロットサイズ：各スコアの重み付き例数
　※縦の点線：全体の平均スコア $\overline{s}= \sum \left( n_{i} s_{i} \right)/N$
　※横の点線：全体の確率 $\overline{p}= A/N$

また、CA検定で扱う偏差は次のもの。これら偏差の平方和を用いて $\chi^{2}$ 統計量を求める。※下図の緑の両矢印が偏差のイメージ。

　確率の偏差　：個々のスコアでの確率 $p_{i}$ - 全体の確率 $\overline{p}$
　スコアの偏差：個々のスコア $s_{i}$ - 全体の平均スコア $\overline{s}$

f:id:cochineal19:20210208093005p:plain

$\quad S= \dfrac{\sum \left( n_{i} s_{i} \right)}{N}=\dfrac{110\times 1 + 100\times 2 + 80\times 3.5+60\times 5.5}{350}=2.629$
$\quad \overline{p}= \dfrac{A}{N}=\dfrac{70}{350}=0.2$
$\quad S_{ss}= \sum \{ n_{i} w \left( s_{i} - \overline{s} \right) ^{2} \}=687.5\left(1-2.629 \right)^{2}+625\left(2-2.629 \right)^{2}+500\left(3.5-2.629 \right)^{2}+375\left(5.5-2.629 \right)^{2}=5541.964$
$\quad S_{pp}= \sum \{ n_{i} w \left( p_{i} - \overline{p} \right) ^{2} \}=687.5\left(0.091-0.2 \right)^{2}+625\left(0.2-0.2 \right)^{2}+500\left(0.25-0.2 \right)^{2}+375\left(0.333-0.2 \right)^{2}=16.098$
$\quad S_{sp}= \sum \{ n_{i} w \left( s_{i} - \overline{s} \right) \left( p_{i} - \overline{p} \right) \}=687.5\left(1-2.629 \right)\left(0.091-0.2 \right)+625\left(2-2.629 \right)\left(0.2-0.2 \right)+500\left(3.5-2.629 \right)\left(0.25-0.2 \right)+375\left(5.5-2.629 \right)\left(0.333-0.2 \right)=287.5$

	$\chi ^{2}$ 統計量	自由度(df)	p値
Slope	$\chi _{Slope}^{2}=S_{sp}^{2}/S_{ss}$ $=287.5^{2}/5541.964$ $=14.9146$	$1$	$0.000112$
Residual	$\chi _{Resid}^{2}=\chi _{Total}^{2} - \chi _{Slope}^{2}$ $= 16.0985-14.9146$ $= 1.1839$	$k-1-1$ $=4-1-1$ $=2$	$0.553255$
Total	$\chi _{Total}^{2}=S_{pp}$ $=16.098$	$k-1$ $=4-1$ $=3$	$0.001082$

上の表のうち、 $\chi_{Slope}^{2}$ が回帰直線の傾きに関連する統計量であり、この統計量を使った $\chi^{2}$ 検定がCA検定。
$\chi_{Resid}^{2}$ は残差項。直線からのズレが有意であるかを示すため、非有意が望まれる。残差項が有意な場合、目的変数と説明変数の関係性を直線で表しきれていないことを意味する。
$\chi_{Total}^{2}$ は通常の $\chi^{2}$ 検定の結果。

今回の架空データでは、病気Xの病状（Score）が悪いほど死亡発生確率が増加する傾向がある（p for trend<0.001）と結論付けられそう。また、残差項は有意でないため、直線モデルが当てはまっていると解釈できる。

Excelでの計算結果は次のとおり。

f:id:cochineal19:20210208092624p:plain

Rでの実行：

> mtx <- matrix(c(10,20,20,20,100,80,60,40), nrow=4, ncol=2)
> chisq.test(mtx, correct = F)
	Pearson's Chi-squared test
data:  mtx
X-squared = 16.098, df = 3, p-value = 0.001082
> prop.trend.test(
+   x=mtx[,1]
+   ,n=mtx[,1]+mtx[,2]
+   ,score = c(1,2,3.5,5.5))
	Chi-squared Test for Trend in Proportions

data:  mtx[, 1] out of mtx[, 1] + mtx[, 2] ,
 using scores: 1 2 3.5 5.5
X-squared = 14.915, df = 1, p-value = 0.0001125

SASでの実行：

data ads;
input SCORE OUTCOME CNT @@;
datalines;
1 1 10 1 2 100 2 1 20 2 2 80 3.5 1 20 3.5 2 60 5.5 1 20 5.5 2 40
;run;
proc freq data=ads;
  table SCORE * OUTCOME / trend nocol nopercent;
  weight CNT / zeros;
  output out=OutputDS trend;
run;

f:id:cochineal19:20201219141532p:plain

$z^{2}=\chi ^{2} \left( df=1 \right)$ の関係があるため、 $z^{2}≒3.8619^{2}≒14.9146$ で等価。

プログラムコード

■ Rのコード

prop.trend.test(x=ADS$x, n=ADS$n, score=ADS$score)

■ SASのコード

proc freq data=[InputDS];
  table Factor * Outcome / trend nocol nopercent;
  weight Count / zeros;
  output out=[OutputDS] trend; /*_TREND_, P2_TRENDを取得。*/
run;

■ Pythonのコード

整備中

その他

その１：

SASの出力結果のところで少し触れたが、統計ソフトによって取り扱う統計量が $\chi^{2}$ のものとz値のものがある。
$\chi^{2}$ 分布（自由度＝１）には $\chi^{2}=z^{2}$ の関係があるため、p値はどちらも同じ。

備忘録がてら z値の導出式も掲載。
平方和の式に重みはなく、z値の導出時に重みを与えている。

$\quad S_{sp}= \sum \{ n_{i} \left( s_{i} - \overline{s} \right) \left( p_{i} - \overline{p} \right) \}$

$\quad S_{ss}= \sum \{ n_{i} \left( s_{i} - \overline{s} \right) ^{2} \}$

$\quad z=\dfrac{\sum n_{i}\left( p_{i}-\overline{p}\right) \left( s_{i}-\overline{s}\right) }{\sqrt{\overline{p}\left( 1-\overline{p}\right) \sum n_{i}\left( s_{i}-\overline{s}\right) }}=\dfrac{S_{sp}}{\sqrt{\overline{p}\left( 1-\overline{p}\right) S_{ss}}}=\dfrac{S_{sp}^2}{\overline{p}\left( 1-\overline{p}\right) S_{ss}}=w\dfrac{S_{sp}^2}{S_{ss}}$

本記事の架空データで試すと次のとおり。

$\quad S_{sp}= \sum \{ n_{i} \left( s_{i} - \overline{s} \right) \left( p_{i} - \overline{p} \right) \}=110\left(1-2.629 \right)\left(0.091-0.2 \right)+100\left(2-2.629 \right)\left(0.2-0.2 \right)+80\left(3.5-2.629 \right)\left(0.25-0.2 \right)+60\left(5.5-2.629 \right)\left(0.333-0.2 \right)=46$

$\quad S_{ss}= \sum \{ n_{i} \left( s_{i} - \overline{s} \right) ^{2} \}=110\left(1-2.629 \right)^{2}+100\left(2-2.629 \right)^{2}+80\left(3.5-2.629 \right)^{2}+60\left(5.5-2.629 \right)^{2}≒886.714$

$\quad w=\dfrac{1}{\overline{p}(1-\overline{p})}=\dfrac{1}{0.2×0.8}=6.25 \verb|　＜再掲＞|$

$\quad z=w\dfrac{S_{sp}^2}{S_{ss}}=6.25\times\dfrac{46^2}{886.714}≒3.86194$

この $z≒3.86194$ の p値を求めると $p≒0.0001124$ になる。
$\quad \verb|= 2 * (1 - NORM.S.DIST(ABS(3.86194), TRUE))　　※Excel関数|$

---------------
その２

Rのlm関数のweightオプションに本記事で取り上げた「各スコアの重み付きN数」を指定すれば、重み付き最小二乗法での回帰分析を実行できる。
この結果から分散分析表を作れば、 $\chi^{2}_{slope}$ 統計量を導出できる。

> library(tidyverse)
> x <- c(10,20,20,20)
> n <- c(110,100,80,60)
> score <- c(1,2,3.5,5.5)
> (p1 <- sum(x)/sum(n))
> (p2 <- x/n)
> (w1 <- 1/(p1*(1-p1)))
> (w2 <- n*w1)
[1] 0.2
[1] 0.09090909 0.20000000 0.25000000 0.33333333
[1] 6.25
[1] 687.5 625.0 500.0 375.0
> (fit <- lm(p2 ~ score, weights=w2))
Call:
lm(formula = p2 ~ score, weights = w2)
Coefficients:
(Intercept)        score  
    0.06364      0.05188  

> anova(fit)
Analysis of Variance Table
Response: p2
          Df  Sum Sq Mean Sq F value  Pr(>F)  
score      1 14.9146 14.9146  25.196 0.03747 *
Residuals  2  1.1839  0.5919                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> #-- 例示した回帰直線のグラフはこちら。
> ads <- data.frame(score=score,p=p2,w=w2)
> ggplot(ads, aes(x=score, y=p)) +
+   geom_vline(xintercept=2.629, lty=2) +
+   geom_hline(yintercept=0.200, lty=2) +
+   geom_point(aes(size=w), colour="blue", alpha=.5) +
+   stat_smooth(method="lm", se=F, aes(weight=w))

参考

統計学（出版：東京図書）, 日本統計学会編

https://www.lexjansen.com/pharmasug/2007/sp/SP05.pdf

統計学入門−第5章