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

log using cda_notes_ch3_bin, replace
version 9.1
set scheme s2manual

//  cda_notes_ch3_bin: ICPSR CDA Chapter 3 - binary regressions
//  18Apr2006 jsl

//  LPM

    use binlfp2, clear
*   1: examine data
    summarize lfp k5 k618 age wc hc lwg inc
    desc lfp k5 k618 age wc hc lwg inc
*   2: estimate lpm
    regress lfp k5 k618 age wc hc lwg inc
*   3: predictions out of range
    prvalue, x(age=35 k5=4 wc=0 hc=0) rest(mean)

//  identifcation

    use science2, clear
*   4: change scale by factor of 100
    gen job100 = 100*job
        label var job100 "job*100"
    sum job job100
*   5: regression on job
    regress job pub1 phd female
        est store regjob
*   6: regression on job*100
    regress job100 pub1 phd female
        est store regjob100
*   7: compare results
    estimates table regjob regjob100, b(%9.4f) t(%6.2f)

//  comparing logit and probit

    use binlfp2, clear
*   8: estimate and save logit
    logit lfp k5 k618 age wc hc lwg inc, nolog
        est store binlgt1
*   9: estimate and save probit
    probit lfp k5 k618 age wc hc lwg inc, nolog
        est store binpbt1
*   10: compare logit and probit estimates
    estimates table bin*, stats(N chi2 bic) b(%9.3f) t(%6.2f)

//  comparing predicted probabilities

    use binlfp2, clear
*   11: lpm
    regress lfp k5-inc
        predict prlpm
        label var prlpm "LPM: predicted probability"
*   12: logit
    logit lfp k5-inc, nolog
        predict prlogit
        label var prlogit "Logit: predicted probability"
*   13: probit
    probit lfp k5-inc, nolog
        predict prprobit
        label var prprobit "Probit: predicted probability"
*   14: compare summary statistics of predictions
    sum pr*
*   15: plot logit against probit predictions
    twoway scatter prlogit prprobit, ysize(1) xsize(1)      ///
        msymbol(Oh) xlabel(0(.1)1,grid) ylabel(0(.1)1,grid) ///
        caption("(`job' `date')", size(tiny))
    graph export `job'_prlgtpbt.emf, replace
*   16: examine distribution of logit predictions
    dotplot prlogit, caption("(`job' `date')", size(tiny))
        graph export `job'_prdotplot.emf, replace

//  ideal types

*   17: base model
    logit lfp k5-inc, nolog
*   18: prediction at average
    prvalue, rest(mean)
*   19: young, mother with little education in family
    prvalue, x(age=35 k5=4 wc=0 hc=0) rest(mean)
*   20: young, educated family without children
    prvalue, x(age=35 k5=0 wc=1 hc=1) rest(mean)
*   21: older, educated women in high income family
    prvalue, x(age=55 k5=0 k618=0 wc=1 hc=1 inc=80) rest(mean)

//  table of probabilities

    probit lfp k5-inc, nolog
*   22: predictions with prtab
    prtab k5 wc, rest(mean) novarlbl

//  plotting probabilities

*   23: predictions while age changes; wc=0 and wc=1
    prgen age, x(wc=1) rest(mean) f(30) t(60) gen(wc1) gap(5) ci
        label var wc1p1 "Attended college"
    prgen age, x(wc=0) rest(mean) f(30) t(60) gen(wc0) gap(5) ci
        label var wc0p1 "Did not attend college"
*   24: graph predictions (with ci) by age by wc
    graph twoway ///
        (rarea wc1p1ub wc1p1lb wc1x, color(gs10))            ///
        (rarea wc0p1ub wc0p1lb wc1x, color(gs10))            ///
        (connected wc1p1 wc1x, lpattern(dash) msize(zero))   ///
        (connected wc0p1 wc1x, lpattern(solid) msize(zero)), ///
        xtitle("Age") legend(on order(3 4))                  ///
        ylabel(0(.2)1,grid) ytitle("Pr(In Labor Force)")     ///
        ysize(4) xsize(6) caption("(`job' `date')", size(tiny))
    graph export `job'_prage.emf, replace

*   25: predictions by inc for various age levels
    prgen inc, from(0) to(100) ncases(13) x(age=30) rest(mean) gen(p30)
        label var p30p1 "Age 30"
    prgen inc, from(0) to(100) ncases(13) x(age=40) rest(mean) gen(p40)
        label var p40p1 "Age 40"
    prgen inc, from(0) to(100) ncases(13) x(age=50) rest(mean) gen(p50)
        label var p50p1 "Age 50"
    prgen inc, from(0) to(100) ncases(13) x(age=60) rest(mean) gen(p60)
        label var p60p1 "Age 60"
*   26: examine predictions created by prgen
    list p30x p30p0 p30p1 p40x p40p1 if p30x<.
*   27: plot predictions
    scatter p30p1 p60x, msymbol(Oh) c(l)          /// at age 30
         ytitle("Pr(In LF)") ylabel(0(.2)1.,grid) ///
         xtitle("Income") xlabel(0(20)100)        ///
     || scatter p40p1 p60x, msymbol(Sh) c(l)      /// at age 40
     || scatter p50p1 p60x, msymbol(Th) c(l)      /// at age 50
     || scatter p60p1 p60x, msymbol(Dh) c(l)      /// at age 60
    ysize(4) xsize(6) caption("(`job' `date')", size(tiny))
    graph export `job'_princ.emf, replace

//  discrete change

    probit lfp k5 k618 age wc hc lwg inc, nolog
*   28: discrete change for wc
    prchange wc
*   28: compute same change with confidence intervals
    quiet prvalue, x(wc=0) rest(mean) save
          prvalue, x(wc=1) rest(mean) dif

log close