*! 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