R smooth.spline(): 平滑样条曲线不平滑但过度拟合我的数据
R smooth.spline(): smoothing spline is not smooth but overfitting my data
我有几个数据点似乎适合通过它们拟合样条曲线。当我这样做的时候，我得到了一个相当颠簸的拟合，就像过度拟合，这不是我理解的平滑。
是否有一个特殊的选项/参数来恢复像 here 这样非常平滑的样条函数的功能。
对 smooth.spline 使用 penalty 参数没有任何可见的效果。也许我做错了？
这里是数据和代码：
results <- structure(
list(
beta = c(
0.983790622281964, 0.645152464354322,
0.924104713597375, 0.657703886566088, 0.788138034115623, 0.801080207252363,
1, 0.858337365965949, 0.999687052533693, 0.666552625121279, 0.717453633245958,
0.621570152961453, 0.964658181346544, 0.65071758770312, 0.788971505000918,
0.980476054183113, 0.670263506919246, 0.600387040967624, 0.759173403408052,
1, 0.986409675965, 0.982996471134736, 1, 0.995340781899163, 0.999855895958986,
1, 0.846179233381267, 0.879226324448832, 0.795820998892035, 0.997586607285667,
0.848036806290156, 0.905320944437968, 0.947709125535428, 0.592172373022407,
0.826847031044922, 0.996916006944244, 0.785967729206612, 0.650346929853076,
0.84206351833549, 0.999043126652724, 0.936879214753098, 0.76674066557003,
0.591431233516217, 1, 0.999833445117791, 0.999606223666537, 0.6224971799303,
1, 0.974537160571494, 0.966717133936379
), inventoryCost = c(
1750702.95138889,
442784.114583333, 1114717.44791667, 472669.357638889, 716895.920138889,
735396.180555556, 3837320.74652778, 872873.4375, 2872414.93055556,
481095.138888889, 538125.520833333, 392199.045138889, 1469500.95486111,
459873.784722222, 656220.486111111, 1654143.83680556, 437511.458333333,
393295.659722222, 630952.170138889, 4920958.85416667, 1723517.10069444,
1633579.86111111, 4639909.89583333, 2167748.35069444, 3062420.65972222,
5132702.34375, 838441.145833333, 937659.288194444, 697767.1875,
2523016.31944444, 800903.819444444, 1054991.49305556, 1266970.92013889,
369537.673611111, 764995.399305556, 2322879.6875, 656021.701388889,
458403.038194444, 844133.420138889, 2430700, 1232256.68402778,
695574.479166667, 351348.524305556, 3827440.71180556, 3687610.41666667,
2950652.51736111, 404550.78125, 4749901.64930556, 1510481.59722222,
1422708.07291667
)
), .Names = c("beta", "inventoryCost"), class = c("data.frame")
)
plot(results$beta,results$inventoryCost)
mySpline <- smooth.spline(results$beta,results$inventoryCost, penalty=999999)
lines(mySpline$x, mySpline$y, col="red", lwd = 2)
我认为你不应该使用/想要 splinefun。我建议改用 GAM：
library(mgcv)
fit <- gam(inventoryCost ~ s(beta, bs = "cr", k = 20), data = results)
summary(fit)
gam.check(fit)
plot(fit)
plot(inventoryCost ~ beta, data = results, col = "dark red", , pch = 16)
curve(predict(fit, newdata = data.frame(beta = x)), add = TRUE,
from = min(results$beta), to = max(results$beta), n = 1e3, lwd = 2)
建模前合理转换数据
根据您results$inventoryCost的规模，对数变换是合适的。为了简单起见，下面我使用x、y。我也在重新排序您的数据，以便 x 升序：
x <- results$beta; y <- log(results$inventoryCost)
reorder <- order(x); x <- x[reorder]; y <- y[reorder]
par(mfrow = c(1,2))
plot(x, y, main = "take log transform")
hist(x, main = "x is skewed")
左图更好看？另外，强烈建议对 x 进一步进行变换，因为它是倾斜的！ （见右图）
以下转换是合适的：
x1 <- -(1-x)^(1/3)
(1-x) 的立方根将使数据在 x = 1 周围更加分散。我添加了一个额外的 -1 以便在 x 和 x1 之间存在正单调关系而不是负单调关系。现在让我们检查一下关系：
par(mfrow = c(1,2))
plot(x1, y, main = expression(y %~% ~ x1))
hist(x1, main = "x1 is well spread out")
拟合样条
现在我们已准备好进行统计建模。尝试以下调用：
fit <- smooth.spline(x1, y, nknots = 10)
pred <- stats:::predict.smooth.spline(fit, x1)$y ## predict at all x1
## or you can simply call: pred <- predict(fit, x1)$y
plot(x1, y) ## scatter plot
lines(x1, pred, lwd = 2, col = 2) ## fitted spline
好看吗？请注意，我已经使用 nknots = 10 告诉 smooth.spline 放置 10 个 interior 节（按分位数）；因此，我们要拟合 惩罚回归样条 而不是平滑样条。事实上，smooth.spline() 函数几乎从不适合平滑样条，除非你输入 all.knots = TRUE（见后面的例子）。
我也放弃了penalty = 999999，因为这与平滑度控制无关。如果你真的想控制平滑度，而不是让 smooth.spline 通过 GCV 找出最优的，你应该使用参数 df 或 spar。后面会举例子
要将适合度转换回原始比例，请执行以下操作：
plot(x, exp(y), main = expression(Inventory %~%~ beta))
lines(x, exp(pred), lwd = 2, col = 2)
如您所见，拟合样条曲线与您预期的一样平滑。
拟合样条的解释
让我们看看您的拟合样条曲线的摘要：
> fit
Smoothing Parameter spar= 0.4549062 lambda= 0.0008657722 (11 iterations)
Equivalent Degrees of Freedom (Df): 6.022959
Penalized Criterion: 0.08517417
GCV: 0.004288539
我们使用了 10 节，最终有 6 个自由度，所以惩罚抑制了大约 4 个参数。 GCV 选择的平滑参数，经过 11 次迭代后，为 lambda= 0.0008657722.
为什么要把x改成x1
样条曲线受到二阶导数的惩罚，但这种惩罚是在所有数据点的 averaged/integrated 二阶导数上。现在，查看您的数据 (x, y)。对于0.98之前的x，关系比较稳定；当 x 接近 1 时，关系会迅速变陡。 "change point",0.98，二阶导数非常高，远高于其他位置的二阶导数。
y0 <- as.numeric(tapply(y, x, mean)) ## remove tied values
x0 <- unique(x) ## remove tied values
dy0 <- diff(y0)/diff(x0) ## 1st order difference
ddy0 <- diff(dy0)/diff(x0[-1]) ## 2nd order difference
plot(x0[1:43], abs(ddy0), pch = 19)
看看那个二阶的巨大尖峰 difference/derivative！现在，如果我们直接拟合样条曲线，围绕这个变化点的样条曲线将受到严重惩罚.
bad <- smooth.spline(x, y, all.knots = TRUE)
bad.pred <- predict(bad, x)$y
plot(x, exp(y), main = expression(Inventory %~% ~ beta))
lines(x, exp(bad.pred), col = 2, lwd = 3)
abline(v = 0.98, lwd = 2, lty = 2)
你可以清楚地看到样条曲线在x = 0.98之后逼近数据有一些困难。
当然有一些方法可以在这个变化点之后实现更好的逼近，例如，通过手动设置更小的平滑参数，或者更高的自由度。但我们正在走向另一个极端。请记住，惩罚和自由度都是 全局度量 。增加模型的复杂度会在x = 0.98之后得到更好的逼近，但也会让其他部分更加颠簸。现在让我们尝试一个自由度为 45 的模型：
worse <- smooth.spline(x, y, all.knots = TRUE, df = 45)
worse.pred <- predict(worse, x)$y
plot(x, exp(y), main = expression(Inventory %~% ~ beta))
lines(x, exp(worse.pred), col = 2, lwd = 2)
如您所见，曲线是颠簸的。当然，我们已经过度拟合了 50 个数据的数据集，具有 45 个自由度。
其实你原来误用smooth.spline()也是在做同样的事情：
> mySpline
Call:
smooth.spline(x = results$beta, y = results$inventoryCost, penalty = 999999)
Smoothing Parameter spar= -0.8074624 lambda= 3.266077e-19 (17 iterations)
Equivalent Degrees of Freedom (Df): 45
Penalized Criterion: 5.598386
GCV: 0.03824885
糟糕，45 自由度，过拟合！
我有几个数据点似乎适合通过它们拟合样条曲线。当我这样做的时候，我得到了一个相当颠簸的拟合，就像过度拟合，这不是我理解的平滑。
是否有一个特殊的选项/参数来恢复像 here 这样非常平滑的样条函数的功能。
对 smooth.spline 使用 penalty 参数没有任何可见的效果。也许我做错了？
这里是数据和代码：
results <- structure(
list(
beta = c(
0.983790622281964, 0.645152464354322,
0.924104713597375, 0.657703886566088, 0.788138034115623, 0.801080207252363,
1, 0.858337365965949, 0.999687052533693, 0.666552625121279, 0.717453633245958,
0.621570152961453, 0.964658181346544, 0.65071758770312, 0.788971505000918,
0.980476054183113, 0.670263506919246, 0.600387040967624, 0.759173403408052,
1, 0.986409675965, 0.982996471134736, 1, 0.995340781899163, 0.999855895958986,
1, 0.846179233381267, 0.879226324448832, 0.795820998892035, 0.997586607285667,
0.848036806290156, 0.905320944437968, 0.947709125535428, 0.592172373022407,
0.826847031044922, 0.996916006944244, 0.785967729206612, 0.650346929853076,
0.84206351833549, 0.999043126652724, 0.936879214753098, 0.76674066557003,
0.591431233516217, 1, 0.999833445117791, 0.999606223666537, 0.6224971799303,
1, 0.974537160571494, 0.966717133936379
), inventoryCost = c(
1750702.95138889,
442784.114583333, 1114717.44791667, 472669.357638889, 716895.920138889,
735396.180555556, 3837320.74652778, 872873.4375, 2872414.93055556,
481095.138888889, 538125.520833333, 392199.045138889, 1469500.95486111,
459873.784722222, 656220.486111111, 1654143.83680556, 437511.458333333,
393295.659722222, 630952.170138889, 4920958.85416667, 1723517.10069444,
1633579.86111111, 4639909.89583333, 2167748.35069444, 3062420.65972222,
5132702.34375, 838441.145833333, 937659.288194444, 697767.1875,
2523016.31944444, 800903.819444444, 1054991.49305556, 1266970.92013889,
369537.673611111, 764995.399305556, 2322879.6875, 656021.701388889,
458403.038194444, 844133.420138889, 2430700, 1232256.68402778,
695574.479166667, 351348.524305556, 3827440.71180556, 3687610.41666667,
2950652.51736111, 404550.78125, 4749901.64930556, 1510481.59722222,
1422708.07291667
)
), .Names = c("beta", "inventoryCost"), class = c("data.frame")
)
plot(results$beta,results$inventoryCost)
mySpline <- smooth.spline(results$beta,results$inventoryCost, penalty=999999)
lines(mySpline$x, mySpline$y, col="red", lwd = 2)
我认为你不应该使用/想要 splinefun。我建议改用 GAM：
library(mgcv)
fit <- gam(inventoryCost ~ s(beta, bs = "cr", k = 20), data = results)
summary(fit)
gam.check(fit)
plot(fit)
plot(inventoryCost ~ beta, data = results, col = "dark red", , pch = 16)
curve(predict(fit, newdata = data.frame(beta = x)), add = TRUE,
from = min(results$beta), to = max(results$beta), n = 1e3, lwd = 2)
建模前合理转换数据
根据您results$inventoryCost的规模，对数变换是合适的。为了简单起见，下面我使用x、y。我也在重新排序您的数据，以便 x 升序：
x <- results$beta; y <- log(results$inventoryCost)
reorder <- order(x); x <- x[reorder]; y <- y[reorder]
par(mfrow = c(1,2))
plot(x, y, main = "take log transform")
hist(x, main = "x is skewed")
左图更好看？另外，强烈建议对 x 进一步进行变换，因为它是倾斜的！ （见右图）
以下转换是合适的：
x1 <- -(1-x)^(1/3)
(1-x) 的立方根将使数据在 x = 1 周围更加分散。我添加了一个额外的 -1 以便在 x 和 x1 之间存在正单调关系而不是负单调关系。现在让我们检查一下关系：
par(mfrow = c(1,2))
plot(x1, y, main = expression(y %~% ~ x1))
hist(x1, main = "x1 is well spread out")
拟合样条
现在我们已准备好进行统计建模。尝试以下调用：
fit <- smooth.spline(x1, y, nknots = 10)
pred <- stats:::predict.smooth.spline(fit, x1)$y ## predict at all x1
## or you can simply call: pred <- predict(fit, x1)$y
plot(x1, y) ## scatter plot
lines(x1, pred, lwd = 2, col = 2) ## fitted spline
好看吗？请注意，我已经使用 nknots = 10 告诉 smooth.spline 放置 10 个 interior 节（按分位数）；因此，我们要拟合 惩罚回归样条 而不是平滑样条。事实上，smooth.spline() 函数几乎从不适合平滑样条，除非你输入 all.knots = TRUE（见后面的例子）。
我也放弃了penalty = 999999，因为这与平滑度控制无关。如果你真的想控制平滑度，而不是让 smooth.spline 通过 GCV 找出最优的，你应该使用参数 df 或 spar。后面会举例子
要将适合度转换回原始比例，请执行以下操作：
plot(x, exp(y), main = expression(Inventory %~%~ beta))
lines(x, exp(pred), lwd = 2, col = 2)
如您所见，拟合样条曲线与您预期的一样平滑。
拟合样条的解释
让我们看看您的拟合样条曲线的摘要：
> fit
Smoothing Parameter spar= 0.4549062 lambda= 0.0008657722 (11 iterations)
Equivalent Degrees of Freedom (Df): 6.022959
Penalized Criterion: 0.08517417
GCV: 0.004288539
我们使用了 10 节，最终有 6 个自由度，所以惩罚抑制了大约 4 个参数。 GCV 选择的平滑参数，经过 11 次迭代后，为 lambda= 0.0008657722.
为什么要把x改成x1
样条曲线受到二阶导数的惩罚，但这种惩罚是在所有数据点的 averaged/integrated 二阶导数上。现在，查看您的数据 (x, y)。对于0.98之前的x，关系比较稳定；当 x 接近 1 时，关系会迅速变陡。 "change point",0.98，二阶导数非常高，远高于其他位置的二阶导数。
y0 <- as.numeric(tapply(y, x, mean)) ## remove tied values
x0 <- unique(x) ## remove tied values
dy0 <- diff(y0)/diff(x0) ## 1st order difference
ddy0 <- diff(dy0)/diff(x0[-1]) ## 2nd order difference
plot(x0[1:43], abs(ddy0), pch = 19)
看看那个二阶的巨大尖峰 difference/derivative！现在，如果我们直接拟合样条曲线，围绕这个变化点的样条曲线将受到严重惩罚.
bad <- smooth.spline(x, y, all.knots = TRUE)
bad.pred <- predict(bad, x)$y
plot(x, exp(y), main = expression(Inventory %~% ~ beta))
lines(x, exp(bad.pred), col = 2, lwd = 3)
abline(v = 0.98, lwd = 2, lty = 2)
你可以清楚地看到样条曲线在x = 0.98之后逼近数据有一些困难。
当然有一些方法可以在这个变化点之后实现更好的逼近，例如，通过手动设置更小的平滑参数，或者更高的自由度。但我们正在走向另一个极端。请记住，惩罚和自由度都是 全局度量 。增加模型的复杂度会在x = 0.98之后得到更好的逼近，但也会让其他部分更加颠簸。现在让我们尝试一个自由度为 45 的模型：
worse <- smooth.spline(x, y, all.knots = TRUE, df = 45)
worse.pred <- predict(worse, x)$y
plot(x, exp(y), main = expression(Inventory %~% ~ beta))
lines(x, exp(worse.pred), col = 2, lwd = 2)
如您所见，曲线是颠簸的。当然，我们已经过度拟合了 50 个数据的数据集，具有 45 个自由度。
其实你原来误用smooth.spline()也是在做同样的事情：
> mySpline
Call:
smooth.spline(x = results$beta, y = results$inventoryCost, penalty = 999999)
Smoothing Parameter spar= -0.8074624 lambda= 3.266077e-19 (17 iterations)
Equivalent Degrees of Freedom (Df): 45
Penalized Criterion: 5.598386
GCV: 0.03824885
糟糕，45 自由度，过拟合！