最近，疫学研究の勉強をしています．疫学研究はR言語で解説されていることが多い印象ですが，Python好きのみーとしては，Pythonで統計処理したいので，Pythonで疫学研究を行うことにしました．
代表的な統計計算をPythonで行えるかをUnittestにしました．問題例と正解はhttp://cse.naro.affrc.go.jp/takezawa/r-tips/r.htmlを参照しています．
検証したバージョン
- Python 3.7.3
- numpy 1.16.2
- scipy 1.3.0
Unittest
最初に統計処理するためのデータを用意します．
import scipy.stats import numpy import unittest class TestStats(unittest.TestCase): def setUp(self): self.a = [0.7, -1.6, -0.2, -1.2, -0.1, 3.4, 3.7, 0.8, 0.0, 2.0] self.b = [1.9, 0.8, 1.1, 0.1, - 0.1, 4.4, 5.5, 1.6, 4.6, 3.4] self.c = [ [14, 8], [4, 17] ] self.d = [70, 72, 62, 64, 71, 76, 60, 65, 74, 72] self.e = [70, 74, 65, 68, 72, 74, 61, 66, 76, 75] # modified from https://www.stat.cmu.edu/~larry/all-of-statistics/=data/faithful.dat self.faithful = [3.6, 3.333, 4.533, 4.7, 3.6, 4.35, 3.917, 4.2, 4.7, 4.8, 4.25, 3.45, 3.067, 4.533, 3.6, 4.083, 3.85, 4.433, 4.3, 4.467, 3.367, 4.033, 3.833, 4.833, 4.783, 4.35, 4.567, 4.533, 3.317, 3.833, 4.633, 4.8, 4.716, 4.833, 4.883, 3.717, 4.567, 4.317, 4.5, 4.8, 4.4, 4.167, 4.7, 4.7, 4.033, 4.5, 4, 5.067, 4.567, 3.883, 3.6, 4.133, 4.333, 4.1, 4.067, 4.933, 3.95, 4.517, 4, 4.333, 4.817, 4.3, 4.667, 3.75, 4.9, 4.367, 4.5, 4.05, 4.7, 4.85, 3.683, 4.733, 4.9, 4.417, 4.633, 4.6, 4.417, 4.067, 4.25, 4.6, 3.767, 4.5, 4.65, 4.167, 4.333, 4.383, 4.933, 3.733, 4.233, 4.533, 4.817, 4.333, 4.633, 5.1, 5.033, 4, 4.6, 3.567, 4, 4.5, 4.083, 3.967, 4.15, 3.833, 3.5, 4.583, 5, 4.617, 4.583, 3.333, 4.167, 4.333, 4.5, 4, 4.167, 4.583, 4.25, 3.767, 4.433, 4.083, 4.417, 4.8, 4.8, 4.1, 3.966, 4.233, 3.5, 4.366, 4.667, 4.35, 4.133, 4.6, 4.367, 3.85, 4.5, 4.7, 3.833, 3.417, 4.233, 4.8, 4.15, 4.267, 4.483, 4, 4.117, 4.083, 4.267, 3.917, 4.55, 4.083, 4.183, 4.45, 4.283, 3.95, 4.15, 4.933, 4.583, 3.833, 4.367, 4.35, 4.45, 3.567, 4.5, 4.15, 3.817, 3.917, 4.45, 4.283, 4.767, 4.533, 4.25, 4.75, 4.117, 4.417, 4.467]
59. 基本統計量の算出
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/59.html
4分位偏差 IQR
def test_quantile(self): # クォンタイル # https://docs.scipy.org/doc/numpy-1.15.1/reference/generated/numpy.quantile.html tests = [ (0, -1.600), (0.25, -0.175), (0.5, 0.350), (0.75, 1.700), (1, 3.700) ] for q, s in tests: with self.subTest(q=q): r = numpy.quantile(self.a, q) self.assertAlmostEquals(r, s)
4分位偏差 IQR
def test_iqr(self): # 4分位偏差 IQR() ： 第3四分位から第1四分位を引いた値 # https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.stats.iqr.html r = scipy.stats.iqr(self.a) self.assertAlmostEquals(r, 1.700-(-0.175))
最小値
def test_min(self): r = min(self.b) self.assertAlmostEquals(r, -0.100)
中央値
def test_median(self): # https://docs.scipy.org/doc/numpy-1.15.1/reference/generated/numpy.median.html r = numpy.median(self.b) self.assertAlmostEquals(r, 1.750)
平均値
def test_mean(self): # https://docs.scipy.org/doc/numpy-1.15.1/reference/generated/numpy.mean.html r = numpy.mean(self.b) self.assertAlmostEquals(r, 2.330)
最大値
def test_max(self): r = max(self.b) self.assertAlmostEquals(r, 5.500)
不偏分散
def test_var(self): # 不偏分散 # https://docs.scipy.org/doc/numpy-1.15.1/reference/generated/numpy.var.html r = numpy.var(self.a, ddof=1) self.assertAlmostEquals(r, 3.200556, places=4)
標本分散
def test_variance(self): # 標本分散 # https://docs.scipy.org/doc/numpy-1.15.1/reference/generated/numpy.var.html r = numpy.var(self.a) self.assertAlmostEquals(r, 2.8805)
不偏標準偏差
def test_sd(self): # 不偏標準偏差 # https://docs.scipy.org/doc/numpy-1.15.0/reference/generated/numpy.std.html r = numpy.std(self.a, ddof=1) self.assertAlmostEquals(r, 1.789010, places=4)
62. 検定事始
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/62.html
一標本 t 検定 (両側検定)
def test_t_test_one_sample_two_sided(self): # 一標本 t 検定 (両側検定) # https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.ttest_1samp.html t, pvalue = scipy.stats.ttest_1samp(self.a, 3) self.assertAlmostEquals(t, -3.9771, places=4) self.assertAlmostEquals(pvalue, 0.00322, places=4)
二標本 t 検定 (対応なし)
def test_t_test_two_samples(self): # 二標本 t 検定 (対応なし) # https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.ttest_ind.html t, pvalue = scipy.stats.ttest_ind(self.a, self.b, equal_var=True) self.assertAlmostEquals(t, -1.8608, places=4) self.assertAlmostEquals(pvalue, 0.07919, places=4)
ウェルチの検定
def test_t_test_two_samples_not_equal_var(self): # 二標本 t 検定 (対応なし) (等分散でない) = Welch Two Sample t-test (ウェルチの検定) # https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.ttest_ind.html t, pvalue = scipy.stats.ttest_ind(self.a, self.b, equal_var=False) self.assertAlmostEquals(t, -1.8608, places=4) self.assertAlmostEquals(pvalue, 0.0794, places=4)
63. 正規性の検定
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/63.html
Shapiro-Wilk 検定
def test_shapiro_test(self): # 正規性の検定 Shapiro-Wilk 検定 # https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.shapiro.html w, pvalue = scipy.stats.shapiro(self.faithful) self.assertAlmostEquals(w, 0.9793, places=4) self.assertAlmostEquals(pvalue, 0.01052, places=4)
64. 一標本検定
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/64.html
一標本 t 検定 (両側検定)
def test_t_test_one_sample_two_sided_with_ci(self): # 一標本 t 検定 (両側検定) # https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.ttest_1samp.html t, pvalue = scipy.stats.ttest_1samp(self.a, 0) self.assertAlmostEquals(t, 1.3257, places=4) self.assertAlmostEquals(pvalue, 0.2176, places=4) ci1, ci2 = scipy.stats.t.interval( 0.95, len(self.a)-1, loc=numpy.mean(self.a), scale=scipy.stats.sem(self.a) ) self.assertAlmostEquals(ci1, -0.5297804, places=4) self.assertAlmostEquals(ci2, 2.0297804, places=4)
二標本 t 検定 (分散等しい)
def test_t_test_two_samples_var_equal_with_ci(self): # 二標本 t 検定 (分散等しい) t, pvalue = scipy.stats.ttest_ind(self.a, self.b, equal_var=True) self.assertAlmostEquals(t, -1.8608, places=4) self.assertAlmostEquals(pvalue, 0.07919, places=4) # CI ?
65. 二標本検定
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/65.html
ウェルチの検定
def test_t_test_two_samples_not_equal_var_with_ci(self): # 二標本 t 検定 (対応なし) (等分散でない) = Welch Two Sample t-test (ウェルチの検定) # https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.ttest_ind.html t, pvalue = scipy.stats.ttest_ind(self.a, self.b, equal_var=False) self.assertAlmostEquals(t, -1.8608, places=4) self.assertAlmostEquals(pvalue, 0.0794, places=4) # CI ?
マン・ホイットニーの U 検定
def test_wilcoxon_rank_sum_test(self): # ウィルコクソンの順位和検定 = マン・ホイットニーの U 検定 # https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ranksums.html w, pvalue = scipy.stats.mannwhitneyu( self.a, self.b, alternative="two-sided") self.assertAlmostEquals(w, 25.5, places=4) self.assertAlmostEquals(pvalue, 0.06933, places=4)
66. カテゴリデータの検定
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/66.html
フィッシャーの直接確率検定
def test_fisher_exact_test(self): # フィッシャーの直接確率検定 # https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.fisher_exact.html odds, pvalue = scipy.stats.fisher_exact( self.c, alternative="two-sided") self.assertAlmostEquals(pvalue, 0.005089, places=4) self.assertNotAlmostEquals(odds, 7.051895, places=4) # The calculated odds ratio is different from the one R uses. # This scipy implementation returns the (more common) “unconditional Maximum Likelihood Estimate”, # while R uses the “conditional Maximum Likelihood Estimate”.
R言語とは異なるオッズ比になるらしい．
67. 相関係数と無相関検定
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/67.html
相関係数
def test_cor(self): # 相関係数 ピアソン 検定 r, pvalue = scipy.stats.pearsonr(self.d, self.e) self.assertAlmostEquals(pvalue, 2.656e-05, places=6) self.assertAlmostEquals(r, 0.9496011, places=6)
とりあえず，みーが必要そうな統計処理は網羅したので，ここまでとします．R言語と違って，scipyは信頼区間を自動で計算してくれないので，不完全なところもあります．
