Using the draw function from the gratia package with a model that contains a smooth s(longitude, latitude) will plot a long-lat plot and effect contours. That's very nice!

I want to add a country shape to the plot

library(giscoR)

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

Plotting the shape with ggplot works

ggplot() +
  geom_sf(data=vatican)

But with gratia not so much.

draw(model,
     select="s(longitude,latitude)") +
  geom_sf(data=vatican)

I get the error message:

Coordinate system already present. Adding a new coordinate system, which will replace the existing one.
Error in `geom_sf()`:
! Problem while computing aesthetics.
ℹ Error occurred in the 5th layer.
Caused by error in `.data[["longitude"]]`:
! Column `longitude` not found in `.data`.

I'd appreciate any help on how to solve this!

EDIT

The model to test is

sound <- sample(0:5, size=960, replace=T)
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=41, max=42)), each=24)
latitude <- rep(c(runif(40, min=12, max=13)), each=24)
pronunciation <- data.frame(sound, word, age, gender, longitude, latitude)

library(mgcv)

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))