十四、线性模型

原文：Linear Models

译者：飞龙

协议：CC BY-NC-SA 4.0

线性（回归）建模是一种方法，将输出值预测为输入值的加权线性组合。

线性模型 - 概述

在最简单的情况下，我们试图拟合一条线，因此我们的模型具有以下形式：

y = ax + b

在上面的这个等式中，我们试图从一些其他数据变量x预测一些数据变量y，其中和b是参数，我们需要通过拟合模型得出，并分别反映模型（直线）的斜率和y`轴截距。

我们需要一些过程来寻找a和b。我们将使用 OLS 这样做 - 我们想要的a和b的值是那些满足 OLS 解的值 - 也就是使模型预测与我们的数据之间距离最小的值。

请注意，你需要x和y已知的数据，以便训练你的模型。

该方法也可以推广，包括例如使用更多特征来预测我们感兴趣的输出。

因此，我们将以一般形式重写我们的模型，如下所示：

y = a0 + a1 x1 + a2 x2 + ... + an xn + ε

在上面的等式中，a0是截距（与上面的b相同），a1到an是我们正在尝试学习的n个参数，因为数据的权重为x1到xn。我们的输出变量（我们试图预测的）仍然是y，我们引入了ε，这是误差，它基本上捕获了无法解释的变动。

线性模型实战

在下文中，我们将生成一些数据，其中两个特征D1和D2是相关的。

鉴于相关性，我们可以尝试从D1预测D2的值，我们将创建一个线性模型来实现。

使用上面第二种表示法，模型将采用以下形式：

D2 = a0 + a1 * D1

# 导入 
%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# Statmodels & patsy
import patsy
import statsmodels.api as sm
# 生成一些互相关数据
# 设置
corr = 0.75
covs = [[1, corr], [corr, 1]]
means = [0, 0]
# 生成数据
dat = np.random.multivariate_normal(means, covs, 1000)
# 查看我们生成的数据
plt.scatter(dat[:, 0], dat[:, 1], alpha=0.5);
#plt.scatter(dat[:, 0], dat[:, 1], alpha=0.5);

# 将数据放入 DataFrame
df = pd.DataFrame(dat, columns=['D1', 'D2'])
# 观察数据
df.head()

	D1	D2
0	0.305994	0.521119
1	0.262456	0.562350
2	-0.289970	0.608496
3	0.025279	0.523315
4	-0.561126	-1.320507

# 检查 D1 和 D2 之间的相关性（它与合成的东西匹配）
df.corr()

	D1	D2
D1	1.000000	0.747094
D2	0.747094	1.000000

Statsmodels & Patsy 的线性模型

Statsmodels 是 Python 中统计分析的模块。 Patsy 是一个有用的包，用于处理和描述统计模型。

这里是 statsmodels 和 patsy 的官方文档。

# Patsy 为我们提供了一种构建设计矩阵的简便方法
# 出于我们的目的，“设计矩阵”只是我们的预测变量和输出变量的有组织矩阵
outcome, predictors = patsy.dmatrices('D2 ~ D1', df)

如果检查“结果”和“预测变量”的类型，你会发现它们是DesignMatrix类型的自定义 patsy 对象。

如果你将它们打印出来，你会看到它们重置了 Pandas Series或DataFrames。

你可以将它们视为自定义的类数据帧对象，以便将其组织成用于建模的矩阵。

# 现在使用 statsmodels 初始化 OLS 线性模型
# 此步骤初始化模型，并提供数据（但实际上不计算模型）
mod = sm.OLS(outcome, predictors)

请注意，statsmodels，就像我们稍后会遇到的 scikit-learn 一样，使用面向对象的方法。

在这种方法中，你初始化将数据和方法存储在一起的复杂对象，为你提供有组织的方式来存储和检查数据和参数，拟合模型，然后甚至使用它们进行预测等等。

# 检查我们刚刚创建的模型对象的类型。
# 你还可以使用 tab 补全，来探索此对象的可用内容
type(mod)
# statsmodels.regression.linear_model.OLS
# 最后，拟合模型
res = mod.fit()
# 检查结果
print(res.summary())
'''
==============================================================================
Dep. Variable:                     D2   R-squared:                       0.558
Model:                            OLS   Adj. R-squared:                  0.558
Method:                 Least Squares   F-statistic:                     1261.
Date:                Tue, 06 Mar 2018   Prob (F-statistic):          3.32e-179
Time:                        00:18:34   Log-Likelihood:                -993.04
No. Observations:                1000   AIC:                             1990.
Df Residuals:                     998   BIC:                             2000.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      0.0011      0.021      0.054      0.957      -0.039       0.042
D1             0.7546      0.021     35.506      0.000       0.713       0.796
==============================================================================
Omnibus:                        1.644   Durbin-Watson:                   1.970
Prob(Omnibus):                  0.440   Jarque-Bera (JB):                1.508
Skew:                          -0.068   Prob(JB):                        0.470
Kurtosis:                       3.133   Cond. No.                         1.05
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
'''

解释输出

Statsmodels 为我们提供了大量信息！

最上面的部分主要是元数据：它包括模型类型，运行它的时间和日期等内容。

它还包括 R 方，它是模型能够捕获的方差量的总体汇总。这个值在 0-1 之间，约为 0.5 ，我们在这里看到，是一个非常高的值，表明模型拟合良好。

中间一列是实际的模型结果。

每行反映一个参数，并给出它的值（coef），误差（std err），统计检验结果，关于该参数是否是输出变量的显着预测值（t，其中相关的 p 值为’P > | t |），以及参数值的置信区间（[0.025 - 0.975]`）。

最后一个模型包括对数据运行的一些其他测试，可以帮助你检查输入数据的某些描述符，以及它们是否满足这种模型拟合的必要条件。

检查我们的模型

就模型本身而言，最有用的组件位于第二行，其中摘要给出了参数值，以及我们的预测变量的 p 值，在这种情况下是Intercept和D2。

从上面的结果中，我们可以获取参数的值，并获得以下模型：

D2 = -0.0284 + 0.7246 * D1

但是，我们还应该记住报告的统计测试，参数值是否显着（显着不同于零）的测试。使用 0.05 的 alpha 值，在这种情况下，D2参数值是显着的，但Intercept值不是。由于Intercept的参数值与零没有显着差异，我们可以决定不将它包含在我们的最终模型中。

因此，我们完成了模型：

D2 = 0.7246 * D1

有了这个模型，a1的值很可能为 0.7246，非常接近数据点的相关值，我们将其设置为 0.75！

## 绘制模型的拟合直线
# 绘制原始数据（像之前一样）
plt.scatter(df['D1'], df['D2'], alpha=0.3, label='Data');
# 生成和绘制模型的拟合直线
xs = np.arange(df['D1'].min(), df['D1'].max())
ys = 0.7246 * xs
plt.plot(xs, ys, '--k', linewidth=4, label='Model')
plt.xlabel('D1')
plt.xlabel('D2')
plt.legend();

使用多个预测值

上面的模型只使用了一个预测值，拟合了一条简单的直线，因此实际上模仿了我们用于拟合直线的先前方法。

我们还可以拟合多个预测变量，这就是 patsy 和 statsmodel 的强大功能，因为这些函数你和更复杂的模型，包括我们想要的许多参数，也处理相关特征的某些方面，等等。

在这里，我们将向数据帧添加一个新变量，并使用两个预测变量拟合 OLS 模型。

# 向 df 添加数据的新列
df['D3'] = pd.Series(np.random.randn(1000), index=df.index)
df.head()

	D1	D2	D3
0	0.305994	0.521119	0.019043
1	0.262456	0.562350	-0.631288
2	-0.289970	0.608496	0.794971
3	0.025279	0.523315	0.447560
4	-0.561126	-1.320507	-1.560199

# 从 D2 和 D3 预测 D1
outcome, predictors = patsy.dmatrices('D1 ~ D2 + D3', df)
mod = sm.OLS(outcome, predictors)
res = mod.fit()
# 检查模型拟合摘要
print(res.summary())
'''
==============================================================================
Dep. Variable:                     D1   R-squared:                       0.558
Model:                            OLS   Adj. R-squared:                  0.557
Method:                 Least Squares   F-statistic:                     630.1
Date:                Tue, 06 Mar 2018   Prob (F-statistic):          1.25e-177
Time:                        00:18:34   Log-Likelihood:                -982.91
No. Observations:                1000   AIC:                             1972.
Df Residuals:                     997   BIC:                             1987.
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     -0.0177      0.020     -0.866      0.387      -0.058       0.022
D2             0.7395      0.021     35.479      0.000       0.699       0.780
D3            -0.0121      0.021     -0.587      0.557      -0.052       0.028
==============================================================================
Omnibus:                        0.273   Durbin-Watson:                   1.993
Prob(Omnibus):                  0.872   Jarque-Bera (JB):                0.357
Skew:                          -0.019   Prob(JB):                        0.836
Kurtosis:                       2.916   Cond. No.                         1.04
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
'''

请注意，如上所述，statsmodels 是一种功能强大的通用 OLS 模型方法。

你可以进一步研究如何包含其他功能，例如输入变量之间的交互，等等。

sklearn 线性回归

Scikit-learn 也具有线性回归模型的实现。

在这里，我们将使用 sklearn 而不是 statsmodels，快速演示运行与上述相同的线性 OLS 模型。

sklearn 中的线性回归。

# sklearn 的线性模型
from sklearn import linear_model
# 转换数据的形状，便于 sklearn 使用
d1 = np.reshape(df.D1.values, [len(df.D1), 1])
d2 = np.reshape(df.D2.values, [len(df.D2), 1])
d3 = np.reshape(df.D3.values, [len(df.D3), 1])
# 初始化线性回归模型
reg = linear_model.LinearRegression()
# 拟合线性回归模型
reg.fit(d2, d1) #d1 = a0 + a1*d2
# LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)
# 检查这个结果
# 如果将这些与我们上面使用的 statsmodel 进行比较，它们确实是相同的
print(reg.intercept_[0])
print(reg.coef_[0][0])
'''
-0.0179497213127
0.739686393719
'''

（在 sklearn 中）使用多个预测变量

# 初始化并拟合线性模型
# d1 = a1*d2 + a2*d3 + a0
reg = linear_model.LinearRegression()
reg.fit(np.hstack([d2, d3]), d1)
# LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)
np.hstack([d2, d3]).shape
# (1000, 2)
# 检查这个结果
# 如果将这些与我们上面使用的 statsmodel 进行比较，它们确实是相同的
print('Intercept: \t', reg.intercept_[0])
print('Theta D2 :\t', reg.coef_[0][0])
print('Theta D3 :\t', reg.coef_[0][1])
'''
Intercept:      -0.0177498825239
Theta D2 :     0.739473377018
Theta D3 :     -0.0120672090985
'''