Как я могу добавить второй затененный доверительный интервал (ie., 95% и 85% ДИ) в визуализацию glm с помощью plot_model ()? - PullRequest
1 голос
/ 16 июня 2020

Используя plot_model(), я визуализировал предельные эффекты для нескольких условий модели (см. Рисунок ниже). Функция plot_model() позволяет вам указать уровень достоверности и автоматически построить заштрихованный доверительный интервал.

Я хотел бы нанести второй доверительный интервал поверх первого, чтобы оба можно увидеть 95% и 85% (чтобы дополнительно продемонстрировать силу ответов). Возможно ли это?

Вот некоторые сокращенные фиктивные данные, которые демонстрируют мой вопрос:

response   per.grs
0          0.2430000
0          0.7142857
0          1.0000000
0          0.7619048
0          0.7619048
0          0.1230000
0          0.6666667
0          0.3560000
0          0.9523810
0          0.1450000
1          0.7619048
1          0.6432000
1          0.6666667
1          0.8571429
1          0.8571429
1          0.5238095
1          0.9523810
1          0.8450000
1          0.5714286
1          0.7619048

Вот код, который я использовал:

grass.cover <- glm(response ~ per.grs + I(per.grs^2), data=data, family=binomial, na.action = "na.fail")
plot_model(grass.cover, type = "eff", terms = "per.grs[all]", ci.lvl = .95)

enter image description here

1 Ответ

0 голосов
/ 16 июня 2020

Этого можно добиться вот так. Я заглянул внутрь plot_model. Данные для графика вычисляются через ggeffects::ggeffect. Поэтому я просто вычисляю данные для 85%, которые затем можно использовать для добавления второй заштрихованной области через geom_ribbon. Однако, поскольку это также перекрывает линию с прогнозами, мы должны добавить geom_line, чтобы вернуть прогнозы в график. Однако, особенно в случае, если вы хотите иметь легенду для обоих CI, я бы рекомендовал создать весь график самостоятельно через ggplot2. (;

data <- structure(list(response = c(
  0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
  0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L
), per.grs = c(
  0.243,
  0.7142857, 1, 0.7619048, 0.7619048, 0.123, 0.6666667, 0.356,
  0.952381, 0.145, 0.7619048, 0.6432, 0.6666667, 0.8571429, 0.8571429,
  0.5238095, 0.952381, 0.845, 0.5714286, 0.7619048
)), class = "data.frame", row.names = c(
  NA,
  -20L
))

grass.cover <- glm(response ~ per.grs + I(per.grs^2), data = data, family = binomial, na.action = "na.fail")

library(sjPlot)
#> Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!
library(ggplot2)

p <- plot_model(grass.cover, type = "eff", terms = "per.grs[all]", ci.lvl = .95)

# Get data for 85% CI
dat85 <- ggeffects::ggeffect(model = grass.cover, terms = "per.grs[all]", ci.lvl = .85)

p +
  # Add 85% CI
  geom_ribbon(data = dat85, aes(x, ymin = conf.low, ymax = conf.high), fill = "grey") +
  # Get the predictions back (;)
  geom_line(data = dat85, aes(x, y = predicted), color = "black")

...