我正在使用R中的glm将SAS PROC GENMOD示例转换为R. 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使用脊稳定的Newton-Raphson算法进行估算,即ML.我读到的关于glmR中函数的信息是结果应该等于ML.我可以做些什么来改变我在R中的估算程序,以便得到SAS中产生的等效系数和标准误差估计？
为了更新,感谢Spacedman的回答,我使用了权重,因为数据来自膳食调查中的个体,并且REPLICATE_VAR是平衡的重复复制权重,即整数(并且非常大,大约1000s或10000s).描述重量的网站在这里.我不知道为什么在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)
Spa*_*man 12
SAS中的FREQ与R的glm中的权重不同.在SAS中,它是该事件的发生次数.对于R,它的"每个响应y_i是w_i单位重量观测值的平均值".这两件事情是不一样的.
如果您希望R提供与SAS相同的输出(无法想出原因),那么您可能需要重复数据帧"权重"次数中的每一行.
这里,数据是10行,所有权重= 2,data2是20行(每行数据2个副本),所有权重= 1:
> 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观测值的平均值少,因此具有重复观测值的SE将具有比平均值更小的SE.
编辑：阅读FREQ 的 SAS 文档以及您在上面和下面的回答，我认为您应该尝试以下操作：weights=REPLICATE_VAR在glm语句中使用来调整组的相对权重（您在上面找到的系数相等表明这是正确的方法去），然后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)