R中的GLM函数，日志链接不起作用

Question

我目前正在使用 Hardin 和 Hilbe 的“广义线性模型和扩展”一书（第二版，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 中，我可以运行这些模型而无需设置起始值且不会出现错误。我尝试了许多不同的模型和不同的数据集，但问题总是相同的（除非我只包含一个自变量）。谁能告诉我为什么会这样，或者我做错了什么，或者为什么书中建议的模型可能不合适？非常感谢您的帮助，谢谢！

编辑：作为重现错误的示例，请考虑可以下载的数据集here。加载此数据集后，我运行以下模型：

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 的 GLM 拟合可能比 R 更稳健（在数值上，而不是统计意义上）。也就是说，拟合这个 特定的数据集并没有好像太难了。

读取数据：

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

（R 对 p 值和 (?) 置信区间使用 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 与概率链接一起使用。请记住，您从未提供过任何数据甚至是对分布的描述。显然二项式错误不适用于@BenBolker 提供的数据集。

如果您有具有对数正态分布误差的非整数值，则应考虑准泊松模型。如果您在 Ben Bolker 提供的数据上运行此模型并比较 gaussian(link="log) 模型a，它们几乎无法区分，并且不需要起始值。

> 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