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

问题描述

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

案例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

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

案例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的结果？

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

推荐答案

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"

因此，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

这篇关于在R公式中使用poly()进行预测的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！