python逐步回归函数,Python做回归

Python的函数都有哪些

【常见的内置函数】

1、enumerate(iterable,start=0)

是python的内置函数，是枚举、列举的意思，对于一个可迭代的(iterable)/可遍历的对象(如列表、字符串)，enumerate将其组成一个索引序列，利用它可以同时获得索引和值。

2、zip(*iterables,strict=False)

用于将可迭代的对象作为参数，将对象中对应的元素打包成一个个元组，然后返回由这些元组组成的列表。如果各个迭代器的元素个数不一致，则返回列表长度与最短的对象相同，利用*号操作符，可以将元组解压为列表。

3、filter(function,iterable)

filter是将一个序列进行过滤，返回迭代器的对象，去除不满足条件的序列。

4、isinstance(object,classinfo)

是用来判断某一个变量或者是对象是不是属于某种类型的一个函数，如果参数object是classinfo的实例，或者object是classinfo类的子类的一个实例，

返回True。如果object不是一个给定类型的的对象， 则返回结果总是False

5、eval(expression[,globals[,locals]])

用来将字符串str当成有效的表达式来求值并返回计算结果，表达式解析参数expression并作为Python表达式进行求值(从技术上说是一个条件列表)，采用globals和locals字典作为全局和局部命名空间。

【常用的句式】

1、format字符串格式化

format把字符串当成一个模板，通过传入的参数进行格式化，非常实用且强大。

2、连接字符串

常使用+连接两个字符串。

3、if...else条件语句

Python条件语句是通过一条或多条语句的执行结果(True或者False)来决定执行的代码块。其中if...else语句用来执行需要判断的情形。

4、for...in、while循环语句

循环语句就是遍历一个序列，循环去执行某个操作，Python中的循环语句有for和while。

5、import导入其他脚本的功能

有时需要使用另一个python文件中的脚本，这其实很简单，就像使用import关键字导入任何模块一样。

相关问答

Q1: 请教怎样用Python做stepwise回归

所说所有的变量都是对象。 对象在python里，其实是一个指针，指向一个数据结构，数据结构里有属性，有方法。

对象通常就是指变量。从面向对象OO的概念来讲，对象是类的一个实例。在python里很简单，对象就是变量。

class A:

myname="class a"

上面就是一个类。不是对象

a=A()

这里变量a就是一个对象。

它有一个属性（类属性），myname，你可以显示出来

print a.myname

所以，你看到一个变量后面跟点一个小数点。那么小数点后面。

Q2: 如何用python作空间自回归模型

基本形式

线性模型(linear model)就是试图通过属性的线性组合来进行预测的函数，基本形式如下：

f(x)=wTx+b

许多非线性模型可在线性模型的基础上通过引入层结构或者高维映射（比如核方法）来解决。线性模型有很好的解释性。

线性回归

线性回归要求均方误差最小:

(w∗,b∗)=argmin∑i=1m(f(xi)−yi)2

均方误差有很好的几何意义，它对应了常用的欧式距离(Euclidean distance)。基于均方误差最小化来进行模型求解称为最小二乘法(least square method)，线性回归中，最小二乘发就是试图找到一条直线，使得所有样本到直线的欧式距离之和最小。

我们把上式写成矩阵的形式：

w∗=argmin(y−Xw)T(y−Xw)

这里我们把b融合到w中，X中最后再加一列1。为了求最小值，我们对w求导并令其为0：

2XT(Xw−y)=0

当XTX为满秩矩阵(full-rank matrix)时是可逆的。此时：

w=(XTX)−1XTy

令xi=(xi,1)，可以得到线性回归模型：

f(xi)=xTi(XTX)−1XTy

Q3: python多元线性回归怎么计算

1、什么是多元线性回归模型？

当y值的影响因素不唯一时,采用多元线性回归模型。

y =y=β0+β1x1+β2x2+...+βnxn

例如商品的销售额可能不电视广告投入,收音机广告投入,报纸广告投入有关系,可以有 sales =β0+β1*TV+β2* radio+β3*newspaper.

2、使用pandas来读取数据

pandas 是一个用于数据探索、数据分析和数据处理的python库

[python] view plain copy

import pandas as pd

[html] view plain copy

pre name="code" class="python"# read csv file directly from a URL and save the results

data = pd.read_csv('/home/lulei/Advertising.csv')

# display the first 5 rows

data.head()

上面代码的运行结果：

    TV  Radio  Newspaper  Sales

0  230.1   37.8       69.2   22.1

1   44.5   39.3       45.1   10.4

2   17.2   45.9       69.3    9.3

3  151.5   41.3       58.5   18.5

4  180.8   10.8       58.4   12.9

上面显示的结果类似一个电子表格，这个结构称为Pandas的数据帧(data frame)，类型全称：pandas.core.frame.DataFrame.

pandas的两个主要数据结构：Series和DataFrame：

Series类似于一维数组，它有一组数据以及一组与之相关的数据标签(即索引)组成。

DataFrame是一个表格型的数据结构，它含有一组有序的列，每列可以是不同的值类型。DataFrame既有行索引也有列索引，它可以被看做由Series组成的字典。

[python] view plain copy

# display the last 5 rows

data.tail()

只显示结果的末尾5行

       TV  Radio  Newspaper  Sales

195   38.2    3.7       13.8    7.6

196   94.2    4.9        8.1    9.7

197  177.0    9.3        6.4   12.8

198  283.6   42.0       66.2   25.5

199  232.1    8.6        8.7   13.4

[html] view plain copy

# check the shape of the DataFrame(rows, colums)

data.shape

查看DataFrame的形状,注意第一列的叫索引，和数据库某个表中的第一列类似。

(200,4) 

3、分析数据

特征：

TV：对于一个给定市场中单一产品，用于电视上的广告费用（以千为单位）

Radio：在广播媒体上投资的广告费用

Newspaper：用于报纸媒体的广告费用

响应：

Sales：对应产品的销量

在这个案例中，我们通过不同的广告投入，预测产品销量。因为响应变量是一个连续的值，所以这个问题是一个回归问题。数据集一共有200个观测值，每一组观测对应一个市场的情况。

注意：这里推荐使用的是seaborn包。网上说这个包的数据可视化效果比较好看。其实seaborn也应该属于matplotlib的内部包。只是需要再次的单独安装。

[python] view plain copy

import seaborn as sns

import matplotlib.pyplot as plt

# visualize the relationship between the features and the response using scatterplots

sns.pairplot(data, x_vars=['TV','Radio','Newspaper'], y_vars='Sales', size=7, aspect=0.8)

plt.show()#注意必须加上这一句，否则无法显示。

[html] view plain copy

这里选择TV、Radio、Newspaper 作为特征，Sales作为观测值

[html] view plain copy

返回的结果：

seaborn的pairplot函数绘制X的每一维度和对应Y的散点图。通过设置size和aspect参数来调节显示的大小和比例。可以从图中看出，TV特征和销量是有比较强的线性关系的，而Radio和Sales线性关系弱一些，Newspaper和Sales线性关系更弱。通过加入一个参数kind='reg'，seaborn可以添加一条最佳拟合直线和95%的置信带。

[python] view plain copy

sns.pairplot(data, x_vars=['TV','Radio','Newspaper'], y_vars='Sales', size=7, aspect=0.8, kind='reg')

plt.show()

结果显示如下：

4、线性回归模型

优点：快速；没有调节参数；可轻易解释；可理解。

缺点：相比其他复杂一些的模型，其预测准确率不是太高，因为它假设特征和响应之间存在确定的线性关系，这种假设对于非线性的关系，线性回归模型显然不能很好的对这种数据建模。

线性模型表达式： y=β0+β1x1+β2x2+...+βnxn 其中

y是响应

β0是截距

β1是x1的系数，以此类推

在这个案例中： y=β0+β1∗TV+β2∗Radio+...+βn∗Newspaper

(1)、使用pandas来构建X(特征向量)和y(标签列)

scikit-learn要求X是一个特征矩阵，y是一个NumPy向量。

pandas构建在NumPy之上。

因此，X可以是pandas的DataFrame，y可以是pandas的Series，scikit-learn可以理解这种结构。

[python] view plain copy

#create a python list of feature names

feature_cols = ['TV', 'Radio', 'Newspaper']

# use the list to select a subset of the original DataFrame

X = data[feature_cols]

# equivalent command to do this in one line

X = data[['TV', 'Radio', 'Newspaper']]

# print the first 5 rows

print X.head()

# check the type and shape of X

print type(X)

print X.shape

输出结果如下：

     TV  Radio  Newspaper

0  230.1   37.8       69.2

1   44.5   39.3       45.1

2   17.2   45.9       69.3

3  151.5   41.3       58.5

4  180.8   10.8       58.4

class 'pandas.core.frame.DataFrame'

(200, 3)

[python] view plain copy

# select a Series from the DataFrame

y = data['Sales']

# equivalent command that works if there are no spaces in the column name

y = data.Sales

# print the first 5 values

print y.head()

输出的结果如下：

0    22.1

1    10.4

2     9.3

3    18.5

4    12.9

Name: Sales

（2）、构建训练集与测试集

[html] view plain copy

pre name="code" class="python"span style="font-size:14px;"##构造训练集和测试集

from sklearn.cross_validation import train_test_split  #这里是引用了交叉验证

X_train,X_test, y_train, y_test = train_test_split(X, y, random_state=1)

#default split is 75% for training and 25% for testing

[html] view plain copy

print X_train.shape

print y_train.shape

print X_test.shape

print y_test.shape

输出结果如下：

(150, 3)

(150,)

(50, 3)

(50,)

注：上面的结果是由train_test_spilit()得到的，但是我不知道为什么我的版本的sklearn包中居然报错：

ImportError                               Traceback (most recent call last)ipython-input-182-3eee51fcba5a in module()      1 ###构造训练集和测试集---- 2 from sklearn.cross_validation import train_test_split      3 #import sklearn.cross_validation      4 X_train,X_test, y_train, y_test = train_test_split(X, y, random_state=1)      5 # default split is 75% for training and 25% for testingImportError: cannot import name train_test_split

处理方法：1、我后来重新安装sklearn包。再一次调用时就没有错误了。

2、自己写函数来认为的随机构造训练集和测试集。(这个代码我会在最后附上。)

（3）sklearn的线性回归

[html] view plain copy

from sklearn.linear_model import LinearRegression

linreg = LinearRegression()

model=linreg.fit(X_train, y_train)

print model

print linreg.intercept_

print linreg.coef_

输出的结果如下：

LinearRegression(copy_X=True, fit_intercept=True, normalize=False)

2.66816623043

[ 0.04641001  0.19272538 -0.00349015]

[html] view plain copy

# pair the feature names with the coefficients

zip(feature_cols, linreg.coef_)

输出如下：

[('TV', 0.046410010869663267),

('Radio', 0.19272538367491721),

('Newspaper', -0.0034901506098328305)]

y=2.668+0.0464∗TV+0.192∗Radio-0.00349∗Newspaper

如何解释各个特征对应的系数的意义？

对于给定了Radio和Newspaper的广告投入，如果在TV广告上每多投入1个单位，对应销量将增加0.0466个单位。就是加入其它两个媒体投入固定，在TV广告上每增加1000美元（因为单位是1000美元），销量将增加46.6（因为单位是1000）。但是大家注意这里的newspaper的系数居然是负数，所以我们可以考虑不使用newspaper这个特征。这是后话，后面会提到的。

（4）、预测

[python] view plain copy

y_pred = linreg.predict(X_test)

print y_pred

[python] view plain copy

print type(y_pred)

输出结果如下：

[ 14.58678373   7.92397999  16.9497993   19.35791038   7.36360284

  7.35359269  16.08342325   9.16533046  20.35507374  12.63160058

 22.83356472   9.66291461   4.18055603  13.70368584  11.4533557

  4.16940565  10.31271413  23.06786868  17.80464565  14.53070132

 15.19656684  14.22969609   7.54691167  13.47210324  15.00625898

 19.28532444  20.7319878   19.70408833  18.21640853   8.50112687

  9.8493781    9.51425763   9.73270043  18.13782015  15.41731544

  5.07416787  12.20575251  14.05507493  10.6699926    7.16006245

 11.80728836  24.79748121  10.40809168  24.05228404  18.44737314

 20.80572631   9.45424805  17.00481708   5.78634105   5.10594849]

type 'numpy.ndarray'

5、回归问题的评价测度

(1) 评价测度

对于分类问题，评价测度是准确率，但这种方法不适用于回归问题。我们使用针对连续数值的评价测度(evaluation metrics)。

这里介绍3种常用的针对线性回归的测度。

1)平均绝对误差(Mean Absolute Error, MAE)

(2)均方误差(Mean Squared Error, MSE)

(3)均方根误差(Root Mean Squared Error, RMSE)

这里我使用RMES。

[python] view plain copy

pre name="code" class="python"#计算Sales预测的RMSE

print type(y_pred),type(y_test)

print len(y_pred),len(y_test)

print y_pred.shape,y_test.shape

from sklearn import metrics

import numpy as np

sum_mean=0

for i in range(len(y_pred)):

sum_mean+=(y_pred[i]-y_test.values[i])**2

sum_erro=np.sqrt(sum_mean/50)

# calculate RMSE by hand

print "RMSE by hand:",sum_erro

最后的结果如下：

type 'numpy.ndarray' class 'pandas.core.series.Series'

50 50

(50,) (50,)

RMSE by hand: 1.42998147691

（2）做ROC曲线

[python] view plain copy

import matplotlib.pyplot as plt

plt.figure()

plt.plot(range(len(y_pred)),y_pred,'b',label="predict")

plt.plot(range(len(y_pred)),y_test,'r',label="test")

plt.legend(loc="upper right") #显示图中的标签

plt.xlabel("the number of sales")

plt.ylabel('value of sales')

plt.show()

显示结果如下：（红色的线是真实的值曲线，蓝色的是预测值曲线）

直到这里整个的一次多元线性回归的预测就结束了。

6、改进特征的选择

在之前展示的数据中，我们看到Newspaper和销量之间的线性关系竟是负关系（不用惊讶，这是随机特征抽样的结果。换一批抽样的数据就可能为正了），现在我们移除这个特征，看看线性回归预测的结果的RMSE如何？

依然使用我上面的代码，但只需修改下面代码中的一句即可：

[python] view plain copy

#create a python list of feature names

feature_cols = ['TV', 'Radio', 'Newspaper']

# use the list to select a subset of the original DataFrame

X = data[feature_cols]

# equivalent command to do this in one line

#X = data[['TV', 'Radio', 'Newspaper']]#只需修改这里即可pre name="code" class="python" style="font-size: 15px; line-height: 35px;"X = data[['TV', 'Radio']]  #去掉newspaper其他的代码不变

# print the first 5 rowsprint X.head()# check the type and shape of Xprint type(X)print X.shape

最后的到的系数与测度如下：

LinearRegression(copy_X=True, fit_intercept=True, normalize=False)

2.81843904823

[ 0.04588771  0.18721008]

RMSE by hand: 1.28208957507

然后再次使用ROC曲线来观测曲线的整体情况。我们在将Newspaper这个特征移除之后，得到RMSE变小了，说明Newspaper特征可能不适合作为预测销量的特征，于是，我们得到了新的模型。我们还可以通过不同的特征组合得到新的模型，看看最终的误差是如何的。

备注：

之前我提到了这种错误：

注：上面的结果是由train_test_spilit()得到的，但是我不知道为什么我的版本的sklearn包中居然报错：

ImportError                               Traceback (most recent call last)ipython-input-182-3eee51fcba5a in module()      1 ###构造训练集和测试集---- 2 from sklearn.cross_validation import train_test_split      3 #import sklearn.cross_validation      4 X_train,X_test, y_train, y_test = train_test_split(X, y, random_state=1)      5 # default split is 75% for training and 25% for testingImportError: cannot import name train_test_split

处理方法：1、我后来重新安装sklearn包。再一次调用时就没有错误了。

2、自己写函数来认为的随机构造训练集和测试集。(这个代码我会在最后附上。)

这里我给出我自己写的函数：

Q4: python函数有哪些

1、print()函数：打印字符串;

2、raw_input()函数：从用户键盘捕获字符;

3、len()函数：计算字符长度;

4、format()函数：实现格式化输出;

5、type()函数：查询对象的类型;

6、int()函数、float()函数、str()函数等：类型的转化函数;

7、id()函数：获取对象的内存地址;

8、help()函数：Python的帮助函数;

9、s.islower()函数：判断字符小写;

10、s.sppace()函数：判断是否为空格;

11、str.replace()函数：替换字符;

12、import()函数：引进库;

13、math.sin()函数：sin()函数;

14、math.pow()函数：计算次方函数;

15、os.getcwd()函数：获取当前工作目录;

16、listdir()函数：显示当前目录下的文件;

17、time.sleep()函数：停止一段时间;

18、random.randint()函数：产生随机数;

19、range()函数：返回一个列表，打印从1到100;

20、file.read()函数：读取文件返回字符串;

21、file.readlines()函数：读取文件返回列表;

22、file.readline()函数：读取一行文件并返回字符串;

23、split()函数：用什么来间隔字符串;

24、isalnum()函数：判断是否为有效数字或字符;

25、isalpha()函数：判断是否全为字符;

26、isdigit()函数：判断是否全为数字;

27、 lower()函数：将数据改成小写;

28、upper()函数：将数据改成大写;

29、startswith(s)函数：判断字符串是否以s开始的;

30、endwith(s)函数：判断字符串是否以s结尾的;

31、file.write()函数：写入函数;

32、file.writeline()函数：写入文件;

33、abs()函数：得到某数的绝对值;

34、file.sort()函数：对书数据排序;

35、tuple()函数：创建一个元组;

36、find()函数：查找 返回的是索引;

37、dict()函数：创建字典;

38、clear()函数：清楚字典中的所有项;

39、copy()函数：复制一个字典，会修改所有的字典;

40、 get()函数：查询字典中的元素。

…………

Q5: 如何用python实现含有虚拟自变量的回归

利用python进行线性回归

理解什么是线性回归

线性回归也被称为最小二乘法回归（Linear Regression, also called Ordinary Least-Squares (OLS) Regression）。它的数学模型是这样的：

y = a+ b* x＋e

其中，a 被称为常数项或截距；b 被称为模型的回归系数或斜率；e 为误差项。a 和 b 是模型的参数。

当然，模型的参数只能从样本数据中估计出来：

y'= a' + b'* x

我们的目标是选择合适的参数，让这一线性模型最好地拟合观测值。拟合程度越高，模型越好。

那么，接下来的问题就是，我们如何判断拟合的质量呢？

这一线性模型可以用二维平面上的一条直线来表示，被称为回归线。

模型的拟合程度越高，也即意味着样本点围绕回归线越紧密。

如何计算样本点与回归线之间的紧密程度呢？

高斯和勒让德找到的方法是：被选择的参数，应该使算出来的回归线与观测值之差的平房和最小。用函数表示为：

这被称为最小二乘法。最小二乘法的原理是这样的：当预测值和实际值距离的平方和最小时，就选定模型中的两个参数（a 和 b）。这一模型并不一定反映解释变量和反应变量真实的关系。但它的计算成本低；相比复杂模型更容易解释。

模型估计出来后，我们要回答的问题是：

我们的模型拟合程度如何？或者说，这个模型对因变量的解释力如何？（R2）

整个模型是否能显著预测因变量的变化？（F 检验）

每个自变量是否能显著预测因变量的变化？（t 检验）

首先回答第一个问题。为了评估模型的拟合程度如何，我们必须有一个可以比较的基线模型。

如果让你预测一个人的体重是多少？在没有任何额外信息的情况下，你可能会用平均值来预测，尽管会存在一定误差，但总比瞎猜好。

现在，如果你知道他的身高信息，你的预测值肯定与平均值不一样。额外信息相比平均值更能准确地预测被预测的变量的能力，就代表模型的解释力大小。

上图中，SSA 代表由自变量 x 引起的 y 的离差平方和，即回归平方和，代表回归模型的解释力；SSE 代表由随机因素引起的 y 的离差平方和，即剩余平方和，代表回归模型未能解释的部分；SST 为总的离差平方和，即我们仅凭 y 的平均值去估计 y 时所产生的误差。

用模型能够解释的变异除以总的变异就是模型的拟合程度：

R2=SSA/SST=1-SSE

R2（R 的平方）也被称为决定系数或判定系数。

第二个问题，我们的模型是否显著预测了 y 的变化？

假设 y 与 x 的线性关系不明显，那么 SSA 相对 SSE 占有较大的比例的概率则越小。换句话说，在 y 与 x 无线性关系的前提下，SSA 相对 SSE 的占比越高的概率是越小的，这会呈现一定的概率分布。统计学家告诉我们它满足 F 分布，就像这样：

如果 SSA 相对 SSE 占比较大的情况出现了，比如根据 F 分布，这个值出现的概率小于 5%。那么，我们最好是拒绝 y 与 x 线性关系不显著的原始假设，认为二者存在显著的线性关系较为合适。

第三个问题，每个自变量是否能显著预测因变量的变化？换句话说，回归系数是否显著？

回归系数的显著性检验是围绕回归系数的抽样分布（t 分布）来进行的，推断过程类似于整个模型的检验过程，不赘言。

实际上，对于只有一个自变量的一元线性模型，模型的显著性检验和回归系数的检验是一致的，但对于多元线性模型来说，二者就不能等价了。

利用 statsmodels 进行最小二乘回归

＃导入相应模块

In [1]: import numpy as np

In [2]: import pandas as pd

In [3]: import statsmodels.api as sm

＃将数据导入 pandas 的 dataframe 对象，第一列（年份）作为行标签

In [4]: df=pd.read_csv('/Users/xiangzhendong/Downloads/vincentarelbundock-Rdatasets-1218370/csv/datasets/longley.csv', index_col=0)

＃查看头部数据

In [5]: df.head()

Out[5]:

GNP.deflator      GNP  Unemployed  Armed.Forces  Population  Year  \

1947          83.0  234.289       235.6         159.0     107.608  1947

1948          88.5  259.426       232.5         145.6     108.632  1948

1949          88.2  258.054       368.2         161.6     109.773  1949

1950          89.5  284.599       335.1         165.0     110.929  1950

1951          96.2  328.975       209.9         309.9     112.075  1951

Employed

1947    60.323

1948    61.122

1949    60.171

1950    61.187

1951    63.221

＃设置预测变量和结果变量，用 GNP 预测 Employed

In [6]: y=df.Employed ＃结果变量

In [7]: X=df.GNP ＃预测变量

＃为模型增加常数项，即回归线在 y 轴上的截距

In [8]: X=sm.add_constant(X)

＃执行最小二乘回归，X 可以是 numpy array 或 pandas dataframe（行数等于数据点个数，列数为预测变量个数），y 可以是一维数组（numpy array）或 pandas series

In [10]: est=sm.OLS(y,X)

使用 OLS 对象的 fit() 方法来进行模型拟合

In [11]: est=est.fit()

＃查看模型拟合的结果

In [12]: est.summary()

Out[12]:

＃查看最终模型的参数

In [13]: est.params

Out[13]:

const    51.843590

GNP       0.034752

dtype: float64

＃选择 100 个从最小值到最大值平均分布（equally spaced）的数据点

In [14]: X_prime=np.linspace(X.GNP.min(), X.GNP.max(),100)[:,np.newaxis]

In [15]: X_prime=sm.add_constant(X_prime)

＃计算预测值

In [16]: y_hat=est.predict(X_prime)

In [17]: plt.scatter(X.GNP, y, alpha=0.3) ＃画出原始数据

＃分别给 x 轴和 y 轴命名

In [18]: plt.xlabel("Gross National Product")

In [19]: plt.ylabel("Total Employment")

In [20]: plt.plot(X_prime[:,1], y_hat, 'r', alpha=0.9) ＃添加回归线，红色

多元线性回归（预测变量不止一个）

我们用一条直线来描述一元线性模型中预测变量和结果变量的关系，而在多元回归中，我们将用一个多维（p）空间来拟合多个预测变量。下面表现了两个预测变量的三维图形：商品的销量以及在电视和广播两种不同媒介的广告预算。

数学模型是：

Sales = beta_0 + beta_1＊TV + beta_2＊Radio

图中，白色的数据点是平面上的点，黑色的数据点事平面下的点。平面的颜色是由对应的商品销量的高低决定的，高是红色，低是蓝色。

利用 statsmodels 进行多元线性回归

In [1]: import pandas as pd

In [2]: import numpy as np

In [3]: import statsmodels.api as sm

In [4]: df_adv=pd.read_csv('g.csv',index_col=0)

In [6]: X=df_adv[['TV','Radio']]

In [7]: y=df_adv['Sales']

In [8]: df_adv.head()

Out[8]:

TV  Radio  Newspaper  Sales

1  230.1   37.8       69.2   22.1

2   44.5   39.3       45.1   10.4

3   17.2   45.9       69.3    9.3

4  151.5   41.3       58.5   18.5

5  180.8   10.8       58.4   12.9

In [9]: X=sm.add_constant(X)

In [10]: est=sm.OLS(y,X).fit()

In [11]: est.summary()

Out[11]:

你也可以使用 statsmodels 的 formula 模块来建立多元回归模型

In [12]: import statsmodels.formula.api as smf

In [13]: est=smf.ols(formula='Sales ~ TV + Radio',data=df_adv).fit()

处理分类变量

性别或地域都属于分类变量。

In [15]: df= pd.read_csv('httd.edu/~tibs/ElemStatLearn/datasets/SAheart.data', index_col=0)

In [16]: X=df.copy()

利用 dataframe 的 pop 方法将 chd 列单独提取出来

In [17]: y=X.pop('chd')

In [18]: df.head()

Out[18]:

sbp  tobacco   ldl  adiposity  famhist  typea  obesity  alcohol  \

row.names

1          160    12.00  5.73      23.11  Present     49    25.30    97.20

2          144     0.01  4.41      28.61   Absent     55    28.87     2.06

3          118     0.08  3.48      32.28  Present     52    29.14     3.81

4          170     7.50  6.41      38.03  Present     51    31.99    24.26

5          134    13.60  3.50      27.78  Present     60    25.99    57.34

age  chd

row.names

1           52    1

2           63    1

3           46    0

4           58    1

5           49    1

In [19]: y.groupby(X.famhist).mean()

Out[19]:

famhist

Absent     0.237037

Present    0.500000

Name: chd, dtype: float64

In [20]: import statsmodels.formula.api as smf

In [21]: df['famhist_ord']=pd.Categorical(df.famhist).labels

In [22]: est=smf.ols(formula="chd ~ famhist_ord", data=df).fit()

分类变量的编码方式有许多，其中一种编码方式是虚拟变量编码（dummy－encoding），就是把一个 k 个水平的分类变量编码成 k－1 个二分变量。在 statsmodels 中使用 C 函数实现。

In [24]: est=smf.ols(formula="chd ~ C(famhist)", data=df).fit()

In [26]: est.summary()

Out[26]:

处理交互作用

随着教育年限（education）的增长，薪酬 (wage) 会增加吗？这种影响对男性和女性而言是一样的吗？

这里的问题就涉及性别与教育年限的交互作用。

换言之，教育年限对薪酬的影响是男女有别的。

＃导入相关模块

In [1]: import pandas as pd

In [2]: import numpy as np

In [4]: import statsmodels.api as sm

＃导入数据，存入 dataframe 对象

In [5]: df=pd.read_csv('/Users/xiangzhendong/Downloads/pydatafromweb/wages.csv')

In [6]: df[['Wage','Education','Sex']].tail()

Out[6]:

Wage  Education  Sex

529  11.36         18    0

530   6.10         12    1

531  23.25         17    1

532  19.88         12    0

533  15.38         16    0

由于性别是一个二分变量，我们可以绘制两条回归线，一条是 sex＝0（男性），一条是 sex＝1（女性）

＃绘制散点图

In [7]: plt.scatter(df.Education,df.Wage, alpha=0.3)

In [9]: plt.xlabel('education')

In [10]: plt.ylabel('wage')

＃linspace 的作用是生成从最小到最大的均匀分布的 n 个数

In [17]: education_linspace=np.linspace(df.Education.min(), df.Education.max(),100)

In [12]: import statsmodels.formula.api as smf

In [13]: est=smf.ols(formula='Wage ~ Education + Sex', data=df).fit()

In [18]: plt.plot(education_linspace, est.params[0]+est.params[1]education_linspace+est.params[2]0, 'r')

In [19]: plt.plot(education_linspace, est.params[0]+est.params[1]education_linspace+est.params[2]1, 'g')

以上两条线是平行的。这是因为分类变量只影响回归线的截距，不影响斜率。

接下来我们可以为回归模型增加交互项来探索交互效应。也就是说，对于两个类别，回归线的斜率是不一样的。

In [32]: plt.scatter(df.Education,df.Wage, alpha=0.3)

In [33]: plt.xlabel('education')

In [34]: plt.ylabel('wage')

＃使用＊代表我们的回归模型中除了交互效应，也包括两个变量的主效应；如果只想看交互效应，可以用：代替，但通常不会只看交互效应

In [35]: est=smf.ols(formula='Wage ~ Sex*Education', data=df).fit()

In [36]: plt.plot(education_linspace, est.params[0]+est.params[1]0+est.params[2]education_linspace+est.params[3]0education_linspace, 'r')

In [37]: plt.plot(education_linspace, est.params[0]+est.params[1]1+est.params[2]education_linspace+est.params[3]1education_linspace, 'g')

参考资料：

DataRobot | Ordinary Least Squares in Python

DataRoboe | Multiple Regression using Statsmodels

AnalyticsVidhya | 7 Types of Regression Techniques you should know!