capture log close
set more off
estimates clear
set scheme s2color
matrix drop _all

log using sug_groups_tenure, replace
version 9.1

//  sug_groups_tenure.do - stata user's group - july 2006
//  logit example of gender differences in tenure
//  Scott Long - 20Jul2006

// #1
//  descriptive statistics

spex tenure01, clear
vardesc tenure female year yearsq select articles prestige presthi

// #2
//  estimation

logit tenure female year yearsq select articles presthi, nolog
listcoef articles female, help

// #3
//  preditions with predict using pseudo-observations

set obs 2946
replace female   =  1    if _n==2946
replace articles =  0    if _n==2946
replace year     =  7    if _n==2946
replace yearsq   = 49    if _n==2946 // year squared
replace select   =  4.97 if _n==2946 // mean level
replace presthi  =  0.05 if _n==2946 // mean level

predict prob in 2946
list prob female articles year yearsq select presthi in 2946, clean
drop in 2946

/*  Or you can use _N to determine the observation to be added.
local newn = _N + 1
set obs `newn'
replace female   =  1    if _n==`newn'
replace articles =  0    if _n==`newn'
replace year     =  7    if _n==`newn'
replace yearsq   = 49    if _n==`newn' // year squared
replace select   =  4.97 if _n==`newn' // mean level
replace presthi  =  0.05 if _n==`newn' // mean level
predict prob in `newn'
list prob female year yearsq select articles presthi in `newn', clean
drop in `newn'
*/

// #4
//  preditions with prvalue

prvalue, x(female=1 articles=0 year=7 yearsq=49) rest(mean)

// #5
//  multiple prvalues

prvalue, x(articles= 0 female=1 year=7 yearsq=49) brief
prvalue, x(articles= 5 female=1 year=7 yearsq=49) brief
prvalue, x(articles=10 female=1 year=7 yearsq=49) brief
prvalue, x(articles=15 female=1 year=7 yearsq=49) brief
prvalue, x(articles=20 female=1 year=7 yearsq=49) brief

// #6
//  returns from prvalue

prvalue, x(articles= 0 female=1 year=7 yearsq=49) brief
return list

// #7
//  automating multiple prvalue's use forvalues

gen plotx = .
  label var plotx "Number of articles"
gen plotp1 = .
  label var plotp1 "Pr(tenure for women in year 7)"
local i = 0
forvalues artval = 0(5)40 {
    local ++i
    prvalue, x(female=1 year=7 yearsq=49 articles=`artval')
    replace plotx = `artval' if _n==`i'
    replace plotp1 = r(p1) if _n==`i'
}
graph twoway connected plotp1 plotx, ///
     xlabel(0(5)40) ylabel(0(.1).8) ///
     ytitle("Pr(tenure for women in year 7)")
graph export sug_groups_07_prvalue.emf, replace

* add marker
gen plotp2 = plotp1 if _n==1

twoway (connected plotp1 plotx) ///
       (connected plotp2 plotx, msym(Sh) msize(large) clcol(red)) ///
     , xlabel(0(5)40) ylabel(0(.1).6) legend(off) ///
     ytitle("Pr(tenure for women in year 7)")
graph export sug_groups_07_prvalue_eg.emf, replace

// #8
//  plot with prgen

prgen articles, from(0) to(40) gap(5) generate(m0) x(female=1 year=7 yearsq=49)
label var m0p1 "Pr(tenure for women in year 7)"
label var m0x  "Number of articles"

graph twoway connected m0p1 m0x, ///
     xlabel(0(5)40) ylabel(0(.1).6)
graph export sug_groups_08_prgen.emf, replace
drop m0*

// #9
//  plot with logit using ci's

prgen articles, from(0) to(40) gap(5) generate(m0) ///
    x(female=1 year=7 yearsq=49) ci
label var m0p1 "Pr(tenure for women in year 7)"
label var m0x  "Number of articles"
label var m0p1lb "95% lower bound"
label var m0p1ub "95% upper bound"

* using shading to show confidence interval
graph twoway ///
    (rarea m0p1lb m0p1ub m0x, color(gs14)) ///
    (connected m0p1 m0x, msymbol(i)), ///
     subtitle("Probability of tenure with 95% confidence interval") ///
     yscale(range(0 .6)) ytitle("Pr(tenure for women in year 7)") ///
     xlabel(0(5)40) legend(off)

graph export sug_groups_09_prgenci.emf, replace

// #10
//  change in predictions

prchange female, x(year=7 yearsq=47)

quietly prvalue, x(female=0 year=7 yearsq=47) save
        prvalue, x(female=1 year=7 yearsq=47) dif

//
//  single predictor of tenure
//

// #11
//  variables in articles only model

spex tenure01, clear
sort male
by male: sum tenure female f_articles male m_articles

by male: vardesc tenure year yearsq select articles presthi

// #12
//  estimate logit model by group

logit tenure female f_articles male m_articles, nolog nocon

// #13
//  incorrect chow test

listcoef f_articles m_articles
test f_articles = m_articles

// #14
//  compute predictions at different #s of articles

* predictions with 5 articles
quietly prvalue, x(fem=1 f_art=5 male=0 m_art=0) save
        prvalue, x(fem=0 f_art=0 male=1 m_art=5) dif

matrix drop _all
foreach art of numlist 0(2)50 {
    quietly prvalue, x(fem=1 f_art=`art' male=0 m_art=0)
    matrix y_fem = nullmat(y_fem) \ pepred[2,2]
    quietly prvalue, x(fem=0 f_art=0 male=1 m_art=`art')
    matrix y_mal = nullmat(y_mal) \ pepred[2,2]
    matrix x_art = nullmat(x_art) \ `art'
}
svmat x_art
label var x_art1 "Number of Articles"
svmat y_fem
label var y_fem1 "Women" // Pr(for women)
svmat y_mal
label var y_mal1 "Men"   // Pr(for men)

// #15
//  plot predictions at different #s of articles

twoway (connected y_fem x_art, msym(i) clcol(red)) ///
       (connected y_mal x_art, msym(i) clcol(blue)) ///
      , ytitle(Pr(tenure)) xlabel(0(10)50) ///
        ylabel(0(.2)1.) yline(0 1) legend(pos(11) ring(0) cols(1))
graph export sug_groups_15_artprob.emf, replace

* cleanup
matrix drop _all
capture drop x_art1 y_fem1 y_mal1

// #16
//  compute discrete change and ci (sb_groups_m1ci.log)

foreach art of numlist 0(2)50 {
    quietly prvalue, save /// for women
       x(fem=1 f_art=`art' male=0 m_art=0)
    quietly prvalue, diff /// for men
       x(fem=0 f_art=0 male=1 m_art=`art')
    matrix y_mal = nullmat(y_mal) \ pepred[2,2]
    matrix y_fem = nullmat(y_fem) \ pepred[4,2]
    matrix y_dc  = nullmat(y_dc)  \ pepred[6,2]
    matrix y_ub  = nullmat(y_ub)  \ peupper[6,2]
    matrix y_lb  = nullmat(y_lb)  \ pelower[6,2]
    matrix x_art = nullmat(x_art) \ `art'
}

// #17
//  plot discrete change and ci (sb_groups_m1ciplot.log)

foreach v in x_art y_dc y_ub y_lb y_fem y_mal {
    svmat `v'
}
label var x_art "Number of Articles"
label var y_fem "Women" // Pr(for women)
label var y_mal "Men"   // Pr(for men)
label var y_dc  "Difference" // Pr(for men) - Pr(for women)
label var y_ub  "95% confidence interval"
label var y_lb  "95% confidence interval"

twoway ///
    (connected y_dc x_art, msym(i) clcol(green)) ///
    (connected y_ub x_art, msym(i) clcol(brown) clpat(dash)) ///
    (connected y_lb x_art, msym(i) clcol(brown) clpat(dash)) ///
    , legend(pos(11) ring(0) cols(1) order(1 2)) ///
    ytitle("Pr(men) - Pr(women)") ylabel(-.1(.1).8) yline(0) ///
    xlabel(0(10)50)

graph export sug_groups_17_artdc.emf, replace

//
//  M2 articles and hiprestige as a predictor
//

// #18
//  logit estimated (sb_groups_m2logit.log)

matrix drop _all
spex tenure01, clear

logit tenure female f_art f_presthi ///
               male m_art m_presthi, nolog nocon

// #19
//  chow tests
    test f_art = m_art
    test f_presthi  = m_presthi

// #20
//  compute predictions at both levels of presthi (sb_groups_m2pred.log)

foreach art of numlist 0(2)50 {
    * lower prestige jobs
    quietly prvalue, save /// for women
        x(fem=1 f_art=`art' male=0 m_art=0     f_presthi=0 m_presthi=0)
    quietly prvalue, diff /// for men
        x(fem=0 f_art=0     male=1 m_art=`art' f_presthi=0 m_presthi=0)
    matrix xlo_art = nullmat(xlo_art) \ `art'        // articles
    matrix ylo_mal = nullmat(ylo_mal) \ pepred[2,2]  // pr men
    matrix ylo_fem = nullmat(ylo_fem) \ pepred[4,2]  // pr women
    matrix ylo_dc = nullmat(ylo_dc)   \ pepred[6,2]  // difference
    matrix ylo_ub = nullmat(ylo_ub)   \ peupper[6,2] // upper limit
    matrix ylo_lb = nullmat(ylo_lb)   \ pelower[6,2] // lower limit

    * high prestige jobs
    quietly prvalue, save /// for women
        x(fem=1 f_art=`art' male=0 m_art=0     f_presthi=1 m_presthi=0)
    quietly prvalue, diff /// for men
        x(fem=0 f_art=0     male=1 m_art=`art' f_presthi=0 m_presthi=1)
    matrix xhi_art  = nullmat(xhi_art) \ `art'
    matrix yhi_mal  = nullmat(yhi_mal) \ pepred[2,2]
    matrix yhi_fem  = nullmat(yhi_fem) \ pepred[4,2]
    matrix yhi_dc   = nullmat(yhi_dc)  \ pepred[6,2]
    matrix yhi_ub = nullmat(yhi_ub)    \ peupper[6,2]
    matrix yhi_lb = nullmat(yhi_lb)    \ pelower[6,2]
}

foreach v in xlo_art ylo_fem ylo_mal ylo_dc ylo_ub ylo_lb ///
             xhi_art yhi_fem yhi_mal yhi_dc yhi_ub yhi_lb {
    svmat `v'
}
label var xlo_art1 "Number of Articles"
label var ylo_fem1 "Women - not distinguished"
label var ylo_mal1 "Men - not distinguished"
label var yhi_fem1 "Women - distinguished"
label var yhi_mal1 "Men - distinguished"
label var ylo_dc1  "Male-Female - not distinguished"
label var ylo_ub1  "95% confidence interval"
label var ylo_lb1  "95% confidence interval"
label var yhi_dc1  "Male-Female - distinguished"
label var yhi_lb1  "95% confidence interval"
label var yhi_ub1  "95% confidence interval"

// #21
//  four probability curves (sb_groups_m2plotprob.log)

twoway ///
    (connected yhi_fem xhi_art, msym(i) clcol(red) clpat(solid)) ///
    (connected ylo_fem xlo_art, msym(i) clcol(red) clpat(shortdash_dot)) ///
    (connected yhi_mal xhi_art, msym(i) clcol(blue) clpat(solid)) ///
    (connected ylo_mal xlo_art, msym(i) clcol(blue) clpat(shortdash_dot)) ///
    , legend(pos(11) order(2 1 4 3) ring(0) cols(1) region(ls(none)))  ///
      ylabel(0(.2)1.) yline(0 1) ytitle(Pr(tenure)) xlabel(0(10)50)
graph export sug_groups_21_prprsthi.emf, replace

// #22
//  two DC curves (sb_groups_m2dcplot.log)

twoway ///
    (connected ylo_dc xlo_art, msym(i) clcol(green))             ///
    (connected yhi_dc xhi_art, msym(i) clcol(orange))            ///
    , legend(pos(11) ring(0) cols(1) order(1 2)) xlabel(0(10)50) ///
      ytitle("Pr(men) - Pr(women)") ylabel(-.1(.2).9) yline(0)
graph export sug_groups_22_dcprst.emf, replace

// #23
//  plot discrete change with ci's (sb_groups_m2dcplotsig.log)

* case as missing if not significant
gen ylo_sigdc = ylo_dc if ylo_lb>=0 & ylo_lb!=.
gen yhi_sigdc = yhi_dc if yhi_lb>=0 & yhi_lb!=.
label var ylo_sigdc "Not distinguished"
label var ylo_dc  "if not significant"
label var yhi_sigdc "Distinguished"
label var yhi_dc  "if not significant"

twoway ///
    (connected ylo_sigdc xlo_art, clpat(solid) msym(i) clcol(green) ) ///
    (connected yhi_sigdc xhi_art, clpat(solid) msym(i) clcol(orange)) ///
    (connected ylo_dc    xlo_art, clpat(dash)  msym(i) clcol(green) ) ///
    (connected yhi_dc    xhi_art, clpat(dash)  msym(i) clcol(orange)) ///
    , legend(pos(11) order(2 4 1 3) ring(0) cols(1) region(ls(none))) ///
      ytitl("Pr(men) - Pr(women)") xlab(0(10)50) ylab(-.1(.2).9) ylin(0)
graph export sug_groups_23_dcsigprsthi.emf, replace

//
//  full model
//

// #25
//  estimate full model (sb_groups_m3logit.log)

matrix drop _all
spex tenure01, clear

logit tenure male m_year m_yearsq m_select m_articles m_prestige ///
             fem  f_year f_yearsq f_select f_articles f_prestige, nolog nocon

// #26
//  compute means for rhs variable (sb_groups_m3means.log)

foreach v in year yearsq select art prestige {
    quietly sum `v'
    local mn_`v' = r(mean)
}
local mn_yr = 7 // year for predictions
local mn_yrsq = `mn_yr' * `mn_yr' // year squared
local f_at_mn "f_year=`mn_yr' f_yearsq=`mn_yrsq' f_select=`mn_select'"
local m_at_mn "m_year=`mn_yr' m_yearsq=`mn_yrsq' m_select=`mn_select'"
local m_at_0  "mal=0 m_art=0 m_year=0 m_yearsq=0 m_select=0 m_prestige=0"
local f_at_0  "fem=0 f_art=0 f_year=0 f_yearsq=0 f_select=0 f_prestige=0"


// #27
//  compute predictions at year 7

foreach p in 1 2 3 4 5 { // prestige level

    di "prestige at `p'"

    foreach art of numlist 0(2)50 { // loop through articles

        quietly prvalue, save ///
          x(fem=1 f_art=`art' f_prestige=`p' `f_at_mn' `m_at_0')
        quietly prvalue, diff ///
          x(mal=1 m_art=`art' m_prestige=`p' `m_at_mn' `f_at_0')

        matrix x`p'_art = nullmat(x`p'_art) \ `art'
        matrix y`p'_mal = nullmat(y`p'_mal) \ pepred[2,2]
        matrix y`p'_fem = nullmat(y`p'_fem) \ pepred[4,2]
        matrix y`p'_dc  = nullmat(y`p'_dc)  \ pepred[6,2]
        matrix y`p'_ub  = nullmat(y`p'_ub)  \ peupper[6,2]
        matrix y`p'_lb  = nullmat(y`p'_lb)  \ pelower[6,2]
    }

    foreach v in x`p'_art y`p'_dc y`p'_ub y`p'_lb y`p'_fem y`p'_mal {
        svmat `v'
    }
    label var x`p'_art1 "Number of Articles"
    label var y`p'_mal1 "Men"
    label var y`p'_fem1 "Women"
    label var y`p'_dc1  "Male-Female difference"
    label var y`p'_ub1  "95% confidence interval"
    label var y`p'_lb1  "95% confidence interval"

} // end of prediction loop

// #28
//  probability plots by level of prestige

foreach p in 1 2 3 4 5 { // prestige level

    twoway ///
      (connected y`p'_fem x`p'_art, msym(i) clcol(red)) ///
      (connected y`p'_mal x`p'_art, msym(i) clcol(blue)) ///
     , subtitle("Plotted at prestige = `p'",pos(11)) ///
       legend(pos(11) order(2 1) ring(0) cols(1) region(ls(none))) ///
       ytitle(Pr(tenure)) xlabel(0(10)50) ///
       ylabel(0(.2)1.) yline(0 1)
    graph export sug_groups_28_probprst`p'.emf, replace

    * discrete change in probability
    twoway ///
      (connected y`p'_dc x`p'_art, msym(i) clcol(orange)) ///
      (connected y`p'_ub x`p'_art, msym(i) clcol(brown) clpat(dash)) ///
      (connected y`p'_lb x`p'_art, msym(i) clcol(brown) clpat(dash)) ///
     , subtitle("Plotted at prestige = `p'",pos(11)) ///
       legend(pos(11) order(2 1) ring(0) cols(1) region(ls(none))) ///
       ytitle("Pr(men) - Pr(women)") xlabel(0(10)50) ///
       ylabel(-.1(.2).7) yline(0)
    graph export sug_groups_28_dcprst`p'.emf, replace

} // loop over prestige

// #29
//  plot of multiple prestige levels

* case as missing if not significant
foreach p in 1 2 3 4 5 {
     gen y`p'_sigdc = y`p'_dc if y`p'_lb>=0 & y`p'_lb!=.
}

label var y1_sigdc "Weak (prestige=1)"
label var y1_dc  "  (not significant)"
label var y2_sigdc "Adequate (prestige=2)"
label var y3_sigdc "Good (prestige=3)"
label var y4_sigdc "Strong (prestige=4)"
label var y5_sigdc "Distinguished (prestige=5)"
label var y5_dc  "  (not significant)"

twoway ///
  (connected y1_sigdc x1_art, clpat(solid) msym(i) clcol(red))    ///
  (connected y2_sigdc x2_art, clpat(solid) msym(i) clcol(orange)) ///
  (connected y3_sigdc x3_art, clpat(solid) msym(i) clcol(green))  ///
  (connected y4_sigdc x4_art, clpat(solid) msym(i) clcol(blue))   ///
  (connected y5_sigdc x5_art, clpat(solid) msym(i) clcol(purple)) ///
  (connected y1_dc    x1_art, clpat(dash)  msym(i) clcol(red))    ///
  (connected y2_dc    x2_art, clpat(dash)  msym(i) clcol(orange)) ///
  (connected y3_dc    x3_art, clpat(dash)  msym(i) clcol(green))  ///
  (connected y4_dc    x4_art, clpat(dash)  msym(i) clcol(blue))   ///
  (connected y5_dc    x5_art, clpat(dash)  msym(i) clcol(purple)) ///
 , legend(pos(11) order(1 2 3 4 5) ring(0) cols(1) region(ls(none))) ///
   ytitle("Pr(men) - Pr(women)") xlab(0(10)50) ylab(.0(.1).5) ylin(0)
graph export sug_groups_29_dcprst.emf, replace

// #30
//  compute predictions at year 7 for various levels of publications

local a = 0
foreach art in 0 10 20 30 40 50 { // articles

    local ++a // to number variables that are created
    di "articles = `art'"

    foreach p of numlist 1(.2)5 { // loop through prestige

        quietly prvalue, save ///
          x(fem=1 f_art=`art' f_prestige=`p' `f_at_mn' `m_at_0')
        quietly prvalue, diff ///
          x(mal=1 m_art=`art' m_prestige=`p' `m_at_mn' `f_at_0')

        matrix x`a'a_prst = nullmat(x`a'a_prst) \ `p'
        matrix y`a'a_mal  = nullmat(y`a'a_mal)  \ pepred[2,2]
        matrix y`a'a_fem  = nullmat(y`a'a_fem)  \ pepred[4,2]
        matrix y`a'a_dc   = nullmat(y`a'a_dc)   \ pepred[6,2]
        matrix y`a'a_ub   = nullmat(y`a'a_ub)   \ peupper[6,2]
        matrix y`a'a_lb   = nullmat(y`a'a_lb)   \ pelower[6,2]

    }  // prestige level

} // end of articles loop

// #31
//  plot of multiple article levels

* create variables to plot
foreach a in 1 2 3 4 5 6 {
    foreach v in x`a'a_prst y`a'a_dc y`a'a_ub y`a'a_lb y`a'a_fem y`a'a_mal {
        svmat `v'
    }
    label var x`a'a_prst "Job Prestige"
    label var y`a'a_mal  "Men"
    label var y`a'a_fem  "Women"
    label var y`a'a_dc   "Male-Female difference"
    label var y`a'a_ub   "95% confidence interval"
    label var y`a'a_lb   "95% confidence interval"

}

* create new variable to indicate significance
foreach a in 1 2 3 4 5 6 {
     gen y`a'a_sigdc = y`a'a_dc if y`a'a_lb>=0 & y`a'a_lb!=.
}

* creaet plot
label var y1a_sigdc "no articles"
label var y1a_dc    "  (not significant)"
label var y2a_sigdc "10 articles"
label var y3a_sigdc "20 articles"
label var y4a_sigdc "30 articles"
label var y5a_sigdc "40 articles"
label var y6a_sigdc "50 articles"
label var y6a_dc    "  (not significant)"
*/
twoway ///
  (connected y1a_sigdc x1a_prst, clpat(solid) msym(i) clcol(red))    ///
  (connected y2a_sigdc x2a_prst, clpat(solid) msym(i) clcol(orange)) ///
  (connected y3a_sigdc x3a_prst, clpat(solid) msym(i) clcol(brown))  ///
  (connected y4a_sigdc x4a_prst, clpat(solid) msym(i) clcol(green))  ///
  (connected y5a_sigdc x5a_prst, clpat(solid) msym(i) clcol(blue))   ///
  (connected y6a_sigdc x6a_prst, clpat(solid) msym(i) clcol(purple)) ///
  (connected y1a_dc    x1a_prst, clpat(dash)  msym(i) clcol(red))    ///
  (connected y2a_dc    x2a_prst, clpat(dash)  msym(i) clcol(orange)) ///
  (connected y3a_dc    x3a_prst, clpat(dash)  msym(i) clcol(brown))  ///
  (connected y4a_dc    x4a_prst, clpat(dash)  msym(i) clcol(green))  ///
  (connected y5a_dc    x5a_prst, clpat(dash)  msym(i) clcol(blue))   ///
  (connected y6a_dc    x6a_prst, clpat(dash)  msym(i) clcol(purple)) ///
 , legend(pos(11) order(1 2 3 4 5 6) ring(0) cols(2) region(ls(none))) ///
   ytitle("Pr(men) - Pr(women)") xlab(1(1)5) ylab(.0(.1).5) ylin(0)

graph export sug_groups_31_dcart.emf, replace

/*
// #0
//  installing spost

net from http://www.indiana.edu/~jslsoc/stata/

//  or:

findit spost9
*/
*/

log close