Rを使用してmatrixクラスのデータから期待値を算出する

はじめに

カイ二乗検定では「期待値が5未満のセルが全体の20％以上になってはいけない」というルールがあります（生物科学研究所・井口研究室のページなどを参考）。このルールは「コクランの規則」と呼ばれます。

この論文を以下に示します（ただし有料）。

Cochran W.G. (1954)
Some methods for strengthening the common x2 tests
Biometrics 10: 417-451.

ここでmatirixクラスのデータがある場合に、その期待値を求める関数があるかを調べたのですが、見つけることができませんでした（もしかしたら調べ方が悪いだけかもしれませんが）。そこで、調べた結果をもとに関数を作成しました。

方法

カイ二乗検定はその計算過程で期待値を求めます。したがって、カイ二乗検定を行うための関数chisq.testのソースコードを調べてみました。

chisq.testのソースコード（抜粋）

Rのソースコード（R-4.0.0）をダウンロードして解凍します。
src/library/stats/Rフォルダにchisq.test.Rがあります。

該当する部分を抜粋します。

n <- sum(x)                  # すべてのマスの数値の合計
sr <- rowSums(x)             # 行ごとの合計
sc <- colSums(x)             # 列ごとの合計
E <- outer(sr, sc, "*") / n  # 期待値の算出（外積を使用）

outer関数を使用することで非常に簡潔な式で期待値を算出していることが分かりました。outer関数についてはRjpWikiのこの記事 や この記事 を参考にして下さい。

このスクリプトをもとにしてもっと簡単にしたら以下のようになります。

E <- outer(rowSums(x), colSums(x), "*") / sum(x)

実例

使用するデータ

以下のデータを例にして期待値を求めます。

	c1	c2	c3	Sum
r1	7	8	5	20
r2	8	12	20	40
Sum	15	20	25	60

x <- matrix(c(7, 8, 5, 8, 12, 20), 2, 3, byrow = TRUE)
rownames(x) <- c("r1", "r2")
colnames(x) <- c("c1", "c2", "c3")

期待値の算出

n <- sum(x)       # すべてのマスの数値の合計
sr <- rowSums(x)  # 行ごとの合計
sc <- colSums(x)  # 列ごとの合計

これらの数値を確認します。

n; sr; sc

## [1] 60

## r1 r2 
## 20 40

## c1 c2 c3 
## 15 20 25

期待値を算出します。

E <- outer(sr, sc, "*") / n  # 期待値の算出
E

##    c1        c2        c3
## r1  5  6.666667  8.333333
## r2 10 13.333333 16.666667

E <- sr %o% sc / n  # %o%演算子でも同じ結果となります
E

##    c1        c2        c3
## r1  5  6.666667  8.333333
## r2 10 13.333333 16.666667

関数化したもの

以上のスクリプトをもとに期待値を返す関数を作成しました。matrixクラス以外の場合にはエラーが出るようにしました。

expected_value <- function(x) {
    if (!is.matrix(x)) {
        stop("x must be 'matrix' class")
    }

    E <- rowSums(x) %o% colSums(x) / sum(x)
    rownames(E) <- rownames(x)
    colnames(E) <- colnames(x)
    return(E)
}

expected_value(x)

##    c1        c2        c3
## r1  5  6.666667  8.333333
## r2 10 13.333333 16.666667

解説

数式を用いた説明

以下のような行列を考えます。例として$2 \times 3$の行列とします。

	c1	c2	c3	合計
r1	$a$	$b$	$c$	$m1$
r2	$d$	$e$	$f$	$m2$
合計	$n1$	$n2$	$n3$	$n$

このとき[r1, c1]のセルの期待値を求めるには

$m1 \times \frac{n1}{n} = \frac{m1 \times n1}{n}$

となります。他のセルでも同様に計算をして求めると、期待値は以下の表のようになります。

	c1	c2	c3	合計
r1	$m1 \times n1 / n$	$m1 \times n2 / n$	$m1 \times n3 / n$	$m1$
r2	$m2 \times n1 / n$	$m2 \times n2 / n$	$m2 \times n3 / n$	$m2$
合計	$n1$	$n2$	$n3$	$n$

Rの関数による説明

• rowSums関数は行ごとの合計（すなわち[m1, m2]）
• colSums関数は列ごとの合計（すなわち[n1, n2, n3]）
• outer関数および%o%演算子は変数同士の外積

[m1, m2]と[n1, n2, n3]の外積による表（rowSums(x) %o% colSums(x)）

	c1	c2	c3
r1	$m1 \times n1$	$m1 \times n2$	$m1 \times n3$
r2	$m2 \times n1$	$m2 \times n2$	$m2 \times n3$

この結果を$n$で割ると（rowSums(x) %o% colSums(x) / n）、上の表で示した期待値となります。

おわりに

Rの関数のソースコードを読むことで新たな発見がありました。
また、outer関数のヘルプを読んで初めて%o%演算子の存在を知りました。