0

按照 Hilbe 的 J. 2011 程序（在此称为“本书”）第 20 页，我没有得到相同的结果。该程序用于使用glm R. Hilbe 中的泰坦尼克数据集计算具有稳健标准误差的泊松回归根据以下链接，源代码在表 2.4 中： Negative Binomial Regression Second edition Errata 2012

我相信自从在此处发布以来， titanic数据集发生了一些变化，这是本书中的程序和Stata结果以及执行的操作，这导致截至 2021 年 7 月 R 中的结果不正确或不同：

library(COUNT)
data("titanic")
attach(titanic)
library(gmodels)

str(titanic)
'data.frame':   1316 obs. of  4 variables:
 $ class   : Factor w/ 3 levels "3rd class","1st class",..: 2 2 2 2 2 2 2 2 2 2 ...
 $ age     : Factor w/ 2 levels "child","adults": 2 2 2 2 2 2 2 2 2 2 ...
  ..- attr(*, "label")= chr "0=child; 1=adult"
  ..- attr(*, "format")= chr "%10.0g"
  ..- attr(*, "value.label.table")= Named int  0 1
  .. ..- attr(*, "names")= chr  "child" "adults"
 $ sex     : Factor w/ 2 levels "women","man": 2 2 2 2 2 2 2 2 2 2 ...
  ..- attr(*, "label")= chr "gender: 0=female; 1=male"
  ..- attr(*, "format")= chr "%8.0g"
  ..- attr(*, "value.label.table")= Named int  0 1
  .. ..- attr(*, "names")= chr  "women" "man"
 $ survived: num  2 2 2 2 2 2 2 2 2 2 ...
 - attr(*, "stata.info")=List of 5
  ..$ datalabel  : chr "Hilbe, Modeling Count Data (CUP, 2014)"
  ..$ version    : int 12
  ..$ time.stamp : chr "14 Jul 2014 15:12"
  ..$ val.labels : chr  "class" "age" "sex" "survived"
  ..$ label.table:List of 4
  .. ..$ class   : Named int  1 2 3 4
  .. .. ..- attr(*, "names")= chr  "1st class" "2nd class" "3rd class" "crew"
  .. ..$ age     : Named int  0 1
  .. .. ..- attr(*, "names")= chr  "child" "adults"
  .. ..$ sex     : Named int  0 1
  .. .. ..- attr(*, "names")= chr  "women" "man"
  .. ..$ survived: Named int  0 1
  .. .. ..- attr(*, "names")= chr  "no" "yes"

这本书重新调整了班级。

titanic$class <- relevel(factor(titanic$class), ref=3)

然而，截至 2021 年，“生存”已成为与我认为曾经是二进制 0="no" 和 1="yes" 整数相反的一个因素，因此，生存被相应地重新编码

titanic$survived <- as.character(titanic$survived)
titanic$survived [which(titanic$survived =="no")] <- "0"
titanic$survived [which(titanic$survived =="yes")] <- "1"
titanic$survived <- as.integer(titanic$survived)

2012 勘误表中的代码：

tit3 <- glm(survived ~ factor(class), family=poisson, data=titanic)
irr <- exp(coef(tit3)) # vector of IRRs
library("sandwich")
rse <- sqrt(diag(vcovHC(tit3, type="HC0"))) # coef robust SEs
irr*rse # IRR robust SEs

R 控制台中的 irr*rse 输出

 (Intercept) factor(class)1st class factor(class)2nd class 
            0.01634255             0.19270871             0.15723303

使用汇总功能

> summary(tit3)

Call:
glm(formula = survived ~ factor(class), family = poisson, data = titanic)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.1177  -0.7101  -0.7101   0.4364   1.1225  

Coefficients:
                       Estimate Std. Error z value Pr(>|z|)    
(Intercept)            -1.37783    0.07495 -18.384  < 2e-16 ***
factor(class)1st class  0.90721    0.10268   8.835  < 2e-16 ***
factor(class)2nd class  0.49603    0.11871   4.179 2.93e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 967.81  on 1315  degrees of freedom
Residual deviance: 889.69  on 1313  degrees of freedom
AIC: 1893.7

Number of Fisher Scoring iterations: 5

以下被认为是正确的估计，因为它是书中的内容。事故率 (IRR) 为：

class2: 1.642184 
class1: 2.477407

和估计和稳健的标准。呃。

class2: 0.1572928
class1: 0.192782

都有P>|z| == 0。有人可以证实这一点吗？谢谢

4

1 回答 1

0

确认的！

data('titanic', package="COUNT")
titanic <- transform(titanic, survived=as.numeric(survived) - 1,
                     class=relevel(class, ref=3))

tit3 <- glm(survived ~ class, family=poisson, data=titanic)

library(sandwich);library(lmtest)
(smy <- coeftest(tit3, vcovHC(tit3, type="HC0")))
# z test of coefficients:
#   
#                 Estimate Std. Error  z value  Pr(>|z|)    
# (Intercept)    -1.377832   0.064819 -21.2565 < 2.2e-16 ***
# class1st class  0.907212   0.077786  11.6629 < 2.2e-16 ***
# class2nd class  0.496027   0.095746   5.1806 2.211e-07 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(irr <- exp(coef(tit3)))
# (Intercept) class1st class class2nd class 
#   0.2521246      2.4774071      1.6421841 

rse <- sqrt(diag(vcovHC(tit3, type="HC0")))
irr*rse
# (Intercept) class1st class class2nd class 
#  0.01634255     0.19270871     0.15723303
于 2021-07-06T04:07:11.873 回答