Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Fitting with ggplot2, geom_smooth and nls

I am trying to fit data on an exponential decay function (RC like system) with equation:

RC

My data are on the following dataframe:

dataset <- data.frame(Exp = c(4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6), t = c(0, 0.33, 0.67, 1, 1.33, 1.67, 2, 4, 6, 8, 10, 0, 33, 0.67, 1, 1.33, 1.67, 2, 4, 6, 8, 10, 0, 0.33, 0.67, 1, 1.33, 1.67, 2, 4, 6, 8, 10), fold = c(1, 0.957066345654286, 1.24139015724819, 1.62889151698633, 1.72008539595879, 1.82725412314402, 1.93164365299958, 1.9722929538061, 2.15842019312484, 1.9200507796933, 1.95804730344453, 1, 0.836176542548747, 1.07077717914707, 1.45471712491441, 1.61069357875771, 1.75576377806756, 1.89280913889538, 2.00219054189937, 1.87795513639311, 1.85242493827193, 1.7409346372629, 1, 0.840498729335292, 0.904130905000499, 1.23116185602517, 1.41897551928886, 1.60167656534099, 1.72389226836308, 1.80635095956481, 1.76640786872057, 1.74327897001172, 1.63581509884482))

I have 3 experiment (Exp: 4, 5 and 6) data I want to fit each experiment on the given equation.

I have managed to do it for of the experiment by subsetting my data and using the parameter calculated by nls

test <- subset(dataset,Exp==4)
fit1 = nls(fold ~ 1+(Vmax*(1-exp(-t/tau))),
  data=test,
  start=c(tau=0.2,Vmax=2))
ggplot(test,aes(t,fold))+
  stat_function(fun=function(t){1+coef(fit1)[[2]]*(1-exp(-t/coef(fit1)[[1]]))})+
  geom_point()

But if I try to use the geom_smooth function directly on the full dataset with this code

d <- ggplot(test,aes(t,fold))+
   geom_point()+
   geom_smooth(method="nls", 
     formula='fold~1+Vmax*(1-exp(-t/tau))',
     start=c(tau=0.2,Fmax=2))
print(d)

I get the following error:

Error in model.frame.default(formula = ~fold, data = data, weights = weight) : 
  variable lengths differ (found for '(weights)')
In addition: Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf

Is there anything wrong with my syntax? I would have this one working in order to use the same function on the dataset and using group to have one fit per Exp level.

like image 575
Oneira Avatar asked Jul 30 '14 07:07

Oneira


People also ask

What does Geom_smooth do in Ggplot?

Key R function: geom_smooth() for adding smoothed conditional means / regression line. Key arguments: color , size and linetype : Change the line color, size and type. fill : Change the fill color of the confidence region.

What does Geom_smooth mean in R?

The curve geom_smooth produces is indeed an estimate of the conditional mean function, i.e. it's an estimate of the mean distance in miles conditional on the number of trips per week (it's a particular kind of estimator called LOESS).

What is the confidence interval in Geom_smooth?

Example 2: Modify Level of Confidence Interval By default, geom_smooth() uses 95% confidence bands but you can use the level argument to specify a different confidence level.


1 Answers

There are several problems:

  1. formula is a parameter of nls and you need to pass a formula object to it and not a character.
  2. ggplot2 passes y and x to nls and not fold and t.
  3. By default, stat_smooth tries to get the confidence interval. That isn't implemented in predict.nls.

In summary:

d <- ggplot(test,aes(x=t, y=fold))+ 
         #to make it obvious I use argument names instead of positional matching
  geom_point()+
  geom_smooth(method="nls", 
              formula=y~1+Vmax*(1-exp(-x/tau)), # this is an nls argument, 
                                                #but stat_smooth passes the parameter along
              start=c(tau=0.2,Vmax=2), # this too
              se=FALSE) # this is an argument to stat_smooth and 
                        # switches off drawing confidence intervals

Edit:

After the major ggplot2 update to version 2, you need:

geom_smooth(method="nls", 
              formula=y~1+Vmax*(1-exp(-x/tau)), # this is an nls argument
              method.args = list(start=c(tau=0.2,Vmax=2)), # this too
              se=FALSE)
like image 95
Roland Avatar answered Oct 23 '22 19:10

Roland