第 12 章 线性回归

线性模型是数据分析中最常用的一种分析方法。最基础的往往最深刻。

12.1 从一个案例开始

这是一份1994年收集1379个对象关于收入、身高、教育水平等信息的数据集,数据在这里下载 Download wages.csv

首先,我们下载后导入数据

earn height sex race ed age
79571 73.89 male white 16 49
96397 66.23 female white 16 62
48711 63.77 female white 16 33
80478 63.22 female other 16 95
82089 63.08 female white 17 43
15313 64.53 female white 15 30

12.1.1 缺失值检查

一般情况下,拿到一份数据,首先要了解数据,知道每个变量的含义,

## [1] "earn"   "height" "sex"    "race"   "ed"    
## [6] "age"

同时检查数据是否有缺失值,这点很重要。在R中 NA(notavailable,不可用)表示缺失值, 比如可以这样检查是否有缺失值。

## # A tibble: 1 x 6
##   earn_na height_na sex_na race_na ed_na age_na
##     <int>     <int>  <int>   <int> <int>  <int>
## 1       0         0      0       0     0      0

程序员都是偷懒的,所以也可以写的简便一点。大家在学习的过程中,也会慢慢的发现tidyverse的函数很贴心,很周到。

## # A tibble: 1 x 6
##    earn height   sex  race    ed   age
##   <int>  <int> <int> <int> <int> <int>
## 1     0      0     0     0     0     0

当然,也可以用purrr::map()的方法。这部分我会在后面的章节中逐步介绍。

## # A tibble: 1 x 6
##    earn height   sex  race    ed   age
##   <int>  <int> <int> <int> <int> <int>
## 1     0      0     0     0     0     0

12.1.2 变量简单统计

然后探索下每个变量的分布。比如调研数据中男女的数量分别是多少?

## # A tibble: 2 x 2
##   sex        n
##   <chr>  <int>
## 1 female   859
## 2 male     520

男女这两组的身高均值分别是多少?收入的均值分别是多少?

## # A tibble: 2 x 4
##   sex        n mean_height mean_earn
##   <chr>  <int>       <dbl>     <dbl>
## 1 female   859        64.5    24246.
## 2 male     520        70.0    45993.

也有可以用可视化的方法,呈现男女收入的分布情况

大家可以自行探索其他变量的情况。现在提出几个问题,希望大家带着这些问题去探索:

  1. 长的越高的人挣钱越多?

  2. 是否男性就比女性挣的多?

  3. 影响收入最大的变量是哪个?

  4. 怎么判定我们建立的模型是不是很好?

12.2 线性回归模型

长的越高的人挣钱越多?

要回答这个问题,我们先介绍线性模型。顾名思义,就是认为\(x\)\(y\)之间有线性关系,数学上可以写为

\[ \begin{aligned} y &= \alpha + \beta x + \epsilon \\ \epsilon &\in \text{Normal}(\mu, \sigma) \end{aligned} \]

$ $ 代表误差项,它与\(x\) 无关,且服从正态分布。 建立线性模型,就是要估计这里的系数\(\hat\alpha\)\(\hat\beta\),即截距项和斜率项。常用的方法是最小二乘法(ordinary least squares (OLS) regression): 就是我们估算的\(\hat\alpha\)\(\hat\beta\), 要使得残差的平方和最小,即\(\sum_i(y_i - \hat y_i)^2\)或者叫\(\sum_i \epsilon_i^2\)最小。当然,数据量很大,手算是不现实的,我们借助R语言代码吧

12.3 使用lm() 函数

用R语言代码(建议大家先?lm看看帮助文档),

lm参数很多, 但很多我们都用不上,所以我们只关注其中重要的两个参数

lm(y ~ x, data) 是最常用的线性模型函数(lm是linear model的缩写)。参数解释说明

  • formula:指定回归模型的公式,对于简单的线性回归模型y ~ x.
  • ~ 符号:代表“预测”,可以读做“y由x预测”。有些学科不同的表述,比如下面都是可以的
    • response ~ explanatory
    • dependent ~ independent
    • outcome ~ predictors
  • data:代表数据框,数据框包含了响应变量和独立变量

在运行lm()之前,先画出身高和收入的散点图(记在我们想干什么,寻找身高和收入的关系)

等不及了,就运行代码吧

这里我们将earn作为响应变量,height为预测变量。lm()返回赋值给mod1. mod1现在是个什么东东呢? mod1是一个叫lm object或者叫的东西,

##  [1] "coefficients"  "residuals"     "effects"      
##  [4] "rank"          "fitted.values" "assign"       
##  [7] "qr"            "df.residual"   "xlevels"      
## [10] "call"          "terms"         "model"

我们打印看看,会发生什么

## 
## Call:
## lm(formula = earn ~ height, data = wages)
## 
## Coefficients:
## (Intercept)       height  
##     -126523         2387

这里有两部分信息。首先第一部分是我们建立的模型;第二部分是R给出了截距(\(\alpha = -126532\))和斜率(\(\beta = 2387\)). 也就是说我们建立的线性回归模型是 \[ \hat y = -126532 + 2387 \; x \]

12.4 模型的解释

建立一个lm模型是简单的,然而最重要的是,我们能解释这个模型。

mod1的解释:

  • 对于斜率\(\beta = 2387\)意味着,当一个人的身高是68英寸时,他的预期收入\(earn = -126532 + 2387 \times 68= 35806\) 美元, 换个方式说,身高\(height\)每增加一个1英寸, 收入\(earn\)会增加2387美元。

  • 对于截距\(\alpha = -126532\),即当身高为0时,期望的收入值-126532。呵呵,人的身高不可能为0,所以这是一种极端的理论情况,现实不可能发生。

12.5 多元线性回归

刚才讨论的单个预测变量height,现在我们增加一个预测变量ed,稍微扩展一下我们的一元线性模型,就是多元回归模型

\[ \begin{aligned} earn &= \alpha + \beta_1 \text{height} + \beta_2 \text{ed} +\epsilon \\ \end{aligned} \]

R语言代码实现也很简单,只需要把变量ed增加在公式的右边

同样,我们打印mod2看看

## 
## Call:
## lm(formula = earn ~ height + ed, data = wages)
## 
## Coefficients:
## (Intercept)       height           ed  
##     -161541         2087         4118

大家试着解释下mod2. 😄

12.7 可能遇到的情形

根据同学们的建议,模型中涉及统计知识,留给统计老师讲,我们这里是R语言课,应该讲代码。 因此,这里再介绍几种线性回归中遇到的几种特殊情况

12.7.2 只有截距

## 
## Call:
## lm(formula = earn ~ 1, data = wages)
## 
## Coefficients:
## (Intercept)  
##       32446
## # A tibble: 1 x 1
##   mean_wages
##        <dbl>
## 1     32446.

12.7.3 分类变量

race变量就是数据框wages的一个分类变量,代表四个不同的种族。用分类变量做回归,本质上是各组之间的进行比较。

## # A tibble: 4 x 1
##   race    
##   <chr>   
## 1 white   
## 2 other   
## 3 hispanic
## 4 black
## Warning: Removed 869 rows containing non-finite values
## (stat_boxplot).

以分类变量作为解释变量,做线性回归

## 
## Call:
## lm(formula = earn ~ race, data = wages)
## 
## Coefficients:
##  (Intercept)  racehispanic     raceother  
##        28372         -2887          3905  
##    racewhite  
##         4993

tidyverse框架下,喜欢数据框的统计结果,因此,可用broom的tidy()函数将模型输出转换为数据框的形式

## # A tibble: 4 x 5
##   term         estimate std.error statistic  p.value
##   <chr>           <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    28372.     2781.    10.2   1.29e-23
## 2 racehispanic   -2887.     4515.    -0.639 5.23e- 1
## 3 raceother       3905.     6428.     0.608 5.44e- 1
## 4 racewhite       4993.     2929.     1.70  8.85e- 2

我们看到输出结果,只有race_hispanic、 race_other和race_white三个系数和Intercept截距,race_black去哪里了呢?

事实上,race变量里有4组,回归时,选择black为基线,hispanic的系数,可以理解为由black切换到hispanic,引起earn收入的变化(效应)

  • 对 black 组的估计,earn = 28372.09 = 28372.09
  • 对 hispanic组的估计,earn = 28372.09 + -2886.79 = 25485.30
  • 对 other 组的估计,earn = 28372.09 + 3905.32 = 32277.41
  • 对 white 组的估计,earn = 28372.09 + 4993.33 = 33365.42

分类变量的线性回归本质上就是方差分析

11 章专题讨论方差分析

12.7.4 因子变量

hispanic组的估计最低,适合做基线,因此可以将race转换为因子变量,这样方便调整因子先后顺序

## # A tibble: 6 x 2
##     earn race 
##    <dbl> <fct>
## 1 79571. white
## 2 96397. white
## 3 48711. white
## 4 80478. other
## 5 82089. white
## 6 15313. white

wages_fct替换wages,然后建立线性模型

## # A tibble: 4 x 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)   25485.     3557.     7.16  1.26e-12
## 2 racewhite      7880.     3674.     2.14  3.22e- 2
## 3 raceblack      2887.     4515.     0.639 5.23e- 1
## 4 raceother      6792.     6800.     0.999 3.18e- 1

以hispanic组作为基线,各组系数也调整了,但加上截距后,实际值是没有变的。

大家可以用sex变量试试看

12.7.5 一个分类变量和一个连续变量

如果预测变量是一个分类变量和一个连续变量

## (Intercept)      height     sexmale 
##    -32479.9       879.4     16874.2
  • height = 879.424 当sex保持不变时,height变化引起的earn变化
  • sexmale = 16874.158 当height保持不变时,sex变化(female变为male)引起的earn变化
## Warning: Removed 50 rows containing missing values
## (geom_point).

12.7.6 偷懒的写法

. is shorthand for “everything else.”

R 语言很多时候都出现了.,不同的场景,含义是不一样的。我会在后面第 17 章专门讨论这个问题, 这是一个非常重要的问题

12.7.7 交互项

##    (Intercept)         height        sexmale 
##       -12167.0          564.5       -30510.4 
## height:sexmale 
##          701.4
  • 对于女性,height增长1个单位,引起earn的增长564.5102
  • 对于男性,height增长1个单位,引起earn的增长564.5102 + 701.4065 = 1265.92
## Warning: Removed 50 rows containing missing values
## (geom_point).

12.7.8 predict vs fit

下次再写

12.7.9 回归和相关的关系

  • 相关,比如求两个变量的相关系数cor(x, y)
  • 回归,也是探寻自变量和因变量的关系,一般用来预测

回归分析中,如果自变量只有一个\(x\),也就是模型lm(y~x),那么回归和相关就有关联了。

比如:计算身高和收入两者的Pearson相关系数的平方

## [1] 0.08503

然后看看,身高和收入的线性模型

## [1] 0.08503

相关系数的平方 和 线性模型的\(R^2\)是相等的

12.8 延伸阅读

一篇极富思考性和启发性的文章《常见统计检验的本质是线性模型》