按组预测拟合逻辑模型的值

Question

尝试将多个逻辑模型拟合到不同县的数据，并希望最终将它们全部放在一个数据框中（所有县、所有预测人口、指定年份）。

这是数据：

county <- structure(list(name = structure(c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 
5L, 5L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 
8L, 9L, 9L, 9L, 9L, 9L), .Label = c("Alachua", "Columbia", "Gilchrist", 
"Lake", "Levy", "Marion", "Orange", "Seminole", "Volusia"), class = 
"factor"), 
year = c(1920L, 1940L, 1970L, 1990L, 2010L, 1920L, 1940L, 
1970L, 1990L, 2010L, 1920L, 1940L, 1970L, 1990L, 2010L, 1920L, 
1940L, 1970L, 1990L, 2010L, 1920L, 1940L, 1970L, 1990L, 2010L, 
1920L, 1940L, 1970L, 1990L, 2010L, 1920L, 1940L, 1970L, 1990L, 
2010L, 1920L, 1940L, 1970L, 1990L, 2010L, 1920L, 1940L, 1970L, 
1990L, 2010L), pop = c(24662.84498, 38518.67335, 105080.0739, 
182378.0527, 247964.4355, 14353.67655, 16988.63031, 25423.53768, 
42636.12851, 67396.52047, 6955.297482, 4331.7027, 3661.621676, 
9835.709676, 16780.95117, 12812.1731, 27202.15681, 65668.28125, 
153585.2153, 297441.8053, 10034.20186, 12707.52359, 12911.58508, 
26370.47373, 41650.51535, 23990.09377, 31340.67059, 69056.41468, 
194358.0547, 334117.7792, 19825.73528, 68559.76913, 337259.2307, 
670422.46, 1140314.083, 11027.52715, 23881.62063, 91628.11201, 
298115.877, 438079.7446, 24526.72497, 55775.68449, 175004.8787, 
382885.1367, 516049.0225)), .Names = c("name", "year", "pop"
), row.names = c(NA, -45L), class = "data.frame")

这就是我最终得到的结果：

library(dplyr) 
county %>% 
    group_by(name) %>%
    (function(x) {
            fm<- nls(pop ~ SSlogis(year, phi1, phi2, phi3), data = x)
            timevalues <- c(1992, 2002, 2007, 2012)
            predict <- predict(fm,list(year=timevalues))
            cbind(predict, predict)
    })

但这只会给我一个包含四个数据点的列表：

out:
  predict  predict
[1,] 226713.5 226713.5
[2,] 293596.4 293596.4
[3,] 326455.5 326455.5
[4,] 357640.8 357640.8

不知道他们在哪个县？如果我单独使用此代码（不使用 groupby），我可以让它工作。但是我必须为每个县单独做，然后自己绑定，一旦我与超过 9 个县合作，这将变得乏味。

Answer 1

正如@Esther 在 cmets 中建议的那样，好的第一步是提取 你的匿名预测函数变成了一个命名的。也会有道理 使函数接受预测年份作为参数，而不是 在函数中修复它们：

predict_pop <- function(data, year) {
  model <- nls(pop ~ SSlogis(year, phi1, phi2, phi3), data = data)

  nd <- data.frame(year)
  pred <- predict(model, nd)

  cbind(nd, pred)
}

让我们检查一下这是否适用于完整数据：

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

years <- c(1992, 2002, 2007, 2012)
predict_pop(county, years)
#>   year     pred
#> 1 1992 226713.5
#> 2 2002 293596.4
#> 3 2007 326455.5
#> 4 2012 357640.8

太棒了！现在有一种方法（如 cmets 中的@eipi10 所建议）来拟合模型 对于每个县将首先将split() 的数据放入数据列表中 每个县的帧，然后使用lapply() 获得每个子集中的预测。

split(county, county$name) %>%
  lapply(predict_pop, years)
#> Error in nls(y ~ 1/(1 + exp((xmid - x)/scal)), data = xy, start = list(xmid = aux[[1L]], : step factor 0.000488281 reduced below 'minFactor' of 0.000976562

但是，这会导致错误：似乎无法拟合模型 对于一些县自己。您可能需要使用模型本身来解决这个问题；但是如果我们想从这个模型中预测那些模型可以适合的县，我们可以将预测函数修改为 处理模型不适合的情况。

一种方法是使用purrr::safely() 制作“安全”版本 nls() 函数的一部分，它不会在出现错误时停止一切，但是 而是总是返回一个包含两个元素的列表：result，其中包含 如果函数执行成功，则返回结果，如果有 错误;和一个error，如果发生错误，则包含错误。

通过安全建模功能，我们可以检查模型是否可以 拟合，如果没有，则返回 NA 作为预测而不是错误。 以下是预测功能的修改版本：

predict_pop <- function(data, year) {
  safe_nls <- function(...) purrr::safely(nls)(...)$result
  model <- safe_nls(pop ~ SSlogis(year, phi1, phi2, phi3), data = data)

  nd <- data.frame(year)
  pred <- NA_real_

  if (!is.null(model))
    pred <- predict(model, nd)

  cbind(nd, pred)
}

现在我们可以使用之前的技术形式来获得预测。我添加了一个 bind_rows() 调用将列表结果合并到一个数据框中：

split(county, county$name) %>%
  lapply(predict_pop, years) %>% 
  bind_rows(.id = "county") %>% 
  head()
#>     county year     pred
#> 1  Alachua 1992 186020.6
#> 2  Alachua 2002 222332.3
#> 3  Alachua 2007 239432.0
#> 4  Alachua 2012 255440.9
#> 5 Columbia 1992       NA
#> 6 Columbia 2002       NA

在这里，我们可以看到哥伦比亚县之一的缺失预测 模型拟合失败。

还有其他几种方法可以预测每个县。一个这样的 @rawr 和 @Esther 在 cmets 中提到的替代方法是使用 do():

county %>% 
  group_by(name) %>% 
  do(predict_pop(., years)) %>% 
  head()
#> # A tibble: 6 x 3
#> # Groups:   name [2]
#>   name      year    pred
#>   <fct>    <dbl>   <dbl>
#> 1 Alachua   1992 186021.
#> 2 Alachua   2002 222332.
#> 3 Alachua   2007 239432.
#> 4 Alachua   2012 255441.
#> 5 Columbia  1992     NA 
#> 6 Columbia  2002     NA

另一种方法是通过分配分组来创建“嵌套”数据框 数据到带有tidyr::nest() 的列表列中。然后我们可以使用lapply() 从模型中获取每个数据子集的预测，最后 tidyr::unnest() 从列表列中获取预测。

county %>% 
  tidyr::nest(-name) %>% 
  tidyr::unnest(lapply(data, predict_pop, years)) %>% 
  head()
#>       name year     pred
#> 1  Alachua 1992 186020.6
#> 2  Alachua 2002 222332.3
#> 3  Alachua 2007 239432.0
#> 4  Alachua 2012 255440.9
#> 5 Columbia 1992       NA
#> 6 Columbia 2002       NA

我们有它：处理许多模型的一整套技术。有关这方面的进一步讨论和示例，您可能对 R for Data Science 一书中的many models chapter 感兴趣。

由reprex package (v0.2.0) 于 2018 年 6 月 4 日创建。