# 非线性数据拟合-nls

## 米氏方程拟合实例

### 数据设定

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))

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

concrate

### 数据拟合

L.minor.m1 <- nls(rate ~ Vm * conc/(K + conc), data = L.minor, #采用M-M动力学方程start = list(K = 20, Vm = 120), #初始值设置为K=20，Vm=120trace = 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#同时以字符串""的形式表示，不能使用expressionparse = 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="残差")

Author：shangfr邮箱：shangfr@foxmail.com