按组为模型添加预测
Add predictions for models by group
我在我的数据集中按组估计回归模型，然后我希望为所有组添加正确的拟合值。
我正在尝试以下操作：
library(dplyr)
library(modelr)
df <- tribble(
~year, ~country, ~value,
2001, "France", 55,
2002, "France", 53,
2003, "France", 31,
2004, "France", 10,
2005, "France", 30,
2006, "France", 37,
2007, "France", 54,
2008, "France", 58,
2009, "France", 50,
2010, "France", 40,
2011, "France", 49,
2001, "USA", 55,
2002, "USA", 53,
2003, "USA", 64,
2004, "USA", 40,
2005, "USA", 30,
2006, "USA", 39,
2007, "USA", 55,
2008, "USA", 53,
2009, "USA", 71,
2010, "USA", 44,
2011, "USA", 40
)
rmod <- df %>%
group_by(country) %>%
do(fitModels = lm("value ~ year", data = .))
df <- df %>%
add_predictions(rmod)
抛出错误：
Error in UseMethod("predict") :
no applicable method for 'predict' applied to an object of class "c('rowwise_df', 'tbl_df', 'tbl', 'data.frame')"
我想得到一列包含国家/地区的每个拟合值或一列包含每个国家/地区的预测值。当模型在 do() 调用后保存为列表时，add_predictions() 函数似乎无法正常工作。
我会用 data.table 执行以下操作：
library(data.table)
setDT(df) # convert to data.table
df[ , value_hat := lm(value ~ year)$fitted.values, by = country]
如果您有 NA，一种选择是：
df[complete.cases(df), value_hat := lm(value ~ year)$fitted.values, by = country]
另一个实际使用 predict:
df[ , value_hat := predict(lm(value ~ year), .SD), by = country]
还有一些其他方法可以解决这个问题。
可能是最直接的，但是你失去了中间模型：
rmod <- df %>%
group_by(country) %>%
mutate(fit = lm(value ~ year)$fitted.values) %>%
ungroup
rmod
# # A tibble: 22 × 4
# year country value fit
# <dbl> <chr> <dbl> <dbl>
# 1 2001 France 55 38.13636
# 2 2002 France 53 39.00000
# 3 2003 France 31 39.86364
# 4 2004 France 10 40.72727
# 5 2005 France 30 41.59091
# 6 2006 France 37 42.45455
# 7 2007 France 54 43.31818
# 8 2008 France 58 44.18182
# 9 2009 France 50 45.04545
# 10 2010 France 40 45.90909
# # ... with 12 more rows
另一种方法是使用 "tidy" 模型将数据、模型和结果封装到框架内的各个单元格中：
rmod <- df %>%
group_by(country) %>%
nest() %>%
mutate(mdl = map(data, ~ lm(value ~ year, data=.))) %>%
mutate(fit = map(mdl, ~ .$fitted.values))
rmod
# # A tibble: 2 × 4
# country data mdl fit
# <chr> <list> <list> <list>
# 1 France <tibble [11 × 2]> <S3: lm> <dbl [11]>
# 2 USA <tibble [11 × 2]> <S3: lm> <dbl [11]>
此方法的优点是您可以根据需要访问模型的其他属性，也许 summary( filter(rmod, country == "France")$mdl[[1]] )。 （[[1]] 是必需的，因为使用 tibbles，$mdl 将始终 return 一个 list。）
您可以 extract/unnest 如下所示：
select(rmod, -mdl) %>% unnest()
# # A tibble: 22 × 4
# country fit year value
# <chr> <dbl> <dbl> <dbl>
# 1 France 38.13636 2001 55
# 2 France 39.00000 2002 53
# 3 France 39.86364 2003 31
# 4 France 40.72727 2004 10
# 5 France 41.59091 2005 30
# 6 France 42.45455 2006 37
# 7 France 43.31818 2007 54
# 8 France 44.18182 2008 58
# 9 France 45.04545 2009 50
# 10 France 45.90909 2010 40
# # ... with 12 more rows
（不幸的是，这些列被重新排序了，但这很美观并且很容易补救。）
编辑
如果您 want/need 在此处使用 modelr- 细节，请尝试：
rmod <- df %>%
group_by(country) %>%
nest() %>%
mutate(mdl = map(data, ~ lm(value ~ year, data=.))) %>%
mutate(fit = map(mdl, ~ .$fitted.values)) %>%
mutate(data = map2(data, mdl, add_predictions))
rmod
# # A tibble: 2 x 4
# country data mdl fit
# <chr> <list> <list> <list>
# 1 France <tibble [11 x 3]> <S3: lm> <dbl [11]>
# 2 USA <tibble [11 x 3]> <S3: lm> <dbl [11]>
select(rmod, -mdl, -fit) %>% unnest()
# # A tibble: 22 x 4
# country year value pred
# <chr> <dbl> <dbl> <dbl>
# 1 France 2001. 55. 38.1
# 2 France 2002. 53. 39.0
# 3 France 2003. 31. 39.9
# 4 France 2004. 10. 40.7
# 5 France 2005. 30. 41.6
# 6 France 2006. 37. 42.5
# 7 France 2007. 54. 43.3
# 8 France 2008. 58. 44.2
# 9 France 2009. 50. 45.0
# 10 France 2010. 40. 45.9
# # ... with 12 more rows
这是使用 broom 包而不是 modelr 的替代方法。 augment 添加拟合值以及其他有用信息，例如原始观测值的残差。它旨在与符合 do 的分组模型的输出完美配合。见下文：
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
library(broom)
df <- tribble(
~year, ~country, ~value,
2001, "France", 55,
2002, "France", 53,
2003, "France", 31,
2004, "France", 10,
2005, "France", 30,
2006, "France", 37,
2007, "France", 54,
2008, "France", 58,
2009, "France", 50,
2010, "France", 40,
2011, "France", 49,
2001, "USA", 55,
2002, "USA", 53,
2003, "USA", 64,
2004, "USA", 40,
2005, "USA", 30,
2006, "USA", 39,
2007, "USA", 55,
2008, "USA", 53,
2009, "USA", 71,
2010, "USA", 44,
2011, "USA", 40
)
rmod <- df %>%
group_by(country) %>%
do(fitModels = lm("value ~ year", data = .))
rmod %>%
augment(fitModels)
#> # A tibble: 22 x 10
#> # Groups: country [2]
#> country value year .fitted .se.fit .resid .hat .sigma .cooksd
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 France 55. 2001. 38.1 8.49 16.9 0.318 14.2 0.430
#> 2 France 53. 2002. 39.0 7.31 14.0 0.236 14.9 0.176
#> 3 France 31. 2003. 39.9 6.25 -8.86 0.173 15.6 0.0438
#> 4 France 10. 2004. 40.7 5.37 -30.7 0.127 10.9 0.349
#> 5 France 30. 2005. 41.6 4.76 -11.6 0.1000 15.4 0.0366
#> 6 France 37. 2006. 42.5 4.54 -5.45 0.0909 15.8 0.00723
#> 7 France 54. 2007. 43.3 4.76 10.7 0.100 15.5 0.0311
#> 8 France 58. 2008. 44.2 5.37 13.8 0.127 15.1 0.0705
#> 9 France 50. 2009. 45.0 6.25 4.95 0.173 15.8 0.0137
#> 10 France 40. 2010. 45.9 7.31 -5.91 0.236 15.8 0.0313
#> # ... with 12 more rows, and 1 more variable: .std.resid <dbl>
由 reprex package (v0.2.0) 创建于 2018-04-19。
我在我的数据集中按组估计回归模型，然后我希望为所有组添加正确的拟合值。
我正在尝试以下操作：
library(dplyr)
library(modelr)
df <- tribble(
~year, ~country, ~value,
2001, "France", 55,
2002, "France", 53,
2003, "France", 31,
2004, "France", 10,
2005, "France", 30,
2006, "France", 37,
2007, "France", 54,
2008, "France", 58,
2009, "France", 50,
2010, "France", 40,
2011, "France", 49,
2001, "USA", 55,
2002, "USA", 53,
2003, "USA", 64,
2004, "USA", 40,
2005, "USA", 30,
2006, "USA", 39,
2007, "USA", 55,
2008, "USA", 53,
2009, "USA", 71,
2010, "USA", 44,
2011, "USA", 40
)
rmod <- df %>%
group_by(country) %>%
do(fitModels = lm("value ~ year", data = .))
df <- df %>%
add_predictions(rmod)
抛出错误：
Error in UseMethod("predict") :
no applicable method for 'predict' applied to an object of class "c('rowwise_df', 'tbl_df', 'tbl', 'data.frame')"
我想得到一列包含国家/地区的每个拟合值或一列包含每个国家/地区的预测值。当模型在 do() 调用后保存为列表时，add_predictions() 函数似乎无法正常工作。
我会用 data.table 执行以下操作：
library(data.table)
setDT(df) # convert to data.table
df[ , value_hat := lm(value ~ year)$fitted.values, by = country]
如果您有 NA，一种选择是：
df[complete.cases(df), value_hat := lm(value ~ year)$fitted.values, by = country]
另一个实际使用 predict:
df[ , value_hat := predict(lm(value ~ year), .SD), by = country]
还有一些其他方法可以解决这个问题。
可能是最直接的，但是你失去了中间模型：
rmod <- df %>%
group_by(country) %>%
mutate(fit = lm(value ~ year)$fitted.values) %>%
ungroup
rmod
# # A tibble: 22 × 4
# year country value fit
# <dbl> <chr> <dbl> <dbl>
# 1 2001 France 55 38.13636
# 2 2002 France 53 39.00000
# 3 2003 France 31 39.86364
# 4 2004 France 10 40.72727
# 5 2005 France 30 41.59091
# 6 2006 France 37 42.45455
# 7 2007 France 54 43.31818
# 8 2008 France 58 44.18182
# 9 2009 France 50 45.04545
# 10 2010 France 40 45.90909
# # ... with 12 more rows
另一种方法是使用 "tidy" 模型将数据、模型和结果封装到框架内的各个单元格中：
rmod <- df %>%
group_by(country) %>%
nest() %>%
mutate(mdl = map(data, ~ lm(value ~ year, data=.))) %>%
mutate(fit = map(mdl, ~ .$fitted.values))
rmod
# # A tibble: 2 × 4
# country data mdl fit
# <chr> <list> <list> <list>
# 1 France <tibble [11 × 2]> <S3: lm> <dbl [11]>
# 2 USA <tibble [11 × 2]> <S3: lm> <dbl [11]>
此方法的优点是您可以根据需要访问模型的其他属性，也许 summary( filter(rmod, country == "France")$mdl[[1]] )。 （[[1]] 是必需的，因为使用 tibbles，$mdl 将始终 return 一个 list。）
您可以 extract/unnest 如下所示：
select(rmod, -mdl) %>% unnest()
# # A tibble: 22 × 4
# country fit year value
# <chr> <dbl> <dbl> <dbl>
# 1 France 38.13636 2001 55
# 2 France 39.00000 2002 53
# 3 France 39.86364 2003 31
# 4 France 40.72727 2004 10
# 5 France 41.59091 2005 30
# 6 France 42.45455 2006 37
# 7 France 43.31818 2007 54
# 8 France 44.18182 2008 58
# 9 France 45.04545 2009 50
# 10 France 45.90909 2010 40
# # ... with 12 more rows
（不幸的是，这些列被重新排序了，但这很美观并且很容易补救。）
编辑
如果您 want/need 在此处使用 modelr- 细节，请尝试：
rmod <- df %>%
group_by(country) %>%
nest() %>%
mutate(mdl = map(data, ~ lm(value ~ year, data=.))) %>%
mutate(fit = map(mdl, ~ .$fitted.values)) %>%
mutate(data = map2(data, mdl, add_predictions))
rmod
# # A tibble: 2 x 4
# country data mdl fit
# <chr> <list> <list> <list>
# 1 France <tibble [11 x 3]> <S3: lm> <dbl [11]>
# 2 USA <tibble [11 x 3]> <S3: lm> <dbl [11]>
select(rmod, -mdl, -fit) %>% unnest()
# # A tibble: 22 x 4
# country year value pred
# <chr> <dbl> <dbl> <dbl>
# 1 France 2001. 55. 38.1
# 2 France 2002. 53. 39.0
# 3 France 2003. 31. 39.9
# 4 France 2004. 10. 40.7
# 5 France 2005. 30. 41.6
# 6 France 2006. 37. 42.5
# 7 France 2007. 54. 43.3
# 8 France 2008. 58. 44.2
# 9 France 2009. 50. 45.0
# 10 France 2010. 40. 45.9
# # ... with 12 more rows
这是使用 broom 包而不是 modelr 的替代方法。 augment 添加拟合值以及其他有用信息，例如原始观测值的残差。它旨在与符合 do 的分组模型的输出完美配合。见下文：
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
library(broom)
df <- tribble(
~year, ~country, ~value,
2001, "France", 55,
2002, "France", 53,
2003, "France", 31,
2004, "France", 10,
2005, "France", 30,
2006, "France", 37,
2007, "France", 54,
2008, "France", 58,
2009, "France", 50,
2010, "France", 40,
2011, "France", 49,
2001, "USA", 55,
2002, "USA", 53,
2003, "USA", 64,
2004, "USA", 40,
2005, "USA", 30,
2006, "USA", 39,
2007, "USA", 55,
2008, "USA", 53,
2009, "USA", 71,
2010, "USA", 44,
2011, "USA", 40
)
rmod <- df %>%
group_by(country) %>%
do(fitModels = lm("value ~ year", data = .))
rmod %>%
augment(fitModels)
#> # A tibble: 22 x 10
#> # Groups: country [2]
#> country value year .fitted .se.fit .resid .hat .sigma .cooksd
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 France 55. 2001. 38.1 8.49 16.9 0.318 14.2 0.430
#> 2 France 53. 2002. 39.0 7.31 14.0 0.236 14.9 0.176
#> 3 France 31. 2003. 39.9 6.25 -8.86 0.173 15.6 0.0438
#> 4 France 10. 2004. 40.7 5.37 -30.7 0.127 10.9 0.349
#> 5 France 30. 2005. 41.6 4.76 -11.6 0.1000 15.4 0.0366
#> 6 France 37. 2006. 42.5 4.54 -5.45 0.0909 15.8 0.00723
#> 7 France 54. 2007. 43.3 4.76 10.7 0.100 15.5 0.0311
#> 8 France 58. 2008. 44.2 5.37 13.8 0.127 15.1 0.0705
#> 9 France 50. 2009. 45.0 6.25 4.95 0.173 15.8 0.0137
#> 10 France 40. 2010. 45.9 7.31 -5.91 0.236 15.8 0.0313
#> # ... with 12 more rows, and 1 more variable: .std.resid <dbl>
由 reprex package (v0.2.0) 创建于 2018-04-19。