将R中glm函数计算的标准误差估计值与SAS PROC GENMOD中计算的进行比较

Question

我正在将一个SAS PROC GENMOD的示例转换为R，使用R中的glm。SAS代码如下：

proc genmod data=data0 namelen=30;
model boxcoxy=boxcoxxy ~ AGEGRP4 + AGEGRP5 + AGEGRP6 + AGEGRP7 + AGEGRP8 + RACE1 + RACE3 + WEEKEND + 
SEQ/dist=normal;
FREQ REPLICATE_VAR;  
run;

我的R代码是：

parmsg2 <- glm(boxcoxxy ~ AGEGRP4 + AGEGRP5 + AGEGRP6 + AGEGRP7 + AGEGRP8 + RACE1 + RACE3 + WEEKEND + 
SEQ , data=data0, family=gaussian, weights = REPLICATE_VAR)

当我使用summary(parmsg2)时，我得到的系数估计和SAS中的相同，但我的标准误差有很大的差异。

来自SAS的总结输出为：

Name         df   Estimate      StdErr    LowerWaldCL  UpperWaldCL      ChiSq   ProbChiSq
Intercept    1   6.5007436    .00078884      6.4991975    6.5022897    67911982 0
agegrp4      1   .64607262    .00105425      .64400633    .64813891   375556.79 0
agegrp5      1    .4191395    .00089722      .41738099    .42089802   218233.76 0
agegrp6      1  -.22518765    .00083118     -.22681672   -.22355857   73401.113 0
agegrp7      1  -1.7445189    .00087569     -1.7462352   -1.7428026   3968762.2 0
agegrp8      1  -2.2908855    .00109766     -2.2930369   -2.2887342   4355849.4 0
race1        1  -.13454883    .00080672     -.13612997   -.13296769    27817.29 0
race3        1  -.20607036    .00070966     -.20746127   -.20467944   84319.131 0
weekend      1    .0327884    .00044731       .0319117    .03366511   5373.1931 0
seq2          1 -.47509583    .00047337     -.47602363   -.47416804   1007291.3 0
Scale         1 2.9328613     .00015586      2.9325559    2.9331668     -127

R的摘要输出如下：

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.50074    0.10354  62.785  < 2e-16 
AGEGRP4      0.64607    0.13838   4.669 3.07e-06 
AGEGRP5      0.41914    0.11776   3.559 0.000374 
AGEGRP6     -0.22519    0.10910  -2.064 0.039031  
AGEGRP7     -1.74452    0.11494 -15.178  < 2e-16
AGEGRP8     -2.29089    0.14407 -15.901  < 2e-16
RACE1       -0.13455    0.10589  -1.271 0.203865    
RACE3       -0.20607    0.09315  -2.212 0.026967 
WEEKEND      0.03279    0.05871   0.558 0.576535 
SEQ         -0.47510    0.06213  -7.646 2.25e-14

重要的区别在于标准误差，SAS系数都是统计显著的，但是R输出中的RACE1和WEEKEND系数不是。我已经找到了一个公式来计算R中的Wald置信区间，但由于标准误差的差异，这是无意义的，因为我不会得到相同的结果。

显然，SAS使用稳定Ridge-Newton-Raphson算法进行其估计，这些估计是ML。我读到的关于R中glm函数的信息是结果应该等价于ML。我该如何更改我的R估计过程，以便获得与SAS产生的等效系数和标准误差估计值？

更新一下，感谢Spacedman的回答，我使用权重，因为数据来自膳食调查中的个人，REPLICATE_VAR是平衡重复复制权重，是一个整数（非常大，约为1000或10000）。描述该权重的网站在此处。我不知道为什么SAS使用FREQ而不是WEIGHT命令。现在，我将通过使用REPLICATE_VAR扩展观测数量并重新运行分析来进行测试。

感谢Ben以下的答案，我现在正在使用的代码是：

parmsg2 <- coef(summary(glm(boxcoxxy ~ AGEGRP4 + AGEGRP5 + AGEGRP6 + AGEGRP7 + AGEGRP8 + RACE1 + RACE3 
+ WEEKEND + SEQ , data=data0, family=gaussian, weights = REPLICATE_VAR)))
#clean up the standard errors
parmsg2[,"Std. Error"] <- parmsg2[,"Std. Error"]/sqrt(mean(data0$REPLICATE_VAR)) 
parmsg2[,"t value"] <- parmsg2[,"Estimate"]/parmsg2[,"Std. Error"] 
#note: using the t-distribution for p-values, correct the t-values
allsummary <- summary.glm(glm(boxcoxxy ~ AGEGRP4 + AGEGRP5 + AGEGRP6 + AGEGRP7 + AGEGRP8 + RACE1 +
RACE3 + WEEKEND + SEQ , data=data0, family=gaussian, weights = REPLICATE_VAR))
parmsg2[,"Pr(>|t|)"] <- 2*pt(-abs(parmsg2[,"t value"]),df=allsummary$df.resid)

- Michelle

---

你的R标准误差大约是SAS标准误差的131.25倍，如果这有任何启示的话。为了更好地理解，我会将问题简化为一个变量并观察结果。 - Spacedman

Answer 1

SAS中的FREQ与R中glm中的weights不同。在SAS中，它是事件发生的次数。对于R来说，则是“每个响应y_i是w_i个单位权观测值的平均值”。这两者并不相同。

如果您希望R给出与SAS相同的输出（想不出为什么），那么您可能需要将数据框中的每一行重复“weight”次。

这里，“data”是所有权重均为2的10行数据，“data2”是所有权重均为1的20行数据（data的每一行都有2个副本）：

> summary(glm(y~x,data=data2,weights=weights))$coef
              Estimate Std. Error   t value   Pr(>|t|)
(Intercept) 0.32859847 0.13413683 2.4497259 0.02475748
x           0.01540002 0.02161811 0.7123667 0.48537003
> summary(glm(y~x,data=data,weights=weights))$coef
              Estimate Std. Error   t value  Pr(>|t|)
(Intercept) 0.32859847 0.20120525 1.6331506 0.1410799
x           0.01540002 0.03242716 0.4749111 0.6475449

简单来说，同一个值的N次观测比将该观测的平均值作为N次观测更加精确，因此重复观测的标准误差会比平均值的标准误差小。

- Spacedman

---

我尝试使用vcdExtra包中的expand.dft()函数，除了在glm中使用的对象外，已删除所有对象，但是R现在崩溃了。我正在使用64位3.14.0版本，在Windows 7 Vmware链接中。行数将为177,050,435，取自REPLICATE_VAR的总和。显然，这将是美国人口中年龄在9岁及以上的男性数量。 - Michelle

Answer 2

编辑：阅读SAS文档中关于FREQ的内容以及您在上下方的回复后，这是我认为您应该尝试的方法：在glm语句中使用weights=REPLICATE_VAR来调整组别的相对权重（您发现的系数相等表明这是正确的方法），然后使用下面建议的N=sum(REPLICATE_VAR)进行调整（我还认为您可以使用lm而不是glm来解决这个问题...虽然差别不大，但速度会更快，更健壮。） 类似这样：

s <- coef(summary(lm(y~x,data=data2, weights=REPLICATE_VAR)))
s[,"Std. Error"] <- s[,"Std. Error"]/sqrt(sum(data2$REPLICATE_VAR))
s[,"t value"] <- s[,"Estimate"]/s[,"Std. Error"]
s[,"Pr(>|t|)"] <- 2*pt(abs(s[,"t value"]),df=g$df.resid)

- Ben Bolker

---

我们现在非常接近，SE估计相差100倍（R估计现在小100倍）。如果它们是相同数量级的话，差异将在第5位小数处，这我不会关心。 - Michelle

是恰好可以被100整除的因子？还是大约可以被100整除的因子？ - Ben Bolker

有点随意尝试的感觉，但是除以 sqrt(data2$REPLICATE_VAR)（而不是 sqrt(sum(data2$REPLICATE_VAR))）会产生什么效果呢...？ - Ben Bolker

奇怪。你能否在某个地方发布你的数据，或者我们可以使用子集进行工作？远程调试很累人...如果你的偏差恰好是100倍，那么我认为按独立的“REPLICATE_VAR”因子除是行不通的... - Ben Bolker

嗯，我很高兴它能够工作，但我希望你能用不同的数据集或这个数据集的子集来复制SAS到R的对应关系，以确保它确实按照你所想的方式工作。 - Ben Bolker