1

我正在使用 splines 包的 bs 函数来创建用于图形目的的 b 样条平滑曲线。（至少有一份报告称 Excel 使用三阶 b 样条曲线作为其平滑线图，我希望能够复制这些曲线。）我无法理解 bs 函数所需的参数。代表性代码如下，改编自 bs 文档：

require(splines)
require(ggplot2)
n <- 10
x <- 1:10
y <- rnorm(n)
d <- data.frame(x=x, y=y)
summary(fm1 <- lm(y ~ bs(x, degree=3)), data=d)
x.spline <- seq(1, 10, length.out=n*10)
spline.data <- data.frame(x=x.spline, y=predict(fm1, data.frame(x=x.spline)))
ggplot(d, aes(x,y)) + geom_point + geom_line(aes(x,y), data=spline.data)

bs 文档中的示例代码在对 bs 的调用中指定了 df=5，并且没有指定度数。我不知道我有多少自由度。我所知道的是我想要一个三阶 b 样条。我已经尝试指定不同的 df 值来代替度数，或者除了度数之外，我得到了截然不同的结果。这就是为什么我怀疑 df 的规范是这里的问题。在这种情况下我将如何计算 df ？

帮助文件建议 df = length(knots) + degree。如果我将内部点视为结，则在此示例中为我提供 df=11，这会生成错误消息和无意义的样条拟合。

先感谢您。

我显然不清楚我的意图。我正在尝试这样做： 如何将 spline() 与 ggplot 一起使用？，但使用 b 样条。

4

2 回答 2

3

你不应该试图适应每一点。目标是找到一个可接受的拟合但取决于有限数量的结的摘要。将多项式的 hte 次数增加到默认值三以上并没有太大价值。只有 10 分，您肯定不希望 df=11。尝试 df=5 ，结果应该相当平坦。rms/Hnisc 包的作者 Frank Harrell 更喜欢限制三次样条，因为在极端情况下的预测是线性的，因此比其他多项式基的预测要少。

我更正了几个拼写错误并添加了一个knots参数以使您的代码正常工作：

require(splines)
require(ggplot2); set.seed(trunc(100000*pi))

n <- 10
x <- 1:10
y <- rnorm(n)
d <- data.frame(x=x, y=y)
summary(fm1 <- lm(y ~ bs(x, degree=3, knots=2)), data=d)
x.spline <- seq(1, 10, length.out=n*10)
spline.data <- data.frame(x=x.spline, y=predict(fm1, data.frame(x=x.spline)))
ggplot(d, aes(x,y)) + geom_point() + geom_line(aes(x,y), data=spline.data)

我从改变随机种子的练习中走出来，认为弗兰克哈雷尔知道他在说什么。在使用他的包时，我不会在极端情况下得到同样的行为。

于 2012-04-27T23:01:29.200 回答
0

我做了更多的工作，并提出了以下建议。首先，道歉。我正在寻找的是平滑样条，而不是回归样条。我没有正确表达问题的词汇。虽然 bs() 帮助文件中的示例似乎提供了这一点，但该函数并未为我的示例数据提供相同的行为。在 stats 包中还有另一个函数，smooth.spline，它提供了我需要的东西。

set.seed(tunc(100000*pi))
n <- 10
x <- 1:n
xx <- seq(1, n, length.out=200)
y <- rnorm(n)
d <- data.frame(x=x, y=y)
spl <- smooth.spline(x,y, spar=0.1)
spline.data <- data.frame(y=predict(spl,xx))
ggplot(d,aes(x,y)) + geom_point() + geom_line(aes(x,y), spline.data)
spl2 <- smooth.spline(x, y, control=
            list(trace=TRUE, tol=1e-6, spar=0.1, low=-1.5, high=0.3))
spline.data2 <- data.frame(predit(spl2,xx))
ggplot(d,aes(x,y)) + geom_point() + geom_line(aes(x,y), spline.data2)

对 smooth.spline 的两次调用代表了两种方法。第一个手动指定平滑参数，第二个迭代到最优解。我发现我必须适当地限制优化以获得我所追求的解决方案类型。

结果旨在与 Excel 线图使用的 b 样条相匹配。我有一些合作者认为 Excel 图形是标准，我至少需要匹配这种性能。

于 2012-05-03T18:52:31.987 回答