Добавление geom_sf к длинному графику с пакетом gratia в R

Разработчик Gratia утверждает, что функция draw():

никогда не предназначался для такой настройки, которая возможна с ggplot() или некоторыми другими пакетами, которые используют ggplot() в качестве слоя построения.

Поэтому я сомневаюсь, что существует способ построить научно-фантастический объект напрямую с помощью draw(). Однако ниже приведены два возможных решения. Обратите внимание: если вы хотели, чтобы данные вашего примера перекрывались с вашим объектом Ватикана, вы неправильно назначили координаты долготы и широты. Это исправлено в примере данных, используемых в этом репрексе.

Вариант 1. Преобразование объекта SF в фрейм данных

• Плюсы: Менее сложный из двух вариантов.
• Минусы: не учитывает относительное соотношение x/y координат CRS.

library(mgcv)
library(dplyr)
library(sf)
library(gratia)
library(giscoR)
library(ggplot2)

set.seed(1) # For reproducibility
sound <- sample(0:5, size = 960, replace = TRUE)
word <- as.factor(rep(c("alpha", "bravo", "charlie", "delta", "echo", "foxtrot"), each = 4))
age <- as.factor(rep(c("young", "old"), times = 480))
gender <- as.factor(rep(c("female", "female", "male", "male"), times = 240))
longitude <- rep(c(runif (40, min = 12, max = 13)), each = 24)
latitude <- rep(c(runif (40, min = 41, max = 42)), each = 24)
pronunciation <- data.frame(sound, word, age, gender, longitude, latitude)

model <- gam(list(sound ~ word + s(longitude, latitude),
                  ~ word + s(longitude, latitude),
                  ~ word + s(longitude, latitude),
                  ~ word + s(longitude, latitude),
                  ~ word + s(longitude, latitude)),
             data = pronunciation,
             family = multinom(K = 5))

vatican <- gisco_get_countries(resolution = "10", country = "VAT") |>
  mutate(res = "10M")

# Convert vatican sf object to df and add longitide and latitude columns (X,Y)
vatican_df <- vatican |>
  select(geometry) |>
  st_cast("POLYGON") |>
  st_coordinates() |>
  as.data.frame()

# Plot
draw(model, select = "s(longitude,latitude)") + 
  geom_polygon(data = vatican_df,
               aes(X, Y),
               fill = "yellow",
               colour = "yellow",
               linewidth = 2) +  # linewidth exaggerated for illustrative purposes
  labs(title = "Title")

Вариант 2. Извлеките данные из объекта ggplot(), созданного draw(), и постройте график с помощью ggplot().

• Плюсы: легче управлять отдельными геометриями, если вы знакомы с ggplot(), можно использовать geom_sf() и, следовательно, поддерживать соотношение координат x/y sf CRS.
• Минусы: более сложный, может потребоваться метод проб и ошибок для извлечения правильных слоев.

# Create plot object
p <- draw(model, select = "s(longitude,latitude)")

# Extract plot data from ggplot2 object
p_df <- ggplot_build(p) 
p_df1 <- p_df$plot$data # Estimates / partial effect
p_df2 <- p_df$data[[2]] # Line object layer
p_df3 <- p_df$data[[3]] # Panel with colour values
p_df4 <- p_df$data[[4]] # Point object layer

# Plot
ggplot() +
  geom_tile(data = p_df1,
            aes(longitude, latitude, fill = .estimate)) +
  geom_line(data = p_df2,
             aes(x, y, group = piece),
             size = 0.3) +
  geom_point(data = p_df4,
             aes(x, y),
             size = 2) +
  geom_sf(data = vatican,
          fill = "yellow",
          colour = "yellow",
          linewidth = 2) +
  labs(title = "Title", caption = paste("Basis:", p_df1[1,".type"])) +
  scale_fill_gradient2(low = p_df3[1,1], midpoint = 0, mid = "white", high = p_df3[2,1],
                       name = "Partial\neffect")

Цветовая шкала не идентична варианту 1, но вы можете установить свой собственный или встроенный градиент, если он должен совпадать.