Как направить приборные модели модели lm в ggplot2?

Question

У меня есть следующий код, который вычисляет для каждого года от 1961:2018 влияния обеих переменных-предикторов: основание на шарах за игру (BB) и хоум-раны за игру (HR) на переменные ответа, запускаемые за игру (R):

rm(list = ls())

library(dbplyr)
library(tidyverse)
library(broom)
library(Lahman)

fit <- Teams %>% 
  filter(yearID %in% 1961:2018) %>% 
  mutate(BB = BB / G, 
         HR = HR / G,
         R = R / G) %>%
  group_by(yearID) %>%
  do(tidy(lm(R ~ BB + HR, data = .), conf.int = TRUE)) %>% filter(term=="BB")
fit

> fit
# A tibble: 58 x 8
# Groups:   yearID [58]
   yearID term  estimate std.error statistic p.value conf.low conf.high
    <int> <chr>    <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
 1   1961 BB      0.0845     0.168     0.502   0.623  -0.274      0.443
 2   1962 BB      0.142      0.273     0.520   0.610  -0.434      0.718
 3   1963 BB      0.339      0.242     1.40    0.178  -0.171      0.849
 4   1964 BB     -0.105      0.302    -0.349   0.731  -0.742      0.532
 5   1965 BB      0.235      0.253     0.928   0.366  -0.299      0.768
 6   1966 BB      0.104      0.216     0.482   0.636  -0.351      0.559
 7   1967 BB      0.0660     0.223     0.296   0.771  -0.405      0.537
 8   1968 BB     -0.199      0.203    -0.983   0.340  -0.627      0.229
 9   1969 BB      0.153      0.163     0.942   0.357  -0.185      0.492
10   1970 BB      0.239      0.157     1.52    0.143  -0.0874     0.566
# ... with 48 more rows

Теперь я хотел бы вывести эту «подгонку», которая на самом деле представляет собой тиббл (или модернизированный фрейм данных), в ggplot, чтобы показать оценки за год в виде точек, а также линию регрессии вместе сCI вычисляется по модели lm, а не просто пересчитывает ее с geom_smooth(method = "lm").

. Я попробовал следующее безуспешно. Я знаю, что augment от метлы должен работать непосредственно с выходом модели lm, и поэтому следующий код неверен, но он иллюстрирует то, чего я пытаюсь достичь:

augment(fit) %>%
  ggplot() +
  geom_point(aes(yearID, estimate)) +
  geom_line(aes(yearID, .fitted), col = "blue")

Как я могу это сделатьэто без "обмана" (дважды вычисляя lm один раз, а затем и на ggplot) и делая:

fit %>% ggplot(aes(yearID,estimate)) + geom_point() + geom_smooth(method = "lm")

Answer 1

Я пошел по аналогичному маршруту к Патрику, используя map() и nest():

library(tidyverse)
library(broom)
library(Lahman)
library(magrittr)

fit <- Teams %>%
  filter(yearID %in% 1961:2018) %>%
  mutate(
    BB = BB / G,
    HR = HR / G,
    R  = R / G
  ) %>%
  nest(data = -yearID) %>%
  mutate(
    model  = map(data, ~ lm(R ~ BB + HR, .x)),  # apply model to all nested groups
    m_tidy = map(model, tidy),                  # tidy up
    est    = map_dbl(m_tidy, ~ .x %>%           # pull BB estimate from each group
      filter(term == "BB") %>%
      pull(estimate)),
  )

Теперь на этом этапе вы можете просто %$% перейти прямо в следующую часть, но я сохранил ихотделитесь здесь, так что говорите о том, как правильно имитировать доверительный интервал. Доверительный интервал geom_smooth() основан на t-распределении, а не на нормальном распределении. Таким образом, нам придется проделать небольшую дополнительную работу, чтобы получить интервалы для работы:

fit %$%
  lm(est ~ yearID) %>%
  augment() %>% 
  mutate(m.se.fit = .se.fit * qt(1 - (1-0.95)/2, nrow(fit))) %>% # 95% conf int calc
  ggplot(aes(yearID, est)) +
  geom_point() +
  geom_line(aes(y = .fitted), col = "blue") + 
  geom_ribbon(aes(ymin = .fitted - m.se.fit, ymax = .fitted + m.se.fit), alpha = .2)

Этот график по существу отражает желаемый график:

fit %>% ggplot(aes(yearID, est)) + 
  geom_point() + 
  geom_smooth(method = "lm")

^{Создано в 2019-10-23 для представительного пакета (v0.3.0)}

Answer 2

Вы можете попробовать map функции из пакета purrr, который включен в tidyverse. Возможный код для описанной проблемы приведен ниже. Также возможно с lapply, если вы не знакомы с пакетом purrr.

library(tidyverse)
library(broom)
library(Lahman)

fit <- Teams %>% 
  filter(yearID %in% 1961:2018) %>% 
  mutate(BB = BB / G, 
         HR = HR / G,
         R = R / G) %>%
  group_by(yearID) %>%
  # consolidate your data
  nest() %>% 
  # creates new nested column with your regression data
  mutate(model = map(data, function(df) 
    tidy(lm(R ~ BB + HR, data = df), conf.int = TRUE) %>%
      filter(term=="BB")
    ),
    # extract the column estimate
    model_est = map_dbl(model, function(df) 
      df %>% pull(estimate)
    ), 
    # extract the column conf.low
    model_conf.low = map_dbl(model, function(df) 
      df %>% pull(conf.low)
    ), 
    # extract the column conf.high
    model_conf.high = map_dbl(model, function(df) 
      df %>% pull(conf.high)
    )
   ) 


fit %>% ggplot(aes(yearID,model_est)) + geom_point() +
  geom_line(aes(yearID, model_conf.low)) + 
  geom_line(aes(yearID, model_conf.high))