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()

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

バイバイ！