RとPythonで重回帰 再追記

#はじめに
R未経験にもかかわらずR対象のデータ分析系の勉強会に出席してしまったので
Pythonだとどうするの？を復習を兼ねて記録に残したいと思います

#改訂2
PythonでもRのように、目的変数と説明変数を指定することができました
備忘録として追記しました

#改訂
勉強会参加メンバーから疑問点について助言いただいたので文末に追記しました

python で決定係数を出力した際、scikit learnとstatsmodelsに差がありました
statsmodelsの設定により差が発生していたことがわかりました

statsmodelsのデフォルトは、切片0とのこと
A intercept is not included by default and should be added by the user.
See statsmodels.tool.add_constant.

↓公式documents
https://www.statsmodels.org/dev/generated/statsmodels.regression.linear_model.OLS.html

別件でscikit-learnの　intercept = False　も結果がおかしくなるとの噂も、、、
この辺は別にまとめたいと思っています

#環境
windows10 Pro 64bit
R :3.6.1
RStudio :1.2.1335

Google Colaboratoryを使用
Python :3.6.8

#何するの？
RのSMCRM内のcustomerAcquisitionデータセットを利用して重回帰分析を行う

```r

library(SMCRM)

#データセットの読み出し
data("customerAcquisition")

#データセットの確認
str(customerAcquisition)

'data.frame':	500 obs. of  16 variables:
 $ customer      : num  1 2 3 4 5 6 7 8 9 10 ...
 $ acquisition   : num  1 0 0 1 1 0 0 0 1 1 ...
 $ first_purchase: num  434 0 0 226 363 ...
 $ clv           : num  0 0 0 5.73 0 ...
 $ duration      : num  384 0 0 730 579 0 0 0 730 730 ...
 $ censor        : num  0 0 0 1 0 0 0 0 1 1 ...
 $ acq_expense   : num  760 148 253 610 672 ...
 $ acq_expense_sq: num  578147 21815 63787 371771 452068 ...
 $ industry      : num  1 1 1 1 1 0 0 0 1 1 ...
 $ revenue       : num  30.2 39.8 54.9 45.8 69 ...
 $ employees     : num  1240 166 1016 122 313 ...
 $ ret_expense   : num  2310 0 0 2193 801 ...
 $ ret_expense_sq: num  5335130 0 0 4807451 641825 ...
 $ crossbuy      : num  5 0 0 2 4 0 0 0 1 1 ...
 $ frequency     : num  2 0 0 12 7 0 0 0 11 14 ...
 $ frequency_sq  : num  4 0 0 144 49 0 0 0 121 196 ...

# acquisition が 1のデータを対象としてdf に代入
df<- customerAcquisition %>% filter(acquisition == 1)

# first_purchase を目的関数、その他すべてのデータを説明関数として重回帰
lm_CA_full <- lm(first_purchase ~ . - acquisition , data = df)

# first_purchase を目的関数、説明関数を4つに絞って重回帰
lm_CA <- lm(first_purchase ~ acq_expense + acq_expense_sq + revenue + employees, data = df)

# lm_CAの回帰計算結果表示
summary(lm_CA)

Call:
lm(formula = first_purchase ~ acq_expense + acq_expense_sq + 
    revenue + employees, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-84.734 -19.765  -0.867  19.986 100.291 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)     7.783e+00  3.150e+01   0.247    0.805    
acq_expense     5.317e-01  9.534e-02   5.577 5.66e-08 ***
acq_expense_sq -6.544e-04  7.192e-05  -9.099  < 2e-16 ***
revenue         2.986e+00  1.086e-01  27.488  < 2e-16 ***
employees       2.483e-01  3.769e-03  65.887  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 30.59 on 287 degrees of freedom
Multiple R-squared:  0.9617,	Adjusted R-squared:  0.9611 
F-statistic:  1800 on 4 and 287 DF,  p-value: < 2.2e-16
```

説明係数の回帰係数は以下
　acq_expense     5.317e-01
　acq_expense_sq -6.544e-04
　revenue         2.986e+00
　employees       2.483e-01
切片は
　(Intercept)     7.783e+00

線形回帰モデルは以下で表せる
first_purchase = 7.783e+00 
                 　　　　　　　+ acq_expense x 5.317e-01 
                　　　　　　　 + acq_expense_sq x -6.544e-04 
                 　　　　　　　+ revenue x 2.986e+00
                　　　　　　　 + employees x 2.483e-01


モデルの評価値は決定係数で表される
決定係数　Multiple R-squared:  0.9617
修正決定係数　Adjusted R-squared:  0.9611

モデル同士の比較はAIC（赤池情報量基準）にて評価

```r

> # AICで比較
> AIC(lm_CA_full)
[1] 2842.858
> AIC(lm_CA)
[1] 2833.373
```
AICの値が小さいほうが理想に近いということらしいので、すべての説明変数を利用したモデルより
4つの説明変数のモデルのほうがよく説明できているモデルということになる。
ではこれがベストか？はstepAICで調べられる。（らしい、、）

```r


> lm_CA_step <- stepAIC(lm_CA_full)
Start:  AIC=2012.2
first_purchase ~ (customer + acquisition + clv + duration + censor + 
    acq_expense + acq_expense_sq + industry + revenue + employees + 
    ret_expense + ret_expense_sq + crossbuy + frequency + frequency_sq) - 
    acquisition

                 Df Sum of Sq     RSS    AIC
- customer        1       238  259372 2010.5
- crossbuy        1       259  259394 2010.5
- duration        1       555  259690 2010.8
- frequency_sq    1       621  259756 2010.9
- ret_expense_sq  1       712  259847 2011.0
- ret_expense     1      1414  260549 2011.8
<none>                         259135 2012.2
- frequency       1      2075  261210 2012.5
- censor          1      2283  261418 2012.8
- clv             1      3111  262246 2013.7
- industry        1      3131  262266 2013.7
- acq_expense     1     23542  282677 2035.6
- acq_expense_sq  1     47453  306588 2059.3
- revenue         1    536393  795528 2337.7
- employees       1   3231713 3490848 2769.6

Step:  AIC=2010.47
first_purchase ~ clv + duration + censor + acq_expense + acq_expense_sq + 
    industry + revenue + employees + ret_expense + ret_expense_sq + 
    crossbuy + frequency + frequency_sq

                 Df Sum of Sq     RSS    AIC
- crossbuy        1       213  259586 2008.7
- duration        1       536  259909 2009.1
- frequency_sq    1       568  259941 2009.1
- ret_expense_sq  1       737  260109 2009.3
- ret_expense     1      1448  260820 2010.1
<none>                         259372 2010.5
- frequency       1      1981  261354 2010.7
- censor          1      2200  261573 2010.9
- industry        1      3009  262381 2011.8
- clv             1      3038  262411 2011.9
- acq_expense     1     23644  283017 2033.9
- acq_expense_sq  1     47614  306987 2057.7
- revenue         1    538374  797746 2336.5
- employees       1   3277066 3536439 2771.3

Step:  AIC=2008.71
first_purchase ~ clv + duration + censor + acq_expense + acq_expense_sq + 
    industry + revenue + employees + ret_expense + ret_expense_sq + 
    frequency + frequency_sq

                 Df Sum of Sq     RSS    AIC
- duration        1       394  259980 2007.2
- frequency_sq    1       496  260082 2007.3
- ret_expense_sq  1       812  260398 2007.6
- ret_expense     1      1428  261014 2008.3
<none>                         259586 2008.7
- frequency       1      1803  261388 2008.7
- censor          1      2000  261586 2009.0
- industry        1      2810  262396 2009.8
- clv             1      2890  262475 2009.9
- acq_expense     1     24239  283825 2032.8
- acq_expense_sq  1     50080  309666 2058.2
- revenue         1    554513  814098 2340.5
- employees       1   3297327 3556912 2771.0

Step:  AIC=2007.15
first_purchase ~ clv + censor + acq_expense + acq_expense_sq + 
    industry + revenue + employees + ret_expense + ret_expense_sq + 
    frequency + frequency_sq

                 Df Sum of Sq     RSS    AIC
- frequency_sq    1       495  260474 2005.7
- ret_expense_sq  1       605  260585 2005.8
- ret_expense     1      1068  261048 2006.3
- censor          1      1606  261586 2007.0
- frequency       1      1609  261589 2007.0
<none>                         259980 2007.2
- industry        1      2419  262399 2007.8
- clv             1      2502  262481 2007.9
- acq_expense     1     31416  291395 2038.5
- acq_expense_sq  1     77789  337769 2081.6
- revenue         1    594439  854419 2352.6
- employees       1   3327401 3587380 2771.5

Step:  AIC=2005.7
first_purchase ~ clv + censor + acq_expense + acq_expense_sq + 
    industry + revenue + employees + ret_expense + ret_expense_sq + 
    frequency

                 Df Sum of Sq     RSS    AIC
- ret_expense_sq  1       651  261125 2004.4
- ret_expense     1      1052  261527 2004.9
- censor          1      1382  261856 2005.2
<none>                         260474 2005.7
- clv             1      2202  262677 2006.2
- industry        1      2341  262815 2006.3
- frequency       1      5221  265695 2009.5
- acq_expense     1     31388  291862 2036.9
- acq_expense_sq  1     77902  338376 2080.1
- revenue         1    612592  873066 2356.9
- employees       1   3335399 3595873 2770.2

Step:  AIC=2004.43
first_purchase ~ clv + censor + acq_expense + acq_expense_sq + 
    industry + revenue + employees + ret_expense + frequency

                 Df Sum of Sq     RSS    AIC
- ret_expense     1       796  261921 2003.3
- censor          1      1429  262554 2004.0
<none>                         261125 2004.4
- clv             1      2225  263349 2004.9
- industry        1      2402  263527 2005.1
- frequency       1      5500  266625 2008.5
- acq_expense     1     30890  292015 2035.1
- acq_expense_sq  1     77273  338398 2078.1
- revenue         1    621020  882145 2357.9
- employees       1   3347586 3608711 2769.3

Step:  AIC=2003.32
first_purchase ~ clv + censor + acq_expense + acq_expense_sq + 
    industry + revenue + employees + frequency

                 Df Sum of Sq     RSS    AIC
- censor          1       929  262850 2002.3
- clv             1      1489  263410 2003.0
<none>                         261921 2003.3
- industry        1      2086  264008 2003.6
- frequency       1      4822  266743 2006.7
- acq_expense     1     30126  292047 2033.1
- acq_expense_sq  1     76673  338594 2076.3
- revenue         1    674215  936136 2373.2
- employees       1   3509176 3771097 2780.1

Step:  AIC=2002.35
first_purchase ~ clv + acq_expense + acq_expense_sq + industry + 
    revenue + employees + frequency

                 Df Sum of Sq     RSS    AIC
- clv             1      1250  264100 2001.7
<none>                         262850 2002.3
- industry        1      1902  264752 2002.5
- frequency       1      4941  267791 2005.8
- acq_expense     1     29200  292049 2031.1
- acq_expense_sq  1     76005  338855 2074.5
- revenue         1    679663  942512 2373.2
- employees       1   3892379 4155229 2806.4

Step:  AIC=2001.74
first_purchase ~ acq_expense + acq_expense_sq + industry + revenue + 
    employees + frequency

                 Df Sum of Sq     RSS    AIC
- industry        1      1233  265333 2001.1
<none>                         264100 2001.7
- frequency       1      3738  267838 2003.8
- acq_expense     1     28352  292452 2029.5
- acq_expense_sq  1     76467  340567 2074.0
- revenue         1    693588  957688 2375.9
- employees       1   3962665 4226765 2809.4

Step:  AIC=2001.1
first_purchase ~ acq_expense + acq_expense_sq + revenue + employees + 
    frequency

                 Df Sum of Sq     RSS    AIC
<none>                         265333 2001.1
- frequency       1      3303  268636 2002.7
- acq_expense     1     28868  294201 2029.3
- acq_expense_sq  1     77127  342460 2073.6
- revenue         1    706100  971433 2378.1
- employees       1   3973948 4239281 2808.3
```
StartのAICが2012.2で単独で評価した2842.858となぜ異なるのかわかりませんが、ベストなモデルは
初めに選択した4つの説明変数に、frequencyを追加したものということが分かりました。

# Pythonで記述
customerAcquisitionデータセットはPythonには用意されていないので、Rからcsvで出力して利用することにします

```python

＃RからcustomerAcquisitionをcsvで出力
write.table(customerAcquisition, "customerAcquisition.csv", sep=",", row.names=F)

```
Pythonでのデータ解析はライブラリーに依存しているところが多々あるため下記のようなライブラリーをインポートする必要がある
pandas  　　　データ分析用ライブラリ。DataFrameを提供
numpy　　　　多次元配列計算用ライブラリ
matplotlib　　グラフ描画ライブラリ
seaborn　　　matplotlibのラッパー、かっこいいグラフが簡単に書ける
sklearn  　　　機械学習用ライブラリー

```python

import numpy as np
import pandas as pd

# 可視化モジュール
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
%matplotlib inline

#重回帰
from sklearn.linear_model import LinearRegression
```

Rから出力したcsvファイルはpandasで読み込み、データ構造も簡単に確認できます

```python

df = pd.read_csv('customerAcquisition.csv')
df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 500 entries, 0 to 499
Data columns (total 16 columns):
customer          500 non-null int64
acquisition       500 non-null int64
first_purchase    500 non-null float64
clv               500 non-null float64
duration          500 non-null int64
censor            500 non-null int64
acq_expense       500 non-null float64
acq_expense_sq    500 non-null float64
industry          500 non-null int64
revenue           500 non-null float64
employees         500 non-null int64
ret_expense       500 non-null float64
ret_expense_sq    500 non-null float64
crossbuy          500 non-null int64
frequency         500 non-null int64
frequency_sq      500 non-null int64
dtypes: float64(7), int64(9)
memory usage: 62.6 KB
```

データフレームは論理式でbool値を取得し、bool値を使ってデータを抽出することができる

```python

df['acquisition'] == 1

0       True
1      False
2      False
3       True
  ... 
497    False
498     True
499    False
Name: acquisition, Length: 500, dtype: bool

# acquisition が 1　のデータのみに抽出
df = df[df['acquisition'] == 1]

#抽出後のデータを確認
df.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 292 entries, 0 to 498
Data columns (total 16 columns):
customer          292 non-null int64
acquisition       292 non-null int64
first_purchase    292 non-null float64
clv               292 non-null float64
duration          292 non-null int64
censor            292 non-null int64
acq_expense       292 non-null float64
acq_expense_sq    292 non-null float64
industry          292 non-null int64
revenue           292 non-null float64
employees         292 non-null int64
ret_expense       292 non-null float64
ret_expense_sq    292 non-null float64
crossbuy          292 non-null int64
frequency         292 non-null int64
frequency_sq      292 non-null int64
dtypes: float64(7), int64(9)
memory usage: 38.8 KB
```
acquisition　= 1　のみのデータ抽出でデータ数が
500 →　292 に減少

このデータフレームをベースに目的変数、説明変数を作りながら
重回帰のclassにデータを入れていきます

```python

#目的変数
y = df['first_purchase']

#説明変数は元のデータフレームから目的変数を削除したもの
X_full = df.drop('first_purchase',axis=1)

#選択した説明変数のデータフレームも作成
df_sel = df[['first_purchase','acq_expense','acq_expense_sq','revenue','employees']]
X_sel = df_sel.drop('first_purchase',axis=1)

#重回帰のインスタンスを作成
model = LinearRegression()

#X_selで学習
model.fit(X_sel,y)

#回帰係数
print('\n回帰係数\n{}'.format(pd.Series(model.coef_,　index=df_sel.drop('first_purchase',axis=1).columns)))

回帰係数
acq_expense       0.531689
acq_expense_sq   -0.000654
revenue           2.986234
employees         0.248334
dtype: float64

#切片
print('切片: {:.3f}'.format(model.intercept_))

切片: 7.783

#決定係数
print(model.score(df_sel.drop('first_purchase',axis=1),df_sel['first_purchase']))

0.962
```
回帰係数、切片ともにRと同様に求められたが
評価値は決定係数のみで修正決定係数は出力できない（と思う、、追記に記載しました）
Pythonの重回帰は機械学習モデルとしてLinearRegressionしか使ったことがなかったが、Rのように統計的な結出力を得たい場合は、statsmodelsで可能となるようです。

```python

from sklearn import datasets
from statsmodels import api as sm

model = sm.OLS(y, X_sel)
result = model.fit()
print (result.summary())

                                 OLS Regression Results                                
=======================================================================================
Dep. Variable:         first_purchase   R-squared (uncentered):                   0.994
Model:                            OLS   Adj. R-squared (uncentered):              0.994
Method:                 Least Squares   F-statistic:                          1.266e+04
Date:                Tue, 13 Aug 2019   Prob (F-statistic):                   3.21e-322
Time:                        14:44:24   Log-Likelihood:                         -1410.7
No. Observations:                 292   AIC:                                      2829.
Df Residuals:                     288   BIC:                                      2844.
Df Model:                           4                                                  
Covariance Type:            nonrobust                                                  
==================================================================================
                     coef    std err          t      P>|t|      [0.025      0.975]
----------------------------------------------------------------------------------
acq_expense        0.5546      0.023     24.391      0.000       0.510       0.599
acq_expense_sq    -0.0007   2.42e-05    -27.743      0.000      -0.001      -0.001
revenue            2.9894      0.108     27.755      0.000       2.777       3.201
employees          0.2486      0.004     69.021      0.000       0.242       0.256
==============================================================================
Omnibus:                        2.317   Durbin-Watson:                   1.891
Prob(Omnibus):                  0.314   Jarque-Bera (JB):                2.013
Skew:                           0.156   Prob(JB):                        0.366
Kurtosis:                       3.260   Cond. No.                     2.97e+04
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.97e+04. This might indicate that there are
strong multicollinearity or other numerical problems.


print(result.aic)

2829.4353329751475

```

修正決定係数まで出力できたが、、、結果が異なる、、
AICも微妙に違う、、なぜ？
Pythonでは、調べた限りではRのようなstepAICはなく自分で組む必要がありそう
Qiitaでも自作している人がちらほら
データが異なる理由、stepAICはじっくり調べてみたいと思います

とりあえず、Rとpythonで同じデータセットを用いて重回帰を行ってみました
Rは初めてでしたが、さすが統計的な扱いは得意ですね、楽そうなのでサクッと計算ができそうです
必要に応じて使い分けですかね？
Rでの勉強会はまだ続くのでpythonとの違いをじっくり調べてみたいと思います。

# 追記

### scikit-learnとstatsmodelsの決定係数差について

statsmodelsで回帰を行う場合は、Xに明示的に切片（statsmodelssでｈconstant)を追加する必要があった
デフォルトではなし

statsmodelsの結果
切片なしの決定係数　0.994
切片ありの决定係数　0.962

scikit-learnの結果
決定係数　0.962

Rの結果
決定係数　Multiple R-squared:  0.9617
修正決定係数　Adjusted R-squared:  0.9611


```python

from sklearn import datasets
from statsmodels import api as sm

X_sel_ = sm.add_constant(X_sel)

model2 = sm.OLS(y, X_sel_)
result2 = model2.fit()
print (result2.summary())

                                                         OLS Regression Results                            
==============================================================================
Dep. Variable:         first_purchase   R-squared:                       0.962
Model:                            OLS   Adj. R-squared:                  0.961
Method:                 Least Squares   F-statistic:                     1800.
Date:                Sat, 31 Aug 2019   Prob (F-statistic):          7.74e-202
Time:                        03:58:27   Log-Likelihood:                -1410.7
No. Observations:                 292   AIC:                             2831.
Df Residuals:                     287   BIC:                             2850.
Df Model:                           4                                         
Covariance Type:            nonrobust                                         
==================================================================================
                     coef    std err          t      P>|t|      [0.025      0.975]
----------------------------------------------------------------------------------
const              7.7830     31.498      0.247      0.805     -54.213      69.779
acq_expense        0.5317      0.095      5.577      0.000       0.344       0.719
acq_expense_sq    -0.0007   7.19e-05     -9.099      0.000      -0.001      -0.001
revenue            2.9862      0.109     27.488      0.000       2.772       3.200
employees          0.2483      0.004     65.887      0.000       0.241       0.256
==============================================================================
Omnibus:                        2.045   Durbin-Watson:                   1.889
Prob(Omnibus):                  0.360   Jarque-Bera (JB):                1.739
Skew:                           0.150   Prob(JB):                        0.419
Kurtosis:                       3.230   Cond. No.                     8.58e+06
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 8.58e+06. This might indicate that there are
strong multicollinearity or other numerical problems.


print(result2.aic)

2831.3732200478385

```

###scikit-learnの自由度修正決定係数の扱いについて
scikit-learnでは決定係数のみのようです
自由度調整決定係数は決定係数から計算できるので試してみました
https://www.statisticshowto.datasciencecentral.com/adjusted-r2/
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/258389/260ca9b1-a733-62f6-e128-748bd3ccc215.png)

R,修正後のstatsmodelsとも数値が合いました！

```python

# 調整決系係数
y= df_sel['first_purchase']
X_sel = df_sel.drop('first_purchase',axis=1).columns

adj_r_sqr = 1 - (1-r_sqr)*(len(y)-1)/(len(y)-len(X_sel)-1)

print('adj_r_sqr = :',adj_r_sqr)

adj_r_sqr = : 0.9611309093220055
```


# 追記2
statsmodels.formula.api を用いると、R-styleのモデル式が使用できます

https://www.statsmodels.org/stable/api.html

statsmodels.formula.api: A convenience interface for specifying models using formula strings and DataFrames. This API directly exposes the from_formula class method of models that support the formula API. Canonically imported using import statsmodels.formula.api as smf

```python
import statsmodels.formula.api as smf

model_Rlike = smf.ols('first_purchase ~ acq_expense + acq_expense_sq + revenue + employees', data=df_sel)
result = model_Rlike.fit()
print(result.summary())

                           OLS Regression Results                            
==============================================================================
Dep. Variable:         first_purchase   R-squared:                       0.962
Model:                            OLS   Adj. R-squared:                  0.961
Method:                 Least Squares   F-statistic:                     1800.
Date:                Sun, 03 May 2020   Prob (F-statistic):          7.74e-202
Time:                        10:19:15   Log-Likelihood:                -1410.7
No. Observations:                 292   AIC:                             2831.
Df Residuals:                     287   BIC:                             2850.
Df Model:                           4                                         
Covariance Type:            nonrobust                                         
==================================================================================
                     coef    std err          t      P>|t|      [0.025      0.975]
----------------------------------------------------------------------------------
Intercept          7.7830     31.498      0.247      0.805     -54.213      69.779
acq_expense        0.5317      0.095      5.577      0.000       0.344       0.719
acq_expense_sq    -0.0007   7.19e-05     -9.099      0.000      -0.001      -0.001
revenue            2.9862      0.109     27.488      0.000       2.772       3.200
employees          0.2483      0.004     65.887      0.000       0.241       0.256
==============================================================================
Omnibus:                        2.045   Durbin-Watson:                   1.889
Prob(Omnibus):                  0.360   Jarque-Bera (JB):                1.739
Skew:                           0.150   Prob(JB):                        0.419
Kurtosis:                       3.230   Cond. No.                     8.58e+06
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 8.58e+06. This might indicate that there are
strong multicollinearity or other numerical problems.
'''