SAS和R中的重复概率单位回归

ix0qys7i 于 2023-01-18 发布在其他

关注(0)|答案(4)|浏览(125)

我试图在R中重复我的SAS工作，但是我得到了稍微不同的结果--这些差异不能用舍入误差来解释。
以下是我的SAS代码：

proc qlim data=mydata;
   model y = x1 x2 x3/ discrete(d=probit);
   output out=outdata marginal;
   title "just ran QLIM model";
run;
quit;

这是我的R代码：

mymodel <- glm(y ~ x1 + x2 + x3, family=binomial(link="probit"), data=mydata)

我真的不知道为什么我会得到不同的结果，并会非常感谢一个解释。

编辑：

以下是我的数据：

2.66  20  0  0
2.89  22  0  0
3.28  24  0  0
2.92  12  0  0
4.00  21  0  1
2.86  17  0  0
2.76  17  0  0
2.87  21  0  0
3.03  25  0  0
3.92  29  0  1
2.63  20  0  0
3.32  23  0  0
3.57  23  0  0
3.26  25  0  1
3.53  26  0  0
2.74  19  0  0
2.75  25  0  0
2.83  19  0  0
3.12  23  1  0
3.16  25  1  1
2.06  22  1  0
3.62  28  1  1
2.89  14  1  0
3.51  26  1  0
3.54  24  1  1
2.83  27  1  1
3.39  17  1  1
2.67  24  1  0
3.65  21  1  1
4.00  23  1  1
3.1   21  1  0
2.39  19  1  1

下面是我估计的系数（括号中为标准误差）：

SAS: -7.452320 (2.542536)
      1.625810 (0.693869)
      0.051729 (0.083891)
      1.426332 (0.595036)
R:   -7.25319  (2.50977)
      1.64888  (0.69427)
      0.03989  (0.07961)
      1.42490  (0.58347)

来源：https://stackoverflow.com/questions/3391644/replicating-probit-regression-in-sas-and-r

4条答案

按热度按时间

iecba09b1#

它可能在默认使用的对比矩阵中。R使用治疗对比，而SAS使用自己的对比。在帮助中查找对比和对照SAS。如果您经常使用SAS对比，您可能需要将选项设置为该值。

options(contrasts=c("contr.SAS", "contr.poly"))

为了了解这是如何影响事情的，观察治疗和SAS对比矩阵的差异

contr.treatment(4)
  2 3 4
1 0 0 0
2 1 0 0
3 0 1 0
4 0 0 1

contr.SAS(4)
  1 2 3
1 1 0 0
2 0 1 0
3 0 0 1
4 0 0 0

赞(0）回复(0）举报 2023-01-18

6pp0gazn2#

当我用你的数据和代码在R中运行它时，我得到的答案（接近）是你为SAS结果显示的：

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept) -7.45231    2.57152  -2.898  0.00376 **
x1           1.62581    0.68973   2.357  0.01841 * 
x2           0.05173    0.08119   0.637  0.52406   
x3           1.42633    0.58695   2.430  0.01510 *

标准误差有百分之几的偏差，但这并不令人惊讶。
我还在glmmADMB（R-forge上提供）中运行了它，这是一个完全不同的实现，得到的估计值与SAS稍有不同，但标准误差更接近SAS --无论如何，与您最初报告的差异都要小得多。

library(glmmADMB)
> mm2 <- glmmadmb(y~x1+x2+x3,family="binomial",link="probit",data=mydata)
["estimated covariance may be non-positive-definite warnings"]
> summary(mm2)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  -7.4519     2.5424   -2.93   0.0034 **
x1            1.6258     0.6939    2.34   0.0191 * 
x2            0.0517     0.0839    0.62   0.5375   
x3            1.4263     0.5950    2.40   0.0165 *

您使用的是哪个版本的R？（虽然glm是非常稳定的代码，但 * 可能 * 版本之间会发生一些变化...）您确定没有搞砸什么吗？

> sessionInfo()
R Under development (unstable) (2011-10-06 r57181)
Platform: i686-pc-linux-gnu (32-bit)

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] glmmADMB_0.6.4

赞(0）回复(0）举报 2023-01-18

thtygnil3#

我是个R新手，但我有个建议。
尝试使用另一个R包运行概率单位...尝试Zelig。

mymodel <- zelig(y ~ x1 + x2 + x3, model="probit", data=mydata)
summary(mymodel)

这个模型中的回归系数是否不同？

赞(0）回复(0）举报 2023-01-18

5sxhfpxr4#

您应该比较哪个软件报告的对数似然最高。这些数字可能不同，只是因为两种算法的终止条件不同。例如，大多数算法使用梯度范数作为停止规则（即：当小于0.0005时），但是每个软件使用其自己的规范。取决于其停止的位置，这些估计的方差将明显不同，因为它们是通过反转Hessian（在其停止的位置评估）获得的。为了100%确定，你可以用R或SAS值来检验，这些值报告了最高的对数似然。或者你可以用这些值手工计算对数似然。如果有人要求你在R和SAS中报告完全相同的值，只需点击两种算法的收敛条件。设置一些非常严格的参数〈0.0000005，在这两种情况下，两个程序应该报告相同的值。
（除非你的可能性有多个极大值，这似乎不是问题所在;在这种情况下，最终估计值将取决于您的初始值）

赞(0）回复(0）举报 2023-01-18

我来回答

SAS和R中的重复概率单位回归

4条答案

相关问题

热门标签

最新问答