logo标签列表

在R中针对每个子集拟合不同的模型

rdata.table
4

我有一个包含多个类别的大型数据集。我的目标是为每个类别拟合一个模型，然后预测结果并在一个区块中为每个类别可视化展示。

为了提供可重复的示例，我使用mtcars创建了一些基本内容。这对于每个类别的简单回归模型非常有效。

mtcars = data.table(mtcars)
model = mtcars[, list(fit = list(lm(mpg~disp+hp+wt))), keyby = cyl]
setkey(mtcars, cyl)
mtcars[model, pred := predict(i.fit[[1]], .SD), by = .EACHI]
ggplot(data = mtcars, aes(x = mpg, y = pred)) + geom_line() + facet_wrap(~cyl)

然而，我想尝试以下类似的内容，但目前还没有生效。这次尝试是使用公式列表，但我也想将不同的模型（一些glm，几棵树）发送到每个数据子集。

mtcars = data.table(mtcars)
factors = list(c("disp","wt"), c("disp"), c("hp"))
form = lapply(factors, function(x) as.formula(paste("mpg~",paste(x,collapse="+"))))
model = mtcars[, list(fit = list(lm(form))), keyby = cyl]
setkey(mtcars, cyl)
mtcars[model, pred := predict(i.fit[[1]], .SD), by = .EACHI]
ggplot(data = mtcars, aes(x = mpg, y = pred)) + geom_line() + facet_wrap(~cyl)
数据表真的必要吗？ - rawr
不，但在大型数据集上速度更快，所以更受欢迎。dplyr也可以。 - Divi
1
我只是在暗示瓶颈将会出现在预测，lm和ggplot上。list(fit = lapply(form, lm, data = .SD)) 这个是你想要的吗？ - rawr
是的，这就是我犯的错误。 - Divi
3个回答
4

以下是一种方法，我们为每个模型设置一个未评估的列表predict，在data.table对象中对它们进行评估，将输出gather并传递给ggplot：

models = quote(list(
      predict(lm(form[[1]], .SD)),
      predict(lm(form[[2]], .SD)), 
      predict(lm(form[[3]], .SD))))

d <- mtcars
d[, c("est1", "est2", "est3") := eval(models), by = cyl]
d <- tidyr::gather(d, key = model, value = pred, est1:est3)

library(ggplot2)
ggplot(d, aes(x = mpg, y = pred)) + geom_line() + facet_grid(cyl ~ model)

输出：

在这里输入图片描述


（注：此文本为中英混合语言，下同）
3

lm()也可以接受字符向量作为公式。因此，我会简单地创建form如下：

form = lapply(factors, function(x) paste("mpg~", paste(x, collapse="+")))

还有，您需要提供正确的数据（使用内置特殊符号 .SD 来对应每个组）：

model = mtcars[, list(fit=lapply(form, lm, data=.SD)), keyby=cyl]

对于每个，
都会被循环遍历，并且相应的公式将作为第一个参数传递给，每次与一起，其中<.SD>代表数据子集，并且本身就是一个data.table。您可以从vignettes中了解更多信息。
如果您还想在结果中获得公式，则：
chform = unlist(form)
model = mtcars[, list(form=chform, fit=lapply(form, lm, data=.SD)), keyby = cyl]

HTH

顺便提一下：如果您计划在使用data.tables中的[...]时使用update()，请阅读此帖子。

这解决了我目前面临的问题。唯一的问题是，我不明白为什么在拟合一个普通模型时，不提供data=.SD也能正常工作？ - Divi
公式对象还捕获它们创建时的环境...这就是将要使用的内容。请查看?lm。 - Arun
1
我目前正在做几乎完全相同的事情，所以时间非常完美。这将是一个以“tidyverse”为主的答案，但我真的很喜欢它的工作方式。 purrr有一些非常方便的map函数，当与tibble中的列表列结合使用时，可以使其非常平滑。使用您的定义（我不试图优化它）
library(data.table)
mtcars = data.table(mtcars)
factors = list(c("disp","wt"), c("disp"), c("hp"))
form = lapply(factors, function(x) as.formula(paste("mpg~",paste(x,collapse="+"))))

这提供了一个函数列表，这些函数可以传递给purrr :: invoke_map，它将一组参数（即您拥有的参数）应用于一组函数（在您的情况下，只有lm，但我怀疑也可以扩展到其他函数），可选参数（在您的示例中，mtcars）。使用tibble，它们被存储为整洁的类似data.frame的list，否则它们将作为lm对象返回。

library(tibble)
library(purrr) 
models <- tibble(fit = invoke_map(lm, form, data = mtcars))
models
#> # A tibble: 3 x 1
#>          fit
#>       <list>
#>   1 <S3: lm>
#>   2 <S3: lm>
#>   3 <S3: lm>

超级有用的部分是当你想对所有这些元素做一些事情时，比如提取拟合系数时：
map(models$fit, coefficients)
#> [[1]]
#> (Intercept)        disp          wt 
#> 34.96055404 -0.01772474 -3.35082533 
#> 
#> [[2]]
#> (Intercept)        disp 
#> 29.59985476 -0.04121512 
#> 
#> [[3]]
#> (Intercept)          hp 
#> 30.09886054 -0.06822828

重新检查所使用的公式。
map(models$fit, formula)
#> [[1]]
#> mpg ~ disp + wt
#> <environment: 0x0000000017ee73a8>
#>   
#>   [[2]]
#> mpg ~ disp
#> <environment: 0x0000000018392c58>
#>   
#>   [[3]]
#> mpg ~ hp
#> <environment: 0x0000000018471d18>

此外，如果您想从模型中添加一些预测结果，可以很容易地使用broom::augment实现。
library(broom)
models_with_predicts <- models %>% mutate(predict = map(fit, augment))
models_with_predicts
#> # A tibble: 3 x 2
#>          fit                predict
#>       <list>                 <list>
#>   1 <S3: lm> <data.frame [32 x 10]>
#>   2 <S3: lm>  <data.frame [32 x 9]>
#>   3 <S3: lm>  <data.frame [32 x 9]>

你可以通过使用 unnest() 返回数据级别（带有预测），但这将合并您的所有数据（添加分组级别以保持拟合分开）。
library(tidyr)
unnest(models_with_predicts, predict)

#> # A tibble: 96 x 11
#> mpg  disp    wt  .fitted   .se.fit     .resid       .hat   .sigma     .cooksd .std.resid    hp
#> <dbl> <dbl> <dbl>    <dbl>     <dbl>      <dbl>      <dbl>    <dbl>       <dbl>      <dbl> <dbl>
#>   1   21.0 160.0 2.620 23.34543 0.6075520 -2.3454326 0.04339369 2.933379 0.010222201 -0.8222164    NA
#> 2   21.0 160.0 2.875 22.49097 0.6221836 -1.4909721 0.04550894 2.954135 0.004351414 -0.5232550    NA
#> 3   22.8 108.0 2.320 25.27237 0.7326015 -2.4723669 0.06309504 2.928665 0.017217431 -0.8757799    NA
#> 4   21.4 258.0 3.215 19.61467 0.5743205  1.7853334 0.03877647 2.948162 0.005241995  0.6243627    NA
#> 5   18.7 360.0 3.440 17.05281 1.0943208  1.6471930 0.14078260 2.949120 0.020275438  0.6092882    NA
#> 6   18.1 225.0 3.460 19.37863 0.6122393 -1.2786309 0.04406584 2.957872 0.003089406 -0.4483953    NA
#> 7   14.3 360.0 3.570 16.61720 0.9897465 -2.3171997 0.11516157 2.931444 0.030948880 -0.8446199    NA
#> 8   24.4 146.7 3.190 21.67120 0.9053245  2.7287988 0.09635365 2.918183 0.034431234  0.9842424    NA
#> 9   22.8 140.8 3.150 21.90981 0.9165259  0.8901898 0.09875274 2.962885 0.003775416  0.3215070    NA
#> 10  19.2 167.6 3.440 20.46305 0.9678618 -1.2630477 0.11012510 2.957375 0.008693734 -0.4590766    NA
#> # ... with 86 more rows
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接