You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
878 lines
31 KiB
Plaintext
878 lines
31 KiB
Plaintext
9 months ago
|
*! 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
|