B样条混乱

Question

我知道这个讨论板块上有关于B样条的话题，但是这些话题更让我感到困惑了，所以我想有人能帮助我。

我的仿真数据的x值范围从0到1。我想用B样条基函数和OLS参数估计（不是惩罚性样条）拟合我的数据，节点为0、0.1、0.2......0.9、1，度数为3。

我认为我需要来自spline软件包的bs函数，但我不太确定，也不知道要输入什么内容。

我还想绘制出结果多项式样条。

谢谢！

- user2249626

---

“B-spline confusion”这个标题似乎很贴切。你怎么可能有10个节点和degree=3的三次样条曲线？ - IRTFM

@DWin那恰好就是bs的功能，不是吗？默认情况下，在结点上用三次立方多项式拟合，并确保各个片段在结点处平滑连接。 - Gavin Simpson

我想为了更好地理解问题。我认为被要求的是在整个数据范围内进行立方多项式拟合。否则，只需对“?bs”页面上的示例代码进行微小修改即可完全回答问题：lm(weight ~ bs(height, df = 3, knots=c(58, 62, 66, 70, 72), ), data = women)。 - IRTFM

Answer 1

## simulate some data - from mgcv::magic
set.seed(1)
n <- 400
x <- 0:(n-1)/(n-1)
f <- 0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10
y <- f + rnorm(n, 0, sd = 2)

## load the splines package - comes with R
require(splines)

如果你想要普通最小二乘估计，你可以在公式中使用bs()函数作为lm。 bs函数提供由节点、多项式次数等确定的基础函数。

mod <- lm(y ~ bs(x, knots = seq(0.1, 0.9, by = 0.1)))

你可以像处理线性模型一样处理它。

> anova(mod)
Analysis of Variance Table

Response: y
                                        Df Sum Sq Mean Sq F value    Pr(>F)    
bs(x, knots = seq(0.1, 0.9, by = 0.1))  12 2997.5 249.792  65.477 < 2.2e-16 ***
Residuals                              387 1476.4   3.815                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

关于结点位置的一些指针。 bs有一个参数Boundary.knots，默认为Boundary.knots = range(x) - 因此当我在上面指定knots参数时，我没有包括边界结点。

阅读?bs获取更多信息。

生成拟合样条的图形

在评论中，我讨论了如何绘制拟合的样条。一种选择是按照协变量的顺序对数据进行排序。这对单个协变量可以正常工作，但对于2个或更多的协变量可能不起作用。另一个问题是您只能评估观察值的x上的拟合样条 - 如果您没有密集采样协变量，则这是可以接受的，但如果没有，样条可能看起来很奇怪，具有长线性部分。

更通用的解决方案是使用predict为协变量的新值生成模型预测。在下面的代码中，我展示了如何为x范围内100个均匀间隔的值进行预测，以对上述模型进行预测。

pdat <- data.frame(x = seq(min(x), max(x), length = 100))
## predict for new `x`
pdat <- transform(pdat, yhat = predict(mod, newdata = pdat))

## now plot
ylim <- range(pdat$y, y) ## not needed, but may be if plotting CIs too
plot(y ~ x)
lines(yhat ~ x, data = pdat, lwd = 2, col = "red")

那会产生：

- Gavin Simpson

---

非常感谢，这帮了我很多！现在，如果我使用points(x,fitted(mod))绘制拟合值，我就能得到我想要的结果。然而，使用lines(x,fitted(mod))无法连接点以显示多项式样条。我该如何绘制结果的样条曲线呢？ - user2249626

抱歉，我是一个新的 R 用户。您所说的“排序”是什么意思？是将 x,y 对放入矩阵中，并按照 x 列对该矩阵进行排序吗？ - user2249626

我稍微考虑了一下。我按照以下方式对数据框进行排序：df <- data.frame(x,y) df_ordered <- df[order(df$x),]. 然后我像这样拟合模型：mod <- lm(df_ordered$y~bs(df_ordered$x,knots=seq(0.1,0.9, by=0.1))). 接着，我将样条线添加到散点图中：lines(df_ordered$x,fitted(mod)) 这样做是正确的吗？ - user2249626

@user2249626 嗯，那是一种方法。更好的方法是在协变量范围内使用 predict。 - Gavin Simpson

@user2249626 我已经添加了一个绘制样条线的示例。 - Gavin Simpson

Answer 2

基于答案中的示例，绘制拟合样条的一种更简单的方法是使用effects包。

## simulate some data - from mgcv::magic
set.seed(1)
n <- 400
x <- 0:(n-1)/(n-1)
f <- 0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10
y <- f + rnorm(n, 0, sd = 2)

## load the splines package - comes with R
require(splines)
require(car)
require(effects)

## estimate model
mod <- lm(y ~ bs(x, knots = seq(0.1, 0.9, by = 0.1)))

然后您可以使用car中的Anova：

> Anova(mod)
Anova Table (Type II tests)

Response: y
                                       Sum Sq  Df F value    Pr(>F)    
bs(x, knots = seq(0.1, 0.9, by = 0.1)) 2997.5  12  65.477 < 2.2e-16 ***
Residuals                              1476.4 387                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

你可以使用effects包轻松绘制拟合的样条曲线。

plot(allEffects(mod))

以下代码将输出：

另请参阅：

• 固定效应下展示logit的边际效应

- landroni