I don't understand what your problem is, but
-nlcom- feeds on one or more formulae, not
the result of a numerical evaluation.
Here's an example of -nlcom- in action.
I'm fitting a negative binomial distribution
to a single response. This isn't necessarily
your problem, but it may help in indicating
some technique.
These are real data: Hilborn and Mangel (1997, p.100)
give data on incidental capture of albatrosses in the New
Zealand subantarctic squid trawl fishery, 1990. (The birds
get trapped accidentally in nets or trawl gear or cables.)
Hilborn, R. and Mangel, M. 1997.
The ecological detective: confronting models with data.
Princeton University Press, Princeton, NJ.
. input x freq
1. 0 807
2. 1 37
3. 2 27
4. 3 8
5. 4 4
6. 5 4
7. 6 1
8. 7 3
9. 8 1
10. 9 0
11. 10 0
12. 11 2
13. 12 1
14. 13 1
15. 14 0
16. 15 0
17. 16 0
18. 17 1
19. end
Suppose I don't like the parameterisation
used by -nbreg-. I can map to a new one,
such as one commonly used in ecology.
Here's the wrapper for -nlcom-:
*! NJC 1.0.0 18 Nov 2003
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']
nlcom (mean: exp(_b[_cons])) (k: exp(-_b[/lnalpha])) ///
(p: exp(_b[_cons])/ (exp(_b[_cons]) + exp(-_b[/lnalpha])))
end
Here are the results:
. mynegbin x [fw=freq]
Fitting negative binomial distribution to x
mean: exp(_b[_cons])
k: exp(-_b[/lnalpha])
p: exp(_b[_cons])/ (exp(_b[_cons]) + exp(-_b[/lnalpha]))
----------------------------------------------------------------------
--------
x | Coef. Std. Err. z P>|z| [95% Conf.
Interval]
-------------+--------------------------------------------------------
--------
mean | 2.777778 .2549476 10.90 0.000 2.27809
3.277466
k | 2.511678 .6363871 3.95 0.000 1.264382
3.758974
p | .5251538 .0672001 7.81 0.000 .393444
.6568636
----------------------------------------------------------------------
--------
-nlcom- does the hard work for you.
Incidentally, this little problem of mapping from one negative
binomial parameterisation to another is also tackled in an ad hoc and
Stata 6 way by -nbfit- on SSC. -nlcom- lets you do the job
properly.
Nick
[email protected]
Jun Xu
>
> Not sure if any solution or ideas about this problem. I try
> to write some
> routine for
> getting the prediction and confidence interval for
> predicted probabilities
> in nbreg. However,
> I got into some trouble. Note the following fourth line with:
>
> .local ai = "1/exp([lnalpha]_b[_cons])" // pay
> attention to this
> line
> I think it is this line giving me trouble. I will get both
> nonsensical
> prediction and for sure incorrect standard errors. And if I have
>
> .local ai "1/exp([lnalpha]_b[_cons])" // pay
> attention to this line
> Again, I get exactly teh same as I did using the line above.
>
> If I take off the quotation mark;that is, if I have
>
> .local ai = 1/exp([lnalpha]_b[_cons]) // pay
> attention to this line
> then the prediction is CORRECT (I checked it), however,
> because the local is
> resolved before getting into nlcom, so the confidence
> interval is incorrect
> for sure.
>
> I got totally confused and don't know what to do. Also note
> that `expxb' has
> been defined to be equal to exp(xb)=mu and `i' refers to
> the number of times
> (y). Because I have tested it for many times and guess the
> line that I just
> talked about giving me trouble, especially how it is
> translated in nlcom
> affects the differences in the output, and I am pretty sure
> it's not the
> length of the expression list giving me trouble. Very
> frustrated, and hope
> someone could help out. In addition, I think the density
> function I specify
> for the nbreg is correct.
>
> *****************************************************************
> local alpha "exp([lnalpha]_b[_cons])"
> *local alpha = "e(alpha)"
> *local ai = 1/`alpha'
> local ai = "1/exp([lnalpha]_b[_cons])" //
> pay attention to
> this line
> local gai "exp(lngamma(`ai'))"
> noi di in y `alpha'
> noi di in y `ai'
> noi di in y `expxb'
> nlcom
> ((exp(lngamma(`i'+`ai'))/(exp(lnfact(`i'))*exp(lngamma(`ai'))))* ///
> ((`ai'/(`ai'+`expxb'))^`ai')* ///
> ((`expxb'/(`ai'+`expxb'))^`i'))
> ******************************************************************
>
*
* 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/
