统计学11-多元线性回归

在统计学10-回归一文中介绍了一元线性回归的概念。假设我们现在有多个解释变量，如何构造多元线性回归模型呢？

第一个实例

在这个例子里，我们不仅考虑房子的大小，还考虑其他的解释变量：
neighborhood：分类变量
bedrooms：数值变量
bathrooms：数值变量
style：分类变量

导入库和数据

import numpy as np
import pandas as pd
import statsmodels.api as sms;

df = pd.read_csv('./house_prices.csv')

添加截距

df['intercept'] = 1

分别用简单线性回归拟合回归模型

lm = sms.OLS(df['price'], df[['intercept','area']])
results = lm.fit()
results.summary()

lm = sms.OLS(df['price'], df[['intercept','bedrooms']])
results = lm.fit()
results.summary()

lm = sms.OLS(df['price'], df[['intercept','bathrooms']])
results = lm.fit()
results.summary()

分析结果

都具有统计显著性。
如果把三个简单线性回归模型的决定系数加起来，值大于1。
在多元线性回归里，我们要如何找到 “正确” 的系数
在简单线性回归那一节中，你知道了我们为什么要最小化每个实际数据点与模型预测值之间的平方距离。

用所有变量建立多元线性回归模型

lm = sms.OLS(df['price'], df[['intercept','area','bedrooms','bathrooms']])
results = lm.fit()
results.summary()

再次分析结果

卧室数量：无统计显著性。
浴室数量：无统计显著性。
面积：具有统计显著性。
决定系数小于1。
不能拟合像“风格”之类的分类变量。

另一种求解方式

但在多元线性回归中，我们要找的数据点所处的空间实际上不仅仅是二维的。
可以以线性代数的方式写出求解公式：
$\bar{β}$ = (X'X)^-X'y

我们将使用 statsmodel 来求得系数，这些系数会与上一概念涉及的系数相似，但你还要用上述等式来解出该系数，证明所得系数不是凭空捏造出来的。

X = df[['intercept', 'bathrooms', 'bedrooms', 'area']]
y = df['price']

np.dot(np.dot(np.linalg(np.dot(X.transpose(), X)), X.transpose()), y)

numpy.matrix.transpose：矩阵转置
np.linalg: 矩阵求逆
np.dot：矩阵相乘

启示

我们可以将各系数解释为： 在模型其它变量不变的情况下，解释变量每增加一个单位，反应变量会随之增加的预测幅度。
当分别单独处理的时候，有些解释变量是统计显著的，而当将他们一起拟合的时候，一些解释变量反而不是统计显著的了。
根据我们的现有知识，无法处理分类变量。

虚拟变量

要往线性模型里添加分类变量，就需要把分类变量转变为 虚拟变量。即为每一个值创建一列：

vv.png

转化后，你需要舍弃一个 虚拟列，才能得到满秩矩阵【1】。

回想一下回归系数闭式解那节内容，我们得用 (X'X)^-X'y来预测 β 值。

为了逆转(X'X)，矩阵X 一定要是满秩的，也就是说，所有X 的列都必须线性独立。

创建虚拟变量时，如果你不舍弃掉一列，那你就得不到稳定解，从 python 里得到的结果也不会可靠到哪去。

本文重点在于 如果你要用 0 、1 编码来创建虚拟变量，你就得舍弃一个虚拟列，确保所得矩阵是满秩的（这样你从 python 里得到的解才会是可靠的。）

之所以要这么做，原因就在于线性代数的本质，更具体地说，要逆转矩阵，你手里的矩阵必须是满秩的（也就是所有列都得线性独立），因此，你得舍弃掉一个虚拟列，方能得到线性独立的各列（和一个满秩矩阵）。

用python演示创建虚拟变量

pd_dummies = pd.get_dummies(df['neighborhood'])
df_new = df.join(pd_dummies)
lm = sms.OLS(df_new['price'], df_new[['intercept','A']])
results = lm.fit()
results.summary()

如果删除的虚拟变量为A，剩下的为‘B’和‘C‘，则可以解释为当neighborhood为“A”时，预测价格为intercept。B和C的斜率解释为和A相比，B和C对反应变量的变化量的影响。
还可以对数据进行可视化以直观的观察分类变量对反应变量的影响：

plt.hist(df_new.query("C == 1")['price'], alpha = 0.3, label = 'C');
plt.hist(df_new.query("A == 1")['price'], alpha = 0.3, label = 'A');
plt.hist(df_new.query("B == 1")['price'], alpha = 0.3, label = 'B');
plt.legend();

其他引入虚拟变量的方法

要用 1、0 编码预测基准类别，你要用到截距，此时斜率是与基准分类比较下的变化；在SAS中，用 1、0、-1 编码，你则得让各分类系数分别乘以 -1，从而填补缺失的分类。

因为如果以 1、0、-1 编码，对于类别为基准类别的数据行来说，所有的虚拟变量都为-1，因此如果假设类别列表为A，B，C，且基准类别为C：
C的平均值 = intercept - B的平均值 - A的平均值。

导入数据和库

import numpy as np
import pandas as pd
import statsmodels.api as sm

df = pd.read_csv('./house_prices.csv')
df2 = df.copy()
df.head(10)

定义创建1，0，-1编码的虚拟变量。

## The below function creates 1, 0, -1 coded dummy variables.

def dummy_cat(df, col):
    '''
    INPUT:
    df - the dataframe where col is stored
    col - the categorical column you want to dummy (as a string)
    OUTPUT:
    df - the dataframe with the added columns
         for dummy variables using 1, 0, -1 coding
    '''
    for idx, val_0 in enumerate(df[col].unique()):
        if idx + 1 < df[col].nunique():            
            df[val_0] = df[col].apply(lambda x: 1 if x == val_0 else 0)
        else:    
            df[val_0] = df[col].apply(lambda x: -1 if x == val_0 else 0)
            for idx, val_1 in enumerate(df[col].unique()):
                if idx + 1 < df[col].nunique():
                    df[val_1] = df[val_0] + df[val_1]
                else:
                    del df[val_1]
    return df

为分类变量调用方法

new_df = dummy_cat(df, 'style') # Use on style
new_df.head(10)

线性回归

new_df['intercept'] = 1

lm = sm.OLS(new_df['price'], new_df[['intercept', 'ranch', 'victorian']])
results = lm.fit()
results.summary()

置信区间不重叠->显著不同？

模型假设及相应的解决方法

之前的视频曾提到统计学习简介一书，其中就提到了如下五个假设：

因变量-自变量关系的非线性
误差项的相关性
非恒定方差和正态分布误差
异常值/高杠杆点
共线性

这里总结了判断上述问题是否存在的方法以及相应的解决办法。这是统计学家面试时经常提出的问题，但该问题是否有实际意义取决于你创建模型的目的。在接下来的概念中，我们会更仔细地研究某些相关知识点。下方列出了各知识点的详尽介绍，我们先来仔细看看下文涉及的每一项。

线性

线性是假设因变量和自变量之间真的存在可用线性模型解释的关系。如果线性假设不为真，那你的预测结果就不会很准确，此外，与系数有关的线性关系也就没什么用了。

为了评估某段线性关系是否合理，一个很实用的方法是做预测值 $\hat{y}$ 的残差(y - $\hat{y}$ )图。如果图中出现多个曲线部分，那就意味着线性模型实际上可能并不拟合数据，自变量和因变量存在其它关系。创建非线性模型的办法有很多（甚至可以线性模型的形式来创建）。

在本页底部的图片里，这些称为偏差模型。理想来说，我们想要的是像图片左上角残差图那样的随机散点图。

非恒定方差和正态分布误差

你预测的值不同，得到的预测值范围也不同，那就意味着方差不恒定。非恒定方差对预测好坏影响不大，但会导致置信区间和 p 值不准确，这种时候，在预测值接近实际值的那部分区域，系数的置信区间会太泛，而在预测值较远离实际值的区域则会太窄。

通常来说，对数函数（或使用其它反应变量的变换方式）能够 “摆脱” 非恒定方差，而要选择合适的变换方式，我们一般会用 Box-Cox。

用预测值的残差图也可以评估非恒定方差。在本页底部的图片中，非恒定方差的标签为 异方差。理想来说，我们要的是一个有异方差残差的无偏模型（其异方差残差在一定数值范围内保持不变）。

虽然本文并不探讨残差的正态性，如果你想创建可靠的置信区间，正态性回归假设就十分重要了，更多相关信息详见这里。

异常值/杠杆点

异常值和杠杆点是远离数据正常趋势的点。这些点会对你的解造成很大的影响，在现实中，这些点甚至可能是错误的。如果从不同来源收集数据，你就可能在记录或收集过程中造成某些数据值出错。

异常值也可能是准确真实的数据点，而不一定是测量或数据输入错误。在这种情况下，'修复'就会变得更为主观。要如何处理这些异常值往往取决于你的分析目的。线性模型，特别是使用最小二乘法的线性模型，比较容易受到影响，也就是说，大异常值可能会大幅度地左右我们的结果。当然，异常值也有一些解决技巧，也就是我们常说的 正则化。本课不会谈及这些技巧，但在机器学习纳米学位免费课程中，我们对这些技巧做了粗略的介绍。

而在宾夕法尼亚州立大学提供的完整回归课程里，就有特别长的篇幅在探讨杠杆点的问题，详见这里。

共线性（多重共线性）

如果我们的自变量彼此相关，就会出现多重共线性。多重共线性的一个主要问题在于：它会导致简单线性回归系数偏离我们想要的方向。

要判断是否有多重共线性，最常见的办法是借助二变量图或 方差膨胀因子 (即 VIFs)。

resid-plots.gif

多重共线性与VIF

多元线性回归模型的重要假设之一是解释变量与反应变量相关，而不是彼此相关。例如前面的例子，bedrooms、bathrooms和房子大小是相关的。

可以用可视化方法查看解释变量彼此之间的相关性：

sb.pairplot(df[['area','bedrooms','bathrooms']]);

除了作图之外，还可以通过方差膨胀因子VIFs 检验解释变量之间的相关性。

如果一个属性的VIF过大，那么可以删除和这个属性corr值较大的属性。

一般情况下，只需了解到多重共线性的原理，以及解决策略即可。其他内容可以作为选修内容。

计算VIF

y, X = dmatrices('price ~ bedrooms + bathrooms + area' , df, return_type='dataframe')

vif = pd.DataFrame()
vif["VIF Factor"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
vif["features"] = X.columns

写成函数：

def vif_calculator(df, response):
    '''
    INPUT:
    df - a dataframe holding the x and y-variables
    response - the column name of the response as a string
    OUTPUT:
    vif - a dataframe of the vifs
    '''
    df2 = df.drop(response, axis = 1, inplace=False)
    features = "+".join(df2.columns)
    y, X = dmatrices(response + ' ~' + features, df, return_type='dataframe')
    vif = pd.DataFrame()
    vif["VIF Factor"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
    vif["features"] = X.columns
    vif = vif.round(1)
    return vif

衡量VIF

一般规则是，如果VIF大于10，那么模式就会出现共线性。VIF_i = 1 / (1- $R_i$ ²)

这里 $R_i$ ²是通过拟合以下线性模型得到的：
x_i = b₀ + b₁x₁ + b₂x₂ + ... + b_nx_n

也就是说如果 $R_i$ ²很大，即 $x_i$ 和其他变量的相关性具有统计显著性，那么VIF的值就会很大。

通常，不会发现只有某一个变量的VIF较大，因为如果VIF较大，意味着有两个或以上的变量具有相关性，这些变量都会有较高的VIF。一般做法是删除最小兴趣的变量。

高阶项

有时候可能需要拟合类似于非线性模型的模型。可以通过在模型中增加高阶项，如交叉项（x₁* x₂），二次式，三次式和更高阶的项。虽然这些高阶项能让我们较好的预测反应变量，但是应该尽量避免使用他们，因为：

高阶项增大了模型的复杂度。
使得模型更加难于解释。（解释变量单位的平方往往没有现实意义）

python 演示

# 增加二次式
# 注意，一定也要把各自的低阶项加入模型
df['bedrooms_squared'] = df['bedrooms'] * df['bed_rooms']
lm = sm.OLS(new_df['price'], new_df[['intercept', 'bedrooms', 'bedrooms_squared']])
results = lm.fit()
results.summary()

此时问题变得复杂：如何解释”卧室数量的平方“？但是决定系数并没有明显的变化！

要怎么确定高阶项？

在创建具备二阶、三阶甚至更高阶变量的模型时，基本上我们要看的是解释变量和反应变量的关系里有多少条曲线。

如果像下图一样，图中曲线折了一个弯，那你会想添加一个二阶项到模型里，因为我们可以明显地发现，只用一条线来拟合这个关系并不能得到最佳效果。

quadraticlinearregression.png

当解释变量和反应变量的关系折了两个弯，我们就会想添加一个三阶的关系，如下图所示。

resid2.jpg

详细的讨论可以参考： https://tamino.wordpress.com/2011/03/31/so-what

解释交叉项

通过上一节的讨论，我们了解了如何通过观察散点图以确定是否采用二次式或者三次式来拟合模型。本节介绍如何确定应该添加交叉项。

当添加交叉项时，我们会认为变量X₁与Y之间的相关性，依赖于另一个变量X₂。

没有交叉项的情况

MR.png

这里我们可以看到，对于街区A和街区B：

得到的直线都能很好的拟合这些点。
得到的直线的斜率是一样的。也就是说，二者是平行的，街区和面积不具有相关性。

有交叉项的情况

MR2.png

这里我们可以看到，对于街区A和街区B：

得到的直线都能很好的拟合这些点。
得到的直线的斜率不同。也就是说，价格和面积的关系取决于房子位于哪个社区，即变量Area与Price之间的相关性，依赖于另一个变量Neighborhood。

特征工程与特征选择

当我们进行线性回归或者机器学习时，模型能给我们带来多大收获取决于特征工程和特征选择。特征选择是指选择应该在模型中使用的解释变量。
本文已经介绍了一些选择变量的方法，如显著性水平（p值）和VIF，另一类方法是交叉验证或正则化。
特征工程包括使用缺失值或删除行，或者用初始特征的不同比例创建新的一行（称为scaling：log，exp，standardizing，quadratics），将文本和图片转化为数字，或者创建虚拟变量。
用初始特征的不同比例创建新的一行有助于对变量进行预测，但是可能会使得回归系数无法解释，应该慎用。

Sebastian 和 Katie 讲授的这个课程是个很好的学习途径，能帮我们更好地理解本文涉及的许多概念，也对我们接下来更深入地理解机器学习很有帮助。

说到这里，我们已经了解了特征工程的多种方式，如虚拟变量和缺失值填充。特别的，如果需要scaling数据，有两个非常著名的库：

Sklearn，若要查看有关特征工程的 Sklearn 文献，请点这里。从scikit-learn中导入pre-processing，可以看到scale方法。然后还可以看到一个叫做“StandardScaler”的东东，他们可以处理不怎么像标准正态分布的数据集，将他们减去mean再除以标准差。此外还有许多其他用途的scaler，MinMax scaler，MaxAbs Scaler，以及其他处理数据的方式quantile transformer，normalization，binarization，encoding categorical features，imputing missing values，custom transformations。
注意到这些技巧大部分都包括拟合转换，即先拟合数据然后按照一定的方式对数据进行转换，这是scikit-learn函数的常见用法。具体来说，preprocessing常依序调用fit、transform对数据进行拟合和转化，或者直接调用fit_transform方法。这样做的目的是使得模型更加健壮，而不会收到scaling的影响。回归是一种非常容易受到离散点影响的一种技巧。这样做可以提高预测的准确度，但是确实会改变解释，

import numpy as np
import pandas as pd
import statsmodels.api as sm;
import sklearn.preprocessing as p

df = pd.DataFrame({'response': [2.4, 3.3, -4.2, 5.6, 1.5, 8.7], 
                         'x1': ['yes','no','yes','maybe','no','yes'],
                         'x2': [-1,-3,np.nan, 0, np.nan, 1],
                         'x3': [2.4, 15, 3.3, 2.4, 1.8, 0.4],
                         'x4': [np.nan, np.nan, 1, 1, 1, 1],
                         'x5': ['A', 'B', np.nan, 'A', 'A', 'A']})

p.scale(df[['x2','x3','x4']])

min_max_scaler = p.MinMaxScaler()
min_max_scaler .fit_transform(df[['x2','x3','x4']])

# 拟合失败
df['intercept'] = 1
lm = sm.OLS(df['response'], df[['intercept','x2','x3','x4']])
results = lm.fit()
results.summary()

# 另一种常用的缩放特征的方法是减去均值并除以标准偏差
norm = p.StandardScaler()
norm.fit(df[['x2','x3','x4']])
norm.transform(df[['x2','x3','x4']])

sklearn.preprocessing包提供了一些通用的工具方法和transformer类，用于将原始特征向量转换成易于评估的形式。

pandas，有关在 pandas 里处理缺失值的其它文献，请点这里。

df.fillna(df.mean(), inplace = True)
lm = sm.OLS(df['response'], df[['intercept','x2','x3','x4']])
results = lm.fit()
results.summary()

模型评估的测量标准

模型评估的测量标准有很多：R-squared，MSE，AIC，BIC和马洛斯的CP值等等。这些标准的值取决于数据的缩放比例和用于预测反应变量的变量。有些时候，人们倾向于在实践中采用R-squared和MSE，但是实际上这些方法在模型拟合中存在误导性，如果我们只在同一个数据集上进行验证，无论我们像模型中添加什么变量，都会改善模型，也就是说R-suqared值会上升，MSE会下降。那么，我们怎么知道增加变量是否能真的改善了模型与数据的拟合度呢？有一个非常有效的技巧，称为交叉验证。

交叉验证（模型评估）

在交叉验证中，我们对数据子集训练回归模型，称为训练数据。我们用称为测试数据的数据子集测试模型如何较好的运行。在交叉验证中，我们通过多次切分数据集，以确保我们的模型能够通用化（generalize）。

k 折交叉验证由 Sebastian 的视频讲解，属于优达学城的免费课程，很好地解释了交叉验证的工作原理。

在交叉验证中，我们希望训练数据和测试数据的两个集合都做到最大化。训练数据集越大，学习效果就越好；测试数据集越大，验证效果越好。显然这里需要折衷。
折衷的要点是将训练数据平分到相同大小的k个容器中，比方说有200个数据，如果k=10，那么每个容器中有20条数据。在k阶交叉验证中，将运行k次单独的学习实验，在每次实验中，从k个子集中挑选一个作为验证集，剩下k-1容器放在一起作为训练集，并在验证集上验证性能。交叉验证的要点是这个操作会运行多次。在这个例子中，k=10，因此运行10次，然后对十个不同的运行结果的表现进行平均，也就是说将这k次实验的测试结果取平均值。这样，学习算法的评估将更加准确，而且可以对全部数据进行验证。

可以在这里找到参考文档。

通过python演示数据集的切分

引入库和数据

import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
df = pd.read_csv('./data/house_prices.csv')

切分数据集
这一步必须在特征工程之前做。因为如果你在切分前得到数据的平均数，那么这就意味这测试集的数据包含了训练集的数据，而这不是我们期待的结果。

X = df[['neighborhood','area','bedroom','bathroom']]
y = df[['price']]
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33, random_state=42)

初始化线性回归模型

lm = LinearRegression()

对训练数据集进行拟合

lm.fit(X_train, y_train)

对测试数据集进行预测

y_preds = lm.predict(X_tests)

对预测结果进行打分
应该用测试数据集而不是训练数据集计算分数，因为测试数据集体现了模型在没有遇到过的数据集上的表现。
根据分数，可以选择在模型中应该保留哪些属性。

# 默认的打分方法为R-squared
lm.score(X_test, y_test)

案例研究

-导入数据集和库

import numpy as np
import pandas as pd
import statsmodels.api as sm
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
from patsy import dmatrices
import matplotlib.pyplot as plt
%matplotlib inline

np.random.seed(42)
boston_data = load_boston()
df = pd.DataFrame()
df['MedianHomePrice'] = boston_data.target
df2 = pd.DataFrame(boston_data.data)
df2.columns = boston_data.feature_names
df = df.join(df2)
df.head()

准备X，Y

y = df['MedianHomePrice']
X = df.drop(['MedianHomePrice'], axis=1)
X.head(10)

split 数据

X_train, X_test, y_train, y_test = train_test_split(X,y,test_size=0.20, random_state=0)
X_train.head()

获取数据集中每个特征的汇总

X_train.isnull().sum()
X_train.info()
X_train['AGE'].median()
X_train['CHAS'].mean()

一般情况下，可以查看：

数据集的缺失值数
数据集的分类变量数
数据集某个变量的中位数
。。。

使用 corr方法对每个变量进行相互比较

df2 = pd.DataFrame()
df2['MedianHomePrice'] = y_train
df2[list(df.columns)[1:]] = X_train
df2.corr()

通过corr的结果可以知道两方面的信息：

和反应变量相关性较强的解释变量。
解释变量之间是否有较强的相关性。

corr的使用可以参照这里。

缩放数据集中的所有 x 变量

# 使用StandardScaler来缩放数据集中的所有 x 变量
norm = StandardScaler()
norm.fit(X_train)
# 将结果存储在 X_scaled_train中
X_scaled_train = norm.transform(X_train)

StandardScaler的使用方法参照这里。

创建一个 pandas 数据帧存储缩放的 x 变量以及training响应变量

X_scaled_train = pd.DataFrame(X_scaled_train)
X_scaled_train.columns = list(df.columns)[1:]
df3 = X_scaled_train.join(y_train.reset_index())

用所有的缩放特征来拟合线性模型，以预测此响应（平均房价)

df3['intercept'] = 1
X_vars_full = df3.drop(['MedianHomePrice','index'] , axis=1, inplace=False)
lm = sm.OLS(df3['MedianHomePrice'], X_vars_full)
results = lm.fit()
results.summary()

观察由线性模型得出的p值，标注变量为具有统计显著性或无显著性。

用于计算vif的函数

def vif_calculator(df, response):
    '''
    INPUT:
    df - a dataframe holding the x and y-variables
    response - the column name of the response as a string
    OUTPUT:
    vif - a dataframe of the vifs
    '''
    df2 = df.drop(response, axis = 1, inplace=False)
    features = "+".join(df2.columns)
    y, X = dmatrices(response + ' ~' + features, df, return_type='dataframe')
    vif = pd.DataFrame()
    vif["VIF Factor"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
    vif["features"] = X.columns
    vif = vif.round(1)
    return vif

计算vif

vif_calculator(df3, 'MedianHomePrice')

可以根据经验设置VIF值的阈值，删除VIF高于该阈值的变量。如果需要的或者在乎的变量的VIF较高，那么就需要考虑删除和VIF较高的变量相关性较强的变量。

为模型打分

X_data = df.drop('MedianHomePrice', axis=1, inplace=False)
X_train, X_test, y_train, y_test = train_test_split(
                 X_data, df['MedianHomePrice'], test_size=0.2, random_state=0)


lm_full = LinearRegression()
lm_full.fit(X_train, y_train)
# 返回预测结果的R-sqaured值
lm_full.score(X_test, y_test)

还可以对不同的模型进行比较(假设RAD的p值较大)

X_train_red = X_train.drop(['AGE','NOX','TAX'] , axis=1, inplace=False)
X_test_red = X_test.drop(['AGE','NOX','TAX'] , axis=1, inplace=False)

X_train_red2 = X_train.drop(['AGE','NOX','TAX','RAD'] , axis=1, inplace=False)
X_test_red2 = X_test.drop(['AGE','NOX','TAX','RAD'] , axis=1, inplace=False)

lm_red = LinearRegression()
lm_red.fit(X_train_red, y_train)
print(lm_red.score(X_test_red, y_test))

lm_red2 = LinearRegression()
lm_red2.fit(X_train_red2, y_train)
print(lm_red2.score(X_test_red2, y_test))

LinearRegression的使用可以参照文档。

决定系数高的模型可能只是因为变量较多，而不意味着这个模型真的能在新（测试）数据集上有更好的效果。

参考书

多元线性回归涉及很多概念和知识，统计学习简介尤其是该书的第三章是很好的资源。
对于线性代数，可观看可汗学院线性代数免费课程。

【1】满秩矩阵是一个很重要的概念, 它是判断一个矩阵是否可逆。