Search code examples
rggplot2tidyversegam

Altering gam (scam) plotting options in geom_smooth across facets


I'm working on something where I want to plot data and put a monotone increasing or decreasing line of 'best fit' over it depending on the sign of Kendall's tau-b correlation. It works like a charm when I don't have to facet using a simple ifelse(). Unfortunately, when I have to facet based on groups (reproducible example below), I can't figure out how to vary the geom_smooth to adjust per facet. I tried (1) inputting a vector of bs (naive, but hopeful); (2) the process described here for switching model types, but neither worked.

The goal is to make this generalizable so I can plug in any dataset and produce a similar plot. While the solution below works, it would be somewhat clunky for that purpose (e.g., a for loop that adds a geom_smooth for each unique level of the facet variable). I hope I'm missing an easy fix!

library(tidyverse)
library(palmerpenguins)
library(scam)

p <- ggplot(data = penguins, aes(x=bill_depth_mm, y=body_mass_g))+
  geom_point()+
  theme_bw()+
  xlab("Bill Depth (mm)")+
  ylab("Body Mass (g)") +
  facet_wrap(~island) 

p <- p + geom_smooth(
  data=subset(penguins, island=="Biscoe"),
  color="black",
  method=scam,
  formula = y ~ s(x, bs="mpd"))

p<- p + geom_smooth(
  data=subset(penguins, island=="Dream"),
  color="black",
  method=scam,
  formula = y ~ s(x, bs="mpi"))

p<- p + geom_smooth(
  data=subset(penguins, island=="Torgersen"),
  color="black",
  method=scam,
  formula = y ~ s(x, bs="mpi"))

p

Interestingly, when I write the loop, it appears the bs method is "overwritten" and I don't get the same plot I do when I brute force it.

p <- ggplot(data = penguins, aes(x=bill_depth_mm, y=body_mass_g))+
  geom_point()+
  theme_bw()+
  xlab("Bill Depth (mm)")+
  ylab("Body Mass (g)") +
  facet_wrap(~island) 

for(facet in unique(penguins$island)){
  subdat<-subset(penguins, island==facet)
  p <- p + 
    geom_smooth(data=subdat,
                color="black",
                method=scam,
                formula = y ~ s(x, 
                                bs=ifelse(cor(subdat$bill_depth_mm, subdat$body_mass_g, method = "kendall")>0, 
                                          "mpi", "mpd")))
  
}

p

Solution

  • Whilst you could write your own wrapper for scam, I would probably do this by mapping a vector of bs as follows:

    library(tidyverse)
    library(palmerpenguins)
    library(scam)
    
    ggplot(data = penguins, aes(x = bill_depth_mm, y = body_mass_g)) +
      geom_point() +
      theme_bw() +
      xlab("Bill Depth (mm)") +
      ylab("Body Mass (g)") +
      facet_wrap(~island) +
      map2(levels(penguins$island), c('mpd', 'mpi', 'mpi'),
           ~ geom_smooth(data = subset(penguins, island == .x),
                         color = 'black', 
                         method = scam,
                         formula = substitute(y ~ s(x, bs = bs), env = list(bs = .y))))
    

    enter image description here