使用`poly(...,2)`拟合二次模型后，使用`nlme::gls`进行预测时出现的问题。

Question

我已经拟合了一个二次模型，并使用方差结构允许每个因素级别有不同的方差水平，但现在只有2个条目的新数据集上预测时遇到了困难。以下是一个可复制的示例：

library(nlme)
set.seed(101)
mtcars$amf <- factor(mtcars$am)
modGLS <- gls(mpg ~ amf*poly(hp, 2),
           weights = varIdent(form = ~ 1|amf), data = mtcars)
minhp <- min(mtcars$hp); maxhp <- max(mtcars$hp)
newdata <- data.frame(amf = as.factor(c(0, 1)),
                        hp = round(runif(2, min = minhp, max = maxhp)))
newdata2 <- data.frame(amf = as.factor(c(0, 0, 1)),
                        hp = round(runif(3, min = minhp, max = maxhp)))
predict(modGLS, newdata = newdata)
# Error in poly(hp, 2) : 'degree' must be less than number of unique points

predict(modGLS, newdata = newdata2)
## [1]  5.973306 13.758955 44.037921
## attr(,"label")
## [1] "Predicted values"

然而，在一个 lm 框架中，预测效果很好：

modLM <- lm(mpg ~ amf*poly(hp, 2), mtcars)
predict(modLM, newdata = newdata)
##        1        2 
## 25.22253 16.83943 

为什么会这样？ emmeans 的一个包维护者似乎认为这可能与attr(, “predvars”)上缺少的信息有关（请参见我们在此处的讨论：https://github.com/rvlenth/emmeans/issues/133）。

我已向Bates博士（nlme的联络点）报告了此事，但我想联系更广泛的社区。提前感谢您的帮助。

- Ceres

---

Doug Bates 不再是 nlme 的维护者; 现在是 R 核心团队。我会看一下这个问题，但它可能值得作为一个缺陷报告（如果还没有的话！） - Ben Bolker

顺便说一下，我怀疑即使你没有收到错误提示，当提供新数据时，预测结果也是错误的（这甚至更糟）。 - Ben Bolker

特别地，在terms(modGLS)中应该有一个“predvars”属性。 - Russ Lenth

刚刚添加了这个：https://bugs.r-project.org/show_bug.cgi?id=18283 - Ben Bolker

你的例子可以大大简化... - Ben Bolker

感谢@BenBolker和@russ-lenth的帮助，并向R-core报告了错误。我简化了示例并提供了一个实际的答案，探讨了GLS模型对象“terms”中缺失的属性。 - Ceres

Answer 1

感谢 @BenBolker 和 @russ-lenth 确认问题与 GLS 对象中缺失的 terms 属性 "predvars" 相关，该属性为 poly 提供了拟合系数。请注意，这在 LM 框架中是有效的（原始帖子），并且该属性存在（另请参见 ?makepredictcall）。请注意，这可能会对预测产生潜在影响。

## e.g. this fails
poly(newdata$hp, 2)
## but this is okay, because the polynomial has been estimated
polyFit <- poly(mtcars$hp, 2)
predict(polyFit, newdata = newdata$hp)

attr(terms(modLM), "predvars")
## list(mpg, amf, poly(hp, 2, coefs = list(alpha = c(146.6875, 198.071514938048), norm2 = c(1, 32, 145726.875, 977168457.086594))))

attr(terms(modGLS), "predvars")
## NULL

- Ceres

---

请注意，自2022年3月中旬以来，“nlme”上的此错误已经修复。 - Ceres