存在哑变量的多重线性回归

转载

信息小飞侠 2024-12-01 22:22:03

文章标签 存在哑变量的多重线性回归线性回归 python 机器学习 Marketing 文章分类 机器学习人工智能

线性回归

statsmodels

# 导入第三方模块
import statsmodels.api as sm

ols(formula, data, subset=None, drop_cols=None)

formula：以字符串的形式指定线性回归模型的公式，如'y~x'就表示简单线性回归模型
data：指定建模的数据集
subset：通过bool类型的数组对象，获取data的子集用于建模
drop_cols：指定需要从data中删除的变量

一元应用

# 导入第三方模块
import statsmodels.api as sm
# 利用收入数据集，构建回归模型
fit = sm.formula.ols('Salary ~ YearsExperience', data = income).fit()
# 返回模型的参数值
fit.params

out:
Intercept 25792.200199
YearsExperience 9449.962321
dtype: float64

多元应用

# 导入模块
from sklearn import model_selection
# 导入数据
Profit = pd.read_excel(r'C:\Users\Administrator\Desktop\Predict to Profit.xlsx')
# 将数据集拆分为训练集和测试集
train, test = model_selection.train_test_split(Profit, test_size = 0.2, random_state=1234)
# 根据train数据集建模
model = sm.formula.ols('Profit ~ RD_Spend+Administration+Marketing_Spend+C(State)', data = train).fit()
print('模型的偏回归系数分别为：\n', model.params)
# 删除test数据集中的Profit变量，用剩下的自变量进行预测
test_X = test.drop(labels = 'Profit', axis = 1)
pred = model.predict(exog = test_X)
print('对比预测值和实际值的差异：\n',pd.DataFrame({'Prediction':pred,'Real':test.Profit}))

通过pandas模块中的get_dummies函数生成哑变量，然后将所需的对照组
对应的哑变量删除即可。为了使读者明白该解决方案，这里不妨重新建模，并以State变量中的New York值作为对照组

# 生成由State变量衍生的哑变量
dummies = pd.get_dummies(Profit.State)
# 将哑变量与原始数据集水平合并
Profit_New = pd.concat([Profit,dummies], axis = 1)
# 删除State变量和California变量（因为State变量已被分解为哑变量，New York变量需要作为参照组）
Profit_New.drop(labels = ['State','New York'], axis = 1, inplace = True)
# 拆分数据集Profit_New
train, test = model_selection.train_test_split(Profit_New, test_size = 0.2, random_state=1234)
# 建模
model2 = sm.formula.ols('Profit~RD_Spend+Administration+Marketing_Spend+Florida+California', data = 
train).fit()
print('模型的偏回归系数分别为：\n', model2.params)

检验

F检验

# 导入第三方模块
import numpy as np
# 计算建模数据中因变量的均值
ybar = train.Profit.mean()
# 统计变量个数和观测个数
p = model2.df_model
n = train.shape[0]
# 计算回归离差平方和
RSS = np.sum((model2.fittedvalues-ybar) ** 2)
# 计算误差平方和
ESS = np.sum(model2.resid ** 2)
# 计算F统计量的值
F = (RSS/p)/(ESS/(n-p-1))
print('F统计量的值：',F)
out:
F统计量的值： 174.6372

# 导入模块
from scipy.stats import f
# 计算F分布的理论值
F_Theroy = f.ppf(q=0.95, dfn = p,dfd = n-p-1)
print('F分布的理论值为：',F_Theroy)
out:
F分布的理论值为： 2.5026

# 有关模型的概览信息
model2.summary()

直方图

# 导入第三方模块
import scipy.stats as stats
# 中文和负号的正常显示
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False
# 绘制商品利润的直方图
sns.distplot(a = Profit_New.Profit, bins = 10, fit = stats.norm, norm_hist = True,
hist_kws = {'color':'steelblue', 'edgecolor':'black'}, 
kde_kws = {'color':'black', 'linestyle':'--', 'label':'核密度曲线'}, 
fit_kws = {'color':'red', 'linestyle':':', 'label':'正态密度曲线'})
# 显示图例
plt.legend()
# 显示图形
plt.show()

pp与qq图

# 残差的正态性检验（PP图和QQ图法）
pp_qq_plot = sm.ProbPlot(Profit_New.Profit)
# 绘制PP图
pp_qq_plot.ppplot(line = '45')
plt.title('P-P图')
# 绘制QQ图
pp_qq_plot.qqplot(line = 'q')
plt.title('Q-Q图')
# 显示图形
plt.show()

正态性检验

# 导入模块
import scipy.stats as stats
# Shapiro检验
stats.shapiro(Profit_New.Profit)
out:
(0.9793398380279541, 0.537902295589447)
Shapiro检验和K-S检验均属于非参数方法，它们的原假设被设定为变量服从正态分布，两者的最
大区别在于适用的数据量不一样，若数据量低于5000，则使用Shapiro检验法比较合理，否则使用
K-S检验法。


# 生成正态分布和均匀分布随机数（样本量超过5000的情况）
rnorm = np.random.normal(loc = 5, scale=2, size = 10000)
runif = np.random.uniform(low = 1, high = 100, size = 10000)
# 正态性检验
KS_Test1 = stats.kstest(rvs = rnorm, args = (rnorm.mean(), rnorm.std()), cdf = 'norm')
KS_Test2 = stats.kstest(rvs = runif, args = (runif.mean(), runif.std()), cdf = 'norm')
print(KS_Test1)
print(KS_Test2)
out:
KstestResult(statistic=0.006035, pvalue=0.8597)
KstestResult(statistic=0.061127, pvalue=7.0185e-33)

多重共线性检验–VIF方法

# 导入statsmodels模块中的函数
from statsmodels.stats.outliers_influence import variance_inflation_factor
# 自变量X(包含RD_Spend、Marketing_Spend和常数列1)
X = sm.add_constant(Profit_New.ix[:,['RD_Spend','Marketing_Spend']])
# 构造空的数据框，用于存储VIF值
vif = pd.DataFrame()
vif["features"] = X.columns
vif["VIF Factor"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
# 返回VIF值
vif

异常值检验

# 异常值检验
outliers = model3.get_influence()
# 高杠杆值点（帽子矩阵）
leverage = outliers.hat_matrix_diag
# dffits值
dffits = outliers.dffits[0]
# 学生化残差
resid_stu = outliers.resid_studentized_external
# cook距离
cook = outliers.cooks_distance[0]

# 合并各种异常值检验的统计量值
contat1 = pd.concat([pd.Series(leverage, name = 'leverage'),pd.Series(dffits, name = 
'dffits'),pd.Series(resid_stu,name = 'resid_stu'),pd.Series(cook, name = 'cook')],axis = 1)
# 重设train数据的行索引
train.index = range(train.shape[0])
# 将上面的统计量与train数据集合并
profit_outliers = pd.concat([train,contat1], axis = 1)
profit_outliers.head()

使用标准化残差，当标准化残差大于2时，即认为对应的数据点为异常值。

# 计算异常值数量的比例
outliers_ratio = sum(np.where((np.abs(profit_outliers.resid_stu)>2),1,0))/ profit_outliers.shape[0]
outliers_ratio
out:
0.025

通过标准化残差监控到了异常值，并且异常比例为2.5%。对于异常值的处
理办法，可以使用两种策略，如果异常样本的比例不高（如小于等于5%），可以考虑将异常点删除；如果异常样本的比例比较高，选择删除会丢失一些重要信息，所以需要衍生哑变量，即对于异常点，设置哑变量的值为1，否则为0。如上可知，建模数据的异常比例只有2.5%，故考虑将其删除。

# 利用学生化残差挑选出非异常的观测点（由于异常比例非常小，故选择删除法）
none_outliers = profit_outliers.ix[np.abs(profit_outliers.resid_stu)<=2,]
# 应用无异常值的数据集重新建模
model4 = sm.formula.ols ('Profit ~ RD_Spend + Marketing_Spend', data = none_outliers).fit()
model4.params

out:
Intercept 51827.416821
RD_Spend 0.797038
Marketing_Spend 0.017740

方差齐次性检验

图形法

方差齐性检验--图形法
# 设置第一张子图的位置
ax1 = plt.subplot2grid(shape = (2,1), loc = (0,0))
# 绘制散点图
ax1.scatter(none_outliers.RD_Spend, (model4.residmodel4.resid.mean())/model4.resid.std())
# 添加水平参考线
ax1.hlines(y = 0 ,xmin = none_outliers.RD_Spend.min(),
xmax = none_outliers.RD_Spend.max(), color = 'red', linestyles = '--')
# 添加x轴和y轴标签
ax1.set_xlabel('RD_Spend')
ax1.set_ylabel('Std_Residual')
# 设置第二张子图的位置
ax2 = plt.subplot2grid(shape = (2,1), loc = (1,0))

# 绘制散点图
ax2.scatter(none_outliers.Marketing_Spend, 
(model4.resid-model4.resid.mean())/model4.resid.std())
# 添加水平参考线
ax2.hlines(y = 0 ,xmin = none_outliers.Marketing_Spend.min(),
xmax = none_outliers.Marketing_Spend.max(), color = 'red', linestyles
= '--')
# 添加x轴和y轴标签
ax2.set_xlabel('Marketing_Spend')
ax2.set_ylabel('Std_Residual')
# 调整子图之间的水平间距和高度间距
plt.subplots_adjust(hspace=0.6, wspace=0.3)
# 显示图形

统计法

# BP检验
sm.stats.diagnostic.het_breushpagan(model4.resid, exog_het = model4.model.exog)
out:
(1.4675103668307794,
0.48010272699007694,
0.70297512371621873,
0.50196597409630139)
如上结果所示，元组中一共包含四个值，第一个值1.468为𝐿𝐿𝐿𝐿统计量；第二个值是统计量对
应的概率𝑝𝑝值，该值大于0.05，说明接受残差方差为常数的原假设；第三个值为𝐹𝐹统计量，用于
检验残差平方项与自变量之间是否独立，如果独立则表明残差方差齐性；第四个值则为𝐹𝐹统计量
的概率𝑝𝑝值，同样大于0.05，则进一步表示残差项满足方差齐性的假设。

本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。