多次元球ベースで2xM表の重み付け検定処理を書く - ryamadaの遺伝学・遺伝統計学メモ

６日目。...一昨昨日、一昨日と昨日(記事はこちらとこちらとこちらとこちらとこちら)の続き
2xM表のM-1自由度の独立性検定カイ自乗値が等しいテーブルが多次元球をなすようなM=df+1次元ユークリッド空間を考える
$\sum_{i=1}^M d_i=0$ ( $d_i$ は第１行第i列のセルの値とその期待数との差)なる制約がある。これは自由度dfに制約された部分空間
1...M列に $w_1,...,w_M$ なる重み付けをした自由度１の検定は、 $\sum_{i=1}^M d_iw_i$ の値が等しいとき、同じ統計量を与える。
その統計量は、 $W=(\sqrt{N_i} w_i)$ (ただし $N_i$ は第i列の周辺度数)というベクトルに垂直な面を等高面とする
この等高面が、自由度dfの部分空間に作る部分面が実際の等高面となる
自由度１の統計量は、この部分空間内の等高面を定める軸への射影となっていることを利用して次のようにして、重み付け自由度１検定が行える
実行はこちらでもできる

SphericalWeightedTest<-function(t=matrix(c(10,20,30,40,50,60),nrow=2,byrow=TRUE),w=1:3){
# What marginal counts determine
# marginal counts
 m1<-apply(t,1,sum)
 m2<-apply(t,2,sum)
 ms<-sum(t)
#expected table
 et<-outer(m1,m2,FUN="*")/ms
# K-vector
 K<-m2*m1[1]*m1[2]/(ms^2)

# What define observed table
# D diff
 d<-t[1,]-et[1,]

# D'
 dpr<-d/sqrt(K)
# Pearson's chi
 C<-sum(dpr^2)
 
# What vector defines
# w to w' to w''
 wpr<-w*sqrt(K)
 shadow<-function(g,f){
  k<-sum(g*f)/sum(f^2)
  g-k*f
 }
 wpr<-shadow(g=wpr,f=sqrt(K))
 wpr<-wpr/(sqrt(sum(wpr^2)))

# What vector and observed table determine
# stat
 st<-sum(wpr*dpr)^2
 
 list(PearsonSt=C,Pearson.p.value=pchisq(C,length(m2)-1,lower.tail=FALSE),statistic=st,p.value=pchisq(abs(st),1,lower.tail=FALSE),wpr=wpr,dpr=dpr)
}



t=matrix(c(10,10,10,10,20,30),nrow=2,byrow=TRUE)
w<-c(1/2,2/3,3/4)
SphericalWeightedTest(t,w)

t=matrix(c(10,10,10,10,20,30),nrow=2,byrow=TRUE)
w<-c(1/2,2/3,3/4)
SphericalWeightedTest(t,w)