回归分析

11、回归分析

1. 回归分析概述

相关分析是研究两个或两个以上的变量之间相关程度及大小的一种统计方法

回归分析是寻找存在相关关系的变量间的数学表达式，并进行统计推断的一种统计方法

在对回归分析进行分类时，主要有两种分类方式：

根据变量的数目，可以分类一元回归、多元回归
根据自变量与因变量的表现形式，分为线性与非线性

所以，回归分析包括四个方向：一元线性回归分析、多元线性回归分析、一元非线性回归分析、多元非线性回归分析。

2. 回归方程定义

一元线性回归分析

因变量(dependent variable)：被预测或被解释的变量，用y表示
自变量(independent variable)：预测或解释因变量的一个或多个变量，用x表示
对于具有线性关系的两个变量，可以用一个方程来表示它们之间的线性关系

描述因变量y如何依赖于自变量x和误差项ε的方程称为回归模型。对于只涉及一个自变量的一元线性回归模型可表示为：

y = β0 + β1*x + ε
- y叫做因变量或被解释变量
- x叫做自变量或解释变量
- β0 表示截距
- β1 表示斜率
- ε表示误差项，反映除x和y之间的线性关系之外的随机因素对y的影响

3. 误差项的定义

4. 最小二乘法的推导和求解

如果回归方程中的参数已知，对于一个给定的x值，利用回归方程就能计算出y的期望值
用样本统计量代替回归方程中的未知参数，就得到估计的回归方程，简称回归直线
对于回归直线，关键在于求解参数，常用高斯提出的最小二乘法，它是使因变量的观察值y与估计值之间的离差平方和达到最小来求解

5. 回归方程求解最小例子

点估计：利用估计的回归方程，对于x的某一个特定的值，求出y的一个估计值就是点估计
区间估计：利用估计的回归方程，对于x的一个特定值，求出y的一个估计值的区间就是区间估计
影响区间宽度的因素：
- 置信水平 (1 - α)，区间宽度随置信水平的增大而增大
- 数据的离散程度Se，区间宽度随离程度的增大而增大
- 样本容量，区间宽度随样本容量的增大而减小
- X0与X均值之间的差异，随着差异程度的增大而增大

6. 回归直线拟合优度

回归直线与各观测点的接近程度称为回归直线对数据的拟合优度

总平方和可以分解为回归平方和、残差平方和两部分：SST＝SSR+SSE

总平方和(SST)，反映因变量的 n 个观察值与其均值的总离差
回归平方和SSR反映了y的总变差中，由于x与y之间的线性关系引起的y的变化部分
残差平方和SSE反映了除了x对y的线性影响之外的其他因素对y变差的作用，是不能由回归直线来解释的y的变差部分

判定系数

回归平方和占总平方和的比例，用R^2表示,其值在0到1之间。

R^2 == 0：说明y的变化与x无关，x完全无助于解释y的变差
R^2 == 1：说明残差平方和为0，拟合是完全的，y的变化只与x有关

显著性检验

显著性检验的主要目的是根据所建立的估计方程用自变量x来估计或预测因变量y的取值。当建立了估计方程后，还不能马上进行估计或预测，因为该估计方程是根据样本数据得到的，它是否真实的反映了变量x和y之间的关系，则需要通过检验后才能证实。

根据样本数据拟合回归方程时，实际上就已经假定变量x与y之间存在着线性关系，并假定误差项是一个服从正态分布的随机变量，且具有相同的方差。但这些假设是否成立需要检验

显著性检验包括两方面：

线性关系检验
回归系数检验

线性关系检验
线性关系检验是检验自变量x和因变量y之间的线性关系是否显著，或者说，它们之间能否用一个线性模型来表示。

将均方回归 (MSR)同均方残差 (MSE)加以比较，应用F检验来分析二者之间的差别是否显著。

均方回归：回归平方和SSR除以相应的自由度(自变量的个数K)
均方残差：残差平方和SSE除以相应的自由度(n-k-1)
H0：β1=0 所有回归系数与零无显著差异，y与全体x的线性关系不显著
计算检验统计量F：

     回归系数显著性检验的目的是通过检验回归系数β的值与0是否有显著性差异，来判断Y与X之间是否有显著的线性关系.若β=0,则总体回归方程中不含X项(即Y不随X变动而变动),因此,变量Y与X之间并不存在线性关系;若β≠0,说明变量Y与X之间存在显著的线性关系。

线性关系的检验是检验自变量与因变量是否可以用线性来表达，而回归系数的检验是对样本数据计算的回归系数检验总体中回归系数是否为0

在一元线性回归中，自变量只有一个，线性关系检验与回归系数检验是等价的
多元回归分析中，这两种检验的意义是不同的。线性关系检验只能用来检验总体回归关系的显著性，而回归系数检验可以对各个回归系数分别进行检验

7. 多元与曲线回归问题

经常会遇到某一现象的发展和变化取决于几个影响因素的情况，也就是一个因变量和几个自变量有依存关系的情况，这时需用多元线性回归分析。

多元线性回归分析预测法，是指通过对两上或两个以上的自变量与一个因变量的相关分析，建立预测模型进行预测和控制的方法
多元线性回归预测模型一般式为：

曲线回归分析

直线关系是两变量间最简单的一种关系，曲线回归分析的基本任务是通过两个相关变量x与y的实际观测数据建立曲线回归方程，以揭示x与y间的曲线联系的形式。

曲线回归分析最困难和首要的工作是确定自变量与因变量间的曲线关系的类型，曲线回归分析的基本过程：

先将x或y进行变量转换
对新变量进行直线回归分析、建立直线回归方程并进行显著性检验和区间估计
将新变量还原为原变量，由新变量的直线回归方程和置信区间得出原变量的曲线回归方程和置信区间

由于曲线回归模型种类繁多，所以没有通用的回归方程可直接使用。但是对于某些特殊的回归模型，可以通过变量代换、取对数等方法将其线性化，然后使用标准方程求解参数，再将参数带回原方程就是所求。

多重共线性

回归模型中两个或两个以上的自变量彼此相关的现象

多重共线性带来的问题有:

回归系数估计值的不稳定性增强
回归系数假设检验的结果不显著等

多重共线性检验的主要方法:

容忍度
方差膨胀因子（VIF）

8. python工具包

9. statsmodels回归分析

python

import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm

nsample = 20
x = np.linspace(0, 10, nsample)
x

# 一元线性回归
# 它可以将一个n行、k列的矩阵X变为一个n行、k+1列的矩阵，其中第一列都是1，代表常数项
X = sm.add_constant(x)
#β0,β1分别设置成2,5
beta = np.array([2, 5])
#误差项
e = np.random.normal(size=nsample)
#实际值y
y = np.dot(X, beta) + e

#最小二乘法
model = sm.OLS(y,X)

#拟合数据
res = model.fit()
#回归系数
res.params
array([ 1.49524076,  5.08701837])

#全部结果
res.summary()
# R-squared:	0.995 R^2值  R-squared:	0.995   调整过的R^2值
# F-statistic:	3668.  F检验 Prob (F-statistic):	2.94e-22 F检验的概率值 比较小相关性强
# const 常数项 x1斜率项

python

#拟合的估计值
y_ = res.fittedvalues
y_
array([  1.49524076,   4.17261885,   6.84999693,   9.52737502,
        12.20475311,  14.8821312 ,  17.55950928,  20.23688737,
        22.91426546,  25.59164354,  28.26902163,  30.94639972,
        33.62377781,  36.30115589,  38.97853398,  41.65591207,
        44.33329015,  47.01066824,  49.68804633,  52.36542442])

python

fig, ax = plt.subplots(figsize=(8,6))
ax.plot(x, y, 'o', label='data')#原始数据
ax.plot(x, y_, 'r--.',label='test')#拟合数据
ax.legend(loc='best')
plt.show()

10. 高阶与分类变量实例

高阶回归

PYTHON

#Y=5+2⋅X+3⋅X^2
 nsample = 50
x = np.linspace(0, 10, nsample)
X = np.column_stack((x, x**2))
X = sm.add_constant(X)
array([[1.00000000e+00, 0.00000000e+00, 0.00000000e+00],

PYTHON

beta = np.array([5, 2, 3])
e = np.random.normal(size=nsample)
y = np.dot(X, beta) + e
model = sm.OLS(y,X)
results = model.fit()
results.params
array([ 4.93210623,  2.16604081,  2.97682135])

results.summary()

PYTHON

y_fitted = results.fittedvalues
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(x, y, 'o', label='data')
ax.plot(x, y_fitted, 'r--.',label='OLS')
ax.legend(loc='best')
plt.show()

分类变量

假设分类变量有4个取值（a,b,c）,比如考试成绩有3个等级。那么a就是（1,0,0），b（0,1,0），c(0,0,1),这个时候就需要3个系数β0,β1,β2，也就是β0x0+β1x1+β2x2

python

nsample = 50
groups = np.zeros(nsample,int)

groups[20:40] = 1
groups[40:] = 2
# 用于将离散变量转化为其它可以用于回归模型的变量
# 参数drop可以设置为True或False，表示是否删除一列
dummy = sm.categorical(groups, drop=True)
[1 2 1 2 3 1 2 1 0 3]
[[0. 1. 0. 0.]
 [0. 0. 1. 0.]
 [0. 1. 0. 0.]
 [0. 0. 1. 0.]
 [0. 0. 0. 1.]
 [0. 1. 0. 0.]
 [0. 0. 1. 0.]
 [0. 1. 0. 0.]
 [1. 0. 0. 0.]
 [0. 0. 0. 1.]]

python

#Y=5+2X+3Z1+6⋅Z2+9⋅Z3.
x = np.linspace(0, 20, nsample)
X = np.column_stack((x, dummy))
X = sm.add_constant(X)
beta = [5, 2, 3, 6, 9]
e = np.random.normal(size=nsample)
y = np.dot(X, beta) + e
result = sm.OLS(y,X).fit()
result.summary()

python

fig, ax = plt.subplots(figsize=(8,6))
ax.plot(x, y, 'o', label="data")
ax.plot(x, result.fittedvalues, 'r--.', label="OLS")
ax.legend(loc='best')
plt.show()

11. 案例: 汽车价格预测任务

任务概述

主要包括3类指标:

汽车的各种特性.
保险风险评级：(-3, -2, -1, 0, 1, 2, 3).
每辆保险车辆年平均相对损失支付.

类别属性

make: 汽车的商标（奥迪，宝马。。。）
fuel-type: 汽油还是天然气
aspiration: 涡轮
num-of-doors: 两门还是四门
body-style: 硬顶车、轿车、掀背车、敞篷车
drive-wheels: 驱动轮
engine-location: 发动机位置
engine-type: 发动机类型
num-of-cylinders: 几个气缸
fuel-system: 燃油系统

连续指标

bore: continuous from 2.54 to 3.94.
stroke: continuous from 2.07 to 4.17.
compression-ratio: continuous from 7 to 23.
horsepower: continuous from 48 to 288.
peak-rpm: continuous from 4150 to 6600.
city-mpg: continuous from 13 to 49.
highway-mpg: continuous from 16 to 54.
price: continuous from 5118 to 45400.

python

# loading packages
import numpy as np
import pandas as pd
from pandas import datetime

# data visualization and missing values
import matplotlib.pyplot as plt
import seaborn as sns # advanced vizs
import missingno as msno # missing values
%matplotlib inline

# stats
from statsmodels.distributions.empirical_distribution import ECDF
from sklearn.metrics import mean_squared_error, r2_score

# machine learning
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import Lasso, LassoCV 
from sklearn.model_selection import train_test_split, cross_val_score 
from sklearn.ensemble import RandomForestRegressor
seed = 123

# importing data ( ? = missing values)
data = pd.read_csv("Auto-Data.csv", na_values = '?')
data.columns

python

data.dtypes

python

print("In total: ",data.shape)
data.head(5)
In total:  (205, 26)
    
data.describe()

缺失值填充

python

# missing values?
sns.set(style = "ticks")

msno.matrix(data)
#https://github.com/ResidentMario/missingno
# 缺失值占比比较小 可以去掉
# 用均值去填充
# 用线性回归去预测缺失值

python

# 查看缺失值比较多的
data[pd.isnull(data['normalized-losses'])].head()

sns.set(style = "ticks")
plt.figure(figsize = (12, 5)) 
c = '#366DE8'

# ECDF 指经验累积分布函数
# 它是在统计学和概率论中，用来描述样本的累积分布情况的一种非参数方法。具体而言，ECDF是对样本的经验分布函数进行逐点估计，它表示所有数值小于某个特定值的样本比例。
plt.subplot(121)
cdf = ECDF(data['normalized-losses'])
plt.plot(cdf.x, cdf.y, label = "statmodels", color = c);
plt.xlabel('normalized losses'); plt.ylabel('ECDF');

# overall distribution
plt.subplot(122)
plt.hist(data['normalized-losses'].dropna(), 
         bins = int(np.sqrt(len(data['normalized-losses']))),
         color = c);

可以发现 80% 的 normalized losses 是低于200 并且绝大多数低于125. 

一个基本的想法就是用中位数来进行填充，但是我们得来想一想，这个特征跟哪些因素可能有关呢？应该是保险的情况吧，所以我们可以分组来进行填充这样会更精确一些。
首先来看一下对于不同保险情况的统计指标:

data.groupby('symboling')['normalized-losses'].describe()

python

# replacing
# 直接删掉缺失值所在的行
data = data.dropna(subset = ['price', 'bore', 'stroke', 'peak-rpm', 'horsepower', 'num-of-doors'])
# 根据分组来填充均值
data['normalized-losses'] = data.groupby('symboling')['normalized-losses'].transform(lambda x: x.fillna(x.mean()))

print('In total:', data.shape)
data.head()
In total: (193, 26)

特征相关性

python

cormatrix = data.corr()
cormatrix *= np.tri(*cormatrix.values.shape, k=-1).T  #返回函数的上三角矩阵，把对角线上的置0，让他们不是最高的。
cormatrix = cormatrix.stack() # 当前指标和其他指标的关系
cormatrix = cormatrix.reindex(cormatrix.abs().sort_values(ascending=False).index).reset_index() # 两两之间的关系
cormatrix.columns = ["FirstVariable", "SecondVariable", "Correlation"]

cormatrix
	FirstVariable	SecondVariable	Correlation
0	city-mpg	highway-mpg	0.971975
1	engine-size	price	0.888778

python

# city_mpg 和 highway-mpg 这哥俩差不多是一个意思了. 对于这个长宽高，他们应该存在某种配对关系，给他们合体吧！
# 相关性差不多的  合体  
data['volume'] = data.length * data.width * data.height
# 删掉不要的列
data.drop(['width', 'length', 'height', 
           'curb-weight', 'city-mpg'], 
          axis = 1, # 1 for columns
          inplace = True) 
# new variables
data.columns

python

# Compute the correlation matrix 
corr_all = data.corr()

# Generate a mask for the upper triangle
mask = np.zeros_like(corr_all, dtype = np.bool)
mask[np.triu_indices_from(mask)] = True

# Set up the matplotlib figure
f, ax = plt.subplots(figsize = (11, 9))

# Draw the heatmap with the mask and correct aspect ratio
sns.heatmap(corr_all, mask = mask,
            square = True, linewidths = .5, ax = ax, cmap = "BuPu")      
plt.show()
# 看起来 price 跟这几个的相关程度比较大 wheel-base,enginine-size, bore,horsepower.

python

sns.pairplot(data, hue = 'fuel-type', palette = 'plasma') # 变量两两之间关系

python

# 让我们仔细看看价格和马力变量之间的关系
# 在fuel-type和num-of-doors条件下，price和horsepower的关系
sns.lmplot('price', 'horsepower', data, 
           hue = 'fuel-type', col = 'fuel-type',  row = 'num-of-doors', 
           palette = 'plasma', 
           fit_reg = True);
# 事实上，对于燃料的类型和数门变量，我们看到，在一辆汽车马力的增加与价格成比例的增加相关的各个层面

预处理问题

如果一个特性的方差比其他的要大得多，那么它可能支配目标函数，使估计者不能像预期的那样正确地从其他特性中学习。这就是为什么我们需要首先对数据进行缩放。

对连续值进行标准化

python

# target and features
target = data.price

regressors = [x for x in data.columns if x not in ['price']]
features = data.loc[:, regressors]
# 只对数字型的进行标准化 其他除外
num = ['symboling', 'normalized-losses', 'volume', 'horsepower', 'wheel-base',
       'bore', 'stroke','compression-ratio', 'peak-rpm']

# scale the data
standard_scaler = StandardScaler() # 标准化
features[num] = standard_scaler.fit_transform(features[num])

# glimpse
features.head()

python

# 对分类属性就行one-hot编码
# categorical vars 对字符串类型的 进行热编码 
classes = ['make', 'fuel-type', 'aspiration', 'num-of-doors', 
           'body-style', 'drive-wheels', 'engine-location',
           'engine-type', 'num-of-cylinders', 'fuel-system']

# create new dataset with only continios vars 
dummies = pd.get_dummies(features[classes])
features = features.join(dummies).drop(classes, 
                                       axis = 1)

# new dataset
print('In total:', features.shape)
features.head()

python

# 划分数据集
# split the data into train/test set
X_train, X_test, y_train, y_test = train_test_split(features, target, 
                                                    test_size = 0.3,
                                                    random_state = seed)
print("Train", X_train.shape, "and test", X_test.shape)
Train (135, 66) and test (58, 66)

回归问题

python

# Lasso回归
# 多加了一个绝对值项来惩罚过大的系数，alphas=0那就是最小二乘了

# logarithmic scale: log base 2
# high values to zero-out more variables
alphas = 2. ** np.arange(2, 12)
scores = np.empty_like(alphas)

for i, a in enumerate(alphas):
    lasso = Lasso(random_state = seed)
    lasso.set_params(alpha = a)
    lasso.fit(X_train, y_train)
    scores[i] = lasso.score(X_test, y_test)

# 交叉验证 平均切成10份
lassocv = LassoCV(cv = 10, random_state = seed)
lassocv.fit(features, target)
lassocv_score = lassocv.score(features, target)
lassocv_alpha = lassocv.alpha_

plt.figure(figsize = (10, 4))
plt.plot(alphas, scores, '-ko')
plt.axhline(lassocv_score, color = c)
plt.xlabel(r'$\alpha$')
plt.ylabel('CV Score')
plt.xscale('log', basex = 2)
sns.despine(offset = 15)
# 选择一个最高点 为惩罚系数
print('CV results:', lassocv_score, lassocv_alpha)

python

# lassocv coefficients
coefs = pd.Series(lassocv.coef_, index = features.columns)

# prints out the number of picked/eliminated features
print("Lasso picked " + str(sum(coefs != 0)) + " features and eliminated the other " +  \
      str(sum(coefs == 0)) + " features.")

# takes first and last 10
coefs = pd.concat([coefs.sort_values().head(5), coefs.sort_values().tail(5)])

plt.figure(figsize = (10, 4))
coefs.plot(kind = "barh", color = c)
plt.title("Coefficients in the Lasso Model")
plt.show()

python

# 用得到的最好的alpha去建模
model_l1 = LassoCV(alphas = alphas, cv = 10, random_state = seed).fit(X_train, y_train)
y_pred_l1 = model_l1.predict(X_test)

model_l1.score(X_test, y_test)
0.83307445226244159

python

# residual plot 残差点画出来  实际值和预测值之间的差距
# 期望在0左右浮动
plt.rcParams['figure.figsize'] = (6.0, 6.0)

preds = pd.DataFrame({"preds": model_l1.predict(X_train), "true": y_train})
preds["residuals"] = preds["true"] - preds["preds"]
preds.plot(x = "preds", y = "residuals", kind = "scatter", color = c)

python

def MSE(y_true,y_pred):
    mse = mean_squared_error(y_true, y_pred)
    print('MSE: %2.3f' % mse)
    return mse

def R2(y_true,y_pred):    
    r2 = r2_score(y_true, y_pred)
    print('R2: %2.3f' % r2)     
    return r2
# 均方误差MSE    r^2  还可以 
MSE(y_test, y_pred_l1); R2(y_test, y_pred_l1);

MSE: 3870543.789
R2: 0.833

python

# predictions
d = {'true' : list(y_test),
     'predicted' : pd.Series(y_pred_l1)
    }

pd.DataFrame(d).head()

回归分析 ​

11、 回归分析 ​

1. 回归分析概述 ​

2. 回归方程定义 ​

3. 误差项的定义 ​

4. 最小二乘法的推导和求解 ​

5. 回归方程求解最小例子 ​

6. 回归直线拟合优度 ​

判定系数 ​

显著性检验 ​

7. 多元与曲线回归问题 ​

曲线回归分析 ​

多重共线性 ​

8. python工具包 ​

9. statsmodels回归分析 ​

10. 高阶与分类变量实例 ​

高阶回归 ​

分类变量 ​

11. 案例: 汽车价格预测任务 ​

任务概述 ​

缺失值填充 ​

特征相关性 ​

预处理问题 ​

回归问题 ​

回归分析

11、回归分析

1. 回归分析概述

2. 回归方程定义

3. 误差项的定义

4. 最小二乘法的推导和求解

5. 回归方程求解最小例子

6. 回归直线拟合优度

判定系数

显著性检验

7. 多元与曲线回归问题

曲线回归分析

多重共线性

8. python工具包

9. statsmodels回归分析

10. 高阶与分类变量实例

高阶回归

分类变量

11. 案例: 汽车价格预测任务

任务概述

缺失值填充

特征相关性

预处理问题

回归问题