R中的GLM功能,日志链接不起作用

Question

我正在使用Hardin和Hilbe的书"Generalized Linear Models and Extension"(第二版,2007).作者提出,"日志链接通常用于响应数据,而不是OLS模型,而这些响应数据只能在连续尺度上采用正值".当然,他们还建议使用残差图来检查是否仍然可以使用使用身份链接的"正常"线性模型.

我试图在R中复制他们在STATA中所做的事情.实际上,我在STATA中没有日志链接的问题.但是,当使用R的glm函数调用相同的模型时,却指定family=gaussian(link="log")我要求提供起始值.当我将它们全部设置为零时,我总是得到算法没有收敛的消息.选择其他值的消息有时是相同的,但我经常得到:

Error in glm.fit(x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  :
     NA/NaN/Inf in 'x'

正如我所说,在STATA中,我可以在不设置起始值且没有错误的情况下运行这些模型.我尝试了很多不同的模型和不同的数据集,但问题总是一样的(除非我只包含一个单独的自变量).谁能告诉我为什么会出现这种情况,或者我做错了什么,或者为什么书中建议的模型可能不合适？我很感激任何帮助,谢谢!

编辑:作为再现错误的示例,请考虑可以在此处下载的数据集.加载此数据集后,我运行以下模型:

mod <- glm(betaplasma ~ age + vituse, family=gaussian(link="log"), data=data2, start=c(0,0,0))

这会产生算法未收敛的警告消息.

Edit2:我被要求提供该模型的STATA输出.这里是:

. glm betaplasma age vituse, link(log)

Iteration 0:   log likelihood = -2162.1385  
Iteration 1:   log likelihood = -2096.4765  
Iteration 2:   log likelihood = -2076.2465  
Iteration 3:   log likelihood = -2076.2244  
Iteration 4:   log likelihood = -2076.2244  

Generalized linear models                          No. of obs      =       315
Optimization     : ML                              Residual df     =       312
                                                   Scale parameter =  31384.51
Deviance         =  9791967.359                    (1/df) Deviance =  31384.51
Pearson          =  9791967.359                    (1/df) Pearson  =  31384.51

Variance function: V(u) = 1                        [Gaussian]
Link function    : g(u) = ln(u)                    [Log]

                                                   AIC             =  13.20142
Log likelihood   = -2076.224437                    BIC             =   9790173

------------------------------------------------------------------------------
             |                 OIM
  betaplasma |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         age |   .0056809   .0032737     1.74   0.083    -.0007354    .0120972
      vituse |   -.273027   .0650773    -4.20   0.000    -.4005762   -.1454779
       _cons |   5.467577   .2131874    25.65   0.000     5.049738    5.885417
------------------------------------------------------------------------------

Answer 1

正如我在评论中所说的那样,Stata可能比R更具有鲁棒性(在数值上,而不是统计意义上)GLM拟合.也就是说,拟合这个特定的数据集似乎并不太难.

读取数据:

data2 <- read.table("http://lib.stat.cmu.edu/datasets/Plasma_Retinol",
         skip=30,nrows=315)
dnames <- c("age","sex","smokstat","quetelet","vituse","calories","fat","fiber",
           "alcohol","cholesterol","betadiet","retdiet","betaplasma","retplasma")
names(data2) <- dnames

绘制数据:

par(mfrow=c(1,2),las=1,bty="l")
with(data2,plot(betaplasma~age))
with(data2,boxplot(betaplasma~vituse))

通过将截距参数的起始值设置为合理的值(即接近对数刻度上数据平均值的东西:这些中的任何一个),可以很容易地使这些变得合适

mod <- glm(betaplasma ~ age + vituse, family=gaussian(link="log"), data=data2,
           start=c(10,0,0))
mod <- glm(betaplasma ~ age + vituse, family=gaussian(link="log"), data=data2,
           start=c(log(mean(data2$betaplasma)),0,0))

后一种情况可能是启动日志链接拟合的合理默认策略.结果(略有缩写)与Stata非常接近:

summary(mod)
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.467575   0.218360  25.039  < 2e-16 ***
## age          0.005681   0.003377   1.682   0.0935 .  
## vituse      -0.273027   0.065552  -4.165 4.03e-05 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
## 
## (Dispersion parameter for gaussian family taken to be 31385.26)
## 
##     Null deviance: 10515638  on 314  degrees of freedom
## Residual deviance:  9791967  on 312  degrees of freedom
## AIC: 4160.4
## 
## Number of Fisher Scoring iterations: 9

confint(mod)
##                     2.5 %      97.5 %
## (Intercept)  5.0364648709  5.87600710
## age         -0.0007913795  0.01211007
## vituse      -0.4075213916 -0.14995759

(对于p值和(？)置信区间,R使用t而不是Z统计量)

但是,有几个原因我可能不适合这些模型这些数据.特别是,恒定方差的假设(与高斯模型相关)不是很合理 - 这些数据似乎更适合对数正态模型(或等效地,仅用于对数变换和用标准高斯模型分析).

在log(1+x)比例上绘图(数据中有一个零条目):

with(data2,plot(log(1+betaplasma)~age))
with(data2,boxplot(log(1+betaplasma)~vituse))

绘图ggplot(这适用于每个值的单独行vituse而不是拟合附加模型)

library(ggplot)
theme_set(theme_bw())
(g1 <- qplot(age,1+betaplasma,colour=factor(vituse),data=data2)+
    geom_smooth(method="lm")+
    scale_y_log10())

没有'异常值'的视图:

g1 %+% subset(data2,betaplasma>0)

另外两点:(1)在这个数据集中有一个值为0的响应有点奇怪 - 不是不可能,而是奇数; (2)它看起来vituse应该被视为一个因素而不是数字("1 =是,经常,2 =是,不经常,3 =否") - 可能是序数.

Answer 2

我想建议可能是非正常的错误.如果您同意(或者更确切地说,如果数据同意),那么请考虑以下结构:

?family
?glm
?binomial
lfit <- glm( dep <- indep1 + indep2, data=dat, family=binomial(link="probit")

这应该提供身份标记模型周围的二项式错误.这样做的好处是您的估算值更容易在变量的原始范围内解释.对早期不正确的建议family = poisson与probit链接一起使用表示歉意.请记住,您从未提供任何数据甚至是分布的描述.显然二项式错误不适合@BenBolker提供的数据集.

如果您有非整数值与数正态分布的错误,你应该考虑quasipoisson模型.如果您运行对本Bolker提供的数据这一模式,比较高斯(链接="日志)MODELA他们几乎没有区别,不需要任何初始值.

> mod2 <- glm(betaplasma ~ age + vituse, family=quasipoisson, data=data2         )
> mod2

Call:  glm(formula = betaplasma ~ age + vituse, family = quasipoisson, 
    data = data2)

Coefficients:
(Intercept)          age       vituse  
   5.452014     0.006096    -0.276679  

Degrees of Freedom: 314 Total (i.e. Null);  312 Residual
Null Deviance:      37270 
Residual Deviance: 33420    AIC: NA 

> glm(betaplasma ~ age + vituse, family=gaussian(link="log"), data=data2,
+            start=c(10,0,0))

Call:  glm(formula = betaplasma ~ age + vituse, family = gaussian(link = "log"), 
    data = data2, start = c(10, 0, 0))

Coefficients:
(Intercept)          age       vituse  
   5.467575     0.005681    -0.273027  

Degrees of Freedom: 314 Total (i.e. Null);  312 Residual
Null Deviance:      10520000 
Residual Deviance: 9792000  AIC: 4160 

您应该使用稍微复杂的模型,因为vituse显然是一个三级因素:

> mod2 <- glm(betaplasma ~ age + factor(vituse), family=quasipoisson, data=data2         )
> mod2

Call:  glm(formula = betaplasma ~ age + factor(vituse), family = quasipoisson, 
    data = data2)

Coefficients:
    (Intercept)              age  factor(vituse)2  factor(vituse)3  
       5.151076         0.006359        -0.224107        -0.562727  

Degrees of Freedom: 314 Total (i.e. Null);  311 Residual
Null Deviance:      37270 
Residual Deviance: 33380    AIC: NA