r语言 - 找到能使结果最大化的模型预测值

Question

如何为模型预测器(线性和非线性的混合)找到一组值，使其产生响应的最高值。

示例模型:

library(lme4); library(splines)
summary(lmer(formula = Solar.R ~ 1 + bs(Ozone) + Wind + Temp + (1 | Month), data = airquality, REML = F))

这里我感兴趣的是什么条件(预测因子)产生最高的太阳辐射(结果)。

这个问题看起来很简单，但是我用谷歌没能找到一个好的答案。

如果模型很简单，我可以求导数来求最大值或最小值。有人建议，如果可以提取模型函数，可以使用stats::optim()函数。作为最后的手段，我可以模拟输入值的所有合理变化，并将其插入predict()函数并寻找最大值。

提到的最后一种方法似乎不是很有效，我想这是一个足够常见的任务(例如，为广告寻找最佳客户)，有人已经构建了一些工具来处理它。

Answer 1

这里有一些概念上的问题。

• 对于简单项(Wind和Temp)，响应是预测因子的线性函数(因此是单调的和无界的)。因此，如果这些项有正的参数估计，将它们的值增加到无穷大(Inf)将得到一个无穷大的响应(Solar.R);如果系数为负，则值应尽可能小(负无穷)。实际上，如果参数估计分别为负或正，那么您希望将这些预测器设置为最小或最大合理的值。

• 对于bs项，我不确定b样条在边界结之外的性质，但我很确定曲线会走向正无穷大或负无穷大，所以你有同样的问题。然而，对于bs，也有可能存在一个或多个内部最大值。对于这种情况，我可能会尝试提取基项并在数据范围内评估样条。

或者，你提到optim让我认为这是一种可能性:

data(airquality)
library(lme4)
library(splines)
m1 <- lmer(formula = Solar.R ~ 1 + bs(Ozone) + Wind + Temp + (1 | Month),
           data = airquality, REML = FALSE)
predval <- function(x) {
    newdata <- data.frame(Ozone=x[1],Wind=x[2],Temp=x[3])
    ## return population-averaged prediction (no Month effect)
    return(predict(m1, newdata=newdata, re.form=~0))
}
aq <- na.omit(airquality)
sval <- with(aq,c(mean(Ozone),mean(Wind),mean(Temp)))
predval(sval)
opt1 <- optim(fn=predval,
      par=sval,
      lower=with(aq,c(min(Ozone),min(Wind),min(Temp))),
      upper=with(aq,c(max(Ozone),max(Wind),max(Temp))),
      method="L-BFGS-B",  ## for constrained opt.
      control=list(fnscale=-1))  ## for maximization
## opt1
## $par
## [1] 70.33851 20.70000 97.00000
## 
## $value
## [1] 282.9784

正如预期的那样，这是臭氧(1-168)范围内的中间值，风(2.3-20.7)和温度(57-97)的最小/最大值。

通过自动选择简单项的最小/最大值，并只对复杂项(多项式/样条等)进行优化，可以使这种蛮力解决方案更有效。