capture log close
clear
set memory 100m
set more off
estimates clear
matrix drop _all

log using sug_groups_sim, replace
version 9

//  sug_groups_sim.do - stata user's group - july 2006
//  simulation data illustrating invariance of predictions
//  Scott Long - 20Jul2006

//  create artifical data
clear
set seed 11020
set obs 5000
local half = _N/2
local half1 = `half'+1
* indep variable is uniform
gen x = uniform()
    label var x "articles"
* divide into two groups
gen f = _n < `half1'
    label var f "is female?"
replace f = 0 if _n > `half'
gen m = 1 - f 
    label var m "is male?"
* group interactions
gen fx = f*x
    label var fx "female*articles"
gen mx = m*x
    label var fx "female*articles"

//  specify parameters - same slopes for men and women, different variances
local alpha = -1
local beta = 2
* sd for men and women differ
local f_sd = 2 // twice as large for women
local m_sd = 1

//  generate errors, ystar and y
gen double e = invnorm(uniform())
    label var e "normal errors"
* rescale errors for each group
replace e = e*`f_sd' if f==1
replace e = e*`m_sd' if f==0
* ystar
gen ystar = `alpha' + `beta'*x + e
    label var ystar "latent y*"
* observed binary y
gen y = ystar>0
    label var y "observed y"

//  estimates based on female only sample
probit y x if f==1 , nolog
    predict f_p
        label var f_p "Pr(y) females only model"
    local f_ll = e(ll)
    estimates store f_only

//  estimates based on male only sample
probit y x if f==0 , nolog
    predict m_p
        label var m_p "Pr(y) males only model"
    local m_ll = e(ll)
    estimates store m_only

//  compare predictions for men and women
sort x
twoway line f_p m_p x

//  estimate model using interactions
probit y m mx f fx , nolog nocon
    mat b = e(b)
    predict all_p
        label var all_p "Pr(y) combined male/female model"
    local all_ll = e(ll)
    local m_f_ll = `f_ll' + `m_ll'
    estimates store mf_inter

//  compare variables for men and women
sum all_p f_p y ystar x e  if f==1
sum all_p m_p y ystar x e  if f==0

//  compare pr(y) from two models
* predictions for women
sum f_p all_p if f==1
pwcorr f_p all_p if f==1
* predictions for men
sum f_p all_p if f==0
pwcorr m_p all_p if f==0

//  compare model results
di "LL in interaction model: `all_ll'"
di "LL-male + LL-female:     `m_f_ll'"

estimates table _all, stats(N ll) eform b(%9.3f) t(%6.2f)

log close