非线性回归模型与平滑回归
10 Dec 2012
开始的话
回归模型regression model是我接触的第一个统计模型。实际上，未学统计之前，高中的课本上就介绍了最小二乘法。本科学计量时，又学了个底朝天。什么一元的，多元的，结构的，滞后的各种假设，各种情形，各种解法。当时就觉得有点没劲，有本事你别整这么多假设，整个放之四海皆为准的出来算你厉害。有种很不屑很愤青的感觉，所以也嗤之以鼻，不怎么学，觉得是一些江湖骗术，那些人的饭碗罢了。后来看到Cox说：All the models are wrong,but some are useful.明白这句话，才知道统计学发展的实质。其实所谓模型没有一个是对的，只是对现实的一种近似，只有上帝才知道事实是长什么摸样的。
这两天看非线性回归和平滑回归。才知道除了经典的计量回归，还有另外一片天。你自己想像符合什么模型，你就设为什么模型，随心所欲，别太离谱就行。另外，平滑回归更是五花八门，诸如样条函数，局部回归，加性模型，甚至还有时髦的神经网络。
对于简单的线性回归模型，主要做两个方面的扩展。其一，由线性到非线性，主要是均值函数形式的变化，包括nls和神经网络。nls比较具体，可以设置固定的非线性函数形式。神经网络比较宽泛，可以近似任意一个连续函数，是一种非常广义的非线性回归。其二，由线性到平滑，包括一维的平滑以及多维的平滑（可加模型），平滑得到结果与参数估计的感觉不一样，它主要是一种探索性工具，用来对数据进行预分析，对数据的关系可以有一个大致的了解。
神经网络模型
神经网络可以看做是一种非线性回归模型。它可以近似任意一个函数。它就相当于一个机制，这个机制在接收数据输入和输出后，就可以进行训练，改变其自身的结构，尤其是其中的一些权重，从而使得它适用于所要考虑的问题。所以，称它为一种自适应系统。一句话，它可以对输入数据和输出数据之间的关系进行建模。
举一个模拟的例子。
p <- 2
N <- 2000
set.seed(1)
x <- matrix(rnorm(p*N, mean = 0, sd = 3), ncol = 2) 以上数据作为输入。
plot(x[, 1], x[, 2], pch = 19, col = rgb(0, 0, 0, 0.1))
输入与输出的关系模式这样定义：
fun.y <- function(x) {
if (x[1] ^2 + x[2]^2 < 9) {
return(1)
} else if (x[1] ^2 + x[2]^2 < 36) {
return(2)
} else if (x[1] ^2 + x[2]^2 >36) {
return(3)
}
}
y <- apply(x, 1, fun.y)
输入与输出关系展示在图中是这样的：
plot(x[, 1], x[, 2], col = c(1, 2, 3)[y], asp = 1, pch = 19)
我们检验利用神经网络算法是否可以发现这样的模式。
mydata <- data.frame(x = x/3, y = y/3)
nn1 <- nnet(y ~ x.1 + x.2, data = mydata, entropy = T,
size =3, maxit = 100000, decay = 0.0000000001, trace = T, skip = T)
nn1记为模型的结果。我们用该模型来对原始数据进行预测，并观察预测结果与真实值之间的差异。
y.predict <- as.numeric(predict(nn1, type = "raw"))
y.predict.hat <- NULL
for (i in 1:2000) {
if (y.predict[i] * 3 < 1.5) {
y.predict.hat <- c(y.predict.hat, 1)
} else if (y.predict[i] * 3 < 2.5) {
y.predict.hat <- c(y.predict.hat, 2)
} else {
y.predict.hat <- c(y.predict.hat, 3)
}
}
oldpar <- par(mfrow = c(1, 2))
plot(x[, 1], x[, 2], col = c(1, 2, 3, 4)[y], asp = 1, pch = 19)
plot(x[, 1], x[, 2], col = c(1, 2, 3, 4)[y.predict.hat],pch = 19, asp = 1)
par(oldpar)
左边是真实数据，右边是预测数据。可见，预测结果与真实结果十分接近。数量上的分析结果如下：
> table(y, y.predict.hat)
y.predict.hat
y 1 2 3
1 735 0 0
2 80 797 65
3 0 2 321
准确率超过90%。
前边提到，神经网络其实是对任意非线性函数的近似。刚才举例，其实是三维空间中的一个阶梯函数，它也属于非线性函数。我们看下，神经网络模型得到的结果是不是与我们最初设计的阶梯函数相似，我们的工具是三维空间中的立体图。
x1 <- seq(-5, 5, length = 50)
x2 <- x1
x <- expand.grid(x1, x2)
colnames(x) <- c("x.1", "x.2")
y.predict <- as.numeric(predict(nn1,x ,type = "raw"))
y.predict.hat <- NULL
for (i in 1:2500) {
if (y.predict[i] * 3 < 1.5) {
y.predict.hat <- c(y.predict.hat, 1)
} else if (y.predict[i] * 3 < 2.5) {
y.predict.hat <- c(y.predict.hat, 2)
} else {
y.predict.hat <- c(y.predict.hat, 3)
}
}
data.grid <- cbind(x, y.predict.hat)
library(lattice)
trellis.device()
wireframe(y.predict.hat ~ x.1 + x.2, data.grid, drape = T,
screen = list(z = 20, x = - 50))
#wireframe(y.predict.hat ~ x.1 + x.2, data.grid, drape = T,
# screen = list(z = 20, x = - 160))
|
|
从函数（经过处理之后的）的立体图像可见，神经网络的结果非常类似于一个阶梯函数。没有经过处理的结果为：
GAM (Generalized Additive Model) 广义可加模型
广义可加模型是多元线性模型的扩展，多元线性模型中，每个相加项都是线性形式，而可加模型中，每个可加项是一个平滑项。该平滑项可以逼近常见的各种函数形式。顾名思义，可加模型处理的是具有可加性质的选择变量。下面举例说明，可加模型的确可以准确的捕获数据本身所具有的一些内在联系。
首先，加载需要的包mgcv:
library(mgcv)
然后生成模拟数据，注意因变量与自变量之间所具有的关系，一个是变量不再是线性变换，而是给出了多样的函数形式；另一方面就是可加性，没有f(x)*f(y)这种乘法形式。把模拟的数据组建为数据框。
x1 <- rnorm(150)
x2 <- rnorm(150)
x3 <- rnorm(150)
y <- x1^2 + x2^3 + 10 * sin(x3) + rnorm(150, 0, 2)
data <- data.frame(x1, x2, x3, y)
然后，就可以用gam来挖掘y与其他几个变量的关系。
m <- gam(y ~ s(x1) + s(x2) + s(x3) - 1, data = data)
提取结果。
> summary(m)
Family: gaussian
Link function: identity
Formula:
y ~ s(x1) + s(x2) + s(x3) - 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(x1) 2.789 3.502 7.90 2.46e-05 ***
s(x2) 6.463 7.612 48.19 < 2e-16 ***
s(x3) 4.992 6.114 161.51 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.931 Deviance explained = 91.8%
GCV score = 6.7232 Scale est. = 6.0848 n = 150
三个平滑可加项s(x1),s(x2),s(x3)都特别显著，另外模型整体解释能力达到Deviance explained = 91.8%.
通过图形观察三个平滑项的的处理结果与开始的设计是否一致。
首先看x1^2与s(x1)的比较。
x <- seq(-3, 3, length = 200)
plot(m, se=F, select = 1, lwd = 2, lty = 1, col = "blue")
lines(x, x^2, col = 2, lwd = 2, lty = 2)
legend(0.2, 15, legend = c("gam处理结果", "数据期初设计"), lty = c(1, 2), lwd = c(2, 2),
bty = "n", col = c("blue", "red"))
gam的处理结果与数据最初的设计比较接近。二次函数的趋势表现的也比较直观，在此基础上，可以设置参数形式进行参数估计。其他两个平滑项与数据开始设计的关系也可以得到：
plot(m, se=F, select = 2, lwd = 2, lty = 1, col = "blue")
lines(x, x^3, col = 2, lwd = 2, lty = 2)
legend(-3, 15, legend = c("gam处理结果", "数据期初设计"), lty = c(1, 2), lwd = c(2, 2),
bty = "n", col = c("blue", "red"))
plot(m, se=F, select = 3, lwd = 2, lty = 1, col = "blue")
lines(x, 10*sin(x), col = 2, lwd = 2, lty = 2)
legend(-3, 15, legend = c("gam处理结果", "数据期初设计"), lty = c(1, 2), lwd = c(2, 2),
bty = "n", col = c("blue", "red"))
|
|
总体来看，gam所使用的平滑方法虽然可以很好地拟合数据，但是，它的解释能力并不尽人意，我们无法通过这个模型来解释一些信息。尽管如此，gam仍可作为一种探索性工具，我们可以得到一些直观信息，例如模型像是一个什么样的参数模型，然后再用参数估计的方法来解决问题。
可以转化为线性回归的非线性回归，是否必要？
学计量时，有一种貌似很牛逼的方法，就是把一些非线性回归的问题转化为线性回归来处理。当时觉得，哈，非线性有什么，照样可以在线性框架内解决。后来发现，非线性回归向线性回归转化，只适用于极少数可以转化的情形。非线性回归自有自己的估计方法。就算某些模型可以转化为线性，线性处理的结果也比非线性直接处理的差。
来举个实际例子。在进行爆破时，需要进行设计，以免爆破对其余建筑物造成损害，理论的依据是萨道夫斯基公式。萨道夫斯基公式描述的是指质点震动速度幅值V与炸药量Q、爆心距R之间的关系。
V = K (Q ^ (1/3)/R) ^ a=K*p ^ a
不同的地形，不同环境下参数K,a并不一样，因此，爆破之前，需要首先进行爆破试验，估计参数K,a.
进行爆破试验时，设置了多组c(Q,R)值，在每一种情况下进行了爆破，并测得震动速度幅值。数据如下：
> explosion <- read.table("explode.txt", header = T, sep = " ")
> head(explosion)
Q R V
1 20 289 0.084
2 30 264 0.094
3 40 338 0.103
4 20 336 0.069
5 30 329 0.065
6 20 213 0.213
将数据转化：
> attach(explosion)
> data <- data.frame(p = Q^(1/3)/R, y = V)
> detach(explosion)
> head(data)
p y
1 0.009392449 0.084
2 0.011769820 0.094
3 0.010118201 0.103
4 0.008078624 0.069
5 0.009444476 0.065
6 0.012743745 0.213
整体数据走势大概为：
> plot(data$p, data$y)
我们直接使用非线性回归估计方法拟合：
> start <- c(K = 100, a = 1.2)
> explode.nls <- nls(y ~ K * p^a, data = data,
+ start = start)
> explode.nls
Nonlinear regression model
model: y ~ K * p^a
data: data
K a
120.254 1.349
residual sum-of-squares: 5.883
Number of iterations to convergence: 5
Achieved convergence tolerance: 7.041e-07
残差平法和为5.883。线性回归的结果：
> lm.model <- lm(y.p ~ x.p)
> sum((data$y - exp(fitted(lm.model)))^2)
[1] 67.762
线性回归的残差平方和（转换回原始数据之后）为67.762，远远超过了非线性回归的拟合结果。初步结论是，非线性回归比线性回归的效果要好。一图胜万言：
> plot(data$p, data$y, ylim = c(0, 20))
> lines(sort(data$p), fitted(explode.nls)[sort(data$p, index.return = T)$ix],
+ col = 3, lwd = 3, lty = 2)
> lines(sort(exp(x.p)), (exp(fitted(lm.model)))[sort(exp(x.p), index.return = T)$ix],
+ col = 4, lwd = 3, lty = 4)
> legend(0.04, 17, legend = c("现场爆破检测数据", "非线性回归拟合曲线", "线性回归拟合曲线"),
+ lty = c(-1, 2, 4), lwd = c(-1, 3, 3), pch = c(1, -1, -1), col = c(1, 3, 4), bty = "n")
由图可见，非线性回归分析法拟合精度更高，尤其是在振动幅值较大的情况下，其优越性更为明显。