小洁详解《R数据科学》--第17章 使用modelr实现基础模型

hijack

建模过程：选定模型族- 拟合模型

1.准备工作

[AppleScript] 纯文本查看 复制代码

library(tidyverse)
## -- Attaching packages -------------------

## √ ggplot2 3.1.0     √ purrr   0.2.5
## √ tibble  1.4.2     √ dplyr   0.7.7
## √ tidyr   0.8.2     √ stringr 1.3.1
## √ readr   1.1.1     √ forcats 0.3.0

## -- Conflicts --- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

[AppleScript] 纯文本查看 复制代码

library(modelr)
options(na.action = na.warn)

2.简单模型

认识一下sim

[AppleScript] 纯文本查看 复制代码

sim1
## # A tibble: 30 x 2
##        x     y
##    <int> <dbl>
##  1     1  4.20
##  2     1  7.51
##  3     1  2.13
##  4     2  8.99
##  5     2 10.2 
##  6     2 11.3 
##  7     3  7.36
##  8     3 10.5 
##  9     3 10.5 
## 10     4 12.4 
## # ... with 20 more rows
class(sim1)
## [1] "tbl_df"     "tbl"        "data.frame"
str(sim1)[mw_shl_code=applescript,true]ggplot(sim1, aes(x, y)) +
       geom_point()

## Classes 'tbl_df', 'tbl' and 'data.frame':    30 obs. of  2 variables:
##  $ x: int  1 1 1 2 2 2 3 3 3 4 ...
##  $ y: num  4.2 7.51 2.13 8.99 10.24 ...[/mw_shl_code]

散点图

确定模型的基本形式

—线性—y= a_0 + a_1 * x

[AppleScript] 纯文本查看 复制代码

models <- tibble(
       a1 = runif(250, -20, 40),
       a2 = runif(250, -5, 5)
)
ggplot(sim1, aes(x, y)) +
       geom_abline(
         aes(intercept = a1, slope = a2),
data = models, alpha = 1/4 )+
       geom_point()

计算出的y值：预测值 y实际值：响应变量

计算预测值与实际值之间的差异

[AppleScript] 纯文本查看 复制代码

model1 <- function(a, data) {
       a[1] + data$x * a[2]
}
model1(c(7, 1.5), sim1)
##  [1]  8.5  8.5  8.5 10.0 10.0 10.0 11.5 11.5 11.5 13.0 13.0 13.0 14.5 14.5
## [15] 14.5 16.0 16.0 16.0 17.5 17.5 17.5 19.0 19.0 19.0 20.5 20.5 20.5 22.0
## [29] 22.0 22.0

这里的a有两个元素，分别是模型的两个参数。

给出了30个距离，实际上应该用均方根误差来表示，量化为一个数值。

[AppleScript] 纯文本查看 复制代码

measure_distance <- function(mod, data) {
       diff <- data$y - model1(mod, data)
       sqrt(mean(diff ^ 2))
}
measure_distance(c(7, 1.5), sim1)
## [1] 2.665212

使用purrr计算所有模型的均方根误差

[AppleScript] 纯文本查看 复制代码

sim1_dist <- function(a1, a2) {
       measure_distance(c(a1, a2), sim1)
}
models <- models %>%
  mutate(dist = purrr::map2_dbl(a1, a2, sim1_dist))
models
## # A tibble: 250 x 3
##        a1     a2  dist
##     <dbl>  <dbl> <dbl>
##  1 -19.1   0.176 34.1 
##  2 -14.1  -3.37  50.6 
##  3   6.48  3.08   8.73
##  4   9.80  2.08   6.12
##  5   5.49 -3.74  34.9 
##  6  36.9  -4.82  20.5 
##  7  34.0   1.90  29.0 
##  8  22.3   1.65  16.1 
##  9   2.46 -4.19  40.3 
## 10  37.0   3.49  40.9 
## # ... with 240 more rows

将最好的（dist最小的）10个模型覆盖到数据上。

[AppleScript] 纯文本查看 复制代码

ggplot(sim1, aes(x, y)) +
       geom_point(size = 2, color = "grey30") +
       geom_abline(
         aes(intercept = a1, slope = a2, color = -dist),
         data = filter(models, rank(dist) <= 10)
       )

intercept是截距，slope是斜率。

a1和a2组成散点图，红色边框显示前10个最佳模型

[AppleScript] 纯文本查看 复制代码

ggplot(models, aes(a1, a2)) +
       geom_point(
         data = filter(models, rank(dist) <= 10),
size = 4, color = "red" )+
       geom_point(aes(colour = -dist))

更加系统化的搜索模型参数的方法：

网格搜索法 这里的“系统”是相对于之前生成随机数的操作说的。

生成网格，显示最佳10个模型

[AppleScript] 纯文本查看 复制代码

grid <- expand.grid(
       a1 = seq(-5, 20, length = 25),
       a2 = seq(1, 3, length = 25)
) %>%
       mutate(dist = purrr::map2_dbl(a1, a2, sim1_dist))

grid %>%
       ggplot(aes(a1, a2)) +
       geom_point(
         data = filter(grid, rank(dist) <= 10),
size = 4, colour = "red" )+
       geom_point(aes(color = -dist))

!](https://upload-images.jianshu.io/upload_images/9475888-bab16e1d583fb073.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)

第三种方法：牛顿—拉夫逊搜索

[AppleScript] 纯文本查看 复制代码

best <- optim(c(0, 0), measure_distance, data = sim1) 
best$par
## [1] 4.222248 2.051204
ggplot(sim1, aes(x, y)) +
       geom_point(size = 2, color = "grey30") +
       geom_abline(intercept = best$par[1], slope = best$par[2])

线性模型拟合函数—lm()

[AppleScript] 纯文本查看 复制代码

sim1_mod <- lm(y ~ x, data = sim1)
coef(sim1_mod)
## (Intercept)           x 
##    4.220822    2.051533

emmm 所以lm这么厉害，我们为啥要学之前那三种奇怪的办法。有一种被大神耍了的感觉，是不是只是为了强调一下lm很厉害，