4

Sorry if this question is trivial, but I'm trying to figure out how to plot a certain type of natural cubic spline (NCS) in R and it's completely eluded me.

In a previous question I learned how to plot the NCS generated by the ns() command in ggplot, but I'm interested in how to plot a slightly different NCS generated the smooth.Pspline command in the pspline package. As far as I know this is the only package that automatically selects the proper smoothing penalty by CV for a given dataset.

Ideally I would be able to provide smooth.Pspline as a method to a stat_smooth layer in ggplot2. My current code is like:

plot <- ggplot(data_plot, aes(x=age, y=wOBA, color=playerID, group=playerID))
plot <- plot + stat_smooth(method = lm, formula = y~ns(x,4),se=FALSE)

I'd like to replace the "lm" formula with smooth.Pspline's functionality. I did a little bit of googling and found a solution to the very similar B-spline function smooth.spline, written by Hadley. But I haven't been able to adapt this to smooth.Pspline perfectly. Does anyone have experience with this?

Thanks so much!

Community
  • 1
  • 1
gogurt
  • 811
  • 1
  • 8
  • 24

1 Answers1

8

You simply need to inspect how predict.smooth.Pspline returns the predicted values.

In the internal workings of stat_smooth, predictdf is called to create the smoothed line. predictdf is an internal (non-exported) function of ggplot2 (it is defined here) it is a standard S3 method.

sm.spline returns an object of class smooth.Pspline, therefore for stat_smooth to work you need to create method for predictdf for class smooth.Pspline.

As such the following will work.

smP <- function(formula,data,...){
  M <- model.frame(formula, data)
  sm.spline(x =M[,2],y =M[,1])

}
# an s3 method for predictdf (called within stat_smooth)
predictdf.smooth.Pspline <- function(model, xseq, se, level) {
  pred <- predict(model, xseq)
  data.frame(x = xseq, y = c(pred))
}

An example (with a pspline fitted using mgcv::gam as comparison). mgcv is awesome and gives great flexibility in fitting methods and smoothing spline choices (although not CV, only GCV/UBRE/REML/ML)

d <- ggplot(mtcars, aes(qsec, wt))
d + geom_point() +  stat_smooth(method = smP, se= FALSE, colour='red', formula = y~x) + 
stat_smooth(method = 'gam', colour = 'blue', formula = y~s(x,bs='ps'))

enter image description here

mnel
  • 113,303
  • 27
  • 265
  • 254
  • 1
    Can you explain what function or class dispatch principles result in `predictdf.smooth.Pspline` being called? And I didn't get it when Hadley wrote it and I don't see the connection here, either. I see that stat_smooth is supposed to return the same values as you are setting for the arguments of that function, but what determines how it is named or gets called? – IRTFM Oct 10 '13 at 01:20
  • I think I get it now. There are several built-in `predictdf` functions that will match a limited number of 'method's, but all you need to do to add one is match the name with the class returned by the 'method' argument. – IRTFM Oct 10 '13 at 02:09
  • Hi @mnel: Can I ask you a somewhat-stupid follow-up question? In the body of the `smP` method, is it possible to make available to `sm.spline(...)` a different column of the original dataset I'm plotting? In particular, I want to specify a weight to each point I'm plotting, and that weight is stored in a separate column of my dataset. It should be accessible from `M<-model.frame(formula,data)`, right? – gogurt Oct 10 '13 at 03:40
  • @gogurt not easily, no. It might be better to use something like `dlply` or `by` to create the fits.... – mnel Oct 10 '13 at 03:45
  • Hmm. That is a problem. Well, at least we're **most** of the way there. Thanks so much again! – gogurt Oct 10 '13 at 04:03