建模过程：选定模型族- 拟合模型
1.准备工作
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()
library(modelr)
options(na.action = na.warn)
2.简单模型
认识一下sim
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)
## 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 ...
散点图
ggplot(sim1, aes(x, y)) +
geom_point()
确定模型的基本形式
--线性--y= a_0 + a_1 * x
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实际值：响应变量
计算预测值与实际值之间的差异
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个距离，实际上应该用均方根误差来表示，量化为一个数值。
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计算所有模型的均方根误差
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个模型覆盖到数据上。
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个最佳模型
ggplot(models, aes(a1, a2)) +
geom_point(
data = filter(models, rank(dist) <= 10),
size = 4, color = "red" )+
geom_point(aes(colour = -dist))
更加系统化的搜索模型参数的方法：
网格搜索法 这里的“系统”是相对于之前生成随机数的操作说的。
生成网格，显示最佳10个模型
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))
第三种方法：牛顿—拉夫逊搜索
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()
sim1_mod <- lm(y ~ x, data = sim1)
coef(sim1_mod)
## (Intercept) x
## 4.220822 2.051533
emmm 所以lm这么厉害，我们为啥要学之前那三种奇怪的办法。有一种被大神耍了的感觉，是不是只是为了强调一下lm很厉害，😢。
3.模型可视化
（1）预测
data_grid() 生成一个数值范围，覆盖数据所在区域。
grid <- sim1 %>%
data_grid(x)
grid
## # A tibble: 10 x 1
## x
## <int>
## 1 1
## 2 2
## 3 3
## 4 4
## 5 5
## 6 6
## 7 7
## 8 8
## 9 9
## 10 10
add_predictions() 添加预测值
grid <- grid %>%
add_predictions(sim1_mod)
grid
## # A tibble: 10 x 2
## x pred
## <int> <dbl>
## 1 1 6.27
## 2 2 8.32
## 3 3 10.4
## 4 4 12.4
## 5 5 14.5
## 6 6 16.5
## 7 7 18.6
## 8 8 20.6
## 9 9 22.7
## 10 10 24.7
ggplot(sim1, aes(x)) +
geom_point(aes(y = y)) +
geom_line(
aes(y = pred),
data = grid,
colour = "red",
size = 1
)
（2）残差
--观测值与预测值之间的距离
--在原始数据种添加残差列
sim1 <- sim1 %>%
add_residuals(sim1_mod)
sim1
## # A tibble: 30 x 3
## x y resid
## <int> <dbl> <dbl>
## 1 1 4.20 -2.07
## 2 1 7.51 1.24
## 3 1 2.13 -4.15
## 4 2 8.99 0.665
## 5 2 10.2 1.92
## 6 2 11.3 2.97
## 7 3 7.36 -3.02
## 8 3 10.5 0.130
## 9 3 10.5 0.136
## 10 4 12.4 0.00763
## # ... with 20 more rows
4.公式和模型族
df <- tribble(
~y, ~x1, ~x2,
4, 2, 5,
5, 1, 6
)
model_matrix(df, y ~ x1)
## # A tibble: 2 x 2
## `(Intercept)` x1
## <dbl> <dbl>
## 1 1 2
## 2 1 1
model_matrix(df, y ~ x1 - 1)
## # A tibble: 2 x 1
## x1
## <dbl>
## 1 2
## 2 1
model_matrix(df, y ~ x1 + x2)
## # A tibble: 2 x 3
## `(Intercept)` x1 x2
## <dbl> <dbl> <dbl>
## 1 1 2 5
## 2 1 1 6
？？？暂时不懂这有啥用。可能是例子不太好。
（1）分类变量
df <- tribble(
~ sex, ~ response,
"male", 1,
"female", 2,
"male", 1
)
model_matrix(df, response ~ sex)
## # A tibble: 3 x 2
## `(Intercept)` sexmale
## <dbl> <dbl>
## 1 1 1
## 2 1 0
## 3 1 1
这里就可以类比到逻辑值和整数的转换，TRUE就是1，FALSE就是0，所以response的数值不管是几，转换后的值都是1和0，至于谁对应1，看结果第二列的列名，sex后面是谁。感觉这个例子很。。。
总结起来这一小节就说了一件事，分类变量模型预测值是观测值得均值。
（2）交互项
有两个预测变量（x1,x2）需要共同预测
mod1 <- lm(y ~ x1 + x2, data = sim3)
mod2 <- lm(y ~ x1 * x2, data = sim3)
grid <- sim3 %>%
data_grid(x1, x2) %>%
gather_predictions(mod1, mod2)
grid
## # A tibble: 80 x 4
## model x1 x2 pred
## <chr> <int> <fct> <dbl>
## 1 mod1 1 a 1.67
## 2 mod1 1 b 4.56
## 3 mod1 1 c 6.48
## 4 mod1 1 d 4.03
## 5 mod1 2 a 1.48
## 6 mod1 2 b 4.37
## 7 mod1 2 c 6.28
## 8 mod1 2 d 3.84
## 9 mod1 3 a 1.28
## 10 mod1 3 b 4.17
## # ... with 70 more rows
两个模型的差别是：x1和x2是相互独立还是相互影响。
比较两个模型哪个更好
sim3 <- sim3 %>%
gather_residuals(mod1, mod2)
ggplot(sim3, aes(x1, resid, color = x2)) +
geom_point() +
facet_grid(model ~ x2)
用残差分布来比较。没有给出如何评价两个模型的方法，大概就是好的模型残差分布应该是比较均匀的（个人理解）。
引自评论soleil：残差那个应该是这样，响应值＝可预测部分＋不可预测的和随机误差。所以在去掉我们fit的模型即可预测的部分后，剩余的残差就应该是不可预测和随机误差，一般情况下是随机分布，即0上下随机，如果能看出有一定模式，就像mod1那个绿色，可能就说明有漏掉可预测部分，缺少交互项或者高价项。
（3） 交互项（两个连续变量）
seq_range()函数生成向量和最小值之间间隔相等的n个值。 参数： pretty 生成整齐一些的序列 trim/expand使数据集中一些或分散一些
mod1 <- lm(y ~ x1 + x2, data = sim4)
mod2 <- lm(y ~ x1 * x2, data = sim4)
grid <- sim4 %>%
data_grid(
x1 = seq_range(x1, 5),
x2 = seq_range(x2, 5)
) %>%
gather_predictions(mod1, mod2)
grid
## # A tibble: 50 x 4
## model x1 x2 pred
## <chr> <dbl> <dbl> <dbl>
## 1 mod1 -1 -1 0.996
## 2 mod1 -1 -0.5 -0.395
## 3 mod1 -1 0 -1.79
## 4 mod1 -1 0.5 -3.18
## 5 mod1 -1 1 -4.57
## 6 mod1 -0.5 -1 1.91
## 7 mod1 -0.5 -0.5 0.516
## 8 mod1 -0.5 0 -0.875
## 9 mod1 -0.5 0.5 -2.27
## 10 mod1 -0.5 1 -3.66
## # ... with 40 more rows
可视化比较两个模型
ggplot(grid, aes(x1, pred, color = x2, group = x2)) +
geom_line() +
facet_wrap(~ model)
ggplot(grid, aes(x2, pred, color = x1, group = x1)) +
geom_line() +
facet_wrap(~ model)
好像并没有太好的办法？大神说还可以逐步完善。（那你讲这干啥）
４．变量转换
大概就是一元多次方程的意思，x2,x3,x^n这样。一元一次方程对应的是直线，多次对应的是各种曲线。
I()包装+，*和^进行变量转换。
用poly简化x^n输入
df <- tribble(
~y, ~x,
1, 1,
2, 2,
3, 3
)
model_matrix(df, y ~ x^2 + x) #不包装时，x^2被删除
## # A tibble: 3 x 2
## `(Intercept)` x
## <dbl> <dbl>
## 1 1 1
## 2 1 2
## 3 1 3
model_matrix(df, y ~ I(x^2) + x) #想要的结果
## # A tibble: 3 x 3
## `(Intercept)` `I(x^2)` x
## <dbl> <dbl> <dbl>
## 1 1 1 1
## 2 1 4 2
## 3 1 9 3
model_matrix(df, y ~ poly(x, 2)) #简化输入
## # A tibble: 3 x 3
## `(Intercept)` `poly(x, 2)1` `poly(x, 2)2`
## <dbl> <dbl> <dbl>
## 1 1 -7.07e- 1 0.408
## 2 1 -7.85e-17 -0.816
## 3 1 7.07e- 1 0.408
library(splines)
model_matrix(df, y ~ ns(x, 2)) #自然样条法 splines::ns()修正poly()超出数据范围的问题（还是不用poly直接用这个ns吧）
## # A tibble: 3 x 3
## `(Intercept)` `ns(x, 2)1` `ns(x, 2)2`
## <dbl> <dbl> <dbl>
## 1 1 0 0
## 2 1 0.566 -0.211
## 3 1 0.344 0.771
非线性函数的模型拟合
#准备数据
sim5 <- tibble(
x = seq(0, 3.5 * pi, length = 50),
y = 4 * sin(x) + rnorm(length(x))
)
ggplot(sim5, aes(x, y)) +
geom_point()
#用5个模型拟合
mod1 <- lm(y ~ ns(x, 1), data = sim5)
mod2 <- lm(y ~ ns(x, 2), data = sim5)
mod3 <- lm(y ~ ns(x, 3), data = sim5)
mod4 <- lm(y ~ ns(x, 4), data = sim5)
mod5 <- lm(y ~ ns(x, 5), data = sim5)
grid <- sim5 %>%
data_grid(x = seq_range(x, n = 50, expand = 0.1)) %>%
gather_predictions(mod1, mod2, mod3, mod4, mod5, .pred = "y")
#可视化、分面对比
ggplot(sim5, aes(x, y)) +
geom_point() +
geom_line(data = grid, color = "red") +
facet_wrap(~ model)
5.缺失值
默认情况是产生警告信息
df <- tribble(
~x, ~y,
1, 2.2,
2, NA,
3, 3.5,
4, 8.3,
NA, 10
)
mod <- lm(y ~ x, data = df)
## Warning: Dropping 2 rows with missing values
不想看警告信息
mod <- lm(y ~ x, data = df, na.action = na.exclude)
看实际使用了多少个有效值（非缺失值）
nobs(mod)
微信公众号生信星球同步更新我的文章
友情链接：
生信技能树公益视频合辑：学习顺序是linux，r，软件安装，geo，小技巧，ngs组学！
B站链接：https://m.bilibili.com/space/338686099
YouTube链接：https://m.youtube.com/channel/UC67sImqK7V8tSWHMG8azIVA/playlists
生信工程师入门最佳指南：https://mp.weixin.qq.com/s/vaX4ttaLIa19MefD86WfUA
学徒培养：https://mp.weixin.qq.com/s/3jw3_PgZXYd7FomxEMxFmw
资料大全：https://mp.weixin.qq.com/s/QcES9u1vYh-l6LMXPgJIlA
网友评论