Проверьте, находятся ли линии в регрессии CI R ggplot

Question

У меня есть набор данных, который я строю с линейной регрессией (всего на geom_smooth()). У меня также есть список уклонов и перехватов для других линий. Мне интересно, как лучше всего проверить, попадают ли эти другие линии в CI регрессии. В реальном наборе данных будет много строк для тестирования, поэтому я надеюсь, что есть какой-то способ массового тестирования.

library(tidyverse)    
# sample data
set.seed(0)
df = tibble(
  x = runif(20,min=5,max=10),
  y = x + rnorm(20,mean=2,sd=3)
)


# df of lines
lines.df <- data.frame(line = c("line 1","line 2","line 3","line 4"),
                       slope = c(.9,1.2,1.4,1.1),
                       intercept = c(2.5,-.5,-.44,0)) 


# plot
ggplot(df) +
  aes(x = x, y = y) +
  geom_point(size = 2) +
  geom_smooth(method = "lm", level = 0.95, fullrange = TRUE, color = "black")+
  geom_abline(data=lines.df,aes(slope=slope,intercept=intercept,color=line))

Мы видим выше, что строки 1, 2 и 4 находятся в пределах CI этой регрессии, а строка 3 - нет. Как я могу проверить это (визуально или иным образом), чтобы добавить столбец ay / n в мой файл lines.df относительно того, будут ли они находиться в этом доверительном интервале или нет?

Единственная идея, которую я имею до сих пор, - это создать регрессионная модель, рассчитайте CI, затем измените столбец в файле lines.df, экстраполируя эту строку и добавляя «если y (line1) при x = 0 находится между y (модель (0) ± CI (0)) & y (line1) при x = .5 находится между y (модель (.5) ± CI (.5)) & y (line1) при x = 1 находится между y (модель (1) ± CI (1)) «затем insideCI ==» да". Но это довольно неуклюже и глупо, поэтому я думаю, что должен быть лучший способ.

Для справки, это не обязательно должно быть графически, это просто самый простой способ объяснить, что я ищу.

Было бы лучше, если бы код был в синтаксисе dplyr / tidyverse, если это возможно.

Спасибо!

Answer 1

Полоса доверия, которую вы получаете из ggplot, является прогнозируемыми значениями подобранных +/- 1,96 * SE. Таким образом, вам нужно проверить для каждого прогнозируемого значения ваших строк, это <1,96 * SE. Чтобы проиллюстрировать это, SE (извините, не очень хорошо с ggplots): </p>

df = df[order(df$x),]
fit = lm(y~x,data=df)
pred=predict(fit,se=TRUE)
plot(df,pch=20)
lines(df$x,fit$fitted.values)
lines(df$x,pred$fit+1.96*pred$se.fit,lty=8)
lines(df$x,pred$fit-1.96*pred$se.fit,lty=8)
for(i in 1:nrow(lines.df)){
with(lines.df,abline(b=slope[i],a=intercept[i],col=terrain.colors(4)[i]))
}

Затем мы go просматриваем ваши данные, сначала сохраняя fit + SE:

library(broom)
res = augment(fit)
# A tibble: 20 x 9
       y     x .fitted .se.fit .resid   .hat .sigma  .cooksd .std.resid
   <dbl> <dbl>   <dbl>   <dbl>  <dbl>  <dbl>  <dbl>    <dbl>      <dbl>
 1  6.64  5.31    7.00   0.981 -0.366 0.204    2.23 0.00455      -0.189
 2  8.28  5.88    7.55   0.816  0.738 0.141    2.23 0.0110        0.366
 3  6.77  6.01    7.66   0.782 -0.890 0.130    2.22 0.0144       -0.439
 4  9.16  6.03    7.68   0.777  1.48  0.128    2.20 0.0388        0.728

А затем мы go через строки, используя purrr и tidyr, (заранее извиняюсь, поскольку вы предпочитаете решение dplyr / tidyverse, и я не так хорошо разбираюсь в те):

library(purrr)
library(tidyr)

lines.df %>% nest(param=c(slope,intercept)) %>%
# calculates values according to slopes 
mutate(pred = map(param,~.x$slope*df$x +.x$intercept),
# calculate the difference between these values and the actual fit 
       deviation_from_lm=map(pred,~abs(.x-res$.fitted)),
#check all of them within 1.96*se
       within=map_lgl(deviation_from_lm,~all(.x<=1.96*res$.se.fit))
)

# A tibble: 4 x 5
  line   param            pred       deviation_from_lm within
  <fct>  <list>           <list>     <list>            <lgl> 
1 line 1 <tibble [1 × 2]> <dbl [20]> <dbl [20]>        TRUE  
2 line 2 <tibble [1 × 2]> <dbl [20]> <dbl [20]>        TRUE  
3 line 3 <tibble [1 × 2]> <dbl [20]> <dbl [20]>        FALSE 
4 line 4 <tibble [1 × 2]> <dbl [20]> <dbl [20]>        TRUE