Search code examples
rggplot2graphsmoothingloess

How do I get a smooth curve from a few data points, in R?


I am trying to plot the rate 1/t as it changes with mue. The code is given below and I have highlighted the relevant lines with input and output.


library("deSolve")
library("reshape")
library("tidyverse")


Fd <- data.frame()
MUES <- c(100, 1000, 2000, 5000, 10000, 20000, 50000, 100000, 100010, 100020, 100050, 100060, 100080, 100090, 100100, 100500) # <------ THIS IS THE INPUT

for (i in 1:length(MUES)){
parameters <- c(tau = 0.005, tau_r = 0.0025, mui=0, Ve=0.06, Vi=-0.01, s=0.015, mue=MUES[i])
state <- c(X = 0.015, Y = 0)
Derivatives <-function(t, state, parameters) {
      #cc <- signal(t)
 with(as.list(c(state, parameters)),{
 # rate of change
 dX <- -(1/tau + mue - mui)*X + (Y-X)/tau_r + mue*Ve - mui*Vi
 dY <- -Y/tau + (X-Y)/tau_r

 # return the rate of change
 list(c(dX, dY))
 }) # end with(as.list ...
 }
times <- seq(0, 0.1, by = 0.0001)
out <- ode(y = state, times = times, func = Derivatives, parms = parameters)

out.1 <- out %>% 
  as.data.frame() %>% summarise(d = min(times[Y >=0.015]))
Time <- out.1$d
    localdf <- data.frame(t=Time, rate= 1/Time, input=MUES[i])
Fd <- rbind.data.frame(Fd, localdf)}. # <----- THIS IS THE DATAFRAME WITH OUTPUT AND INPUT

spline_int <- as.data.frame(spline(Fd$input, Fd$rate))
ggplot(Fd) + 
  geom_point(aes(x = input, y = rate), size = 3) +
  geom_line(data = spline_int, aes(x = x, y = y))

The rate 1/t has a limiting value at 1276 and thats why I have taken quite a few values of mue in the end, to highlight this. I get a graph like this:

first

What I want is something like below, so I can highlight the fact that the rate 1/t doesn't grow to infinity and infact has a limiting value. The below figure is from the Python question.

second

How do I accomplish this in R? I have tried loess and splines and geom_smooth (but just with changing span), perhaps I am missing something obvious.


Solution

  • Splines are polynomials with multiple inflection points. It sounds like you instead want to fit a logarithmic curve:

    # fit a logarithmic curve with your data
    logEstimate <- lm(rate~log(input),data=Fd)
    
    # create a series of x values for which to predict y 
    xvec <- seq(0,max(Fd$input),length=1000)
    
    # predict y based on the log curve fitted to your data
    logpred <- predict(logEstimate,newdata=data.frame(input=xvec))
    
    # save the result in a data frame
    # these values will be used to plot the log curve 
    pred <- data.frame(x = xvec, y = logpred)
    
    ggplot() + 
      geom_point(data = Fd, size = 3, aes(x=input, y=rate)) +
      geom_line(data = pred, aes(x=x, y=y))
    

    Result: enter image description here

    I borrowed some of the code from this answer.