![]() |
ノート/Python/統計/質的変数の連関解析https://pepper.is.sci.toho-u.ac.jp:443/pepper/index.php?%A5%CE%A1%BC%A5%C8%2FPython%2F%C5%FD%B7%D7%2F%BC%C1%C5%AA%CA%D1%BF%F4%A4%CE%CF%A2%B4%D8%B2%F2%C0%CF |
![]() |
分割表(contingency table)とは、2つ以上の質的変数(特に名義尺度〜量的大小が無い場合も含めて)の間の関係を記録し分析するための表。
例: Rで得られるTitanic遭難データから、Class=2nd、Age=Adultのデータを抽出した。
生存 | 死亡 | |
男性 | 154.0 | 14.0 |
女性 | 13.0 | 80.0 |
この例は2x2の表であるが、一般的にmxnの表が可能。
抽出のプログラムは
import pandas as pd from rpy2.robjects import pandas2ri pandas2ri.activate() from rpy2.robjects import r # RのTitanicのデータを使いたい。 # https://stat.ethz.ch/R-manual/R-devel/library/datasets/html/Titanic.html #4次元の配列で、内容は # No Name Levels # 1 Class 1st, 2nd, 3rd, Crew # 2 Sex Male, Female # 3 Age Child, Adult # 4 Survived No, Yes #データそのものは、 # Class Sex Age Survived Freq #1 1st Male Child No 0 #2 2nd Male Child No 0 #3 3rd Male Child No 35 #4 Crew Male Child No 0 #5 1st Female Child No 0 #6 2nd Female Child No 0 #7 3rd Female Child No 17 #8 Crew Female Child No 0 #9 1st Male Adult No 118 #10 2nd Male Adult No 154 #11 3rd Male Adult No 387 #12 Crew Male Adult No 670 #13 1st Female Adult No 4 #14 2nd Female Adult No 13 #15 3rd Female Adult No 89 #16 Crew Female Adult No 3 #17 1st Male Child Yes 5 #18 2nd Male Child Yes 11 #19 3rd Male Child Yes 13 #20 Crew Male Child Yes 0 #21 1st Female Child Yes 1 #22 2nd Female Child Yes 13 #23 3rd Female Child Yes 14 #24 Crew Female Child Yes 0 #25 1st Male Adult Yes 57 #26 2nd Male Adult Yes 14 #27 3rd Male Adult Yes 75 #28 Crew Male Adult Yes 192 #29 1st Female Adult Yes 140 #30 2nd Female Adult Yes 80 #32 Crew Female Adult Yes 20 #データはnumpyのndarrayで返ってくる # titanic = r["Titanic"] #pandasのDataFrameに移したいが4次元なのでいきなりは無理(DFは2次元)。 #適当な2つの次元で切り出す。 #たとえば、Class=2nd、Age=Adultに固定して、SexとSurvived/Noを切り出す。 Sex_2nd_Adult=pd.DataFrame(r["Titanic"][1,:,1,:]) print(Sex_2nd_Adult) #結果は、 # 0 1 # 0 154.0 14.0 # 1 13.0 80.0 #pandasで列名、行名を変更するには # df.index = [1, 2, 3] # df.columns = ['a', 'b', 'c'] Sex_2nd_Adult.index = ['男性', '女性'] Sex_2nd_Adult.columns = ['生存', '死亡'] print(Sex_2nd_Adult)
このような(カテゴリカル)データに対して、相関(連関と呼ぶ)つまり男女と生存・死亡に 関連があるかを調べたい。
a | b | |
c | d |
$Q=\frac{ad-bc}{ad+bc}$
2$\times$2の分割表の連関を示す係数。2つの導入がありそれぞれ同値。
(1) それぞれの軸の値を0/1として(連続値と同じ)ピアソンの積算相関係数を計算する方法
$Cor = \sum_i (X_i – \mu_x)(Y_i – \mu_y) / ( N\times SD_x \times SD_y )$
(2) 計算式 $P = \frac{(ad - bc)}{\sqrt{(a + b)(c + d)(a + b)(c + d)}}$
これが(1)と等しいことは、東京大学 上田博人先生の言語データ分析の授業資料『言語研究のための数値データ分析法』の第6章 関係[PDF]の17ページにある式変形で証明される。
オッズ比 = a/bのオッズ $\frac{a}{b}$ と c/dのオッズ $\frac{c}{d}$ の比率
$\frac{a/b}{c/d}$ = $\frac{ad}{bc}$
観測値を$O$、理論値(期待値)を$E$とするとき、カイ二乗は
$\chi^2 = \sum \frac{(O-E)^2}{E}$
である。これを分割表で適合度検定・独立性検定に用いる。
適合度検定は、観測された度数分布が理論分布と同じかどうかを検定する。例題として100人の標本中の男女の人数の観測度数を、理論値(男女が等しい)に対して比較する。男女が等しいとする帰無仮説を有意水準$\alpha$に対して検定する。
観測値が、男性40、女性60であったとすると、
$\chi^2=\sum \frac{(v_i-np_i)^2}{np_i} = \frac{(40-50)^2}{50}+\frac{(60-50)^2}{50} = 4$
この場合の自由度は1であるので、自由度1の$\chi^2$分布表から$\alpha =0.05$のパーセント点を求めると3.84であり、$\chi^2$値が越えてしまうので帰無仮説「男女は等しい」は棄却され、等しいと言えないという結論になる。またもし、男性45、女性55であれば、$\chi^2=1$になるので、仮説は棄却されない(等しいか等しくないか言えない)。
なお、$\chi^2$分布で近似できるのは、期待度数が$np_i>=10$をすべての属性に対して満たすときとされている。
独立性の検定は、2x2の分割表で示される2つの変数が互いに独立かどうかを検定する。
1 | 2 | ... | c | 計 | |
1 | n11 | n12 | ... | n1c | n1* |
2 | n21 | n22 | ... | n2c | n2* |
... | |||||
r | nr1 | nr2 | ... | nrc | nr* |
計 | n*1 | n*2 | ... | n*c | N |
この時の期待度数Eは次のように計算する。独立性の検定 期待度数の最尤推定量の導出
$E_{ij} = \frac{n_{i*} \cdot n_{*j}}{N}$
自由度は、
$自由度 = (r-1)(c-1)$
これを用いて、$\chi^2$を求める。
2x2の表の場合、Titanicの例では、
生存 | 死亡 | 計 | |
男性 | 154 | 14 | 168 |
女性 | 13 | 80 | 93 |
計 | 167 | 94 | 261 |
であり、期待度数はそれぞれ
生存 | 死亡 | 計 | |
男性 | 168x167/281=107.5 | 168x94/281=56.2 | 168 |
女性 | 93x167/281=55.3 | 93x94/281=31.1 | 93 |
計 | 167 | 94 | 261 |
検定統計量は
生存 | 死亡 | |
男性 | (154-107.5)^2/107.5=20.1 | (14-56.2)^2/56.2=31.7 |
女性 | (13-55.3)^2/55.3=32.4 | (80-31.1)^2/31.1=76.9 |
検定統計量はこの4者の和なので、161.1となる。
自由度=(2-)(2-1)=1のカイ二乗分布に従うので、有意水準$\alpha=0.05$の値を見ると3.84であり、帰無仮説は棄却され、独立ではないと結論できる。
なお、2x2の場合、
B1 | B2 | |
A1 | a | b |
A2 | c | d |
式を整理すると、 $\chi^2=\frac{N(ad-bc)^2}{(a+b)(c+d)(a+c)(b+d)}$ となる。
クラメールの連関係数は、$\chi^2$をピアソンのカイ二乗検定で得られる値、$N$をデータの総数、$k$を行数・列数の少ない方とすると
$V=\sqrt{\frac{\chi^2}{N(k-1)}}$
2$\times$2の分割表の場合、$k=2$となるので、
$V=\sqrt{\frac{\chi^2}{N}}$
但し、$\chi^2$検定による独立性の検定は、すべてのセルの期待値が10未満程度であると正しくない。フィッシャーの正確検定はセルの値が小さい場合にも検定できる。
なお、2$\times$2の分割表の場合に限り、クラメールの連関係数$V$とファイ係数$\phi$の絶対値は一致する。
(南風原朝和「心理統計学の基礎」(有斐閣アルマ)p188、カテゴリカルデータ分析入門 p8)
つまり、$V=\sqrt{\frac{\chi^2}{N}}=|\phi|$ である。
クラメールの連関係数 V プログラム例 2020-01-17追加
pandas - Categorical features correlation - Stack Overflow
python - Using pandas, calculate Cramér's coefficient matrix - Stack Overflow
import pandas as pd import numpy as np import scipy.stats as ss import seaborn as sns tips = sns.load_dataset("tips") print(tips) tips["total_bill_cut"] = pd.cut(tips["total_bill"], np.arange(0, 55, 5), include_lowest=True, right=False) def cramers_v(confusion_matrix): """ calculate Cramers V statistic for categorial-categorial association. uses correction from Bergsma and Wicher, Journal of the Korean Statistical Society 42 (2013): 323-328 """ chi2 = ss.chi2_contingency(confusion_matrix)[0] n = confusion_matrix.sum() phi2 = chi2 / n r, k = confusion_matrix.shape phi2corr = max(0, phi2 - ((k-1)*(r-1))/(n-1)) rcorr = r - ((r-1)**2)/(n-1) kcorr = k - ((k-1)**2)/(n-1) return np.sqrt(phi2corr / min((kcorr-1), (rcorr-1))) # day: Sun...Sat(7選択肢), time: Dinner or Lunch # confusion_matrix = pd.crosstab(tips["day"], tips["time"]).as_matrix() # FutureWarning: Method .as_matrix will be removed in a future version. Use .values instead. confusion_matrix = pd.crosstab(tips["day"], tips["time"]).values print(cramers_v(confusion_matrix)) # Out[10]: 0.93866193407222209 # total_bill_cut: # total_bill=金額, # total_bill_cut=5ずつに区間分けしたもの pd.cut(tips["total_bill"], np.arange(0, 55, 5), include_lowest=True, right=False) #onfusion_matrix = pd.crosstab(tips["total_bill_cut"], tips["time"]).as_matrix() confusion_matrix = pd.crosstab(tips["total_bill_cut"], tips["time"]).values print(cramers_v(confusion_matrix)) # Out[24]: 0.16498707494988371
上記のインターフェースを変えたバージョン
The Search for Categorical Correlation - Towards Data Science
# 上記stackoverflow版を改良したといっている import pandas as pd import numpy as np import scipy.stats as ss import seaborn as sns def cramers_v(x, y): confusion_matrix = pd.crosstab(x,y) chi2 = ss.chi2_contingency(confusion_matrix)[0] n = confusion_matrix.sum().sum() phi2 = chi2/n r,k = confusion_matrix.shape phi2corr = max(0, phi2-((k-1)*(r-1))/(n-1)) rcorr = r-((r-1)**2)/(n-1) kcorr = k-((k-1)**2)/(n-1) return np.sqrt(phi2corr/min((kcorr-1),(rcorr-1))) tips = sns.load_dataset("tips") tips["total_bill_cut"] = pd.cut(tips["total_bill"], np.arange(0, 55, 5), include_lowest=True, right=False) print(cramers_v(tips["day"], tips["time"])) # 元DFの列だけ(=Series)与えればいい
連続変数に対するピアソンの(積率)相関係数の算出法を、順序尺度の場合に適用したもの。順位を数値としてピアソンの相関係数を計算する。
$\rho = 1-\frac{6\sum_i{d_i}^2}{n(n^2-1)}$
但し$n$はデータの個数、$d_i$は順位の差。同順位は存在しないとする。
この式は、ピアソンの相関係数の定義から単純に式変形で導出される。(スピアマンの順位相関係数)
$r_{xy} = \frac{s_{xy}}{s_x\cdot s_y}$
$= \frac{(1/n)\sum (x_i-\overline{x})(y_i-\overline{y})}{\sqrt{(1/n)\sum (x_i-\overline{x})^2\times (1/n)\sum (y_i-\overline{y})^2}} $
$= \frac{\sum x_i\cdot y_i - n\overline{x}\cdot \overline{y}}{\sqrt{\sum ({x_i}^2-n\overline{x}^2) \times \sum ({y_i}^2 - n\overline{y}^2) }} $
ここで、$x$の個数=$y$の個数=$n$で同順位がないとすると、
$\sum x_i=\sum y_i = n(n+1)/2$
$\sum {x_i}^2=\sum {y_i}^2 = n(n+1)(2n+1)/6 $
$\overline{x}=(n+1)/2$、 $\overline{y} = (n+1)/2$
であるので、
$r_{xy}=\frac{\sum x_i y_i -n((n+1)^2/4)}{n(n+1)(n-1)/12}$
$=\frac{12\cdot \sum x_i y_i -3n\cdot (n^2+2n+1)}{n(n+1)(n-1)}$
$=1 - \frac{-12\cdot \sum x_i y_i +3n\cdot (n^2+2n+1) + n(n+1)(n-1)}{n(n+1)(n-1)}$
$=1 - \frac{1}{n(n+1)(n-1)} \cdot (-12\cdot \sum x_i y_i +3n\cdot (n+1)^2 + n(n+1)(n-1))$
$=1 - \frac{6}{n(n+1)(n-1)} \cdot (-2\cdot \sum x_i y_i + (1/6)n(n+1)(3(n+1)+(n-1)))$
$=1 - \frac{6}{n(n+1)(n-1)}\cdot (-2\cdot \sum x_i y_i + (1/6)n(n+1)(2n+1))$
$=1 - \frac{6}{n(n+1)(n-1)}\cdot (-2\cdot \sum x_i y_i +2\frac{n(n+1)(2n+1)}{6}))$
$=1 - \frac{6}{n(n+1)(n-1)}\cdot (-2\cdot \sum x_i y_i + \sum {x_i}^2 + \sum {y_i}^2))$
$=1 - \frac{6}{n(n+1)(n-1)}\cdot (\sum {x_i}^2 + \sum {y_i}^2) + -2\sum x_i y_i)$
$=1 - \frac{6}{n(n^2-1)}\cdot \sum(x_i-y_i)^2 $
標本から求めた相関係数$r$から、母集団の相関係数(母相関係数)$\rho$ を区間推定したい。 そのとき、母相関係数$\rho$が与えられた時の標本相関係数$r$の分布は、$\rho$の値によって著しく異なる。(図は相関分析にある。)
一般に相関の対象となる2つの変数が2変量正規分布(bivariate normal distribution)であるとき、標本相関係数$r$の分布は、
$\rho=0$のとき、統計量 $t=\frac{r\sqrt{n-2}}{\sqrt{1-r^2}}$ は自由度$n-2$の$t$分布に従う
$\rho\neq 0$のとき、統計量 $Z=\frac{1}{2}\ln \frac{1+r}{1-r}$ は近似的に $N(\frac{1}{2}\ln \frac{1+\rho}{1-\rho}, \frac{1}{n-3})$ に従う
ことが分かっている。
この関係を使って、母相関係数を区間推定できる。
点推定の場合: $\rho$の推定値 $\hat\rho = r$
区間推定の場合: 正規分布で推定
また、相関の検定は、通常は無相関の検定(母集団に相関が無い$\rho=0$という帰無仮説を検定)として行われる。従って、帰無仮説$\rho=0$が成り立つとして、統計量$t$
$t=\frac{r\sqrt{n-2}}{\sqrt{1-r^2}}$
が自由度$n-2$の$t$分布に従うことを用いて、$t$検定する。
相関係数の検定に関わるとき、ピアソンの積率相関係数の検定は対象となる2変数が(2変量)正規分布を仮定するが(「なぜか」へのMatt Williamsの返答)、スピアマンの順序相関係数は仮定しない。
フィッシャーの正確検定は、$\chi^2$を用いた独立性の検定が各要素の値が10未満程度の場合には向いていないが、この場合でも適用できる。考え方は、縦計a+c, b+dと横計a+b, c+dを固定したときに上記の分割表を得られる確率$p$は、すべての組合せの計算から、
$p = \frac{(a+b)!(c+d)!(a+c)!(b+d)!}{n!a!b!c!s!}$
となる。この値はいわゆる$p$値であり、$p$が小さいとこのような組み合わせの表が生じる確率は$p$程度に少ないということになる。
また、