r - 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具有比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))

通过将intercept参数的起始值设置为合理的值(即接近于对数刻度上的数据均值的值)，可以很容易地使这些参数适合:

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())

没有“异常值”的 View :

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

另外两点:(1)在此数据集中有一个值为0的响应，这有点奇怪-不是不可能，而是奇怪； (2)看来 vituse应该被视为一个因素而不是数字(“1 =是，经常是2 =是，不经常是，3 =否”)-可能是序数。