从 lm 对象计算变量的临界点

Question

下面的二次模型示例（包括数据）来自 Steve 的数据提示和技巧。

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

我可以使用

在 R 中提取多项式模型对象的临界点？

中的解决方案计算

study_hours

的值，其中 exam_score 最大化（拐点或临界值）。

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

另一个数据集和模型（有错误）

现在，我在下面有另一个数据集和模型。 我想获取变量

ceoten

的值，其中

log(salary)

最大化。但是，我得到以下错误。

我不确定我是否正确理解了

optimize

函数的工作原理。

(1) 谁能分享一些见解吗？

(2) 还有其他方法可以完成任务吗？

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) 

f(arg, ...) 中的错误：缺少参数“comten”，没有默认值

Answer 1

最简单的方法是使用一些微积分。假设

y = b0 + b1*x1 + b2*x2 + b3*x3 + b4*x3^2

并且我们想要找到关于

x3

的临界点。然后：

dy/dx3 = b3 + 2*b4*x3

我们将 LHS 设置为零并得到

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)

或者，您可以使用

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)

（这适用于您不能/不想进行微积分的情况）