R vs Python：統計するならどっちいいの？

データ解析をする上で、Rを使うべきかPythonを使うべきか、この議論は多くの人が色々な意見を持っています。最近はPythonユーザーが増えていますが、Rをメインで使う人が少なからずいるのもまた事実です。

今回は統計解析をするならどっち？という観点からRとPythonを比較したいと思います。ちなみ私はRユーザーなので、その辺りは加味した上でご覧ください^^;

Python統計に関してはこちらを参考にさせて頂きました。

Pythonで学ぶあたらしい統計学の教科書 (AI & TECHNOLOGY)

RやPythonでの統計解析の辞書代わりにもどうぞ！

　平均の比較：t検定

まずはt検定から

R

library(tidyverse)
#データの生成
a <- c(3.9,4.5,8.1,5.1,5.3,8.4,5.9,2.5,3.6,4.1)
b <- c(12.4,10.7,10.8,10.2,8.9,13.6,11,6.1,11.4,9.1)

>t.test(a,b,paired = F,var.equal = F,conf.level = 0.95,alternative = "two.sided")

     Welch Two Sample t-test

data:  a and b
t = -5.9581, df = 17.873, p-value = 1.265e-05
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -7.142764 -3.417236
sample estimates:
mean of x mean of y 
     5.14     10.42 

Rでt検定を行うにはt.test{stats}を使います。指定できるオプションとしては、pairedで対応の有無、var.equalで等分散性の仮定、conf.levelで信頼区間、alternativeで両側、右片測、左片測が指定できます。ちなみにデフォルトは、対応の無い、等分散性を仮定しないwelchのt検定となります。出力には、t統計量や自由度、p値、95%信頼区間、それぞれのグループの平均値などが出力されます。

お次はPythonです。

Python

import numpy as np
from scipy import stats

a = np.array([3.9,4.5,8.1,5.1,5.3,8.4,5.9,2.5,3.6,4.1])
b = np.array([12.4,10.7,10.8,10.2,8.9,13.6,11,6.1,11.4,9.1])

stats.ttest_ind(a,b,equal_var = False)
>Ttest_indResult(statistic=-5.958087910142942, pvalue=1.2647201884359637e-05)

Pythonではscipyのstatsをインポートして行います。指定できるのは等分散性の仮定のみで、デフォルトはTRUE、つまり等分散性を仮定するt検定です。ちなみに対応ありの場合は関数が異なり

#対応有り
stats.ttest_rel(a,b)

となります。信頼区間の算出はややこしく、t分布に平均と分散、自由度を指定して求める必要があります。

bottom,up=stats.t.interval(
    alpha=0.95, # 信頼区間
    loc=a.mean()-b.mean(),  # 標本平均
    scale=0.8861903, #不偏分散
    df=17.873 # 自由度
)
print('95%信頼区間: {:.2f} < μ < {:.2f}'.format(bottom, up))

>95%信頼区間: -7.14 < μ < -3.42

対応の有無、等分散の仮定の違いによって分散と自由度が異なるため、Pythonで正確な信頼区間を算出するのは、骨が折れそうです。

一般化線形モデルによるt検定

少しAdvanceなやり方ですが、一般化線形モデルを使って対応の無いt検定を行うこともできます。線形モデルでの平均値の比較については、こちらがわかり易いです。

統計学が最強の学問である[実践編]――データ分析のための思想と方法

もしくはこちらを

www.medi-08-data-06.work

R

X <-  rep(c(1,0),c(length(a),length(b)))
y <- c(a,b)

lm_mod <- lm(y~X)

>lm_mod 
Call:
lm(formula = y ~ X)

Coefficients:
(Intercept)            X  
      10.42        -5.28  


>confint(lm_mod)
                2.5 %    97.5 %
(Intercept)  9.103497 11.736503
X           -7.141817 -3.418183

CoefficientのXが2群の平均値の差、confint()の出力のXが95%信頼区間です。先ほどのt検定と比較しても、同じような結果になっていますね。

Python

import statsmodels.api as sm

y = np.r_[a,b]
X = np.c_[np.ones(len(a)+len(b)),np.repeat([1,0],[len(a),len(b)])]
lm_mod = sm.OLS(y,X).fit()
lm_mod.summary().tables[1]

	coef	std err	t	P>|t|	[0.025	0.975]
const	10.4200	0.627	16.629	0.000	9.103	11.737
x1	-5.2800	0.886	-5.958	0.000	-7.142	-3.418

Pythonで線形モデルを作るにはstatsmodels.apiを使います。また、Xには切片項の1を加えることを忘れないようにしてください。Pythonでも結果は同じですね。

もし、切片項を追加するのが面倒であれば、from_formulaメソッドを使って、Rライクに書くこともできます。

import pandas as pd

#データフレーム にする必要がある
data = pd.DataFrame({"y":np.r_[a,b],
                     "X" : np.repeat([1,0],[len(a),len(b)])})

lm_mod = sm.OLS.from_formula("y~X",data=data).fit()
lm_mod.summary().tables[1]

	coef	std err	t	P>|t|	[0.025	0.975]
Intercept	10.4200	0.627	16.629	0.000	9.103	11.737
X	-5.2800	0.886	-5.958	0.000	-7.142	-3.418

どちらでも結果は同じです。pandasを使ってデータフレーム にする必要がありますが、Rに慣れている私はこちらの方が好きです^^;

　t検定まとめ

以上から、結果の出力や指定できるオプションの使いやすさなどは、Rが優勢なイメージです。ただ、点推定値や大まかな検定のみであればPythonでも行うことができますね。

平均値の比較：一元配置分散分析

t検定ときたら、やはり次は３群以上の比較を行う分散分析です。まずはRからです。

#分散分析~oneway~
#データの作成
a <- c(3.9,4.5,8.1,5.1,5.3,8.4,5.9,2.5,3.6,4.1)
b <- c(12.4,10.7,10.8,10.2,8.9,13.6,11,6.1,11.4,9.1)
c <- c(5.9,7.6,5.9,6.5,6.7,4.6,9.7,8.3,5.7,10.5)
id = rep(c(0,1,2),c(length(a),length(b),length(c)))

data <-  tibble(y=c(a,b,c),ID=as.character(id))

>data
# A tibble: 6 x 2
      y ID   
  <dbl> <chr>
1   3.9 0    
2   4.5 0    
3   8.1 0    
4   5.1 0    
5   5.3 0    
6   8.4 0    

分散分析はデータフレーム になっていた方が使いやすいので、 データフレーム にしてあります。IDがそれぞれの群を表していて、1~3の3群間の平均値の比較です。Rで分散分析を行う一般的な方法はaov{stats}で分散分析表を作り、Anova{car}で検定を行うやり方です。

R

library(car)
#分散分析表を作る
> aov(data = data,y ~ ID)
Call:
   aov(formula = y ~ ID, data = data)

Terms:
                      ID Residuals
Sum of Squares  142.1227  102.2840
Deg. of Freedom        2        27

Residual standard error: 1.946355
Estimated effects may be unbalanced

#検定
> Anova(aov(data = data,y ~ ID))
Anova Table (Type II tests)

Response: y
          Sum Sq Df F value    Pr(>F)    
ID        142.12  2  18.758 7.814e-06 ***
Residuals 102.28 27                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Anova{car}を使わなくてもデフォルトのanova{stats}でもできるのですが、anova{stats}はtype1Anovaにしか対応していないため、type2またはtype3のAnovaをするためには必要になります。

多重比較をする場合はいくつかやり方がありますが、有名なボンフェローニの検定を行うにはpairwise.t.test{stats}を使います。

> pairwise.t.test(x = data$y, g=data$ID, p.adj = "bonf")

    Pairwise comparisons using t tests with pooled SD 

data:  data$y and data$ID 

  0       1     
1 5.3e-06 -     
2 0.0887  0.0024

P value adjustment method: bonferroni 

0群と1群、1群と2群の間に有意差があるようです。ボンフェローニの他にも、p.adjには様々な多重比較の補正方法を指定できます。

また、t検定同様分散分析を線形モデルで表現することもできます。

> Anova(lm(data = data,y~ID))
Anova Table (Type II tests)

Response: y
          Sum Sq Df F value    Pr(>F)    
ID        142.12  2  18.758 7.814e-06 ***
Residuals 102.28 27                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

結果は同じですね！

Pyhthonではこちらの表現に近い記述になります。Pythonで分散分析をやってみます。

Python

anova = sm.OLS.from_formula("y~ID",data=data).fit()
sm.stats.anova_lm(anova,type=2)

	df	sum_sq	mean_sq	F	PR(>F)
ID	2.0	142.122667	71.061333	18.758124	0.000008
Residual	27.0	102.284000	3.788296	NaN	NaN

まずは、statsmodels.apiのOLSで線形モデルを作り、anova_lmメソッドを使って検定を行うことができます。typeには1,2,3とあるのですが、その違いは上記参考書をご覧ください。結果はRと同じですね。

Pythonでの多重比較は、Turkeyの多重比較法になります。(ボンフェローニは情報がありませんでした...、p値を水準数で割るだけのですが)

print(sm.stats.multicomp.pairwise_tukeyhsd(data.y,data.ID))

Multiple Comparison of Means - Tukey HSD,FWER=0.05

---

group1 group2 meandiff  lower   upper  reject

---------------------------------------------
  0      1      5.28    3.1222  7.4378  True 
  0      2      2.0    -0.1578  4.1578 False 
  1      2     -3.28   -5.4378 -1.1222  True 
---------------------------------------------

結果はRと同じように0群と1群、1群と2群の間に有意差があるようです。

平均値の比較：二元配置分散分析

お次は二元配置分散分析です。やり方は簡単で、先ほどの一元配置分散分析の式に交互作用項を加えるだけです。

R

#もう一水準追加
id2 = rep(c("A","B"),15)

data <- data %>% 
  mutate(ID2 = id2)

> Anova(aov(data=data,y~ID*ID2))
#もしくは線形モデルでの書き方
#Anova(lm(data=data,y~ID*ID2))

Anova Table (Type II tests)

Response: y
           Sum Sq Df F value    Pr(>F)    
ID        142.123  2 17.3673 2.167e-05 ***
ID2         0.385  1  0.0942    0.7616    
ID:ID2      3.699  2  0.4520    0.6417    
Residuals  98.200 24                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

ID*ID2はID+ID2+ID:ID2を表し、交互作用項を含めたモデルになります。

Pythonでも記述方法はほとんど同じです。

Python

#もう一水準追加
ID2 = np.repeat(["A","B"],15)
data["ID2"] = ID2

#2元配置分散分析
anova_two = sm.OLS.from_formula("y~ID+ID2+ID*ID2",data=data).fit()
sm.stats.anova_lm(anova_two,type=2)

	df	sum_sq	mean_sq	F	PR(>F)
ID	2.0	142.122667	71.061333	18.120779	0.000012
ID2	1.0	0.324000	0.324000	0.082621	0.776056
ID*ID2	2.0	10.726240	5.363120	1.367606	0.272437
Residual	26.0	101.960000	3.921538	NaN	NaN

こちらもRと結果は同じですね！

分散分析まとめ

単に分散分析をするだけであれば、RもPythonもそれほど違いはなさそうです。ただ、今回は全て対応がないことを前提としましたが、対応の有無や、多重比較の細やかな指定などが入ってくるとPythonでは苦しい印象です。

単回帰分析・重回帰分析

続いて、統計では最もよく使われる単回帰分析・重回帰分析です。どちらも先ほどの線形モデルによる分散分析と記述法はほとんど同じです。

R

#データの作成

x1 <- c(5,4.6,2.3,3.8,5.6,5.8,2.6,4.3,1.7,4.5)
x2 <- c(13.5,11.8,6.8,12.9,11.7,8.4,3.2,7,12.6,10.4)
x3 <- c(8.2,3.4,8,3.6,6.2,8.9,7.9,7.4,9.7,9.2)
error <- c(-1.9,-0.1,1,0.2,-1.4,-1.4,0.4,-1.8,-0.3,-0.7)

y <- 5*x1+10*x2+x3+error

data <- tibble(y=y,
               x1 = x1,
               x2 = x2,
               x3 = x3)


#単回帰
>lm(y~x1,data=data) %>% 
  summary()

Call:
lm(formula = y ~ x1, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-56.110 -22.517   6.006  25.095  44.402 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   80.775     33.772   2.392   0.0437 *
x1            11.014      7.976   1.381   0.2047  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 33.53 on 8 degrees of freedom
Multiple R-squared:  0.1925, Adjusted R-squared:  0.09153 
F-statistic: 1.907 on 1 and 8 DF,  p-value: 0.2047


#重回帰
>lm(y~x1+x2+x3,data=data) %>% 
  summary()

Call:
lm(formula = y ~ x1 + x2 + x3, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.1898 -0.2333  0.0449  0.3491  0.6329 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.64616    1.30930   2.785 0.031794 *  
x1           4.45684    0.16512  26.992 1.71e-07 ***
x2           9.93916    0.07004 141.897 8.26e-12 ***
x3           0.79799    0.10485   7.611 0.000268 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6572 on 6 degrees of freedom
Multiple R-squared:  0.9998, Adjusted R-squared:  0.9997 
F-statistic:  8596 on 3 and 6 DF,  p-value: 2.753e-11

Python

#データの作成
x1 = np.array([5,4.6,2.3,3.8,5.6,5.8,2.6,4.3,1.7,4.5])
x2 = np.array([13.5,11.8,6.8,12.9,11.7,8.4,3.2,7,12.6,10.4])
x3 = np.array([8.2,3.4,8,3.6,6.2,8.9,7.9,7.4,9.7,9.2])
error = np.array([-1.9,-0.1,1,0.2,-1.4,-1.4,0.4,-1.8,-0.3,-0.7])

y = 5*x1+10*x2+x3+error

data = pd.DataFrame({"y":y,
               "x1" : x1,
               "x2" : x2,
               "x3" : x3})

#単回帰
lm_mod = sm.OLS.from_formula("y~x1",data=data).fit()
lm_mod.summary().tables[1]

	coef	std err	t	P>|t|	[0.025	0.975]
Intercept	80.7745	33.772	2.392	0.044	2.895	158.654
x1	11.0138	7.976	1.381	0.205	-7.379	29.407

#重回帰
lm_mod = sm.OLS.from_formula("y~x1+x2+x3",data=data).fit()
lm_mod.summary().tables[1]

	coef	std err	t	P>|t|	[0.025	0.975]
Intercept	3.6462	1.309	2.785	0.032	0.442	6.850
x1	4.4568	0.165	26.992	0.000	4.053	4.861
x2	9.9392	0.070	141.897	0.000	9.768	10.111
x3	0.7980	0.105	7.611	0.000	0.541	1.055

ちなみに重回帰分析を　Pythonでよく使われる記法で書くと

X_t = data.drop("y",axis=1)
X_t = np.c_[np.ones(X_t.shape[0]),X_t]
y_t = data.y

lm_mod = sm.OLS(y_t,X_t,data=data).fit()

と書くことができます。こうすると何が良いかというと説明変数が膨大にあるときにはRではy~x1+x2+.....としなければいけないところをPythonでは一発で書くことができます。この辺りはビックデータ解析を得意とするPythonの強みですね。下記追記参照

　単回帰・重回帰まとめ

単回帰・重回帰はモデルが単純なだけにRもPythonも使いやすさに違いは感じませんでした。単純な解析のみであれば、どちらを使っても良いでしょう。

---

補足：Rでたくさんの説明変数がある場合の小技

Rで重回帰分析を行う際にはy~x1+x2+...と書く必要があります。説明変数の数が少なければいいですが、膨大になったときには大変です。これを簡単に書く方法として私がよく使うのは、列名を+で繋げるという小技です。実際にやってみましょう。

data %>% 
  colnames() %>%
  str_subset("x.") %>% 
  str_c(collapse = "+")

> "x1+x2+x3"

colnames()で列名を抜き出し、str_subset()で目的の説明変数をぬきだす、そして最後はstr_c()を使って+で繋げるというものです。後は出力結果をコピペして完了です。

少し無理矢理感はありますが、結構便利です笑。もっとスマートなやり方があれば教えて頂きたいですm(__)m

2019/07/08 追記

上記に関してコメントを頂きました。

Rではモデル式の右辺に `.` を使うと左辺に指定した変数以外のx1+x2+x3+...+xnになります．さらに - x1などとすることで，使いたくない変数を外せます．
lm(Sepal.Length ~ . - Species, data = iris)

— atusy (@Atsushi776) 2019年7月7日

やってみます。

> lm(data=data,y~.) %>% 
+   summary()

Call:
lm(formula = y ~ ., data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.1898 -0.2333  0.0449  0.3491  0.6329 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.64616    1.30930   2.785 0.031794 *  
x1           4.45684    0.16512  26.992 1.71e-07 ***
x2           9.93916    0.07004 141.897 8.26e-12 ***
x3           0.79799    0.10485   7.611 0.000268 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6572 on 6 degrees of freedom
Multiple R-squared:  0.9998, Adjusted R-squared:  0.9997 
F-statistic:  8596 on 3 and 6 DF,  p-value: 2.753e-11

おお！これで私の小技も必要なくなりました（笑）

---

　ロジスティック回帰

最後はロジスティック回帰です。これは一般化線形モデルの確率分布を変化させることで実装できます。

R

#データの作成
x <- c(62,53,46,58,47,56,52,27,36,39,45,62,56,40,57,63,63,74,87,77,79,75,69,81,64,78,52,63,93,84)
y <- c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)

data <- tibble(y=y,
               x = x)


#ロジスティック回帰

>glm(data=data,y~x,family = "binomial") %>% 
  summary()

Call:
glm(formula = y ~ x, family = "binomial", data = data)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.35197  -0.31486   0.00094   0.14479   2.34600  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) -18.7406     8.1897  -2.288   0.0221 *
data$x        0.3087     0.1353   2.282   0.0225 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 41.589  on 29  degrees of freedom
Residual deviance: 14.915  on 28  degrees of freedom
AIC: 18.915

Number of Fisher Scoring iterations: 7

Python

#データの生成
x = np.array([62,53,46,58,47,56,52,27,36,39,45,62,56,40,57,63,63,74,87,77,79,75,69,81,64,78,52,63,93,84])
y = np.array([0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1])

data = pd.DataFrame({"y":y,
               "x" : x})


logit_mod = sm.GLM.from_formula("y~x",
                         data=data,
                         family=sm.families.Binomial()).fit()
logit_mod.summary().tables[1]

	coef	std err	z	P>|z|	[0.025	0.975]
Intercept	-18.7406	8.190	-2.288	0.022	-34.792	-2.689
x	0.3087	0.135	2.282	0.022	0.044	0.574

Rではglm{stats}を使って、引数familyにbinomialを指定する、PythonではOLSではなく、GLMを使ってこちらもfamilyにBinomial指定することでロジスティック回帰分析を行うことができます。

どちらも違いはありませんね。RでもPythonでも、familyやリンク関数を指定して、ポワソン回帰などその他のモデルを作ることもできます。

ロジスティック回帰まとめ

ロジスティック回帰も単回帰・重回帰同様簡単なモデルであれば、R、Python共に同じような記述量と出力結果となるようです。また確率分布やリンク関数の設定などもPythonでできることは少し驚きでした....

で、結局どっち？

実は私は統計なら断然Rが優勢と思っていたのですが、Pythonでも基本的な検定は問題なく行えるということを実感しました 。研究のような統計解析の細かなプロセスが求められる場合はやはりRを使うべきだと思いますが、そうでない場合はPythonでもほとんど問題なく統計解析ができそうです。

で結局どっちが良いかというと、ズバリ時と場合による！が正直なところです。ただ、統計解析に関するレファレンスの多さや、少し高度な統計モデルとなるとやはりRが優勢と言ったところでしょうか。またRではデフォルトで読み込まれている{stats}でほとんどのことができてしまいますが、Pythonでは色々な部品を読み込む必要があります。そもそもPythonは解析そのものよりも、何かに組み込んで使うプログラミング言語であるというルーツの違いもありますね。

RとPythonのどちらを使うべきかという議論よりも、目的によってどちらが合っているのか？というところから考えた方が良さそうです^^;

それを考える上で、この記事がお役に立てれば本望です^^

※本記事は筆者が個人的に学んだことをまとめた記事なります。所属する組織の意見・見解とは無関係です。また、数学の記法や詳細な理論、用語等で誤りがあった際はご指摘頂けると幸いです。

参考

Pythonで統計を行う記述法だけでなく、シミュレーションを通して、統計学の基本事項から学べると言ったとても豪華な内容です。Pythonで統計をしたいという方にはオススメです！

Pythonで学ぶあたらしい統計学の教科書 (AI & TECHNOLOGY)