capture log close
set matsize 500
set more off
log using rm4ch8, text replace

***
*** RM4CLDV Ch 8: Count Outcomes - using SPost - verified 3/9/01
***

*** Page 220: Figure 8.2 - fit a Poisson distribution to art
use couart2, clear
poisson art, nolog
prcounts null, plot max(9)
label var nullobeq "Observed"
label var nullpreq "Predicted: Poisson Dist"
graph nullobeq nullpreq nullval, c(ll) gap(3) l2(Probability) /*
*/ yscale(0,.4) ylabel(0,.1,.2,.3,.4) xlabel(0,1,2,3,4,5,6,7,8,9)
translate @Graph rm4ch8fig2.wmf, replace
drop null*

*** Page 227: Table 8.1
gen lnart = log(art + .5)
summarize

*** Page 228: Table 8.2 - OLS with lnart
regress lnart fem mar kid5 phd ment
listcoef, help

*** Page 228: Table 8.2 - Poison regression
poisson art fem mar kid5 phd ment, nolog
listcoef, help

*** Pages 228-229: Poisson regression - predicted values
prvalue, x(fem 1) rest(mean)
prvalue, x(fem 0) rest(mean)

*** Page 228: Table 8.2 - NB regression
nbreg art fem mar kid5 phd ment, nolog
listcoef, help

*** Page 229: Figure 8.4 - predictions from PRM and NBRM
quietly poisson art fem mar kid5 phd ment, nolog
prcounts preg, plot max(9)
sum preg*
quietly nbreg art fem mar kid5 phd ment, nolog
prcounts nbrg, plot max(9)
sum nbrg*
label var pregpreq "Predicted: PRM"
label var pregobeq "Observed"
label var nbrgpreq "Predicted: NBRM"
graph pregobeq pregpreq nbrgpreq pregval, c(ll) gap(3) l2(Probability) /*
*/ yscale(0,.4) ylabel(0,.1,.2,.3,.4) xlabel(0,1,2,3,4,5,6,7,8,9)
translate @Graph rm4ch8fig4.wmf, replace

*** Page 237: LR test of dispersion
** -lrtest- will not work, but see output from nbreg. Or, try this:
quietly poisson art fem mar kid5 phd ment, nolog
scalar llp = e(ll)
quietly nbreg art fem mar kid5 phd ment, nolog
scalar lln = e(ll)
di "LR test of dispersion: " -2*(llp-lln)

*** Page 238: plot of prob of 0
quietly nbreg art fem mar kid5 phd ment, nolog
prgen ment, rest(mean) from(0) to(50) gen(nb) n(20)
save tmpnbpr, replace
use couart2,clear
quietly poisson art fem mar kid5 phd ment
prgen ment, rest(mean) from(0) to(50) gen(psn) n(20)
save tmppsnpr, replace
merge using tmpnbpr
label var psnp0 "Pr(0) for PRM"
label var nbp0 "Pr(0) for NBRM"
graph psnp0 nbp0 nbx, c(ll) gap(3) l2(Probability) /*
*/ yscale(0,.4) ylabel(0,.1,.2,.3,.4) xlabel(0,1,2,3,4,5,6,7,8,9)
translate @Graph rm4ch8fig8.wmf, replace

*** Page 246: Table 8.3 - zip and zinb models with Vuong statistics
zip art fem mar kid5 phd ment, inf(fem mar kid5 phd ment) vuong nolog
listcoef, help
scalar llzip = e(ll)
zinb art fem mar kid5 phd ment, inf(fem mar kid5 phd ment) vuong nolog
listcoef, help
scalar llzinb = e(ll)

*** Page 247: Comparing models
di "LR between zip and zinb: " -2*(llzip-llzinb)

*** Page 248: Figure 8.9 - from PRM, NBRM, ZIP and ZINB
use couart2, clear
quietly poisson art fem mar kid5 phd ment, nolog
prcounts prm, plot max(9)
label var prmpreq "Predicted: PRM"
label var prmobeq "Observed"
quietly nbreg art fem mar kid5 phd ment, nolog
prcounts nbrm, plot max(9)
label var nbrmpreq "Predicted: NBRM"
quietly zip art fem mar kid5 phd ment, /*
*/ inf(fem mar kid5 phd ment) vuong nolog
prcounts zip, plot max(9)
label var zippreq "Predicted: ZIP"
quietly zinb art fem mar kid5 phd ment, /*
*/ inf(fem mar kid5 phd ment) vuong nolog
prcounts zinb, plot max(9)
label var zinbpreq "Predicted: ZINB"
** create deviations
gen obs = prmobeq
gen dprm = obs - prmpreq
label var dprm "PRM"
gen dnbrm = obs - nbrmpreq
label var dnbrm "NBRM"
gen dzip = obs - zippreq
label var dzip "ZIP"
gen dzinb = obs - zinbpreq
label var dzinb "ZINB"
graph dprm dnbrm dzip dzinb prmval, c(llll) gap(3) /*
*/ l2(Observed-Predicted) ylabel(-.1,0,.1) /*
*/ xlabel(0,1,2,3,4,5,6,7,8) yline(0)
translate @Graph rm4ch8fig9.wmf, replace

log close