搜索
查看: 688|回复: 0

[R] 第17章

[复制链接]

31

主题

36

帖子

1130

积分

金牌会员

Rank: 6Rank: 6

积分
1130
发表于 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)



本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



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

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-10-22 23:04 , Processed in 0.037574 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.