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)
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()
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)
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()