Extract critical points of a polynomial model object in R?

Question

I am trying to solve for the inflection points of a cubic polynomial function which has been fitted to data, i.e. values of x where the first derivative is zero.

I also need a way to find the values of y at the critical points of x.

It is easy enough to fit the model using lm() and to view the model quality with summary(). And I can plot the function easily enough by adding predictions and using geom_line().

There must be a package or a base R function dedicated to this problem. Can anyone suggest a method?

Below is a reprex to depict the problem. Needless to say, the arrows are drawn only to illustrate the question; they are not mapped to the true inflection points or I would not be asking this question...

library(tidyverse)
library(modelr)

set.seed(0)
#generate random data and plot the values
df <- tibble(x= sample(x= c(-100:200), size= 50),
             y= -0.5*(x^3) + 50*(x^2) + 7*(x)  + rnorm(n=50, mean=10000, sd=50000) )

df %>% ggplot(aes(x, y)) +
  geom_point()

# fit a model to the data
cubic_poly_model <- lm(data= df, formula = y~poly(x, 3))

# plot the fitted model
df %>%
  add_predictions(model = cubic_poly_model) %>%
  ggplot(aes(x, y))+
  geom_point(alpha=1/3)+
  geom_line(aes(x, y=pred))+
  annotate('text', label= 'critical point A', x=-50, y=-250000)+
  geom_segment(x=-50, xend=-10, y=-200000, yend=-5000, arrow = arrow(length=unit(3, 'mm'), type = 'closed'))+
  annotate('text', label= 'critical point B', x=140, y=400000)+
  geom_segment(x=110, xend=90, y=300000, yend=100000, arrow = arrow(length=unit(3, 'mm'), type = 'closed'))

# But how can I get the critical values of x and the y values they produce?

^{Created on 2020-09-03 by the reprex package (v0.3.0)}

Answer 1

I devised a solution using the mosaic package . The makeFun() function allows a model object to be converted to a function. You can then use base R optimize()to find the max or min value of that function over a specified interval (in this case, the range of x values). Specify the "maximum" argument in optimize() to state whether you want the local maximum or local minimum.

See code below:

library(magrittr) 
 set.seed(0)

    #generate random data and plot the values
  
  df <- tibble::tibble(x= sample(x= c(-100:200), size= 50),
               y= -0.5*(x^3) + 50*(x^2) + 7*(x)  + rnorm(n=50, mean=10000, sd=50000) )
  
  
  
  cubic_poly_model <- lm(data= df, formula = y~poly(x, 3))
  
  crit_values <- cubic_poly_model %>% 
    mosaic::makeFun() %>% 
    optimize(interval = c(min(df$x), max(df$x)), maximum = TRUE) 
  
  funct_crit_x <- crit_values[['maximum']][[1]]
  funct_max <- crit_values[['objective']]
  
  funct_crit_x
  funct_max