从lm对象计算变量临界点：解决ceoten最大化log(salary)报错问题

解答 #

(1) 报错原因 #

你的model2包含多个自变量：log(sales)、comten、ceoten及ceoten的平方项。mosaic::makeFun()生成的函数需要所有自变量的取值才能计算预测值，但调用optimize时仅指定了ceoten的区间，未给log(sales)、comten传入固定值，函数无法确定这些变量的取值来计算log(salary)，因此抛出参数缺失的报错。

而第一个示例的模型仅包含study_hours一个自变量（含平方项），生成的函数只需这一个变量的值，因此optimize可正常运行。

(2) 替代方法 #

方法1：固定其他自变量取值（如均值）后求最优ceoten

先计算其他自变量的均值，在生成的函数中固定这些值后再调用optimize：

library(mosaic)

# 计算其他自变量的均值
mean_log_sales <- mean(log(ceosalary$sales), na.rm = TRUE)
mean_comten <- mean(ceosalary$comten, na.rm = TRUE)

# 生成模型函数并固定其他参数
model_fun <- mosaic::makeFun(model2)
optimize(
  f = function(x) model_fun(log(sales) = mean_log_sales, comten = mean_comten, ceoten = x),
  interval = range(ceosalary$ceoten, na.rm = TRUE),
  maximum = TRUE
)

方法2：利用二次函数解析解直接计算

模型中ceoten的项是二次形式（β1*ceoten + β2*ceoten²），对于二次函数y = a*x² + b*x + c，极值点为x = -b/(2a)，直接提取系数计算即可：

# 提取模型中ceoten相关项的系数
coefs <- coef(model2)
b <- coefs["ceoten"]
a <- coefs["I(ceoten^2)"]

# 计算最优ceoten值
optimal_ceoten <- -b/(2*a)
optimal_ceoten

注：该结果对应其他自变量取均值时的最优值，和方法1结果一致；若需其他自变量取特定值，优先选择方法1。

内容的提问来源于stack exchange，提问作者J.K.