我经常使用 B 样条进行回归。到目前为止,我从未需要了解 bs
的输出详细说明:我会选择我感兴趣的型号,然后搭配 lm
.但是,我现在需要在外部(非 R)代码中重现 b 样条模型。那么,bs
生成的矩阵是什么意思? ?示例:
x <- c(0.0, 11.0, 17.9, 49.3, 77.4)
bs(x, df = 3, degree = 1) # generate degree 1 (linear) B-splines with 2 internal knots
# 1 2 3
# [1,] 0.0000000 0.0000000 0.0000000
# [2,] 0.8270677 0.0000000 0.0000000
# [3,] 0.8198433 0.1801567 0.0000000
# [4,] 0.0000000 0.7286085 0.2713915
# [5,] 0.0000000 0.0000000 1.0000000
# attr(,"degree")
# [1] 1
# attr(,"knots")
# 33.33333% 66.66667%
# 13.30000 38.83333
# attr(,"Boundary.knots")
# [1] 0.0 77.4
# attr(,"intercept")
# [1] FALSE
# attr(,"class")
# [1] "bs" "basis" "matrix"
好的,所以
degree
是 1,正如我在输入中指定的那样。 knots
告诉我两个内部结分别位于 x = 13.3000 和 x = 38.8333。看到结在固定分位数上有点惊讶,我希望 R 会为我的数据找到最好的分位数,但这当然会使模型不是线性的,而且如果不知道响应数据也是不可能的。 intercept = FALSE
意味着没有在基础中包含截距(这是一件好事吗?我一直被教导不要在没有截距的情况下拟合线性模型......好吧猜测 lm
无论如何只是添加一个)。然而,矩阵呢?我真的不明白如何解释它。三列,我认为这意味着基函数是三个。这是有道理的:如果我有两个内部结
K1
和 K2
,我将在左边界结B1
之间有一个样条线和 K1
, K1
之间的另一个样条和 K2
,以及 K2
之间的最后一个和 B2
,所以...三个基函数,好的。但到底哪些是基函数呢?例如,此列是什么意思?# 1
# [1,] 0.0000000
# [2,] 0.8270677
# [3,] 0.8198433
# [4,] 0.0000000
# [5,] 0.0000000
编辑:这与 this question 相似但不完全相同.该问题询问回归系数的解释,但我在此之前的一步:我想了解模型矩阵系数的含义。如果我尝试制作与 the first answer 中建议的相同的图,我得到了一个困惑的情节:
b <- bs(x, df = 3, degree = 1)
b1 <- b[, 1] ## basis 1
b2 <- b[, 2] ## basis 2
b3 <- b[,3]
par(mfrow = c(1, 3))
plot(x, b1, type = "l", main = "basis 1: b1")
plot(x, b2, type = "l", main = "basis 2: b2")
plot(x, b3, type = "l", main = "basis 3: b3")
这些不能是 B 样条基函数,因为它们有太多节点(每个函数应该只有一个)。
The second answer实际上可以让我在 R 之外重建我的模型,所以我想我可以这样做。然而,这个答案也没有完全解释
b
的元素是什么。矩阵是:它处理线性回归的系数,这里我还没有介绍。这确实是我的最终目标,但我也想了解这个中间步骤。
最佳答案
矩阵b
# 1 2 3
# [1,] 0.0000000 0.0000000 0.0000000
# [2,] 0.8270677 0.0000000 0.0000000
# [3,] 0.8198433 0.1801567 0.0000000
# [4,] 0.0000000 0.7286085 0.2713915
# [5,] 0.0000000 0.0000000 1.0000000
实际上只是
x
每个点的三个基函数的值的矩阵,这对我来说应该是显而易见的,因为它与多项式线性模型的解释完全相同。事实上,由于边界节点是bknots <- attr(b,"Boundary.knots")
# [1] 0.0 77.4
内部结是
iknots <- attr(b,"knots")
# 33.33333% 66.66667%
# 13.30000 38.83333
然后是三个基函数,如图here , 是:
knots <- c(bknots[1],iknots,bknots[2])
y1 <- c(0,1,0,0)
y2 <- c(0,0,1,0)
y3 <- c(0,0,0,1)
par(mfrow = c(1, 3))
plot(knots, y1, type = "l", main = "basis 1: b1")
plot(knots, y2, type = "l", main = "basis 2: b2")
plot(knots, b3, type = "l", main = "basis 3: b3")
现在,考虑
b[,1]
# 1
# [1,] 0.0000000
# [2,] 0.8270677
# [3,] 0.8198433
# [4,] 0.0000000
# [5,] 0.0000000
这些必须是
b1
的值在 x <- c(0.0, 11.0, 17.9, 49.3, 77.4)
.事实上,b1
在 knots[1] = 0
中为 0和 1 英寸 knots[2] = 13.3000
,这意味着在 x[2]
(11.0) 值必须是 11/13.3 = 0.8270677
,正如预期的那样。同样,由于 b1
knots[3] = 38.83333
为 0 , x[3]
中的值(17.9) 必须是 (38.83333-13.3)/17.9 = 0.8198433
.自 x[4], x[5] > knots[3] = 38.83333
, b1
那里是 0。可以对其他两列给出类似的解释。
关于r - R 函数 bs() 输出的解释(B 样条基矩阵),我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/39446842/