смещение в glm () и предиката.glm () в R

Question

Меня смущает (кажущаяся) несогласованность в определении смещения в glm (), которое должно быть преобразовано в лог (как указано выше), и вgnett.glm () с использованием newdata, в котором не используется преобразование в лог.

Я строю набор данных с 5 «обработками» каждый с лямбда = 2. Единственная разница заключается во времени, в течение которого производятся подсчеты, т. Е. Смещение.

    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, одна без смещения, другая с учетом смещения. Поскольку ссылка - log, переменная смещения должна регистрироваться при вводе в модель. Из документации glm (): смещение равно «1006 * априорному известному компоненту, который будет включен в линейный предиктор во время подгонки»

    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)

Когда я вызываю Foregit.glm для модели noOffset с использованием newdata = preddata, я получаю прогнозы по исходной шкале подсчета, которая ожидается, поскольку в модели не было смещения.

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

Когда я вызываю Foregit.glm для модели смещения, не передавая ей новые данные = Я получаю прогнозы по исходной шкале подсчета, что опять-таки следует ожидать.

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

Когда я вызываю Forext.glm с новыми данными = предданные (т. е. со смещением = 1) я получаю прогнозы, которые учитывают смещение - значения приблизительно равны 2 для всех обработок. Этого не ожидается, потому что смещение в предданных не было преобразовано в лог, даже если смещение в исходном вызове glm было.

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

Когда я преобразовываю в журнале свое смещение в предданных0, в соответствии с тем, как я Мне нужно ввести offset () в исходном glm (), я получаю мусор

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

Кажется очень непоследовательным (и склонным к ошибкам) ​​требовать смещения по исходной шкале счетчика вgnett.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)

Таким образом, когда вы выполняете Forex.Glm, он проходит через eval(object$call$offset, newdata) и добавляет журнал (смещение), чтобы сделать ваш прогноз. Например, вы можете попробовать:

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

Вы можете прочитать оставшуюся часть кода для Foregnet. Смещение добавляется как предиктор.

Нижняя строка, если вы вызвали offet = log () .., предикат.glm также будет обрабатывать ваш столбец в newdata