平滑样条
在第二章中，我们学习了如何利用节点和基函数来拟合回归样条。分析发现，节点数量过多将导致较低的 MSE ，这意味着过度拟合了曲线特征。
以使用 25 个节点的自然样条曲线为例：
显然，这条曲线过度拟合了数据。因此，我们将使用类似正则化的方法来解决这个问题。我们可以选择许多节点，但是需要限制（惩罚）曲线的粗糙度。
注意，一阶导数表示曲线的斜率，二阶导数表示斜率的变化速度。因此，曲线的二阶导数与曲线的粗糙度有关，可以将其作为残差平方和中的惩罚项：
这就是平滑样条。在平滑样条中，惩罚项将处理曲线的粗糙度，因此节点的数量不是那么重要。
下面的动图，显示了不同 构建的平滑样条对曲线的拟合效果，每个模型的节点个数都是由smooth.spline()函数自动选择的， 并给出了对应模型的均方误差。
越大，曲线越光滑； 越小，曲线越粗糙。利用交叉验证可以选择最优 。此外，平滑样条可以被纳入广义线性模型框架，形成的模型通常被称为广义加性模型 （generalised additive model, GAM） 。我们可以用平滑样条曲线来模拟自变量与因变量之间的非线性关系:
其中，可以是一个线性函数。如果所有的 都是线性函数，那么该模型就是一个广义线性模型（generalised linear model, GLM）。
练习题
Q1：拟合平滑样条
我们将继续使用MultiKink 包中的triceps数据集，使用平滑三次样条模拟三头肌皮褶厚度（triceps）与年龄（age）的关系。
下面使用smooth.spline()函数拟合平滑三次样条：
library(splines) library(MultiKink) #for the data library(ggplot2) #for the plots set.seed(1974) #fix the random generator seed data("triceps") #load the dataset triceps #notice that the variable of interest #it is also called triceps. Don't get #confused! #smooth spline with automatic number of knots chosen #and penalisation chosen by leave-one-out CV (this is the #option cv=T, otherwise generalized’ cross-validation is used) sspline <- smooth.spline(triceps$age, triceps$triceps, cv=TRUE)
可视化拟合情况：
plot(triceps$age, triceps$triceps) lines(sspline, col="blue")
【注意】:smooth.spline()函数自带的交叉验证默认选择了一个非常低的 值，这导致欠拟合。可以通过自定义 值来避免这种情况。例如：
sspline <- smooth.spline(triceps$age, triceps$triceps, lambda=.0001) plot(triceps$age, triceps$triceps) lines(sspline, col="blue")
也可以改变节点的数量来调整平滑度。从下图可以看到，节点数量越多曲线越粗糙：
sspline19 <- smooth.spline(triceps$age, triceps$triceps, df=19) sspline30 <- smooth.spline(triceps$age, triceps$triceps, df=30) plot(triceps$age, triceps$triceps) lines(sspline19, col="red") lines(sspline30, col="green")
下面预测女性在 10 岁和 30 岁时三头肌皮褶厚度：
predict(sspline, x=c(10,30)) ## $x ## [1] 10 30 ## ## $y ## [1] 6.470573 14.290032
Q2: 拟合一个可加模型
bmd.csv 数据集包含 169 个不同年龄男性和女性的骨密度（bmd）测量值。我们想用年龄（age）、性别（sex）和 bmi 作为预测指标来拟合骨密度模型。
#read the data and compute BMI bmd.data <- read.csv("bmd.csv", stringsAsFactors = TRUE) bmd.data$bmi <- bmd.data$weight_kg / (bmd.data$height_cm/100)^2
Q3: 拟合自然三次样条
我们将继续使用triceps数据集来拟合一个自然三次样条，来模拟女性的三头肌皮褶厚度（triceps）与年龄（age）的关系。
我们将使用mgcv包中的gam()函数。注意，gam包中有一个同名gam()函数。两个函数的区别在于：mgcv中gam()函数可以使用广义交叉验证来自动选择参数。
模型中的s()函数表明为相应的变量创建平滑样条，我们将使用与三次回归样条相同的基样条：bs="cr"。
#libraries that we will need library(mgcv) #package for gam set.seed(1974) #fix the random generator seed bmd.gam <- gam(bmd ~ s(age, bs="cr")+ s(bmi, bs="cr") + sex, data=bmd.data) summary(bmd.gam) #libraries that we will need library(mgcv) #package for gam set.seed(1974) #fix the random generator seed bmd.gam <- gam(bmd ~ s(age, bs="cr")+ s(bmi, bs="cr") + sex, data=bmd.data) summary(bmd.gam) ## ## Family: gaussian ## Link function: identity ## ## Formula: ## bmd ~ s(age, bs = "cr") + s(bmi, bs = "cr") + sex ## ## Parametric coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 0.74193 0.01456 50.97 < 2e-16 *** ## sexM 0.08092 0.02064 3.92 0.000131 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Approximate significance of smooth terms: ## edf Ref.df F p-value ## s(age) 1.035 1.070 17.65 3.02e-05 *** ## s(bmi) 5.687 6.611 10.09 < 2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## R-sq.(adj) = 0.381 Deviance explained = 40.9% ## GCV = 0.018099 Scale est. = 0.017165 n = 169
age变量的有效自由度数约为 1，这表明年龄对骨密度的影响是线性的。我们可以为每个预测变量绘制相应的样条曲线：
plot(bmd.gam)
我们还可以绘制出每个性别对应的拟合曲面。以女性为例：
# Let's create a grid to be used in persp() steps <- 60 age <- with(bmd.data, seq(min(age), max(age), length = steps) ) bmi <- with(bmd.data, seq(min(bmi), max(bmi), length = steps) ) newdat <- expand.grid(age = age, #grid bmi = bmi, sex="F") bmdpred <- matrix(predict(bmd.gam, newdat), steps, steps) #predictions # Now plot it p <- persp(age, bmi, bmdpred, theta = 65, #angle of the perspective col = "green") # To add the points, you need the same 3d transformation obs <- with(bmd.data[bmd.data$sex=="F", ], trans3d(age, bmi, bmd, p)) pred <- with(bmd.data[bmd.data$sex=="F", ], trans3d(age, bmi, fitted(bmd.gam)[bmd.data$sex=="F"], p)) # Add segments to show the points and where they are in 3d points(obs, col = "red", pch = 16) segments(obs$x, obs$y, pred$x, pred$y)
参考资料
[1]
原文地址: https://bookdown.org/tpinto_home/Beyond-Linearity/polynomial-regression.html
