Decay functions and ggplot - How do I fit a decay function to my data, with conf intervals?

Question

I am working with concentration over time data. The intent is to find a model for the decay rate of the phenol content of algae over time. Not much info is available on the literature about the shape of decay of these phenols, but my data showed a linear decay for one species and one non-linear for another species (see scatterplots below).

I want to use a decay function of the type: y = ae^kt, where y is the final concentration at time t, a is the initial concentration, and k is the decay constant. This is what I believe is theoretically the correct function to use for the 'Carpo' species, since phenols are degrading over time in an exponential fashion (see graph)

Here is my data for Carpo spp:

structure(list(Time = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 24, 
24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 72, 72, 72, 72, 72, 
72, 72, 72, 72, 72, 72, 72, 168, 168, 168, 168, 168, 168, 168, 
168, 168, 168, 168, 168), Phenols = c(155.05396358031, 292.448819046919, 
392.09254588543, 373.777907851168, 158.870919593217, 227.401160002983, 
158.545124752453, 348.038178016855, 196.004621459891, 198.711112888521, 
214.843845072805, 125.476444729792, 115.269817901552, 50.5754768174938, 
182.514953753076, 413.728300560148, 90.3366033152076, 116.823625331981, 
145.874402114001, 18.4955652722213, 52.868657898569, 130.266434542432, 
164.41413201994, 73.7144602576443, 15.0853972561612, 37.200244360702, 
32.0793538786566, 27.2060777211984, 26.4422227157876, 25.2930201395159, 
30.7405411687273, 9.67833536375294, 19.8968659319271, 28.3441845349511, 
18.9431364122308, 25.34674890162, 3.72857179653475, 10.5561871026197, 
14.1864805110413, 0.521827126052154, 7.4465801315315, 6.77057861607754, 
2.26713632376336, -1.90266594684443, 3.61534055996961, 3.42587396256444, 
3.12017043262619, 2.30896861408144)), row.names = c(NA, -48L), class = c("tbl_df", 
"tbl", "data.frame"))

That gives this scatterplot:

Carpo spp phenol concentration over time

The Macro spp has this scatterplot- just for reference:

Macro spp phenol concentration over time, with fitten linear regression (y ~ x)

So what I am trying to do, is to find line of best fit + 95% CI for the Carpo data, with a line of the type y = ae^kt. The only way I have thought of to do this so far, is by log-transforming the data and fitting a linear regression (with abs values of phenols as some negative values/zeros exist):

ggplot(Carpo.exp.data, aes(x=log(Time), y=log(abs(Phenols))))+
  geom_point()+
  theme_minimal()+
  geom_smooth(method = "lm")

which gives this: log-transformed concentration of phenols over time (Carpo spp)

To look at the coefficient estimates I fit a linear model and look at the summary (excluding the first 12 rows of data, as Time = 0, which produces an error for the log conversion):

> m1<-lm(log(abs(Phenols))~log(Time), Carpo.exp.data[13:46,])
> summary(m1)

Call:
lm(formula = log(abs(Phenols)) ~ log(Time), data = Carpo.exp.data[13:46, 
    ])

Residuals:
    Min      1Q  Median      3Q     Max 
-2.1706 -0.2753  0.1180  0.4235  1.3116 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   9.9291     0.6810   14.58 1.10e-15 ***
log(Time)    -1.6411     0.1617  -10.15 1.57e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7411 on 32 degrees of freedom
Multiple R-squared:  0.763, Adjusted R-squared:  0.7556 
F-statistic:   103 on 1 and 32 DF,  p-value: 1.571e-11

But now I feel like I am stuck at this stage and cannot figure out how to go from this linear model back to the exponential decay using the non-transformed data, as I want:

the estimated decay rates (k) from my data,
a graph with the non-transformed data overlayed with the exponential decay curve derived from the log trasnformed lm coefficients.

Does this make sense? Can someone help me out? Thank you - and sorry in advance if I have forgotten to include something.

Vero Rodriguez-j · Accepted Answer · 2023-02-22 05:03:55Z

Update!! I think I have figured out a way to do this. Posting it in case someone else has a similar issue in the future. Please let me know if something seems off.

I figured that by taking the ln of my original phenol data (y), the decay equation becomes: ln(y) = ln(a) + k * x

So by fitting a linear regression on my log-transformed data, I can estimate the values for a and k so:

lm.carpo<-lm(log(abs(Phenols))~Time, Carpo.exp.data)
summary(lm.carpo)

we get the results:

Call: lm(formula = log(abs(Phenols)) ~ Time, data = Carpo.exp.data)

Residuals: Min 1Q Median 3Q Max -1.80722 -0.32138 -0.01005 0.27216 1.49549

Coefficients: Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.203863 0.137204 37.93 <2e-16 ***
Time -0.024090 0.001489 -16.18 <2e-16 *** --- Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6636 on 46 degrees of freedom Multiple R-squared: 0.8506, Adjusted R-squared: 0.8474 F-statistic: 261.9 on 1 and 46 DF, p-value: < 2.2e-16

so, the estimated coefficients are: ln(a) = 5.203863 +- 0.137204 k = -0.024090 +- 0.001489

this means that the estimated decay rate curve for this species is: y = exp(5.2039)exp(-0.0241t), where t=x

Then, I figured out how to get a graph of this curve with 95% Confidence Intervals with ggplot, by making a ribbon under a scatterplot of the original data. Here is how I did it:

1- Get Confidence Interval values for both coefficients (a and k)

confint(lm.carpo, 'Time', level = 0.95)
confint(lm.carpo, "(Intercept)", level = 0.95)

confint(lm.carpo, 'Time', level = 0.95)

2.5 % 97.5 % Time -0.02708595 -0.02109341

confint(lm.carpo, "(Intercept)", level = 0.95)

2.5 % 97.5 % (Intercept) 4.927686 5.480041

Then, I make functions representing the estimated decay curve, and an upper and lower CI curve:

curveCarpo<- function(v) {
  exp(5.2039)*exp(-0.0241*v)}
  
curveCarpo.upper <- function(v){
  exp(5.4800)*exp(-0.0211*v)}

curveCarpo.lower <- function(v){
  exp(4.9277)*exp(-0.0271*v)}

And lastly, I use a vector of x and y to generate a smooth area for the ribbon. Then I make a ggplot with the ribbon first (using the functions I made above to define the y-min and y-max), and then add the other layers: the estimated decay curve (as a stat_function) and the datapoints (as geom_points).

ribbon.x.carpo<- c(0:168, 0.05)
ribbon.y.carpo<- curveCarpo(ribbon.x.carpo)
ribbon.vals<-data.frame(ribbon.x.carpo, ribbon.y.carpo)

p.carpo.curves<- ggplot(Carpo.exp.data, aes(x=Time, y=Phenols))+ 
  geom_ribbon(data=ribbon.vals, mapping=aes(x=ribbon.x.carpo, 
                                            y=ribbon.y.carpo, 
                                            ymin = curveCarpo.lower(ribbon.x.carpo), 
                                            ymax = curveCarpo.upper(ribbon.x.carpo)),
                                            fill = "#FDB863",alpha = 0.4)+
  stat_function(fun = "curveCarpo", color="#E66101", lwd=1.6)+
  geom_point(color="black",lwd=2.5) 
p.carpo.curves

Which yields this graph:

Decay rate curve estimate with 95% CI

Hope this helps someone, as it took me several days of work to figure it out! Cheers

Collectives™ on Stack Overflow

Decay functions and ggplot - How do I fit a decay function to my data, with conf intervals?

1 Answer 1

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Comments

Your Answer

Sign up or log in

Post as a guest

Related