非线性数据拟合-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="残差")

通过散点图中点的大小可以判断出拟合数据的优劣。