目的

コクラン・アーミテージ検定を行う
R には prop.trend.test という関数として用意されている。
上に示した結果のうち，傾向検定の結果(Trend の項目）のみが得られる。

注：コクラン・アーミテージ検定のトレンドを表すカイ二乗値(CA)と
Linear-by-Linear 検定（Mantel の傾向検定）のカイ二乗値(M)には，全サンプルサイズを n としたとき，
M = CA*(n-1)/n の関係がある。

使用法

Cochran.Armitage(r.i, n.i, x.i=1:length(r.i))

引数

r.i	反応数のベクトル
n.i	ケース数のベクトル
X.i	外的基準。省略されたときは1から始まる整数値が仮定される

ソース

インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/Cochran_Armitage.R", encoding="euc-jp")

# コクラン・アーミテージ検定
Cochran.Armitage <- function(        r.i,                    # 反応数のベクトル
                                n.i,                    # ケース数のベクトル
                                x.i=seq(along=r.i))     # 外的基準。省略されたときは1から始まる整数値が仮定される
{
        k <- length(r.i)                             # 群の数
        stopifnot(length(n.i) == k, length(x.i) == k)

        p.i <- r.i/n.i                                       # 標本比率のベクトル
        r <- sum(r.i)                                        # 反応数の合計
        n <- sum(n.i)                                        # ケース数の合計
        p.m <- r/n                                   # プールした標本比率（標本比率の平均値）
        x.m <- sum(n.i*x.i)/n                                # 外的基準の平均値
        xx <- x.i-x.m                                        # 平均偏差
        b <- sum(n.i*(p.i-p.m)*xx)/sum(n.i*xx^2)     # 傾き
        a <- p.m-b*x.m                                       # 切片
        xt <- b^2*sum(n.i*xx^2)/(p.m*(1-p.m))                # 傾き（トレンド）に対するカイ二乗値
        xh <- n^2*(sum(r.i^2/n.i)-r^2/n)/r/(n-r)     # 非一様性
        xq <- xh-xt                                  # 直線からの乖離に対するカイ二乗値
        chisq <- c(xt, xq, xh)                               # カイ二乗値
        df <- c(1, k-2, k-1)                         # 自由度
        P <- pchisq(chisq, df, lower.tail=FALSE)     # P 値
        res <- data.frame(chisq=chisq, df=df, P=P)
        colnames(res) <- c("カイ二乗値", "自由度", "P 値")
        rownames(res) <- c("トレンド", "直線からの乖離", "非一様性")
        return(res)
}


インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/wilcoxon_score.R", encoding="euc-jp")

# ウィルコクソン・スコア（平均順位）
wilcoxon.score <- function(n)                # カテゴリーの度数ベクトル
{
        a <- cumsum(n)
        a <- c(0, a[-length(a)])
        return(a+(n+1)/2)               # カテゴリーのウィルコクソン・スコア
}


使用例

x.i <- c(10, 20, 30, 40, 50)   # 各群の外的基準変数の値
n.i <- c(30, 35, 47, 21, 45)   # 各群のケース数
r.i <- c(2, 4, 14, 13, 39)     # 各群の反応ケース数

Cochran.Armitage(r.i, n.i, x.i)	  # 例題の場合は，第三引数を省略しても同じ結果になる
( score <- wilcoxon.score(n.i) )
Cochran.Armitage(r.i, n.i, score) # 外的基準としてウィルコクソン・スコア（平均順位）を使う

出力結果例

> Cochran.Armitage(r.i, n.i, x.i)
               カイ二乗値 自由度         P 値
トレンド        68.572732      1 1.222834e-16
直線からの乖離   4.014484      3 2.599043e-01
非一様性        72.587216      4 6.449422e-15

> score <- wilcoxon.score(n.i)
[1]  15.5  48.0  89.0 123.0 156.0 # ウィルコクソン・スコア

> Cochran.Armitage(r.i, n.i, score)
               カイ二乗値 自由度         P 値
トレンド        68.085779      1 1.565355e-16
直線からの乖離   4.501437      3 2.121622e-01
非一様性        72.587216      4 6.449422e-15
注

R には prop.trend.test という関数として用意されている。
上に示した結果のうち，傾向検定の結果(Trend の項目）のみが得られる。

> prop.trend.test(r.i, n.i) # デフォルトのスコアを使う

	Chi-squared Test for Trend in Proportions

data:  r.i out of n.i ,
 using scores: 1 2 3 4 5 # デフォルトのスコア
X-squared = 68.5727, df = 1, p-value < 2.2e-16

> prop.trend.test(r.i, n.i, score) # 使用するスコアを指定する

	Chi-squared Test for Trend in Proportions

data:  r.i out of n.i ,
 using scores: 15.5 48 89 123 156 # 使用されたスコア（ウィルコクソン・スコア）
X-squared = 68.0858, df = 1, p-value < 2.2e-16

 解説ページ