我正在将模型的 R 和 SAS 的 GLM 输出与 Gamma 分布进行比较。点估计是相同的,但它们对标准误差的估计不同,因此 p 值也不同。

有谁知道为什么?我想知道 R 和 SAS 是否使用不同的方法来估计标准误差?也许 MLE 与时刻方法?

R 示例代码

set.seed(2) 
test = data.table(y = rnorm(100, 1000, 100), x1 = rnorm(100, 50, 20), x2 = rgamma(100, 0.01)) 
model = summary(glm(formula = y ~ x1+x2 , family = Gamma(link = "log"), data = test)) 

使用此处生成的相同数据,我使用以下代码在 SAS 中运行模型:
proc genmod data= test_data; 
                model y =  x1 x2 /link= log dist= gamma; 
    run; 

R 的输出如下:
Call: 
glm(formula = y ~ x1 + x2, family = Gamma(link = "log"), data = test) 
 
Deviance Residuals:  
     Min        1Q    Median        3Q       Max   
-0.26213  -0.08456  -0.01033   0.08364   0.20878   
 
Coefficients: 
              Estimate Std. Error t value Pr(>|t|)     
(Intercept)  6.9210757  0.0324674 213.170   <2e-16 *** 
x1          -0.0003371  0.0005985  -0.563    0.575     
x2           0.0234097  0.0627251   0.373    0.710     
--- 
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
 
(Dispersion parameter for Gamma family taken to be 0.01376099) 
 
    Null deviance: 1.3498  on 99  degrees of freedom 
Residual deviance: 1.3436  on 97  degrees of freedom 
AIC: 1240.6 
 
Number of Fisher Scoring iterations: 4 

SAS 的输出:

请您参考如下方法:

默认情况下,R 以与 sas/genmod/model 选项 scale=pearson 相同的方式计算 scale=1/dispersion 参数。尺度参数的选择会影响 SE。请参阅此处的文档:
https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_genmod_sect022.htm

默认情况下,SAS/genmod 提供形状参数的 MLE。假设拟合的 gamma 模型存储在列表“fit”中。要在 R 中加载 MASS 库,然后键入

gamma.shape(fit) 

这给出了形状参数 alpha 的 MLE。
如果你输入
summary(fit, dispersion=1/gamma.shape(fit)$alpha)) 

汇总函数在计算 SE 时将使用 alpha 的 MLE,并且它们将完全匹配 SAS/genmod。

我将就此单独发表一篇文章。虽然summary.glm 给出了正确的SE(使用指定的离散值),但它没有打印正确的AIC 值(它不使用指定的离散,而是使用用Pearson 残差计算的值)。差异很小,但我称之为错误。


评论关闭
IT干货网

微信公众号号:IT虾米 (左侧二维码扫一扫)欢迎添加!