質的変数の連関分析 2017-12-28 †

↑

分割表 †

分割表（contingency table）とは、2つ以上の質的変数（特に名義尺度〜量的大小が無い場合も含めて）の間の関係を記録し分析するための表。

例：　Rで得られるTitanic遭難データから、Class=2nd、Age=Adultのデータを抽出した。

生存死亡
男性154.014.0
女性13.080.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
#４次元の配列で、内容は
# 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に移したいが４次元なのでいきなりは無理(DFは２次元)。
#適当な２つの次元で切り出す。
#たとえば、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)

このような（カテゴリカル）データに対して、相関（連関と呼ぶ）つまり男女と生存・死亡に 関連があるかを調べたい。

↑

連関分析の方法・係数 †

ab
cd
↑

ユールのQ †

$Q=\frac{ad-bc}{ad+bc}$

↑

ファイ係数（$\phi$係数） †

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}$

↑

ピアソンの$\chi^2$ と クラメールの連関係数$V$ †

観測値を$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の分割表で示される２つの変数が互いに独立かどうかを検定する。　

12...c計
1n11n12...n1cn1*
2n21n22...n2cn2*
...
rnr1nr2...nrcnr*
計n*1n*2...n*cN

この時の期待度数Eは次のように計算する。独立性の検定 期待度数の最尤推定量の導出
　　$E_{ij} = \frac{n_{i*}　\cdot 　n_{*j}}{N}$
自由度は、 　　$自由度 = (r-1)(c-1)$
これを用いて、$\chi^2$を求める。

2x2の表の場合、Titanicの例では、

生存死亡計
男性15414168
女性138093
計16794261

であり、期待度数はそれぞれ

生存死亡計
男性168x167/281=107.5168x94/281=56.2168
女性93x167/281=55.393x94/281=31.193
計16794261

検定統計量は

生存死亡
男性(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

検定統計量はこの４者の和なので、161.1となる。

自由度=(2-)(2-1)=1のカイ二乗分布に従うので、有意水準$\alpha=0.05$の値を見ると3.84であり、帰無仮説は棄却され、独立ではないと結論できる。

なお、2x2の場合、

B1B2
A1ab
A2cd

式を整理すると、 $\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$の値によって著しく異なる。（図は相関分析にある。）

一般に相関の対象となる２つの変数が２変量正規分布（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$検定する。

相関係数の検定に関わるとき、ピアソンの積率相関係数の検定は対象となる２変数が（２変量）正規分布を仮定するが（「なぜか」への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$程度に少ないということになる。

↑

無相関検定 †

また、

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2020-01-17 (金) 08:00:42 (1048d)