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

st: RE: local macro in nlcom


From   "Nick Cox" <n.j.cox@durham.ac.uk>
To   <statalist@hsphsun2.harvard.edu>
Subject   st: RE: local macro in nlcom
Date   Tue, 18 Nov 2003 11:50:04 -0000

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
n.j.cox@durham.ac.uk

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/



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