Прогнозное моделирование данных подсчета рождений (с линейной регрессией)
Мне дали эти данные (ниже), и я хочу попробовать прогностическое моделирование с этими данными. В частности, я хочу построить линейную регрессию для прогнозирования количества рождений (по месяцам, если возможно) до 2032 года. Я предполагаю, что мне придется использовать функцию predict().
Вот мои данные
В широком формате:
> dput(births_monthly_cross_022822_clean)
structure(list(month = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12
), `2010` = c(323816, 302551, 339219, 325582, 328960, 335180,
345875, 350473, 351439, 337477, 326868, 339665), `2011` = c(321104,
298537, 330752, 313875, 327242, 337890, 346192, 360079, 346254,
329174, 322135, 327986), `2012` = c(316959, 305060, 324944, 307262,
330779, 327811, 348157, 361902, 340901, 346251, 325810, 324960
), `2013` = c(324314, 292360, 321166, 312274, 330193, 320475,
349920, 354293, 338863, 341353, 319233, 336320), `2014` = c(327154,
299087, 324317, 319788, 335444, 326757, 356446, 355201, 349403,
344265, 318819, 341494), `2015` = c(326747, 298815, 329714, 321618,
328709, 331400, 354384, 352782, 348479, 339904, 319605, 336576
), `2016` = c(317445, 306750, 329341, 314312, 328434, 333166,
343643, 358737, 346525, 331939, 320258, 325562), `2017` = c(314597,
289694, 320327, 300801, 323169, 324633, 335708, 353033, 338012,
330762, 317280, 316738), `2018` = c(315593, 284940, 316824, 299125,
321448, 315585, 329851, 345687, 323685, 327740, 309475, 311581
), `2019` = c(311678, 280679, 304999, 299755, 317160, 304843,
334518, 342617, 326662, 326165, 298899, 309607), `2020` = c(305536,
283385, 302331, 290940, 301902, 302574, 322096, 320072, 312152,
305548, 283010, 290280)), class = c("spec_tbl_df", "tbl_df",
"tbl", "data.frame"), row.names = c(NA, -12L), spec = structure(list(
cols = list(month = structure(list(), class = c("collector_double",
"collector")), `2010` = structure(list(), class = c("collector_number",
"collector")), `2011` = structure(list(), class = c("collector_number",
"collector")), `2012` = structure(list(), class = c("collector_number",
"collector")), `2013` = structure(list(), class = c("collector_number",
"collector")), `2014` = structure(list(), class = c("collector_number",
"collector")), `2015` = structure(list(), class = c("collector_number",
"collector")), `2016` = structure(list(), class = c("collector_number",
"collector")), `2017` = structure(list(), class = c("collector_number",
"collector")), `2018` = structure(list(), class = c("collector_number",
"collector")), `2019` = structure(list(), class = c("collector_number",
"collector")), `2020` = structure(list(), class = c("collector_number",
"collector"))), default = structure(list(), class = c("collector_guess",
"collector")), skip = 1), class = "col_spec"))
Те же данные, но в длинном формате:
> dput(Birthcount)
structure(list(month = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7,
7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9,
9, 9, 9, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 11, 11,
11, 11, 11, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 12, 12,
12, 12, 12, 12), name = c("2010", "2011", "2012", "2013", "2014",
"2015", "2016", "2017", "2018", "2019", "2020", "2010", "2011",
"2012", "2013", "2014", "2015", "2016", "2017", "2018", "2019",
"2020", "2010", "2011", "2012", "2013", "2014", "2015", "2016",
"2017", "2018", "2019", "2020", "2010", "2011", "2012", "2013",
"2014", "2015", "2016", "2017", "2018", "2019", "2020", "2010",
"2011", "2012", "2013", "2014", "2015", "2016", "2017", "2018",
"2019", "2020", "2010", "2011", "2012", "2013", "2014", "2015",
"2016", "2017", "2018", "2019", "2020", "2010", "2011", "2012",
"2013", "2014", "2015", "2016", "2017", "2018", "2019", "2020",
"2010", "2011", "2012", "2013", "2014", "2015", "2016", "2017",
"2018", "2019", "2020", "2010", "2011", "2012", "2013", "2014",
"2015", "2016", "2017", "2018", "2019", "2020", "2010", "2011",
"2012", "2013", "2014", "2015", "2016", "2017", "2018", "2019",
"2020", "2010", "2011", "2012", "2013", "2014", "2015", "2016",
"2017", "2018", "2019", "2020", "2010", "2011", "2012", "2013",
"2014", "2015", "2016", "2017", "2018", "2019", "2020"), value = c(323816,
321104, 316959, 324314, 327154, 326747, 317445, 314597, 315593,
311678, 305536, 302551, 298537, 305060, 292360, 299087, 298815,
306750, 289694, 284940, 280679, 283385, 339219, 330752, 324944,
321166, 324317, 329714, 329341, 320327, 316824, 304999, 302331,
325582, 313875, 307262, 312274, 319788, 321618, 314312, 300801,
299125, 299755, 290940, 328960, 327242, 330779, 330193, 335444,
328709, 328434, 323169, 321448, 317160, 301902, 335180, 337890,
327811, 320475, 326757, 331400, 333166, 324633, 315585, 304843,
302574, 345875, 346192, 348157, 349920, 356446, 354384, 343643,
335708, 329851, 334518, 322096, 350473, 360079, 361902, 354293,
355201, 352782, 358737, 353033, 345687, 342617, 320072, 351439,
346254, 340901, 338863, 349403, 348479, 346525, 338012, 323685,
326662, 312152, 337477, 329174, 346251, 341353, 344265, 339904,
331939, 330762, 327740, 326165, 305548, 326868, 322135, 325810,
319233, 318819, 319605, 320258, 317280, 309475, 298899, 283010,
339665, 327986, 324960, 336320, 341494, 336576, 325562, 316738,
311581, 309607, 290280)), class = c("tbl_df", "tbl", "data.frame"
), row.names = c(NA, -132L))
Пока что моя грубая модель выглядит следующим образом:
mod.fert.natl <- lm( value ~ name, data = Birthcount)
Мне интересно, как лучше всего это сделать. Огромное спасибо заранее!
И используйте теги time-series и forecasting на Cross Validated.
Также имейте в виду, что в Cross Validated стандарт задавать хорошие вопросы отличается от стандарта в ОС. stats.stackexchange.com/help/how-to-ask
Ответы 1
Взяв ваши данные в длинном формате, нам нужно кое-что изменить. Годы должны быть числовыми, если мы хотим смоделировать долгосрочные тенденции. Месяцы будет легче отслеживать, если мы сделаем их факторами в правильном порядке. Мы также должны создать столбец даты, который позволит нам комбинировать год и месяц для построения графика, и мы должны дать всем этим столбцам репрезентативные имена:
Birthcount <- dplyr::transmute(Birthcount,
date = as.Date(paste(name, month, 1, sep = '-')),
month = factor(month.name[month], month.name),
year = as.numeric(name),
count = value)
Теперь создать модель несложно. Когда вы имеете дело с данными подсчета, вам в идеале следует использовать glm с family = poisson, но когда вы имеете дело с данными подсчета в сотнях тысяч, а минимальное значение также исчисляется сотнями тысяч, это действительно не имеет значения. практическая разница, используете ли вы glm или lm (см. Приложение). Чтобы быть строгим здесь, давайте использовать glm.
Предполагая, что вы хотите смоделировать сезонные колебания числа рождений по месяцам и тенденцию по годам, вы можете сделать следующее:
mod <- glm(count ~ month + year, data = Birthcount, family = poisson)
Чтобы получить прогнозы до 2032 года, нам нужна новая структура данных о годах и месяцах, которые мы хотим предсказать.
new_df <- expand.grid(month = factor(month.name, month.name), year = 2010:2032)
new_df$date <- seq(as.Date("2010-01-01"), as.Date("2032-12-01"), by = "month")
Теперь мы можем добавить наши прогнозируемые значения:
new_df$count <- predict(mod, newdata = new_df)
Теперь new_df содержит все наши прогнозы. Чтобы увидеть, как это выглядит, мы можем построить прогнозы и фактические данные:
library(ggplot2)
ggplot(Birthcount, aes(date, count)) +
geom_line(aes(colour = "Predicted"), data = new_df, linetype = 2) +
geom_point(aes(colour = "Actual")) +
scale_colour_manual(values = c("red4", "gray80")) +
ggtitle("Monthly births (Actual and Predicted) 2010 - 2032") +
labs(colour = "") +
theme_bw()
Приложение
Если мы рассмотрим прогнозы между glm и lm, мы увидим, что практическая разница для целей прогнозирования невелика.
mod_glm <- mod
mod_lm <- lm(count ~ month + year, data = Birthcount)
plot(predict(mod_glm, type = "response"), predict(mod_lm))
abline(a = 0, b = 1)
Вы уверены, что использовать линейную модель для подсчета данных — это хорошая идея?
@dipetkov в общем нет, вы бы использовали glm с семейством Пуассона и функцию ссылки на журнал, но в случаях, когда вы имеете дело с подсчетами в сотни тысяч, это различие не имеет практического значения. Я попробую это с glm, чтобы увидеть, как только вернусь к своему компьютеру, но я предполагаю, что разница будет незначительной, если вообще будет.
ОП и многие люди, которые натыкаются на этот вопрос и ответ, могут этого не знать. Я думаю, что ваш ответ выиграет от наличия баннера типа «Есть более подходящие модели для данных подсчета». И лучшие способы моделирования временных рядов в этом отношении.
@dipetkov У меня складывается впечатление, что ОП - новичок, и я не думаю, что сейчас самое подходящее время, чтобы провести различие там, где это не имеет значения - это сложно сделать, поскольку вы не хотите делать что-то более сложное, чем должно быть, но не хочу быть неточным. Я подозреваю, что у вас не возникнет проблем с использованием линейных моделей при работе, скажем, с долларовыми ценностями от 3000 до 4000 долларов, хотя технически это целые числа центов в диапазоне от 300 000 до 400 000.
@dipetkov, но вы правы, и я добавлю сноску, когда вернусь к своему компьютеру.
Я согласен! Какое-то время я думал о том, чтобы показать, как подогнать модель временных рядов к этим данным, но потом подумал: нет, это слишком сложно. Но я все же думаю, что лучше указать на это, чтобы сделать это правильно, анализ должен быть более сложным. А затем укажите OP на Cross Validated. Затем им решать, сколько сорняков они хотят получить.
Интересная дискуссия, всем спасибо. Я на самом деле знаком с моделями Пуассона и рассматривал модель. Будет ли правильная спецификация модели summary(m1 <- glm(count ~ month + year, family = "poisson", data=Birthcount)) ? Кроме того, спасибо всем за вашу помощь.
да, это было бы более «технически» правильным, но я не думаю, что предсказания были бы другими. Возможно, вы также захотите изучить моделирование этого как временного ряда, как предлагает @dipetkov.
Так в чем именно ваш вопрос? Если вам нужен совет по созданию модели прогнозирования, вам действительно следует обратиться за помощью в Перекрестная проверка, где обсуждаются статистические вопросы. Это не похоже на конкретный вопрос программирования, который подходит для переполнения стека.