# st: Re: statalist-digest V4 #2986

 From Dirk Enzmann To statalist@hsphsun2.harvard.edu Subject st: Re: statalist-digest V4 #2986 Date Thu, 28 Feb 2008 14:52:25 +0100

Nick, you wished to know what I'm doing with -fitnb-. Its similarity to -nbfit- is no accident: I try to identify outlying cases of a presumable negative binomial distributed count variable.

Thanks for your hints, I tried to modify the program accordingly. See below:

=================================================
capture program drop fitnb
program define fitnb, rclass
* -----------------------------------------------------------------
* fitnb identifies likely outlying cases of a count variable `x'. A
* case is defined as outlying if its value is bigger than the value
* for which less than half a case is expected assuming a negative
* binomial distribution of x with parameters 'size' and 'mu'.
* -----------------------------------------------------------------
version 9
syntax varname(numeric) [if] [in]
marksample touse
quietly count if `touse'
if r(N) == 0 {
display "(no observations)"
exit 2000
}
local x `varlist'
tempname x_max size mu N prob i nbdat maxfit pmisfit
sum `x' if `touse'
display ""
scalar `x_max' = r(max)
if `x_max'==0 {
display "(`x' is zero for all observations)"
return scalar N=r(N)
return scalar max=`x_max'
return scalar maxfit=`x_max'
exit 498
}
quietly nbreg(`x') if `touse'
scalar `size' = 1/exp(_b[/lnalpha])
scalar `mu' = exp(_b[_cons])
scalar `N' = e(N)
scalar `prob' = `size'/(`size'+`mu')
scalar `i'=0
scalar `nbdat'=1
while (`i'<=`x_max' & `nbdat'>0) {
if `i' > 0 {
scalar `nbdat' = round((ibeta(`size',`i'+1,`prob')-ibeta(`size',`i',`prob'))*`N')
}
else {
scalar `nbdat' = round(ibeta(`size',`i'+1,`prob')*`N')
}
if `nbdat' > 0 {
scalar `i'=`i'+1
}
}
scalar `maxfit' = `i'-1
quietly count if `touse' & `x'<=`maxfit'
if `N'==r(N) {
display "(no outlying cases)"
}
else {
display "without outlying cases:"
sum `x' if `touse' & `x'<=`maxfit'
scalar `pmisfit' = 100*(1-r(N)/`N')
display ""
display "outlying cases > " `maxfit' ": n = " `N'-r(N) " (" round(`pmisfit',0.001) "%)"
display ""
display "outlying values > " `maxfit' ":"
tab `x' if `touse' & `x' > `maxfit'
}
return scalar N=`N'
return scalar max=`x_max'
return scalar maxfit=`maxfit'
end
=================================================

Dirk

Nick wrote:

On a different note, your program reminds me of -nbfit- on SSC, written
by Roberto Gutierrez and friend. As the friend in question, I am interested to know what you are doing
that is different. Perhaps -nbfit- should support it.
Your program starts something like this. I am omitting details
irrelevant to the points I am going
to make.
program fitnb
version 9
syntax [varlist] [, maxval(integer 0)]
tokenize `varlist'
args x if `maxval' == 0 local selection "if `x' < ."
else local selection "if `x' < . & `x' <= `maxval'"
sum `x' `selection'

First, it seems that you want to feed this a single variable.
program fitnb
version 9
syntax varname(numeric) [, maxval(integer 0)]
local x `varlist'
if `maxval'==0 local selection "if `x'<."
else local selection "if `x' < . & `x' <= `maxval'"
sum `x' `selection'

is then a cleaner way to start that will trap some kinds of
inappropriate input.
Second, stipulating within your program -if `x' < .- is redundant in the
case of commands like -summarize- or -nbreg- which automatically omit
observations that are missing.
What -maxval()- does, as far as your visible code is concerned, can be
done easily and directly if your program supports -if-. (Adding -in- does no harm either.)
program fitnb
version 9
syntax varname [if] [in]

marksample touse qui count if `touse' if r(N) == 0 exit 2000
local x `varlist'
sum `x' if `touse'
Then you can call with
fitnb foobar if foobar <= 20
Or whatever.
Nick n.j.cox@durham.ac.uk
*************************************************
Dr. Dirk Enzmann
Institute of Criminal Sciences
Dept. of Criminology
Schlueterstr. 28
D-20146 Hamburg
Germany

phone: +49-(0)40-42838.7498 (office)
+49-(0)40-42838.4591 (Mrs Billon)
fax: +49-(0)40-42838.2344
email: dirk.enzmann@uni-hamburg.de
www: http://www2.jura.uni-hamburg.de/instkrim/kriminologie/Mitarbeiter/Enzmann/Enzmann.html
*************************************************
*
* 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/