r - ggplot : fit a curve (geom_smooth method ="nls") with CI95% bands

Question

数据发布于this网站，关于 enzyme 动力学:

Enz <- c("WT","WT","WT","WT","WT",
         "WT","WT","WT","WT","WT",
         "WT","WT","WT",
         "H297F","H297F","H297F",
         "H297F","H297F","H297F",
         "H297F","H297F")
S <- c(2.00, 1.00, 0.60, 0.50, 0.40, 
         0.30, 0.20, 0.10, 0.09, 0.08, 
         0.06, 0.04, 0.02, 
         0.05, 0.10, 0.20, 
         0.30, 0.40, 0.50, 
         1.00, 2.00)
v <- c(59.01, 58.29, 54.17, 51.82, 49.76, 
         45.15, 36.88, 26.10, 23.50, 22.26, 
         16.45, 13.67, 6.14, 
         11.8, 19.9, 30.3, 
         36.6, 40.2, 42.1, 
         47.8, 50.0)

以及绘图代码:

ggplot(data=enzdata,         
            aes(x=S,            
            y=v,            
            colour = Enz)) +  
            geom_point() +            
            xlab("Substrate (mM)") +  
            ylab("Velocity (uM/min/mg.enzyme)") +    
            ggtitle("Glucose Dehydrogenase \n wild type and mutant") +  
            geom_smooth(method = "nls", 
                        formula = y ~ Vmax * x / (Km + x), 
                        start = list(Vmax = 50, Km = 0.2),
                        se = F, size = 0.5, 
                        data = subset(enzdata, Enz=="WT")) +
            geom_smooth(method = "nls", 
                        formula = y ~ Vmax * x / (Km + x), 
                        start = list(Vmax = 50, Km = 0.2),
                        se = F, size = 0.5, 
                        data = subset(enzdata, Enz=="H297F"))

是否可以为两条独立曲线添加 CI95% strip ？

使用 @adiana 解决方案，它会产生下一个图:

Answer 1

不幸的是，由于 nls 的预测有点棘手，因此无法使用 ggplot2 自动进行预测，而是需要手动拟合模型。

首先拟合模型:

library(nls2)
nsmodel1<-nls(formula = v ~ Vmax * S / (Km + S),data=subset(enzdata, Enz=="WT"),start = list(Vmax = 50, Km = 0.2))
nsmodel2<-nls(formula = v ~ Vmax * S / (Km + S),data=subset(enzdata, Enz=="H297F"),start = list(Vmax = 50, Km = 0.2))

然后预测两个区间。在此处查找 as.lm.nls 的代码

http://www.leg.ufpr.br/~walmes/cursoR/ciaeear/as.lm.R

fit1<-predict(as.lm.nls(nsmodel1), interval = "confidence") 
fit2<-predict(as.lm.nls(nsmodel2), interval = "confidence") 
enzdata$lowerfit[enzdata$Enz=="WT"]<-fit1[,2]
enzdata$upperfit[enzdata$Enz=="WT"]<-fit1[,3]
enzdata$lowerfit[enzdata$Enz=="H297F"]<-fit2[,2]
enzdata$upperfit[enzdata$Enz=="H297F"]<-fit2[,3]

最后使用 geom_ribbon 绘制间隔，我假设 p 是您之前的拟合

p+geom_ribbon(aes(x=S,ymin=lowerfit,ymax=upperfit),data=subset(enzdata, Enz=="WT"),alpha=0.5)+
  geom_ribbon(aes(x=S,ymin=lowerfit,ymax=upperfit),data=subset(enzdata, Enz=="H297F"),alpha=0.5)