みーのぺーじ

みーが趣味でやっているPCやソフトウェアについて.Python, Javascript, Processing, Unityなど.

疫学研究の統計に関するunittest

最近,疫学研究の勉強をしています.疫学研究は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は信頼区間を自動で計算してくれないので,不完全なところもあります.