Custom color gradient for multiple GAM plots using Gratia

Question

I fitted a GAM model which includes two tensor smooths, using mgcv package. All variables except 'season' are numeric in nature.

library(mgcv)
# Fit the GAM model
gam_model <- gam(
  score ~ te(duration, difference) + te(duration, temperature) + season,
  data = my_data
)

Then I plotted the model to see the partial effects, using draw() function from gratia package and it plots fine. However I want to change the color gradient such that it is red for low values and green for high values. I tried the below code for the same.

gam_plot <- gratia::draw(gam_model, contour = FALSE)
gam_plot + scale_fill_gradient(low = "red", high = "green")

However, the color gradient is getting applied only for the second tensor product and the first one has the default blue to red gradient.

I then tried splitting & plotting them based on the interaction name,

# Draw the first interaction term with custom color
plot1 <- gratia::draw(gam_model, terms = "te(duration, difference)", contour = FALSE)
plot1 <- plot1 + ggplot2::scale_fill_gradient(low = "red", high = "green")
print(plot1)

# Draw the second interaction term with custom color
plot2 <- gratia::draw(gam_model, terms = "te(duration, temperature)", contour = FALSE)
plot2 <- plot2 + ggplot2::scale_fill_gradient(low = "red", high = "green")
print(plot2)

But still, both print the two plots together as I had earlier i.e. one with color gradient changed and other without.

Moreover all the above plots come with a warning

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

Appreciate your help in getting both plots with red -> green gradient. Thank you

Answer 1

library("mgcv")
library("dplyr")
library("readr")
library("gratia")

data_url <- "https://github.com/gavinsimpson/intro-gam-webinar-2020/raw/master/data/gbtemp.csv"

galveston <- read_csv(data_url) |>
    mutate(datetime = as.POSIXct(paste(DATE, TIME), format = '%m/%d/%y %H:%M', tz = "CDT"),
           STATION_ID = factor(STATION_ID), DoY = as.numeric(format(datetime, format = '%j')),
           ToD = as.numeric(format(datetime, format = '%H')) +
               (as.numeric(format(datetime, format = '%M')) / 60))

knots <- list(DoY = c(0.5, 366.5))
m <- bam(MEASUREMENT ~
             s(ToD, k = 10) +
             s(DoY, k = 12, bs = 'cc') +
             s(YEAR, k = 30) +
             s(LONGITUDE, LATITUDE, k = 100, bs = 'ds', m = c(1, 0.5)) +
             ti(DoY, YEAR, bs = c('cc', 'tp'), k = c(12, 15)) +
             ti(LONGITUDE, LATITUDE, ToD, d = c(2,1), bs = c('ds','tp'),
                m = list(c(1, 0.5), NA), k = c(20, 10)) +
             ti(LONGITUDE, LATITUDE, DoY, d = c(2,1), bs = c('ds','cc'),
                m = list(c(1, 0.5), NA), k = c(25, 12)) +
             ti(LONGITUDE, LATITUDE, YEAR, d = c(2,1), bs = c('ds','tp'),
                m = list(c(1, 0.5), NA), k = c(25, 15)),
         data = galveston, method = 'fREML', knots = knots,
         nthreads = c(4, 1), discrete = TRUE)

One way in general to change all the plots in the patchwork is to use the & operator:

draw(m, n = 25, rug = FALSE) &
  ggplot2::scale_fill_gradient(low = "red", high = "green")

but this is noisy (note all the messages from ggplot) as you are forcing the replacement of the scales. Instead, draw.gam() has arguments to change these aesthetics. In this case use continuous_fill:

draw(m, n = 25, rug = FALSE,
  continuous_fill = ggplot2::scale_fill_gradient(low = "red", high = "green") )

That said, I would rethink doing this. You are plotting a surface with a natural 0 point and as such a diverging colour palette would be most appropriate. With your red-green gradient it is much harder to judge the 0 point (which is the overall average of the data). Then there's the problems such a scale would cause for many individuals with certain colour vision deficiencies.