ggplot: подобрать кривую (метод geom_smooth=nls) с полосами CI95%

Данные опубликованы на этом сайте. , о кинетике ферментов:

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% для обеих независимых кривых?

С решением @adiana получается следующий график:

введите здесь описание изображения


person Juanchi    schedule 13.03.2016    source источник


Ответы (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) 
person adaien    schedule 15.03.2016
comment
в функции as.lm.nls мне пришлось изменить это: fo <- as.formula(fo, env = proto:::as.proto.list(L)) - person c0bra; 28.02.2017