我在使用R语言。我有25个变量,跨越15个时间点,每个变量在每个时间点有3个或更多的重复测量。我已经使用melt
函数将这些数据转换为data.frame
格式,我可以使用ggplot的facet_wrap()
命令等来愉快地绘制这些数据。我的转换后的数据框名为lis
;以下是它的头部和尾部数据,以便你了解数据的结构:
> head(lis) time variable value1 10 SELL 8.1704682 10 SELL 8.2158923 10 SELL 8.2142464 15 SELL 8.9106545 15 SELL 7.9285376 15 SELL 8.777784> tail(lis) time variable value145 1 GAS5 10.92248146 1 GAS5 11.37983147 1 GAS5 10.95310148 1 GAS5 11.60476149 1 GAS5 11.69092150 1 GAS5 11.70777
我可以使用以下ggplot2命令绘制所有时间序列的美丽图表,并添加拟合的样条曲线和95%的置信区间:
p <- ggplot(lis, aes(x=time, y=value)) + facet_wrap(~variable)p <- p + geom_point() + stat_smooth(method = "lm", formula = y ~ ns(x,3))
问题是这个平滑器不符合我的要求——95%的置信区间偏差很大。我希望使用高斯过程(GP)来获得更好的回归和时间序列的协方差估计。
我可以使用类似下面的代码来拟合一个高斯过程:
library(tgp) out <- bgp(X, Y, XX = seq(0, 200, length = 100))
这个函数接受时间X
、观测值Y
,并在XX
中的每个点进行预测。对象out
包含了关于这些预测的许多信息,包括一个协方差矩阵,我可以用它来替代从ns()
得到的95%置信区间(我认为是这样?)。
问题是我不知道如何包装这个函数以使其与ggplot2::stat_smooth()
接口。任何关于如何进行的建议或指导将不胜感激!
回答:
如文档中所述,stat_smooth
有y
、ymin
和ymax
美学属性,你可以使用这些属性来应用自定义平滑器,文档链接在这里:http://had.co.nz/ggplot2/stat_smooth.html。你可以创建一个包含自定义平滑器预测和置信区间的数据框,并直接在stat_smooth
中使用它(指定一个新的数据参数)。你可能可以使用stat_smooth(method="tgp::bgp",XX=seq(0,200,length=100))
,但我还没有尝试过。