*! version 0.8.1 - 2009-07-31
*   work around S11 bug in estimates table

*  To do:
*   1) add weights
*   2) use returns rather than scalars

// compare count models

capture program drop countfit
program define countfit, rclass
    version 9.0

    * 2009-07-31
    local vwidth = 32 // causes error in Stata 11
    local vwidth = 30

    syntax varlist [if] [in] ///
    , [Generate(string)  replace ///
      INFlate(string) /// inflation variables
      MAXcount(integer 9) ///
      NOGraph /// suppress graph of differences in predicted probabilities
      NODifferences /// suppress table of differences in predicted probabilities
      NOEstimates /// suppress table of estimated parameters
      NOFit /// suppress table of fit statistics
      Prm Nbreg ZIP ZINb ///
      nodash ///
      NOConstant ///
      NOPRTable ///
      note(string) /// note for graph
      GRAPHexport(string) /// options passed to graph export
      NOIsily ///
      ]

// trap weights

    if "`e(wtype)'"!="" {
        di in r "-countfix- does not work with weights."
        exit
    }

//  define variable lists

    * inflate is same as rhs if not specified
    if "`inflate'"=="" {
        ** better way to do this?
        local nvars : word count `varlist'
        foreach i of numlist 2(1)`nvars' {
            local var : word `i' of `varlist'
            local inflate "`inflate' `var'"
        }
    }

// set up printing

    local f83 "%8.3f"
    local f93 "%9.3f"
    local f93l "%-9.3f"
    local f103 "%10.3f"
    local dash ""
    if "`nodash'"=="" {
        local dash ///
        "--------------------------------------------------------------------------"
    }

// construct information on models being assessed

    * defalult name is CF for variables created
    if "`generate'"=="" local set "CF"
    local set "`generate'" // label for sets of models
    local modellist "" // list of types of models estimated
    local mdlnmlist "" // like modellist but begin with generate prefix _
    local pltsym "" // symbols used in the plot
    local plty "" // y variables to plot
    local pltx "" // x variables to plot
    local mdlnum = 1 // number associated with each model

    * set up information for all models specified
    if "`prm'"=="prm" {
        local mdlnmlist "`mdlnmlist' `set'PRM"
        local modellist "PRM "
        local pltx "`set'PRMval"
        local plty "`set'PRMdif "
        local pltsym "Th "
        local pltleg1 "`set'PRM"
        local mdlnum = `mdlnum' + 1
    }
    if "`nbreg'"=="nbreg" {
        local mdlnmlist "`mdlnmlist' `set'NBRM"
        local modellist "`modellist'NBRM "
        local pltx "`set'NBRMval"
        local plty "`plty' `set'NBRMdif "
        local pltsym "`pltsym' Sh "
        local pltleg`mdlnum' "`set'NBRM"
        local mdlnum = `mdlnum' + 1
    }
    if "`zip'"=="zip" {
        local mdlnmlist "`mdlnmlist' `set'ZIP"
        local modellist "`modellist'ZIP "
        local pltx "`set'ZIPval"
        local plty "`plty' `set'ZIPdif "
        local pltsym "`pltsym' T "
        local pltleg`mdlnum' "`set'ZIP"
        local mdlnum = `mdlnum' + 1
    }
    if "`zinb'"=="zinb" {
        local mdlnmlist "`mdlnmlist' `set'ZINB"
        local modellist "`modellist'ZINB "
        local pltx "`set'ZINBval"
        local plty "`plty' `set'ZINBdif "
        local pltsym "`pltsym' S "
        local pltleg`mdlnum' "`set'ZINB"
        local mdlnum = `mdlnum' + 1
    }
    * if none, then all
    if "`prm'"=="" & "`nbreg'"=="" & "`zip'"=="" & "`zinb'"=="" {
        local pltx "`set'PRMval"
        local plty "`set'PRMdif `set'NBRMdif `set'ZIPdif `set'ZINBdif"
        local pltsym "Th Sh T S"
        local pltleg1 "`set'PRM"
        local pltleg2 "`set'NBRM"
        local pltleg3 "`set'ZIP"
        local pltleg4 "`set'ZINB"
        local mdlnmlist "`set'PRM `set'NBRM `set'ZIP `set'ZINB"
        local modellist "PRM NBRM ZIP ZINB"
        local prm "prm"
        local nbreg "nbreg"
        local zip "zip"
        local zinb "zinb"
        * drop? local alphaopt "drop(lnalpha:_cons)"
    }

// estimate models

    * 2009-07-31
    _count_estimate `varlist' `if' `in' [`fweight' `pweight' `iweight'], ///
        max(`maxcount') gen(`generate') `replace' inflate(`inflate') ///
        `prm' `nbreg' `zip' `zinb' `noconstant' `noisily'

// table of estimates


    if "`noestimates'"=="" {
        estimates table `mdlnmlist', eform ///
            b(%9.3f) t(%7.2f) label varwidth(`vwidth') ///
            stats(alpha N ll bic aic)
    }

// summary table of mean observed - predicted probabilities

    local c1 = 10
    local c2 = 25
    local c3 = 31
    local c4 = 47

    if "`nodifferences'"=="" {

        di in g "Comparison of Mean Observed and Predicted Count"
        di
        di in g "     " _col(`c1') "   Maximum"    ///
            _col(`c2') "  At" ///
            _col(`c3') "    Mean"
        di in g "Model" _col(`c1') " Difference" ///
            _col(`c2') "Value" ///
            _col(`c3') "   |Diff|"
        di in g "---------------------------------------------"
    }

    * loop through models and compute differences
    foreach m in `modellist' {

        local modelnm "`set'`m'"

        * stats on difference
        qui sum `modelnm'dif
        local toplot "`toplot' `modelnm'dif"
        local `m'difsd = r(sd) // sd of obs-pred
        local `m'difmin = r(min) // min
        local `m'difmax = r(max) // max

        * stats on mean absolute difference
        capture drop `modelnm'difabs
        qui gen `modelnm'absdif = abs(`modelnm'dif)
        qui sum `modelnm'absdif
        local `m'absdifmax = r(max) // max
        scalar absdifsd`m' = r(sd) // sd
        scalar absdifmn`m' = r(mean) // mean
        * find values for largest difference
        tempname difval
        qui gen `difval' = (`modelnm'absdif>``m'absdifmax'-.00001)*_n ///
                if `modelnm'absdif!=.
        qui sum `difval'
        local `m'absdifmaxval = r(max) - 1

        * sign of deviation
        tempname maxdif // 0.8.0
        if ``m'absdifmax' == abs(``m'difmin') {
        *local maxdif = ``m'difmin' // 0.8.0
            scalar `maxdif' = ``m'difmin'
        }
        if ``m'absdifmax' == abs(``m'difmax') {
        *local maxdif = ``m'difmax' // 0.8.0
            scalar `maxdif' = ``m'difmax'
        }

        if "`nodifferences'"=="" {
            * print summary of differences from predictions
            di in y "`modelnm'" _col(`c1') in y `f83' `maxdif' ///
                _col(`c2') `f83' "  ``m'absdifmaxval'" ///
                _col(`c3') `f83' absdifmn`m'
        }

    } // loop through models

// TABLES OF OBSERVED AND PREDICTED COUNTS

    if "`noprtable'" == "" {

        foreach t in `mdlnmlist' {

            qui {
                sum `t'absdif
                local max = r(N) - 1 // largest count
                tempname sumde sumpe sumob sumpr // 0.8.0
                scalar `sumde' = r(sum) // sum of abs dif
                sum `t'pearson
                scalar `sumpe' = r(sum) // sum of pearson dif
                sum `t'obeq
                scalar `sumob' = r(sum) // sum of pr observed
                sum `t'preq
                scalar `sumpr' = r(sum) // sum of pr predicted
            } // qui
            local c1 = 7
            local c2 = 19
            local c3 = 30
            local c4 = 40
            local c5 = 50
            di
            di in y "`t'" in g ": Predicted and actual probabilities"
            di
            di in g "Count" _col(`c1') "  Actual" ///
                _col(`c2') "Predicted" ///
                _col(`c3') "  |Diff|" ///
                _col(`c4') " Pearson"

            local dash ""
            if "`nodash'"=="" {
                local dash ///
                "------------------------------------------------"
            }

            di in g "`dash'"

            foreach c of numlist 0(1)`max' {
                local i = `c' + 1
                tempname ob pr de pe // 0.8.0
                scalar `ob' = `t'obeq[`i']
                scalar `pr' = `t'preq[`i']
                scalar `de' = abs(`t'dif[`i'])
                scalar `pe' = `t'pearson[`i']

                di in y "`c'" ///
                    _col(`c1') `f83' `ob' ///
                    _col(`c2') `f83' `pr' ///
                    _col(`c3') `f83' `de' ///
                    _col(`c4') `f83' `pe'

            } // loop through counts

            di in g "`dash'"
            di in y "Sum" ///
                _col(`c1') `f83' `sumob' ///
                _col(`c2') `f83' `sumpr' ///
                _col(`c3') `f83' `sumde' ///
                _col(`c4') `f83' `sumpe'

        } // loop through modellist

    } // print table of predictions

// PLOT DIFFERENCES

    if "`nograph'"=="" {

        twoway (connected `plty' `pltx', ///
            msymbol(`pltsym') ///
            clpat(tight_dot tight_dot tight_dot tight_dot ) ///
            ytitle("Observed-Predicted") ///
            subtitle("Note: positive deviations show underpredictions.",  ///
            pos(11) size(small)) ///
            ylabel(-.10(.05).10, grid gmax gmin) ///
            xlabel(0(1)`maxcount') ///
            note(`note') ///
            ysize(3.5) xsize(4.5) ///
            legend(order(1 "`pltleg1'" 2 "`pltleg2'" ///
            3 "`pltleg3'" 4 "`pltleg4'")) ///
        )

        * export graph
        if "`graphexport'" != "" {
            qui graph export `graphexport'
        }

    }

// COMPARE FIT STATISTICS

    if "`nofit'"=="" {

        local dash ""
        if "`nodash'"=="" {
        local dash ///
    "-------------------------------------------------------------------------"
        }
        local c1 = 16
        local c2 = 32
        local c3 = 48
        local c4 = 56
        local c5 = 62

        di _n in g "Tests and Fit Statistics"
        di
        tempname aicd  // 0.8.0

        * base PRM
        if "`prm'"!="" {
            local mdl1 "`set'PRM"
            local mdl1type "PRM"
            di in y "`mdl1'" ///
                _col(`c1') in g "BIC=" in y `f103' bic`mdl1' ///
                _col(`c2') in g "AIC=" in y `f103' aic`mdl1' ///
                _col(`c3') in g "Prefer" ///
                _col(`c4') in g "Over" ///
                _col(`c5') in g "Evidence"

            * prm vs nbreg
            if "`nbreg'"!="" {
                local mdl2 "`set'NBRM"
                local mdl2type "NBRM"
                global bic1 = bic`mdl1'
                global bic2 = bic`mdl2'
                _count_bic_dif
                local pref = $bicpref
                local nopref = $bicnopref
                di "`dash'"
                * bic
                di in y " " in g " vs " in y "`mdl2'" ///
                    _col(`c1') in g "BIC=" in y `f103' bic`mdl2' ///
                    _col(`c2') in g "dif=" in y `f103' $bicdif ///
                    _col(`c3') in y "`mdl`pref'type'" ///
                    _col(`c4') in y "`mdl`nopref'type'" ///
                    _col(`c5') in y "$bicsup"
                * aic
                scalar `aicd' = aic`mdl1' - aic`mdl2'
                local aicfav "`mdl2type'"
                local aicnofav "`mdl1type'"
                if aic`mdl1' < aic`mdl2' {
                    local aicfav "`mdl1type'"
                    local aicnofav "`mdl2type'"
                }
                di  _col(`c1') in g "AIC=" in y `f103' aic`mdl2' ///
                    _col(`c2') in g "dif=" in y `f103' `aicd' ///
                    _col(`c3') in y "`aicfav'" ///
                    _col(`c4') in y "`aicnofav'"
                * lr test
                local lrfav "`mdl1type'"
                local lrnofav "`mdl2type'"
                if lrnb_prmp < .05 {
                    local lrfav "`mdl2type'"
                    local lrnofav "`mdl1type'"
                }
                di _col(`c1') in g "LRX2=" in y `f93' lrnb_prm ///
                    _col(`c2') in g "prob=" in y `f93' lrnb_prmp in g ///
                    _col(`c3') in y "`lrfav'" ///
                    _col(`c4') in y "`lrnofav'" ///
                    _col(`c5') in y "p=" `f93l' lrnb_prmp
            } // no nbreg vs prm

            * prm vs zip
            if "`zip'"!="" {
                local mdl2 "`set'ZIP"
                local mdl2type "ZIP"
                * bic
                global bic1 = bic`mdl1'
                global bic2 = bic`mdl2'
                _count_bic_dif
                local pref = $bicpref
                local nopref = $bicnopref
                di in g "`dash'"
                di in y " " in g " vs " in y "`mdl2'" ///
                    _col(`c1') in g "BIC=" in y `f103' bic`mdl2' ///
                    _col(`c2') in g "dif=" in y `f103' $bicdif ///
                    _col(`c3') in y "`mdl`pref'type'" ///
                    _col(`c4') in y "`mdl`nopref'type'" ///
                    _col(`c5') in y "$bicsup"
                * aic
                tempname aicd // 0.8.0
                scalar `aicd' = aic`mdl1' - aic`mdl2'
                local aicfav "`mdl2type'"
                local aicnofav "`mdl1type'"
                if aic`mdl1' < aic`mdl2' {
                    local aicfav "`mdl1type'"
                    local aicnofav "`mdl2type'"
                }
                di _col(`c1') in g "AIC=" in y `f103' aic`mdl2' ///
                        _col(`c2') in g "dif=" in y `f103' `aicd' ///
                        _col(`c3') in y "`aicfav'" ///
                        _col(`c4') in y "`aicnofav'"
                * vuong test
                local vufav "`mdl1type'"
                local vunofav "`mdl2type'"
                if vu`mdl2'>0  {
                    local vufav "`mdl2type'"
                    local vunofav "`mdl1type'"
                }
                di _col(`c1') in g "Vuong=" in y `f83' vu`mdl2' ///
                    _col(`c2') in g "prob=" in y `f93' vu`mdl2'p in g ///
                    _col(`c3') in y "`vufav'" ///
                    _col(`c4') in y "`vunofav'" ///
                    _col(`c5') in y "p=" `f93l' vu`mdl2'p
            } // no zip vs prm

            * prm vs zinb
            if "`zinb'"!="" {
                local mdl2 "`set'ZINB"
                local mdl2type "ZINB"
                * bic
                global bic1 = bic`mdl1'
                global bic2 = bic`mdl2'
                _count_bic_dif
                local pref = $bicpref
                local nopref = $bicnopref
                di in g "`dash'"
                di in y " " in g " vs " in y "`mdl2'" ///
                    _col(`c1') in g "BIC=" in y `f103' bic`mdl2' ///
                    _col(`c2') in g "dif=" in y `f103' $bicdif ///
                    _col(`c3') in y "`mdl`pref'type'" ///
                    _col(`c4') in y "`mdl`nopref'type'" ///
                    _col(`c5') in y "$bicsup"
                * aic
                scalar `aicd' = aic`mdl1' - aic`mdl2'
                local aicfav "`mdl2type'"
                local aicnofav "`mdl1type'"

                if aic`mdl1' < aic`mdl2' {
                    local aicfav "`mdl1type'"
                    local aicnofav "`mdl2type'"
                }
                di _col(`c1') in g "AIC=" in y `f103' aic`mdl2' ///
                    _col(`c2') in g "dif=" in y `f103' `aicd' ///
                    _col(`c3') in y "`aicfav'" ///
                    _col(`c4') in y "`aicnofav'"
            } // no zinb vs prm

        } // if no prm

        * base nbreg
        if "`nbreg'"!="" {
            local mdl1 "`set'NBRM"
            local mdl1type "NBRM"
            di in g "`dash'"
            di in y "`mdl1'" ///
                _col(`c1') in g "BIC=" in y `f103' bic`mdl1' ///
                _col(`c2') in g "AIC=" in y `f103' aic`mdl1' ///
                _col(`c3') in g "Prefer" ///
                _col(`c4') in g "Over" ///
                _col(`c5') in g "Evidence"

            * nbreg vs zip
            if "`zip'"!="" {
                local mdl2 "`set'ZIP"
                local mdl2type "ZIP"
                * bic
                global bic1 = bic`mdl1'
                global bic2 = bic`mdl2'
                _count_bic_dif
                local pref = $bicpref
                local nopref = $bicnopref
                di in g "`dash'"
                di in y " " in g " vs " in y "`mdl2'" ///
                    _col(`c1') in g "BIC=" in y `f103' bic`mdl2' ///
                    _col(`c2') in g "dif=" in y `f103' $bicdif ///
                    _col(`c3') in y "`mdl`pref'type'" ///
                    _col(`c4') in y "`mdl`nopref'type'" ///
                    _col(`c5') in y "$bicsup"
                * aic
                scalar `aicd' = aic`mdl1' - aic`mdl2'
                local aicfav "`mdl2type'"
                local aicnofav "`mdl1type'"

                if aic`mdl1' < aic`mdl2' {
                    local aicfav "`mdl1type'"
                    local aicnofav "`mdl2type'"
                }
                di _col(`c1') in g "AIC=" in y `f103' aic`mdl2' ///
                    _col(`c2') in g "dif=" in y `f103' `aicd' ///
                    _col(`c3') in y "`aicfav'" ///
                    _col(`c4') in y "`aicnofav'"
             } // no zip vs nbreg

            * nbreg vs zinb
            if "`zinb'"!="" {
                local mdl2 "`set'ZINB"
                local mdl2type "ZINB"
                * bic
                global bic1 = bic`mdl1'
                global bic2 = bic`mdl2'
                _count_bic_dif
                local pref = $bicpref
                local nopref = $bicnopref
                di in g "`dash'"
                di in y " " in g " vs " in y "`mdl2'" ///
                    _col(`c1') in g "BIC=" in y `f103' bic`mdl2' ///
                    _col(`c2') in g "dif=" in y `f103' $bicdif ///
                    _col(`c3') in y "`mdl`pref'type'" ///
                    _col(`c4') in y "`mdl`nopref'type'" ///
                    _col(`c5') in y "$bicsup"
                * aic
                scalar `aicd' = aic`mdl1' - aic`mdl2'
                local aicfav "`mdl2type'"
                local aicnofav "`mdl1type'"
                if aic`mdl1' < aic`mdl2' {
                    local aicfav "`mdl1type'"
                    local aicnofav "`mdl2type'"
                }
                di _col(`c1') in g "AIC=" in y `f103' aic`mdl2' ///
                    _col(`c2') in g "dif=" in y `f103' `aicd' ///
                    _col(`c3') in y "`aicfav'" ///
                    _col(`c4') in y "`aicnofav'"
                * vuong test
                local vufav "`mdl1type'"
                local vunofav "`mdl2type'"
                if vu`mdl2'>0  {
                    local vufav "`mdl2type'"
                    local vunofav "`mdl1type'"
                }
                di _col(`c1') in g "Vuong=" in y `f83' vu`mdl2' ///
                    _col(`c2') in g "prob=" in y `f93' vu`mdl2'p in g ///
                    _col(`c3') in y "`vufav'" ///
                    _col(`c4') in y "`vunofav'" ///
                    _col(`c5') in y "p=" `f93l' vu`mdl2'p

            } // no zinb vs nbreg

        } // if no nbreg

        if "`zip'"!="" {

            local mdl1 "`set'ZIP"
            local mdl1type "ZIP"
            di in g "`dash'"
            di in y "`mdl1'" ///
                _col(`c1') in g "BIC=" in y `f103' bic`mdl1' ///
                _col(`c2') in g "AIC=" in y `f103' aic`mdl1' ///
                _col(`c3') in g "Prefer" ///
                _col(`c4') in g "Over" ///
                _col(`c5') in g "Evidence"

            * zip vs zinb
            if "`zinb'"!="" {
                local mdl2 "`set'ZINB"
                local mdl2type "ZINB"
                * bic
                global bic1 = bic`mdl1'
                global bic2 = bic`mdl2'
                _count_bic_dif
                local pref = $bicpref
                local nopref = $bicnopref
                di in g "`dash'"
                di in y " " in g " vs " in y "`mdl2'" ///
                    _col(`c1') in g "BIC=" in y `f103' bic`mdl2' ///
                    _col(`c2') in g "dif=" in y `f103' $bicdif ///
                    _col(`c3') in y "`mdl`pref'type'" ///
                    _col(`c4') in y "`mdl`nopref'type'" ///
                    _col(`c5') in y "$bicsup"
                * aic
                scalar `aicd' = aic`mdl1' - aic`mdl2'
                local aicfav "`mdl2type'"
                local aicnofav "`mdl1type'"
                if aic`mdl1' < aic`mdl2' {
                    local aicfav "`mdl1type'"
                    local aicnofav "`mdl2type'"
                }
                di _col(`c1') in g "AIC=" in y `f103' aic`mdl2' ///
                    _col(`c2') in g "dif=" in y `f103' `aicd' ///
                    _col(`c3') in y "`aicfav'" ///
                    _col(`c4') in y "`aicnofav'"
                * lr test
                local lrfav "`mdl1type'"
                local lrnofav "`mdl2type'"
                if lrnb_prmp < .05 {
                    local lrfav "`mdl2type'"
                    local lrnofav "`mdl1type'"
                }
                di _col(`c1') in g "LRX2=" in y `f93' lrzip_zinb ///
                    _col(`c2') in g "prob=" in y `f93' lrzip_zinbp in g ///
                    _col(`c3') in y "`lrfav'" ///
                    _col(`c4') in y "`lrnofav'" ///
                    _col(`c5') in y "p=" `f93l' lrnb_prmp
                di in g "`dash'"

            } // no zinb vs zip

        } // no zip

}  // no fit

end

// compute strength of bic difference?

capture program drop _count_bic_dif
program define _count_bic_dif

    * compute bin diff from model 1 and model 2
    global bicdif = $bic1 - $bic2
    tempname bicdabs // 0.8.0
    scalar `bicdabs' = abs($bicdif)

    * evaluate support based on bic difference
    if `bicdabs'~=. {
        global bicsup "Very strong"
        if `bicdabs'<= .0000000001 {
            global bicsup "no"
            }
        else if `bicdabs' <=2      {
            global bicsup "Weak"
        }
        else if `bicdabs' <=6      {
            global bicsup "Positive"
        }
        else if `bicdabs' <=10     {
            global bicsup "Strong"
        }
        global bicpref = 2
        global bicnopref = 1
        if $bicdif < 0 {
            global bicpref = 1
            global bicnopref = 2
        }
        if `bicdabs'< .0000000001 & `bicdabs'>-.0000000001 {
            global bicpref = 0
        }

    }
end

// ESTIMATE MODELS AND STORE RESULTS

capture program drop _count_estimate
program define _count_estimate, rclass

    syntax varlist [if] [in] [,     ///
      inflate(varlist) ///
      MAXcount(integer 9) Generate(string) replace ///
      Prm Nbreg ZIP ZINb NOConstant ///
      NOIsily ]

// estimate models & create globals with stats

    local set "`generate'"
    tempname n

    local noise ""
    if "`noisily'"=="noisily" {
        local noise "noisily"
    }

// prm

    if "`prm'"!="" {
        qui {
            local modelnm "PRM"
            local fullnm "`set'`modelnm'"
            `noise' poisson `varlist' `if' `in', `noconstant'
            estimates store `fullnm'
            scalar ll`fullnm' = e(ll) // log lik
            fitstat, bic
            scalar bicp`fullnm' = r(bic_p) // bic'
            return scalar bicp`fullnm' = r(bic_p) // 0.8.0
            scalar aic`fullnm' = r(aic) // aic
            scalar x2`fullnm' = r(lrx2) // lrx2 all b=0
            scalar x2p`fullnm' = r(lrx2_p)
            scalar bic`fullnm' =  r(bic) // bic
            if "`replace'"=="replace" {
                capture drop `fullnm'rate
                capture drop `fullnm'prgt
                capture drop `fullnm'val
                capture drop `fullnm'obeq
                capture drop `fullnm'preq
                capture drop `fullnm'prle
                capture drop `fullnm'oble
                capture drop `fullnm'absdif
                capture drop `fullnm'dif
                capture drop `fullnm'pearson
                local i = 0
                while `i'<=`maxcount' {
                    capture drop `fullnm'pr`i'
                    capture drop `fullnm'cu`i'
                    local i = `i' + 1
                }
            }
            prcounts `fullnm', plot max(`maxcount') // predicted counts
            label var `fullnm'preq "`modelnm' predicted" // predicted Pr(y)
            gen `fullnm'dif = `fullnm'obeq - `fullnm'preq // obs - predicted
            label var `fullnm'dif  "`modelnm' obs - pred"
            * 2004-10-29 add CT 5.34
            scalar `n' = e(N) // sample size
            gen `fullnm'pearson = ///
                (`n' * `fullnm'dif * `fullnm'dif) / `fullnm'preq
            label var `fullnm'pearson "`modelnm' contribution to Pearson X2"
            sum `fullnm'pearson
        } // qui
    }

// nbreg

    if "`nbreg'"!="" {
        qui {
            local modelnm "NBRM"
            local fullnm "`set'`modelnm'"
            `noise' nbreg `varlist' `if' `in', `noconstant'
            estimates store `fullnm'
            scalar ll`fullnm' = e(ll) // log lik

            scalar lrnb_prm = e(chi2_c) // lrx2 of nb vs prm
            scalar lrnb_prmp = chiprob(1, e(chi2_c))*0.5

            fitstat, bic
            scalar bicp`fullnm' = r(bic_p) // bic'
            scalar aic`fullnm' = r(aic) // aic
            scalar x2`fullnm' = r(lrx2) // lrx2 all b=0
            scalar x2p`fullnm' = r(lrx2_p)
            scalar bic`fullnm' =  r(bic) // bic

            if "`replace'"=="replace" {
                capture drop `fullnm'rate
                capture drop `fullnm'prgt
                capture drop `fullnm'val
                capture drop `fullnm'obeq
                capture drop `fullnm'preq
                capture drop `fullnm'prle
                capture drop `fullnm'oble
                capture drop `fullnm'dif
                capture drop `fullnm'absdif
                capture drop `fullnm'pearson
                local i = 0
                while `i'<=`maxcount' {
                    capture drop `fullnm'pr`i'
                    capture drop `fullnm'cu`i'
                    local i = `i' + 1
                }
            }

            prcounts `fullnm', plot max(`maxcount') // predicted counts
                label var `fullnm'preq "`modelnm' predicted"
            gen `fullnm'dif = `fullnm'obeq - `fullnm'preq
                label var `fullnm'dif  "`modelnm' obs - pred"
            * 2004-10-29 add CT 5.34
            scalar `n' = e(N) // sample size
            gen `fullnm'pearson = ///
                (`n' * `fullnm'dif * `fullnm'dif) / `fullnm'preq
            label var `fullnm'pearson "`modelnm' contribution to Pearson X2"
            sum `fullnm'pearson
        } // qui
    }

// zip

    if "`zip'"!="" {
        qui {
            local modelnm "ZIP"
            local fullnm "`set'`modelnm'"
            `noise' zip `varlist' `if' `in', ///
                inf(`inflate') vuong `noconstant'
            estimates store `fullnm'
            scalar vu`fullnm' = e(vuong) // vuong vs prm
            scalar vu`fullnm'p = 1-norm(abs(e(vuong)))
            scalar ll`fullnm' = e(ll) // log lik
            fitstat, bic
            scalar bicp`fullnm' = r(bic_p) // bic'
            scalar aic`fullnm' = r(aic) // aic
            scalar x2`fullnm' = r(lrx2) // lrx2 all b=0
            scalar x2p`fullnm' = r(lrx2_p)
            scalar bic`fullnm' =  r(bic) // bic
            if "`replace'"=="replace" {
                capture drop `fullnm'rate
                capture drop `fullnm'prgt
                capture drop `fullnm'val
                capture drop `fullnm'obeq
                capture drop `fullnm'preq
                capture drop `fullnm'prle
                capture drop `fullnm'oble
                capture drop `fullnm'all0
                capture drop `fullnm'dif
                capture drop `fullnm'absdif
                capture drop `fullnm'pearson
                local i = 0
                while `i'<=`maxcount' {
                    capture drop `fullnm'pr`i'
                    capture drop `fullnm'cu`i'
                    local i = `i' + 1
                }
            }

            prcounts `fullnm', plot max(`maxcount') // predicted counts
            label var `fullnm'preq "`modelnm' predicted"
            gen `fullnm'dif = `fullnm'obeq - `fullnm'preq
            label var `fullnm'dif  "`modelnm' obs - pred"
            * 2004-10-29 add CT 5.34
            scalar `n' = e(N) // sample size
            gen `fullnm'pearson = ///
                (`n' * `fullnm'dif * `fullnm'dif) / `fullnm'preq
            label var `fullnm'pearson "`modelnm' contribution to Pearson X2"
            qui sum `fullnm'pearson
        } // qui
    }

// zinb

    if "`zinb'"!="" {
        qui {
            local modelnm "ZINB"
            local fullnm "`set'`modelnm'"
            `noise' zinb `varlist'  `if' `in', ///
                inf(`inflate') vuong zip `noconstant'
            estimates store `fullnm'
            scalar vu`fullnm' = e(vuong) // vuong vs nbreg
            scalar vu`fullnm'p = 1-norm(abs(e(vuong)))
            scalar ll`fullnm' = e(ll) // loglik
            scalar lrzip_zinb = e(chi2_cp) // lrx2 zinb vs zip
            scalar lrzip_zinbp = chiprob(1, e(chi2_cp))*0.5
            fitstat, bic
            scalar bicp`fullnm' = r(bic_p) // bic'
            scalar aic`fullnm' = r(aic) // aic
            scalar x2`fullnm' = r(lrx2) // lrx2 all b=0
            scalar x2p`fullnm' = r(lrx2_p)
            scalar bic`fullnm' =  r(bic) // bic
            if "`replace'"=="replace" {
                capture drop `fullnm'rate
                capture drop `fullnm'prgt
                capture drop `fullnm'val
                capture drop `fullnm'obeq
                capture drop `fullnm'preq
                capture drop `fullnm'prle
                capture drop `fullnm'oble
                capture drop `fullnm'absdif
                capture drop `fullnm'dif
                capture drop `fullnm'all0
                capture drop `fullnm'pearson
                local i = 0
                while `i'<=`maxcount' {
                    capture drop `fullnm'pr`i'
                    capture drop `fullnm'cu`i'
                    local i = `i' + 1
                }
            }

            prcounts `fullnm', plot max(`maxcount') // predicted counts
            label var `fullnm'preq "`modelnm' predicted"
            gen `fullnm'dif = `fullnm'obeq - `fullnm'preq
            label var `fullnm'dif  "`modelnm' obs - pred"
            * 2004-10-29 add CT 5.34
            scalar `n' = e(N) // sample size
            gen `fullnm'pearson = ///
                (`n' * `fullnm'dif * `fullnm'dif) / `fullnm'preq
            label var `fullnm'pearson "`modelnm' contribution to Pearson X2"
            qui sum `fullnm'pearson
        } // qui
    }

end // _count_estimate

exit
* version 0.8.0 fix pr bug
* version 0.2.1 13Apr2005 add replace; trap weights
* version 0.2.0 27Feb2005 first documented version