capture log close
set more off
estimates clear
local job cda_notes_ch7_cou
local date "18Apr2006"

log using cda_notes_ch7_cou, replace
version 9.1
set scheme s2manual
set scheme lean1

//  cda_notes_ch7_cou: ICPSR CDA Chapter 7 - count models
//  18Apr2006 jsl

//  poisson pdf

    use couart2, clear

*   1: fit pdf
    poisson art, nolog
*   2: compute expected values
    prcounts null, plot
        label var nullobeq "Observed Probability"
        label var nullpreq "Predicted Poisson"
*   3: plot predictions
    twoway connected nullobeq nullpreq nullval,         ///
        ytitle("Pr(y=k)") xtitle("# of articles")       ///
        ylabel(0(.1).4, grid gmax gmin gstyle(dot))     ///
        xlabel(0(1)9, nogrid)                           ///
        caption("(`job' `date')", size(tiny))
    graph export `job'_psnpdf.emf, replace

//  prm

*   4: fit prm
    sum art fem mar kid5 phd ment
    poisson art fem mar kid5 phd ment
*   5: discrete change
    prchange, rest(mean)
*   6: factor change
    listcoef, help
*   7: list data for cases 1 and 2
    list fem mar kid5 phd ment in 1/2, nolabel clean
*   8: predictions for case 1
    prvalue, x(fem=0 mar=1 kid5=0 phd=2.52 ment=7)
*   9: predictions for case 2
    prvalue, x(fem=1 mar=0 kid5=0 phd=2.05 ment=6)
*   10: compute predictions to plot
    prcounts prm, plot
        label var prmobeq "Observed"
        label var prmpreq "PRM"
*   11: examine predictions
    desc prmrate prmpr0 prmpr1 prmval prmpreq prmobeq
    sum  prmrate prmpr0 prmpr1 prmval prmpreq prmobeq
    list prmrate prmpr0-prmpr4 in 1/2 // see prvalue above
*   12: observed frequences
    tab art
*   13: list predictions
    list art prmrate prmpr0 prmpr1 prmval prmpreq prmobeq ///
        in 1/12, nodisplay clean
*   14: plot observed and predicted
    twoway ///
        (connected prmobeq prmval, msymbol(circle)                  ///
            msize(medium) mcolor(black) clcolor(black) clpat(dot))  ///
        (connected prmpreq prmval, msymbol(square_hollow)           ///
            msize(small) mcolor(black) clcolor(black) clpat(dot))   ///
        , ytitle("Pr(y=k)") xtitle("# of articles")                 ///
          ylabel(0(.1).4, grid gmax gmin gstyle(dot))               ///
          xlabel(0(1)9, nogrid)                                     ///
          caption("(`job' `date')", size(tiny))
    graph export `job'_prmprobs.emf, replace

//  negative binomial model

*   15: estimate model
    nbreg art fem mar kid5 phd ment, nolog
*   16: compute predictions over values of ment
    prgen ment, from(0) to(50) generate(nb) rest(mean)
        label var nbp0 "Pr(0) for Nebgin Regression"
*   17: compute corresponding values from PRM
    qui poisson art fem mar kid5 phd ment, nolog
    prgen ment, from(0) to(50) generate(psn) rest(mean)
        label var psnp0 "Pr(0) for Poisson Regression"
*   18: compare prm and nbrm
    graph twoway connected psnp0 nbp0 nbx,                  ///
        ytitle("Pr(Zero Articles)") ylabel(0(.1).4, grid)   ///
        msymbol(Th Sh) caption("(`job' `date')", size(tiny))
    graph export `job'_prmnbrm.emf, replace

//  zero inflated models

*   19: zip
    zip art fem mar kid5 phd ment, inflate(fem mar kid5 phd ment) nolog
*   20: zinb
    zinb art fem mar kid5 phd ment, inf(fem mar kid5 phd ment) nolog
*   21: coefficients from zinb
    listcoef, help
*   22: predictions
    prvalue

//  comparing models the easy way

    use couart2, clear
*   23: compare base models
    countfit art fem mar kid5 phd ment, gen(Base_)   ///
        inflate(fem mar kid5 phd ment) maxcount(6)  ///
        prm nbreg zip zinb nodash
    graph export `job'_countfit.emf, replace

log close