version 8
capture log close
set matsize 500
set more off
set scheme spost // only used for generating graphs for publication
log using st8ch7, text replace

// *
// * RM4STATA Ch 7: Models for Count Outcomes - 5/26/2003
// *

// * section 7.1: the poisson distribution

* generate the graph from this section
clear
set obs 25
gen ya = .8
poisson ya, nolog
prcounts pya, plot max(20)
gen yb = 1.5
poisson yb, nolog
prcounts pyb, plot max(20)
gen yc = 2.9
poisson yc, nolog
prcounts pyc, plot max(20)
gen yd = 10.5
poisson yd, nolog
prcounts pyd, plot max(20)
label var pyapreq "mu=0.8"
label var pybpreq "mu=1.5"
label var pycpreq "mu=2.9"
label var pydpreq "mu=10.5"
label var pyaval "y=# of Events"
graph twoway connected pyapreq pybpreq pycpreq pydpreq pyaval, ///
    ytitle("Probability") ylabel(0(.1).5) xlabel(0(2)20)
graph  export 07poisson.eps, replace

// * section 7.1.1: fitting the poisson distribution with -poisson-

use couart2, clear
describe
summarize
poisson art, nolog

// * section 7.1.3: plotting the poisson distribution with -prcounts- 
prcounts psn, plot max(9)
label var psnobeq "Observed Proportion"
label var psnpreq "Poisson Prediction"
label var psnval "# of Articles"
list psnval psnobeq psnpreq in 1/10
graph twoway connected psnobeq psnpreq psnval, ///
    ytitle("Probability") ylabel(0(.1).4) ///
    xlabel(0(1)9)   ///
    ysize(2.7051) xsize(4.0421)
graph export 07psnpred.eps, replace

// * section 7.2.1: prm - estimation

* non-integer outcomes - note warning!
use couart2, clear
gen artx = art + .3
poisson artx fem mar kid5 phd ment, nolog

* if and in conditions
poisson art mar kid5 phd ment if fem==1, nolog

// * section 7.2.2: prm - example

use couart2, clear
poisson art fem mar kid5 phd ment, nolog

// * section 7.2.3: prm - interpretation using the rate mu

* factor changes
listcoef fem ment, help

* percent change
listcoef fem ment, percent help

* marginal change with -prchange-
prchange phd ment, rest(mean)

* marginal change with -mfx compute-
mfx compute

* discrete change
prchange fem ment, rest(mean)
prchange kid5, uncentered x(kid5=1)
prchange kid5, uncentered x(kid5=1) delta(2)

// * section 7.2.4: prm - interpretation using predicted probabilities

* -prvalue-  
* single women without children
quietly prvalue, x(mar=0 fem=1 kid5=0) rest(mean) save
* compared to married women without children
prvalue, x(mar=1 fem=1 kid5=0) rest(mean) dif
* effects of number of children
prvalue, x(mar=1 fem=1 kid5=0) rest(mean) brief
prvalue, x(mar=1 fem=1 kid5=1) rest(mean) brief
prvalue, x(mar=1 fem=1 kid5=2) rest(mean) brief
prvalue, x(mar=1 fem=1 kid5=3) rest(mean) brief

* -prgen-
* for women
prgen kid5, x(fem=1 mar=1) rest(mean) from(0) to(3) gen(fprm) n(4)
* for men
prgen kid5, x(fem=0 mar=1) rest(mean) from(0) to(3) gen(mprm) n(4)
label var fprmp0 "Married Women"
label var mprmp0 "Married Men"
label var mprmx  "Number of Children"
* graph predictions
graph twoway connected fprmp0 mprmp0 mprmx, ///
    ylabel(0(.1).4) yline(.1 .2 .3) ///
    xlabel(0(1)3) ///
    ytitle("Probability of No Articles") ///
    ysize(2.5051) xsize(4.0421)
graph export 07prmprob0.eps, replace

*  list the predictions
list fprmp0 mprmp0 mprmx in 1/4

* -prcounts-
* for univariate poisson distribution
poisson art, nolog
prcounts psn, plot max(9)
label var psnpreq "Univariate Poisson Dist."
poisson art fem mar kid5 phd ment, nolog
prcounts prm, plot max(9)
label var prmpreq "PRM"
label var prmobeq "Observed"
graph twoway connected prmobeq psnpreq prmpreq prmval, ///
    ytitle("Probability of Count") ylabel(0(.1).4) ///
    xlabel(0(1)9)   ///
    ysize(2.7051) xsize(4.0413)
graph export 07prmpred.eps, replace

// * section 7.2.5: prm - exposure time

* simulate data to illustrate exposure
* construct random integers
gen profage = round((1+10*uniform()),1)
label var profage "Professional Age"
gen lnage = ln(profage)
label var lnage "Log of Professional Age"
* let total articles equal artiles*exposure time
gen totalarts = art*profage
label var totalarts "Total Articles during Career"

* PRM without including time
poisson totalarts fem mar kid5 phd ment, nolog

* PRM exposure time
poisson totalarts fem mar kid5 phd ment, nolog exposure(profage)

* PRM with exposure time as an independent variable with a constraint
constraint define 1 lnage=1
poisson totalarts fem mar kid5 phd ment lnage, nolog constraint(1)

* PRM with offset()
poisson totalarts fem mar kid5 phd ment, nolog offset(lnage)

// * section 7.3.2: nbrm - estimation

use couart2, clear
nbreg art fem mar kid5 phd ment, nolog

* combining poisson and nbreg
poisson art fem mar kid5 phd ment, nolog
estimates store PRM
nbreg art fem mar kid5 phd ment, nolog
estimates store NBRM
estimates table PRM NBRM, b(%9.3f) t label varwidth(32) drop(lnalpha:_cons) stats(alpha N)

// * section 7.3.4: nbrm - interpretation using the rate mu

nbreg art fem mar kid5 phd ment, nolog
listcoef fem ment, help
listcoef fem ment, help percent

// * section 7.3.5: nbrm - interpretation using predicted probabilities

* -prvalue- for prm and nbrm
quietly poisson art fem mar kid5 phd ment
prvalue
quietly nbreg art fem mar kid5 phd ment
prvalue

* pr(y=0) for prm vs nbrm
quietly nbreg art fem mar kid5 phd ment, nolog
prgen ment, rest(mean) f(0) t(50) gen(nb) n(20)
quietly poisson art fem mar kid5 phd ment
prgen ment, rest(mean) f(0) t(50) gen(psn) n(20)
label var psnp0 "Pr(0) for PRM"
label var nbp0 "Pr(0) for NBRM"
graph twoway connected psnp0 nbp0 nbx, ///
    ylabel(0(.1).4) yline(.1 .2 .3) ///
    ytitle("Probability of a Zero Count") ///
    ysize(2.7051) xsize(4.0421)
graph export 07nbprm0.eps, replace

// * section 7.4.2: zip/zinb - example of estimating the zip and zinb models

zip art fem mar kid5 phd ment, inf(fem mar kid5 phd ment) nolog
estimates store ZIP
zinb art fem mar kid5 phd ment, inf(fem mar kid5 phd ment) nolog
estimates store ZINB
estimates table ZIP ZINB, b(%9.3f) t label varwidth(35)

// * section 7.4.3: zip/zinb - interpretation of coefficients

zip art fem mar kid5 phd ment, inf(fem mar kid5 phd ment) nolog
listcoef, help
zinb art fem mar kid5 phd ment, inf(fem mar kid5 phd ment) nolog
listcoef, help

// * section 7.4.4: zip/zinb - interpretation of predicted probabilities

* -prvalue-
zinb art fem mar kid5 phd ment, inf(fem mar kid5 phd ment) nolog
quietly prvalue, x(fem=0 mar=1 kid5=3 phd=3 ment=10) save
prvalue, x(fem=1 mar=1 kid5=3 phd=1 ment=0) dif

* -prgen-
prgen ment, rest(mean) f(0) t(20) gen(zinb) n(21)
gen zinbnb0 = zinbp0 - zinball0
label var zinbp0 "0s from Both Equations"
label var zinball0 "0s from Binary Equation"
label var zinbnb0  "0s from Count Equation"
label var zinbx"Mentor's publications"
graph twoway connected zinball0 zinbnb0 zinbp0 zinbx, ///
    xlabel(0(5)20) ylabel(0(.1).7) ///
    ytitle(Probability of zero) msymbol(Oh Sh O) ///
    ysize(2.6541) xsize(3.9678) 
graph export 07zinbpr0.eps, replace

// * section 7.5.1: comparing mean probabilities

* estimate various models and save predictions
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"
* plot deviations
graph twoway connected dprm dnbrm dzip dzinb prmval, ///
    ytitle(Observed-Predicted) ylabel(-.10(.05).10) ///
    xlabel(0(1)9) msymbol(Oh Sh O S)  ///
    ysize(2.7051) xsize(4.0413)
graph export 07compare.eps, replace

// * section 6.5.2: tests to compare count models

* LR test for -zip- and -zinb-
use couart2, clear
nbreg art fem mar kid5 phd ment, nolog
quietly zip art fem mar kid5 phd ment, inf(fem mar kid5 phd ment) vuong nolog
scalar llzip = e(ll)
quietly zinb art fem mar kid5 phd ment, inf(fem mar kid5 phd ment) vuong nolog
scalar llzinb = e(ll)
scalar lr = -2*(llzip-llzinb)
scalar pvalue = chiprob(1,lr)/2
scalar lnalpha = -.9763565
if (lnalpha < -20) { 
    scalar pvalue= 1 
}
di as text "Likelihood-ratio test comparing ZIP to ZINB: " as res %8.3f ///
    lr as text " Prob>=" as res %5.3f pvalue

* vuong test of non-nested models
zip art fem mar kid5 phd ment, inf(fem mar kid5 phd ment) vuong nolog
listcoef, help
zinb art fem mar kid5 phd ment, inf(fem mar kid5 phd ment) vuong nolog
listcoef, help

log close