Как построить вектор загруженных склонов в ggplot2?

Question

Я использовал ggplot2 для построения графиков результатов начальной загрузки различных статистических данных, таких как коэффициенты корреляции.Совсем недавно я запустил наклон модели линейной регрессии.Вот как это выглядит при использовании функции plot() из пакета graphics:

plot(main="Relationship Between Eruption Length at Wait Time at \n 
 Old Faithful With Bootstrapped Regression Lines", 
 xlab = "Eruption Length (minutes)", 
 ylab = "Wait Time (minutes)", 
 waiting ~ eruptions, 
 data = faithful, 
 col = spot_color, 
 pch = 19)

index <- 1:nrow(faithful)
for (i in 1:10000) {
    index_boot <- sample(index, replace = TRUE) #getting a boostrap sample (of indices) 
    faithful_boot <- faithful[index_boot, ]
    # Fitting the linear model to the bootstrapped data:
    fit.boot <- lm(waiting ~ eruptions, data = faithful_boot)
    abline(fit.boot, lwd = 0.1, col = rgb(0, 0.1, 0.25, alpha = 0.05)) # Add line to plot
}
fit <- lm(waiting ~ eruptions, data=faithful)
abline(fit, lwd = 2.5, col = "blue")

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

Например, функция может выглядеть примерно так:

set.seed(777) # included so the following output is reproducible
n_resample <- 10000 # set the number of times to resample the data

# First argument is the data; second is the number of resampled datasets
bootstrap <- function(df, n_resample) {
    slope_resample <- matrix(NA, nrow = n_resample) # initialize vector 
    index <- 1:nrow(df) # create an index for supplied table

    for (i in 1:n_resample) {
        index_boot <- sample(index, replace = TRUE) # sample row numbers, with replacement
        df_boot <- df[index_boot, ] # create a bootstrap sample from original data
        a <- lm(waiting ~ eruptions, data=df_boot) # compute linear model
        slope_resample[i] <- slope <- a$coefficients[2] # take the slope
    }
    return(slope_resample) # Return a vector of differences of proportion
}

bootstrapped_slopes <- bootstrap(faithful, 10000)

Но как получить geom_line() или geom_smooth() для получения данных из bootstrapped_slopes?Любая помощь очень ценится.

Answer 1

РЕДАКТИРОВАТЬ: Более прямая адаптация из OP

Для построения графика я предполагаю, что вам нужны и уклоны, и пересечения, так что вот модифицированная функция bootstrap:

bootstrap <- function(df, n_resample) {
  # Note 2 dimensions here, for slope and intercept
  slope_resample <- matrix(NA, 2, nrow = n_resample) # initialize vector 
  index <- 1:nrow(df) # create an index for supplied table

  for (i in 1:n_resample) {
    index_boot <- sample(index, replace = TRUE) # sample row numbers, with replacement
    df_boot <- df[index_boot, ] # create a bootstrap sample from original data
    a <- lm(waiting ~ eruptions, data=df_boot) # compute linear model
    slope_resample[i, 1] <- slope <- a$coefficients[1] # take the slope
    slope_resample[i, 2] <- intercept <- a$coefficients[2] # take the intercept
  }
  # Return a data frame with all the slopes and intercepts 
  return(as.data.frame(slope_resample))      
}

Тогдазапустите его и построите линии из этого фрейма данных:

bootstrapped_slopes <- bootstrap(faithful, 10000)

library(dplyr); library(ggplot2)
ggplot(faithful, aes(eruptions, waiting)) +
  geom_abline(data = bootstrapped_slopes %>% 
                sample_n(1000), # 10k lines look about the same as 1k, just darker and slower
              aes(slope =  V2, intercept = V1), #, group = id), 
              alpha = 0.01) +
  geom_point(shape = 19, color = "red")

---

Альтернативное решение

Это также можно сделать с помощью modelr и broom, чтобы упростить некоторые из процессов начальной загрузки,Основываясь на примере основной справки для modelr::bootstrap, мы можем сделать следующее:

library(purrr); library(modelr); library(broom); library(dplyr)
set.seed(777) 

# Creates bootstrap object with 10k extracts from faithful
boot <- modelr::bootstrap(faithful, 10000)
# Applies the linear regression to each
models <- map(boot$strap, ~ lm(waiting ~ eruptions, data = .))
# Extracts the model results into a tidy format
tidied <- map_df(models, broom::tidy, .id = "id")
# We just need the slope and intercept here
tidied_wide <- tidied %>% select(id, term, estimate) %>% spread(term, estimate) 

ggplot(faithful, aes(eruptions, waiting)) +
  geom_abline(data = tidied_wide %>%
     sample_n(1000), # 10k lines look about the same as 1k, just darker and slower
     aes(slope =  eruptions, intercept = `(Intercept)`, group = id), alpha = 0.05) +
  geom_point(shape = 19, color = "red") # spot_color wasn't provided in OP