Search code examples
rggplot2nls

nls and log scale in ggplot2


I'm trying to plot 3 non-linear models in ggplot2. It's working in automatic scale but not in log10 scale where I get the "singular gradient error". What could be the problem(s)?

The data I'm trying to fit (df1):

x   y
4.17    0.55
10.08   0.48
40.25   0.38
101.17  0.32
400.33  0.24

The code I tried:

plot <- ggplot(df1, aes(x=x, y=y))+
  stat_smooth(method="nls",
              formula=y~I(a*x^(-n)),
              data=df1,
              start=list(a=1,n=1),
              se=FALSE,
              colour="red")+
  stat_smooth(method="nls",
              formula=y~m*(x+m^(1/n))^(-n),
              data=df1,
              start=list(m=0.7, n=0.17),
              se=FALSE,
              colour="blue")+
  stat_smooth(method="nls",
              formula=y~m*(x+m^(1/n))^(-n)+b,
              data=df1,
              start=list(m=0.7, n=0.17, b=1),
              se=FALSE,
              colour="green")+
  geom_point()+
  scale_x_log10()+
  scale_y_log10()+
  theme_bw()
plot

Solution

  • The problem seems to be that when you specify scale_x_log10 or scale_y_log10, the values of your data are transformed before being passed along to the different stats or geoms. This means while your nls may work on the untransformed data, it does not work on the log-transformed data.

    #OK
    nls(y~m*(x+m^(1/n))^(-n), df1, start=list(m=0.7, n=0.17))
    #NOT OK
    nls(y~m*(x+m^(1/n))^(-n), log10(df1), start=list(m=0.7, n=0.17))
    

    There doesn't seem to be much you can do in ggplot2 to fix this. Instead, you could fit the NLS models ahead of time on the untransformed scale and just plot the results with ggplot2. For example

    mods<-list(
        list(y~I(a*x^(-n)), list(a=1,n=1)),
        list(y~m*(x+m^(1/n))^(-n), list(m=0.7, n=0.17)),
        list(y~m*(x+m^(1/n))^(-n)+b, list(m=0.7, n=0.17, b=1))
    )
    
    fits<-lapply(mods, function(x, xr) {
        mod<-nls(x[[1]], data=df1, start=x[[2]])
        xx<-seq(xr[1], xr[2], length.out=100)
        yy<-predict(mod, newdata=data.frame(x=xx))
        data.frame(x=xx, y=yy)
    }, xr=range(df1$x))
    
    library(ggplot2)
    ggplot(df1, aes(x=x, y=y))+
      geom_line(data=fits[[1]], color="red") +
      geom_line(data=fits[[2]], color="blue") +
      geom_line(data=fits[[3]], color="green") +
      geom_point()+
      scale_x_log10()+
      scale_y_log10()+
      theme_bw()
    

    will produce

    enter image description here