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

 From "Jun Xu" 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

_________________________________________________________________