書記が数学やるだけ#463 重回帰分析の計算

# ライブラリの読み込み
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import scipy as sp
import statsmodels
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.stats.outliers_influence import variance_inflation_factor

#データセットの読み込み（空白に_を追加）
wine = pd.read_csv("winequality-red.csv", sep=";", header=0, 
                  names=["fixed_acidity", "volatile_acidity", "citric_acid", "residual_sugar", "chlorides", 
                        "free_sulfur_dioxide", "total_sulfur_dioxide", "density", "pH", "sulphates", "alcohol", 
                        "quality"])
wine.head()

#Qualityを応答変数，sulphatesとalcoholを説明変数
wine2 = wine[["sulphates", "alcohol", "quality"]]
wine2

#Qualityを応答変数，sulphatesとalcoholを説明変数
model2 = smf.ols('quality ~ sulphates+alcohol', data=wine2).fit()
model2.summary()

#回帰式の導出
y=wine2["quality"]
x1=wine2["sulphates"]
x2=wine2["alcohol"]
n=len(wine2)

y_sum = y.sum()
y_sum2 = (y**2).sum()
y_bar = y_sum/n
s_yy=y_sum2-((y_sum)**2/n)

x1_sum = x1.sum()
x1_sum2 = (x1**2).sum()
x1_bar = x1_sum/n
s_11=x1_sum2-((x1_sum)**2)/n

x2_sum = x2.sum()
x2_sum2 = (x2**2).sum()
x2_bar = x2_sum/n
s_22=x2_sum2-((x2_sum)**2)/n

x1y_sum=(x1*y).sum()
s_1y=x1y_sum-(x1_sum*y_sum)/n
x2y_sum=(x2*y).sum()
s_2y=x2y_sum-(x2_sum*y_sum)/n
x1x2_sum=(x1*x2).sum()
s_12=x1x2_sum-(x1_sum*x2_sum)/n

delta=1/(s_11*s_22-(s_12)**2)
beta1=delta*(s_22*s_1y - s_12*s_2y)
beta2=delta*(-s_12*s_1y + s_11*s_2y)
beta0=y_bar-(beta1*x1_bar)-(beta2*x2_bar)

print("回帰式：y=", beta0, "+", beta1, "x1+", beta2, "x2")

#寄与率，自由度調整済寄与率の導出
s_R=(beta1*s_1y)+(beta2*s_2y)
s_e=s_yy-s_R
V_e=s_e/(n-3)

R2= s_R/s_yy
R2a= 1-(s_e/(n-3))/(s_yy/(n-1))

print("寄与率：", R2, ", 自由度調整済寄与率：", R2a)

#予測値，標準化残差，テコ比の導出
wine2["y_hat"]=beta0+beta1*x1+beta2*x2
wine2["e"]=y-wine2["y_hat"]
wine2["es"]=wine2["e"]/np.sqrt(V_e)

delta=1/(s_11*s_22-(s_12)**2)
S_11=delta*s_22
S_12=delta*(-s_12)
S_22=delta*s_11
D2=(n-1)* (((x1-x1_bar)**2)*S_11 +2*(x1-x1_bar)*(x2-x2_bar)*S_12 + ((x2-x2_bar)**2)*S_22)

wine2["h"]=1/n + D2/(n-1)
wine2

wine2.plot.scatter(x='sulphates', y="es")

wine2.plot.scatter(x='alcohol', y="es")

wine2.plot.scatter(x="es", y="h")

#信頼区間の導出
x10=np.linspace(0, 2, 100)
x20=np.linspace(5, 15, 100)
y0=beta0+beta1*x10+beta2*x20
D20=(n-1)* (((x10-x1_bar)**2)*S_11 +2*(x10-x1_bar)*(x20-x2_bar)*S_12 + ((x20-x2_bar)**2)*S_22)

from scipy.stats import t
t = t.ppf(0.975, n-3)

y_CI_pm=t*np.sqrt( (1/n + D20/(n-1)) *V_e )
y_CI_m=y0-y_CI_pm
y_CI_p=y0+y_CI_pm

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
fig, ax = plt.subplots(figsize=(6, 6), subplot_kw={'projection' : '3d'})
ax.plot(x10, x20, y0, c='tomato')
ax.plot(x10, x20, y_CI_m, c='orange')
ax.plot(x10, x20, y_CI_p, c='orange')
plt.show()

#予測区間の導出
y_PI_pm=t*np.sqrt(1+ (1/n + D20/(n-1)) *V_e )
y_PI_m=y0-y_PI_pm
y_PI_p=y0+y_PI_pm

fig, ax = plt.subplots(figsize=(6, 6), subplot_kw={'projection' : '3d'})
ax.plot(x10, x20, y0, c='tomato')
ax.plot(x10, x20, y_PI_m, c='orange')
ax.plot(x10, x20, y_PI_p, c='orange')
plt.show()