r - 为什么无论我如何将偏移量输入到模型中，预测都不会忽略 R 中泊松模型的偏移量？

Question

我正在 R 中工作，但一直在 Stata 中验证我的结果，并通过这样做观察到 R 中的预测并没有忽略我与泊松模型的偏移。让我解释一下:

我在 R 中拟合了以下模型 - 模拟超额死亡率而不是简单的死亡率(ExpDeaths 是基于总人口的每个受试者年龄、性别和时期的预期死亡人数，接下来显示的 Stata 代码中的 logExpDeaths 是只是 ExpDeaths 的自然对数):

model <- glm(Event ~ relevel( as.factor(Period), ref=2) + relevel( as.factor(AgeCat), ref="50-59") + relevel( as.factor(Sex), ref="Female") relevel( as.factor(AlcCombo), ref="0") + relevel( as.factor(ScoreSurv), ref="0") + relevel( as.factor(DrugCombo), ref="0"), offset = (log(ExpDeaths)), data=data, family = poisson)

并使用 Stata 验证结果:

poisson Event ib1.Period ib1.Age i.Sex ib1.AlcCombo ib0.ScoreSurv ib0.DrugCombo, 
   offset(logExpDeaths)

使用上述代码行在 R 和 Stata 中的模型结果完全相同。

但是，当我尝试从模型中获取每个主题的线性预测变量时:

在 R 中，使用代码 predict(model, type="link") 我得到前五个值: -3.812156 -2.472995 -2.499536 -2.299561 -2.217279

但是，当我在 Stata 中使用代码 predict lp, xb nooffset 时，我得到了前五个值: 0.6458265 0.8994361 0.8994361 0.8588267 1.338368

这些是我想在 R 中生成的值，但我意识到问题是 R 没有忽略偏移量，就像我在 Stata 中预测 lb, xb 时一样，即保持根据预期死亡人数进行偏移，我得到的值与在 R 中得到的值相同: -3.812156 -2.472995 -2.499536 -2.299561 -2.217279

glm 的 R 文档(请参阅 https://www.math.ucla.edu/~anderson/rw1001/library/base/html/glm.html )指出“由 offset 指定的偏移量将不会包含在 Predict.glm 的预测中，而由公式中的偏移项指定的偏移量将包含在预测中”，即如果我使用像我一样模型，应该忽略偏移量:

model <- glm(Event ~ relevel( as.factor(Period), ref=2) + relevel( as.factor(AgeCat), ref="50-59") + relevel( as.factor(Sex), ref="Female") + relevel( as.factor(AlcCombo), ref="0") + relevel( as.factor(ScoreSurv), ref="0") + relevel( as.factor(DrugCombo), ref="0"), offset = (log(ExpDeaths)), data=data, family = poisson)

与使用以下内容相反，这意味着根据文档使用预测时不会忽略偏移量:

model <- glm(Event ~ relevel( as.factor(Period), ref=2) + relevel( as.factor(AgeCat), ref="50-59") + relevel( as.factor(Sex), ref="Female") + relevel( as.factor(AlcCombo), ref="0") + relevel( as.factor(ScoreSurv), ref="0") + relevel( as.factor(DrugCombo), ref="0") + offset(log(ExpDeaths)), data=data, family = poisson)

但是，我使用两者得到了完全相同的模型(这是我所期望的)和线性预测器(应该是不同的)，这使我得出结论，在 R 中编写模型的两种方式都不会导致偏移量被忽略使用预测时。

我知道我可以使用 Stata 来获得所需的结果，但我真的想知道如何使用 R 来获得 Stata 结果，只是为了我自己的理智，即如何使用 R 进行预测以忽略偏移量。

Answer 1

当您调用nooffset时，您只需从线性预测器中减去偏移量。

Stata

use https://data.princeton.edu/wws509/datasets/ceb.dta,clear
gen y=round(mean*n,1)
gen os=log(n)
poisson y i.res, offset(os)
predict xb, xb
predict lp, xb nooffset
list in 1/6,clean


       i   dur     res            educ   mean    var    n    y         os         xb         lp  
  1.   1   0-4    Suva            None     .5   1.14    8    4   2.079442   3.284039   1.204598  
  2.   2   0-4    Suva   Lower primary   1.14    .73   21   24   3.044523    4.24912   1.204598  
  3.   3   0-4    Suva   Upper primary     .9    .67   42   38    3.73767   4.942267   1.204598  
  4.   4   0-4    Suva      Secondary+    .73    .48   51   37   3.931826   5.136423   1.204598  
  5.   5   0-4   Urban            None   1.17   1.06   12   14   2.484907   3.833794   1.348887  
  6.   6   0-4   Urban   Lower primary    .85   1.59   27   23   3.295837   4.644724   1.348887

R

在这里，请注意，我可以复制 stata 调用predict lp, xb nooffset，只需从 xb 中减去 os(请参阅 >ceb$lp=ceb$xb-ceb$os)

library(foreign)
ceb<- read.dta("http://data.princeton.edu/wws509/datasets/ceb.dta")
ceb$y  <- round(ceb$mean*ceb$n, 0)
ceb$os <- log(ceb$n)
m1 = glm(y~res, offset=os,data=ceb,family="poisson")
ceb$xb=predict(m1, type="link")
ceb$lp=ceb$xb-ceb$os 
head(ceb)

  i dur   res          educ mean  var  n  y       os       xb       lp
1 1 0-4  Suva          None 0.50 1.14  8  4 2.079442 3.284039 1.204598
2 2 0-4  Suva Lower primary 1.14 0.73 21 24 3.044522 4.249120 1.204598
3 3 0-4  Suva Upper primary 0.90 0.67 42 38 3.737670 4.942267 1.204598
4 4 0-4  Suva    Secondary+ 0.73 0.48 51 37 3.931826 5.136423 1.204598
5 5 0-4 Urban          None 1.17 1.06 12 14 2.484907 3.833794 1.348887
6 6 0-4 Urban Lower primary 0.85 1.59 27 23 3.295837 4.644724 1.348887