做数据分析为何要学统计学（10）——如何进行回归分析-编程知识

做数据分析为何要学统计学（10）——如何进行回归分析

回归分析（regression analysis)是量化两种或两种以上因素/变量间相互依赖关系的统计分析方法。回归分析根据因素的数量，分为一元回归和多元回归分析；按因素之间依赖关系的复杂程度，可分为线性回归分析和非线性回归分析。我们通过一下两个例子介绍如何使用python完成回归分析。

在python中有多个软件包可以用于回归分析，在这里我们选择 sklearn软件包中的LinearRegression训练算法，之所以选择该算法是因为它支持多元回归，还可以用于非线性回归分析（多项式回归）。

1.线性回归分析

某调查公司采集了多人健康数据，试图建立体重与身高和每天运动时长的量化关系。

人员	体重	身高	运动时长
1	52	1.65	2
2	63	1.68	2
3	71	1.75	1
4	82	1.78	1
5	90	1.92	2
6	108	1.8	1
7	85	1.78	2
8	72	1.75	1
9	63	1.62	1

代码如下：

import numpy as np
from sklearn.linear_model import LinearRegression
#构造样本数据之因变量（体重）
y=np.array([52,63,71,82,90,108,85,72,63])
#构造样本数据之因变量（身高和运动时长）
X=np.array([[1.65,2],[1.68,2],[1.75,1],[1.78,1],[1.92,2],[1.8,1],[1.78,2],[1.75,1],[1.62,1]])
#创建回归模型
model=LinearRegression().fit(X,y)
#查看模型，第一个是截距，后面是自变量的系数（身高变量和运动时长变量）
model.intercept_,model.coef_
#(-173.9154414624661, array([150.83625049,  -9.33963438]))

公式为： $weight=150.83*height-9.34*time-173.92$ 。拟合值及其残差如下表所示。

人员	体重	身高	运动时长	拟合值	残差
1	52	1.65	2	56.2695	-4.2695
2	63	1.68	2	60.7944	2.2056
3	71	1.75	1	80.6925	-9.6925
4	82	1.78	1	85.2174	-3.2174
5	90	1.92	2	96.9936	-6.9936
6	108	1.8	1	88.234	19.766
7	85	1.78	2	75.8774	9.1226
8	72	1.75	1	80.6925	-8.6925
9	63	1.62	1	61.0846	1.9154

如果残差均值接近0且服务正态分布，也就是说残差是白噪声，则模型通过质量评价。以下为模型评价代码。

#可决系数，自变量对因变量变化的影响程度，越接近1越好，但无经验阈值
r2 = model.score(X, y)
print(r2)
#求预测值
y1=model.predict(X)
#查看模型准确率
print(1-(abs(y1-y)/(y+0.00001)).mean())
#计算残差，残差是白噪声（均值为0的正态分布）系列说明拟合公式已经提取完全部有用信息
res=y-model.predict(X)
#绘制QQ图
from statsmodels.graphics.api import qqplot
%matplotlib inline
ax=qqplot(res,line="s")
#或使用Ljung-Box检验，p值>0.05可认定为白噪声
from statsmodels.stats.diagnostic import acorr_ljungbox
print(acorr_ljungbox(res))

输出结果为：

0.6816823621107787
0.9097391886649957lb_stat  lb_pvalue
1  0.022422   0.880971

上述模型的可决系数 $R^2=0.6816823621107787$ , 准确率为0.9097，Ljung-Box检验 p值=0.880971。模型通过评价。其残差QQ图如下（没有按对角线分布，残差质量不是很高，说明该模型因素间不是质量很高的线性回归关系，因此拟合出来的线性回归模型准确率也不算高）。

2.非线性回归分析

我们对上例进行多项式回归分析，试图用多项式（单项式的线性组合）来拟合其非线性关系。所谓单项式就是数字与字母的乘积，如 $a,5a,ab,ab^2$ 。

代码如下：

#引入高阶单项式构造函数
from sklearn.preprocessing import  PolynomialFeatures
#构造2阶单项式
new_X=PolynomialFeatures(degree=2).fit_transform(X)

构造出的新的单项式为

人员	常数项	$height$	$time$	$height^2$	$height*time$	$time^2$
1	1	1.65	2	2.7225	3.3	4
2	1	1.68	2	2.8224	3.36	4
3	1	1.75	1	3.0625	1.75	1
4	1	1.78	1	3.1684	1.78	1
5	1	1.92	2	3.6864	3.84	4
6	1	1.8	1	3.24	1.8	1
7	1	1.78	2	3.1684	3.56	4
8	1	1.75	1	3.0625	1.75	1
9	1	1.62	1	2.6244	1.62	1

然后以上述数据为自变量，与y建立线性回归关系，就构造出现二元二次多项式回归公式。代码如下。

model=LinearRegression().fit(new_X,y)
model.intercept_,model.coef_
#(-706.4617834336859,array([   0.        ,  739.84211301,    2.80313388, -157.23758678,-20.80508628,    8.40940164]))
r2 = model.score(new_X, y)

公式为 $weight=739.84*height+2.8*time-157.24*height^2-20.81*height*time+8.41*time^2-706.46$