R + ggplot : how to use a custom smoother (Gaussian Process)
I'm using R. I have 25 variables over 15 time points, with 3 or more replicates per variable per time point. I've melt
ed 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 thin开发者_StackOverflow社区gs 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!
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.
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.
精彩评论