How to plot the output from an nls model fit in gg

Question

I have some data where I would like to fit a nonlinear model to each subset of the data using nls, then superimpose the fitted models onto a graph of the data points using ggplot2. Specifically the model is of the form

y~V*x/(K+x)

which you may recognize as Michaelis-Menten. One way to do this is using geom_smooth, but if I use geom_smooth I don't have any way to retrieve the coefficients for the model fit. Alternatively I could fit the data using nls then plot lines fitted using geom_smooth, but then how do I know that the curves which geom_smooth plotted are the same as those given by my nls fit? I can't pass the coefficients from my nls fit to geom_smooth and tell it to use them unless I can get geom_smooth to only use a subset of the data, then I can specify the starting parameters so that would work, but... Every time I've tried that I get an error reading as follows:

Aesthetics must be either length 1 or the same as the data (8): x, y, colour

Here's some sample made-up data I have been using:

Concentration <- c(500.0,250.0,100.0,62.5,50.0,25.0,12.5,5.0,
                   500.0,250.0,100.0,62.5,50.0,25.0,12.5,5.0)

drug <- c(1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2)

rate <- c(1.889220,1.426500,0.864720,0.662210,0.564340,0.343140,0.181120,0.077170,
          3.995055,3.011800,1.824505,1.397237,1.190078,0.723637,0.381865,0.162771)

file<-data.frame(Concentration,drug,rate)

where Concentration will be x in my plot and rate will be y; drug will be the color variable. If I write the following I get that error:

plot <- ggplot(file,aes(x=file[,1],y=file[,3],color=Compound))+geom_point()

plot<-plot+geom_smooth(data=subset(file,file[,2]==drugNames[i]),method.args=list(formula=y~Vmax*x/(Km+x),start=list(Vmax=coef(models[[i]])[1],Km=coef(models[[i]])[2])),se=FALSE,size=0.5)

where models[[]] is a list of model parameters returned by nls.

Any ideas on how I can subset a data frame in geom_smooth so I can individually plot curves using starting parameters from my nls fit?

Answer 1

The ideal solution would plot the results of nls() using ggplot, but here's a "quick and dirty" solution based on a couple of observations.

First, you can be sure that if you use the same formula for nls() and geom_smooth(method = "nls"), you will get the same coefficients. That's because the latter is calling the former.

Second, using your example data, nls() converges to the same values of Vmax and Km (different for each drug), regardless of start value. In other words, there's no need to build models using start values in the range for each individual drug. Any of the following give the same result for drug 1 (and similarly for drug 2):

library(dplyr)
# use maximum as start
df1 %>% 
  filter(drug == 1) %>% 
  nls(rate ~ Vm * Concentration/(K + Concentration), 
      data = ., 
      start = list(K = max(.$Concentration), Vm = max(.$rate)))

# use minimum as start
df1 %>% 
  filter(drug == 1) %>% 
  nls(rate ~ Vm * Concentration/(K + Concentration), 
      data = ., 
      start = list(K = min(.$Concentration), Vm = min(.$rate)))

# use arbitrary values as start
df1 %>% 
  filter(drug == 1) %>% 
  nls(rate ~ Vm * Concentration/(K + Concentration), 
      data = ., 
      start = list(K = 50, Vm = 2))

So the quickest way to plot the curves is simply to map the drug to a ggplot aesthetic, such as color. This will construct separate nls curves from the same start values and you can then go back to nls() if required to get the coefficients, knowing that the models should be the same as the plot.

Using your example data file (but don't call it file, I used df1):

library(ggplot2)
df1 <- structure(list(Concentration = c(500, 250, 100, 62.5, 50, 25, 12.5, 5, 
                                        500, 250, 100, 62.5, 50, 25, 12.5, 5), 
                      drug = c(1, 1, 1, 1, 1, 1, 1, 1, 
                               2, 2, 2, 2, 2, 2, 2, 2), 
                      rate = c(1.88922, 1.4265, 0.86472, 0.66221, 0.56434, 0.34314, 
                               0.18112, 0.07717, 3.995055, 3.0118, 1.824505, 1.397237, 
                               1.190078, 0.723637, 0.381865, 0.162771)),
                      .Names = c("Concentration", "drug", "rate"), 
                      row.names = c(NA, -16L), 
                      class = "data.frame")

# could use e.g. Km = min(df1$Concentration) for start
# but here we use arbitrary values
ggplot(df1, aes(Concentration, rate)) + 
  geom_point() + 
  geom_smooth(method = "nls", 
              method.args = list(formula = y ~ Vmax * x / (Km + x),
                                 start = list(Km = 50, Vmax = 2)), 
              data = df1,
              se = FALSE,
              aes(color = factor(drug)))