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

log using cda_notes_ch4_fit, replace
version 9.1
set scheme s2manual

//  cda_notes_ch4_fit: ICPSR CDA Chapter 4 - testing and fit
//  18Apr2006 jsl

//  testing

    use binlfp2, clear
*   1: simple z test
    logit lfp k5 k618 age wc hc lwg inc, nolog
*   2: wald test
    test k5
    test wc hc
    test wc=hc
*   3: LR test
   *a: full model
    logit lfp k5 k618 age wc hc lwg inc, nolog
        estimates store full
   *b: k5=0
    logit lfp k618 age wc hc lwg inc, nolog
        estimates store dropk5
    lrtest full dropk5
   *c: wc=hc=0
    logit lfp k5 k618 age lwg inc, nolog
        estimates store dropwchc
    lrtest full dropwchc
   *d: all coefficients zero
    logit lfp, nolog
        estimates store constant
    lrtest full constant

//  residuals and influence

*   4: create index number
    sort inc, stable
    generate index = _n
        label var index "Observation Number"
*   5: estimate model
    logit lfp k5 k618 age wc hc lwg inc, nolog
*   6: standardized residuals
    predict jslrstd, rs
        label var jslrstd "Standardized residual"
*   7: index plot of residuals
    graph twoway scatter jslrstd index, msymbol(Oh) mcolor(black) ///
        xtitle("Observation Number") xlabel(0(200)800)            ///
        ylabel(-4(2)4) yline(0, lpattern(solid))                  ///
        yline(2 -2, lpattern(dot)) caption("(`job' `date')", size(tiny))
    graph export `job'_rstd.emf, replace
*   8: graph labeling points with index number
    graph twoway scatter jslrstd index, msymbol(none) ///
        mlabel(index) mcolor(black) mlabposition(0) ///
        xtitle("Observation Number") xlabel(0(200)800) ///
        ylabel(-4(2)4) yline(0, lpattern(solid)) ///
        yline(2 -2, lpattern(dot)) caption("(`job' `date')", size(tiny))
    graph export `job'_rstdnum.emf, replace
*   9: list large residuals
    sort jslrstd
    list jslrstd lfp k5 wc inc lwg index if jslrstd>2.25 | jslrstd<-2.25
*   10: label residuals with value of k5
    graph twoway scatter jslrstd index, msymbol(none) ///
        mlabel(k5) mcolor(black) mlabposition(0) ///
        xtitle("Observation Number") xlabel(0(200)800) ///
        ylabel(-4(2)4) yline(0, lpattern(solid)) ///
        yline(2 -2, lpattern(dot)) caption("(`job' `date')", size(tiny))
    graph export `job'_rstdk5.emf, replace
*   11: list only large residuals with value of k5
    graph twoway scatter jslrstd index if k5!=0, msymbol(none) ///
        mlabel(k5) mcolor(black) mlabposition(0) ///
        xtitle("Observation Number") xlabel(0(200)800) ///
        ylabel(-4(2)4) yline(0, lpattern(solid)) ///
        yline(2 -2, lpattern(dot)) caption("(`job' `date')", size(tiny))
    graph export `job'_rstdk5not0.emf, replace
*   12: try using dummies for k5
    foreach k in 1 2 3 {
        gen dumk5`k' = k5==`k'
        label var dumk5`k' "k5==`k'"
    }
    logit lfp dumk5* k618 age wc hc lwg inc, nolog
    predict jslrstd2, rs
        label var jslrstd2 "Std resid for dummy k5"
    sort jslrstd2
*   13: graph new residuals
graph twoway scatter jslrstd2 index, msymbol(none) ///
        mlabel(k5) mcolor(black) mlabposition(0) ///
        xtitle("Observation Number") xlabel(0(200)800) ///
        ylabel(-4(2)4) yline(0, lpattern(solid)) ///
        yline(2 -2, lpattern(dot)) caption("(`job' `date')", size(tiny))
    graph export `job'_rstdk5dum.emf, replace

//  influence

*   14: compute influence
    predict jslcook, dbeta
        label var jslcook "Cook's Statistic"
*   15: plot influence statistics
    graph twoway scatter jslcook index, ///
        mlabel(index) msymbol(none) mlabposition(0) ///
        xtitle("Observation Number") xlabel(0(200)800) ///
        ylabel(0(.1).3) yline(.1 .2, lpattern(dot)) ///
        caption("(`job' `date')", size(tiny))
    graph export `job'_cook.emf, replace

//  compare LRM models

*   16: base model
    regress lfp k5 k618 age wc hc lwg inc
        estimates store fitlrm1, title(LRM-base)
        quietly fitstat, save
*   17: model adding age squared and removing other variables
    gen age2 = age*age
        label var age2 "Age-squared"
    regress lfp k5 age age2 wc inc
        estimates store fitlrm2, title(LRM-age2 w/o hc)
*   18: compare estimates
    estimates table fitlrm*, b(%8.3f) ///
       stats(N ll chi2 df_m bic aic) star
*   19: compare fit statistics
    fitstat, dif

//  compare logits with pseudo-R2's based on LRM

*   20: base model
    logit lfp k5 k618 age wc hc lwg inc, nolog
        estimates store fitbrm1, title(BRM-base)
        quietly fitstat, save
*   21: add age squared and remove others
    logit lfp k5 age age2 wc inc, nolog
        estimates store fitbrm2, title(BRM-age2 w/o hc)
*   22: compare estimates
    estimates table fitbrm*, b(%8.3f) ///
        stats(N ll chi2 df_m bic aid) star
*   23: compare fit statistics
    fitstat, dif

//  pseudo-R2s comparing observed and predicted values

*   24: estimate model and predict probabilities
    qui logit lfp k5 k618 age wc hc lwg inc, nolog
    predict prlfp
        label var prlfp "Estimated Pr(lfp)"
        sum prlfp
*   25: are predictions above .5?
    gen prlfp1 = prlfp>.5
        label var prlfp1 "Predict lfp=1"
        label def prlfp1 0 "PredinLF" 1 "PredNotLF"
        label val prlfp1 prlfp1
*   26: compute count R2
    tab lfp prlfp1, row col cell
    di "Count R2: " 100*(180+342)/(753)
    di "Adjusted Count R2: " 100*(180+342-428)/(753-428)

log close