ROSALI-SIM/Modules/ado/plus/p/prcounts.ado

*! version 1.6.8 13Apr2005
*  version 1.6.7 18Feb2005 ztp and ztnb
*  version 1.6.6 8/20/03
*  version 1.6.5 1/10/01

capture program drop prcounts
program define prcounts
    version 6.0
    tempname b alpha ai gai mu
    syntax newvarname [if] [in] [, Max(integer 9) Plot]

*=> classify each valid type of model
    local cmd = e(cmd)
    *zt 18Feb2005
    if "`cmd'"=="nbreg" | "`cmd'"=="poisson" | /*
    */ "`cmd'"=="zinb"  | "`cmd'"=="zip"     | /*
    */ "`cmd'"=="ztp"   | "`cmd'"=="ztnb"    {
    }
    else {
        di
        di in r "prcounts does not work for the last type of model estimated."
        exit
    }
    if `max' > 99 | `max' < 0 {
        di in r "max() must be in range from 0 to 99"
        exit 198
    }
    local stem "`varlist'"
    local stem = substr("`stem'", 1, 28)
    cap version 7
        if _rc!=0 { local stem = substr("`stem'", 1, 4) }
    version 6.0
    local modelis "from `cmd'"

*-> GENERATE PREDICTED RATE
    quietly predict `stem'rate `if' `in', n
    label variable `stem'rate "Predicted rate `modelis'"
    gen `mu' = `stem'rate

    *18Feb2005
    * compute conditional mu
    if  "`cmd'"=="ztp"   | "`cmd'"=="ztnb"  {
        quietly predict `stem'Crate `if' `in', cm
        label variable `stem'Crate "Predicted conditional rate `modelis'"
    }

*-> GENERATE PREDICTED FREQUENCY OF ALWAYS ZERO
    if "`cmd'"=="zip" | "`cmd'"=="zinb" {
        quietly predict `stem'all0 `if' `in', p
        label variable `stem'all0 "Pr(always 0) `modelis'"
        * predict n, n computes: exp(x*beta)*(1-all0).
        * The last term needs to be removed.
        quietly replace `mu' = `mu'/(1-`stem'all0) `if' `in'
    }

*-> TAKE CARE OF ALPHA
    *18Feb2005
    if "`cmd'"=="nbreg" | "`cmd'"=="ztnb" {
        sca `alpha' = e(alpha)
    }
    if "`cmd'"=="zinb" {
        mat `b' = e(b)
        local temp = colsof(`b')
        sca `alpha' = `b'[1, `temp']
        sca `alpha' = exp(`alpha')
    }
    *18Feb2005
    if "`cmd'"=="nbreg" | "`cmd'"=="zinb" | "`cmd'"=="ztnb" {
        sca `ai' = 1/`alpha'
        sca `gai' = exp(lngamma(`ai'))
        if `gai'==. {
            di in r "problem with alpha prevents" /*
            */ " estimation of predicted probabilities."
            exit 198
        }
    }

*-> GENERATE PREDICTED PROBABILITIES
    local i 0
    while `i' <= `max' {
        local newvar "`stem'pr`i'"
        *18Feb2005
        if "`cmd'"=="poisson" | "`cmd'"=="ztp" {
            quietly gen `newvar' = /*
            */ ((exp(-`mu'))*(`mu'^`i'))/(round(exp(lnfact(`i'))), 1) /*
            */ `if' `in'
        }
        *18Feb2005
        if "`cmd'"=="nbreg" | "`cmd'"=="ztnb" {
            quietly gen `newvar' = /*
            */ (exp(lngamma(`i'+`ai')) / /*
            */ (round(exp(lnfact(`i')),1) * exp(lngamma(`ai')))) /*
            */ * ((`ai'/(`ai'+`mu'))^`ai') * ((`mu'/(`ai'+`mu'))^`i') /*
            */ `if' `in'
        }
        if "`cmd'"=="zip" {
            quietly gen `newvar' = (1-`stem'all0)*((exp(-`mu')) /*
            */ * (`mu'^`i'))/(round(exp(lnfact(`i'))), 1) `if' `in'
            if `i'==0 {
                quietly replace `newvar' = `newvar' + `stem'all0 `if' `in'
            }
        }
        if "`cmd'"=="zinb" {
            quietly gen `newvar' = (1-`stem'all0)*(exp(lngamma(`i'+`ai')) /*
            */ / ( round(exp(lnfact(`i')),1) * exp(lngamma(`ai')) )  ) /*
            */ * ((`ai'/(`ai'+`mu'))^`ai') * ((`mu'/(`ai'+`mu'))^`i') /*
            */ `if' `in'
            if `i'==0 {
                quietly replace `newvar' = `newvar' + `stem'all0 `if' `in'
            }
        }
        label variable `newvar' "Pr(y=`i') `modelis'"
        local i = `i' + 1
    }

*spost9 zt needs to compute conditional probabilities 19feb2005
*-> GENERATE CONDITIONAL PREDICTED PROBABILITIES
if "`cmd'"=="ztp" | "`cmd'"=="ztnb" {
    local i 1
    while `i' <= `max' {
        local newvar "`stem'Cpr`i'"
        quietly gen `newvar' = `stem'pr`i'/(1-`stem'pr0)
        label variable `newvar' "Pr(y=`i'|y>0) `modelis'"
        local i = `i' + 1
    }
} /// zt

*-> GENERATE CUMULATIVE PROBABILITIES
    quietly gen `stem'cu0=`stem'pr0
    label variable `stem'cu0 "Pr(y=0) `modelis'"
    local i 1
    while `i' <= `max' {
        quietly egen `stem'cu`i' = rsum(`stem'pr0-`stem'pr`i') if `mu'~=.
        label variable `stem'cu`i' "Pr(y<=`i') `modelis'"
*spost9 zt cumulative probabilities 18Feb2005
        if "`cmd'"=="ztp" | "`cmd'"=="ztnb" {
            quietly egen `stem'Ccu`i' = rsum(`stem'Cpr1-`stem'Cpr`i') if `mu'~=.
            label variable `stem'Ccu`i' "Pr(y<=`i'|y<0) `modelis'"
        }
        local i = `i' + 1
    }

*-> GENERATE GREATER THAN VARIABLE
    quietly gen `stem'prgt = 1-`stem'cu`max' if `mu'~=.
    label variable `stem'prgt "Pr(y>`max') `modelis'"
    *18Feb2005
    if "`cmd'"=="ztp" | "`cmd'"=="ztnb" {
        quietly gen `stem'Cprgt = 1-`stem'Ccu`max' if `mu'~=.
        label variable `stem'Cprgt "Pr(y>`max'|y>0) `modelis'"
    }

*-> IF PLOT OPTION SPECIFIED
    if "`plot'"=="plot" {
        quietly gen `stem'val = .
            label variable `stem'val "Count"
        quietly gen `stem'obeq = .
            label variable `stem'obeq "Observed Pr(y=k) `modelis'"
        quietly gen `stem'preq = .
            label variable `stem'preq "Predicted Pr(y=k) `modelis'"
        quietly gen `stem'oble = .
            label variable `stem'oble "Observed Pr(y<=k) `modelis'"
        quietly gen `stem'prle = .
            label variable `stem'prle "Predicted Pr(y<=k) `modelis'"
    *18Feb2005
        if "`cmd'"=="ztp" | "`cmd'"=="ztnb" {
            quietly gen `stem'Cpreq = .
            label variable `stem'Cpreq "Predicted Pr(y=k|y>0) `modelis'"
            quietly gen `stem'Cprle = .
            label variable `stem'Cprle "Predicted Pr(y<=k|y>0) `modelis'"
        }
        ** bug fix -- makes sure observed probabilities are
        **      computed on estimation sample
        if "`if'" == "" {
            local if "if e(sample)==1"
        }
        else {
            local if "`if' & e(sample)==1"
        }

        local i 0
        while `i' <= `max' {
            quietly {
                local obs = `i' + 1
                replace `stem'val = `i' in `obs'
                tempvar count1 count2
                * 1 if outcomes equal to i
                gen `count1' = (`e(depvar)'==`i') `if' `in'
                sum `count1' `if' `in'
                replace `stem'obeq = r(mean) in `obs'
                * 1 if outcome lt i
                gen `count2' = (`e(depvar)'<=`i') `if' `in'
                sum `count2' `if' `in'
                replace `stem'oble = r(mean) in `obs'
                * compute average predicted prob
                sum `stem'pr`i' `if' `in'
                replace `stem'preq = r(mean) in `obs'
                * compute average cumulative predicted prob
                sum `stem'cu`i' `if' `in'
                replace `stem'prle = r(mean) in `obs'

                *zt only compute if count > 0 18Feb2005*ZT
                if ("`cmd'"=="ztp" | "`cmd'"=="ztnb") & `i'>0 {
                    sum `stem'Cpr`i' `if' `in'
                    replace `stem'Cpreq = r(mean) in `obs'
                    sum `stem'Ccu`i' `if' `in'
                    replace `stem'Cprle = r(mean) in `obs'
                }
            }
            local i = `i' + 1
        }
    }

end