轮船事故的计数数据模型

摘要
数据集
模型方法
- Poisson模型

摘要

计数变量主要用于描述某一事件发生的次数，它仅取整数值。例如，每户家庭的子女数、某人在一年中看医生的次数、两国之间的贸易纠纷案件诉讼数等。因变量为计数变量的模型称为计数模型（Count Model）。

模型相关数据是4个制造期、2个服务期的五种类型的轮船的事故发生次数的数据。TYPE表示轮船类型；TA、TB、TC、TD和TE是表示轮船类型的虚拟变量；T6064、T6569、T7074和T7579是制造期间虚拟变量；O6074和O7579是运营期间的虚拟变量；Mon是服务量的测量；Acc是发生事故的次数。由于在1975-1979年制造的轮船不可能在1960-1964年期间运营，因此每一类型的轮船都有一个缺失数据。此外，由于原数据引用处未交待的原因，第五类型的轮船还有一个缺失数据，因此实际样本容量只有34个。

我们感兴趣的是轮船发生事故的次数以及其影响因素，轮船发生事故的次数是正整数，不同于普通回归的因变量是连续的实数，该如何分析？？？

当因变量为计数数据时，一般是以泊松分布而非正态分布来描述它的概率的。

计数数据模型的估计：对计数模型公式进行参数估计，虽然可以对对数方程直接进行OLS估计，但是会丢失那些因变量取值为0的数据，因为0是无法取对数的。另一种选择是直接采用非线性的迭代方法，但往往无法克服模型中的异方差特征，因此更多的是采用最大似然法进行估计。

R中泊松回归的估计也用glm()函数，只是对应的family改为poission()即可。

\(Ｅ(ACC|x)\)=\(exp({\beta}_1 + {\beta}_2 TB + {\beta}_3 TC + {\beta}_4 TD + {\beta}_5 TE + {\beta}_6 T6569 + {\beta}_7T7074 + {\beta}_8T7579 + {\beta}_9O7579 + {\beta}_{10}ln{(Mon)}) + \xi\)

其中，TA、TB、TC、TD和TE是表示轮船类型的虚拟变量；T6064、T6569、T7074和T7579是制造期间的虚拟变量；O6064和O7579是运营期间的虚拟变量；Mon是服务量的测量；Acc是发生事故的次数。

数据集

mdata <- read.csv('http://data.galaxystatistics.com/blog_data/regression/data_transform_ship_accident.csv',header=T)

modeldata <- mdata
dim(modeldata)

## [1] 34 10

head(modeldata)

##   Acc TB TC TD TE T6569 T7074 T7579 O7579 log.Mon.
## 1   0  0  0  0  0     0     0     0     0 4.844187
## 2   0  0  0  0  0     0     0     0     1 4.143135
## 3   3  0  0  0  0     1     0     0     0 6.998510
## 4   4  0  0  0  0     1     0     0     1 6.998510
## 5   6  0  0  0  0     0     1     0     0 7.321189
## 6  18  0  0  0  0     0     1     0     1 8.117611

class(modeldata)

## [1] "data.frame"

typeof(modeldata)

## [1] "list"

mode(modeldata)

## [1] "list"

#数据集字段名称转换
renames <- function(x) {
  ocol <- ncol(x)
  if (ocol==2) {
    xnam <- paste0("X")
    names(x) <- c('Y', xnam)
  }else if(ocol>2) {
    xnam <- paste0("X", 1:(ocol-1))
    names(x) <- c('Y', xnam)
  }
  return(x)
}
modeldata <- renames(modeldata)
head(modeldata)

##    Y X1 X2 X3 X4 X5 X6 X7 X8       X9
## 1  0  0  0  0  0  0  0  0  0 4.844187
## 2  0  0  0  0  0  0  0  0  1 4.143135
## 3  3  0  0  0  0  1  0  0  0 6.998510
## 4  4  0  0  0  0  1  0  0  1 6.998510
## 5  6  0  0  0  0  0  1  0  0 7.321189
## 6 18  0  0  0  0  0  1  0  1 8.117611

class(modeldata)

## [1] "data.frame"

dim(modeldata)

## [1] 34 10

regFormula <- as.formula(paste("Y ~ ", paste(names(modeldata)[-1], collapse= "+")))
regFormula

## Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + X9

模型方法

Poisson模型

#Poisson模型
Poisson.fit <- glm(regFormula, family=poisson(), data=modeldata)
#展示拟合模型的详细结果
summary(Poisson.fit)

## 
## Call:
## glm(formula = regFormula, family = poisson(), data = modeldata)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6569  -0.8872  -0.4843   0.4789   2.7436  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -5.6185     0.8737  -6.431 1.27e-10 ***
## X1           -0.3586     0.2697  -1.330  0.18360    
## X2           -0.7621     0.3383  -2.253  0.02426 *  
## X3           -0.1314     0.2970  -0.442  0.65820    
## X4            0.2697     0.2419   1.115  0.26492    
## X5            0.6618     0.1539   4.301 1.70e-05 ***
## X6            0.7602     0.1781   4.268 1.97e-05 ***
## X7            0.3615     0.2473   1.462  0.14379    
## X8            0.3699     0.1182   3.129  0.00175 ** 
## X9            0.9062     0.1017   8.906  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 614.539  on 33  degrees of freedom
## Residual deviance:  38.132  on 24  degrees of freedom
## AIC: 156
## 
## Number of Fisher Scoring iterations: 5

#列出拟合模型的模型参数
coefficients(Poisson.fit)

## (Intercept)          X1          X2          X3          X4          X5 
##  -5.6185211  -0.3585995  -0.7621295  -0.1314032   0.2696666   0.6618199 
##          X6          X7          X8          X9 
##   0.7602049   0.3614621   0.3698627   0.9061702

# 预测
predict(Poisson.fit, modeldata[-1])

##           1           2           3           4           5           6 
## -1.22886322 -1.49427324  1.38513957  1.75500229  1.77592657  2.86748334 
##           7           8           9          10          11          12 
##  2.10482667  3.72958596  3.22905070  3.98334345  4.04541812  2.81425927 
##          13          14          15          16          17          18 
##  3.74370421  2.79215263  0.02816723 -0.28963890  0.31678393  0.55581192 
##          19          20          21          22          23          24 
##  0.41748650  1.61325609 -0.56287644 -0.74292364 -1.16278150  0.04350153 
##          25          26          27          28          29          30 
##  0.04594387  0.31597219  1.81098063  1.89192930 -1.16804567  1.35781490 
##          31          32          33          34 
##  1.19228227  1.80309943  2.73908348  1.08705261

#在泊松回归中，因变量以条件均值的对数形式来建模。例如，log(Mon)增加一个单位时，轮船事故发生次数的对数平均值将增加0.9062。截距项表示当所有自变量都为0时，轮船事故发生次数的对数均值。通常在因变量的初始尺度（轮船事故发生数而非发生数的对数）上解释回归系数比较容易。因此，经常需要对回归系数进行指数化。
#首先利用coef()提取回归模型的系数。
#回归系数指数化
round(exp(coefficients(Poisson.fit)),3)

## (Intercept)          X1          X2          X3          X4          X5 
##       0.004       0.699       0.467       0.877       1.310       1.938 
##          X6          X7          X8          X9 
##       2.139       1.435       1.448       2.475

轮船事故的计数数据模型

baorui

摘要

数据集

模型方法

Poisson模型