LinearRegressionを使い重回帰分析でt値を出したいのですが、一項目だけ正しくない値になってしまう。

Question

前提

djangoでゴルフのwebアプリを作っています。
その中で重回帰分析をしているのですが、一つの項目だけやたら大きい値になってしまいます。
エクセルやスプレッドシートでも重回帰分析をしてみたのですが、そちらでは正しい値になります。
Statモデルにスコア（目的変数）とパット数など（説明変数）があり、データベースに登録されたそれらの情報をもとに重回帰分析で影響度を計算しようとしています。

該当のソースコード

df = pd.DataFrame(Stat.objects.filter(player_id=pk).values())
        df.columns = ["id", "player_id", "date", "total_score", "ob", "penalty", "fw", "par_on", "putt", "stat_number"]
        x = df.drop(['id','player_id','date','total_score','stat_number'], axis=1)
        y = df['total_score']

        reg = LinearRegression()
        results = reg.fit(x,y)
        coef = reg.coef_.round(4)
        n = x.shape[0]
        p = x.shape[1]

        y_hat = reg.predict(x)
        sse = np.sum((y - y_hat) **2, axis=0)
        sse = sse / (n - p - 1)
        s = np.linalg.inv(np.dot(x.T, x))
        std_err = np.sqrt(np.diagonal(sse * s)).round(4)

        t_values = (coef/ std_err).round(4)

total_score ob penalty fw par_on putt bunker
83	　　　       1	1	 57	55	 34	 1
85	         1	3	 85	55	 38	 0
81	         0	0	 64	50	 36	 0
90	         1	1	 53	27	 37	 1
86	         1	1	 42	27	 34	 1
86	         0	0	 57	33	 31	 0
85	         2	2	 78	55	 36	 1
86	         3	0	 50	22	 31	 0
87	         0	0  	 71	27 	 32	 1
86	         0	2	 42	22	 33	 0
82	         0	2	 57	66	 37	 1
93	         2	1	 28	11	 33	 0
88	         1	0	 50	16	 31	 0
90	         2	0	 46	27	 33	 1
81	         0	1	 50	27	 31	 0
80	         0	0	 57	38	 29	 0
85	         1	0	 85	33	 31	 0
85	         0	1	 71	22	 33	 0
89	         0	0	 57	33	 37	 0
83	         1	0	 64	50	 32	 0

#スプレッドシートでの計算結果
ob         2.0533
penalty  -0.822
fw         0.1249
par_on   -4.1637
putt       3.1677
bunker   0.7899

#現在作っているアプリ上での計算結果
ob         2.1573   -
penalty  -0.9684
fw         0.1144
par_on   -4.7625
putt       12.6963
bunker   0.8931

#Djangoを使用しない、google colab上での計算結果
ob         2.1573   -
penalty  -0.9684
fw         0.1144
par_on   -4.7625
putt       12.6963
bunker   0.8931
(アプリ上と完全に同じになりました)

google colabではcsvファイルを読み込む方法と、csvと同じデータを手書きでリストを作る方法の両方を試しましたが、どちらも同じ値になりました。

puttの値が大きくなり過ぎてしまいます。他の項目("total_score", "ob", "penalty", "fw")は全て正しく計算できています。
計算部分のみ載せています。

試したこと

puttの値だけ上手く変数に入れられていないのではと思い、xをテンプレートに渡し値を確認しましたが問題ありませんでした。
同じくcoefとstd_errの値を確認しましたが、どちらともおかしい値でした（エクセルなどで確認した値と比較して）。

補足情報（FW/ツールのバージョンなど）

django3
python3
ここにより詳細な情報を記載してください。

Comments

[退会済みユーザー]

何がどうと言う前に問題点をちゃんと「自分で」切り分けましょう。

[8524ba23]

django処理はおいといて
- dfデータ内容とPython上でのT値
- Excelで同データを検証したときの（できれば手順と）T値などの結果。
といった具体的なものを提示すると回答得られやすくなるかと思います。
とくにdfデータが数十行程度であればCSV形式のテキストとして質問に追記なりすると
第三者が実行、検証しやすくなります。

[meg_]

・django上でDBから値を取得する処理に問題はありませんか？
・djangoを介さずにDBとPythonで重回帰分析で実行した場合は正しい結果が得られますか？

[yoshiringoo]

dameo様　ありがとうございます。

can110様　ご返信ありがとうございます。質問するにあたり情報不足で大変失礼いたしました！csvデータを追記しましたので、もしお時間ございましたら確認頂ければ大変助かります。よろしくお願いします。

meg_様　ご返信ありがとうございます。
views.pyのxの値をテンプレートに渡して確認してみたのですが、csvファイル通りの値がしっかり入っていました。
シンプルにpythonでまず試してみるという基本的なことを失念しておりました・・・ありがとうございます。本文に追記させていただきましたが、それでも間違った値になってしまいました。

Answer 1

回帰分析については詳しくありませんが、t値の求め方元のデータが標準化されていないことが原因だと思われます。
提示コードでは以下などを参考に係数(coef)のみから求めていると思います。

• 回帰分析のt値の求め方:Pythonで実装
• [Docker, Python] (重回帰分析編, sk-learn, t値算出, 忘備録)

上記の場合、事前にデータを標準化することで切片項がなくなって正しい結果が得られるようです。
参考：【重回帰分析】標準化と切片の扱いについて (表現の違いについて)

いっぽう以下では切片(intercept)も含めて求めるようなコードになっています。

• Find p-value (significance) in scikit-learn LinearRegression
• 回帰係数のp値とt値が不明なときに、p値とt値を求める方法

この方法であれば元のデータのままでもstatsmodels.regression.linear_model.OLSでの結果と一致します。

import numpy as np
import pandas as pd
from io import StringIO
from sklearn.linear_model import LinearRegression
import statsmodels.api as sm

# データ
s = """total_score,ob,penalty,fw,par_on,putt,bunker
83,1,1,57,55,34,1
85,1,3,85,55,38,0
81,0,0,64,50,36,0
90,1,1,53,27,37,1
86,1,1,42,27,34,1
86,0,0,57,33,31,0
85,2,2,78,55,36,1
86,3,0,50,22,31,0
87,0,0,71,27,32,1
86,0,2,42,22,33,0
82,0,2,57,66,37,1
93,2,1,28,11,33,0
88,1,0,50,16,31,0
90,2,0,46,27,33,1
81,0,1,50,27,31,0
80,0,0,57,38,29,0
85,1,0,85,33,31,0
85,0,1,71,22,33,0
89,0,0,57,33,37,0
83,1,0,64,50,32,0"""
df = pd.read_csv(StringIO(s))
x = df.drop(['total_score'], axis=1)
y = df['total_score']


#-------------------------------------------------------------------------------
# 質問のコード

reg = LinearRegression()
results = reg.fit(x,y)
coef = reg.coef_.round(4)
n = x.shape[0]
p = x.shape[1]

y_hat = reg.predict(x)
sse = np.sum((y - y_hat) **2, axis=0)
sse = sse / (n - p - 1)
s = np.linalg.inv(np.dot(x.T, x))
std_err = np.sqrt(np.diagonal(sse * s)).round(4)

t_values = (coef/ std_err).round(4)
print(t_values)
# [ 2.1573 -0.9684  0.1144 -4.7625 12.6963  0.8931]


# 以下とコードを合わせるために
X = x
params = np.append(reg.intercept_,reg.coef_)
predictions = y_hat


#-------------------------------------------------------------------------------
# Find p-value (significance) in scikit-learn LinearRegression
# https://stackoverflow.com/a/42677750

newX = pd.DataFrame({"Constant":np.ones(len(X))}).join(pd.DataFrame(X))
MSE = (sum((y-predictions)**2))/(len(newX)-len(newX.columns))

# Note if you don't want to use a DataFrame replace the two lines above with
# newX = np.append(np.ones((len(X),1)), X, axis=1)
# MSE = (sum((y-predictions)**2))/(len(newX)-len(newX[0]))

var_b = MSE*(np.linalg.inv(np.dot(newX.T,newX)).diagonal())
sd_b = np.sqrt(var_b)
ts_b = params/ sd_b
print(ts_b)
# [ 8.48895225  2.14428475 -0.85496786  0.11136259 -4.60776964  3.33811352 0.87347851]


#-------------------------------------------------------------------------------

X2 = sm.add_constant(X)
est = sm.OLS(y, X2)
est2 = est.fit()
print(est2.summary())
"""
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         64.4012      7.586      8.489      0.000      48.012      80.791
ob             1.1246      0.524      2.144      0.051      -0.008       2.258
penalty       -0.5263      0.616     -0.855      0.408      -1.856       0.804
fw             0.0042      0.037      0.111      0.913      -0.077       0.085
par_on        -0.1905      0.041     -4.608      0.000      -0.280      -0.101
putt           0.7986      0.239      3.338      0.005       0.282       1.315
bunker         0.9099      1.042      0.873      0.398      -1.340       3.160
"""

Comments

[yoshiringoo]

can110様　こんなご丁寧にありがとうございます！こちらの方法で正しいt値を求めることができました。
また本日回帰分析について調べていたところ、私はデータの標準化を飛ばしていてしまったようでした。標準化をしましたら、正しい値になりましたので下記に載せておきます。
df.apply(lambda x: (x-x.mean())/x.std(), axis=0)
とすることで、dfを行毎に標準化できました。
ありがとうございました。

[8524ba23]

なるほど。切片がなくなるので正しい結果が得られるようですね。
勉強になります。回答にも追記しておきました。