python 分散分析を行う


統計検定の問題で、一元配置の分散分析に対する問題がありました。
python で分散分析を行なったことがなく、ライブラリと実装方法について調べてみましたので、結果を記載します。


分散分析について

史上最強 図解 これならわかる!統計学 | 涌井 良幸, 涌井 貞美 |本 | 通販 | Amazon
には、以下記載されています。引用します。
* バラツキの中から、有効な結果が得られたかどうかを判定するのが分散分析。
* 統計資料の1要因だけに着目して、資料への影響を分析する分析術が一元配置の分散分析。
* 統計資料の2要因の影響の有無を調べるのが二元配置の分散分析。


分散分析のできるライブラリ

調べたところ、以下で分散分析ができるようです。

  • spipy
    scipy.stats.f_oneway(*args) で、一元配置分散分析が可能です。

  • statsmodels
    statsmodels.api.stats.anova_lm() で、一元配置分散分析、二限配置分散分析が可能です。

  • pyvttbl
    python3.x では動作しなそうですが、pyvttbl というライブラリがあり、このライブラリでも分散分析が可能です。


参考

以下、スクリプト作成時に参考にした記事になります。
* pythonで一元配置分散分析(one way ANOVA) - 技術メモ
* Four ways to conduct one-way ANOVAs with Python - Erik Marsja
* https://www.marsja.se/three-ways-to-carry-out-2-way-anova-with-python/
* 私のための統計処理 ー基礎解説ーANOVA
* Anova in Python | plotly


実際にスクリプトを書いて試す

spipy、statsmodels で実際にスクリプトを書いて分散分析を行います。


spipy

spipy で 一元配置分散分析を行います。
使用するデータは、Rdatasets から拝借しました。

import pandas as pd
from scipy import stats
url = "https://vincentarelbundock.github.io/Rdatasets/csv/datasets/PlantGrowth.csv"
data = pd.read_csv(url)
data
Unnamed: 0weightgroup
014.17ctrl
125.58ctrl
235.18ctrl
346.11ctrl
454.50ctrl
564.61ctrl
675.17ctrl
784.53ctrl
895.33ctrl
9105.14ctrl
10114.81trt1
11124.17trt1
12134.41trt1
13143.59trt1
14155.87trt1
15163.83trt1
16176.03trt1
17184.89trt1
18194.32trt1
19204.69trt1
20216.31trt2
21225.12trt2
22235.54trt2
23245.50trt2
24255.37trt2
25265.29trt2
26274.92trt2
27286.15trt2
28295.80trt2
29305.26trt2
# CSVデータは縦持ちなので、横持ちのデータに変換する
ctrl = data['weight'][data.group == 'ctrl']
grps = pd.unique(data.group.values)
d_data = {grp:data['weight'][data.group == grp] for grp in grps}
d_data
{'ctrl': 0    4.17
 1    5.58
 2    5.18
 3    6.11
 4    4.50
 5    4.61
 6    5.17
 7    4.53
 8    5.33
 9    5.14
 Name: weight, dtype: float64, 'trt1': 10    4.81
 11    4.17
 12    4.41
 13    3.59
 14    5.87
 15    3.83
 16    6.03
 17    4.89
 18    4.32
 19    4.69
 Name: weight, dtype: float64, 'trt2': 20    6.31
 21    5.12
 22    5.54
 23    5.50
 24    5.37
 25    5.29
 26    4.92
 27    6.15
 28    5.80
 29    5.26
 Name: weight, dtype: float64}
# stats.f_onewayを実行して結果を出力する
result = stats.f_oneway(d_data['ctrl'], d_data['trt1'], d_data['trt2'])
print(result)
F_onewayResult(statistic=4.846087862380136, pvalue=0.015909958325622899)

statistic は、F値で、pvalue は P値です。
有意水準5%とすると、0.0159 <= 0.05 なので、集団間の平均値に差異があることが言えそうです。


statsmodels

statsmodels で、一元配置分散分析、二元配置分散分析を行います。

一元配置分散分析

import pandas as pd
from scipy import stats
url = "https://vincentarelbundock.github.io/Rdatasets/csv/datasets/PlantGrowth.csv"
data = pd.read_csv(url)
data
Unnamed: 0weightgroup
014.17ctrl
125.58ctrl
235.18ctrl
346.11ctrl
454.50ctrl
564.61ctrl
675.17ctrl
784.53ctrl
895.33ctrl
9105.14ctrl
10114.81trt1
11124.17trt1
12134.41trt1
13143.59trt1
14155.87trt1
15163.83trt1
16176.03trt1
17184.89trt1
18194.32trt1
19204.69trt1
20216.31trt2
21225.12trt2
22235.54trt2
23245.50trt2
24255.37trt2
25265.29trt2
26274.92trt2
27286.15trt2
28295.80trt2
29305.26trt2
import statsmodels.api as sm
from statsmodels.formula.api import ols

# 最小二乗法 ols 実行。 正直あまり、ここで何故登場するのか意味がわかっておりません。   
mod = ols('weight ~ group', data=data).fit()                
print(mod.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                 weight   R-squared:                       0.264
Model:                            OLS   Adj. R-squared:                  0.210
Method:                 Least Squares   F-statistic:                     4.846
Date:                Sun, 08 Jul 2018   Prob (F-statistic):             0.0159
Time:                        04:18:39   Log-Likelihood:                -26.810
No. Observations:                  30   AIC:                             59.62
Df Residuals:                      27   BIC:                             63.82
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
=================================================================================
                    coef    std err          t      P>|t|      [0.025      0.975]
---------------------------------------------------------------------------------
Intercept         5.0320      0.197     25.527      0.000       4.628       5.436
group[T.trt1]    -0.3710      0.279     -1.331      0.194      -0.943       0.201
group[T.trt2]     0.4940      0.279      1.772      0.088      -0.078       1.066
==============================================================================
Omnibus:                        1.835   Durbin-Watson:                   2.704
Prob(Omnibus):                  0.400   Jarque-Bera (JB):                1.406
Skew:                           0.524   Prob(JB):                        0.495
Kurtosis:                       2.835   Cond. No.                         3.73
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
# 分散分析 を実行
aov_table = sm.stats.anova_lm(mod, typ=2) # The type of Anova test to perform. See notes.
print(aov_table)
            sum_sq    df         F   PR(>F)
group      3.76634   2.0  4.846088  0.01591
Residual  10.49209  27.0       NaN      NaN

skipy と 同様の、F値、P値が取得できました。

二元配置分散分析

続いて、二元配置分散分析を実行します。
データは、繰り返しありのデータですので、繰り返しありの二元配置分散分析になります。

# 歯の成長データを取得
data = pd.read_csv('https://vincentarelbundock.github.io/Rdatasets/csv/datasets/ToothGrowth.csv')
data
Unnamed: 0lensuppdose
014.2VC0.5
1211.5VC0.5
237.3VC0.5
345.8VC0.5
456.4VC0.5
5610.0VC0.5
6711.2VC0.5
7811.2VC0.5
895.2VC0.5
9107.0VC0.5
101116.5VC1.0
111216.5VC1.0
121315.2VC1.0
131417.3VC1.0
141522.5VC1.0
151617.3VC1.0
161713.6VC1.0
171814.5VC1.0
181918.8VC1.0
192015.5VC1.0
202123.6VC2.0
212218.5VC2.0
222333.9VC2.0
232425.5VC2.0
242526.4VC2.0
252632.5VC2.0
262726.7VC2.0
272821.5VC2.0
282923.3VC2.0
293029.5VC2.0
303115.2OJ0.5
313221.5OJ0.5
323317.6OJ0.5
33349.7OJ0.5
343514.5OJ0.5
353610.0OJ0.5
36378.2OJ0.5
37389.4OJ0.5
383916.5OJ0.5
39409.7OJ0.5
404119.7OJ1.0
414223.3OJ1.0
424323.6OJ1.0
434426.4OJ1.0
444520.0OJ1.0
454625.2OJ1.0
464725.8OJ1.0
474821.2OJ1.0
484914.5OJ1.0
495027.3OJ1.0
505125.5OJ2.0
515226.4OJ2.0
525322.4OJ2.0
535424.5OJ2.0
545524.8OJ2.0
555630.9OJ2.0
565726.4OJ2.0
575827.3OJ2.0
585929.4OJ2.0
596023.0OJ2.0
import statsmodels.api as sm
from statsmodels.formula.api import ols
formula = 'len ~ C(supp) + C(dose) + C(supp):C(dose)'
model = ols(formula, data).fit()
print(model.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                    len   R-squared:                       0.794
Model:                            OLS   Adj. R-squared:                  0.775
Method:                 Least Squares   F-statistic:                     41.56
Date:                Sun, 08 Jul 2018   Prob (F-statistic):           2.50e-17
Time:                        04:18:54   Log-Likelihood:                -159.35
No. Observations:                  60   AIC:                             330.7
Df Residuals:                      54   BIC:                             343.3
Df Model:                           5                                         
Covariance Type:            nonrobust                                         
================================================================================================
                                   coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------------------------
Intercept                       13.2300      1.148     11.521      0.000      10.928      15.532
C(supp)[T.VC]                   -5.2500      1.624     -3.233      0.002      -8.506      -1.994
C(dose)[T.1.0]                   9.4700      1.624      5.831      0.000       6.214      12.726
C(dose)[T.2.0]                  12.8300      1.624      7.900      0.000       9.574      16.086
C(supp)[T.VC]:C(dose)[T.1.0]    -0.6800      2.297     -0.296      0.768      -5.285       3.925
C(supp)[T.VC]:C(dose)[T.2.0]     5.3300      2.297      2.321      0.024       0.725       9.935
==============================================================================
Omnibus:                        0.336   Durbin-Watson:                   2.025
Prob(Omnibus):                  0.846   Jarque-Bera (JB):                0.324
Skew:                           0.164   Prob(JB):                        0.850
Kurtosis:                       2.852   Cond. No.                         9.77
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
aov_table = sm.stats.anova_lm(model, typ=2)
print(aov_table)
                      sum_sq    df          F        PR(>F)
C(supp)           205.350000   1.0  15.571979  2.311828e-04
C(dose)          2426.434333   2.0  91.999965  4.046291e-18
C(supp):C(dose)   108.319000   2.0   4.106991  2.186027e-02
Residual          712.106000  54.0        NaN           NaN
import math
print("C(supp) pvalue", 2.311828 * pow(math.e, -4))
print("C(dose) pvalue", 4.046291 * pow(math.e, -18))
print("C(supp):C(dose) pvalue", 2.186027 * pow(math.e, -2))
C(supp) pvalue 0.042342606820864576
C(dose) pvalue 6.162492997121306e-08
C(supp):C(dose) pvalue 0.29584658320788276

C(supp) の p値は、0.04234 <= 5 で有意差があると言えそうで、
C(dose) の p値は、6.1625 > 5 で、有意差があるとは言えなそうです。
おそらく、suppと、dose の組み合わせを示すp値は、0.2958 <= 5 で有意差があると言えそうです。
supp のみが、有意差があると考えるのがいいのかもしれません。


まとめ

skipy と、statsmodel を使って、分散分析を行いました。 以下、感想まとめます。

  • 「統計用語 + python」 で検索すると、statsmodel 記事がよく引っかかる。
  • skipy と、statsmodel だったら、一元配置分散分析も、二元配置分散分析も、同じ手順で実装できる statsmodel の方が良いと思った。
  • 因子分析と、分散分析の違いがよくわからなくなったが、因子が決まっていてどれくらいの影響があるのかを調べるのが分散分析で、因子が何かわからない状況で探すのが因子分析。 と、個人的には理解した。

以上です。

コメント