Statalist


[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

RE: st: gologit2


From   "Martin Weiss" <martin.weiss@uni-tuebingen.de>
To   <statalist@hsphsun2.harvard.edu>
Subject   RE: st: gologit2
Date   Fri, 18 Apr 2008 15:33:40 +0200

Very much apart from the heated discussion in this thread, could anyone tell
me what - list retokenize- does. Cannot find it anywhere. -findit list
retokenize- does not throw up anything, either.

Furthermore, I would be interested how long Maarten`s original code (10,000
reps) a couple of days ago took to execute for you. I ran it last night on
2GB Ram, Vista SP1 32-bit, with AMD 64 2.2 GhZ dual core under Stata 10.0
MP/2, and it took 7.3 hours (26212 seconds). 

Martin Weiss
_________________________________________________________________

Diplom-Kaufmann Martin Weiss
Mohlstrasse 36
Room 415
72074 Tuebingen
Germany

Fon: 0049-7071-2978184

Home: http://www.wiwi.uni-tuebingen.de/cms/index.php?id=1130

Publications: http://www.wiwi.uni-tuebingen.de/cms/index.php?id=1131

SSRN: http://papers.ssrn.com/sol3/cf_dev/AbsByAuth.cfm?per_id=669945

-----Original Message-----
From: owner-statalist@hsphsun2.harvard.edu
[mailto:owner-statalist@hsphsun2.harvard.edu] On Behalf Of Richard Williams
Sent: Friday, April 18, 2008 4:18 PM
To: statalist@hsphsun2.harvard.edu
Subject: RE: st: gologit2

For anyone who is still following this thread - Maarten and I have 
been conversing offline and I think we both agree the original 
simulations of the brant test and its counterparts were flawed and 
the test works much better than we initially thought.  Here are new 
simulations.  For each simulation, the expected value for reject H0 
is .05, i.e. 50 false rejections in 1000 tries using the .05 level of 
significance.

Brant Simulations, 1000 reps

# of Xs   | reject H0
----------+-----------
         1 |      .051
         2 |      .042
         3 |      .054
         4 |      .051
         5 |      .057
         6 |      .054
         7 |      .048
         8 |      .056
         9 |      .050
        10 |      .056
        12 |      .067
        14 |      .038
        16 |      .048
        18 |      .041
        20 |      .036

Altogether, the observed number of false rejections of the null 
totaled 749; the expected number of rejections was 750, so I'd say 
the simulations worked pretty good.  Most critically, you don't see 
an increased tendency to reject the null as the number of x variables 
increases.

For anyone who is interested as to what has changed - in Maarten's 
original code, as the number of Xs increased, the distribution of Y 
also changed.  Cases were less likely to be in the middle categories 
of 2 and 3 and more likely to be in the extreme categories of 1 and 
4.  I don't think this would necessarily be bad, except that it does 
result in some very thin cell counts in some simulations which 
probably results in poor estimation.

I tweaked Maarten's code so that the distribution of Y was always the 
same no matter how many Xs were in the model.  This seems more 
realistic anyway, i.e. the distribution of the observed y shouldn't 
change just because you have more explanatory variables.  That 
produced the results above.  Also, Maarten reports that he tweaked 
his original code, increasing the sample size, and the brant test 
also performed better then.

Here is the modified code.

set more off
set seed 12345


capture program drop sim


program define sim, rclass
         syntax, [nx(integer 1)]
         drop _all
         set obs 500
         forvalues i = 1/`nx' {
                 gen x`i' = invnorm(uniform())
                 local x `x' x`i'
         }
         local x : list retokenize x
         local xsum : subinstr local x " " " + ", all
         gen u = uniform()
         gen ystar = `xsum' + ln(u/(1-u))
         * original code for y.  This code causes
         * the distribution of y to change as the
         * number of x vars increases, probably
         * resulting in very small cell counts in
         * some simulations.
         *gen y = cond(ystar < -2, 1,     ///
         *  cond(ystar <  0, 2,     ///
         *        cond(ystar <  2, 3, 4)))

         * New code for y.  Distribution for y is
         * the same no matter how many Xs there are.
         egen y = cut(ystar), group(4)
         replace y = y + 1

         * Uncomment the section you want
         * brant test code
         ologit y `x'
         brant
         return scalar p = r(p)

         * omodel test code
         *omodel logit y `x'
         *return scalar p = chi2tail($S_2, $S_1)

         * gologit2 test code
         *ologit y `x'
         *est store m1
         *gologit2 y `x', npl store(m2)
         *lrtest m1 m2, force
         *return scalar p = r(p)
end


simulate p=r(p), reps(1000): sim, nx(1)
count if p < .05
matrix res = 1, r(N)/1000


foreach i of numlist 2/10 12(2)20 {
         simulate p=r(p), reps(1000): sim, nx(`i')
         count if p < .05
         matrix res = res \ `i', r(N)/1000
}
matlist res


*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/

*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index