polr probit versus stata oprobit
Thank you Thomas for your answer. It was the weights that are giving me problems and I still have no idea why. i.e. when I try your example, everything work fine. However when I do not include the weights=Freq and [fw=Freq] in both softwares, I do get verry different results. Jean,
On Thu, 11 Nov 2004, Thomas Lumley wrote:
On Wed, 10 Nov 2004, Jean Eid wrote:
Dear Thomas, Where you also able to replicate the second example? (the exaample that I turned the housing data into numerical variables) That is the one that my estimates differ.
I don't have your second example, but I get the same results from
polr(formula = Sat ~ as.numeric(Infl) + as.numeric(Type) +
as.numeric(Cont), data = housing, weights = Freq, method = "probit")
and
oprobit Sat Infl Type Cont [fw=Freq]
for example.
-thomas
On Wed, 10 Nov 2004, Thomas Lumley wrote:
On Wed, 10 Nov 2004, Jean Eid wrote:
Dear All, I have been struggling to understand why for the housing data in MASS library R and stata give coef. estimates that are really different. I also tried to come up with many many examples myself (see below, of course I did not have the set.seed command included) and all of my `random' examples seem to give verry similar output. For the housing data, I have changed the data into numeric vectors instead of factors/ordered factors. I did so to try and get the same results as in STATA and to have the housing example as close as possible to the one I constructed. I run a debian sid, kernel 2.4, R 2.0.0, and STATA version 8.2, MASS version 7.2-8. here's the example ( I assume that you have STATA installed and can run in batch mode, if not the output is also given below)
That example shows the same results with Stata and polr() from MASS.
For the housing data, I also get the same coefficients in Stata as with
polr():
In R:
library(MASS)
library(foreign)
write.dta(housing, file="housing.dta")
house.probit<-polr(Sat ~ Infl + Type + Cont, data = housing, weights =
Freq, method = "probit")
summary(house.probit)
-------------------------
Re-fitting to get Hessian
Call:
polr(formula = Sat ~ Infl + Type + Cont, data = housing, weights = Freq,
method = "probit")
Coefficients:
Value Std. Error t value
InflMedium 0.3464233 0.06413706 5.401297
InflHigh 0.7829149 0.07642620 10.244063
TypeApartment -0.3475372 0.07229093 -4.807480
TypeAtrium -0.2178874 0.09476607 -2.299213
TypeTerrace -0.6641737 0.09180004 -7.235005
ContHigh 0.2223862 0.05812267 3.826153
Intercepts:
Value Std. Error t value
Low|Medium -0.2998 0.0762 -3.9371
Medium|High 0.4267 0.0764 5.5850
Residual Deviance: 3479.689
AIC: 3495.689
------------------------
In Stata
-----------------
. use housing.dta
. xi: oprobit Sat i.Infl i.Type i.Cont [fw=Freq]
i.Infl _IInfl_1-3 (naturally coded; _IInfl_1 omitted)
i.Type _IType_1-4 (naturally coded; _IType_1 omitted)
i.Cont _ICont_1-2 (naturally coded; _ICont_1 omitted)
Iteration 0: log likelihood = -1824.4388
Iteration 1: log likelihood = -1739.9254
Iteration 2: log likelihood = -1739.8444
Ordered probit estimates Number of obs = 1681
LR chi2(6) = 169.19
Prob > chi2 = 0.0000
Log likelihood = -1739.8444 Pseudo R2 = 0.0464
------------------------------------------------------------------------------
Sat | Coef. Std. Err. z P>|z| [95% Conf.
Interval]
-------------+----------------------------------------------------------------
_IInfl_2 | .3464228 .064137 5.40 0.000 .2207165 .472129
_IInfl_3 | .7829146 .076426 10.24 0.000 .6331224 .9327069
_IType_2 | -.3475367 .0722908 -4.81 0.000 -.4892241 -.2058493
_IType_3 | -.2178875 .094766 -2.30 0.021 -.4036254 -.0321497
_IType_4 | -.6641735 .0917999 -7.24 0.000 -.844098 -.484249
_ICont_2 | .2223858 .0581226 3.83 0.000 .1084676 .336304
-------------+----------------------------------------------------------------
_cut1 | -.2998279 .0761537 (Ancillary parameters)
_cut2 | .4267208 .0764043
------------------------------------------------------------------------------
-thomas
______________________________________________ R-help at stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Thomas Lumley Assoc. Professor, Biostatistics tlumley at u.washington.edu University of Washington, Seattle