pairwise_tukeyhsdのp値の算出方法

Question

Pythonで統計学を学ぶ(6)を参考にさせていただきました。

from statsmodels.stats.multicomp import pairwise_tukeyhsd
import numpy as np

def tukey_hsd( lst, ind, n ):
    data_arr = np.hstack( lst )
    ind_arr = np.repeat(ind, n)
    print(pairwise_tukeyhsd(data_arr,ind_arr))

上記のようにTukeyの多重比較比較検定を行うような関数を定義します。
以下、A,B,C,Dを以下のように定義し多重比較検定を実行します。

A = np.array([15,9,18,14,18])
B = np.array([13,8,8,12,7])
C = np.array([10,6,11,7,12])
D = np.array([10,7,3,5,7])

tukey_hsd( (A,B,C,D), list('ABCD') , 5) 

出力

Multiple Comparison of Means - Tukey HSD,FWER=0.05
==============================================
group1 group2 meandiff  lower    upper  reject
----------------------------------------------
  A      B      -5.2   -10.5145  0.1145 False 
  A      C      -5.6   -10.9145 -0.2855  True 
  A      D      -8.4   -13.7145 -3.0855  True 
  B      C      -0.4   -5.7145   4.9145 False 
  B      D      -3.2   -8.5145   2.1145 False 
  C      D      -2.8   -8.1145   2.5145 False 
----------------------------------------------

と出力されます。
ここでrejectをFWER=0.05の値を用いて判定していることが分かります。
つまり、group1とgroup2を比較するとき必ずp値を算出していることが分かります。
しかし、私はその算出されているであろうp値を調べる方法が分かりません。

statsmodels.stats.multicomp.pairwise_tukeyhsd

を見る限りはstatsmodels.sandbox.stats.multicomp.TukeyHSDResultsのインスタンスが出力(return)されていて
statsmodels.sandbox.stats.multicomp.TukeyHSDResults
を見る限るにはAttributesに(pvaluesadjusted p-values from the HSD test)と書かれているので、
そのTukeyHSDResultsインスタンスのAttributesの中身を見る方法があればこの問題は解決すると思うのですが。
どのようにすればいいのかここで詰まってしまってしまっています。
分かる方がいれば教えていただけると大変恐縮です。

Answer 1

内容を全面改訂しました。

print(pairwise_tukeyhsd(data_arr,ind_arr))
を実行するとpairwise_tukeyhsd(data_arr,ind_arr)はTukeyHSDResultsインスタンスを返してきますがメソッド__str__()で文字列に変換されてしまいます。

一旦、変数で受けて、vars()で表示すればよいです。

def tukey_hsd( lst, ind, n ):
    data_arr = np.hstack( lst )
    ind_arr = np.repeat(ind, n)
    res = pairwise_tukeyhsd(data_arr, ind_arr)
    print(vars(res))


A = np.array([15,9,18,14,18])
B = np.array([13,8,8,12,7])
C = np.array([10,6,11,7,12])
D = np.array([10,7,3,5,7])

tukey_hsd( (A,B,C,D), list('ABCD') , 5) 

{'_multicomp': <statsmodels.sandbox.stats.multicomp.MultiComparison object at 0x0000015F51131C88>, '_results_table': <class 'statsmodels.iolib.table.SimpleTable'>, 'q_crit': 4.046412438282385, 'reject': array([False,  True,  True, False, False, False]), 'meandiffs': array([-5.2, -5.6, -8.4, -0.4, -3.2, -2.8]), 'std_pairs': array([1.31339255, 1.31339255, 1.31339255, 1.31339255, 1.31339255,
       1.31339255]), 'confint': array([[-10.51452797,   0.11452797],
       [-10.91452797,  -0.28547203],
       [-13.71452797,  -3.08547203],
       [ -5.71452797,   4.91452797],
       [ -8.51452797,   2.11452797],
       [ -8.11452797,   2.51452797]]), 'df_total': 16, 'reject2': array([False,  True,  True, False, False, False]), 'variance': 8.625000000000002, 'pvalues': array([0.0562591 , 0.03714849, 0.00177409, 0.9       , 0.34502168,
       0.45735327]), 'data': array([15,  9, 18, 14, 18, 13,  8,  8, 12,  7, 10,  6, 11,  7, 12, 10,  7,
        3,  5,  7]), 'groups': array(['A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'C', 'C', 'C',
       'C', 'C', 'D', 'D', 'D', 'D', 'D'], dtype='<U1'), 'groupsunique': array(['A', 'B', 'C', 'D'], dtype='<U1')}

'pvalues': array([0.0562591 , 0.03714849, 0.00177409, 0.9 , 0.34502168, 0.45735327])が所望のデータだと思います。

Comments

[yu__]

vars()の引数にインスタンスを入れるとそのクラスの__init__()の情報を確認することができるということですね。初めて知りました。ものすごく勉強になりました。ありがとうございました。
ベストアンサーにさせていただきました。

しかしながら、__init__()の
self.pvalues
Taken out of _multicomp for ease of access for unknowledgeable users
と書かれていて
pvaluesだけ結果が算出されません。
nomukenさんが表示しているresultを出したとき何か特別設定はしましたか？
回答していただけると幸いです。

[nomuken]

> nomukenさんが表示しているresultを出したとき何か特別設定はしましたか？
特にしてません。`print(res.pvalues)`と書いてpvaluesだけ出力することもできます。

> Taken out of _multicomp for ease of access for unknowledgeable users
> pvaluesだけ結果が算出されません。
この文章はそのコメントの後ろ３行の命令にかかっている文章だと思います。
```
self.data = self._multicomp.data
self.groups = self._multicomp.groups
self.groupsunique = self._multicomp.groupsunique
```
self.dataとself.groupsとself.groupsuniqueはself._multicompに入っているがデータにアクセスしやすくするために行っているという意味だと思います。

pvaluesの算出は
https://github.com/statsmodels/statsmodels/blob/bc5680db6265d275d89505815a5cec9e9f632239/statsmodels/sandbox/stats/multicomp.py#L1298
の計算でやっていると思います。