Maximum rate of curve change on a smoothed curve

Question

I have a time series dataset and I am trying to calculate the maximum rate of change to estimate greenup date on NDVI data. My data is below:

    date       NDVI
 1 2021-01-01 0.307
 2 2021-01-17 0.291
 3 2021-02-02 0.287
 4 2021-02-18 0.247
 5 2021-03-06 0.286
 6 2021-03-22 0.268
 7 2021-04-07 0.299
 8 2021-04-23 0.295
 9 2021-05-09 0.349
10 2021-05-25 0.402
11 2021-06-10 0.359
12 2021-06-26 0.432
13 2021-07-12 0.564
14 2021-07-28 0.654
15 2021-08-13 0.699
16 2021-08-29 0.614
17 2021-09-14 0.588
18 2021-09-30 0.553
19 2021-10-16 0.450
20 2021-11-01 0.377
21 2021-11-17 0.341
22 2021-12-03 0.331
23 2021-12-19 0.323

#I plot my dataset and fit a curve

p1 <- ggplot(data,aes(x = date, y = NDVI)) + stat_smooth(method = "lm", formula = y ~ ns(x,3), color="blue") + geom_point() 

p1

I now want to be able to calculate the maximum rate of curve change to identify when the vegetation is starting to greenup (my guess based on the figure is sometime in May).

Thanks for your help.

Answer 1

The maximum rate of increase in this model occurs at 9:13pm on 11th June 2021. A simple way to get this result is to create your model outside of ggplot:

library(splines)

mod <- lm(NDVI ~ ns(date, 3), data = data)

Now use this to predict NDVI for each day in the year:

dates <- seq(as.Date("2021-01-01"), as.Date("2021-12-31"), by = "day")

predictions <- predict(mod, newdata = list(date = dates))

The maximum rate of change occurs at the maximum difference between consecutive days' predictions:

dates[which.max(diff(predictions))]
# [1] "2021-06-11"

If you need the precision down to the hour and minute (which isn't justified based on the number of data points), you would do:

data$date <- as.POSIXct(data$date)
mod <- lm(NDVI ~ ns(date, 3), data = data)

dates <- seq(as.POSIXct("2021-01-01"), as.POSIXct("2021-12-31"), by = "min")

predictions <- predict(mod, newdata = list(date = dates))

dates[which.max(diff(predictions))]
#> [1] "2021-06-11 21:13:00 GMT"

Answer 2

One possible approach could be to build the ggplot, extract the spline, and then calculate on its control points:

library(dplyr); library(ggplot2)
a <- mtcars %>%
  ggplot(aes(drat, mpg)) +
  stat_smooth(method = "lm", formula = y ~ splines::ns(x,3), color="blue") +
  geom_point()

b <- ggplot_build(a)
spline_x <- b[["data"]][[1]][["x"]]
spline_y <- b[["data"]][[1]][["y"]]
change_rate <- c(NA, diff(spline_y)/ diff(spline_x))

a +
  geom_line(data = data.frame(spline_x, change_rate),
            aes(spline_x, change_rate), color = "red")

Then we could extract the value with max change_rate:

data.frame(spline_x, change_rate) %>%
  filter(change_rate == max(change_rate, na.rm = TRUE))

  spline_x change_rate
1 3.748861    10.68438