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

["Jun Xu" <mystata@hotmail.com>]

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/