搜索
查看: 1237|回复: 0

[R] 第17章

[复制链接]

31

主题

36

帖子

1550

积分

金牌会员

Rank: 6Rank: 6

积分
1550
发表于 2018-11-11 15:04:18 | 显示全部楼层 |阅读模式
# chapter 17
# 使用modelr实现基础模型
#1
# 建模分两个阶段
# (1)定义一个模型族表示一种精确但一般性的模式
# (2)生成一个拟合模型,方位为从模型族中找出最接近数据的一个模型
# 准备工作

[AppleScript] 纯文本查看 复制代码
library(tidyverse)

library(modelr)

options(na.action = na.warn)


#2 一个简单模型
# 以sim1数据集为例

[AppleScript] 纯文本查看 复制代码
ggplot(sim1, aes(x,y)) +

  geom_point()



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 =.25

  ) +

  geom_point()


# 以上代码的图片中有250个模型,大多很糟糕,找出最佳的模型

[AppleScript] 纯文本查看 复制代码
model1 <- function(a, data){

  a[1] + data$x * a[2]

}

model1(c(7, 1.5), sim1)


# 计算“均方根误差”

[AppleScript] 纯文本查看 复制代码
measure_distance <- function(mod, data){

  diff <- data$y - model1(mod, data)

  sqrt(mean(diff^2))

}

measure_distance(c(7, 1.5), sim1)


# 现在,使用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


# 接下来,我们将最好的10个模型覆盖到数据上,使用-dist为模型上色

[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)

  )


# 还可将这些模型看做观测,并使用由a1和a2组成的一张散点图表示它们

[AppleScript] 纯文本查看 复制代码
ggplot(models, aes(a1, a2)) + 

  geom_point(

    data = filter(models, rank(dist) <= 10),

    size = 4, color = "red"

  ) + 

  geom_point(aes(color = -dist))


# 相较于检查多个随机模型,可使用更系统化的方法找出模型参数(网格搜索法)

[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))


# 看一下10个最佳模型重新覆盖到原始数据的效果

[AppleScript] 纯文本查看 复制代码
ggplot(sim1, aes(x, y)) + 

  geom_point(size = 2, color = "grey30") + 

  geom_abline(

    aes(intercept = a1, slope = a2, color = -dist),

    data = filter(grid ,rank(dist) <= 10)

  )


# 更佳的一种方式,“牛顿-拉夫逊搜索”的数值化最小化工具,
# 先选择一个起点,寻找最陡的斜坡,向下滑行一小段,不断重复该过程,直到不能下滑为止
# 使用optim()函数完成该任务

[AppleScript] 纯文本查看 复制代码
best <- optim(c(0, 0), measure_distance, data = sim1)

best$par



ggplot(sim1, aes(x, y)) + 

  geom_point(size = 2, color = "grey30") + 

  geom_abline(intercept = best$par[1], slope = best$par[2])


# 另一种方法,线性模型法

[AppleScript] 纯文本查看 复制代码
sim1_mod <- lm(y ~ x, data = sim1)

coef(sim1_mod)



#3 模型可视化
#预测——告知模型捕获的模式
# 第一步——要对模型的预测进行可视化表示,首先需要生成一个分部均匀的数值网络,以覆盖数据所在区域

[AppleScript] 纯文本查看 复制代码
grid <- sim1 %>%

  data_grid(x)

grid


# 接着使用modelr::add_prediction()函数添加预测值,
# 该函数可将模型的预测值作为一个新列添加到数据框中,其参数是一个数据框和一个模型

[AppleScript] 纯文本查看 复制代码
grid <- grid %>%

  modelr::add_predictions(sim1_mod)

grid


# 第二步——绘制预测值

[AppleScript] 纯文本查看 复制代码
ggplot(sim1, aes(x)) + 

  geom_point(aes(y = y)) + 

  geom_line(

    aes(y = pred),

    data = grid,

    color = "red",

    size =1

  )


#残差——与预测值相对应,表示模型漏掉的部分
# 使用add_residuals()将残差添加到数据中,
# 但是使用的是原始数据,不是生成的网络,因为计算残差需要使用实际的y值

[AppleScript] 纯文本查看 复制代码
sim1 <- sim1 %>%

  add_residuals(sim1_mod)


# 对于残差可以反映出模型的信息的不同方法
# 1 简单的绘制频率多边形图

[AppleScript] 纯文本查看 复制代码
ggplot(sim1, aes(resid)) + 

  geom_freqpoly(binwidth = .5)


#4 公式和模型族

[AppleScript] 纯文本查看 复制代码
df <- tribble(

  ~y, ~x1, ~x2,

  4, 2, 5,

  5, 1, 6

)


# R中增加截距项的方法是加入一个值全是1的列,默认情况下,R总是添加这一列。
# 如不想要截距项,必须使用-1来明确丢弃它

[AppleScript] 纯文本查看 复制代码
model_matrix(df, y~ x1-1)


# 如想添加更多变量,那么模型矩阵也会随之增长

[AppleScript] 纯文本查看 复制代码
model_matrix(df, y ~ x1 +x2)


#4.1 分类变量

[AppleScript] 纯文本查看 复制代码
df <- tribble(

  ~sex, ~ response,

  "male", 1,

  "female", 2,

  "male", 1

)

model_matrix(df, response ~ sex)

ggplot(sim2) + 

  geom_point(aes(x, y))


# 拟合一个模型,生成预测

[AppleScript] 纯文本查看 复制代码
mod2 <- lm(y ~ x, data = sim2)

grid <- sim2 %>%

  data_grid(x) %>%

  add_predictions(mod2)

grid


#带有分类变量x的模型会为每个分类预测出均值

[AppleScript] 纯文本查看 复制代码
ggplot(sim2, aes(x)) + 

  geom_point(aes(y = y)) +

  geom_point(

    data = grid, 

    aes(y = pred),

    color = "red",

    size = 4

  )


# 不能对未观测到的水平进行预测

#4.2 交互项(连续变量与分类变量)
# 将连续变量和分类变量组合使用

[AppleScript] 纯文本查看 复制代码
ggplot(sim3, aes(x1, y)) + 

  geom_point(aes(color = x2))



mod1 <- lm(y ~ x1 + x2, data = sim3)

mod2 <- lm(y ~ x1 * x2, data = sim3)


# 因为有两个预测变量, 所以需要将这两个变量都传给data_grid ()函数,
# 该函数会找出x1和x2中的所有唯一值,并生成所有组合
# 要想为以上的两个模型同时生成预测,可以使用gather_predictions()函数,
# 可以将每个预测作为一行加入数据框。与gather_predictions()互补的函数是spread_predictions()
# spread_predictions()可将每个预测作为一列加入数据框

[AppleScript] 纯文本查看 复制代码
grid <- sim3 %>%

  data_grid(x1, x2) %>%

  gather_predictions(mod1, mod2)


# 可以使用分面将两个模型的可视化结果放在一张图中

[AppleScript] 纯文本查看 复制代码
ggplot(sim3, aes(x1, y, color =x2)) +

  geom_point() + 

  geom_line(data = grid, aes(y = pred)) +

  facet_wrap(~model)



sim3 <- sim3 %>%

  gather_residuals(mod1, mod2)

ggplot(sim3, aes(x1, resid, color = x2)) +

  geom_point() + 

  facet_grid(model ~ x2)


#4.3 交互项(两个连续变量)

[AppleScript] 纯文本查看 复制代码
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


# pretty = T 会生成一个看起来比较舒服的序列,对于生成输出表格,该参数非常有用

[AppleScript] 纯文本查看 复制代码
seq_range(c(.0123, .923423), n = 5 )

seq_range(c(.0123, .923423), n = 5 ,pretty = T)


# trim = .1 会阶段10%的尾部值,如变量具有长尾分布, 而希望生成中心附近的值,就可使用该参数

[AppleScript] 纯文本查看 复制代码
x1 <- rcauchy(100)

seq_range(x1, n = 5)

seq_range(x1, n = 5, trim = .1)

seq_range(x1, n = 5, trim = .25)


# expand = .1 某种程度是trim()的反函数,可将取值范围扩大10%

[AppleScript] 纯文本查看 复制代码
x2 <- c(0, 1)

seq_range(x2, n = 5)


# 使用geom_tile()显示模型

[AppleScript] 纯文本查看 复制代码
ggplot(grid, aes(x1, x2)) +

  geom_tile(aes(fill = pred)) +

  facet_wrap(~ model)


# 以上代码不易看出差异,使用geom_line()显示

[AppleScript] 纯文本查看 复制代码
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)


# 交叉项说明两变量是相互影响的,如果要预测y值,必须同时考虑x1的值和x2的值

#4.4 变量转换

[AppleScript] 纯文本查看 复制代码
df <- tribble(

  ~y,~x,

  1,1,

  2,2,

  3,3

)

model_matrix(df, y ~ x^2 + x)

model_matrix(df, y ~ I(x^2) + x)


# 变量转换非常有用,可使用它们来近似表示非线性函数,poly()辅助输入拟合方程

[AppleScript] 纯文本查看 复制代码
model_matrix(df, y ~ poly(x, 2))


# 使用poly()时,多项式的值会超出数据范围,很容易接近正无穷或负无穷

[AppleScript] 纯文本查看 复制代码
library(splines)

model_matrix(df, y ~ ns(x, 2))


# 近似非线性函数的情况

[AppleScript] 纯文本查看 复制代码
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 缺失值

[AppleScript] 纯文本查看 复制代码
df <- tribble(

  ~x,~y,

  1,2.2,

  2,NA,

  3,3.5,

  4,8.3,

  NA,10

)

mod <- lm(y ~ x,data = df)


# 要想防止错误信息,可以设置na.action = na.exclude

[AppleScript] 纯文本查看 复制代码
mod <- lm(y ~ x, data = df, na.action = na.exclude)


# 通过nobs()函数,可知道模型实际使用了多少个观测
[AppleScript] 纯文本查看 复制代码
nobs(mod)

sin.png sim1 models.png resid frepoly.png prediction 1.png net search plot.png netsearch fitting 1.png interaction 1.png interaction 2.png interaction 3.png interaction 4.png facet.png facet 2.png dot plot.png category 1.png category 2.png best$par.png abline plot.png




上一篇:想用Deseq2来做表达差异分析,但是老师执着于FPKM
下一篇:慧美——R for data science 第九章 使用dplyr处理关系数据
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2020-6-4 07:25 , Processed in 0.026834 second(s), 28 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.