R + ggplot : how to use a custom smoother (Gaussian Process)

Question

I'm using R. I have 25 variables over 15 time points, with 3 or more replicates per variable per time point. I've melted this into a data.frame, which I can plot happily using (amongst other things) ggplot's facet_wrap() command. My melted data frame is called lis; here's its head and tail, so you get an idea of the data:

> head(lis)
  time variable    value
1   10     SELL 8.170468
2   10     SELL 8.215892
3   10     SELL 8.214246
4   15     SELL 8.910654
5   15     SELL 7.928537
6   15     SELL 8.777784
> tail(lis)
    time variable    value
145    1     GAS5 10.92248
146    1     GAS5 11.37983
147    1     GAS5 10.95310
148    1     GAS5 11.60476
149    1     GAS5 11.69092
150    1     GAS5 11.70777

I can get a beautiful plot of all the time series, along with a fitted spline and 95% confidence intervals using the following ggplot2 commands:

p <- ggplot(lis, aes(x=time, y=value)) + facet_wrap(~variable)
p <- p + geom_point() + stat_smooth(method = "lm", formula = y ~ ns(x,3))

The trouble is that the smoother is not to my liking - the 95% confidence intervals are way off. I would like to use Gaussian Processes (GP) to get a better regression and estimate of covariance for my time series.

I can fit a GP using something like

library(tgp) 
out <- bgp(X, Y, XX = seq(0, 200, length = 100))

which takes time X, observations Y and makes predictions at each point in XX. The object out contains a bunch of things about those predictions, including a covariance matrix I can use in place of the 95% confidence interval I get (I think?) from ns().

The trouble is I'm not how to wrap this function to make it interface with ggplot2::stat_smooth(). Any ideas or pointers as to how to proceed would be greatly appreciated!

Answer 1

Stat_smooth has y, ymin, and ymax aesthetics that you can use with a custom smoother, as documented here: http://had.co.nz/ggplot2/stat_smooth.html. You create a data frame with the predictions and CI from your custom smoother and use that directly in stat_smooth (specifying a new data argument). You might be able to use stat_smooth(method="tgp::bgp",XX=seq(0,200,length=100)) but I haven't tried it.

Answer 2

It looks like bgp doesn't follow the standard R style for modelling functions. This means that you can't use it inside geom_smooth and you'll need to fit the model outside of the ggplot2 call. You might also want to email the tgp package author and encourage them to follow R standards.

Answer 3

Dear Mike,

How are you planning to use your replicated observations with bgp? It only allows a single vector of observations and I cannot figure out how to do it as I have replicated vectors.