3

我已经拟合了一个具有方差结构的二次模型，该结构允许每个因子水平的不同方差水平，并且我无法预测仅包含 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联系人），但我想我也会接触到更广泛的社区。提前致谢

4

1 回答 1

0

感谢@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
于 2022-01-18T20:13:52.630 回答