为什么 `ns` 和 `rcs` 在 R 中生成不同的预测?
Why does `ns` and `rcs` generate different predictions in R?
我的理解是 rcs()
(来自 rms
包)使用截断幂基来表示自然(受限)三次样条。或者,我可以使用使用 B 样条基础的 ns()
(来自 splines
包)。
但是,我注意到训练拟合和测试预测可能非常不同(尤其是在外推 x
时)。我试图了解 rcs()
和 ns()
之间的区别,以及我是否可以互换使用这些函数。
伪造非线性数据。
library(tidyverse)
library(splines)
library(rms)
set.seed(100)
xx <- rnorm(1000)
yy <- 10 + 5*xx - 0.5*xx^2 - 2*xx^3 + rnorm(1000, 0, 4)
df <- data.frame(x=xx, y=yy)
用 ns
和另一个 rcs
用相同的结拟合一个模型。
ns_mod <- lm(y ~ ns(x, knots=c(-2, 0, 2)), data=df)
ddist <- datadist(df)
options("datadist" = "ddist")
trunc_power_mod <- ols(y ~ rcs(x, knots=c(-2, 0, 2)), data=df)
检查它们的拟合度 (MSE)。
mean(ns_mod$residuals^2)
mean(trunc_power_mod$residuals^2)
df$pred_ns <- ns_mod$fitted.values
df$pred_trunc_power <- trunc_power_mod$fitted.values
df_melt <- df %>%
gather(key="model", value="predictions", -x, -y)
ggplot(df_melt, aes(x=x, y=y)) +
geom_point(alpha=0.1) +
geom_line(aes(x=x, y=predictions, group=model, linetype=model))
生成测试数据集并绘制两个模型之间的预测。
newdata <- data.frame(x=seq(-10, 10, 0.1))
pred_ns_new <- predict(ns_mod, newdata=newdata)
pred_trunc_new <- predict(trunc_power_mod, newdata=newdata)
newdata$pred_ns_new <- pred_ns_new
newdata$pred_trunc_new <- pred_trunc_new
newdata_melted <- newdata %>%
gather(key="model", value="predictions", -x)
ggplot(newdata_melted, aes(x=x, y=predictions, group=model, linetype=model)) +
geom_line()
有一个相当简单的解释:knots
不是 rcs()
的参数。它希望使用参数 parms
指定节点。另一个问题是 ns()
的 knots
参数没有指定“边界结”,默认为 range(x)
。所以要得到相同的预测,你需要
trunc_power_mod <- ols(y ~ rcs(x, parms=c(min(x), -2, 0, 2, max(x))), data=df)
我的理解是 rcs()
(来自 rms
包)使用截断幂基来表示自然(受限)三次样条。或者,我可以使用使用 B 样条基础的 ns()
(来自 splines
包)。
但是,我注意到训练拟合和测试预测可能非常不同(尤其是在外推 x
时)。我试图了解 rcs()
和 ns()
之间的区别,以及我是否可以互换使用这些函数。
伪造非线性数据。
library(tidyverse)
library(splines)
library(rms)
set.seed(100)
xx <- rnorm(1000)
yy <- 10 + 5*xx - 0.5*xx^2 - 2*xx^3 + rnorm(1000, 0, 4)
df <- data.frame(x=xx, y=yy)
用 ns
和另一个 rcs
用相同的结拟合一个模型。
ns_mod <- lm(y ~ ns(x, knots=c(-2, 0, 2)), data=df)
ddist <- datadist(df)
options("datadist" = "ddist")
trunc_power_mod <- ols(y ~ rcs(x, knots=c(-2, 0, 2)), data=df)
检查它们的拟合度 (MSE)。
mean(ns_mod$residuals^2)
mean(trunc_power_mod$residuals^2)
df$pred_ns <- ns_mod$fitted.values
df$pred_trunc_power <- trunc_power_mod$fitted.values
df_melt <- df %>%
gather(key="model", value="predictions", -x, -y)
ggplot(df_melt, aes(x=x, y=y)) +
geom_point(alpha=0.1) +
geom_line(aes(x=x, y=predictions, group=model, linetype=model))
生成测试数据集并绘制两个模型之间的预测。
newdata <- data.frame(x=seq(-10, 10, 0.1))
pred_ns_new <- predict(ns_mod, newdata=newdata)
pred_trunc_new <- predict(trunc_power_mod, newdata=newdata)
newdata$pred_ns_new <- pred_ns_new
newdata$pred_trunc_new <- pred_trunc_new
newdata_melted <- newdata %>%
gather(key="model", value="predictions", -x)
ggplot(newdata_melted, aes(x=x, y=predictions, group=model, linetype=model)) +
geom_line()
有一个相当简单的解释:knots
不是 rcs()
的参数。它希望使用参数 parms
指定节点。另一个问题是 ns()
的 knots
参数没有指定“边界结”,默认为 range(x)
。所以要得到相同的预测,你需要
trunc_power_mod <- ols(y ~ rcs(x, parms=c(min(x), -2, 0, 2, max(x))), data=df)