使用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()处理公式中的用户定义函数.我的问题是如何才能得到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