Search code examples
rggplot2confidence-intervalnls

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

With data published in this site, about enzyme kinetics:

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)

and the code for the plot:

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

enter image description here

Is it possible to add the CI95% bands for both independent curves?

With @adiana solution it yields the next plot:

enter image description here

Solution

  • Unfortunately, as the predict for nls is a little bit tricky, there is no way to do it automatically with ggplot2 but you need to fit the model manually.

    First fit the models:

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

    Then predict the two intervals. Find the code for as.lm.nls here

    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]

    Finally use geom_ribbon to plot the intervals, I assume p is your previous fit

    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)