R語(yǔ)言非線性回歸入門(mén)
作者 Lionel Hertzog
在一簇散點(diǎn)中擬合一條回歸線（即線性回歸）是數(shù)據(jù)分析的基本方法之一。有時(shí)，線性模型能很好地?cái)M合數(shù)據(jù)，但在某些（很多）情形下，變量間的關(guān)系未必是線性的。這時(shí)，一般有三類方法解決這個(gè)問(wèn)題: (1) 通過(guò)變換數(shù)據(jù)使得其關(guān)系線性化, (2) 用多項(xiàng)式或者比較復(fù)雜的樣條來(lái)擬合數(shù)據(jù), (3) 用非線性函數(shù)來(lái)擬合數(shù)據(jù)
從標(biāo)題你應(yīng)該已經(jīng)猜到非線性回歸是本文的重點(diǎn)
什么是非線性回歸
在非線性回歸中，分析師通常采用一個(gè)確定的函數(shù)形式和相應(yīng)的參數(shù)來(lái)擬合數(shù)據(jù)。最常用的參數(shù)估計(jì)方法是利用非線性最小二乘法（R中的nls函數(shù)）。該方法使用線性函數(shù)來(lái)逼近非線性函數(shù)，并且通過(guò)不斷迭代這個(gè)過(guò)程來(lái)得到參數(shù)的最優(yōu)解（本段來(lái)自維基百科）。非線性回歸的良好性質(zhì)之一是估計(jì)出的參數(shù)都有清晰的解釋（如Michaelis-Menten模型的Vmax是指最大速率），而變換數(shù)據(jù)后得到的線性模型其參數(shù)往往難以解釋。
非線性最小二乘擬合
首先，我們以Michaelis-Menten方程為例。
# 生成一些仿真數(shù)據(jù)
set.seed(20160227)
x <- seq(0, 50, 1)
y <- ((runif(1, 10, 20)*x)/(runif(1, 0, 10)+x)) + rnorm(51, 0, 1)
# 對(duì)于一些簡(jiǎn)單的模型，nls函數(shù)可以自動(dòng)找到合適的參數(shù)初值
m <- nls(y ~ a*x/(b+x))
# 計(jì)算模型的擬合優(yōu)度
cor(y, predict(m))
[1] 0.9496598
# 將結(jié)果可視化
plot(x, y)
lines(x, predict(m), lty = 2, col = "red", lwd = 3)
輸出的圖片如下：
選擇適宜的迭代初值
在非線性回歸中，找到合適的迭代初值對(duì)于整個(gè)模型算法的收斂性而言至關(guān)重要。假如你設(shè)定的參數(shù)初值完全脫離了其潛在的取值范圍，迭代算法可能不收斂或者返回一些沒(méi)有意義的參數(shù)值。比如返回一個(gè)大小為1000的增長(zhǎng)率，但其真值卻是1.04。尋找合適初值的最好辦法是“緊盯著”數(shù)據(jù)，繪制相應(yīng)圖表并結(jié)合你對(duì)方程的理解來(lái)確定參數(shù)的合適初值。
# 生成仿真數(shù)據(jù)，并且此次對(duì)于參數(shù)沒(méi)有先驗(yàn)信息
y <- runif(1, 5, 15)*exp(-runif(1, 0.01, 0.05)*x)+rnorm(51, 0, 0.5)
# 可視化數(shù)據(jù)并選擇一些參數(shù)初值
plot(x, y)
# 通過(guò)這個(gè)散點(diǎn)圖確定參數(shù)a, b的初值
a_start <- 8 # 參數(shù)a是x = 0時(shí)y的取值
b_start<- 2*log(2)/a_start # b 是衰減速率
# 擬合模型
m <- nls(y ~ a*exp(-b*x), start = list(a = a_start, b = b_start))
# 計(jì)算擬合優(yōu)度
cor(y, predict(m))
[1] 0.9811831
# 將結(jié)果可視化
lines(x, predict(m), col = "red", lty = 2, lwd = 3)
輸出的圖片如下
使用自啟動(dòng)函數(shù)
不同的科學(xué)研究領(lǐng)域會(huì)對(duì)同一個(gè)模型設(shè)定不同的參數(shù)形式（即不同的方程），比如研究人口增長(zhǎng)的邏輯斯蒂模型，在生態(tài)學(xué)中一般采用如下形式：
等式中的 代表時(shí)間t時(shí)的個(gè)體數(shù)， 是個(gè)體增長(zhǎng)速率， 是環(huán)境承載能力。我們可以將這個(gè)等式改寫(xiě)為微分方程的形式：
library(deSolve)
# 利用邏輯斯蒂模型生成人口增長(zhǎng)的仿真數(shù)據(jù)，并用nls估計(jì)參數(shù)
log_growth <- function(Time, State, Pars) {
with(as.list(c(State, Pars)), {
dN <- R*N*(1-N/K)
return(list(c(dN)))
})
}
# 邏輯斯蒂增長(zhǎng)的參數(shù)
pars <- c(R = 0.2, K = 1000)
# 設(shè)定初值
N_ini <- c(N = 1)
# 常微分方程的時(shí)間階段（下標(biāo)t）
times <- seq(0, 50, by = 1)
# 常微分方程
out <- ode(N_ini, times, log_growth, pars)
# 添加一些隨機(jī)波動(dòng)
N_obs <- out[, 2]+rnorm(51, 0, 50)
# 個(gè)體數(shù)值不能小于1
N_obs <- ifelse(N_obs<1, 1, N_obs)
# 畫(huà)圖
plot(times, N_obs)
這部分代碼只是生成了帶有隨機(jī)誤差的仿真數(shù)據(jù)，接下來(lái)的部分會(huì)展現(xiàn)估計(jì)參數(shù)初值的技巧。R語(yǔ)言中有一個(gè)估計(jì)邏輯斯蒂方程參數(shù)的內(nèi)建函數(shù)（SSlogis），但它使用的是如下方程：
# 尋找方程的參數(shù)
SS <- getInitial(N_obs ~ SSlogis(times, alpha, xmid, scale), data = data.frame(N_obs = N_obs, times = times))
我們可使用getInitial函數(shù)來(lái)對(duì)模型參數(shù)做一個(gè)基于數(shù)據(jù)的初步估計(jì)。然后把該函數(shù)的輸出作為一個(gè)向量化參數(shù)傳遞給自啟動(dòng)函數(shù)（SSlogis），同時(shí)也將無(wú)引號(hào)的三個(gè)參數(shù)名賦值給邏輯斯蒂方程（譯者注：即alpha，xmid，scale三個(gè)參數(shù)）。
然而，由于SSlogis的參數(shù)設(shè)定有些不同，我們需要對(duì)SSlogis的輸出值做一些處理，使得其與邏輯斯蒂方程中的形式一致。
# 改變參數(shù)形式
K_start <- SS["alpha"]
R_start <- 1/SS["scale"]
N0_start <- SS["alpha"]/(exp(SS["xmid"]/SS["scale"])+1)
# 構(gòu)建模型的公式
log_formula <- formula(N_obs ~ K*N0*exp(R*times)/(K + N0*(exp(R*times) - 1)))
# 擬合模型
m <- nls(log_formula, start = list(K = K_start, R = R_start, N0 = N0_start))
# 估計(jì)參數(shù)
summary(m)
Formula: N_obs ~ K * N0 * exp(R * times)/(K + N0 * (exp(R * times) - 1))
Parameters:
Estimate Std. Error t value Pr(>|t|)
K 1.012e+03 3.446e+01 29.366 <2e-16 ***
R 2.010e-01 1.504e-02 13.360 <2e-16 ***
N0 9.600e-01 4.582e-01 2.095 0.0415 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 49.01 on 48 degrees of freedom
Number of iterations to convergence: 1
Achieved convergence tolerance: 1.537e-06
# 計(jì)算擬合優(yōu)度
cor(N_obs,predict(m))
[1] 0.9910316
# 結(jié)果可視化
lines(times, predict(m), col = "red", lty = 2, lwd = 3)
輸出圖形如下:
修改SSlogis輸出的參數(shù)結(jié)構(gòu)確實(shí)有些繁瑣，不過(guò)值得一試。
在后續(xù)的文章中，除了非線性最小二乘法，我們將利用更為可靠和強(qiáng)大的極大似然估計(jì)法來(lái)擬合模型。它能使你構(gòu)建你能想到的任何模型。
原文刊載于datascience+網(wǎng)站
|
