Прогнозировать значения из подобранной логистической модели по группам - PullRequest
Новая
       66

Прогнозировать значения из подобранной логистической модели по группам

0 голосов
/ 04 июня 2018

Попытка согласовать несколько логистических моделей с данными для разных округов, и в конце они хотели бы объединить все это в один фрейм данных (все округа, все прогнозируемые группы населения за указанные годы).

Вот данные: 

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 округами.

1 Ответ

0 голосов
/ 04 июня 2018

Как подсказывает @Esther в комментариях, хорошим первым шагом будет извлечение вашей функции анонимного предсказания в именованную.Также имеет смысл заставить функцию принимать прогнозируемые годы в качестве аргумента, а не фиксировать их внутри функции: 

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

Отлично!Теперь один из способов (как предложил @ 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, который содержитрезультат, если функция выполнена успешно, и NULL, если произошла ошибка;и 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, заключается в использовании 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 для науки о данных.

Создана в 2018-06-04 представпакет (v0.2.0).

...