0

I have a non-linear model ExpDec2RC fitted to the following data.

data <- read.table("https://gitlab.com/-/snippets/1999226/raw", sep = "\t")
colnames(data) <- c("irf", "decay")
data$t <- 1:nrow(data)

data <- data[367:1279,]
data$t <- data$t - 366

The model uses fft() for convolution and is as follows.

ExpDec2RC <- function(t, irf, i0, shift, alpha1, tau1, alpha2, tau2) {
  
  irfbins <- seq_along(t) - 1 # start from 0
  irfcounts <- irf
  n <- length(irf)
  
  irfshifted <- ((1 - shift + floor(shift)) *
                   irfcounts[((((irfbins - floor(shift) - 1) %% n) + n) %% n) + 1]) +
    ((shift - floor(shift)) *
       irfcounts[((((irfbins - ceiling(shift) - 1) %% n) + n) %% n) + 1])
  irfnormalised <- irfshifted / sum(irfshifted)
  
  ft_ymodel = fft((alpha1 * exp(-t / tau1)) + (alpha2 * exp(-t / tau2)))
  ft_irfnormalised = fft(irfnormalised)
  
  ftcp = ft_ymodel*ft_irfnormalised
  
  i0 + Re(fft(z = ftcp, inverse = TRUE)/length(ftcp))
}

With a set of fitted parameters, I can get the fitted values as follows.

data$fit <- ExpDec2RC(t = data$t,
                       irf = data$irf,
                       i0 = 11.7764066,
                       shift = 2.00000973,
                       alpha1 =  58414.9896,
                       tau1 = 1.85542518,
                       alpha2 = 1042.19136,
                       tau2 = 75.4146475)

But, when I use curve() to plot the same, I am getting a curve completely different from the fitted values. Why is that and how to get the correct fitted values with curve()?

plot(data$t, data$decay, col="red", pch="o", log="y")
points(data$t, data$irf, col="blue", pch="o")
# plot fitted values
lines(data$t, data$fit, col="cyan", lwd = 5)
# plot with curve()
curve(ExpDec2RC(t = x,
                irf = data$irf,
                i0 = 11.7764066,
                shift = 2.00000973,
                alpha1 =  58414.9896,
                tau1 = 1.85542518,
                alpha2 = 1042.19136,
                tau2 = 75.4146475), col = "yellow", add = TRUE, lwd = 5)
legend(600, 9000, legend=c("Fitted values", "With curve()"),
       col=c("cyan", "yellow"), lty = 1, lwd = 5)

enter image description here

The same is also true with stat_function() in ggplot2

ggplot(data, aes(t, decay)) +
  geom_point(alpha =0.25, colour = "red") +
  geom_point(aes(t, irf), alpha =0.25, colour = "blue") +
  geom_point(aes(t, fit), alpha =0.25, colour = "cyan") +
  scale_y_log10() +
  stat_function(fun = function (x) ExpDec2RC(t = x,
                                             irf = data$irf,
                                             i0 = 11.7764066,
                                             shift = 2.00000973,
                                             alpha1 =  58414.9896,
                                             tau1 = 1.85542518,
                                             alpha2 = 1042.19136,
                                             tau2 = 75.4146475),
                colour = "yellow", inherit.aes = FALSE, alpha = 0.5, size= 2) +
  theme_bw()

enter image description here

With a simple model, the fitted values as well as the curve overlap perfectly.

ExpDec2 <- function(t, i0, alpha1, tau1, alpha2, tau2) {
  
  i0 + (alpha1 * exp(-t / tau1)) + (alpha2 * exp(-t / tau2))
}

data$fit2 <- ExpDec2(t = data$t,
                     i0 = 11.63431,
                     alpha1 = 6618.36,
                     tau1 = 8.475745201,
                     alpha2 = 644.1506,
                     tau2 = 103.1719055)

plot(data$t, data$decay, col="red", pch="o", log="y")
# plot fitted values
lines(data$t, data$fit2, col="cyan", lwd = 8)
# plot with curve()
curve(ExpDec2(t = x,
                i0 = 11.63431,
                alpha1 = 6618.36,
                tau1 = 8.475745201,
                alpha2 = 644.1506,
                tau2 = 103.1719055), col = "yellow", add = TRUE, lwd = 3)
legend(600, 9000, legend=c("Fitted values", "With curve()"),
       col=c("cyan", "yellow"), lty = 1, lwd = 5)

enter image description here

ggplot(data, aes(t, decay)) +
  geom_point(alpha =0.25, colour = "red") +
  geom_point(aes(t, fit2), alpha =0.25, colour = "cyan") +
  scale_y_log10() +
  stat_function(fun = function (x) ExpDec2(t = x,
                                             i0 = 11.63431,
                                             alpha1 = 6618.36,
                                             tau1 = 8.475745201,
                                             alpha2 = 644.1506,
                                             tau2 = 103.1719055),
                colour = "yellow", inherit.aes = FALSE, alpha = 0.5, size= 2) +
theme_bw()

enter image description here

Crops
  • 5,024
  • 5
  • 38
  • 65
  • On running your code, I get: `Error in strwidth(legend, units = "user", cex = cex, font = text.font) : plot.new has not been called yet`. In `curve(ExpDec2RC(t = x,`, where is `x` defined? – Limey Mar 20 '21 at 09:49
  • You have to run all the code chunks before the one with `plot`. In `curve(ExpDec2RC(t = x`, as `add=TRUE` `x` would be the same as the one the previous plot ie. `data$t`. – Crops Mar 20 '21 at 10:22
  • That’s what did. – Limey Mar 20 '21 at 14:25
  • @Limey I am not getting the error. So alternatively I have added equivalent `ggplot2` code. – Crops Mar 20 '21 at 17:28

0 Answers0