Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

st: RE: RE: local macro in nlcom [nlcom bug?] Response to Nick


From   "Jun Xu" <mystata@hotmail.com>
To   statalist@hsphsun2.harvard.edu
Subject   st: RE: RE: local macro in nlcom [nlcom bug?] Response to Nick
Date   Tue, 18 Nov 2003 10:46:35 -0600

Nick,

I have to say, thank you very much. To reduce the information load, I was reluctant to respond with a thank-you note for your valuable and many quick responses to my stupid questions and through the process, I have improved my stata programming skill. I think there should be some best Stata helper of the year award for your effort to help out both novice and experienced users. Though I think SAS is more powerful than Stata in some respects (well vice versa), but the main reason why I converted to use Stata for most of my research projects is because we have a wonderful Stata user community and people like you. That said, I am able to reprodce the problem I talked about using the data you send to me, and I think it is a Stata bug in nlcom (well, with 95% subjective confidence interval):

1) here is your data, just in case any mistakes in the process of my transfer:
******
x freq
0 807
1 37
2 27
3 8
4 4
5 4
6 1
7 3
8 1
9 0
10 0
11 2
12 1
13 1
14 0
15 0
16 0
17 1
******

Here is a slightly modified program of yours
******
cap program drop mynegbin
program mynegbin
version 8.0
syntax varname(numeric) [if] [in] [fweight]
marksample touse
di _n as txt "Fitting negative binomial distribution to `varlist'"
qui nbreg `varlist' if `touse' [`weight' `exp']
local expxb "exp(_b[_cons])"
local alpha "exp(-_b[/lnalpha])"
local ai = 1/exp(-_b[/lnalpha])
local i=0
nlcom ((exp(lngamma(`i'+`ai'))/(exp(lnfact(`i'))*exp(lngamma(`ai'))))* ///
((`ai'/(`ai'+`expxb'))^`ai')* ///
((`expxb'/(`ai'+`expxb'))^`i'))
end

mynegbin x [fw=freq]
******

In this case, note that "local ai=1/exp(-_b[/lnalpha])", and the results:

Fitting negative binomial distribution to x

_nl_1: (exp(lngamma(0+16.089697466578))/(exp(lnfact(0))*exp(lngamma(16.089697466578))))*

((16.089697466578/(16.089697466578+exp(_b[_cons])))^16.089697466578)* ((exp(_b[_cons])/(16.0896
97466578+exp(_b[_cons])))^0)
------------------------------------------------------------------------------
x | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_nl_1 | .7585698 .0307805 24.64 0.000 .6982412 .8188985
------------------------------------------------------------------------------

This time, let's try a slightly different version, let's make "local ai = 1/exp(-_b[/lnalpha])" (and the same for using: local ai "1/exp(-_b[/lnalpha])" ) and the program now becomes:

******
cap program drop mynegbin
program mynegbin
version 8.0
syntax varname(numeric) [if] [in] [fweight]
marksample touse
di _n as txt "Fitting negative binomial distribution to `varlist'"
qui nbreg `varlist' if `touse' [`weight' `exp']
local expxb "exp(_b[_cons])"
local alpha "exp(-_b[/lnalpha])"
local ai = "1/exp(-_b[/lnalpha])"
local i=0
nlcom ((exp(lngamma(`i'+`ai'))/(exp(lnfact(`i'))*exp(lngamma(`ai'))))* ///
((`ai'/(`ai'+`expxb'))^`ai')* ///
((`expxb'/(`ai'+`expxb'))^`i'))
end

mynegbin x [fw=freq]
******
Now the results become:


_nl_1: (exp(lngamma(0+1/exp(-_b[/lnalpha])))/(exp(lnfact(0))*exp(lngamma(1/exp(-_b[/lnalp

ha])))))* ((1/exp(-_b[/lnalpha])/(1/exp(-_b[/lnalpha])+exp(_b[_cons])))^1/exp(-_b[/lnalpha]))* ((exp(_b[_cons])/(1/exp(-_b[/lnalpha])+exp(_b[_cons])))^0)
------------------------------------------------------------------------------
x | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_nl_1 | 15.81574 2.407928 6.57 0.000 11.09628 20.53519
------------------------------------------------------------------------------

The striking differences between the two predictions puzzled me and I really appreicate your example because I am very confident it is how that local gets resolved affected the results in my own example and I suspect this might be a bug in nlcom. Probably it won't be able to deal with local string with certain characteristics?

Jun Xu
Department of Sociology
Indiana University

_________________________________________________________________
Send a QuickGreet with MSN Messenger http://www.msnmessenger-download.com/tracking/cdp_games

*
* 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