はじめに

昨今はp値よりも効果量を検討することが重要視されている。 特に\(2 \times 2\)行列に集計されたデータの場合には、\(\phi\)係数や odds比は以前から良く知られており、 \(\phi\) 係数と log(odds比)の散布図は特有の境界で句切られ、それぞれに特徴があり、疫学などでは odds比 が重視される。

一般の行列の場合には、Cramer’s Vが効果量として広く使われている。 定義からCramer’s V は非負値をとるが、実は符号を付与することもできる。そして、Cramer’s V は\(\phi\)係数を一般の行列に拡張したものになっている。

それでは、正方行列の場合に限って、odds比の拡張は存在するのだろうか。また、\(\phi\)係数の拡張として、Cramer’s V とは異なるものは存在するのだろうか。性質としては、\(2 \times 2\)の場合と似ていることが望まれる。

\(OR, \Phi\)

命題

複素 \(n \times n\)正則行列からリーマン球への有理写像\(F\) で、

\(F: GL(n,{C}) \rightarrow {C} \cup \{\infty\} \backslash\{1\}\)

次の条件を満たすものが存在する。

条件

  1. 行列\(M\)のある行（又は列）を（0でない）\(c\)倍した行列を\(M_c\)とすれば、
    \(F(M) = F(M_c).\)
  2. 行列\(M\)の2つの行（又は列）を互換した行列を\(M_t\)とすれば、
    \(F(M_t) = 1/F(M_t).\)

実際、後に定義する \(OR\) はこの条件をみたす。
非負実数値行列 \(GL(n,{R_0})\)に制限した場合、

\(OR: GL(n,{R_0}) \rightarrow {R_0} \cup \{\infty\}\backslash\{1\}.\)

であり、その対数をとって \(logOR = log(OR)\)とすれば、

\(logOR: Mat(n,{R_0}) \rightarrow {R} \cup \{\infty\}\backslash\{0\}.\)

  1. \(logOR(M) = logOR(M_c).\)
  2. \(logOR(M_t) = -logOR(M_t).\)

さらに \(\Phi\)も同じく

\(\Phi: Mat(n,{R_0}) \rightarrow {R} \cup \{\infty\}.\)

となるように、次式で定義する。ここで\(det\)は行列式である。

\(\Phi(M) = {det}(M)/\sqrt{prod(rowSums(M),colSums(M))}\)

すると、\(\Phi\) もまた同じ性質を満たす。

  1. \(\Phi(M) = \Phi(M_c).\)
  2. \(\Phi(M_t) = -\Phi(M_t).\)

さらに、修正Cramer’s V として、符号を付与することができる。これは、\(V\)が0に近い 値をとるときに、有効な場合がある。

　\(CrV(M) = sign(det(M))*CramerV(M).\)

library(vcd)
## Loading required package: grid
PhiV <- function(x, signed=FALSE){
   dt <- det(x)
   c(dt/sqrt(prod(rowSums(x),colSums(x))),
   if(signed) assocstats(x)$cramer * sign(dt)
   else assocstats(x)$cramer)
}

OR の定義

対称式、交代式、行列式

\(K\)を体とする。\(K\)上の\(n\)次元ベクトル空間\(V=K^n\)のaffine座標環は\(K[V]=K[x_1,x_2,...,x_n]\)である。 この添字の集合 \(I=1:n\) には、\(n\)文字の置換群 \(S(n)\)が作用する。 その作用が誘導する\(K[V]\)上の作用による不変式環 \(S[V]=K[V]^{S(n)}\) は対称式からなる。 \(S(n)\)には index 2 の部分群として、偶置換からなる\(A(n)\)が存在し、その不変式環\(A[V]=K[V]^{A(n)}\)は、Sの２次拡大となる（なぜなら[S(n):A(n)]=2）。 その生成元は、\(K\)の標数が２でないならば、基本差積

\[\Delta = \prod_{i<j}(x_j - x_i)\]

である。\(K\)の標数が２のときには、

\[\Delta' = \sum_{\sigma \in A(n)} x_{\sigma(1)}^{n-1}x_{\sigma(2)}^{n-2}...x_{\sigma(n)}^{2}x_{\sigma(n)}\]

を採用する。いずれも\(S[V]=A[V]+\Delta A[V]\), \(S[V]=A[V]+\Delta' A[V]\) となる。

なお、置換の符号とよばれる関数\(sign: S(n) \rightarrow \{1,-1\}\)があり、乗法群としての準同型である。このとき \(A(n) = Ker(sign)\)　となる。

行列環

次に、行列環\(M(n,K)\)を\(K\)上の\(n^2\)次元ベクトル空間として考察する。座標環は \(K[M]=K[x_{ij}], i,j \in I\) である。\(f \in K[M]\)に対して、\(\sigma \in S(n)\)の作用を、 \(f^\sigma(x_{ij})=f(x_{i\sigma(j)})\) と定める。

正方行列\(X=(x_{ij})\)の行列式\(det(X)\)を、

\(det(X) = \sum_{\sigma \in S(n)} sign(\sigma) x_{1\sigma(1)}x_{2\sigma(2)}...x_{n\sigma(n)}\)

により定める。これを２つに分割すれば、

\(detp(X) = \sum_{\sigma \in A(n)} x_{1\sigma(1)}x_{2\sigma(2)}...x_{n\sigma(n)},\)

\(detn(X) = \sum_{\sigma \in (12)A(n)} x_{1\sigma(1)}x_{2\sigma(2)}...x_{n\sigma(n)},\)

\(det(X) = detp(X) - detn(X).\)

\(n=2: det(X)=x_{11}x_{22} - x_{12}x_{21}\)

\(n=3: det(X)=x_{11}x_{22}x_{33} +x_{12}x_{23}x_{31} + x_{13}x_{21}x_{32} - (x_{13}x_{22}x_{31} + x_{12}x_{21}x_{33} + x_{11}x_{23}x_{32})\)

以上の事実を背景として、つぎのように\(OR\) を定める。

定義 & 命題　拡張oddsratio

\(OR(M) = detp(M) / detn(M), M \in GL(n,R_0)\)

\(OR\) は 条件 1,2 を満たす。（証明はstraightforward）

\(OR\)の計算

\(OR\) は通常の行列式の関数では計算できず、\(detp\) と \(detn\)を別個に求めねばならない。この各々の直接計算は\(n!/2\)個の,\(n\)個の積からなり、計算時間がかかるのみならず、巨大な数値となるので、多倍長計算も必要となる。スクリプトは再帰的に構成するのがよい。たとえば rcsv.det のとおりである。これによって

detp(x)=pn[1], detn(x)=pn[2], det(x)=pn[1]-pn[2]

となる。ただし、有効数字15桁までで近似計算が行われる。
正確に行う場合は、Rmpfrを使うこととなり、psはRmpfr object となる。下記の例では、k=100(bits)を宣言しており、10進で30桁まで中途計算を保証する。それらの割り算である、odds.ratioもRmpfr objectである。それを、通常の数値 num.odds.ratio に変換している。
行列のサイズが小さい場合は、rcsv.det で十分である。 なお、Rでは 1/0 はInf, 0/0 は NaN を与える。特に try()で error 処理をする必要はない。
package:Rmpfr をloadすると下記のメッセージがでる。以下では loadしていないので、参考のため記す。 ——-
Loading required package: gmp Attaching package: ‘gmp’
以下のオブジェクトは ‘package:base’ からマスクされています: %*%, apply, crossprod, matrix, tcrossprod
C code of R package ‘Rmpfr’: GMP using 64 bits per limb
Attaching package: ‘Rmpfr’
以下のオブジェクトは ‘package:stats’ からマスクされています: dbinom, dnorm, dpois, pnorm
以下のオブジェクトは ‘package:base’ からマスクされています: cbind, pmax, pmin, rbind
——-

rcsv.det <- function(x)  # x: matrix
{
  n <- nrow(x)
  pn <- c(0,0)
  if (n==2) {
    pn[1] <- x[1,1]*x[2,2]
    pn[2] <- x[1,2]*x[2,1]
    return(pn)
  } else {
    indx <- rep(c(1,2),n)
    for (cnt in 1:n){
      pn[1] <- pn[1]+x[1,cnt]*rcsv.det(x[-1,-cnt])[indx[cnt]]
      pn[2] <- pn[2]+x[1,cnt]*rcsv.det(x[-1,-cnt])[indx[cnt+1]]
    }
    return(pn)
  }
}

odds.ratio <- function(x)
{
  pn <- rcsv.det(x)
  try(pn[1]/pn[2])
}


# library(Rmpfr)
k <- 100L
rcsvm.det <- function(x)  # x: matrix
{
  n <- nrow(x)
  pn <- c(mpfr(0,k),mpfr(0,k))
  if (n==2) {
    pn[1] <- mpfr(x[1,1],k)*mpfr(x[2,2],k)
    pn[2] <- mpfr(x[1,2],k)*mpfr(x[2,1],k)
    return(pn)
  } else {
    indx <- rep(c(1,2),n)
    for (cnt in 1:n){
      pn[1] <- pn[1] + mpfr(x[1,cnt],k)*rcsvm.det(x[-1, -cnt])[indx[cnt]]
      pn[2] <- pn[2] + mpfr(x[1,cnt],k)*rcsvm.det(x[-1, -cnt])[indx[cnt+1]]
    }
    return(pn)
  }
}

num.odds.ratio <- function(x)
{
  pn <- rcsvm.det(x)
  odds.ratio <- try(pn[1]/pn[2])
  as.numeric(odds.ratio)
}

実例

実例をみれば、計算の簡単な \(\Phi\)と面倒な\(logOR\)とがよく相関することがわかる。

実験　marginal を固定したランダム行列

rs <- c(50,50,60,60)
cs <- c(60,60,60,40)
m <- 200
matr <- r2dtable(m, rs, cs)
ans <- matrix(0, nrow=m, ncol=3)
for (cnt in 1:m){
  PV <- PhiV(matr[[cnt]])
  ans[cnt,1] <- PV[1]    # Phi 
  ans[cnt,2] <- PV[2]    # V
  ans[cnt,3] <- log(odds.ratio(matr[[cnt]]))
}
par(mfrow=c(1,3))
plot(ans[,2], ans[,1],  col=1,main=expression(Phi), xlab="V")
plot(ans[,2], ans[,3], col=2, pch=0 , main="log odds.ratio", xlab="V")
plot(ans[,1], ans[,3], col=4, pch=2 , main="log odds.ratio", xlab=expression(Phi))

par(mfrow=c(1,1))

実験　固定した行列の近傍からなるランダム行列

中心は次の通り centre <- matrix(c(30,10,15,10,10,25,10,10, 10,10,30,10,15,15,15,20), 4)

## Loading required package: manipulate

対角行列のランダム近傍

diag(c(30,30,30,30)) を中心にして 行列成分に[0,1]のランダム一様分布の20倍を加える

成分を1~500 からランダムに