R中glm()对predict.glm()的偏移量

Question

我感到困惑的是，在glm()中定义一个偏移量(应该是上述的)和在predict.glm()中使用不使用日志转换偏移量的newdata定义偏移量(显然)不一致。

我用lambda=2构造了一个数据集，每个数据集有5个“处理”，唯一的区别是计数的时间，即偏移量。

    pmean <- 2
    Offset <- rep(c(10,31.6,100,316,1000),5)

    set.seed(7919)
    data <- data.frame(cnt = rpois(length(Offset),pmean*Offset),
                       trt = factor(rep(c("A","B","C","D","E"),5)),
                       Offset = Offset)

我安装了两个glm模型，一个没有偏移量，另一个包括偏移量。因为链接是日志，所以当输入到模型时，必须记录偏移量变量。从glm()文档中:偏移量是“在拟合过程中包含在线性预测器中的一个先验已知分量”。

    noOffsetModel <- glm(cnt ~ trt, family = poisson(log), data = data)
    noOffsetModel

    # Call:  glm(formula = cnt ~ trt, family = poisson(log), data = data)
    # 
    # Coefficients:
    # (Intercept)         trtB         trtC         trtD         trtE  
    #       2.845        1.395        2.435        3.629        4.746  
    # 
    # Degrees of Freedom: 24 Total (i.e. Null);  20 Residual
    # Null Deviance:        20740 
    # Residual Deviance: 21.29  AIC: 209.3

    OffsetModel <- glm(cnt ~ trt, family = poisson(log), offset = log(Offset), data = data)
    OffsetModel

    # Call:  glm(formula = cnt ~ trt, family = poisson(log), data = data, 
    #     offset = log(Offset))
    # 
    # Coefficients:
    # (Intercept)         trtB         trtC         trtD         trtE  
    #      0.5423       0.2444       0.1327       0.1761       0.1409  
    # 
    # Degrees of Freedom: 24 Total (i.e. Null);  20 Residual
    # Null Deviance:        29.64 
    # Residual Deviance: 21.29  AIC: 209.3

我为预测数据的结构创建了两个数据框架，即对单个时间单位的预测，即偏移量= 1，以及偏移量被日志转换的另一个数据帧，即log(1) = 0。

    preddata  <- data.frame(trt = unique(data$trt), Offset = 1)
    preddata0 <- data.frame(trt = unique(data$trt), Offset = 0)

当我使用newdata = preddata为noOffset模型调用noOffset时，就会得到原始计数尺度上的预测，这是因为模型中没有偏移量。

    exp(predict.glm(noOffsetModel, newdata = preddata, type = "link"))
    # [1]   17.2   69.4  196.4  648.2 1980.2

当我为偏移模型调用predict.glm而不给它一个newdata=时，我得到了原始计数尺度上的预测，这也是预期的。

    unique(exp(predict.glm(OffsetModel, type = "link")))
    # [1]   17.2   69.4  196.4  648.2 1980.2

当我用newdata = preddata (即使用偏移量= 1)调用predict.glm时，会得到考虑偏移量的预测值--所有处理的值大约为2。这是预料不到的，因为先前数据中的偏移量没有被日志转换，即使原始调用glm中的偏移量是。

    exp(predict.glm(OffsetModel, newdata = preddata, type = "link"))
    # [1] 1.720000 2.196203 1.964000 2.051266 1.980200

当我在preddata0中对偏移量进行日志转换时，与在原始glm()中输入偏移量()的方式一致，我会得到垃圾。

    exp(predict.glm(OffsetModel, newdata = preddata0, type = "link"))
    # 1 2 3 4 5 
    # 0 0 0 0 0 

在predict.glm(newdata=)中要求在原始计数比例尺上设置偏移量似乎非常不一致(而且容易出错)。在glm()中，偏移量中的值必须对日志链接进行日志转换。

谢谢你澄清。

Answer 1

可能不知何故忘记了这个问题，但下面是您不需要提供日志的原因。

在predict.glm下面，调用了predict.lm，这些行与您的问题相关：

predict.lm
[...]
offset <- rep(0, nrow(X))
if (!is.null(off.num <- attr(tt, "offset"))) 
   for (i in off.num) offset <- offset + eval(attr(tt,"variables")[[i + 1]],newdata)
if (!is.null(object$call$offset)) 
offset <- offset + eval(object$call$offset, newdata)

我们看看你的目标：

OffsetModel$call$offset
log(Offset)

因此，当您执行predict.glm时，它会通过eval(object$call$offset, newdata)并添加日志(偏移量)来进行预测。例如，您可以尝试：

eval(OffsetModel$call$offset, preddata)
[1] 0 0 0 0 0

您可以阅读predict.lm的其余代码，偏移量是作为预测器添加的。

底线是，如果调用offet=log().，predict.glm也将在新数据中以相同的方式对待列