Репликация пробит регрессии в SAS и R

Question

Я пытаюсь повторить свою работу SAS в R, но я получаю немного другие результаты - различия, которые нельзя объяснить ошибкой округления.

Вот мой код 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)

Я не совсем уверен, почему я получил бы другие результаты, и был бы очень благодарен за объяснение.

EDIT Вот мои данные:

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)

Answer 1

Это возможно в контрастной матрице, используемой по умолчанию.R использует контрасты лечения, в то время как SAS использует свои собственные.Посмотрите контрасты и contr.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

Answer 2

Когда я запускаю его в 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 - намного меньшеразличия, чем вы первоначально сообщили в любом случае.

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 

Answer 3

Вы должны сравнить, какое программное обеспечение сообщает о наибольшем правдоподобии. Эти числа могут отличаться только потому, что критерий завершения отличается в обоих алгоритмах. Например, большинство алгоритмов используют норму градиента в качестве правила остановки (т. Е. Когда меньше 0,0005), но каждое программное обеспечение использует свою собственную спецификацию. В зависимости от того, где он останавливается, дисперсия этих оценок будет, очевидно, отличаться, так как они получены путем инвертирования гессиана (оцениваемого там, где он останавливается). Просто чтобы быть на 100% уверенным, вы можете проверить, используя значения R или SAS, которые сообщают о самой высокой вероятности записи. Или вы можете вручную рассчитать логарифмическую вероятность, используя эти значения. Если кто-то требует от вас сообщать одинаковые значения в R и SAS, просто коснитесь критерия сходимости обоих алгоритмов. Установите какой-то очень жесткий параметр <0,00000005, в обоих случаях обе программы должны сообщать одно и то же значение. </p>

(хорошо, если ваша вероятность не имеет нескольких максимумов, что, похоже, не является проблемой здесь; в этом случае окончательные оценки будут зависеть от ваших начальных значений)

Answer 4

Это отличный источник http://sas -and-r.blogspot.com /

Answer 5

Я новичок в R, но у меня есть предложение.

Попробуйте запустить пробит, используя другой пакет R ... попробуйте Zelig.

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

Отличаются ли коэффициенты регрессии в этой модели?