*! version 0.0.2 2011-02-28 jx * add probit and cloglog * version 0.0.1 2009-12-28 jx // Compute averaged marginal effects between groups in binary regression models // TO DO LIST: capture program drop grdiame program define grdiame, rclass version 10.0 tempname z /// b /// I /// xvec /// xb /// sxb /// prob pdf dfdl grad /// g1gradsum g0gradsum avggraddif /// marg avgmarg g1margsum g0margsum avgmargdif avgmargdifse avgmargdiflo avgmargdifhi /// V tobase tempvar touse tempfile syntax [if] [in], /// [x(passthru) rest(passthru) level(integer $S_level) all /// noBAse /// Group(string)] // // MODELS APPLIED if "`e(cmd)'"=="logistic" | /// "`e(cmd)'"=="logit" | /// "`e(cmd)'"=="probit" | /// "`e(cmd)'"=="cloglog" { local io = "typical binary" } if "`io'"=="" { di in r "this grcompare command does not work for the last type of model estimated." exit } matrix `b' = e(b) // command from SPost, get base values _pebase `if' `in', `x' `rest' `all' mat `tobase' = r(pebase) loc rhsnms "`r(rhsnms)'" loc nrhs `r(nrhs)' mat `I' = I(`nrhs'+1) sca `z' = -invnorm((1-`level'/100)/2) mark `touse' if e(sample) qui count loc totobs = r(N) // decompose groupvar into three gettoken grpvar group: group, parse(" ") // group variable name gettoken grp0 group: group, parse(" ") // group 0 group value gettoken grp1: group, parse(" ") // group 1 group value if "`grp0'"=="" | "`grp1'"=="" { di in r "You need to specify values for the two comparisons in the group opion!" exit } // group 1 group value: group 1 - group 0 /* di in y "grpvar=`grpvar'" di in g "grp0=`grp0'" di in r "grp1=`grp1'" */ // command from SPost for group size counts _pecats `grpvar' loc numcats = r(numcats) loc catvals "`r(catvals)'" forvalues i = 1 / `numcats' { loc crtval: word `i' of `catvals' tempname catcnt_`crtval' qui count if `grpvar'==`crtval' & `touse' // counts of each category for averaging late sca `catcnt_`crtval'' = r(N) } // establish the null matrices for marginal sum, gradient sum for two groups mat `g1margsum' = J(1, `nrhs'+1, 0) mat `g1gradsum' = J(`nrhs'+1, `nrhs'+1, 0) mat `g0margsum' = J(1, `nrhs'+1, 0) mat `g0gradsum' = J(`nrhs'+1, `nrhs'+1, 0) forvalues k= 1/`totobs' { // check if an observation is in the estimation sample or not loc sampchk = `touse'[`k'] if `sampchk'==0 { continue } // otherwise, take that into computation else { cap mat drop `xvec' forvalues m = 1 / `nrhs'{ tempname varval`m' loc crtvar: word `m' of `rhsnms' sca `varval`m'' = `crtvar'[`k'] mat `xvec' = nullmat(`xvec'), `varval`m'' } // add constant one mat `xvec' = `xvec',J(1,1,1) *mat list `xvec' *mat list `b' // get matrix linear predictor mat `xb' = `xvec' * `b'' *local nb = colsof(`b') - 1 *matrix `b' = `b'[1,1..`nb'] /* get rid of _con */ // construct marginal effect and gradient for each observation sca `sxb' = `xb'[1,1] if "`e(cmd)'"=="logit" | "`e(cmd)'"=="logistic" { sca `prob' = exp(`sxb')/(1+exp(`sxb')) sca `pdf' = `prob'*(1-`prob') sca `dfdl' = (1-2*`prob')*`prob'*(1-`prob') } else if "`e(cmd)'"=="probit" { sca `prob' = normprob(`sxb') sca `pdf' = exp(-`sxb'*`sxb'/2)/sqrt(2*_pi) sca `dfdl' = -`sxb'*`pdf' } else if "`e(cmd)'"=="cloglog" { sca `prob' = 1- exp(-exp(`sxb')) sca `pdf' = (`prob'-1)*ln(1-`prob') sca `dfdl'= `pdf'*(ln(1-`prob')+1) } mat `marg' = `pdf'*`b' mat `grad' = `pdf'*`I'+`dfdl'*(`b')'*(`xvec') // grab the group value from the groupvar loc yesno= `grpvar'[`k'] /* di in r "yesno=`yesno' `grpvar'[`k']" mat list `marg' mat list `grad' di in y "========================" */ //computation for sum of marginal effects and gradient for group 1 and 0 if "`yesno'" == "`grp1'" { mat `g1margsum' = `g1margsum' + `marg' mat `g1gradsum' = `g1gradsum' + `grad' } if "`yesno'" == "`grp0'" { mat `g0margsum' = `g0margsum' + `marg' mat `g0gradsum' = `g0gradsum' + `grad' } /* di in y "matrice marg and grad" mat list `marg' mat list `grad' */ /* di in y "matricex g1margsum and g0margsum" mat list `g1margsum' mat list `g0margsum' */ } } // sum and average mat `avgmargdif' = `g1margsum'/`catcnt_`grp1'' - `g0margsum'/`catcnt_`grp0'' mat `avggraddif' = `g1gradsum'/`catcnt_`grp1'' - `g0gradsum'/`catcnt_`grp0'' /* mat list `g1margsum' mat list `g0margsum' mat list `avggraddif' */ // variance of difference in average marginal effect mat `V' = (`avggraddif') * e(V) * (`avggraddif')' // get s.e. for marginal effect mat `avgmargdifse' = vecdiag(cholesky(diag(vecdiag(`V')))) mat `avgmargdiflo' = `avgmargdif' - `z'*`avgmargdifse' mat `avgmargdifhi' = `avgmargdif' + `z'*`avgmargdifse' mat `avgmargdif'= `avgmargdiflo' \ `avgmargdif' \ `avgmargdifhi' mat `avgmarg' = `g0margsum' / `catcnt_`grp0'' \ `g1margsum' / `catcnt_`grp1'' mat rownames `avgmarg' = group0 group1 mat rownames `avgmargdif' = diflo g1-g0 difhi mat list `avgmargdif', noheader ret mat avgmargdif = `avgmargdif' ret mat avgmarg = `avgmarg' end