* Alternate Count Data Model Chapter
drop _all

set obs 21
scalar mu     = 4
scalar alpha1 = 0.5
scalar alpha2 = 1.5
scalar p1     = 1/(1+alpha1*mu)
scalar p2     = 1/(1+alpha2*mu)
gen y = _n-1
gen double fy0p0 = exp(-mu)*mu^y/exp(lngamma(y+1))
gen double fy0p5 = exp(lngamma(y+1/alpha1)-lngamma(y+1)-lngamma(1/alpha1))*p1^(1/alpha1)*(1-p1)^y
gen double fy1p5 = exp(lngamma(y+1/alpha2)-lngamma(y+1)-lngamma(1/alpha2))*p2^(1/alpha2)*(1-p2)^y
label var fy0p0 "NB-alpha=0.0"
label var fy0p5 "NB-alpha=0.5"
label var fy1p5 "NB-alpha=1.5"
twoway (scatter fy0p0 y, c(l)) (scatter fy0p5 y, c(l)) (scatter fy1p5 y, c(l)), title("Mean=4: Count Data Models") 
graph export nbplot.eps, replace

drop _all
set seed 234987
set obs 50000
generate x1 = invnormal(uniform())
generate x2 = invnormal(uniform())
generate xb = -1.0 + .5*x1 + 1.5*x2
genpoisson yp, xbeta(xb)

poisson yp x1 x2

poisson yp x1 x2 if y>0, nolog

ztp yp x1 x2 if y>0, nolog

tab yp

zip yp x1 x2, inflate(x1 x2) nolog

hplogit yp x1 x2, nolog

gen yp0 = yp > 0
logit yp0 x1 x2, nolog

ztp yp x1 x2 if yp>0, nolog

hplogit yp x1 x2, nolog
predict xb1, eq(#1)       /* logit prediction                  */
predict xb2, eq(#2)       /* zero-truncated Poisson prediction */
gen mu1 = 1/(1+exp(-xb1)) /* predicted Pr(yp=0)                */
gen mu2 = exp(xb2)        /* predicted ztp count               */

list yp0 mu1 yp mu2 in 1000/1020

use medpar, clear
nbreg los hmo white type2 type3 died, nolog

gnbreg los hmo white type2 type3 died, lnalpha(type2 type3 died) nolog

glm los hmo white type2 type3 died, family(nb .4339959) notable nolog

glm los hmo white type2 type3 died, family(poisson) notable nolog

gnbreg los hmo white type2 type3, lnalpha(died) irr nolog

gnbreg 

gpoisson los hmo white type2 type3, irr nolog

use doll, clear
glm deaths smokes a2-a5, lnoffset(pyears) family(poi) irls nolog 
predict double pearson, pearson
predict double yhatp

gen lpy = ln(pyears)
gpoisson deaths smokes a2-a5, offset(lpy) irr nolog 

local delta = `e(delta)'
predict double yhat
gen double ehat = (deaths-yhat)/sqrt(yhat/(1-`delta'^2))
twoway (scatter ehat yhat, yline(0)) (scatter pearson yhatp, yline(0))

use medparc, clear
describe

tab cenx
 
list los cenx if cenx == -1

tab los if los > 62

tab los if los > 20 & los < 30

cpoisson los hmo white type2 type3, cen(cenx) nolog

censornb los hmo white type2 type3, cen(cen1) nolog
censornb los hmo white type2 type3, cen(cenx) nolog

censornb los hmo white type2 type3, cen(cen1) nolog

censornb los hmo white type2 type3, cen(cenx) cluster(provnum) irr nolog