在 R 公式中使用 poly() 进行预测

Question

我有一个关于公式和用户定义函数的问题：

案例一：

 clotting <- data.frame(
     u = c(5,10,15,20,30,40,60,80,100),
     lot1 = c(118,58,42,35,27,25,21,19,18),
     lot2 = c(69,35,26,21,18,16,13,12,12))

 g1 = glm(lot1 ~ log(u) + poly(u,1), data = clotting, family = Gamma)
 dc = clotting
 dc$u = 1
 predict(g1, dc)

      1           2           3           4           5           6           7           8           9
 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929

但是，如果我只是简单地将 poly 包装为用户定义的函数（实际上我会有自己的更复杂的函数），那么我会得到错误：

案例2：

 xpoly <- function(x, degree=1){poly(x,degree)}
 g2 = glm(lot1 ~ log(u) + xpoly(u,1), data = clotting, family = Gamma)
 predict(g2, dc)
       Error in poly(x, degree) :
      'degree' must be less than number of unique points

似乎预测用 I() 处理公式中的用户定义函数。我的问题是如何使 Case2 的结果与 case1 相同？

任何人都可以对此有所了解吗？

Answer 1

poly 在这里有点独特的功能。默认情况下，它返回一组正交多项式，因此它正在对数据进行一些居中和重新缩放。如果您希望能够使用拟合模型中的系数进行预测，则需要以与原始数据相同的方式转换新数据。这意味着必须传递一些额外的数据。

首先我要指出，如果您使用原始的非正交值，您不会遇到这个问题。

g1 <- glm(lot1 ~ log(u) + poly(u,1, raw=T), data = clotting, family = Gamma)
xpoly<-function(x,degree=1){poly(x,degree, raw=T)}
g2 <- glm(lot1 ~ log(u) + xpoly(u,1), data = clotting, family = Gamma)

dc=clotting
dc$u=1
predict(g1,dc)
#       1           2           3           4           5           6           7           8           9 
#-0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 
predict(g2,dc)
#       1           2           3           4           5           6           7           8           9 
#-0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929

但让我们进一步探讨poly 如何将缩放信息传递给predict。这项工作实际上发生在model.frame 函数中。比较这两个结果

attr(terms(model.frame(lot1 ~ log(u) + poly(u,1), clotting)), "predvar")
# list(lot1, log(u), poly(u, 1, coefs = list(alpha = 40, norm2 = c(1, 
9, 8850))))
attr(terms(model.frame(lot1 ~ log(u) + xpoly(u,1), clotting)), "predvar")
# list(lot1, log(u), xpoly(u, 1))

您可以看到第一个公式中对poly() 的调用已在返回的公式的predvar 属性中进行了调整。这是在model.frame 代码中完成的

...
if (is.null(attr(formula, "predvars"))) {
    for (i in seq_along(varnames)) predvars[[i + 1L]] <- makepredictcall(variables[[i]], 
        vars[[i + 1L]])
    attr(formula, "predvars") <- predvars
}
...

请注意，它调用makepredictcall() 函数，这是一个通用函数，它根据返回对象的类进行调度。碰巧poly 返回了一个“poly”类的对象

class(poly(1:5, 1))
# [1] "poly"   "matrix"

所以“poly”数据调用的就是这个函数

stats:::makepredictcall.poly
function (var, call) 
{
    if (as.character(call)[1L] != "poly") 
        return(call)
    call$coefs <- attr(var, "coefs")
    call
}
<bytecode: 0x123262178>
<environment: namespace:stats>

这是添加coef= 属性的地方。但还要注意，它会检查调用是否来自“poly”函数本身。由于您的函数名为“xpoly”但返回一个“poly”对象，因此不会返回系数信息。一种解决方法是更改​​对象的返回类并创建自己的makepredictcall 函数。例如你可以做

xpoly <- function(...){p<-poly(...); class(p)[1]<-"xpoly"; p}
makepredictcall.xpoly <- function(var, call) {
    call$coefs <- attr(var, "coefs")
    call
}

请注意，这个新版本的xpoly 还将接受coef= 参数并通过... 参数将其传递给poly()。然后就可以运行了

g1 <- glm(lot1 ~ log(u) + poly(u,1), data = clotting, family = Gamma)
g2 <- glm(lot1 ~ log(u) + xpoly(u,1), data = clotting, family = Gamma)
predict(g1,dc)
#          1           2           3           4           5           6           7           8           9 
#-0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929
predict(g2,dc)
#          1           2           3           4           5           6           7           8           9 
#-0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929