The below quadratic model example (including data) is from Steve's Data Tips and Tricks.

# Sample Data and Quadratic Model
study_hours <- c(6, 9, 12, 14, 30, 35, 40, 47, 51, 55, 60)
exam_scores <- c(14, 28, 50, 70, 89, 94, 90, 75, 59, 44, 27)
data <- data.frame(study_hours, exam_scores)
model <- lm(exam_scores ~ study_hours + I(study_hours^2),
  data = data)

I can calculate the value of study_hours where the exam_score is maximized (inflection point or critical value) using a solution in Extract critical points of a polynomial model object in R?.

model %>% 
  mosaic::makeFun() %>%     
  optimize(interval = c(min(df$study_hours), max(df$study_hours)), maximum = TRUE) 

Another dataset and model (with error)

Now, I have another dataset and model below. I would like to get the value of variable ceoten where the log(salary) is maximized. However, I get the error below.

I am not sure if I correctly understood how the optimize function works.

(1) Can anyone share some insights?

(2) Would there be any other way to achieve the task?

library(wooldridge)
data("ceosal2")
ceosalary <- ceosal2
model2 <- lm(log(salary) ~ log(sales) + comten + 
     ceoten + I(ceoten^2),
                 data=ceosalary) 
    
library(mosaic)
model2 %>% 
# This function is used to extract the formula of the model as a function
mosaic::makeFun() %>%     
  optimize(interval = c(min(ceosalary$ceoten),
   max(ceosalary$ceoten)), maximum = TRUE) 

Error in f(arg, ...) : argument "comten" is missing, with no default

Your recipe doesn't work when there are predictor variables in the model other than the focal (quadratic) predictor.

The easiest way to do this is with some calculus. Suppose y = b0 + b1*x1 + b2*x2 + b3*x3 + b4*x3^2 and we want to find a critical point with respect to x3. Then:

dy/dx3 = b3 + 2*b4*x3

We set the LHS to zero and get 0 = b3 + 2*b4*x3 ⇒ x3 = -b3/(2*b4).

b3 <- coef(model2)[["ceoten"]]
b4 <- coef(model2)[["I(ceoten^2)"]]
maxval <- -b3/(2*b4)  ## 21.416

pframe <- with(ceosalary,
   data.frame(sales = mean(sales), comten = mean(comten),
              ceoten = seq(min(ceoten), max(ceoten), length.out = 51)))
plot(pframe$ceoten, predict(model2, newdata = pframe), type = "l")
abline(v = maxval, lty = 2, col = 2)

Alternatively, you can set up your function using predict:

pfun <- function(x) {
   pframe <- with(ceosalary,
        data.frame(sales = mean(sales), comten = mean(comten),
                   ceoten = x))
   predict(model2, newdata = pframe)
}
optimize(pfun, interval = c(min(ceosalary$ceoten), max(ceosalary$ceoten)), maximum = TRUE)

(this will work for cases where you can't/don't feel like doing the calculus)