专栏名称: 数据派THU

本订阅号是“THU数据派”的姊妹账号，致力于传播大数据价值、培养数据思维。

【独家】手把手教线性回归分析（附R语言实例）

数据派THU · 公众号 · 大数据 · 2017-06-11 20:18

正文

本文长度为8619字，建议阅读15分钟

本文为你介绍线性回归分析。

通常在现实应用中，我们需要去理解一个变量是如何被一些其他变量所决定的。

回答这样的问题，需要我们去建立一个模型。一个模型就是一个公式之中，一个因变量（dependent variable）(需要预测的值)会随着一个或多个数值型的自变量（independent variable）（预测变量）而改变的。我们能够构建的最简单的模型之一就是线性模型，我们可以假设因变量和自变量间是线性的关系。回归分方法可用于预测数值型数据以及量化预测结果与其预测变量之间关系的大小及强度。本文将介绍如何将回归方法应用到你自己的数据中，主要介绍学习内容：

用线性回归方法来拟合数据方程的基本统计原则和它们如何描述数据元素之间的关系。
如何使用R准备数据进行回归分析，定义一个线性方程并估计回归模型。

一、理解回归

回归主要关注确定一个唯一的因变量（dependent variable）(需要预测的值)和一个或多个数值型的自变量（independent variable）（预测变量）之间的关系。我们首先假设因变量和自变量之间的关系遵循一条直线，即线性关系。

你可能还记得数学中是以类似于Y=aX + b的斜截式来定义直线的，其中，y是因变量，x是自变量。在这个公式中，斜率（slope）a表示每增加一个单位的x，直接会上升的高度；变量b表示X=0时y的值，它称为截距，因为它指定了直线穿过y轴时的位置。

回归方程使用类似于斜截式的形式对数据建立模型。该机器的工作就是确定a和b的值，从而使指定的直线最适合用来反映所提供的x值和y值之间的关系，这可能不是完美的匹配，所以该机器也需要有一些方法来量化误差范围，很快我们就会讨论这个问题。

回归分析通常用来对数据元素之间的复杂关系建立模型，用来估计一种处理方法对结果的影响和推断未来。一些具体应用案例包括：

根据种群和个体测得的特征，研究他们之间如何不同（差异性），从而用于不同领域的科学研究，如经济学、社会学、心理学、物理学和生态学；
量化事件及其相应的因果关系，比如可应用于药物临床试验、工程安全检测、销售研究等。
给定已知的规则，确定可用来预测未来行为的模型，比如用来预测保险赔偿、自然灾害的损失、选举的结果和犯罪率等。

回归方法也可用于假设检验，其中包括数据是否能够表明原假设更可能是真还是假。回归模型对关系强度和一致性的估计提供了信息用于评估结果是否是由于偶然性造成的。回归分析是大量方法的一个综合体，几乎可以应用于所有的机器学习任务。如果被限制只能选择单一的分析方法，那么回归方法将是一个不错的选择。

本文只关注最基本的回归模型，即那些使用直线回归的模型，这叫做线性回归（linearregression）。如果只有一个单一的自变量，那就是所谓的简单线性回归（simple linear regression），否则，称为多元回归（multiple regression），这两个模型都假定因变量是连续的。对其他类型的因变量，即使是分类任务，使用回归方法也是可能的。逻辑回归（logistic regression）可以用来对二元分类的结果建模；泊松分布（Possion regression）可以用来对整型的计数数据建模。相同的基本原则适用于所有的回归方法，所以一旦理解了线性情况下的回归方法，就可以研究其他的回归方法。

二、简单线性回归

让我们从基础开始。记得高中时学过的直线方程吗？

Y=aX + b

a就是斜率，b就是y轴截距。简单而言，线性回归就是一系列技术用于找出拟合一系列数据点的直线。这也可以被认为是从数据之中反推出一个公式。我们会从最基础的一些规则开始，慢慢增加数学复杂度，增进对这个概念了解的深入程度。但是在此之前，也许你会很好奇这里的a和b的值分别是多少。接下来，我们通过一个例子，使用软件R来为我们计算，我们的数据来源于一组真实的关于儿童的身高和年龄，记录的数据。首先我们先直观地显示年龄与身高之间的关系，画出一张散点图，以年龄age为横坐标，身高height为纵坐标，R的代码如下：

> age=18:29 #年龄从18到29岁

> height=c(76.1,77,78.1,78.2,78.8,79.7,79.9,81.1,81.2,81.8,82.8,83.5)

> plot(age,height,main = "身高与年龄散点图")

散点图结果如图1所示。

图1 身高与年龄散点图

从图中可以观察到，年龄与身高基本在一条直线附近，可以认为两者具有线性关系，接下来建立回归模型，R代码如下：

> lm.reg #建立回归方程

> lm.reg

> abline(lm.reg) #画出拟合的线性回归线

产生以下的输出：

Call:

lm(formula = height ~ age)

cients:

(Intercept) age

64.928 0.635

图2 身高与年龄拟合直线

我们可以看到两个数值，“截距”和“斜率”。无论我们用什么软件来做线性回归（本文中的例子统一采用R语言），它都会用某种形式来报告这两个数值。截距就是我们的公式中的b，斜率就是Y和自变量之间的倾斜程度。

总结起来，我们有一个数据集（观测值）和一个模型（我们猜测可以拟合数据的一个公式），我们还要去找出模型的参数（我们的最佳拟合模型中的参数a和b），这样，模型就可以“最佳”拟合数据了。我们希望用我们的数据来找出一个公式的参数，这样，这个公式也能够“最佳”拟合数据了。

1. 用模型来做预测

一旦你目测出最佳拟合直线并且读出a和b，也许你大概会说大意是这样的话：“这些数据遵循一个形式为Y = aX + b 的线性方程，其中a（斜率）= 某个数，b（y轴截距）= 另外某个数”。你也许记得，这个等式并非是确切的表述，因为很有可能你的数据并不是所有的都在一条完美的直线上，所以数据点之间可能会有不同程度的误差。你的目测是主观地尝试把一些直觉上的总的误差给降到最低。你做的事情就是直觉上的“线性回归”。你按照“我看起来很顺眼”的算法来估计a和b。我们会以这个直觉的行为为开端，并迅速带入一些重量级的机械，使得我们能够解决一些相当复杂的问题。

在这个节点，你的实验室练习也许会要求你为不在你的观测值集合以内的，某个给定数值的X，给出Y的估值。然后你就会用上面的等式，比如说a是2.1，b是0.3的等式Y = 2.1 X + 0.3作为你的模型，你把X输入进去，就会得到一个Y。这时候你就是在用你的模型去预测一个值，换句话说，你正在陈述这样的事实：我在实验之中并没有用这个X值，并且我的数据里也没有它，但是我想要知道这个X值是怎样投射到Y轴上的。你也许会想要能够说出：“我的误差会是某个数，所以我相信实际上的值会在[Y-误差,Y+误差]之间”。在这样的情况下，我们把变量X叫做“预测变量”，而Y的值是基于X的一个值来预测的，所以变量Y是“反应”。

2. 一个概念：总误差

我们想要创造一个有关误差，或者说我们的直线所给出的Y值与我们数据集中的真实Y值之间的差异的简单公式。除非我们的直线正好穿过一个特定的点，否则，这个误差是非零的。它可能是正值，也可能是负值。我们可以取这个误差的平方（或者我们也可以取它的绝对值），然后我们把每个点的误差项相加起来，然后得到直线和这个数据集的总误差。在同一个实验的不同的样例集合中，我们会得到一个不同的数据集，很有可能一条不同的直线，并且几乎可以肯定一个不同的总误差。我们所用的误差的平方值是一个非常常用的总误差形式，它就是“方差”。它用同样的方式来处理正值的误差和负值的误差，所以方差总是正值的。从现在开始我们会使用“方差”作为我们误差的代表。总的来说，回归是我们所使用的任意的数据，最小化方差，来估测模型系数的手段。统计软件用多变量微积分这些专业技术来最小化误差，并且给我们提供系数的估测值。对于回归方程回归系数的检验，检验一般用方差分析或t检验，两者的检验结果是等价的。方差分析主要是针对整个模型的，而t检验是关于回归系数的。

对于上例中的回归方程，我们对模型进行检验，方差分析的R代码如下：

> anova(lm.reg) #模型方差分析

产生以下的输出：

Analysis of Variance Table

Response: height

Df Sum Sq Mean Sq F value Pr(>F)

age 1 57.655 57.655 879.99 4.428e-11 ***

Residuals 10 0.655 0.066

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

由于P<0.05，于是在α=0.05水平下，本例的回归系数有统计学意义，身高和年龄存在直线回归关系。

同理，对于上例中的回归方程，我们对模型进行回归系数的t检验，t检验的R代码如下：

> summary(lm.reg) #回归系数的t检验

产生以下的输出：

Call:

lm(formula = height ~ age)

Residuals:

Min 1Q Median 3Q Max

-0.27238 -0.24248 -0.02762 0.16014 0.47238

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 64.9283 0.5084 127.71 < 2e-16 ***

age 0.6350 0.0214 29.66 4.43e-11 ***

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.256 on 10 degrees of freedom

Multiple R-squared: 0.9888,Adjusted R-squared: 0.9876

F-statistic: 880 on 1 and 10 DF, p-value: 4.428e-11

同方差分析，由于P<0.05，于是在α=0.05水平下，本例的回归系数有统计学意义，身高和年龄存在回归关系。

三、多元线性回归

很自然地，如果Y依赖于多于一个变量的时候，会发生什么。这时候，数学上的普遍性的力量就显现了。这个简单的原理在多维的情况下同样适用。不仅仅是两个或者三个，还有更多更大的数值。如果我们想要构建现实数据的模型的话，二十个，三十个，甚至百来个自变量都毫无问题。但现在我们就看看Y，作为一个含有多个自变量的函数，例如含有来个自变量X1和X2的形式：

Y = a₀ + a₁X₁ + a₂X₂

这里的a₀就是截距项，a₁，a₂就是自变量X₁，X₂ 的系数。为了看看具有多个潜在的自变量的真实数据集，我们会在下一步使用具体的数据——医疗费用的数据。

多元线性回归优缺点
优点	缺点
迄今为止，它是数值型数据建模最常用的方法	对数据做出了很强的假设
可适用于几乎所有的数据	该模型的形式必须由使用者事先指定
提供了特征（变量）与结果之间关系的强度和大小的估计	不能很好地处理缺失数据
	只能处理数值特征，所以分类数据需要额外的处理
	需要一些统计知识来理解模型

四、线性回归数据探索：医疗费用

医疗费用很难估计，因为花费最高的情况是罕见的而且似乎是随机的。但是有些情况对于部分特定的群体还是比较普遍存在的。例如，吸烟者比不吸烟者得肺癌的可能性更大，肥胖的人更有可能得心脏病。此分析的目的是利用病人的数据，来预测这部分群体的平均医疗费用。这些估计可以用来创建一个精算表，根据预期的治疗费用来设定年度保费价格是高一点还是低一点。我们想要去探索这些数据，尝试从中获得一些对于建立线性回归模型有用处的见解。

我们会依据以下步骤完成，将来，我们还会对其他数据集座同样的事情。

1) 收集/观察数据；

2) 探索和准备数据；

3) 基于数据训练模型；

4) 评估模型的性能；

5) 提高模型的性能。

第1步——收集/观察数据

为了便于分析，我们使用一个模拟数据集，该数据集包含了美国病人的医疗费用。而本文创建的这些数据使用了来自美国人口普查局（U.S. Census Bureau）的人口统计资料，因此可以大致反映现实世界的情况。

注：如果你想一起学习这个例子，那么你需要从Packt出版社的网站（https://github.com/stedy/Machine-Learning-with-R-datasets/find/master）下载insurance.csv文件，并将该文件保存到R的工作文件夹中。

该文件（insurance.csv）包含1338个案例，即目前已经登记过的保险计划受益者以及表示病人特点和历年计划计入的总的医疗费用的特征。这些特征是：

age: 这是一个整数，表示主要受益者的年龄（不包括超过64岁的人，因为他们一般由政府支付）。
sex: 这是保单持有人的性别，要么是male，要么是female。
bmi: 这是身体质量指数（Body Mass Index,BMI），它提供了一个判断人的体重相对于身高是过重还是偏轻的方法，BMI指数等于体重（公斤）除以身高（米）的平方。一个理想的BMI指数在18.5~24.9的范围内。
children: 这是一个整数，表示保险计划中所包括的孩子/受抚养者的数量。
smoker: 根据被保险人是否吸烟判断yes或者no。
region: 根据受益人在美国的居住地，分为4个地理区域：northeast、southeast、southwest和northwest。

如何将这些变量与已结算的医疗费用联系在一起是非常重要的。例如，我们可能认为老年人和吸烟者在大额医疗费用上是有较高的风险。与许多其他的方法不同，在回归分析中，特征之间的关系通常由使用者指定而不是自动检测出来。

第2步——探索和准备数据

在R中，我们将使用read.csv()函数来加载用于分析的数据。我们可以使用stringAsFactors = TRUE，因为将名义变量转换成因子变量是恰当的：

> insurance

函数str()确认该数据转换了我们之前所期望的形式：

> str(insurance)

产生以下的输出：

'data.frame':1338 obs. of 7 variables:

$ age : int 19 18 28 33 32 31 46 37 37 60 ...

$ sex : Factor w/ 2 levels "female","male": 1 2 2 2 2 1 1 1 2 1 ...

$ bmi : num 27.9 33.8 33 22.7 28.9 ...

$ children: int 0 1 3 0 0 0 1 3 2 0 ...

$ smoker : Factor w/ 2 levels "no","yes": 2 1 1 1 1 1 1 1 1 1 ...

$ region : Factor w/ 4 levels "northeast","northwest",..: 4 3 3 2 2 3 3 2 1 2 ...

$ charges : num 16885 1726 4449 21984 3867 ...

既然因变量是changes，那么让我们一起来看一下它是如何分布的：

> summary(insurance$charges)

产生以下的输出：

Min. 1st Qu. Median Mean 3rd Qu. Max.

1122 4740 9382 13270 16640 63770

因为平均数远大于中位数，表明保险费用的分布是右偏的，我们可以用直方图证实这一点。

> hist(insurance$charges)

图3 charges直方图

在我们的数据中，绝大多数的个人每年的费用都在0~15000美元，尽管分布的尾部经过直方图的峰部后延伸得很远。即将面临的另一个问题就是回归模型需要每一个特征都是数值型的，而在我们的数据框中，我们有3个因子类型的特征。很快，我们会看到R中的线性回归函数如何处理我们的变量。

1.探索特征之间的关系——相关系数矩阵

在使用回归模型拟合数据之前，有必要确定自变量与因变量之间以及自变量之间是如何相关的。相关系数矩阵（correlation matrix）提供了这些关系的快速概览。给定一组变量，它可以为每一对变量之间的关系提供一个相关系数。

为insurance数据框中的4个数值型变量创建一个相关系数矩阵，可以使用cor()命令：

> cor(insurance[c("age","bmi","children","charges")])

产生以下的输出：

age bmi children charges

age 1.0000000 0.1092719 0.04246900 0.29900819

bmi 0.1092719 1.0000000 0.01275890 0.19834097

children 0.0424690 0.0127589 1.00000000 0.06799823

charges 0.2990082 0.1983410 0.06799823 1.00000000

该矩阵中中的相关系数不是强相关的，但还是存在一些显著的关联。例如，age和bmi显示出中度相关，这意味着随着年龄（age）的增长，身体质量指数（bmi）也会增加。此外，age和charges，bmi和charges，以及children和charges也都呈现处中度相关。当我们建立最终的回归模型时，我们会尽量更加清晰地梳理出这些关系。

2.可视化特征之间的关系——散点图矩阵

或许通过使用散点图，可视化特征之间的关系更有帮助。虽然我们可以为每个可能的关系创建一个散点图，但对于大量的特征，这样做可能会变得比较繁琐。

另一种方法就是创建一个散点图矩阵（scatterplot matrix)，就是简单地将一个散点图集合排列在网格中，里边包含着相互紧邻在一起的多种因素的图表。它显示了每个因素相互之间的关系。斜对角线上的图并不符合这个形式。为何不符合呢？在这个语境下，这意味着找到某个事物和自身的关系，而我们正在尝试确定某些变量对于另一个变量的影响。默认的R中提供了函数pairs()，该函数产生散点图矩阵提供了基本的功能。对医疗费用数据之中的四个变量的散点图矩阵如下图所示。R代码如下：

pairs(insurance[c("age","bmi","children","charges")])

图4 散点图矩阵

与相关系数矩阵一样，每个行与列的交叉点所在的散点图表示其所在的行与列的两个变量的相关关系。由于对角线上方和下方的x轴和y轴是交换的，所以对角线上方的图和下方的图是互为转置的。

你注意到这些散点图中的一些图案了吗？尽管有一些看上去像是随机密布的点，但还是有一些似乎呈现了某种趋势。age和charges之间的关系呈现出几条相对的直线，而bmi和charges的散点图构成了两个不同的群体。

如果我们对散点图添加更多的信息，那么它就会更加有用。一个改进后的散点图矩阵可以用psych包中的pairs.panels()函数来创建。R中如果你还没有安装这个包，那么可以输入install.packages("psych")命令将其安装到你的系统中，并使用library(psych)命令加载它。R代码及散点图矩阵如下：

pairs.panels(insurance[c("age","bmi","children","charges")])

图5 散点图矩阵

在对角线的上方，散点图被相关系数矩阵所取代。在对角线上，直方图描绘了每个特征的数值分布。最后，对角线下方的散点图带有额外的可视化信息。

每个散点图中呈椭圆形的对象称为相关椭圆（correlation ellipse），它提供了一种变量之间是如何密切相关的可视化信息。位于椭圆中心的点表示x轴变量的均值和y轴变量的均值所确定的点。两个变量的相关性由椭圆的形状所表示，椭圆越被拉伸，其相关性越强。一个几乎类似于圆的完美的椭圆形，如bmi和children，表示一种非常弱的相关性。

散点图中绘制的曲线称为局部回归平滑（loess smooth），它表示x轴和y轴变量之间的一般关系。最好通过例子来理解。散点图中age和childr的曲线是一个倒置的U，峰值在中年附近，这意味着案例中年龄最大的人和年龄最小的人比年龄大约在中年附近的人拥有的孩子更少。因为这种趋势是非线性的，所以这一发现已经不能单独从相关性推断出来。另一方面，对于age和bmi，局部回归光滑是一条倾斜的逐步上升的线，这表明BMI会随着年龄（age）的增长而增加，从相关系数矩阵中我们也可推断出该结论。

第3步——基于数据训练模型

用R对数据拟合一个线性回归模型时，可以使用lm()函数。该函数包含在stats添加包中，当安装R时，该包已经被默认安装并在R启动时自动加载好。使用R拟合称为ins_model的线性回归模型，该模型将6个自变量与总的医疗费用联系在一起。代码如下：

ins_model

建立模型后，只需输入该模型对象的名称，就可以看到估计的a系数：

> ins_model

产生以下的输出：

Call:

lm(formula = charges ~ ., data = insurance)

Coefficients:

(Intercept) age sexmale

-11938.5 256.9 -131.3

bmi children smokeryes

339.2 475.5 23848.5

regionnorthwest regionsoutheast regionsouthwest

-353.0 -1035.0 -960.1

你可能注意到，在我们的模型公式中，我们仅指定了6个变量，但是输出时，除了截距项外，却输出了8个系数。之所以发生这种情况，是因为lm()函数自动将一种称为虚拟编码(dummy coding)的技术应用于模型所包含的每一个因子类型的变量中。当添加一个虚拟编码的变量到回归模型中时，一个类别总是被排除在外作为参照类别。然后，估计的系数就是相对于参照类别解释的。在我们的模型中，R自动保留sexfemale、smokerno和regionnortheast变量，使东北地区的女性非吸烟者作为参照组。因此，相对于女性来说，男性每年的医疗费用要少$131.30；吸烟者平均多花费$23848.50，远超过非吸烟者。此外，模型中另外3个地区的系数是负的，这意味着东北地区倾向于具有最高的平均医疗费用。

线性回归模型的结果是合乎逻辑的。高龄、吸烟和肥胖往往与其他健康问题联系在一起，而额外的家庭成员或者受抚养者可能会导致就诊次数增加和预防保健（比如接种疫苗、每年体检）费用的增加。然而，我们目前并不知道该模型对数据的拟合有多好？我们将在下一部分回答这个问题。

第4步——评估模型的性能

通过在R命令行输入ins_model，可以获得参数的估计值，它们告诉我们关于自变量是如何与因变量相关联的。但是它们根本没有告诉我们用该模型来拟合数据有多好。为了评估模型的性能，可以使用summary()命令来分析所存储的回归模型：

> summary(ins_model)

产生以下的输出：

Call:

lm(formula = charges ~ ., data = insurance)

Residuals:

Min 1Q Median 3Q Max

-11304.9 -2848.1 -982.1 1393.9 29992.8

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) -11938.5 987.8 -12.086 < 2e-16 ***

age 256.9 11.9 21.587 < 2e-16 ***

sexmale -131.3 332.9 -0.394 0.693348

bmi 339.2 28.6 11.860 < 2e-16 ***

children 475.5 137.8 3.451 0.000577 ***

smokeryes 23848.5 413.1 57.723 < 2e-16 ***

regionnorthwest -353.0 476.3 -0.741 0.458769

regionsoutheast -1035.0 478.7 -2.162 0.030782 *

regionsouthwest -960.0 477.9 -2.009 0.044765 *

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6062 on 1329 degrees of freedom

Multiple R-squared: 0.7509,Adjusted R-squared: 0.7494

F-statistic: 500.8 on 8 and 1329 DF, p-value: < 2.2e-16

开始时，summary()的输出可能看起来令人费解，但基本原理是很容易掌握的。与上述输出中用标签编号所表示的一样，该输出为评估模型的性能提供了3个关键的方面：

1） Residuals（残差）部分提供了预测误差的主要统计量；

2）星号（例如，***）表示模型中每个特征的预测能力；

3）多元R方值（也称为判定系数）提供度量模型性能的方式，即从整体上，模型能多大程度解释因变量的值。

给定前面3个性能指标，我们的模型表现得相当好。对于现实世界数据的回归模型，R方值相当低的情况并不少见，因此0.75的R方值实际上是相当不错的。考虑到医疗费用的性质，其中有些误差的大小是需要关注的，但并不令人吃惊。如下节所述，我们会以略微不同的方式来指定模型，从而提高模型的性能。

第5步——提高模型的性能

正如前面所提到的，回归模型通常会让使用者来选择特征和设定模型。因此，如果我们有关于一个特征是如何与结果相关的学科知识，我们就可以使用该信息来对模型进行设定，并可能提高模型的性能。

1. 模型的设定——添加非线性关系

在线性回归中，自变量和因变量之间的关系被假定为是线性的，然而这不一定是正确的。例如，对于所有的年龄值来讲，年龄对医疗费用的影响可能不是恒定的；对于最老的人群，治疗可能会过于昂贵。

2. 转换——将一个数值型变量转换为一个二进制指标

假设我们有一种预感，一个特征的影响不是累积的，而是当特征的取值达到一个给定的阈值后才产生影响。例如，对于在正常体重范围内的个人来说，BMI对医疗费用的影响可能为0，但是对于肥胖者（即BMI不低于30）来说，它可能与较高的费用密切相关。我们可以通过创建一个二进制指标变量来建立这种关系，即如果BMI大于等于30，那么设定为1，否则设定为0。

注：如果你在决定是否要包含一个变量时遇到困难，一种常见的做法就是包含它并检验其显著性水平。然后，如果该变量在统计上不显著，那么就有证据支持在将来排除该变量。

3. 模型的设定——加入相互作用的影响

到目前为止，我们只考虑了每个特征对结果的单独影响（贡献）。如果某些特征对因变量有综合影响，那么该怎么办呢？例如，吸烟和肥胖可能分别都有有害的影响，但是假设它们的共同影响可能会比它们每一个单独影响之和更糟糕是合理的。

当两个特征存在共同的影响时，这称为相互作用（interaction）。如果怀疑两个变量相互作用，那么可以通过在模型中添加它们的相互作用来检验这一假设，可以使用R中的公式语法来指定相互作用的影响。为了体现肥胖指标（bmi30）和吸烟指标（smoker）的相互作用，可以这样的形式写一个公式：charge~bmi30*smoker。

4. 全部放在一起——一个改进的回归模型

基于医疗费用如何与患者特点联系在一起的一点学科知识，我们采用一个我们认为更加精确的专用的回归公式。下面就总结一下我们的改进：

增加一个非线性年龄项
为肥胖创建一个指标
指定肥胖与吸烟之间的相互作用

我们将像之前一样使用lm()函数来训练模型，但是这一次，我们将添加新构造的变量和相互作用项：

> ins_model2

接下来，我们概述结果：

> summary(ins_model2)

产生以下的输出：

Call:

lm(formula = charges ~ age + age2 + children + bmi + sex + bmi30 *

smoker + region, data = insurance)

Residuals:

Min 1Q Median 3Q Max

-17296.4 -1656.0 -1263.3 -722.1 24160.2

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 134.2509 1362.7511 0.099 0.921539

age -32.6851 59.8242 -0.546 0.584915

age2 3.7316 0.7463 5.000 6.50e-07 ***

children 678.5612 105.8831 6.409 2.04e-10 ***

bmi 120.0196 34.2660 3.503 0.000476 ***

sexmale -496.8245 244.3659 -2.033 0.042240 *

bmi30 -1000.1403 422.8402 -2.365 0.018159 *

smokeryes 13404.6866 439.9491 30.469 < 2e-16 ***

regionnorthwest -279.2038 349.2746 -0.799 0.424212

regionsoutheast -828.5467 351.6352 -2.356 0.018604 *

regionsouthwest -1222.6437 350.5285 -3.488 0.000503 ***

bmi30:smokeryes 19810.7533 604.6567 32.764 < 2e-16 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4445 on 1326 degrees of freedom

Multiple R-squared: 0.8664,Adjusted R-squared: 0.8653

F-statistic: 781.7 on 11 and 1326 DF, p-value: < 2.2e-16

分析该模型的拟合统计量有助于确定我们的改变是否提高了回归模型的性能。相对于我们的第一个模型，R方值从0.75提高到约0.87，我们的模型现在能解释医疗费用变化的87%。此外，我们关于模型函数形式的理论似乎得到了验证，高阶项age2在在统计上是显著的，肥胖指标bmi30也是显著的。肥胖和吸烟之间的相互作用表明了一个巨大的影响，除了单独吸烟增加的超过$13404的费用外，肥胖的吸烟者每年要另外花费$19810，这可能表明吸烟会加剧（恶化）与肥胖有关的疾病。

作者简介

慕生鹏，数据派研究部志愿者。北京林业大学计算数学专业硕士在读学生。

日常喜欢长跑，健身等活动。对数据的分析、学习很感兴趣。日常会借助网络等资源，自助学习各类数据的分析方法。希望在数据的分析算法方面，不断地加强功底。

数据派研究部介绍

数据派研究部成立于2017年初，志于打造一流的结构化知识分享平台、活跃的数据科学爱好者社群，致力于传播数据思维、提升数据能力、探索数据价值、实现产学研结合！

研究部的逻辑在于知识结构化、实践出真知：梳理打造结构化基础知识网络；原创手把手教以及实践经验等文章；形成专业兴趣社群，交流学习、组队实践、追踪前沿。
兴趣组是研究部的核心，各组既遵循研究部整体的知识分享和实践项目规划，又会各具特色：
算法模型组：积极组队参加kaggle等比赛，原创手把手教系列文章；
调研分析组：通过专访等方式调研大数据的应用，探索数据产品之美；
系统平台组：追踪大数据&人工智能系统平台技术前沿，对话专家；
自然语言处理组：重于实践，积极参加比赛及策划各类文本分析项目；
制造业大数据组：秉工业强国之梦，产学研政结合，挖掘数据价值；
数据可视化组：将信息与艺术融合，探索数据之美，学用可视化讲故事；
网络爬虫组：爬取网络信息，配合其他各组开发创意项目。