非線性數據擬合-nls

在許多實際問題中，回歸模型中響應變量和預測變量之間的關系可能是復雜的非線性函數。這時就需要采取非線性回歸方法來建立模型。非線性回歸是在對變量的非線性關系有一定認識前提下，對非線性函數的參數進行最優化的過程，最優化后的參數會使得模型的RSS（殘差平方和）達到最小。 在R語言中最為常用的非線性回歸建模函數是nls，下面以米氏方程為例，介紹一下這個函數。

米氏方程（Michaelis-Menten equation）表示一個酶促反應的起始速度與底物濃度關系的速度方程。在酶促反應中，在低濃度底物情況下，反應相對於底物是一級反應（first order reaction）；而當底物濃度處於中間范圍時，反應（相對於底物）是混合級反應（mixed order reaction）。當底物濃度增加時，反應由一級反應向零級反應（zero order reaction）過渡。\(v_0 = \frac{V_{max}[S]}{K_m+[S]}\)這個方程稱為Michaelis-Menten方程，是在假定存在一個穩態反應條件下推導出來的，其中\(K_m\)值稱為米氏常數，\(V_{max}\)是酶被底物飽和時的反應速度，\([S]\)為底物濃度。

米氏方程擬合實例

例如：已知底物濃度數據和相應的速率，求米氏常數\(K_m\)和酶被底物飽和時的反應速度\(V_{max}\)。使用nls函數時，需要指定函數形式，並且指定參數的初始值；此外，還有一種更為簡便的方法就是采用內置自啟動模型(self-starting Models)， 此時只需要指定函數形式，而不需要指定參數初始值。

數據設定

conc <- c(2.856829, 5.005303, 7.519473, 22.101664, 27.769976, 39.198025, 45.483269, 203.784238) #底物濃度
rate <- c(14.58342, 24.74123, 31.34551, 72.96985, 77.50099, 96.08794, 96.96624, 108.88374) #速率
L.minor <- data.frame(conc, rate)
knitr::kable(head(L.minor,8))

conc	rate
2.856829	14.58342
5.005303	24.74123
7.519473	31.34551
22.101664	72.96985
27.769976	77.50099
39.198025	96.08794
45.483269	96.96624
203.784238	108.88374

數據擬合

使用nls函數擬合實驗數據，指定函數形式：M-M動力學方程，初始值設置為K=20，Vm=120。

L.minor.m1 <- nls(rate ~ Vm * conc/(K + conc), data = L.minor, #采用M-M動力學方程
                  start = list(K = 20, Vm = 120), #初始值設置為K=20，Vm=120
                  trace = TRUE) #占線擬合過程

## 624.3282 :   20 120
## 244.5458 :   15.92382 124.57149
## 234.5196 :   17.25299 126.43878
## 234.3593 :   17.04442 125.96181
## 234.3531 :   17.08574 126.04671
## 234.3529 :   17.07774 126.03017
## 234.3529 :   17.0793 126.0334
## 234.3529 :   17.07899 126.03276

summary(L.minor.m1)

## 
## Formula: rate ~ Vm * conc/(K + conc)
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## K    17.079      2.953   5.784  0.00117 ** 
## Vm  126.033      7.173  17.570 2.18e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.25 on 6 degrees of freedom
## 
## Number of iterations to convergence: 7 
## Achieved convergence tolerance: 8.137e-06

#確定x軸范圍並構建數據集
min <- range(L.minor$conc)[1]
max <- range(L.minor$conc)[2]
line.data <- data.frame(conc = seq(min, max, length.out = 1000))
#用模型預測數據構建數據集
line.data$p.predict <- predict(L.minor.m1, newdata = line.data)

可視化-判斷擬合效果

非線性回歸模型建立后需要判斷擬合效果，因為有時候參數最優化過程會捕捉到局部極值點而非全局極值點。最直觀的方法是在原始數據點上繪制擬合曲線。

require(ggplot2)

## Loading required package: ggplot2

ggplot() +
  geom_point(aes(x = conc, y = rate), data = L.minor,
             alpha = 0.5, size = 5, color = "red") +
  geom_line(aes(x = conc, y = p.predict), data = line.data,
            size = 1, color = "blue") +
  scale_x_continuous(
    name = expression(Substrate ~~ concentration(mmol ~~ m^3)),#采用expression來表示數學公式
    breaks = seq(0, 200, by = 25)) +
  scale_y_continuous(
    name = "Uptake rate (weight/h)",
    breaks = seq(0, 120, by = 10)) +
  geom_text(aes(x = 100, y = 60),
            label = "bolditalic(f(list(x, (list(K, V[m])))) == frac(V[m]%.%x, K+x))",
            #注意 geom_text中如果用expression()來進行表達，必須開啟parse = TRUE
            #同時以字符串""的形式表示，不能使用expression
            parse = TRUE, 
            size = 5, family = "times"
            ) +
  theme_bw() +
  theme(
        axis.title.x=element_text(size=16),
        axis.title.y=element_text(size=16),
        axis.text.x=element_text(size=12),
        axis.text.y=element_text(size=12))

殘差診斷

和線性回歸類似，非線性回歸假設誤差是正態、獨立和同方差性。為了檢測這些假設是否成立我們用擬合模型的殘差來代替誤差進行判斷。

plot(fitted(L.minor.m1),resid(L.minor.m1),pch=20,type="b",cex=0.6*round(abs(resid(L.minor.m1))),col="#00EEEE",xlab = "擬合數", ylab="殘差")

通過散點圖中點的大小可以判斷出擬合數據的優劣。

<embed src="http://www.zzsky.cn/flash/flash/200951029516379.swf"width="150" height="110">

Author：shangfr
郵箱：[e-mail]