「統計学関連なんでもあり」の過去ログ--- 048

No.22612　ゼロでの割り算　　【青木繁伸】　2018/09/09(Sun) 11:10

中澤さんが，2x2 より大きい行列で R の mcnemar.test の結果が NaN になる場合があることについて，書いている。

http://minato.sip21c.org/bulbul2/20180908.html
http://minato.sip21c.org/bulbul2/20160913.html
http://minato.sip21c.org/bulbul2/20160911.html

「対角成分の和がゼロの組合せを計算からスキップできる拡張マクネマー検定の関数」として

mct <- function(x) {
  L <- NROW(x)
  X <- 0
  N <- 0
  for (i in 1:(L - 1)) {
    for (j in (i + 1):L) {
      s <- (x[i, j] + x[j, i])
      if (s > 0) {
        N <- N + 1
        X <- X + (x[i, j] - x[j, i]) ^ 2 / s
      }
    }
  }
  return(c( X2 = X, df = N, p = (1 - pchisq(X, N))))
}

を挙げている。
しかし，「このやり方が正しいという理論的根拠を見つけられないので悩ましいところだが」とあるのだが。

「N23とN32がゼロ，という状態だと，(N23+N32)が分母になる項で"Division by zero"エラーがでてしまって，カイ二乗値の計算結果自体がNaNになってしまう」そこで「分母がゼロになる組合せをスキップする」ということを書いている。

この記述は正しくない。N23=N32=0のとき，つまり(0-0)/(0+0) = 0/0 の計算結果は，NaN である（"Division by zero"ではない。Python だと，ZeroDivisionError: division by zero になる）。何かに NaN を足すと NaN になり，それ以降 NaN に何を足しても NaN のまま。ということで，R の mcnemar.test は NaN を返すと，中澤さんも書いている。

N23=N32=0 のときの計算結果が NaN になった場合は足し算をスキップすればよいわけで，N23=N32=0 の場合は計算をスキップしても同じではある。またその場合には自由度のカウントはしない。

以下のような関数を書いてみた。

mct2 <- function(x) {
  i <- col(x)[lower.tri(x)]
  j <- row(x)[lower.tri(x)]
  a <- mapply(function(i, j) (x[i, j] - x[j, i]) ^ 2 / (x[i, j] + x[j, i]), i, j)
  chisq <- sum(a, na.rm = TRUE) # NaN も除かれる
  df <- sum(!is.nan(a)) # !is.na(a) でもよい
  return(c(X2 = chisq, df = df, p.value = pchisq(chisq, df, lower.tail = FALSE)))
}

「ゼロ割」は避けるのではなく，積極的に正しく対処するべきということか？（同じようなことではあるが）

蛇足ではあるが，自由度を減ずるのは，カイ二乗分布の再生性から明らか（自由度 1 のカイ二乗分布にしたがう確率変数 k 個の和は自由度 k のカイ二乗分布にしたがう）。

R の mcnemar.test を修正するとすれば，

    STATISTIC <- sum(y[upper.tri(x)]^2/x[upper.tri(x)])

を

    STATISTIC <- sum(y[upper.tri(x)]^2/x[upper.tri(x)], na.rm=TRUE)
    PARAMETER <- sum(!is.nan(y[upper.tri(x)]^2/x[upper.tri(x)]))

にする（前の方の PARAMETER <- r * (r - 1)/2 は不要）。
プログラムにおいて，特殊な条件の発生を見逃すと命取りというところ。
しかしさすがに，sum(y[upper.tri(x)]^2/x[upper.tri(x)]) とは，にくいねえ。

> x = matrix(c(c(7, 2, 5, 3, 11, 2, 0, 3, 4, 0, 5, 9, 3, 15, 14, 10)), 4, 4)
> mct(x)
          X2           df            p 
15.428836864  5.000000000  0.008678852 
> mct2(x)
          X2           df      p.value 
15.428836864  5.000000000  0.008678852 

> mcnemar.test(x) # 元のまま

	McNemar's Chi-squared test

data:  x
McNemar's chi-squared = NaN, df = 6, p-value = NA

> mcnemar.test2(x) # 修正したもの

	McNemar's Chi-squared test

data:  x
McNemar's chi-squared = 15.429, df = 5, p-value = 0.008679

なお，「"Division by zero"を避けるために全部のセルに0.01を足す」というのは明らかに間違い（自由度が減らない）。

No.22613　Re: ゼロでの割り算　　【青木繁伸】　2018/09/09(Sun) 12:49

2 0
0 5

というのと

2 1
1 5

というのは，0-0 = 1-1 = 0 という情報を与えるという意味では同じなので，自由度は減らさなくてもよいか
STATISTIC <- sum(y[upper.tri(x)]^2/x[upper.tri(x)], na.rm=TRUE)
だけの変更で，これは実質的には NaN の場合に 0 を加えるというのと同じであるから
PARAMETER はそのまま
PARAMETER <- r * (r - 1)/2
とすれば，mcnemar.test(x+0.01) もありか。ただし，加える数は他のまともなセルのことも考えるとできるだけ小さい方が良いだろうとは思う。例えば +1e-15 とか

No.22614　Re: ゼロでの割り算　　【中澤】　2018/09/09(Sun) 21:30

記事ありがとうございます。
確かにNaNが正しいですね。不正確な記述でした。
自由度の扱いは悩ましいですね。スキップでなくて分子がゼロなのでゼロとみなすということですね。勉強になりました。

No.22615　Re: ゼロでの割り算　　【中澤】　2018/09/19(Wed) 17:55

http://minato.sip21c.org/bulbul2/20180915.html
に書きましたが，カテゴリ数が3以上の場合の拡張マクネマー検定と同様の帰無仮説を検定する方法としてBhapkarの検定というのがあって，irrパッケージのbhapkar()で実行できました。いま書いている論文ではこれを使おうと思います。