Plotting two smooth terms from different GAM's on a single plot in R-CodePudding

Is it possible to display smooth terms from different GAM's in R if those terms are estimated using the same form of data?

I have two ecological datasets with values of species diversity along an elevational gradient. One site ranges from 1700-2800 m a.s.l. while the other ranges from 2500-3800 m a.s.l. I've modelled the relationship of species diversity against elevation using gam in R and I would like to display the smooth term from each site's GAM in the same plot area. Something like this:

I know you can use compare_smooths to compare smooths between two GAMs but I was wondering if there was a more versatile approach.

Hopefully this question is fine without reproducible data.

CodePudding user response：

The easiest way would be to create predictions from the two models covering the elevational ranges you want, and then plot the fitted values.

For example (code not checked as I am AFK right now)

library("gratia")
library("dplyr")

newd <- data.frame(Elevation = c(seq(1700, 2800, length = 100), 
                                 seq(2500, 3800, length = 100)),
                   Site = rep(c("A", "B"), each = 100))
fv1 <- fitted_values(gam1, data = filter(newd, Site == "A"),
                     type = "response")
fv2 <- fitted_values(gam2, data = filter(newd, Site == "B"),
                     type = "response")
fv <- bind_rows(fv1, fv2)

Then you should have enough data in fv to plot, say

library("ggplot2")

fv %>% ggplot(aes(x = Elevation, colour = Site))  
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2)  
  geom_line(aes(y = fitted))  
  geom_point(data = foo, aes(x = Elevation, y = Diversity))

(where the last line is for adding the raw data, but I'm not sure how you have this arranged, so figure that bit out for yourself.)

CodePudding user response：

You can pass the subset of each site to a couple of different geom_smooth(). Something like this:

your_plot  
    geom_smooth(data = subset1) 
    geom_smooth(data = subset2)

I could be more precise if you provided a reprex.