【问题标题】:Multivariate regression splines in RR中的多元回归样条
【发布时间】:2017-01-12 21:20:52
【问题描述】:

大多数人可能都熟悉样条曲线中的bs:

library(splines)

workingModel <- lm(mpg ~ factor(gear) + bs(wt, knots = 5) + hp, data = mtcars)

bs(mtcars$wt, knots = 4)

这对单变量权重使用 b 样条,但您也可以使用多元样条:

bs(cbind(mtcars$wt,mtcars$hp), knots = 4)

但这会产生一个行数是mtcars 两倍的矩阵,所以当我尝试时:

brokenModel <- lm(mpg ~ bs(cbind(mtcars$wt,mtcars$hp), knots = 4), data = mtcars)

我收到关于不同长度的错误。

我的问题是:如果模型的行数与结果变量不同,如何在模型中使用多元样条?我是否将结果变量堆叠在自身之上 y &lt;- c(y, y)?为什么多元样条会产生额外的行?

谢谢。

【问题讨论】:

标签: r regression lm spline


【解决方案1】:

在这种情况下,您不能使用splines::bs,因为它严格用于构造单变量样条。如果你做bs(mat) 其中mat 是一个矩阵,它只是在做bs(c(mat))。例如,

mat <- matrix(runif(8), 4, 2)
identical(bs(mat), bs(c(mat)))
# [1] TRUE

这解释了为什么在执行bs(cbind(mtcars$wt,mtcars$hp) 时会得到双倍的行数。


要创建 2D 样条,最简单的方法是创建附加样条:

lm(mpg ~ factor(gear) + bs(wt, knots = 5) + bs(hp, knots = 4), mtcars)

但这可能不是您想要的。然后考虑交互:

model <- lm(mpg ~ factor(gear) + bs(wt, knots = 5):bs(hp, knots = 4), mtcars)

bs(wt, knots = 5):bs(hp, knots = 4) 在两个设计矩阵之间形成逐行克罗内克积。由于bs(wt, knots = 5) 是4 列矩阵,bs(hp, knots = 4) 是3 列矩阵,因此交互有4 * 3 = 12 列。


或者,考虑使用mgcv 包。在mgcv 中,可以通过两种方式构造多变量样条:

  • 各向同性薄板样条;
  • 尺度不变张量积样条。

显然你想要第二个,因为wt 和hp 有不同的单位。要构造张量积样条,我们可以使用:

library(mgcv)
fit <- gam(mpg ~ factor(gear)
                 + s(wt, bs = 'cr', k = 4, fx = TRUE)
                 + s(hp, bs = 'cr', k = 4, fx = TRUE)
                 + ti(wt, hp, bs = 'cr', k = c(4, 4), d = c(1, 1), fx = TRUE),
                 data = mtcars)

这里我特意设置fx = TRUE禁用惩罚回归。

我不想写一个广泛的答案来介绍mgcv。关于s、ti 和gam 的工作原理,请阅读文档。如果您需要弥补理论上的差距,请阅读 Simon Wood 于 2006 年出版的书:Generalized Additive Models: an Introduction with R。


mgcv 用法的实际示例?

我有一个答案Cubic spline method for longitudinal series data,它可能会帮助您熟悉mgcv。但作为一个介绍性示例,它仅展示了如何使用单变量样条。幸运的是,这也是关键。张量积样条是由单变量样条构造的。

我与mgcv 相关的其他答案更多是理论方面的;虽然并非所有与spline 相关的答案都参考了mgcv。所以这个问题和答案是我在这个阶段能给你的最好的。

尺度不变张量积样条是否等同于径向平滑,还是各向同性的薄位样条?

径向平滑等效于薄板样条,因为薄板样条的基函数是径向的。这就是为什么它是各向同性的,可以用于空间回归。

张量积样条是尺度不变的,因为它被构造为单变量样条基的(成对)乘法。

【讨论】:

    猜你喜欢
    • 2019-07-02
    • 1970-01-01
    • 2016-11-07
    • 2013-07-11
    • 2014-06-26
    • 2017-07-30
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多