2016-01-10

matplotlibで、箱ひげ図と日本語の使い方

Python

matplotlibを使うとお手軽に箱ひげ図が書けるようです。

リファレンス：pyplot — Matplotlib 1.5.0 documentation
例：http://matplotlib.org/examples/pylab_examples/boxplot_demo.html

#coding: utf-8

import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties

fp = FontProperties(fname=r'C:\WINDOWS\Fonts\YuGothic.ttf', size=14)
#日本語表示用のフォント

data = [0.09, 0.10, 0.12, 0.04, 0.06, 0.14, 0.13, 0.28 ,0.29, 0.17]
data = [data[:], data[:]] #リストのコピーにはスライス
del data[1][7:9] #2つ目のデータから外れ値を削除

label = (u'なんとか', u'かんとか') #x軸のラベル
fig = plt.figure()
ax = fig.add_subplot(111)
ax.grid() #グリッド線
bp = ax.boxplot(data) #箱ひげ図を描画
ax.set_xticklabels(label, fontproperties=fp)
plt.ylim(0,0.3)
plt.show()

f:id:shinep:20160110052900p:plain

私が使いそうだと思ったboxplotのオプションを書いておきます。

sym:外れ値のマーカー

　'gD'(ひし形)、'rs'(正方形)、'o'(丸)

vert:バーの縦横

　True(縦)、False(横)

whis:ひげの長さ

　デフォルトは1.5、'range'で外れ値がなくなります。

matplotlibで日本語を使うには、fontproperties（凡例はprop）で和文フォント指定する必要があるようです。

Python - matplotlibで日本語 - Qiita

from matplotlib.font_manager import FontProperties
fp = FontProperties(fname=r'C:\WINDOWS\Fonts\YuGothic.ttf', size=14)
#使いたいフォントファイルを指定

が必要になります。

文字列の前にuをつけてUnicodeオブジェクトにしないと、エラーが出るみたいですね。

2016-01-05

多重比較法（Williams）

統計

Tukey-Kramer・Bonferroniについて書きましたが、他の方法が解説されているページを見つけたので、その中でも自分が使いそうなWilliamsを書いておきます。

http://www012.upp.so-net.ne.jp/doi/biostat/biostat.html
下の方の（多重比較）のところです。

Williamsは、群間に順位のあるときに採用される検定です。
全ての群をまとめた仮説を立て、それから徐々に群を狭めていく方法をとります。

群が4つある場合
$H_{1234},\ H_{123},\ H_{12}$
と仮説を立て、左から検定していきます。

仮説が棄却されたときは、さらに右の仮説を検定します。
仮説が棄却されなかったときは、その仮説より右の仮説は保留されます。

群のプール

$\mu_1\leq\mu_2\leq\mu_3\leq\mu_4$
を示したいときに、

$\bar{x}_1\leq\bar{x}_2\leq\bar{x}_3>\bar{x}_4$
のように標本平均の不等号が逆になっている箇所がある場合は、群をプールする必要があります。

やり方は簡単で、群1と群2の平均を出すだけです。
$\bar{x}_1\leq\bar{x}_2\leq\bar{x}_{34}$

それでも不等号が逆になっている場合、
$\bar{x}_1\leq\bar{x}_2>\bar{x}_{34}$

不等号の順番が全て同じになるまで、群をプールしていきます。
$\bar{x}_1\leq\bar{x}_{234}$

統計量y

標本が4群あって、群のプールが必要ないときを例にして計算方法を説明していきます。

各群のデータの合計を $T_i$ として、統計量yを計算します。
$\displaystyle y_{24}=\frac{T_2+T_3+T_4}{n_2+n_3+n_4}\\ \displaystyle y_{34}=\frac{T_3+T_4}{n_3+n_4}\\ \displaystyle y_{44}=\frac{T_4}{n_4}$

$\mu_1\leq\mu_2\leq\mu_3\leq\mu_4$ を示したいときは、計算したyの中での最大値Mを、
$\mu_1\geq\mu_2\geq\mu_3\geq\mu_4$ を示したいときは、計算したyの中での最小値mを導きます。

検定

検定は、t検定で行います。

$\mu_1\leq\mu_2\leq\mu_3\leq\mu_4$ を示したいときは、
$\displaystyle \mathrm{t}値 = \frac{M-\overline{x}_{1}}{\sqrt{{\hat\sigma}^2(\frac{1}{n_1}+\frac{1}{n_4})}}$

$\mu_1\geq\mu_2\geq\mu_3\geq\mu_4$ を示したいときは、
$\displaystyle \mathrm{t}値 = \frac{\overline{x}_1-m}{\sqrt{{\hat\sigma}^2(\frac{1}{n_1}+\frac{1}{n_4})}}$

$\hat{\sigma}^2$ は、誤差分散(郡内分散)です。
群数をaとすると、

自由度 = (標本サイズ-1) × a

$\displaystyle \hat{\sigma}^2 = \left\{{\sum_{i=1}^{a} (n_i-1)V_i}\right\}\div{自由度}$

小さい側の群は、常にプールしていない平均を使用します（ $\overline{x}_{12}$ とプールしていても、 $\overline{x}_{1}$ とする）。

t値を求めたら、TDIST関数でp値を求め、棄却されるか確かめます。
棄却されたならば、右の検定を保留されるまでやっていきます。

$H_{1234},\ H_{123}$ が棄却され、 $H_{12}$ が棄却されなければ、
$\bar{x}_1\leq\bar{x}_3,\ \bar{x}_4$
と言うことができます。

2016/1/20：方法の説明が全く合ってなかったので、下記の本を参考にして、手直ししました。

統計的多重比較法の基礎

作者: 永田靖,吉田道弘
出版社/メーカー: サイエンティスト社
発売日: 1997/12
メディア: 単行本
購入: 3人クリック: 31回
この商品を含むブログ (10件) を見る

2016-01-04

分散分析と多重比較法（Tukey-Kramer・Bonferroni）後編

統計

さて、昨日の記事の続きを書いていきます。
後編では、多重比較法を扱っていきます。
分散分析を扱った前編はこちら。
turtlewalk.hatenablog.jp

多重比較法

分散分析では、どれか1つ以上の群間に差があることしか分からず、どの群間にあるのかは分かりませんでした。
それを調べるには、多重比較が必要となります。

単純に群同士を検定していくと、有意水準が乗算されていき、第1種の過誤を犯しやすくなってしまいます。

比較する群の組み合わせ数	1	2	3	4
第1種の過誤を犯さない確率	0.95	0.90	0.86	0.81

（有意水準α=0.05のとき）

説明に昨日と同じ標本集団を使っていきます。

	群1	群2	群3	群4
標本1	10	15	24	28
標本2	11	12	26	30
標本3	8	13	23	31

有意水準は0.05とします。

t検定

まず試しに、それぞれをt検定で単純に検定していきます。
TTEST関数を使えば、簡単にp値を算出できます。

TTEST(配列1,配列2,尾部,検定の種類)

配列1と配列2にそれぞれ群の標本を入れます。
尾部（1:片側検定, 2:両側検定）
検定の種類(1:対応のある標本, 2:分散の等しい標本, 3:分散の等しくない標本)
にしたがって、該当する数値を入れます。（今回は尾部が2、検定の種類も2）

p値	群2	群3	群4
群1	0.042	<0.001	<0.001
群2		<0.001	<0.001
群3			0.013

全ての群間で有意だと検定されました。

t検定で分析ツールを使う場合は、
http://www1.tcue.ac.jp/home1/abek/htdocs/stat/Excel/t-test/t-test.html
が参考になります。

t検定での結果は、先に書いたように「甘い」検定となっていて、適当だとは言えません。
そこで多重比較法を使って、第1種の過誤の起こす確率を減らしていきます。
方法は多数ありますが、分析ツールを使わずに計算できそうだったのが、Tukey-Kramer・Bonferroniの2種類だったので、その2つを書いていきます。

2つの方法

どちらの方法にもメリット・デメリットがあります。

Tukey-Kramer

メリット

比較する標本サイズが異なっていてもよい

使えるソフトウェアが多い？

デメリット

データが正規分布に従っていなければならない

群間が等分散である必要がある

Bonferroni

メリット

ノンパラメトリックでも使える

操作が簡便

デメリット

比較する標本サイズが同じでなければならない

5群以上では検出力が落ちる

Tukey-Kramer

まずはTukey-Kramer法からやってみます。

ステューデント化された範囲のq値というものを計算し、それをq境界値と比べて検定します。
（関数での、p値の算出方法は見つけられませんでした）

$\displaystyle \mathrm{q}値=\frac{\overline{x}_1-\overline{x}_2}{\sqrt{{\hat{\sigma}_e}^2(\frac{1}{n_1}+\frac{1}{n_2})}}$

$\overline{x}$ が比較する群の平均、 $\hat{\sigma_e}$ が群内変動（誤差）の分散(前編より、2.33)、 $n$ が標本サイズ(すべて3)となります。

それぞれの群間でq値を計算すると、

q値	群2	群3	群4
群1	2.94	11.8	16.0
群2		8.82	13.1
群3			4.28

q境界値は分布表から求めます。
Deus ex machinaな日々: スチューデント化された範囲の表の補間

vは群内変動の自由度（自由度=(標本サイズ-1) × 群数）なので、8になります。
よって4.53、これを $\sqrt{2}$ で割ると、Tukey-Kramer法でのq境界値になります。

q境界値

3.20

q境界値よりもq値が大きければ、有意だと検定されます。
Tukey-Kramerでは、単純なt検定と異なり、群1-群2間では帰無仮説が有意ではありませんでした。

Bonferroni

次に、Bonferroni法です。

Bonferroni法は、t検定での結果を利用して、そのときの有意水準を検定する組み合わせの数で割った数にするだけです。
また代わりに、p値を検定する組み合わせの数で掛けることもできます。
今回は、後者の方法でやっていきます。

p値	群2	群3	群4
群1	0.254	0.002	<0.001
群2		0.005	0.001
群3			0.077

よって、群1-群2間・群3-群4間では有意と判定されませんでした。

Bonferroni法は、第2種の過誤が大きくなりやすい方法であり、その点を解消するための方法も開発されています。
http://www.med.osaka-u.ac.jp/pub/kid/clinicaljournalclub1.html

このように多重比較法では、検定が厳しくなっています。

この他にも、多重比較法はたくさんあり、用途によってどれを用いるか検討していく必要があります。

亀の歩み

備忘録として

matplotlibで、箱ひげ図と日本語の使い方

多重比較法（Williams）

群のプール

統計量y

検定

分散分析と多重比較法（Tukey-Kramer・Bonferroni）後編

多重比較法

t検定

2つの方法

Tukey-Kramer

Bonferroni

Tukey-Kramer

Bonferroni