在 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 包装为用户定义的函数（实际上我会有自己的更复杂的函数），那么我会得到错误：

案例二：

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

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

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() 的调用已在 return 编辑的公式的 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() 函数，该函数是基于 returned 对象的 class 进行调度的通用函数。碰巧 poly return 是 class "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" 但 return 是一个 "poly" 对象，因此系数信息未被 return 编辑。一种解决方法是更改对象的 return class 并创建自己的 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

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

Use poly() in R formula to predict

r

function

formula

predict

案例一：

案例二：