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

Question

我有一个关于公式和用户定义函数的问题:
情况1:

 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() 处理公式中的用户定义函数。我的问题是如何获得与案例 1 相同的案例 2 的结果？
任何人都可以对此有任何想法吗？

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() function 是一个通用函数，它根据返回对象的类进行调度。碰巧 poly返回类“poly”的对象

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

所以 get 对“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