R語言用nls做非線性回歸以及函數模型的參數估計
非線性回歸是在對變量的非線性關系有一定認識前提下，對非線性函數的參數進行最優化的過程，最優化後的參數會使得模型的RSS（殘差平方和）達到最小。在R語言中最為常用的非線性回歸建模函數是nls，下面以car包中的USPop數據集為例來講解其用法。數據中population表示人口數，year表示年份。如果將二者繪制散點圖可以發現它們之間的非線性關系。在建立非線性回歸模型時需要事先確定兩件事，一個是非線性函數形式，另一個是參數初始值。
一、模型擬合
對於人口模型可以采用Logistic增長函數形式，它考慮了初期的指數增長以及總資源的限制。其函數形式如下。
首先載入car包以便讀取數據，然後使用nls函數進行建模，其中theta1、theta2、theta3表示三個待估計參數，start設置了參數初始值，設定trace為真以顯示叠代過程。nls函數默認采用Gauss-Newton方法尋找極值，叠代過程中第一列為RSS值，後面三列是各參數估計值。然後用summary返回回歸結果。
library(car) pop.mod1 <- nls(population ~ theta1/(1+exp(-(theta2+theta3*year))),start=list(theta1 = 400, theta2 = -49, theta3 = 0.025), data=USPop, trace=T) summary(pop.mod)
還有一種更為簡便的方法就是采用內置自啟動模型(self-starting Models)，此時我們只需要指定函數形式，而不需要指定參數初始值。本例的logistic函數所對應的selfstarting函數名為SSlogis
pop.mod2 <- nls(population ~ SSlogis(year,phi1,phi2,phi3),data=USPop)
二、判斷擬合效果
非線性回歸模型建立後需要判斷擬合效果，因為有時候參數最優化過程會捕捉到局部極值點而非全局極值點。最直觀的方法是在原始數據點上繪制擬合曲線。
library(ggplot2) p <- ggplot(USPop,aes(year, population)) p+geom_point(size=3)+geom_line(aes(year,fitted(pop.mod1)),col=‘red‘)
附註：關於fitted詳細講解轉——一個不錯的博客
若比較多個模型的擬合效果可使用AIC函數，取最小值為佳。（AIC是赤池系數用於比較模型的好壞，類似的BIC是貝葉斯系數）
三、殘差診斷
為了檢測這些假設是否成立我們用擬合模型的殘差來代替誤差進行判斷。
plot(fitted(pop.mod1) , resid(pop.mod1),type=‘b‘)
fitted是擬合值(predict是預測值) resid和residuals表示殘差
四、函數模型的參數估計
關於函數估計至少有這麽幾個問題是需要關心的：
1、知道函數的一個大概的模型，需要估計函數的參數；
2、不知道模型，但想用一個不壞的模型刻畫它；
3、不知道模型，不太關心其顯式表達是什麽，只想知道它在沒觀測到的點的取值。
這三個問題第一個是擬合或者叫參數估計，第二個叫函數逼近，第三個叫函數插值。從統計的角度來看，第一個是參數問題，剩下的是非參數的問題
以含常數項的指數函數為例
模擬模型( y=x^beta+mu +varepsilon ),這裏假設( beta=3,mu=5.2 )
產生仿真數據
len <- 24 x <- runif(len) y <- x^3 + 5.2 + rnorm(len, 0, 0.06) ds <- data.frame(x = x, y = y) str(ds)
‘data.frame‘: 24 obs. of 2 variables:
$ x: num 0.3961 0.2004 0.0407 0.9873 0.83 ...
$ y: num 5.37 5.15 5.21 6.06 5.75 ...
plot(y ~ x, main = "指數模型") s <- seq(0, 1, length = 100) lines(s, s^3, lty = 2, col = "red")
使用nls函數估計如下：
rhs <- function(x, b0, b1) {
b0 + x^b1
}
m.1 <- nls(y ~ rhs(x, intercept, power), data = ds, start = list(intercept = 0,
power = 2), trace = T)
回顯如下：
629.9495 : 0 2
0.08918652 : 5.174334 2.526742
0.08346069 : 5.184992 2.786072
0.08303884 : 5.188992 2.870896
0.08301127 : 5.190252 2.894080
0.08300947 : 5.190594 2.900075
0.08300935 : 5.190683 2.901602
0.08300934 : 5.190705 2.901989
0.08300934 : 5.190711 2.902088
0.08300934 : 5.190713 2.902112
summary(m.1)
回顯如下：
Formula: y ~ rhs(x, intercept, power)
Parameters:
Estimate Std. Error t value Pr(>|t|)
intercept 5.19071 0.02189 237.150 < 2e-16
power 2.90211 0.30825 9.415 3.58e-09
intercept ***
power ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.06143 on 22 degrees of freedom
Number of iterations to convergence: 9
Achieved convergence tolerance: 4.296e-06
如果采用最小二乘估計方法，得到的結果是：
model <- lm(I(log(y)) ~ I(log(x))) summary(model)
回顯如下：
Call:
lm(formula = I(log(y)) ~ I(log(x)))
Residuals:
Min 1Q Median 3Q Max
-0.054991 -0.029944 -0.008994 0.037530 0.073063
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.755422 0.011128 157.755 < 2e-16
I(log(x)) 0.055445 0.009502 5.835 7.17e-06
(Intercept) ***
I(log(x)) ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.03851 on 22 degrees of freedom
Multiple R-squared: 0.6075, Adjusted R-squared: 0.5897
F-statistic: 34.05 on 1 and 22 DF, p-value: 7.166e-06
我們可以將估計數據、真實模型、nls估計模型、最小二乘模型得到的結果展示在下圖中，來擬合好壞有個直觀的判斷：
plot(ds$y ~ ds$x, main = "Fitted power model, with intercept", sub = "Blue: fit; magenta（洋紅）: fit LSE ; green: known") lines(s, s^3 + 5.2, lty = 2, col = "green") lines(s, predict(m.2, list(x = s)), lty = 1, col = "blue") lines(s, exp(predict(model, list(x = s))), lty = 2, col = "magenta") segments(x, y, x, fitted(m.2), lty = 2, col = "red")
legend("topleft",c("nsl擬合","最小二乘法(洋紅)","真實","nls估計"),col=c("blue","magenta","green","red"),pch=15:15,cex = 0.7)
R語言用nls做非線性回歸以及函數模型的參數估計