ただのメモ

他の人に見せても良い方のメモ

Cochran-Mantel-Haenszel test

Cochran Mantel Haenszel検定(以下、CMH検定)について調べて、お絵かきしてみる。

始めましょうか。

  • まずこちら。
  •  CMH検定は交絡因子を考えた上で、曝露と結果の関係を推定する手法。
    • 交絡因子で複数の群に分けて、それぞれに2×2の表を作り、その平均のリスク比を考える。

お絵かきしてみる。

  • χ二乗分布に漸近的に近づくことを言うために、普通は分布を表すグラフを用いるのだろう。
  • しかし、今回は違う表現の仕方が無いかと考えた。(頭の体操)
  • 取り出したい特徴
    • グラフ
      • 高さで分布におけるそこでの多さを示す。
    • 散布図
      • 密度で分布におけるそこでの多さを示す。
      • →散布図が使えそう。
  • 散布図を使えると思ったが、1次元圧縮されすぎてわかんないので、2次元の表にする。
  • コードはこちら。
#%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import math
import random

#initial settings
N=3000;dt=0.01;T=N*dt
#a=0.1;b=0.1;c=20;d=0.3;K0=5;K=20;


#function settings
def T(a,b,c,d):
    return a+b+c+d
def N1(a,b):
    return a+b
def N2(c,d):
    return c+d
def M1(a,c):
    return a+c
def M2(b,d):
    return b+d
def partialCMHm(a,b,c,d):
    return a-N1(a,b)*M1(a,c)/T(a,b,c,d)
def partialCMHc(a,b,c,d):
    return N1(a,b)*N2(c,d)*M1(a,c)*M2(b,d)/(T(a,b,c,d)**3-T(a,b,c,d)**2)

def CMH(data):
    CMHm=0;CMHc=0;K=len(data)
    for i in range(K):
        CMHm+=partialCMHm(data[i][0],data[i][1],data[i][2],data[i][3])
        CMHc+=partialCMHc(data[i][0],data[i][1],data[i][2],data[i][3])
    CMHm=CMHm**2
    return CMHm/CMHc

def main():
    s=10
    data = [[random.randint(1,100) for i in range(s)]]
    plt.scatter(CMH(data),random.random(),s=10)

for i in range(N):
    main()

plt.xlabel("x")
plt.ylabel("y")
plt.title("Example 408")

plt.show()

 

f:id:medical-science:20220408095221p:plain

  • CMH検定がχ二乗検定の塊みたいなことをしていて、それ故にχ二乗分布との関わりが出てくる、ということか。
  • 交絡因子は名義尺度でも順序尺度でも良いみたい。
    • BMI25以上か以下か、は、肥満か非肥満か、という言い換えが出来る。

バイバイ!